/[PAMELA software]/anticounter/flight/macros/ACEXE.cra
ViewVC logotype

Contents of /anticounter/flight/macros/ACEXE.cra

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Thu Mar 9 10:13:31 2006 UTC (18 years, 8 months ago) by pam-se
Branch point for: MAIN, v1r00
Initial revision

1 // ------------------------------------------------------------
2 // | |
3 // | AC LEVEL2 Software - P.Hofverberg, E.Mocchiutti, S.Orsi |
4 // | |
5 // | ACEXE.cra version 1.00 (2006-01-19) |
6 // | |
7 // ------------------------------------------------------------
8 //
9 // Based on:
10 // - template TEMPLATEEXE.cra (v2.00 2006-01-13) by E. Mocchiutti
11 // - routine AcLEVEL1.c (v1.01 2005-05-17) by P.Hofverberg
12 //
13 //
14 // Changelog:
15 //
16 //
17 // Comments (TEMPLATEEXE.cra
18 // (2006-01-13): linked to latest external packages (root2paw and calocommon), uses ROOT classes to connect to DB, assume you always want to create
19 // a rootple (no ntuple in intermediate level data) and that you will force processing.
20 //
21 // Include ROOT header files for macro compilation.
22 //
23 #include <TTree.h>
24 #include <TClassEdit.h>
25 #include <TObject.h>
26 #include <TList.h>
27 #include <TSystem.h>
28 #include <TSystemDirectory.h>
29 #include <TString.h>
30 #include <TFile.h>
31 #include <TClass.h>
32 #include <TCanvas.h>
33 #include <TH1.h>
34 #include <TH1F.h>
35 #include <TH2D.h>
36 #include <TLatex.h>
37 #include <TPad.h>
38 //
39 // other general header files
40 //
41 #include <fstream>
42 #include <TSQLServer.h>
43 #include <TSQLRow.h>
44 #include <TSQLResult.h>
45 //
46 // files in the inc directory
47 //
48 #include <acl2class.h>
49 //
50 // these files are the YODA header files needed here (you could need to add more of them)
51 //
52 #include <event/PamelaRun.h>
53 #include <event/RegistryEvent.h>
54 #include <event/physics/anticounter/AnticounterEvent.h>
55 #include <event/CalibCalPedEvent.h>
56 //
57 // this file comes from the calorimeter COMMON package (must be installed)
58 //
59 #include <CaloFunctions.h>
60 //
61 //
62 // MAIN ROUTINE
63 //
64 // Input variables:
65 //
66 // - run : the run ID we want to process.
67 //
68 // - dbc : the pointer to a DB connection to be able to query the DB.
69 //
70 //
71 short int ACEXE(Int_t run, TSQLServer *dbc){
72 //
73 // x Silvio: se non ti piace ACEXE e cambi nome devi poi cambiarlo anche in ACL2.cc
74 //
75 // swcode contain the version of this software. It is an integer of this type: subsubversion*10000+subversion*100+version
76 // For example, version 3.45.99 gives swcode = 994503
77 //
78 Int_t swcode = 2;
79 //
80 // when set debug = 1 you will have a verbose printout on the screen.
81 //
82 int debug = 0;
83 //
84 // filename is the name of the file that contains the run we are going to process; filename will be recovered from the DB.
85 //
86 TString filename;
87 //
88 // where to store results: for the moment we will leave here the starting directory
89 //
90 TString outDir ="./";
91 //
92 // SQL variables
93 //
94 TSQLResult *pResult;
95 TSQLRow *Row;
96 int t;
97 int r;
98 int procev = 0;
99 stringstream myquery;
100 //
101 // Here we query the DB to obtain the informations needed to process the run number "run"
102 //
103 printf("\n Querying DB for data files informations...\n");
104 //
105 // define variables where to store the absolute run header and run trailer times (unsigned long long integers, when set to a number use ULL to store the correct number).
106 //
107 ULong64_t runheadtime = 0ULL;
108 ULong64_t runtrailtime = 0ULL;
109 //
110 // first query: which is the ID of the file which contains the run number "run"?
111 // and at what absolute time the run header and run trailer time were recorded?
112 //
113 myquery.str("");
114 myquery << "select ID_REG_RUNHEADER,RUNHEADER_TIME,RUNTRAILER_TIME from GL_RUN where ID=" << run;
115 myquery << ";";
116 //
117 if ( debug ) printf("The query is:\n \"%s\" \n",myquery.str().c_str());
118 //
119 // query the db and store results in pResult
120 //
121 pResult = dbc->Query(myquery.str().c_str());
122 //
123 Int_t id_reg_runheader=-1;
124 //
125 for( r=0; r < 1000; r++){
126 Row = pResult->Next();
127 if( Row == NULL ) break;
128 for( t = 0; t < pResult->GetFieldCount(); t++){
129 char *pEnd;
130 //
131 // store in id_reg_runheader the answer about the file ID
132 //
133 if ( t == 0 ) id_reg_runheader=atoi(Row->GetField(t));
134 //
135 // store in runheadertime the run header time
136 //
137 if ( t == 1 ) runheadtime = strtoull(Row->GetField(t),&pEnd,0);
138 //
139 // store in runtrailertime the run trailer time
140 //
141 if ( t == 2 ) runtrailtime = strtoull(Row->GetField(t),&pEnd,0);
142 //
143 if ( debug ) printf(" Row %i value = %s \n",t,Row->GetField(t));
144 };
145 };
146 delete pResult;
147 //
148 if ( id_reg_runheader == -1 ){
149 printf("\n ERROR: no run with ID_RUN = %i \n\n Exiting... \n\n",run);
150 return(102);
151 };
152 //
153 // where can I find the file with the id I have just recovered?
154 //
155 myquery.str("");
156 myquery << "select PATH from GL_ROOT where ID=" << id_reg_runheader;
157 myquery << ";";
158 //
159 if ( debug ) printf("The query is:\n \"%s\" \n",myquery.str().c_str());
160 //
161 pResult = dbc->Query(myquery.str().c_str());
162 stringstream fpath;
163 fpath.str("");
164 for( r=0; r < 1000; r++){
165 Row = pResult->Next();
166 if( Row == NULL ) break;
167 for( t = 0; t < pResult->GetFieldCount(); t++){
168 //
169 // put in fpath the path to that file
170 //
171 fpath << Row->GetField(t);
172 if ( debug ) printf(" Row %i value = %s \n",t,Row->GetField(t));
173 };
174 };
175 delete pResult;
176 //
177 // this path is of the form /home/pamela/filesfromyoda/dw_YYMMDD_NNNMM/RunHeader
178 // while we want something of the form of /home/pamela/filesfromyoda/dw_YYMMDD_NNNMM/
179 // here the fpath string is manipulated and the output is put in filename
180 //
181 const string ufpath = fpath.str().c_str();
182 Int_t runheaderpos = ufpath.find("RunHeader");
183 TString runhpath=(fpath.str().c_str());
184 stringcopy(filename,runhpath,0,runheaderpos);
185 //
186 // print out informations
187 //
188 printf("\n LEVEL0 data directory: %s \n",filename.Data());
189 printf(" RUN HEADER absolute time is: %llu \n",runheadtime);
190 printf(" RUN TRAILER absolute time is: %llu \n\n",runtrailtime);
191 //
192 // finished the first set of queries
193 //
194 printf(" ...finished querying DB!\n");
195 //
196 //
197 // define pointers to the root file we want to open
198 //
199 // registry file
200 TFile *regFile;
201 // header file
202 TFile *headerFile;
203 // trigger file
204 TFile *acFile;
205 //
206 TString filety;
207 stringstream file2;
208 stringstream file3;
209 const char *file4;
210 stringstream qy;
211 qy.str("");
212 //
213 // store the information about the directory where the program has been started
214 //
215 gDirectory->GetList()->Delete();
216 const char* startingdir = gSystem->WorkingDirectory();
217 TString path;
218 stringcopy(path,startingdir);
219 //
220 // Open files:
221 //
222 // the registry file
223 //
224 regFile = emigetFile(filename, "Physics", "Registry");
225 if ( !regFile ){
226 printf("\n\n ERROR: no registry file, exiting...\n\n");
227 return(103);
228 };
229 TTree *regTree = (TTree*)regFile->Get("Registry");
230 //
231 // the header file
232 //
233 headerFile = emigetFile(filename, "Physics", "Header");
234 if ( !headerFile ){
235 printf("\n\n ERROR: no header file, exiting...\n\n");
236 return(1);
237 };
238 TTree *otr = (TTree*)headerFile->Get("Pscu");
239 //
240 // the trigger file
241 //
242 acFile = emigetFile(filename, "Anticounter");
243 if ( !acFile ){
244 printf("\n\n ERROR: No Anticounter file! \n");
245 return(1);
246 } else {
247 otr->AddFriend("Anticounter",acFile);
248 };
249 //
250 // x Silvio: controlla che i vari branch si chiamino proprio Anticounter, non ho controllato...
251 //
252 // Define some basic variables
253 //
254 // determine the number of entries in the rootples
255 //
256 Long64_t nevents = otr->GetEntries();
257 //
258 // check if the file contains any event
259 //
260 if ( nevents < 1 ) {
261 printf(" ERROR: the file is empty!\n\n");
262 return(7);
263 };
264 //
265 // from the YODA filename determine the YODA processed level and check if we have a DW or dw unpacked file
266 //
267 const string myfil = (const char *)filename;
268 Int_t myposiz = myfil.find("dw_");
269 Int_t DW = 0;
270 if ( myposiz == -1 ) {
271 myposiz = myfil.find("DW_");
272 DW = 1;
273 };
274 if ( myposiz == -1 ) return(6);
275 TString basepath;
276 TString yodalev;
277 stringcopy(basepath,filename,0,myposiz);
278 stringcopy(yodalev,filename,myposiz+13,myposiz+15);
279 //
280 // since we are going to create an intermediate data output we create a rootple which is smaller respect to ntuple
281 //
282 filety = "root";
283 //
284 gSystem->ChangeDirectory(path);
285 //
286 // if no output is given put the output in the YODA directory structure under /Physics/Level2
287 //
288 if ( outDir == "" ){
289 file4 = filename;
290 file3.str("");
291 file3 << file4 << "/Physics/Level2";
292 } else {
293 file4 = outDir;
294 file3.str("");
295 file3 << file4;
296 }
297 //
298 file2.str("");
299 file2 << file3.str().c_str() << "/";
300 //
301 // the filename must contain the run ID number and the detector name. In the future it will be perhaps necessary to add also the software version which generated the file (?).
302 //
303 file2 << "run_id-" << run;
304 file2 << "_template."<< filety;
305 //
306 printf("\n Filename will be: \n %s \n\n",file2.str().c_str());
307 //
308 // check if Level2 directory exists, if not create it.
309 //
310 Int_t ERR = gSystem->MakeDirectory(file3.str().c_str());
311 //
312 if ( !ERR ) {
313 printf(" LEVEL2 directory doesn't exist! Creating LEVEL2 directory \n\n");
314 };
315 //
316 // Open file to write
317 //
318 //
319 // In any case we must create a tree
320 //
321 TTree *tree = 0;
322 //
323 gSystem->ChangeDirectory(path);
324 stringstream name;
325 TString tmpfile2;
326 stringcopy(tmpfile2,file2.str().c_str());
327 TFile *hfile = 0;
328 //
329 // To pass data from C to FORTRAN we need a structure, while to save data into a rootple
330 // the best is to have a class. Here the definition of both of them
331 // NOTICE that the name of the struct MUST be the same of the name of the fortran COMMON
332 // plus the underscore at the end.
333 //
334 AcLevel2 *acl2 = new AcLevel2();
335 //
336 // x Silvio: AcLevel2 e` il nome della classe, quello lo setti in ../inc/acl2class.h
337 // x Silvio: cambia tutti i ctempl2 in acl2 da qui in poi (OK!)
338 //
339
340 //
341 // here we open the rootple file
342 //
343 char *type;
344 type = "RECREATE";
345 hfile = new TFile(file2.str().c_str(),type,"AC data");
346 //
347 // and we book the rootple
348 //
349 tree = new TTree("AcLevel2","AC Level2 data");
350 //
351 // cambia tutte le variabili !! (Silvio)
352 //
353 tree->Branch("status",acl2->status,"status[2]/s");
354 tree->Branch("hitmap",acl2->hitmap,"hitmap[2]/s");
355 tree->Branch("hitstatus",acl2->hitstatus,"hitstatus[2]/s");
356 tree->Branch("trigger",acl2->trigger,"trigger[2]/s");
357 tree->Branch("OBT",&acl2->OBT,"OBT/l");
358 tree->Branch("pro_num",&acl2->pro_num,"pro_num/I");
359 tree->Branch("pkt_num",&acl2->pkt_num,"pkt_num/l");
360 //
361 // x Silvio qui ci metti le variabili della tua classe. In teoria puoi puntare direttamente alla classe senza elencare tutte le variabili ma il riempimento e` molto piu` lento
362 //
363 //
364 // here we have the software code branch, filled once...
365 //
366 TTree *software = 0;
367 software = new TTree("Software","Software used to generate data");
368 software->Branch("swcode",&swcode,"swcode/I");
369 //
370 // ...here:
371 //
372 software->Fill();
373 //
374 // run over all the events
375 //
376 printf("\n Processed events: \n\n");
377 //
378 // define pointers to the data branches
379 //
380 pamela::PscuHeader *ph = 0;
381 pamela::EventHeader *eh = 0;
382 pamela::RegistryEvent *reg=0;
383 pamela::anticounter::AnticounterEvent *ac = 0;
384 //
385 // x Silvio: controlla che sia giusta questa riga sopra e da qui in poi sostituisci tutti i trig con ac (OK, level 0)
386 //
387 regTree->SetBranchAddress("Registry", &reg);
388 otr->SetBranchAddress("Header", &eh);
389 otr->SetBranchAddress("Anticounter.Event", &ac);
390 //
391 // x Silvio: come sopra (OK)
392 //
393 // here we define the number of entry of the registry file
394 // notice that this number can be different from the number of entries of the physics files
395 // due to the possibility of multiple records of the same event in that files
396 //
397 Long64_t regevents = regTree->GetEntries();
398 //
399 // definition of the first event minus one, this is done this way for synchronization purpose
400 //
401 Int_t i = -1;
402 //
403 // some variable definition
404 //
405 Int_t procerr = 0;
406 Int_t retval = 0;
407 Int_t pktnum = 0;
408 Int_t obt = 0;
409 ULong64_t atime = 0ULL;
410 //
411 //
412 // start looping on events cointained in the registry file
413 //
414 for (Int_t re=0; re<regevents; re++){
415 //
416 if ( re%1000 == 0 && re > 0 ) printf(" %iK \n",re/1000);
417 //
418 // get the first registry entry
419 //
420 regTree->GetEntry(re);
421 //
422 // read the event absolute time
423 //
424 atime = reg->absTime;
425 //
426 // introduce here declaration to avoid cross-inizialitation (see commented declaration ~35 lines below)
427 UShort_t cnt = 1;
428 //
429 // if this event does not belong to the processed run jump to the next event (goto ev_jump)
430 //
431 if ( atime > runtrailtime || atime < runheadtime ) {
432 goto ev_jump;
433 };
434 //
435 // ok, this event belongs to the processed run
436 //
437 i = reg->event;
438 //
439 // get data
440 //
441 otr->GetEntry(i);
442 //
443 // retrieve obt and packet number for this entry
444 //
445 ph = eh->GetPscuHeader();
446 pktnum = ph->GetCounter();
447 obt = ph->GetOrbitalTime();
448 //
449 procev++;
450 //
451 // ok, here we must do our data processing
452 //
453 //
454 // x Silvio: qui ci metti il codice in c++ scritto da Petter per il software a terra
455 // ti riferisci alle variabili che vuoi salvare sempre con acl2->
456 //
457 // --------------------------------------------------------------
458 //
459 // declare and initialize variables ...
460 //
461 //UShort_t cnt = 1; // move to line 428 before goto and it works
462 //
463 UShort_t CRCcheck[2];
464 UShort_t Status[2];
465 //UShort_t Reg[2][6];
466 //UShort_t Header[2][2];
467 UShort_t Hitmap[2];
468 UShort_t Hitstatus[2];
469 UShort_t Trigger[2];
470 UShort_t Trigg_old[2];
471 //UShort_t Temp[2][4];
472
473 //UShort_t timeLSB[2];
474 //UShort_t timeMSB[2];
475 //Double_t time[2];
476
477 UShort_t Counter[2][16];
478 UShort_t Counter_old[2][16];
479
480 UShort_t Shift[2][16];
481 Int_t vec[2][16][16];
482 //UShort_t pack[2][16];
483 Int_t tmp;
484 //
485 for(Int_t gg = 0; gg < 2; gg++)
486 {
487 Hitstatus[gg] = 0;
488 }
489 //
490 //fetch data
491 //
492 for(Int_t j = 0; j < 2; j++)
493 {
494 if(re>0)
495 {
496 Trigg_old[j] = Trigger[j];
497 for(Int_t jj = 0; jj < 16; jj++)
498 Counter_old[j][jj] = Counter[j][jj];
499 }
500
501 for(Int_t jj = 0; jj < 16; jj++)
502 {
503 Counter[j][jj] = ac->counters[j][jj];
504 Shift[j][jj] = ac->shift[j][jj];
505 }
506 Status[j] = ac->status[j];
507 Hitmap[j] = ac->hitmap[j];
508 Trigger[j] = ac->trigg[j];
509 CRCcheck[j] = ac->CRCcheck[j];
510 }
511
512 /***********************************************/
513
514
515 //process data
516 /***********************************************/
517 //shiftregisters
518 cnt=1;
519 for(Int_t b = 0; b < 2; b++){ //card
520 for(Int_t k = 0; k < 16; k++){ //shift register
521 for(Int_t l = 0; l < 16; l++){ //bin
522 tmp = ((Shift[b][k] & cnt) > 0 ? 1 : 0);
523 vec[b][k][l]=tmp;
524 cnt=cnt<<1;
525 }
526 cnt=1;
527 }
528 }
529
530 //fill Level1 file
531 /************************************************/
532 for(Int_t s = 0; s < 2; s++)
533 {
534 acl2->hitmap[s] = Hitmap[s];
535 acl2->trigger[s] = Trigger[s];
536 cnt=1;
537 for(Int_t k = 0; k < 16; k++)
538 {
539 for(Int_t bin = 5; bin < 9; bin++) //acceptance window
540 {
541 if(vec[s][bin][k]==1)
542 Hitstatus[s] = Hitstatus[s] | cnt;
543 }
544 cnt=cnt<<1;
545 }
546 acl2->hitstatus[s] = Hitstatus[s];
547
548 //Status
549 /****************************************/
550 if(s==0){
551 if(Trigger[0] != (Trigg_old[0]+1))
552 acl2->status[0] = acl2->status[0] | 0x1;
553 if(Status[0] & 0x001F < 0x001F)
554 acl2->status[0] = acl2->status[0] | 0x2;
555 if(CRCcheck[0] == 0)
556 acl2->status[0] = acl2->status[0] | 0x4;
557 for(Int_t gg = 0; gg < 16; gg++){
558 if((Counter[0][gg] == Counter_old[0][gg]) && i>0){
559 acl2->status[0] = acl2->status[0] | 0x8;
560 }
561 }
562 }
563 else{
564 if(Trigger[1] != (Trigg_old[1]+2))
565 acl2->status[1] = acl2->status[1] | 0x1;
566 if(Status[1] & 0x001F < 0x001F)
567 acl2->status[1] = acl2->status[1] | 0x2;
568 if(CRCcheck[1] == 0)
569 acl2->status[1] = acl2->status[1] | 0x4;
570 for(Int_t gg = 0; gg < 16; gg++){
571 if((Counter[1][gg] == Counter_old[1][gg]) && i>0)
572 acl2->status[1] = acl2->status[1] | 0x8;
573 }
574 }
575
576 /****************************************/
577 }
578 acl2->OBT = ph->GetOrbitalTime();
579 acl2->pkt_num = re+1;
580 acl2->pro_num = ph->GetCounter();
581 /************************************************/
582
583
584 // --------------------------------------------------------------
585
586 //
587 // fill the rootple
588 //
589 tree->Fill();
590 //
591 // finished processing, go to the next event.
592 //
593 ev_jump:
594 //
595 // if any error occurred during processing try to save and close the files anyway.
596 //
597 if ( procerr ){
598 if ( i > 0 ){
599 printf("\nTrying to close the file anyway...\n");
600 hfile->Write();
601 hfile->Close();
602 printf("...done!\n");
603 } else {
604 //
605 // if no data were processed remove the file from disk
606 //
607 printf("\nERROR: output file is empty! \n");
608 stringstream remfile;
609 remfile.str("");
610 remfile << "rm -f " << file2.str().c_str();
611 gSystem->Exec(remfile.str().c_str());
612 };
613 printf("\nERROR: processed %i out of %i events\n",i,(int)nevents);
614 printf("\nERROR: an error occurred during processing!\n\n Exiting...\n\n");
615 goto theend;
616 };
617 };
618 //
619 // close rootple files and exit
620 //
621 hfile->Write();
622 hfile->Close();
623 printf("\n");
624 printf(" Finished! Processed events: %i \n\n Exiting... \n",procev);
625 printf("\n");
626 theend:
627 //
628 // go back to the starting path and unload fortran libraries
629 //
630 gSystem->ChangeDirectory(path);
631 return(retval);
632 }
633

  ViewVC Help
Powered by ViewVC 1.1.23