8f95f3c288bd44568c9803d232c3ed4df42d62d0
[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 //////////////////////////////////////////////////////////////////////////////
57
58 #include "AliESDtrack.h"
59 #include "AliESDEvent.h"
60 #include "AliTOFT0v1.h"
61
62 ClassImp(AliTOFT0v1)
63            
64 //____________________________________________________________________________ 
65 AliTOFT0v1::AliTOFT0v1():
66   TObject(),
67   fLowerMomBound(0.5),
68   fUpperMomBound(3),  
69   fTimeResolution(0.80e-10), 
70   fTimeCorr(0.), 
71   fEvent(0x0)
72 //   fCalib(0x0)
73 {
74   //
75   // default constructor
76   //
77
78   Init(NULL);
79     
80 }
81
82            
83 //____________________________________________________________________________ 
84 AliTOFT0v1::AliTOFT0v1(AliESDEvent* event): 
85   TObject(),
86   fLowerMomBound(0.5),
87   fUpperMomBound(3.0),  
88   fTimeResolution(0.80e-10), 
89   fTimeCorr(0.), 
90   fEvent(event)
91 //   fCalib(0x0)
92 {
93   //
94   // real constructor
95   //
96   
97   Init(event);
98
99 }
100
101 /* copy-constructor and operator= suppresed 
102
103 //____________________________________________________________________________ 
104 AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero):
105   TObject(),
106   fLowerMomBound(tzero.fLowerMomBound),
107   fUpperMomBound(tzero.fUpperMomBound),  
108   fTimeResolution(tzero.fTimeResolution), 
109   fTimeCorr(tzero.fTimeCorr), 
110   fEvent(tzero.fEvent)
111 //   fCalib(tzero.fCalib)
112 {
113   //
114   // copy constructor
115   //
116     
117   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
118   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
119   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
120   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
121
122 }
123
124 //____________________________________________________________________________ 
125 AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
126 {
127  //
128   // assign. operator
129   //
130
131   if (this == &tzero)
132     return *this;
133   
134   fLowerMomBound=tzero.fLowerMomBound;
135   fUpperMomBound=tzero.fUpperMomBound;  
136   fTimeResolution=tzero.fTimeResolution; 
137   fTimeCorr=tzero.fTimeCorr; 
138   fEvent=tzero.fEvent;
139 //   fCalib=tzero.fCalib;
140   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
141   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
142   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
143   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
144
145   return *this;
146 }
147
148 */
149 //____________________________________________________________________________ 
150 AliTOFT0v1::~AliTOFT0v1()
151 {
152   // dtor
153 //   fCalib=NULL;
154   fEvent=NULL;
155
156 }
157 //____________________________________________________________________________ 
158
159 void
160 AliTOFT0v1::Init(AliESDEvent *event) 
161 {
162
163   /* 
164    * init
165    */
166
167   fEvent = event;
168   fT0SigmaT0def[0]=0.;
169   fT0SigmaT0def[1]=0.6;
170   fT0SigmaT0def[2]=0.;
171   fT0SigmaT0def[3]=0.;
172
173 }
174
175 //____________________________________________________________________________ 
176 void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
177   // Set the TOF time resolution
178   fTimeResolution=timeresolution;
179 }
180 //____________________________________________________________________________
181 //____________________________________________________________________________
182 Double_t * AliTOFT0v1::DefineT0(Option_t *option) 
183
184   // Caluclate the Event Time using the ESD TOF time
185
186   fT0SigmaT0def[0]=0.;
187   fT0SigmaT0def[1]=0.600;
188   fT0SigmaT0def[2]=0.;
189   fT0SigmaT0def[3]=0.;
190
191  Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
192   
193   const Int_t nmaxtracksinset=10;
194 //   if(strstr(option,"all")){
195 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
196 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
197 //   }
198   
199   Int_t nsets=0;
200   Int_t nUsedTracks=0;
201   Int_t ngoodsetsSel= 0;
202   Float_t t0bestSel[300];
203   Float_t eT0bestSel[300];
204   Float_t chiSquarebestSel[300];
205   Float_t confLevelbestSel[300];
206   Float_t t0bestallSel=0.;
207   Float_t eT0bestallSel=0.;
208   Float_t sumWt0bestallSel=0.;
209   Float_t eMeanTzeroPi=0.;
210   Float_t meantzeropi=0.;
211   Float_t sumAllweightspi=0.;
212   Double_t t0def=-999;
213   Double_t deltat0def=999;
214   Int_t ngoodtrktrulyused=0;
215   Int_t ntracksinsetmyCut = 0;
216
217   Int_t ntrk=fEvent->GetNumberOfTracks();
218   
219   AliESDtrack **tracks=new AliESDtrack*[ntrk];
220   Int_t ngoodtrk=0;
221   Int_t ngoodtrkt0 =0;
222   Float_t mintime =1E6;
223   
224   // First Track loop, Selection of good tracks
225
226   for (Int_t itrk=0; itrk<ntrk; itrk++) {
227     AliESDtrack *t=fEvent->GetTrack(itrk);
228     Double_t momOld=t->GetP();
229     Double_t mom=momOld-0.0036*momOld;
230     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
231     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
232     Double_t time=t->GetTOFsignal();
233     
234     time*=1.E-3; // tof given in nanoseconds       
235     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
236    
237     if (!AcceptTrack(t)) continue;
238
239     if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
240     if(time <= mintime) mintime=time;
241     tracks[ngoodtrk]=t;
242     ngoodtrk++;
243   }
244   
245   
246 //    cout << " N. of ESD tracks                    : " << ntrk << endl;
247 //    cout << " N. of preselected tracks            : " << ngoodtrk << endl;
248 //    cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
249   
250   AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
251   
252   for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
253     AliESDtrack *t=tracks[jtrk];
254     Double_t time=t->GetTOFsignal();
255
256     if((time-mintime*1.E3)<50.E3){ // For pp and per 
257       gtracks[ngoodtrkt0]=t;
258       ngoodtrkt0++;
259     }
260   }
261   
262
263   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
264   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
265   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
266
267   if(ngoodtrkt0<2){
268 //     cout << "less than 2 tracks, skip event " << endl;
269     t0def=-999;
270     deltat0def=0.600;
271     fT0SigmaT0def[0]=t0def;
272     fT0SigmaT0def[1]=deltat0def;
273     fT0SigmaT0def[2]=ngoodtrkt0;
274     fT0SigmaT0def[3]=ngoodtrkt0;
275     //goto finish;
276   }
277   if(ngoodtrkt0>=2){
278   // Decide how many tracks in set 
279     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
280     Int_t nset=1;
281
282     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} 
283         
284     // Loop over selected sets
285     
286     if(nset>=1){
287       for (Int_t i=0; i< nset; i++) {   
288         
289         Float_t t0best=999.;
290         Float_t eT0best=999.;
291         Float_t chisquarebest=99999.;
292         Int_t npionbest=0;
293         
294         Int_t ntracksinsetmy=0;      
295         AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
296         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
297           Int_t index = itrk+i*ntracksinset;
298           if(index < ngoodtrkt0){
299             AliESDtrack *t=gtracks[index];
300             tracksT0[itrk]=t;
301             ntracksinsetmy++;
302           }
303         }
304         
305         // Analyse it
306         
307         Int_t   assparticle[nmaxtracksinset];
308         Float_t exptof[nmaxtracksinset][3];
309         Float_t timeofflight[nmaxtracksinset];
310         Float_t momentum[nmaxtracksinset];
311         Float_t timezero[nmaxtracksinset];
312         Float_t weightedtimezero[nmaxtracksinset];
313         Float_t beta[nmaxtracksinset];
314         Float_t texp[nmaxtracksinset];
315         Float_t dtexp[nmaxtracksinset];
316         Float_t sqMomError[nmaxtracksinset];
317         Float_t sqTrackError[nmaxtracksinset];
318         Float_t massarray[3]={0.13957,0.493677,0.9382723};
319         Float_t tracktoflen[nmaxtracksinset];
320         Float_t besttimezero[nmaxtracksinset];
321         Float_t besttexp[nmaxtracksinset];
322         Float_t besttimeofflight[nmaxtracksinset];
323         Float_t bestmomentum[nmaxtracksinset];
324         Float_t bestchisquare[nmaxtracksinset];
325         Float_t bestweightedtimezero[nmaxtracksinset];
326         Float_t bestsqTrackError[nmaxtracksinset];
327         Int_t imass[nmaxtracksinset];
328         
329         for (Int_t j=0; j<ntracksinset; j++) {
330           assparticle[j] = 3;
331           timeofflight[j] = 0;
332           momentum[j] = 0;
333           timezero[j] = 0;
334           weightedtimezero[j] = 0;
335           beta[j] = 0;
336           texp[j] = 0;
337           dtexp[j] = 0;
338           sqMomError[j] = 0;
339           sqTrackError[j] = 0;
340           tracktoflen[j] = 0;
341           besttimezero[j] = 0;
342           besttexp[j] = 0;
343           besttimeofflight[j] = 0;
344           bestmomentum[j] = 0;
345           bestchisquare[j] = 0;
346           bestweightedtimezero[j] = 0;
347           bestsqTrackError[j] = 0;
348           imass[j] = 1;
349         }
350         
351         for (Int_t j=0; j<ntracksinsetmy; j++) {
352           AliESDtrack *t=tracksT0[j];
353           Double_t momOld=t->GetP();
354           Double_t mom=momOld-0.0036*momOld;
355           Double_t time=t->GetTOFsignal();
356           
357           time*=1.E-3; // tof given in nanoseconds         
358           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
359           Double_t toflen=t->GetIntegratedLength();
360           toflen=toflen/100.; // toflen given in m 
361           
362           timeofflight[j]=time;
363           tracktoflen[j]=toflen;
364           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
365           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
366           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
367           momentum[j]=mom;
368           assparticle[j]=3;
369           
370         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
371         
372         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
373           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
374                                        +momentum[itz]*momentum[itz]);
375           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])); 
376           sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
377           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
378           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
379           sumAllweightspi+=1./sqTrackError[itz];
380           meantzeropi+=weightedtimezero[itz];   
381         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
382         
383         
384         // Then, Combinatorial Algorithm
385         
386         if(ntracksinsetmy<2 )break;
387         
388         for (Int_t j=0; j<ntracksinsetmy; j++) {
389           imass[j] = 3;
390         }
391         
392         Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
393         
394         // Loop on mass hypotheses
395         for (Int_t k=0; k < ncombinatorial;k++) {
396           for (Int_t j=0; j<ntracksinsetmy; j++) {
397             imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
398             texp[j]=exptof[j][imass[j]];
399             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
400           }
401           Float_t sumAllweights=0.;
402           Float_t meantzero=0.;
403           Float_t eMeanTzero=0.;
404           
405           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
406             sqTrackError[itz]=
407               (timeresolutioninns*
408                timeresolutioninns
409                +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
410             
411             timezero[itz]=texp[itz]-timeofflight[itz];// in ns                    
412             
413             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
414             sumAllweights+=1./sqTrackError[itz];
415             meantzero+=weightedtimezero[itz];
416             
417           } // end loop for (Int_t itz=0; itz<15;itz++)
418           
419           meantzero=meantzero/sumAllweights; // it is given in [ns]
420           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
421           
422           // calculate chisquare
423           
424           Float_t chisquare=0.;         
425           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
426             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
427             
428           } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
429           
430           if(chisquare<=chisquarebest){
431             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
432               
433               bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
434               besttimezero[iqsq]=timezero[iqsq]; 
435               bestmomentum[iqsq]=momentum[iqsq]; 
436               besttimeofflight[iqsq]=timeofflight[iqsq]; 
437               besttexp[iqsq]=texp[iqsq]; 
438               bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
439               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; 
440             }
441             
442             Int_t npion=0;
443             for (Int_t j=0; j<ntracksinsetmy; j++) {
444               assparticle[j]=imass[j];
445               if(imass[j] == 0) npion++;
446             }
447             npionbest=npion;
448             chisquarebest=chisquare;          
449             t0best=meantzero;
450             eT0best=eMeanTzero;
451           } // close if(dummychisquare<=chisquare)
452           
453         }
454         
455         Double_t chi2cut[nmaxtracksinset];
456         chi2cut[0] = 0;
457         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
458         for (Int_t j=2; j<ntracksinset; j++) {
459           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
460         }
461         
462         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
463         
464 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
465         
466         Bool_t kRedoT0 = kFALSE;
467         ntracksinsetmyCut = ntracksinsetmy;
468         Bool_t usetrack[nmaxtracksinset];
469         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
470           usetrack[icsq] = kTRUE;
471           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
472             kRedoT0 = kTRUE;
473             ntracksinsetmyCut--;
474             usetrack[icsq] = kFALSE;
475           }
476         } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
477         
478         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
479         
480         // Loop on mass hypotheses Redo
481         if(kRedoT0 && ntracksinsetmyCut > 1){
482           //      printf("Redo T0\n");
483           for (Int_t k=0; k < ncombinatorial;k++) {
484             for (Int_t j=0; j<ntracksinsetmy; j++) {
485               imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
486               texp[j]=exptof[j][imass[j]];
487               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
488             }
489             
490             Float_t sumAllweights=0.;
491             Float_t meantzero=0.;
492             Float_t eMeanTzero=0.;
493             
494             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
495               if(! usetrack[itz]) continue;
496               sqTrackError[itz]=
497                 (timeresolutioninns*
498                  timeresolutioninns
499                  +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
500               
501               timezero[itz]=texp[itz]-timeofflight[itz];// in ns                          
502               
503               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
504               sumAllweights+=1./sqTrackError[itz];
505               meantzero+=weightedtimezero[itz];
506               
507             } // end loop for (Int_t itz=0; itz<15;itz++)
508             
509             meantzero=meantzero/sumAllweights; // it is given in [ns]
510             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
511             
512             // calculate chisquare
513             
514             Float_t chisquare=0.;               
515             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
516               if(! usetrack[icsq]) continue;
517               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
518               
519             } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
520             
521             Int_t npion=0;
522             for (Int_t j=0; j<ntracksinsetmy; j++) {
523               assparticle[j]=imass[j];
524               if(imass[j] == 0) npion++;
525             }
526             
527             if(chisquare<=chisquarebest){
528               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
529                 if(! usetrack[iqsq]) continue;
530                 bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
531                 besttimezero[iqsq]=timezero[iqsq]; 
532                 bestmomentum[iqsq]=momentum[iqsq]; 
533                 besttimeofflight[iqsq]=timeofflight[iqsq]; 
534                 besttexp[iqsq]=texp[iqsq]; 
535                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
536                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; 
537               }
538               
539               npionbest=npion;
540               chisquarebest=chisquare;        
541               t0best=meantzero;
542               eT0best=eMeanTzero;
543             } // close if(dummychisquare<=chisquare)
544             
545           }
546         }
547                 
548         // filling histos
549         Float_t confLevel=999;
550         
551         // Sets with decent chisquares
552         
553         if(chisquarebest<999.){
554           Double_t dblechisquare=(Double_t)chisquarebest;
555           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); 
556 //        cout << " Set Number " << nsets << endl;      
557 //        cout << "Best Assignment, selection " << assparticle[0] << 
558 //          assparticle[1] << assparticle[2] << 
559 //          assparticle[3] << assparticle[4] << 
560 //          assparticle[5] << endl;
561 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
562 //        cout << " C.L. of the set "<< confLevel <<endl;
563 //        cout << " T0 for this set (in ns)  " << t0best << endl;
564
565           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
566
567             if(! usetrack[icsq]) continue;
568             
569 //          cout << "Track # " << icsq  << " T0 offsets = " 
570 //               << besttimezero[icsq]-t0best << 
571 //            " track error = "  << bestsqTrackError[icsq]
572 //               << " Chisquare = " << bestchisquare[icsq] 
573 //               << " Momentum  = " << bestmomentum[icsq] 
574 //               << " TOF   = "     << besttimeofflight[icsq] 
575 //               << " TOF tracking  = " << besttexp[icsq]
576 //               << " is used = " << usetrack[icsq] << endl;
577           }
578           
579           // Pick up only those with C.L. >1%
580           //      if(confLevel>0.01 && ngoodsetsSel<200){
581           if(confLevel>0.01 && ngoodsetsSel<200){
582             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
583             confLevelbestSel[ngoodsetsSel]=confLevel;
584             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
585             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
586             t0bestallSel += t0best/eT0best/eT0best;
587             sumWt0bestallSel += 1./eT0best/eT0best;
588             ngoodsetsSel++;
589             ngoodtrktrulyused+=ntracksinsetmyCut;           
590           }
591           else{
592             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
593           }
594         }       
595         delete[] tracksT0;
596         nsets++;
597         
598       } // end for the current set
599       
600       nUsedTracks =  ngoodtrkt0;  
601       if(strstr(option,"all")){
602         if(sumAllweightspi>0.){
603           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
604           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
605         }      
606         
607         if(sumWt0bestallSel>0){
608           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
609           eT0bestallSel = sqrt(1./sumWt0bestallSel);
610           
611         }// end of if(sumWt0bestallSel>0){
612         
613 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
614       }
615       
616       t0def=t0bestallSel;
617       deltat0def=eT0bestallSel;
618       if ((TMath::Abs(t0bestallSel) < 0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
619         t0def=-999; deltat0def=0.600;
620       }
621       
622       fT0SigmaT0def[0]=t0def;
623       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*ngoodtrktrulyused/(ngoodtrktrulyused-1));
624       fT0SigmaT0def[2]=ngoodtrkt0;
625       fT0SigmaT0def[3]=ngoodtrktrulyused;
626     }
627   }
628   
629   //   if(strstr(option,"tim") || strstr(option,"all")){
630   //     cout << "AliTOFT0v1:" << endl ;
631   //}
632   
633   if(fT0SigmaT0def[1] < 0.01) fT0SigmaT0def[1] = 0.6;
634
635   return fT0SigmaT0def;
636   }
637 //__________________________________________________________________
638 Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut) 
639
640   // Caluclate the Event Time using the ESD TOF time
641
642   fT0SigmaT0def[0]=0.;
643   fT0SigmaT0def[1]=0.600;
644   fT0SigmaT0def[2]=0.;
645   fT0SigmaT0def[3]=0.;
646
647  Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
648   
649   const Int_t nmaxtracksinset=10;
650 //   if(strstr(option,"all")){
651 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
652 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
653 //   }
654   
655   
656   Int_t nsets=0;
657   Int_t nUsedTracks=0;
658   Int_t ngoodsetsSel= 0;
659   Float_t t0bestSel[300];
660   Float_t eT0bestSel[300];
661   Float_t chiSquarebestSel[300];
662   Float_t confLevelbestSel[300];
663   Float_t t0bestallSel=0.;
664   Float_t eT0bestallSel=0.;
665   Float_t sumWt0bestallSel=0.;
666   Float_t eMeanTzeroPi=0.;
667   Float_t meantzeropi=0.;
668   Float_t sumAllweightspi=0.;
669   Double_t t0def=-999;
670   Double_t deltat0def=999;
671   Int_t ngoodtrktrulyused=0;
672   Int_t ntracksinsetmyCut = 0;
673
674   Int_t ntrk=fEvent->GetNumberOfTracks();
675   
676   AliESDtrack **tracks=new AliESDtrack*[ntrk];
677   Int_t ngoodtrk=0;
678   Int_t ngoodtrkt0 =0;
679   Float_t mintime =1E6;
680   
681   // First Track loop, Selection of good tracks
682
683   for (Int_t itrk=0; itrk<ntrk; itrk++) {
684     AliESDtrack *t=fEvent->GetTrack(itrk);
685     Double_t momOld=t->GetP();
686     Double_t mom=momOld-0.0036*momOld;
687     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
688     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
689     Double_t time=t->GetTOFsignal();
690     
691     time*=1.E-3; // tof given in nanoseconds       
692     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
693    
694     if (!AcceptTrack(t)) continue;
695
696     if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
697     if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
698     if(time <= mintime) mintime=time;
699     tracks[ngoodtrk]=t;
700     ngoodtrk++;
701   }
702   
703   
704 //    cout << " N. of ESD tracks                    : " << ntrk << endl;
705 //    cout << " N. of preselected tracks            : " << ngoodtrk << endl;
706 //    cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
707   
708   AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
709   
710   for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
711     AliESDtrack *t=tracks[jtrk];
712     Double_t time=t->GetTOFsignal();
713
714     if((time-mintime*1.E3)<50.E3){ // For pp and per 
715       gtracks[ngoodtrkt0]=t;
716       ngoodtrkt0++;
717     }
718   }
719   
720
721   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
722   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
723   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
724
725   if(ngoodtrkt0<2){
726 //     cout << "less than 2 tracks, skip event " << endl;
727     t0def=-999;
728     deltat0def=0.600;
729     fT0SigmaT0def[0]=t0def;
730     fT0SigmaT0def[1]=deltat0def;
731     fT0SigmaT0def[2]=ngoodtrkt0;
732     fT0SigmaT0def[3]=ngoodtrkt0;
733     //goto finish;
734   }
735   if(ngoodtrkt0>=2){
736   // Decide how many tracks in set 
737     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
738     Int_t nset=1;
739
740     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} 
741         
742     // Loop over selected sets
743     
744     if(nset>=1){
745       for (Int_t i=0; i< nset; i++) {   
746         
747         Float_t t0best=999.;
748         Float_t eT0best=999.;
749         Float_t chisquarebest=99999.;
750         Int_t npionbest=0;
751         
752         Int_t ntracksinsetmy=0;      
753         AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
754         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
755           Int_t index = itrk+i*ntracksinset;
756           if(index < ngoodtrkt0){
757             AliESDtrack *t=gtracks[index];
758             tracksT0[itrk]=t;
759             ntracksinsetmy++;
760           }
761         }
762         
763         // Analyse it
764         
765         Int_t   assparticle[nmaxtracksinset];
766         Float_t exptof[nmaxtracksinset][3];
767         Float_t timeofflight[nmaxtracksinset];
768         Float_t momentum[nmaxtracksinset];
769         Float_t timezero[nmaxtracksinset];
770         Float_t weightedtimezero[nmaxtracksinset];
771         Float_t beta[nmaxtracksinset];
772         Float_t texp[nmaxtracksinset];
773         Float_t dtexp[nmaxtracksinset];
774         Float_t sqMomError[nmaxtracksinset];
775         Float_t sqTrackError[nmaxtracksinset];
776         Float_t massarray[3]={0.13957,0.493677,0.9382723};
777         Float_t tracktoflen[nmaxtracksinset];
778         Float_t besttimezero[nmaxtracksinset];
779         Float_t besttexp[nmaxtracksinset];
780         Float_t besttimeofflight[nmaxtracksinset];
781         Float_t bestmomentum[nmaxtracksinset];
782         Float_t bestchisquare[nmaxtracksinset];
783         Float_t bestweightedtimezero[nmaxtracksinset];
784         Float_t bestsqTrackError[nmaxtracksinset];
785         Int_t imass[nmaxtracksinset];
786         
787         for (Int_t j=0; j<ntracksinset; j++) {
788           assparticle[j] = 3;
789           timeofflight[j] = 0;
790           momentum[j] = 0;
791           timezero[j] = 0;
792           weightedtimezero[j] = 0;
793           beta[j] = 0;
794           texp[j] = 0;
795           dtexp[j] = 0;
796           sqMomError[j] = 0;
797           sqTrackError[j] = 0;
798           tracktoflen[j] = 0;
799           besttimezero[j] = 0;
800           besttexp[j] = 0;
801           besttimeofflight[j] = 0;
802           bestmomentum[j] = 0;
803           bestchisquare[j] = 0;
804           bestweightedtimezero[j] = 0;
805           bestsqTrackError[j] = 0;
806           imass[j] = 1;
807         }
808         
809         for (Int_t j=0; j<ntracksinsetmy; j++) {
810           AliESDtrack *t=tracksT0[j];
811           Double_t momOld=t->GetP();
812           Double_t mom=momOld-0.0036*momOld;
813           Double_t time=t->GetTOFsignal();
814           
815           time*=1.E-3; // tof given in nanoseconds         
816           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
817           Double_t toflen=t->GetIntegratedLength();
818           toflen=toflen/100.; // toflen given in m 
819           
820           timeofflight[j]=time;
821           tracktoflen[j]=toflen;
822           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
823           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
824           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
825           momentum[j]=mom;
826           assparticle[j]=3;
827           
828         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
829         
830         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
831           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
832                                        +momentum[itz]*momentum[itz]);
833           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])); 
834           sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
835           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
836           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
837           sumAllweightspi+=1./sqTrackError[itz];
838           meantzeropi+=weightedtimezero[itz];   
839         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
840         
841         
842         // Then, Combinatorial Algorithm
843         
844         if(ntracksinsetmy<2 )break;
845         
846         for (Int_t j=0; j<ntracksinsetmy; j++) {
847           imass[j] = 3;
848         }
849         
850         Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
851         
852         // Loop on mass hypotheses
853         for (Int_t k=0; k < ncombinatorial;k++) {
854           for (Int_t j=0; j<ntracksinsetmy; j++) {
855             imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
856             texp[j]=exptof[j][imass[j]];
857             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
858           }
859           Float_t sumAllweights=0.;
860           Float_t meantzero=0.;
861           Float_t eMeanTzero=0.;
862           
863           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
864             sqTrackError[itz]=
865               (timeresolutioninns*
866                timeresolutioninns
867                +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
868             
869             timezero[itz]=texp[itz]-timeofflight[itz];// in ns                    
870             
871             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
872             sumAllweights+=1./sqTrackError[itz];
873             meantzero+=weightedtimezero[itz];
874             
875           } // end loop for (Int_t itz=0; itz<15;itz++)
876           
877           meantzero=meantzero/sumAllweights; // it is given in [ns]
878           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
879           
880           // calculate chisquare
881           
882           Float_t chisquare=0.;         
883           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
884             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
885             
886           } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
887           
888           if(chisquare<=chisquarebest){
889             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
890               
891               bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
892               besttimezero[iqsq]=timezero[iqsq]; 
893               bestmomentum[iqsq]=momentum[iqsq]; 
894               besttimeofflight[iqsq]=timeofflight[iqsq]; 
895               besttexp[iqsq]=texp[iqsq]; 
896               bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
897               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; 
898             }
899             
900             Int_t npion=0;
901             for (Int_t j=0; j<ntracksinsetmy; j++) {
902               assparticle[j]=imass[j];
903               if(imass[j] == 0) npion++;
904             }
905             npionbest=npion;
906             chisquarebest=chisquare;          
907             t0best=meantzero;
908             eT0best=eMeanTzero;
909           } // close if(dummychisquare<=chisquare)
910           
911         }
912         
913         Double_t chi2cut[nmaxtracksinset];
914         chi2cut[0] = 0;
915         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
916         for (Int_t j=2; j<ntracksinset; j++) {
917           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
918         }
919         
920         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
921         
922 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
923         
924         Bool_t kRedoT0 = kFALSE;
925         ntracksinsetmyCut = ntracksinsetmy;
926         Bool_t usetrack[nmaxtracksinset];
927         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
928           usetrack[icsq] = kTRUE;
929           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
930             kRedoT0 = kTRUE;
931             ntracksinsetmyCut--;
932             usetrack[icsq] = kFALSE;
933           }
934         } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
935         
936         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
937         
938         // Loop on mass hypotheses Redo
939         if(kRedoT0 && ntracksinsetmyCut > 1){
940           //      printf("Redo T0\n");
941           for (Int_t k=0; k < ncombinatorial;k++) {
942             for (Int_t j=0; j<ntracksinsetmy; j++) {
943               imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
944               texp[j]=exptof[j][imass[j]];
945               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
946             }
947             
948             Float_t sumAllweights=0.;
949             Float_t meantzero=0.;
950             Float_t eMeanTzero=0.;
951             
952             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
953               if(! usetrack[itz]) continue;
954               sqTrackError[itz]=
955                 (timeresolutioninns*
956                  timeresolutioninns
957                  +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
958               
959               timezero[itz]=texp[itz]-timeofflight[itz];// in ns                          
960               
961               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
962               sumAllweights+=1./sqTrackError[itz];
963               meantzero+=weightedtimezero[itz];
964               
965             } // end loop for (Int_t itz=0; itz<15;itz++)
966             
967             meantzero=meantzero/sumAllweights; // it is given in [ns]
968             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
969             
970             // calculate chisquare
971             
972             Float_t chisquare=0.;               
973             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
974               if(! usetrack[icsq]) continue;
975               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
976               
977             } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
978             
979             Int_t npion=0;
980             for (Int_t j=0; j<ntracksinsetmy; j++) {
981               assparticle[j]=imass[j];
982               if(imass[j] == 0) npion++;
983             }
984             
985             if(chisquare<=chisquarebest){
986               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
987                 if(! usetrack[iqsq]) continue;
988                 bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
989                 besttimezero[iqsq]=timezero[iqsq]; 
990                 bestmomentum[iqsq]=momentum[iqsq]; 
991                 besttimeofflight[iqsq]=timeofflight[iqsq]; 
992                 besttexp[iqsq]=texp[iqsq]; 
993                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
994                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; 
995               }
996               
997               npionbest=npion;
998               chisquarebest=chisquare;        
999               t0best=meantzero;
1000               eT0best=eMeanTzero;
1001             } // close if(dummychisquare<=chisquare)
1002             
1003           }
1004         }
1005                 
1006         // filling histos
1007         Float_t confLevel=999;
1008         
1009         // Sets with decent chisquares
1010         
1011         if(chisquarebest<999.){
1012           Double_t dblechisquare=(Double_t)chisquarebest;
1013           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); 
1014 //        cout << " Set Number " << nsets << endl;      
1015 //        cout << "Best Assignment, selection " << assparticle[0] << 
1016 //          assparticle[1] << assparticle[2] << 
1017 //          assparticle[3] << assparticle[4] << 
1018 //          assparticle[5] << endl;
1019 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
1020 //        cout << " C.L. of the set "<< confLevel <<endl;
1021 //        cout << " T0 for this set (in ns)  " << t0best << endl;
1022
1023           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1024
1025             if(! usetrack[icsq]) continue;
1026             
1027 //          cout << "Track # " << icsq  << " T0 offsets = " 
1028 //               << besttimezero[icsq]-t0best << 
1029 //            " track error = "  << bestsqTrackError[icsq]
1030 //               << " Chisquare = " << bestchisquare[icsq] 
1031 //               << " Momentum  = " << bestmomentum[icsq] 
1032 //               << " TOF   = "     << besttimeofflight[icsq] 
1033 //               << " TOF tracking  = " << besttexp[icsq]
1034 //               << " is used = " << usetrack[icsq] << endl;
1035           }
1036           
1037           // Pick up only those with C.L. >1%
1038           //      if(confLevel>0.01 && ngoodsetsSel<200){
1039           if(confLevel>0.01 && ngoodsetsSel<200){
1040             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
1041             confLevelbestSel[ngoodsetsSel]=confLevel;
1042             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
1043             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
1044             t0bestallSel += t0best/eT0best/eT0best;
1045             sumWt0bestallSel += 1./eT0best/eT0best;
1046             ngoodsetsSel++;
1047             ngoodtrktrulyused+=ntracksinsetmyCut;           
1048           }
1049           else{
1050             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
1051           }
1052         }       
1053         delete[] tracksT0;
1054         nsets++;
1055         
1056       } // end for the current set
1057       
1058       nUsedTracks =  ngoodtrkt0;  
1059       if(strstr(option,"all")){
1060         if(sumAllweightspi>0.){
1061           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
1062           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
1063         }      
1064         
1065         if(sumWt0bestallSel>0){
1066           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
1067           eT0bestallSel = sqrt(1./sumWt0bestallSel);
1068           
1069         }// end of if(sumWt0bestallSel>0){
1070         
1071 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
1072       }
1073       
1074       t0def=t0bestallSel;
1075       deltat0def=eT0bestallSel;
1076       if ((TMath::Abs(t0bestallSel) < 0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
1077         t0def=-999; deltat0def=0.600;
1078       }
1079       
1080       fT0SigmaT0def[0]=t0def;
1081       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*ngoodtrktrulyused/(ngoodtrktrulyused-1));
1082       fT0SigmaT0def[2]=ngoodtrkt0;
1083       fT0SigmaT0def[3]=ngoodtrktrulyused;
1084     }
1085   }
1086   
1087   //   if(strstr(option,"tim") || strstr(option,"all")){
1088   //     cout << "AliTOFT0v1:" << endl ;
1089   //}
1090
1091   if(fT0SigmaT0def[1] < 0.01) fT0SigmaT0def[1] = 0.6;
1092
1093   return fT0SigmaT0def;
1094   }
1095 //__________________________________________________________________
1096 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
1097 {
1098   // Take the error extimate for the TOF time in the track reconstruction
1099
1100   static const Double_t kMasses[]={
1101     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
1102   };
1103
1104   Double_t mass=kMasses[index+2];
1105   Double_t dpp=0.02;      //mean relative pt resolution;
1106   //  if(mom > 1) dpp = 0.02*mom;
1107   Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
1108
1109   sigma =TMath::Sqrt(sigma*sigma);
1110
1111   return sigma;
1112 }
1113
1114 //__________________________________________________________________
1115 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
1116 {
1117
1118   /* TPC refit */
1119   if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1120   /* do not accept kink daughters */
1121   if (track->GetKinkIndex(0)>0) return kFALSE;
1122   /* N clusters TPC */
1123   if (track->GetTPCclusters(0) < 50) return kFALSE;
1124   /* chi2 TPC */
1125   if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
1126   /* sigma to vertex */
1127   if (GetSigmaToVertex(track) > 4.) return kFALSE;
1128   
1129   /* accept track */
1130   return kTRUE;
1131
1132 }
1133
1134 //____________________________________________________________________
1135 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
1136 {
1137   // Calculates the number of sigma to the vertex.
1138
1139   Float_t b[2];
1140   Float_t bRes[2];
1141   Float_t bCov[3];
1142   esdTrack->GetImpactParameters(b,bCov);
1143   
1144   if (bCov[0]<=0 || bCov[2]<=0) {
1145     bCov[0]=0; bCov[2]=0;
1146   }
1147   bRes[0] = TMath::Sqrt(bCov[0]);
1148   bRes[1] = TMath::Sqrt(bCov[2]);
1149
1150   // -----------------------------------
1151   // How to get to a n-sigma cut?
1152   //
1153   // The accumulated statistics from 0 to d is
1154   //
1155   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1156   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1157   //
1158   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1159   // Can this be expressed in a different way?
1160
1161   if (bRes[0] == 0 || bRes[1] ==0)
1162     return -1;
1163
1164   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1165
1166   // work around precision problem
1167   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1168   // 1e-15 corresponds to nsigma ~ 7.7
1169   if (TMath::Exp(-d * d / 2) < 1e-15)
1170     return 1000;
1171
1172   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1173   return nSigma;
1174 }