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