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