1 |
// |
2 |
// Calorimeter useful functions and subroutines - Emiliano Mocchiutti |
3 |
// |
4 |
// FCaloFUNCTIONS.cxx version 1.00 (2006-03-03) |
5 |
// |
6 |
// Programs in this file are called by other programs and cannot be run by hand. |
7 |
// |
8 |
// Changelog: |
9 |
// |
10 |
// 1.00 - 1.01 (2006-03-03): read unique YODA file and compile correctly with ROOT classes etc. Small changes in the fetchpreviousfile/whatnamewith routines. |
11 |
// |
12 |
// 0.00 - 1.00 (2006-03-03) Clone of CaloFunctions.h |
13 |
// |
14 |
#include <iostream> |
15 |
#include <sstream> |
16 |
#include <fstream> |
17 |
#include <TSystem.h> |
18 |
#include <TH1.h> |
19 |
#include <TF1.h> |
20 |
#include <TROOT.h> |
21 |
#include <TStyle.h> |
22 |
#include <TClassEdit.h> |
23 |
#include <TObject.h> |
24 |
#include <TList.h> |
25 |
#include <TSystem.h> |
26 |
#include <TSystemDirectory.h> |
27 |
#include <TString.h> |
28 |
#include <TTree.h> |
29 |
#include <TChain.h> |
30 |
#include <TFile.h> |
31 |
#include <TClass.h> |
32 |
#include <TCanvas.h> |
33 |
#include <TKey.h> |
34 |
#include <TClassTable.h> |
35 |
#include <TThread.h> |
36 |
// |
37 |
#include <PamelaRun.h> |
38 |
#include <physics/calorimeter/CalorimeterEvent.h> |
39 |
#include <CalibCalPedEvent.h> |
40 |
// |
41 |
#include <fcalostructs.h> |
42 |
#include <FCaloFUNCTIONSfun.h> |
43 |
#include <caloclassesfun.h> |
44 |
// |
45 |
using namespace std; |
46 |
|
47 |
void stringcopy(TString& s1, const TString& s2, Int_t from=0, Int_t to=0){ |
48 |
if ( to == 0 ){ |
49 |
Int_t t2length = s2.Length(); |
50 |
s1 = ""; |
51 |
to = t2length; |
52 |
}; |
53 |
for (Int_t i = from; i<to; i++){ |
54 |
s1.Append(s2[i],1); |
55 |
}; |
56 |
} |
57 |
|
58 |
void stringappend(TString& s1, const TString& s2){ |
59 |
Int_t t2length = s2.Length(); |
60 |
for (Int_t i = 0; i<t2length; i++){ |
61 |
s1.Append(s2[i],1); |
62 |
}; |
63 |
} |
64 |
|
65 |
void CaloCompressData(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ |
66 |
Float_t minstrip = 100000.; |
67 |
Float_t rms = 0.; |
68 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
69 |
if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) { |
70 |
minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; |
71 |
rms = calib.calthr[l][m][pre]; |
72 |
}; |
73 |
}; |
74 |
// |
75 |
// compression |
76 |
// |
77 |
if ( minstrip < evento.base[l][m][pre] && minstrip != 100000.){ |
78 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
79 |
if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.; |
80 |
}; |
81 |
}; |
82 |
} |
83 |
|
84 |
void CaloFindBase(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ |
85 |
Float_t ominstrip = 100000.; |
86 |
Float_t orms = 0.; |
87 |
Float_t minstrip = 100000.; |
88 |
Float_t rms = 0.; |
89 |
evento.base[l][m][pre] = 0.; |
90 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
91 |
if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) { |
92 |
ominstrip = minstrip; |
93 |
minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; |
94 |
orms = rms; |
95 |
rms = calib.calthr[l][m][pre]; |
96 |
}; |
97 |
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 ) { |
98 |
ominstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; |
99 |
orms = calib.calthr[l][m][pre]; |
100 |
}; |
101 |
}; |
102 |
if ( minstrip != 100000. ) { |
103 |
if ( ominstrip != 10000.) { |
104 |
minstrip = ominstrip; |
105 |
rms = orms; |
106 |
}; |
107 |
Float_t strip6s = 0.; |
108 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
109 |
if ( (evento.dexyc[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexyc[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { |
110 |
strip6s += 1.; |
111 |
evento.base[l][m][pre] += (evento.dexyc[l][m][e] - calib.calped[l][m][e]); |
112 |
}; |
113 |
// |
114 |
// compression |
115 |
// |
116 |
if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.; |
117 |
}; |
118 |
if ( strip6s >= 9. ){ |
119 |
Double_t arro = evento.base[l][m][pre]/strip6s; |
120 |
Float_t deci = 1000.*((float)arro - float(int(arro))); |
121 |
if ( deci < 500. ) { |
122 |
arro = double(int(arro)); |
123 |
} else { |
124 |
arro = 1. + double(int(arro)); |
125 |
}; |
126 |
evento.base[l][m][pre] = arro; |
127 |
} else { |
128 |
evento.base[l][m][pre] = 31000.; |
129 |
}; |
130 |
} else { |
131 |
evento.base[l][m][pre] = 31000.; |
132 |
}; |
133 |
} |
134 |
|
135 |
void CaloFindBaseRaw(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ |
136 |
Float_t minstrip = 100000.; |
137 |
Float_t rms = 0.; |
138 |
evento.base[l][m][pre] = 0.; |
139 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
140 |
if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) { |
141 |
minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e]; |
142 |
rms = calib.calthr[l][m][pre]; |
143 |
}; |
144 |
}; |
145 |
if ( minstrip != 100000. ) { |
146 |
Float_t strip6s = 0.; |
147 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
148 |
if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { |
149 |
strip6s += 1.; |
150 |
evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]); |
151 |
}; |
152 |
// |
153 |
// compression |
154 |
// |
155 |
if ( abs((int)(evento.dexy[l][m][e]-calib.calped[l][m][e])) <= (minstrip+rms) ) { |
156 |
evento.dexyc[l][m][e] = 0.; |
157 |
} else { |
158 |
evento.dexyc[l][m][e] = evento.dexy[l][m][e]; |
159 |
}; |
160 |
}; |
161 |
if ( strip6s >= 9. ){ |
162 |
Double_t arro = evento.base[l][m][pre]/strip6s; |
163 |
Float_t deci = 1000.*((float)arro - float(int(arro))); |
164 |
if ( deci < 500. ) { |
165 |
arro = double(int(arro)); |
166 |
} else { |
167 |
arro = 1. + double(int(arro)); |
168 |
}; |
169 |
evento.base[l][m][pre] = arro; |
170 |
} else { |
171 |
evento.base[l][m][pre] = 31000.; |
172 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
173 |
evento.dexyc[l][m][e] = evento.dexy[l][m][e]; |
174 |
}; |
175 |
}; |
176 |
} else { |
177 |
evento.base[l][m][pre] = 31000.; |
178 |
}; |
179 |
} |
180 |
|
181 |
void CaloFindBaseRawNC(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ |
182 |
Float_t minstrip = 100000.; |
183 |
Float_t rms = 0.; |
184 |
evento.base[l][m][pre] = 0.; |
185 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
186 |
if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) { |
187 |
minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e]; |
188 |
rms = calib.calthr[l][m][pre]; |
189 |
}; |
190 |
}; |
191 |
if ( minstrip != 100000. ) { |
192 |
Float_t strip6s = 0.; |
193 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
194 |
if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { |
195 |
strip6s += 1.; |
196 |
evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]); |
197 |
}; |
198 |
// |
199 |
// compression |
200 |
// |
201 |
evento.dexyc[l][m][e] = evento.dexy[l][m][e]; |
202 |
}; |
203 |
if ( strip6s >= 9. ){ |
204 |
Double_t arro = evento.base[l][m][pre]/strip6s; |
205 |
Float_t deci = 1000.*((float)arro - float(int(arro))); |
206 |
if ( deci < 500. ) { |
207 |
arro = double(int(arro)); |
208 |
} else { |
209 |
arro = 1. + double(int(arro)); |
210 |
}; |
211 |
evento.base[l][m][pre] = arro; |
212 |
} else { |
213 |
evento.base[l][m][pre] = 31000.; |
214 |
for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ |
215 |
evento.dexyc[l][m][e] = evento.dexy[l][m][e]; |
216 |
}; |
217 |
}; |
218 |
} else { |
219 |
evento.base[l][m][pre] = 31000.; |
220 |
}; |
221 |
} |
222 |
|
223 |
Bool_t existfile(TString filename){ |
224 |
ifstream myfile; |
225 |
myfile.open(filename.Data()); |
226 |
if ( !myfile ){ |
227 |
// printf("\n\nERROR: not able to open file: %s \n",filename.Data()); |
228 |
return(false); |
229 |
}; |
230 |
myfile.close(); |
231 |
return(true); |
232 |
} |
233 |
|
234 |
int CaloPede(TString filename, Int_t s, Int_t atime, Calib & calib){ |
235 |
if ( !existfile(filename) ) return(3); |
236 |
// |
237 |
TFile *File = new TFile(filename.Data()); |
238 |
TTree *tr = (TTree*)File->Get("CalibCalPed"); |
239 |
// |
240 |
pamela::CalibCalPedEvent *ce = 0; |
241 |
pamela::PscuHeader *cph = 0; |
242 |
pamela::EventHeader *ceh = 0; |
243 |
tr->SetBranchAddress("Header", &ceh); |
244 |
tr->SetBranchAddress("CalibCalPed", &ce); |
245 |
// |
246 |
Long64_t ncalibs = tr->GetEntries(); |
247 |
for (Int_t ci = 0; ci < ncalibs ; ci++){ |
248 |
tr->GetEntry(ci); |
249 |
cph = ceh->GetPscuHeader(); |
250 |
if ( atime == cph->GetOrbitalTime()){ |
251 |
calib.iev = ce->iev; |
252 |
if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { |
253 |
for ( Int_t d=0 ; d<11 ;d++ ){ |
254 |
Int_t pre = -1; |
255 |
for ( Int_t j=0; j<96 ;j++){ |
256 |
if ( j%16 == 0 ) pre++; |
257 |
if ( s == 2 ){ |
258 |
calib.calped[0][2*d+1][j] = ce->calped[3][d][j]; |
259 |
calib.cstwerr[3] = ce->cstwerr[3]; |
260 |
calib.cperror[3] = ce->cperror[3]; |
261 |
calib.calgood[0][2*d+1][j] = ce->calgood[3][d][j]; |
262 |
calib.calthr[0][2*d+1][pre] = ce->calthr[3][d][pre]; |
263 |
calib.calrms[0][2*d+1][j] = ce->calrms[3][d][j]; |
264 |
calib.calbase[0][2*d+1][pre] = ce->calbase[3][d][pre]; |
265 |
calib.calvar[0][2*d+1][pre] = ce->calvar[3][d][pre]; |
266 |
}; |
267 |
if ( s == 3 ){ |
268 |
calib.calped[0][2*d][j] = ce->calped[1][d][j]; |
269 |
calib.cstwerr[1] = ce->cstwerr[1]; |
270 |
calib.cperror[1] = ce->cperror[1]; |
271 |
calib.calgood[0][2*d][j] = ce->calgood[1][d][j]; |
272 |
calib.calthr[0][2*d][pre] = ce->calthr[1][d][pre]; |
273 |
calib.calrms[0][2*d][j] = ce->calrms[1][d][j]; |
274 |
calib.calbase[0][2*d][pre] = ce->calbase[1][d][pre]; |
275 |
calib.calvar[0][2*d][pre] = ce->calvar[1][d][pre]; |
276 |
}; |
277 |
if ( s == 0 ){ |
278 |
calib.calped[1][2*d][j] = ce->calped[0][d][j]; |
279 |
calib.cstwerr[0] = ce->cstwerr[0]; |
280 |
calib.cperror[0] = ce->cperror[0]; |
281 |
calib.calgood[1][2*d][j] = ce->calgood[0][d][j]; |
282 |
calib.calthr[1][2*d][pre] = ce->calthr[0][d][pre]; |
283 |
calib.calrms[1][2*d][j] = ce->calrms[0][d][j]; |
284 |
calib.calbase[1][2*d][pre] = ce->calbase[0][d][pre]; |
285 |
calib.calvar[1][2*d][pre] = ce->calvar[0][d][pre]; |
286 |
}; |
287 |
if ( s == 1 ){ |
288 |
calib.calped[1][2*d+1][j] = ce->calped[2][d][j]; |
289 |
calib.cstwerr[2] = ce->cstwerr[2]; |
290 |
calib.cperror[2] = ce->cperror[2]; |
291 |
calib.calgood[1][2*d+1][j] = ce->calgood[2][d][j]; |
292 |
calib.calthr[1][2*d+1][pre] = ce->calthr[2][d][pre]; |
293 |
calib.calrms[1][2*d+1][j] = ce->calrms[2][d][j]; |
294 |
calib.calbase[1][2*d+1][pre] = ce->calbase[2][d][pre]; |
295 |
calib.calvar[1][2*d+1][pre] = ce->calvar[2][d][pre]; |
296 |
}; |
297 |
}; |
298 |
}; |
299 |
}; |
300 |
}; |
301 |
}; |
302 |
tr->Delete(); |
303 |
return(0); |
304 |
} |
305 |
|
306 |
TString getFilename(const TString filename){ |
307 |
const string fil = (const char*)filename; |
308 |
Int_t posiz = fil.find("dw_"); |
309 |
if ( posiz == -1 ) posiz = fil.find("DW_"); |
310 |
if ( posiz == -1 ) return 0; |
311 |
Int_t posiz2 = posiz+13; |
312 |
TString file2; |
313 |
stringcopy(file2,filename,posiz,posiz2); |
314 |
TString pdat(".dat"); |
315 |
stringappend(file2,pdat); |
316 |
return file2; |
317 |
} |
318 |
|
319 |
Int_t fetchpreviousfile(TString ffile, Int_t s){ |
320 |
const char *file = ffile; |
321 |
stringstream fil; |
322 |
fil.str(""); |
323 |
fil << file ; |
324 |
Int_t posiz = fil.str().find("dw_"); |
325 |
Int_t upper = 0; |
326 |
if ( posiz == -1 ) { |
327 |
posiz = fil.str().find("DW_"); |
328 |
upper = 1; |
329 |
}; |
330 |
if ( posiz == -1 ) return(1); |
331 |
// |
332 |
TString nomefile = getFilename(ffile); |
333 |
const char *ufnome = nomefile; |
334 |
stringstream nyear; |
335 |
stringstream nmonth; |
336 |
stringstream nday; |
337 |
stringstream nfno; |
338 |
TString year; |
339 |
TString month; |
340 |
TString day; |
341 |
TString fno; |
342 |
TString yno; |
343 |
stringcopy(year, ufnome, 3, 5); |
344 |
stringcopy(month, ufnome, 5, 7); |
345 |
stringcopy(day, ufnome, 7, 9); |
346 |
stringcopy(fno, ufnome, 10, 13); |
347 |
stringcopy(yno, ffile.Data(), posiz+14, posiz+16); |
348 |
const char *cyear = year; |
349 |
const char *cmonth = month; |
350 |
const char *cday = day; |
351 |
const char *cfno = fno; |
352 |
const char *cyno = yno; |
353 |
Int_t jy = atoi(cyear); |
354 |
Int_t jm = atoi(cmonth); |
355 |
Int_t jd = atoi(cday); |
356 |
Int_t jn = atoi(cfno); |
357 |
Int_t jynst = atoi(cyno); |
358 |
Int_t jyn = 0; |
359 |
Int_t inter = 0; |
360 |
char *zy; |
361 |
char *zm; |
362 |
char *zd; |
363 |
char *zn = ""; |
364 |
char *zyn; |
365 |
Int_t goodthis = 0; |
366 |
Int_t onegood = 0; |
367 |
Long64_t ncalibs; |
368 |
if ( jn > 1 ) { |
369 |
jn--; |
370 |
} else { |
371 |
jn = 999; |
372 |
if ( jd > 1 ) { |
373 |
jd--; |
374 |
} else { |
375 |
if ( jm > 1 ) { |
376 |
jm--; |
377 |
} else { |
378 |
jm = 12; |
379 |
if ( jy > 3 ) { |
380 |
jy--; |
381 |
} else { |
382 |
return(1); |
383 |
}; |
384 |
}; |
385 |
if ( jm==4 || jm==6 || jm==9 || jm==11 ){ |
386 |
jd = 30; |
387 |
} else { |
388 |
jd = 31; |
389 |
}; |
390 |
if ( jm==2 ) jd = 29; |
391 |
}; |
392 |
}; |
393 |
Int_t gday = 0; |
394 |
Int_t retvar = 1; |
395 |
stringstream nyd; |
396 |
while ( jy>3 ) { |
397 |
while ( jm>0 ){ |
398 |
while ( jd > 0 ){ |
399 |
while ( jn > 0 ){ |
400 |
retvar--; |
401 |
nyear.str(""); |
402 |
nyear << jy; |
403 |
if ( jy<10 ){ |
404 |
zy = "0"; |
405 |
} else { |
406 |
zy = ""; |
407 |
}; |
408 |
nmonth.str(""); |
409 |
nmonth << jm; |
410 |
if ( jm<10 ){ |
411 |
zm = "0"; |
412 |
} else { |
413 |
zm = ""; |
414 |
}; |
415 |
nday.str(""); |
416 |
nday << jd; |
417 |
if ( jd<10 ){ |
418 |
zd = "0"; |
419 |
} else { |
420 |
zd = ""; |
421 |
}; |
422 |
nfno.str(""); |
423 |
nfno << jn; |
424 |
if ( jn>99 ) zn = ""; |
425 |
if ( jn<100 && jn>9 ) zn = "0"; |
426 |
if ( jn<10 ) zn = "00"; |
427 |
jyn = jynst; |
428 |
goodthis = 0; |
429 |
onegood = 0; |
430 |
while ( jyn < 100 ){ |
431 |
begwhile: |
432 |
nyd.str(""); |
433 |
nyd << jyn; |
434 |
if ( jyn<10 ){ |
435 |
zyn = "0"; |
436 |
} else { |
437 |
zyn = ""; |
438 |
}; |
439 |
stringstream nfile; |
440 |
nfile.str(""); |
441 |
if ( upper ){ |
442 |
nfile << "DW_" << zy; |
443 |
} else { |
444 |
nfile << "dw_" << zy; |
445 |
}; |
446 |
nfile << nyear.str().c_str(); |
447 |
nfile << zm; |
448 |
nfile << nmonth.str().c_str(); |
449 |
nfile << zd; |
450 |
nfile << nday.str().c_str() << "_"; |
451 |
nfile << zn; |
452 |
nfile << nfno.str().c_str(); |
453 |
nfile << zyn; |
454 |
nfile << nyd.str().c_str(); |
455 |
nfile << ".root"; // addded |
456 |
// |
457 |
TString fgfile = ""; |
458 |
stringcopy(fgfile,ffile,0,posiz); |
459 |
stringappend(fgfile,nfile.str().c_str()); |
460 |
// |
461 |
TFile *File = 0; |
462 |
TTree *tr = 0; |
463 |
pamela::CalibCalPedEvent *ce = 0; |
464 |
pamela::PscuHeader *cph = 0; |
465 |
pamela::EventHeader *ceh = 0; |
466 |
if ( !existfile(fgfile) ){ |
467 |
if ( jyn > jynst && onegood ) { |
468 |
goodthis = 1; |
469 |
jyn--; |
470 |
goto begwhile; |
471 |
}; |
472 |
if ( jyn == jynst ) { |
473 |
onegood = 0; |
474 |
jyn = 100; |
475 |
goto salta; |
476 |
}; |
477 |
if ( jyn > jynst && !onegood ) { |
478 |
onegood = 0; |
479 |
jyn = 100; |
480 |
goto salta; |
481 |
}; |
482 |
goto salta; |
483 |
} else { |
484 |
onegood = 1; |
485 |
if ( !goodthis ) { |
486 |
goto salta; |
487 |
}; |
488 |
}; |
489 |
File = new TFile(fgfile.Data()); |
490 |
tr = (TTree*)File->Get("CalibCalPed"); |
491 |
//Takes the tree of the header file |
492 |
tr->SetBranchAddress("Header", &ceh); |
493 |
tr->SetBranchAddress("CalibCalPed", &ce); |
494 |
ncalibs = tr->GetEntries(); |
495 |
if ( ncalibs == 0 ) { |
496 |
jyn = 100; |
497 |
onegood = 0; |
498 |
goodthis = 0; |
499 |
File->Close(); |
500 |
goto salta; |
501 |
}; |
502 |
inter = 0; |
503 |
for (Int_t c = 0; c < ncalibs; c++){ |
504 |
tr->GetEntry(c); |
505 |
cph = ceh->GetPscuHeader(); |
506 |
if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { |
507 |
inter++; |
508 |
}; |
509 |
}; |
510 |
if ( inter ){ |
511 |
jyn = 100; |
512 |
File->Close(); |
513 |
return(retvar); |
514 |
}; |
515 |
File->Close(); |
516 |
salta: |
517 |
jyn++; |
518 |
}; |
519 |
jn--; |
520 |
}; |
521 |
jn = 999; |
522 |
jd--; |
523 |
gday++; |
524 |
if ( gday>1 ) return(1); |
525 |
}; |
526 |
jm--; |
527 |
if ( jm==4 || jm==6 || jm==9 || jm==11 ){ |
528 |
jd = 30; |
529 |
} else { |
530 |
jd = 31; |
531 |
}; |
532 |
if ( jm==2 ) jd = 29; |
533 |
}; |
534 |
jm = 12; |
535 |
jy--; |
536 |
}; |
537 |
return(retvar); |
538 |
} |
539 |
|
540 |
TString whatnamewith(TString ffile, Int_t retvar){ |
541 |
Int_t debug = 0; |
542 |
const char *file = ffile; |
543 |
stringstream fil; |
544 |
fil.str(""); |
545 |
fil << file ; |
546 |
if ( debug ) printf("ffile %s retvar %i \n",file,retvar); |
547 |
Int_t posiz = fil.str().find("dw_"); |
548 |
Int_t upper = 0; |
549 |
if ( posiz == -1 ) { |
550 |
posiz = fil.str().find("DW_"); |
551 |
upper = 1; |
552 |
}; |
553 |
if ( posiz == -1 ) return("wrong"); |
554 |
// |
555 |
TString nomefile = getFilename(ffile); |
556 |
const char *ufnome = nomefile; |
557 |
stringstream nyear; |
558 |
stringstream nmonth; |
559 |
stringstream nday; |
560 |
stringstream nfno; |
561 |
TString year; |
562 |
TString month; |
563 |
TString day; |
564 |
TString fno; |
565 |
TString yno; |
566 |
stringcopy(year, ufnome, 3, 5); |
567 |
stringcopy(month, ufnome, 5, 7); |
568 |
stringcopy(day, ufnome, 7, 9); |
569 |
stringcopy(fno, ufnome, 10, 13); |
570 |
stringcopy(yno, ffile.Data(), posiz+14, posiz+16); |
571 |
const char *cyear = year; |
572 |
const char *cmonth = month; |
573 |
const char *cday = day; |
574 |
const char *cfno = fno; |
575 |
const char *cyno = yno; |
576 |
Int_t jy = atoi(cyear); |
577 |
Int_t jm = atoi(cmonth); |
578 |
Int_t jd = atoi(cday); |
579 |
Int_t jn = atoi(cfno); |
580 |
Int_t jynst = atoi(cyno); |
581 |
Int_t jyn = 0; |
582 |
Int_t goodthis = 0; |
583 |
Int_t onegood = 0; |
584 |
char *zy; |
585 |
char *zm; |
586 |
char *zd; |
587 |
char *zn = ""; |
588 |
char *zyn; |
589 |
Long64_t ncalibs; |
590 |
TString fgfile; |
591 |
if ( jn > 1 ) { |
592 |
jn--; |
593 |
} else { |
594 |
jn = 999; |
595 |
if ( jd > 1 ) { |
596 |
jd--; |
597 |
} else { |
598 |
if ( jm > 1 ) { |
599 |
jm--; |
600 |
} else { |
601 |
jm = 12; |
602 |
if ( jy > 3 ) { |
603 |
jy--; |
604 |
} else { |
605 |
return("wrong"); |
606 |
}; |
607 |
}; |
608 |
if ( jm==4 || jm==6 || jm==9 || jm==11 ){ |
609 |
jd = 30; |
610 |
} else { |
611 |
jd = 31; |
612 |
}; |
613 |
if ( jm==2 ) jd = 29; |
614 |
}; |
615 |
}; |
616 |
Int_t retval = 1; |
617 |
Int_t gday = 0; |
618 |
stringstream nyd; |
619 |
while ( jy>3 ) { |
620 |
while ( jm>0 ){ |
621 |
while ( jd > 0 ){ |
622 |
while ( jn > 0 ){ |
623 |
retval--; |
624 |
if ( retval < retvar ) return("wrong"); |
625 |
nyear.str(""); |
626 |
nyear << jy; |
627 |
if ( jy<10 ){ |
628 |
zy = "0"; |
629 |
} else { |
630 |
zy = ""; |
631 |
}; |
632 |
nmonth.str(""); |
633 |
nmonth << jm; |
634 |
if ( jm<10 ){ |
635 |
zm = "0"; |
636 |
} else { |
637 |
zm = ""; |
638 |
}; |
639 |
nday.str(""); |
640 |
nday << jd; |
641 |
if ( jd<10 ){ |
642 |
zd = "0"; |
643 |
} else { |
644 |
zd = ""; |
645 |
}; |
646 |
nfno.str(""); |
647 |
nfno << jn; |
648 |
if ( jn>99 ) zn = ""; |
649 |
if ( jn<100 && jn>9 ) zn = "0"; |
650 |
if ( jn<10 ) zn = "00"; |
651 |
jyn = jynst; |
652 |
goodthis = 0; |
653 |
onegood = 0; |
654 |
while ( jyn < 100 ){ |
655 |
begwhile: |
656 |
nyd.str(""); |
657 |
nyd << jyn; |
658 |
if ( jyn<10 ){ |
659 |
zyn = "0"; |
660 |
} else { |
661 |
zyn = ""; |
662 |
}; |
663 |
stringstream nfile; |
664 |
nfile.str(""); |
665 |
if ( upper ){ |
666 |
nfile << "DW_" << zy; |
667 |
} else { |
668 |
nfile << "dw_" << zy; |
669 |
}; |
670 |
nfile << nyear.str().c_str(); |
671 |
nfile << zm; |
672 |
nfile << nmonth.str().c_str(); |
673 |
nfile << zd; |
674 |
nfile << nday.str().c_str() << "_"; |
675 |
nfile << zn; |
676 |
nfile << nfno.str().c_str(); |
677 |
nfile << zyn; |
678 |
nfile << nyd.str().c_str(); |
679 |
nfile << ".root"; // added |
680 |
// |
681 |
fgfile = ""; |
682 |
stringcopy(fgfile,ffile,0,posiz); |
683 |
stringappend(fgfile,nfile.str().c_str()); |
684 |
const char *tosee = fgfile; |
685 |
if ( debug ) printf("file: %s retval %i retvar %i \n",tosee,retval,retvar); |
686 |
TFile *File = 0; |
687 |
TTree *tr = 0; |
688 |
// |
689 |
if ( !existfile(fgfile) ){ |
690 |
if ( jyn > jynst && onegood ) { |
691 |
goodthis = 1; |
692 |
jyn--; |
693 |
goto begwhile; |
694 |
}; |
695 |
if ( jyn == jynst ) { |
696 |
onegood = 0; |
697 |
jyn = 100; |
698 |
goto salta; |
699 |
}; |
700 |
if ( jyn > jynst && !onegood ) { |
701 |
onegood = 0; |
702 |
jyn = 100; |
703 |
goto salta; |
704 |
}; |
705 |
goto salta; |
706 |
} else { |
707 |
onegood = 1; |
708 |
if ( !goodthis ){ |
709 |
goto salta; |
710 |
}; |
711 |
}; |
712 |
File = new TFile(fgfile.Data()); |
713 |
tr = (TTree*)File->Get("CalibCalPed"); |
714 |
ncalibs = tr->GetEntries(); |
715 |
if ( ncalibs == 0 ) { |
716 |
jyn = 100; |
717 |
onegood = 0; |
718 |
goodthis = 0; |
719 |
File->Close(); |
720 |
goto salta; |
721 |
}; |
722 |
File->Close(); |
723 |
if ( retval == retvar ) return(fgfile); |
724 |
salta: |
725 |
jyn++; |
726 |
}; |
727 |
jn--; |
728 |
}; |
729 |
jn = 999; |
730 |
jd--; |
731 |
gday++; |
732 |
if ( gday>7 ) return("wrong"); |
733 |
}; |
734 |
jm--; |
735 |
if ( jm==4 || jm==6 || jm==9 || jm==11 ){ |
736 |
jd = 30; |
737 |
} else { |
738 |
jd = 31; |
739 |
}; |
740 |
if ( jm==2 ) jd = 29; |
741 |
}; |
742 |
jm = 12; |
743 |
jy--; |
744 |
}; |
745 |
return("wrong"); |
746 |
} |
747 |
|
748 |
char *getLEVname(TString filename, TString detc, TString numb, TString frame="root"){ |
749 |
// char *file; |
750 |
stringstream file; |
751 |
file.str(""); |
752 |
const char *num = numb; |
753 |
const char *det = detc; |
754 |
const char *fra = frame; |
755 |
string fil = (const char *)filename; |
756 |
Int_t posiz = fil.find("dw_"); |
757 |
if ( posiz == -1 ) posiz = fil.find("DW_"); |
758 |
if ( posiz == -1 ) return("unkwnown.root"); |
759 |
Int_t spos = posiz; |
760 |
Int_t epos = posiz+15; |
761 |
TString tmptempf; |
762 |
stringcopy(tmptempf,filename,spos,epos); |
763 |
const char *tempf = tmptempf; |
764 |
file << tempf << ".Physics.Level"; |
765 |
file << num << "."; |
766 |
file << det << ".Event."; |
767 |
file << fra; |
768 |
const char *rfile = file.str().c_str(); |
769 |
return((char*)rfile); |
770 |
}; |
771 |
|
772 |
char *getfileLEVname(TString filename, TString detc, TString numb, TString frame="root"){ |
773 |
// char *file; |
774 |
stringstream file; |
775 |
const char *num = numb; |
776 |
const char *det = detc; |
777 |
const char *fra = frame; |
778 |
string fil = (const char *)filename; |
779 |
Int_t posiz = fil.find("dw_"); |
780 |
if ( posiz == -1 ) posiz = fil.find("DW_"); |
781 |
if ( posiz == -1 ) return("unknown.root"); |
782 |
Int_t spos = posiz; |
783 |
Int_t epos = posiz+13; |
784 |
TString tmptempf; |
785 |
stringcopy(tmptempf,filename,spos,epos); |
786 |
const char *tempf = tmptempf; |
787 |
file.str(""); |
788 |
file << tempf << "00.Physics.Level"; |
789 |
// file << "00.Physics.Level"; |
790 |
file << num << "."; |
791 |
file << det << ".Event."; |
792 |
file << fra; |
793 |
const char *rfile = file.str().c_str(); |
794 |
return((char*)rfile); |
795 |
}; |
796 |
|
797 |
int OLDCaloFindCalibs(TString &filename, Calib & calib){ |
798 |
for (Int_t s = 0; s < 4; s++){ |
799 |
for (Int_t d = 1; d<50; d++){ |
800 |
calib.ttime[s][d] = 0; |
801 |
if ( d < 49 ) calib.time[s][d] = 0; |
802 |
}; |
803 |
}; |
804 |
if ( !existfile(filename) ) return(1); |
805 |
// |
806 |
TFile *File = new TFile(filename.Data()); |
807 |
TTree *tr = (TTree*)File->Get("CalibCalPed"); |
808 |
pamela::CalibCalPedEvent *ce = 0; |
809 |
pamela::PscuHeader *cph = 0; |
810 |
pamela::EventHeader *ceh = 0; |
811 |
tr->SetBranchAddress("Header", &ceh); |
812 |
tr->SetBranchAddress("CalibCalPed", &ce); |
813 |
Long64_t ncalibs = tr->GetEntries(); |
814 |
Int_t inter; |
815 |
for (Int_t s = 0; s < 4; s++){ |
816 |
for (Int_t d = 1; d<50; d++){ |
817 |
calib.ttime[s][d] = 0; |
818 |
if ( d < 49 ) calib.time[s][d] = 0; |
819 |
}; |
820 |
}; |
821 |
for (Int_t s = 0; s < 4; s++){ |
822 |
inter = 0; |
823 |
for (Int_t c = 0; c < ncalibs; c++){ |
824 |
tr->GetEntry(c); |
825 |
cph = ceh->GetPscuHeader(); |
826 |
calib.ttime[s][inter] = 0; |
827 |
if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { |
828 |
calib.ttime[s][inter] = cph->GetOrbitalTime(); |
829 |
inter++; |
830 |
} else { |
831 |
if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) { |
832 |
printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]); |
833 |
}; |
834 |
}; |
835 |
}; |
836 |
if ( inter == 0 ){ |
837 |
printf(" ERROR: no suitable calibration for section %i in this file!\n",s); |
838 |
}; |
839 |
for (Int_t d = 1; d<50; d++){ |
840 |
if ( calib.ttime[s][d] != 0 ) { |
841 |
calib.time[s][d-1] = calib.ttime[s][d-1] + (int)((calib.ttime[s][d] - calib.ttime[s][d-1])/2.); |
842 |
} else { |
843 |
if ( d == 1 ) { |
844 |
calib.time[s][d-1] = 0; |
845 |
}; |
846 |
}; |
847 |
}; |
848 |
}; |
849 |
File->Close(); |
850 |
return(0); |
851 |
} |
852 |
|
853 |
int CaloFindCalibs(TString &filename, TString &calibfile, Int_t &wused, Calib & calib){ |
854 |
// |
855 |
// First of all check if the file has a monotone growing OBT time function. |
856 |
// |
857 |
Int_t debug = 0; |
858 |
Int_t obtjump = 0; |
859 |
Int_t evjump[50]; |
860 |
Int_t OBT = 0; |
861 |
Int_t OOBT = 0; |
862 |
Int_t OBT1 = 0; |
863 |
Int_t LOBT = 0; |
864 |
pamela::PscuHeader *ph = 0; |
865 |
pamela::EventHeader *eh = 0; |
866 |
// |
867 |
if ( !existfile(filename) ) return(1); |
868 |
TFile *File = new TFile(filename.Data()); |
869 |
TTree *ctr = (TTree*)File->Get("Physics"); |
870 |
ctr->SetBranchStatus("*",0); |
871 |
ctr->SetBranchStatus("Header*",1); |
872 |
ctr->SetBranchStatus("Pscu*",1); |
873 |
Long64_t nevents = ctr->GetEntries(); |
874 |
if ( nevents == 0 ) { |
875 |
printf(" The file is empty! \n "); |
876 |
return(1); |
877 |
}; |
878 |
ctr->SetBranchAddress("Header", &eh); |
879 |
for (Int_t i = 0; i < nevents; i++){ |
880 |
ctr->GetEntry(i); |
881 |
ph = eh->GetPscuHeader(); |
882 |
OBT = ph->GetOrbitalTime(); |
883 |
if ( !i ) OBT1 = OBT; |
884 |
if ( i == nevents-1 ) LOBT = OBT; |
885 |
if ( OOBT > OBT ) { |
886 |
evjump[obtjump] = i; |
887 |
obtjump++; |
888 |
}; |
889 |
OOBT = OBT; |
890 |
}; |
891 |
printf("\n Obtjump = %i - FIRST OBT %i - LAST OBT %i \n\n",obtjump,OBT1,LOBT); |
892 |
// |
893 |
TTree *tr; |
894 |
pamela::CalibCalPedEvent *ce = 0; |
895 |
pamela::PscuHeader *cph = 0; |
896 |
pamela::EventHeader *ceh = 0; |
897 |
Long64_t ncalibs; |
898 |
// |
899 |
for (Int_t s = 0; s < 4; s++){ |
900 |
for (Int_t d = 1; d<50; d++){ |
901 |
calib.ttime[s][d] = 0; |
902 |
if ( d < 49 ) calib.time[s][d] = 0; |
903 |
}; |
904 |
}; |
905 |
for (Int_t s = 0; s < 4; s++){ |
906 |
// |
907 |
Int_t firstlook = 0; |
908 |
repeatsearch: |
909 |
// |
910 |
Int_t inter = 0; |
911 |
// |
912 |
tr = (TTree*)File->Get("CalibCalPed"); |
913 |
tr->SetBranchAddress("Header", &ceh); |
914 |
tr->SetBranchAddress("CalibCalPed", &ce); |
915 |
if ( tr->GetEntries() == 0 ){ |
916 |
printf(" No calibration in this file! \n"); |
917 |
if ( !firstlook ){ |
918 |
wused = 2; |
919 |
filename = calibfile; |
920 |
firstlook = 1; |
921 |
goto repeatsearch; |
922 |
} else { |
923 |
wused = 0; |
924 |
goto jte; |
925 |
}; |
926 |
}; |
927 |
wused = 1; |
928 |
// |
929 |
ncalibs = tr->GetEntries(); |
930 |
calib.obtjump = 0; |
931 |
inter = 0; |
932 |
for (Int_t c = 0; c < ncalibs; c++){ |
933 |
tr->GetEntry(c); |
934 |
cph = ceh->GetPscuHeader(); |
935 |
calib.ttime[s][inter] = 0; |
936 |
if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { |
937 |
calib.ttime[s][inter] = cph->GetOrbitalTime(); |
938 |
inter++; |
939 |
} else { |
940 |
if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) { |
941 |
printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]); |
942 |
}; |
943 |
}; |
944 |
}; |
945 |
jte: |
946 |
if ( inter == 0 ){ |
947 |
printf(" WARNING: no suitable calibration for section %i in this file!\n",s); |
948 |
printf(" I WILL SEARCH IN PREVIOUS FILES\n"); |
949 |
if ( !firstlook ){ |
950 |
wused = 2; |
951 |
filename = calibfile; |
952 |
firstlook = 1; |
953 |
goto repeatsearch; |
954 |
}; |
955 |
}; |
956 |
if ( inter > 50 ){ |
957 |
printf(" WARNING: cannot handle more than 50 calibrations for file!\n"); |
958 |
printf(" I WILL SEARCH IN PREVIOUS FILES\n"); |
959 |
inter = 0; |
960 |
if ( !firstlook ){ |
961 |
wused = 2; |
962 |
filename = calibfile; |
963 |
firstlook = 1; |
964 |
goto repeatsearch; |
965 |
}; |
966 |
}; |
967 |
if ( obtjump ){ |
968 |
calib.obtjump = 1; |
969 |
printf(" ERROR: jump in OBT! CPU restarded? \n"); |
970 |
printf(" WARNING: I will use only the first calibration (if any)\n"); |
971 |
}; |
972 |
// |
973 |
Int_t lastcal = 0; |
974 |
Int_t of = 0; |
975 |
if ( OBT1 < calib.ttime[s][0] || !inter ) { |
976 |
// |
977 |
// here we must look for a calibration in the previous file... |
978 |
// |
979 |
TString ffile = ""; |
980 |
stringcopy(ffile,filename,0,0); |
981 |
Int_t fcodep = fetchpreviousfile(ffile,s); |
982 |
if ( fcodep == 1 ) { |
983 |
if ( inter ){ |
984 |
of = -3; |
985 |
printf(" WARNING: Problems to find a good calibration for section %i \n",s); |
986 |
} else { |
987 |
of = -2; |
988 |
printf(" ERROR: I was not able to find any good calibration for section %i \n",s); |
989 |
}; |
990 |
} else { |
991 |
stringcopy(ffile,filename,0,0); |
992 |
TString pfile = whatnamewith(ffile,fcodep); |
993 |
struct Calib tcal; |
994 |
OLDCaloFindCalibs(pfile,tcal); |
995 |
for ( Int_t d = 0; d<51; d++ ){ |
996 |
if ( tcal.ttime[s][d] == 0 ){ |
997 |
lastcal = d-1; |
998 |
break; |
999 |
}; |
1000 |
}; |
1001 |
if ( !inter ){ |
1002 |
calib.ttime[s][0] = tcal.ttime[s][lastcal]; |
1003 |
calib.time[s][0] = OBT1; |
1004 |
calib.time[s][1] = LOBT; |
1005 |
calib.fcode[s][0] = fcodep; |
1006 |
if ( debug ) printf("boh! %i \n",fcodep); |
1007 |
if ( debug ) getchar(); |
1008 |
of = -1; |
1009 |
} else { |
1010 |
if ( obtjump ){ |
1011 |
if ( debug ) printf("boh2! %i \n",fcodep); |
1012 |
if ( debug ) getchar(); |
1013 |
inter = 1; |
1014 |
calib.time[s][0] = OBT1; |
1015 |
calib.time[s][1] = LOBT; |
1016 |
calib.fcode[s][0] = 10; |
1017 |
of = -5; |
1018 |
for (Int_t d = 1; d<50; d++){ |
1019 |
calib.ttime[s][d] = 0; |
1020 |
if ( d>1 && d<49 ) calib.time[s][d] = 0; |
1021 |
}; |
1022 |
} else { |
1023 |
if ( lastcal >= 0 && tcal.ttime[s][lastcal]<OBT1 && (OBT1-tcal.ttime[s][lastcal])<7200000 && (OBT1-tcal.ttime[s][lastcal])>0 ){ |
1024 |
for (Int_t i = inter+1; i>0; i--){ |
1025 |
calib.ttime[s][i] = calib.ttime[s][i-1]; |
1026 |
}; |
1027 |
calib.time[s][0] = OBT1; |
1028 |
calib.ttime[s][0] = tcal.ttime[s][lastcal]; |
1029 |
calib.fcode[s][0] = fcodep; |
1030 |
of = 2; |
1031 |
} else { |
1032 |
calib.time[s][0] = OBT1; |
1033 |
calib.fcode[s][0] = 10; |
1034 |
of = 1; |
1035 |
}; |
1036 |
}; |
1037 |
}; |
1038 |
}; |
1039 |
} else { |
1040 |
of = 0; |
1041 |
}; |
1042 |
if ( debug ) printf("boh3! %i \n",of); |
1043 |
Int_t of1 = 0; |
1044 |
if ( inter && of != 0 && of != 2 && of != -2 && of != -5) { |
1045 |
for (Int_t i = inter+1; i>0; i--){ |
1046 |
calib.ttime[s][i] = calib.ttime[s][i-1]; |
1047 |
}; |
1048 |
of1 = 0; |
1049 |
}; |
1050 |
if ( of == -2 ) of = -1; |
1051 |
if ( of == -3 ) { |
1052 |
of = 1; |
1053 |
calib.fcode[s][0] = 10; |
1054 |
}; |
1055 |
// |
1056 |
if ( of != -1 && of != -5 ){ |
1057 |
if ( of == 2 ) of = 1; |
1058 |
for (Int_t d = 0; d<inter+1; d++){ |
1059 |
calib.time[s][d+of] = calib.ttime[s][d+of]; |
1060 |
calib.fcode[s][d+of] = 10; |
1061 |
if ( d == inter ){ |
1062 |
if ( debug ) printf("boh4! %i \n",of); |
1063 |
calib.time[s][d+of1+of] = LOBT; |
1064 |
calib.fcode[s][d+of] = 0; |
1065 |
}; |
1066 |
}; |
1067 |
}; |
1068 |
// if ( calheadFile ) calheadFile->Close(); |
1069 |
// if ( calcalibFile ) calcalibFile->Close(); |
1070 |
//if ( File ) File->Close(); |
1071 |
}; |
1072 |
return(0); |
1073 |
} |
1074 |
|
1075 |
void Calo1stcalib(TString filename, Calib & calib, Int_t b[4]){ |
1076 |
Float_t estrip[2][22][96]; |
1077 |
// |
1078 |
// this is the value of the mip for each strip. To be changed when we will have the real values |
1079 |
// |
1080 |
for (Int_t s=0; s<4;s++){ |
1081 |
for (Int_t d = 0; d<50; d++){ |
1082 |
calib.ttime[s][d] = 0 ; |
1083 |
if ( d < 49 ) calib.time[s][d] = 0 ; |
1084 |
}; |
1085 |
}; |
1086 |
for (Int_t m = 0; m < 2 ; m++ ){ |
1087 |
for (Int_t k = 0; k < 22; k++ ){ |
1088 |
for (Int_t l = 0; l < 96; l++ ){ |
1089 |
calib.calped[m][k][l] = 0. ; |
1090 |
estrip[m][k][l] = 0.; |
1091 |
}; |
1092 |
}; |
1093 |
} |
1094 |
// |
1095 |
// first of all find the calibrations in the file |
1096 |
// |
1097 |
OLDCaloFindCalibs(filename, calib); |
1098 |
// |
1099 |
// print on the screen the results: |
1100 |
// |
1101 |
printf(" ---------------------------------------------------------- \n \n"); |
1102 |
Int_t calibex = 0; |
1103 |
for (Int_t s=0; s<4;s++){ |
1104 |
Int_t stop = 0; |
1105 |
for (Int_t d = 0; d<48; d++){ |
1106 |
if ( calib.ttime[s][d] != 0 ) calibex++; |
1107 |
if ( calib.time[s][0] != 0 ){ |
1108 |
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]); |
1109 |
if ( calib.time[s][d+1] != 0 ) { |
1110 |
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]); |
1111 |
} else { |
1112 |
if ( !stop ){ |
1113 |
printf(" Section %i from time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d+1]); |
1114 |
stop = 1; |
1115 |
}; |
1116 |
}; |
1117 |
} else { |
1118 |
if ( calib.ttime[s][d] != 0 ) printf(" Section %i from time 0 use calibration at time %i \n",s,calib.ttime[s][d]); |
1119 |
}; |
1120 |
}; |
1121 |
printf("\n"); |
1122 |
}; |
1123 |
printf(" ---------------------------------------------------------- \n"); |
1124 |
if ( calibex < 1 ) { |
1125 |
printf("No full calibration data in this file, sorry!\n"); |
1126 |
} else { |
1127 |
// |
1128 |
// calibrate before starting |
1129 |
// |
1130 |
for (Int_t s = 0; s < 4; s++){ |
1131 |
b[s]=0; |
1132 |
CaloPede(filename,s,calib.ttime[s][b[s]],calib); |
1133 |
}; |
1134 |
}; |
1135 |
} |
1136 |
|
1137 |
|
1138 |
void ColorMIP(Float_t mip, int& colo){ |
1139 |
// printf("mip = %f \n",mip); |
1140 |
if ( colo > 0 ){ |
1141 |
colo = 10; |
1142 |
if ( mip > 0.7 ) colo = 38; |
1143 |
if ( mip > 2. ) colo = 4; |
1144 |
if ( mip > 10. ) colo = 3; |
1145 |
if ( mip > 100. ) colo = 2; |
1146 |
if ( mip > 500. ) colo = 6; |
1147 |
} else { |
1148 |
colo = 10; |
1149 |
if ( mip > 0.7 ) colo = 17; |
1150 |
if ( mip > 2. ) colo = 15; |
1151 |
if ( mip > 10. ) colo = 14; |
1152 |
if ( mip > 100. ) colo = 13; |
1153 |
if ( mip > 500. ) colo = 12; |
1154 |
}; |
1155 |
} |
1156 |
|
1157 |
void ColorTOFMIP(Float_t mip, int& colo){ |
1158 |
// printf("mip = %f \n",mip); |
1159 |
if ( colo > 0 ){ |
1160 |
colo = 10; |
1161 |
if ( mip > 0. ) colo = 38; |
1162 |
if ( mip > 2. ) colo = 4; |
1163 |
if ( mip > 10. ) colo = 3; |
1164 |
if ( mip > 100. ) colo = 2; |
1165 |
if ( mip > 500. ) colo = 6; |
1166 |
} else { |
1167 |
colo = 10; |
1168 |
if ( mip > 0. ) colo = 17; |
1169 |
if ( mip > 2. ) colo = 15; |
1170 |
if ( mip > 10. ) colo = 14; |
1171 |
if ( mip > 100. ) colo = 13; |
1172 |
if ( mip > 500. ) colo = 12; |
1173 |
}; |
1174 |
} |
1175 |
|
1176 |
void *processevents(void *ptr){ |
1177 |
char *tellme = (char *)ptr; |
1178 |
cin.getline(tellme,256); |
1179 |
delete TThread::Self(); |
1180 |
return(NULL); |
1181 |
} |
1182 |
|
1183 |
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) { |
1184 |
char *bw; |
1185 |
if ( jumpto == -1 ){ |
1186 |
bw = "bw"; |
1187 |
} else { |
1188 |
bw = ""; |
1189 |
}; |
1190 |
jumpto = 0; |
1191 |
char tellme[256]; |
1192 |
stringstream input; |
1193 |
char tellme2[256]; |
1194 |
stringstream input2; |
1195 |
char tellme3[256]; |
1196 |
stringstream input3; |
1197 |
stringstream out; |
1198 |
out.str("x"); |
1199 |
input2.str("zzzzzzzzzzzzzz"); |
1200 |
input3.str("z"); |
1201 |
input << out.str().c_str(); |
1202 |
stringstream stc; |
1203 |
stringstream stc2; |
1204 |
while ( !strcmp(input.str().c_str(),out.str().c_str()) ) { |
1205 |
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"); |
1206 |
// |
1207 |
input.str(""); |
1208 |
if ( !gROOT->GetListOfClasses()->FindObject("TRint") ){ |
1209 |
TThread *thread = new TThread((TThread::VoidFunc_t)processevents,(void *)tellme); |
1210 |
thread->Run(); |
1211 |
while ( thread->GetState() == TThread::kRunningState ){ |
1212 |
gSystem->ProcessEvents(); |
1213 |
}; |
1214 |
thread->Kill(); |
1215 |
} else { |
1216 |
cin.getline(tellme,256); |
1217 |
}; |
1218 |
input << tellme; |
1219 |
// |
1220 |
out << "y"; |
1221 |
// |
1222 |
stc.str("b"); |
1223 |
// |
1224 |
if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { |
1225 |
if ( i > 0 ) { |
1226 |
printf("WARNING: going backward!\n\n"); |
1227 |
doflag = 2; |
1228 |
} else { |
1229 |
printf("This is the first event, you can't go backward! \n"); |
1230 |
out.str(""); |
1231 |
out << input.str().c_str(); |
1232 |
}; |
1233 |
}; |
1234 |
// |
1235 |
stc.str("j"); |
1236 |
if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { |
1237 |
printf("\n Do you want to jump to a progressive number [p] or to an event number [e]? "); |
1238 |
cin.getline(tellme3,256); |
1239 |
input3.str(""); |
1240 |
input3 << tellme3; |
1241 |
stc.str("p"); |
1242 |
if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) { |
1243 |
printf("\n Enter the progressive number you want to jump to: "); |
1244 |
cin.getline(tellme2,256); |
1245 |
input2.str(""); |
1246 |
input2 << tellme2; |
1247 |
Int_t j; |
1248 |
j = atoi(input2.str().c_str()); |
1249 |
if ( j < 1 || j > nevents+1 ) { |
1250 |
printf("\n You can choose between 1 and %i \n",(int)nevents+1); |
1251 |
out.str(""); |
1252 |
out << input.str().c_str(); |
1253 |
} else { |
1254 |
printf("\n Jumping to progressive number %i\n\n",j); |
1255 |
i = j-2; |
1256 |
}; |
1257 |
}; |
1258 |
stc.str("e"); |
1259 |
if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) { |
1260 |
printf("\n Enter the event number you want to jump to: "); |
1261 |
cin.getline(tellme2,256); |
1262 |
input2.str(""); |
1263 |
input2 << tellme2; |
1264 |
Int_t j; |
1265 |
j = atoi(input2.str().c_str()); |
1266 |
if ( j < firstevno || j > lastevno ) { |
1267 |
printf("\n You can choose between %i and %i \n",(int)firstevno,(int)lastevno); |
1268 |
out.str(""); |
1269 |
out << input.str().c_str(); |
1270 |
} else { |
1271 |
printf("\n Jumping to event number %i\n\n",j); |
1272 |
jumpto = j; |
1273 |
i = -1; |
1274 |
}; |
1275 |
}; |
1276 |
stc.str("e"); |
1277 |
stc2.str("p"); |
1278 |
if ( strcmp(input3.str().c_str(),stc.str().c_str()) && strcmp(input3.str().c_str(),stc2.str().c_str()) ) { |
1279 |
printf(" You must type or \"p\" or \"e\"\n"); |
1280 |
out.str(""); |
1281 |
out << input.str().c_str(); |
1282 |
}; |
1283 |
}; |
1284 |
// |
1285 |
stc.str("q"); |
1286 |
if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { |
1287 |
printf("Exiting...\n"); |
1288 |
return(1); |
1289 |
}; |
1290 |
// |
1291 |
stc.str("p"); |
1292 |
if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { |
1293 |
printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n"); |
1294 |
cin.getline(tellme2,256); |
1295 |
input2.str(""); |
1296 |
input2 << tellme2; |
1297 |
out.str(""); |
1298 |
out << input.str().c_str(); |
1299 |
// |
1300 |
TString filename = file; |
1301 |
const string fil = (const char*)filename; |
1302 |
Int_t posiz = fil.find("dw_"); |
1303 |
if ( posiz == -1 ) posiz = fil.find("DW_"); |
1304 |
Int_t posiz2 = posiz+13; |
1305 |
TString file2; |
1306 |
stringcopy(file2,filename,posiz,posiz2); |
1307 |
const char *figrec = file2; |
1308 |
const char *outdir = outDir; |
1309 |
stringstream figsave; |
1310 |
const char *ty = figty; |
1311 |
figsave.str(""); |
1312 |
figsave << outdir << "/"; |
1313 |
figsave << ty << "_"; |
1314 |
figsave << (i+1) << "_"; |
1315 |
figsave << figrec; |
1316 |
figsave << bw << "."; |
1317 |
figsave << input2.str().c_str(); |
1318 |
figure.SaveAs(figsave.str().c_str()); |
1319 |
printf("\n"); |
1320 |
}; |
1321 |
}; |
1322 |
return(0); |
1323 |
} |
1324 |
|
1325 |
void PrintFigure(TString filename, TString outDir, TString figty, TString saveas, TCanvas& figure) { |
1326 |
const string fil = (const char*)filename; |
1327 |
Int_t posiz = fil.find("dw_"); |
1328 |
if ( posiz == -1 ) posiz = fil.find("DW_"); |
1329 |
Int_t posiz2 = posiz+13; |
1330 |
TString file2; |
1331 |
stringcopy(file2,filename,posiz,posiz2); |
1332 |
// |
1333 |
const char *figrec = file2; |
1334 |
stringstream figsave; |
1335 |
figsave.str(""); |
1336 |
figsave << outDir.Data() << "/"; |
1337 |
figsave << figrec << "_"; |
1338 |
figsave << figty.Data() << "."; |
1339 |
figsave << saveas.Data(); |
1340 |
figure.SaveAs(figsave.str().c_str()); |
1341 |
return; |
1342 |
} |
1343 |
|
1344 |
Double_t langaufun(Double_t *x, Double_t *par) { |
1345 |
// Numeric constants |
1346 |
Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) |
1347 |
Double_t mpshift = -0.22278298; // Landau maximum location |
1348 |
// Control constants |
1349 |
Double_t np = 100.0; // number of convolution steps |
1350 |
Double_t sc = 3.; // convolution extends to +-sc Gaussian sigmas |
1351 |
// Variables |
1352 |
Double_t xx; |
1353 |
Double_t mpc; |
1354 |
Double_t fland; |
1355 |
Double_t sum = 0.0; |
1356 |
Double_t xlow,xupp; |
1357 |
Double_t step; |
1358 |
Double_t i; |
1359 |
// MP shift correction |
1360 |
mpc = par[1] - mpshift * par[0]; |
1361 |
// Range of convolution integral |
1362 |
xlow = x[0] - sc * par[3]; |
1363 |
xupp = x[0] + sc * par[3]; |
1364 |
step = (xupp-xlow) / np; |
1365 |
// Convolution integral of Landau and Gaussian by sum |
1366 |
for(i=1.0; i<=np/2; i++) { |
1367 |
xx = xlow + (i-.5) * step; |
1368 |
fland = TMath::Landau(xx,mpc,par[0]) / par[0]; |
1369 |
sum += fland * TMath::Gaus(x[0],xx,par[3]); |
1370 |
xx = xupp - (i-.5) * step; |
1371 |
fland = TMath::Landau(xx,mpc,par[0]) / par[0]; |
1372 |
sum += fland * TMath::Gaus(x[0],xx,par[3]); |
1373 |
}; |
1374 |
return (par[2] * step * sum * invsq2pi / par[3]); |
1375 |
} |
1376 |
|
1377 |
|
1378 |
|
1379 |
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"){ |
1380 |
// Once again, here are the Landau * Gaussian parameters: |
1381 |
// par[0]=Width (scale) parameter of Landau density |
1382 |
// par[1]=Most Probable (MP, location) parameter of Landau density |
1383 |
// par[2]=Total area (integral -inf to inf, normalization constant) |
1384 |
// par[3]=Width (sigma) of convoluted Gaussian function |
1385 |
// |
1386 |
// Variables for langaufit call: |
1387 |
// his histogram to fit |
1388 |
// fitrange[2] lo and hi boundaries of fit range |
1389 |
// startvalues[4] reasonable start values for the fit |
1390 |
// parlimitslo[4] lower parameter limits |
1391 |
// parlimitshi[4] upper parameter limits |
1392 |
// fitparams[4] returns the final fit parameters |
1393 |
// fiterrors[4] returns the final fit errors |
1394 |
// ChiSqr returns the chi square |
1395 |
// NDF returns ndf |
1396 |
Int_t i; |
1397 |
Char_t FunName[100]; |
1398 |
sprintf(FunName,"Fitfcn_%s",his->GetName()); |
1399 |
// |
1400 |
TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); |
1401 |
if (ffitold) delete ffitold; |
1402 |
// |
1403 |
TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); |
1404 |
ffit->SetParameters(startvalues); |
1405 |
ffit->SetParNames("Width","MP","Area","GSigma"); |
1406 |
// |
1407 |
for (i=0; i<4; i++) { |
1408 |
ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); |
1409 |
}; |
1410 |
his->Fit(FunName,drw); // fit within specified range, use ParLimits, do not plot |
1411 |
// |
1412 |
ffit->GetParameters(fitparams); // obtain fit parameters |
1413 |
for (i=0; i<4; i++) { |
1414 |
fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors |
1415 |
} |
1416 |
ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 |
1417 |
NDF[0] = ffit->GetNDF(); // obtain ndf |
1418 |
// |
1419 |
return (ffit); // return fit function |
1420 |
} |
1421 |
|
1422 |
|
1423 |
//Double_t langaupro(Double_t *params, Double_t &maxx) { |
1424 |
Double_t langaupro(Double_t *params) { |
1425 |
// Seaches for the location (x value) at the maximum of the |
1426 |
// Landau-Gaussian convolute and its full width at half-maximum. |
1427 |
// |
1428 |
Double_t p,x; |
1429 |
Double_t step; |
1430 |
Double_t l,lold; |
1431 |
Int_t i = 0; |
1432 |
Int_t MAXCALLS = 5000; |
1433 |
// Search for maximum |
1434 |
p = params[1] - 0.1 * params[0]; |
1435 |
step = 0.05 * params[0]; |
1436 |
lold = -2.0; |
1437 |
l = -1.0; |
1438 |
// while ( (l != lold) && (i < MAXCALLS) ) { |
1439 |
while ( fabs(l-lold)>10e-12 && (i < MAXCALLS) ) { |
1440 |
i++; |
1441 |
lold = l; |
1442 |
x = p + step; |
1443 |
l = langaufun(&x,params); |
1444 |
if ( l < lold ) step = -step/10.; |
1445 |
p += step; |
1446 |
// printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step); |
1447 |
// printf(" l %.16g \n",l); |
1448 |
}; |
1449 |
// printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step); |
1450 |
//printf(" l %.16g \n",l); |
1451 |
if (i == MAXCALLS) return (-1.); |
1452 |
// |
1453 |
return(x); |
1454 |
} |
1455 |
|
1456 |
|
1457 |
void delay(int selftrig, int &del){ |
1458 |
int i; |
1459 |
i = 0; |
1460 |
del = 300; |
1461 |
// |
1462 |
while ((i < 16) && ((selftrig & (0x8000>>i)) == 0)) { |
1463 |
del += 50; |
1464 |
i++; |
1465 |
}; |
1466 |
if (i > 15) del = 1100; |
1467 |
} |
1468 |
|
1469 |
Double_t fitraw(Double_t *x, Double_t *par){ |
1470 |
// |
1471 |
Double_t fitval = par[3] + par[0]*TMath::Exp(-par[2]*(x[0]-par[1])); |
1472 |
// |
1473 |
return fitval; |
1474 |
} |
1475 |
|
1476 |
|
1477 |
|
1478 |
int checkmerging(Chmerge & chmerge){ |
1479 |
for (Int_t i = chmerge.fromflag; i<1000; i++){ |
1480 |
if ( chmerge.pscuOBT == 0 && chmerge.pscucount == 0 ) return(1); |
1481 |
if ( chmerge.pscuOBT == chmerge.lastOBT[i] && chmerge.pscucount == chmerge.lastcount[i] ){ |
1482 |
chmerge.fromflag = i; |
1483 |
printf("\n WARNING! %i event number %i at OBT %i already processed, skipped!\n",i,chmerge.pscucount,chmerge.pscuOBT); |
1484 |
return(1); |
1485 |
}; |
1486 |
}; |
1487 |
return(0); |
1488 |
} |
1489 |
|
1490 |
int fillmerging(Chmerge & chmerge){ |
1491 |
chmerge.lastOBT[chmerge.scanflag] = chmerge.pscuOBT; |
1492 |
chmerge.lastcount[chmerge.scanflag] = chmerge.pscucount; |
1493 |
chmerge.scanflag++; |
1494 |
if ( chmerge.scanflag > 999 ) return(0); |
1495 |
chmerge.lastOBT[chmerge.scanflag] = 0; |
1496 |
chmerge.lastcount[chmerge.scanflag] = 0; |
1497 |
return(0); |
1498 |
} |
1499 |
|
1500 |
//short int checkifpretyodaversion(TString filename){ |
1501 |
// string fil = (const char *)filename; |
1502 |
// Int_t posiz = fil.find("dw_"); |
1503 |
// if ( posiz == -1 ) posiz = fil.find("DW_"); |
1504 |
// if ( posiz == -1 ) return (-1); |
1505 |
// TString year; |
1506 |
// stringcopy(year,filename,posiz+3,posiz+5); |
1507 |
// const char *ye = year; |
1508 |
// Int_t y = atoi(ye); |
1509 |
// TString month; |
1510 |
// stringcopy(month,filename,posiz+5,posiz+7); |
1511 |
// const char *mo = month; |
1512 |
// Int_t m = atoi(mo); |
1513 |
// TString day; |
1514 |
// stringcopy(day,filename,posiz+7,posiz+9); |
1515 |
// const char *da = day; |
1516 |
// Int_t d = atoi(da); |
1517 |
// TString fno; |
1518 |
// stringcopy(fno,filename,posiz+10,posiz+13); |
1519 |
// const char *fn = fno; |
1520 |
// Int_t f = atoi(fn); |
1521 |
// if ( y >= 5 && m >= 5 && d < 15 ) return(1); |
1522 |
// if ( y == 5 && m == 5 && d == 15 && f < 7 ) return(1); |
1523 |
// return(0); |
1524 |
//} |
1525 |
|
1526 |
|