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