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