ProcessOutput update: deal with both TFileMerged merged data (histos inside TList...
[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 Int_t AliTOFAnalysisTaskCalibPass0::fgkMaxNumberOfPoints = 10000; // max number of points
55 const Double_t AliTOFAnalysisTaskCalibPass0::fgkMinVertexIntegral = 1000.;
56 const Double_t AliTOFAnalysisTaskCalibPass0::fgkMinDeltatIntegral = 20000.;
57
58 //_______________________________________________________
59   
60 AliTOFAnalysisTaskCalibPass0::AliTOFAnalysisTaskCalibPass0() :
61   AliAnalysisTaskSE("TOFCalib-Pass0"),
62   fInitFlag(kFALSE),
63   fEventSelectionFlag(kFALSE),
64   fVertexSelectionFlag(kFALSE),
65   fVertexCut(10.),
66   fRunNumber(0),
67   fESDEvent(NULL),
68   fEventCuts(new AliPhysicsSelection()),
69   fTrackCuts(new AliESDtrackCuts()),
70   fStartTime(0),
71   fEndTime(0),
72   fEventTime(0),
73   fElapsedTime(0),
74   fkVertex(NULL),
75   fGRPManager(new AliGRPManager()),
76   fkGRPObject(NULL),
77   fTOFcalib(new AliTOFcalib()),
78   fHistoList(new TList()),
79   fHistoVertexTimestamp(NULL),
80   fHistoDeltatTimestamp(NULL),
81   fHistoDeltazEta(NULL),
82   fHistoDeltazCosTheta(NULL),
83   fHistoAcceptedTracksEtaPt(NULL),
84   fHistoMatchedTracksEtaPt(NULL)
85 {
86   /* 
87    * default constructor 
88    */
89
90   /* define output */
91   fHistoList->SetOwner(kTRUE);
92   DefineOutput(1, TList::Class());
93 }
94
95 //_______________________________________________________
96
97 AliTOFAnalysisTaskCalibPass0::~AliTOFAnalysisTaskCalibPass0()
98 {
99   /*
100    * default destructor
101    */
102
103   if (fHistoList) delete fHistoList;
104 }
105
106 //_______________________________________________________
107
108 Bool_t
109 AliTOFAnalysisTaskCalibPass0::InitRun()
110 {
111   /*
112    * init run
113    */
114
115   /* get ESD event */
116   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
117   if (!fESDEvent) {
118     AliError("cannot get ESD event");
119     return kFALSE;
120   }
121   /* get run number */
122   Int_t runNb = fESDEvent->GetRunNumber();
123   /* check run already initialized */
124   if (fInitFlag && fRunNumber == runNb) return kTRUE;
125   if (fInitFlag && fRunNumber != runNb) {
126     AliFatal(Form("run number has changed within same job and this is not allowed: %d -> %d", fRunNumber, runNb));
127     return kTRUE;
128   }
129   fInitFlag = kFALSE;
130   /* init cdb if not yet done */
131   AliCDBManager *cdb = AliCDBManager::Instance();
132   if (!cdb->IsDefaultStorageSet())
133     cdb->SetDefaultStorage("raw://");
134   if (cdb->GetRun() < 0)
135     cdb->SetRun(runNb);
136     /* init TOF calib */
137   if (!fTOFcalib->Init()) {
138     AliError("cannot init TOF calib");
139     return kFALSE;
140   }
141   /* get GRP data */
142   if (!fGRPManager->ReadGRPEntry()) {
143     AliError("error while reading \"GRPEntry\" from OCDB");
144     return kFALSE;
145   }
146   fkGRPObject = fGRPManager->GetGRPData();
147   if (!fkGRPObject) {
148     AliError("cannot get \"GRPData\" from GRP manager");
149     return kFALSE;
150   }
151   fStartTime = fkGRPObject->GetTimeStart();
152   fEndTime = fkGRPObject->GetTimeEnd();
153   AliInfo(Form("got \"GRPData\": startTime=%d, endTime=%d", fStartTime, fEndTime));
154   
155   AliInfo(Form("initialized for run %d", runNb));
156   fInitFlag = kTRUE;
157   fRunNumber = runNb;
158   
159   /* set histo title with run-number and start-time */
160   fHistoVertexTimestamp->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
161   fHistoDeltatTimestamp->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
162   fHistoDeltazEta->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
163   fHistoDeltazCosTheta->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
164   fHistoAcceptedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
165   fHistoMatchedTracksEtaPt->SetTitle(Form("run: %d, startTimestamp: %u", fRunNumber, fStartTime));
166   
167   return kTRUE;
168 }
169
170 //_______________________________________________________
171
172 Bool_t
173 AliTOFAnalysisTaskCalibPass0::InitEvent()
174 {
175   /*
176    * init event
177    */
178
179   /* get ESD event */
180   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
181   if (!fESDEvent) return kFALSE;
182   /* set event time and elapsed time */
183   fEventTime = fESDEvent->GetTimeStamp();
184   fElapsedTime = fESDEvent->GetTimeStamp() - fStartTime;
185   /* event selection */
186   if (fEventSelectionFlag && !fEventCuts->IsCollisionCandidate(fESDEvent)) return kFALSE;
187   /* vertex selection */
188   fkVertex = fESDEvent->GetPrimaryVertexTracks();
189   if (fVertexSelectionFlag && fkVertex->GetNContributors() < 1) {
190     fkVertex = fESDEvent->GetPrimaryVertexSPD();
191     if (fkVertex->GetNContributors() < 1) return kFALSE;
192   }
193   if (fVertexSelectionFlag && TMath::Abs(fkVertex->GetZ()) > fVertexCut) return kFALSE;
194   /* calibrate ESD if requested */
195   fTOFcalib->CalibrateESD(fESDEvent);
196
197   return kTRUE;
198 }
199
200 //_______________________________________________________
201
202 Bool_t
203 AliTOFAnalysisTaskCalibPass0::HasTOFMeasurement(const AliESDtrack *track) const
204 {
205   /*
206    * has TOF measurement
207    */
208
209   /* check TOF status flags */
210   if (!(track->GetStatus() & AliESDtrack::kTOFout) ||
211       !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE;
212   /* check integrated length */
213   if (track->GetIntegratedLength() < 350.) return kFALSE;
214
215   /* TOF measurement ok */
216   return kTRUE;
217 }
218
219 //_______________________________________________________
220
221 void
222 AliTOFAnalysisTaskCalibPass0::UserCreateOutputObjects()
223 {
224   /*
225    * user create output objects
226    */
227
228   /* time binning */
229   Int_t timeBins = 288;
230   Float_t timeMin = 0.;
231   Float_t timeMax = 24. * 3600.;
232   /* vertex binning */
233   Int_t vertexBins = 200;
234   Float_t vertexMin = -50.;
235   Float_t vertexMax = 50.;
236   /* deltat binning */
237   Int_t deltatBins = 2000;
238   Float_t deltatMin = -24400.;
239   Float_t deltatMax = 24400.;
240   /* deltaz binning */
241   Int_t deltazBins = 200;
242   Float_t deltazMin = -10.;
243   Float_t deltazMax = 10.;
244   /* eta binning */
245   Int_t etaBins = 200;
246   Float_t etaMin = -1.;
247   Float_t etaMax = 1.;
248   /* p binning */
249   Int_t pBins = 500;
250   Float_t pMin = 0.;
251   Float_t pMax = 5.;
252   
253   fHistoVertexTimestamp = new TH2F("hHistoVertexTimestamp", "Vertex position;elapsed time (s);z (cm);", timeBins, timeMin, timeMax, vertexBins, vertexMin, vertexMax);
254   fHistoList->Add(fHistoVertexTimestamp);
255
256   fHistoDeltatTimestamp = new TH2F("hHistoDeltatTimestamp", "Global time shift (T0-fill);elapsed time (s);t - t_{exp}^{(#pi)} (ps);", timeBins, timeMin, timeMax, deltatBins, deltatMin, deltatMax);
257   fHistoList->Add(fHistoDeltatTimestamp);
258
259   fHistoDeltazEta = new TH2F("hHistoDeltazEta", "Matching residuals (longitudinal);#eta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
260   fHistoList->Add(fHistoDeltazEta);
261
262   fHistoDeltazCosTheta = new TH2F("hHistoDeltazCosTheta", "Matching residuals (longitudinal);cos #theta;#Deltaz (cm);", etaBins, etaMin, etaMax, deltazBins, deltazMin, deltazMax);
263   fHistoList->Add(fHistoDeltazCosTheta);
264
265   fHistoAcceptedTracksEtaPt = new TH2F("hHistoAcceptedTracksEtaPt", "Accepted tracks;#eta;#p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
266   fHistoList->Add(fHistoAcceptedTracksEtaPt);
267
268   fHistoMatchedTracksEtaPt = new TH2F("hHistoMatchedTracksEtaPt", "Matched tracks;#eta;p_{T} (GeV/c);", etaBins, etaMin, etaMax, pBins, pMin, pMax);
269   fHistoList->Add(fHistoMatchedTracksEtaPt);
270
271   /* post data */
272   PostData(1, fHistoList);
273 }
274
275 //_______________________________________________________
276
277 void
278 AliTOFAnalysisTaskCalibPass0::UserExec(Option_t *)
279 {
280   /*
281    * user exec
282    */
283
284   /* init run */
285   if (!InitRun()) return;
286   /* init event */
287   if (!InitEvent()) return;
288
289   /*** ACCEPTED EVENT ***/
290   
291   /* fill vertex histo */
292   fHistoVertexTimestamp->Fill(fElapsedTime, fkVertex->GetZ());
293
294   /* loop over ESD tracks */
295   Int_t nTracks = fESDEvent->GetNumberOfTracks();
296   AliESDtrack *track;
297   Double_t eta, costheta, pt, time, timei[AliPID::kSPECIES], deltat, deltaz;
298   for (Int_t itrk = 0; itrk < nTracks; itrk++) {
299     /* get track */
300     track = fESDEvent->GetTrack(itrk);
301     if (!track) continue;
302     /* check accept track */
303     if (!fTrackCuts->AcceptTrack(track)) continue;
304     /* get track info */
305     eta = track->Eta();
306     costheta = TMath::Cos(track->Theta());
307     pt = track->Pt();
308     /* fill accepted tracks histo */
309     fHistoAcceptedTracksEtaPt->Fill(eta, pt);
310     /* check TOF measurement */
311     if (!HasTOFMeasurement(track)) continue;
312
313     /*** ACCEPTED TRACK WITH TOF MEASUREMENT ***/
314
315     /* fill matched tracks histo */
316     fHistoMatchedTracksEtaPt->Fill(eta, pt);
317     /* get TOF info */
318     time = track->GetTOFsignal();
319     track->GetIntegratedTimes(timei);
320     deltat = time - timei[AliPID::kPion];
321     deltaz = track->GetTOFsignalDz();
322     
323     /* fill histos */
324     fHistoDeltatTimestamp->Fill(fElapsedTime, deltat);
325     fHistoDeltazEta->Fill(eta, deltaz);
326     fHistoDeltazCosTheta->Fill(costheta, deltaz);
327     
328   } /* end of loop over ESD tracks */
329
330   /* post data */
331   PostData(1, fHistoList);
332   
333 }
334
335 //_______________________________________________________
336
337 Bool_t
338 AliTOFAnalysisTaskCalibPass0::ProcessOutput(const Char_t *filename, const Char_t *dbString)
339 {
340   /*
341    * process output
342    */
343
344   /* open file */
345   TFile *file = TFile::Open(filename);
346   if (!file || !file->IsOpen()) {
347     AliError(Form("cannot open output file %s", filename));
348     return kFALSE;
349   }
350   /* get histograms */
351   TList *list = (TList *)file->Get("Histos");
352   TH2F *histoVertexTimestamp = NULL;
353   TH2F *histoDeltatTimestamp = NULL;
354   TH2F *histoDeltazEta = NULL;
355   TH2F *histoDeltazCosTheta = NULL;
356   TH2F *histoAcceptedTracksEtaPt = NULL;
357   TH2F *histoMatchedTracksEtaPt = NULL;
358   if (list) {
359     AliInfo(Form("getting histograms from \"Histos\" list from file %s", filename));
360     histoVertexTimestamp = (TH2F *)list->FindObject("hHistoVertexTimestamp");
361     histoDeltatTimestamp = (TH2F *)list->FindObject("hHistoDeltatTimestamp");
362     histoDeltazEta = (TH2F *)list->FindObject("hHistoDeltazEta");
363     histoDeltazCosTheta = (TH2F *)list->FindObject("hHistoDeltazCosTheta");
364     histoAcceptedTracksEtaPt = (TH2F *)list->FindObject("hHistoAcceptedTracksEtaPt");
365     histoMatchedTracksEtaPt = (TH2F *)list->FindObject("hHistoMatchedTracksEtaPt");
366   }
367   else {
368     AliInfo(Form("getting histograms directly from file %s", filename));
369     histoVertexTimestamp = (TH2F *)file->Get("hHistoVertexTimestamp");
370     histoDeltatTimestamp = (TH2F *)file->Get("hHistoDeltatTimestamp");
371     histoDeltazEta = (TH2F *)file->Get("hHistoDeltazEta");
372     histoDeltazCosTheta = (TH2F *)file->Get("hHistoDeltazCosTheta");
373     histoAcceptedTracksEtaPt = (TH2F *)file->Get("hHistoAcceptedTracksEtaPt");
374     histoMatchedTracksEtaPt = (TH2F *)file->Get("hHistoMatchedTracksEtaPt");
375   }
376   /* check histos */ 
377   if (!histoVertexTimestamp) {
378     AliError(Form("cannot get \"hHistoVertexTimestamp\" object from file %s", filename));
379     return kFALSE;
380   }
381   if (!histoDeltatTimestamp) {
382     AliError(Form("cannot get \"hHistoDeltatTimestamp\" object from file %s", filename));
383     return kFALSE;
384   }
385   if (!histoDeltazEta) {
386     AliError(Form("cannot get \"hHistoDeltazEta\" object from file %s", filename));
387     return kFALSE;
388   }
389   if (!histoDeltazCosTheta) {
390     AliError(Form("cannot get \"hHistoDeltazCosTheta\" object from file %s", filename));
391     return kFALSE;
392   }
393   if (!histoAcceptedTracksEtaPt) {
394     AliError(Form("cannot get \"hHistoAccptedTracksEtaPt\" object from file %s", filename));
395     return kFALSE;
396   }
397   if (!histoMatchedTracksEtaPt) {
398     AliError(Form("cannot get \"hHistoMatchedTracksEtaPt\" object from file %s", filename));
399     return kFALSE;
400   }
401   /* check matching performance */
402   if (!CheckMatchingPerformance(histoDeltazEta, histoAcceptedTracksEtaPt, histoMatchedTracksEtaPt)) {
403     AliError("error while checking matching efficiency");
404     return kFALSE;
405   }
406   /* calibrate and store */
407   if (!CalibrateAndStore(histoVertexTimestamp, histoDeltatTimestamp, dbString)) {
408     AliError("error while calibrating and storing");
409     return kFALSE;
410   }
411
412   /* success */
413   return kTRUE;
414 }
415
416 //_______________________________________________________
417
418 Bool_t
419 AliTOFAnalysisTaskCalibPass0::CheckMatchingPerformance(const TH2F *histoDeltazCosTheta, const TH2F *histoAcceptedTracksEtaPt, const TH2F *histoMatchedTracksEtaPt) const
420 {
421   /*
422    * check matching performance
423    */
424
425   /* check pointers */
426   if (!histoDeltazCosTheta || !histoAcceptedTracksEtaPt || !histoMatchedTracksEtaPt)
427     return kFALSE;
428   /* dummy for the time being */
429   return kTRUE;
430 }
431
432 //_______________________________________________________
433
434 Bool_t
435 AliTOFAnalysisTaskCalibPass0::CalibrateAndStore(TH2F *histoVertexTimestamp, TH2F *histoDeltatTimestamp, const Char_t *dbString)
436 {
437   /*
438    * calibrate and store
439    */
440
441   /* check pointers */
442   if (!histoVertexTimestamp || !histoDeltatTimestamp)
443     return kFALSE;
444
445   /*** GET RUN-NUMBER AND START-TIMESTAMP ***/
446
447   TString str;
448   TObjArray *strarr = NULL;
449   TObjString *ostr = NULL;
450   str = histoVertexTimestamp->GetTitle();
451   strarr = str.Tokenize(",");
452   if (!strarr) {
453     AliError("problems whith tokenize histogram title");
454     return kFALSE;
455   }
456   
457   /* run number */
458   ostr = (TObjString *)strarr->At(0);
459   if (!ostr) {
460     AliError("problems while getting run number from histogram title");
461     return kFALSE;
462   }
463   str = ostr->GetString();
464   if (!str.BeginsWith("run:")) {
465     AliError("problems while getting run number from histogram title");
466     return kFALSE;
467   }
468   str.Remove(0, 5);
469   Int_t runNb = atoi(str.Data());
470   if (runNb <= 0) {
471     AliError(Form("bad run number: %d", runNb));
472     return kFALSE;
473   }
474   AliInfo(Form("got run number: %d", runNb));
475
476   /* start timestamp */
477   ostr = (TObjString *)strarr->At(1);
478   if (!ostr) {
479     AliError("problems while getting start timestamp from histogram title");
480     return kFALSE;
481   }
482   str = ostr->GetString();
483   str.Remove(0, 1); /* remove empty space at the beginning */
484   if (!str.BeginsWith("startTimestamp:")) {
485     AliError("problems while getting start timestamp from histogram title");
486     return kFALSE;
487   }
488   str.Remove(0, 16);
489   UInt_t startTimestamp = atoi(str.Data());
490   if (startTimestamp <= 0) {
491     AliError(Form("bad start timestamp: %d", startTimestamp));
492     return kFALSE;
493   }
494   TTimeStamp ts = startTimestamp;
495   AliInfo(Form("got start timestamp: %d (%s)", startTimestamp, ts.AsString()));
496
497   /*** CALIBRATION STAGE ***/
498
499   /* get fit function */
500   TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
501
502   /* projection-x */
503   TH1D *histoVertexTimestamppx = histoVertexTimestamp->ProjectionX("histoVertexTimestamppx");
504   TH1D *histoDeltatTimestamppx = histoDeltatTimestamp->ProjectionX("histoDeltatTimestamppx");
505
506   /* define mix and max time bin */
507   Int_t minBin = histoVertexTimestamppx->FindFirstBinAbove(0);
508   Int_t maxBin = histoVertexTimestamppx->FindLastBinAbove(0);
509   Float_t minTime = histoVertexTimestamppx->GetBinLowEdge(minBin);
510   Float_t maxTime = histoVertexTimestamppx->GetBinLowEdge(maxBin + 1);
511   AliInfo(Form("min/max time defined: %d < t < %d s [%d, %d]", (Int_t)minTime, (Int_t)maxTime, minBin, maxBin));
512
513   /* loop over time bins */
514   Int_t nPoints = 0;
515   Float_t time[fgkMaxNumberOfPoints], timeerr[fgkMaxNumberOfPoints];
516   Float_t vertexMean[fgkMaxNumberOfPoints], vertexMeanerr[fgkMaxNumberOfPoints];
517   Float_t vertexSigma[fgkMaxNumberOfPoints], vertexSigmaerr[fgkMaxNumberOfPoints];
518   Float_t timeZeroMean[fgkMaxNumberOfPoints], timeZeroMeanerr[fgkMaxNumberOfPoints];
519   Float_t timeZeroSigma[fgkMaxNumberOfPoints], timeZeroSigmaerr[fgkMaxNumberOfPoints];
520   Float_t averageTimeZero = 0.;
521   Float_t averageTimeSigma = 0.;
522   Float_t averageVertexSpread = 0.;
523   for (Int_t ibin = minBin; ibin <= maxBin; ibin++) {
524
525     /* define time window */
526     Int_t startBin = ibin;
527     Int_t endBin = ibin;
528     while(histoVertexTimestamppx->Integral(startBin, endBin) < fgkMinVertexIntegral ||
529           histoDeltatTimestamppx->Integral(startBin, endBin) < fgkMinDeltatIntegral) {
530       if (endBin < maxBin) endBin++;
531       else if (startBin > minBin) startBin--;
532       else break;
533     }
534     if (histoVertexTimestamppx->Integral(startBin, endBin) <= 0 ||
535         histoDeltatTimestamppx->Integral(startBin, endBin) <= 0) continue;
536     Float_t startTime = histoVertexTimestamppx->GetBinLowEdge(startBin);
537     Float_t endTime = histoVertexTimestamppx->GetBinLowEdge(endBin + 1);
538     Float_t vertexIntegral = histoVertexTimestamppx->Integral(startBin, endBin);
539     Float_t deltatIntegral = histoDeltatTimestamppx->Integral(startBin, endBin);
540     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));
541
542     /* projection-y */
543     TH1D *histoVertexTimestamppy = histoVertexTimestamp->ProjectionY("histoVertexTimestamppy", startBin, endBin);
544     TH1D *histoDeltatTimestamppy = histoDeltatTimestamp->ProjectionY("histoDeltatTimestamppy", startBin, endBin);
545
546     /* average time */
547     histoVertexTimestamppx->GetXaxis()->SetRange(startBin, endBin);
548     time[nPoints] = histoVertexTimestamppx->GetMean();
549     timeerr[nPoints] = histoVertexTimestamppx->GetMeanError();
550
551     /* fit vertex */
552     if (FitPeak(fitFunc, histoVertexTimestamppy, 10., 3., 3.) != 0) {
553       AliError("troubles fitting vertex, skip");
554       delete histoVertexTimestamppy;
555       delete histoDeltatTimestamppy;
556       continue;
557     }
558     vertexMean[nPoints] = fitFunc->GetParameter(1);
559     vertexMeanerr[nPoints] = fitFunc->GetParError(1);
560     vertexSigma[nPoints] = fitFunc->GetParameter(2);
561     vertexSigmaerr[nPoints] = fitFunc->GetParError(2);
562     averageVertexSpread += fitFunc->GetParameter(2);
563
564     /* fit time-zero */
565     if (FitPeak(fitFunc, histoDeltatTimestamppy, 500., 2., 1.) != 0) {
566       AliError("troubles fitting time-zero TRACKS, skip");
567       delete histoVertexTimestamppy;
568       delete histoDeltatTimestamppy;
569       continue;
570     }
571     timeZeroMean[nPoints] = fitFunc->GetParameter(1);
572     timeZeroMeanerr[nPoints] = fitFunc->GetParError(1);
573     timeZeroSigma[nPoints] = fitFunc->GetParameter(2);
574     timeZeroSigmaerr[nPoints] = fitFunc->GetParError(2);
575     averageTimeZero += fitFunc->GetParameter(1);
576     averageTimeSigma += fitFunc->GetParameter(2);
577
578     /* delete projection-y */
579     delete histoVertexTimestamppy;
580     delete histoDeltatTimestamppy;
581
582     /* increment n points */
583     nPoints++;
584
585     /* set current bin */
586     ibin = endBin;
587
588   } /* end of loop over time bins */
589
590   /* delete projection-x */
591   delete histoVertexTimestamppx;
592   delete histoDeltatTimestamppx;
593
594   /* check points */
595   if (nPoints <= 0) {
596     AliError("no measurement available, quit");
597     return kFALSE;
598   }
599   AliInfo(Form("average time-zero  = %f", averageTimeZero / nPoints));
600   AliInfo(Form("average time-sigma = %f", averageTimeSigma / nPoints));
601   AliInfo(Form("average v-spread   = %f", averageVertexSpread / nPoints));
602
603   /*** CREATE RUN PARAMS OBJECT ***/
604
605 #if 0
606   /* get start time from GRP */
607   TGrid::Connect("alien://", 0, 0, "t");
608   AliCDBManager *cdb = AliCDBManager::Instance();
609   cdb->SetDefaultStorage("raw://");
610   cdb->SetRun(runNb);
611   AliGRPManager grp;
612   if (!grp.ReadGRPEntry()) {
613     AliError("error while reading GRP entry");
614     return kFALSE;
615   }
616   UInt_t startTimestamp = grp.GetGRPData()->GetTimeStart();
617   TTimeStamp ts;
618   ts = startTimestamp;
619   AliInfo(Form("got start time from GRP: %s", ts.AsString()));
620 #endif
621     
622   /* create arrays */
623   UInt_t timestamp[fgkMaxNumberOfPoints];
624   Float_t t0[fgkMaxNumberOfPoints];
625   Float_t tofReso[fgkMaxNumberOfPoints];
626   Float_t t0Spread[fgkMaxNumberOfPoints];
627   for (Int_t ipoint = 0; ipoint < nPoints; ipoint++) {
628     timestamp[ipoint] = (UInt_t)time[ipoint] + startTimestamp;
629     t0[ipoint] = timeZeroMean[ipoint];
630     tofReso[ipoint] = timeZeroSigma[ipoint];
631     t0Spread[ipoint] = vertexSigma[ipoint] / 2.99792457999999984e-02;
632   }
633   UInt_t run[1] = {runNb};
634   UInt_t runFirstPoint[1] = {0};
635   UInt_t runLastPoint[1] = {nPoints - 1};
636   
637   /* create run params object */
638   AliTOFRunParams obj(nPoints, 1);
639   AliInfo(Form("create run params object for run %d with %d points", runNb, nPoints));
640   obj.SetTimestamp(timestamp);
641   obj.SetT0(t0);
642   obj.SetTOFResolution(tofReso);
643   obj.SetT0Spread(t0Spread);
644   obj.SetRunNb(run);
645   obj.SetRunFirstPoint(runFirstPoint);
646   obj.SetRunLastPoint(runLastPoint);
647   
648   /*** CREATE OCDB ENTRY ***/
649
650   if (!dbString) {
651     AliError("cannot store object because of NULL string");
652     return kFALSE;
653   }
654
655   /* install run params object in OCDB */
656   AliCDBManager *cdb = AliCDBManager::Instance();
657   AliCDBStorage *sto = cdb->GetStorage(dbString);
658   if (!sto) {
659     AliError(Form("cannot get storage %s", dbString));
660     return kFALSE;
661   }
662   AliCDBId id("TOF/Calib/RunParams", runNb, runNb);
663   AliCDBMetaData md;
664   md.SetResponsible("Roberto Preghenella");
665   md.SetComment("offline TOF run parameters");
666   md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
667   md.SetBeamPeriod(0);
668   if (!sto->Put(&obj, id, &md)) {
669     AliError(Form("error while putting object in storage %s", dbString));
670   }
671
672   /* success */
673   return kTRUE;
674 }
675
676 //_____________________________________________________________
677
678 Int_t
679 AliTOFAnalysisTaskCalibPass0::FitPeak(TF1 *fitFunc, TH1D *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
680 {
681   /*
682    * fit peak
683    */
684
685   Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
686   Double_t fitMin = fitCent - nSigmaMin * startSigma;
687   Double_t fitMax = fitCent + nSigmaMax * startSigma;
688   if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
689   if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
690   fitFunc->SetParameter(1, fitCent);
691   fitFunc->SetParameter(2, startSigma);
692   Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
693   if (fitres != 0) return fitres;
694   /* refit with better range */
695   for (Int_t i = 0; i < 3; i++) {
696     fitCent = fitFunc->GetParameter(1);
697     fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
698     fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
699     if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
700     if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
701     fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
702     if (fitres != 0) return fitres;
703   }
704   return fitres;
705 }