1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 /* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
17 /////////////////////////////////////////////////////////////////////////////
19 // This class contains the basic functions for the time zero //
20 // evaluation with TOF detector informations. //
21 // Use case in an analysis task: //
23 // Create the object in the task constructor (fTOFmaker is a private var) //
24 // AliESDpid *extPID=new AliESDpid(); //
25 // fTOFmaker = new AliTOFT0maker(extPID); //
26 // fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res //
27 // 115 ps is the TOF default resolution value //
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 //
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 //
41 // Let consider that: //
42 // - the PIF is automatically recalculated with the event time subtrction //
44 /////////////////////////////////////////////////////////////////////////////
46 #include "AliTOFT0v1.h"
47 #include "AliTOFT0maker.h"
50 #include "AliESDpid.h"
51 #include "AliESDEvent.h"
54 #include "AliTOFcalib.h"
55 #include "AliTOFRunParams.h"
57 #include "AliTOFHeader.h"
59 ClassImp(AliTOFT0maker)
61 //____________________________________________________________________________
62 AliTOFT0maker::AliTOFT0maker():
66 fExternalPIDFlag(kFALSE),
88 if(AliPID::ParticleMass(0) == 0) new AliPID();
90 fPIDesd = new AliESDpid();
92 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
95 fT0TOF = new AliTOFT0v1(fPIDesd);
98 //____________________________________________________________________________
99 AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
102 fPIDesd(externalPID),
103 fExternalPIDFlag(kTRUE),
107 fTimeResolution(100),
125 if(AliPID::ParticleMass(0) == 0) new AliPID();
128 fPIDesd = new AliESDpid();
129 fExternalPIDFlag = kFALSE;
132 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
135 fT0TOF = new AliTOFT0v1(fPIDesd);
139 //____________________________________________________________________________
140 AliTOFT0maker::~AliTOFT0maker()
144 if (!fExternalPIDFlag) delete fPIDesd;
146 //____________________________________________________________________________
147 Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
149 // Remake TOF PID probabilities
153 if(fKmask) ApplyMask(esd);
155 Double_t t0fill = 0.;
157 fPIDesd->GetTOFResponse().ResetT0info();
159 /* get T0 spread from TOFcalib if available otherwise use default value */
160 if (fTOFcalib && esd) {
161 AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
162 if (runParams && runParams->GetTimestamp(0) != 0) {
163 Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
164 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
166 SetT0FillWidth(t0spread);
171 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
176 Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
177 if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;
179 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
181 SetT0FillWidth(t0spread);
186 Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
190 AliTOFT0v1* t0maker = fT0TOF;
191 if (fTOFT0algorithm==2) t0maker->SetOptimization(kTRUE);
192 t0maker->DefineT0("all",1.5,3.0);
193 t0tof[0] = t0maker->GetResult(0);
194 t0tof[1] = t0maker->GetResult(1);
195 t0tof[2] = t0maker->GetResult(2);
196 t0tof[3] = t0maker->GetResult(3);
197 t0tof[4] = t0maker->GetResult(4);
198 t0tof[5] = t0maker->GetResult(5);
200 Float_t lT0Current=0.;
203 // Int_t nrun = esd->GetRunNumber();
207 Float_t sigmaFill = fT0width;
209 if(sigmaFill < 20) sigmaFill = 140;
211 fCalculated[0]=-1000*t0tof[0]; // best t0
212 fCalculated[1]=1000*t0tof[1]; // sigma best t0
213 fCalculated[2] = t0fill; //t0 fill
214 fCalculated[3] = t0tof[2]; // n TOF tracks
215 fCalculated[4]=-1000*t0tof[0]; // TOF t0
216 fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
217 fCalculated[6]=sigmaFill; // sigma t0 fill
218 fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
220 if(fCalculated[7] > 30) thrGood = 10000000;
223 fCalculated[8] = t0tof[4]; // real time in s
224 fCalculated[9] = t0tof[5]; // cpu time in s
226 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
227 fT0sigma=fCalculated[1];
228 lT0Current=fCalculated[0];
231 fCalculated[4] = t0fill;
232 fCalculated[5] = sigmaFill;
235 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
237 fCalculated[4] = t0fill;
238 fCalculated[5] = sigmaFill;
243 Double_t w1 = 1./t0sigma/t0sigma;
244 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
246 Double_t wtot = w1+w2;
248 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
249 fT0sigma = TMath::Sqrt(1./wtot);
257 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
258 fCalculated[1]=fT0sigma;
259 fCalculated[0]=lT0Current;
262 if(fT0sigma >= 1000 || fNoTOFT0){
264 fT0sigma = sigmaFill;
266 fCalculated[0] = t0fill;
267 fCalculated[1] = sigmaFill;
271 Float_t *t0values = new Float_t[fNmomBins];
272 Float_t *t0resolution = new Float_t[fNmomBins];
273 if(fCalculated[7] < 100){
274 for(Int_t i=0;i<fNmomBins;i++){
275 t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
276 t0tof[0] = t0maker->GetResult(0);
277 t0tof[1] = t0maker->GetResult(1);
278 t0tof[2] = t0maker->GetResult(2);
279 t0tof[3] = t0maker->GetResult(3);
281 Float_t t0bin =-1000*t0tof[0]; // best t0
282 Float_t t0binRes =1000*t0tof[1]; // sigma best t0
284 if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
287 Double_t w1 = 1./t0sigma/t0sigma;
288 Double_t w2 = 1./t0binRes/t0binRes;
290 Double_t wtot = w1+w2;
292 t0bin = (w1*t0time + w2*t0bin) / wtot;
293 t0binRes = TMath::Sqrt(1./wtot);
298 t0binRes = sigmaFill;
305 t0resolution[i] = t0binRes;
309 for(Int_t i=0;i<fNmomBins;i++){
310 t0values[i] = lT0Current;
311 t0resolution[i] = fT0sigma;
314 for(Int_t i=0;i<fNmomBins;i++){
315 fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
316 fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
320 delete[] t0resolution;
324 //____________________________________________________________________________
325 Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
326 Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
328 fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
329 fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
333 //____________________________________________________________________________
334 void AliTOFT0maker::SetTOFResponse(){
335 fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
337 //____________________________________________________________________________
338 Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
339 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
343 //____________________________________________________________________________
344 void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
346 // Recalculate TOF PID probabilities
349 // subtruct t0 for each track
350 Int_t ntracks = esd->GetNumberOfTracks();
353 AliESDtrack *t=esd->GetTrack(ntracks);
355 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
357 Double_t time=t->GetTOFsignal();
358 Float_t p = t->GetP();
360 Double_t *t0=GetT0p(p);
362 t->SetTOFsignal(time);
365 for(Int_t i=0;i<fNmomBins;i++){
366 fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
371 //____________________________________________________________________________
372 void AliTOFT0maker::LoadChannelMap(const char *filename){
373 // Load the histo with the channel off map
374 TFile *f= new TFile(filename);
376 printf("Cannot open the channel map file (%s)\n",filename);
380 fHmapChannel = (TH1F *) f->Get("hChEnabled");
383 printf("Cannot laod the channel map histo (from %s)\n",filename);
388 //____________________________________________________________________________
389 void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
390 // Switch off the disable channel
392 printf("Channel Map is not available\n");
396 Int_t ntracks = esd->GetNumberOfTracks();
399 AliESDtrack *t=esd->GetTrack(ntracks);
401 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
403 Int_t chan = t->GetTOFCalChannel();
405 if(fHmapChannel->GetBinContent(chan) < 0.01){
406 t->ResetStatus(AliESDtrack::kTOFout);
412 AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
417 Float_t TOFtimeResolutionDefault=80;
419 Float_t t0 = gRandom->Gaus(0.,fT0width);
421 Float_t extraSmearing = 0;
423 if(fTimeResolution > TOFtimeResolutionDefault){
424 extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
427 // subtruct t0 for each track
428 Int_t ntracks = esd->GetNumberOfTracks();
431 AliESDtrack *t=esd->GetTrack(ntracks);
433 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
435 /* check if channel is enabled */
437 if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
438 /* reset TOF status */
439 t->ResetStatus(AliESDtrack::kTOFin);
440 t->ResetStatus(AliESDtrack::kTOFout);
441 t->ResetStatus(AliESDtrack::kTOFmismatch);
442 t->ResetStatus(AliESDtrack::kTOFpid);
446 Double_t time=t->GetTOFsignal();
451 Float_t smearing = gRandom->Gaus(0.,extraSmearing);
455 t->SetTOFsignal(time);
460 //_________________________________________________________________________
461 void AliTOFT0maker::WriteInESD(AliESDEvent *esd){
463 // Write t0Gen, t0ResGen, nt0;
464 // t0resESD[0:nt0], it0ESD[0:nt0]
465 // in the AliESDEvent
467 Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();
470 Float_t *t0 = new Float_t[nMomBins];
471 Float_t *t0res = new Float_t[nMomBins];
472 Int_t *multT0 = new Int_t[nMomBins];
474 for(Int_t i=0;i<nMomBins;i++){
475 // printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
478 for(Int_t j=0;j < nt0;j++){
479 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
486 t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
487 t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
492 Int_t iMultT0=0,nmult=0;
493 for(Int_t j=0;j < nt0;j++){ // find the most frequent
494 if(multT0[j] > nmult){
500 Float_t *t0ESD = new Float_t[nMomBins];
501 Float_t *t0resESD = new Float_t[nMomBins];
502 Int_t *it0ESD = new Int_t[nMomBins];
504 Float_t t0Gen,t0ResGen;
506 t0ResGen = t0res[iMultT0];
508 // printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
509 for(Int_t i=0;i<nMomBins;i++){
510 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
511 t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
512 t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
514 // printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
519 // Write t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent
521 AliTOFHeader *tofHeader =
522 new AliTOFHeader(t0Gen,t0ResGen,nt0,
523 t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
525 esd->SetTOFHeader(tofHeader);
529 AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
530 AliDebug(1,Form("%d ",nt0));
531 for (Int_t ii=0; ii<nt0; ii++)
532 AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));