TOF/AliTOFAnalysisTaskCalibPass0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14 ***************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////
  17 //                                                           //
  18 //   This class provides TOF pass0/passX calibration tools   //
  19 //                                                           //
  20 ///////////////////////////////////////////////////////////////
  21
  22 #include "AliTOFAnalysisTaskCalibPass0.h"
  23 #include "AliESDtrack.h"
  24 #include "AliESDEvent.h"
  25 #include "AliPhysicsSelection.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliCDBManager.h"
  28 #include "AliGRPManager.h"
  29 #include "AliGRPObject.h"
  30 #include "AliTOFcalib.h"
  31 #include "TList.h"
  32 #include "TH2F.h"
  33 #include "TFile.h"
  34 #include "TH1D.h"
  35 #include "TObjArray.h"
  36 #include "TString.h"
  37 #include "TObjString.h"
  38 #include "TROOT.h"
  39 #include "TF1.h"
  40 #include "TTime.h"
  41 #include "TGrid.h"
  42 #include "TTimeStamp.h"
  43 #include "AliTOFRunParams.h"
  44 #include "AliCDBStorage.h"
  45 #include "AliCDBId.h"
  46 #include "AliCDBMetaData.h"
  47 #include "TSystem.h"
  48 #include "AliLog.h"
  49
  50 ClassImp(AliTOFAnalysisTaskCalibPass0)
  51
  52 //_______________________________________________________
  53
  54 const Char_t *AliTOFAnalysisTaskCalibPass0::fgkStatusCodeName[AliTOFAnalysisTaskCalibPass0::kNStatusCodes] = {
  55   "ok",
  56   "input error",
  57   "data error",
  58   "not active",
  59   "low statistics",
  60   "no measurement",
  61   "store error"
  62 };
  63
  64 const Int_t AliTOFAnalysisTaskCalibPass0::fgkMaxNumberOfPoints = 10000; // max number of points
  65 Double_t AliTOFAnalysisTaskCalibPass0::fgMinVertexIntegral = 100.;
  66 Double_t AliTOFAnalysisTaskCalibPass0::fgMinDeltatIntegral = 2000.;
  67 Double_t AliTOFAnalysisTaskCalibPass0::fgMinVertexIntegralSample = 1000.;
  68 Double_t AliTOFAnalysisTaskCalibPass0::fgMinDeltatIntegralSample = 20000.;
  69
  70 //_______________________________________________________
  71
  72 AliTOFAnalysisTaskCalibPass0::AliTOFAnalysisTaskCalibPass0() :
  73   AliAnalysisTaskSE("TOFCalib-Pass0"),
  74   fStatus(kOk),
  75   fInitFlag(kFALSE),
  76   fEventSelectionFlag(kFALSE),
  77   fVertexSelectionFlag(kFALSE),
  78   fVertexCut(10.),
  79   fRunNumber(0),
  80   fESDEvent(NULL),
  81   fEventCuts(new AliPhysicsSelection()),
  82   fTrackCuts(new AliESDtrackCuts()),
  83   fStartTime(0),
  84   fEndTime(0),
  85   fEventTime(0),
  86   fElapsedTime(0),
  87   fkVertex(NULL),
  88   fGRPManager(new AliGRPManager()),
  89   fkGRPObject(NULL),
  90   fTOFcalib(new AliTOFcalib()),
  91   fHistoList(new TList()),
  92   fHistoVertexTimestamp(NULL),
  93   fHistoDeltatTimestamp(NULL),
  94   fHistoDeltazEta(NULL),
  95   fHistoDeltazCosTheta(NULL),
  96   fHistoAcceptedTracksEtaPt(NULL),
  97   fHistoMatchedTracksEtaPt(NULL)
  98 {
  99   /*
 100    * default constructor
 101    */
 102
 103   /* define output */
 104   fHistoList->SetOwner(kTRUE);
 105   DefineOutput(1, TList::Class());
 106 }
 107
 108 //_______________________________________________________
 109
 110 AliTOFAnalysisTaskCalibPass0::~AliTOFAnalysisTaskCalibPass0()
 111 {
 112   /*
 113    * default destructor
 114    */
 115
 116   if (fHistoList) delete fHistoList;
 117 }
 118
 119 //_______________________________________________________
 120
 121 Bool_t
 122 AliTOFAnalysisTaskCalibPass0::InitRun()
 123 {
 124   /*
 125    * init run
 126    */
 127
 128   /* get ESD event */
 129   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 130   if (!fESDEvent) {
 131     AliError("cannot get ESD event");
 132     return kFALSE;
 133   }
 134   /* get run number */
 135   Int_t runNb = fESDEvent->GetRunNumber();
 136   /* check run already initialized */
 137   if (fInitFlag && fRunNumber == runNb) return kTRUE;
 138   if (fInitFlag && fRunNumber != runNb) {
 139     AliFatal(Form("run number has changed within same job and this is not allowed: %d -> %d", fRunNumber, runNb));
 140     return kTRUE;
 141   }
 142   fInitFlag = kFALSE;
 143   /* init cdb if not yet done */
 144   AliCDBManager *cdb = AliCDBManager::Instance();
 145   if (!cdb->IsDefaultStorageSet())
 146     cdb->SetDefaultStorage("raw://");
 147   if (cdb->GetRun() < 0)
 148     cdb->SetRun(runNb);
 149     /* init TOF calib */
 150   if (!fTOFcalib->Init()) {
 151     AliError("cannot init TOF calib");
 152     return kFALSE;
 153   }
 154   /* get GRP data */
 155   if (!fGRPManager->ReadGRPEntry()) {
 156     AliError("error while reading \"GRPEntry\" from OCDB");
 157     return kFALSE;
 158   }
 159   fkGRPObject = fGRPManager->GetGRPData();
 160   if (!fkGRPObject) {
 161     AliError("cannot get \"GRPData\" from GRP manager");
 162     return kFALSE;
 163   }
 164   fStartTime = fkGRPObject->GetTimeStart();
 165   fEndTime = fkGRPObject->GetTimeEnd();
 166   AliInfo(Form("got \"GRPData\": startTime=%d, endTime=%d", fStartTime, fEndTime));
 167
 168   AliInfo(Form("initialized for run %d", runNb));
 169   fInitFlag = kTRUE;
 170   fRunNumber = runNb;
 171
 172   Int_t useBPTX = fTOFcalib->GetUseLHCClockPhase() ? 1 : 0;
 173
 174   /* set histo title with run-number and start-time */
 175   fHistoVertexTimestamp->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 176   fHistoDeltatTimestamp->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 177   fHistoDeltazEta->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 178   fHistoDeltazCosTheta->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 179   fHistoAcceptedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 180   fHistoMatchedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u, BPTX: %d", fRunNumber, fStartTime, useBPTX));
 181
 182   return kTRUE;
 183 }
 184
 185 //_______________________________________________________
 186
 187 Bool_t
 188 AliTOFAnalysisTaskCalibPass0::InitEvent()
 189 {
 190   /*
 191    * init event
 192    */
 193
 194   /* get ESD event */
 195   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 196   if (!fESDEvent) return kFALSE;
 197   /* set event time and elapsed time */
 198   fEventTime = fESDEvent->GetTimeStamp();
 199   fElapsedTime = fESDEvent->GetTimeStamp() - fStartTime;
 200   /* event selection */
 201   if (fEventSelectionFlag && !fEventCuts->IsCollisionCandidate(fESDEvent)) return kFALSE;
 202   /* vertex selection */
 203   fkVertex = fESDEvent->GetPrimaryVertexTracks();
 204   if (fVertexSelectionFlag && fkVertex->GetNContributors() < 1) {
 205     fkVertex = fESDEvent->GetPrimaryVertexSPD();
 206     if (fkVertex->GetNContributors() < 1) return kFALSE;
 207   }
 208   if (fVertexSelectionFlag && TMath::Abs(fkVertex->GetZ()) > fVertexCut) return kFALSE;
 209   /* calibrate ESD if requested */
 210   fTOFcalib->CalibrateESD(fESDEvent);
 211
 212   return kTRUE;
 213 }
 214
 215 //_______________________________________________________
 216
 217 Bool_t
 218 AliTOFAnalysisTaskCalibPass0::HasTOFMeasurement(const AliESDtrack *track) const
 219 {
 220   /*
 221    * has TOF measurement
 222    */
 223
 224   /* check TOF status flags */
 225   if (!(track->GetStatus() & AliESDtrack::kTOFout) ||
 226       !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE;
 227   /* check integrated length */
 228   if (track->GetIntegratedLength() < 350.) return kFALSE;
 229
 230   /* TOF measurement ok */
 231   return kTRUE;
 232 }
 233
 234 //_______________________________________________________
 235
 236 void
 237 AliTOFAnalysisTaskCalibPass0::UserCreateOutputObjects()
 238 {
 239   /*
 240    * user create output objects
 241    */
 242
 243   /* time binning */
 244   Int_t timeBins = 288;
 245   Float_t timeMin = 0.;
 246   Float_t timeMax = 24. * 3600.;
 247   /* vertex binning */
 248   Int_t vertexBins = 200;
 249   Float_t vertexMin = -50.;
 250   Float_t vertexMax = 50.;
 251   /* deltat binning */
 252   Int_t deltatBins = 2000;
 253   Float_t deltatMin = -24400.;
 254   Float_t deltatMax = 24400.;
 255   /* deltaz binning */
 256   Int_t deltazBins = 200;
 257   Float_t deltazMin = -10.;
 258   Float_t deltazMax = 10.;
 259   /* eta binning */
 260   Int_t etaBins = 200;
 261   Float_t etaMin = -1.;
 262   Float_t etaMax = 1.;
 263   /* p binning */
 264   Int_t pBins = 500;
 265   Float_t pMin = 0.;
 266   Float_t pMax = 5.;
 267
 268   fHistoVertexTimestamp = new TH2F("hHistoVertexTimestamp", "Vertex position;elapsed time (s);z (cm);", timeBins, timeMin, timeMax, vertexBins, vertexMin, vertexMax);
 269   fHistoList->Add(fHistoVertexTimestamp);
 270
 271   fHistoDeltatTimestamp = new TH2F("hHistoDeltatTimestamp", "Global time shift (T0-fill);elapsed time (s);t - t_{exp}^{(#pi)} (ps);", timeBins, timeMin, timeMax, deltatBins, deltatMin, deltatMax);
 272   fHistoList->Add(fHistoDeltatTimestamp);
 273
 274   fHistoDeltazEta = new TH2F("hHistoDeltazEta", "Matching residuals (longitudinal);#eta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
 275   fHistoList->Add(fHistoDeltazEta);
 276
 277   fHistoDeltazCosTheta = new TH2F("hHistoDeltazCosTheta", "Matching residuals (longitudinal);cos #theta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
 278   fHistoList->Add(fHistoDeltazCosTheta);
 279
 280   fHistoAcceptedTracksEtaPt = new TH2F("hHistoAcceptedTracksEtaPt", "Accepted tracks;#eta;#p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
 281   fHistoList->Add(fHistoAcceptedTracksEtaPt);
 282
 283   fHistoMatchedTracksEtaPt = new TH2F("hHistoMatchedTracksEtaPt", "Matched tracks;#eta;p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
 284   fHistoList->Add(fHistoMatchedTracksEtaPt);
 285
 286   /* post data */
 287   PostData(1, fHistoList);
 288 }
 289
 290 //_______________________________________________________
 291
 292 void
 293 AliTOFAnalysisTaskCalibPass0::UserExec(Option_t *)
 294 {
 295   /*
 296    * user exec
 297    */
 298
 299   /* init run */
 300   if (!InitRun()) return;
 301   /* init event */
 302   if (!InitEvent()) return;
 303
 304   /*** ACCEPTED EVENT ***/
 305
 306   /* fill vertex histo */
 307   fHistoVertexTimestamp->Fill(fElapsedTime, fkVertex->GetZ());
 308
 309   /* loop over ESD tracks */
 310   Int_t nTracks = fESDEvent->GetNumberOfTracks();
 311   AliESDtrack *track;
 312   Double_t eta, costheta, pt, time, timei[AliPID::kSPECIES], deltat, deltaz;
 313   for (Int_t itrk = 0; itrk < nTracks; itrk++) {
 314     /* get track */
 315     track = fESDEvent->GetTrack(itrk);
 316     if (!track) continue;
 317     /* check accept track */
 318     if (!fTrackCuts->AcceptTrack(track)) continue;
 319     /* get track info */
 320     eta = track->Eta();
 321     costheta = TMath::Cos(track->Theta());
 322     pt = track->Pt();
 323     /* fill accepted tracks histo */
 324     fHistoAcceptedTracksEtaPt->Fill(eta, pt);
 325     /* check TOF measurement */
 326     if (!HasTOFMeasurement(track)) continue;
 327
 328     /*** ACCEPTED TRACK WITH TOF MEASUREMENT ***/
 329
 330     /* fill matched tracks histo */
 331     fHistoMatchedTracksEtaPt->Fill(eta, pt);
 332     /* get TOF info */
 333     time = track->GetTOFsignal();
 334     track->GetIntegratedTimes(timei);
 335     deltat = time - timei[AliPID::kPion];
 336     deltaz = track->GetTOFsignalDz();
 337
 338     /* fill histos */
 339     fHistoDeltatTimestamp->Fill(fElapsedTime, deltat);
 340     fHistoDeltazEta->Fill(eta, deltaz);
 341     fHistoDeltazCosTheta->Fill(costheta, deltaz);
 342
 343   } /* end of loop over ESD tracks */
 344
 345   /* post data */
 346   PostData(1, fHistoList);
 347
 348 }
 349
 350 //_______________________________________________________
 351
 352 Int_t
 353 AliTOFAnalysisTaskCalibPass0::GetStatus()
 354 {
 355   /*
 356    * get status
 357    */
 358
 359   switch (fStatus) {
 360
 361     /* OK, return zero */
 362   case kOk:
 363     return 0;
 364     break;
 365
 366     /* non-fatal error, return negative status */
 367   case kNotActive:
 368   case kLowStatistics:
 369   case kNoMeasurement:
 370     return -fStatus;
 371     break;
 372
 373     /* fatal error, return positive status */
 374   case kInputError:
 375   case kDataError:
 376   case kStoreError:
 377     return fStatus;
 378     break;
 379
 380     /* anything else, return negative large number */
 381   default:
 382     return -999;
 383     break;
 384   }
 385
 386   /* should never arrive here, anyway return negative large number */
 387   return -999;
 388 }
 389
 390 //_______________________________________________________
 391
 392 void
 393 AliTOFAnalysisTaskCalibPass0::PrintStatus()
 394 {
 395   /*
 396    * print status
 397    */
 398
 399   Int_t status = GetStatus();
 400   if (status == 0) {
 401     AliInfo(Form("TOF calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 402   }
 403   else if (status > 0) {
 404     AliInfo(Form("TOF calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 405   }
 406   else if (status < 0) {
 407     AliInfo(Form("TOF calibration failed (expected): %s (status=%d)", fgkStatusCodeName[fStatus], status));
 408   }
 409
 410 }
 411
 412 //_______________________________________________________
 413
 414 Bool_t
 415 AliTOFAnalysisTaskCalibPass0::ProcessOutput(const Char_t *filename, const Char_t *dbString)
 416 {
 417   /*
 418    * process output
 419    */
 420
 421   Int_t ret = DoProcessOutput(filename, dbString);
 422   PrintStatus();
 423   return ret;
 424 }
 425
 426 //_______________________________________________________
 427
 428 Bool_t
 429 AliTOFAnalysisTaskCalibPass0::DoProcessOutput(const Char_t *filename, const Char_t *dbString)
 430 {
 431   /*
 432    * do process output
 433    */
 434
 435   /* reset status to OK */
 436   fStatus = kOk;
 437
 438   /* open file */
 439   TFile *file = TFile::Open(filename);
 440   if (!file || !file->IsOpen()) {
 441     AliError(Form("cannot open output file %s", filename));
 442     fStatus = kInputError;
 443     return kFALSE;
 444   }
 445   /* get histograms */
 446   TList *list = (TList *)file->Get("TOFHistos");
 447   TH2F *histoVertexTimestamp = NULL;
 448   TH2F *histoDeltatTimestamp = NULL;
 449   TH2F *histoDeltazEta = NULL;
 450   TH2F *histoDeltazCosTheta = NULL;
 451   TH2F *histoAcceptedTracksEtaPt = NULL;
 452   TH2F *histoMatchedTracksEtaPt = NULL;
 453   if (list) {
 454     AliInfo(Form("getting histograms from \"Histos\" list from file %s", filename));
 455     histoVertexTimestamp = (TH2F *)list->FindObject("hHistoVertexTimestamp");
 456     histoDeltatTimestamp = (TH2F *)list->FindObject("hHistoDeltatTimestamp");
 457     histoDeltazEta = (TH2F *)list->FindObject("hHistoDeltazEta");
 458     histoDeltazCosTheta = (TH2F *)list->FindObject("hHistoDeltazCosTheta");
 459     histoAcceptedTracksEtaPt = (TH2F *)list->FindObject("hHistoAcceptedTracksEtaPt");
 460     histoMatchedTracksEtaPt = (TH2F *)list->FindObject("hHistoMatchedTracksEtaPt");
 461   }
 462   else {
 463     AliInfo(Form("getting histograms directly from file %s", filename));
 464     histoVertexTimestamp = (TH2F *)file->Get("hHistoVertexTimestamp");
 465     histoDeltatTimestamp = (TH2F *)file->Get("hHistoDeltatTimestamp");
 466     histoDeltazEta = (TH2F *)file->Get("hHistoDeltazEta");
 467     histoDeltazCosTheta = (TH2F *)file->Get("hHistoDeltazCosTheta");
 468     histoAcceptedTracksEtaPt = (TH2F *)file->Get("hHistoAcceptedTracksEtaPt");
 469     histoMatchedTracksEtaPt = (TH2F *)file->Get("hHistoMatchedTracksEtaPt");
 470   }
 471   /* check histos */
 472   if (!histoVertexTimestamp) {
 473     AliError(Form("cannot get \"hHistoVertexTimestamp\" object from file %s", filename));
 474     fStatus = kInputError;
 475     return kFALSE;
 476   }
 477   if (!histoDeltatTimestamp) {
 478     AliError(Form("cannot get \"hHistoDeltatTimestamp\" object from file %s", filename));
 479     fStatus = kInputError;
 480     return kFALSE;
 481   }
 482   if (!histoDeltazEta) {
 483     AliError(Form("cannot get \"hHistoDeltazEta\" object from file %s", filename));
 484     fStatus = kInputError;
 485     return kFALSE;
 486   }
 487   if (!histoDeltazCosTheta) {
 488     AliError(Form("cannot get \"hHistoDeltazCosTheta\" object from file %s", filename));
 489     fStatus = kInputError;
 490     return kFALSE;
 491   }
 492   if (!histoAcceptedTracksEtaPt) {
 493     AliError(Form("cannot get \"hHistoAccptedTracksEtaPt\" object from file %s", filename));
 494     fStatus = kInputError;
 495     return kFALSE;
 496   }
 497   if (!histoMatchedTracksEtaPt) {
 498     AliError(Form("cannot get \"hHistoMatchedTracksEtaPt\" object from file %s", filename));
 499     fStatus = kInputError;
 500     return kFALSE;
 501   }
 502
 503   /* check matching performance */
 504   if (!CheckMatchingPerformance(histoDeltazEta, histoAcceptedTracksEtaPt, histoMatchedTracksEtaPt)) {
 505     AliError("error while checking matching efficiency");
 506     return kFALSE;
 507   }
 508   /* calibrate and store */
 509   if (!CalibrateAndStore(histoVertexTimestamp, histoDeltatTimestamp, dbString)) {
 510     AliError("error while calibrating and storing");
 511     return kFALSE;
 512   }
 513
 514   /* success */
 515   return kTRUE;
 516 }
 517
 518 //_______________________________________________________
 519
 520 Bool_t
 521 AliTOFAnalysisTaskCalibPass0::CheckMatchingPerformance(const TH2F *histoDeltazCosTheta, const TH2F *histoAcceptedTracksEtaPt, const TH2F *histoMatchedTracksEtaPt) const
 522 {
 523   /*
 524    * check matching performance
 525    */
 526
 527   /* check pointers */
 528   if (!histoDeltazCosTheta || !histoAcceptedTracksEtaPt || !histoMatchedTracksEtaPt)
 529     return kFALSE;
 530   /* dummy for the time being */
 531   return kTRUE;
 532 }
 533
 534 //_______________________________________________________
 535
 536 Bool_t
 537 AliTOFAnalysisTaskCalibPass0::CalibrateAndStore(TH2F *histoVertexTimestamp, TH2F *histoDeltatTimestamp, const Char_t *dbString)
 538 {
 539   /*
 540    * calibrate and store
 541    */
 542
 543   /* check pointers */
 544   if (!histoVertexTimestamp || !histoDeltatTimestamp)
 545     return kFALSE;
 546
 547   /*** GET RUN-NUMBER AND START-TIMESTAMP ***/
 548
 549   TString str;
 550   TObjArray *strarr = NULL;
 551   TObjString *ostr = NULL;
 552   str = histoVertexTimestamp->GetTitle();
 553   strarr = str.Tokenize(",");
 554   if (!strarr) {
 555     AliError("problems whith tokenize histogram title");
 556     fStatus = kDataError;
 557     return kFALSE;
 558   }
 559
 560   /* run number */
 561   ostr = (TObjString *)strarr->At(0);
 562   if (!ostr) {
 563     AliError("problems while getting run number from histogram title");
 564     fStatus = kDataError;
 565     return kFALSE;
 566   }
 567   str = ostr->GetString();
 568   if (!str.BeginsWith("run:")) {
 569     AliError("problems while getting run number from histogram title");
 570     fStatus = kDataError;
 571     return kFALSE;
 572   }
 573   str.Remove(0, 5);
 574   Int_t runNb = atoi(str.Data());
 575   if (runNb <= 0) {
 576     AliError(Form("bad run number: %d", runNb));
 577     fStatus = kDataError;
 578     return kFALSE;
 579   }
 580   AliInfo(Form("got run number: %d", runNb));
 581
 582   /* start timestamp */
 583   ostr = (TObjString *)strarr->At(1);
 584   if (!ostr) {
 585     AliError("problems while getting start timestamp from histogram title");
 586     fStatus = kDataError;
 587     return kFALSE;
 588   }
 589   str = ostr->GetString();
 590   str.Remove(0, 1); /* remove empty space at the beginning */
 591   if (!str.BeginsWith("startTimestamp:")) {
 592     AliError("problems while getting start timestamp from histogram title");
 593     fStatus = kDataError;
 594     return kFALSE;
 595   }
 596   str.Remove(0, 16);
 597   UInt_t startTimestamp = atoi(str.Data());
 598   if (startTimestamp <= 0) {
 599     AliError(Form("bad start timestamp: %d", startTimestamp));
 600     fStatus = kDataError;
 601     return kFALSE;
 602   }
 603   TTimeStamp ts = startTimestamp;
 604   AliInfo(Form("got start timestamp: %d (%s)", startTimestamp, ts.AsString()));
 605
 606   /* BPTX */
 607   ostr = (TObjString *)strarr->At(2);
 608   if (!ostr) {
 609     AliError("problems while getting BPTX from histogram title");
 610     fStatus = kDataError;
 611     return kFALSE;
 612   }
 613   str = ostr->GetString();
 614   str.Remove(0, 1); /* remove empty space at the beginning */
 615   if (!str.BeginsWith("BPTX:")) {
 616     AliError("problems while getting BPTX from histogram title");
 617     fStatus = kDataError;
 618     return kFALSE;
 619   }
 620   str.Remove(0, 6);
 621   Bool_t useBPTX = atoi(str.Data());
 622   AliInfo(Form("got BPTX: %d", useBPTX));
 623
 624   /*** CALIBRATION STAGE ***/
 625
 626   /* get fit function */
 627   TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
 628
 629   /* projection-x */
 630   TH1D *histoVertexTimestamppx = histoVertexTimestamp->ProjectionX("histoVertexTimestamppx");
 631   TH1D *histoDeltatTimestamppx = histoDeltatTimestamp->ProjectionX("histoDeltatTimestamppx");
 632
 633   /* check statistics */
 634   if (histoVertexTimestamppx->Integral() < fgMinVertexIntegral ||
 635       histoDeltatTimestamppx->Integral() < fgMinDeltatIntegral) {
 636     fStatus = kLowStatistics;
 637     return kFALSE;
 638   }
 639
 640   /* define mix and max time bin */
 641   Int_t minBin = histoVertexTimestamppx->FindFirstBinAbove(0);
 642   Int_t maxBin = histoVertexTimestamppx->FindLastBinAbove(0);
 643   Float_t minTime = histoVertexTimestamppx->GetBinLowEdge(minBin);
 644   Float_t maxTime = histoVertexTimestamppx->GetBinLowEdge(maxBin + 1);
 645   AliInfo(Form("min/max time defined: %d < t < %d s [%d, %d]", (Int_t)minTime, (Int_t)maxTime, minBin, maxBin));
 646
 647   /* loop over time bins */
 648   Int_t nPoints = 0;
 649   Float_t time[fgkMaxNumberOfPoints], timeerr[fgkMaxNumberOfPoints];
 650   Float_t vertexMean[fgkMaxNumberOfPoints], vertexMeanerr[fgkMaxNumberOfPoints];
 651   Float_t vertexSigma[fgkMaxNumberOfPoints], vertexSigmaerr[fgkMaxNumberOfPoints];
 652   Float_t timeZeroMean[fgkMaxNumberOfPoints], timeZeroMeanerr[fgkMaxNumberOfPoints];
 653   Float_t timeZeroSigma[fgkMaxNumberOfPoints], timeZeroSigmaerr[fgkMaxNumberOfPoints];
 654   Float_t averageTimeZero = 0.;
 655   Float_t averageTimeSigma = 0.;
 656   Float_t averageVertexSpread = 0.;
 657   for (Int_t ibin = minBin; ibin <= maxBin; ibin++) {
 658
 659     /* define time window */
 660     Int_t startBin = ibin;
 661     Int_t endBin = ibin;
 662     while(histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegralSample ||
 663           histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegralSample) {
 664       if (endBin < maxBin) endBin++;
 665       else if (startBin > minBin) startBin--;
 666       else break;
 667     }
 668     if (histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegral ||
 669         histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegral) continue;
 670     Float_t startTime = histoVertexTimestamppx->GetBinLowEdge(startBin);
 671     Float_t endTime = histoVertexTimestamppx->GetBinLowEdge(endBin + 1);
 672     Float_t vertexIntegral = histoVertexTimestamppx->Integral(startBin, endBin);
 673     Float_t deltatIntegral = histoDeltatTimestamppx->Integral(startBin, endBin);
 674     AliInfo(Form("time window defined: %d < t < %d s [%d, %d]: %d vertices, %d tracks", (Int_t)startTime, (Int_t)endTime, startBin, endBin, (Int_t)vertexIntegral, (Int_t)deltatIntegral));
 675
 676     /* projection-y */
 677     TH1D *histoVertexTimestamppy = histoVertexTimestamp->ProjectionY("histoVertexTimestamppy", startBin, endBin);
 678     TH1D *histoDeltatTimestamppy = histoDeltatTimestamp->ProjectionY("histoDeltatTimestamppy", startBin, endBin);
 679
 680     /* average time */
 681     histoVertexTimestamppx->GetXaxis()->SetRange(startBin, endBin);
 682     time[nPoints] = histoVertexTimestamppx->GetMean();
 683     timeerr[nPoints] = histoVertexTimestamppx->GetMeanError();
 684
 685     /* fit vertex */
 686     if (FitPeak(fitFunc, histoVertexTimestamppy, 10., 3., 3.) != 0) {
 687       AliError("troubles fitting vertex, skip");
 688       delete histoVertexTimestamppy;
 689       delete histoDeltatTimestamppy;
 690       continue;
 691     }
 692     vertexMean[nPoints] = fitFunc->GetParameter(1);
 693     vertexMeanerr[nPoints] = fitFunc->GetParError(1);
 694     vertexSigma[nPoints] = fitFunc->GetParameter(2);
 695     vertexSigmaerr[nPoints] = fitFunc->GetParError(2);
 696     averageVertexSpread += fitFunc->GetParameter(2);
 697
 698     /* fit time-zero */
 699     if (FitPeak(fitFunc, histoDeltatTimestamppy, 500., 2., 1.) != 0) {
 700       AliError("troubles fitting time-zero TRACKS, skip");
 701       delete histoVertexTimestamppy;
 702       delete histoDeltatTimestamppy;
 703       continue;
 704     }
 705     timeZeroMean[nPoints] = fitFunc->GetParameter(1);
 706     timeZeroMeanerr[nPoints] = fitFunc->GetParError(1);
 707     timeZeroSigma[nPoints] = fitFunc->GetParameter(2);
 708     timeZeroSigmaerr[nPoints] = fitFunc->GetParError(2);
 709     averageTimeZero += fitFunc->GetParameter(1);
 710     averageTimeSigma += fitFunc->GetParameter(2);
 711
 712     /* delete projection-y */
 713     delete histoVertexTimestamppy;
 714     delete histoDeltatTimestamppy;
 715
 716     /* increment n points */
 717     nPoints++;
 718
 719     /* set current bin */
 720     ibin = endBin;
 721
 722   } /* end of loop over time bins */
 723
 724   /* delete projection-x */
 725   delete histoVertexTimestamppx;
 726   delete histoDeltatTimestamppx;
 727
 728   /* check points */
 729   if (nPoints <= 0) {
 730     AliError("no measurement available, quit");
 731     fStatus = kNoMeasurement;
 732     return kFALSE;
 733   }
 734   AliInfo(Form("average time-zero  = %f", averageTimeZero / nPoints));
 735   AliInfo(Form("average time-sigma = %f", averageTimeSigma / nPoints));
 736   AliInfo(Form("average v-spread   = %f", averageVertexSpread / nPoints));
 737
 738   /*** CREATE RUN PARAMS OBJECT ***/
 739
 740 #if 0
 741   /* get start time from GRP */
 742   TGrid::Connect("alien://", 0, 0, "t");
 743   AliCDBManager *cdb = AliCDBManager::Instance();
 744   cdb->SetDefaultStorage("raw://");
 745   cdb->SetRun(runNb);
 746   AliGRPManager grp;
 747   if (!grp.ReadGRPEntry()) {
 748     AliError("error while reading GRP entry");
 749     return kFALSE;
 750   }
 751   UInt_t startTimestamp = grp.GetGRPData()->GetTimeStart();
 752   TTimeStamp ts;
 753   ts = startTimestamp;
 754   AliInfo(Form("got start time from GRP: %s", ts.AsString()));
 755 #endif
 756
 757   /* create arrays */
 758   UInt_t timestamp[fgkMaxNumberOfPoints];
 759   Float_t t0[fgkMaxNumberOfPoints];
 760   Float_t tofReso[fgkMaxNumberOfPoints];
 761   Float_t t0Spread[fgkMaxNumberOfPoints];
 762   for (Int_t ipoint = 0; ipoint < nPoints; ipoint++) {
 763     timestamp[ipoint] = (UInt_t)time[ipoint] + startTimestamp;
 764     t0[ipoint] = timeZeroMean[ipoint];
 765     tofReso[ipoint] = timeZeroSigma[ipoint];
 766     t0Spread[ipoint] = vertexSigma[ipoint] / 2.99792457999999984e-02;
 767   }
 768   UInt_t run[1] = {runNb};
 769   UInt_t runFirstPoint[1] = {0};
 770   UInt_t runLastPoint[1] = {nPoints - 1};
 771
 772   /* create run params object */
 773   AliTOFRunParams obj(nPoints, 1);
 774   AliInfo(Form("create run params object for run %d with %d points", runNb, nPoints));
 775   obj.SetTimestamp(timestamp);
 776   obj.SetT0(t0);
 777   obj.SetTOFResolution(tofReso);
 778   obj.SetT0Spread(t0Spread);
 779   obj.SetRunNb(run);
 780   obj.SetRunFirstPoint(runFirstPoint);
 781   obj.SetRunLastPoint(runLastPoint);
 782   obj.SetUseLHCClockPhase(useBPTX);
 783
 784   /*** CREATE OCDB ENTRY ***/
 785
 786   if (!dbString) {
 787     AliError("cannot store object because of NULL string");
 788     fStatus = kStoreError;
 789     return kFALSE;
 790   }
 791
 792   /* install run params object in OCDB */
 793   AliCDBManager *cdb = AliCDBManager::Instance();
 794   AliCDBStorage *sto = cdb->GetStorage(dbString);
 795   if (!sto) {
 796     AliError(Form("cannot get storage %s", dbString));
 797     fStatus = kStoreError;
 798     return kFALSE;
 799   }
 800   AliCDBId id("TOF/Calib/RunParams", runNb, runNb);
 801   AliCDBMetaData md;
 802   md.SetResponsible("Roberto Preghenella");
 803   md.SetComment("offline TOF run parameters");
 804   md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
 805   md.SetBeamPeriod(0);
 806   if (!sto->Put(&obj, id, &md)) {
 807     fStatus = kStoreError;
 808     AliError(Form("error while putting object in storage %s", dbString));
 809     return kFALSE;
 810   }
 811
 812   /* success */
 813   return kTRUE;
 814 }
 815
 816 //_____________________________________________________________
 817
 818 Int_t
 819 AliTOFAnalysisTaskCalibPass0::FitPeak(TF1 *fitFunc, TH1D *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
 820 {
 821   /*
 822    * fit peak
 823    */
 824
 825   Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
 826   Double_t fitMin = fitCent - nSigmaMin * startSigma;
 827   Double_t fitMax = fitCent + nSigmaMax * startSigma;
 828   if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 829   if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 830   fitFunc->SetParameter(1, fitCent);
 831   fitFunc->SetParameter(2, startSigma);
 832   Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
 833   if (fitres != 0) return fitres;
 834   /* refit with better range */
 835   for (Int_t i = 0; i < 3; i++) {
 836     fitCent = fitFunc->GetParameter(1);
 837     fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
 838     fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
 839     if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 840     if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 841     fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
 842     if (fitres != 0) return fitres;
 843   }
 844   return fitres;
 845 }