/[PAMELA software]/DarthVader/S4Level2/src/S4Core.cpp
ViewVC logotype

Contents of /DarthVader/S4Level2/src/S4Core.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.16 - (show annotations) (download)
Wed May 30 11:45:11 2007 UTC (17 years, 7 months ago) by mocchiut
Branch: MAIN
Changes since 1.15: +2 -1 lines
S4Core crashes if no calibration are found, fixed

1 //
2 // Given a calibration and a data file this program create an ntuple with LEVEL2 S4 variables
3 //
4 //
5 // C/C++ headers
6 //
7 #include <fstream>
8 #include <string.h>
9 #include <iostream>
10 #include <cstring>
11 #include <stdio.h>
12 //
13 // ROOT headers
14 //
15 #include <TGraph.h>
16 #include <TF1.h>
17 #include <TTree.h>
18 #include <TClassEdit.h>
19 #include <TObject.h>
20 #include <TList.h>
21 #include <TArrayI.h>
22 #include <TArrayD.h>
23 #include <TSystem.h>
24 #include <TSystemDirectory.h>
25 #include <TString.h>
26 #include <TFile.h>
27 #include <TClass.h>
28 #include <TSQLServer.h>
29 #include <TSQLRow.h>
30 #include <TSQLResult.h>
31 #include <TClonesArray.h>
32 #include <stdlib.h>
33 #include <math.h>
34 //
35 // RunInfo header
36 //
37 #include <RunInfo.h>
38 #include <GLTables.h>
39 //
40 // YODA headers
41 //
42 #include <PamelaRun.h>
43 #include <PscuHeader.h>
44 #include <PscuEvent.h>
45 #include <EventHeader.h>
46 #include <CalibS4Event.h>
47 #include <physics/S4/S4Event.h>
48 //
49 // This program headers
50 //
51 #include <S4Level2.h>
52 #include <S4Core.h>
53 #include <S4Verl2.h>
54 //
55 using namespace std;
56 //
57 /*
58 * Fit function Received from Valeria Malvezzi 06/02/2006
59 */
60 Double_t fitf(Double_t *x, Double_t *par){
61 Double_t fitval =(par[0]*x[0])+par[1];
62 return fitval;
63 }
64
65 /*
66 * Fit the S4 calibration with a straight line - Received from Valeria Malvezzi 06/02/2006
67 */
68 TArrayD *S4_paramfit(UInt_t atime, TSQLServer *dbc){
69 //
70 TArrayD *parametri = new TArrayD(2);
71 //
72 GL_S4_CALIB *glS4calib = new GL_S4_CALIB();
73 //
74 if ( !dbc->IsConnected() ) throw -504;
75 Int_t s4sig = glS4calib->Query_GL_S4_CALIB(atime, dbc);
76 if ( s4sig < 0 ) throw s4sig;
77 //
78 GL_ROOT *glroot = new GL_ROOT();
79 if ( !dbc->IsConnected() ) throw -504;
80 glroot->Query_GL_ROOT(glS4calib->ID_ROOT_L0,dbc);
81 //
82 stringstream ftmpname;
83 ftmpname.str("");
84 ftmpname << glroot->PATH.Data() << "/";
85 ftmpname << glroot->NAME.Data();
86 //
87 TFile *file = 0;
88 file = new TFile(ftmpname.str().c_str());
89 //
90 if ( !file ) return(NULL);
91 //
92 TTree *tr = 0;
93 tr = (TTree*)file->Get("CalibS4");
94 if ( !tr ){
95 file->Close();
96 return(NULL);
97 };
98 //
99 pamela::CalibS4Event *S4CalibEvent = 0;
100 tr->SetBranchAddress("CalibS4", &S4CalibEvent);
101 if ( tr->GetEntries() < glS4calib->EV_ROOT ) return(NULL);
102 //
103 tr->GetEntry(glS4calib->EV_ROOT);
104 //
105 // Variables initialization
106 //
107 Double_t mip[3]={1, 30, 300};
108 Double_t adc[3] = {0.,0.,0.};
109 //
110 // Fit calibrations and find parameters to calibrate data
111 //
112 pamela::S4::S4Event *s4Record = 0;
113 //
114 for (Int_t j = 0; j < 4; j++){
115 for (Int_t i = 0; i < 128; i++){
116 s4Record = (pamela::S4::S4Event*)S4CalibEvent->Records->At((j*128 + i));
117 switch (j) {
118 case 0 :{
119 adc[0]=adc[0]+((s4Record->S4_DATA)-32);
120 break;
121 };
122 case 1 :{
123 adc[1]=adc[1]+((s4Record->S4_DATA)-32);
124 break;
125 };
126 case 3 :{
127 adc[2]=adc[2]+((s4Record->S4_DATA)-32);
128 break;
129 };
130 };
131 };
132 };
133 //
134 adc[0]=adc[0]/128;
135 adc[1]=adc[1]/128;
136 adc[2]=adc[2]/128;
137 // if ( IsDebug() ) printf(" adc1 = %g adc2 = %g adc3 = %g\n",adc[0],adc[1],adc[2]);
138 TGraph *fitpar = new TGraph (3, adc, mip);
139 TF1 *func = new TF1("fitf", fitf, -1., 400., 2); // function definition with 2 parameters
140 //
141 func->SetParameters(0.3,1.); // parameters initialization to 1
142 func->SetParNames("m","q"); //parameter's name
143 fitpar->Fit(func,"qr"); //function fit
144 parametri->AddAt(func->GetParameter(0),0);
145 parametri->AddAt(func->GetParameter(1),1);
146 // if ( parametri[0] < 0. || parametri[0] > 0.5 || parametri[1] < 0.7 || parametri[1] > 1.1 ) valid = 0;
147 //
148 delete glroot;
149 delete glS4calib;
150 delete fitpar;
151 delete func;
152 file->Close();
153 //
154 return parametri;
155 };
156
157 //
158 // CORE ROUTINE
159 //
160 //
161 int S4Core(UInt_t run, TFile *file, TSQLServer *dbc, Int_t S4argc, char *S4argv[]){
162 //
163 // Set these to true to have a verbose output.
164 //
165 Bool_t debug = false;
166 Bool_t verbose = false;
167 //
168 //
169 // Output directory is the working directoy.
170 //
171 const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));
172 //
173 Int_t i = 0;
174 TString processFolder = Form("S4Folder_%u",run);
175 if ( S4argc > 0 ){
176 i = 0;
177 while ( i < S4argc ){
178 if ( !strcmp(S4argv[i],"-processFolder") ) {
179 if ( S4argc < i+1 ){
180 throw -3;
181 };
182 processFolder = (TString)S4argv[i+1];
183 i++;
184 };
185 if ( (!strcmp(S4argv[i],"--debug")) || (!strcmp(S4argv[i],"-g"))) {
186 verbose = true;
187 debug = true;
188 };
189 if ( (!strcmp(S4argv[i],"--verbose")) || (!strcmp(S4argv[i],"-v"))) {
190 verbose = true;
191 };
192 i++;
193 };
194 };
195 // Variables for level2
196 //
197 TTree *S4tr = 0;
198 UInt_t nevents = 0;
199 //
200 // Variables needed to reprocess data
201 //
202 Long64_t maxsize = 10000000000LL;
203 TTree::SetMaxTreeSize(maxsize);
204 //
205 TString S4version;
206 ItoRunInfo *runinfo = 0;
207 TArrayI *runlist = 0;
208 TTree *S4trclone = 0;
209 Bool_t reproc = false;
210 Bool_t reprocall = false;
211 UInt_t nobefrun = 0;
212 UInt_t noaftrun = 0;
213 UInt_t numbofrun = 0;
214 stringstream ftmpname;
215 TString fname;
216 UInt_t totfileentries = 0;
217 UInt_t idRun = 0;
218 Double_t ParamFit0 = 0.;
219 Double_t ParamFit1 = 0.;
220 //
221 // variables needed to handle error signals
222 //
223 Int_t code = 0;
224 Int_t sgnl;
225 //
226 // S4 level2 classes
227 //
228 S4Level2 *s4 = new S4Level2();
229 S4Level2 *s4clone = new S4Level2();
230 //
231 // define variables for opening and reading level0 file
232 //
233 TFile *l0File = 0;
234 TTree *l0tr = 0;
235 TBranch *l0head = 0;
236 TBranch *l0S4 =0;
237 pamela::S4::S4Event *l0s4e = 0;
238 pamela::EventHeader *eh = 0;
239 pamela::PscuHeader *ph = 0;
240 //
241 // Define other basic variables
242 //
243 UInt_t procev = 0;
244 stringstream file2;
245 stringstream file3;
246 stringstream qy;
247 Int_t totevent = 0;
248 UInt_t atime = 0;
249 // Int_t ei = 0;
250 UInt_t re = 0;
251 //
252 // Working filename
253 //
254 TString outputfile;
255 stringstream name;
256 name.str("");
257 name << outDir << "/";
258 //
259 // temporary file and folder
260 //
261 TFile *tempfile = 0;
262 TTree *tempS4 = 0;
263 stringstream tempname;
264 stringstream S4folder;
265 tempname.str("");
266 tempname << outDir;
267 tempname << "/" << processFolder.Data();
268 S4folder.str("");
269 S4folder << tempname.str().c_str();
270 gSystem->MakeDirectory(S4folder.str().c_str());
271 tempname << "/S4tree_run";
272 tempname << run << ".root";
273 //
274 // DB classes
275 //
276 GL_ROOT *glroot = new GL_ROOT();
277 GL_TIMESYNC *dbtime = 0;
278 //
279 // Let's start!
280 //
281 // 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
282 // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file
283 // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
284 //
285 if ( run == 0 ) reproc = true;
286 //
287 //
288 // Output file is "outputfile"
289 //
290 if ( !file->IsOpen() ){
291 //printf(" S4 - ERROR: cannot open file for writing\n");
292 throw -501;
293 };
294 //
295 // Retrieve GL_RUN variables from the level2 file
296 //
297 S4version = S4Info(false); // we should decide how to handle versioning system
298 //
299 // create an interface to RunInfo called "runinfo"
300 //
301 runinfo = new ItoRunInfo(file);
302 //
303 // open "Run" tree in level2 file, if not existing return an error (sngl != 0)
304 //
305 sgnl = 0;
306 sgnl = runinfo->Update(run,"S4",S4version);
307 if ( sgnl ){
308 //printf(" S4 - ERROR: RunInfo exited with non-zero status\n");
309 code = sgnl;
310 goto closeandexit;
311 } else {
312 sgnl = 0;
313 };
314 //
315 // number of events in the file BEFORE the first event of our run
316 //
317 nobefrun = runinfo->GetFirstEntry();
318 //
319 // total number of events in the file
320 //
321 totfileentries = runinfo->GetFileEntries();
322 //
323 // first file entry AFTER the last event of our run
324 //
325 noaftrun = runinfo->GetLastEntry() + 1;
326 //
327 // number of run to be processed
328 //
329 numbofrun = runinfo->GetNoRun();
330 UInt_t totnorun = runinfo->GetRunEntries();
331 //
332 // Try to access the S4 tree in the file, if it exists we are reprocessing data if not we are processing a new run
333 //
334 S4trclone = (TTree*)file->Get("S4");
335 //
336 if ( !S4trclone ){
337 //
338 // tree does not exist, we are not reprocessing
339 //
340 reproc = false;
341 if ( run == 0 ){
342 if (verbose) printf(" S4 - WARNING: you are reprocessing data but S4 tree does not exist!\n");
343 }
344 if ( runinfo->IsReprocessing() && run != 0 ) {
345 if (verbose) printf(" S4 - WARNING: it seems you are not reprocessing data but S4\n versioning information already exists in RunInfo.\n");
346 }
347 } else {
348 //
349 // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?
350 //
351 S4trclone->SetAutoSave(900000000000000LL);
352 reproc = true;
353 //
354 // update versioning information
355 //
356 if (verbose) printf("\n Preparing the pre-processing...\n");
357 //
358 if ( run == 0 || totnorun == 1 ){
359 //
360 // we are reprocessing all the file
361 // 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
362 //
363 reprocall = true;
364 //
365 if (verbose) printf("\n S4 - WARNING: Reprocessing all runs\n");
366 //
367 } else {
368 //
369 // 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
370 //
371 reprocall = false;
372 //
373 if (verbose) printf("\n S4 - WARNING: Reprocessing run number %u \n",run);
374 //
375 // copying old tree to a new file
376 //
377 tempfile = new TFile(tempname.str().c_str(),"RECREATE");
378 tempS4 = S4trclone->CloneTree(-1,"fast");
379 tempS4->SetName("S4-old");
380 tempfile->Write();
381 tempfile->Close();
382 }
383 //
384 // Delete the old tree from old file and memory
385 //
386 S4trclone->Delete("all");
387 //
388 if (verbose) printf(" ...done!\n");
389 //
390 };
391 //
392 // create mydetector tree mydect
393 //
394 file->cd();
395 S4tr = new TTree("S4-new","PAMELA Level2 S4 data");
396 S4tr->SetAutoSave(900000000000000LL);
397 S4tr->Branch("S4Level2","S4Level2",&s4);
398 //
399 if ( reproc && !reprocall ){
400 //
401 // open new file and retrieve also tree informations
402 //
403 tempfile = new TFile(tempname.str().c_str(),"READ");
404 S4trclone = (TTree*)tempfile->Get("S4-old");
405 S4trclone->SetAutoSave(900000000000000LL);
406 S4trclone->SetBranchAddress("S4Level2",&s4clone);
407 //
408 if ( nobefrun > 0 ){
409 if (verbose){
410 printf("\n Pre-processing: copying events from the old tree before the processed run\n");
411 printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
412 printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
413 }
414 for (UInt_t j = 0; j < nobefrun; j++){
415 //
416 S4trclone->GetEntry(j);
417 //
418 // copy s4clone to mydec
419 //
420 s4->Clear();
421 //
422 memcpy(&s4,&s4clone,sizeof(s4clone));
423 //
424 // Fill entry in the new tree
425 //
426 S4tr->Fill();
427 //
428 };
429 if (verbose) printf(" Finished successful copying!\n");
430 };
431 };
432 //
433 // Get the list of run to be processed, if only one run has to be processed the list will contain one entry only.
434 //
435 runlist = runinfo->GetRunList();
436 //
437 // Loop over the run to be processed
438 //
439 for (UInt_t irun=0; irun < numbofrun; irun++){
440 //
441 // retrieve the first run ID to be processed using the RunInfo list
442 //
443 idRun = runlist->At(irun);
444 if (verbose){
445 printf("\n\n\n ####################################################################### \n");
446 printf(" PROCESSING RUN NUMBER %u \n",idRun);
447 printf(" ####################################################################### \n\n\n");
448 }
449 //
450 runinfo->ID_ROOT_L0 = 0;
451 //
452 // store in the runinfo class the GL_RUN variables for our run
453 //
454 sgnl = 0;
455 sgnl = runinfo->GetRunInfo(idRun);
456 if ( sgnl ){
457 if ( debug ) printf(" S4 - ERROR: RunInfo exited with non-zero status\n");
458 code = sgnl;
459 goto closeandexit;
460 } else {
461 sgnl = 0;
462 };
463 //
464 // now you can access that variables using the RunInfo class this way runinfo->ID_ROOT_L0
465 //
466 if ( runinfo->ID_ROOT_L0 == 0 ){
467 if ( debug ) printf("\n S4 - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
468 code = -5;
469 goto closeandexit;
470 };
471 //
472 // prepare the timesync for the db
473 //
474 if ( !dbc->IsConnected() ) throw -504;
475 dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
476 //
477 // Search in the DB the path and name of the LEVEL0 file to be processed.
478 //
479 if ( !dbc->IsConnected() ) throw -504;
480 glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
481 //
482 ftmpname.str("");
483 ftmpname << glroot->PATH.Data() << "/";
484 ftmpname << glroot->NAME.Data();
485 fname = ftmpname.str().c_str();
486 //
487 // print out informations
488 //
489 totevent = runinfo->NEVENTS;
490 if (verbose){
491 printf("\n LEVEL0 data file: %s \n",fname.Data());
492 printf(" RUN HEADER absolute time is: %u \n",runinfo->RUNHEADER_TIME);
493 printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
494 printf(" %i events to be processed for run %u: from %u to %u \n\n",totevent,idRun,runinfo->EV_FROM,(runinfo->EV_FROM+totevent));
495 }
496 //
497 // Open Level0 file
498 //
499 l0File = new TFile(fname.Data());
500 if ( !l0File ) {
501 if ( debug ) printf(" S4 - ERROR: problems opening Level0 file\n");
502 code = -6;
503 goto closeandexit;
504 };
505 //
506 l0tr = (TTree*)l0File->Get("Physics");
507 if ( !l0tr ) {
508 if ( debug ) printf(" S4 - ERROR: no Physics tree in Level0 file\n");
509 l0File->Close();
510 code = -7;
511 goto closeandexit;
512 };
513 l0head = l0tr->GetBranch("Header");
514 if ( !l0head ) {
515 if ( debug ) printf(" S4 - ERROR: no Header branch in Level0 tree\n");
516 l0File->Close();
517 code = -8;
518 goto closeandexit;
519 };
520 l0S4 = l0tr->GetBranch("S4");
521 if ( !l0S4 ) {
522 if ( debug ) printf(" S4 - ERROR: no S4 branch in Level0 tree\n");
523 l0File->Close();
524 code = -503;
525 goto closeandexit;
526 };
527 //
528 l0tr->SetBranchAddress("S4", &l0s4e);
529 l0tr->SetBranchAddress("Header", &eh);
530 //
531 nevents = l0S4->GetEntries();
532 //
533 if ( nevents < 1 ) {
534 if ( debug ) printf(" S4 - ERROR: Level0 file is empty\n\n");
535 l0File->Close();
536 code = -11;
537 goto closeandexit;
538 };
539 //
540 if ( runinfo->EV_TO > nevents-1 ) {
541 if ( debug ) printf(" S4 - ERROR: too few entries in the S4 tree\n");
542 l0File->Close();
543 code = -12;
544 goto closeandexit;
545 };
546 //
547 // Check if we have to load parameter files (or calibration associated to runs and not to events)
548 // second query: which is the value of paramfit relative to the calibration?
549 //
550 TArrayD *params = 0;
551 params = S4_paramfit(runinfo->RUNHEADER_TIME, dbc);
552 if ( !params ){
553 code = -13;
554 goto closeandexit;
555 };
556 //
557 ParamFit0 = params->At(0);
558 ParamFit1 = params->At(1);
559 //
560 // run over all the events of the run
561 //
562 if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
563 //
564 for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
565 //
566 if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);
567 //
568 l0head->GetEntry(re);
569 //
570 // absolute time of this event
571 //
572 ph = eh->GetPscuHeader();
573 atime = dbtime->DBabsTime(ph->GetOrbitalTime());
574 //
575 // paranoid check
576 //
577 if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME) ) {
578 if (verbose) printf(" S4 - WARNING: event at time outside the run time window, skipping it\n");
579 goto jumpev;
580 };
581 //
582 procev++;
583 //
584 // start processing
585 //
586 s4->Clear();
587 l0S4->GetEntry(re);
588 if (l0s4e->unpackError == 0){
589 s4->S4adc = l0s4e->S4_DATA;
590 //
591 if ((l0s4e->S4_DATA) > 31 ){
592 s4->S4calibrated = ParamFit0*((l0s4e->S4_DATA)-32)+ParamFit1;
593 }else{
594 s4->S4calibrated = 0;
595 }
596 };
597 //
598 s4->unpackError = l0s4e->unpackError;
599 //
600 S4tr->Fill();
601 //
602 //
603 jumpev:
604 debug = false;
605 //
606 };
607 //
608 // Here you may want to clear some variables before processing another run
609 //
610 delete dbtime;
611 if ( params ) delete params;
612 //
613 }; // process all the runs
614 //
615 if (verbose) printf("\n Finished processing data \n");
616 //
617 closeandexit:
618 //
619 // 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.
620 //
621 if ( !reprocall && reproc && code >= 0 ){
622 if ( totfileentries > noaftrun ){
623 if (verbose){
624 printf("\n Post-processing: copying events from the old tree after the processed run\n");
625 printf(" Copying %i events in the file which are after the end of the run %i \n",(int)(totfileentries-noaftrun),(int)run);
626 printf(" Start copying at event number %i end copying at event number %i \n",(int)noaftrun,(int)totfileentries);
627 }
628 for (UInt_t j = noaftrun; j < totfileentries; j++ ){
629 //
630 // Get entry from old tree
631 //
632 S4trclone->GetEntry(j);
633 //
634 // copy s4clone to s4
635 //
636 // s4 = new S4Level2();
637 s4->Clear();
638 memcpy(&s4,&s4clone,sizeof(s4clone));
639 //
640 // Fill entry in the new tree
641 //
642 S4tr->Fill();
643 };
644 if (verbose) printf(" Finished successful copying!\n");
645 };
646 };
647 //
648 // Close files, delete old tree(s), write and close level2 file
649 //
650 if ( l0File ) l0File->Close();
651 if ( tempfile ) tempfile->Close();
652 gSystem->Unlink(tempname.str().c_str());
653 //
654 if ( runinfo ) runinfo->Close();
655 if ( S4tr ) S4tr->SetName("S4");
656 if ( file ){
657 file->cd();
658 file->Write();
659 };
660 //
661 gSystem->Unlink(S4folder.str().c_str());
662 //
663 // the end
664 //
665 if (verbose) printf("\n Exiting...\n");
666 if ( S4tr ) S4tr->Delete();
667 //
668 if ( s4 ) delete s4;
669 if ( s4clone ) delete s4clone;
670 if ( glroot ) delete glroot;
671 if ( runinfo ) delete runinfo;
672 //
673 if(code < 0) throw code;
674 return(code);
675 }
676
677

  ViewVC Help
Powered by ViewVC 1.1.23