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