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