/[PAMELA software]/DarthVader/TrackerLevel2/src/TrkCore.cpp
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/TrkCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (hide annotations) (download)
Fri Aug 4 08:18:06 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.4: +13 -12 lines
some memory leak bugs fixed + CN computation modified

1 mocchiut 1.1 /**
2     * \file TrkCore.cpp
3     * \author Elena Vannuccini
4     */
5     // .........................................................
6     // ROOT header files
7     // .........................................................
8     #include <TTree.h>
9     #include <TClassEdit.h>
10     #include <TObject.h>
11     #include <TList.h>
12     #include <TSystem.h>
13     #include <TSystemDirectory.h>
14     #include <TString.h>
15     #include <TFile.h>
16     #include <TClass.h>
17     #include <TCanvas.h>
18     #include <TH1.h>
19     #include <TH1F.h>
20     #include <TH2D.h>
21     #include <TLatex.h>
22     #include <TPad.h>
23     #include <TSQLServer.h>
24     #include <TSQLRow.h>
25     #include <TSQLResult.h>
26     // .........................................................
27     // other general header files
28     // .........................................................
29     #include <fstream>
30     #include <iostream>
31     // .........................................................
32     // files in the inc directory
33     // .........................................................
34     #include <RunInfo.h>
35     #include <GLTables.h>
36     #include <TrkVerl2.h>
37     #include <TrkProcess.h>
38 pam-fi 1.2 //#include <TrkStruct.h>
39 mocchiut 1.1 #include <TrkLevel2.h>
40     //#include <TrkLevel1.h>
41     #include <TrkLevel0.h>
42     // .........................................................
43     // YODA header files
44     // .........................................................
45     #include <PamelaRun.h>
46     #include <RegistryEvent.h>
47     #include <physics/tracker/TrackerEvent.h>
48     #include <CalibTrk1Event.h>
49     #include <CalibTrk2Event.h>
50 pam-fi 1.2
51 mocchiut 1.1 // .........................................................
52     // namespaces
53     // .........................................................
54     using namespace std;
55     using namespace pamela;
56     // ================================================================================
57     //
58     //
59     // ================================================================================
60     /**
61     * \brief Tracker data reduction routine.
62     *
63     * It performs data reduction LEVEL0->LEVEL1->LEVEL2, producing LEVEL1 and or LEVEL2 output, in ROOT or HBOOK format.
64     * Input parameters:
65     * @param run id of the run to be processed (if run=0 a whole file is reprocessed)
66     * @param dbc pointer to BD server
67     * @param file1 output LEVEL1 file name (null if no LEVEL1 output is required)
68     * @param file2 output LEVEL2 file name (null if no LEVEL2 output is required)
69     * @param standalone (flag to run the program in standalone mode, that is without reading RunInfo)
70     */
71     //short int TrkCore(Int_t run, TSQLServer *dbc, TString file1, TString file2, Bool_t standalone)
72     int TrkCore(ULong64_t run, TFile *f2, TSQLServer *dbc, int ncustom, char*vcustom[])
73     {
74     // ---------------------------
75     // Define some basic variables
76     // ---------------------------
77    
78     TString filename = 0;
79     TString trkcalibfile = 0;
80     Int_t last_trk_calib_used = 0;
81 pam-fi 1.2 Long64_t nevents = 0LL;
82    
83 mocchiut 1.1 // LEVEL0 input
84     TFile *f0 = 0;
85 pam-fi 1.2 TFile *f0_c = 0;
86 mocchiut 1.1 TTree *physicsTree = 0;
87     TBranch *b_registry = 0;
88     TBranch *b_trk = 0;
89     TBranch *b_header = 0;
90     RegistryEvent *reg = 0;
91     EventHeader *header = 0;
92     PscuHeader *pscu = 0;
93     TrkLevel0 *l0_event = 0;
94     // RunInfo
95     ItoRunInfo *runinfo = 0;
96     TArrayL *runlist = 0;
97     ULong64_t from_run;
98     ULong64_t to_run;
99     Bool_t reprocessing = false;
100     //LEVEL1 output - ROOT
101     // TFile *f1;
102     // TTree *t_level1;
103     // TrkLevel1 *l1_event;
104     //LEVEL2 output - ROOT
105     TTree *t_level2 = 0;
106     TTree *t_clone = 0;
107     TrkLevel2 *l2_event = new TrkLevel2();
108     TrkLevel2 *l2_clone = new TrkLevel2();
109     // -----------------------
110     // -----------------------
111     // Handle input parameters
112     // (data)
113     //
114     // - open/create output files, determining the environment root/hbook from the estension
115     // - create the level1/level2 tree+branch/nt-uple
116     // -----------------------
117     // -----------------------
118     TrkProcess *p = new TrkProcess(run,f2);
119     p->HandleCustomPar(ncustom,vcustom);
120 pam-fi 1.2 if( p->DebugMode() )p->Dump();
121 mocchiut 1.1
122    
123    
124     // ===================================
125     // Open/Create level1 output file
126     // (NB output files other than LEVEL2 are
127     // created in a separate folder)
128     // ===================================
129     if(p->get1){
130     if(p->ifroot1){
131     throw -299;
132     }else{
133     throw -299;
134     };
135     };
136     // ===================================
137     // Open/Create level2 output file
138     // ===================================
139     if(p->get2){
140     // if(p->ifroot2){ // root
141     //-------------------------------------------
142     // read RunInfo
143     //-------------------------------------------
144     if(!(p->standalone)){
145     // Open "RunInfo" tree and get list of run
146     runinfo = new ItoRunInfo(f2);
147 pam-fi 1.3 char *trkversion = TrkInfo(false);
148     int runinfo_error = runinfo->Update(run,"TRK",trkversion);
149     if( runinfo_error ) throw runinfo_error;
150 mocchiut 1.1 runlist = runinfo->GetRunList();
151     reprocessing = runinfo->IsReprocessing();//????
152 pam-fi 1.2 if(p->VerboseMode()){
153 mocchiut 1.1 cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file
154     cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file
155     cout << "1-st run "<< runlist->At(0) << endl;
156     cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl;
157     cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run
158     cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run
159     cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl;
160     };
161     };
162     //-------------------------------------------
163     //
164     //-------------------------------------------
165     // Take (if present) the old Tracker tree
166     t_clone = (TTree*)f2->Get("Tracker");
167     if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone);
168     // Create NEW Tracker tree
169     t_level2 = new TTree("Tracker","PAMELA tracker level2 data ");
170     t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
171     // }else{ // hbook
172     // throw -299;
173     // }
174     };
175    
176    
177     // -------------------------------------------
178     // define runs to be processed/reprocessed
179     // -------------------------------------------
180     if(run == 0){
181     // reprocessing ALL runs
182     if(p->standalone)throw -298; // reprocessing not implemented
183     from_run = runlist->At(0);
184     to_run = runlist->At(runinfo->GetRunEntries()-1);
185     }else{
186     // processing/reprocessing ONE single run
187     from_run = run;
188     to_run = run;
189     };
190    
191     //
192     // init "expiration date" of calibrations and parameters
193     //
194     ULong64_t pedsig_time =0ULL;
195     ULong64_t B_time =0ULL;
196     ULong64_t mip_time =0ULL;
197     ULong64_t charge_time =0ULL;
198     ULong64_t eta_time =0ULL;
199     ULong64_t align_time =0ULL;
200     ULong64_t mask_time =0ULL;
201     //
202     // init event counter
203     //
204     Int_t ev_count =0;
205     //
206     // create query-results objects
207     //
208     GL_RUN q1 = GL_RUN();
209     GL_TRK_CALIB q2 = GL_TRK_CALIB();
210     GL_ROOT q3 = GL_ROOT();
211     GL_PARAM q4 = GL_PARAM();
212    
213     // ------------------------------------------------------------
214     // if reprocessing one run, copy all the events BEFORE the run
215     // ------------------------------------------------------------
216     if( !(p->standalone) ){
217     for(UInt_t i=0; i<runinfo->GetFirstEntry(); i++){
218     t_clone->GetEntry(i);
219     *l2_event = *l2_clone;
220     t_level2->Fill();
221     l2_event->Clear();
222     // COPY COPY COPY
223     };
224     };
225     // ------------------------------------------------------------
226     // ------------------------------------------------------------
227     // START LOOP OVER RUNS TO PROCESSED/REPROCESSED
228     // ------------------------------------------------------------
229     // ------------------------------------------------------------
230     for(ULong64_t idRun=from_run; idRun <= to_run; idRun++){
231    
232 pam-fi 1.2 if(p->VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
233 mocchiut 1.1 ULong64_t runheadtime = 0ULL;
234     ULong64_t runtrailtime = 0ULL;
235     UInt_t evfrom = 0;
236     UInt_t evto = 0;
237     Int_t id_reg_run =-1;
238     Int_t trk_calib_used = 0;
239    
240     if(p->standalone){
241     // ==============================
242     // first query: retrieve run info
243     // ==============================
244     if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
245     id_reg_run = q1.ID_REG_RUN;
246     runheadtime = q1.RUNHEADER_TIME;
247     runtrailtime = q1.RUNTRAILER_TIME;
248     evfrom = q1.EV_REG_PHYS_FROM;
249     evto = q1.EV_REG_PHYS_TO;
250     trk_calib_used = q1.TRK_CALIB_USED;
251     }else{
252     // ==============================
253     // get run info from RunInfo tree
254     // ==============================
255 pam-fi 1.3 int runinfo_error = runinfo->GetRunInfo(idRun);
256     if( runinfo_error ) throw runinfo_error;
257 mocchiut 1.1 id_reg_run = runinfo->ID_REG_RUN;
258     runheadtime = runinfo->RUNHEADER_TIME;
259     runtrailtime = runinfo->RUNTRAILER_TIME;
260     evfrom = runinfo->EV_REG_PHYS_FROM;
261     evto = runinfo->EV_REG_PHYS_TO;
262     trk_calib_used = runinfo->TRK_CALIB_USED;
263     };
264     //
265 pam-fi 1.2 if(p->VerboseMode()){
266 mocchiut 1.1 cout << "ROOT file ID "<< id_reg_run << endl;
267     cout << "RunHeader time "<< runheadtime << endl;
268     cout << "RunTrailer time "<< runtrailtime << endl;
269     cout << " from event "<< evfrom << endl;
270     cout << " to event "<< evto << endl;
271     cout << "trk-calibration used "<< trk_calib_used << endl;
272     };
273     // ========================================================
274     // second query: search the LEVEL0 file that contains idRun
275     // ========================================================
276     TString lastfilename = filename;
277     if( q3.Query_GL_ROOT(id_reg_run,dbc) )throw -51;
278     filename = q3.PATH + q3.NAME;
279     // ========================================================
280     // Open the input LEVEL0 data file
281     // ========================================================
282     if(filename.CompareTo(lastfilename)){
283     if(!lastfilename.IsNull())f0->Close();
284     //if( debug ) cout << "Opening LEVEL0 file: "<< filename << endl;
285 pam-fi 1.2 if(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
286     FileStat_t t;
287     if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
288 mocchiut 1.1 f0 = new TFile(filename);
289     if ( !f0 ) throw -6;
290     physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7;
291     b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8;
292     b_registry = physicsTree ->GetBranch("Registry"); if(!b_registry ) throw -9;
293     b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
294     physicsTree->SetBranchAddress("Tracker" ,&l0_event);
295     physicsTree->SetBranchAddress("Registry",&reg);
296     physicsTree->SetBranchAddress("Header",&header);
297    
298     nevents = physicsTree->GetEntries();
299     if ( nevents < 1 ) throw -11;
300 pam-fi 1.3 if ( nevents < (evto+1)) throw -12;
301 mocchiut 1.1 };
302    
303     // =============================================================
304 pam-fi 1.2 // retrieve calibration file needed to reduce data
305 mocchiut 1.1 // =============================================================
306 pam-fi 1.2 // if run OBT is > last calibration "expiration date"
307     // - search for new calibration packet
308     // - load calibration parameters (full + truncated)
309 pam-fi 1.5 if(p->VerboseMode() )cout << "Full pedestals for cluster finding:";
310 pam-fi 1.2 if(runheadtime > pedsig_time){
311     if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
312     pedsig_time = q2.TO_TIME;
313     if(q2.EV_REG_CALIBTRK1 != q2.EV_REG_CALIBTRK2)
314     printf("WARNING!! ---> EV_REG_CALIBTRK1=%d it's different from EV_REG_CALIBTRK2=%d \n\n",q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
315     if( q3.Query_GL_ROOT(q2.ID_REG_CALIBTRK,dbc) )throw -51;
316 pam-fi 1.5 trkcalibfile = q3.PATH + q3.NAME;
317     if(p->VerboseMode()) cout << " >> Loading from LEVEL0 file: "<< trkcalibfile << endl;
318 pam-fi 1.2 if(trkcalibfile.CompareTo(filename)){
319 pam-fi 1.5 FileStat_t t;
320     if( gSystem->GetPathInfo(trkcalibfile.Data(),t) )throw -6;
321     f0_c=new TFile(trkcalibfile);
322     if ( !f0_c ) throw -6;
323 pam-fi 1.2 }else{
324 pam-fi 1.5 f0_c = f0;
325 pam-fi 1.2 };
326     if(p->VerboseMode()) {
327 pam-fi 1.5 cout << " calibration entries "<< q2.EV_REG_CALIBTRK1 << " " << q2.EV_REG_CALIBTRK2;
328 pam-fi 1.2 cout << " (from time "<< q2.FROM_TIME <<" to time "<< q2.TO_TIME <<")"<<endl;
329     };
330     pedsigbad_.FillACalibFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
331 pam-fi 1.5 }else{
332     if(p->VerboseMode() )cout<< " already loaded. "<< endl;
333 pam-fi 1.2 };
334     if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal[64][12][i] << endl;
335     // =============================================================
336     // retrieve calibration file needed to uncompress data
337     // =============================================================
338     // if the run was compressed using default calib
339     // load truncated pedestals from default
340     // otherwise reload them from on-line calibration
341     if(p->VerboseMode() )cout << "Truncated pedestals for uncompression:";
342 mocchiut 1.1 if(trk_calib_used==104 && last_trk_calib_used!=104){
343    
344 pam-fi 1.5 if(p->VerboseMode() )cout << " >> Loading default calibration " << endl;
345 pam-fi 1.2 if (q4.Query_GL_PARAM(runheadtime,"default calibration",dbc) )throw -52;
346     pedsigbad_.FillTCalibFrom(q4.PATH+q4.NAME);
347    
348     }else if( (trk_calib_used==1 | trk_calib_used==2) && last_trk_calib_used==104 ){
349    
350     if(p->VerboseMode()) cout << ">> Re-loading on-line calibration " << endl;
351     pedsigbad_.FillTCalibFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
352     }else{
353     if(p->VerboseMode()) cout << ">> Using on-line calibration " << endl;
354     };
355     if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal_t[64][12][i] << endl;
356     last_trk_calib_used = trk_calib_used;
357    
358     /* // =============================================================
359     // retrieve calibration file needed to uncompress data
360     // =============================================================
361     if(trk_calib_used==104 && last_trk_calib_used!=104){
362    
363     if(p->VerboseMode() )cout << "** Using default calibration **" << endl;
364 mocchiut 1.1
365     if (q4.Query_GL_PARAM(runheadtime,"default calibration",dbc) )throw -52;
366     path_.FillWith(q4.PATH+q4.NAME);
367     fillpedsigfromdefault_();
368     if(path_.error) throw -216;
369    
370     }else if( (trk_calib_used==1 | trk_calib_used==2) && runheadtime > pedsig_time){
371    
372     if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
373     pedsig_time = q2.TO_TIME;
374     if(q2.EV_REG_CALIBTRK1 != q2.EV_REG_CALIBTRK2)
375     printf("WARNING!! ---> EV_REG_CALIBTRK1=%d it's different from EV_REG_CALIBTRK2=%d \n\n",q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
376    
377     if( q3.Query_GL_ROOT(q2.ID_REG_CALIBTRK,dbc) )throw -51;
378    
379     trkcalibfile = q3.PATH + q3.NAME;
380     // store calibration informations
381 pam-fi 1.2 if(p->VerboseMode()) cout << "Online calibration: "<< trkcalibfile << endl;
382 mocchiut 1.1 // pedsigbad_ = TrkFlightPede(trkcalibfile,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
383     TFile* f0_c= f0;
384     if(trkcalibfile.CompareTo(filename))f0_c = new TFile(trkcalibfile);
385    
386 pam-fi 1.2 pedsigbad_.FillAllFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
387 mocchiut 1.1
388     if(trkcalibfile.CompareTo(filename))f0_c->Close();
389     };
390 pam-fi 1.2 last_trk_calib_used = trk_calib_used;*/
391 mocchiut 1.1 // =============================================================
392     // retrieve information about parameters to process LEVEL2
393     // =============================================================
394     //
395     // magnetic field
396     //
397 pam-fi 1.2 if(p->VerboseMode()) cout << "Loading parameter files : "<< endl;
398    
399 mocchiut 1.1 if(runheadtime > B_time){
400     if( q4.Query_GL_PARAM(runheadtime,"field",dbc) )throw -52;
401     B_time = q4.TO_TIME;
402     l2_event->LoadField(q4.PATH+q4.NAME);
403 pam-fi 1.2 // l2_event->LoadField(q4.PATH+q4.NAME);
404 mocchiut 1.1 if(path_.error) throw -215;
405     };
406     //
407     // mip conversion parameters
408     //
409     if(runheadtime > mip_time){
410     if( q4.Query_GL_PARAM(runheadtime,"trk mip",dbc) )throw -52;
411     mip_time = q4.TO_TIME;
412     path_.FillWith(q4.PATH+q4.NAME);
413     readmipparam_();
414     if(path_.error) throw -212;
415     };
416     //
417     // charge correlation parameters
418     //
419     if(runheadtime > charge_time){
420     if( q4.Query_GL_PARAM(runheadtime,"trk charge",dbc) )throw -52;
421     charge_time = q4.TO_TIME;
422     path_.FillWith(q4.PATH+q4.NAME);
423     readchargeparam_();
424     if(path_.error) throw -213;
425     };
426     //
427     // eta p.f.a. parameters
428     //
429     if(runheadtime > eta_time){
430     if( q4.Query_GL_PARAM(runheadtime,"trk pfa",dbc) )throw -52;
431     eta_time = q4.TO_TIME;
432     path_.FillWith(q4.PATH+q4.NAME);
433     readetaparam_();
434     if(path_.error) throw -214;
435     };
436     //
437     // alignment parameters
438     //
439     if(runheadtime > align_time){
440     if( q4.Query_GL_PARAM(runheadtime,"trk alignment",dbc) )throw -52;
441     align_time = q4.TO_TIME;
442     path_.FillWith(q4.PATH+q4.NAME);
443     readalignparam_();
444     if(path_.error) throw -211;
445     };
446     //
447     // viking mask
448     //
449     if(runheadtime > mask_time){
450     if( q4.Query_GL_PARAM(runheadtime,"trk mask",dbc) )throw -52;
451     mask_time = q4.TO_TIME;
452     path_.FillWith(q4.PATH+q4.NAME);
453     readvkmask_();
454     if(path_.error) throw -210;
455     };
456    
457     //
458     // ====================================================
459     // start looping on events cointained in the data file
460     // ====================================================
461     for (UInt_t re = evfrom; re <= evto; re++){
462    
463     ev_count++;
464     // if ( p->DEBUG && re%100 == 0 && re > 0 ) printf(" %iK \n",re/1000);
465 pam-fi 1.2 if ( p->DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
466 mocchiut 1.1
467     b_registry->GetEntry(re);
468     b_trk->GetEntry(reg->event);
469     b_header->GetEntry(reg->event);
470     pscu = header->GetPscuHeader();
471     // pscu->Print();
472     // cout <<endl<< " *************** Pkt ID "<< (UINT32)pscu->Counter << endl;
473    
474 pam-fi 1.2 p->ProcessEvent(l0_event);
475 pam-fi 1.5
476 mocchiut 1.1 // ----------------
477     // LEVEL1 output
478     // ----------------
479     if(p->get1){
480     if(p->ifroot1){ // root
481     }else{ // hbook
482     };
483     };
484     // ----------------
485     // LEVEL2 output
486     // ----------------
487     if(p->get2){
488     // if(p->ifroot2){ // root
489 pam-fi 1.4 l2_event->SetFromLevel2Struct(&level2event_);
490 mocchiut 1.1 // l2_event->Dump();
491     t_level2->Fill();
492     l2_event->Clear();
493     // }else{ // hbook
494     // }
495     };
496    
497     }; // end loop on events
498 pam-fi 1.2 if(p->VerboseMode())cout << " >>> processed "<< ev_count <<" events"<< endl;
499 mocchiut 1.1 }; // end loop on runs
500     // ------------------------------------------------------------
501     // if reprocessing one run, copy all the events AFTER the run
502     // ------------------------------------------------------------
503     if( !(p->standalone) ){
504     for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
505     t_clone->GetEntry(i);
506     *l2_event = *l2_clone;
507     t_level2->Fill();
508     l2_event->Clear();
509     // COPY COPY COPY
510     };
511     };
512     // ---------------
513     // close the files
514     // ---------------
515     if(p->get1){
516     if(p->ifroot1){ // root
517     }else{ // hbook
518     }
519     };
520     if(p->get2){
521     // if(p->ifroot2){ // root
522     if( t_clone )t_clone->Delete("all");//delete old tree from file
523     if( !(p->standalone) )runinfo->Close();
524     f2->Write("Tracker");
525     if( t_level2 )t_level2->Delete(); //delete new tree from memory
526    
527     // }else{ // hbook
528     // }
529     };
530 pam-fi 1.2 if( f0->IsOpen() ) f0->Close();
531     if( f0_c->IsOpen() )f0_c->Close();
532    
533     return(p->ostatus);
534 mocchiut 1.1 }
535    

  ViewVC Help
Powered by ViewVC 1.1.23