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