/[PAMELA software]/calo/ground/LEVEL2/macros/CaloLEVEL2.cra
ViewVC logotype

Annotation of /calo/ground/LEVEL2/macros/CaloLEVEL2.cra

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Mon Dec 5 16:13:53 2005 UTC (19 years, 2 months ago) by mocchiut
Branch: LEVEL2
CVS Tags: v4r00, start
Changes since 1.1: +0 -0 lines
Imported sources

1 mocchiut 1.1 //
2     // Given a calibration and a data file this program create an ntuple with LEVEL2 calorimeter variables - Emiliano Mocchiutti
3     //
4     // CaloLEVEL2.c version 4.00 (2005-11-29)
5     //
6     // The only input needed is the path to the directory created by YODA for the data file you want to analyze.
7     //
8     // Changelog:
9     //
10     // 3.17 - 4.00 (2005-11-29): preparing for the final release, small changes in the makefile.
11     //
12     // 3.16 - 3.17 (2005-10-13): do not exclude all strips marked "bad" but only the ones with RMS>6.
13     //
14     // 3.15 - 3.16 (2005-10-04): now it reads Tracker Data version 2.00. Process self-trigger events assuming R=1000 GV and using the calorimeter fit of the track;
15     // self-trigger events are still marked as "bad" and there is no energy info, nor any angular correction on the track.
16     //
17     // 3.14 - 3.15 (2005-09-20): QX22 and NX22 are not related to the last plane but the 11th! fixed!
18     //
19     // 3.13 - 3.14 (2005-09-09): Still some bugs in 64 bit arch, fixed.
20     //
21     // 3.12 - 3.13 (2005-09-08): Bug: variable cbar is filled incorrectly. Fixed.
22     //
23     // 3.11 - 3.12 (2005-08-31): Added two variables in the ntuple (thex, they sin of x/y angles as measured by the calorimeter)
24     //
25     // 3.10 - 3.11 (2005-08-30): Don't run in MySQL standalone mode even if in FORCE mode.
26     //
27     // 3.09 - 3.10 (2005-08-05): 64 bit arch bugs fixed.
28     //
29     // 3.08 - 3.09 (2005-07-27): x and y view were inverted when fitting the track using calorimeter information, fixed. (affect variables cbar, cibar, npcfit and varcfit).
30     //
31     // 3.07 - 3.08 (2005-07-25): planetot and qmean always empty in rootple, fixed. Added variables on goodness of the calorimeter fit (npcfit, varcfit).
32     //
33     // 3.06 - 3.07 (2005-07-21): If not in FORCE mode abort processing if no calibration file is found.
34     //
35     // 3.05 - 3.06 (2005-07-19): Infinite loop in case of OBT jumps, fixed. NaN in the ntuple and rootple in some cases, fixed (lowered optimization level in
36     // the compilation of the ntuple and check for NaN from the tracking routine). Added variables in the ntuple/rootple: planetot, qmean.
37     //
38     // 3.04 - 3.05 (2005-07-18): No more class CalorimeterLevel2 since it makes very slow the filling of the rootple (unknown reason). Now the tree is filled
39     // with variables from structure CaLevel2.
40     //
41     // 3.03 - 3.04 (2005-07-15): More tuning, now it can compile and it works! Not loaded anymore yodaUtility.c, use clone routines instead.
42     //
43     // 3.02 - 3.03 (2005-07-13): Start the tuning for compilation and speed optimization.
44     //
45     // 3.01 - 3.02 (2005-07-11): Correction for high energy deposit. Correction in the calorimeter planes position.
46     //
47     // 3.00 - 3.01 (2005-07-11): Infinite loop in some cases of lost of sync, fixed. Cannot write ntuples bigger than 3 MB, fixed.
48     //
49     // 2.10 - 3.00 (2005-07-05): Include track information.
50     //
51     // 2.09 - 2.10 (2005-07-05): Started the tuning for the compilation of the macro. Added include files of ROOT subroutines.
52     //
53     // 2.08 - 2.09 (2005-06-29): Simplification in the check of goodness for the events. Integer convention for software code.
54     //
55     // 2.07 - 2.08 (2005-06-24): Changed the structure of the release and added the enviromental variables.
56     //
57     // 2.06 - 2.07 (2005-06-23): Added boolean variables to ntuple: good, crc, perr, swerr.
58     //
59     // 2.05 - 2.06 (2005-06-22): Trigger flag set as default to -1. It will remain -1 if the program is not able to determine who gave the trigger.
60     // Added: new tree (Software)/ntuple (number 2) with one entry only which contains the information about the software used to
61     // produce the LEVEL2 data. Changed names in commons/struct/ntuples/rootples: evfile -> pro_num, headc -> pkt_num .
62     // Disabled MySQL input variable (still there for testing purpouse).
63     //
64     // 2.04 - 2.05 (2005-06-17): Trigger flag is not always correct for S4 triggers, fixed with a workaround.
65     //
66     // 2.03 - 2.04 (2005-06-16): Bug, when casting from integer to float the resulting float is WRONG! changed in the rootple OBT, HEADC, EVFILE from real to integer.
67     // Changes also in the fortran libraries and passing structures.
68     //
69     // 2.02 - 2.03 (2005-06-15): Bug, in PAW and FORCE mode the existence of file is not checked, finally fixed.
70     //
71     // 2.01 - 2.02 (2005-06-07): Bug, the calibration used could be incorrect for some events. Fixed: check that the event time is in the calibration window,
72     // if not search for another calibration.
73     //
74     // 2.00 - 2.01 (2005-06-07): Check if file is older than 050515_007.
75     //
76     // 1.04 - 2.00 (2005-06-06): Added support to MySQL Trieste database.
77     //
78     // 1.03 - 1.04 (2005-05-31): Now it can unload all the libraries and you can run this script as many times you want without quitting ROOT.
79     //
80     // 1.02 - 1.03 (2005-05-27): cannot find libraries and calibration file when running twice in the same ROOT section, fixed. KNOWN BUG: cannot unload cfillcalol2_C!!
81     //
82     // 1.01 - 1.02 (2005-05-18): filename corruption in some cases. Workaround for the moment...
83     //
84     // 1.00 - 1.01 (2005-05-17): problems when starting with a possible backup calibration file. Fixed.
85     //
86     // 0.01 - 1.00 (2005-05-06): first time it works. evfile must be i+1 not i.
87     //
88     // 0.01 (2005-05-04): created.
89     //
90     #include <fstream>
91     #include <TTree.h>
92     #include <TClassEdit.h>
93     #include <TObject.h>
94     #include <TList.h>
95     #include <TSystem.h>
96     #include <TSystemDirectory.h>
97     #include <TString.h>
98     #include <TFile.h>
99     #include <TClass.h>
100     #include <TCanvas.h>
101     #include <TH1.h>
102     #include <TH1F.h>
103     #include <TH2D.h>
104     #include <TLatex.h>
105     #include <TPad.h>
106     //
107     #include <calomysqlstruct.h>
108     #include <ccal2struct.h>
109     #include <ctrkstruct.h>
110     #include <ctrkinclude.h>
111     //
112     #if !defined (__CINT__)
113     extern short int creadmysql(char *, const char *, struct Mystruct &);
114     extern void cfillcalol2(struct CaLevel1 &, struct CaLevel2 &);
115     extern void creadB(const char []);
116     extern void coptrklev2(const char [], struct Tracklev2 &, int &);
117     extern void copencalol2(char [], Int_t &, Int_t &, int &);
118     extern void cretrklev2(int &, Tracklev2 &);
119     extern void cclosecalol2(char []);
120     extern void ccltrklev2(struct Tracklev2 &);
121     extern const char *pathtocalibration();
122     #include <event/PamelaRun.h>
123     #include <event/physics/calorimeter/CalorimeterEvent.h>
124     #include <event/physics/trigger/TriggerEvent.h>
125     #include <event/CalibCalPedEvent.h>
126     #endif
127     //
128     #include <caloclasses.h>
129     #include <CaloFunctions.h>
130     //
131     #if !defined (__CINT__)
132     void unloadf77lib(TString Framework, Int_t MySQL){}
133     #endif
134     #if defined (__CINT__)
135     void unloadf77lib(TString Framework, Int_t MySQL){
136     const char *franame = Framework;
137     //
138     char *pam_lib=gSystem->Getenv("PAM_LIB");
139     if ( MySQL == 1 || MySQL == -1 ){
140     stringstream libload;
141     libload.str("");
142     libload << pam_lib << "/libreadmy";
143     gSystem->Unload(libload.str().c_str());
144     libload.str("");
145     libload << pam_lib << "/libreadmy_C";
146     gSystem->Unload(libload.str().c_str());
147     };
148     if ( !strcmp(franame,"paw") ) {
149     stringstream libload;
150     libload.str("");
151     libload << pam_lib << "/libclcalol2_C";
152     gSystem->Unload(libload.str().c_str());
153     libload.str("");
154     libload << pam_lib << "/libclcalol2";
155     gSystem->Unload(libload.str().c_str());
156     libload.str("");
157     libload << pam_lib << "/libopcalol2_C";
158     gSystem->Unload(libload.str().c_str());
159     libload.str("");
160     libload << pam_lib << "/libopcalol2";
161     gSystem->Unload(libload.str().c_str());
162     };
163     stringstream libload;
164     libload.str("");
165     libload << pam_lib << "/libcrcalol2_C";
166     gSystem->Unload(libload.str().c_str());
167     char *pam_lib=gSystem->Getenv("PAM_LIB");
168     stringstream libload;
169     libload.str("");
170     libload << pam_lib << "/libcrcalol2";
171     gSystem->Unload(libload.str().c_str());
172     char *pam_lib=gSystem->Getenv("PAM_LIB");
173     stringstream libload;
174     libload.str("");
175     libload << pam_lib << "/libreadb2maps_C";
176     gSystem->Unload(libload.str().c_str());
177     char *pam_lib=gSystem->Getenv("PAM_LIB");
178     stringstream libload;
179     libload.str("");
180     libload << pam_lib << "/libreadb2maps";
181     gSystem->Unload(libload.str().c_str());
182     char *pam_lib=gSystem->Getenv("PAM_LIB");
183     stringstream libload;
184     libload.str("");
185     libload << pam_lib << "/libtrack";
186     gSystem->Unload(libload.str().c_str());
187     }
188     #endif
189    
190     short int calolevel2core(Int_t ei, Int_t b[4] , TTree *otr, TTree *tree, CaLevel2 & clevel2, Evento & evento, Calib & calib, TString calcalibfile){
191     //
192     // Define variables
193     //
194     struct CaLevel1 clevel1;
195     clevel1.xalig = 120.4;
196     clevel1.yalig = 118.6;
197     clevel1.zalig = -260.;
198     struct Mystruct mystruct;
199     clevel1.paw = calib.ispaw;
200     Int_t upnn ;
201     Int_t done = 0;
202     Int_t pre = -1;
203     Int_t etime;
204     Int_t S3 = 0;
205     Int_t S2 = 0;
206     Int_t S12 = 0;
207     Int_t S11 = 0;
208     Int_t se = 5;
209     bool isCOMP = 0;
210     bool isFULL = 0;
211     bool isRAW = 0;
212     Float_t ener, basel;
213     Float_t estripnull[96][22][2];
214     //
215     Int_t doneb = 0;
216     Int_t nn;
217     Int_t ck = 0;
218     Int_t ipre = 0;
219     Int_t ip[3] = {0};
220     Float_t base0, base1, base2;
221     base0 = 0.;
222     base1 = 0.;
223     base2 = 0.;
224     Float_t qpre[6] = {0};
225     Float_t ene[96] = {0};
226     //
227     pamela::calorimeter::CalorimeterEvent *de = 0;
228     pamela::PscuHeader *ph = 0;
229     pamela::EventHeader *eh = 0;
230     pamela::trigger::TriggerEvent *trig = 0;
231     otr->SetBranchAddress("Header", &eh);
232     otr->SetBranchAddress("Calorimeter.Event", &de);
233     otr->SetBranchAddress("Trigger.Event", &trig);
234     //
235     // other variables
236     //
237     char *yodala2;
238     Int_t chdone[4] = {0,0,0,0};
239     Int_t pe = 0;
240     Float_t tmptrigty = -1.;
241     //
242     // start processing
243     //
244     //
245     // read from the header of the event the OBT and event number
246     //
247     otr->GetEntry(ei);
248     ph = eh->GetPscuHeader();
249     clevel2.pkt_num = ph->GetCounter();
250     etime = ph->GetOrbitalTime();
251     clevel2.OBT = etime;
252     clevel2.pro_num = (ei+1);
253     //
254     // determine who generate the trigger for this event (TOF, S4/PULSER, CALO)
255     //
256     S3 = 0;
257     S2 = 0;
258     S12 = 0;
259     S11 = 0;
260     S3 = trig->patterntrig[2];
261     S2 = trig->patterntrig[3];
262     S12 = trig->patterntrig[4];
263     S11 = trig->patterntrig[5];
264     if ( trig->patterntrig[0] ) tmptrigty = 2.;
265     if ( S3 || S2 || S12 || S11 ) tmptrigty = 0.;
266     if ( trig->patterntrig[1] & (1<<0) || (!trig->patterntrig[0] && !S3 && !S2 && !S12 && !S11) ) tmptrigty = 1.;
267     clevel2.trigty = tmptrigty;
268     //
269     yodala2 = "";
270     if ( calib.mysql ){
271     if ( !calib.obtjump ){
272     for (Int_t s = 0; s < 4; s++){
273     if ( etime > calib.time[s][0] ){
274     tryagain:
275     printf(" ** SECTION %i **\n",s);
276     stringstream qy;
277     qy.str("");
278     Int_t chkmysql = 0;
279     stringstream myfile;
280     const char *tabula = calib.tablename;
281     qy << "select * from " << tabula;
282     qy << " where section = " << s;
283     qy << " and fromtime >= " << calib.time[s][0];
284     qy << " limit 1";
285     chkmysql = creadmysql(calib.db,(char *)qy.str().c_str(),mystruct);
286     if ( chkmysql ) {
287     printf("- ERROR: problems querying MySQL database! -\n- ERROR: Empty table? -\n\n");
288     return(1);
289     } else {
290     if ( mystruct.totime == -1 ){
291     printf("- ERROR: problems querying MySQL database! -\n- ERROR: Corrupted table? -\n\n");
292     return(2);
293     };
294     calib.time[s][0] = mystruct.totime;
295     calib.time[s][1] = mystruct.fromtime;
296     calib.ttime[s][0] = mystruct.calibtime;
297     const char *caliba = calib.basepath;
298     const char *yodala = calib.yodalev;
299     yodala2 = (char *)yodala;
300     const char *fila = mystruct.calcalibfile;
301     myfile.str("");
302     myfile << caliba;
303     if ( !calib.DW ) {
304     myfile << fila;
305     } else {
306     myfile << "DW";
307     string myfil = (const char*)fila;
308     Int_t myposiz = myfil.find("dw_");
309     TString mytmp;
310     Int_t mspos = myposiz+2;
311     Int_t mepos = myposiz+13;
312     stringcopy(mytmp,fila,mspos,mepos);
313     const char *myf = mytmp;
314     myfile << mytmp;
315     };
316     myfile << yodala << "/";
317     calib.obtjump = mystruct.obtjump;
318     printf(" - event at time %i. From time %i to time %i \n use calibration at time %i, file %s\n",etime,calib.time[s][1],calib.time[s][0],calib.ttime[s][0],myfile.str().c_str());
319     if ( calib.obtjump ) printf("\n WARNING: jump in the OBT! this calibration might not be the best! \n");
320     if ( etime > calib.time[s][0] && !calib.obtjump ){
321     printf("\n WARNING: event at time greater than the upper limit for this calibration!\nSearch for the correct calibration\n");
322     goto tryagain;
323     };
324     };
325     Int_t pedeerr = 0;
326     pedeerr = CaloPede(myfile.str().c_str(),s,calib.ttime[s][0],calib);
327     if ( pedeerr ) {
328     printf("\n\nERROR: problems opening calibration file! \n\nERROR: YODA version of the calibration file is not %s? \n",yodala2);
329     return(1);
330     };
331     printf("\n");
332     };
333     };
334     };
335     } else {
336     //
337     // for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
338     //
339     if ( !calib.obtjump ) {
340     for (Int_t s = 0; s < 4; s++){
341     if ( calib.ttime[s][b[s]] ){
342     while ( etime > calib.time[s][b[s]+1] ){
343     printf(" CALORIMETER: \n" );
344     printf(" - Section %i, event at time %i while old calibration time limit at %i. Use new calibration at time %i -\n",s,etime,calib.time[s][b[s]+1],calib.ttime[s][b[s]+1]);
345     printf(" END CALORIMETER. \n\n" );
346     b[s]++;
347     TString pfile;
348     if ( calib.fcode[s][b[s]] != 10 ){
349     TString file2f = "";
350     stringcopy(file2f,calcalibfile);
351     TString pfile = whatnamewith(file2f,calib.fcode[s][b[s]]);
352     } else {
353     pfile = (TString)calcalibfile;
354     };
355     CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
356     };
357     };
358     };
359     };
360     };
361     //
362     // run over views and planes
363     //
364     for (Int_t l = 0; l < 2; l++){
365     for (Int_t m = 0; m < 22; m++){
366     //
367     // determine the section number
368     //
369     se = 5;
370     if (l == 0 && m%2 == 0) se = 3;
371     if (l == 0 && m%2 != 0) se = 2;
372     if (l == 1 && m%2 == 0) se = 1;
373     if (l == 1 && m%2 != 0) se = 0;
374     //
375     // determine what kind of event we are going to analyze
376     //
377     isCOMP = 0;
378     isFULL = 0;
379     isRAW = 0;
380     if ( de->stwerr[se] & (1 << 16) ) isCOMP = 1;
381     if ( de->stwerr[se] & (1 << 17) ) isFULL = 1;
382     if ( de->stwerr[se] & (1 << 3) ) isRAW = 1;
383     if ( !chdone[se] ){
384     //
385     // check for any error in the event
386     //
387     clevel2.crc[se] = 0;
388     if ( de->perror[se] == 132 ){
389     clevel2.crc[se] = 1;
390     pe++;
391     };
392     clevel2.perr[se] = 0;
393     if ( de->perror[se] != 0 ){
394     clevel2.perr[se] = 1;
395     pe++;
396     };
397     clevel2.swerr[se] = 0;
398     for (Int_t j = 0; j < 7 ; j++){
399     if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
400     clevel2.swerr[se] = 1;
401     pe++;
402     };
403     };
404     chdone[se] = 1;
405     };
406     if ( clevel2.crc[se] == 0 && calib.good2 == 1 ){
407     pre = -1;
408     if ( isRAW ){
409     for (Int_t nn = 0; nn < 96; nn++){
410     evento.dexy[l][m][nn] = de->dexy[l][m][nn] ;
411     evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
412     };
413     };
414     //
415     // run over preamplifiers
416     //
417     pre = -1;
418     for (Int_t i = 0; i < 3; i++){
419     for (Int_t j = 0; j < 2; j++){
420     pre = j + i*2;
421     //
422     // baseline check and calculation
423     //
424     if ( !isRAW ) {
425     evento.base[l][m][pre] = de->base[l][m][pre] ;
426     } else {
427     //
428     // if it is a raw event and we haven't checked
429     // yet, calculate the baseline.
430     //
431     CaloFindBaseRaw(calib,evento,l,m,pre);
432     };
433     };
434     };
435     //
436     // run over strips
437     //
438     pre = -1;
439     for (Int_t i = 0 ; i < 3 ; i++){
440     ck = 0;
441     ip[i] = 0;
442     for (Int_t n = i*32 ; n < (i+1)*32 ; n++){
443     if (n%16 == 0) {
444     done = 0;
445     doneb = 0;
446     pre++;
447     qpre[pre] = 0;
448     };
449     //
450     // baseline check and calculation
451     //
452     if ( !isRAW ) {
453     evento.dexy[l][m][n] = de->dexy[l][m][n] ;
454     evento.dexyc[l][m][n] = de->dexyc[l][m][n] ;
455     };
456     //
457     // no suitable new baseline, use old ones and compress data!
458     //
459     if ( !done ){
460     if ( (evento.base[l][m][pre] == 31000. || evento.base[l][m][pre] == 0.) ){
461     ck = 1;
462     if (pre%2 == 0)
463     ip[i] = pre + 1;
464     else
465     ip[i] = pre - 1;
466     if ( (evento.base[l][m][ip[i]] == 31000. || evento.base[l][m][ip[i]] == 0.) ){
467     //
468     evento.base[l][m][pre] = calib.sbase[l][m][pre];
469     //
470     ck = 2;
471     upnn = n+16;
472     if ( upnn > 96 ) nn = 96;
473     for ( Int_t nn = n; nn<upnn; nn++ ){
474     evento.dexyc[l][m][nn] = de->dexyc[l][m][nn];
475     };
476     CaloCompressData(calib,evento,l,m,pre);
477     done = 1;
478     };
479     };
480     };
481     //
482     // CALIBRATION ALGORITHM
483     //
484     basel = evento.base[l][m][pre];
485     if ( !doneb ){
486     switch (ck) {
487     case 0:
488     base0 = evento.base[l][m][pre];
489     base2 = calib.calbase[l][m][pre];
490     break;
491     case 1:
492     base0 = evento.base[l][m][ip[i]];
493     base2 = calib.calbase[l][m][ip[i]];
494     break;
495     case 2:
496     base0 = calib.sbase[l][m][pre];
497     base2 = calib.calbase[l][m][pre];
498     break;
499     };
500     base1 = calib.calbase[l][m][pre];
501     doneb = 1;
502     };
503     ener = evento.dexyc[l][m][n];
504     clevel1.estrip[n][m][l] = 0.;
505     if ( base0>0 && base0 < 30000. && ener > 0. ){
506     clevel1.estrip[n][m][l] = (ener - calib.calped[l][m][n] - base0 - base1 + base2)/calib.mip[l][m][n] ;
507     //
508     // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)
509     //
510     // if ( clevel1.estrip[n][m][l] < evento.emin || calib.calgood[l][m][n] != 0 ) {
511     if ( clevel1.estrip[n][m][l] < evento.emin || calib.calrms[l][m][n] > 26 ) {
512     clevel1.estrip[n][m][l] = 0.;
513     };
514     qpre[pre] += clevel1.estrip[n][m][l];
515     };
516     calib.sbase[l][m][pre] = evento.base[l][m][pre];
517     };
518     if (ck == 1){
519     if (ip[i]%2 == 0)
520     ipre = ip[i] + 1;
521     else
522     ipre = ip[i] - 1;
523     for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
524     clevel1.estrip[j][m][l] = clevel1.estrip[j][m][l] + (qpre[ipre] - qpre[ip[i]]) * 0.00478;
525     };
526     };
527     if (ck == 2){
528     for (Int_t j = i*32 ; j < (i+1)*32 ; j++){
529     ipre = j/16 + 1;
530     clevel1.estrip[j][m][l] = clevel1.estrip[j][m][l] + qpre[ipre] * 0.00478;
531     };
532     };
533     };
534     //
535     Int_t jj = -2;
536     Int_t jjj = -3;
537     for (Int_t j = 0 ; j < 98 ; j++){
538     jj++;
539     jjj++;
540     if ( j < 96 ) ene[j] = clevel1.estrip[j][m][l];
541     if ( jj >= 0 && jj < 96 ){
542     if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] = ene[jj-1] - clevel1.estrip[jj][m][l] * 0.01581;
543     if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] = ene[jj+1] - clevel1.estrip[jj][m][l] * 0.01581;
544     };
545     if ( jjj >= 0 && jjj < 96 ){
546     if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1.estrip[jjj-1][m][l] = clevel1.estrip[jjj-1][m][l] - ene[jjj] * 0.01581;
547     if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1.estrip[jjj+1][m][l] = clevel1.estrip[jjj+1][m][l] - ene[jjj] * 0.01581;
548     };
549     };
550     //
551     } else {
552     for (Int_t nn = 0; nn < 96; nn++){
553     clevel1.estrip[nn][m][l] = 0.;
554     };
555     };
556     };
557     };
558     //
559     //
560     //
561     if ( !pe && calib.good2 ){
562     clevel2.good = 1;
563     } else {
564     clevel2.good = 0;
565     };
566     //
567     // all we need from the tracker is the state vector al_p and the boolean for good events:
568     //
569     if ( calib.good2 ){
570     for (Int_t e = 0; e<5 ; e++){
571     clevel1.al_p[e][0] = calib.al_p[e][0];
572     clevel1.al_p[e][1] = calib.al_p[e][1];
573     };
574     clevel1.good2 = 1;
575     } else {
576     for (Int_t e = 0; e<5 ; e++){
577     clevel1.al_p[e][0] = 0.;
578     clevel1.al_p[e][1] = 0.;
579     };
580     clevel1.good2 = 0;
581     };
582     //
583     // for each event generate level2 data (clevel2) calling the fortran routine and passing calibrated data (clevel1):
584     //
585     cfillcalol2(clevel1,clevel2);
586     //
587     // save in the class the level2 calorimeter informations
588     //
589     //
590     // fill the rootple if not in paw mode
591     //
592     if ( clevel1.paw == 0. ) tree->Fill();
593     //
594     // clear estrip values.
595     //
596     memcpy(clevel1.estrip, estripnull, sizeof(estripnull));
597     //
598     return(0);
599     }
600     //
601     short int CaloLEVEL2(TString filename, TString TrackerDir="", TString outDir ="", TString Framework = "root", Int_t FORCE = 0){
602     //
603     // define variables
604     //
605     Int_t swcode = 04;
606     Int_t swtrkcode = 2;
607     Int_t MySQL = 1;
608     Int_t debug = 0;
609     //
610     #if defined(__CINT__)
611     if ( !gSystem->Getenv("PAM_INC") || !gSystem->Getenv("PAM_LIB") || !gSystem->Getenv("PAM_CALIB") ){
612     printf("\n ERROR: you first must set the PAMELA enviromental variables\n read the README file \n");
613     printf("\n Exiting... \n");
614     return(11);
615     };
616     emicheckLib();
617     #endif
618     Int_t trkdir = 1;
619     TString trackerdir;
620     if ( TrackerDir=="" ){
621     trackerdir = filename;
622     trkdir = 0;
623     } else {
624     trackerdir = TrackerDir;
625     };
626     TString calcalibfile;
627     //
628     // other variables
629     //
630     TFile *headerFile;
631     TFile *caloFile;
632     TFile *triggerFile;
633     TString filety;
634     stringstream file2;
635     stringstream file3;
636     const char *file4;
637     Int_t okcalo = 0;
638     TTree *ctree = 0;
639     TFile *chfile;
640     Long64_t cnevents=0;
641     stringstream qy;
642     qy.str("");
643     Int_t chkmysql = 0;
644     Int_t imtrack;
645     //
646     // load libraries to determine the track information.
647     //
648     #if defined (__CINT__)
649     const char *pam_lib=gSystem->Getenv("PAM_LIB");
650     stringstream libload;
651     libload.str("");
652     libload << pam_lib << "/libtrack";
653     gSystem->Load(libload.str().c_str());
654     libload.str("");
655     libload << pam_lib << "/libreadb2maps";
656     gSystem->Load(libload.str().c_str());
657     libload.str("");
658     libload << pam_lib << "/libreadb2maps_C";
659     gSystem->Load(libload.str().c_str());
660     //
661     // open libraries to extract calorimeter level2 data
662     //
663     libload.str("");
664     libload << pam_lib << "/libcrcalol2";
665     gSystem->Load(libload.str().c_str());
666     libload.str("");
667     libload << pam_lib << "/libcrcalol2_C";
668     gSystem->Load(libload.str().c_str());
669     #endif
670     const char *franame = Framework;
671     #if defined (__CINT__)
672     if ( !strcmp(franame,"paw") ) {
673     libload.str("");
674     libload << pam_lib << "/libopcalol2";
675     gSystem->Load(libload.str().c_str());
676     libload.str("");
677     libload << pam_lib << "/libopcalol2_C";
678     gSystem->Load(libload.str().c_str());
679     libload.str("");
680     libload << pam_lib << "/libclcalol2";
681     gSystem->Load(libload.str().c_str());
682     libload.str("");
683     libload << pam_lib << "/libclcalol2_C";
684     gSystem->Load(libload.str().c_str());
685     };
686     #endif
687     if ( MySQL != 1 ){
688     printf("\nWARNING: you MUST run this script using the MySQL database. Reverting MySQL to 1 \n\n");
689     MySQL = 1;
690     };
691     //MySQL = 0;
692     gDirectory->GetList()->Delete();
693     const char* startingdir = gSystem->WorkingDirectory();
694     TString path;
695     stringcopy(path,startingdir);
696     //
697     #if defined (__CINT__)
698     Int_t MySQLload = 0;
699     #endif
700     //
701     // Open files
702     //
703     headerFile = emigetFile(filename, "Physics", "Header");
704     if ( !headerFile ){
705     printf("\n\n ERROR: no header file, exiting...\n\n");
706     #if defined (__CINT__)
707     gSystem->ChangeDirectory(path);
708     unloadf77lib(Framework,MySQLload);
709     #endif
710     return(1);
711     };
712     caloFile = emigetFile(filename, "Calorimeter");
713     if ( !caloFile ){
714     printf("\n\n ERROR: no calorimeter file! \n\n");
715     #if defined (__CINT__)
716     gSystem->ChangeDirectory(path);
717     unloadf77lib(Framework,MySQLload);
718     #endif
719     return(2);
720     };
721     //
722     TTree *otr = (TTree*)headerFile->Get("Pscu");
723     otr->AddFriend("Calorimeter",caloFile);
724     triggerFile = emigetFile(filename, "Trigger");
725     if ( !triggerFile ){
726     printf("\n\n ERROR: No trigger file! \n");
727     #if defined (__CINT__)
728     gSystem->ChangeDirectory(path);
729     unloadf77lib(Framework,MySQLload);
730     #endif
731     return(1);
732     } else {
733     otr->AddFriend("Trigger",triggerFile);
734     };
735     //
736     // Define some basic variables
737     //
738     struct Evento evento;
739     struct Calib calib;
740     struct Mystruct mystruct;
741     evento.emin = 0.7;
742     Long64_t nevents = otr->GetEntries();
743     Int_t snevents = (Int_t)nevents;
744     if ( nevents < 1 ) {
745     printf(" ERROR: the file is empty!\n\n");
746     #if defined (__CINT__)
747     gSystem->ChangeDirectory(path);
748     unloadf77lib(Framework,MySQLload);
749     #endif
750     return(7);
751     };
752     Int_t b[4];
753     Int_t paw = 0;
754     //
755     const string myfil = (const char *)filename;
756     Int_t myposiz = myfil.find("dw_");
757     calib.DW = 0;
758     if ( myposiz == -1 ) {
759     myposiz = myfil.find("DW_");
760     calib.DW = 1;
761     };
762     if ( myposiz == -1 ) return(6);
763     stringcopy(calib.basepath,filename,0,myposiz);
764     stringcopy(calib.yodalev,filename,myposiz+13,myposiz+15);
765     //
766     // calibration filename
767     //
768     calcalibfile = filename;
769     //
770     // create the ntuple level2 name depending on framework and starting from the YODA name
771     //
772     if ( !strcmp(franame,"paw") ) paw = 1;
773     //
774     if ( paw ){
775     calib.ispaw = 1.;
776     filety = "rz";
777     } else {
778     calib.ispaw = 0.;
779     filety = "root";
780     };
781     //
782     TString numb("2"); // level
783     TString detc("Calorimeter"); // detector
784     char *file = getLEVname(filename,detc,numb,filety);
785     stringstream truefile;
786     truefile.str("");
787     truefile << file;
788     //
789     if ( outDir == "" ){
790     file4 = filename;
791     file3.str("");
792     file3 << file4 << "/Physics/Level2";
793     } else {
794     file4 = outDir;
795     file3.str("");
796     file3 << file4;
797     }
798     //
799     file2.str("");
800     file2 << file3.str().c_str() << "/";
801     file2 << truefile.str().c_str();
802     printf("\n Filename will be: \n %s \n\n",file2.str().c_str());
803     //
804     // check if Level2 directory exists, if not create it.
805     //
806     Int_t ERR = gSystem->MakeDirectory(file3.str().c_str());
807     //
808     if ( !ERR ) {
809     printf(" LEVEL2 directory doesn't exist! Creating LEVEL2 directory \n\n");
810     } else {
811     //
812     // directory exists, check if file exists (if we are not in the force mode)
813     //
814     if ( !FORCE ) {
815     printf(" Not in FORCE mode, check the existence of LEVEL2 data: \n\n");
816     Int_t nofile = 0;
817     if ( paw ){
818     ifstream mypawfile;
819     mypawfile.open(file2.str().c_str());
820     if ( mypawfile ){
821     nofile = 1;
822     } else {
823     printf("Error in opening file: file %s does not exist \n",file2.str().c_str());
824     };
825     } else {
826     TFile tfile(file2.str().c_str());
827     if ( !tfile.IsZombie() ) nofile = 1;
828     };
829     if ( nofile ){
830     printf(" ERROR: file already exists! Use FORCE = 1 to override \n\n");
831     #if defined (__CINT__)
832     gSystem->ChangeDirectory(path);
833     unloadf77lib(Framework,MySQLload);
834     #endif
835     return(3);
836     } else {
837     printf("\n OK, I will create it!\n\n");
838     };
839     };
840     };
841     //
842     // clear calibration time intervals
843     //
844     for (Int_t s=0; s<4;s++){
845     for (Int_t d = 0; d<51; d++){
846     calib.ttime[s][d] = 0 ;
847     calib.time[s][d] = 0 ;
848     };
849     };
850     //
851     // open and read the rootple containing the ADC to MIP conversion values
852     //
853     stringstream calfile;
854     calfile.str("");
855     const char *pam_calib = pathtocalibration();
856     calfile << pam_calib << "/CaloADC2MIP.root";
857     //
858     printf(" Using calorimeter calibration file: \n %s \n",calfile.str().c_str());
859     chfile = new TFile(calfile.str().c_str(),"READ","Calorimeter CALIBRATION data");
860     CalorimeterCalibration *ccalo = 0;
861     if ( chfile->IsZombie() ){
862     if ( !FORCE ){
863     printf("\n ERROR: no calorimeter calibration file! \n If you want to proceed using 26 as conversion factor for all strips use FORCE = 1. \n\n");
864     return(10);
865     } else {
866     printf("\n WARNING: no calorimeter calibration file! \n WARNING: using 26 as conversion factor for all strips. \n");
867     okcalo = 0;
868     };
869     } else {
870     okcalo = 1;
871     ctree = (TTree*)chfile->Get("CaloADC");
872     ctree->SetBranchAddress("Event", &ccalo);
873     //
874     cnevents = ctree->GetEntries();
875     ctree->GetEntry(cnevents-1);
876     };
877     //
878     if ( okcalo ) {
879     for (Int_t m = 0; m < 2 ; m++ ){
880     for (Int_t k = 0; k < 22; k++ ){
881     for (Int_t l = 0; l < 96; l++ ){
882     if ( (ccalo->fp[1][m][k][l] > 20. && ccalo->fp[1][m][k][l] < 32.) || ccalo->mask[m][k][l] == 1. ) {
883     calib.mip[m][k][l] = ccalo->mip[m][k][l];
884     if ( calib.mip[m][k][l] == 0 ) calib.mip[m][k][l] = 26. ;
885     } else {
886     calib.mip[m][k][l] = 26. ;
887     };
888     calib.calped[m][k][l] = 0. ;
889     evento.dexy[m][k][l] = 0. ;
890     evento.dexyc[m][k][l] = 0. ;
891     };
892     };
893     };
894     } else {
895     for (Int_t m = 0; m < 2 ; m++ ){
896     for (Int_t k = 0; k < 22; k++ ){
897     for (Int_t l = 0; l < 96; l++ ){
898     calib.mip[m][k][l] = 26. ;
899     calib.calped[m][k][l] = 0. ;
900     evento.dexy[m][k][l] = 0. ;
901     evento.dexyc[m][k][l] = 0. ;
902     };
903     };
904     };
905     };
906     chfile->Close();
907     //
908     // Are we able to connect to the MySQL server?
909     //
910     if ( MySQL ){
911     calib.mysql = 1;
912     gSystem->ChangeDirectory(path);
913     #if defined (__CINT__)
914     libload.str("");
915     libload << pam_lib << "/libreadmy";
916     gSystem->Load(libload.str().c_str());
917     libload.str("");
918     libload << pam_lib << "/libreadmy_C";
919     gSystem->Load(libload.str().c_str());
920     #endif
921     calib.db = "romemuons";
922     string myfil = (const char*)filename;
923     Int_t myposiz = myfil.find("dw_");
924     if ( myposiz == -1 ) myposiz = myfil.find("DW_");
925     if ( myposiz == -1 ) return(6);
926     //
927     TString mytmp;
928     Int_t mspos = myposiz+2;
929     Int_t mepos = myposiz+13;
930     stringcopy(mytmp,filename,mspos,mepos);
931     const char *myfile = mytmp;
932     //
933     stringstream tablename;
934     tablename.str("");
935     tablename << "calocalib_dw" << myfile;
936     TString tabula;
937     tabula = tablename.str().c_str();
938     //
939     printf("\n Try the connection to the MySQL database in Trieste...\n");
940     qy.str("");
941     qy << "select * from " << tablename.str().c_str();
942     chkmysql = creadmysql(calib.db,(char *)qy.str().c_str(),mystruct);
943     if ( chkmysql || mystruct.totime == -1 ) {
944     printf("\n- ERROR: problems querying MySQL database! -\n\n- ERROR: Cannot use the MySQL database! -\n\n");
945     if ( mystruct.totime == -1 ) printf("- ERROR: it seems there are no data in table called %s -\n\n",tablename.str().c_str());
946     if ( chkmysql == 1 ) printf("- ERROR: problems in MySQL initialization -\n\n");
947     if ( chkmysql == 2 ) printf("- ERROR: problems in MySQL login, check username and password -\n\n");
948     if ( chkmysql == 3 ) printf("- ERROR: it seems there is no table called %s -\n\n",tablename.str().c_str());
949     // if ( !FORCE ){
950     // printf(" Exiting! Use FORCE = 1 or MySQL = 0 to run in standalone mode \n\n");
951     printf(" Exiting! Check your installation or contact calorimeter developers \n\n");
952     #if defined (__CINT__)
953     gSystem->ChangeDirectory(path);
954     unloadf77lib(Framework,MySQL);
955     #endif
956     return(5);
957     //};
958     // printf("WARNING: running in standalone mode! \n\n");
959     // printf("WARNING: I will search for the best calibration\n\n");
960     // calib.mysql = 0;
961     // MySQL = -1;
962     } else {
963     printf(" ...OK, the connection is fine! \n\n Using database \"%s\", table \"%s\"\n\n",calib.db,tablename.str().c_str());
964     };
965     calib.tablename = tablename.str().c_str();
966     } else {
967     calib.tablename = "";
968     calib.mysql = 0;
969     };
970     if ( MySQL < 1 ){
971     //
972     // first of all find the best calibrations for this file
973     //
974     Int_t wused = 0;
975     TString nfilen;
976     stringcopy(nfilen,file2.str().c_str());
977     CaloFindCalibs(filename, calcalibfile, wused, calib);
978     if ( wused == 1 ) calcalibfile = filename;
979     //
980     // print on the screen the results:
981     //
982     const char *ffile = filename;
983     printf(" ------ %s ------- \n \n",ffile);
984     Int_t calibex = 0;
985     for (Int_t s=0; s<4;s++){
986     printf(" ** SECTION %i **\n",s);
987     for (Int_t d = 0; d<51; d++){
988     TString pfile;
989     if ( calib.ttime[s][d] != 0 ) {
990     calibex++;
991     if ( calib.fcode[s][d] != 10 ){
992     TString file2f = "";
993     stringcopy(file2f,calcalibfile);
994     pfile = whatnamewith(file2f,calib.fcode[s][d]);
995     } else {
996     pfile = (TString)calcalibfile;
997     };
998     const char *ffile = pfile;
999     printf(" - from time %i to time %i use calibration at\n time %i, file: %s \n",calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d],ffile);
1000     if ( !strcmp(ffile,"wrong") ) calibex--;
1001     };
1002     };
1003     printf("\n");
1004     };
1005     printf(" ----------------------------------------------------------------------- \n \n");
1006    
1007     //
1008     // const char *file2 = nfilen; // the full path and name of the level1 ntuple
1009     file2.str("");
1010     file2 << nfilen;
1011     if ( calibex < 4 ) {
1012     printf("No full calibration data in this file, sorry!\n\n ");
1013     //
1014     // remove the empty file!
1015     //
1016     stringstream remfile;
1017     remfile.str("");
1018     remfile << "rm -f " << file2.str().c_str();
1019     gSystem->Exec(remfile.str().c_str());
1020     #if defined (__CINT__)
1021     gSystem->ChangeDirectory(path);
1022     unloadf77lib(Framework,MySQL);
1023     #endif
1024     return(4);
1025     };
1026    
1027     //
1028     // fill with the first calibration values the common containing calibration informations
1029     //
1030     for (Int_t s = 0; s < 4; s++){
1031     b[s]=0;
1032     TString pfile;
1033     if ( calib.fcode[s][b[s]] != 10 ){
1034     TString file2f = "";
1035     stringcopy(file2f,calcalibfile);
1036     pfile = whatnamewith(file2f,calib.fcode[s][b[s]]);
1037     } else {
1038     pfile = (TString)calcalibfile;
1039     };
1040     CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
1041     };
1042     } else {
1043     calib.obtjump = 0;
1044     calib.time[0][0] = -1;
1045     calib.time[1][0] = -1;
1046     calib.time[2][0] = -1;
1047     calib.time[3][0] = -1;
1048     };
1049    
1050     //
1051     // Open file to write depending on the framework
1052     //
1053     //
1054     // In any case we must create a tree
1055     //
1056     TTree *tree = 0;
1057     //
1058     // Open tracker ntuple or rootples and load magnetic field maps.
1059     //
1060     //
1061     printf(" TRACKER: loading the magnetic field maps...\n\n\n");
1062     pam_calib = pathtocalibration();
1063     stringstream bdir;
1064     bdir.str("");
1065     bdir << pam_calib << ".";
1066     creadB(bdir.str().c_str());
1067     //
1068     printf(" ...done! \n");
1069     //
1070     printf("\n Check the existence of tracker data... \n");
1071     Int_t isrootple = 1;
1072     struct Tracklev2 trk;
1073     TFile *trfile = emigetFile(trackerdir,"Physics.Level2","Tracker");
1074     TTree *tr = 0;
1075     Int_t trkpaw = 0;
1076     Int_t trnevents = 0;
1077     stringstream filetrk;
1078     if ( !trfile ) {
1079     const char *file0 = filename;
1080     string fil = (const char *)filename;
1081     Int_t posiz = fil.find("dw_");
1082     if ( posiz == -1 ) posiz = fil.find("DW_");
1083     if ( posiz == -1 ) return(5);
1084     //
1085     TString trntmp;
1086     Int_t spos = posiz+3;
1087     Int_t epos = posiz+13;
1088     stringcopy(trntmp,filename,spos,epos);
1089     const char *trname = trntmp;
1090     if ( !trkdir ){
1091     filetrk.str("");
1092     filetrk << file0 << "/Physics/Level2/DW_";
1093     filetrk << trname << "_level2.rz";
1094     } else {
1095     const char *trackdir = trackerdir;
1096     filetrk.str("");
1097     filetrk << trackdir << "/DW_";
1098     filetrk << trname << "_level2.rz";
1099     };
1100     //
1101     ifstream mypawfile;
1102     mypawfile.open(filetrk.str().c_str());
1103     if ( mypawfile ){
1104     trkpaw = 1;
1105     isrootple = 0;
1106     printf(" ...found tracker level2 NTUPLE:\n %s \n\n",filetrk.str().c_str());
1107     #if defined (__CINT__)
1108     const char *sdir=gSystem->Getenv("PAM_LIB");
1109     bdir.str("");
1110     bdir << sdir << "/liboptrklev2.so";
1111     gSystem->Load(bdir.str().c_str());
1112     bdir.str("");
1113     bdir << sdir << "/liboptrklev2_C.so";
1114     gSystem->Load(bdir.str().c_str());
1115     bdir.str("");
1116     bdir << sdir << "/libretrklev2.so";
1117     gSystem->Load(bdir.str().c_str());
1118     bdir.str("");
1119     bdir << sdir << "/libretrklev2_C.so";
1120     gSystem->Load(bdir.str().c_str());
1121     bdir.str("");
1122     bdir << sdir << "/libcltrklev2.so";
1123     gSystem->Load(bdir.str().c_str());
1124     bdir.str("");
1125     bdir << sdir << "/libcltrklev2_C.so";
1126     gSystem->Load(bdir.str().c_str());
1127     #endif
1128     trnevents = 0;
1129     coptrklev2((char *)filetrk.str().c_str(),trk,trnevents);
1130     } else {
1131     printf("Error in opening file: file %s does not exist \n",file2.str().c_str());
1132     printf("\n No tracker data! You must first process data with tracker programs\n and then you can run CaloLEVEL2. \n Exiting... \n\n");
1133     return(8);
1134     };
1135     } else {
1136     //
1137     const char *file0 = filename;
1138     string fil = (const char *)filename;
1139     Int_t posiz = fil.find("dw_");
1140     if ( posiz == -1 ) posiz = fil.find("DW_");
1141     if ( posiz == -1 ) return(5);
1142     //
1143     TString trntmp;
1144     Int_t spos = posiz;
1145     Int_t epos = posiz+15;
1146     stringcopy(trntmp,filename,spos,epos);
1147     const char *trname = trntmp;
1148     //
1149     if ( !trkdir ){
1150     filetrk.str("");
1151     filetrk << file0 << "/Physics/Level2/";
1152     filetrk << trname << ".Physics.Level2.Tracker.Event.root";
1153     } else {
1154     const char *trkdirname = trackerdir;
1155     filetrk.str("");
1156     filetrk << trkdirname << "/";
1157     filetrk << trname << ".Physics.Level2.Tracker.Event.root";
1158     };
1159     //
1160     printf(" ...found tracker level2 ROOTPLE:\n %s \n",filetrk.str().c_str());
1161     tr = (TTree*) trfile->Get("TrkLevel2");
1162     settrklev2(tr,trk);
1163     trnevents = tr->GetEntries();
1164     };
1165     if ( trnevents != nevents ){
1166     printf("\n WARNING: different length for tracker and calorimeter ntuples!\n");
1167     printf(" calorimeter: %i events \n",(int)nevents);
1168     printf(" tracker : %i events \n\n",(int)trnevents);
1169     };
1170     //
1171     stringstream name;
1172     TString tmpfile2;
1173     stringcopy(tmpfile2,file2.str().c_str());
1174     const char *tmpname = tmpfile2;
1175     TFile *hfile = 0;
1176     struct CaLevel2 clevel2;
1177     if ( paw ){
1178     name.str("");
1179     name << tmpname;
1180     copencalol2((char *)name.str().c_str(),trkpaw,swcode,swtrkcode);
1181     } else {
1182     char *type;
1183     type = "NEW";
1184     if ( FORCE ) type = "RECREATE";
1185     hfile = new TFile(file2.str().c_str(),type,"Calorimeter LEVEL2 data");
1186     //
1187     // hfile = new TFile(file2,type,"Calorimeter LEVEL2 data",0);
1188     //
1189     tree = new TTree("CaloLevel2","PAMELA Level2 calorimeter data");
1190     tree->Branch("OBT",&clevel2.OBT,"OBT/I");
1191     tree->Branch("pkt_num",&clevel2.pkt_num,"pkt_num/I");
1192     tree->Branch("pro_num",&clevel2.pro_num,"pro_num/I");
1193     tree->Branch("trigty",&clevel2.trigty,"trigty/F");
1194     tree->Branch("good",&clevel2.good,"good/I");
1195     tree->Branch("perr",clevel2.perr,"perr[4]/I");
1196     tree->Branch("swerr",clevel2.swerr,"swerr[4]/I");
1197     tree->Branch("crc",clevel2.crc,"crc[4]/I");
1198     tree->Branch("nstrip",&clevel2.nstrip,"nstrip/F");
1199     tree->Branch("qtot",&clevel2.qtot,"qtot/F");
1200     tree->Branch("ncore",&clevel2.ncore,"ncore/F");
1201     tree->Branch("qcore",&clevel2.qcore,"qcore/F");
1202     tree->Branch("impx",&clevel2.impx,"impx/F");
1203     tree->Branch("impy",&clevel2.impy,"impy/F");
1204     tree->Branch("tanx",&clevel2.tanx,"tanx/F");
1205     tree->Branch("tany",&clevel2.tany,"tany/F");
1206     tree->Branch("nint",&clevel2.nint,"nint/F");
1207     tree->Branch("ncyl",&clevel2.ncyl,"ncyl/F");
1208     tree->Branch("qcyl",&clevel2.qcyl,"qcyl/F");
1209     tree->Branch("qtrack",&clevel2.qtrack,"qtrack/F");
1210     tree->Branch("qmax",&clevel2.qmax,"qmax/F");
1211     tree->Branch("nx22",&clevel2.nx22,"nx22/F");
1212     tree->Branch("qx22",&clevel2.nx22,"qx22/F");
1213     tree->Branch("qq",clevel2.qq,"qq[4]/F");
1214     tree->Branch("qtrackx",&clevel2.qtrackx,"qtrackx/F");
1215     tree->Branch("qtracky",&clevel2.qtracky,"qtracky/F");
1216     tree->Branch("dxtrack",&clevel2.dxtrack,"dxtrack/F");
1217     tree->Branch("dytrack",&clevel2.dytrack,"dytrack/F");
1218     tree->Branch("qlast",&clevel2.qlast,"qlast/F");
1219     tree->Branch("nlast",&clevel2.nlast,"nlast/F");
1220     tree->Branch("qpre",&clevel2.qpre,"qpre/F");
1221     tree->Branch("npre",&clevel2.npre,"npre/F");
1222     tree->Branch("qpresh",&clevel2.qpresh,"qpresh/F");
1223     tree->Branch("npresh",&clevel2.npresh,"npresh/F");
1224     tree->Branch("qlow",&clevel2.qlow,"qlow/F");
1225     tree->Branch("nlow",&clevel2.nlow,"nlow/F");
1226     tree->Branch("qtr",&clevel2.qtr,"qtr/F");
1227     tree->Branch("ntr",&clevel2.ntr,"ntr/F");
1228     tree->Branch("planetot",&clevel2.planetot,"planetot/F");
1229     tree->Branch("qmean",&clevel2.qmean,"qmean/F");
1230     tree->Branch("varcfit",clevel2.varcfit,"varcfit[2]/F");
1231     tree->Branch("npcfit",clevel2.npcfit,"npcfit[2]/I");
1232     tree->Branch("thex",&clevel2.thex,"thex/F");
1233     tree->Branch("they",&clevel2.they,"they/F");
1234     tree->Branch("cibar",clevel2.cibar,"cibar[22][2]/I");
1235     tree->Branch("tibar",clevel2.tibar,"tibar[22][2]/I");
1236     tree->Branch("cbar",clevel2.cbar,"cbar[22][2]/F");
1237     tree->Branch("tbar",clevel2.tbar,"tbar[22][2]/F");
1238     //
1239     TTree *software = 0;
1240     software = new TTree("Software","Software used to generate data");
1241     software->Branch("swcode",&swcode,"swcode/I");
1242     software->Branch("swtrkcode",&swtrkcode,"swtrkcode/I");
1243     software->Fill();
1244     };
1245     //
1246     // run over all the events
1247     //
1248     printf("\n Processed events: \n\n");
1249     //
1250     pamela::PscuHeader *ph = 0;
1251     pamela::EventHeader *eh = 0;
1252     Int_t caloerr = 0;
1253     Int_t i = -1;
1254     Int_t itr;
1255     if ( isrootple ){
1256     itr = -1;
1257     } else {
1258     itr = 0;
1259     };
1260     Int_t retval = 0;
1261     Int_t syncro = 1;
1262     Int_t trklost = 0;
1263     Int_t calolost = 0;
1264     Int_t pktnum = 0;
1265     Int_t obt = 0;
1266     //
1267     while ( i < (nevents-1) ){
1268     //
1269     if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);
1270     // printf(" %i \n",i);
1271     //
1272     // look for tracker data
1273     //
1274     itr++;
1275     i++;
1276     syncro = 1;
1277     #if !defined (__CINT__)
1278     trkcalosync:
1279     #endif
1280     #if defined (__CINT__)
1281     trkcalosync: printf("");
1282     #endif
1283     //
1284     // check if we have tracker data
1285     //
1286     if ( i >= nevents ){
1287     printf(" WARNING: no more calorimeter data.\n");
1288     if ( nevents > trnevents ) calolost += (snevents-trnevents);
1289     if ( nevents < trnevents ) trklost += (-snevents+trnevents);
1290     retval = 9;
1291     goto closeandexit;
1292     };
1293     if ( itr > trnevents && isrootple ){
1294     printf(" WARNING: no more tracker data.\n");
1295     if ( nevents > trnevents ) calolost += (snevents-trnevents);
1296     if ( nevents < trnevents ) trklost += (-snevents+trnevents);
1297     retval = 9;
1298     goto closeandexit;
1299     };
1300     if ( itr >= trnevents+1 && !isrootple ){
1301     printf(" WARNING: no more tracker data.\n");
1302     if ( nevents > trnevents ) calolost += (snevents-trnevents);
1303     if ( nevents < trnevents ) trklost += (-snevents+trnevents);
1304     retval = 9;
1305     goto closeandexit;
1306     };
1307     //
1308     // retrieve tracker informations
1309     //
1310     if ( isrootple ){
1311     if ( debug ) printf(" itr %i i %i nevents %i trnevents %i \n",itr,i,(int)nevents,trnevents);
1312     tr->GetEntry(itr);
1313     } else {
1314     if ( debug ) printf(" itr %i i %i nevents %i trnevents %i \n",itr,i,(int)nevents,trnevents);
1315     cretrklev2(itr,trk);
1316     };
1317     //
1318     // check synchronization tracker and calorimeter informations:
1319     //
1320     otr->SetBranchAddress("Header", &eh);
1321     otr->GetEntry(i);
1322     ph = eh->GetPscuHeader();
1323     pktnum = ph->GetCounter();
1324     obt = ph->GetOrbitalTime();
1325     if ( pktnum != trk.pkt_num || obt != trk.obt ){
1326     if ( pktnum > trk.pkt_num || obt > trk.obt ){
1327     if ( debug ) printf("itr++ %i pktnum calo %i trk %i obt calo %i trk %i \n",itr,pktnum,trk.pkt_num,obt,trk.obt);
1328     itr++;
1329     trklost++;
1330     if ( syncro ) printf(" WARNING: lost sync! try to recover... \n");
1331     syncro = 0;
1332     goto trkcalosync;
1333     };
1334     if ( pktnum < trk.pkt_num || obt < trk.obt ){
1335     if ( debug ) printf("i++ %i pktnum calo %i trk %i obt calo %i trk %i \n",i,pktnum,trk.pkt_num,obt,trk.obt);
1336     i++;
1337     calolost++;
1338     if ( syncro ) printf(" WARNING: lost sync! try to recover... \n");
1339     syncro = 0;
1340     goto trkcalosync;
1341     };
1342     };
1343     //
1344     // here we have synchronized data
1345     //
1346     if ( !syncro ) {
1347     printf(" ...synchronization recovered! \n");
1348     printf(" Sync info: \n - tracker packets without calorimeter %i\n - calorimeter packets without tracker %i\n\n",trklost,calolost);
1349     syncro = 1;
1350     if ( debug ) printf("pktnum calo %i trk %i obt calo %i trk %i \n",pktnum,trk.pkt_num,obt,trk.obt);
1351     };
1352     //
1353     // store track information in the calib structure
1354     //
1355     if ( trk.ntrk > 0 ){
1356     imtrack = 0;
1357     for (Int_t e = 0; e<trk.ntrk ; e++){
1358     if ( trk.image[e] != 0 ) imtrack++;
1359     };
1360     if ( (trk.ntrk-(imtrack/2)) == 1 ){
1361     calib.good2 = trk.good2;
1362     for (Int_t e = 0; e<5 ; e++){
1363     calib.al_p[e][0] = trk.al[0][e];
1364     if ( imtrack != 0 ) calib.al_p[e][1] = trk.al[1][e];
1365     };
1366     if ( imtrack != 0 ) {
1367     if ( trk.chi2[0]>trk.chi2[1] ) {
1368     calib.trkchi2 = 1;
1369     } else {
1370     calib.trkchi2 = 2;
1371     };
1372     } else {
1373     calib.trkchi2 = 1;
1374     };
1375     } else {
1376     calib.good2 = false;
1377     for (Int_t e = 0; e<5 ; e++){
1378     calib.al_p[e][0] = 0.;
1379     calib.al_p[e][1] = 0.;
1380     };
1381     };
1382     } else {
1383     calib.good2 = false;
1384     for (Int_t e = 0; e<5 ; e++){
1385     calib.al_p[e][0] = 0.;
1386     calib.al_p[e][1] = 0.;
1387     };
1388     };
1389     //
1390     // calibrate calorimeter data and retrieve level2 informations
1391     //
1392     // ====> NOTICE: in the case of no tracks or not good tracker events the program will fill the ntuple with zero values; <====
1393     // ====> in case of multiple tracks the program will calculate variables using the state vector and rigidity of the first track stored <====
1394     //
1395     caloerr = calolevel2core(i,b,otr,tree,clevel2,evento,calib,calcalibfile);
1396     //
1397     //
1398     //
1399     if ( caloerr ){
1400     if ( i > 0 ){
1401     printf("\nTrying to close the file anyway...\n");
1402     if ( paw ){
1403     name.str("");
1404     name << tmpname;
1405     cclosecalol2((char *)name.str().c_str());
1406     } else {
1407     hfile->Write();
1408     hfile->Close();
1409     };
1410     printf("...done!\n");
1411     } else {
1412     printf("\nERROR: output file is empty! \n");
1413     stringstream remfile;
1414     remfile.str("");
1415     remfile << "rm -f " << file2.str().c_str();
1416     gSystem->Exec(remfile.str().c_str());
1417     };
1418     printf("\nERROR: processed %i out of %i events\n",i,(int)nevents);
1419     printf("\nERROR: an error occurred during processing!\n\n Exiting...\n\n");
1420     goto theend;
1421     };
1422     Int_t ciccio;
1423     ciccio = 1;
1424     };
1425     if ( (i+1) < trnevents ){
1426     printf(" WARNING: no more calorimeter data. \n");
1427     if ( nevents > trnevents ) calolost += (snevents-trnevents);
1428     if ( nevents < trnevents ) trklost += (-snevents+trnevents);
1429     printf(" Sync info: \n - tracker packets without calorimeter %i\n - calorimeter packets without tracker %i\n",trklost,calolost);
1430     retval = 9;
1431     };
1432     //
1433     // close rootple/ntuple files
1434     //
1435     closeandexit:
1436     if ( paw ){
1437     name.str("");
1438     name << tmpname;
1439     cclosecalol2((char *)name.str().c_str());
1440     } else {
1441     hfile->Write();
1442     hfile->Close();
1443     };
1444     if ( isrootple ){
1445     trfile->Close();
1446     } else {
1447     ccltrklev2(trk);
1448     };
1449     printf("\n");
1450     printf(" Finished, exiting... \n");
1451     printf("\n");
1452     if ( trklost || calolost ){
1453     printf("\n Sync info over the whole file: \n - tracker packets without calorimeter %i\n - calorimeter packets without tracker %i\n",trklost,calolost);
1454     };
1455     printf("\n");
1456     theend:
1457     //
1458     // go back to the starting path and unload fortran libraries
1459     //
1460     gSystem->ChangeDirectory(path);
1461     #if defined (__CINT__)
1462     if ( !isrootple ){
1463     const char *sdir=gSystem->Getenv("PAM_LIB");
1464     bdir.str("");
1465     bdir << sdir << "/libcltrklev2_C.so";
1466     gSystem->Unload(bdir.str().c_str());
1467     bdir.str("");
1468     bdir << sdir << "/libcltrklev2.so";
1469     gSystem->Unload(bdir.str().c_str());
1470     bdir.str("");
1471     bdir << sdir << "/libretrklev2_C.so";
1472     gSystem->Unload(bdir.str().c_str());
1473     bdir.str("");
1474     bdir << sdir << "/libretrklev2.so";
1475     gSystem->Unload(bdir.str().c_str());
1476     bdir.str("");
1477     bdir << sdir << "/liboptrklev2_C.so";
1478     gSystem->Unload(bdir.str().c_str());
1479     bdir.str("");
1480     bdir << sdir << "/liboptrklev2.so";
1481     gSystem->Unload(bdir.str().c_str());
1482     };
1483     // unloadf77lib(Framework,MySQL);
1484     #endif
1485     return(retval);
1486     }
1487    

  ViewVC Help
Powered by ViewVC 1.1.23