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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show 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 //
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.04 (2006-03-20)
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 // 4.03 - 4.04 (2006-03-20): Fixed erroneous distance between x/y planes in the core routine
11 //
12 // 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 // 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 // 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 const char *franame = Framework;
146 //
147 char *pam_lib=gSystem->Getenv("PAM_LIB");
148 if ( MySQL == 1 || MySQL == -1 ){
149 stringstream libload;
150 libload.str("");
151 libload << pam_lib << "/libreadmy";
152 gSystem->Unload(libload.str().c_str());
153 libload.str("");
154 libload << pam_lib << "/libreadmy_C";
155 gSystem->Unload(libload.str().c_str());
156 };
157 if ( !strcmp(franame,"paw") ) {
158 stringstream libload;
159 libload.str("");
160 libload << pam_lib << "/libclcalol2_C";
161 gSystem->Unload(libload.str().c_str());
162 libload.str("");
163 libload << pam_lib << "/libclcalol2";
164 gSystem->Unload(libload.str().c_str());
165 libload.str("");
166 libload << pam_lib << "/libopcalol2_C";
167 gSystem->Unload(libload.str().c_str());
168 libload.str("");
169 libload << pam_lib << "/libopcalol2";
170 gSystem->Unload(libload.str().c_str());
171 };
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 }
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 //
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 //
345 // for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
346 //
347 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 };
363 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 };
422 //
423 // run over preamplifiers
424 //
425 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 };
441 };
442 };
443 //
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 //
462 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 //
467 // no suitable new baseline, use old ones and compress data!
468 //
469 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 };
477 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 };
484 // 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 };
518 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 //
528 // 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 //
530 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 };
571 };
572 };
573 //
574 } else {
575 for (Int_t nn = 0; nn < 96; nn++){
576 clevel1.estrip[nn][m][l] = 0.;
577 };
578 };
579 };
580 };
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 }
624 //
625 short int CaloLEVEL2(TString filename, TString TrackerDir="", TString outDir ="", TString Framework = "root", Int_t FORCE = 0){
626 //
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 #if defined(__CINT__)
635 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 #endif
694 const char *franame = Framework;
695 #if defined (__CINT__)
696 if ( !strcmp(franame,"paw") ) {
697 libload.str("");
698 libload << pam_lib << "/libopcalol2";
699 gSystem->Load(libload.str().c_str());
700 libload.str("");
701 libload << pam_lib << "/libopcalol2_C";
702 gSystem->Load(libload.str().c_str());
703 libload.str("");
704 libload << pam_lib << "/libclcalol2";
705 gSystem->Load(libload.str().c_str());
706 libload.str("");
707 libload << pam_lib << "/libclcalol2_C";
708 gSystem->Load(libload.str().c_str());
709 };
710 #endif
711 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 #endif
724 //
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 #endif
734 return(1);
735 };
736 caloFile = emigetFile(filename, "Calorimeter");
737 if ( !caloFile ){
738 printf("\n\n ERROR: no calorimeter file! \n\n");
739 #if defined (__CINT__)
740 gSystem->ChangeDirectory(path);
741 unloadf77lib(Framework,MySQLload);
742 #endif
743 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 #if defined (__CINT__)
752 gSystem->ChangeDirectory(path);
753 unloadf77lib(Framework,MySQLload);
754 #endif
755 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 #if defined (__CINT__)
771 gSystem->ChangeDirectory(path);
772 unloadf77lib(Framework,MySQLload);
773 #endif
774 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 #if defined (__CINT__)
856 gSystem->ChangeDirectory(path);
857 unloadf77lib(Framework,MySQLload);
858 #endif
859 return(3);
860 } else {
861 printf("\n OK, I will create it!\n\n");
862 };
863 };
864 };
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 } else {
890 printf("\n WARNING: no calorimeter calibration file! \n WARNING: using 26 as conversion factor for all strips. \n");
891 okcalo = 0;
892 };
893 } 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 };
916 };
917 };
918 } 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 };
927 };
928 };
929 };
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 //
951 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 //
963 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 } else {
987 printf(" ...OK, the connection is fine! \n\n Using database \"%s\", table \"%s\"\n\n",calib.db,tablename.str().c_str());
988 };
989 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 };
1026 };
1027 printf("\n");
1028 };
1029 printf(" ----------------------------------------------------------------------- \n \n");
1030
1031 //
1032 // 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 #if defined (__CINT__)
1045 gSystem->ChangeDirectory(path);
1046 unloadf77lib(Framework,MySQL);
1047 #endif
1048 return(4);
1049 };
1050
1051 //
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 } else {
1119 const char *trackdir = trackerdir;
1120 filetrk.str("");
1121 filetrk << trackdir << "/DW_";
1122 filetrk << trname << "_level2.rz";
1123 };
1124 //
1125 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 #if defined (__CINT__)
1132 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 #endif
1152 trnevents = 0;
1153 coptrklev2((char *)filetrk.str().c_str(),trk,trnevents);
1154 } else {
1155 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 } else {
1178 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 //
1313 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 //
1335 // retrieve tracker informations
1336 //
1337 if ( isrootple ){
1338 if ( debug ) printf(" itr %i i %i nevents %i trnevents %i \n",itr,i,(int)nevents,trnevents);
1339 tr->GetEntry(itr);
1340 } else {
1341 if ( debug ) printf(" itr %i i %i nevents %i trnevents %i \n",itr,i,(int)nevents,trnevents);
1342 cretrklev2(itr,trk);
1343 };
1344 //
1345 // 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 itr++;
1356 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 i++;
1364 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 };
1393 if ( imtrack != 0 ) {
1394 if ( trk.chi2[0]>trk.chi2[1] ) {
1395 calib.trkchi2 = 1;
1396 } else {
1397 calib.trkchi2 = 2;
1398 };
1399 } else {
1400 calib.trkchi2 = 1;
1401 };
1402 } 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 };
1416 //
1417 // 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 //
1426 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 closeandexit:
1463 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 theend:
1484 //
1485 // go back to the starting path and unload fortran libraries
1486 //
1487 gSystem->ChangeDirectory(path);
1488 #if defined (__CINT__)
1489 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 #endif
1512 return(retval);
1513 }
1514

  ViewVC Help
Powered by ViewVC 1.1.23