7f857aac2b31a5346a0d1c29c4adfa4f91a40a85
[u/mrichter/AliRoot.git] / 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 }