Fix for report #72089 changes to AliTOFPIDResponse to be commited in STEER (Roberto)
[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.cxx
CommitLineData
8f589502 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 **************************************************************************/
6a4e212e 15/* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
8f589502 16
17/////////////////////////////////////////////////////////////////////////////
18// //
19// This class contains the basic functions for the time zero //
6a4e212e 20// evaluation with TOF detector informations. //
21// Use case in an analysis task: //
22// //
23// Create the object in the task constructor (fTOFmaker is a private var) //
5b4ed716 24// AliESDpid *extPID=new AliESDpid(); //
25// fTOFmaker = new AliTOFT0maker(extPID); //
26// fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res //
6a4e212e 27// 115 ps is the TOF default resolution value //
28// //
29// Use the RemakePID method in the task::Exec //
30// Double_t* calcolot0; //
31// calcolot0=fTOFmaker->RemakePID(fESD); //
32// //calcolot0[0] = calculated event time //
33// //calcolot0[1] = event time time resolution //
34// //calcolot0[2] = average event time for the current fill //
03bd764d 35// //calcolot0[3] = tracks at TOF //
36// //calcolot0[4] = calculated event time (only TOF) //
37// //calcolot0[5] = event time time resolution (only TOF) //
38// //calcolot0[6] = sigma t0 fill //
39// //calcolot0[7] = tracks at TOF really used in tht algorithm //
6a4e212e 40// //
41// Let consider that: //
42// - the PIF is automatically recalculated with the event time subtrction //
8f589502 43// //
44/////////////////////////////////////////////////////////////////////////////
45
536031f2 46#include "AliTOFT0v1.h"
47#include "AliTOFT0maker.h"
536031f2 48#include "AliPID.h"
49#include "AliESDpid.h"
03bd764d 50#include "AliESDEvent.h"
51#include "TFile.h"
52#include "TH1F.h"
5b4ed716 53#include "AliTOFcalib.h"
54#include "AliTOFRunParams.h"
55#include "TRandom.h"
536031f2 56
57ClassImp(AliTOFT0maker)
58
59//____________________________________________________________________________
5b4ed716 60AliTOFT0maker::AliTOFT0maker():
61 TObject(),
62 fT0TOF(NULL),
63 fPIDesd(NULL),
64 fExternalPIDFlag(kFALSE),
65 fTOFcalib(NULL),
66 fNoTOFT0(0),
67 fTimeResolution(100),
68 fT0sigma(1000),
69 fHmapChannel(0),
70 fKmask(0),
71 fT0width(150.)
72{
73 // ctr
74 fCalculated[0] = 0;
75 fCalculated[1] = 0;
76 fCalculated[2] = 0;
77 fCalculated[3] = 0;
78
79 fT0TOF = new AliTOFT0v1();
80 if(AliPID::ParticleMass(0) == 0) new AliPID();
81
82 fPIDesd = new AliESDpid();
83
84 fPtCutMin[0] = 0.3;
85 fPtCutMin[1] = 0.5;
86 fPtCutMin[2] = 0.6;
87 fPtCutMin[3] = 0.7;
88 fPtCutMin[4] = 0.8;
89 fPtCutMin[5] = 0.9;
90 fPtCutMin[6] = 1;
91 fPtCutMin[7] = 1.2;
92 fPtCutMin[8] = 1.5;
93 fPtCutMin[9] = 2;
94
95 fPtCutMax[0] = 0.5;
96 fPtCutMax[1] = 0.6;
97 fPtCutMax[2] = 0.7;
98 fPtCutMax[3] = 0.8;
99 fPtCutMax[4] = 0.9;
100 fPtCutMax[5] = 1;
101 fPtCutMax[6] = 1.2;
102 fPtCutMax[7] = 1.5;
103 fPtCutMax[8] = 2;
104 fPtCutMax[9] = 3;
105
106 /* init arrays */
107 for (Int_t i = 0; i < 10; i++) {
108 fT0pt[i] = 0.;
109 fT0ptSigma[i] = 0.;
110 }
111}
112//____________________________________________________________________________
113AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
114 TObject(),
115 fT0TOF(NULL),
116 fPIDesd(externalPID),
117 fExternalPIDFlag(kTRUE),
118 fTOFcalib(tofCalib),
03bd764d 119 fNoTOFT0(0),
5b4ed716 120 fTimeResolution(100),
03bd764d 121 fT0sigma(1000),
122 fHmapChannel(0),
5b4ed716 123 fKmask(0),
124 fT0width(150.)
536031f2 125{
6a4e212e 126 // ctr
6a4e212e 127 fCalculated[0] = 0;
128 fCalculated[1] = 0;
129 fCalculated[2] = 0;
03bd764d 130 fCalculated[3] = 0;
6a4e212e 131
5b4ed716 132 fT0TOF = new AliTOFT0v1();
536031f2 133 if(AliPID::ParticleMass(0) == 0) new AliPID();
03bd764d 134
5b4ed716 135 if(!fPIDesd){
136 fPIDesd = new AliESDpid();
137 fExternalPIDFlag = kFALSE;
138 }
139
140 fPtCutMin[0] = 0.3;
141 fPtCutMin[1] = 0.5;
142 fPtCutMin[2] = 0.6;
143 fPtCutMin[3] = 0.7;
144 fPtCutMin[4] = 0.8;
145 fPtCutMin[5] = 0.9;
146 fPtCutMin[6] = 1;
147 fPtCutMin[7] = 1.2;
148 fPtCutMin[8] = 1.5;
149 fPtCutMin[9] = 2;
150
151 fPtCutMax[0] = 0.5;
152 fPtCutMax[1] = 0.6;
153 fPtCutMax[2] = 0.7;
154 fPtCutMax[3] = 0.8;
155 fPtCutMax[4] = 0.9;
156 fPtCutMax[5] = 1;
157 fPtCutMax[6] = 1.2;
158 fPtCutMax[7] = 1.5;
159 fPtCutMax[8] = 2;
160 fPtCutMax[9] = 3;
161
162 /* init arrays */
163 for (Int_t i = 0; i < 10; i++) {
164 fT0pt[i] = 0.;
165 fT0ptSigma[i] = 0.;
166 }
536031f2 167}
5b4ed716 168
169/* copy-constructor and operator= suppressed
170
171
536031f2 172//____________________________________________________________________________
8f589502 173AliTOFT0maker::AliTOFT0maker(const AliTOFT0maker & t) :
5b4ed716 174 TObject(t),
175 fT0TOF(t.fT0TOF),
176 fPIDESD(t.fPIDESD),
03bd764d 177 fNoTOFT0(t.fNoTOFT0),
8f589502 178 fTimeResolution(t.fTimeResolution),
03bd764d 179 fT0sigma(t.fT0sigma),
180 fHmapChannel(t.fHmapChannel),
181 fKmask(t.fKmask)
5b4ed716 182 {
03bd764d 183 // copy ctr
8f589502 184}
185
186//____________________________________________________________________________
187AliTOFT0maker& AliTOFT0maker::operator=(const AliTOFT0maker &t)
188{
6a4e212e 189 //
8f589502 190 // assign. operator
191 //
192
193 if (this == &t)
194 return *this;
5b4ed716 195
196 TObject::operator=(t);
8f589502 197 fTimeResolution = t.fTimeResolution;
198 fT0sigma = t.fT0sigma;
199
200 return *this;
201}
5b4ed716 202
203*/
204
8f589502 205//____________________________________________________________________________
536031f2 206AliTOFT0maker::~AliTOFT0maker()
207{
208 // dtor
5b4ed716 209
210 delete fT0TOF;
211 if (!fExternalPIDFlag) delete fPIDesd;
536031f2 212}
213//____________________________________________________________________________
5b4ed716 214Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
8f589502 215 //
216 // Remake TOF PID probabilities
217 //
218
5b4ed716 219 Double_t t0tof[4];
536031f2 220
03bd764d 221 if(fKmask) ApplyMask(esd);
222
5b4ed716 223 /* get T0 spread from TOFcalib if available otherwise use default value */
224 if (fTOFcalib && esd) {
225 AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
226 if (runParams && runParams->GetTimestamp(0) != 0) {
227 Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
228 SetT0FillWidth(t0spread);
229 }
536031f2 230 }
03bd764d 231
5b4ed716 232 fT0TOF->Init(esd);
233 AliTOFT0v1* t0maker= fT0TOF;
234 t0maker->SetTimeResolution(fTimeResolution*1e-12);
235
236 t0maker->DefineT0("all",1.5,3.0);
237 t0tof[0] = t0maker->GetResult(0);
238 t0tof[1] = t0maker->GetResult(1);
239 t0tof[2] = t0maker->GetResult(2);
240 t0tof[3] = t0maker->GetResult(3);
536031f2 241
8f589502 242 Float_t lT0Current=0.;
536031f2 243 fT0sigma=1000;
244
03bd764d 245// Int_t nrun = esd->GetRunNumber();
5b4ed716 246 Double_t t0fill = 0.;
247
03bd764d 248 t0time += t0fill;
8f589502 249
5b4ed716 250 Float_t sigmaFill = fT0width;
251
252 if(sigmaFill < 20) sigmaFill = 140;
8f589502 253
03bd764d 254 fCalculated[0]=-1000*t0tof[0]; // best t0
255 fCalculated[1]=1000*t0tof[1]; // sigma best t0
256 fCalculated[2] = t0fill; //t0 fill
257 fCalculated[3] = t0tof[2]; // n TOF tracks
258 fCalculated[4]=-1000*t0tof[0]; // TOF t0
259 fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
260 fCalculated[6]=sigmaFill; // sigma t0 fill
261 fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
262
5b4ed716 263 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500 && fCalculated[1] < fTimeResolution*1.2){
6a4e212e 264 fT0sigma=fCalculated[1];
265 lT0Current=fCalculated[0];
536031f2 266 }
03bd764d 267 else{
268 fCalculated[4] = t0fill;
269 fCalculated[5] = sigmaFill;
270 }
271
5b4ed716 272 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
03bd764d 273 fT0sigma =1000;
274 fCalculated[4] = t0fill;
275 fCalculated[5] = sigmaFill;
276 }
536031f2 277
278 if(t0sigma < 1000){
279 if(fT0sigma < 1000){
280 Double_t w1 = 1./t0sigma/t0sigma;
6a4e212e 281 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
536031f2 282
283 Double_t wtot = w1+w2;
284
6a4e212e 285 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
536031f2 286 fT0sigma = TMath::Sqrt(1./wtot);
287 }
288 else{
8f589502 289 lT0Current=t0time;
536031f2 290 fT0sigma=t0sigma;
291 }
292 }
293
03bd764d 294 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < 500){
295 fCalculated[1]=fT0sigma;
296 fCalculated[0]=lT0Current;
297 }
298
299 if(fT0sigma >= 1000 || fNoTOFT0){
8f589502 300 lT0Current = t0fill;
03bd764d 301 fT0sigma = sigmaFill;
8f589502 302
6a4e212e 303 fCalculated[0] = t0fill;
03bd764d 304 fCalculated[1] = sigmaFill;
536031f2 305 }
306
5b4ed716 307 // T0 pt bin
308 for(Int_t i=0;i<10;i++){
309 t0maker->DefineT0("all",fPtCutMin[i],fPtCutMax[i]);
310 t0tof[0] = t0maker->GetResult(0);
311 t0tof[1] = t0maker->GetResult(1);
312 t0tof[2] = t0maker->GetResult(2);
313 t0tof[3] = t0maker->GetResult(3);
314 fT0pt[i] =-1000*t0tof[0]; // best t0
315 fT0ptSigma[i] =1000*t0tof[1]; // sigma best t0
316
317 if(fT0ptSigma[i] < sigmaFill && fT0ptSigma[i] < fTimeResolution * 1.2 && TMath::Abs(fT0pt[i] - t0fill) < 500){
318 // Ok T0
319 }
320 else{
321 fT0pt[i] = t0fill;
322 fT0ptSigma[i] = sigmaFill;
323 }
324 }
325 //----
326 SetTOFResponse();
536031f2 327
6a4e212e 328 return fCalculated;
536031f2 329}
330//____________________________________________________________________________
5b4ed716 331Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
332 Int_t i=0;
333 while(p > fPtCutMin[i] && i < 10) i++;
334 if(i > 0) i--;
03bd764d 335
5b4ed716 336 fT0cur[0] = fT0pt[i];
337 fT0cur[1] = fT0ptSigma[i];
338 return fT0cur;
536031f2 339}
340//____________________________________________________________________________
5b4ed716 341void AliTOFT0maker::SetTOFResponse(){
342 fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
343}
344//____________________________________________________________________________
345Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
346 Double_t *sigmaT0 = GetT0p(mom);
347 fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(sigmaT0[1]*sigmaT0[1] + fTimeResolution*fTimeResolution));
348 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
349 fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
8f589502 350
5b4ed716 351 return sigma;
536031f2 352}
353//____________________________________________________________________________
5b4ed716 354void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
8f589502 355 //
5b4ed716 356 // Recalculate TOF PID probabilities
8f589502 357 //
358
5b4ed716 359 // subtruct t0 for each track
360 Int_t ntracks = esd->GetNumberOfTracks();
536031f2 361
5b4ed716 362 while (ntracks--) {
363 AliESDtrack *t=esd->GetTrack(ntracks);
03bd764d 364
5b4ed716 365 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
366
367 Double_t time=t->GetTOFsignal();
368 Float_t p = t->GetP();
03bd764d 369
5b4ed716 370 Double_t *t0=GetT0p(p);
371 time -= t0[0];
372 t->SetTOFsignal(time);
03bd764d 373 }
5b4ed716 374 //
536031f2 375}
03bd764d 376//____________________________________________________________________________
377void AliTOFT0maker::LoadChannelMap(char *filename){
378 // Load the histo with the channel off map
379 TFile *f= new TFile(filename);
380 if(!f){
381 printf("Cannot open the channel map file (%s)\n",filename);
382 return;
383 }
384
385 fHmapChannel = (TH1F *) f->Get("hChEnabled");
386
387 if(!fHmapChannel){
388 printf("Cannot laod the channel map histo (from %s)\n",filename);
389 return;
390 }
391
392}
393//____________________________________________________________________________
394void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
395 // Switch off the disable channel
396 if(!fHmapChannel){
397 printf("Channel Map is not available\n");
398 return;
399 }
400
401 Int_t ntracks = esd->GetNumberOfTracks();
402
403 while (ntracks--) {
404 AliESDtrack *t=esd->GetTrack(ntracks);
405
406 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
407
408 Int_t chan = t->GetTOFCalChannel();
409
410 if(fHmapChannel->GetBinContent(chan) < 0.01){
411 t->ResetStatus(AliESDtrack::kTOFout);
412 }
413 }
414}
415
5b4ed716 416Float_t
417AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
418 //
419 // tune for MC data
420 //
421
422 Float_t TOFtimeResolutionDefault=80;
423
424 Float_t t0 = gRandom->Gaus(0.,fT0width);
425
426 Float_t extraSmearing = 0;
427
428 if(fTimeResolution > TOFtimeResolutionDefault){
429 extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
430 }
431
432 // subtruct t0 for each track
433 Int_t ntracks = esd->GetNumberOfTracks();
434
435 while (ntracks--) {
436 AliESDtrack *t=esd->GetTrack(ntracks);
437
438 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
439
440 /* check if channel is enabled */
441 if (fTOFcalib && !fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
442 /* reset TOF status */
443 t->ResetStatus(AliESDtrack::kTOFin);
444 t->ResetStatus(AliESDtrack::kTOFout);
445 t->ResetStatus(AliESDtrack::kTOFrefit);
446 t->ResetStatus(AliESDtrack::kTOFpid);
447 }
448
449 Double_t time=t->GetTOFsignal();
450
451 time += t0;
452
453 if(extraSmearing>0){
454 Float_t smearing = gRandom->Gaus(0.,extraSmearing);
455 time += smearing;
456 }
457
458 t->SetTOFsignal(time);
459 }
460 //
461 return t0;
462}