Optional scoring for background studies in HALL.
[u/mrichter/AliRoot.git] / TOF / AliTOFAnalysisTaskCalibPass0.cxx
CommitLineData
3c8efc07 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
50ClassImp(AliTOFAnalysisTaskCalibPass0)
51
52//_______________________________________________________
53
07901ecf 54const 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
3c8efc07 64const Int_t AliTOFAnalysisTaskCalibPass0::fgkMaxNumberOfPoints = 10000; // max number of points
d501ffec 65Double_t AliTOFAnalysisTaskCalibPass0::fgMinVertexIntegral = 100.;
66Double_t AliTOFAnalysisTaskCalibPass0::fgMinDeltatIntegral = 2000.;
67Double_t AliTOFAnalysisTaskCalibPass0::fgMinVertexIntegralSample = 1000.;
68Double_t AliTOFAnalysisTaskCalibPass0::fgMinDeltatIntegralSample = 20000.;
3c8efc07 69
70//_______________________________________________________
71
72AliTOFAnalysisTaskCalibPass0::AliTOFAnalysisTaskCalibPass0() :
73 AliAnalysisTaskSE("TOFCalib-Pass0"),
07901ecf 74 fStatus(kOk),
3c8efc07 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
110AliTOFAnalysisTaskCalibPass0::~AliTOFAnalysisTaskCalibPass0()
111{
112 /*
113 * default destructor
114 */
115
116 if (fHistoList) delete fHistoList;
117}
118
119//_______________________________________________________
120
121Bool_t
122AliTOFAnalysisTaskCalibPass0::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
3fd7a5cd 172 Int_t useBPTX = fTOFcalib->GetUseLHCClockPhase() ? 1 : 0;
173
3c8efc07 174 /* set histo title with run-number and start-time */
3fd7a5cd 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));
3c8efc07 181
182 return kTRUE;
183}
184
185//_______________________________________________________
186
187Bool_t
188AliTOFAnalysisTaskCalibPass0::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
217Bool_t
218AliTOFAnalysisTaskCalibPass0::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
236void
237AliTOFAnalysisTaskCalibPass0::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
292void
293AliTOFAnalysisTaskCalibPass0::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
07901ecf 352Int_t
353AliTOFAnalysisTaskCalibPass0::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
3c8efc07 392Bool_t
393AliTOFAnalysisTaskCalibPass0::ProcessOutput(const Char_t *filename, const Char_t *dbString)
394{
395 /*
396 * process output
397 */
398
07901ecf 399 Int_t ret = DoProcessOutput(filename, dbString);
400 Int_t status = GetStatus();
401 if (status == 0) {
402 AliInfo(Form("TOF calibration successful: %s (status=%d)", fgkStatusCodeName[fStatus], status));
403 }
404 else if (status > 0) {
405 AliInfo(Form("TOF calibration failed: %s (status=%d)", fgkStatusCodeName[fStatus], status));
406 }
407 else if (status < 0) {
408 AliInfo(Form("TOF calibration failed (expected): %s (status=%d)", fgkStatusCodeName[fStatus], status));
409 }
410
411 return ret;
412}
413
414//_______________________________________________________
415
416Bool_t
417AliTOFAnalysisTaskCalibPass0::DoProcessOutput(const Char_t *filename, const Char_t *dbString)
418{
419 /*
420 * do process output
421 */
422
423 /* reset status to OK */
424 fStatus = kOk;
425
3c8efc07 426 /* open file */
427 TFile *file = TFile::Open(filename);
428 if (!file || !file->IsOpen()) {
429 AliError(Form("cannot open output file %s", filename));
07901ecf 430 fStatus = kInputError;
3c8efc07 431 return kFALSE;
432 }
36c9ca5c 433 /* get histograms */
87ddb812 434 TList *list = (TList *)file->Get("TOFHistos");
36c9ca5c 435 TH2F *histoVertexTimestamp = NULL;
436 TH2F *histoDeltatTimestamp = NULL;
437 TH2F *histoDeltazEta = NULL;
438 TH2F *histoDeltazCosTheta = NULL;
439 TH2F *histoAcceptedTracksEtaPt = NULL;
440 TH2F *histoMatchedTracksEtaPt = NULL;
441 if (list) {
442 AliInfo(Form("getting histograms from \"Histos\" list from file %s", filename));
443 histoVertexTimestamp = (TH2F *)list->FindObject("hHistoVertexTimestamp");
444 histoDeltatTimestamp = (TH2F *)list->FindObject("hHistoDeltatTimestamp");
445 histoDeltazEta = (TH2F *)list->FindObject("hHistoDeltazEta");
446 histoDeltazCosTheta = (TH2F *)list->FindObject("hHistoDeltazCosTheta");
447 histoAcceptedTracksEtaPt = (TH2F *)list->FindObject("hHistoAcceptedTracksEtaPt");
448 histoMatchedTracksEtaPt = (TH2F *)list->FindObject("hHistoMatchedTracksEtaPt");
449 }
450 else {
451 AliInfo(Form("getting histograms directly from file %s", filename));
452 histoVertexTimestamp = (TH2F *)file->Get("hHistoVertexTimestamp");
453 histoDeltatTimestamp = (TH2F *)file->Get("hHistoDeltatTimestamp");
454 histoDeltazEta = (TH2F *)file->Get("hHistoDeltazEta");
455 histoDeltazCosTheta = (TH2F *)file->Get("hHistoDeltazCosTheta");
456 histoAcceptedTracksEtaPt = (TH2F *)file->Get("hHistoAcceptedTracksEtaPt");
457 histoMatchedTracksEtaPt = (TH2F *)file->Get("hHistoMatchedTracksEtaPt");
3c8efc07 458 }
36c9ca5c 459 /* check histos */
3c8efc07 460 if (!histoVertexTimestamp) {
461 AliError(Form("cannot get \"hHistoVertexTimestamp\" object from file %s", filename));
07901ecf 462 fStatus = kInputError;
3c8efc07 463 return kFALSE;
464 }
3c8efc07 465 if (!histoDeltatTimestamp) {
466 AliError(Form("cannot get \"hHistoDeltatTimestamp\" object from file %s", filename));
07901ecf 467 fStatus = kInputError;
3c8efc07 468 return kFALSE;
469 }
3c8efc07 470 if (!histoDeltazEta) {
471 AliError(Form("cannot get \"hHistoDeltazEta\" object from file %s", filename));
07901ecf 472 fStatus = kInputError;
3c8efc07 473 return kFALSE;
474 }
3c8efc07 475 if (!histoDeltazCosTheta) {
476 AliError(Form("cannot get \"hHistoDeltazCosTheta\" object from file %s", filename));
07901ecf 477 fStatus = kInputError;
3c8efc07 478 return kFALSE;
479 }
3c8efc07 480 if (!histoAcceptedTracksEtaPt) {
481 AliError(Form("cannot get \"hHistoAccptedTracksEtaPt\" object from file %s", filename));
07901ecf 482 fStatus = kInputError;
3c8efc07 483 return kFALSE;
484 }
3c8efc07 485 if (!histoMatchedTracksEtaPt) {
486 AliError(Form("cannot get \"hHistoMatchedTracksEtaPt\" object from file %s", filename));
07901ecf 487 fStatus = kInputError;
3c8efc07 488 return kFALSE;
489 }
07901ecf 490
3c8efc07 491 /* check matching performance */
492 if (!CheckMatchingPerformance(histoDeltazEta, histoAcceptedTracksEtaPt, histoMatchedTracksEtaPt)) {
493 AliError("error while checking matching efficiency");
494 return kFALSE;
495 }
496 /* calibrate and store */
497 if (!CalibrateAndStore(histoVertexTimestamp, histoDeltatTimestamp, dbString)) {
498 AliError("error while calibrating and storing");
499 return kFALSE;
500 }
501
502 /* success */
503 return kTRUE;
504}
505
506//_______________________________________________________
507
508Bool_t
509AliTOFAnalysisTaskCalibPass0::CheckMatchingPerformance(const TH2F *histoDeltazCosTheta, const TH2F *histoAcceptedTracksEtaPt, const TH2F *histoMatchedTracksEtaPt) const
510{
511 /*
512 * check matching performance
513 */
514
515 /* check pointers */
516 if (!histoDeltazCosTheta || !histoAcceptedTracksEtaPt || !histoMatchedTracksEtaPt)
517 return kFALSE;
518 /* dummy for the time being */
519 return kTRUE;
520}
521
522//_______________________________________________________
523
524Bool_t
525AliTOFAnalysisTaskCalibPass0::CalibrateAndStore(TH2F *histoVertexTimestamp, TH2F *histoDeltatTimestamp, const Char_t *dbString)
526{
527 /*
528 * calibrate and store
529 */
530
531 /* check pointers */
532 if (!histoVertexTimestamp || !histoDeltatTimestamp)
533 return kFALSE;
534
535 /*** GET RUN-NUMBER AND START-TIMESTAMP ***/
536
537 TString str;
538 TObjArray *strarr = NULL;
539 TObjString *ostr = NULL;
540 str = histoVertexTimestamp->GetTitle();
541 strarr = str.Tokenize(",");
542 if (!strarr) {
543 AliError("problems whith tokenize histogram title");
07901ecf 544 fStatus = kDataError;
3c8efc07 545 return kFALSE;
546 }
547
548 /* run number */
549 ostr = (TObjString *)strarr->At(0);
550 if (!ostr) {
551 AliError("problems while getting run number from histogram title");
07901ecf 552 fStatus = kDataError;
3c8efc07 553 return kFALSE;
554 }
555 str = ostr->GetString();
556 if (!str.BeginsWith("run:")) {
557 AliError("problems while getting run number from histogram title");
07901ecf 558 fStatus = kDataError;
3c8efc07 559 return kFALSE;
560 }
561 str.Remove(0, 5);
562 Int_t runNb = atoi(str.Data());
563 if (runNb <= 0) {
564 AliError(Form("bad run number: %d", runNb));
07901ecf 565 fStatus = kDataError;
3c8efc07 566 return kFALSE;
567 }
568 AliInfo(Form("got run number: %d", runNb));
569
570 /* start timestamp */
571 ostr = (TObjString *)strarr->At(1);
572 if (!ostr) {
573 AliError("problems while getting start timestamp from histogram title");
07901ecf 574 fStatus = kDataError;
3c8efc07 575 return kFALSE;
576 }
577 str = ostr->GetString();
578 str.Remove(0, 1); /* remove empty space at the beginning */
579 if (!str.BeginsWith("startTimestamp:")) {
580 AliError("problems while getting start timestamp from histogram title");
07901ecf 581 fStatus = kDataError;
3c8efc07 582 return kFALSE;
583 }
584 str.Remove(0, 16);
585 UInt_t startTimestamp = atoi(str.Data());
586 if (startTimestamp <= 0) {
587 AliError(Form("bad start timestamp: %d", startTimestamp));
07901ecf 588 fStatus = kDataError;
3c8efc07 589 return kFALSE;
590 }
591 TTimeStamp ts = startTimestamp;
592 AliInfo(Form("got start timestamp: %d (%s)", startTimestamp, ts.AsString()));
593
3fd7a5cd 594 /* BPTX */
595 ostr = (TObjString *)strarr->At(2);
596 if (!ostr) {
597 AliError("problems while getting BPTX from histogram title");
598 fStatus = kDataError;
599 return kFALSE;
600 }
601 str = ostr->GetString();
602 str.Remove(0, 1); /* remove empty space at the beginning */
603 if (!str.BeginsWith("BPTX:")) {
604 AliError("problems while getting BPTX from histogram title");
605 fStatus = kDataError;
606 return kFALSE;
607 }
608 str.Remove(0, 6);
609 Bool_t useBPTX = atoi(str.Data());
610 AliInfo(Form("got BPTX: %d", useBPTX));
611
3c8efc07 612 /*** CALIBRATION STAGE ***/
613
614 /* get fit function */
615 TF1 *fitFunc = (TF1 *)gROOT->GetFunction("gaus");
616
617 /* projection-x */
618 TH1D *histoVertexTimestamppx = histoVertexTimestamp->ProjectionX("histoVertexTimestamppx");
619 TH1D *histoDeltatTimestamppx = histoDeltatTimestamp->ProjectionX("histoDeltatTimestamppx");
620
07901ecf 621 /* check statistics */
622 if (histoVertexTimestamppx->Integral() < fgMinVertexIntegral ||
623 histoDeltatTimestamppx->Integral() < fgMinDeltatIntegral) {
624 fStatus = kLowStatistics;
625 return kFALSE;
626 }
627
3c8efc07 628 /* define mix and max time bin */
629 Int_t minBin = histoVertexTimestamppx->FindFirstBinAbove(0);
630 Int_t maxBin = histoVertexTimestamppx->FindLastBinAbove(0);
631 Float_t minTime = histoVertexTimestamppx->GetBinLowEdge(minBin);
632 Float_t maxTime = histoVertexTimestamppx->GetBinLowEdge(maxBin + 1);
633 AliInfo(Form("min/max time defined: %d < t < %d s [%d, %d]", (Int_t)minTime, (Int_t)maxTime, minBin, maxBin));
634
635 /* loop over time bins */
636 Int_t nPoints = 0;
637 Float_t time[fgkMaxNumberOfPoints], timeerr[fgkMaxNumberOfPoints];
638 Float_t vertexMean[fgkMaxNumberOfPoints], vertexMeanerr[fgkMaxNumberOfPoints];
639 Float_t vertexSigma[fgkMaxNumberOfPoints], vertexSigmaerr[fgkMaxNumberOfPoints];
640 Float_t timeZeroMean[fgkMaxNumberOfPoints], timeZeroMeanerr[fgkMaxNumberOfPoints];
641 Float_t timeZeroSigma[fgkMaxNumberOfPoints], timeZeroSigmaerr[fgkMaxNumberOfPoints];
642 Float_t averageTimeZero = 0.;
643 Float_t averageTimeSigma = 0.;
644 Float_t averageVertexSpread = 0.;
645 for (Int_t ibin = minBin; ibin <= maxBin; ibin++) {
646
647 /* define time window */
648 Int_t startBin = ibin;
649 Int_t endBin = ibin;
d501ffec 650 while(histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegralSample ||
651 histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegralSample) {
3c8efc07 652 if (endBin < maxBin) endBin++;
653 else if (startBin > minBin) startBin--;
654 else break;
655 }
d501ffec 656 if (histoVertexTimestamppx->Integral(startBin, endBin) < fgMinVertexIntegral ||
657 histoDeltatTimestamppx->Integral(startBin, endBin) < fgMinDeltatIntegral) continue;
3c8efc07 658 Float_t startTime = histoVertexTimestamppx->GetBinLowEdge(startBin);
659 Float_t endTime = histoVertexTimestamppx->GetBinLowEdge(endBin + 1);
660 Float_t vertexIntegral = histoVertexTimestamppx->Integral(startBin, endBin);
661 Float_t deltatIntegral = histoDeltatTimestamppx->Integral(startBin, endBin);
662 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));
663
664 /* projection-y */
665 TH1D *histoVertexTimestamppy = histoVertexTimestamp->ProjectionY("histoVertexTimestamppy", startBin, endBin);
666 TH1D *histoDeltatTimestamppy = histoDeltatTimestamp->ProjectionY("histoDeltatTimestamppy", startBin, endBin);
667
668 /* average time */
669 histoVertexTimestamppx->GetXaxis()->SetRange(startBin, endBin);
670 time[nPoints] = histoVertexTimestamppx->GetMean();
671 timeerr[nPoints] = histoVertexTimestamppx->GetMeanError();
672
673 /* fit vertex */
674 if (FitPeak(fitFunc, histoVertexTimestamppy, 10., 3., 3.) != 0) {
675 AliError("troubles fitting vertex, skip");
676 delete histoVertexTimestamppy;
677 delete histoDeltatTimestamppy;
678 continue;
679 }
680 vertexMean[nPoints] = fitFunc->GetParameter(1);
681 vertexMeanerr[nPoints] = fitFunc->GetParError(1);
682 vertexSigma[nPoints] = fitFunc->GetParameter(2);
683 vertexSigmaerr[nPoints] = fitFunc->GetParError(2);
684 averageVertexSpread += fitFunc->GetParameter(2);
685
686 /* fit time-zero */
687 if (FitPeak(fitFunc, histoDeltatTimestamppy, 500., 2., 1.) != 0) {
688 AliError("troubles fitting time-zero TRACKS, skip");
689 delete histoVertexTimestamppy;
690 delete histoDeltatTimestamppy;
691 continue;
692 }
693 timeZeroMean[nPoints] = fitFunc->GetParameter(1);
694 timeZeroMeanerr[nPoints] = fitFunc->GetParError(1);
695 timeZeroSigma[nPoints] = fitFunc->GetParameter(2);
696 timeZeroSigmaerr[nPoints] = fitFunc->GetParError(2);
697 averageTimeZero += fitFunc->GetParameter(1);
698 averageTimeSigma += fitFunc->GetParameter(2);
699
700 /* delete projection-y */
701 delete histoVertexTimestamppy;
702 delete histoDeltatTimestamppy;
703
704 /* increment n points */
705 nPoints++;
706
707 /* set current bin */
708 ibin = endBin;
709
710 } /* end of loop over time bins */
711
712 /* delete projection-x */
713 delete histoVertexTimestamppx;
714 delete histoDeltatTimestamppx;
715
716 /* check points */
717 if (nPoints <= 0) {
718 AliError("no measurement available, quit");
07901ecf 719 fStatus = kNoMeasurement;
3c8efc07 720 return kFALSE;
721 }
722 AliInfo(Form("average time-zero = %f", averageTimeZero / nPoints));
723 AliInfo(Form("average time-sigma = %f", averageTimeSigma / nPoints));
724 AliInfo(Form("average v-spread = %f", averageVertexSpread / nPoints));
725
726 /*** CREATE RUN PARAMS OBJECT ***/
727
728#if 0
729 /* get start time from GRP */
730 TGrid::Connect("alien://", 0, 0, "t");
731 AliCDBManager *cdb = AliCDBManager::Instance();
732 cdb->SetDefaultStorage("raw://");
733 cdb->SetRun(runNb);
734 AliGRPManager grp;
735 if (!grp.ReadGRPEntry()) {
736 AliError("error while reading GRP entry");
737 return kFALSE;
738 }
739 UInt_t startTimestamp = grp.GetGRPData()->GetTimeStart();
740 TTimeStamp ts;
741 ts = startTimestamp;
742 AliInfo(Form("got start time from GRP: %s", ts.AsString()));
743#endif
744
745 /* create arrays */
746 UInt_t timestamp[fgkMaxNumberOfPoints];
747 Float_t t0[fgkMaxNumberOfPoints];
748 Float_t tofReso[fgkMaxNumberOfPoints];
749 Float_t t0Spread[fgkMaxNumberOfPoints];
750 for (Int_t ipoint = 0; ipoint < nPoints; ipoint++) {
7ca3d6d6 751 timestamp[ipoint] = (UInt_t)time[ipoint] + startTimestamp;
3c8efc07 752 t0[ipoint] = timeZeroMean[ipoint];
753 tofReso[ipoint] = timeZeroSigma[ipoint];
754 t0Spread[ipoint] = vertexSigma[ipoint] / 2.99792457999999984e-02;
755 }
756 UInt_t run[1] = {runNb};
757 UInt_t runFirstPoint[1] = {0};
758 UInt_t runLastPoint[1] = {nPoints - 1};
759
760 /* create run params object */
761 AliTOFRunParams obj(nPoints, 1);
762 AliInfo(Form("create run params object for run %d with %d points", runNb, nPoints));
763 obj.SetTimestamp(timestamp);
764 obj.SetT0(t0);
765 obj.SetTOFResolution(tofReso);
766 obj.SetT0Spread(t0Spread);
767 obj.SetRunNb(run);
768 obj.SetRunFirstPoint(runFirstPoint);
769 obj.SetRunLastPoint(runLastPoint);
3fd7a5cd 770 obj.SetUseLHCClockPhase(useBPTX);
3c8efc07 771
772 /*** CREATE OCDB ENTRY ***/
773
774 if (!dbString) {
775 AliError("cannot store object because of NULL string");
07901ecf 776 fStatus = kStoreError;
3c8efc07 777 return kFALSE;
778 }
779
780 /* install run params object in OCDB */
781 AliCDBManager *cdb = AliCDBManager::Instance();
782 AliCDBStorage *sto = cdb->GetStorage(dbString);
783 if (!sto) {
784 AliError(Form("cannot get storage %s", dbString));
07901ecf 785 fStatus = kStoreError;
3c8efc07 786 return kFALSE;
787 }
788 AliCDBId id("TOF/Calib/RunParams", runNb, runNb);
789 AliCDBMetaData md;
790 md.SetResponsible("Roberto Preghenella");
791 md.SetComment("offline TOF run parameters");
792 md.SetAliRootVersion(gSystem->Getenv("ARVERSION"));
793 md.SetBeamPeriod(0);
794 if (!sto->Put(&obj, id, &md)) {
07901ecf 795 fStatus = kStoreError;
3c8efc07 796 AliError(Form("error while putting object in storage %s", dbString));
07901ecf 797 return kFALSE;
3c8efc07 798 }
799
800 /* success */
801 return kTRUE;
802}
803
804//_____________________________________________________________
805
806Int_t
807AliTOFAnalysisTaskCalibPass0::FitPeak(TF1 *fitFunc, TH1D *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
808{
809 /*
810 * fit peak
811 */
812
813 Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
814 Double_t fitMin = fitCent - nSigmaMin * startSigma;
815 Double_t fitMax = fitCent + nSigmaMax * startSigma;
816 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
817 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
818 fitFunc->SetParameter(1, fitCent);
819 fitFunc->SetParameter(2, startSigma);
820 Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
821 if (fitres != 0) return fitres;
822 /* refit with better range */
823 for (Int_t i = 0; i < 3; i++) {
824 fitCent = fitFunc->GetParameter(1);
825 fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
826 fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
827 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
828 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
829 fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
830 if (fitres != 0) return fitres;
831 }
832 return fitres;
833}