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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.24 - (show annotations) (download)
Fri Apr 27 10:39:57 2007 UTC (17 years, 7 months ago) by pam-fi
Branch: MAIN
CVS Tags: v3r04, v3r05, v3r06, v3r03
Changes since 1.23: +139 -133 lines
v3r00: new hough parameters, new variables, and other things...

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 #include <TBenchmark.h>
27
28 // .........................................................
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 //#include <TrkStruct.h>
41 #include <TrkLevel2.h>
42 #include <TrkLevel1.h>
43 #include <TrkLevel0.h>
44 #include <TrkHough.h>
45 // .........................................................
46 // YODA header files
47 // .........................................................
48 #include <PamelaRun.h>
49 //#include <RegistryEvent.h>
50 #include <physics/tracker/TrackerEvent.h>
51 #include <CalibTrk1Event.h>
52 #include <CalibTrk2Event.h>
53
54 // .........................................................
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 int TrkCore(UInt_t run, TFile *f2, TSQLServer *dbc, int ncustom, char*vcustom[])
76 {
77 // ---------------------------
78 // Define some basic variables
79 // ---------------------------
80
81 TString filename = 0;
82 // TString trkcalibfile = 0;
83 // Int_t last_trk_calib_used = 0;
84 Long64_t nentries = 0LL;
85
86 // LEVEL0 input
87 TFile *f0 = 0;
88 // TFile *f0_c = 0;
89 TTree *physicsTree = 0;
90 TBranch *b_trk = 0;
91 TBranch *b_header = 0;
92 EventHeader *header = 0;
93 PscuHeader *pscu = 0;
94 TrkLevel0 *l0_event = new TrkLevel0();
95 // RunInfo
96 ItoRunInfo *runinfo = 0;
97 TArrayI *runlist = 0;
98 UInt_t from_run;
99 UInt_t to_run;
100 Bool_t reprocessing = false;
101 //LEVEL2 output - ROOT
102 TTree *t_level2 = 0;
103 TTree *t_clone = 0;
104 TrkLevel2 *l2_event = new TrkLevel2();
105 TrkLevel2 *l2_clone = new TrkLevel2();
106 //LEVEL1 output - ROOT
107 TrkLevel1 *l1_event = new TrkLevel1();
108 TrkHough *lh_event = new TrkHough();
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 if( p->DebugMode() )p->Dump();
121
122 // ===================================
123 // Open/Create level2 output file
124 // ===================================
125 if(p->get2){
126 //-------------------------------------------
127 // read RunInfo
128 //-------------------------------------------
129 if(!(p->standalone)){
130 // Open "RunInfo" tree and get list of run
131 runinfo = new ItoRunInfo(f2);
132 char *trkversion = TrkInfo(false);
133 int runinfo_error = runinfo->Update(run,"TRK",trkversion);
134 if( runinfo_error ) throw runinfo_error;
135 runlist = runinfo->GetRunList();
136 reprocessing = runinfo->IsReprocessing();//????
137 if(p->VerboseMode()){
138 cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file
139 cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file
140 cout << "1-st run "<< runlist->At(0) << endl;
141 cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl;
142 cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run
143 cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run
144 cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl;
145 };
146 };
147 //-------------------------------------------
148 //
149 //-------------------------------------------
150 // Take (if present) the old Tracker tree
151 t_clone = (TTree*)f2->Get("Tracker");
152 if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone);
153 // Create NEW Tracker tree
154 t_level2 = new TTree("Tracker","PAMELA tracker level2 data ");
155 l2_event->Set(); // ****NBNBNBN*****
156 t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
157 if(p->get1){
158 if(p->VerboseMode())cout << endl << "Requested LEVEL1 output" << endl;
159 l1_event->Set(); // ****NBNBNBN*****
160 t_level2->Branch("TrkLevel1","TrkLevel1",&l1_event);
161 }
162 if(p->geth){
163 if(p->VerboseMode())cout << endl << "Requested Hough-Transform output" << endl;
164 t_level2->Branch("TrkHough","TrkHough",&lh_event);
165 };
166 };
167
168
169 // -------------------------------------------
170 // define runs to be processed/reprocessed
171 // -------------------------------------------
172 if(run == 0){
173 // reprocessing ALL runs
174 if(p->standalone)throw -298; // reprocessing not implemented
175 from_run = runlist->At(0);
176 to_run = runlist->At(runinfo->GetRunEntries()-1);
177 }else{
178 // processing/reprocessing ONE single run
179 from_run = run;
180 to_run = run;
181 };
182
183 //
184 // init "expiration date" of calibrations and parameters
185 //
186 // UInt_t pedsig_time =0;
187 // UInt_t B_time =0;
188 // UInt_t mip_time =0;
189 // UInt_t charge_time =0;
190 // UInt_t eta_time =0;
191 // UInt_t align_time =0;
192 // UInt_t mask_time =0;
193 //
194 // init event counter
195 //
196 Int_t ev_count =0;
197 //
198 // create query-results objects
199 //
200 GL_RUN q1 = GL_RUN();
201 GL_TRK_CALIB q2 = GL_TRK_CALIB();
202 GL_ROOT q3 = GL_ROOT();
203 GL_PARAM q4 = GL_PARAM();
204
205 // ------------------------------------------------------------
206 // if reprocessing one run, copy all the events BEFORE the run
207 // ------------------------------------------------------------
208 if( !(p->standalone) ){
209 for(UInt_t i=0; i<runinfo->GetFirstEntry(); i++){
210 t_clone->GetEntry(i);
211 *l2_event = *l2_clone;
212 t_level2->Fill();
213 l2_event->Clear();
214 // COPY COPY COPY
215 };
216 };
217 // ------------------------------------------------------------
218 // ------------------------------------------------------------
219 // START LOOP OVER RUNS TO PROCESSED/REPROCESSED
220 // ------------------------------------------------------------
221 // ------------------------------------------------------------
222 for(UInt_t idRun = from_run; idRun <= to_run; idRun++){
223
224 if(p->VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
225 UInt_t runheadtime = 0;
226 UInt_t runtrailtime = 0;
227 UInt_t evfrom = 0;
228 UInt_t evto = 0;
229 UInt_t nevents = 0;
230 Int_t id_root_l0 =-1;
231 Int_t trk_calib_used = 0;
232
233 if(p->standalone){
234 // ==============================
235 // first query: retrieve run info
236 // ==============================
237 if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
238 id_root_l0 = q1.ID_ROOT_L0;
239 runheadtime = q1.RUNHEADER_TIME;
240 runtrailtime = q1.RUNTRAILER_TIME;
241 evfrom = q1.EV_FROM;
242 evto = q1.EV_TO;
243 nevents = q1.NEVENTS;
244 trk_calib_used = q1.TRK_CALIB_USED;
245 }else{
246 // ==============================
247 // get run info from RunInfo tree
248 // ==============================
249 int runinfo_error = runinfo->GetRunInfo(idRun);
250 if( runinfo_error ) throw runinfo_error;
251 id_root_l0 = runinfo->ID_ROOT_L0;
252 runheadtime = runinfo->RUNHEADER_TIME;
253 runtrailtime = runinfo->RUNTRAILER_TIME;
254 evfrom = runinfo->EV_FROM;
255 evto = runinfo->EV_TO;
256 nevents = runinfo->NEVENTS;
257 trk_calib_used = runinfo->TRK_CALIB_USED;
258
259 };
260
261 //
262 if(p->VerboseMode()){
263 cout << "ROOT file ID "<< id_root_l0 << endl;
264 cout << "RunHeader time "<< runheadtime << endl;
265 cout << "RunTrailer time "<< runtrailtime << endl;
266 cout << " from event "<< evfrom << endl;
267 cout << " to event "<< evto << endl;
268 cout << " num. events "<< nevents << 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_root_l0,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 if(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
284 FileStat_t t;
285 if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
286 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_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
291 l0_event->Set();
292 physicsTree->SetBranchAddress("Tracker" ,l0_event->GetPointerToTrackerEvent());
293 physicsTree->SetBranchAddress("Header",&header);
294
295 nentries = physicsTree->GetEntries();
296 if ( nentries < 1 ) throw -11;
297 if ( nentries < (evto+1)) throw -12;
298
299 };
300
301 GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc);
302
303
304 // =============================================================
305 // retrieve information about parameters to process LEVEL2
306 // =============================================================
307
308 TrkParams::Set(runinfo->GetGL_RUN(),dbc);
309 TrkParams::Load();
310 if( !TrkParams::IsLoaded() )throw -52;
311
312 // =============================================================
313 // retrieve calibration file needed to reduce data
314 // =============================================================
315
316 TrkParams::SetCalib(runinfo->GetGL_RUN(),dbc);
317 TrkParams::LoadCalib( );
318 if( !TrkParams::CalibIsLoaded() )throw -53;
319
320 // =============================================================
321 // retrieve calibration file needed to reduce data
322 // =============================================================
323 // if run OBT is > last calibration "expiration date"
324 // - search for new calibration packet
325 // - load calibration parameters (full + truncated)
326 // =============================================================
327 // if(p->VerboseMode() )cout << "Full pedestals for cluster finding:";
328 // if(runheadtime > pedsig_time){
329 // if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
330 // pedsig_time = q2.TO_TIME;
331 // if(q2.EV_ROOT_CALIBTRK1 != q2.EV_ROOT_CALIBTRK2)
332 // printf("WARNING!! ---> EV_ROOT_CALIBTRK1=%d it's different from EV_ROOT_CALIBTRK2=%d \n\n",q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
333 // if( q3.Query_GL_ROOT(q2.ID_ROOT_L0,dbc) )throw -51;
334 // trkcalibfile = q3.PATH + q3.NAME;
335 // if(p->VerboseMode()) cout << " >> Loading from LEVEL0 file: "<< trkcalibfile << endl;
336
337 // if(trkcalibfile.CompareTo(filename)){
338 // FileStat_t t;
339 // if( gSystem->GetPathInfo(trkcalibfile.Data(),t) )throw -6;
340 // f0_c=new TFile(trkcalibfile);
341 // if ( !f0_c ) throw -6;
342 // }else{
343 // f0_c = f0;
344 // };
345 // if(p->VerboseMode()) {
346 // cout << " calibration entries "<< q2.EV_ROOT_CALIBTRK1 << " " << q2.EV_ROOT_CALIBTRK2;
347 // cout << " (from time "<< q2.FROM_TIME <<" to time "<< q2.TO_TIME <<")"<<endl;
348 // };
349 // TrkParams::FillACalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
350 // TrkParams::FillMask(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
351
352 // }else{
353 // if(p->VerboseMode() )cout<< " already loaded. "<< endl;
354 // };
355 // if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal[64][12][i] << endl;
356 // =============================================================
357 // retrieve calibration file needed to uncompress data
358 // =============================================================
359 // if the run was compressed using default calib
360 // load truncated pedestals from default
361 // otherwise reload them from on-line calibration
362 // =============================================================
363 // if(p->VerboseMode() )cout << "Truncated pedestals for uncompression:";
364 // if(trk_calib_used==104 && last_trk_calib_used!=104){
365
366 // if(p->VerboseMode() )cout << " >> Loading default calibration " << endl;
367 // if (q4.Query_GL_PARAM(runheadtime,7,dbc) )throw -52;
368 // TrkParams::FillTCalibFrom(q4.PATH+q4.NAME);
369
370 // }else if( (trk_calib_used==1 | trk_calib_used==2) && last_trk_calib_used==104 ){
371
372 // if(p->VerboseMode()) cout << ">> Re-loading on-line calibration " << endl;
373 // TrkParams::FillTCalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
374 // }else{
375 // if(p->VerboseMode()) cout << ">> already loaded " << endl;
376 // };
377 // if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal_t[64][12][i] << endl;
378 // last_trk_calib_used = trk_calib_used;
379
380
381
382 // //
383 // // magnetic field
384 // //
385 // if(p->VerboseMode()) cout << "Loading parameter files : "<< endl;
386
387 // if(runheadtime > B_time){
388 // if( q4.Query_GL_PARAM(runheadtime,1,dbc) )throw -52;
389 // B_time = q4.TO_TIME;
390 // l2_event->LoadField(q4.PATH+q4.NAME);
391 // if(path_.error) throw -215;
392 // };
393
394 // //
395 // // mip conversion parameters
396 // //
397 // if(runheadtime > mip_time){
398 // if( q4.Query_GL_PARAM(runheadtime,2,dbc) )throw -52;
399 // mip_time = q4.TO_TIME;
400 // path_.FillWith(q4.PATH+q4.NAME);
401 // readmipparam_();
402 // if(path_.error) throw -212;
403 // };
404 // //
405 // // charge correlation parameters
406 // //
407 // if(runheadtime > charge_time){
408 // if( q4.Query_GL_PARAM(runheadtime,3,dbc) )throw -52;
409 // charge_time = q4.TO_TIME;
410 // path_.FillWith(q4.PATH+q4.NAME);
411 // readchargeparam_();
412 // if(path_.error) throw -213;
413 // };
414 // //
415 // // eta p.f.a. parameters
416 // //
417 // if(runheadtime > eta_time){
418 // if( q4.Query_GL_PARAM(runheadtime,4,dbc) )throw -52;
419 // eta_time = q4.TO_TIME;
420 // path_.FillWith(q4.PATH+q4.NAME);
421 // readetaparam_();
422 // if(path_.error) throw -214;
423 // };
424 // //
425 // // alignment parameters
426 // //
427 // if(runheadtime > align_time){
428 // if( q4.Query_GL_PARAM(runheadtime,5,dbc) )throw -52;
429 // align_time = q4.TO_TIME;
430 // path_.FillWith(q4.PATH+q4.NAME);
431 // readalignparam_();
432 // if(path_.error) throw -211;
433 // };
434 // //
435 // // viking mask
436 // //
437 // if(runheadtime > mask_time){
438 // if( q4.Query_GL_PARAM(runheadtime,6,dbc) )throw -52;
439 // mask_time = q4.TO_TIME;
440 // path_.FillWith(q4.PATH+q4.NAME);
441 // readvkmask_();
442 // if(path_.error) throw -210;
443 // };
444
445 TBenchmark *reduction = new TBenchmark();
446 if(p->VerboseMode())reduction->Start("reduction");
447 Int_t ntrk=0;
448 // ====================================================
449 // start looping on events cointained in the data file
450 // ====================================================
451 for (UInt_t re = evfrom; re < evfrom+nevents; re++){
452
453 ev_count++;
454
455 // if ( p->DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
456
457 b_trk->GetEntry(re);
458 b_header->GetEntry(re);
459 pscu = header->GetPscuHeader();
460
461 if( p->DebugMode() )cout << ">>> "<<ev_count-1<<" @ OBT "<<pscu->GetOrbitalTime()<<endl;
462
463 if ( dbtime->DBabsTime(pscu->GetOrbitalTime()) > runtrailtime || dbtime->DBabsTime(pscu->GetOrbitalTime()) < runheadtime) {
464
465 if (p->VerboseMode()){
466 printf(" TrkCore - WARNING: event outside the run time window, skipping it\n");
467 cout << " OBT "<<pscu->GetOrbitalTime()<<" ABS-time "<<dbtime->DBabsTime(pscu->GetOrbitalTime())<<" RUN "<<runheadtime<<"-"<<runtrailtime<<endl;
468 };
469 }else{
470 if ( p->DebugMode() )
471 printf("\n-------------------------------\nevent %d\n",re-evfrom);
472
473 p->ProcessEvent(l0_event);
474
475 // ----------------
476 // LEVEL1 output
477 // ----------------
478 if(p->get1){
479 if(p->ifroot1){ // root
480 l1_event->Clear();
481 // l1_event->SetFromLevel1Struct(&level1event_,p->full1);
482 l1_event->SetFromLevel1Struct(p->full1);
483 // t_level1->Fill();
484 }else{ // hbook
485 throw -299;
486 };
487 };
488 // ----------------
489 // HOUGH output
490 // ----------------
491 if(p->geth){
492 if(p->ifrooth){ // root
493 lh_event->Delete();
494 lh_event->SetFromHoughStruct(&houghevent_);
495 }else{ // hbook
496 throw -299;
497 };
498 };
499 // ----------------
500 // LEVEL2 output
501 // ----------------
502 if(p->get2){
503 l2_event->Clear();
504 if(p->get1){
505 l2_event->SetFromLevel2Struct(&level2event_,l1_event);//set references to level1
506 }else{
507 l2_event->SetFromLevel2Struct(&level2event_);
508 }
509 // l2_event->Dump();
510 t_level2->Fill();
511 if( l2_event->ntrk()>0 )ntrk++;
512 };
513 };
514 }; // end loop on events
515 if(p->VerboseMode()){
516 cout << " >>> processed "<< ev_count <<" events"<< endl;
517 if(ev_count)cout << ntrk << " events with at least one track ("<<(Int_t)(100*ntrk/ev_count)<<"%)\n";
518 reduction->Show("reduction");
519 }
520 delete reduction;
521
522 delete dbtime;
523
524 }; // end loop on runs
525
526
527 // ------------------------------------------------------------
528 // if reprocessing one run, copy all the events AFTER the run
529 // ------------------------------------------------------------
530 if( !(p->standalone) ){
531 for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
532 t_clone->GetEntry(i);
533 *l2_event = *l2_clone;
534 t_level2->Fill();
535 l2_event->Clear();
536 // COPY COPY COPY
537 };
538 };
539 // ---------------
540 // close the files
541 // ---------------
542 if(p->get2){
543 if( t_clone )t_clone->Delete("all");//delete old tree from file
544 if( !(p->standalone) )runinfo->Close();
545 f2->Write("Tracker");
546 if( t_level2 )t_level2->Delete(); //delete new tree from memory
547
548 };
549
550 if(f0) if(f0->IsOpen()) f0->Close();
551 // if( f0_c->IsOpen() )f0_c->Close();
552
553 // lh_event->Delete();
554 l1_event->Delete();
555 l2_event->Delete();
556 l2_clone->Delete();
557
558 return(p->ostatus);
559 }
560

  ViewVC Help
Powered by ViewVC 1.1.23