Fixed compilator warning: "conversion from string literal to 'char *' is deprecated...
[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 "AliTOFT0maker.h"
536031f2 47#include "AliPID.h"
f858b00e 48#include "AliLog.h"
536031f2 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"
f858b00e 56#include "AliTOFHeader.h"
536031f2 57
58ClassImp(AliTOFT0maker)
59
60//____________________________________________________________________________
5b4ed716 61AliTOFT0maker::AliTOFT0maker():
62 TObject(),
63 fT0TOF(NULL),
64 fPIDesd(NULL),
65 fExternalPIDFlag(kFALSE),
66 fTOFcalib(NULL),
67 fNoTOFT0(0),
2a258f40 68 fNmomBins(0),
5b4ed716 69 fTimeResolution(100),
70 fT0sigma(1000),
71 fHmapChannel(0),
72 fKmask(0),
2a258f40 73 fT0width(150.),
74 fT0spreadExt(-1.),
75 fT0fillExt(0)
5b4ed716 76{
77 // ctr
78 fCalculated[0] = 0;
79 fCalculated[1] = 0;
80 fCalculated[2] = 0;
81 fCalculated[3] = 0;
82
020dbfc2 83 fT0cur[0]=0.;
84 fT0cur[1]=0.;
85
5b4ed716 86 if(AliPID::ParticleMass(0) == 0) new AliPID();
87
88 fPIDesd = new AliESDpid();
89
2a258f40 90 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
91 SetTOFResponse();
92
93 fT0TOF = new AliTOFT0v1(fPIDesd);
94
5b4ed716 95}
96//____________________________________________________________________________
97AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
98 TObject(),
99 fT0TOF(NULL),
100 fPIDesd(externalPID),
101 fExternalPIDFlag(kTRUE),
102 fTOFcalib(tofCalib),
03bd764d 103 fNoTOFT0(0),
2a258f40 104 fNmomBins(0),
5b4ed716 105 fTimeResolution(100),
03bd764d 106 fT0sigma(1000),
107 fHmapChannel(0),
5b4ed716 108 fKmask(0),
2a258f40 109 fT0width(150.),
110 fT0spreadExt(-1.),
111 fT0fillExt(0)
536031f2 112{
6a4e212e 113 // ctr
6a4e212e 114 fCalculated[0] = 0;
115 fCalculated[1] = 0;
116 fCalculated[2] = 0;
03bd764d 117 fCalculated[3] = 0;
6a4e212e 118
020dbfc2 119 fT0cur[0]=0.;
120 fT0cur[1]=0.;
121
536031f2 122 if(AliPID::ParticleMass(0) == 0) new AliPID();
03bd764d 123
5b4ed716 124 if(!fPIDesd){
125 fPIDesd = new AliESDpid();
126 fExternalPIDFlag = kFALSE;
127 }
128
2a258f40 129 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
130 SetTOFResponse();
5b4ed716 131
2a258f40 132 fT0TOF = new AliTOFT0v1(fPIDesd);
8f589502 133
8f589502 134}
5b4ed716 135
8f589502 136//____________________________________________________________________________
536031f2 137AliTOFT0maker::~AliTOFT0maker()
138{
139 // dtor
5b4ed716 140
141 delete fT0TOF;
142 if (!fExternalPIDFlag) delete fPIDesd;
536031f2 143}
144//____________________________________________________________________________
5b4ed716 145Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
8f589502 146 //
147 // Remake TOF PID probabilities
148 //
2a258f40 149 Double_t t0tof[6];
536031f2 150
03bd764d 151 if(fKmask) ApplyMask(esd);
152
2a258f40 153 Double_t t0fill = 0.;
154
155 fPIDesd->GetTOFResponse().ResetT0info();
156
5b4ed716 157 /* get T0 spread from TOFcalib if available otherwise use default value */
158 if (fTOFcalib && esd) {
159 AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
160 if (runParams && runParams->GetTimestamp(0) != 0) {
161 Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
2a258f40 162 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
163 else{
164 SetT0FillWidth(t0spread);
165 t0fill = fT0fillExt;
166 }
167 }
168 else{
90beec49 169 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
2a258f40 170 t0fill = fT0fillExt;
5b4ed716 171 }
536031f2 172 }
bc4fdef6 173 else if(esd){
174 Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
175 if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;
176
177 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
178 else{
179 SetT0FillWidth(t0spread);
180 t0fill = fT0fillExt;
181 }
182 }
03bd764d 183
bf40b7fa 184 Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
185
5b4ed716 186 fT0TOF->Init(esd);
187 AliTOFT0v1* t0maker= fT0TOF;
5b4ed716 188
189 t0maker->DefineT0("all",1.5,3.0);
190 t0tof[0] = t0maker->GetResult(0);
191 t0tof[1] = t0maker->GetResult(1);
192 t0tof[2] = t0maker->GetResult(2);
193 t0tof[3] = t0maker->GetResult(3);
2a258f40 194 t0tof[4] = t0maker->GetResult(4);
195 t0tof[5] = t0maker->GetResult(5);
536031f2 196
8f589502 197 Float_t lT0Current=0.;
536031f2 198 fT0sigma=1000;
199
03bd764d 200// Int_t nrun = esd->GetRunNumber();
5b4ed716 201
03bd764d 202 t0time += t0fill;
8f589502 203
5b4ed716 204 Float_t sigmaFill = fT0width;
205
206 if(sigmaFill < 20) sigmaFill = 140;
8f589502 207
03bd764d 208 fCalculated[0]=-1000*t0tof[0]; // best t0
209 fCalculated[1]=1000*t0tof[1]; // sigma best t0
210 fCalculated[2] = t0fill; //t0 fill
211 fCalculated[3] = t0tof[2]; // n TOF tracks
212 fCalculated[4]=-1000*t0tof[0]; // TOF t0
213 fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
214 fCalculated[6]=sigmaFill; // sigma t0 fill
215 fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
216
bf40b7fa 217 if(fCalculated[7] > 30) thrGood = 10000000;
218
2a258f40 219 //statistics
220 fCalculated[8] = t0tof[4]; // real time in s
221 fCalculated[9] = t0tof[5]; // cpu time in s
222
bf40b7fa 223 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
6a4e212e 224 fT0sigma=fCalculated[1];
225 lT0Current=fCalculated[0];
536031f2 226 }
03bd764d 227 else{
228 fCalculated[4] = t0fill;
229 fCalculated[5] = sigmaFill;
230 }
231
5b4ed716 232 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
03bd764d 233 fT0sigma =1000;
234 fCalculated[4] = t0fill;
235 fCalculated[5] = sigmaFill;
236 }
536031f2 237
238 if(t0sigma < 1000){
239 if(fT0sigma < 1000){
240 Double_t w1 = 1./t0sigma/t0sigma;
6a4e212e 241 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
536031f2 242
243 Double_t wtot = w1+w2;
244
6a4e212e 245 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
536031f2 246 fT0sigma = TMath::Sqrt(1./wtot);
247 }
248 else{
8f589502 249 lT0Current=t0time;
536031f2 250 fT0sigma=t0sigma;
251 }
252 }
253
bf40b7fa 254 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
03bd764d 255 fCalculated[1]=fT0sigma;
256 fCalculated[0]=lT0Current;
257 }
258
259 if(fT0sigma >= 1000 || fNoTOFT0){
8f589502 260 lT0Current = t0fill;
03bd764d 261 fT0sigma = sigmaFill;
8f589502 262
6a4e212e 263 fCalculated[0] = t0fill;
03bd764d 264 fCalculated[1] = sigmaFill;
536031f2 265 }
266
5b4ed716 267 // T0 pt bin
721ed687 268 Float_t *t0values = new Float_t[fNmomBins];
269 Float_t *t0resolution = new Float_t[fNmomBins];
2a258f40 270 if(fCalculated[7] < 100){
271 for(Int_t i=0;i<fNmomBins;i++){
272 t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
273 t0tof[0] = t0maker->GetResult(0);
274 t0tof[1] = t0maker->GetResult(1);
275 t0tof[2] = t0maker->GetResult(2);
276 t0tof[3] = t0maker->GetResult(3);
2a258f40 277
278 Float_t t0bin =-1000*t0tof[0]; // best t0
279 Float_t t0binRes =1000*t0tof[1]; // sigma best t0
280
bf40b7fa 281 if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
2a258f40 282 // Ok T0
283 if(t0sigma < 1000){
284 Double_t w1 = 1./t0sigma/t0sigma;
285 Double_t w2 = 1./t0binRes/t0binRes;
286
287 Double_t wtot = w1+w2;
288
289 t0bin = (w1*t0time + w2*t0bin) / wtot;
290 t0binRes = TMath::Sqrt(1./wtot);
291 }
292 }
293 else{
294 t0bin = t0fill;
295 t0binRes = sigmaFill;
296 if(t0sigma < 1000){
297 t0bin = t0time;
298 t0binRes= t0sigma;
299 }
300 }
721ed687 301 t0values[i] = t0bin;
302 t0resolution[i] = t0binRes;
5b4ed716 303 }
2a258f40 304 }
305 else{
306 for(Int_t i=0;i<fNmomBins;i++){
721ed687 307 t0values[i] = lT0Current;
308 t0resolution[i] = fT0sigma;
5b4ed716 309 }
310 }
721ed687 311 for(Int_t i=0;i<fNmomBins;i++){
312 fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
313 fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
314 }
315
316 delete[] t0values;
317 delete[] t0resolution;
536031f2 318
6a4e212e 319 return fCalculated;
536031f2 320}
321//____________________________________________________________________________
5b4ed716 322Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
2a258f40 323 Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
03bd764d 324
2a258f40 325 fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
326 fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
327
5b4ed716 328 return fT0cur;
536031f2 329}
330//____________________________________________________________________________
5b4ed716 331void AliTOFT0maker::SetTOFResponse(){
2a258f40 332 fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
5b4ed716 333}
334//____________________________________________________________________________
335Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
5b4ed716 336 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
8f589502 337
5b4ed716 338 return sigma;
536031f2 339}
340//____________________________________________________________________________
5b4ed716 341void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
8f589502 342 //
5b4ed716 343 // Recalculate TOF PID probabilities
8f589502 344 //
345
90beec49 346// subtruct t0 for each track
347 Int_t ntracks = esd->GetNumberOfTracks();
536031f2 348
90beec49 349 while (ntracks--) {
350 AliESDtrack *t=esd->GetTrack(ntracks);
03bd764d 351
90beec49 352 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
5b4ed716 353
90beec49 354 Double_t time=t->GetTOFsignal();
355 Float_t p = t->GetP();
356
357 Double_t *t0=GetT0p(p);
358 time -= t0[0];
359 t->SetTOFsignal(time);
360 }
361
362 for(Int_t i=0;i<fNmomBins;i++){
363 fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
364 }
365
5b4ed716 366 //
536031f2 367}
03bd764d 368//____________________________________________________________________________
777fcfde 369void AliTOFT0maker::LoadChannelMap(const char *filename){
03bd764d 370 // Load the histo with the channel off map
371 TFile *f= new TFile(filename);
372 if(!f){
373 printf("Cannot open the channel map file (%s)\n",filename);
374 return;
375 }
376
377 fHmapChannel = (TH1F *) f->Get("hChEnabled");
378
379 if(!fHmapChannel){
380 printf("Cannot laod the channel map histo (from %s)\n",filename);
381 return;
382 }
383
384}
385//____________________________________________________________________________
386void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
387 // Switch off the disable channel
388 if(!fHmapChannel){
389 printf("Channel Map is not available\n");
390 return;
391 }
392
393 Int_t ntracks = esd->GetNumberOfTracks();
394
395 while (ntracks--) {
396 AliESDtrack *t=esd->GetTrack(ntracks);
397
398 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
399
400 Int_t chan = t->GetTOFCalChannel();
401
402 if(fHmapChannel->GetBinContent(chan) < 0.01){
403 t->ResetStatus(AliESDtrack::kTOFout);
404 }
405 }
406}
407
5b4ed716 408Float_t
409AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
410 //
411 // tune for MC data
412 //
413
414 Float_t TOFtimeResolutionDefault=80;
415
416 Float_t t0 = gRandom->Gaus(0.,fT0width);
417
418 Float_t extraSmearing = 0;
419
420 if(fTimeResolution > TOFtimeResolutionDefault){
421 extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
422 }
423
424 // subtruct t0 for each track
425 Int_t ntracks = esd->GetNumberOfTracks();
426
427 while (ntracks--) {
428 AliESDtrack *t=esd->GetTrack(ntracks);
429
430 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
431
432 /* check if channel is enabled */
bc4fdef6 433 if (fTOFcalib){
434 if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
435 /* reset TOF status */
436 t->ResetStatus(AliESDtrack::kTOFin);
437 t->ResetStatus(AliESDtrack::kTOFout);
438 t->ResetStatus(AliESDtrack::kTOFmismatch);
439 t->ResetStatus(AliESDtrack::kTOFpid);
440 }
5b4ed716 441 }
442
443 Double_t time=t->GetTOFsignal();
444
445 time += t0;
446
447 if(extraSmearing>0){
448 Float_t smearing = gRandom->Gaus(0.,extraSmearing);
449 time += smearing;
450 }
451
452 t->SetTOFsignal(time);
453 }
454 //
455 return t0;
456}
f858b00e 457//_________________________________________________________________________
458void AliTOFT0maker::WriteInESD(AliESDEvent *esd){
459 //
460 // Write t0Gen, t0ResGen, nt0;
461 // t0resESD[0:nt0], it0ESD[0:nt0]
462 // in the AliESDEvent
463 //
464 Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();
465
466 Int_t nt0=0;
467 Float_t *t0 = new Float_t[nMomBins];
468 Float_t *t0res = new Float_t[nMomBins];
469 Int_t *multT0 = new Int_t[nMomBins];
470
471 for(Int_t i=0;i<nMomBins;i++){
472// printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
473 multT0[i]=0;
474 Bool_t kNewT0=kTRUE;
475 for(Int_t j=0;j < nt0;j++){
476 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
477 kNewT0=kFALSE;
478 multT0[j]++;
479 j=nMomBins*10;
480 }
481 }
482 if(kNewT0){
483 t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
484 t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
485 nt0++;
486 }
487 }
488
489 Int_t iMultT0=0,nmult=0;
490 for(Int_t j=0;j < nt0;j++){ // find the most frequent
491 if(multT0[j] > nmult){
492 iMultT0 = j;
493 nmult = multT0[j];
494 }
495 }
496
497 Float_t *t0ESD = new Float_t[nMomBins];
498 Float_t *t0resESD = new Float_t[nMomBins];
499 Int_t *it0ESD = new Int_t[nMomBins];
500
501 Float_t t0Gen,t0ResGen;
502 t0Gen = t0[iMultT0];
503 t0ResGen = t0res[iMultT0];
504 nt0=0;
505 // printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
506 for(Int_t i=0;i<nMomBins;i++){
507 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
508 t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
509 t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
510 it0ESD[nt0]=i;
511// printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
512 nt0++;
513 }
514 }
515
516 // Write t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent
517
518 AliTOFHeader *tofHeader =
519 new AliTOFHeader(t0Gen,t0ResGen,nt0,
c85a1b24 520 t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
f858b00e 521
522 esd->SetTOFHeader(tofHeader);
523
c85a1b24 524 AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
525 AliDebug(1,Form("%d ",nt0));
f858b00e 526 for (Int_t ii=0; ii<nt0; ii++)
c85a1b24 527 AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
f858b00e 528
529 delete[] t0;
530 t0 = NULL;
531 delete[] t0res;
532 t0res = NULL;
533 delete[] multT0;
534 multT0 = NULL;
535 delete[] t0ESD;
536 t0ESD = NULL;
537 delete[] t0resESD;
538 t0resESD = NULL;
539 delete[] it0ESD;
540 it0ESD = NULL;
541
542}