/[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.32 - (hide annotations) (download)
Wed Mar 5 17:00:19 2008 UTC (16 years, 9 months ago) by pam-fi
Branch: MAIN
CVS Tags: v5r00
Changes since 1.31: +10 -3 lines
modified TrkSinglet, optimized DoTrack2, fixed bug in evaluation of effective angle

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     TString filename = 0;
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     t_clone->GetEntry(i);
199     *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     if ( nentries < 1 ) throw -11;
311     if ( nentries < (evto+1)) throw -12;
312    
313 mocchiut 1.1 };
314    
315 pam-fi 1.7 GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc);
316    
317 pam-fi 1.24
318     // =============================================================
319     // retrieve information about parameters to process LEVEL2
320     // =============================================================
321    
322     TrkParams::Set(runinfo->GetGL_RUN(),dbc);
323 pam-fi 1.25 for(int i=0; i<p->npar; i++){
324     if(TrkParams::VerboseMode())cout<<" ((( force parameters from input path )))"<<endl;
325     TrkParams::Set(p->parpath[i],p->partype[i]);
326     }
327    
328 pam-fi 1.24 TrkParams::Load();
329     if( !TrkParams::IsLoaded() )throw -52;
330    
331     // =============================================================
332     // retrieve calibration file needed to reduce data
333     // =============================================================
334    
335     TrkParams::SetCalib(runinfo->GetGL_RUN(),dbc);
336     TrkParams::LoadCalib( );
337     if( !TrkParams::CalibIsLoaded() )throw -53;
338    
339 pam-fi 1.11 TBenchmark *reduction = new TBenchmark();
340 pam-fi 1.25 if(TrkParams::VerboseMode())reduction->Start("reduction");
341 pam-fi 1.11 Int_t ntrk=0;
342 mocchiut 1.1 // ====================================================
343     // start looping on events cointained in the data file
344     // ====================================================
345 mocchiut 1.28 if(dbc){
346 mocchiut 1.29 dbc->Close();
347     delete dbc;};
348 mocchiut 1.28 //
349 pam-fi 1.32 for (UInt_t re = evfrom+min(p->nskip,nevents); re < evfrom+nevents; re++){
350 mocchiut 1.1
351     ev_count++;
352 pam-fi 1.32
353    
354 pam-fi 1.25 // if ( TrkParams::DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
355 mocchiut 1.1
356 pam-fi 1.13 b_trk->GetEntry(re);
357     b_header->GetEntry(re);
358     pscu = header->GetPscuHeader();
359 pam-fi 1.23
360 pam-fi 1.25 if( TrkParams::DebugMode() )cout << ">>> "<<ev_count-1<<" @ OBT "<<pscu->GetOrbitalTime()<<endl;
361 pam-fi 1.23
362 pam-fi 1.13 if ( dbtime->DBabsTime(pscu->GetOrbitalTime()) > runtrailtime || dbtime->DBabsTime(pscu->GetOrbitalTime()) < runheadtime) {
363 pam-fi 1.7
364 pam-fi 1.25 if (TrkParams::VerboseMode()){
365 pam-fi 1.14 printf(" TrkCore - WARNING: event outside the run time window, skipping it\n");
366 pam-fi 1.32 cout << " OBT "<<pscu->GetOrbitalTime()<<" RUN "<<runheadobt<<"-"<<runtrailobt<<" ABS-time "<<dbtime->DBabsTime(pscu->GetOrbitalTime())<<" RUN "<<runheadtime<<"-"<<runtrailtime<<endl;
367 pam-fi 1.14 };
368 pam-fi 1.13 }else{
369 pam-fi 1.25 if ( TrkParams::DebugMode() )
370 pam-fi 1.23 printf("\n-------------------------------\nevent %d\n",re-evfrom);
371 pam-fi 1.16
372 pam-fi 1.13 p->ProcessEvent(l0_event);
373 pam-fi 1.7
374 pam-fi 1.13 // ----------------
375     // LEVEL1 output
376     // ----------------
377     if(p->get1){
378     if(p->ifroot1){ // root
379     l1_event->Clear();
380 pam-fi 1.22 // l1_event->SetFromLevel1Struct(&level1event_,p->full1);
381     l1_event->SetFromLevel1Struct(p->full1);
382 pam-fi 1.13 // t_level1->Fill();
383 pam-fi 1.23 }else{ // hbook
384 pam-fi 1.13 throw -299;
385     };
386     };
387     // ----------------
388     // HOUGH output
389     // ----------------
390     if(p->geth){
391     if(p->ifrooth){ // root
392 pam-fi 1.23 lh_event->Delete();
393     lh_event->SetFromHoughStruct(&houghevent_);
394 pam-fi 1.13 }else{ // hbook
395     throw -299;
396     };
397     };
398     // ----------------
399     // LEVEL2 output
400     // ----------------
401     if(p->get2){
402 pam-fi 1.23 l2_event->Clear();
403     if(p->get1){
404     l2_event->SetFromLevel2Struct(&level2event_,l1_event);//set references to level1
405     }else{
406     l2_event->SetFromLevel2Struct(&level2event_);
407     }
408 pam-fi 1.10 // l2_event->Dump();
409 pam-fi 1.23 t_level2->Fill();
410     if( l2_event->ntrk()>0 )ntrk++;
411 pam-fi 1.27 if(TrkParams::VerboseMode())l2_event->Dump();
412 mocchiut 1.1 };
413 pam-fi 1.13 };
414 mocchiut 1.1 }; // end loop on events
415 pam-fi 1.25 if(TrkParams::VerboseMode()){
416 pam-fi 1.13 cout << " >>> processed "<< ev_count <<" events"<< endl;
417     if(ev_count)cout << ntrk << " events with at least one track ("<<(Int_t)(100*ntrk/ev_count)<<"%)\n";
418     reduction->Show("reduction");
419 pam-fi 1.11 }
420     delete reduction;
421 pam-fi 1.7
422     delete dbtime;
423    
424 pam-fi 1.13 }; // end loop on runs
425    
426    
427 mocchiut 1.1 // ------------------------------------------------------------
428     // if reprocessing one run, copy all the events AFTER the run
429     // ------------------------------------------------------------
430     if( !(p->standalone) ){
431     for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
432     t_clone->GetEntry(i);
433     *l2_event = *l2_clone;
434     t_level2->Fill();
435     l2_event->Clear();
436     // COPY COPY COPY
437     };
438     };
439     // ---------------
440     // close the files
441     // ---------------
442     if(p->get2){
443 pam-fi 1.13 if( t_clone )t_clone->Delete("all");//delete old tree from file
444     if( !(p->standalone) )runinfo->Close();
445     f2->Write("Tracker");
446     if( t_level2 )t_level2->Delete(); //delete new tree from memory
447    
448     };
449 pam-fi 1.6
450 pam-fi 1.24 if(f0) if(f0->IsOpen()) f0->Close();
451     // if( f0_c->IsOpen() )f0_c->Close();
452 pam-fi 1.13
453     // lh_event->Delete();
454     l1_event->Delete();
455     l2_event->Delete();
456     l2_clone->Delete();
457 pam-fi 1.2
458     return(p->ostatus);
459 mocchiut 1.1 }
460    

  ViewVC Help
Powered by ViewVC 1.1.23