/[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.4 - (hide annotations) (download)
Mon Mar 20 10:12:13 2006 UTC (18 years, 11 months ago) by mocchiut
Branch: MAIN
CVS Tags: v4r04, HEAD
Changes since 1.3: +3 -1 lines
Error in x/y plane distance in the core routine fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23