]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/AliTOFT0v1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Remake TOF-PID with event-time resolution calculated event-by-event (T0 fill or T0...
[u/mrichter/AliRoot.git] / TOF / AliTOFT0v1.cxx
CommitLineData
536031f2 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
6a4e212e 16/* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
536031f2 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.
6a4e212e 22// We select in the ESD some "primary" particles - or tracks in the following -
536031f2 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
6a4e212e 25// of the track when the TOF is reached,
26// the momentum and the time of flight
536031f2 27// given by the TOF detector.
536031f2 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
6a4e212e 55//-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
536031f2 56//////////////////////////////////////////////////////////////////////////////
57
536031f2 58#include "AliESDtrack.h"
8f589502 59#include "AliESDEvent.h"
03bd764d 60#include "AliTOFT0v1.h"
536031f2 61
62ClassImp(AliTOFT0v1)
63
64//____________________________________________________________________________
8f589502 65AliTOFT0v1::AliTOFT0v1():
03bd764d 66 fLowerMomBound(0.5),
67 fUpperMomBound(1.5),
8f589502 68 fTimeResolution(0.80e-10),
69 fTimeCorr(0.),
03bd764d 70 fEvent(0x0)
71// fCalib(0x0)
536031f2 72{
8f589502 73 //
74 // default constructor
75 //
76
536031f2 77 fT0SigmaT0def[0]=-999.;
78 fT0SigmaT0def[1]=999.;
79 fT0SigmaT0def[2]=-999.;
80 fT0SigmaT0def[3]=-999.;
81
8f589502 82}
83
84
85//____________________________________________________________________________
86AliTOFT0v1::AliTOFT0v1(AliESDEvent* event):
03bd764d 87 fLowerMomBound(0.5),
88 fUpperMomBound(1.5),
8f589502 89 fTimeResolution(0.80e-10),
90 fTimeCorr(0.),
03bd764d 91 fEvent(event)
92// fCalib(0x0)
8f589502 93{
94 //
95 // real constructor
96 //
97
98 fT0SigmaT0def[0]=-999.;
99 fT0SigmaT0def[1]= 999.;
100 fT0SigmaT0def[2]=-999.;
101 fT0SigmaT0def[3]=-999.;
102
536031f2 103}
104
105//____________________________________________________________________________
8f589502 106AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero):
107 TObject(),
108 fLowerMomBound(tzero.fLowerMomBound),
109 fUpperMomBound(tzero.fUpperMomBound),
110 fTimeResolution(tzero.fTimeResolution),
111 fTimeCorr(tzero.fTimeCorr),
03bd764d 112 fEvent(tzero.fEvent)
113// fCalib(tzero.fCalib)
536031f2 114{
8f589502 115 //
116 // copy constructor
117 //
118
119 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
120 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
121 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
122 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
123
536031f2 124}
125
8f589502 126//____________________________________________________________________________
127AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
128{
129 //
130 // assign. operator
131 //
132
133 if (this == &tzero)
134 return *this;
135
136 fLowerMomBound=tzero.fLowerMomBound;
137 fUpperMomBound=tzero.fUpperMomBound;
138 fTimeResolution=tzero.fTimeResolution;
139 fTimeCorr=tzero.fTimeCorr;
140 fEvent=tzero.fEvent;
03bd764d 141// fCalib=tzero.fCalib;
8f589502 142 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
143 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
144 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
145 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
146
147 return *this;
148}
536031f2 149//____________________________________________________________________________
150AliTOFT0v1::~AliTOFT0v1()
151{
152 // dtor
03bd764d 153// fCalib=NULL;
8f589502 154 fEvent=NULL;
155
536031f2 156}
8f589502 157//____________________________________________________________________________
536031f2 158void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
8f589502 159 // Set the TOF time resolution
536031f2 160 fTimeResolution=timeresolution;
161}
162//____________________________________________________________________________
163//____________________________________________________________________________
164Double_t * AliTOFT0v1::DefineT0(Option_t *option)
165{
8f589502 166 // Caluclate the Event Time using the ESD TOF time
167
03bd764d 168 Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
536031f2 169
170 const Int_t nmaxtracksinset=10;
6a4e212e 171// if(strstr(option,"all")){
172// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
173// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
174// }
536031f2 175
176
177 Int_t nsets=0;
8f589502 178 Int_t nUsedTracks=0;
536031f2 179 Int_t ngoodsetsSel= 0;
180 Float_t t0bestSel[300];
8f589502 181 Float_t eT0bestSel[300];
182 Float_t chiSquarebestSel[300];
183 Float_t confLevelbestSel[300];
536031f2 184 Float_t t0bestallSel=0.;
8f589502 185 Float_t eT0bestallSel=0.;
536031f2 186 Float_t sumWt0bestallSel=0.;
8f589502 187 Float_t eMeanTzeroPi=0.;
536031f2 188 Float_t meantzeropi=0.;
189 Float_t sumAllweightspi=0.;
190 Double_t t0def=-999;
191 Double_t deltat0def=999;
192 Int_t ngoodtrktrulyused=0;
193 Int_t ntracksinsetmyCut = 0;
194
195 Int_t ntrk=fEvent->GetNumberOfTracks();
196
197 AliESDtrack **tracks=new AliESDtrack*[ntrk];
198 Int_t ngoodtrk=0;
199 Int_t ngoodtrkt0 =0;
200 Float_t mintime =1E6;
201
202 // First Track loop, Selection of good tracks
203
204 for (Int_t itrk=0; itrk<ntrk; itrk++) {
205 AliESDtrack *t=fEvent->GetTrack(itrk);
206 Double_t momOld=t->GetP();
207 Double_t mom=momOld-0.0036*momOld;
208 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
03bd764d 209 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
536031f2 210 Double_t time=t->GetTOFsignal();
211
212 time*=1.E-3; // tof given in nanoseconds
213 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
14b2cbea 214
215 if (!AcceptTrack(t)) continue;
536031f2 216
217 if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
218 if(time <= mintime) mintime=time;
219 tracks[ngoodtrk]=t;
220 ngoodtrk++;
221 }
222
223
03bd764d 224// cout << " N. of ESD tracks : " << ntrk << endl;
225// cout << " N. of preselected tracks : " << ngoodtrk << endl;
226// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
536031f2 227
228 AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
229
230 for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
231 AliESDtrack *t=tracks[jtrk];
232 Double_t time=t->GetTOFsignal();
233
234 if((time-mintime*1.E3)<50.E3){ // For pp and per
235 gtracks[ngoodtrkt0]=t;
236 ngoodtrkt0++;
237 }
238 }
239
240
241 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
242 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
243 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
244
245 if(ngoodtrkt0<2){
6a4e212e 246// cout << "less than 2 tracks, skip event " << endl;
536031f2 247 t0def=-999;
248 deltat0def=0.600;
249 fT0SigmaT0def[0]=t0def;
250 fT0SigmaT0def[1]=deltat0def;
251 fT0SigmaT0def[2]=ngoodtrkt0;
252 fT0SigmaT0def[3]=ngoodtrkt0;
253 //goto finish;
254 }
255 if(ngoodtrkt0>=2){
256 // Decide how many tracks in set
03bd764d 257 Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
536031f2 258 Int_t nset=1;
259
260 if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
6a4e212e 261
536031f2 262 // Loop over selected sets
263
264 if(nset>=1){
265 for (Int_t i=0; i< nset; i++) {
266
267 Float_t t0best=999.;
8f589502 268 Float_t eT0best=999.;
536031f2 269 Float_t chisquarebest=99999.;
270 Int_t npionbest=0;
271
272 Int_t ntracksinsetmy=0;
273 AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
274 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
275 Int_t index = itrk+i*ntracksinset;
276 if(index < ngoodtrkt0){
277 AliESDtrack *t=gtracks[index];
278 tracksT0[itrk]=t;
279 ntracksinsetmy++;
280 }
281 }
282
283 // Analyse it
284
285 Int_t assparticle[nmaxtracksinset];
286 Float_t exptof[nmaxtracksinset][3];
287 Float_t timeofflight[nmaxtracksinset];
288 Float_t momentum[nmaxtracksinset];
289 Float_t timezero[nmaxtracksinset];
290 Float_t weightedtimezero[nmaxtracksinset];
291 Float_t beta[nmaxtracksinset];
292 Float_t texp[nmaxtracksinset];
293 Float_t dtexp[nmaxtracksinset];
294 Float_t sqMomError[nmaxtracksinset];
295 Float_t sqTrackError[nmaxtracksinset];
296 Float_t massarray[3]={0.13957,0.493677,0.9382723};
297 Float_t tracktoflen[nmaxtracksinset];
298 Float_t besttimezero[nmaxtracksinset];
299 Float_t besttexp[nmaxtracksinset];
300 Float_t besttimeofflight[nmaxtracksinset];
301 Float_t bestmomentum[nmaxtracksinset];
302 Float_t bestchisquare[nmaxtracksinset];
303 Float_t bestweightedtimezero[nmaxtracksinset];
304 Float_t bestsqTrackError[nmaxtracksinset];
305 Int_t imass[nmaxtracksinset];
306
307 for (Int_t j=0; j<ntracksinset; j++) {
308 assparticle[j] = 3;
309 timeofflight[j] = 0;
310 momentum[j] = 0;
311 timezero[j] = 0;
312 weightedtimezero[j] = 0;
313 beta[j] = 0;
314 texp[j] = 0;
315 dtexp[j] = 0;
316 sqMomError[j] = 0;
317 sqTrackError[j] = 0;
318 tracktoflen[j] = 0;
319 besttimezero[j] = 0;
320 besttexp[j] = 0;
321 besttimeofflight[j] = 0;
322 bestmomentum[j] = 0;
323 bestchisquare[j] = 0;
324 bestweightedtimezero[j] = 0;
325 bestsqTrackError[j] = 0;
326 imass[j] = 1;
327 }
328
329 for (Int_t j=0; j<ntracksinsetmy; j++) {
330 AliESDtrack *t=tracksT0[j];
331 Double_t momOld=t->GetP();
332 Double_t mom=momOld-0.0036*momOld;
333 Double_t time=t->GetTOFsignal();
334
335 time*=1.E-3; // tof given in nanoseconds
336 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
337 Double_t toflen=t->GetIntegratedLength();
338 toflen=toflen/100.; // toflen given in m
339
8f589502 340 timeofflight[j]=time;
341 tracktoflen[j]=toflen;
536031f2 342 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
343 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
344 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
345 momentum[j]=mom;
346 assparticle[j]=3;
347
348 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
349
536031f2 350 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
351 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
352 +momentum[itz]*momentum[itz]);
353 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]));
354 sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
355 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
356 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
357 sumAllweightspi+=1./sqTrackError[itz];
358 meantzeropi+=weightedtimezero[itz];
359 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
360
361
362 // Then, Combinatorial Algorithm
363
364 if(ntracksinsetmy<2 )break;
365
366 for (Int_t j=0; j<ntracksinsetmy; j++) {
367 imass[j] = 3;
368 }
369
370 Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
371
372 // Loop on mass hypotheses
373 for (Int_t k=0; k < ncombinatorial;k++) {
374 for (Int_t j=0; j<ntracksinsetmy; j++) {
375 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
376 texp[j]=exptof[j][imass[j]];
377 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
378 }
379 Float_t sumAllweights=0.;
380 Float_t meantzero=0.;
8f589502 381 Float_t eMeanTzero=0.;
536031f2 382
383 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
384 sqTrackError[itz]=
385 (timeresolutioninns*
386 timeresolutioninns
387 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
388
389 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
390
391 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
392 sumAllweights+=1./sqTrackError[itz];
393 meantzero+=weightedtimezero[itz];
394
395 } // end loop for (Int_t itz=0; itz<15;itz++)
396
397 meantzero=meantzero/sumAllweights; // it is given in [ns]
8f589502 398 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
536031f2 399
400 // calculate chisquare
401
402 Float_t chisquare=0.;
403 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
404 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
405
406 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
407
408 if(chisquare<=chisquarebest){
409 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
410
411 bestsqTrackError[iqsq]=sqTrackError[iqsq];
412 besttimezero[iqsq]=timezero[iqsq];
413 bestmomentum[iqsq]=momentum[iqsq];
414 besttimeofflight[iqsq]=timeofflight[iqsq];
415 besttexp[iqsq]=texp[iqsq];
416 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
417 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
418 }
419
420 Int_t npion=0;
421 for (Int_t j=0; j<ntracksinsetmy; j++) {
422 assparticle[j]=imass[j];
423 if(imass[j] == 0) npion++;
424 }
425 npionbest=npion;
426 chisquarebest=chisquare;
427 t0best=meantzero;
8f589502 428 eT0best=eMeanTzero;
536031f2 429 } // close if(dummychisquare<=chisquare)
430
431 }
432
433 Double_t chi2cut[nmaxtracksinset];
434 chi2cut[0] = 0;
435 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
436 for (Int_t j=2; j<ntracksinset; j++) {
437 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
438 }
439
440 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
441
03bd764d 442// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
536031f2 443
444 Bool_t kRedoT0 = kFALSE;
445 ntracksinsetmyCut = ntracksinsetmy;
446 Bool_t usetrack[nmaxtracksinset];
447 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
448 usetrack[icsq] = kTRUE;
449 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
450 kRedoT0 = kTRUE;
451 ntracksinsetmyCut--;
452 usetrack[icsq] = kFALSE;
453 }
454 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
455
6a4e212e 456 // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
536031f2 457
458 // Loop on mass hypotheses Redo
459 if(kRedoT0 && ntracksinsetmyCut > 1){
6a4e212e 460 // printf("Redo T0\n");
536031f2 461 for (Int_t k=0; k < ncombinatorial;k++) {
462 for (Int_t j=0; j<ntracksinsetmy; j++) {
463 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
464 texp[j]=exptof[j][imass[j]];
465 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
466 }
467
468 Float_t sumAllweights=0.;
469 Float_t meantzero=0.;
8f589502 470 Float_t eMeanTzero=0.;
536031f2 471
472 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
473 if(! usetrack[itz]) continue;
474 sqTrackError[itz]=
475 (timeresolutioninns*
476 timeresolutioninns
477 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
478
479 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
480
481 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
482 sumAllweights+=1./sqTrackError[itz];
483 meantzero+=weightedtimezero[itz];
484
485 } // end loop for (Int_t itz=0; itz<15;itz++)
486
487 meantzero=meantzero/sumAllweights; // it is given in [ns]
8f589502 488 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
536031f2 489
490 // calculate chisquare
491
492 Float_t chisquare=0.;
493 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
494 if(! usetrack[icsq]) continue;
495 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
496
497 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
498
499 Int_t npion=0;
500 for (Int_t j=0; j<ntracksinsetmy; j++) {
501 assparticle[j]=imass[j];
502 if(imass[j] == 0) npion++;
503 }
504
505 if(chisquare<=chisquarebest){
506 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
507 if(! usetrack[iqsq]) continue;
508 bestsqTrackError[iqsq]=sqTrackError[iqsq];
509 besttimezero[iqsq]=timezero[iqsq];
510 bestmomentum[iqsq]=momentum[iqsq];
511 besttimeofflight[iqsq]=timeofflight[iqsq];
512 besttexp[iqsq]=texp[iqsq];
513 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
514 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
515 }
516
517 npionbest=npion;
518 chisquarebest=chisquare;
519 t0best=meantzero;
8f589502 520 eT0best=eMeanTzero;
536031f2 521 } // close if(dummychisquare<=chisquare)
522
523 }
524 }
6a4e212e 525
536031f2 526 // filling histos
527 Float_t confLevel=999;
528
529 // Sets with decent chisquares
530
531 if(chisquarebest<999.){
532 Double_t dblechisquare=(Double_t)chisquarebest;
533 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
03bd764d 534// cout << " Set Number " << nsets << endl;
535// cout << "Best Assignment, selection " << assparticle[0] <<
536// assparticle[1] << assparticle[2] <<
537// assparticle[3] << assparticle[4] <<
538// assparticle[5] << endl;
539// cout << " Chisquare of the set "<< chisquarebest <<endl;
540// cout << " C.L. of the set "<< confLevel <<endl;
541// cout << " T0 for this set (in ns) " << t0best << endl;
536031f2 542
543 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
544
545 if(! usetrack[icsq]) continue;
546
03bd764d 547// cout << "Track # " << icsq << " T0 offsets = "
548// << besttimezero[icsq]-t0best <<
549// " track error = " << bestsqTrackError[icsq]
550// << " Chisquare = " << bestchisquare[icsq]
551// << " Momentum = " << bestmomentum[icsq]
552// << " TOF = " << besttimeofflight[icsq]
553// << " TOF tracking = " << besttexp[icsq]
554// << " is used = " << usetrack[icsq] << endl;
536031f2 555 }
556
557 // Pick up only those with C.L. >1%
558 // if(confLevel>0.01 && ngoodsetsSel<200){
559 if(confLevel>0.01 && ngoodsetsSel<200){
8f589502 560 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
561 confLevelbestSel[ngoodsetsSel]=confLevel;
562 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
563 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
564 t0bestallSel += t0best/eT0best/eT0best;
565 sumWt0bestallSel += 1./eT0best/eT0best;
536031f2 566 ngoodsetsSel++;
567 ngoodtrktrulyused+=ntracksinsetmyCut;
568 }
569 else{
6a4e212e 570 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
536031f2 571 }
572 }
573 delete[] tracksT0;
574 nsets++;
575
576 } // end for the current set
577
8f589502 578 nUsedTracks = ngoodtrkt0;
536031f2 579 if(strstr(option,"all")){
580 if(sumAllweightspi>0.){
581 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
8f589502 582 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
536031f2 583 }
584
585 if(sumWt0bestallSel>0){
586 t0bestallSel = t0bestallSel/sumWt0bestallSel;
8f589502 587 eT0bestallSel = sqrt(1./sumWt0bestallSel);
536031f2 588
589 }// end of if(sumWt0bestallSel>0){
590
6a4e212e 591// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
536031f2 592 }
593
594 t0def=t0bestallSel;
8f589502 595 deltat0def=eT0bestallSel;
596 if ((TMath::Abs(t0bestallSel) < 0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
536031f2 597 t0def=-999; deltat0def=0.600;
598 }
599
600 fT0SigmaT0def[0]=t0def;
03bd764d 601 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
6a4e212e 602 fT0SigmaT0def[2]=ngoodtrkt0;
536031f2 603 fT0SigmaT0def[3]=ngoodtrktrulyused;
604 }
605 }
606
6a4e212e 607 // if(strstr(option,"tim") || strstr(option,"all")){
608 // cout << "AliTOFT0v1:" << endl ;
609 //}
536031f2 610
611 return fT0SigmaT0def;
612 }
613//__________________________________________________________________
8f589502 614Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
536031f2 615{
8f589502 616 // Take the error extimate for the TOF time in the track reconstruction
536031f2 617
536031f2 618 static const Double_t kMasses[]={
619 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
620 };
621
622 Double_t mass=kMasses[index+2];
623 Double_t dpp=0.01; //mean relative pt resolution;
03bd764d 624 if(mom > 1) dpp = 0.01*mom;
536031f2 625 Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
626
627 sigma =TMath::Sqrt(sigma*sigma);
628
629 return sigma;
630}
14b2cbea 631
632//__________________________________________________________________
633Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
634{
635
636 /* TPC refit */
637 if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
638 /* do not accept kink daughters */
639 if (track->GetKinkIndex(0)>0) return kFALSE;
640 /* N clusters TPC */
641 if (track->GetTPCclusters(0) < 50) return kFALSE;
642 /* chi2 TPC */
643 if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
644 /* sigma to vertex */
645 if (GetSigmaToVertex(track) > 4.) return kFALSE;
646
647 /* accept track */
648 return kTRUE;
649
650}
651
652//____________________________________________________________________
8f589502 653Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
14b2cbea 654{
655 // Calculates the number of sigma to the vertex.
656
657 Float_t b[2];
658 Float_t bRes[2];
659 Float_t bCov[3];
660 esdTrack->GetImpactParameters(b,bCov);
661
662 if (bCov[0]<=0 || bCov[2]<=0) {
663 bCov[0]=0; bCov[2]=0;
664 }
665 bRes[0] = TMath::Sqrt(bCov[0]);
666 bRes[1] = TMath::Sqrt(bCov[2]);
667
668 // -----------------------------------
669 // How to get to a n-sigma cut?
670 //
671 // The accumulated statistics from 0 to d is
672 //
673 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
674 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
675 //
676 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
677 // Can this be expressed in a different way?
678
679 if (bRes[0] == 0 || bRes[1] ==0)
680 return -1;
681
682 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
683
684 // work around precision problem
685 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
686 // 1e-15 corresponds to nsigma ~ 7.7
687 if (TMath::Exp(-d * d / 2) < 1e-15)
688 return 1000;
689
690 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
691 return nSigma;
692}