/[PAMELA software]/calo/flight/MyDetector2Level2/src/MyDect2Core.cpp
ViewVC logotype

Contents of /calo/flight/MyDetector2Level2/src/MyDect2Core.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (show annotations) (download)
Wed May 3 10:03:40 2006 UTC (18 years, 8 months ago) by mocchiut
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +101 -195 lines
Major update, requires GLTables and RunInfo

1 //
2 // C/C++ headers
3 //
4 #include <fstream>
5 #include <string.h>
6 //
7 // ROOT headers
8 //
9 #include <TTree.h>
10 #include <TClassEdit.h>
11 #include <TObject.h>
12 #include <TList.h>
13 #include <TArrayL.h>
14 #include <TSystem.h>
15 #include <TSystemDirectory.h>
16 #include <TString.h>
17 #include <TFile.h>
18 #include <TClass.h>
19 #include <TCanvas.h>
20 #include <TH1.h>
21 #include <TH1F.h>
22 #include <TH2D.h>
23 #include <TLatex.h>
24 #include <TPad.h>
25 #include <TSQLServer.h>
26 #include <TSQLRow.h>
27 #include <TSQLResult.h>
28 #include <TClonesArray.h>
29 //
30 // YODA headers
31 //
32 #include <PamelaRun.h>
33 #include <RegistryEvent.h>
34 #include <physics/trigger/TriggerEvent.h>
35 //
36 // RunInfo header
37 //
38 #include <RunInfo.h>
39 //
40 // This program headers
41 //
42 #include <MyDect2Core.h>
43 #include <MyDect2Level2.h>
44 #include <MyDect2Verl2.h>
45 //
46 //
47 // Declaration of the core fortran routines
48 //
49 #define mydeccom mydeccom_
50 extern "C" int mydeccom();
51 #define mydectrk mydectrk_
52 extern "C" int mydectrk();
53 //
54 // Tracker classes headers and definitions
55 //
56 #include <TrkLevel2.h>
57 //
58 using namespace std;
59 //
60 Bool_t existfile(TString filename){
61 ifstream myfile;
62 myfile.open(filename.Data());
63 if ( !myfile ){
64 return(false);
65 };
66 myfile.close();
67 return(true);
68 }
69 //
70 //
71 // CORE ROUTINE
72 //
73 //
74 int MyDect2Core(ULong64_t run, TString filename, TString processFolder, TSQLServer *dbc){
75 //
76 // Set these to true to have a very verbose output.
77 //
78 // Bool_t debug = true;
79 Bool_t debug = false;
80 //
81 // Output directory is the working directoy.
82 //
83 const char* outdir = gSystem->WorkingDirectory();
84 //
85 // Variables for level2
86 //
87 TFile *file = 0;
88 TTree *tracker = 0;
89 TTree *mydect = 0;
90 Long64_t nevents = 0LL;
91 //
92 // variables needed to reprocess data
93 //
94 TString mydectversion;
95 ItoRunInfo *runinfo = 0;
96 TArrayL *runlist = 0;
97 TTree *mydectclone = 0;
98 Bool_t reproc = false;
99 Bool_t reprocall = false;
100 UInt_t nobefrun = 0;
101 UInt_t noaftrun = 0;
102 UInt_t numbofrun = 0;
103 stringstream ftmpname;
104 TString fname;
105 Long64_t totfileentries = 0ULL;
106 Long64_t idRun = 0LL;
107 //
108 // variables needed to handle error signals
109 //
110 Int_t code = 0;
111 Int_t sgnl;
112 //
113 // mydetector2 level2 classes
114 //
115 MyDect2Level2 *mydec = new MyDect2Level2();
116 MyDect2Level2 *mydecclone = new MyDect2Level2();
117 //
118 // tracker level2 variables
119 //
120 TrkLevel2 *trk = new TrkLevel2();
121 Long64_t nevtrkl2 = 0;
122 //
123 // define variables for opening and reading level0 file
124 //
125 TFile *l0File = 0;
126 TTree *l0tr = 0;
127 TBranch *l0head = 0;
128 TBranch *l0registry = 0;
129 TBranch *l0trig = 0;
130 pamela::RegistryEvent *l0reg=0;
131 pamela::trigger::TriggerEvent *trig = 0;
132 //
133 // Define other basic variables
134 //
135 UInt_t procev = 0;
136 stringstream file2;
137 stringstream file3;
138 stringstream qy;
139 Int_t itr = -1;
140 Int_t totevent = 0;
141 ULong64_t atime = 0ULL;
142 Int_t ei = 0;
143 Int_t re = 0;
144 //
145 // Working filename
146 //
147 TString outputfile;
148 stringstream name;
149 name.str("");
150 name << outdir << "/";
151 //
152 // temporary file and folder
153 //
154 TFile *tempfile = 0;
155 TTree *tempmydect = 0;
156 stringstream tempname;
157 stringstream mydectfolder;
158 tempname.str("");
159 tempname << outdir;
160 tempname << "/" << processFolder.Data();
161 mydectfolder.str("");
162 mydectfolder << tempname.str().c_str();
163 gSystem->MakeDirectory(mydectfolder.str().c_str());
164 tempname << "/mydect2tree_run";
165 tempname << run << ".root";
166 //
167 // variables needed to load magnetic field maps
168 //
169 Int_t ntrkentry = 0;
170 ULong64_t tttrkpar1 = 0ULL;
171 Bool_t trkpar1 = true;
172 //
173 // DB classes
174 //
175 GL_ROOT *glroot = new GL_ROOT();
176 GL_PARAM *glparam = new GL_PARAM();
177 //
178 // declaring external output and input structures
179 //
180 extern struct InputStruct cinput_;
181 extern struct OutputStruct coutput_;
182 //
183 // Let's start!
184 //
185 //
186 // As a first thing we must check what we have to do: if run = 0 we must process all events in the file has been passed
187 // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file
188 // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
189 //
190 if ( run != 0ULL ){
191 //
192 // Run is not 0, we must process only one run .
193 //
194 //
195 // If no filename is given, use the deafult one run_id.Level2.root.
196 //
197 if ( !strcmp(filename.Data(),"") ){
198 name << run << ".Level2.root";
199 } else {
200 name << filename.Data();
201 };
202 //
203 } else {
204 //
205 // Run is 0, we must process all entries in the given file.
206 //
207 // OUCH! no filename has been given, run out!
208 //
209 if ( !strcmp(filename.Data(),"") ) {
210 printf(" MYDETECTOR2 - ERROR: processing all runs but no filename is given\n");
211 return(-4);
212 };
213 name << filename.Data();
214 reproc = true;
215 //
216 };
217 //
218 // Output file is "outputfile"
219 //
220 outputfile = name.str().c_str();
221 printf("\n Output filename is: \n %s \n\n",outputfile.Data());
222 //
223 // Check if file exists, if not exit with error (we need tracker data).
224 //
225 if ( !existfile(outputfile.Data()) ) {
226 printf(" MYDETECTOR2 - ERROR: no Level2 file\n");
227 return(-100);
228 };
229 //
230 // OK, file exists, open it in "update" mode.
231 //
232 file = new TFile(outputfile.Data(),"UPDATE");
233 if ( !file->IsOpen() ){
234 printf(" MYDETECTOR2 - ERROR: cannot open file for writing\n");
235 return(-101);
236 };
237 //
238 // Does it contain the Tracker tree?
239 //
240 tracker = (TTree*)file->Get("Tracker");
241 if ( !tracker ) {
242 printf(" MYDETECTOR2 - ERROR: no tracker tree\n");
243 code = -102;
244 goto closeandexit;
245 };
246 //
247 // get tracker level2 data pointer
248 //
249 tracker->SetBranchAddress("TrkLevel2",&trk);
250 nevtrkl2 = tracker->GetEntries();
251 //
252 // Retrieve GL_RUN variables from the level2 file
253 //
254 // mydectversion = MyDect2Info(false); // we should decide how to handle versioning system
255 //
256 // create an interface to RunInfo called "runinfo"
257 //
258 runinfo = new ItoRunInfo(file);
259 //
260 // open "Run" tree in level2 file, if not existing return an error (sngl != 0)
261 //
262 sgnl = 0;
263 sgnl = runinfo->Read(run);
264 if ( sgnl ){
265 printf(" MYDETECTOR2 - ERROR: RunInfo exited with non-zero status\n");
266 code = sgnl;
267 goto closeandexit;
268 } else {
269 sgnl = 0;
270 };
271 //
272 // number of events in the file BEFORE the first event of our run
273 //
274 nobefrun = runinfo->GetFirstFileEntry();
275 //
276 // total number of events in the file
277 //
278 totfileentries = runinfo->GetFileEntries();
279 //
280 // first file entry AFTER the last event of our run
281 //
282 noaftrun = runinfo->GetLastFileEntry();
283 //
284 // number of run to be processed
285 //
286 numbofrun = runinfo->GetNoRun();
287 //
288 // Try to access the MyDetector2 tree in the file, if it exists we are reprocessing data if not we are processing a new run
289 //
290 mydectclone = (TTree*)file->Get("MyDetector2");
291 //
292 if ( !mydectclone ){
293 //
294 // tree does not exist, we are not reprocessing
295 //
296 reproc = false;
297 if ( run == 0ULL ) printf(" MYDETECTOR2 - WARNING: you are reprocessing data but MyDetector2 tree does not exist!\n");
298 if ( runinfo->IsReprocessing() && run != 0ULL ) printf(" MYDETECTOR2 - WARNING: it seems you are not reprocessing data but mydetector2\n versioning information already exists in RunInfo.\n");
299
300 } else {
301 //
302 // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?
303 //
304 reproc = true;
305 //
306 // update versioning information
307 //
308 // sgnl = runinfo->Update(run,"MYDECT2");
309 //
310 printf("\n Preparing the pre-processing...\n");
311 //
312 if ( run == 0ULL ){
313 //
314 // we are reprocessing all the file
315 // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately
316 //
317 reprocall = true;
318 //
319 printf("\n MYDETECTOR2 - WARNING: Reprocessing all runs\n");
320 //
321 } else {
322 //
323 // we are reprocessing a single run, we must copy to the new tree the events in the file which preceed the first event of the run
324 //
325 reprocall = false;
326 //
327 printf("\n MYDETECTOR2 - WARNING: Reprocessing run number %llu \n",run);
328 //
329 // copying old tree to a new file
330 //
331 tempfile = new TFile(tempname.str().c_str(),"RECREATE");
332 tempmydect = mydectclone->CloneTree(-1,"fast");
333 tempmydect->SetName("MyDetector2-old");
334 tempfile->Write();
335 tempfile->Close();
336 }
337 //
338 // Delete the old tree from old file and memory
339 //
340 mydectclone->Delete("all");
341 //
342 printf(" ...done!\n");
343 //
344 };
345 //
346 // create mydetector tree mydect
347 //
348 file->cd();
349 mydect = new TTree("MyDetector2-new","PAMELA Level2 MyDetector2 data");
350 mydect->Branch("MyDect2Level2","MyDect2Level2",&mydec);
351 //
352 if ( reproc && !reprocall ){
353 //
354 // open new file and retrieve alo tree informations
355 //
356 tempfile = new TFile(tempname.str().c_str(),"READ");
357 mydectclone = (TTree*)tempfile->Get("MyDetector2-old");
358 mydectclone->SetBranchAddress("MyDect2Level2",&mydecclone);
359 //
360 if ( nobefrun > 0 ){
361 printf("\n Pre-processing: copying events from the old tree before the processed run\n");
362 printf(" Copying %u events in the file which are before the beginning of the run %llu \n",nobefrun,run);
363 printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
364 for (UInt_t j = 0; j < nobefrun; j++){
365 //
366 mydectclone->GetEntry(j);
367 //
368 // copy mydecclone to mydec
369 //
370 mydec = new MyDect2Level2();
371 memcpy(&mydec,&mydecclone,sizeof(mydecclone));
372 //
373 // Fill entry in the new tree
374 //
375 mydect->Fill();
376 //
377 };
378 printf(" Finished successful copying!\n");
379 };
380 };
381 //
382 // Get the list of run to be processed, if only one run has to be processed the list will contain one entry only.
383 //
384 runlist = runinfo->GetRunList();
385 //
386 // Loop over the run to be processed
387 //
388 for (UInt_t irun=0; irun < numbofrun; irun++){
389 //
390 // retrieve the first run ID to be processed using the RunInfo list
391 //
392 idRun = runlist->At(irun);
393 printf("\n\n\n ####################################################################### \n");
394 printf(" PROCESSING RUN NUMBER %i \n",(int)idRun);
395 printf(" ####################################################################### \n\n\n");
396 //
397 runinfo->ID_REG_RUN = 0ULL;
398 //
399 // store in the runinfo class the GL_RUN variables for our run
400 //
401 sgnl = 0;
402 sgnl = runinfo->GetRunInfo(idRun);
403 if ( sgnl ){
404 printf(" MYDETECTOR2 - ERROR: RunInfo exited with non-zero status\n");
405 code = sgnl;
406 goto closeandexit;
407 } else {
408 sgnl = 0;
409 };
410 //
411 // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN
412 //
413 if ( runinfo->ID_REG_RUN == 0 ){
414 printf("\n MYDETECTOR2 - ERROR: no run with ID_RUN = %i \n\n Exiting... \n\n",(int)idRun);
415 code = -5;
416 goto closeandexit;
417 };
418 //
419 // Search in the DB the path and name of the LEVEL0 file to be processed.
420 //
421 glroot->Query_GL_ROOT(runinfo->ID_REG_RUN,dbc);
422 //
423 ftmpname.str("");
424 ftmpname << glroot->PATH.Data() << "/";
425 ftmpname << glroot->NAME.Data();
426 fname = ftmpname.str().c_str();
427 //
428 // print out informations
429 //
430 totevent = runinfo->EV_REG_PHYS_TO - runinfo->EV_REG_PHYS_FROM + 1;
431 printf("\n LEVEL0 data file: %s \n",fname.Data());
432 printf(" RUN HEADER absolute time is: %llu \n",runinfo->RUNHEADER_TIME);
433 printf(" RUN TRAILER absolute time is: %llu \n",runinfo->RUNTRAILER_TIME);
434 printf(" %i events to be processed for run %llu: from %i to %i (reg entries)\n\n",totevent,idRun,runinfo->EV_REG_PHYS_FROM,runinfo->EV_REG_PHYS_TO);
435 //
436 // Open Level0 file
437 //
438 l0File = new TFile(fname.Data());
439 if ( !l0File ) {
440 printf(" MYDETECTOR2 - ERROR: problems opening Level0 file\n");
441 code = -6;
442 goto closeandexit;
443 };
444 l0tr = (TTree*)l0File->Get("Physics");
445 if ( !l0tr ) {
446 printf(" MYDETECTOR2 - ERROR: no Physics tree in Level0 file\n");
447 l0File->Close();
448 code = -7;
449 goto closeandexit;
450 };
451 l0head = l0tr->GetBranch("Header");
452 if ( !l0head ) {
453 printf(" MYDETECTOR2 - ERROR: no Header branch in Level0 tree\n");
454 l0File->Close();
455 code = -8;
456 goto closeandexit;
457 };
458 l0registry = l0tr->GetBranch("Registry");
459 if ( !l0registry ) {
460 printf(" MYDETECTOR2 - ERROR: no Registry branch in Level0 tree\n");
461 l0File->Close();
462 code = -9;
463 goto closeandexit;
464 };
465 l0trig = l0tr->GetBranch("Trigger");
466 if ( !l0trig ) {
467 printf(" MYDETECTOR2 - ERROR: no Trigger branch in Level0 tree\n");
468 l0File->Close();
469 code = -104;
470 goto closeandexit;
471 };
472 //
473 l0tr->SetBranchAddress("Trigger", &trig);
474 l0tr->SetBranchAddress("Registry", &l0reg);
475 //
476 nevents = l0registry->GetEntries();
477 //
478 if ( nevents < 1 ) {
479 printf(" MYDETECTOR2 - ERROR: Level0 file is empty\n\n");
480 l0File->Close();
481 code = -11;
482 goto closeandexit;
483 };
484 //
485 if ( runinfo->EV_REG_PHYS_TO > nevents-1 ) {
486 printf(" MYDETECTOR2 - ERROR: too few entries in the registry tree\n");
487 l0File->Close();
488 code = -12;
489 goto closeandexit;
490 };
491 //
492 // Check if we have to load parameter files (or calibration associated to runs and not to events)
493 //
494 // for example let's assume that we could have different magnetic field maps for different runs:
495 //
496 if ( trkpar1 || ( tttrkpar1 != 0 && tttrkpar1 < runinfo->RUNHEADER_TIME ) ){
497 trkpar1 = false;
498 glparam->Query_GL_PARAM(runinfo->RUNHEADER_TIME,"field",dbc);
499 tttrkpar1 = glparam->TO_TIME;
500 // ----------------------------
501 // Read the magnetic field
502 // ----------------------------
503 printf(" Reading magnetic field maps: \n");
504 trk->LoadField(glparam->PATH+glparam->NAME);
505 printf("\n");
506 };
507 //
508 // run over all the events of the run
509 //
510 printf("\n Ready to start! \n\n Processed events: \n\n");
511 //
512 for ( re = runinfo->EV_REG_PHYS_FROM; re <= runinfo->EV_REG_PHYS_TO; re++){
513 //
514 if ( procev%1000 == 0 && procev > 0 ) printf(" %iK \n",procev/1000);
515 //
516 l0registry->GetEntry(re);
517 //
518 // absolute time of this event
519 //
520 atime = l0reg->absTime;
521 //
522 // physics events is at entry number ei where
523 //
524 ei = l0reg->event;
525 //
526 // paranoid check
527 //
528 if ( atime > runinfo->RUNTRAILER_TIME || atime < runinfo->RUNHEADER_TIME ) {
529 printf(" MYDETECTOR2 - WARNING: event at time outside the run time window, skipping it\n");
530 goto jumpev;
531 };
532 //
533 // retrieve tracker informations, the LEVEL2 entry which correspond to our event will be "itr"
534 //
535 if ( !reprocall ){
536 itr = nobefrun + (re - runinfo->EV_REG_PHYS_FROM);
537 } else {
538 itr = runinfo->GetFirstFileEntry() + (re - runinfo->EV_REG_PHYS_FROM);
539 };
540 if ( itr > nevtrkl2 ){
541 printf(" MYDETECTOR2 - ERROR: no tracker events with entry = %i in Level2 file\n",itr);
542 l0File->Close();
543 code = -113;
544 goto closeandexit;
545 };
546 tracker->GetEntry(itr);
547 //
548 procev++;
549 //
550 // start processing
551 //
552 mydec = new MyDect2Level2();
553
554 //
555 // Here we will use some procedure to calibrate our data and put some kind of informations in the cinput structure
556 //
557 cinput_.myinputvar = (int)ei;
558
559 //
560 // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related variables.
561 //
562 //
563 // Calculate variables common to all tracks
564 //
565 mydeccom();
566 //
567 // and now we must copy from the output structure to the level2 class:
568 //
569 mydec->myvar1 = coutput_.myvar1;
570 mydec->myvar2 = coutput_.myvar2;
571 mydec->myvar3 = coutput_.myvar3;
572 //
573 // Calculate track-related variables
574 //
575 //
576 // Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
577 //
578 ntrkentry = 0;
579 //
580 if ( trk->GetNTracks() > 0 ){
581 //
582 // We have at least one track
583 //
584 //
585 // Run over tracks
586 //
587 for(Int_t nt=0; nt < trk->ntrk(); nt++){
588 //
589 TrkTrack *ptt = trk->GetStoredTrack(nt);
590 //
591 // Copy the alpha vector in the input structure
592 //
593 for (Int_t e = 0; e < 5 ; e++){
594 cinput_.al_pp[e] = ptt->al[e];
595 };
596 //
597 // Get tracker related variables for this track
598 //
599 mydectrk();
600 //
601 // Copy values in the class from the structure (we need to use a temporary class to store variables).
602 //
603 MyDect2TrkVar *t_mydec = new MyDect2TrkVar();
604 t_mydec->mytrkvar = coutput_.mytrkvar;
605 //
606 // Store the tracker track number in order to be sure to have shyncronized data during analysis
607 //
608 t_mydec->trkseqno = nt;
609 //
610 // create a new object for this event with track-related variables
611 //
612 TClonesArray &t = *mydec->myTrk;
613 new(t[ntrkentry]) MyDect2TrkVar(*t_mydec);
614 //
615 ntrkentry++;
616 //
617 }; // loop on all the tracks
618 };
619 //
620 // Here you may want to clear structures used to communicate with fortran
621 //
622 cinput_.myinputvar = 0;
623 //
624 // Fill the rootple
625 //
626 mydect->Fill();
627 //
628 //
629 jumpev:
630 debug = false;
631 //
632 };
633 //
634 // Here you may want to clear some variables before processing another run
635 //
636 ei = 0;
637 }; // process all the runs
638 //
639 printf("\n Finished processing data \n");
640 //
641 closeandexit:
642 //
643 // we have finished processing the run(s). If we processed a single run now we must copy all the events after our run from the old tree to the new one and delete the old tree.
644 //
645 if ( !reprocall && reproc && code >= 0 ){
646 if ( totfileentries > noaftrun ){
647 printf("\n Post-processing: copying events from the old tree after the processed run\n");
648 printf(" Copying %i events in the file which are after the end of the run %i \n",(int)(totfileentries-noaftrun),(int)run);
649 printf(" Start copying at event number %i end copying at event number %i \n",(int)noaftrun,(int)totfileentries);
650 for (UInt_t j = noaftrun; j < totfileentries; j++ ){
651 //
652 // Get entry from old tree
653 //
654 mydectclone->GetEntry(j);
655 //
656 // copy mydecclone to mydec
657 //
658 mydec = new MyDect2Level2();
659 memcpy(&mydec,&mydecclone,sizeof(mydecclone));
660 //
661 // Fill entry in the new tree
662 //
663 mydect->Fill();
664 };
665 printf(" Finished successful copying!\n");
666 };
667 };
668 //
669 // Close files, delete old tree(s), write and close level2 file
670 //
671 if ( l0File ) l0File->Close();
672 if ( tempfile ) tempfile->Close();
673 gSystem->Unlink(tempname.str().c_str());
674 if ( tracker ) tracker->Delete(); // delete tracker tree from memory only to avoid writing a copy to file!
675 //
676 if ( code < 0 ) printf("\n MYDETECTOR2 - ERROR: an error occurred, try to save anyway...\n");
677 printf("\n Writing and closing rootple\n");
678 if ( runinfo ) runinfo->Close();
679 if ( mydect ) mydect->SetName("MyDetector2");
680 if ( file ){
681 file->cd();
682 file->Write();
683 file->Close();
684 };
685 //
686 gSystem->Unlink(mydectfolder.str().c_str());
687 //
688 // the end
689 //
690 printf("\n Exiting...\n");
691 return(code);
692 }
693
694
695

  ViewVC Help
Powered by ViewVC 1.1.23