/[PAMELA software]/DarthVader/CalorimeterLevel2/src/CaloCore.cpp
ViewVC logotype

Annotation of /DarthVader/CalorimeterLevel2/src/CaloCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.14 - (hide annotations) (download)
Tue Nov 28 10:20:27 2006 UTC (18 years, 3 months ago) by mocchiut
Branch: MAIN
CVS Tags: v2r01
Changes since 1.13: +3 -3 lines
Implemented CaloStrip::Set(Float_t,Float_t,Float_t) method, masked calo plane 18X

1 mocchiut 1.1 //
2     // Given a calibration and a data file this program create an ntuple with LEVEL2 calorimeter variables - Emiliano Mocchiutti
3     //
4 mocchiut 1.3 // CaloCore.cxx version 3.05 (2006-05-30)
5 mocchiut 1.1 //
6     // The only input needed is the path to the directory created by YODA for the data file you want to analyze.
7     //
8     // Changelog:
9 mocchiut 1.3 //
10 mocchiut 1.11 // 3.08 (2006-11-13): Added high energy nuclei capability and "process all events" capability.
11     //
12 mocchiut 1.4 // 3.04 - 3.05 (2006-05-30): Qlast and nlast are now calculated using 4 (not 8) strips aournd the shower axis. Small bug fixed.
13 mocchiut 1.2 //
14     // 3.03 - 3.04 (2006-05-23): Forgot to put impx and impy in the PAMELA reference system, fixed.
15 mocchiut 1.1 //
16     // 3.02 - 3.03 (2006-05-18): updated to be called in DarthVader. Output dimension are now in cm and in the PAMELA reference system.
17     //
18     // 3.01 - 3.02 (2006-04-21): when copying entries get size of caclone and not of ca... (shouldn't matter). Fixed increasing file dimension bug when reprocessing.
19     // Added variable planemax[2], plane of maximum energy release (x and y) in final output. Use ItoRunInfo instead of RunInfo.
20     //
21     // 3.00 - 3.01 (2006-04-14): fixed small bug in tagging the track used to determine track-related variables, put in caloprocessing the opening of parameters files, fixed
22     // small bug in fortran routines
23     //
24     // 2.01 - 3.00 (2006-04-14): almost everything has changed. Now it can process one, all or some runs, introduced the final CaloLevel2 class+methods and the
25     // working class "CaloProcessing", linked to the preliminary tracker flight software v0r00, reads YODA unique output files,
26     // reduced the number of installed libraries, F77 programs splitted depending on the function contained, introduced the readout and
27     // processing of self-trigger events, if the tracker provides more than one track all calorimeter track-related variables are saved
28     // as many times as the number of tracks (via the TClonesArray object in the level2 rootple) and many other small changes.
29     //
30     // 2.00 - 2.01 (2006-01-26): bug: wrong calculation of baselines in some cases, fixed.
31     //
32     // 1.00 - 2.00 (2006-01-11): use TSQL ROOT classes instead of directly calling MySQL.
33     //
34     // 0.00 - 1.00 (2005-09-14): seems working.
35     //
36     // 0.00 (2005-09-09): clone of CaloLEVEL2.c .
37     //
38     // C/C++ headers
39     //
40     #include <fstream>
41     #include <string.h>
42     //
43     // ROOT headers
44     //
45     #include <TTree.h>
46     #include <TClassEdit.h>
47     #include <TObject.h>
48     #include <TList.h>
49 mocchiut 1.8 #include <TArrayI.h>
50 mocchiut 1.1 #include <TSystem.h>
51     #include <TSystemDirectory.h>
52     #include <TString.h>
53     #include <TFile.h>
54     #include <TClass.h>
55     #include <TCanvas.h>
56     #include <TH1.h>
57     #include <TH1F.h>
58     #include <TH2D.h>
59     #include <TLatex.h>
60     #include <TPad.h>
61     #include <TSQLServer.h>
62     #include <TSQLRow.h>
63     #include <TSQLResult.h>
64     #include <TClonesArray.h>
65     #include <TStreamerInfo.h>
66     //
67     // YODA headers
68     //
69     #include <PamelaRun.h>
70     #include <physics/trigger/TriggerEvent.h>
71     //
72     // This program headers
73     //
74     #include <RunInfo.h>
75     #include <CaloCore.h>
76 mocchiut 1.13 #include <CaloLevel1.h>
77 mocchiut 1.1 #include <CaloLevel2.h>
78     #include <CaloProcessing.h>
79     #include <CaloVerl2.h>
80     //
81     // Tracker classes headers and definitions
82     //
83     #include <TrkLevel2.h>
84     //
85     using namespace std;
86     //
87     // CORE ROUTINE
88     //
89 mocchiut 1.8 int CaloCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t calargc, char *calargv[]){
90 mocchiut 1.1 //
91     // Set these to true to have a very verbose output.
92     //
93     Bool_t verbose = false;
94     Bool_t debug = false;
95     //
96 mocchiut 1.11 Bool_t trackanyway = true;
97     //
98     Float_t rigdefault = 50.;
99     //
100     Bool_t hZn = true;
101     //
102     Bool_t withtrk = true;
103     //
104     Bool_t st = true;
105     //
106 mocchiut 1.13 Bool_t getl1 = true;
107     //
108 mocchiut 1.1 // Output directory is the working directoy.
109     //
110     const char* outdir = gSystem->DirName(gSystem->DirName(file->GetPath()));
111     //
112     Int_t ri = 0;
113     TString processFolder = "calorimeterFolder";
114     if ( calargc > 0 ){
115     ri = 0;
116     while ( ri < calargc ){
117     if ( !strcmp(calargv[ri],"-processFolder") ) {
118     if ( calargc < ri+1 ){
119     throw -3;
120     };
121     processFolder = (TString)calargv[ri+1];
122     ri++;
123     };
124     if ( !strcmp(calargv[ri],"-v") || !strcmp(calargv[ri],"--verbose") ) {
125     verbose = true;
126     };
127     if ( !strcmp(calargv[ri],"-g") || !strcmp(calargv[ri],"--debug") ) {
128 mocchiut 1.13 verbose = true;
129 mocchiut 1.1 debug = true;
130     };
131 mocchiut 1.11 if ( !strcmp(calargv[ri],"--alltracks") ) {
132     trackanyway = true;
133     };
134     if ( !strcmp(calargv[ri],"--no-tracker") ) {
135     withtrk = false;
136     };
137     if ( !strcmp(calargv[ri],"--with-tracker") ) {
138     withtrk = true;
139     };
140     if ( !strcmp(calargv[ri],"--defrig") ) {
141     if ( calargc < ri+1 ){
142     throw -3;
143     };
144     rigdefault = atof(calargv[ri+1]);
145     ri++;
146     };
147     if ( !strcmp(calargv[ri],"--no-alltracks") ) {
148     trackanyway = false;
149     };
150     if ( !strcmp(calargv[ri],"--highZnuclei") ) {
151     hZn = true;
152     };
153     if ( !strcmp(calargv[ri],"--no-highZnuclei") ) {
154     hZn = false;
155     };
156     if ( !strcmp(calargv[ri],"--selftrigger") ) {
157     st = true;
158     };
159     if ( !strcmp(calargv[ri],"--no-selftrigger") ) {
160     st = false;
161     };
162 mocchiut 1.13 if ( !strcmp(calargv[ri],"--no-level1") ) {
163     getl1 = false;
164     };
165 mocchiut 1.11 if ( !strcmp(calargv[ri],"--help") ) {
166     printf("\n\n CALORIMETER HELP CALLED\n\n");
167     printf(" CaloCore options: \n");
168     printf(" -v | --verbose be verbose\n");
169     printf(" -g | --debug be really verbose\n");
170     printf(" --defrig rig rig is the default rigidity in GV to be used to\n");
171     printf(" obtain calorimeter variables in the routines\n");
172     printf(" \"alltracks\" and \"higZnuclei\" [default = 50]\n");
173     printf(" --alltracks fill the track related variables even in the case\n");
174     printf(" of no tracks from tracker and no selftrigger event\n");
175 mocchiut 1.14 printf(" when we have a calorimeter fit for both views [default]\n");
176     printf(" --no-alltracks fill the track related variables only in the case\n");
177     printf(" of a good track from tracker or selftrigger\n");
178 mocchiut 1.11 printf(" --highZnuclei call the routine to analyze high Z nuclei\n");
179     printf(" selftrigger events [default]\n");
180     printf(" --no-highZnuclei do not call the routine to analyze high Z nuclei\n");
181     printf(" selftrigger events\n");
182     printf(" --no-tracker do not use tracker level2\n");
183     printf(" --with-tracker use tracker level2 [default]\n");
184     printf(" --selftrigger process selftrigger events [default]\n");
185     printf(" --no-selftrigger skip selftrigger events\n");
186 mocchiut 1.13 printf(" --no-level1 do not save Level1 TBranch\n");
187 mocchiut 1.11 throw -114;
188     };
189 mocchiut 1.1 ri++;
190     };
191     };
192     //
193 mocchiut 1.13 if ( verbose ){
194     printf("\n");
195     if ( getl1 ) printf(" Saving calorimeter level1 data \n");
196     if ( st ) printf(" Calling selftrigger subroutine \n");
197     if ( hZn ) printf(" Calling high energy nuclei subroutine \n");
198     if ( trackanyway ) printf(" Filling track related variables for all the possible tracks \n");
199     if ( hZn || trackanyway ) printf(" => default assumed rigidity %f \n",rigdefault);
200     if ( withtrk ) printf(" Using tracker level2 data \n");
201     printf("\n");
202     };
203 mocchiut 1.11 //
204 mocchiut 1.1 // Working filename
205     //
206     TString outputfile;
207     stringstream name;
208     name.str("");
209     name << outdir << "/";
210     //
211     // Variables.
212     //
213     TTree *tracker = 0;
214     TTree *calo = 0;
215     TTree *caloclone = 0;
216     Bool_t reproc = false;
217     Bool_t reprocall = false;
218 mocchiut 1.8 UInt_t nevents = 0;
219 mocchiut 1.1 UInt_t nobefrun = 0;
220     UInt_t noaftrun = 0;
221     UInt_t numbofrun = 0;
222 mocchiut 1.8 UInt_t totnorun = 0;
223 mocchiut 1.1 //
224     Int_t code = 0;
225     Int_t sgnl;
226     //
227     // calorimeter level2 classes
228     //
229 mocchiut 1.13 CaloLevel1 *c1 = 0;
230     CaloLevel1 *c1clone = 0;
231     if ( getl1 ){
232     c1 = new CaloLevel1();
233     c1clone = new CaloLevel1();
234     };
235     //
236     // calorimeter level2 classes
237     //
238 mocchiut 1.1 CaloLevel2 *ca = new CaloLevel2();
239     CaloLevel2 *caclone = new CaloLevel2();
240     //
241     TrkLevel2 *trk = new TrkLevel2();
242 mocchiut 1.8 Int_t nevtrkl2 = 0;
243 mocchiut 1.1 //
244     UInt_t procev = 0;
245     //
246 mocchiut 1.8 // define variables where to store the absolute run header and run trailer times (unsigned long long integers, when set to a number use to store the correct number).
247 mocchiut 1.1 //
248 mocchiut 1.8 UInt_t runheadtime = 0;
249     UInt_t runtrailtime = 0;
250 mocchiut 1.1 UInt_t evfrom = 0;
251     UInt_t evto = 0;
252 mocchiut 1.8 UInt_t totfileentries = 0;
253     UInt_t idRun = 0;
254 mocchiut 1.1 Int_t id_reg_run=-1;
255     stringstream ftmpname;
256     TString fname;
257     //
258     // define variables for opening and reading level0 file
259     //
260     TFile *l0File = 0;
261     TTree *l0tr = 0;
262     TBranch *l0head = 0;
263     TBranch *l0calo = 0;
264     TBranch *l0trig = 0;
265 mocchiut 1.8 pamela::EventHeader *eh = 0;
266     pamela::PscuHeader *ph = 0;
267 mocchiut 1.1 pamela::trigger::TriggerEvent *trig = 0;
268     //
269     // Define some basic variables
270     //
271     CaloProcessing *event = new CaloProcessing(); // NOTICE: very important to call here the constructor!
272     stringstream file2;
273     stringstream file3;
274     stringstream qy;
275     // Bool_t imtrack = false;
276     Bool_t filled = false;
277     //
278 mocchiut 1.8 UInt_t caloevents = 0;
279 mocchiut 1.1 stringstream calfile;
280     stringstream aligfile;
281     //
282     Int_t i = -1;
283     Int_t itr = -1;
284     Int_t badevent = 0;
285     Int_t totevent = 0;
286     //
287 mocchiut 1.8 UInt_t atime = 0;
288 mocchiut 1.1 //
289     Int_t S3 = 0;
290     Int_t S2 = 0;
291     Int_t S12 = 0;
292     Int_t S11 = 0;
293     UInt_t re = 0;
294 mocchiut 1.9 UInt_t jumped = 0;
295 mocchiut 1.1 //
296     TString caloversion;
297     ItoRunInfo *runinfo = 0;
298 mocchiut 1.8 TArrayI *runlist = 0;
299 mocchiut 1.1 //
300     Float_t tmptrigty = -1.;
301     Int_t ntrkentry = 0;
302     GL_PARAM *q4 = new GL_PARAM();
303 mocchiut 1.8 UInt_t tttrkpar1 = 0;
304 mocchiut 1.1 Bool_t trkpar1 = true;
305     GL_ROOT *glroot = new GL_ROOT();
306 mocchiut 1.8 GL_TIMESYNC *dbtime = 0;
307 mocchiut 1.1 //
308     //
309     //
310     TFile *tempfile = 0;
311     TTree *tempcalo = 0;
312     stringstream tempname;
313     stringstream calofolder;
314     tempname.str("");
315     tempname << outdir;
316     tempname << "/" << processFolder.Data();
317     calofolder.str("");
318     calofolder << tempname.str().c_str();
319     gSystem->MakeDirectory(calofolder.str().c_str());
320     tempname << "/calotree_run";
321     tempname << run << ".root";
322     //
323     // As a first thing we must check what we have to do: if run = -1 we must process all events in the file has been passed
324     // if run != -1 we must process only that run but first we have to check if the branch calorimeter already exist in the file
325     // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
326     //
327 mocchiut 1.8 if ( run == 0 ) reproc = true;
328 mocchiut 1.1 //
329     //
330     //
331     if ( !file->IsOpen() ){
332     if ( verbose ) printf(" CALORIMETER - ERROR: cannot open file for writing\n");
333     throw -101;
334     };
335     //
336 mocchiut 1.11 if ( withtrk ){
337     //
338     // Does it contain the Tracker tree?
339     //
340     tracker = (TTree*)file->Get("Tracker");
341     if ( !tracker ) {
342     if ( verbose ) printf(" CALORIMETER - ERROR: no tracker tree\n");
343     code = -102;
344     goto closeandexit;
345     };
346     tracker->SetBranchAddress("TrkLevel2",&trk);
347     nevtrkl2 = tracker->GetEntries();
348 mocchiut 1.1 };
349     //
350     // Call runinfo
351     //
352     sgnl = 0;
353     runinfo = new ItoRunInfo(file);
354     //
355     // update versioning information and retrieve informations about the run to be processed
356     //
357     caloversion = CaloInfo(false);
358     sgnl = runinfo->Update(run,"CALO",caloversion);
359     //
360     if ( sgnl ){
361     if ( verbose ) printf(" CALORIMETER - ERROR: RunInfo exited with non-zero status\n");
362     code = sgnl;
363     goto closeandexit;
364     } else {
365     sgnl = 0;
366     };
367     //
368     // number of events in the file BEFORE the first event of our run
369     //
370     nobefrun = runinfo->GetFirstEntry();
371     //
372     // total number of events in the file
373     //
374     totfileentries = runinfo->GetFileEntries();
375     //
376     // first file entry AFTER the last event of our run
377     //
378     noaftrun = runinfo->GetLastEntry() + 1;
379     //
380     // number of run to be processed
381     //
382     numbofrun = runinfo->GetNoRun();
383     //
384     // number of runs in the file
385     //
386     totnorun = runinfo->GetRunEntries();
387     //
388     // Does it contain already a Calorimeter branch? if so we are reprocessing data, if not we must create it.
389     //
390     caloclone = (TTree*)file->Get("Calorimeter");
391     //
392     if ( !caloclone ){
393     reproc = false;
394 mocchiut 1.8 if ( run == 0 && verbose ) printf(" CALORIMETER - WARNING: you are reprocessing data but calorimeter tree does not exist!\n");
395     if ( runinfo->IsReprocessing() && run != 0 && verbose ) printf(" CALORIMETER - WARNING: it seems you are not reprocessing data but calorimeter\n versioning information already exists in RunInfo.\n");
396 mocchiut 1.1 //
397     } else {
398     //
399     reproc = true;
400     //
401     if ( verbose ) printf("\n Preparing the pre-processing...\n");
402     //
403 mocchiut 1.8 if ( run == 0 ){
404 mocchiut 1.1 //
405     // if we are reprocessing everything we don't need to copy any old event and we can just create a new branch in the clone tree and jump steps 4/5/7.
406     //
407     if ( verbose ) printf("\n CALORIMETER - WARNING: Reprocessing all runs in the file\n");
408     reprocall = true;
409     //
410     } else {
411     //
412     // we are reprocessing a single run
413     //
414 mocchiut 1.8 if ( verbose ) printf("\n CALORIMETER - WARNING: Reprocessing run number %u \n",run);
415 mocchiut 1.1 reprocall = false;
416     //
417     //
418     //
419     tempfile = new TFile(tempname.str().c_str(),"RECREATE");
420     tempcalo = caloclone->CloneTree(-1,"fast");
421     tempcalo->SetName("Calorimeter-old");
422     tempfile->Write();
423     tempfile->Close();
424     };
425     //
426     // delete old tree
427     //
428     caloclone->Delete("all");
429     //
430     if ( verbose ) printf("\n ...done!\n");
431     //
432     };
433     //
434     // create calorimeter tree calo
435     //
436     file->cd();
437     calo = new TTree("Calorimeter-new","PAMELA Level2 calorimeter data");
438     calo->Branch("CaloLevel2","CaloLevel2",&ca);
439 mocchiut 1.13 if ( getl1 ) calo->Branch("CaloLevel1","CaloLevel1",&c1);
440 mocchiut 1.1 //
441     if ( reproc && !reprocall ){
442     //
443     //
444     //
445     tempfile = new TFile(tempname.str().c_str(),"READ");
446     caloclone = (TTree*)tempfile->Get("Calorimeter-old");
447     caloclone->SetBranchAddress("CaloLevel2",&caclone);
448 mocchiut 1.13 TBranch *havel1 = caloclone->GetBranch("CaloLevel1");
449     if ( !havel1 && getl1 ) throw -117;
450     if ( havel1 && !getl1 ) throw -118;
451     if ( getl1 ) caloclone->SetBranchAddress("CaloLevel1",&c1clone);
452 mocchiut 1.1 //
453     if ( nobefrun > 0 ){
454     if ( verbose ) printf("\n Pre-processing: copying events from the old tree before the processed run\n");
455 mocchiut 1.8 if ( verbose ) printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
456 mocchiut 1.1 if ( verbose ) printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
457     for (UInt_t j = 0; j < nobefrun; j++){
458     //
459     caloclone->GetEntry(j);
460     //
461     // copy caclone to ca
462     //
463 mocchiut 1.13 if ( getl1 ) memcpy(&c1,&c1clone,sizeof(c1clone));
464 mocchiut 1.1 memcpy(&ca,&caclone,sizeof(caclone));
465     //
466     // Fill entry in the new tree
467     //
468     calo->Fill();
469     //
470 mocchiut 1.13 if ( getl1 ) c1->Clear();
471 mocchiut 1.6 ca->Clear();
472     //
473 mocchiut 1.1 };
474     if ( verbose ) printf(" Finished successful copying!\n");
475     };
476     };
477     //
478     // Get the list of run to be processed
479     //
480     runlist = runinfo->GetRunList();
481     //
482     // Loop over the run to be processed
483     //
484     for (UInt_t irun=0; irun < numbofrun; irun++){
485     //
486     badevent = 0;
487     //
488     idRun = runlist->At(irun);
489     if ( verbose ) printf("\n\n\n ####################################################################### \n");
490 mocchiut 1.8 if ( verbose ) printf(" PROCESSING RUN NUMBER %u \n",idRun);
491 mocchiut 1.1 if ( verbose ) printf(" ####################################################################### \n\n\n");
492     //
493     sgnl = runinfo->GetRunInfo(idRun);
494     if ( sgnl ){
495     if ( verbose ) printf(" CALORIMETER - ERROR: RunInfo exited with non-zero status\n");
496     code = sgnl;
497     goto closeandexit;
498     } else {
499     sgnl = 0;
500     };
501 mocchiut 1.8 id_reg_run = runinfo->ID_ROOT_L0;
502 mocchiut 1.1 runheadtime = runinfo->RUNHEADER_TIME;
503     runtrailtime = runinfo->RUNTRAILER_TIME;
504 mocchiut 1.8 evfrom = runinfo->EV_FROM;
505     evto = runinfo->EV_TO;
506 mocchiut 1.1 //
507     if ( id_reg_run == -1 ){
508 mocchiut 1.8 if ( verbose ) printf("\n CALORIMETER - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
509 mocchiut 1.1 code = -5;
510     goto closeandexit;
511     };
512     //
513 mocchiut 1.8 // prepare the timesync for the db
514     //
515 mocchiut 1.12 if ( !dbc->IsConnected() ) throw -116;
516 mocchiut 1.8 dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
517     //
518 mocchiut 1.1 // Search in the DB the path and name of the LEVEL0 file to be processed.
519     //
520 mocchiut 1.12 if ( !dbc->IsConnected() ) throw -116;
521 mocchiut 1.8 glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
522 mocchiut 1.1 //
523     ftmpname.str("");
524     ftmpname << glroot->PATH.Data() << "/";
525     ftmpname << glroot->NAME.Data();
526     fname = ftmpname.str().c_str();
527     //
528     // print out informations
529     //
530 mocchiut 1.8 totevent = runinfo->NEVENTS;
531 mocchiut 1.1 if ( verbose ) printf("\n LEVEL0 data file: %s \n",fname.Data());
532 mocchiut 1.8 if ( verbose ) printf(" RUN HEADER absolute time is: %u \n",runheadtime);
533     if ( verbose ) printf(" RUN TRAILER absolute time is: %u \n",runtrailtime);
534     if ( verbose ) printf(" %i events to be processed for run %u: from %i to %i (reg entries)\n\n",totevent,idRun,evfrom,evfrom+totevent);
535 mocchiut 1.1 //
536     // Open Level0 file
537     //
538     l0File = new TFile(fname.Data());
539     if ( !l0File ) {
540     if ( verbose ) printf(" CALORIMETER - ERROR: problems opening Level0 file\n");
541     code = -6;
542     goto closeandexit;
543     };
544     l0tr = (TTree*)l0File->Get("Physics");
545     if ( !l0tr ) {
546     if ( verbose ) printf(" CALORIMETER - ERROR: no Physics tree in Level0 file\n");
547     l0File->Close();
548     code = -7;
549     goto closeandexit;
550     };
551     l0head = l0tr->GetBranch("Header");
552     if ( !l0head ) {
553     if ( verbose ) printf(" CALORIMETER - ERROR: no Header branch in Level0 tree\n");
554     l0File->Close();
555     code = -8;
556     goto closeandexit;
557     };
558     l0calo = l0tr->GetBranch("Calorimeter");
559     if ( !l0calo ) {
560     if ( verbose ) printf(" CALORIMETER - ERROR: no Calorimeter branch in Level0 tree\n");
561     l0File->Close();
562     code = -103;
563     goto closeandexit;
564     };
565     l0trig = l0tr->GetBranch("Trigger");
566     if ( !l0trig ) {
567     if ( verbose ) printf(" CALORIMETER - ERROR: no Trigger branch in Level0 tree\n");
568     l0File->Close();
569     code = -104;
570     goto closeandexit;
571     };
572     //
573     l0tr->SetBranchAddress("Trigger", &trig);
574 mocchiut 1.8 l0tr->SetBranchAddress("Header", &eh);
575 mocchiut 1.1 //
576     // Construct the event object, look for the calibration which include the first header
577     //
578     sgnl = 0;
579     if ( verbose ) printf(" Check for calorimeter calibrations and initialize event object \n");
580 mocchiut 1.12 if ( !dbc->IsConnected() ) throw -116;
581 mocchiut 1.1 event->ProcessingInit(dbc,runheadtime,sgnl,l0tr,debug,verbose);
582     if ( verbose ) printf("\n");
583     if ( sgnl == 100 ) {
584     code = sgnl;
585     if ( verbose ) printf(" CALORIMETER - WARNING: run header not included in any calibration interval\n");
586     sgnl = 0;
587     };
588     if ( sgnl ){
589     l0File->Close();
590     code = sgnl;
591     goto closeandexit;
592     };
593     //
594     qy.str("");
595     //
596 mocchiut 1.8 nevents = l0head->GetEntries();
597 mocchiut 1.1 caloevents = l0calo->GetEntries();
598     //
599     if ( nevents < 1 ) {
600     if ( verbose ) printf(" CALORIMETER - ERROR: Level0 file is empty\n\n");
601     l0File->Close();
602     code = -11;
603     goto closeandexit;
604     };
605     //
606     if ( evto > nevents-1 ) {
607     if ( verbose ) printf(" CALORIMETER - ERROR: too few entries in the registry tree\n");
608     l0File->Close();
609     code = -12;
610     goto closeandexit;
611     };
612     //
613     // Check if we have to load parameter files
614     //
615     sgnl = 0;
616 mocchiut 1.12 if ( !dbc->IsConnected() ) throw -116;
617 mocchiut 1.1 sgnl = event->ChkParam(dbc,runheadtime); // calorimeter parameter files
618     if ( sgnl < 0 ){
619     code = sgnl;
620     l0File->Close();
621     goto closeandexit;
622     };
623     //
624 mocchiut 1.11 if ( withtrk ){
625     if ( trkpar1 || ( tttrkpar1 != 0 && tttrkpar1 < runheadtime ) ){
626     trkpar1 = false;
627 mocchiut 1.12 if ( !dbc->IsConnected() ) throw -116;
628 mocchiut 1.11 Int_t glpar = q4->Query_GL_PARAM(runinfo->RUNHEADER_TIME,1,dbc);
629     if ( glpar < 0 ){
630     code = glpar;
631     goto closeandexit;
632     };
633     tttrkpar1 = q4->TO_TIME;
634     // ----------------------------
635     // Read the magnetic field
636     // ----------------------------
637     if ( verbose ) printf(" Reading magnetic field maps at %s\n",(q4->PATH+q4->NAME).Data());
638     trk->LoadField(q4->PATH+q4->NAME);
639     if ( verbose ) printf("\n");
640 mocchiut 1.4 };
641 mocchiut 1.1 };
642     //
643     // run over all the events of the run
644     //
645     if ( verbose ) printf("\n Ready to start! \n\n Processed events: \n\n");
646     //
647 mocchiut 1.8 //
648     for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
649 mocchiut 1.1 //
650     if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);
651     //
652 mocchiut 1.8 l0head->GetEntry(re);
653 mocchiut 1.1 //
654     // absolute time of this event
655     //
656 mocchiut 1.8 ph = eh->GetPscuHeader();
657     atime = dbtime->DBabsTime(ph->GetOrbitalTime());
658 mocchiut 1.1 //
659     //
660     //
661 mocchiut 1.8 if ( re > caloevents-1 ){
662 mocchiut 1.1 if ( verbose ) printf(" CALORIMETER - ERROR: no physics events with entry = %i in Level0 file\n",i);
663     l0File->Close();
664     code = -112;
665     goto closeandexit;
666     };
667     //
668     if ( atime > runtrailtime || atime < runheadtime ) {
669     if ( verbose ) printf(" CALORIMETER - WARNING: event at time outside the run time window, skipping it\n");
670 mocchiut 1.9 jumped++;
671 mocchiut 1.1 goto jumpev;
672     };
673     //
674     // retrieve tracker informations
675     //
676     if ( !reprocall ){
677 mocchiut 1.9 itr = nobefrun + (re - evfrom - jumped);
678 mocchiut 1.1 } else {
679 mocchiut 1.9 itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
680 mocchiut 1.1 };
681 mocchiut 1.11 //
682     if ( withtrk ){
683     if ( itr > nevtrkl2 ){
684     if ( verbose ) printf(" CALORIMETER - ERROR: no tracker events with entry = %i in Level2 file\n",itr);
685     l0File->Close();
686     code = -113;
687     goto closeandexit;
688     };
689     //
690     trk->Clear();
691     //
692     tracker->GetEntry(itr);
693     //
694 mocchiut 1.1 };
695 mocchiut 1.10 //
696 mocchiut 1.1 procev++;
697     //
698     // start processing
699     //
700 mocchiut 1.13 if ( getl1 ) c1->Clear();
701 mocchiut 1.6 ca->Clear();
702 mocchiut 1.1 //
703     // determine who generate the trigger for this event (TOF, S4/PULSER, CALO)
704     //
705 mocchiut 1.8 l0trig->GetEntry(re);
706 mocchiut 1.1 S3 = 0;
707     S2 = 0;
708     S12 = 0;
709     S11 = 0;
710     S3 = trig->patterntrig[2];
711     S2 = trig->patterntrig[3];
712     S12 = trig->patterntrig[4];
713     S11 = trig->patterntrig[5];
714     if ( trig->patterntrig[0] ) tmptrigty = 2.;
715     if ( S3 || S2 || S12 || S11 ) tmptrigty = 0.;
716     if ( trig->patterntrig[1] & (1<<0) || (!trig->patterntrig[0] && !S3 && !S2 && !S12 && !S11) ) tmptrigty = 1.;
717     event->clevel2->trigty = tmptrigty;
718     //
719     // check if the calibration we are using is still good, if not load another calibration
720     //
721     sgnl = 0;
722     sgnl = event->ChkCalib(dbc,atime);
723     if ( sgnl < 0 ){
724     code = sgnl;
725     goto closeandexit;
726     };
727     if ( sgnl == 100 ){
728     code = sgnl;
729     if ( verbose ) printf(" CALORIMETER - WARNING: data not associated to any calibration interval\n");
730     badevent++;
731     sgnl = 0;
732     };
733     //
734 mocchiut 1.11 // do we have at least one track from the tracker? this check has been disabled
735 mocchiut 1.1 //
736     event->clevel1->good2 = 1;
737 mocchiut 1.8 //
738 mocchiut 1.1 //
739 mocchiut 1.8 // Calibrate calorimeter event "re" and store output in the two structures that will be passed to fortran routine
740 mocchiut 1.1 //
741 mocchiut 1.11 event->Calibrate(re);
742 mocchiut 1.1 //
743     // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract topological variables.
744     //
745     //
746     // Calculate variables common to all tracks (qtot, nstrip, etc.)
747     //
748     event->GetCommonVar();
749     //
750     // Fill common variables
751     //
752 mocchiut 1.13 event->FillCommonVar(c1,ca);
753 mocchiut 1.1 //
754     // Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
755     //
756     ntrkentry = 0;
757     //
758     filled = false;
759     //
760 mocchiut 1.11 // Run over tracks (tracker or calorimeter )
761     //
762     if ( withtrk ){
763     for(Int_t nt=0; nt < trk->ntrk(); nt++){
764     //
765     event->clevel1->good2 = 1;
766     //
767     TrkTrack *ptt = trk->GetStoredTrack(nt);
768     //
769     event->clevel1->trkchi2 = 0;
770     //
771     // Copy the alpha vector in the input structure
772     //
773     for (Int_t e = 0; e < 5 ; e++){
774     event->clevel1->al_p[e][0] = ptt->al[e];
775     };
776     //
777     // Get tracker related variables for this track
778     //
779     event->GetTrkVar();
780     //
781     // Save tracker track sequence number
782     //
783     event->trkseqno = nt;
784     //
785     // Copy values in the class ca from the structure clevel2
786     //
787     event->FillTrkVar(ca,ntrkentry);
788     ntrkentry++;
789     filled = true;
790     //
791     }; // loop on all the tracks
792     };
793     //
794     // if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter
795     // fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off).
796     // here we make an event selection so it must be done very carefully...
797     //
798     // conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event
799     //
800     if ( trackanyway && !filled && event->clevel2->npcfit[0] >= 2 && event->clevel2->npcfit[1] >= 2 && event->clevel2->good != 0 && event->clevel2->trigty < 2. ){
801     if ( debug ) printf(" Event with a track not fitted by the tracker at entry %i \n",itr);
802     //
803     // Disable "track mode" in the fortran routine
804     //
805     event->clevel1->good2 = 0;
806     event->clevel1->riginput = rigdefault;
807     if ( debug ) printf(" Using as default rigidity: %f \n",event->clevel1->riginput);
808     //
809     // We have a selftrigger event to analyze.
810     //
811     for (Int_t e = 0; e < 5 ; e++){
812     event->clevel1->al_p[e][0] = 0.;
813     event->clevel1->al_p[e][1] = 0.;
814     };
815     event->clevel1->trkchi2 = 0;
816     //
817     event->GetTrkVar();
818     //
819     // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
820 mocchiut 1.1 //
821 mocchiut 1.11 if ( event->clevel1->good2 == 0 ) {
822 mocchiut 1.1 //
823 mocchiut 1.11 // In selftrigger mode the trkentry variable is set to -1
824 mocchiut 1.1 //
825 mocchiut 1.11 event->trkseqno = -3;
826 mocchiut 1.1 //
827 mocchiut 1.11 // Copy values in the class ca from the structure clevel2
828 mocchiut 1.1 //
829 mocchiut 1.11 event->FillTrkVar(ca,ntrkentry);
830     ntrkentry++;
831     filled = true;
832     //
833     } else {
834     if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
835     };
836     //
837     };
838     //
839     // Call high energy nuclei routine
840     //
841     if ( hZn && event->clevel2->trigty >= 2. ){
842     if ( debug ) printf(" Calling selftrigger high energy nuclei routine for entry %i \n",itr);
843     //
844     // Disable "track mode" in the fortran routine
845     //
846     event->clevel1->good2 = 0;
847     //
848     // Set high energy nuclei flag to one
849     //
850     event->clevel1->hzn = 1;
851     event->clevel1->riginput = rigdefault;
852     //
853     // We have a selftrigger event to analyze.
854     //
855     for (Int_t e = 0; e < 5 ; e++){
856     event->clevel1->al_p[e][0] = 0.;
857     event->clevel1->al_p[e][1] = 0.;
858 mocchiut 1.1 };
859 mocchiut 1.11 event->clevel1->trkchi2 = 0;
860     //
861     event->GetTrkVar();
862     //
863     // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
864     //
865     if ( event->clevel1->good2 == 0 ) {
866     //
867     // In selftrigger mode the trkentry variable is set to -1
868     //
869     event->trkseqno = -2;
870     //
871     // Copy values in the class ca from the structure clevel2
872 mocchiut 1.1 //
873 mocchiut 1.11 event->FillTrkVar(ca,ntrkentry);
874     ntrkentry++;
875     filled = true;
876 mocchiut 1.1 //
877 mocchiut 1.11 } else {
878     if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
879     };
880     //
881     };
882     //
883     // self trigger event
884     //
885     if ( st && event->clevel2->trigty >= 2. ){
886     if ( verbose ) printf(" Selftrigger event at entry %i \n",itr);
887     //
888     // Disable "track mode" in the fortran routine
889     //
890     event->clevel1->good2 = 0;
891     //
892     // disable high enery nuclei flag;
893     //
894     event->clevel1->hzn = 0;
895     //
896     // We have a selftrigger event to analyze.
897     //
898     for (Int_t e = 0; e < 5 ; e++){
899     event->clevel1->al_p[e][0] = 0.;
900     event->clevel1->al_p[e][1] = 0.;
901     };
902     event->clevel1->trkchi2 = 0;
903     //
904     event->GetTrkVar();
905     //
906     // if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on
907     //
908     if ( event->clevel1->good2 == 0 ) {
909 mocchiut 1.1 //
910 mocchiut 1.11 // In selftrigger mode the trkentry variable is set to -1
911 mocchiut 1.1 //
912 mocchiut 1.11 event->trkseqno = -1;
913 mocchiut 1.1 //
914 mocchiut 1.11 // Copy values in the class ca from the structure clevel2
915 mocchiut 1.1 //
916 mocchiut 1.11 event->FillTrkVar(ca,ntrkentry);
917     ntrkentry++;
918     filled = true;
919 mocchiut 1.1 //
920 mocchiut 1.11 } else {
921     if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
922 mocchiut 1.1 };
923     };
924     //
925     if ( !filled ) badevent++;
926     //
927     // Clear structures used to communicate with fortran
928     //
929     event->ClearStructs();
930     //
931     // Fill the rootple
932     //
933     calo->Fill();
934     //
935 mocchiut 1.6 // delete ca;
936 mocchiut 1.1 //
937     jumpev:
938     debug = false;
939     //
940     };
941     //
942     if ( verbose ) printf("\n SUMMARY:\n");
943     if ( verbose ) printf(" Total number of events: %i \n",totevent);
944     if ( verbose ) printf(" Events with at least one track: %i \n",totevent-badevent);
945     if ( verbose ) printf(" Events without tracks: %i \n",badevent);
946     //
947     if ( badevent == totevent ){
948 mocchiut 1.8 if ( verbose ) printf("\n CALORIMETER - WARNING no tracks or good events in run %u \n",idRun);
949 mocchiut 1.1 code = 101;
950     };
951     //
952 mocchiut 1.8 delete dbtime;
953 mocchiut 1.1 //
954     // Clear variables before processing another run (needed to have always the same result when reprocessing data with the same software).
955     //
956     event->RunClose();
957     if ( l0File ) l0File->Close();
958     //
959     }; // process all the runs
960     //
961     if ( verbose ) printf("\n Finished processing data \n");
962     //
963     closeandexit:
964     //
965     if ( !reprocall && reproc && code >= 0 ){
966     if ( totfileentries > noaftrun ){
967     if ( verbose ) printf("\n Post-processing: copying events from the old tree after the processed run\n");
968 mocchiut 1.8 if ( verbose ) printf(" Copying %u events in the file which are after the end of the run %u \n",(totfileentries-noaftrun),run);
969     if ( verbose ) printf(" Start copying at event number %u end copying at event number %u \n",noaftrun,totfileentries);
970 mocchiut 1.1 for (UInt_t j = noaftrun; j < totfileentries; j++ ){
971     //
972     // Get entry from old tree
973     //
974     caloclone->GetEntry(j);
975     //
976     // copy caclone to ca
977     //
978 mocchiut 1.13 if ( getl1 ) c1->Clear();
979 mocchiut 1.6 ca->Clear();
980 mocchiut 1.8 //
981 mocchiut 1.13 if ( getl1 ) memcpy(&c1,&c1clone,sizeof(c1clone));
982 mocchiut 1.1 memcpy(&ca,&caclone,sizeof(caclone));
983     //
984     // Fill entry in the new tree
985     //
986     calo->Fill();
987     //
988     };
989     if ( verbose ) printf(" Finished successful copying!\n");
990     };
991     };
992     //
993     // Case of no errors: close files, delete old tree(s), write and close level2 file
994     //
995     if ( l0File ) l0File->Close();
996     if ( tempfile ) tempfile->Close();
997     gSystem->Unlink(tempname.str().c_str());
998     if ( tracker ) tracker->Delete();
999     //
1000     if ( code < 0 ) printf("\n CALORIMETER - ERROR: an error occurred, try to save anyway...\n");
1001     if ( verbose ) printf("\n Writing and closing rootple\n");
1002     if ( runinfo ) runinfo->Close();
1003     if ( calo ) calo->SetName("Calorimeter");
1004     if ( file ){
1005     file->cd();
1006     file->Write();
1007     };
1008     if ( calo ) calo->Delete();
1009     //
1010     gSystem->Unlink(calofolder.str().c_str());
1011     //
1012 mocchiut 1.7 if ( ca ) delete ca;
1013     if ( caclone ) delete caclone;
1014 mocchiut 1.13 if ( c1 ) delete c1;
1015     if ( c1clone ) delete c1clone;
1016 mocchiut 1.7 if ( trk ) delete trk;
1017     if ( q4 ) delete q4;
1018     if ( glroot ) delete glroot;
1019     if ( runinfo ) delete runinfo;
1020     //
1021 mocchiut 1.1 // the end
1022     //
1023     if ( verbose ) printf("\n Exiting...\n");
1024     if ( code < 0 ) throw code;
1025     return(code);
1026     }
1027    
1028    
1029    

  ViewVC Help
Powered by ViewVC 1.1.23