]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/TOFcalibAPI.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update timestamp for new data points simulation
[u/mrichter/AliRoot.git] / TOF / TOFcalibAPI.C
CommitLineData
9d67299b 1\13/*
2 * calibapi.cxx
3 * set of functions to read TOFcompactCalib data and
4 * deal with TOF calibration in a centralized way
5 */
6
7#include <stdio.h>
8#include <signal.h>
9
10#define MAXHITS 100000
11#define MAXPOINTS 100000
12
13/* input data */
14TFile *datafile = NULL;
15TTree *datatree = NULL;
16Int_t nevents = 0;
17Int_t curev = 0;
18Int_t runNb = 0;
19UInt_t timestamp = 0;
20Float_t vertex = 0.;
21Float_t timezero = 0.;
22Int_t nhits = 0;
23Float_t momentum[MAXHITS];
24Float_t length[MAXHITS];
25Int_t index[MAXHITS];
26Float_t time[MAXHITS];
27Float_t tot[MAXHITS];
28Float_t texp[MAXHITS];
29
30/* calib histos */
31enum EParam_t {
32 kTRM,
33 kFEA,
34 kChannel,
35 kNParams
36};
37const Char_t *paramName[kNParams] = {"trm", "fea", "channel"};
38Int_t paramBins[kNParams] = {720, 6552, 157248};
39Float_t paramMin[kNParams] = {0., 0., 0.};
40Float_t paramMax[kNParams] = {720., 6552., 157248.};
41Int_t deltatBins = 201;
42Float_t deltatMin = -2440. - 12.2;
43Float_t deltatMax = 2440. + 12.2;
44Int_t totBins = 400;
45Float_t totMin = 4.88;
46Float_t totMax = 24.4;
47UInt_t timeZeroSampling = 600; /* seconds */
48UInt_t timeMin = 0;
49UInt_t timeMax = 0;
50UInt_t timeBins = 0;
51TH2F *hFEAHitMap = NULL;
52TH2F *hChannelHitMap = NULL;
53TH2F *hTimeZeroFillHisto = NULL;
54TH1F *hTimeZeroFill_mean = NULL;
55TH1F *hTimeZeroFill_sigma = NULL;
56TProfile *hTimePressureHisto = NULL;
57TH2F *hCalibHisto[kNParams] = {NULL, NULL, NULL};
58TH1F *hCalibParam_mean[kNParams] = {NULL, NULL, NULL};
59TH1F *hCalibParam_sigma[kNParams] = {NULL, NULL, NULL};
60TH1F *hChannelParam_mean = NULL;
61TH1F *hChannelParam_sigma = NULL;
62
63TH2F *hPerfHistoDeltaT = NULL;
64TH2F *hPerfHistoBeta = NULL;
65
66AliTOFcalibHisto *calibHisto = NULL;
67
68/* other */
69#define MAXRUNS 1000000
70AliLHCClockPhase *lhcClockPhase[MAXRUNS];
71AliDCSSensor *cavernPressure[MAXRUNS];
72
73/* monitoring */
74TStopwatch stopwatch;
75
76//_____________________________________________________________
77
78Bool_t init()
79{
80
81 AliCDBManager *cdb = AliCDBManager::Instance();
82 cdb->SetDefaultStorage("raw://");
83
84 /* loop over events */
85 printf("init: loop over events\n");
86 for (curev = 0; curev < nevents; curev++) {
87 /* get and check event */
88 datatree->GetEvent(curev);
89 if (cdb->GetRun() == runNb) continue;
90 /* set run and read GRP */
91 cdb->SetRun(runNb);
92 /* get LHC clock-phase */
93 AliCDBEntry *cdbe = (AliCDBEntry *)cdb->Get("GRP/Calib/LHCClockPhase");
94 lhcClockPhase[runNb] = (AliLHCClockPhase *)cdbe->GetObject();
95 }
96
97 return kFALSE;
98
99}
100
101//_____________________________________________________________
102
103Bool_t
104runCalib(const Char_t *filename, const Char_t *paramfilename = NULL, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE, Bool_t runFix = kTRUE, Int_t nSteps = 2, Int_t evMax = kMaxInt)
105{
106
107 /* open data */
108 if (openData(filename, evMax))
109 return kTRUE;
110 /* open calib params */
111 if (paramfilename)
112 if (openCalibParams(paramfilename))
113 return kTRUE;
114
115 /* init if needed */
116 if (useLHCClockPhase)
117 init();
118
119 /* fill and fit time-zero fill */
120 deltatMin = -24400. - 122.;
121 deltatMax = +24400. + 122.;
122 fillTimeZeroFillHisto(useTimeZeroTOF, useLHCClockPhase);
123 fitTimeZeroFillHisto();
124 deltatMin = -2440. - 12.2;
125 deltatMax = +2440. + 12.2;
126 fillTimeZeroFillHisto(useTimeZeroTOF, useLHCClockPhase);
127 fitTimeZeroFillHisto();
128 writeCalibHistos("calibhistos.step-timezero.root");
129 writeCalibParams("calibparams.step-timezero.root");
130
131 /* fill with extended range and fix what is far away */
132 if (runFix) {
133 /* fix up to 250 ns shift */
134 deltatMin = -244000. - 1220.;
135 deltatMax = +244000. + 1220.;
136 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
137 writeCalibHistos("calibhistos.step-fix-250ns.root");
138 fitCalibHistos(10., 10., 100, kTRUE);
139 writeCalibParams("calibparams.step-fix-250ns.root");
140 /* fix up to 25 ns shift */
141 deltatMin = -24400. - 122.;
142 deltatMax = +24400. + 122.;
143 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
144 writeCalibHistos("calibhistos.step-fix-25ns.root");
145 fitCalibHistos(10., 10., 100, kTRUE);
146 writeCalibParams("calibparams.step-fix-25ns.root");
147 }
148
149 /* fill and fit with limited range */
150 deltatMin = -2440. - 12.2;
151 deltatMax = +2440. + 12.2;
152 for (Int_t istep = 0; istep < nSteps; istep++) {
153 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
154 writeCalibHistos(Form("calibhistos.step-%d.root", istep));
155 fitCalibHistos(3., 2., 100, kFALSE);
156 writeCalibParams(Form("calibparams.step-%d.root", istep));
157 }
158
159 /* fill and fit common shift */
160 // fillCalibHistos(kFALSE, useLHCClockPhase);
161 // writeCalibHistos("calibhistos.step-common.root");
162 // fitCommonShift();
163 // writeCalibParams("calibparams.step-common.root");
164
165 /* final resuls */
166 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
167 writeCalibHistos("calibhistos.root");
168 writeCalibParams("calibparams.root");
169
170 return kFALSE;
171}
172
173//_____________________________________________________________
174
175Bool_t
176checkCalib(const Char_t *filename, const Char_t *paramfilename = NULL, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE, Bool_t runExtended = kFALSE, Int_t evMax = kMaxInt)
177{
178
179 /* open data */
180 if (openData(filename, evMax))
181 return kTRUE;
182 /* open calib params */
183 if (paramfilename)
184 if (openCalibParams(paramfilename))
185 return kTRUE;
186
187 /* fill and fit time-zero fill */
188 deltatMin = -24400.;
189 deltatMax = +24400.;
190 fillTimeZeroFillHisto(useTimeZeroTOF, useLHCClockPhase);
191 fitTimeZeroFillHisto();
192 deltatMin = -2440.;
193 deltatMax = +2440.;
194 fillTimeZeroFillHisto(useTimeZeroTOF, useLHCClockPhase);
195 fitTimeZeroFillHisto();
196
197 /* fill with extended range */
198 if (runExtended) {
199 deltatMin = -24400.;
200 deltatMax = +24400.;
201 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
202 writeCalibHistos("checkcalib.extended.root");
203 }
204
205 /* fill with limited range */
206 deltatMin = -2440.;
207 deltatMax = +2440.;
208 fillCalibHistos(useTimeZeroTOF, useLHCClockPhase);
209 writeCalibHistos("checkcalib.root");
210 writeCalibParams("checkcalib.calibparams.root");
211
212 return kFALSE;
213}
214
215//_____________________________________________________________
216
217Bool_t
218checkPerf(const Char_t *filename, const Char_t *paramfilename = NULL, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE, Int_t evMax = kMaxInt)
219{
220
221 /* open data */
222 if (openData(filename, evMax))
223 return kTRUE;
224 /* open calib params */
225 if (paramfilename)
226 if (openCalibParams(paramfilename))
227 return kTRUE;
228
229 fillTimeZeroFillHisto(useTimeZeroTOF, useLHCClockPhase);
230 fitTimeZeroFillHisto();
231 fillPerfHistos(useTimeZeroTOF, useLHCClockPhase);
232 writePerfHistos("checkperf.root");
233
234 return kFALSE;
235}
236
237//_____________________________________________________________
238
239Bool_t
240openData(const Char_t *filename, Int_t evMax = kMaxInt)
241{
242 /*
243 * open data
244 */
245
246 /* open file */
247 printf("openData: opening file: %s\n", filename);
248 datafile = TFile::Open(filename);
249 if (!datafile || !datafile->IsOpen()) {
250 printf("openData: cannot open file %s\n", filename)
251 return kTRUE;
252 }
253 /* get tree */
254 datatree = (TTree *)datafile->Get("aodTree");
255 if (!datatree) {
256 printf("openData: cannot find \'aodTree\' tree in %s\n", filename);
257 return kTRUE;
258 }
259 nevents = datatree->GetEntries();
260 printf("openData: got \'aodTree\' tree: %d entries\n", nevents);
261 if (nevents > evMax) {
262 printf("openData: setting max event to %d as requested\n", evMax);
263 nevents = evMax;
264 }
265 /* connect inputs */
266 datatree->SetBranchAddress("run", &runNb);
267 datatree->SetBranchAddress("timestamp", &timestamp);
268 datatree->SetBranchAddress("vertex", &vertex);
269 datatree->SetBranchAddress("timezero", &timezero);
270 datatree->SetBranchAddress("nhits", &nhits);
271 datatree->SetBranchAddress("momentum", &momentum);
272 datatree->SetBranchAddress("length", &length);
273 datatree->SetBranchAddress("index", &index);
274 datatree->SetBranchAddress("time", &time);
275 datatree->SetBranchAddress("tot", &tot);
276 datatree->SetBranchAddress("texp", &texp);
277
278 /* get first event, and retrieve the year */
279 datatree->GetEvent(0);
280 TTimeStamp ts(timestamp);
281 UInt_t year;
282 ts.GetDate(kTRUE, 0, &year);
283 TTimeStamp firstTimestamp(year, 1, 1, 0, 0, 0, 0, kTRUE);
284 TTimeStamp lastTimestamp(year + 1, 1, 1, 0, 0, 0, 0, kTRUE);
285 timeMin = firstTimestamp.GetTimeSpec().tv_sec;
286 timeMax = lastTimestamp.GetTimeSpec().tv_sec;
287 timeBins = (timeMax - timeMin) / timeZeroSampling;
288 printf("openData: calibration running on %d data\n", year);
289
290 return kFALSE;
291}
292
293//_____________________________________________________________
294
295Bool_t
296acceptEvent(Bool_t useTimeZeroTOF = kFALSE)
297{
298 /*
299 * acceptEvent
300 */
301
302 if (useTimeZeroTOF && timezero == 999999.) return kFALSE;
303 return kTRUE;
304}
305
306//_____________________________________________________________
307
308Float_t
309getTimeZeroFill(UInt_t ts)
310{
311 /*
312 * getTimeZeroFill
313 */
314
315 if (!hTimeZeroFill_mean) return 0.;
316 Int_t tsbin = hTimeZeroFill_mean->FindBin(ts);
317 return hTimeZeroFill_mean->GetBinContent(tsbin);
318}
319
320//_____________________________________________________________
321
322Float_t
323getLHCClockPhase(UInt_t ts)
324{
325 /*
326 * getLHCClockPhase
327 */
328
329 return lhcClockPhase[runNb] ? 1.e3 * lhcClockPhase[runNb]->GetPhase(ts) : 0.;
330}
331
332//_____________________________________________________________
333
334Float_t
335getCalib(Int_t idx)
336{
337 /*
338 * getCalib
339 */
340
341 if (!hChannelParam_mean) return 0.;
342 return hChannelParam_mean->GetBinContent(idx + 1);
343}
344
345//_____________________________________________________________
346
347Float_t
348getDeltaT(Int_t ihit, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
349{
350 /*
351 * getDeltaT
352 */
353
354 Float_t val = time[ihit] - getCalib(index[ihit]) - getTimeZeroFill(timestamp) - texp[ihit];
355 if (useTimeZeroTOF) val -= timezero ;
356 if (useLHCClockPhase) val += getLHCClockPhase(timestamp);
357 return val;
358}
359
360//_____________________________________________________________
361
362Float_t
363getBeta(Int_t ihit, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
364{
365 /*
366 * getBeta
367 */
368
369 Float_t tof = time[ihit] - getCalib(index[ihit]) - getTimeZeroFill(timestamp);
370 if (useTimeZeroTOF) tof -= timezero ;
371 if (useLHCClockPhase) tof += getLHCClockPhase(timestamp);
372 Float_t beta = length[ihit] / tof / 2.99792458000000000e-2;
373 return beta;
374}
375
376//_____________________________________________________________
377
378Float_t
379getMass(Int_t ihit, Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
380{
381 /*
382 * getMass
383 */
384
385 Float_t tof = time[ihit] - getCalib(index[ihit]) - getTimeZeroFill(timestamp);
386 if (useTimeZeroTOF) tof -= timezero ;
387 if (useLHCClockPhase) tof += getLHCClockPhase(timestamp);
388 Float_t beta = length[ihit] / tof / 2.99792458000000000e-2;
389 if (beta > 1.) beta *= -1.
390 Float_t mass = momentum[ihit] * TMath::Sqrt(1. / (beta * beta) - 1.)
391 return mass;
392}
393
394//_____________________________________________________________
395
396Int_t
397getIndex(Int_t param, Int_t idx)
398{
399 /*
400 * getIndex
401 */
402
403 if (!calibHisto) {
404 calibHisto = new AliTOFcalibHisto();
405 calibHisto->LoadCalibHisto();
406 }
407
408 /* switch param */
409 Int_t sector, ddl, trm, strip, padx, paramIndex;
410 switch (param) {
411 case kTRM:
412 ddl = calibHisto->GetCalibMap(AliTOFcalibHisto::kDDL, idx);
413 trm = calibHisto->GetCalibMap(AliTOFcalibHisto::kTRM, idx);
414 paramIndex = trm + 10 * ddl;
415 break;
416 case kFEA:
417 sector = calibHisto->GetCalibMap(AliTOFcalibHisto::kSector, idx);
418 strip = calibHisto->GetCalibMap(AliTOFcalibHisto::kSectorStrip, idx);
419 padx = calibHisto->GetCalibMap(AliTOFcalibHisto::kPadX, idx);
420 paramIndex = padx / 12 + 4 * strip + 364 * sector;
421 break;
422 case kChannel:
423 paramIndex = idx;
424 break;
425 }
426
427 return paramIndex;
428}
429
430//_____________________________________________________________
431
432Int_t
433getHitMapXY(Int_t param, Int_t idx, Float_t *hitmap)
434{
435 /*
436 * getHitMapXY
437 */
438
439 if (!calibHisto) {
440 calibHisto = new AliTOFcalibHisto();
441 calibHisto->LoadCalibHisto();
442 }
443
444 /* get from calib histo */
445 Int_t sector, strip, padx, padz, fea;
446 sector = calibHisto->GetCalibMap(AliTOFcalibHisto::kSector, idx);
447 strip = calibHisto->GetCalibMap(AliTOFcalibHisto::kSectorStrip, idx);
448 padx = calibHisto->GetCalibMap(AliTOFcalibHisto::kPadX, idx);
449 padz = calibHisto->GetCalibMap(AliTOFcalibHisto::kPadZ, idx);
450 fea = padx / 12;
451 /* switch param */
452 switch (param) {
453 case kFEA:
454 hitmap[0] = sector + ((Double_t)(3 - fea) + 0.5) / 4.;
455 hitmap[1] = strip;
456 break;
457 case kChannel:
458 hitmap[0] = sector + ((Double_t)(47 - padx) + 0.5) / 48.;
459 hitmap[1] = strip + ((Double_t)(padz) + 0.5) / 2.;
460 break;
461 default:
462 hitmap[0] = 0.;
463 hitmap[1] = 0.;
464 break;
465 }
466
467}
468
469//_____________________________________________________________
470
471void
472fillTimeZeroFillHisto(Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
473{
474 /*
475 * fillTimeZeroFillHisto
476 */
477
478 /* create/reset histos */
479 if (hTimeZeroFillHisto) delete hTimeZeroFillHisto;
480 hTimeZeroFillHisto = new TH2F("hTimeZeroFillHisto", "", timeBins, timeMin, timeMax, deltatBins, deltatMin, deltatMax);
481 if (hTimePressureHisto) delete hTimePressureHisto;
482 hTimePressureHisto = new TProfile("hTimePressureHisto", "", timeBins, timeMin, timeMax);
483
484 /* reset and start stopwatch */
485 stopwatch.Reset();
486 stopwatch.Start();
487
488 /* loop over events */
489 printf("fillTimeZeroFillHisto: loop over events\n");
490 if (useTimeZeroTOF)
491 printf("fillTimeZeroFillHisto: time-zero TOF requested\n");
492 else
493 printf("fillTimeZeroFillHisto: not using time-zero TOF\n");
494 if (useLHCClockPhase)
495 printf("fillTimeZeroFillHisto: BPTX clock-phase requested\n");
496 else
497 printf("fillTimeZeroFillHisto: not using BPTX clock-phase\n");
498 Float_t hitmap[2], deltat;
499 for (curev = 0; curev < nevents; curev++) {
500 /* get and check event */
501 datatree->GetEvent(curev);
502 if (!acceptEvent(useTimeZeroTOF)) continue;
503 /* loop over hits */
504 for (Int_t ihit = 0; ihit < nhits; ihit++) {
505 deltat = getDeltaT(ihit, useTimeZeroTOF, useLHCClockPhase);
506 /* fill histos */
507 hTimeZeroFillHisto->Fill(timestamp, deltat);
508 } /* end of loop over hits */
509 } /* end of loop over events */
510
511 /* print monitor */
512 monitor();
513
514}
515
516//_____________________________________________________________
517
518void
519fillCalibHistos(Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
520{
521 /*
522 * fillCalibHistos
523 */
524
525 /* create/reset histos */
526 if (!hFEAHitMap)
527 hFEAHitMap = new TH2F("hFEAHitMap", "", 72, 0., 18., 91, 0., 91.);
528 else
529 hFEAHitMap->Reset();
530
531 if (!hChannelHitMap)
532 hChannelHitMap = new TH2F("hChannelHitMap", "", 864, 0., 18., 91, 0., 91.);
533 else
534 hChannelHitMap->Reset();
535
536 for (Int_t iparam = 0; iparam < kNParams; iparam++) {
537 if (hCalibHisto[iparam]) delete hCalibHisto[iparam];
538 hCalibHisto[iparam] = new TH2F(Form("hCalibHisto_%s", paramName[iparam]), "", paramBins[iparam], paramMin[iparam], paramMax[iparam], deltatBins, deltatMin, deltatMax);
539 }
540
541 /* reset and start stopwatch */
542 stopwatch.Reset();
543 stopwatch.Start();
544
545 /* loop over events */
546 printf("fillCalibHistos: loop over events\n");
547 if (useTimeZeroTOF)
548 printf("fillCalibHistos: time-zero TOF requested\n");
549 else
550 printf("fillCalibHistos: not using time-zero TOF\n");
551 if (useLHCClockPhase)
552 printf("fillTimeZeroFillHisto: BPTX clock-phase requested\n");
553 else
554 printf("fillTimeZeroFillHisto: not using BPTX clock-phase\n");
555 Float_t hitmap[2], deltat;
556 for (curev = 0; curev < nevents; curev++) {
557 /* get and check event */
558 datatree->GetEvent(curev);
559 if (!acceptEvent(useTimeZeroTOF)) continue;
560 /* loop over hits */
561 for (Int_t ihit = 0; ihit < nhits; ihit++) {
562 deltat = getDeltaT(ihit, useTimeZeroTOF, useLHCClockPhase);
563 /* fill histos */
564 getHitMapXY(kFEA, index[ihit], hitmap);
565 hFEAHitMap->Fill(hitmap[0], hitmap[1]);
566 getHitMapXY(kChannel, index[ihit], hitmap);
567 hChannelHitMap->Fill(hitmap[0], hitmap[1]);
568 for (Int_t iparam = 0; iparam < kNParams; iparam++) {
569 hCalibHisto[iparam]->Fill(getIndex(iparam, index[ihit]), deltat);
570 }
571 } /* end of loop over hits */
572 } /* end of loop over events */
573
574 /* print monitor */
575 monitor();
576
577}
578
579//_____________________________________________________________
580
581void
582fillPerfHistos(Bool_t useTimeZeroTOF = kFALSE, Bool_t useLHCClockPhase = kFALSE)
583{
584 /*
585 * fillPerfHistos
586 */
587
588 /* create/reset histos */
589 if (!hFEAHitMap)
590 hFEAHitMap = new TH2F("hFEAHitMap", "", 72, 0., 18., 91, 0., 91.);
591 else
592 hFEAHitMap->Reset();
593
594 if (!hChannelHitMap)
595 hChannelHitMap = new TH2F("hChannelHitMap", "", 864, 0., 18., 91, 0., 91.);
596 else
597 hChannelHitMap->Reset();
598
599 if (!hPerfHistoDeltaT)
600 hPerfHistoDeltaT = new TH2F("hPerfHistoDeltaT", "", 1000, 0., 10., 2000, -24400., 24400.);
601 else
602 hPerfHistoDeltaT->Reset();
603
604 if (!hPerfHistoBeta)
605 hPerfHistoBeta = new TH2F("hPerfHistoBeta", "", 1000, 0., 10., 2200, 0., 1.1);
606 else
607 hPerfHistoBeta->Reset();
608
609
610 /* reset and start stopwatch */
611 stopwatch.Reset();
612 stopwatch.Start();
613
614 /* loop over events */
615 printf("fillPerfHistos: loop over events\n");
616 if (useTimeZeroTOF)
617 printf("fillPerfHistos: time-zero TOF requested\n");
618 else
619 printf("fillPerfHistos: not using time-zero TOF\n");
620 if (useLHCClockPhase)
621 printf("fillPerfHisto: BPTX clock-phase requested\n");
622 else
623 printf("fillPerfHisto: not using BPTX clock-phase\n");
624 Float_t hitmap[2], deltat, beta, mass;
625 for (curev = 0; curev < nevents; curev++) {
626 /* get and check event */
627 datatree->GetEvent(curev);
628 if (!acceptEvent(useTimeZeroTOF)) continue;
629 /* loop over hits */
630 for (Int_t ihit = 0; ihit < nhits; ihit++) {
631 deltat = getDeltaT(ihit, useTimeZeroTOF, useLHCClockPhase);
632 beta = getBeta(ihit, useTimeZeroTOF, useLHCClockPhase);
633
634 // mass = getMass(ihit, useTimeZeroTOF, useLHCClockPhase);
635 /* fill histos */
636 getHitMapXY(kFEA, index[ihit], hitmap);
637 hFEAHitMap->Fill(hitmap[0], hitmap[1]);
638 getHitMapXY(kChannel, index[ihit], hitmap);
639 hChannelHitMap->Fill(hitmap[0], hitmap[1]);
640 hPerfHistoDeltaT->Fill(momentum[ihit], deltat);
641 hPerfHistoBeta->Fill(momentum[ihit], beta);
642 // hPerfHistoMass->Fill(p, mass);
643 } /* end of loop over hits */
644 } /* end of loop over events */
645
646 /* print monitor */
647 monitor();
648
649}
650
651//_____________________________________________________________
652
653Bool_t
654writeCalibHistos(const Char_t *filename)
655{
656 /*
657 * writeCalibHistos
658 */
659
660 /* open output file and write */
661 TFile *fileout = TFile::Open(filename, "RECREATE");
662 if (!fileout || !fileout->IsOpen()) {
663 printf("writeCalibHistos: error while opening output file %s\n", filename);
664 return kTRUE;
665 }
666 if (hFEAHitMap) hFEAHitMap->Write();
667 if (hChannelHitMap) hChannelHitMap->Write();
668 if (hTimeZeroFillHisto) hTimeZeroFillHisto->Write();
669 if (hTimePressureHisto) hTimePressureHisto->Write();
670 for (Int_t iparam = 0; iparam < kNParams; iparam++) {
671 if (hCalibHisto[iparam]) hCalibHisto[iparam]->Write();
672 }
673 fileout->Close();
674 printf("writeCalibHistos: output written on %s\n", filename);
675 return kFALSE;
676}
677
678//_____________________________________________________________
679
680Bool_t
681writePerfHistos(const Char_t *filename)
682{
683 /*
684 * writePerfHistos
685 */
686
687 /* open output file and write */
688 TFile *fileout = TFile::Open(filename, "RECREATE");
689 if (!fileout || !fileout->IsOpen()) {
690 printf("writePerfHistos: error while opening output file %s\n", filename);
691 return kTRUE;
692 }
693 if (hFEAHitMap) hFEAHitMap->Write();
694 if (hChannelHitMap) hChannelHitMap->Write();
695 if (hPerfHistoDeltaT) hPerfHistoDeltaT->Write();
696 if (hPerfHistoBeta) hPerfHistoBeta->Write();
697 fileout->Close();
698 printf("writePerfHistos: output written on %s\n", filename);
699 return kFALSE;
700}
701
702//_____________________________________________________________
703
704Bool_t
705writeCalibParams(const Char_t *filename)
706{
707 /*
708 * writeCalibParams
709 */
710
711 /* open output file and write */
712 TFile *fileout = TFile::Open(filename, "RECREATE");
713 if (!fileout || !fileout->IsOpen()) {
714 printf("writeCalibParams: error while opening output file %s\n", filename);
715 return kTRUE;
716 }
717 if (hTimeZeroFill_mean) hTimeZeroFill_mean->Write();
718 if (hTimeZeroFill_sigma) hTimeZeroFill_sigma->Write();
719 for (Int_t iparam = 0; iparam < kNParams; iparam++) {
720 if (hCalibParam_mean[iparam]) hCalibParam_mean[iparam]->Write();
721 if (hCalibParam_sigma[iparam]) hCalibParam_sigma[iparam]->Write();
722 }
723 if (hChannelParam_mean) hChannelParam_mean->Write();
724 if (hChannelParam_sigma) hChannelParam_sigma->Write();
725 fileout->Close();
726 printf("writeCalibParams: output written on %s\n", filename);
727 return kFALSE;
728}
729
730//_____________________________________________________________
731
732Bool_t
733openCalibHistos(const Char_t *filename)
734{
735 /*
736 * openCalibHistos
737 */
738
739 /* open file and write */
740 printf("openCalibHistos: opening file: %s\n", filename);
741 TFile *filein = TFile::Open(filename);
742 if (!filein || !filein->IsOpen()) {
743 printf("openCalibHistos: cannot open file %s\n", filename);
744 return kTRUE;
745 }
746 /* time-zero fill */
747 if (hTimeZeroFillHisto) {
748 printf("openCalibHistos: histo \'hTimeZeroFillHisto\' will be replaced\n");
749 delete hTimeZeroFillHisto;
750 }
751 hTimeZeroFillHisto = (TH2F *)filein->Get("hTimeZeroFillHisto");
752 if (!hTimeZeroFillHisto) {
753 printf("openCalibHistos: cannot get \'hTimeZeroFillHisto\' from file %s\n", filename);
754 return kTRUE;
755 }
756 printf("openCalibHistos: got \'hTimeZeroFillHisto\' from file\n");
757 for (Int_t iparam = 0; iparam < kNParams; iparam++) {
758 /* calib */
759 if (hCalibHisto[iparam]) {
760 printf("openCalibHistos: histo \'hCalibHisto_%s\' will be replaced\n", paramName[iparam]);
761 delete hCalibHisto[iparam];
762 }
763 hCalibHisto[iparam] = (TH2F *)filein->Get(Form("hCalibHisto_%s", paramName[iparam]));
764 if (!hCalibHisto[iparam]) {
765 printf("openCalibHistos: cannot get \'hCalibHisto_%s\' from file %s\n", paramName[iparam], filename);
766 return kTRUE;
767 }
768 printf("openCalibHistos: got \'hCalibHisto_%s\' from file\n", paramName[iparam]);
769 }
770
771 return kFALSE;
772}
773
774//_____________________________________________________________
775
776Bool_t
777openCalibParams(const Char_t *filename)
778{
779 /*
780 * openCalibParams
781 */
782
783 /* open file */
784 printf("openCalibParams: opening file: %s\n", filename);
785 TFile *filein = TFile::Open(filename);
786 if (!filein || !filein->IsOpen()) {
787 printf("openCalibParams: cannot open file %s\n", filename);
788 return kTRUE;
789 }
790 /* mean */
791 if (hChannelParam_mean) {
792 printf("openCalibParams: histo \'hChannelParam_mean\' will be replaced\n");
793 delete hChannelParam_mean;
794 }
795 hChannelParam_mean = (TH1F *)filein->Get("hChannelParam_mean");
796 if (!hChannelParam_mean) {
797 printf("openCalibParams: cannot get \'hChannelParam_mean\' from file %s\n", filename);
798 return kTRUE;
799 }
800 printf("openCalibParams: got \'hChannelParam_mean\' from file\n");
801 /* sigma */
802 if (hChannelParam_sigma) {
803 printf("openCalibParams: histo \'hChannelParam_sigma\' will be replaced\n");
804 delete hChannelParam_sigma;
805 }
806 hChannelParam_sigma = (TH1F *)filein->Get("hChannelParam_sigma");
807 if (!hChannelParam_sigma) {
808 printf("openCalibParams: cannot get \'hChannelParam_sigma\' from file %s\n", filename);
809 return kTRUE;
810 }
811 printf("openCalibParams: got \'hChannelParam_sigma\' from file\n");
812
813 return kFALSE;
814}
815
816//_____________________________________________________________
817
818Bool_t
819fitTimeZeroFillHisto(Float_t nSigmaMin = 2., Float_t nSigmaMax = 1., Int_t minIntegral = 1000, Bool_t useMaxBin = kFALSE)
820{
821 /*
822 * fitTimeZeroFillHisto
823 */
824
825 // if (!hChannelParam_mean || !hCalibHisto[kChannel]) return kTRUE;
826
827 if (!hTimeZeroFill_mean) hTimeZeroFill_mean = new TH1F("hTimeZeroFill_mean", "", timeBins, timeMin, timeMax);
828 if (!hTimeZeroFill_sigma) hTimeZeroFill_sigma = new TH1F("hTimeZeroFill_sigma", "", timeBins, timeMin, timeMax);
829
830 if (useMaxBin)
831 printf("fitTimeZeroFill: fitting time-zero fill (using max bin, intMin=%d)\n", minIntegral);
832 else
833 printf("fitTimeZeroFill: fitting time-zero fill (sigmaMin=%.1f, sigmaMax=%.1f, intMin=%d)\n", nSigmaMin, nSigmaMax, minIntegral);
834
835 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
836 for (Int_t ibin = 0; ibin < hTimeZeroFillHisto->GetNbinsX(); ibin++) {
837 TH1D *hpy = hTimeZeroFillHisto->ProjectionY("hpy", ibin + 1, ibin +1);
838 if (hpy->GetEntries() <= minIntegral) {
839 delete hpy;
840 continue;
841 }
842 if (fitPeak(fitFunc, hpy, nSigmaMin, nSigmaMax) != 0) {
843 delete hpy;
844 continue;
845 }
846 hTimeZeroFill_mean->AddBinContent(ibin + 1, fitFunc->GetParameter(1));
847 hTimeZeroFill_mean->SetBinError(ibin + 1, fitFunc->GetParError(1));
848 hTimeZeroFill_sigma->SetBinContent(ibin + 1, fitFunc->GetParameter(2));
849 hTimeZeroFill_sigma->SetBinError(ibin + 1, fitFunc->GetParError(2));
850 delete hpy;
851 }
852
853 return kFALSE;
854
855}
856
857//_____________________________________________________________
858
859Bool_t
860fitCommonShift(Float_t nSigmaMin = 2., Float_t nSigmaMax = 1., Int_t minIntegral = 1000, Bool_t useMaxBin = kFALSE)
861{
862 /*
863 * fitCommonShift
864 */
865
866 if (!hChannelParam_mean || !hCalibHisto[kChannel]) return kTRUE;
867
868 if (useMaxBin)
869 printf("fitCommonShift: fitting common shift (using max bin, intMin=%d)\n", minIntegral);
870 else
871 printf("fitCommonShift: fitting common shift (sigmaMin=%.1f, sigmaMax=%.1f, intMin=%d)\n", nSigmaMin, nSigmaMax, minIntegral);
872
873 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
874 TH1D *hpy = hCalibHisto[kChannel]->ProjectionY("hpy");
875 if (hpy->GetEntries() <= minIntegral) {
876 delete hpy;
877 continue;
878 }
879 if (useMaxBin) {
880 for (Int_t idx = 0; idx < 157248; idx++)
881 hChannelParam_mean->AddBinContent(idx + 1, hpy->GetBinCenter(hpy->GetMaximumBin()));
882 delete hpy;
883 return kFALSE;
884 }
885 if (fitPeak(fitFunc, hpy, nSigmaMin, nSigmaMax) != 0) {
886 delete hpy;
887 return kTRUE;
888 }
889 for (Int_t idx = 0; idx < 157248; idx++)
890 hChannelParam_mean->AddBinContent(idx + 1, fitFunc->GetParameter(1));
891 delete hpy;
892
893 return kFALSE;
894
895}
896
897//_____________________________________________________________
898
899Bool_t
900fitCalibHistos(Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral, Bool_t useMaxBin)
901{
902 /*
903 * fitCalibHistos
904 */
905
906 for (Int_t iparam = 0; iparam < kNParams; iparam++)
907 fitCalibHisto(iparam, nSigmaMin, nSigmaMax, minIntegral, useMaxBin);
908 updateChannelParams();
909
910}
911
912//_____________________________________________________________
913
914Bool_t
915updateChannelParams()
916{
917
918 if (!hChannelParam_mean) hChannelParam_mean = new TH1F("hChannelParam_mean", "", paramBins[kChannel], paramMin[kChannel], paramMax[kChannel]);
919 if (!hChannelParam_sigma) hChannelParam_sigma = new TH1F("hChannelParam_sigma", "", paramBins[kChannel], paramMin[kChannel], paramMax[kChannel]);
920
921 Double_t mean, meanerr, sigma, sigmaerr;
922 for (Int_t idx = 0; idx < 157248; idx++) {
923 for (Int_t iparam = kNParams - 1; iparam >= 0; iparam--) {
924 if (!hCalibParam_mean[iparam]) continue;
925 if (hCalibParam_mean[iparam]->GetBinError(getIndex(iparam, idx) + 1) == 0.) continue;
926 mean = hCalibParam_mean[iparam]->GetBinContent(getIndex(iparam, idx) + 1);
927 meanerr = hCalibParam_mean[iparam]->GetBinError(getIndex(iparam, idx) + 1);
928 sigma = hCalibParam_sigma[iparam]->GetBinContent(getIndex(iparam, idx) + 1);
929 sigmaerr = hCalibParam_sigma[iparam]->GetBinError(getIndex(iparam, idx) + 1);
930 hChannelParam_mean->AddBinContent(idx + 1, mean);
931 hChannelParam_mean->SetBinError(idx + 1, meanerr);
932 hChannelParam_sigma->SetBinContent(idx + 1, sigma);
933 hChannelParam_sigma->SetBinError(idx + 1, sigmaerr);
934 break;
935 }
936 }
937
938}
939
940//_____________________________________________________________
941
942Bool_t
943fitCalibHisto(Int_t param, Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral, Bool_t useMaxBin)
944{
945 /*
946 * fitCalibHisto
947 */
948
949 /* check data histo */
950 if (!hCalibHisto[param]) {
951 printf("fitCalibHisto: cannot get \'hCalibHisto_%s\'\n", paramName[param]);
952 return kTRUE;
953 }
954 /* check mean histo */
955 if (!hCalibParam_mean[param]) hCalibParam_mean[param] = new TH1F(Form("hCalibParam_%s_mean", paramName[param]), "", paramBins[param], paramMin[param], paramMax[param]);
956 else hCalibParam_mean[param]->Reset();
957 /* check sigma histo */
958 if (!hCalibParam_sigma[param]) hCalibParam_sigma[param] = new TH1F(Form("hCalibParam_%s_sigma", paramName[param]), "", paramBins[param], paramMin[param], paramMax[param]);
959 else hCalibParam_sigma[param]->Reset();
960 /* fit */
961 return fitCalibHisto(hCalibHisto[param], hCalibParam_mean[param], hCalibParam_sigma[param], nSigmaMin, nSigmaMax, minIntegral, useMaxBin);
962}
963
964//_____________________________________________________________
965
966Bool_t
967fitCalibHisto(TH2F *hdata, TH1F *hmean, TH1F *hsigma, Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral, Bool_t useMaxBin)
968{
969 /*
970 * fitCalibHisto
971 */
972
973 if (!hdata || !hmean || !hsigma) return kTRUE;
974 if (useMaxBin)
975 printf("fitCalibHisto: fitting \'%s\' (using max bin, intMin=%d)\n", hdata->GetName(), minIntegral);
976 else
977 printf("fitCalibHisto: fitting \'%s\' (sigmaMin=%.1f, sigmaMax=%.1f, intMin=%d)\n", hdata->GetName(), nSigmaMin, nSigmaMax, minIntegral);
978 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
979 Int_t nDone = 0;
980 for (Int_t ibin = 0; ibin < hdata->GetNbinsX(); ibin++) {
981 TH1D *hpy = hdata->ProjectionY("hpy", ibin + 1, ibin + 1);
982 if (hpy->GetEntries() <= minIntegral) {
983 delete hpy;
984 continue;
985 }
986 if (useMaxBin) {
987 hmean->SetBinContent(ibin + 1, hpy->GetBinCenter(hpy->GetMaximumBin()));
988 hmean->SetBinError(ibin + 1, hpy->GetBinWidth(hpy->GetMaximumBin()));
989 hsigma->SetBinContent(ibin + 1, hpy->GetBinWidth(hpy->GetMaximumBin()));
990 hsigma->SetBinError(ibin + 1, hpy->GetBinWidth(hpy->GetMaximumBin()));
991 // hmean->SetBinContent(ibin + 1, hpy->GetMean());
992 // hmean->SetBinError(ibin + 1, hpy->GetMeanError());
993 // hsigma->SetBinContent(ibin + 1, hpy->GetRMS());
994 // hsigma->SetBinError(ibin + 1, hpy->GetRMSError());
995 nDone++;
996 delete hpy;
997 continue;
998 }
999 /* fit and check result */
1000 if (fitPeak(fitFunc, hpy, nSigmaMin, nSigmaMax) != 0) {
1001 delete hpy;
1002 continue;
1003 }
1004#if 0
1005 /* check mean value larger than error */
1006 if (TMath::Abs(fitFunc->GetParameter(1)) < fitFunc->GetParameter(2)) {
1007 delete hpy;
1008 continue;
1009 }
1010#endif
1011 hmean->SetBinContent(ibin + 1, fitFunc->GetParameter(1));
1012 hmean->SetBinError(ibin + 1, fitFunc->GetParError(1));
1013 hsigma->SetBinContent(ibin + 1, fitFunc->GetParameter(2));
1014 hsigma->SetBinError(ibin + 1, fitFunc->GetParError(2));
1015 nDone++;
1016 delete hpy;
1017 }
1018 printf("fitCalibHisto: %d done\n", nDone);
1019 return kFALSE;
1020}
1021
1022//_____________________________________________________________
1023
1024problematicChannels(const Char_t *filename, Float_t nsigma = 5., Int_t minHits = 0, Int_t startRun = 0, Int_t endRun = AliCDBRunRange::Infinity(), const Char_t *dbString = "local://$HOME/OCDB")
1025{
1026
1027 TFile *fin = TFile::Open(filename);
1028 TH2 *hin = (TH2 *)fin->Get("hCalibHisto_channel");
1029 TH2 *hingood = (TH2 *)hin->Clone("hingood");
1030 TH2 *hinbad = (TH2 *)hin->Clone("hinbad");
1031
1032 TH1D *hpy;
1033
1034 TH1F *hEntries = new TH1F("hEntries", "", 10. * hin->GetEntries() / 157248, 0., 5. * hin->GetEntries() / 157248);
1035 for (Int_t ich = 0; ich < 157248; ich++) {
1036 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1037 hEntries->Fill(hpy->GetEntries());
1038 }
1039
1040 new TCanvas("cEntries");
1041 hEntries->Draw();
1042
1043 TH1D *hinpy_all = hin->ProjectionY("hinpy_all");
1044 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
1045 fitPeak(fitFunc, hinpy_all, 2., 1.);
1046 Double_t mean = fitFunc->GetParameter(1);
1047 Double_t sigma = fitFunc->GetParameter(2);
1048 Double_t intmin = mean - nsigma * sigma;
1049 Double_t intmax = mean + nsigma * sigma;
1050 Int_t binmin = hinpy_all->FindBin(intmin);
1051 Int_t binmax = hinpy_all->FindBin(intmax);
1052
1053 Double_t integral_all = hinpy_all->Integral(binmin, binmax);
1054 Double_t entries_all = hinpy_all->GetEntries();
1055 Double_t fraction_all = integral_all / entries_all;
1056
1057 printf("found %f hits in +- %d sigma\n", fraction_all, nsigma);
1058
1059 TH1D *hFraction = new TH1D("hFraction", "", 2000, 0., 2.);
1060
1061 TH1C *obj = new TH1C("hProblematic", "", 157248, 0., 157248.);
1062
1063 TH2F *hFEAHitMap = new TH2F("hFEAHitMap", "", 72, 0., 18., 91, 0., 91.);
1064
1065 Float_t hitmap[2];
1066
1067 for (Int_t ich = 0; ich < 157248; ich++) {
1068 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1069 if (hpy->GetEntries() <= minHits) {
1070 delete hpy;
1071 continue;
1072 }
1073 Double_t integral_ch = hpy->Integral(binmin, binmax);
1074 Double_t entries_ch = hpy->GetEntries();
1075 Double_t fraction_ch = integral_ch / entries_ch / fraction_all;
1076 hFraction->Fill(fraction_ch);
1077
1078 }
1079
1080 TF1 *dgaus = new TF1("dgaus", "[0] * ((x < [1]) * TMath::Gaus(x, [1], [2]) + (x > [1]) * TMath::Gaus(x, [1], [3]))", 0., 2.);
1081 dgaus->SetParameter(0, 1.);
1082 dgaus->SetParameter(1, 1.);
1083 dgaus->SetParameter(2, 1.);
1084 dgaus->SetParameter(3, 1.);
1085 dgaus->SetParLimits(2, 0., 1.);
1086 dgaus->SetParLimits(3, 0., 1.);
1087
1088 hFraction->Fit(dgaus, "0", "");
1089 hFraction->Fit(dgaus, "0", "", dgaus->GetParameter(1) - 3. * dgaus->GetParameter(2), dgaus->GetParameter(1) + 3. * dgaus->GetParameter(3));
1090 Double_t cutLimit = dgaus->GetParameter(1) - 3. * dgaus->GetParameter(2);
1091 Double_t cutLimit2 = dgaus->GetParameter(1) - 3. * dgaus->GetParameter(3);
1092
1093 TLine *lcutfrac = new TLine(cutLimit, 0, cutLimit, hFraction->GetMaximum());
1094 TLine *lcutfrac2 = new TLine(cutLimit2, 0, cutLimit2, hFraction->GetMaximum());
1095
1096 Int_t ntotch = 0, nprobch = 0;
1097
1098 for (Int_t ich = 0; ich < 157248; ich++) {
1099 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1100 if (hpy->GetEntries() <= minHits) {
1101 delete hpy;
1102 continue;
1103 }
1104 Double_t integral_ch = hpy->Integral(binmin, binmax);
1105 Double_t entries_ch = hpy->GetEntries();
1106 Double_t fraction_ch = integral_ch / entries_ch / fraction_all;
1107
1108 ntotch++;
1109 if (fraction_ch < cutLimit) {
1110 nprobch++;
1111 obj->SetBinContent(ich + 1, 0x1);
1112 getHitMapXY(kFEA, ich, hitmap);
1113 hFEAHitMap->Fill(hitmap[0], hitmap[1]);
1114 for (Int_t i = 0; i < hingood->GetNbinsY(); i++)
1115 hingood->SetBinContent(ich + 1, i + 1, 0);
1116 }
1117 else {
1118 for (Int_t i = 0; i < hinbad->GetNbinsY(); i++)
1119 hinbad->SetBinContent(ich + 1, i + 1, 0);
1120 }
1121 }
1122
1123 printf("flagged %d problematic over %d checked channels\n", nprobch, ntotch);
1124
1125 new TCanvas("cFraction");
1126 hFraction->Draw();
1127 dgaus->Draw("same");
1128 lcutfrac->Draw("same");
1129 lcutfrac2->Draw("same");
1130 gPad->Update();
1131 new TCanvas("cGood");
1132 hingood->Draw("colz");
1133 gPad->Update();
1134 new TCanvas("cBad");
1135 hinbad->Draw("colz");
1136 gPad->Update();
1137 new TCanvas;
1138 hin->Draw("colz");
1139 new TCanvas("hFEA");
1140 hFEAHitMap->Draw("colz");
1141 gPad->Update();
1142
1143 TFile *fout = TFile::Open("problematicChannels.root", "RECREATE");
1144 hFraction->Write();
1145 hingood->Write();
1146 hinbad->Write();
1147 hFEAHitMap->Write("hBadFEAMap");
1148 fout->Close();
1149
1150 /* create cdb info */
1151 AliCDBId id("TOF/Calib/Problematic", startRun, endRun);
1152 AliCDBMetaData *md = new AliCDBMetaData();
1153 md->SetResponsible("Roberto Preghenella");
1154 md->SetComment("Problematic");
1155 md->SetAliRootVersion(gSystem->Getenv("ARVERSION"));
1156 md->SetBeamPeriod(0);
1157
1158 /* put object in cdb */
1159 AliCDBManager *cdb = AliCDBManager::Instance();
1160 cdb->SetDefaultStorage(dbString);
1161 cdb->GetDefaultStorage()->Put(obj, id, md);
1162
1163}
1164
1165//_____________________________________________________________
1166
1167problematicFEAs(const Char_t *filename, Float_t nsigma = 3.)
1168{
1169
1170 /* map FEA-channels */
1171 Float_t hitmap[6552][2];
1172 for (Int_t ich = 0; ich < 157248; ich++) {
1173 Int_t ifea = getIndex(kFEA, ich);
1174 getHitMapXY(kFEA, ich, hitmap[ifea]);
1175 }
1176
1177 TFile *fin = TFile::Open(filename);
1178 TH2 *hin = (TH2 *)fin->Get("hCalibHisto_fea");
1179 TH2 *hingood = (TH2 *)hin->Clone("hingood");
1180 TH2 *hinbad = (TH2 *)hin->Clone("hinbad");
1181
1182 TH1D *hinpy_all = hin->ProjectionY("hinpy_all");
1183 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
1184 fitPeak(fitFunc, hinpy_all, 2., 1.);
1185 Double_t mean = fitFunc->GetParameter(1);
1186 Double_t sigma = fitFunc->GetParameter(2);
1187 Double_t intmin = mean - nsigma * sigma;
1188 Double_t intmax = mean + nsigma * sigma;
1189 Int_t binmin = hinpy_all->FindBin(intmin);
1190 Int_t binmax = hinpy_all->FindBin(intmax);
1191
1192 Double_t integral_all = hinpy_all->Integral(binmin, binmax);
1193 Double_t entries_all = hinpy_all->GetEntries();
1194 Double_t fraction_all = integral_all / entries_all;
1195
1196 printf("found %f hits in +- 3 sigma\n", fraction_all);
1197
1198 TH1D *hFraction = new TH1D("hFraction", "", 2000, 0., 2.);
1199
1200 TH1C *obj = new TH1C("hProblematic", "", 157248, 0., 157248.);
1201
1202 TH2F *hFEAHitMap = new TH2F("hFEAHitMap", "", 72, 0., 18., 91, 0., 91.);
1203 TH2F *hBadFEAHitMap = new TH2F("hBadFEAHitMap", "", 72, 0., 18., 91, 0., 91.);
1204
1205 Float_t hitmap[2];
1206
1207 TH1D *hpy;
1208 for (Int_t ich = 0; ich < 6552; ich++) {
1209 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1210 if (hpy->GetEntries() <= 0) {
1211 delete hpy;
1212 continue;
1213 }
1214 Double_t integral_ch = hpy->Integral(binmin, binmax);
1215 Double_t entries_ch = hpy->GetEntries();
1216 Double_t fraction_ch = integral_ch / entries_ch;// / fraction_all;
1217 hFraction->Fill(fraction_ch);
1218
1219 }
1220 TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus");
1221 hFraction->Fit(gaus, "WW");
1222 hFraction->Fit(gaus, "WW", "", gaus->GetParameter(1) - gaus->GetParameter(2), 2.);
1223 fraction_all = gaus->GetParameter(1);
1224 Float_t cutLimit = gaus->GetParameter(1) - 3. * gaus->GetParameter(2);
1225 cutLimit /= fraction_all;
1226
1227 hFraction->Reset();
1228 for (Int_t ich = 0; ich < 6552; ich++) {
1229 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1230 if (hpy->GetEntries() <= 0) {
1231 delete hpy;
1232 continue;
1233 }
1234 Double_t integral_ch = hpy->Integral(binmin, binmax);
1235 Double_t entries_ch = hpy->GetEntries();
1236 Double_t fraction_ch = integral_ch / entries_ch / fraction_all;
1237 hFraction->Fill(fraction_ch);
1238
1239 hFEAHitMap->Fill(hitmap[ich][0], hitmap[ich][1], fraction_ch);
1240 if (fraction_ch < cutLimit)
1241 hBadFEAHitMap->Fill(hitmap[ich][0], hitmap[ich][1]);
1242 }
1243
1244 new TCanvas;
1245 hFraction->Draw();
1246 new TCanvas;
1247 hFEAHitMap->Draw("colz");
1248 new TCanvas;
1249 hBadFEAHitMap->Draw("colz");
1250
1251 TFile *fout = TFile::Open("problematicFEAs.root", "RECREATE");
1252 hFraction->Write();
1253 hFEAHitMap->Write("hFEAFractionMap");
1254 hBadFEAHitMap->Write("hBadFEAMap");
1255 fout->Close();
1256
1257 return;
1258
1259 TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus");
1260 hFraction->Fit(gaus, "WW");
1261 Double_t cutLimit = gaus->GetParameter(1) - 3. * gaus->GetParameter(2);
1262
1263 for (Int_t ich = 0; ich < 157248; ich++) {
1264 hpy = hin->ProjectionY("hpy", ich + 1, ich + 1);
1265 if (hpy->GetEntries() <= 0) {
1266 delete hpy;
1267 continue;
1268 }
1269 Double_t integral_ch = hpy->Integral(binmin, binmax);
1270 Double_t entries_ch = hpy->GetEntries();
1271 Double_t fraction_ch = integral_ch / entries_ch / fraction_all;
1272
1273 if (fraction_ch < cutLimit) {
1274 obj->SetBinContent(ich + 1, 0x1);
1275 getHitMapXY(kFEA, ich, hitmap);
1276 hFEAHitMap->Fill(hitmap[0], hitmap[1]);
1277 for (Int_t i = 0; i < hingood->GetNbinsY(); i++)
1278 hingood->SetBinContent(ich + 1, i + 1, 0);
1279 }
1280 else {
1281 for (Int_t i = 0; i < hinbad->GetNbinsY(); i++)
1282 hinbad->SetBinContent(ich + 1, i + 1, 0);
1283 }
1284 }
1285
1286 hFraction->Draw();
1287 new TCanvas;
1288 hingood->Draw("colz");
1289 new TCanvas;
1290 hinbad->Draw("colz");
1291 new TCanvas;
1292 hin->Draw("colz");
1293 new TCanvas;
1294 hFEAHitMap->Draw("colz");
1295
1296 TFile *fout = TFile::Open("problematicChannels.root", "RECREATE");
1297 hFraction->Write();
1298 hingood->Write();
1299 hinbad->Write();
1300 hFEAHitMap->Write("hBadFEAMap");
1301 fout->Close();
1302
1303 /* create cdb info */
1304 AliCDBId id("TOF/Calib/Problematic", startRun, endRun);
1305 AliCDBMetaData *md = new AliCDBMetaData();
1306 md->SetResponsible("Roberto Preghenella");
1307 md->SetComment("Problematic");
1308 md->SetAliRootVersion(gSystem->Getenv("ARVERSION"));
1309 md->SetBeamPeriod(0);
1310
1311 /* put object in cdb */
1312 AliCDBManager *cdb = AliCDBManager::Instance();
1313 cdb->SetDefaultStorage(dbString);
1314 cdb->GetDefaultStorage()->Put(obj, id, md);
1315
1316}
1317
1318//_____________________________________________________________
1319
1320mergeProblematicMaps(const Char_t *mapfilename1, const Char_t *mapfilename2, Int_t startRun = 0, Int_t endRun = AliCDBRunRange::Infinity(), const Char_t *dbString = "local://$HOME/OCDB")
1321{
1322 TFile *file1 = TFile::Open(mapfilename1);
1323 TFile *file2 = TFile::Open(mapfilename2);
1324
1325 AliCDBEntry *cdbe1 = (AliCDBEntry *)file1->Get("AliCDBEntry");
1326 AliCDBEntry *cdbe2 = (AliCDBEntry *)file2->Get("AliCDBEntry");
1327
1328 TH1C *h1 = (TH1C *)cdbe1->GetObject();
1329 TH1C *h2 = (TH1C *)cdbe2->GetObject();
1330 TH1C *obj = new TH1C("hProblematic", "", 157248, 0., 157248.);
1331
1332 for (Int_t ich = 0; ich < 157248; ich++) {
1333 if (h1->GetBinContent(ich + 1) || h2->GetBinContent(ich + 1))
1334 obj->SetBinContent(ich + 1, 0x1);
1335 }
1336
1337 /* create cdb info */
1338 AliCDBId id("TOF/Calib/Problematic", startRun, endRun);
1339 AliCDBMetaData *md = new AliCDBMetaData();
1340 md->SetResponsible("Roberto Preghenella");
1341 md->SetComment("Problematic");
1342 md->SetAliRootVersion(gSystem->Getenv("ARVERSION"));
1343 md->SetBeamPeriod(0);
1344
1345 /* put object in cdb */
1346 AliCDBManager *cdb = AliCDBManager::Instance();
1347 cdb->SetDefaultStorage(dbString);
1348 cdb->GetDefaultStorage()->Put(obj, id, md);
1349}
1350
1351//_____________________________________________________________
1352
1353Int_t
1354fitPeak(TF1 *f, TH1D *h, Float_t nSigmaMin, Float_t nSigmaMax)
1355{
1356 /*
1357 * fitPeak
1358 */
1359
1360 /* rebin if needed */
1361 Double_t ent = h->GetEntries();
1362 // while (h->GetBinContent(h->GetMaximumBin()) < 0.1 * h->GetEntries())
1363 // h->Rebin(2);
1364 Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
1365 Double_t fitMin = fitCent - nSigmaMin * 300.;
1366 Double_t fitMax = fitCent + nSigmaMax * 300.;
1367 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
1368 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
1369 f->SetParameter(1, fitCent);
1370 f->SetParameter(2, 150.);
1371 Int_t fitres = h->Fit(f, "WWq0", "", fitMin, fitMax);
1372 if (fitres != 0) return fitres;
1373 /* refit with better range */
1374 for (Int_t i = 0; i < 3; i++) {
1375 fitCent = f->GetParameter(1);
1376 fitMin = fitCent - nSigmaMin * f->GetParameter(2);
1377 fitMax = fitCent + nSigmaMax * f->GetParameter(2);
1378 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
1379 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
1380 fitres = h->Fit(f, "q0", "", fitMin, fitMax);
1381 if (fitres != 0) return fitres;
1382 }
1383 return fitres;
1384}
1385
1386//_____________________________________________________________
1387
1388void
1389monitor()
1390{
1391 /*
1392 * monitor
1393 */
1394
1395 Float_t evfrac = (Float_t)curev / (Float_t)nevents;
1396 Float_t eltime = stopwatch.RealTime();
1397 stopwatch.Continue();
1398 Float_t tottime = evfrac > 0. ? eltime / evfrac : eltime;
1399 Float_t speed = eltime > 0. ? (Float_t)curev / eltime : 0.;
1400 printf("%4.2f \% | elapsed: %d s | total: %d s | remaining: %d s | speed: %d ev/s\n", 100. * evfrac, eltime, tottime, tottime - eltime, speed);
1401
1402}
1403
1404//_____________________________________________________________
1405
1406Bool_t
1407updateParOfflineEntry(const Char_t *paramfilename, const Char_t *ocdbfilename, Int_t startrun, Int_t endrun = AliCDBRunRange::Infinity(), const Char_t *dbString = "alien://?folder=/alice/cern.ch/user/r/rpreghen/OCDB?user=rpreghen?se=ALICE::CERN::SE")
1408{
1409
1410 TFile *paramfile = TFile::Open(paramfilename);
1411 if (!paramfile || !paramfile->IsOpen()) {
1412 return kTRUE;
1413 }
1414 TH1D *hshiftparams = (TH1D *)paramfile->Get("hChannelParam_mean");
1415 if (!hshiftparams) {
1416 return kTRUE;
1417 }
1418 TFile *ocdbfile = TFile::Open(ocdbfilename);
1419 if (!ocdbfile || !ocdbfile->IsOpen()) {
1420 return kTRUE;
1421 }
1422 AliCDBEntry *cdbe = (AliCDBEntry *)ocdbfile->Get("AliCDBEntry");
1423 if (!cdbe) {
1424 return kTRUE;
1425 }
1426 TObjArray *currentpararray = (TObjArray *)cdbe->GetObject();
1427 if (!currentpararray) {
1428 return kTRUE;
1429 }
1430
1431 /* create calib structure */
1432 AliTOFcalib *tofcalib = new AliTOFcalib();
1433 tofcalib->CreateCalArrays();
1434 TObjArray *newpararray = (TObjArray *) tofcalib->GetTOFCalArrayOffline();
1435 AliTOFChannelOffline *currentparchannel, *newparchannel;
1436 Float_t currentslewpar, newslewpar;
1437 /* loop over channels */
1438 for (Int_t idx = 0; idx < 157248; idx++) {
1439 currentparchannel = (AliTOFChannelOffline *)currentpararray->At(idx);
1440 newparchannel = (AliTOFChannelOffline *)newpararray->At(idx);
1441 /* loop over slew params */
1442 for (Int_t ipar = 0; ipar < 6; ipar++) {
1443 currentslewpar = currentparchannel->GetSlewPar(ipar);
1444 if (ipar == 0) /* add shift to par[0] */
1445 newslewpar = currentslewpar + 1.e-3 * hshiftparams->GetBinContent(idx + 1);
1446 else /* leave other params unchanged */
1447 newslewpar = currentslewpar;
1448 newparchannel->SetSlewPar(ipar, newslewpar);
1449 }
1450 }
1451
1452 AliCDBManager *man = AliCDBManager::Instance();
1453 man->SetDefaultStorage("local://$HOME/OCDB");
1454 man->SetSpecificStorage("TOF/Calib/ParOffline", dbString);
1455 tofcalib->WriteParOfflineOnCDB("TOF/Calib", "valid", startrun, endrun);
1456
1457}
1458
1459//_____________________________________________________________
1460
1461Bool_t
1462updateRunParamsEntry(const Char_t *paramfilename, Int_t startrun, Int_t endrun = AliCDBRunRange::Infinity(), const Char_t *dbString = "alien://?folder=/alice/cern.ch/user/r/rpreghen/OCDB?user=rpreghen?se=ALICE::CERN::SE")
1463{
1464
1465 TFile *paramfile = TFile::Open(paramfilename);
1466 if (!paramfile || !paramfile->IsOpen()) {
1467 return kTRUE;
1468 }
1469 TH1D *htimezerofill = (TH1D *)paramfile->Get("hTimeZeroFill_mean");
1470 if (!htimezerofill) {
1471 return kTRUE;
1472 }
1473
1474 /* get t0-fill data */
1475 UInt_t timestamp[MAXPOINTS];
1476 Float_t t0[MAXPOINTS];
1477 Float_t t0spread[MAXPOINTS];
1478 Float_t tofreso[MAXPOINTS];
1479 Int_t npoints = 0;
1480 for (Int_t ibin = 0; ibin < htimezerofill->GetNbinsX(); ibin++) {
1481 if (htimezerofill->GetBinError(ibin + 1) == 0.) continue;
1482 timestamp[npoints] = (UInt_t)htimezerofill->GetBinCenter(ibin + 1);
1483 t0[npoints] = htimezerofill->GetBinContent(ibin + 1);
1484 t0spread[npoints] = -1.;
1485 tofreso[npoints] = -1.;
1486 npoints++;
1487 }
1488
1489 /* setup dummy run */
1490 UInt_t run[1] = {-1};
1491 UInt_t runFirstPoint[1] = {0};
1492 UInt_t runLastPoint[1] = {npoints - 1};
1493
1494 /* create runParams object */
1495 AliTOFRunParams obj(npoints, 1);
1496 obj.SetTimestamp(timestamp);
1497 obj.SetT0(t0);
1498 obj.SetTOFResolution(tofreso);
1499 obj.SetT0Spread(t0spread);
1500 obj.SetRunNb(run);
1501 obj.SetRunFirstPoint(runFirstPoint);
1502 obj.SetRunLastPoint(runLastPoint);
1503
1504 /* install run params object in OCDB */
1505 AliCDBManager *cdb = AliCDBManager::Instance();
1506 AliCDBStorage *sto = cdb->GetStorage(dbString);
1507 if (!sto) {
1508 printf("cannot get storage %s", dbString);
1509 return kTRUE;
1510 }
1511 AliCDBId id("TOF/Calib/RunParams", startrun, endrun);
1512 AliCDBMetaData md;
1513 md.SetResponsible("Roberto Preghenella");
1514 md.SetComment("offline TOF run parameters (emergency entry, no distinction between runs)");
1515 md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
1516 md.SetBeamPeriod(0);
1517 if (!sto->Put(&obj, id, &md)) {
1518 printf("error while putting object in storage %s", dbString);
1519 return kTRUE;
1520 }
1521
1522 return kFALSE;
1523}
1524
1525//_____________________________________________________________
1526
1527inspectTRM(const Char_t *filename, Int_t trmIndex)
1528{
1529
1530 TFile *fin = TFile::Open(filename);
1531 TFile *fout = TFile::Open(Form("inspectTRM_%d.%s", trmIndex, filename), "RECREATE");
1532 TH2 *hin = (TH2 *)fin->Get("hCalibHisto_channel");
1533 TH1 *hpy;
1534
1535 for (Int_t i = 0; i < 157248; i++) {
1536 if (getIndex(kTRM, i) != trmIndex) continue;
1537 hpy = hin->ProjectionY(Form("hChannel_%d", i), i + 1, i + 1);
1538 fout->cd();
1539 hpy->Write();
1540 }
1541
1542 fout->Close();
1543}