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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide 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 pam-se 1.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