/[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.1.1.1 - (show annotations) (download) (vendor branch)
Mon Dec 5 16:13:53 2005 UTC (19 years, 2 months ago) by mocchiut
Branch: LEVEL2
CVS Tags: v4r00, start
Changes since 1.1: +0 -0 lines
Imported sources

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

  ViewVC Help
Powered by ViewVC 1.1.23