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