/[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.2 - (hide annotations) (download)
Tue May 30 16:30:37 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.1: +81 -58 lines
Error handling from F77 routine / Fixed some bugs with default calibration

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     char *trkversion = TrkInfo(false);
148     if( runinfo->Update(run,"TRK",trkversion) ) throw -205;
149     runlist = runinfo->GetRunList();
150     reprocessing = runinfo->IsReprocessing();//????
151 pam-fi 1.2 if(p->VerboseMode()){
152 mocchiut 1.1 cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file
153     cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file
154     cout << "1-st run "<< runlist->At(0) << endl;
155     cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl;
156     cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run
157     cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run
158     cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl;
159     };
160     };
161     //-------------------------------------------
162     //
163     //-------------------------------------------
164     // Take (if present) the old Tracker tree
165     t_clone = (TTree*)f2->Get("Tracker");
166     if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone);
167     // Create NEW Tracker tree
168     t_level2 = new TTree("Tracker","PAMELA tracker level2 data ");
169     t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
170     // }else{ // hbook
171     // throw -299;
172     // }
173     };
174    
175    
176     // -------------------------------------------
177     // define runs to be processed/reprocessed
178     // -------------------------------------------
179     if(run == 0){
180     // reprocessing ALL runs
181     if(p->standalone)throw -298; // reprocessing not implemented
182     from_run = runlist->At(0);
183     to_run = runlist->At(runinfo->GetRunEntries()-1);
184     }else{
185     // processing/reprocessing ONE single run
186     from_run = run;
187     to_run = run;
188     };
189    
190     //
191     // init "expiration date" of calibrations and parameters
192     //
193     ULong64_t pedsig_time =0ULL;
194     ULong64_t B_time =0ULL;
195     ULong64_t mip_time =0ULL;
196     ULong64_t charge_time =0ULL;
197     ULong64_t eta_time =0ULL;
198     ULong64_t align_time =0ULL;
199     ULong64_t mask_time =0ULL;
200     //
201     // init event counter
202     //
203     Int_t ev_count =0;
204     //
205     // create query-results objects
206     //
207     GL_RUN q1 = GL_RUN();
208     GL_TRK_CALIB q2 = GL_TRK_CALIB();
209     GL_ROOT q3 = GL_ROOT();
210     GL_PARAM q4 = GL_PARAM();
211    
212     // ------------------------------------------------------------
213     // if reprocessing one run, copy all the events BEFORE the run
214     // ------------------------------------------------------------
215     if( !(p->standalone) ){
216     for(UInt_t i=0; i<runinfo->GetFirstEntry(); i++){
217     t_clone->GetEntry(i);
218     *l2_event = *l2_clone;
219     t_level2->Fill();
220     l2_event->Clear();
221     // COPY COPY COPY
222     };
223     };
224     // ------------------------------------------------------------
225     // ------------------------------------------------------------
226     // START LOOP OVER RUNS TO PROCESSED/REPROCESSED
227     // ------------------------------------------------------------
228     // ------------------------------------------------------------
229     for(ULong64_t idRun=from_run; idRun <= to_run; idRun++){
230    
231 pam-fi 1.2 if(p->VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
232 mocchiut 1.1 ULong64_t runheadtime = 0ULL;
233     ULong64_t runtrailtime = 0ULL;
234     UInt_t evfrom = 0;
235     UInt_t evto = 0;
236     Int_t id_reg_run =-1;
237     Int_t trk_calib_used = 0;
238    
239     if(p->standalone){
240     // ==============================
241     // first query: retrieve run info
242     // ==============================
243     if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
244     id_reg_run = q1.ID_REG_RUN;
245     runheadtime = q1.RUNHEADER_TIME;
246     runtrailtime = q1.RUNTRAILER_TIME;
247     evfrom = q1.EV_REG_PHYS_FROM;
248     evto = q1.EV_REG_PHYS_TO;
249     trk_calib_used = q1.TRK_CALIB_USED;
250     }else{
251     // ==============================
252     // get run info from RunInfo tree
253     // ==============================
254     if( runinfo->GetRunInfo(idRun) ) throw -205;
255     id_reg_run = runinfo->ID_REG_RUN;
256     runheadtime = runinfo->RUNHEADER_TIME;
257     runtrailtime = runinfo->RUNTRAILER_TIME;
258     evfrom = runinfo->EV_REG_PHYS_FROM;
259     evto = runinfo->EV_REG_PHYS_TO;
260     trk_calib_used = runinfo->TRK_CALIB_USED;
261     };
262     //
263 pam-fi 1.2 if(p->VerboseMode()){
264 mocchiut 1.1 cout << "ROOT file ID "<< id_reg_run << endl;
265     cout << "RunHeader time "<< runheadtime << endl;
266     cout << "RunTrailer time "<< runtrailtime << endl;
267     cout << " from event "<< evfrom << endl;
268     cout << " to event "<< evto << endl;
269     cout << "trk-calibration used "<< trk_calib_used << endl;
270     };
271     // ========================================================
272     // second query: search the LEVEL0 file that contains idRun
273     // ========================================================
274     TString lastfilename = filename;
275     if( q3.Query_GL_ROOT(id_reg_run,dbc) )throw -51;
276     filename = q3.PATH + q3.NAME;
277     // ========================================================
278     // Open the input LEVEL0 data file
279     // ========================================================
280     if(filename.CompareTo(lastfilename)){
281     if(!lastfilename.IsNull())f0->Close();
282     //if( debug ) cout << "Opening LEVEL0 file: "<< filename << endl;
283 pam-fi 1.2 if(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
284     FileStat_t t;
285     if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
286 mocchiut 1.1 f0 = new TFile(filename);
287     if ( !f0 ) throw -6;
288     physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7;
289     b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8;
290     b_registry = physicsTree ->GetBranch("Registry"); if(!b_registry ) throw -9;
291     b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
292     physicsTree->SetBranchAddress("Tracker" ,&l0_event);
293     physicsTree->SetBranchAddress("Registry",&reg);
294     physicsTree->SetBranchAddress("Header",&header);
295    
296     nevents = physicsTree->GetEntries();
297     if ( nevents < 1 ) throw -11;
298     };
299    
300     // =============================================================
301 pam-fi 1.2 // retrieve calibration file needed to reduce data
302 mocchiut 1.1 // =============================================================
303 pam-fi 1.2 // if run OBT is > last calibration "expiration date"
304     // - search for new calibration packet
305     // - load calibration parameters (full + truncated)
306     if(p->VerboseMode() )cout << "Full pedestals for cluster finding:" << endl;
307     if(runheadtime > pedsig_time){
308     if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
309     pedsig_time = q2.TO_TIME;
310     if(q2.EV_REG_CALIBTRK1 != q2.EV_REG_CALIBTRK2)
311     printf("WARNING!! ---> EV_REG_CALIBTRK1=%d it's different from EV_REG_CALIBTRK2=%d \n\n",q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
312     if( q3.Query_GL_ROOT(q2.ID_REG_CALIBTRK,dbc) )throw -51;
313     trkcalibfile = q3.PATH + q3.NAME;
314     if(trkcalibfile.CompareTo(filename)){
315     // if( f0_c->IsOpen() )f0_c->Close();
316     if(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< trkcalibfile << endl;
317     FileStat_t t;
318     if( gSystem->GetPathInfo(trkcalibfile.Data(),t) )throw -6;
319     f0_c=new TFile(trkcalibfile);
320     if ( !f0_c ) throw -6;
321     }else{
322     f0_c = f0;
323     };
324     if(p->VerboseMode()) {
325     cout << "Loading calibration entries "<< q2.EV_REG_CALIBTRK1 << " " << q2.EV_REG_CALIBTRK2;
326     cout << " (from time "<< q2.FROM_TIME <<" to time "<< q2.TO_TIME <<")"<<endl;
327     };
328     pedsigbad_.FillACalibFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
329     };
330     if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal[64][12][i] << endl;
331     // =============================================================
332     // retrieve calibration file needed to uncompress data
333     // =============================================================
334     // if the run was compressed using default calib
335     // load truncated pedestals from default
336     // otherwise reload them from on-line calibration
337     if(p->VerboseMode() )cout << "Truncated pedestals for uncompression:";
338 mocchiut 1.1 if(trk_calib_used==104 && last_trk_calib_used!=104){
339    
340 pam-fi 1.2 if(p->VerboseMode() )cout << ">> Loading default calibration " << endl;
341     if (q4.Query_GL_PARAM(runheadtime,"default calibration",dbc) )throw -52;
342     pedsigbad_.FillTCalibFrom(q4.PATH+q4.NAME);
343    
344     }else if( (trk_calib_used==1 | trk_calib_used==2) && last_trk_calib_used==104 ){
345    
346     if(p->VerboseMode()) cout << ">> Re-loading on-line calibration " << endl;
347     pedsigbad_.FillTCalibFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
348     }else{
349     if(p->VerboseMode()) cout << ">> Using on-line calibration " << endl;
350     };
351     if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal_t[64][12][i] << endl;
352     last_trk_calib_used = trk_calib_used;
353    
354     /* // =============================================================
355     // retrieve calibration file needed to uncompress data
356     // =============================================================
357     if(trk_calib_used==104 && last_trk_calib_used!=104){
358    
359     if(p->VerboseMode() )cout << "** Using default calibration **" << endl;
360 mocchiut 1.1
361     if (q4.Query_GL_PARAM(runheadtime,"default calibration",dbc) )throw -52;
362     path_.FillWith(q4.PATH+q4.NAME);
363     fillpedsigfromdefault_();
364     if(path_.error) throw -216;
365    
366     }else if( (trk_calib_used==1 | trk_calib_used==2) && runheadtime > pedsig_time){
367    
368     if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
369     pedsig_time = q2.TO_TIME;
370     if(q2.EV_REG_CALIBTRK1 != q2.EV_REG_CALIBTRK2)
371     printf("WARNING!! ---> EV_REG_CALIBTRK1=%d it's different from EV_REG_CALIBTRK2=%d \n\n",q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
372    
373     if( q3.Query_GL_ROOT(q2.ID_REG_CALIBTRK,dbc) )throw -51;
374    
375     trkcalibfile = q3.PATH + q3.NAME;
376     // store calibration informations
377 pam-fi 1.2 if(p->VerboseMode()) cout << "Online calibration: "<< trkcalibfile << endl;
378 mocchiut 1.1 // pedsigbad_ = TrkFlightPede(trkcalibfile,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
379     TFile* f0_c= f0;
380     if(trkcalibfile.CompareTo(filename))f0_c = new TFile(trkcalibfile);
381    
382 pam-fi 1.2 pedsigbad_.FillAllFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
383 mocchiut 1.1
384     if(trkcalibfile.CompareTo(filename))f0_c->Close();
385     };
386 pam-fi 1.2 last_trk_calib_used = trk_calib_used;*/
387 mocchiut 1.1 // =============================================================
388     // retrieve information about parameters to process LEVEL2
389     // =============================================================
390     //
391     // magnetic field
392     //
393 pam-fi 1.2 if(p->VerboseMode()) cout << "Loading parameter files : "<< endl;
394    
395 mocchiut 1.1 if(runheadtime > B_time){
396     if( q4.Query_GL_PARAM(runheadtime,"field",dbc) )throw -52;
397     B_time = q4.TO_TIME;
398     l2_event->LoadField(q4.PATH+q4.NAME);
399 pam-fi 1.2 // l2_event->LoadField(q4.PATH+q4.NAME);
400 mocchiut 1.1 if(path_.error) throw -215;
401     };
402     //
403     // mip conversion parameters
404     //
405     if(runheadtime > mip_time){
406     if( q4.Query_GL_PARAM(runheadtime,"trk mip",dbc) )throw -52;
407     mip_time = q4.TO_TIME;
408     path_.FillWith(q4.PATH+q4.NAME);
409     readmipparam_();
410     if(path_.error) throw -212;
411     };
412     //
413     // charge correlation parameters
414     //
415     if(runheadtime > charge_time){
416     if( q4.Query_GL_PARAM(runheadtime,"trk charge",dbc) )throw -52;
417     charge_time = q4.TO_TIME;
418     path_.FillWith(q4.PATH+q4.NAME);
419     readchargeparam_();
420     if(path_.error) throw -213;
421     };
422     //
423     // eta p.f.a. parameters
424     //
425     if(runheadtime > eta_time){
426     if( q4.Query_GL_PARAM(runheadtime,"trk pfa",dbc) )throw -52;
427     eta_time = q4.TO_TIME;
428     path_.FillWith(q4.PATH+q4.NAME);
429     readetaparam_();
430     if(path_.error) throw -214;
431     };
432     //
433     // alignment parameters
434     //
435     if(runheadtime > align_time){
436     if( q4.Query_GL_PARAM(runheadtime,"trk alignment",dbc) )throw -52;
437     align_time = q4.TO_TIME;
438     path_.FillWith(q4.PATH+q4.NAME);
439     readalignparam_();
440     if(path_.error) throw -211;
441     };
442     //
443     // viking mask
444     //
445     if(runheadtime > mask_time){
446     if( q4.Query_GL_PARAM(runheadtime,"trk mask",dbc) )throw -52;
447     mask_time = q4.TO_TIME;
448     path_.FillWith(q4.PATH+q4.NAME);
449     readvkmask_();
450     if(path_.error) throw -210;
451     };
452    
453     //
454     // ====================================================
455     // start looping on events cointained in the data file
456     // ====================================================
457     for (UInt_t re = evfrom; re <= evto; re++){
458    
459     ev_count++;
460     // if ( p->DEBUG && re%100 == 0 && re > 0 ) printf(" %iK \n",re/1000);
461 pam-fi 1.2 if ( p->DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
462 mocchiut 1.1
463     b_registry->GetEntry(re);
464     b_trk->GetEntry(reg->event);
465     b_header->GetEntry(reg->event);
466     pscu = header->GetPscuHeader();
467     // pscu->Print();
468     // cout <<endl<< " *************** Pkt ID "<< (UINT32)pscu->Counter << endl;
469    
470 pam-fi 1.2 p->ProcessEvent(l0_event);
471 mocchiut 1.1
472     // ----------------
473     // LEVEL1 output
474     // ----------------
475     if(p->get1){
476     if(p->ifroot1){ // root
477     }else{ // hbook
478     };
479     };
480     // ----------------
481     // LEVEL2 output
482     // ----------------
483     if(p->get2){
484     // if(p->ifroot2){ // root
485     l2_event->FillCommonVar(&level2event_);
486     // l2_event->Dump();
487     t_level2->Fill();
488     l2_event->Clear();
489     // }else{ // hbook
490     // }
491     };
492    
493     }; // end loop on events
494 pam-fi 1.2 if(p->VerboseMode())cout << " >>> processed "<< ev_count <<" events"<< endl;
495 mocchiut 1.1 }; // end loop on runs
496     // ------------------------------------------------------------
497     // if reprocessing one run, copy all the events AFTER the run
498     // ------------------------------------------------------------
499     if( !(p->standalone) ){
500     for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
501     t_clone->GetEntry(i);
502     *l2_event = *l2_clone;
503     t_level2->Fill();
504     l2_event->Clear();
505     // COPY COPY COPY
506     };
507     };
508     // ---------------
509     // close the files
510     // ---------------
511     if(p->get1){
512     if(p->ifroot1){ // root
513     }else{ // hbook
514     }
515     };
516     if(p->get2){
517     // if(p->ifroot2){ // root
518     if( t_clone )t_clone->Delete("all");//delete old tree from file
519     if( !(p->standalone) )runinfo->Close();
520     f2->Write("Tracker");
521     if( t_level2 )t_level2->Delete(); //delete new tree from memory
522    
523     // }else{ // hbook
524     // }
525     };
526 pam-fi 1.2 if( f0->IsOpen() ) f0->Close();
527     if( f0_c->IsOpen() )f0_c->Close();
528    
529     return(p->ostatus);
530 mocchiut 1.1 }
531    

  ViewVC Help
Powered by ViewVC 1.1.23