warning fix
[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.cxx
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 /* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
16
17 /////////////////////////////////////////////////////////////////////////////
18 //                                                                         //
19 //  This class contains the basic functions for the time zero              //
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)  //
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                              //
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                //
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             // 
40 //                                                                         //
41 // Let consider that:                                                      //
42 // - the PIF is automatically recalculated with the event time subtrction  //
43 //                                                                         //
44 /////////////////////////////////////////////////////////////////////////////
45
46 #include "AliTOFT0maker.h"
47 #include "AliPID.h"
48 #include "AliLog.h"
49 #include "AliESDpid.h"
50 #include "AliESDEvent.h"
51 #include "TFile.h"
52 #include "TH1F.h"
53 #include "AliTOFcalib.h"
54 #include "AliTOFRunParams.h"
55 #include "TRandom.h"
56 #include "AliTOFHeader.h"
57
58 ClassImp(AliTOFT0maker)
59            
60 //____________________________________________________________________________ 
61 AliTOFT0maker::AliTOFT0maker():
62   TObject(),
63   fT0TOF(NULL),
64   fPIDesd(NULL),
65   fExternalPIDFlag(kFALSE),
66   fTOFcalib(NULL),
67   fNoTOFT0(0),
68   fNmomBins(0),
69   fTimeResolution(100),
70   fT0sigma(1000),
71   fHmapChannel(0),
72   fKmask(0),
73   fT0width(150.),
74   fT0spreadExt(-1.),
75   fT0fillExt(0)
76 {
77   // ctr
78   fCalculated[0] = 0;
79   fCalculated[1] = 0;
80   fCalculated[2] = 0;
81   fCalculated[3] = 0;
82
83   fT0cur[0]=0.;
84   fT0cur[1]=0.;
85
86   if(AliPID::ParticleMass(0) == 0) new AliPID();
87
88   fPIDesd = new AliESDpid();
89
90   fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
91   SetTOFResponse();
92
93   fT0TOF = new AliTOFT0v1(fPIDesd);
94
95 }
96 //____________________________________________________________________________ 
97 AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
98     TObject(),
99     fT0TOF(NULL),
100     fPIDesd(externalPID),
101     fExternalPIDFlag(kTRUE),
102     fTOFcalib(tofCalib),
103     fNoTOFT0(0),
104     fNmomBins(0),
105     fTimeResolution(100),
106     fT0sigma(1000),
107     fHmapChannel(0),
108     fKmask(0),
109     fT0width(150.),
110     fT0spreadExt(-1.),
111     fT0fillExt(0)
112 {
113   // ctr
114   fCalculated[0] = 0;
115   fCalculated[1] = 0;
116   fCalculated[2] = 0;
117   fCalculated[3] = 0;
118
119   fT0cur[0]=0.;
120   fT0cur[1]=0.;
121
122   if(AliPID::ParticleMass(0) == 0) new AliPID();
123
124   if(!fPIDesd){
125     fPIDesd = new AliESDpid();
126     fExternalPIDFlag = kFALSE;
127   }
128
129   fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
130   SetTOFResponse();
131
132   fT0TOF = new AliTOFT0v1(fPIDesd);
133
134 }
135
136 //____________________________________________________________________________ 
137 AliTOFT0maker::~AliTOFT0maker()
138 {
139   // dtor
140
141   delete fT0TOF;
142   if (!fExternalPIDFlag) delete fPIDesd;
143 }
144 //____________________________________________________________________________ 
145 Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
146   //
147   // Remake TOF PID probabilities
148   //
149   Double_t t0tof[6];
150
151   if(fKmask) ApplyMask(esd);
152
153   Double_t t0fill = 0.;
154
155   fPIDesd->GetTOFResponse().ResetT0info();
156
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());
162       if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
163       else{
164         SetT0FillWidth(t0spread);
165         t0fill = fT0fillExt;
166       }
167     }
168     else{
169       if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
170       t0fill = fT0fillExt;
171     }
172   }
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   }
183
184   Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
185
186   fT0TOF->Init(esd);
187   AliTOFT0v1* t0maker= fT0TOF;
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);
194   t0tof[4] = t0maker->GetResult(4);
195   t0tof[5] = t0maker->GetResult(5);
196
197   Float_t lT0Current=0.;
198   fT0sigma=1000;
199
200 //   Int_t nrun = esd->GetRunNumber();
201
202   t0time += t0fill;
203
204   Float_t sigmaFill = fT0width;
205
206   if(sigmaFill < 20) sigmaFill = 140;
207
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
217   if(fCalculated[7] > 30) thrGood = 10000000;
218
219   //statistics
220   fCalculated[8] = t0tof[4]; // real time in s
221   fCalculated[9] = t0tof[5]; // cpu time in s
222
223   if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
224     fT0sigma=fCalculated[1];
225     lT0Current=fCalculated[0];
226   }
227   else{
228     fCalculated[4] = t0fill;
229     fCalculated[5] = sigmaFill;
230   }
231
232   if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
233     fT0sigma =1000;
234     fCalculated[4] = t0fill;
235     fCalculated[5] = sigmaFill;
236   }
237
238   if(t0sigma < 1000){
239     if(fT0sigma < 1000){
240       Double_t w1 = 1./t0sigma/t0sigma;
241       Double_t w2 = 1./fCalculated[1]/fCalculated[1];
242
243       Double_t wtot = w1+w2;
244
245       lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
246       fT0sigma = TMath::Sqrt(1./wtot);
247     }
248     else{
249       lT0Current=t0time;
250       fT0sigma=t0sigma;
251     }
252   }
253
254   if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
255     fCalculated[1]=fT0sigma;
256     fCalculated[0]=lT0Current;
257   }
258
259   if(fT0sigma >= 1000 || fNoTOFT0){
260     lT0Current = t0fill;
261     fT0sigma = sigmaFill;
262
263     fCalculated[0] = t0fill;
264     fCalculated[1] = sigmaFill;
265   }
266
267   // T0 pt bin
268   Float_t *t0values = new Float_t[fNmomBins];
269   Float_t *t0resolution = new Float_t[fNmomBins];
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);
277
278       Float_t t0bin =-1000*t0tof[0]; // best t0
279       Float_t t0binRes =1000*t0tof[1]; // sigma best t0
280       
281       if(t0binRes < sigmaFill  && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
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       }
301       t0values[i] = t0bin;
302       t0resolution[i] = t0binRes;
303     }
304   }
305   else{
306     for(Int_t i=0;i<fNmomBins;i++){
307       t0values[i] = lT0Current;
308       t0resolution[i] = fT0sigma;
309     }
310   }
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;
318
319   return fCalculated;
320 }
321 //____________________________________________________________________________ 
322 Double_t  *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
323   Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
324   
325   fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
326   fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
327
328   return fT0cur;
329 }
330 //____________________________________________________________________________ 
331 void AliTOFT0maker::SetTOFResponse(){
332   fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
333 }
334 //____________________________________________________________________________ 
335 Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
336   Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
337
338   return sigma;
339 }
340 //____________________________________________________________________________ 
341 void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
342   //
343   // Recalculate TOF PID probabilities
344   //
345
346 // subtruct t0 for each track
347    Int_t ntracks = esd->GetNumberOfTracks();
348   
349    while (ntracks--) {
350        AliESDtrack *t=esd->GetTrack(ntracks);
351     
352        if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
353     
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    
366   //
367 }
368 //____________________________________________________________________________ 
369 void  AliTOFT0maker::LoadChannelMap(const char *filename){
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 //____________________________________________________________________________ 
386 void 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
408 Float_t  
409 AliTOFT0maker::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 */
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         }
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 }
457 //_________________________________________________________________________
458 void     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,
520                          t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
521   
522     esd->SetTOFHeader(tofHeader);
523
524     AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
525     AliDebug(1,Form("%d ",nt0));
526     for (Int_t ii=0; ii<nt0; ii++)
527       AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
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 }