Fix for Savannah bug 74286
[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),
2a258f40 67 fNmomBins(0),
5b4ed716 68 fTimeResolution(100),
69 fT0sigma(1000),
70 fHmapChannel(0),
71 fKmask(0),
2a258f40 72 fT0width(150.),
73 fT0spreadExt(-1.),
74 fT0fillExt(0)
5b4ed716 75{
76 // ctr
77 fCalculated[0] = 0;
78 fCalculated[1] = 0;
79 fCalculated[2] = 0;
80 fCalculated[3] = 0;
81
5b4ed716 82 if(AliPID::ParticleMass(0) == 0) new AliPID();
83
84 fPIDesd = new AliESDpid();
85
2a258f40 86 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
87 SetTOFResponse();
88
89 fT0TOF = new AliTOFT0v1(fPIDesd);
90
5b4ed716 91}
92//____________________________________________________________________________
93AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
94 TObject(),
95 fT0TOF(NULL),
96 fPIDesd(externalPID),
97 fExternalPIDFlag(kTRUE),
98 fTOFcalib(tofCalib),
03bd764d 99 fNoTOFT0(0),
2a258f40 100 fNmomBins(0),
5b4ed716 101 fTimeResolution(100),
03bd764d 102 fT0sigma(1000),
103 fHmapChannel(0),
5b4ed716 104 fKmask(0),
2a258f40 105 fT0width(150.),
106 fT0spreadExt(-1.),
107 fT0fillExt(0)
536031f2 108{
6a4e212e 109 // ctr
6a4e212e 110 fCalculated[0] = 0;
111 fCalculated[1] = 0;
112 fCalculated[2] = 0;
03bd764d 113 fCalculated[3] = 0;
6a4e212e 114
536031f2 115 if(AliPID::ParticleMass(0) == 0) new AliPID();
03bd764d 116
5b4ed716 117 if(!fPIDesd){
118 fPIDesd = new AliESDpid();
119 fExternalPIDFlag = kFALSE;
120 }
121
2a258f40 122 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
123 SetTOFResponse();
5b4ed716 124
2a258f40 125 fT0TOF = new AliTOFT0v1(fPIDesd);
8f589502 126
8f589502 127}
5b4ed716 128
8f589502 129//____________________________________________________________________________
536031f2 130AliTOFT0maker::~AliTOFT0maker()
131{
132 // dtor
5b4ed716 133
134 delete fT0TOF;
135 if (!fExternalPIDFlag) delete fPIDesd;
536031f2 136}
137//____________________________________________________________________________
5b4ed716 138Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
8f589502 139 //
140 // Remake TOF PID probabilities
141 //
142
2a258f40 143 Double_t t0tof[6];
536031f2 144
03bd764d 145 if(fKmask) ApplyMask(esd);
146
2a258f40 147 Double_t t0fill = 0.;
148
149 fPIDesd->GetTOFResponse().ResetT0info();
150
5b4ed716 151 /* get T0 spread from TOFcalib if available otherwise use default value */
152 if (fTOFcalib && esd) {
153 AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
154 if (runParams && runParams->GetTimestamp(0) != 0) {
155 Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
2a258f40 156 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
157 else{
158 SetT0FillWidth(t0spread);
159 t0fill = fT0fillExt;
160 }
161 }
162 else{
90beec49 163 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
2a258f40 164 t0fill = fT0fillExt;
5b4ed716 165 }
536031f2 166 }
03bd764d 167
5b4ed716 168 fT0TOF->Init(esd);
169 AliTOFT0v1* t0maker= fT0TOF;
5b4ed716 170
171 t0maker->DefineT0("all",1.5,3.0);
172 t0tof[0] = t0maker->GetResult(0);
173 t0tof[1] = t0maker->GetResult(1);
174 t0tof[2] = t0maker->GetResult(2);
175 t0tof[3] = t0maker->GetResult(3);
2a258f40 176 t0tof[4] = t0maker->GetResult(4);
177 t0tof[5] = t0maker->GetResult(5);
536031f2 178
8f589502 179 Float_t lT0Current=0.;
536031f2 180 fT0sigma=1000;
181
03bd764d 182// Int_t nrun = esd->GetRunNumber();
5b4ed716 183
03bd764d 184 t0time += t0fill;
8f589502 185
5b4ed716 186 Float_t sigmaFill = fT0width;
187
188 if(sigmaFill < 20) sigmaFill = 140;
8f589502 189
03bd764d 190 fCalculated[0]=-1000*t0tof[0]; // best t0
191 fCalculated[1]=1000*t0tof[1]; // sigma best t0
192 fCalculated[2] = t0fill; //t0 fill
193 fCalculated[3] = t0tof[2]; // n TOF tracks
194 fCalculated[4]=-1000*t0tof[0]; // TOF t0
195 fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
196 fCalculated[6]=sigmaFill; // sigma t0 fill
197 fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
198
2a258f40 199 //statistics
200 fCalculated[8] = t0tof[4]; // real time in s
201 fCalculated[9] = t0tof[5]; // cpu time in s
202
5b4ed716 203 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500 && fCalculated[1] < fTimeResolution*1.2){
6a4e212e 204 fT0sigma=fCalculated[1];
205 lT0Current=fCalculated[0];
536031f2 206 }
03bd764d 207 else{
208 fCalculated[4] = t0fill;
209 fCalculated[5] = sigmaFill;
210 }
211
5b4ed716 212 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
03bd764d 213 fT0sigma =1000;
214 fCalculated[4] = t0fill;
215 fCalculated[5] = sigmaFill;
216 }
536031f2 217
218 if(t0sigma < 1000){
219 if(fT0sigma < 1000){
220 Double_t w1 = 1./t0sigma/t0sigma;
6a4e212e 221 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
536031f2 222
223 Double_t wtot = w1+w2;
224
6a4e212e 225 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
536031f2 226 fT0sigma = TMath::Sqrt(1./wtot);
227 }
228 else{
8f589502 229 lT0Current=t0time;
536031f2 230 fT0sigma=t0sigma;
231 }
232 }
233
03bd764d 234 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < 500){
235 fCalculated[1]=fT0sigma;
236 fCalculated[0]=lT0Current;
237 }
238
239 if(fT0sigma >= 1000 || fNoTOFT0){
8f589502 240 lT0Current = t0fill;
03bd764d 241 fT0sigma = sigmaFill;
8f589502 242
6a4e212e 243 fCalculated[0] = t0fill;
03bd764d 244 fCalculated[1] = sigmaFill;
536031f2 245 }
246
5b4ed716 247 // T0 pt bin
2a258f40 248 if(fCalculated[7] < 100){
249 for(Int_t i=0;i<fNmomBins;i++){
250 t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
251 t0tof[0] = t0maker->GetResult(0);
252 t0tof[1] = t0maker->GetResult(1);
253 t0tof[2] = t0maker->GetResult(2);
254 t0tof[3] = t0maker->GetResult(3);
255
256
257 Float_t t0bin =-1000*t0tof[0]; // best t0
258 Float_t t0binRes =1000*t0tof[1]; // sigma best t0
259
260 if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < 500){
261 // Ok T0
262 if(t0sigma < 1000){
263 Double_t w1 = 1./t0sigma/t0sigma;
264 Double_t w2 = 1./t0binRes/t0binRes;
265
266 Double_t wtot = w1+w2;
267
268 t0bin = (w1*t0time + w2*t0bin) / wtot;
269 t0binRes = TMath::Sqrt(1./wtot);
270 }
271 }
272 else{
273 t0bin = t0fill;
274 t0binRes = sigmaFill;
275 if(t0sigma < 1000){
276 t0bin = t0time;
277 t0binRes= t0sigma;
278 }
279 }
280 fPIDesd->GetTOFResponse().SetT0bin(i,t0bin);
281 fPIDesd->GetTOFResponse().SetT0binRes(i,t0binRes);
5b4ed716 282 }
2a258f40 283 }
284 else{
285 for(Int_t i=0;i<fNmomBins;i++){
286 fPIDesd->GetTOFResponse().SetT0bin(i,lT0Current);
287 fPIDesd->GetTOFResponse().SetT0binRes(i,fT0sigma);
5b4ed716 288 }
289 }
536031f2 290
6a4e212e 291 return fCalculated;
536031f2 292}
293//____________________________________________________________________________
5b4ed716 294Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
2a258f40 295 Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
03bd764d 296
2a258f40 297 fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
298 fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
299
5b4ed716 300 return fT0cur;
536031f2 301}
302//____________________________________________________________________________
5b4ed716 303void AliTOFT0maker::SetTOFResponse(){
2a258f40 304 fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
5b4ed716 305}
306//____________________________________________________________________________
307Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
5b4ed716 308 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
8f589502 309
5b4ed716 310 return sigma;
536031f2 311}
312//____________________________________________________________________________
5b4ed716 313void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
8f589502 314 //
5b4ed716 315 // Recalculate TOF PID probabilities
8f589502 316 //
317
90beec49 318// subtruct t0 for each track
319 Int_t ntracks = esd->GetNumberOfTracks();
536031f2 320
90beec49 321 while (ntracks--) {
322 AliESDtrack *t=esd->GetTrack(ntracks);
03bd764d 323
90beec49 324 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
5b4ed716 325
90beec49 326 Double_t time=t->GetTOFsignal();
327 Float_t p = t->GetP();
328
329 Double_t *t0=GetT0p(p);
330 time -= t0[0];
331 t->SetTOFsignal(time);
332 }
333
334 for(Int_t i=0;i<fNmomBins;i++){
335 fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
336 }
337
5b4ed716 338 //
536031f2 339}
03bd764d 340//____________________________________________________________________________
341void AliTOFT0maker::LoadChannelMap(char *filename){
342 // Load the histo with the channel off map
343 TFile *f= new TFile(filename);
344 if(!f){
345 printf("Cannot open the channel map file (%s)\n",filename);
346 return;
347 }
348
349 fHmapChannel = (TH1F *) f->Get("hChEnabled");
350
351 if(!fHmapChannel){
352 printf("Cannot laod the channel map histo (from %s)\n",filename);
353 return;
354 }
355
356}
357//____________________________________________________________________________
358void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
359 // Switch off the disable channel
360 if(!fHmapChannel){
361 printf("Channel Map is not available\n");
362 return;
363 }
364
365 Int_t ntracks = esd->GetNumberOfTracks();
366
367 while (ntracks--) {
368 AliESDtrack *t=esd->GetTrack(ntracks);
369
370 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
371
372 Int_t chan = t->GetTOFCalChannel();
373
374 if(fHmapChannel->GetBinContent(chan) < 0.01){
375 t->ResetStatus(AliESDtrack::kTOFout);
376 }
377 }
378}
379
5b4ed716 380Float_t
381AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
382 //
383 // tune for MC data
384 //
385
386 Float_t TOFtimeResolutionDefault=80;
387
388 Float_t t0 = gRandom->Gaus(0.,fT0width);
389
390 Float_t extraSmearing = 0;
391
392 if(fTimeResolution > TOFtimeResolutionDefault){
393 extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
394 }
395
396 // subtruct t0 for each track
397 Int_t ntracks = esd->GetNumberOfTracks();
398
399 while (ntracks--) {
400 AliESDtrack *t=esd->GetTrack(ntracks);
401
402 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
403
404 /* check if channel is enabled */
405 if (fTOFcalib && !fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
406 /* reset TOF status */
407 t->ResetStatus(AliESDtrack::kTOFin);
408 t->ResetStatus(AliESDtrack::kTOFout);
90beec49 409 t->ResetStatus(AliESDtrack::kTOFmismatch);
5b4ed716 410 t->ResetStatus(AliESDtrack::kTOFpid);
411 }
412
413 Double_t time=t->GetTOFsignal();
414
415 time += t0;
416
417 if(extraSmearing>0){
418 Float_t smearing = gRandom->Gaus(0.,extraSmearing);
419 time += smearing;
420 }
421
422 t->SetTOFsignal(time);
423 }
424 //
425 return t0;
426}