/[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.20 - (hide annotations) (download)
Tue Jan 23 15:36:42 2007 UTC (17 years, 10 months ago) by mocchiut
Branch: MAIN
CVS Tags: v3r00
Changes since 1.19: +4 -1 lines
Compiling bug under fedora4 fixed, Syntax error inc/TrkLevelX bug fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23