/[PAMELA software]/DarthVader/NDLevel2/src/NDCore.cpp
ViewVC logotype

Contents of /DarthVader/NDLevel2/src/NDCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.16 - (show annotations) (download)
Thu Sep 6 12:49:28 2007 UTC (17 years, 5 months ago) by mocchiut
Branch: MAIN
Changes since 1.15: +5 -3 lines
Small improvements of *Core routines, added simulation identification feature in CaloCore

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

  ViewVC Help
Powered by ViewVC 1.1.23