free arrays before of return in PropagateBack
[u/mrichter/AliRoot.git] / TOF / AliTOFT0v1.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
16 /* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
17
18 //_________________________________________________________________________
19 // This is a TTask that made the calculation of the Time zero using TOF.
20 // Description: The algorithm used to calculate the time zero of interaction
21 // using TOF detector is the following.
22 // We select in the ESD some "primary" particles - or tracks in the following - 
23 // that strike the TOF detector (the larger part are pions, kaons or protons). 
24 // We choose a set of 10 selected tracks, for each track You have the length
25 // of the track when the TOF is reached, 
26 // the momentum and the time of flight
27 // given by the TOF detector.
28 // Let consider now only one set of 10 tracks (the algorithm is the same for all sets).
29 // Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton,
30 // we consider all the 3 at 10 possible cases. 
31 // For each track in each (mass) configuration
32 // (a configuration can be e.g. pion/pion/kaon/proton/pion/proton/kaon/kaon/pion/pion)
33 // we calculate the time zero (we know in fact the velocity of the track after 
34 // the assumption about its mass, the time of flight given by the TOF, and the 
35 // corresponding path travelled till the TOF detector). Then for each mass configuration we have
36 // 10 time zero and we can calculate the ChiSquare for the current configuration using the 
37 // weighted mean over all 10 time zero.
38 // We call the best assignment the mass configuration that gives the minimum value of the ChiSquare. 
39 // We plot the weighted mean over all 10 time zero for the best assignment, 
40 // the ChiSquare for the best assignment and the corresponding confidence level.
41 // The strong assumption is the MC selection of primary particles. It will be introduced
42 // in the future also some more realistic simulation about this point. 
43 // Use case:
44 // root [0] AliTOFT0v1 * tzero = new AliTOFT0v1("galice.root")
45 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
46 // root [1] tzero->ExecuteTask()
47 // root [2] tzero->ExecuteTask("tim")
48 //             // available parameters:
49 //             tim - print benchmarking information
50 //             all - print usefull informations about the number of misidentified tracks 
51 //                   and a comparison about the true configuration (known from MC) and the best
52 //                   assignment
53 // Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions 
54 //-- Author: F. Pierella
55 //-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
56 //-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table
57 //   for Power3)
58 //////////////////////////////////////////////////////////////////////////////
59
60 #include "AliESDtrack.h"
61 #include "AliESDEvent.h"
62 #include "AliTOFT0v1.h"
63 #include "TBenchmark.h"
64 #include "AliPID.h"
65 #include "AliESDpid.h"
66
67 ClassImp(AliTOFT0v1)
68            
69 //____________________________________________________________________________ 
70 AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
71   TObject(),
72   fLowerMomBound(0.5),
73   fUpperMomBound(3),  
74   fTimeCorr(0.), 
75   fEvent(0x0),
76   fPIDesd(extPID),
77   fTracks(new TObjArray(10)),
78   fGTracks(new TObjArray(10)),
79   fTracksT0(new TObjArray(10)),
80   fOptFlag(kFALSE)
81 {
82   //
83   // default constructor
84   //
85   if(AliPID::ParticleMass(0) == 0) new AliPID();
86
87   if(!fPIDesd){
88     fPIDesd = new AliESDpid();
89   }
90
91   Init(NULL);
92
93   //initialise lookup table for power 3
94   // a set should only have 10 tracks a t maximum
95   // so up to 15 should be more than enough
96   for (Int_t i=0; i<15; ++i) {
97     fLookupPowerThree[i]=ToCalculatePower(3,i);
98   }
99 }
100
101 //____________________________________________________________________________ 
102 AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID): 
103   TObject(),
104   fLowerMomBound(0.5),
105   fUpperMomBound(3.0),  
106   fTimeCorr(0.), 
107   fEvent(event),
108   fPIDesd(extPID),
109   fTracks(new TObjArray(10)),
110   fGTracks(new TObjArray(10)),
111   fTracksT0(new TObjArray(10)),
112   fOptFlag(kFALSE)
113 {
114   //
115   // real constructor
116   //
117   if(AliPID::ParticleMass(0) == 0) new AliPID();
118
119   if(!fPIDesd){
120     fPIDesd = new AliESDpid();
121   }
122
123   Init(event);
124   //initialise lookup table for power 3
125   for (Int_t i=0; i<15; ++i) {
126     fLookupPowerThree[i]=Int_t(TMath::Power(3,i));
127   }
128 }
129 //____________________________________________________________________________ 
130 AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
131 {
132  //
133   // assign. operator
134   //
135
136   if (this == &tzero)
137     return *this;
138   
139   fLowerMomBound=tzero.fLowerMomBound;
140   fUpperMomBound=tzero.fUpperMomBound;  
141   fTimeCorr=tzero.fTimeCorr; 
142   fEvent=tzero.fEvent;
143   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
144   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
145   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
146   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
147
148   fTracks=tzero.fTracks;
149   fGTracks=tzero.fGTracks;
150   fTracksT0=tzero.fTracksT0;
151
152   for (Int_t ii=0; ii<tzero.fTracks->GetEntries(); ii++)
153     fTracks->AddLast(tzero.fTracks->At(ii));
154
155   for (Int_t ii=0; ii<tzero.fGTracks->GetEntries(); ii++)
156     fGTracks->AddLast(tzero.fGTracks->At(ii));
157
158   for (Int_t ii=0; ii<tzero.fTracksT0->GetEntries(); ii++)
159     fTracksT0->AddLast(tzero.fTracksT0->At(ii));
160   
161   fOptFlag=tzero.fOptFlag;
162
163   return *this;
164 }
165
166 //____________________________________________________________________________ 
167 AliTOFT0v1::~AliTOFT0v1()
168 {
169   // dtor
170   fEvent=NULL;
171   
172   if (fTracks) {
173     fTracks->Clear();
174     delete fTracks;
175     fTracks=0x0;
176   }
177
178   if (fGTracks) {
179     fGTracks->Clear();
180     delete fGTracks;
181     fGTracks=0x0;
182   }
183
184   if (fTracksT0) {
185     fTracksT0->Clear();
186     delete fTracksT0;
187     fTracksT0=0x0;
188   }
189
190
191 }
192 //____________________________________________________________________________ 
193
194 void
195 AliTOFT0v1::Init(AliESDEvent *event) 
196 {
197
198   /* 
199    * init
200    */
201
202   fEvent = event;
203   fT0SigmaT0def[0]=0.;
204   fT0SigmaT0def[1]=0.6;
205   fT0SigmaT0def[2]=0.;
206   fT0SigmaT0def[3]=0.;
207
208 }
209 //____________________________________________________________________________ 
210 Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut) 
211
212   TBenchmark *bench=new TBenchmark();
213   bench->Start("t0computation");
214
215   // Caluclate the Event Time using the ESD TOF time
216
217   fT0SigmaT0def[0]=0.;
218   fT0SigmaT0def[1]=0.600;
219   fT0SigmaT0def[2]=0.;
220   fT0SigmaT0def[3]=0.;
221   
222   const Int_t nmaxtracksinsetMax=10;
223   Int_t nmaxtracksinset=10;
224 //   if(strstr(option,"all")){
225 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
226 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
227 //   }
228   
229   
230   Int_t nsets=0;
231   Int_t nUsedTracks=0;
232   Int_t ngoodsetsSel= 0;
233   Float_t t0bestSel[300];
234   Float_t eT0bestSel[300];
235   Float_t chiSquarebestSel[300];
236   Float_t confLevelbestSel[300];
237   Float_t t0bestallSel=0.;
238   Float_t eT0bestallSel=0.;
239   Float_t sumWt0bestallSel=0.;
240   Float_t eMeanTzeroPi=0.;
241   Float_t meantzeropi=0.;
242   Float_t sumAllweightspi=0.;
243   Double_t t0def=-999;
244   Double_t deltat0def=999;
245   Int_t ngoodtrktrulyused=0;
246   Int_t ntracksinsetmyCut = 0;
247
248   Int_t ntrk=fEvent->GetNumberOfTracks();
249   
250   Int_t ngoodtrk=0;
251   Int_t ngoodtrkt0 =0;
252   Float_t meantime =0;
253   
254   // First Track loop, Selection of good tracks
255
256   fTracks->Clear();
257   for (Int_t itrk=0; itrk<ntrk; itrk++) {
258     AliESDtrack *t=fEvent->GetTrack(itrk);
259     Double_t momOld=t->GetP();
260     Double_t mom=momOld-0.0036*momOld;
261     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
262     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
263     Double_t time=t->GetTOFsignal();
264     
265     time*=1.E-3; // tof given in nanoseconds       
266     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
267
268     if (!AcceptTrack(t)) continue;
269
270     if(t->GetIntegratedLength() < 350)continue; //skip decays
271     if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
272
273     meantime+=time;
274     fTracks->AddLast(t);
275     ngoodtrk++;
276   }
277
278   if(ngoodtrk > 1) meantime /= ngoodtrk;
279
280   if(ngoodtrk>22) nmaxtracksinset = 6;
281
282   fGTracks->Clear();
283   for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
284     AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
285     //    Double_t time=t->GetTOFsignal();
286     //    if((time-meantime*1.E3)<50.E3){ // For pp and per 
287     fGTracks->AddLast(t);
288     ngoodtrkt0++;
289       //    }
290   }
291
292   fTracks->Clear();
293
294   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
295   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
296   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
297
298
299   if(ngoodtrkt0<2){
300     t0def=-999;
301     deltat0def=0.600;
302     fT0SigmaT0def[0]=t0def;
303     fT0SigmaT0def[1]=deltat0def;
304     fT0SigmaT0def[2]=ngoodtrkt0;
305     fT0SigmaT0def[3]=ngoodtrkt0;
306     //goto finish;
307   }
308   if(ngoodtrkt0>=2){
309   // Decide how many tracks in set 
310     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
311     Int_t nset=1;
312
313     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} 
314         
315     // Loop over selected sets
316     
317     if(nset>=1){
318       for (Int_t i=0; i< nset; i++) {   
319         //      printf("Set %i of %i\n",i+1,nset);
320         Float_t t0best=999.;
321         Float_t eT0best=999.;
322         Float_t chisquarebest=99999.;
323         Int_t npionbest=0;
324         
325         fTracksT0->Clear();
326         Int_t ntracksinsetmy=0;      
327         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
328           Int_t index = itrk+i*ntracksinset;
329           if(index < fGTracks->GetEntries()){
330             AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
331             fTracksT0->AddLast(t);
332             ntracksinsetmy++;
333           }
334         }
335
336         // Analyse it
337         
338         Int_t   assparticle[nmaxtracksinsetMax];
339         Float_t exptof[nmaxtracksinsetMax][3];
340         Float_t momErr[nmaxtracksinsetMax][3];
341         Float_t timeofflight[nmaxtracksinsetMax];
342         Float_t momentum[nmaxtracksinsetMax];
343         Float_t timezero[nmaxtracksinsetMax];
344         Float_t weightedtimezero[nmaxtracksinsetMax];
345         Float_t beta[nmaxtracksinsetMax];
346         Float_t texp[nmaxtracksinsetMax];
347         Float_t dtexp[nmaxtracksinsetMax];
348         Float_t sqMomError[nmaxtracksinsetMax];
349         Float_t sqTrackError[nmaxtracksinsetMax];
350         Float_t massarray[3]={0.13957,0.493677,0.9382723};
351         Float_t tracktoflen[nmaxtracksinsetMax];
352         Float_t besttimezero[nmaxtracksinsetMax];
353         Float_t besttexp[nmaxtracksinsetMax];
354         Float_t besttimeofflight[nmaxtracksinsetMax];
355         Float_t bestmomentum[nmaxtracksinsetMax];
356         Float_t bestchisquare[nmaxtracksinsetMax];
357         Float_t bestweightedtimezero[nmaxtracksinsetMax];
358         Float_t bestsqTrackError[nmaxtracksinsetMax];
359         Int_t imass[nmaxtracksinsetMax];
360         
361         for (Int_t j=0; j<ntracksinset; j++) {
362           assparticle[j] = 3;
363           timeofflight[j] = 0;
364           momentum[j] = 0;
365           timezero[j] = 0;
366           weightedtimezero[j] = 0;
367           beta[j] = 0;
368           texp[j] = 0;
369           dtexp[j] = 0;
370           sqMomError[j] = 0;
371           sqTrackError[j] = 0;
372           tracktoflen[j] = 0;
373           besttimezero[j] = 0;
374           besttexp[j] = 0;
375           besttimeofflight[j] = 0;
376           bestmomentum[j] = 0;
377           bestchisquare[j] = 0;
378           bestweightedtimezero[j] = 0;
379           bestsqTrackError[j] = 0;
380           imass[j] = 1;
381         }
382         
383         for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
384           AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
385           Double_t momOld=t->GetP();
386           Double_t mom=momOld-0.0036*momOld;
387           Double_t time=t->GetTOFsignal();
388           
389           time*=1.E-3; // tof given in nanoseconds         
390           Double_t exptime[AliPID::kSPECIESC]; 
391           t->GetIntegratedTimes(exptime,AliPID::kSPECIESC);
392           Double_t toflen=t->GetIntegratedLength();
393           toflen=toflen/100.; // toflen given in m 
394           
395           timeofflight[j]=time;
396           tracktoflen[j]=toflen;
397           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
398           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
399           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
400           momentum[j]=mom;
401           assparticle[j]=3;
402
403           // in principle GetMomError only depends on two indices k=imass[j] and j itslef (see blow in the ncombinatorial loop)
404           // so it should be possible to make a lookup in order to speed up the code:
405           if (fOptFlag) {
406             momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]);
407             momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]);
408             momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]);
409           }
410   } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
411         
412         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
413           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
414                                        +momentum[itz]*momentum[itz]);
415           sqMomError[itz]= ((1.-beta[itz]*beta[itz])*0.01)*((1.-beta[itz]*beta[itz])*0.01)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz])); 
416           sqTrackError[itz]=sqMomError[itz]; //in ns
417           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
418           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
419           sumAllweightspi+=1./sqTrackError[itz];
420           meantzeropi+=weightedtimezero[itz];   
421         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
422         
423         
424         // Then, Combinatorial Algorithm
425         
426         if(ntracksinsetmy<2 )break;
427         
428         for (Int_t j=0; j<ntracksinsetmy; j++) {
429           imass[j] = 3;
430         }
431
432         Int_t ncombinatorial;
433         if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
434         else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
435
436
437         // Loop on mass hypotheses
438         for (Int_t k=0; k < ncombinatorial;k++) {
439           for (Int_t j=0; j<ntracksinsetmy; j++) {
440             imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
441             texp[j]=exptof[j][imass[j]];
442             if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
443             else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
444           }
445
446           Float_t sumAllweights=0.;
447           Float_t meantzero=0.;
448           Float_t eMeanTzero=0.;
449           
450           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
451             sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
452             
453             timezero[itz]=texp[itz]-timeofflight[itz];// in ns                    
454             
455             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
456             sumAllweights+=1./sqTrackError[itz];
457             meantzero+=weightedtimezero[itz];
458             
459           } // end loop for (Int_t itz=0; itz<15;itz++)
460           
461           meantzero=meantzero/sumAllweights; // it is given in [ns]
462           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
463           
464           // calculate chisquare
465           Float_t chisquare=0.;         
466           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
467             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
468           } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
469           
470           if(chisquare<=chisquarebest){
471             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
472               
473               bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
474               besttimezero[iqsq]=timezero[iqsq]; 
475               bestmomentum[iqsq]=momentum[iqsq]; 
476               besttimeofflight[iqsq]=timeofflight[iqsq]; 
477               besttexp[iqsq]=texp[iqsq]; 
478               bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
479               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
480             }
481             
482             Int_t npion=0;
483             for (Int_t j=0; j<ntracksinsetmy; j++) {
484               assparticle[j]=imass[j];
485               if(imass[j] == 0) npion++;
486             }
487             npionbest=npion;
488             chisquarebest=chisquare;          
489             t0best=meantzero;
490             eT0best=eMeanTzero;
491           } // close if(dummychisquare<=chisquare)
492         }
493         
494         Double_t chi2cut[nmaxtracksinsetMax];
495         chi2cut[0] = 0;
496         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
497         for (Int_t j=2; j<ntracksinset; j++) {
498           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
499         }
500         
501         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
502         
503         //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
504         
505         Bool_t kRedoT0 = kFALSE;
506         ntracksinsetmyCut = ntracksinsetmy;
507         Bool_t usetrack[nmaxtracksinsetMax];
508         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
509           usetrack[icsq] = kTRUE;
510           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
511               kRedoT0 = kTRUE;
512               ntracksinsetmyCut--;
513               usetrack[icsq] = kFALSE;
514               //              printf("tracks chi2 = %f\n",bestchisquare[icsq]);
515           }
516         } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
517         
518         // Loop on mass hypotheses Redo
519         if(kRedoT0 && ntracksinsetmyCut > 1){
520           //      printf("Redo T0\n");
521           for (Int_t k=0; k < ncombinatorial;k++) {
522             for (Int_t j=0; j<ntracksinsetmy; j++) {
523               imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
524               texp[j]=exptof[j][imass[j]];
525               if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
526               else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
527             }
528             
529             Float_t sumAllweights=0.;
530             Float_t meantzero=0.;
531             Float_t eMeanTzero=0.;
532             
533             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
534               if(! usetrack[itz]) continue;
535               sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
536               
537               timezero[itz]=texp[itz]-timeofflight[itz];// in ns                          
538               
539               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
540               sumAllweights+=1./sqTrackError[itz];
541               meantzero+=weightedtimezero[itz];
542               
543             } // end loop for (Int_t itz=0; itz<15;itz++)
544             
545             meantzero=meantzero/sumAllweights; // it is given in [ns]
546             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
547             
548             // calculate chisquare
549             
550             Float_t chisquare=0.;               
551             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
552               if(! usetrack[icsq]) continue;
553               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
554             } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
555             
556             Int_t npion=0;
557             for (Int_t j=0; j<ntracksinsetmy; j++) {
558               assparticle[j]=imass[j];
559               if(imass[j] == 0) npion++;
560             }
561             
562             if(chisquare<=chisquarebest && npion>0){
563               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
564                 if(! usetrack[iqsq]) continue;
565                 bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
566                 besttimezero[iqsq]=timezero[iqsq]; 
567                 bestmomentum[iqsq]=momentum[iqsq]; 
568                 besttimeofflight[iqsq]=timeofflight[iqsq]; 
569                 besttexp[iqsq]=texp[iqsq]; 
570                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
571                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
572               }
573               
574               npionbest=npion;
575               chisquarebest=chisquare;        
576               t0best=meantzero;
577               eT0best=eMeanTzero;
578             } // close if(dummychisquare<=chisquare)
579             
580           }
581         }
582                 
583         // filling histos
584         Float_t confLevel=999;
585         
586         // Sets with decent chisquares
587         //      printf("Chi2best of the set = %f \n",chisquarebest);
588         
589         if(chisquarebest<999.){
590           Double_t dblechisquare=(Double_t)chisquarebest;
591           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); 
592
593           Int_t ntrackincurrentsel=0;
594           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
595
596             if(! usetrack[icsq]) continue;
597             
598             ntrackincurrentsel++;
599           }
600           
601           //      printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
602
603           // Pick up only those with C.L. >1%
604           if(confLevel>0.01 && ngoodsetsSel<200){
605             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
606             confLevelbestSel[ngoodsetsSel]=confLevel;
607             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
608             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
609             t0bestallSel += t0best/eT0best/eT0best;
610             sumWt0bestallSel += 1./eT0best/eT0best;
611             ngoodsetsSel++;
612             ngoodtrktrulyused+=ntracksinsetmyCut;
613             //      printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
614           }
615           else{
616             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
617           }
618         }       
619         fTracksT0->Clear();
620         nsets++;
621         
622       } // end for the current set
623       
624       //Redo the computation of the best in order to esclude very bad samples
625         if(ngoodsetsSel > 1){
626             Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
627             Int_t nsamples=ngoodsetsSel;
628             ngoodsetsSel=0;
629             t0bestallSel=0;
630             sumWt0bestallSel=0;
631             for (Int_t itz=0; itz<nsamples;itz++) {
632                 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
633                     t0bestallSel += t0bestSel[itz];
634                     sumWt0bestallSel += eT0bestSel[itz];              
635                     ngoodsetsSel++;   
636                     //        printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
637                 }
638                 else{
639                   //          printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
640                 }
641             }
642         }
643         if(ngoodsetsSel < 1){
644             sumWt0bestallSel = 0.0;
645         }
646       //--------------------------------End recomputation
647
648       nUsedTracks =  ngoodtrkt0;  
649       if(strstr(option,"all")){
650         if(sumAllweightspi>0.){
651           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
652           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
653         }      
654         
655         //      printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
656
657         if(sumWt0bestallSel>0){
658           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
659           eT0bestallSel = sqrt(1./sumWt0bestallSel);
660           //      printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);          
661         }// end of if(sumWt0bestallSel>0){
662         
663       }
664       
665       t0def=t0bestallSel;
666       deltat0def=eT0bestallSel;
667       
668       fT0SigmaT0def[0]=t0def;
669       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
670       fT0SigmaT0def[2]=ngoodtrkt0;
671       fT0SigmaT0def[3]=ngoodtrktrulyused;
672     }
673   }
674
675   fGTracks->Clear();
676
677   if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
678
679   bench->Stop("t0computation");
680
681   fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
682   fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
683
684 //   bench->Print("t0computation");
685 //   printf("(%4.1f < p < %4.1f GeV/c) T0-TOF =%9.1f +/- %5.1f ps (n_track = %i)\n\n",pMinCut,pMaxCut,-fT0SigmaT0def[0]*1000,fT0SigmaT0def[1]*1000,Int_t(fT0SigmaT0def[3]));
686
687   delete bench;
688   bench=NULL;
689
690   return fT0SigmaT0def;
691   }
692 //__________________________________________________________________
693 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
694 {
695   // Take the error extimate for the TOF time in the track reconstruction
696
697   static const Double_t kMasses[]={
698     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
699   };
700
701   Double_t mass=kMasses[index+2];
702
703   Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
704
705   return sigma;
706 }
707
708 //__________________________________________________________________
709 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
710 {
711
712   /* TPC refit */
713   if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
714   /* do not accept kink daughters */
715   if (track->GetKinkIndex(0)>0) return kFALSE;
716   /* N clusters TPC */
717   if (track->GetTPCclusters(0) < 50) return kFALSE;
718   /* chi2 TPC */
719   if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
720   /* sigma to vertex */
721   if (GetSigmaToVertex(track) > 4.) return kFALSE;
722   
723   /* accept track */
724   return kTRUE;
725
726 }
727
728 //____________________________________________________________________
729 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
730 {
731   // Calculates the number of sigma to the vertex.
732
733   Float_t b[2];
734   Float_t bRes[2];
735   Float_t bCov[3];
736   esdTrack->GetImpactParameters(b,bCov);
737   
738   if (bCov[0]<=0 || bCov[2]<=0) {
739     bCov[0]=0; bCov[2]=0;
740   }
741   bRes[0] = TMath::Sqrt(bCov[0]);
742   bRes[1] = TMath::Sqrt(bCov[2]);
743
744   // -----------------------------------
745   // How to get to a n-sigma cut?
746   //
747   // The accumulated statistics from 0 to d is
748   //
749   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
750   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
751   //
752   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
753   // Can this be expressed in a different way?
754
755   if (bRes[0] == 0 || bRes[1] ==0)
756     return -1;
757
758   //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
759   Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
760
761   // work around precision problem
762   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
763   // 1e-15 corresponds to nsigma ~ 7.7
764   if (TMath::Exp(-d * d / 2) < 1e-15)
765     return 1000;
766
767   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
768   return nSigma;
769 }
770 //____________________________________________________________________
771
772 Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
773     Bool_t status = kFALSE;
774     
775     Double_t exptimes[AliPID::kSPECIESC];
776     track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC);
777     
778     Float_t dedx = track->GetTPCsignal();
779     
780     Double_t ptpc[3];
781     track->GetInnerPxPyPz(ptpc);
782     Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
783
784     if(imass > 2 || imass < 0) return status;
785     Int_t i = imass+2;
786     
787     AliPID::EParticleType type=AliPID::EParticleType(i);
788     
789     Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
790     Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
791         
792     if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
793         status = kTRUE;
794     }
795     
796     return status;
797 }
798
799 //____________________________________________________________________
800 Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const
801 {
802   //
803   // Returns base^exponent
804   //
805
806   Float_t power=1.;
807
808   for (Int_t ii=exponent; ii>0; ii--)
809       power=power*base;
810
811   return power;
812
813 }
814 //____________________________________________________________________
815 Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const
816 {
817   //
818   // Returns base^exponent
819   //
820
821   Int_t power=1;
822
823   for (Int_t ii=exponent; ii>0; ii--)
824       power=power*base;
825
826   return power;
827
828 }