/[PAMELA software]/calo/ground/COMMON/inc/CaloFunctions.h
ViewVC logotype

Contents of /calo/ground/COMMON/inc/CaloFunctions.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Dec 6 10:18:53 2005 UTC (19 years, 2 months ago) by mocchiut
Branch: COMMON
CVS Tags: start, v2r00
Changes since 1.1: +0 -0 lines
File MIME type: text/plain
Imported sources

1 //
2 // Calorimeter useful functions and subroutines - Emiliano Mocchiutti
3 //
4 // CaloFunctions.h version 4.00 (2005-11-29)
5 //
6 // Programs in this file are called by other programs and cannot be run by hand.
7 //
8 // Changelog:
9 //
10 // 3.10 - 4.00 (2005-11-29): preparing for the final release, changed checkemilib (PAM_YODALIB -> PAM_YODA).
11 //
12 // 3.09 - 3.10 (2005-11-15): changed langaupro since it did not work in the compiled version; changed WhatToDo to work in compiled version.
13 //
14 // 3.08 - 3.09 (2005-11-11): small changes to PrintFigure subroutine
15 //
16 // 3.07 - 3.08 (2005-10-13): classes are now in another file (caloclasses.h)
17 //
18 // 3.06 - 3.07 (2005-10-10): added subroutine ColorTOFMIP for TOF colors.
19 //
20 // 3.05 - 3.06 (2005-10-04): added subroutine getfileLEVname.
21 //
22 // 3.04 - 3.05 (2005-08-04): Fixed bugs when running on a 64 bit arch.
23 //
24 // 3.03 - 3.04 (2005-07-14): Bug in getLEVname, fixed.
25 //
26 // 3.02 - 3.03 (2005-07-13): Finally CaloFunctions.h can be compiled under ROOT! Added clones of routines in the yodautility file no needed anymore.
27 //
28 // 3.01 - 3.02 (2005-07-11): Some more tuning for the macro compilation.
29 //
30 // 3.00 - 3.01 (2005-07-06): changed getLEVname since it gave some problems when called more than once in a ROOT sesssion. Tracker structurs now placed in a tracker
31 // include file.
32 //
33 // 2.27 - 3.00 (2005-07-05): Started the tuning for the compilation of the macros. Changed order of functions, added include files for ROOT routines, fixed some
34 // implicit casting. Fixed redefinition of class constructors. Small change in the calib structure.
35 //
36 // 2.26 - 2.27 (2005-06-23): Changed class CalorimeterLevel2 and structures to add good, perr,swerr and crc infos.
37 //
38 // 2.25 - 2.26 (2005-06-16): Changed class CalorimeterLevel2 and structures to convert OBT, headc and evfile to integers.
39 //
40 // 2.24 - 2.25 (2005-06-07): Added subroutine checkifpretyodaversion to check if file is older than 050515_007 (since that file there is a new format in TMTC data).
41 //
42 // 2.23 - 2.24 (2005-06-06): Small changes in the struct "calib". Change type of CaloPede subroutine, now it is an integer and not void: return 3 if it doesn't find the
43 // calibration file.
44 //
45 // 2.22 - 2.23 (2005-05-30): Added checkmerging and fillmerging subroutines to avoid processing twice the same event.
46 //
47 // 2.21 - 2.22 (2005-05-19): Speed-up and better behaviour of CaloFindcalibs. No need to change calls to this subroutine, all internal changes.
48 //
49 // 2.20 - 2.21 (2005-05-17): Changed CaloFindCalibs in order to check also the calibration file provided if any but first check the current file.
50 //
51 // 2.19 - 2.20 (2005-05-06): Added mip to structure calib.
52 //
53 // 2.18 - 2.19 (2005-05-05): Added calorimeter class CalorimeterLevel2 and structures CaLevel1 CaLevel2. Changed getLEVname subroutines (still compatible with older version).
54 //
55 // 2.17 - 2.18 (2005-04-19): Added the old subroutine Calo1stcalib with the name Calo1stcalibOLD (to be used in ShowEvent.c).
56 //
57 // 2.16 - 2.17 (2005-04-14): Added classes for calorimeter ADC to MIP conversion. Bug fixed in CalorimeterCalibration class. CaloFindBaseRawNC subroutine added to
58 // calculate baseline from raw data without any compression. Added getEmiFile subroutine to find YODA files (equivalent of getFile in
59 // yodautility.c); sintax is similar to getFile but it doesn't crash after 1111 or 2222 files.
60 //
61 // 2.15 - 2.16 (2005-04-11): Bugs fixes in the CaloFindCalibs, fetchpreviousfile and whatnamewith subroutines. Consider the case of an OBT jump in data.
62 //
63 // 2.14 - 2.15 (2005-04-07): Changes in the CaloFindCalibs subroutine. Added CaloCompressData subroutine, some small bugs fixed.
64 //
65 // 2.13 - 2.14 (2005-04-06): Changes in the CaloFindCalibs subroutine. Added fetchpreviousfile and whatnamewith subroutines. Now
66 // the program looks for calorimeter calibrations also in previous files and associates them to data in
67 // a correct way.
68 //
69 // 2.12 - 2.13 (2005-03-17): Changes in the WhatToDo subroutine.
70 //
71 // 2.11 - 2.12 (2005-03-09): PrintFigure subroutine changed to accomplish the new "saveas" script style.
72 //
73 // 2.10 - 2.11 (2005-03-04): Changed structure CalorimeterCalibration to save also fit parameters.
74 //
75 // 2.09 - 2.10 (2005-03-02): Show calibration CRC errors (in CaloFindCalibs). Added langaupro subroutines.
76 //
77 // 2.08 - 2.09 (2005-02-28): Calorimeter calibration classes changed (required by the new CaloADC2MIP).
78 //
79 // 2.07 - 2.08 (2005-02-25): check the existence of calibration header file in CaloFindCalibs!
80 //
81 // 2.06 - 2.07 (2005-02-24): bug in CaloFindBase! DO NOT put to zero the ADC values if you can't determine the baseline!
82 //
83 // 2.05 - 2.06 (2005-02-24): CalorimeterLevel1 class must be cleared during initialization to avoid the so-called "gundam-bug".
84 //
85 // 2.04 - 2.05 (2005-02-24): Changed variable definition from C/C++ style to ROOT style (int->Int_t).
86 //
87 // 2.03 - 2.04 (2005-02-23): Mispelling in the WhatToDo subroutine, fixed.
88 //
89 // 2.02 - 2.03 (2005-02-18): Calo1stCalib call CaloPede also if we don't have all four the calorimeter calibrations.
90 //
91 // 2.01 - 2.02 (2005-02-10): "progressive number" instead of "event number" in WhatToDo subroutine.
92 //
93 // 2.00 - 2.01 (2005-02-07): Added langau fit subroutines and the PrintFigure subroutine (used in CaloPLANES mainly). Some bug fixed,
94 // added user interaction routine and getfilename and getlevname subroutines to retrieve file name and to give
95 // the level name correctly.
96 //
97 // 1.00 - 2.00 (2005-01-17): Cleanup of the code
98 //
99 #include <iostream>
100 #include <TSystem.h>
101 #include <TH1.h>
102 #include <TF1.h>
103 #include <TROOT.h>
104 #include <TStyle.h>
105 #include <TClassEdit.h>
106 #include <TObject.h>
107 #include <TList.h>
108 #include <TSystem.h>
109 #include <TSystemDirectory.h>
110 #include <TString.h>
111 #include <TFile.h>
112 #include <TClass.h>
113 #include <TCanvas.h>
114 #include <TKey.h>
115 #include <TClassTable.h>
116
117 typedef struct Chmerge {
118 Int_t fromflag;
119 Int_t scanflag;
120 Int_t pscuOBT;
121 Int_t pscucount;
122 Int_t lastOBT[1000];
123 Int_t lastcount[1000];
124 } chmerge;
125
126 typedef struct Evento {
127 Int_t iev;
128 Int_t stwerr[4];
129 Float_t emin;
130 Float_t perror[4];
131 Float_t dexy[2][22][96];
132 Float_t dexyc[2][22][96];
133 Double_t base[2][22][6];
134 Float_t calselftrig[4][7];
135 Float_t calIItrig[4];
136 Float_t calstriphit[4];
137 Float_t calDSPtaberr[4];
138 Float_t calevnum[4];
139 Float_t shift;
140 } evento;
141
142 typedef struct Calib {
143 TString tablename;
144 TString basepath;
145 TString yodalev;
146 char *db;
147 Int_t DW;
148 Int_t mysql;
149 Int_t obtjump;
150 Int_t time[4][51];
151 Int_t ttime[4][51];
152 Int_t fcode[4][51];
153 Int_t iev;
154 Int_t cstwerr[4];
155 Float_t ispaw;
156 Float_t cperror[4];
157 Float_t mip[2][22][96];
158 Float_t calped[2][22][96];
159 Float_t calgood[2][22][96];
160 Float_t calthr[2][22][6];
161 Float_t calrms[2][22][96];
162 Float_t calbase[2][22][6];
163 Float_t calvar[2][22][6];
164 Float_t calpuls[2][22][96];
165 Double_t sbase[2][22][6];
166 Double_t al_p[5][2];
167 Int_t trkchi2;
168 Bool_t good2;
169 } calib;
170
171 void emicheckLib(){
172 if (!TClassTable::GetDict("pamela::TmtcEvent")){
173 char *pamyodalib = gSystem->ExpandPathName("$PAM_YODA");
174 stringstream libYoda;
175 libYoda.str("");
176 libYoda << pamyodalib << "/lib/libyoda.so";
177 gSystem->Load(libYoda.str().c_str());
178 }
179 }
180
181 void emidigForFiles(TList& out, TSystemDirectory *tsd, string defin){
182 const char *initialDir = gSystem->pwd();
183 TList *lnk = tsd->GetListOfFiles();
184 if (lnk==0) return;
185
186 TSystemFile *file = (TSystemFile*)lnk->First();
187 if (file->IsZombie()) {
188 lnk->Delete();
189 gSystem->cd(initialDir);
190 return;
191 }
192
193 string *fileDes = new string(file->GetName());
194 unsigned int loc = 0;
195 while(file){
196 fileDes = new string(file->GetName());
197 if(!((fileDes->compare(".") == 0) || (fileDes->compare("..") == 0))){
198 if (file->IsDirectory()){
199 emidigForFiles(out, (TSystemDirectory*)file, defin);
200 } else {
201 loc = fileDes->find(defin.c_str(), 0);
202 if (loc != string::npos){
203 out.AddLast((TObject*)file);
204 }
205 }
206
207 }
208 file = (TSystemFile*)lnk->After(file);
209 }
210 gSystem->cd(initialDir);
211 }
212
213 void emimakeAllFriend(TTree* out, TList* input){
214 TList *keys;
215 TKey *key;
216 TTree *tr;
217 string *base;
218 string *className;
219 unsigned int loc;
220 bool cont = false;
221 TSystemFile *tsf = (TSystemFile*)input->First();
222 if (tsf->IsZombie()) return;
223 while(tsf){
224 base = new string(tsf->GetTitle());
225 base->append("/");
226 base->append(tsf->GetName());
227 TFile *tf = new TFile(base->c_str());
228 if (tf->IsZombie()) {
229 tsf = (TSystemFile*)input->After(tsf);
230 continue;
231 }
232 keys = tf->GetListOfKeys();
233 key = (TKey*)keys->First();
234 while(key){
235 className = new string(key->GetClassName());
236 loc = className->find("Tree", 0);
237 if(loc != string::npos){
238 tr = (TTree*)key->ReadObj();
239 if (tr->GetEntries() > 0){
240 if(cont){
241 tr = (TTree*)key->ReadObj();
242 out->AddFriend(tr->GetName(), tf);
243 } else {
244 out = ((TTree*)key->ReadObj());
245 cont = true;
246 }
247 }
248 }
249 key = (TKey*)keys->After(key);
250 }
251
252 tsf = (TSystemFile*)input->After(tsf);
253 }
254 }
255
256 void emimakeAllChained(TChain& out, TList* input){
257 TList *keys;
258 TKey *key;
259 TTree *tr;
260 string *base;
261 string *className;
262 unsigned int loc;
263
264 TSystemFile *tsf = (TSystemFile*)input->First();
265 if (tsf->IsZombie()) return;
266
267 while(tsf){
268 base = new string(tsf->GetTitle());
269 base->append("/");
270 base->append(tsf->GetName());
271 TFile *tf = new TFile(base->c_str());
272 if (tf->IsZombie()) {
273 tsf = (TSystemFile*)input->After(tsf);
274 continue;
275 }
276
277 keys = tf->GetListOfKeys();
278 key = (TKey*)keys->First();
279 while(key){
280 className = new string(key->GetClassName());
281 loc = className->find("Tree", 0);
282 if(loc != string::npos){
283 tr = (TTree*)key->ReadObj();
284 if (tr->GetEntries() > 0){
285 out.Add(base->c_str());
286 }
287 }
288 key = (TKey*)keys->After(key);
289 }
290 tsf = (TSystemFile*)input->After(tsf);
291 }
292 }
293
294 TFile* emigetFile(TString base, TString packetType, TString subType = "Event"){
295 TSystemDirectory *targetDir = new TSystemDirectory("", base);
296 TList *filesList = new TList;
297 TSystemFile *tsf = 0;
298 string *tmpString = 0;
299 tmpString = new string(".");
300 tmpString->append(packetType);
301 tmpString->append(".");
302 tmpString->append(subType);
303 tmpString->append(".");
304 emidigForFiles(*filesList, targetDir, tmpString->c_str());
305 if ( filesList->IsEmpty() ) {
306 filesList->Delete();
307 return 0;
308 };
309 targetDir = new TSystemDirectory("", base);
310 tsf = (TSystemFile*)filesList->First();
311 tmpString = new string(tsf->GetTitle());
312 tmpString->append("/");
313 tmpString->append(tsf->GetName());
314 targetDir->Delete();
315 filesList->Delete();
316 return new TFile(tmpString->c_str());
317 }
318
319 Int_t emigetLastNotZeroBin(TH1 *histo){
320 Int_t lastBin = 0;
321 Int_t tempBin = 0;
322 Stat_t minVal = 0;
323 Stat_t tempVal = 0;
324 Int_t range = histo->GetNbinsX();
325
326 tempBin = histo->GetMaximumBin();
327 tempVal = histo->GetBinContent(tempBin);
328 minVal = tempVal;
329 lastBin = tempBin;
330
331 for (Int_t i = histo->GetMaximumBin() + 1; i < range; i++){
332 tempVal = histo->GetBinContent(i);
333 if ((tempVal != 0) && (tempVal < minVal)) {
334 minVal = tempVal;
335 lastBin = i;
336 }
337 }
338 return (Int_t)(lastBin*1.10);
339 }
340
341 Int_t emigetFirstNotZeroBin(TH1 *histo){
342 Int_t lastBin = 0;
343 Int_t tempBin = 0;
344 Stat_t minVal = 0;
345 Stat_t tempVal = 0;
346
347 tempBin = histo->GetMaximumBin();
348 tempVal = histo->GetBinContent(tempBin);
349 minVal = tempVal;
350 lastBin = tempBin;
351
352 for (Int_t i = histo->GetMaximumBin() - 1; i > 0; i--){
353 tempVal = histo->GetBinContent(i);
354 if ((tempVal != 0) && (tempVal < minVal)) {
355 minVal = tempVal;
356 lastBin = i;
357 }
358 }
359 return (Int_t)(lastBin*0.90);
360 }
361
362 void stringcopy(TString& s1, const TString& s2, Int_t from=0, Int_t to=0){
363 if ( to == 0 ){
364 Int_t t2length = s2.Length();
365 s1 = "";
366 to = t2length;
367 };
368 for (Int_t i = from; i<to; i++){
369 s1.Append(s2[i],1);
370 };
371 }
372
373 void stringappend(TString& s1, const TString& s2){
374 Int_t t2length = s2.Length();
375 for (Int_t i = 0; i<t2length; i++){
376 s1.Append(s2[i],1);
377 };
378 }
379
380 TFile* getEmiFile(TString base, TString packetType, TString subType="Event"){
381 //
382 TString newfile;
383 stringcopy(newfile,base);
384 //
385 const char *base2 = base;
386 stringstream fil;
387 fil.str("");
388 fil << base2;
389 Int_t posiz = fil.str().find("dw_");
390 Int_t upper = 0;
391 if ( posiz == -1 ) {
392 posiz = fil.str().find("DW_");
393 upper = 1;
394 };
395 if ( posiz == -1 ) return 0;
396 TString fullname;
397 stringcopy(fullname,base,posiz,posiz+15);
398 stringappend(fullname,".");
399 stringappend(fullname,packetType);
400 stringappend(fullname,".");
401 stringappend(fullname,subType);
402 stringappend(fullname,".root");
403 //
404 stringappend(newfile,"/");
405 stringappend(newfile,packetType);
406 stringappend(newfile,"/");
407 stringappend(newfile,fullname);
408 //
409 if ( !gSystem->IsFileInIncludePath(newfile) ){
410 return 0;
411 } else {
412 return new TFile(newfile);
413 };
414 }
415
416 void CaloCompressData(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
417 Float_t minstrip = 100000.;
418 Float_t rms = 0.;
419 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
420 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) {
421 minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
422 rms = calib.calthr[l][m][pre];
423 };
424 };
425 //
426 // compression
427 //
428 if ( minstrip < evento.base[l][m][pre] && minstrip != 100000.){
429 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
430 if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
431 };
432 };
433 }
434
435 void CaloFindBase(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
436 Float_t ominstrip = 100000.;
437 Float_t orms = 0.;
438 Float_t minstrip = 100000.;
439 Float_t rms = 0.;
440 evento.base[l][m][pre] = 0.;
441 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
442 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) {
443 ominstrip = minstrip;
444 minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
445 orms = rms;
446 rms = calib.calthr[l][m][pre];
447 };
448 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < ominstrip && evento.dexyc[l][m][e]-calib.calped[l][m][e] >minstrip ) {
449 ominstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
450 orms = calib.calthr[l][m][pre];
451 };
452 };
453 if ( minstrip != 100000. ) {
454 if ( ominstrip != 10000.) {
455 minstrip = ominstrip;
456 rms = orms;
457 };
458 Float_t strip6s = 0.;
459 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
460 if ( (evento.dexyc[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexyc[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
461 strip6s += 1.;
462 evento.base[l][m][pre] += (evento.dexyc[l][m][e] - calib.calped[l][m][e]);
463 };
464 //
465 // compression
466 //
467 if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
468 };
469 if ( strip6s >= 9. ){
470 Double_t arro = evento.base[l][m][pre]/strip6s;
471 Float_t deci = 1000.*((float)arro - float(int(arro)));
472 if ( deci < 500. ) {
473 arro = double(int(arro));
474 } else {
475 arro = 1. + double(int(arro));
476 };
477 evento.base[l][m][pre] = arro;
478 } else {
479 evento.base[l][m][pre] = 31000.;
480 };
481 } else {
482 evento.base[l][m][pre] = 31000.;
483 };
484 }
485
486 void CaloFindBaseRaw(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
487 Float_t minstrip = 100000.;
488 Float_t rms = 0.;
489 evento.base[l][m][pre] = 0.;
490 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
491 if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) {
492 minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
493 rms = calib.calthr[l][m][pre];
494 };
495 };
496 if ( minstrip != 100000. ) {
497 Float_t strip6s = 0.;
498 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
499 if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
500 strip6s += 1.;
501 evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
502 };
503 //
504 // compression
505 //
506 if ( abs((int)(evento.dexy[l][m][e]-calib.calped[l][m][e])) <= (minstrip+rms) ) {
507 evento.dexyc[l][m][e] = 0.;
508 } else {
509 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
510 };
511 };
512 if ( strip6s >= 9. ){
513 Double_t arro = evento.base[l][m][pre]/strip6s;
514 Float_t deci = 1000.*((float)arro - float(int(arro)));
515 if ( deci < 500. ) {
516 arro = double(int(arro));
517 } else {
518 arro = 1. + double(int(arro));
519 };
520 evento.base[l][m][pre] = arro;
521 } else {
522 evento.base[l][m][pre] = 31000.;
523 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
524 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
525 };
526 };
527 } else {
528 evento.base[l][m][pre] = 31000.;
529 };
530 }
531
532 void CaloFindBaseRawNC(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
533 Float_t minstrip = 100000.;
534 Float_t rms = 0.;
535 evento.base[l][m][pre] = 0.;
536 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
537 if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) {
538 minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
539 rms = calib.calthr[l][m][pre];
540 };
541 };
542 if ( minstrip != 100000. ) {
543 Float_t strip6s = 0.;
544 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
545 if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
546 strip6s += 1.;
547 evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
548 };
549 //
550 // compression
551 //
552 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
553 };
554 if ( strip6s >= 9. ){
555 Double_t arro = evento.base[l][m][pre]/strip6s;
556 Float_t deci = 1000.*((float)arro - float(int(arro)));
557 if ( deci < 500. ) {
558 arro = double(int(arro));
559 } else {
560 arro = 1. + double(int(arro));
561 };
562 evento.base[l][m][pre] = arro;
563 } else {
564 evento.base[l][m][pre] = 31000.;
565 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
566 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
567 };
568 };
569 } else {
570 evento.base[l][m][pre] = 31000.;
571 };
572 }
573
574 int CaloPede(TString filename, Int_t s, Int_t atime, Calib & calib){
575 TFile *calheadFile = 0;
576 calheadFile = getEmiFile(filename, "CalibCalPed", "Header");
577 TFile *calcalibFile = 0;
578 calcalibFile = getEmiFile(filename, "CalibCalPed","Event");
579 //
580 const char *filen = 0;
581 filen = filename;
582 if ( !calheadFile || !calcalibFile ) {
583 if ( calheadFile ) calheadFile->Close();
584 if ( calcalibFile ) calcalibFile->Close();
585 printf("\n\nERROR: not able to open file: %s \n",filen);
586 return(3);
587 };
588 TTree *tr = (TTree*)calheadFile->Get("Pscu");
589 tr->AddFriend("CalibCalPed", calcalibFile);
590 pamela::CalibCalPedEvent *ce = 0;
591 pamela::PscuHeader *cph = 0;
592 pamela::EventHeader *ceh = 0;
593 tr->SetBranchAddress("Header", &ceh);
594 tr->SetBranchAddress("Event", &ce);
595 Long64_t ncalibs = tr->GetEntries();
596 for (Int_t ci = 0; ci < ncalibs ; ci++){
597 tr->GetEntry(ci);
598 cph = ceh->GetPscuHeader();
599 if ( atime == cph->GetOrbitalTime()){
600 calib.iev = ce->iev;
601 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
602 for ( Int_t d=0 ; d<11 ;d++ ){
603 Int_t pre = -1;
604 for ( Int_t j=0; j<96 ;j++){
605 if ( j%16 == 0 ) pre++;
606 if ( s == 2 ){
607 calib.calped[0][2*d+1][j] = ce->calped[3][d][j];
608 calib.cstwerr[3] = ce->cstwerr[3];
609 calib.cperror[3] = ce->cperror[3];
610 calib.calgood[0][2*d+1][j] = ce->calgood[3][d][j];
611 calib.calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
612 calib.calrms[0][2*d+1][j] = ce->calrms[3][d][j];
613 calib.calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
614 calib.calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
615 };
616 if ( s == 3 ){
617 calib.calped[0][2*d][j] = ce->calped[1][d][j];
618 calib.cstwerr[1] = ce->cstwerr[1];
619 calib.cperror[1] = ce->cperror[1];
620 calib.calgood[0][2*d][j] = ce->calgood[1][d][j];
621 calib.calthr[0][2*d][pre] = ce->calthr[1][d][pre];
622 calib.calrms[0][2*d][j] = ce->calrms[1][d][j];
623 calib.calbase[0][2*d][pre] = ce->calbase[1][d][pre];
624 calib.calvar[0][2*d][pre] = ce->calvar[1][d][pre];
625 };
626 if ( s == 0 ){
627 calib.calped[1][2*d][j] = ce->calped[0][d][j];
628 calib.cstwerr[0] = ce->cstwerr[0];
629 calib.cperror[0] = ce->cperror[0];
630 calib.calgood[1][2*d][j] = ce->calgood[0][d][j];
631 calib.calthr[1][2*d][pre] = ce->calthr[0][d][pre];
632 calib.calrms[1][2*d][j] = ce->calrms[0][d][j];
633 calib.calbase[1][2*d][pre] = ce->calbase[0][d][pre];
634 calib.calvar[1][2*d][pre] = ce->calvar[0][d][pre];
635 };
636 if ( s == 1 ){
637 calib.calped[1][2*d+1][j] = ce->calped[2][d][j];
638 calib.cstwerr[2] = ce->cstwerr[2];
639 calib.cperror[2] = ce->cperror[2];
640 calib.calgood[1][2*d+1][j] = ce->calgood[2][d][j];
641 calib.calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
642 calib.calrms[1][2*d+1][j] = ce->calrms[2][d][j];
643 calib.calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
644 calib.calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
645 };
646 };
647 };
648 };
649 };
650 };
651 tr->Delete();
652 if ( calheadFile ) calheadFile->Close();
653 if ( calcalibFile ) calcalibFile->Close();
654 return(0);
655 }
656
657 TString getFilename(const TString filename){
658 const string fil = (const char*)filename;
659 Int_t posiz = fil.find("dw_");
660 if ( posiz == -1 ) posiz = fil.find("DW_");
661 if ( posiz == -1 ) return 0;
662 Int_t posiz2 = posiz+13;
663 TString file2;
664 stringcopy(file2,filename,posiz,posiz2);
665 TString pdat(".dat");
666 stringappend(file2,pdat);
667 return file2;
668 }
669
670 Int_t fetchpreviousfile(TString ffile, Int_t s){
671 const char *file = ffile;
672 stringstream fil;
673 fil.str("");
674 fil << file ;
675 Int_t posiz = fil.str().find("dw_");
676 Int_t upper = 0;
677 if ( posiz == -1 ) {
678 posiz = fil.str().find("DW_");
679 upper = 1;
680 };
681 if ( posiz == -1 ) return(1);
682 //
683 TString nomefile = getFilename(ffile);
684 const char *ufnome = nomefile;
685 stringstream nyear;
686 stringstream nmonth;
687 stringstream nday;
688 stringstream nfno;
689 TString year;
690 TString month;
691 TString day;
692 TString fno;
693 stringcopy(year, ufnome, 3, 5);
694 stringcopy(month, ufnome, 5, 7);
695 stringcopy(day, ufnome, 7, 9);
696 stringcopy(fno, ufnome, 10, 13);
697 const char *cyear = year;
698 const char *cmonth = month;
699 const char *cday = day;
700 const char *cfno = fno;
701 Int_t jy = atoi(cyear);
702 Int_t jm = atoi(cmonth);
703 Int_t jd = atoi(cday);
704 Int_t jn = atoi(cfno);
705 Int_t inter = 0;
706 char *zy;
707 char *zm;
708 char *zd;
709 char *zn = "";
710 char *zyn;
711 Int_t jyn = 0;
712 Int_t goodthis = 0;
713 Int_t onegood = 0;
714 Long64_t ncalibs;
715 TTree *tr;
716 pamela::CalibCalPedEvent *ce = 0;
717 pamela::PscuHeader *cph = 0;
718 pamela::EventHeader *ceh = 0;
719 if ( jn > 1 ) {
720 jn--;
721 } else {
722 jn = 999;
723 if ( jd > 1 ) {
724 jd--;
725 } else {
726 if ( jm > 1 ) {
727 jm--;
728 } else {
729 jm = 12;
730 if ( jy > 3 ) {
731 jy--;
732 } else {
733 return(1);
734 };
735 };
736 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
737 jd = 30;
738 } else {
739 jd = 31;
740 };
741 if ( jm==2 ) jd = 29;
742 };
743 };
744 Int_t gday = 0;
745 Int_t retvar = 1;
746 stringstream nyd;
747 while ( jy>3 ) {
748 while ( jm>0 ){
749 while ( jd > 0 ){
750 while ( jn > 0 ){
751 retvar--;
752 nyear.str("");
753 nyear << jy;
754 if ( jy<10 ){
755 zy = "0";
756 } else {
757 zy = "";
758 };
759 nmonth.str("");
760 nmonth << jm;
761 if ( jm<10 ){
762 zm = "0";
763 } else {
764 zm = "";
765 };
766 nday.str("");
767 nday << jd;
768 if ( jd<10 ){
769 zd = "0";
770 } else {
771 zd = "";
772 };
773 nfno.str("");
774 nfno << jn;
775 if ( jn>99 ) zn = "";
776 if ( jn<100 && jn>9 ) zn = "0";
777 if ( jn<10 ) zn = "00";
778 jyn = 0;
779 goodthis = 0;
780 onegood = 0;
781 while ( jyn < 100 ){
782 begwhile:
783 nyd.str("");
784 nyd << jyn;
785 if ( jyn<10 ){
786 zyn = "0";
787 } else {
788 zyn = "";
789 };
790 stringstream nfile;
791 nfile.str("");
792 if ( upper ){
793 nfile << "DW_" << zy;
794 } else {
795 nfile << "dw_" << zy;
796 };
797 nfile << nyear.str().c_str();
798 nfile << zm;
799 nfile << nmonth.str().c_str();
800 nfile << zd;
801 nfile << nday.str().c_str() << "_";
802 nfile << zn;
803 nfile << nfno.str().c_str();
804 nfile << zyn;
805 nfile << nyd.str().c_str();
806 //
807 TString fgfile = "";
808 stringcopy(fgfile,ffile,0,posiz);
809 stringappend(fgfile,nfile.str().c_str());
810 TFile *calheadFile = 0;
811 calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header");
812 TFile *calcalibFile = 0;
813 calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event");
814 //
815 if ( !calheadFile || !calcalibFile ) {
816 if ( calheadFile ) calheadFile->Close();
817 if ( calcalibFile ) calcalibFile->Close();
818 if ( jyn > 0 && onegood ) {
819 goodthis = 1;
820 jyn--;
821 goto begwhile;
822 };
823 if ( jyn == 0 ) {
824 onegood = 0;
825 jyn = 100;
826 goto salta;
827 };
828 if ( jyn > 0 && !onegood ) {
829 onegood = 0;
830 jyn = 100;
831 goto salta;
832 };
833 goto salta;
834 } else {
835 onegood = 1;
836 if ( !goodthis ) {
837 calheadFile->Close();
838 calcalibFile->Close();
839 goto salta;
840 };
841 };
842 //Takes the tree of the header file
843 tr = (TTree*)calheadFile->Get("Pscu");
844 tr->AddFriend("CalibCalPed", calcalibFile);
845 tr->SetBranchAddress("Header", &ceh);
846 tr->SetBranchAddress("Event", &ce);
847 ncalibs = tr->GetEntries();
848 if ( ncalibs == 0 ) {
849 jyn = 100;
850 onegood = 0;
851 goodthis = 0;
852 calheadFile->Close();
853 calcalibFile->Close();
854 goto salta;
855 };
856 inter = 0;
857 for (Int_t c = 0; c < ncalibs; c++){
858 tr->GetEntry(c);
859 cph = ceh->GetPscuHeader();
860 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
861 inter++;
862 };
863 };
864 if ( inter ){
865 jyn = 100;
866 calheadFile->Close();
867 calcalibFile->Close();
868 return(retvar);
869 };
870 calheadFile->Close();
871 calcalibFile->Close();
872 salta:
873 jyn++;
874 };
875 jn--;
876 };
877 jn = 999;
878 jd--;
879 gday++;
880 if ( gday>1 ) return(1);
881 };
882 jm--;
883 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
884 jd = 30;
885 } else {
886 jd = 31;
887 };
888 if ( jm==2 ) jd = 29;
889 };
890 jm = 12;
891 jy--;
892 };
893 return(retvar);
894 }
895
896 TString whatnamewith(TString ffile, Int_t retvar){
897 Int_t debug = 0;
898 const char *file = ffile;
899 stringstream fil;
900 fil.str("");
901 fil << file ;
902 if ( debug ) printf("ffile %s retvar %i \n",file,retvar);
903 Int_t posiz = fil.str().find("dw_");
904 Int_t upper = 0;
905 if ( posiz == -1 ) {
906 posiz = fil.str().find("DW_");
907 upper = 1;
908 };
909 if ( posiz == -1 ) return("wrong");
910 //
911 // printf("posiz %i \n",posiz);
912 TString nomefile = getFilename(ffile);
913 const char *ufnome = nomefile;
914 stringstream nyear;
915 stringstream nmonth;
916 stringstream nday;
917 stringstream nfno;
918 TString year;
919 TString month;
920 TString day;
921 TString fno;
922 stringcopy(year, ufnome, 3, 5);
923 stringcopy(month, ufnome, 5, 7);
924 stringcopy(day, ufnome, 7, 9);
925 stringcopy(fno, ufnome, 10, 13);
926 const char *cyear = year;
927 const char *cmonth = month;
928 const char *cday = day;
929 const char *cfno = fno;
930 Int_t jy = atoi(cyear);
931 Int_t jm = atoi(cmonth);
932 Int_t jd = atoi(cday);
933 Int_t jn = atoi(cfno);
934 Int_t jyn = 0;
935 Int_t goodthis = 0;
936 Int_t onegood = 0;
937 char *zy;
938 char *zm;
939 char *zd;
940 char *zn = "";
941 char *zyn;
942 Long64_t ncalibs;
943 TString fgfile;
944 TTree *tr;
945 TFile *calheadFile = 0;
946 TFile *calcalibFile = 0;
947 if ( jn > 1 ) {
948 jn--;
949 } else {
950 jn = 999;
951 if ( jd > 1 ) {
952 jd--;
953 } else {
954 if ( jm > 1 ) {
955 jm--;
956 } else {
957 jm = 12;
958 if ( jy > 3 ) {
959 jy--;
960 } else {
961 return("wrong");
962 };
963 };
964 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
965 jd = 30;
966 } else {
967 jd = 31;
968 };
969 if ( jm==2 ) jd = 29;
970 };
971 };
972 Int_t retval = 1;
973 Int_t gday = 0;
974 stringstream nyd;
975 while ( jy>3 ) {
976 while ( jm>0 ){
977 while ( jd > 0 ){
978 while ( jn > 0 ){
979 retval--;
980 if ( retval < retvar ) return("wrong");
981 nyear.str("");
982 nyear << jy;
983 if ( jy<10 ){
984 zy = "0";
985 } else {
986 zy = "";
987 };
988 nmonth.str("");
989 nmonth << jm;
990 if ( jm<10 ){
991 zm = "0";
992 } else {
993 zm = "";
994 };
995 nday.str("");
996 nday << jd;
997 if ( jd<10 ){
998 zd = "0";
999 } else {
1000 zd = "";
1001 };
1002 nfno.str("");
1003 nfno << jn;
1004 if ( jn>99 ) zn = "";
1005 if ( jn<100 && jn>9 ) zn = "0";
1006 if ( jn<10 ) zn = "00";
1007 jyn = 0;
1008 goodthis = 0;
1009 onegood = 0;
1010 while ( jyn < 100 ){
1011 begwhile:
1012 nyd.str("");
1013 nyd << jyn;
1014 if ( jyn<10 ){
1015 zyn = "0";
1016 } else {
1017 zyn = "";
1018 };
1019 stringstream nfile;
1020 nfile.str("");
1021 if ( upper ){
1022 nfile << "DW_" << zy;
1023 } else {
1024 nfile << "dw_" << zy;
1025 };
1026 nfile << nyear.str().c_str();
1027 nfile << zm;
1028 nfile << nmonth.str().c_str();
1029 nfile << zd;
1030 nfile << nday.str().c_str() << "_";
1031 nfile << zn;
1032 nfile << nfno.str().c_str();
1033 nfile << zyn;
1034 nfile << nyd.str().c_str();
1035 //
1036 fgfile = "";
1037 stringcopy(fgfile,ffile,0,posiz);
1038 stringappend(fgfile,nfile.str().c_str());
1039 const char *tosee = fgfile;
1040 if ( debug ) printf("file: %s retval %i retvar %i \n",tosee,retval,retvar);
1041 calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header");
1042 calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event");
1043 //
1044 if ( !calheadFile || !calcalibFile ) {
1045 if ( calheadFile ) calheadFile->Close();
1046 if ( calcalibFile ) calcalibFile->Close();
1047 if ( jyn > 0 && onegood ) {
1048 goodthis = 1;
1049 jyn--;
1050 goto begwhile;
1051 };
1052 if ( jyn == 0 ) {
1053 onegood = 0;
1054 jyn = 100;
1055 goto salta;
1056 };
1057 if ( jyn > 0 && !onegood ) {
1058 onegood = 0;
1059 jyn = 100;
1060 goto salta;
1061 };
1062 goto salta;
1063 } else {
1064 onegood = 1;
1065 if ( !goodthis ){
1066 calheadFile->Close();
1067 calcalibFile->Close();
1068 goto salta;
1069 };
1070 };
1071 tr = (TTree*)calheadFile->Get("Pscu");
1072 tr->AddFriend("CalibCalPed", calcalibFile);
1073 ncalibs = tr->GetEntries();
1074 if ( ncalibs == 0 ) {
1075 jyn = 100;
1076 onegood = 0;
1077 goodthis = 0;
1078 calheadFile->Close();
1079 calcalibFile->Close();
1080 goto salta;
1081 };
1082 calheadFile->Close();
1083 calcalibFile->Close();
1084 if ( retval == retvar ) return(fgfile);
1085 salta:
1086 jyn++;
1087 };
1088 jn--;
1089 };
1090 jn = 999;
1091 jd--;
1092 gday++;
1093 if ( gday>7 ) return("wrong");
1094 };
1095 jm--;
1096 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
1097 jd = 30;
1098 } else {
1099 jd = 31;
1100 };
1101 if ( jm==2 ) jd = 29;
1102 };
1103 jm = 12;
1104 jy--;
1105 };
1106 return("wrong");
1107 }
1108
1109 char *getLEVname(TString filename, TString detc, TString numb, TString frame="root"){
1110 // char *file;
1111 stringstream file;
1112 file.str("");
1113 const char *num = numb;
1114 const char *det = detc;
1115 const char *fra = frame;
1116 string fil = (const char *)filename;
1117 Int_t posiz = fil.find("dw_");
1118 if ( posiz == -1 ) posiz = fil.find("DW_");
1119 if ( posiz == -1 ) return(0);
1120 Int_t spos = posiz;
1121 Int_t epos = posiz+15;
1122 TString tmptempf;
1123 stringcopy(tmptempf,filename,spos,epos);
1124 const char *tempf = tmptempf;
1125 file << tempf << ".Physics.Level";
1126 file << num << ".";
1127 file << det << ".Event.";
1128 file << fra;
1129 const char *rfile = file.str().c_str();
1130 return((char*)rfile);
1131 };
1132
1133 char *getfileLEVname(TString filename, TString detc, TString numb, TString frame="root"){
1134 // char *file;
1135 stringstream file;
1136 const char *num = numb;
1137 const char *det = detc;
1138 const char *fra = frame;
1139 string fil = (const char *)filename;
1140 Int_t posiz = fil.find("dw_");
1141 if ( posiz == -1 ) posiz = fil.find("DW_");
1142 if ( posiz == -1 ) return(0);
1143 Int_t spos = posiz;
1144 Int_t epos = posiz+13;
1145 TString tmptempf;
1146 stringcopy(tmptempf,filename,spos,epos);
1147 const char *tempf = tmptempf;
1148 file.str("");
1149 file << tempf << ".Physics.Level";
1150 // file << "00.Physics.Level";
1151 file << num << ".";
1152 file << det << ".Event.";
1153 file << fra;
1154 const char *rfile = file.str().c_str();
1155 return((char*)rfile);
1156 };
1157
1158 void OLDCaloFindCalibs(TString &filename, Calib & calib){
1159 for (Int_t s = 0; s < 4; s++){
1160 for (Int_t d = 1; d<50; d++){
1161 calib.ttime[s][d] = 0;
1162 if ( d < 49 ) calib.time[s][d] = 0;
1163 };
1164 };
1165 TFile *calheadFile = 0;
1166 calheadFile = getEmiFile(filename, "CalibCalPed", "Header");
1167 TFile *calcalibFile = 0;
1168 calcalibFile = getEmiFile(filename, "CalibCalPed","Event");
1169 //
1170 if ( !calheadFile || !calcalibFile ) {
1171 if ( calheadFile ) calheadFile->Close();
1172 if ( calcalibFile ) calcalibFile->Close();
1173 printf("No calibration header file! Exiting...\n");
1174 return;
1175 };
1176 //Takes the tree of the header file
1177 TTree *tr = (TTree*)calheadFile->Get("Pscu");
1178 tr->AddFriend("CalibCalPed", calcalibFile);
1179 pamela::CalibCalPedEvent *ce = 0;
1180 pamela::PscuHeader *cph = 0;
1181 pamela::EventHeader *ceh = 0;
1182 tr->SetBranchAddress("Header", &ceh);
1183 tr->SetBranchAddress("Event", &ce);
1184 Long64_t ncalibs = tr->GetEntries();
1185 Int_t inter;
1186 for (Int_t s = 0; s < 4; s++){
1187 for (Int_t d = 1; d<50; d++){
1188 calib.ttime[s][d] = 0;
1189 if ( d < 49 ) calib.time[s][d] = 0;
1190 };
1191 };
1192 for (Int_t s = 0; s < 4; s++){
1193 inter = 0;
1194 for (Int_t c = 0; c < ncalibs; c++){
1195 tr->GetEntry(c);
1196 cph = ceh->GetPscuHeader();
1197 calib.ttime[s][inter] = 0;
1198 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
1199 calib.ttime[s][inter] = cph->GetOrbitalTime();
1200 inter++;
1201 } else {
1202 if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
1203 printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
1204 };
1205 };
1206 };
1207 if ( inter == 0 ){
1208 printf(" ERROR: no suitable calibration for section %i in this file!\n",s);
1209 };
1210 for (Int_t d = 1; d<50; d++){
1211 if ( calib.ttime[s][d] != 0 ) {
1212 calib.time[s][d-1] = calib.ttime[s][d-1] + (int)((calib.ttime[s][d] - calib.ttime[s][d-1])/2.);
1213 } else {
1214 if ( d == 1 ) {
1215 calib.time[s][d-1] = 0;
1216 };
1217 };
1218 };
1219 };
1220 calheadFile->Close();
1221 calcalibFile->Close();
1222 }
1223
1224 void CaloFindCalibs(TString &filename, TString &calibfile, Int_t &wused, Calib & calib){
1225 //
1226 // First of all check if the file has a monotone growing OBT time function.
1227 //
1228 Int_t debug = 0;
1229 Int_t obtjump = 0;
1230 Int_t evjump[50];
1231 Int_t OBT = 0;
1232 Int_t OOBT = 0;
1233 Int_t OBT1 = 0;
1234 Int_t LOBT = 0;
1235 pamela::PscuHeader *ph = 0;
1236 pamela::EventHeader *eh = 0;
1237 TFile *headerF = 0;
1238 headerF = emigetFile(filename, "Physics", "Header");
1239 if ( !headerF ) {
1240 printf(" No physics file! \n ");
1241 return;
1242 };
1243 TTree *ctr = (TTree*)headerF->Get("Pscu");
1244 Long64_t nevents = ctr->GetEntries();
1245 if ( nevents == 0 ) {
1246 printf(" The file is empty! \n ");
1247 return;
1248 };
1249 for (Int_t i = 0; i < nevents; i++){
1250 ctr->SetBranchAddress("Header", &eh);
1251 ctr->GetEntry(i);
1252 ph = eh->GetPscuHeader();
1253 OBT = ph->GetOrbitalTime();
1254 if ( !i ) OBT1 = OBT;
1255 if ( i == nevents-1 ) LOBT = OBT;
1256 if ( OOBT > OBT ) {
1257 evjump[obtjump] = i;
1258 obtjump++;
1259 };
1260 OOBT = OBT;
1261 };
1262 printf("\n Obtjump = %i - FIRST OBT %i - LAST OBT %i \n\n",obtjump,OBT1,LOBT);
1263 //
1264 TTree *tr;
1265 pamela::CalibCalPedEvent *ce = 0;
1266 pamela::PscuHeader *cph = 0;
1267 pamela::EventHeader *ceh = 0;
1268 Long64_t ncalibs;
1269 //
1270 for (Int_t s = 0; s < 4; s++){
1271 for (Int_t d = 1; d<50; d++){
1272 calib.ttime[s][d] = 0;
1273 if ( d < 49 ) calib.time[s][d] = 0;
1274 };
1275 };
1276 for (Int_t s = 0; s < 4; s++){
1277 //
1278 Int_t firstlook = 0;
1279 repeatsearch:
1280 //
1281 Int_t inter = 0;
1282 TFile *calheadFile = 0;
1283 calheadFile = getEmiFile(filename, "CalibCalPed", "Header");
1284 TFile *calcalibFile = 0;
1285 calcalibFile = getEmiFile(filename, "CalibCalPed","Event");
1286 //
1287 if ( !calheadFile || !calcalibFile ) {
1288 if ( calheadFile ) calheadFile->Close();
1289 if ( calcalibFile ) calcalibFile->Close();
1290 printf(" No calibration in this file! \n");
1291 if ( !firstlook ){
1292 wused = 2;
1293 filename = calibfile;
1294 firstlook = 1;
1295 goto repeatsearch;
1296 } else {
1297 wused = 0;
1298 goto jte;
1299 };
1300 };
1301 wused = 1;
1302 //Takes the tree of the header file
1303 tr = (TTree*)calheadFile->Get("Pscu");
1304 tr->AddFriend("CalibCalPed", calcalibFile);
1305 tr->SetBranchAddress("Header", &ceh);
1306 tr->SetBranchAddress("Event", &ce);
1307 //
1308 ncalibs = tr->GetEntries();
1309 calib.obtjump = 0;
1310 inter = 0;
1311 for (Int_t c = 0; c < ncalibs; c++){
1312 tr->GetEntry(c);
1313 cph = ceh->GetPscuHeader();
1314 calib.ttime[s][inter] = 0;
1315 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
1316 calib.ttime[s][inter] = cph->GetOrbitalTime();
1317 inter++;
1318 } else {
1319 if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
1320 printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
1321 };
1322 };
1323 };
1324 jte:
1325 if ( inter == 0 ){
1326 printf(" WARNING: no suitable calibration for section %i in this file!\n",s);
1327 printf(" I WILL SEARCH IN PREVIOUS FILES\n");
1328 };
1329 if ( inter > 50 ){
1330 printf(" WARNING: cannot handle more than 50 calibrations for file!\n");
1331 printf(" I WILL SEARCH IN PREVIOUS FILES\n");
1332 inter = 0;
1333 };
1334 if ( obtjump ){
1335 calib.obtjump = 1;
1336 printf(" ERROR: jump in OBT! CPU restarded? \n");
1337 printf(" WARNING: I will use only the first calibration (if any)\n");
1338 };
1339 //
1340 Int_t lastcal = 0;
1341 Int_t of = 0;
1342 if ( OBT1 < calib.ttime[s][0] || !inter ) {
1343 //
1344 // here we must look for a calibration in the previous file...
1345 //
1346 TString ffile = "";
1347 stringcopy(ffile,filename);
1348 Int_t fcodep = fetchpreviousfile(ffile,s);
1349 if ( fcodep == 1 ) {
1350 if ( inter ){
1351 of = -3;
1352 printf(" WARNING: Problems to find a good calibration for section %i \n",s);
1353 } else {
1354 of = -2;
1355 printf(" ERROR: I was not able to find any good calibration for section %i \n",s);
1356 };
1357 } else {
1358 stringcopy(ffile,filename);
1359 TString pfile = whatnamewith(ffile,fcodep);
1360 struct Calib tcal;
1361 OLDCaloFindCalibs(pfile,tcal);
1362 for ( Int_t d = 0; d<51; d++ ){
1363 if ( tcal.ttime[s][d] == 0 ){
1364 lastcal = d-1;
1365 break;
1366 };
1367 };
1368 if ( !inter ){
1369 calib.ttime[s][0] = tcal.ttime[s][lastcal];
1370 calib.time[s][0] = OBT1;
1371 calib.time[s][1] = LOBT;
1372 calib.fcode[s][0] = fcodep;
1373 if ( debug ) printf("boh! %i \n",fcodep);
1374 if ( debug ) getchar();
1375 of = -1;
1376 } else {
1377 if ( obtjump ){
1378 if ( debug ) printf("boh2! %i \n",fcodep);
1379 if ( debug ) getchar();
1380 inter = 1;
1381 calib.time[s][0] = OBT1;
1382 calib.time[s][1] = LOBT;
1383 calib.fcode[s][0] = 10;
1384 of = -5;
1385 for (Int_t d = 1; d<50; d++){
1386 calib.ttime[s][d] = 0;
1387 if ( d>1 && d<49 ) calib.time[s][d] = 0;
1388 };
1389 } else {
1390 if ( lastcal >= 0 && tcal.ttime[s][lastcal]<OBT1 && (OBT1-tcal.ttime[s][lastcal])<7200000 && (OBT1-tcal.ttime[s][lastcal])>0 ){
1391 for (Int_t i = inter+1; i>0; i--){
1392 calib.ttime[s][i] = calib.ttime[s][i-1];
1393 };
1394 calib.time[s][0] = OBT1;
1395 calib.ttime[s][0] = tcal.ttime[s][lastcal];
1396 calib.fcode[s][0] = fcodep;
1397 of = 2;
1398 } else {
1399 calib.time[s][0] = OBT1;
1400 calib.fcode[s][0] = 10;
1401 of = 1;
1402 };
1403 };
1404 };
1405 };
1406 } else {
1407 of = 0;
1408 };
1409 if ( debug ) printf("boh3! %i \n",of);
1410 Int_t of1 = 0;
1411 if ( inter && of != 0 && of != 2 && of != -2 && of != -5) {
1412 for (Int_t i = inter+1; i>0; i--){
1413 calib.ttime[s][i] = calib.ttime[s][i-1];
1414 };
1415 of1 = 0;
1416 };
1417 if ( of == -2 ) of = -1;
1418 if ( of == -3 ) {
1419 of = 1;
1420 calib.fcode[s][0] = 10;
1421 };
1422 //
1423 if ( of != -1 && of != -5 ){
1424 if ( of == 2 ) of = 1;
1425 for (Int_t d = 0; d<inter+1; d++){
1426 calib.time[s][d+of] = calib.ttime[s][d+of];
1427 calib.fcode[s][d+of] = 10;
1428 if ( d == inter ){
1429 if ( debug ) printf("boh4! %i \n",of);
1430 calib.time[s][d+of1+of] = LOBT;
1431 calib.fcode[s][d+of] = 0;
1432 };
1433 };
1434 };
1435 if ( calheadFile ) calheadFile->Close();
1436 if ( calcalibFile ) calcalibFile->Close();
1437 };
1438 }
1439
1440 void Calo1stcalib(TString filename, Calib & calib, Int_t b[4]){
1441 Float_t estrip[2][22][96];
1442 //
1443 // this is the value of the mip for each strip. To be changed when we will have the real values
1444 //
1445 for (Int_t s=0; s<4;s++){
1446 for (Int_t d = 0; d<50; d++){
1447 calib.ttime[s][d] = 0 ;
1448 if ( d < 49 ) calib.time[s][d] = 0 ;
1449 };
1450 };
1451 for (Int_t m = 0; m < 2 ; m++ ){
1452 for (Int_t k = 0; k < 22; k++ ){
1453 for (Int_t l = 0; l < 96; l++ ){
1454 calib.calped[m][k][l] = 0. ;
1455 estrip[m][k][l] = 0.;
1456 };
1457 };
1458 }
1459 //
1460 // first of all find the calibrations in the file
1461 //
1462 OLDCaloFindCalibs(filename, calib);
1463 //
1464 // print on the screen the results:
1465 //
1466 printf(" ---------------------------------------------------------- \n \n");
1467 Int_t calibex = 0;
1468 for (Int_t s=0; s<4;s++){
1469 Int_t stop = 0;
1470 for (Int_t d = 0; d<48; d++){
1471 if ( calib.ttime[s][d] != 0 ) calibex++;
1472 if ( calib.time[s][0] != 0 ){
1473 if ( d == 0 ) printf(" Section %i from time 0 to time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d]);
1474 if ( calib.time[s][d+1] != 0 ) {
1475 printf(" Section %i from time %i to time %i use calibration at time %i \n",s,calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d+1]);
1476 } else {
1477 if ( !stop ){
1478 printf(" Section %i from time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d+1]);
1479 stop = 1;
1480 };
1481 };
1482 } else {
1483 if ( calib.ttime[s][d] != 0 ) printf(" Section %i from time 0 use calibration at time %i \n",s,calib.ttime[s][d]);
1484 };
1485 };
1486 printf("\n");
1487 };
1488 printf(" ---------------------------------------------------------- \n");
1489 if ( calibex < 1 ) {
1490 printf("No full calibration data in this file, sorry!\n");
1491 } else {
1492 //
1493 // calibrate before starting
1494 //
1495 for (Int_t s = 0; s < 4; s++){
1496 b[s]=0;
1497 CaloPede(filename,s,calib.ttime[s][b[s]],calib);
1498 };
1499 };
1500 }
1501
1502
1503 void ColorMIP(Float_t mip, int& colo){
1504 // printf("mip = %f \n",mip);
1505 if ( colo > 0 ){
1506 colo = 10;
1507 if ( mip > 0.7 ) colo = 38;
1508 if ( mip > 2. ) colo = 4;
1509 if ( mip > 10. ) colo = 3;
1510 if ( mip > 100. ) colo = 2;
1511 if ( mip > 500. ) colo = 6;
1512 } else {
1513 colo = 10;
1514 if ( mip > 0.7 ) colo = 17;
1515 if ( mip > 2. ) colo = 15;
1516 if ( mip > 10. ) colo = 14;
1517 if ( mip > 100. ) colo = 13;
1518 if ( mip > 500. ) colo = 12;
1519 };
1520 }
1521
1522 void ColorTOFMIP(Float_t mip, int& colo){
1523 // printf("mip = %f \n",mip);
1524 if ( colo > 0 ){
1525 colo = 10;
1526 if ( mip > 0. ) colo = 38;
1527 if ( mip > 2. ) colo = 4;
1528 if ( mip > 10. ) colo = 3;
1529 if ( mip > 100. ) colo = 2;
1530 if ( mip > 500. ) colo = 6;
1531 } else {
1532 colo = 10;
1533 if ( mip > 0. ) colo = 17;
1534 if ( mip > 2. ) colo = 15;
1535 if ( mip > 10. ) colo = 14;
1536 if ( mip > 100. ) colo = 13;
1537 if ( mip > 500. ) colo = 12;
1538 };
1539 }
1540
1541 Int_t WhatToDo(int& i, int& doflag, int& jumpto, Long64_t nevents, Float_t lastevno, Float_t firstevno, const char *file, TString outDir, TString figty, TCanvas &figure) {
1542 char *bw;
1543 if ( jumpto == -1 ){
1544 bw = "bw";
1545 } else {
1546 bw = "";
1547 };
1548 jumpto = 0;
1549 char tellme[256];
1550 stringstream input;
1551 char tellme2[256];
1552 stringstream input2;
1553 char tellme3[256];
1554 stringstream input3;
1555 stringstream out;
1556 out.str("x");
1557 input2.str("zzzzzzzzzzzzzz");
1558 input3.str("z");
1559 input << out.str().c_str();
1560 stringstream stc;
1561 stringstream stc2;
1562 while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
1563 printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nparticular event number, p<enter> to save the figure, q<enter> to quit: \n");
1564 cin.getline(tellme,256);
1565 //
1566 input.str("");
1567 input << tellme;
1568 out << "y";
1569 //
1570 stc.str("b");
1571 //
1572 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1573 if ( i > 0 ) {
1574 printf("WARNING: going backward!\n\n");
1575 doflag = 2;
1576 } else {
1577 printf("This is the first event, you can't go backward! \n");
1578 out.str("");
1579 out << input.str().c_str();
1580 };
1581 };
1582 stc.str("j");
1583 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1584 printf("\n Do you want to jump to a progressive number [p] or to an event number [e]? ");
1585 cin.getline(tellme3,256);
1586 input3.str("");
1587 input3 << tellme3;
1588 stc.str("p");
1589 if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
1590 printf("\n Enter the progressive number you want to jump to: ");
1591 cin.getline(tellme2,256);
1592 input2.str("");
1593 input2 << tellme2;
1594 Int_t j;
1595 j = atoi(input2.str().c_str());
1596 if ( j < 1 || j > nevents+1 ) {
1597 printf("\n You can choose between 1 and %i \n",(int)nevents+1);
1598 out.str("");
1599 out << input.str().c_str();
1600 } else {
1601 printf("\n Jumping to progressive number %i\n\n",j);
1602 i = j-2;
1603 };
1604 };
1605 stc.str("e");
1606 if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
1607 printf("\n Enter the event number you want to jump to: ");
1608 cin.getline(tellme2,256);
1609 input2.str("");
1610 input2 << tellme2;
1611 Int_t j;
1612 j = atoi(input2.str().c_str());
1613 if ( j < firstevno || j > lastevno ) {
1614 printf("\n You can choose between %i and %i \n",(int)firstevno,(int)lastevno);
1615 out.str("");
1616 out << input.str().c_str();
1617 } else {
1618 printf("\n Jumping to event number %i\n\n",j);
1619 jumpto = j;
1620 i = -1;
1621 };
1622 };
1623 stc.str("e");
1624 stc2.str("p");
1625 if ( strcmp(input3.str().c_str(),stc.str().c_str()) && strcmp(input3.str().c_str(),stc2.str().c_str()) ) {
1626 printf(" You must type or \"p\" or \"e\"\n");
1627 out.str("");
1628 out << input.str().c_str();
1629 };
1630 };
1631 //
1632 stc.str("q");
1633 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1634 printf("Exiting...\n");
1635 return(1);
1636 };
1637 //
1638 stc.str("p");
1639 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1640 printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
1641 cin.getline(tellme2,256);
1642 input2.str("");
1643 input2 << tellme2;
1644 out.str("");
1645 out << input.str().c_str();
1646 //
1647 TString filename = file;
1648 const string fil = (const char*)filename;
1649 Int_t posiz = fil.find("dw_");
1650 if ( posiz == -1 ) posiz = fil.find("DW_");
1651 Int_t posiz2 = posiz+13;
1652 TString file2;
1653 stringcopy(file2,filename,posiz,posiz2);
1654 const char *figrec = file2;
1655 const char *outdir = outDir;
1656 stringstream figsave;
1657 const char *ty = figty;
1658 figsave.str("");
1659 figsave << outdir << "/";
1660 figsave << ty << "_";
1661 figsave << (i+1) << "_";
1662 figsave << figrec;
1663 figsave << bw << ".";
1664 figsave << input2.str().c_str();
1665 figure.SaveAs(figsave.str().c_str());
1666 printf("\n");
1667 };
1668 };
1669 return(0);
1670 }
1671
1672 void PrintFigure(TString filename, TString outDir, TString figty, TString saveas, TCanvas& figure) {
1673 const string fil = (const char*)filename;
1674 Int_t posiz = fil.find("dw_");
1675 if ( posiz == -1 ) posiz = fil.find("DW_");
1676 Int_t posiz2 = posiz+13;
1677 TString file2;
1678 stringcopy(file2,filename,posiz,posiz2);
1679 //
1680 const char *figrec = file2;
1681 stringstream figsave;
1682 figsave.str("");
1683 figsave << outDir.Data() << "/";
1684 figsave << figrec << "_";
1685 figsave << figty.Data() << ".";
1686 figsave << saveas.Data();
1687 figure.SaveAs(figsave.str().c_str());
1688 return;
1689 }
1690
1691 Double_t langaufun(Double_t *x, Double_t *par) {
1692 // Numeric constants
1693 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
1694 Double_t mpshift = -0.22278298; // Landau maximum location
1695 // Control constants
1696 Double_t np = 100.0; // number of convolution steps
1697 Double_t sc = 3.; // convolution extends to +-sc Gaussian sigmas
1698 // Variables
1699 Double_t xx;
1700 Double_t mpc;
1701 Double_t fland;
1702 Double_t sum = 0.0;
1703 Double_t xlow,xupp;
1704 Double_t step;
1705 Double_t i;
1706 // MP shift correction
1707 mpc = par[1] - mpshift * par[0];
1708 // Range of convolution integral
1709 xlow = x[0] - sc * par[3];
1710 xupp = x[0] + sc * par[3];
1711 step = (xupp-xlow) / np;
1712 // Convolution integral of Landau and Gaussian by sum
1713 for(i=1.0; i<=np/2; i++) {
1714 xx = xlow + (i-.5) * step;
1715 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
1716 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1717 xx = xupp - (i-.5) * step;
1718 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
1719 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1720 };
1721 return (par[2] * step * sum * invsq2pi / par[3]);
1722 }
1723
1724
1725
1726 TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF, TString drw="QR0"){
1727 // Once again, here are the Landau * Gaussian parameters:
1728 // par[0]=Width (scale) parameter of Landau density
1729 // par[1]=Most Probable (MP, location) parameter of Landau density
1730 // par[2]=Total area (integral -inf to inf, normalization constant)
1731 // par[3]=Width (sigma) of convoluted Gaussian function
1732 //
1733 // Variables for langaufit call:
1734 // his histogram to fit
1735 // fitrange[2] lo and hi boundaries of fit range
1736 // startvalues[4] reasonable start values for the fit
1737 // parlimitslo[4] lower parameter limits
1738 // parlimitshi[4] upper parameter limits
1739 // fitparams[4] returns the final fit parameters
1740 // fiterrors[4] returns the final fit errors
1741 // ChiSqr returns the chi square
1742 // NDF returns ndf
1743 Int_t i;
1744 Char_t FunName[100];
1745 sprintf(FunName,"Fitfcn_%s",his->GetName());
1746 //
1747 TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
1748 if (ffitold) delete ffitold;
1749 //
1750 TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
1751 ffit->SetParameters(startvalues);
1752 ffit->SetParNames("Width","MP","Area","GSigma");
1753 //
1754 for (i=0; i<4; i++) {
1755 ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
1756 };
1757 his->Fit(FunName,drw); // fit within specified range, use ParLimits, do not plot
1758 //
1759 ffit->GetParameters(fitparams); // obtain fit parameters
1760 for (i=0; i<4; i++) {
1761 fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
1762 }
1763 ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
1764 NDF[0] = ffit->GetNDF(); // obtain ndf
1765 //
1766 return (ffit); // return fit function
1767 }
1768
1769
1770 //Double_t langaupro(Double_t *params, Double_t &maxx) {
1771 Double_t langaupro(Double_t *params) {
1772 // Seaches for the location (x value) at the maximum of the
1773 // Landau-Gaussian convolute and its full width at half-maximum.
1774 //
1775 Double_t p,x;
1776 Double_t step;
1777 Double_t l,lold;
1778 Int_t i = 0;
1779 Int_t MAXCALLS = 5000;
1780 // Search for maximum
1781 p = params[1] - 0.1 * params[0];
1782 step = 0.05 * params[0];
1783 lold = -2.0;
1784 l = -1.0;
1785 // while ( (l != lold) && (i < MAXCALLS) ) {
1786 while ( fabs(l-lold)>10e-12 && (i < MAXCALLS) ) {
1787 i++;
1788 lold = l;
1789 x = p + step;
1790 l = langaufun(&x,params);
1791 if ( l < lold ) step = -step/10.;
1792 p += step;
1793 // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
1794 // printf(" l %.16g \n",l);
1795 };
1796 // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
1797 //printf(" l %.16g \n",l);
1798 if (i == MAXCALLS) return (-1.);
1799 //
1800 return(x);
1801 }
1802
1803
1804 void delay(int selftrig, int &del){
1805 int i;
1806 i = 0;
1807 del = 300;
1808 //
1809 while ((i < 16) && ((selftrig & (0x8000>>i)) == 0)) {
1810 del += 50;
1811 i++;
1812 };
1813 if (i > 15) del = 1100;
1814 }
1815
1816 Double_t fitraw(Double_t *x, Double_t *par){
1817 //
1818 Double_t fitval = par[3] + par[0]*TMath::Exp(-par[2]*(x[0]-par[1]));
1819 //
1820 return fitval;
1821 }
1822
1823
1824
1825 int checkmerging(Chmerge & chmerge){
1826 for (Int_t i = chmerge.fromflag; i<1000; i++){
1827 if ( chmerge.pscuOBT == 0 && chmerge.pscucount == 0 ) return(1);
1828 if ( chmerge.pscuOBT == chmerge.lastOBT[i] && chmerge.pscucount == chmerge.lastcount[i] ){
1829 chmerge.fromflag = i;
1830 printf("\n WARNING! %i event number %i at OBT %i already processed, skipped!\n",i,chmerge.pscucount,chmerge.pscuOBT);
1831 return(1);
1832 };
1833 };
1834 return(0);
1835 }
1836
1837 int fillmerging(Chmerge & chmerge){
1838 chmerge.lastOBT[chmerge.scanflag] = chmerge.pscuOBT;
1839 chmerge.lastcount[chmerge.scanflag] = chmerge.pscucount;
1840 chmerge.scanflag++;
1841 if ( chmerge.scanflag > 999 ) return(0);
1842 chmerge.lastOBT[chmerge.scanflag] = 0;
1843 chmerge.lastcount[chmerge.scanflag] = 0;
1844 return(0);
1845 }
1846
1847 //short int checkifpretyodaversion(TString filename){
1848 // string fil = (const char *)filename;
1849 // Int_t posiz = fil.find("dw_");
1850 // if ( posiz == -1 ) posiz = fil.find("DW_");
1851 // if ( posiz == -1 ) return (-1);
1852 // TString year;
1853 // stringcopy(year,filename,posiz+3,posiz+5);
1854 // const char *ye = year;
1855 // Int_t y = atoi(ye);
1856 // TString month;
1857 // stringcopy(month,filename,posiz+5,posiz+7);
1858 // const char *mo = month;
1859 // Int_t m = atoi(mo);
1860 // TString day;
1861 // stringcopy(day,filename,posiz+7,posiz+9);
1862 // const char *da = day;
1863 // Int_t d = atoi(da);
1864 // TString fno;
1865 // stringcopy(fno,filename,posiz+10,posiz+13);
1866 // const char *fn = fno;
1867 // Int_t f = atoi(fn);
1868 // if ( y >= 5 && m >= 5 && d < 15 ) return(1);
1869 // if ( y == 5 && m == 5 && d == 15 && f < 7 ) return(1);
1870 // return(0);
1871 //}
1872
1873

  ViewVC Help
Powered by ViewVC 1.1.23