/[PAMELA software]/calo/flight/FCaloREPROC/src/CaloCore.cpp
ViewVC logotype

Annotation of /calo/flight/FCaloREPROC/src/CaloCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Fri Oct 26 11:41:32 2007 UTC (17 years, 1 month ago) by mocchiut
Branch point for: MAIN, FCaloREPROC
Initial revision

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

  ViewVC Help
Powered by ViewVC 1.1.23