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

  ViewVC Help
Powered by ViewVC 1.1.23