/[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.1 - (hide annotations) (download)
Fri May 19 13:15:54 2006 UTC (18 years, 6 months ago) by mocchiut
Branch: MAIN
Branch point for: DarthVader
Initial revision

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

  ViewVC Help
Powered by ViewVC 1.1.23