TOF/TOFcalibAPI.C

   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 */
  14 TFile *datafile = NULL;
  15 TTree *datatree = NULL;
  16 Int_t nevents = 0;
  17 Int_t curev = 0;
  18 Int_t runNb = 0;
  19 UInt_t timestamp = 0;
  20 Float_t vertex = 0.;
  21 Float_t timezero = 0.;
  22 Int_t nhits = 0;
  23 Float_t momentum[MAXHITS];
  24 Float_t length[MAXHITS];
  25 Int_t index[MAXHITS];
  26 Float_t time[MAXHITS];
  27 Float_t tot[MAXHITS];
  28 Float_t texp[MAXHITS];
  29
  30 /* calib histos */
  31 enum EParam_t {
  32     kTRM,
  33     kFEA,
  34     kChannel,
  35     kNParams
  36 };
  37 const Char_t *paramName[kNParams] = {"trm", "fea", "channel"};
  38 Int_t paramBins[kNParams] = {720, 6552, 157248};
  39 Float_t paramMin[kNParams] = {0., 0., 0.};
  40 Float_t paramMax[kNParams] = {720., 6552., 157248.};
  41 Int_t deltatBins = 201;
  42 Float_t deltatMin = -2440. - 12.2;
  43 Float_t deltatMax = 2440. + 12.2;
  44 Int_t totBins = 400;
  45 Float_t totMin = 4.88;
  46 Float_t totMax = 24.4;
  47 UInt_t timeZeroSampling = 600; /* seconds */
  48 UInt_t timeMin = 0;
  49 UInt_t timeMax = 0;
  50 UInt_t timeBins = 0;
  51 TH2F *hFEAHitMap = NULL;
  52 TH2F *hChannelHitMap = NULL;
  53 TH2F *hTimeZeroFillHisto = NULL;
  54 TH1F *hTimeZeroFill_mean = NULL;
  55 TH1F *hTimeZeroFill_sigma = NULL;
  56 TProfile *hTimePressureHisto = NULL;
  57 TH2F *hCalibHisto[kNParams] = {NULL, NULL, NULL};
  58 TH1F *hCalibParam_mean[kNParams] = {NULL, NULL, NULL};
  59 TH1F *hCalibParam_sigma[kNParams] = {NULL, NULL, NULL};
  60 TH1F *hChannelParam_mean = NULL;
  61 TH1F *hChannelParam_sigma = NULL;
  62
  63 TH2F *hPerfHistoDeltaT = NULL;
  64 TH2F *hPerfHistoBeta = NULL;
  65
  66 AliTOFcalibHisto *calibHisto = NULL;
  67
  68 /* other */
  69 #define MAXRUNS 1000000
  70 AliLHCClockPhase *lhcClockPhase[MAXRUNS];
  71 AliDCSSensor *cavernPressure[MAXRUNS];
  72
  73 /* monitoring */
  74 TStopwatch stopwatch;
  75
  76 //_____________________________________________________________
  77
  78 Bool_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
 103 Bool_t
 104 runCalib(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
 175 Bool_t
 176 checkCalib(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
 217 Bool_t
 218 checkPerf(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
 239 Bool_t
 240 openData(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
 295 Bool_t
 296 acceptEvent(Bool_t useTimeZeroTOF = kFALSE)
 297 {
 298     /*
 299      * acceptEvent
 300      */
 301
 302     if (useTimeZeroTOF && timezero == 999999.) return kFALSE;
 303     return kTRUE;
 304 }
 305
 306 //_____________________________________________________________
 307
 308 Float_t
 309 getTimeZeroFill(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
 322 Float_t
 323 getLHCClockPhase(UInt_t ts)
 324 {
 325     /*
 326      * getLHCClockPhase
 327      */
 328
 329     return lhcClockPhase[runNb] ? 1.e3 * lhcClockPhase[runNb]->GetPhase(ts) : 0.;
 330 }
 331
 332 //_____________________________________________________________
 333
 334 Float_t
 335 getCalib(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
 347 Float_t
 348 getDeltaT(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
 362 Float_t
 363 getBeta(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
 378 Float_t
 379 getMass(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
 396 Int_t
 397 getIndex(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
 432 Int_t
 433 getHitMapXY(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
 471 void
 472 fillTimeZeroFillHisto(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
 518 void
 519 fillCalibHistos(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
 581 void
 582 fillPerfHistos(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
 653 Bool_t
 654 writeCalibHistos(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
 680 Bool_t
 681 writePerfHistos(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
 704 Bool_t
 705 writeCalibParams(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
 732 Bool_t
 733 openCalibHistos(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
 776 Bool_t
 777 openCalibParams(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
 818 Bool_t
 819 fitTimeZeroFillHisto(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
 859 Bool_t
 860 fitCommonShift(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
 899 Bool_t
 900 fitCalibHistos(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
 914 Bool_t
 915 updateChannelParams()
 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
 942 Bool_t
 943 fitCalibHisto(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
 966 Bool_t
 967 fitCalibHisto(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
1024 problematicChannels(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
1167 problematicFEAs(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
1320 mergeProblematicMaps(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
1353 Int_t
1354 fitPeak(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
1388 void
1389 monitor()
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
1406 Bool_t
1407 updateParOfflineEntry(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
1461 Bool_t
1462 updateRunParamsEntry(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
1527 inspectTRM(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 }