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   /* set histo title with run-number and start-time */
 173   fHistoVertexTimestamp->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 174   fHistoDeltatTimestamp->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 175   fHistoDeltazEta->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 176   fHistoDeltazCosTheta->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 177   fHistoAcceptedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 178   fHistoMatchedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
 179
 180   return kTRUE;
 181 }
 182
 183 //_______________________________________________________
 184
 185 Bool_t
 186 AliTOFAnalysisTaskCalibPass0::InitEvent()
 187 {
 188   /*
 189    * init event
 190    */
 191
 192   /* get ESD event */
 193   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 194   if (!fESDEvent) return kFALSE;
 195   /* set event time and elapsed time */
 196   fEventTime = fESDEvent->GetTimeStamp();
 197   fElapsedTime = fESDEvent->GetTimeStamp() - fStartTime;
 198   /* event selection */
 199   if (fEventSelectionFlag && !fEventCuts->IsCollisionCandidate(fESDEvent)) return kFALSE;
 200   /* vertex selection */
 201   fkVertex = fESDEvent->GetPrimaryVertexTracks();
 202   if (fVertexSelectionFlag && fkVertex->GetNContributors() < 1) {
 203     fkVertex = fESDEvent->GetPrimaryVertexSPD();
 204     if (fkVertex->GetNContributors() < 1) return kFALSE;
 205   }
 206   if (fVertexSelectionFlag && TMath::Abs(fkVertex->GetZ()) > fVertexCut) return kFALSE;
 207   /* calibrate ESD if requested */
 208   fTOFcalib->CalibrateESD(fESDEvent);
 209
 210   return kTRUE;
 211 }
 212
 213 //_______________________________________________________
 214
 215 Bool_t
 216 AliTOFAnalysisTaskCalibPass0::HasTOFMeasurement(const AliESDtrack *track) const
 217 {
 218   /*
 219    * has TOF measurement
 220    */
 221
 222   /* check TOF status flags */
 223   if (!(track->GetStatus() & AliESDtrack::kTOFout) ||
 224       !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE;
 225   /* check integrated length */
 226   if (track->GetIntegratedLength() < 350.) return kFALSE;
 227
 228   /* TOF measurement ok */
 229   return kTRUE;
 230 }
 231
 232 //_______________________________________________________
 233
 234 void
 235 AliTOFAnalysisTaskCalibPass0::UserCreateOutputObjects()
 236 {
 237   /*
 238    * user create output objects
 239    */
 240
 241   /* time binning */
 242   Int_t timeBins = 288;
 243   Float_t timeMin = 0.;
 244   Float_t timeMax = 24. * 3600.;
 245   /* vertex binning */
 246   Int_t vertexBins = 200;
 247   Float_t vertexMin = -50.;
 248   Float_t vertexMax = 50.;
 249   /* deltat binning */
 250   Int_t deltatBins = 2000;
 251   Float_t deltatMin = -24400.;
 252   Float_t deltatMax = 24400.;
 253   /* deltaz binning */
 254   Int_t deltazBins = 200;
 255   Float_t deltazMin = -10.;
 256   Float_t deltazMax = 10.;
 257   /* eta binning */
 258   Int_t etaBins = 200;
 259   Float_t etaMin = -1.;
 260   Float_t etaMax = 1.;
 261   /* p binning */
 262   Int_t pBins = 500;
 263   Float_t pMin = 0.;
 264   Float_t pMax = 5.;
 265
 266   fHistoVertexTimestamp = new TH2F("hHistoVertexTimestamp", "Vertex position;elapsed time (s);z (cm);", timeBins, timeMin, timeMax, vertexBins, vertexMin, vertexMax);
 267   fHistoList->Add(fHistoVertexTimestamp);
 268
 269   fHistoDeltatTimestamp = new TH2F("hHistoDeltatTimestamp", "Global time shift (T0-fill);elapsed time (s);t - t_{exp}^{(#pi)} (ps);", timeBins, timeMin, timeMax, deltatBins, deltatMin, deltatMax);
 270   fHistoList->Add(fHistoDeltatTimestamp);
 271
 272   fHistoDeltazEta = new TH2F("hHistoDeltazEta", "Matching residuals (longitudinal);#eta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
 273   fHistoList->Add(fHistoDeltazEta);
 274
 275   fHistoDeltazCosTheta = new TH2F("hHistoDeltazCosTheta", "Matching residuals (longitudinal);cos #theta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
 276   fHistoList->Add(fHistoDeltazCosTheta);
 277
 278   fHistoAcceptedTracksEtaPt = new TH2F("hHistoAcceptedTracksEtaPt", "Accepted tracks;#eta;#p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
 279   fHistoList->Add(fHistoAcceptedTracksEtaPt);
 280
 281   fHistoMatchedTracksEtaPt = new TH2F("hHistoMatchedTracksEtaPt", "Matched tracks;#eta;p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
 282   fHistoList->Add(fHistoMatchedTracksEtaPt);
 283
 284   /* post data */
 285   PostData(1, fHistoList);
 286 }
 287
 288 //_______________________________________________________
 289
 290 void
 291 AliTOFAnalysisTaskCalibPass0::UserExec(Option_t *)
 292 {
 293   /*
 294    * user exec
 295    */
 296
 297   /* init run */
 298   if (!InitRun()) return;
 299   /* init event */
 300   if (!InitEvent()) return;
 301
 302   /*** ACCEPTED EVENT ***/
 303
 304   /* fill vertex histo */
 305   fHistoVertexTimestamp->Fill(fElapsedTime, fkVertex->GetZ());
 306
 307   /* loop over ESD tracks */
 308   Int_t nTracks = fESDEvent->GetNumberOfTracks();
 309   AliESDtrack *track;
 310   Double_t eta, costheta, pt, time, timei[AliPID::kSPECIES], deltat, deltaz;
 311   for (Int_t itrk = 0; itrk < nTracks; itrk++) {
 312     /* get track */
 313     track = fESDEvent->GetTrack(itrk);
 314     if (!track) continue;
 315     /* check accept track */
 316     if (!fTrackCuts->AcceptTrack(track)) continue;
 317     /* get track info */
 318     eta = track->Eta();
 319     costheta = TMath::Cos(track->Theta());
 320     pt = track->Pt();
 321     /* fill accepted tracks histo */
 322     fHistoAcceptedTracksEtaPt->Fill(eta, pt);
 323     /* check TOF measurement */
 324     if (!HasTOFMeasurement(track)) continue;
 325
 326     /*** ACCEPTED TRACK WITH TOF MEASUREMENT ***/
 327
 328     /* fill matched tracks histo */
 329     fHistoMatchedTracksEtaPt->Fill(eta, pt);
 330     /* get TOF info */
 331     time = track->GetTOFsignal();
 332     track->GetIntegratedTimes(timei);
 333     deltat = time - timei[AliPID::kPion];
 334     deltaz = track->GetTOFsignalDz();
 335
 336     /* fill histos */
 337     fHistoDeltatTimestamp->Fill(fElapsedTime, deltat);
 338     fHistoDeltazEta->Fill(eta, deltaz);
 339     fHistoDeltazCosTheta->Fill(costheta, deltaz);
 340
 341   } /* end of loop over ESD tracks */
 342
 343   /* post data */
 344   PostData(1, fHistoList);
 345
 346 }
 347
 348 //_______________________________________________________
 349
 350 Int_t
 351 AliTOFAnalysisTaskCalibPass0::GetStatus()
 352 {
 353   /*
 354    * get status
 355    */
 356
 357   switch (fStatus) {
 358
 359     /* OK, return zero */
 360   case kOk:
 361     return 0;
 362     break;
 363
 364     /* non-fatal error, return negative status */
 365   case kNotActive:
 366   case kLowStatistics:
 367   case kNoMeasurement:
 368     return -fStatus;
 369     break;
 370
 371     /* fatal error, return positive status */
 372   case kInputError:
 373   case kDataError:
 374   case kStoreError:
 375     return fStatus;
 376     break;
 377
 378     /* anything else, return negative large number */
 379   default:
 380     return -999;
 381     break;
 382   }
 383
 384   /* should never arrive here, anyway return negative large number */
 385   return -999;
 386 }
 387
 388 //_______________________________________________________
 389
 390 Bool_t
 391 AliTOFAnalysisTaskCalibPass0::ProcessOutput(const Char_t *filename, const Char_t *dbString)
 392 {
 393   /*
 394    * process output
 395    */
 396
 397   Int_t ret = DoProcessOutput(filename, dbString);
 398   Int_t status = GetStatus();
 399   if (status == 0) {
 400     AliInfo(Form("TOF calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 401   }
 402   else if (status > 0) {
 403     AliInfo(Form("TOF calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
 404   }
 405   else if (status < 0) {
 406     AliInfo(Form("TOF calibration failed (expected): %s (status=%d)", fgkStatusCodeName[fStatus], status));
 407   }
 408
 409   return ret;
 410 }
 411
 412 //_______________________________________________________
 413
 414 Bool_t
 415 AliTOFAnalysisTaskCalibPass0::DoProcessOutput(const Char_t *filename, const Char_t *dbString)
 416 {
 417   /*
 418    * do process output
 419    */
 420
 421   /* reset status to OK */
 422   fStatus = kOk;
 423
 424   /* open file */
 425   TFile *file = TFile::Open(filename);
 426   if (!file || !file->IsOpen()) {
 427     AliError(Form("cannot open output file %s", filename));
 428     fStatus = kInputError;
 429     return kFALSE;
 430   }
 431   /* get histograms */
 432   TList *list = (TList *)file->Get("TOFHistos");
 433   TH2F *histoVertexTimestamp = NULL;
 434   TH2F *histoDeltatTimestamp = NULL;
 435   TH2F *histoDeltazEta = NULL;
 436   TH2F *histoDeltazCosTheta = NULL;
 437   TH2F *histoAcceptedTracksEtaPt = NULL;
 438   TH2F *histoMatchedTracksEtaPt = NULL;
 439   if (list) {
 440     AliInfo(Form("getting histograms from \"Histos\" list from file %s", filename));
 441     histoVertexTimestamp = (TH2F *)list->FindObject("hHistoVertexTimestamp");
 442     histoDeltatTimestamp = (TH2F *)list->FindObject("hHistoDeltatTimestamp");
 443     histoDeltazEta = (TH2F *)list->FindObject("hHistoDeltazEta");
 444     histoDeltazCosTheta = (TH2F *)list->FindObject("hHistoDeltazCosTheta");
 445     histoAcceptedTracksEtaPt = (TH2F *)list->FindObject("hHistoAcceptedTracksEtaPt");
 446     histoMatchedTracksEtaPt = (TH2F *)list->FindObject("hHistoMatchedTracksEtaPt");
 447   }
 448   else {
 449     AliInfo(Form("getting histograms directly from file %s", filename));
 450     histoVertexTimestamp = (TH2F *)file->Get("hHistoVertexTimestamp");
 451     histoDeltatTimestamp = (TH2F *)file->Get("hHistoDeltatTimestamp");
 452     histoDeltazEta = (TH2F *)file->Get("hHistoDeltazEta");
 453     histoDeltazCosTheta = (TH2F *)file->Get("hHistoDeltazCosTheta");
 454     histoAcceptedTracksEtaPt = (TH2F *)file->Get("hHistoAcceptedTracksEtaPt");
 455     histoMatchedTracksEtaPt = (TH2F *)file->Get("hHistoMatchedTracksEtaPt");
 456   }
 457   /* check histos */
 458   if (!histoVertexTimestamp) {
 459     AliError(Form("cannot get \"hHistoVertexTimestamp\" object from file %s", filename));
 460     fStatus = kInputError;
 461     return kFALSE;
 462   }
 463   if (!histoDeltatTimestamp) {
 464     AliError(Form("cannot get \"hHistoDeltatTimestamp\" object from file %s", filename));
 465     fStatus = kInputError;
 466     return kFALSE;
 467   }
 468   if (!histoDeltazEta) {
 469     AliError(Form("cannot get \"hHistoDeltazEta\" object from file %s", filename));
 470     fStatus = kInputError;
 471     return kFALSE;
 472   }
 473   if (!histoDeltazCosTheta) {
 474     AliError(Form("cannot get \"hHistoDeltazCosTheta\" object from file %s", filename));
 475     fStatus = kInputError;
 476     return kFALSE;
 477   }
 478   if (!histoAcceptedTracksEtaPt) {
 479     AliError(Form("cannot get \"hHistoAccptedTracksEtaPt\" object from file %s", filename));
 480     fStatus = kInputError;
 481     return kFALSE;
 482   }
 483   if (!histoMatchedTracksEtaPt) {
 484     AliError(Form("cannot get \"hHistoMatchedTracksEtaPt\" object from file %s", filename));
 485     fStatus = kInputError;
 486     return kFALSE;
 487   }
 488
 489   /* check matching performance */
 490   if (!CheckMatchingPerformance(histoDeltazEta, histoAcceptedTracksEtaPt, histoMatchedTracksEtaPt)) {
 491     AliError("error while checking matching efficiency");
 492     return kFALSE;
 493   }
 494   /* calibrate and store */
 495   if (!CalibrateAndStore(histoVertexTimestamp, histoDeltatTimestamp, dbString)) {
 496     AliError("error while calibrating and storing");
 497     return kFALSE;
 498   }
 499
 500   /* success */
 501   return kTRUE;
 502 }
 503
 504 //_______________________________________________________
 505
 506 Bool_t
 507 AliTOFAnalysisTaskCalibPass0::CheckMatchingPerformance(const TH2F *histoDeltazCosTheta, const TH2F *histoAcceptedTracksEtaPt, const TH2F *histoMatchedTracksEtaPt) const
 508 {
 509   /*
 510    * check matching performance
 511    */
 512
 513   /* check pointers */
 514   if (!histoDeltazCosTheta || !histoAcceptedTracksEtaPt || !histoMatchedTracksEtaPt)
 515     return kFALSE;
 516   /* dummy for the time being */
 517   return kTRUE;
 518 }
 519
 520 //_______________________________________________________
 521
 522 Bool_t
 523 AliTOFAnalysisTaskCalibPass0::CalibrateAndStore(TH2F *histoVertexTimestamp, TH2F *histoDeltatTimestamp, const Char_t *dbString)
 524 {
 525   /*
 526    * calibrate and store
 527    */
 528
 529   /* check pointers */
 530   if (!histoVertexTimestamp || !histoDeltatTimestamp)
 531     return kFALSE;
 532
 533   /*** GET RUN-NUMBER AND START-TIMESTAMP ***/
 534
 535   TString str;
 536   TObjArray *strarr = NULL;
 537   TObjString *ostr = NULL;
 538   str = histoVertexTimestamp->GetTitle();
 539   strarr = str.Tokenize(",");
 540   if (!strarr) {
 541     AliError("problems whith tokenize histogram title");
 542     fStatus = kDataError;
 543     return kFALSE;
 544   }
 545
 546   /* run number */
 547   ostr = (TObjString *)strarr->At(0);
 548   if (!ostr) {
 549     AliError("problems while getting run number from histogram title");
 550     fStatus = kDataError;
 551     return kFALSE;
 552   }
 553   str = ostr->GetString();
 554   if (!str.BeginsWith("run:")) {
 555     AliError("problems while getting run number from histogram title");
 556     fStatus = kDataError;
 557     return kFALSE;
 558   }
 559   str.Remove(0, 5);
 560   Int_t runNb = atoi(str.Data());
 561   if (runNb <= 0) {
 562     AliError(Form("bad run number: %d", runNb));
 563     fStatus = kDataError;
 564     return kFALSE;
 565   }
 566   AliInfo(Form("got run number: %d", runNb));
 567
 568   /* start timestamp */
 569   ostr = (TObjString *)strarr->At(1);
 570   if (!ostr) {
 571     AliError("problems while getting start timestamp from histogram title");
 572     fStatus = kDataError;
 573     return kFALSE;
 574   }
 575   str = ostr->GetString();
 576   str.Remove(0, 1); /* remove empty space at the beginning */
 577   if (!str.BeginsWith("startTimestamp:")) {
 578     AliError("problems while getting start timestamp from histogram title");
 579     fStatus = kDataError;
 580     return kFALSE;
 581   }
 582   str.Remove(0, 16);
 583   UInt_t startTimestamp = atoi(str.Data());
 584   if (startTimestamp <= 0) {
 585     AliError(Form("bad start timestamp: %d", startTimestamp));
 586     fStatus = kDataError;
 587     return kFALSE;
 588   }
 589   TTimeStamp ts = startTimestamp;
 590   AliInfo(Form("got start timestamp: %d (%s)", startTimestamp, ts.AsString()));
 591
 592   /*** CALIBRATION STAGE ***/
 593
 594   /* get fit function */
 595   TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
 596
 597   /* projection-x */
 598   TH1D *histoVertexTimestamppx = histoVertexTimestamp->ProjectionX("histoVertexTimestamppx");
 599   TH1D *histoDeltatTimestamppx = histoDeltatTimestamp->ProjectionX("histoDeltatTimestamppx");
 600
 601   /* check statistics */
 602   if (histoVertexTimestamppx->Integral() < fgMinVertexIntegral ||
 603       histoDeltatTimestamppx->Integral() < fgMinDeltatIntegral) {
 604     fStatus = kLowStatistics;
 605     return kFALSE;
 606   }
 607
 608   /* define mix and max time bin */
 609   Int_t minBin = histoVertexTimestamppx->FindFirstBinAbove(0);
 610   Int_t maxBin = histoVertexTimestamppx->FindLastBinAbove(0);
 611   Float_t minTime = histoVertexTimestamppx->GetBinLowEdge(minBin);
 612   Float_t maxTime = histoVertexTimestamppx->GetBinLowEdge(maxBin + 1);
 613   AliInfo(Form("min/max time defined: %d < t < %d s [%d, %d]", (Int_t)minTime, (Int_t)maxTime, minBin, maxBin));
 614
 615   /* loop over time bins */
 616   Int_t nPoints = 0;
 617   Float_t time[fgkMaxNumberOfPoints], timeerr[fgkMaxNumberOfPoints];
 618   Float_t vertexMean[fgkMaxNumberOfPoints], vertexMeanerr[fgkMaxNumberOfPoints];
 619   Float_t vertexSigma[fgkMaxNumberOfPoints], vertexSigmaerr[fgkMaxNumberOfPoints];
 620   Float_t timeZeroMean[fgkMaxNumberOfPoints], timeZeroMeanerr[fgkMaxNumberOfPoints];
 621   Float_t timeZeroSigma[fgkMaxNumberOfPoints], timeZeroSigmaerr[fgkMaxNumberOfPoints];
 622   Float_t averageTimeZero = 0.;
 623   Float_t averageTimeSigma = 0.;
 624   Float_t averageVertexSpread = 0.;
 625   for (Int_t ibin = minBin; ibin <= maxBin; ibin++) {
 626
 627     /* define time window */
 628     Int_t startBin = ibin;
 629     Int_t endBin = ibin;
 630     while(histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegralSample ||
 631           histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegralSample) {
 632       if (endBin < maxBin) endBin++;
 633       else if (startBin > minBin) startBin--;
 634       else break;
 635     }
 636     if (histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegral ||
 637         histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegral) continue;
 638     Float_t startTime = histoVertexTimestamppx->GetBinLowEdge(startBin);
 639     Float_t endTime = histoVertexTimestamppx->GetBinLowEdge(endBin + 1);
 640     Float_t vertexIntegral = histoVertexTimestamppx->Integral(startBin, endBin);
 641     Float_t deltatIntegral = histoDeltatTimestamppx->Integral(startBin, endBin);
 642     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));
 643
 644     /* projection-y */
 645     TH1D *histoVertexTimestamppy = histoVertexTimestamp->ProjectionY("histoVertexTimestamppy", startBin, endBin);
 646     TH1D *histoDeltatTimestamppy = histoDeltatTimestamp->ProjectionY("histoDeltatTimestamppy", startBin, endBin);
 647
 648     /* average time */
 649     histoVertexTimestamppx->GetXaxis()->SetRange(startBin, endBin);
 650     time[nPoints] = histoVertexTimestamppx->GetMean();
 651     timeerr[nPoints] = histoVertexTimestamppx->GetMeanError();
 652
 653     /* fit vertex */
 654     if (FitPeak(fitFunc, histoVertexTimestamppy, 10., 3., 3.) != 0) {
 655       AliError("troubles fitting vertex, skip");
 656       delete histoVertexTimestamppy;
 657       delete histoDeltatTimestamppy;
 658       continue;
 659     }
 660     vertexMean[nPoints] = fitFunc->GetParameter(1);
 661     vertexMeanerr[nPoints] = fitFunc->GetParError(1);
 662     vertexSigma[nPoints] = fitFunc->GetParameter(2);
 663     vertexSigmaerr[nPoints] = fitFunc->GetParError(2);
 664     averageVertexSpread += fitFunc->GetParameter(2);
 665
 666     /* fit time-zero */
 667     if (FitPeak(fitFunc, histoDeltatTimestamppy, 500., 2., 1.) != 0) {
 668       AliError("troubles fitting time-zero TRACKS, skip");
 669       delete histoVertexTimestamppy;
 670       delete histoDeltatTimestamppy;
 671       continue;
 672     }
 673     timeZeroMean[nPoints] = fitFunc->GetParameter(1);
 674     timeZeroMeanerr[nPoints] = fitFunc->GetParError(1);
 675     timeZeroSigma[nPoints] = fitFunc->GetParameter(2);
 676     timeZeroSigmaerr[nPoints] = fitFunc->GetParError(2);
 677     averageTimeZero += fitFunc->GetParameter(1);
 678     averageTimeSigma += fitFunc->GetParameter(2);
 679
 680     /* delete projection-y */
 681     delete histoVertexTimestamppy;
 682     delete histoDeltatTimestamppy;
 683
 684     /* increment n points */
 685     nPoints++;
 686
 687     /* set current bin */
 688     ibin = endBin;
 689
 690   } /* end of loop over time bins */
 691
 692   /* delete projection-x */
 693   delete histoVertexTimestamppx;
 694   delete histoDeltatTimestamppx;
 695
 696   /* check points */
 697   if (nPoints <= 0) {
 698     AliError("no measurement available, quit");
 699     fStatus = kNoMeasurement;
 700     return kFALSE;
 701   }
 702   AliInfo(Form("average time-zero  = %f", averageTimeZero / nPoints));
 703   AliInfo(Form("average time-sigma = %f", averageTimeSigma / nPoints));
 704   AliInfo(Form("average v-spread   = %f", averageVertexSpread / nPoints));
 705
 706   /*** CREATE RUN PARAMS OBJECT ***/
 707
 708 #if 0
 709   /* get start time from GRP */
 710   TGrid::Connect("alien://", 0, 0, "t");
 711   AliCDBManager *cdb = AliCDBManager::Instance();
 712   cdb->SetDefaultStorage("raw://");
 713   cdb->SetRun(runNb);
 714   AliGRPManager grp;
 715   if (!grp.ReadGRPEntry()) {
 716     AliError("error while reading GRP entry");
 717     return kFALSE;
 718   }
 719   UInt_t startTimestamp = grp.GetGRPData()->GetTimeStart();
 720   TTimeStamp ts;
 721   ts = startTimestamp;
 722   AliInfo(Form("got start time from GRP: %s", ts.AsString()));
 723 #endif
 724
 725   /* create arrays */
 726   UInt_t timestamp[fgkMaxNumberOfPoints];
 727   Float_t t0[fgkMaxNumberOfPoints];
 728   Float_t tofReso[fgkMaxNumberOfPoints];
 729   Float_t t0Spread[fgkMaxNumberOfPoints];
 730   for (Int_t ipoint = 0; ipoint < nPoints; ipoint++) {
 731     timestamp[ipoint] = (UInt_t)time[ipoint] + startTimestamp;
 732     t0[ipoint] = timeZeroMean[ipoint];
 733     tofReso[ipoint] = timeZeroSigma[ipoint];
 734     t0Spread[ipoint] = vertexSigma[ipoint] / 2.99792457999999984e-02;
 735   }
 736   UInt_t run[1] = {runNb};
 737   UInt_t runFirstPoint[1] = {0};
 738   UInt_t runLastPoint[1] = {nPoints - 1};
 739
 740   /* create run params object */
 741   AliTOFRunParams obj(nPoints, 1);
 742   AliInfo(Form("create run params object for run %d with %d points", runNb, nPoints));
 743   obj.SetTimestamp(timestamp);
 744   obj.SetT0(t0);
 745   obj.SetTOFResolution(tofReso);
 746   obj.SetT0Spread(t0Spread);
 747   obj.SetRunNb(run);
 748   obj.SetRunFirstPoint(runFirstPoint);
 749   obj.SetRunLastPoint(runLastPoint);
 750
 751   /*** CREATE OCDB ENTRY ***/
 752
 753   if (!dbString) {
 754     AliError("cannot store object because of NULL string");
 755     fStatus = kStoreError;
 756     return kFALSE;
 757   }
 758
 759   /* install run params object in OCDB */
 760   AliCDBManager *cdb = AliCDBManager::Instance();
 761   AliCDBStorage *sto = cdb->GetStorage(dbString);
 762   if (!sto) {
 763     AliError(Form("cannot get storage %s", dbString));
 764     fStatus = kStoreError;
 765     return kFALSE;
 766   }
 767   AliCDBId id("TOF/Calib/RunParams", runNb, runNb);
 768   AliCDBMetaData md;
 769   md.SetResponsible("Roberto Preghenella");
 770   md.SetComment("offline TOF run parameters");
 771   md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
 772   md.SetBeamPeriod(0);
 773   if (!sto->Put(&obj, id, &md)) {
 774     fStatus = kStoreError;
 775     AliError(Form("error while putting object in storage %s", dbString));
 776     return kFALSE;
 777   }
 778
 779   /* success */
 780   return kTRUE;
 781 }
 782
 783 //_____________________________________________________________
 784
 785 Int_t
 786 AliTOFAnalysisTaskCalibPass0::FitPeak(TF1 *fitFunc, TH1D *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
 787 {
 788   /*
 789    * fit peak
 790    */
 791
 792   Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
 793   Double_t fitMin = fitCent - nSigmaMin * startSigma;
 794   Double_t fitMax = fitCent + nSigmaMax * startSigma;
 795   if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 796   if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 797   fitFunc->SetParameter(1, fitCent);
 798   fitFunc->SetParameter(2, startSigma);
 799   Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
 800   if (fitres != 0) return fitres;
 801   /* refit with better range */
 802   for (Int_t i = 0; i < 3; i++) {
 803     fitCent = fitFunc->GetParameter(1);
 804     fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
 805     fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
 806     if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
 807     if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
 808     fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
 809     if (fitres != 0) return fitres;
 810   }
 811   return fitres;
 812 }