1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
18 //_________________________________________________________________________
19 // This is a TTask that made the calculation of the Time zero using TOF.
20 // Description: The algorithm used to calculate the time zero of interaction
21 // using TOF detector is the following.
22 // We select in the ESD some "primary" particles - or tracks in the following -
23 // that strike the TOF detector (the larger part are pions, kaons or protons).
24 // We choose a set of 10 selected tracks, for each track You have the length
25 // of the track when the TOF is reached,
26 // the momentum and the time of flight
27 // given by the TOF detector.
28 // Let consider now only one set of 10 tracks (the algorithm is the same for all sets).
29 // Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton,
30 // we consider all the 3 at 10 possible cases.
31 // For each track in each (mass) configuration
32 // (a configuration can be e.g. pion/pion/kaon/proton/pion/proton/kaon/kaon/pion/pion)
33 // we calculate the time zero (we know in fact the velocity of the track after
34 // the assumption about its mass, the time of flight given by the TOF, and the
35 // corresponding path travelled till the TOF detector). Then for each mass configuration we have
36 // 10 time zero and we can calculate the ChiSquare for the current configuration using the
37 // weighted mean over all 10 time zero.
38 // We call the best assignment the mass configuration that gives the minimum value of the ChiSquare.
39 // We plot the weighted mean over all 10 time zero for the best assignment,
40 // the ChiSquare for the best assignment and the corresponding confidence level.
41 // The strong assumption is the MC selection of primary particles. It will be introduced
42 // in the future also some more realistic simulation about this point.
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
53 // Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions
54 //-- Author: F. Pierella
55 //-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
56 //////////////////////////////////////////////////////////////////////////////
58 #include "Riostream.h"
60 #include "AliTOFT0v1.h"
61 #include "AliESDtrack.h"
62 #include "AliTOFcalibHisto.h"
63 #include "AliESDEvent.h"
67 //____________________________________________________________________________
68 AliTOFT0v1::AliTOFT0v1():
71 fTimeResolution(0.80e-10),
77 // default constructor
80 fT0SigmaT0def[0]=-999.;
81 fT0SigmaT0def[1]=999.;
82 fT0SigmaT0def[2]=-999.;
83 fT0SigmaT0def[3]=-999.;
88 //____________________________________________________________________________
89 AliTOFT0v1::AliTOFT0v1(AliESDEvent* event):
92 fTimeResolution(0.80e-10),
101 fT0SigmaT0def[0]=-999.;
102 fT0SigmaT0def[1]= 999.;
103 fT0SigmaT0def[2]=-999.;
104 fT0SigmaT0def[3]=-999.;
108 //____________________________________________________________________________
109 AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero):
111 fLowerMomBound(tzero.fLowerMomBound),
112 fUpperMomBound(tzero.fUpperMomBound),
113 fTimeResolution(tzero.fTimeResolution),
114 fTimeCorr(tzero.fTimeCorr),
115 fEvent(tzero.fEvent),
122 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
123 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
124 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
125 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
129 //____________________________________________________________________________
130 AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
139 fLowerMomBound=tzero.fLowerMomBound;
140 fUpperMomBound=tzero.fUpperMomBound;
141 fTimeResolution=tzero.fTimeResolution;
142 fTimeCorr=tzero.fTimeCorr;
145 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
146 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
147 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
148 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
152 //____________________________________________________________________________
153 AliTOFT0v1::~AliTOFT0v1()
160 //____________________________________________________________________________
161 void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
162 // Set the TOF time resolution
163 fTimeResolution=timeresolution;
165 //____________________________________________________________________________
166 //____________________________________________________________________________
167 Double_t * AliTOFT0v1::DefineT0(Option_t *option)
169 // Caluclate the Event Time using the ESD TOF time
171 Float_t timeresolutioninns=fTimeResolution*(1.e+9) * TMath::Sqrt(2.); // convert in [ns]
173 const Int_t nmaxtracksinset=10;
174 // if(strstr(option,"all")){
175 // cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
176 // cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
182 Int_t ngoodsetsSel= 0;
183 Float_t t0bestSel[300];
184 Float_t eT0bestSel[300];
185 Float_t chiSquarebestSel[300];
186 Float_t confLevelbestSel[300];
187 Float_t t0bestallSel=0.;
188 Float_t eT0bestallSel=0.;
189 Float_t sumWt0bestallSel=0.;
190 Float_t eMeanTzeroPi=0.;
191 Float_t meantzeropi=0.;
192 Float_t sumAllweightspi=0.;
194 Double_t deltat0def=999;
195 Int_t ngoodtrktrulyused=0;
196 Int_t ntracksinsetmyCut = 0;
198 Int_t ntrk=fEvent->GetNumberOfTracks();
200 AliESDtrack **tracks=new AliESDtrack*[ntrk];
203 Float_t mintime =1E6;
205 // First Track loop, Selection of good tracks
207 for (Int_t itrk=0; itrk<ntrk; itrk++) {
208 AliESDtrack *t=fEvent->GetTrack(itrk);
209 Double_t momOld=t->GetP();
210 Double_t mom=momOld-0.0036*momOld;
211 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
212 Double_t time=t->GetTOFsignal();
214 time*=1.E-3; // tof given in nanoseconds
215 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
217 if (!AcceptTrack(t)) continue;
219 if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
220 if(time <= mintime) mintime=time;
226 // cout << " N. of ESD tracks : " << ntrk << endl;
227 // cout << " N. of preselected tracks : " << ngoodtrk << endl;
228 // cout << " Minimum tof time in set (in ns) : " << mintime << endl;
230 AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
232 for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
233 AliESDtrack *t=tracks[jtrk];
234 Double_t time=t->GetTOFsignal();
236 if((time-mintime*1.E3)<50.E3){ // For pp and per
237 gtracks[ngoodtrkt0]=t;
243 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
244 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
245 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
248 // cout << "less than 2 tracks, skip event " << endl;
251 fT0SigmaT0def[0]=t0def;
252 fT0SigmaT0def[1]=deltat0def;
253 fT0SigmaT0def[2]=ngoodtrkt0;
254 fT0SigmaT0def[3]=ngoodtrkt0;
258 // Decide how many tracks in set
259 Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent);
262 if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
264 // Loop over selected sets
267 for (Int_t i=0; i< nset; i++) {
270 Float_t eT0best=999.;
271 Float_t chisquarebest=99999.;
274 Int_t ntracksinsetmy=0;
275 AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
276 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
277 Int_t index = itrk+i*ntracksinset;
278 if(index < ngoodtrkt0){
279 AliESDtrack *t=gtracks[index];
287 Int_t assparticle[nmaxtracksinset];
288 Float_t exptof[nmaxtracksinset][3];
289 Float_t timeofflight[nmaxtracksinset];
290 Float_t momentum[nmaxtracksinset];
291 Float_t timezero[nmaxtracksinset];
292 Float_t weightedtimezero[nmaxtracksinset];
293 Float_t beta[nmaxtracksinset];
294 Float_t texp[nmaxtracksinset];
295 Float_t dtexp[nmaxtracksinset];
296 Float_t sqMomError[nmaxtracksinset];
297 Float_t sqTrackError[nmaxtracksinset];
298 Float_t massarray[3]={0.13957,0.493677,0.9382723};
299 Float_t tracktoflen[nmaxtracksinset];
300 Float_t besttimezero[nmaxtracksinset];
301 Float_t besttexp[nmaxtracksinset];
302 Float_t besttimeofflight[nmaxtracksinset];
303 Float_t bestmomentum[nmaxtracksinset];
304 Float_t bestchisquare[nmaxtracksinset];
305 Float_t bestweightedtimezero[nmaxtracksinset];
306 Float_t bestsqTrackError[nmaxtracksinset];
307 Int_t imass[nmaxtracksinset];
309 for (Int_t j=0; j<ntracksinset; j++) {
314 weightedtimezero[j] = 0;
323 besttimeofflight[j] = 0;
325 bestchisquare[j] = 0;
326 bestweightedtimezero[j] = 0;
327 bestsqTrackError[j] = 0;
331 for (Int_t j=0; j<ntracksinsetmy; j++) {
332 AliESDtrack *t=tracksT0[j];
333 Double_t momOld=t->GetP();
334 Double_t mom=momOld-0.0036*momOld;
335 Double_t time=t->GetTOFsignal();
337 time*=1.E-3; // tof given in nanoseconds
338 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
339 Double_t toflen=t->GetIntegratedLength();
340 toflen=toflen/100.; // toflen given in m
342 timeofflight[j]=time;
343 tracktoflen[j]=toflen;
344 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
345 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
346 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
350 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
352 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
353 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
354 +momentum[itz]*momentum[itz]);
355 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]));
356 sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
357 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
358 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
359 sumAllweightspi+=1./sqTrackError[itz];
360 meantzeropi+=weightedtimezero[itz];
361 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
364 // Then, Combinatorial Algorithm
366 if(ntracksinsetmy<2 )break;
368 for (Int_t j=0; j<ntracksinsetmy; j++) {
372 Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
374 // Loop on mass hypotheses
375 for (Int_t k=0; k < ncombinatorial;k++) {
376 for (Int_t j=0; j<ntracksinsetmy; j++) {
377 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
378 texp[j]=exptof[j][imass[j]];
379 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
381 Float_t sumAllweights=0.;
382 Float_t meantzero=0.;
383 Float_t eMeanTzero=0.;
385 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
389 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
391 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
393 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
394 sumAllweights+=1./sqTrackError[itz];
395 meantzero+=weightedtimezero[itz];
397 } // end loop for (Int_t itz=0; itz<15;itz++)
399 meantzero=meantzero/sumAllweights; // it is given in [ns]
400 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
402 // calculate chisquare
404 Float_t chisquare=0.;
405 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
406 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
408 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
410 if(chisquare<=chisquarebest){
411 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
413 bestsqTrackError[iqsq]=sqTrackError[iqsq];
414 besttimezero[iqsq]=timezero[iqsq];
415 bestmomentum[iqsq]=momentum[iqsq];
416 besttimeofflight[iqsq]=timeofflight[iqsq];
417 besttexp[iqsq]=texp[iqsq];
418 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
419 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
423 for (Int_t j=0; j<ntracksinsetmy; j++) {
424 assparticle[j]=imass[j];
425 if(imass[j] == 0) npion++;
428 chisquarebest=chisquare;
431 } // close if(dummychisquare<=chisquare)
435 Double_t chi2cut[nmaxtracksinset];
437 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
438 for (Int_t j=2; j<ntracksinset; j++) {
439 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
442 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
444 // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
446 Bool_t kRedoT0 = kFALSE;
447 ntracksinsetmyCut = ntracksinsetmy;
448 Bool_t usetrack[nmaxtracksinset];
449 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
450 usetrack[icsq] = kTRUE;
451 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
454 usetrack[icsq] = kFALSE;
456 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
458 // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
460 // Loop on mass hypotheses Redo
461 if(kRedoT0 && ntracksinsetmyCut > 1){
462 // printf("Redo T0\n");
463 for (Int_t k=0; k < ncombinatorial;k++) {
464 for (Int_t j=0; j<ntracksinsetmy; j++) {
465 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
466 texp[j]=exptof[j][imass[j]];
467 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
470 Float_t sumAllweights=0.;
471 Float_t meantzero=0.;
472 Float_t eMeanTzero=0.;
474 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
475 if(! usetrack[itz]) continue;
479 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
481 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
483 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
484 sumAllweights+=1./sqTrackError[itz];
485 meantzero+=weightedtimezero[itz];
487 } // end loop for (Int_t itz=0; itz<15;itz++)
489 meantzero=meantzero/sumAllweights; // it is given in [ns]
490 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
492 // calculate chisquare
494 Float_t chisquare=0.;
495 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
496 if(! usetrack[icsq]) continue;
497 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
499 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
502 for (Int_t j=0; j<ntracksinsetmy; j++) {
503 assparticle[j]=imass[j];
504 if(imass[j] == 0) npion++;
507 if(chisquare<=chisquarebest){
508 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
509 if(! usetrack[iqsq]) continue;
510 bestsqTrackError[iqsq]=sqTrackError[iqsq];
511 besttimezero[iqsq]=timezero[iqsq];
512 bestmomentum[iqsq]=momentum[iqsq];
513 besttimeofflight[iqsq]=timeofflight[iqsq];
514 besttexp[iqsq]=texp[iqsq];
515 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
516 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
520 chisquarebest=chisquare;
523 } // close if(dummychisquare<=chisquare)
529 Float_t confLevel=999;
531 // Sets with decent chisquares
533 if(chisquarebest<999.){
534 Double_t dblechisquare=(Double_t)chisquarebest;
535 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
536 // cout << " Set Number " << nsets << endl;
537 // cout << "Best Assignment, selection " << assparticle[0] <<
538 // assparticle[1] << assparticle[2] <<
539 // assparticle[3] << assparticle[4] <<
540 // assparticle[5] << endl;
541 // cout << " Chisquare of the set "<< chisquarebest <<endl;
542 // cout << " C.L. of the set "<< confLevel <<endl;
543 // cout << " T0 for this set (in ns) " << t0best << endl;
545 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
547 if(! usetrack[icsq]) continue;
549 // cout << "Track # " << icsq << " T0 offsets = "
550 // << besttimezero[icsq]-t0best <<
551 // " track error = " << bestsqTrackError[icsq]
552 // << " Chisquare = " << bestchisquare[icsq]
553 // << " Momentum = " << bestmomentum[icsq]
554 // << " TOF = " << besttimeofflight[icsq]
555 // << " TOF tracking = " << besttexp[icsq]
556 // << " is used = " << usetrack[icsq] << endl;
559 // Pick up only those with C.L. >1%
560 // if(confLevel>0.01 && ngoodsetsSel<200){
561 if(confLevel>0.01 && ngoodsetsSel<200){
562 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
563 confLevelbestSel[ngoodsetsSel]=confLevel;
564 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
565 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
566 t0bestallSel += t0best/eT0best/eT0best;
567 sumWt0bestallSel += 1./eT0best/eT0best;
569 ngoodtrktrulyused+=ntracksinsetmyCut;
572 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
578 } // end for the current set
580 nUsedTracks = ngoodtrkt0;
581 if(strstr(option,"all")){
582 if(sumAllweightspi>0.){
583 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
584 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
587 if(sumWt0bestallSel>0){
588 t0bestallSel = t0bestallSel/sumWt0bestallSel;
589 eT0bestallSel = sqrt(1./sumWt0bestallSel);
591 }// end of if(sumWt0bestallSel>0){
593 // cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
597 deltat0def=eT0bestallSel;
598 if ((TMath::Abs(t0bestallSel) < 0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
599 t0def=-999; deltat0def=0.600;
602 fT0SigmaT0def[0]=t0def;
603 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
604 fT0SigmaT0def[2]=ngoodtrkt0;
605 fT0SigmaT0def[3]=ngoodtrktrulyused;
609 // if(strstr(option,"tim") || strstr(option,"all")){
610 // cout << "AliTOFT0v1:" << endl ;
613 return fT0SigmaT0def;
615 //__________________________________________________________________
616 Double_t * AliTOFT0v1::DefineT0RawCorrection(Option_t *option)
618 // Caluclate the Event Time using the RAW+correction TOF time
620 Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
622 const Int_t nmaxtracksinset=10;
623 // if(strstr(option,"all")){
624 // cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
625 // cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
628 Float_t stripmean = 0;
632 Int_t ngoodsetsSel= 0;
633 Float_t t0bestSel[300];
634 Float_t eT0bestSel[300];
635 Float_t chiSquarebestSel[300];
636 Float_t confLevelbestSel[300];
637 Float_t t0bestallSel=0.;
638 Float_t eT0bestallSel=0.;
639 Float_t sumWt0bestallSel=0.;
640 Float_t eMeanTzeroPi=0.;
641 Float_t meantzeropi=0.;
642 Float_t sumAllweightspi=0.;
644 Double_t deltat0def=999;
645 Int_t ngoodtrktrulyused=0;
646 Int_t ntracksinsetmyCut = 0;
648 Int_t ntrk=fEvent->GetNumberOfTracks();
650 AliESDtrack **tracks=new AliESDtrack*[ntrk];
653 Float_t mintime =1E6;
655 // First Track loop, Selection of good tracks
657 for (Int_t itrk=0; itrk<ntrk; itrk++) {
658 AliESDtrack *t=fEvent->GetTrack(itrk);
659 Double_t momOld=t->GetP();
660 Double_t mom=momOld-0.0036*momOld;
661 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
662 Double_t tot = t->GetTOFsignalToT();
663 Double_t time=t->GetTOFsignalRaw();
664 Int_t chan = t->GetTOFCalChannel();
665 Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
668 Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
669 if(crate == 63 || crate == 62){
673 Int_t strip = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kSectorStrip,chan);
675 time*=1.E-3; // tof given in nanoseconds
676 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
678 if (!AcceptTrack(t)) continue;
680 if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
681 if(time <= mintime) mintime=time;
686 if(ngoodtrk) stripmean /= ngoodtrk;
688 // cout << " N. of ESD tracks : " << ntrk << endl;
689 // cout << " N. of preselected tracks : " << ngoodtrk << endl;
690 // cout << " Minimum tof time in set (in ns) : " << mintime << endl;
692 AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
694 for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
695 AliESDtrack *t=tracks[jtrk];
696 Double_t tot = t->GetTOFsignalToT();
697 Double_t time=t->GetTOFsignalRaw();
698 Int_t chan = t->GetTOFCalChannel();
699 Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
702 Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
703 if(crate == 63 || crate == 62){
707 if((time-mintime*1.E3)<50.E3){ // For pp and per
708 gtracks[ngoodtrkt0]=t;
713 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
714 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
715 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
718 Int_t nlastset = (ngoodtrkt0 % nmaxtracksinsetCurrent);
720 while(nlastset-nseteq+1 > 2 ){
721 nmaxtracksinsetCurrent++;
722 nlastset -= nseteq-1;
724 if(nmaxtracksinsetCurrent > nmaxtracksinset) nmaxtracksinsetCurrent = nmaxtracksinset;
728 // cout << "less than 2 tracks, skip event " << endl;
731 fT0SigmaT0def[0]=t0def;
732 fT0SigmaT0def[1]=deltat0def;
733 fT0SigmaT0def[2]=ngoodtrkt0;
734 fT0SigmaT0def[3]=ngoodtrkt0;
738 // Decide how many tracks in set
739 Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent);
741 if(ngoodtrkt0>nmaxtracksinset) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
743 // Loop over selected sets
746 for (Int_t i=0; i< nset; i++) {
749 Float_t eT0best=999.;
750 Float_t chisquarebest=99999.;
753 Int_t ntracksinsetmy=0;
754 AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
755 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
756 Int_t index = itrk+i*ntracksinset;
757 if(index < ngoodtrkt0){
758 AliESDtrack *t=gtracks[index];
766 Int_t assparticle[nmaxtracksinset];
767 Float_t exptof[nmaxtracksinset][3];
768 Float_t timeofflight[nmaxtracksinset];
769 Float_t momentum[nmaxtracksinset];
770 Float_t timezero[nmaxtracksinset];
771 Float_t weightedtimezero[nmaxtracksinset];
772 Float_t beta[nmaxtracksinset];
773 Float_t texp[nmaxtracksinset];
774 Float_t dtexp[nmaxtracksinset];
775 Float_t sqMomError[nmaxtracksinset];
776 Float_t sqTrackError[nmaxtracksinset];
777 Float_t massarray[3]={0.13957,0.493677,0.9382723};
778 Float_t tracktoflen[nmaxtracksinset];
779 Float_t besttimezero[nmaxtracksinset];
780 Float_t besttexp[nmaxtracksinset];
781 Float_t besttimeofflight[nmaxtracksinset];
782 Float_t bestmomentum[nmaxtracksinset];
783 Float_t bestchisquare[nmaxtracksinset];
784 Float_t bestweightedtimezero[nmaxtracksinset];
785 Float_t bestsqTrackError[nmaxtracksinset];
786 Int_t imass[nmaxtracksinset];
788 for (Int_t j=0; j<ntracksinset; j++) {
793 weightedtimezero[j] = 0;
802 besttimeofflight[j] = 0;
804 bestchisquare[j] = 0;
805 bestweightedtimezero[j] = 0;
806 bestsqTrackError[j] = 0;
810 for (Int_t j=0; j<ntracksinsetmy; j++) {
811 AliESDtrack *t=tracksT0[j];
812 Double_t momOld=t->GetP();
813 Double_t mom=momOld-0.0036*momOld;
814 Double_t tot = t->GetTOFsignalToT();
815 Double_t time=t->GetTOFsignalRaw();
816 Int_t chan = t->GetTOFCalChannel();
817 Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
820 Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
821 if(crate == 63 || crate == 62){
825 time*=1.E-3; // tof given in nanoseconds
826 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
827 Double_t toflen=t->GetIntegratedLength();
828 toflen=toflen/100.; // toflen given in m
830 timeofflight[j]=time;
831 tracktoflen[j]=toflen;
832 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
833 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
834 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
838 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
840 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
841 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
842 +momentum[itz]*momentum[itz]);
843 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]));
844 sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
845 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
846 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
847 sumAllweightspi+=1./sqTrackError[itz];
848 meantzeropi+=weightedtimezero[itz];
849 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
852 // Then, Combinatorial Algorithm
854 if(ntracksinsetmy<2 )break;
856 for (Int_t j=0; j<ntracksinsetmy; j++) {
860 Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
862 // Loop on mass hypotheses
863 for (Int_t k=0; k < ncombinatorial;k++) {
864 for (Int_t j=0; j<ntracksinsetmy; j++) {
865 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
866 texp[j]=exptof[j][imass[j]];
867 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
869 Float_t sumAllweights=0.;
870 Float_t meantzero=0.;
871 Float_t eMeanTzero=0.;
872 Double_t sumAllSquare=0.;
874 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
878 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
880 // printf("pt=%f -- TOF res=%f -- track res=%f -- all res=%f\n",momentum[itz],timeresolutioninns,dtexp[itz],TMath::Sqrt(sqTrackError[itz]));
882 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
884 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
885 sumAllSquare += timezero[itz]*timezero[itz];
886 sumAllweights+=1./sqTrackError[itz];
887 meantzero+=weightedtimezero[itz];
889 } // end loop for (Int_t itz=0; itz<15;itz++)
891 meantzero=meantzero/sumAllweights; // it is given in [ns]
892 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
895 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
896 sumAllSquare+= (timezero[itz] - meantzero)*(timezero[itz] - meantzero);
898 // eMeanTzero = TMath::Sqrt(sumAllSquare/ntracksinsetmy);
900 // calculate chisquare
902 Float_t chisquare=0.;
903 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
904 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
906 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
908 if(chisquare<=chisquarebest){
909 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
911 bestsqTrackError[iqsq]=sqTrackError[iqsq];
912 besttimezero[iqsq]=timezero[iqsq];
913 bestmomentum[iqsq]=momentum[iqsq];
914 besttimeofflight[iqsq]=timeofflight[iqsq];
915 besttexp[iqsq]=texp[iqsq];
916 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
917 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
921 for (Int_t j=0; j<ntracksinsetmy; j++) {
922 assparticle[j]=imass[j];
923 if(imass[j] == 0) npion++;
926 chisquarebest=chisquare;
929 } // close if(dummychisquare<=chisquare)
933 Double_t chi2cut[nmaxtracksinset];
935 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
936 for (Int_t j=2; j<ntracksinset; j++) {
937 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
940 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
942 // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
944 Bool_t kRedoT0 = kFALSE;
945 ntracksinsetmyCut = ntracksinsetmy;
946 Bool_t usetrack[nmaxtracksinset];
947 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
948 usetrack[icsq] = kTRUE;
949 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
952 usetrack[icsq] = kFALSE;
954 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
956 // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
958 // Loop on mass hypotheses Redo
959 if(kRedoT0 && ntracksinsetmyCut > 1){
960 // printf("Redo T0\n");
961 for (Int_t k=0; k < ncombinatorial;k++) {
962 for (Int_t j=0; j<ntracksinsetmy; j++) {
963 imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
964 texp[j]=exptof[j][imass[j]];
965 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
968 Float_t sumAllweights=0.;
969 Float_t meantzero=0.;
970 Float_t eMeanTzero=0.;
971 Double_t sumAllSquare=0;
973 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
974 if(! usetrack[itz]) continue;
978 +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
980 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
982 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
983 sumAllweights+=1./sqTrackError[itz];
984 meantzero+=weightedtimezero[itz];
985 } // end loop for (Int_t itz=0; itz<15;itz++)
987 meantzero=meantzero/sumAllweights; // it is given in [ns]
988 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
991 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
992 if(! usetrack[itz]) continue;
993 sumAllSquare+= (timezero[itz] - meantzero)*(timezero[itz] - meantzero);
995 // eMeanTzero = TMath::Sqrt(sumAllSquare/ntracksinsetmyCut);
997 // calculate chisquare
999 Float_t chisquare=0.;
1000 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1001 if(! usetrack[icsq]) continue;
1002 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
1004 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1007 for (Int_t j=0; j<ntracksinsetmy; j++) {
1008 assparticle[j]=imass[j];
1009 if(imass[j] == 0) npion++;
1012 if(chisquare<=chisquarebest){
1013 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
1014 if(! usetrack[iqsq]) continue;
1015 bestsqTrackError[iqsq]=sqTrackError[iqsq];
1016 besttimezero[iqsq]=timezero[iqsq];
1017 bestmomentum[iqsq]=momentum[iqsq];
1018 besttimeofflight[iqsq]=timeofflight[iqsq];
1019 besttexp[iqsq]=texp[iqsq];
1020 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1021 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
1025 chisquarebest=chisquare;
1028 } // close if(dummychisquare<=chisquare)
1033 if(chisquarebest >= 999){
1034 printf("How is it possible (chi2 = %f)? T0best = %f\n",chisquarebest,t0best);
1036 // for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1037 // cout << "Track # " << icsq << " T0 offsets = "
1038 // << besttimezero[icsq]-t0best <<
1039 // " track error = " << bestsqTrackError[icsq]
1040 // << " Chisquare = " << bestchisquare[icsq]
1041 // << " Momentum = " << bestmomentum[icsq]
1042 // << " TOF = " << besttimeofflight[icsq]
1043 // << " TOF tracking = " << besttexp[icsq]
1044 // << " is used = " << usetrack[icsq] << endl;
1049 Float_t confLevel=999;
1051 // Sets with decent chisquares
1053 if(chisquarebest<999.){
1054 Double_t dblechisquare=(Double_t)chisquarebest;
1055 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
1056 // cout << " Set Number " << nsets << endl;
1057 // cout << "Best Assignment, selection " << assparticle[0] <<
1058 // assparticle[1] << assparticle[2] <<
1059 // assparticle[3] << assparticle[4] <<
1060 // assparticle[5] << endl;
1061 // cout << " Chisquare of the set "<< chisquarebest <<endl;
1062 // cout << " C.L. of the set "<< confLevel <<endl;
1063 // cout << " T0 for this set (in ns) " << t0best << endl;
1065 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1067 if(! usetrack[icsq]) continue;
1069 // cout << "Track # " << icsq << " T0 offsets = "
1070 // << besttimezero[icsq]-t0best <<
1071 // " track error = " << bestsqTrackError[icsq]
1072 // << " Chisquare = " << bestchisquare[icsq]
1073 // << " Momentum = " << bestmomentum[icsq]
1074 // << " TOF = " << besttimeofflight[icsq]
1075 // << " TOF tracking = " << besttexp[icsq]
1076 // << " is used = " << usetrack[icsq] << endl;
1079 // Pick up only those with C.L. >1%
1080 // if(confLevel>0.01 && ngoodsetsSel<200){
1081 if(confLevel>0.01 && ngoodsetsSel<200){
1082 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
1083 confLevelbestSel[ngoodsetsSel]=confLevel;
1084 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
1085 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
1086 t0bestallSel += t0best/eT0best/eT0best;
1087 sumWt0bestallSel += 1./eT0best/eT0best;
1090 ngoodtrktrulyused+=ntracksinsetmyCut;
1093 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
1099 } // end for the current set
1101 nUsedTracks = ngoodtrkt0;
1102 if(strstr(option,"all")){
1103 if(sumAllweightspi>0.){
1104 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
1105 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
1108 if(sumWt0bestallSel>0){
1109 t0bestallSel = t0bestallSel/sumWt0bestallSel;
1110 eT0bestallSel = sqrt(1./sumWt0bestallSel);
1111 }// end of if(sumWt0bestallSel>0){
1116 deltat0def=eT0bestallSel;
1117 if ((TMath::Abs(t0bestallSel)<0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
1118 t0def=-999; deltat0def=0.600;
1121 fT0SigmaT0def[0]=t0def;
1122 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
1123 fT0SigmaT0def[2]=ngoodtrkt0;
1124 fT0SigmaT0def[3]=ngoodtrktrulyused;
1128 return fT0SigmaT0def;
1131 //__________________________________________________________________
1132 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
1134 // Take the error extimate for the TOF time in the track reconstruction
1136 static const Double_t kMasses[]={
1137 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
1140 Double_t mass=kMasses[index+2];
1141 Double_t dpp=0.01; //mean relative pt resolution;
1142 Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
1144 sigma =TMath::Sqrt(sigma*sigma);
1149 //__________________________________________________________________
1150 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
1154 if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1155 /* do not accept kink daughters */
1156 if (track->GetKinkIndex(0)>0) return kFALSE;
1157 /* N clusters TPC */
1158 if (track->GetTPCclusters(0) < 50) return kFALSE;
1160 if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
1161 /* sigma to vertex */
1162 if (GetSigmaToVertex(track) > 4.) return kFALSE;
1169 //____________________________________________________________________
1170 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
1172 // Calculates the number of sigma to the vertex.
1177 esdTrack->GetImpactParameters(b,bCov);
1179 if (bCov[0]<=0 || bCov[2]<=0) {
1180 bCov[0]=0; bCov[2]=0;
1182 bRes[0] = TMath::Sqrt(bCov[0]);
1183 bRes[1] = TMath::Sqrt(bCov[2]);
1185 // -----------------------------------
1186 // How to get to a n-sigma cut?
1188 // The accumulated statistics from 0 to d is
1190 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1191 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1193 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1194 // Can this be expressed in a different way?
1196 if (bRes[0] == 0 || bRes[1] ==0)
1199 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1201 // work around precision problem
1202 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1203 // 1e-15 corresponds to nsigma ~ 7.7
1204 if (TMath::Exp(-d * d / 2) < 1e-15)
1207 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);