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 Bool_t
 393 AliTOFAnalysisTaskCalibPass0::ProcessOutput(const Char_t *filename, const Char_t *dbString)
 394 {
 395   /*
 396    * process output
 397    */
 398
 399   Int_t ret = DoProcessOutput(filename, dbString);
 400   Int_t status = GetStatus();
 401   if (status == 0) {
 402     AliInfo(Form("TOF calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 403   }
 404   else if (status > 0) {
 405     AliInfo(Form("TOF calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 406   }
 407   else if (status < 0) {
 408     AliInfo(Form("TOF calibration failed (expected): %s (status=%d)", fgkStatusCodeName[fStatus], status));
 409   }
 410
 411   return ret;
 412 }
 413
 414 //_______________________________________________________
 415
 416 Bool_t
 417 AliTOFAnalysisTaskCalibPass0::DoProcessOutput(const Char_t *filename, const Char_t *dbString)
 418 {
 419   /*
 420    * do process output
 421    */
 422
 423   /* reset status to OK */
 424   fStatus = kOk;
 425
 426   /* open file */
 427   TFile *file = TFile::Open(filename);
 428   if (!file || !file->IsOpen()) {
 429     AliError(Form("cannot open output file %s", filename));
 430     fStatus = kInputError;
 431     return kFALSE;
 432   }
 433   /* get histograms */
 434   TList *list = (TList *)file->Get("TOFHistos");
 435   TH2F *histoVertexTimestamp = NULL;
 436   TH2F *histoDeltatTimestamp = NULL;
 437   TH2F *histoDeltazEta = NULL;
 438   TH2F *histoDeltazCosTheta = NULL;
 439   TH2F *histoAcceptedTracksEtaPt = NULL;
 440   TH2F *histoMatchedTracksEtaPt = NULL;
 441   if (list) {
 442     AliInfo(Form("getting histograms from \"Histos\" list from file %s", filename));
 443     histoVertexTimestamp = (TH2F *)list->FindObject("hHistoVertexTimestamp");
 444     histoDeltatTimestamp = (TH2F *)list->FindObject("hHistoDeltatTimestamp");
 445     histoDeltazEta = (TH2F *)list->FindObject("hHistoDeltazEta");
 446     histoDeltazCosTheta = (TH2F *)list->FindObject("hHistoDeltazCosTheta");
 447     histoAcceptedTracksEtaPt = (TH2F *)list->FindObject("hHistoAcceptedTracksEtaPt");
 448     histoMatchedTracksEtaPt = (TH2F *)list->FindObject("hHistoMatchedTracksEtaPt");
 449   }
 450   else {
 451     AliInfo(Form("getting histograms directly from file %s", filename));
 452     histoVertexTimestamp = (TH2F *)file->Get("hHistoVertexTimestamp");
 453     histoDeltatTimestamp = (TH2F *)file->Get("hHistoDeltatTimestamp");
 454     histoDeltazEta = (TH2F *)file->Get("hHistoDeltazEta");
 455     histoDeltazCosTheta = (TH2F *)file->Get("hHistoDeltazCosTheta");
 456     histoAcceptedTracksEtaPt = (TH2F *)file->Get("hHistoAcceptedTracksEtaPt");
 457     histoMatchedTracksEtaPt = (TH2F *)file->Get("hHistoMatchedTracksEtaPt");
 458   }
 459   /* check histos */
 460   if (!histoVertexTimestamp) {
 461     AliError(Form("cannot get \"hHistoVertexTimestamp\" object from file %s", filename));
 462     fStatus = kInputError;
 463     return kFALSE;
 464   }
 465   if (!histoDeltatTimestamp) {
 466     AliError(Form("cannot get \"hHistoDeltatTimestamp\" object from file %s", filename));
 467     fStatus = kInputError;
 468     return kFALSE;
 469   }
 470   if (!histoDeltazEta) {
 471     AliError(Form("cannot get \"hHistoDeltazEta\" object from file %s", filename));
 472     fStatus = kInputError;
 473     return kFALSE;
 474   }
 475   if (!histoDeltazCosTheta) {
 476     AliError(Form("cannot get \"hHistoDeltazCosTheta\" object from file %s", filename));
 477     fStatus = kInputError;
 478     return kFALSE;
 479   }
 480   if (!histoAcceptedTracksEtaPt) {
 481     AliError(Form("cannot get \"hHistoAccptedTracksEtaPt\" object from file %s", filename));
 482     fStatus = kInputError;
 483     return kFALSE;
 484   }
 485   if (!histoMatchedTracksEtaPt) {
 486     AliError(Form("cannot get \"hHistoMatchedTracksEtaPt\" object from file %s", filename));
 487     fStatus = kInputError;
 488     return kFALSE;
 489   }
 490
 491   /* check matching performance */
 492   if (!CheckMatchingPerformance(histoDeltazEta, histoAcceptedTracksEtaPt, histoMatchedTracksEtaPt)) {
 493     AliError("error while checking matching efficiency");
 494     return kFALSE;
 495   }
 496   /* calibrate and store */
 497   if (!CalibrateAndStore(histoVertexTimestamp, histoDeltatTimestamp, dbString)) {
 498     AliError("error while calibrating and storing");
 499     return kFALSE;
 500   }
 501
 502   /* success */
 503   return kTRUE;
 504 }
 505
 506 //_______________________________________________________
 507
 508 Bool_t
 509 AliTOFAnalysisTaskCalibPass0::CheckMatchingPerformance(const TH2F *histoDeltazCosTheta, const TH2F *histoAcceptedTracksEtaPt, const TH2F *histoMatchedTracksEtaPt) const
 510 {
 511   /*
 512    * check matching performance
 513    */
 514
 515   /* check pointers */
 516   if (!histoDeltazCosTheta || !histoAcceptedTracksEtaPt || !histoMatchedTracksEtaPt)
 517     return kFALSE;
 518   /* dummy for the time being */
 519   return kTRUE;
 520 }
 521
 522 //_______________________________________________________
 523
 524 Bool_t
 525 AliTOFAnalysisTaskCalibPass0::CalibrateAndStore(TH2F *histoVertexTimestamp, TH2F *histoDeltatTimestamp, const Char_t *dbString)
 526 {
 527   /*
 528    * calibrate and store
 529    */
 530
 531   /* check pointers */
 532   if (!histoVertexTimestamp || !histoDeltatTimestamp)
 533     return kFALSE;
 534
 535   /*** GET RUN-NUMBER AND START-TIMESTAMP ***/
 536
 537   TString str;
 538   TObjArray *strarr = NULL;
 539   TObjString *ostr = NULL;
 540   str = histoVertexTimestamp->GetTitle();
 541   strarr = str.Tokenize(",");
 542   if (!strarr) {
 543     AliError("problems whith tokenize histogram title");
 544     fStatus = kDataError;
 545     return kFALSE;
 546   }
 547
 548   /* run number */
 549   ostr = (TObjString *)strarr->At(0);
 550   if (!ostr) {
 551     AliError("problems while getting run number from histogram title");
 552     fStatus = kDataError;
 553     return kFALSE;
 554   }
 555   str = ostr->GetString();
 556   if (!str.BeginsWith("run:")) {
 557     AliError("problems while getting run number from histogram title");
 558     fStatus = kDataError;
 559     return kFALSE;
 560   }
 561   str.Remove(0, 5);
 562   Int_t runNb = atoi(str.Data());
 563   if (runNb <= 0) {
 564     AliError(Form("bad run number: %d", runNb));
 565     fStatus = kDataError;
 566     return kFALSE;
 567   }
 568   AliInfo(Form("got run number: %d", runNb));
 569
 570   /* start timestamp */
 571   ostr = (TObjString *)strarr->At(1);
 572   if (!ostr) {
 573     AliError("problems while getting start timestamp from histogram title");
 574     fStatus = kDataError;
 575     return kFALSE;
 576   }
 577   str = ostr->GetString();
 578   str.Remove(0, 1); /* remove empty space at the beginning */
 579   if (!str.BeginsWith("startTimestamp:")) {
 580     AliError("problems while getting start timestamp from histogram title");
 581     fStatus = kDataError;
 582     return kFALSE;
 583   }
 584   str.Remove(0, 16);
 585   UInt_t startTimestamp = atoi(str.Data());
 586   if (startTimestamp <= 0) {
 587     AliError(Form("bad start timestamp: %d", startTimestamp));
 588     fStatus = kDataError;
 589     return kFALSE;
 590   }
 591   TTimeStamp ts = startTimestamp;
 592   AliInfo(Form("got start timestamp: %d (%s)", startTimestamp, ts.AsString()));
 593
 594   /* BPTX */
 595   ostr = (TObjString *)strarr->At(2);
 596   if (!ostr) {
 597     AliError("problems while getting BPTX from histogram title");
 598     fStatus = kDataError;
 599     return kFALSE;
 600   }
 601   str = ostr->GetString();
 602   str.Remove(0, 1); /* remove empty space at the beginning */
 603   if (!str.BeginsWith("BPTX:")) {
 604     AliError("problems while getting BPTX from histogram title");
 605     fStatus = kDataError;
 606     return kFALSE;
 607   }
 608   str.Remove(0, 6);
 609   Bool_t useBPTX = atoi(str.Data());
 610   AliInfo(Form("got BPTX: %d", useBPTX));
 611
 612   /*** CALIBRATION STAGE ***/
 613
 614   /* get fit function */
 615   TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
 616
 617   /* projection-x */
 618   TH1D *histoVertexTimestamppx = histoVertexTimestamp->ProjectionX("histoVertexTimestamppx");
 619   TH1D *histoDeltatTimestamppx = histoDeltatTimestamp->ProjectionX("histoDeltatTimestamppx");
 620
 621   /* check statistics */
 622   if (histoVertexTimestamppx->Integral() < fgMinVertexIntegral ||
 623       histoDeltatTimestamppx->Integral() < fgMinDeltatIntegral) {
 624     fStatus = kLowStatistics;
 625     return kFALSE;
 626   }
 627
 628   /* define mix and max time bin */
 629   Int_t minBin = histoVertexTimestamppx->FindFirstBinAbove(0);
 630   Int_t maxBin = histoVertexTimestamppx->FindLastBinAbove(0);
 631   Float_t minTime = histoVertexTimestamppx->GetBinLowEdge(minBin);
 632   Float_t maxTime = histoVertexTimestamppx->GetBinLowEdge(maxBin + 1);
 633   AliInfo(Form("min/max time defined: %d < t < %d s [%d, %d]", (Int_t)minTime, (Int_t)maxTime, minBin, maxBin));
 634
 635   /* loop over time bins */
 636   Int_t nPoints = 0;
 637   Float_t time[fgkMaxNumberOfPoints], timeerr[fgkMaxNumberOfPoints];
 638   Float_t vertexMean[fgkMaxNumberOfPoints], vertexMeanerr[fgkMaxNumberOfPoints];
 639   Float_t vertexSigma[fgkMaxNumberOfPoints], vertexSigmaerr[fgkMaxNumberOfPoints];
 640   Float_t timeZeroMean[fgkMaxNumberOfPoints], timeZeroMeanerr[fgkMaxNumberOfPoints];
 641   Float_t timeZeroSigma[fgkMaxNumberOfPoints], timeZeroSigmaerr[fgkMaxNumberOfPoints];
 642   Float_t averageTimeZero = 0.;
 643   Float_t averageTimeSigma = 0.;
 644   Float_t averageVertexSpread = 0.;
 645   for (Int_t ibin = minBin; ibin <= maxBin; ibin++) {
 646
 647     /* define time window */
 648     Int_t startBin = ibin;
 649     Int_t endBin = ibin;
 650     while(histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegralSample ||
 651           histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegralSample) {
 652       if (endBin < maxBin) endBin++;
 653       else if (startBin > minBin) startBin--;
 654       else break;
 655     }
 656     if (histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegral ||
 657         histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegral) continue;
 658     Float_t startTime = histoVertexTimestamppx->GetBinLowEdge(startBin);
 659     Float_t endTime = histoVertexTimestamppx->GetBinLowEdge(endBin + 1);
 660     Float_t vertexIntegral = histoVertexTimestamppx->Integral(startBin, endBin);
 661     Float_t deltatIntegral = histoDeltatTimestamppx->Integral(startBin, endBin);
 662     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));
 663
 664     /* projection-y */
 665     TH1D *histoVertexTimestamppy = histoVertexTimestamp->ProjectionY("histoVertexTimestamppy", startBin, endBin);
 666     TH1D *histoDeltatTimestamppy = histoDeltatTimestamp->ProjectionY("histoDeltatTimestamppy", startBin, endBin);
 667
 668     /* average time */
 669     histoVertexTimestamppx->GetXaxis()->SetRange(startBin, endBin);
 670     time[nPoints] = histoVertexTimestamppx->GetMean();
 671     timeerr[nPoints] = histoVertexTimestamppx->GetMeanError();
 672
 673     /* fit vertex */
 674     if (FitPeak(fitFunc, histoVertexTimestamppy, 10., 3., 3.) != 0) {
 675       AliError("troubles fitting vertex, skip");
 676       delete histoVertexTimestamppy;
 677       delete histoDeltatTimestamppy;
 678       continue;
 679     }
 680     vertexMean[nPoints] = fitFunc->GetParameter(1);
 681     vertexMeanerr[nPoints] = fitFunc->GetParError(1);
 682     vertexSigma[nPoints] = fitFunc->GetParameter(2);
 683     vertexSigmaerr[nPoints] = fitFunc->GetParError(2);
 684     averageVertexSpread += fitFunc->GetParameter(2);
 685
 686     /* fit time-zero */
 687     if (FitPeak(fitFunc, histoDeltatTimestamppy, 500., 2., 1.) != 0) {
 688       AliError("troubles fitting time-zero TRACKS, skip");
 689       delete histoVertexTimestamppy;
 690       delete histoDeltatTimestamppy;
 691       continue;
 692     }
 693     timeZeroMean[nPoints] = fitFunc->GetParameter(1);
 694     timeZeroMeanerr[nPoints] = fitFunc->GetParError(1);
 695     timeZeroSigma[nPoints] = fitFunc->GetParameter(2);
 696     timeZeroSigmaerr[nPoints] = fitFunc->GetParError(2);
 697     averageTimeZero += fitFunc->GetParameter(1);
 698     averageTimeSigma += fitFunc->GetParameter(2);
 699
 700     /* delete projection-y */
 701     delete histoVertexTimestamppy;
 702     delete histoDeltatTimestamppy;
 703
 704     /* increment n points */
 705     nPoints++;
 706
 707     /* set current bin */
 708     ibin = endBin;
 709
 710   } /* end of loop over time bins */
 711
 712   /* delete projection-x */
 713   delete histoVertexTimestamppx;
 714   delete histoDeltatTimestamppx;
 715
 716   /* check points */
 717   if (nPoints <= 0) {
 718     AliError("no measurement available, quit");
 719     fStatus = kNoMeasurement;
 720     return kFALSE;
 721   }
 722   AliInfo(Form("average time-zero  = %f", averageTimeZero / nPoints));
 723   AliInfo(Form("average time-sigma = %f", averageTimeSigma / nPoints));
 724   AliInfo(Form("average v-spread   = %f", averageVertexSpread / nPoints));
 725
 726   /*** CREATE RUN PARAMS OBJECT ***/
 727
 728 #if 0
 729   /* get start time from GRP */
 730   TGrid::Connect("alien://", 0, 0, "t");
 731   AliCDBManager *cdb = AliCDBManager::Instance();
 732   cdb->SetDefaultStorage("raw://");
 733   cdb->SetRun(runNb);
 734   AliGRPManager grp;
 735   if (!grp.ReadGRPEntry()) {
 736     AliError("error while reading GRP entry");
 737     return kFALSE;
 738   }
 739   UInt_t startTimestamp = grp.GetGRPData()->GetTimeStart();
 740   TTimeStamp ts;
 741   ts = startTimestamp;
 742   AliInfo(Form("got start time from GRP: %s", ts.AsString()));
 743 #endif
 744
 745   /* create arrays */
 746   UInt_t timestamp[fgkMaxNumberOfPoints];
 747   Float_t t0[fgkMaxNumberOfPoints];
 748   Float_t tofReso[fgkMaxNumberOfPoints];
 749   Float_t t0Spread[fgkMaxNumberOfPoints];
 750   for (Int_t ipoint = 0; ipoint < nPoints; ipoint++) {
 751     timestamp[ipoint] = (UInt_t)time[ipoint] + startTimestamp;
 752     t0[ipoint] = timeZeroMean[ipoint];
 753     tofReso[ipoint] = timeZeroSigma[ipoint];
 754     t0Spread[ipoint] = vertexSigma[ipoint] / 2.99792457999999984e-02;
 755   }
 756   UInt_t run[1] = {runNb};
 757   UInt_t runFirstPoint[1] = {0};
 758   UInt_t runLastPoint[1] = {nPoints - 1};
 759
 760   /* create run params object */
 761   AliTOFRunParams obj(nPoints, 1);
 762   AliInfo(Form("create run params object for run %d with %d points", runNb, nPoints));
 763   obj.SetTimestamp(timestamp);
 764   obj.SetT0(t0);
 765   obj.SetTOFResolution(tofReso);
 766   obj.SetT0Spread(t0Spread);
 767   obj.SetRunNb(run);
 768   obj.SetRunFirstPoint(runFirstPoint);
 769   obj.SetRunLastPoint(runLastPoint);
 770   obj.SetUseLHCClockPhase(useBPTX);
 771
 772   /*** CREATE OCDB ENTRY ***/
 773
 774   if (!dbString) {
 775     AliError("cannot store object because of NULL string");
 776     fStatus = kStoreError;
 777     return kFALSE;
 778   }
 779
 780   /* install run params object in OCDB */
 781   AliCDBManager *cdb = AliCDBManager::Instance();
 782   AliCDBStorage *sto = cdb->GetStorage(dbString);
 783   if (!sto) {
 784     AliError(Form("cannot get storage %s", dbString));
 785     fStatus = kStoreError;
 786     return kFALSE;
 787   }
 788   AliCDBId id("TOF/Calib/RunParams", runNb, runNb);
 789   AliCDBMetaData md;
 790   md.SetResponsible("Roberto Preghenella");
 791   md.SetComment("offline TOF run parameters");
 792   md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
 793   md.SetBeamPeriod(0);
 794   if (!sto->Put(&obj, id, &md)) {
 795     fStatus = kStoreError;
 796     AliError(Form("error while putting object in storage %s", dbString));
 797     return kFALSE;
 798   }
 799
 800   /* success */
 801   return kTRUE;
 802 }
 803
 804 //_____________________________________________________________
 805
 806 Int_t
 807 AliTOFAnalysisTaskCalibPass0::FitPeak(TF1 *fitFunc, TH1D *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
 808 {
 809   /*
 810    * fit peak
 811    */
 812
 813   Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
 814   Double_t fitMin = fitCent - nSigmaMin * startSigma;
 815   Double_t fitMax = fitCent + nSigmaMax * startSigma;
 816   if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 817   if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 818   fitFunc->SetParameter(1, fitCent);
 819   fitFunc->SetParameter(2, startSigma);
 820   Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
 821   if (fitres != 0) return fitres;
 822   /* refit with better range */
 823   for (Int_t i = 0; i < 3; i++) {
 824     fitCent = fitFunc->GetParameter(1);
 825     fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
 826     fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
 827     if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 828     if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 829     fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
 830     if (fitres != 0) return fitres;
 831   }
 832   return fitres;
 833 }