/[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.38 - (hide annotations) (download)
Tue May 15 14:31:41 2012 UTC (12 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.37: +3 -3 lines
Reprocessing bugs fixed, ToF bugs fixed, new software versions, new quaternions, IGRF bug fixed, code cleanup

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 pam-fi 1.11 #include <TBenchmark.h>
27    
28 mocchiut 1.1 // .........................................................
29     // other general header files
30     // .........................................................
31     #include <fstream>
32     #include <iostream>
33     // .........................................................
34     // files in the inc directory
35     // .........................................................
36     #include <RunInfo.h>
37     #include <GLTables.h>
38     #include <TrkVerl2.h>
39     #include <TrkProcess.h>
40 pam-fi 1.2 //#include <TrkStruct.h>
41 mocchiut 1.1 #include <TrkLevel2.h>
42 pam-fi 1.6 #include <TrkLevel1.h>
43 mocchiut 1.1 #include <TrkLevel0.h>
44 pam-fi 1.15 #include <TrkHough.h>
45 mocchiut 1.1 // .........................................................
46     // YODA header files
47     // .........................................................
48     #include <PamelaRun.h>
49 pam-fi 1.7 //#include <RegistryEvent.h>
50 mocchiut 1.1 #include <physics/tracker/TrackerEvent.h>
51     #include <CalibTrk1Event.h>
52     #include <CalibTrk2Event.h>
53 pam-fi 1.2
54 mocchiut 1.1 // .........................................................
55     // namespaces
56     // .........................................................
57     using namespace std;
58     using namespace pamela;
59     // ================================================================================
60     //
61     //
62     // ================================================================================
63     /**
64     * \brief Tracker data reduction routine.
65     *
66     * It performs data reduction LEVEL0->LEVEL1->LEVEL2, producing LEVEL1 and or LEVEL2 output, in ROOT or HBOOK format.
67     * Input parameters:
68     * @param run id of the run to be processed (if run=0 a whole file is reprocessed)
69     * @param dbc pointer to BD server
70     * @param file1 output LEVEL1 file name (null if no LEVEL1 output is required)
71     * @param file2 output LEVEL2 file name (null if no LEVEL2 output is required)
72     * @param standalone (flag to run the program in standalone mode, that is without reading RunInfo)
73     */
74     //short int TrkCore(Int_t run, TSQLServer *dbc, TString file1, TString file2, Bool_t standalone)
75 mocchiut 1.31 //int TrkCore(UInt_t run, TFile *f2, TSQLServer *dbc, int ncustom, char*vcustom[]) // EMILIANO
76     int TrkCore(UInt_t run, TFile *f2, GL_TABLES *glt, int ncustom, char*vcustom[]) // EMILIANO
77 mocchiut 1.1 {
78     // ---------------------------
79     // Define some basic variables
80     // ---------------------------
81    
82 pizzolot 1.36 TString filename = "";
83 pam-fi 1.7 Long64_t nentries = 0LL;
84 pam-fi 1.2
85 mocchiut 1.1 // LEVEL0 input
86     TFile *f0 = 0;
87     TTree *physicsTree = 0;
88     TBranch *b_trk = 0;
89     TBranch *b_header = 0;
90     EventHeader *header = 0;
91     PscuHeader *pscu = 0;
92 pam-fi 1.22 TrkLevel0 *l0_event = new TrkLevel0();
93 mocchiut 1.1 // RunInfo
94     ItoRunInfo *runinfo = 0;
95 pam-fi 1.7 TArrayI *runlist = 0;
96     UInt_t from_run;
97     UInt_t to_run;
98 mocchiut 1.1 Bool_t reprocessing = false;
99     //LEVEL2 output - ROOT
100     TTree *t_level2 = 0;
101     TTree *t_clone = 0;
102     TrkLevel2 *l2_event = new TrkLevel2();
103     TrkLevel2 *l2_clone = new TrkLevel2();
104 pam-fi 1.13 //LEVEL1 output - ROOT
105     TrkLevel1 *l1_event = new TrkLevel1();
106 pam-fi 1.15 TrkHough *lh_event = new TrkHough();
107 mocchiut 1.1 // -----------------------
108     // -----------------------
109     // Handle input parameters
110     // (data)
111     //
112     // - open/create output files, determining the environment root/hbook from the estension
113     // - create the level1/level2 tree+branch/nt-uple
114     // -----------------------
115     // -----------------------
116     TrkProcess *p = new TrkProcess(run,f2);
117 pam-fi 1.26 if( p->HandleCustomPar(ncustom,vcustom) )return(p->ostatus);;
118 pam-fi 1.25 if( TrkParams::VerboseMode() )p->Dump();
119 mocchiut 1.1
120     // ===================================
121     // Open/Create level2 output file
122     // ===================================
123     if(p->get2){
124 pam-fi 1.13 //-------------------------------------------
125     // read RunInfo
126     //-------------------------------------------
127     if(!(p->standalone)){
128     // Open "RunInfo" tree and get list of run
129     runinfo = new ItoRunInfo(f2);
130     char *trkversion = TrkInfo(false);
131     int runinfo_error = runinfo->Update(run,"TRK",trkversion);
132     if( runinfo_error ) throw runinfo_error;
133     runlist = runinfo->GetRunList();
134     reprocessing = runinfo->IsReprocessing();//????
135 pam-fi 1.25 if(TrkParams::VerboseMode()){
136 pam-fi 1.13 cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file
137     cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file
138     cout << "1-st run "<< runlist->At(0) << endl;
139     cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl;
140     cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run
141     cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run
142     cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl;
143     };
144     };
145     //-------------------------------------------
146     //
147     //-------------------------------------------
148     // Take (if present) the old Tracker tree
149     t_clone = (TTree*)f2->Get("Tracker");
150     if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone);
151     // Create NEW Tracker tree
152     t_level2 = new TTree("Tracker","PAMELA tracker level2 data ");
153 pam-fi 1.21 l2_event->Set(); // ****NBNBNBN*****
154 pam-fi 1.13 t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
155     if(p->get1){
156 pam-fi 1.25 if(TrkParams::VerboseMode())cout << endl << "Requested LEVEL1 output" << endl;
157 pam-fi 1.21 l1_event->Set(); // ****NBNBNBN*****
158 pam-fi 1.13 t_level2->Branch("TrkLevel1","TrkLevel1",&l1_event);
159 pam-fi 1.16 }
160     if(p->geth){
161 pam-fi 1.25 if(TrkParams::VerboseMode())cout << endl << "Requested Hough-Transform output" << endl;
162 pam-fi 1.16 t_level2->Branch("TrkHough","TrkHough",&lh_event);
163 pam-fi 1.7 };
164 pam-fi 1.13 };
165 mocchiut 1.1
166    
167     // -------------------------------------------
168     // define runs to be processed/reprocessed
169     // -------------------------------------------
170     if(run == 0){
171 pam-fi 1.13 // reprocessing ALL runs
172     if(p->standalone)throw -298; // reprocessing not implemented
173     from_run = runlist->At(0);
174     to_run = runlist->At(runinfo->GetRunEntries()-1);
175 mocchiut 1.1 }else{
176     // processing/reprocessing ONE single run
177 pam-fi 1.13 from_run = run;
178     to_run = run;
179 mocchiut 1.1 };
180    
181     //
182     // init event counter
183     //
184     Int_t ev_count =0;
185     //
186     // create query-results objects
187     //
188 pam-fi 1.7 GL_RUN q1 = GL_RUN();
189 mocchiut 1.1 GL_TRK_CALIB q2 = GL_TRK_CALIB();
190     GL_ROOT q3 = GL_ROOT();
191     GL_PARAM q4 = GL_PARAM();
192    
193     // ------------------------------------------------------------
194     // if reprocessing one run, copy all the events BEFORE the run
195     // ------------------------------------------------------------
196     if( !(p->standalone) ){
197 pam-fi 1.13 for(UInt_t i=0; i<runinfo->GetFirstEntry(); i++){
198 mocchiut 1.37 if ( t_clone->GetEntry(i) <= 0 ) throw -36;// EM
199 pam-fi 1.13 *l2_event = *l2_clone;
200     t_level2->Fill();
201     l2_event->Clear();
202     // COPY COPY COPY
203     };
204 mocchiut 1.1 };
205     // ------------------------------------------------------------
206     // ------------------------------------------------------------
207 pam-fi 1.25 // START LOOP OVER RUNS TO PROCESS/REPROCESS
208 mocchiut 1.1 // ------------------------------------------------------------
209     // ------------------------------------------------------------
210 pam-fi 1.7 for(UInt_t idRun = from_run; idRun <= to_run; idRun++){
211 mocchiut 1.1
212 pam-fi 1.25 if(TrkParams::VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
213 pam-fi 1.7 UInt_t runheadtime = 0;
214     UInt_t runtrailtime = 0;
215 pam-fi 1.32 UInt_t runheadobt = 0;
216     UInt_t runtrailobt = 0;
217 mocchiut 1.1 UInt_t evfrom = 0;
218     UInt_t evto = 0;
219 pam-fi 1.7 UInt_t nevents = 0;
220     Int_t id_root_l0 =-1;
221 mocchiut 1.1 Int_t trk_calib_used = 0;
222 pam-fi 1.7
223 mocchiut 1.31 TString host = glt->CGetHost(); // EMILIANO
224     TString user = glt->CGetUser(); // EMILIANO
225     TString psw = glt->CGetPsw(); // EMILIANO
226     // TString host = "mysql://localhost/pamelaprod";
227     // TString user = "anonymous";
228     // TString psw = "";
229     // const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
230     // const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
231     // const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
232     // if ( !pamdbhost ) pamdbhost = "";
233     // if ( !pamdbuser ) pamdbuser = "";
234     // if ( !pamdbpsw ) pamdbpsw = "";
235     // if ( strcmp(pamdbhost,"") ) host = pamdbhost;
236     // if ( strcmp(pamdbuser,"") ) user = pamdbuser;
237     // if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
238 mocchiut 1.29 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
239 mocchiut 1.31 stringstream myquery; // EMILIANO
240     myquery.str(""); // EMILIANO
241     myquery << "SET time_zone='+0:00'"; // EMILIANO
242     dbc->Query(myquery.str().c_str()); // EMILIANO
243 mocchiut 1.1 if(p->standalone){
244     // ==============================
245     // first query: retrieve run info
246     // ==============================
247 pam-fi 1.7 if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
248     id_root_l0 = q1.ID_ROOT_L0;
249 mocchiut 1.1 runheadtime = q1.RUNHEADER_TIME;
250     runtrailtime = q1.RUNTRAILER_TIME;
251 pam-fi 1.32 runheadobt = q1.RUNHEADER_OBT;
252     runtrailobt = q1.RUNTRAILER_OBT;
253 pam-fi 1.7 evfrom = q1.EV_FROM;
254     evto = q1.EV_TO;
255 pam-fi 1.13 nevents = q1.NEVENTS;
256     trk_calib_used = q1.TRK_CALIB_USED;
257 mocchiut 1.1 }else{
258     // ==============================
259     // get run info from RunInfo tree
260     // ==============================
261 pam-fi 1.3 int runinfo_error = runinfo->GetRunInfo(idRun);
262 pam-fi 1.7 if( runinfo_error ) throw runinfo_error;
263     id_root_l0 = runinfo->ID_ROOT_L0;
264 mocchiut 1.1 runheadtime = runinfo->RUNHEADER_TIME;
265     runtrailtime = runinfo->RUNTRAILER_TIME;
266 pam-fi 1.32 runheadobt = runinfo->RUNHEADER_OBT;
267     runtrailobt = runinfo->RUNTRAILER_OBT;
268 pam-fi 1.7 evfrom = runinfo->EV_FROM;
269     evto = runinfo->EV_TO;
270 pam-fi 1.13 nevents = runinfo->NEVENTS;
271     trk_calib_used = runinfo->TRK_CALIB_USED;
272 pam-fi 1.24
273 mocchiut 1.1 };
274 pam-fi 1.24
275 mocchiut 1.1 //
276 pam-fi 1.25 if(TrkParams::VerboseMode()){
277 pam-fi 1.7 cout << "ROOT file ID "<< id_root_l0 << endl;
278 mocchiut 1.1 cout << "RunHeader time "<< runheadtime << endl;
279     cout << "RunTrailer time "<< runtrailtime << endl;
280 pam-fi 1.7 cout << " from event "<< evfrom << endl;
281     cout << " to event "<< evto << endl;
282 pam-fi 1.24 cout << " num. events "<< nevents << endl;
283     cout << "trk-calibration used "<< trk_calib_used << endl;
284 mocchiut 1.1 };
285     // ========================================================
286     // second query: search the LEVEL0 file that contains idRun
287     // ========================================================
288     TString lastfilename = filename;
289 pam-fi 1.7 if( q3.Query_GL_ROOT(id_root_l0,dbc) )throw -51;
290 mocchiut 1.1 filename = q3.PATH + q3.NAME;
291     // ========================================================
292     // Open the input LEVEL0 data file
293     // ========================================================
294     if(filename.CompareTo(lastfilename)){
295     if(!lastfilename.IsNull())f0->Close();
296     //if( debug ) cout << "Opening LEVEL0 file: "<< filename << endl;
297 pam-fi 1.25 if(TrkParams::VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
298 pam-fi 1.2 FileStat_t t;
299     if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
300 mocchiut 1.1 f0 = new TFile(filename);
301     if ( !f0 ) throw -6;
302     physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7;
303     b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8;
304     b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
305 pam-fi 1.22 l0_event->Set();
306     physicsTree->SetBranchAddress("Tracker" ,l0_event->GetPointerToTrackerEvent());
307 mocchiut 1.1 physicsTree->SetBranchAddress("Header",&header);
308    
309 pam-fi 1.7 nentries = physicsTree->GetEntries();
310 mocchiut 1.34 if ( nentries < 1 && nevents ) throw -11; // EMILIANO: go on if the file is empty, why not? in principle we should not have any event to process but if the case the next line will take care of exiting; exit only if the file is empty but nevents is not zero
311 mocchiut 1.33 if ( nentries < (evto+1) && nevents > 0 ) throw -12; // EMILIANO: if NEVENTS = 0 and the file is empty everything is ok but due to a mistake in the
312     // variable EV_TO type (UInt_t instead of Int_t) that we don't want to change to avoid changing a lot of code, evto becomes inf hence
313     // the condition ( nentries < (evto+1) ) is satisfied and DV exit with error even if the error is only in the DB.
314 pam-fi 1.7
315 mocchiut 1.1 };
316    
317 pam-fi 1.7 GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc);
318    
319 pam-fi 1.24
320     // =============================================================
321     // retrieve information about parameters to process LEVEL2
322     // =============================================================
323    
324     TrkParams::Set(runinfo->GetGL_RUN(),dbc);
325 pam-fi 1.25 for(int i=0; i<p->npar; i++){
326     if(TrkParams::VerboseMode())cout<<" ((( force parameters from input path )))"<<endl;
327     TrkParams::Set(p->parpath[i],p->partype[i]);
328     }
329    
330 pam-fi 1.24 TrkParams::Load();
331     if( !TrkParams::IsLoaded() )throw -52;
332    
333     // =============================================================
334     // retrieve calibration file needed to reduce data
335     // =============================================================
336    
337     TrkParams::SetCalib(runinfo->GetGL_RUN(),dbc);
338     TrkParams::LoadCalib( );
339     if( !TrkParams::CalibIsLoaded() )throw -53;
340    
341 pam-fi 1.11 TBenchmark *reduction = new TBenchmark();
342 pam-fi 1.25 if(TrkParams::VerboseMode())reduction->Start("reduction");
343 pam-fi 1.11 Int_t ntrk=0;
344 mocchiut 1.1 // ====================================================
345     // start looping on events cointained in the data file
346     // ====================================================
347 mocchiut 1.28 if(dbc){
348 mocchiut 1.29 dbc->Close();
349     delete dbc;};
350 mocchiut 1.28 //
351 pam-fi 1.32 for (UInt_t re = evfrom+min(p->nskip,nevents); re < evfrom+nevents; re++){
352 mocchiut 1.1
353     ev_count++;
354 pam-fi 1.32
355    
356 pam-fi 1.25 // if ( TrkParams::DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
357 mocchiut 1.1
358 mocchiut 1.37 if ( b_trk->GetEntry(re) <= 0 ) throw -36;//EM
359     if ( b_header->GetEntry(re) <= 0 ) throw -36;//EM
360 pam-fi 1.13 pscu = header->GetPscuHeader();
361 pam-fi 1.23
362 pam-fi 1.25 if( TrkParams::DebugMode() )cout << ">>> "<<ev_count-1<<" @ OBT "<<pscu->GetOrbitalTime()<<endl;
363 pam-fi 1.23
364 mocchiut 1.35 if ( dbtime->DBabsTime(pscu->GetOrbitalTime()) > (runtrailtime+1) || dbtime->DBabsTime(pscu->GetOrbitalTime()) < (runheadtime-1)) {
365 pam-fi 1.7
366 pam-fi 1.25 if (TrkParams::VerboseMode()){
367 pam-fi 1.14 printf(" TrkCore - WARNING: event outside the run time window, skipping it\n");
368 pam-fi 1.32 cout << " OBT "<<pscu->GetOrbitalTime()<<" RUN "<<runheadobt<<"-"<<runtrailobt<<" ABS-time "<<dbtime->DBabsTime(pscu->GetOrbitalTime())<<" RUN "<<runheadtime<<"-"<<runtrailtime<<endl;
369 pam-fi 1.14 };
370 pam-fi 1.13 }else{
371 pam-fi 1.25 if ( TrkParams::DebugMode() )
372 pam-fi 1.23 printf("\n-------------------------------\nevent %d\n",re-evfrom);
373 pam-fi 1.16
374 pam-fi 1.13 p->ProcessEvent(l0_event);
375 pam-fi 1.7
376 pam-fi 1.13 // ----------------
377     // LEVEL1 output
378     // ----------------
379     if(p->get1){
380     if(p->ifroot1){ // root
381     l1_event->Clear();
382 pam-fi 1.22 // l1_event->SetFromLevel1Struct(&level1event_,p->full1);
383     l1_event->SetFromLevel1Struct(p->full1);
384 pam-fi 1.13 // t_level1->Fill();
385 pam-fi 1.23 }else{ // hbook
386 pam-fi 1.13 throw -299;
387     };
388     };
389     // ----------------
390     // HOUGH output
391     // ----------------
392     if(p->geth){
393     if(p->ifrooth){ // root
394 pam-fi 1.23 lh_event->Delete();
395     lh_event->SetFromHoughStruct(&houghevent_);
396 pam-fi 1.13 }else{ // hbook
397     throw -299;
398     };
399     };
400     // ----------------
401     // LEVEL2 output
402     // ----------------
403     if(p->get2){
404 pam-fi 1.23 l2_event->Clear();
405     if(p->get1){
406     l2_event->SetFromLevel2Struct(&level2event_,l1_event);//set references to level1
407     }else{
408     l2_event->SetFromLevel2Struct(&level2event_);
409     }
410 pam-fi 1.10 // l2_event->Dump();
411 pam-fi 1.23 t_level2->Fill();
412     if( l2_event->ntrk()>0 )ntrk++;
413 pam-fi 1.27 if(TrkParams::VerboseMode())l2_event->Dump();
414 mocchiut 1.1 };
415 pam-fi 1.13 };
416 mocchiut 1.1 }; // end loop on events
417 pam-fi 1.25 if(TrkParams::VerboseMode()){
418 pam-fi 1.13 cout << " >>> processed "<< ev_count <<" events"<< endl;
419     if(ev_count)cout << ntrk << " events with at least one track ("<<(Int_t)(100*ntrk/ev_count)<<"%)\n";
420     reduction->Show("reduction");
421 pam-fi 1.11 }
422     delete reduction;
423 pam-fi 1.7
424     delete dbtime;
425    
426 pam-fi 1.13 }; // end loop on runs
427    
428    
429 mocchiut 1.1 // ------------------------------------------------------------
430     // if reprocessing one run, copy all the events AFTER the run
431     // ------------------------------------------------------------
432     if( !(p->standalone) ){
433     for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
434 mocchiut 1.37 if ( t_clone->GetEntry(i) <= 0 ) throw -36;//EM
435 mocchiut 1.1 *l2_event = *l2_clone;
436     t_level2->Fill();
437     l2_event->Clear();
438     // COPY COPY COPY
439     };
440     };
441     // ---------------
442     // close the files
443     // ---------------
444     if(p->get2){
445 mocchiut 1.38 // if( t_clone )t_clone->Delete("all");//delete old tree from file
446 pam-fi 1.13 if( !(p->standalone) )runinfo->Close();
447 mocchiut 1.38 t_level2->Write("Tracker", TObject::kOverwrite);
448     // if( t_level2 )t_level2->Delete(); //delete new tree from memory
449 pam-fi 1.13
450     };
451 pam-fi 1.6
452 pam-fi 1.24 if(f0) if(f0->IsOpen()) f0->Close();
453     // if( f0_c->IsOpen() )f0_c->Close();
454 pam-fi 1.13
455     // lh_event->Delete();
456     l1_event->Delete();
457     l2_event->Delete();
458     l2_clone->Delete();
459 pam-fi 1.2
460     return(p->ostatus);
461 mocchiut 1.1 }
462    

  ViewVC Help
Powered by ViewVC 1.1.23