]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TOF/AliTOFT0v1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
#101318: Patch for various problems in AliROOT
[u/mrichter/AliRoot.git] / TOF / AliTOFT0v1.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
17
18//_________________________________________________________________________
19// This is a TTask that made the calculation of the Time zero using TOF.
20// Description: The algorithm used to calculate the time zero of interaction
21// using TOF detector is the following.
22// We select in the ESD some "primary" particles - or tracks in the following -
23// that strike the TOF detector (the larger part are pions, kaons or protons).
24// We choose a set of 10 selected tracks, for each track You have the length
25// of the track when the TOF is reached,
26// the momentum and the time of flight
27// given by the TOF detector.
28// Let consider now only one set of 10 tracks (the algorithm is the same for all sets).
29// Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton,
30// we consider all the 3 at 10 possible cases.
31// For each track in each (mass) configuration
32// (a configuration can be e.g. pion/pion/kaon/proton/pion/proton/kaon/kaon/pion/pion)
33// we calculate the time zero (we know in fact the velocity of the track after
34// the assumption about its mass, the time of flight given by the TOF, and the
35// corresponding path travelled till the TOF detector). Then for each mass configuration we have
36// 10 time zero and we can calculate the ChiSquare for the current configuration using the
37// weighted mean over all 10 time zero.
38// We call the best assignment the mass configuration that gives the minimum value of the ChiSquare.
39// We plot the weighted mean over all 10 time zero for the best assignment,
40// the ChiSquare for the best assignment and the corresponding confidence level.
41// The strong assumption is the MC selection of primary particles. It will be introduced
42// in the future also some more realistic simulation about this point.
43// Use case:
44// root [0] AliTOFT0v1 * tzero = new AliTOFT0v1("galice.root")
45// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
46// root [1] tzero->ExecuteTask()
47// root [2] tzero->ExecuteTask("tim")
48// // available parameters:
49// tim - print benchmarking information
50// all - print usefull informations about the number of misidentified tracks
51// and a comparison about the true configuration (known from MC) and the best
52// assignment
53// Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions
54//-- Author: F. Pierella
55//-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
56//-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table
57// for Power3)
58//////////////////////////////////////////////////////////////////////////////
59
60#include "AliESDtrack.h"
61#include "AliESDEvent.h"
62#include "AliTOFT0v1.h"
63#include "TBenchmark.h"
64#include "AliPID.h"
65#include "AliESDpid.h"
66
67ClassImp(AliTOFT0v1)
68
69//____________________________________________________________________________
70AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
71 TObject(),
72 fLowerMomBound(0.5),
73 fUpperMomBound(3),
74 fTimeCorr(0.),
75 fEvent(0x0),
76 fPIDesd(extPID),
77 fTracks(new TObjArray(10)),
78 fGTracks(new TObjArray(10)),
79 fTracksT0(new TObjArray(10)),
80 fOptFlag(kFALSE)
81{
82 //
83 // default constructor
84 //
85 if(AliPID::ParticleMass(0) == 0) new AliPID();
86
87 if(!fPIDesd){
88 fPIDesd = new AliESDpid();
89 }
90
91 Init(NULL);
92
93 //initialise lookup table for power 3
94 // a set should only have 10 tracks a t maximum
95 // so up to 15 should be more than enough
96 for (Int_t i=0; i<15; ++i) {
97 fLookupPowerThree[i]=ToCalculatePower(3,i);
98 }
99}
100
101//____________________________________________________________________________
102AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
103 TObject(),
104 fLowerMomBound(0.5),
105 fUpperMomBound(3.0),
106 fTimeCorr(0.),
107 fEvent(event),
108 fPIDesd(extPID),
109 fTracks(new TObjArray(10)),
110 fGTracks(new TObjArray(10)),
111 fTracksT0(new TObjArray(10)),
112 fOptFlag(kFALSE)
113{
114 //
115 // real constructor
116 //
117 if(AliPID::ParticleMass(0) == 0) new AliPID();
118
119 if(!fPIDesd){
120 fPIDesd = new AliESDpid();
121 }
122
123 Init(event);
124 //initialise lookup table for power 3
125 for (Int_t i=0; i<15; ++i) {
126 fLookupPowerThree[i]=Int_t(TMath::Power(3,i));
127 }
128}
129//____________________________________________________________________________
130AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
131{
132 //
133 // assign. operator
134 //
135
136 if (this == &tzero)
137 return *this;
138
139 fLowerMomBound=tzero.fLowerMomBound;
140 fUpperMomBound=tzero.fUpperMomBound;
141 fTimeCorr=tzero.fTimeCorr;
142 fEvent=tzero.fEvent;
143 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
144 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
145 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
146 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
147
148 fTracks=tzero.fTracks;
149 fGTracks=tzero.fGTracks;
150 fTracksT0=tzero.fTracksT0;
151
152 for (Int_t ii=0; ii<tzero.fTracks->GetEntries(); ii++)
153 fTracks->AddLast(tzero.fTracks->At(ii));
154
155 for (Int_t ii=0; ii<tzero.fGTracks->GetEntries(); ii++)
156 fGTracks->AddLast(tzero.fGTracks->At(ii));
157
158 for (Int_t ii=0; ii<tzero.fTracksT0->GetEntries(); ii++)
159 fTracksT0->AddLast(tzero.fTracksT0->At(ii));
160
161 fOptFlag=tzero.fOptFlag;
162
163 return *this;
164}
165
166//____________________________________________________________________________
167AliTOFT0v1::~AliTOFT0v1()
168{
169 // dtor
170 fEvent=NULL;
171
172 if (fTracks) {
173 fTracks->Clear();
174 delete fTracks;
175 fTracks=0x0;
176 }
177
178 if (fGTracks) {
179 fGTracks->Clear();
180 delete fGTracks;
181 fGTracks=0x0;
182 }
183
184 if (fTracksT0) {
185 fTracksT0->Clear();
186 delete fTracksT0;
187 fTracksT0=0x0;
188 }
189
190
191}
192//____________________________________________________________________________
193
194void
195AliTOFT0v1::Init(AliESDEvent *event)
196{
197
198 /*
199 * init
200 */
201
202 fEvent = event;
203 fT0SigmaT0def[0]=0.;
204 fT0SigmaT0def[1]=0.6;
205 fT0SigmaT0def[2]=0.;
206 fT0SigmaT0def[3]=0.;
207
208}
209//____________________________________________________________________________
210Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
211{
212 TBenchmark *bench=new TBenchmark();
213 bench->Start("t0computation");
214
215 // Caluclate the Event Time using the ESD TOF time
216
217 fT0SigmaT0def[0]=0.;
218 fT0SigmaT0def[1]=0.600;
219 fT0SigmaT0def[2]=0.;
220 fT0SigmaT0def[3]=0.;
221
222 const Int_t nmaxtracksinsetMax=10;
223 Int_t nmaxtracksinset=10;
224// if(strstr(option,"all")){
225// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
226// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
227// }
228
229
230 Int_t nsets=0;
231 Int_t nUsedTracks=0;
232 Int_t ngoodsetsSel= 0;
233 Float_t t0bestSel[300];
234 Float_t eT0bestSel[300];
235 Float_t chiSquarebestSel[300];
236 Float_t confLevelbestSel[300];
237 Float_t t0bestallSel=0.;
238 Float_t eT0bestallSel=0.;
239 Float_t sumWt0bestallSel=0.;
240 Float_t eMeanTzeroPi=0.;
241 Float_t meantzeropi=0.;
242 Float_t sumAllweightspi=0.;
243 Double_t t0def=-999;
244 Double_t deltat0def=999;
245 Int_t ngoodtrktrulyused=0;
246 Int_t ntracksinsetmyCut = 0;
247
248 Int_t ntrk=fEvent->GetNumberOfTracks();
249
250 Int_t ngoodtrk=0;
251 Int_t ngoodtrkt0 =0;
252 Float_t meantime =0;
253
254 // First Track loop, Selection of good tracks
255
256 fTracks->Clear();
257 for (Int_t itrk=0; itrk<ntrk; itrk++) {
258 AliESDtrack *t=fEvent->GetTrack(itrk);
259 Double_t momOld=t->GetP();
260 Double_t mom=momOld-0.0036*momOld;
261 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
262 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
263 Double_t time=t->GetTOFsignal();
264
265 time*=1.E-3; // tof given in nanoseconds
266 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
267
268 if (!AcceptTrack(t)) continue;
269
270 if(t->GetIntegratedLength() < 350)continue; //skip decays
271 if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
272
273 meantime+=time;
274 fTracks->AddLast(t);
275 ngoodtrk++;
276 }
277
278 if(ngoodtrk > 1) meantime /= ngoodtrk;
279
280 if(ngoodtrk>22) nmaxtracksinset = 6;
281
282 fGTracks->Clear();
283 for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
284 AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
285 // Double_t time=t->GetTOFsignal();
286 // if((time-meantime*1.E3)<50.E3){ // For pp and per
287 fGTracks->AddLast(t);
288 ngoodtrkt0++;
289 // }
290 }
291
292 fTracks->Clear();
293
294 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
295 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
296 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
297
298
299 if(ngoodtrkt0<2){
300 t0def=-999;
301 deltat0def=0.600;
302 fT0SigmaT0def[0]=t0def;
303 fT0SigmaT0def[1]=deltat0def;
304 fT0SigmaT0def[2]=ngoodtrkt0;
305 fT0SigmaT0def[3]=ngoodtrkt0;
306 //goto finish;
307 }
308 if(ngoodtrkt0>=2){
309 // Decide how many tracks in set
310 Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
311 Int_t nset=1;
312
313 if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
314
315 // Loop over selected sets
316
317 if(nset>=1){
318 for (Int_t i=0; i< nset; i++) {
319 // printf("Set %i of %i\n",i+1,nset);
320 Float_t t0best=999.;
321 Float_t eT0best=999.;
322 Float_t chisquarebest=99999.;
323 Int_t npionbest=0;
324
325 fTracksT0->Clear();
326 Int_t ntracksinsetmy=0;
327 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
328 Int_t index = itrk+i*ntracksinset;
329 if(index < fGTracks->GetEntries()){
330 AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
331 fTracksT0->AddLast(t);
332 ntracksinsetmy++;
333 }
334 }
335
336 // Analyse it
337
338 Int_t assparticle[nmaxtracksinsetMax];
339 Float_t exptof[nmaxtracksinsetMax][3];
340 Float_t momErr[nmaxtracksinsetMax][3];
341 Float_t timeofflight[nmaxtracksinsetMax];
342 Float_t momentum[nmaxtracksinsetMax];
343 Float_t timezero[nmaxtracksinsetMax];
344 Float_t weightedtimezero[nmaxtracksinsetMax];
345 Float_t beta[nmaxtracksinsetMax];
346 Float_t texp[nmaxtracksinsetMax];
347 Float_t dtexp[nmaxtracksinsetMax];
348 Float_t sqMomError[nmaxtracksinsetMax];
349 Float_t sqTrackError[nmaxtracksinsetMax];
350 Float_t massarray[3]={0.13957,0.493677,0.9382723};
351 Float_t tracktoflen[nmaxtracksinsetMax];
352 Float_t besttimezero[nmaxtracksinsetMax];
353 Float_t besttexp[nmaxtracksinsetMax];
354 Float_t besttimeofflight[nmaxtracksinsetMax];
355 Float_t bestmomentum[nmaxtracksinsetMax];
356 Float_t bestchisquare[nmaxtracksinsetMax];
357 Float_t bestweightedtimezero[nmaxtracksinsetMax];
358 Float_t bestsqTrackError[nmaxtracksinsetMax];
359 Int_t imass[nmaxtracksinsetMax];
360
361 for (Int_t j=0; j<ntracksinset; j++) {
362 assparticle[j] = 3;
363 timeofflight[j] = 0;
364 momentum[j] = 0;
365 timezero[j] = 0;
366 weightedtimezero[j] = 0;
367 beta[j] = 0;
368 texp[j] = 0;
369 dtexp[j] = 0;
370 sqMomError[j] = 0;
371 sqTrackError[j] = 0;
372 tracktoflen[j] = 0;
373 besttimezero[j] = 0;
374 besttexp[j] = 0;
375 besttimeofflight[j] = 0;
376 bestmomentum[j] = 0;
377 bestchisquare[j] = 0;
378 bestweightedtimezero[j] = 0;
379 bestsqTrackError[j] = 0;
380 imass[j] = 1;
381 }
382
383 for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
384 AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
385 Double_t momOld=t->GetP();
386 Double_t mom=momOld-0.0036*momOld;
387 Double_t time=t->GetTOFsignal();
388
389 time*=1.E-3; // tof given in nanoseconds
390 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
391 Double_t toflen=t->GetIntegratedLength();
392 toflen=toflen/100.; // toflen given in m
393
394 timeofflight[j]=time;
395 tracktoflen[j]=toflen;
396 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
397 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
398 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
399 momentum[j]=mom;
400 assparticle[j]=3;
401
402 // in principle GetMomError only depends on two indices k=imass[j] and j itslef (see blow in the ncombinatorial loop)
403 // so it should be possible to make a lookup in order to speed up the code:
404 if (fOptFlag) {
405 momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]);
406 momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]);
407 momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]);
408 }
409 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
410
411 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
412 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
413 +momentum[itz]*momentum[itz]);
414 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]));
415 sqTrackError[itz]=sqMomError[itz]; //in ns
416 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
417 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
418 sumAllweightspi+=1./sqTrackError[itz];
419 meantzeropi+=weightedtimezero[itz];
420 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
421
422
423 // Then, Combinatorial Algorithm
424
425 if(ntracksinsetmy<2 )break;
426
427 for (Int_t j=0; j<ntracksinsetmy; j++) {
428 imass[j] = 3;
429 }
430
431 Int_t ncombinatorial;
432 if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
433 else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
434
435
436 // Loop on mass hypotheses
437 for (Int_t k=0; k < ncombinatorial;k++) {
438 for (Int_t j=0; j<ntracksinsetmy; j++) {
439 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
440 texp[j]=exptof[j][imass[j]];
441 if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
442 else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
443 }
444
445 Float_t sumAllweights=0.;
446 Float_t meantzero=0.;
447 Float_t eMeanTzero=0.;
448
449 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
450 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
451
452 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
453
454 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
455 sumAllweights+=1./sqTrackError[itz];
456 meantzero+=weightedtimezero[itz];
457
458 } // end loop for (Int_t itz=0; itz<15;itz++)
459
460 meantzero=meantzero/sumAllweights; // it is given in [ns]
461 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
462
463 // calculate chisquare
464 Float_t chisquare=0.;
465 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
466 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
467 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
468
469 if(chisquare<=chisquarebest){
470 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
471
472 bestsqTrackError[iqsq]=sqTrackError[iqsq];
473 besttimezero[iqsq]=timezero[iqsq];
474 bestmomentum[iqsq]=momentum[iqsq];
475 besttimeofflight[iqsq]=timeofflight[iqsq];
476 besttexp[iqsq]=texp[iqsq];
477 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
478 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
479 }
480
481 Int_t npion=0;
482 for (Int_t j=0; j<ntracksinsetmy; j++) {
483 assparticle[j]=imass[j];
484 if(imass[j] == 0) npion++;
485 }
486 npionbest=npion;
487 chisquarebest=chisquare;
488 t0best=meantzero;
489 eT0best=eMeanTzero;
490 } // close if(dummychisquare<=chisquare)
491 }
492
493 Double_t chi2cut[nmaxtracksinsetMax];
494 chi2cut[0] = 0;
495 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
496 for (Int_t j=2; j<ntracksinset; j++) {
497 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
498 }
499
500 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
501
502 // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
503
504 Bool_t kRedoT0 = kFALSE;
505 ntracksinsetmyCut = ntracksinsetmy;
506 Bool_t usetrack[nmaxtracksinsetMax];
507 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
508 usetrack[icsq] = kTRUE;
509 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
510 kRedoT0 = kTRUE;
511 ntracksinsetmyCut--;
512 usetrack[icsq] = kFALSE;
513 // printf("tracks chi2 = %f\n",bestchisquare[icsq]);
514 }
515 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
516
517 // Loop on mass hypotheses Redo
518 if(kRedoT0 && ntracksinsetmyCut > 1){
519 // printf("Redo T0\n");
520 for (Int_t k=0; k < ncombinatorial;k++) {
521 for (Int_t j=0; j<ntracksinsetmy; j++) {
522 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
523 texp[j]=exptof[j][imass[j]];
524 if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
525 else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
526 }
527
528 Float_t sumAllweights=0.;
529 Float_t meantzero=0.;
530 Float_t eMeanTzero=0.;
531
532 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
533 if(! usetrack[itz]) continue;
534 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
535
536 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
537
538 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
539 sumAllweights+=1./sqTrackError[itz];
540 meantzero+=weightedtimezero[itz];
541
542 } // end loop for (Int_t itz=0; itz<15;itz++)
543
544 meantzero=meantzero/sumAllweights; // it is given in [ns]
545 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
546
547 // calculate chisquare
548
549 Float_t chisquare=0.;
550 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
551 if(! usetrack[icsq]) continue;
552 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
553 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
554
555 Int_t npion=0;
556 for (Int_t j=0; j<ntracksinsetmy; j++) {
557 assparticle[j]=imass[j];
558 if(imass[j] == 0) npion++;
559 }
560
561 if(chisquare<=chisquarebest && npion>0){
562 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
563 if(! usetrack[iqsq]) continue;
564 bestsqTrackError[iqsq]=sqTrackError[iqsq];
565 besttimezero[iqsq]=timezero[iqsq];
566 bestmomentum[iqsq]=momentum[iqsq];
567 besttimeofflight[iqsq]=timeofflight[iqsq];
568 besttexp[iqsq]=texp[iqsq];
569 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
570 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
571 }
572
573 npionbest=npion;
574 chisquarebest=chisquare;
575 t0best=meantzero;
576 eT0best=eMeanTzero;
577 } // close if(dummychisquare<=chisquare)
578
579 }
580 }
581
582 // filling histos
583 Float_t confLevel=999;
584
585 // Sets with decent chisquares
586 // printf("Chi2best of the set = %f \n",chisquarebest);
587
588 if(chisquarebest<999.){
589 Double_t dblechisquare=(Double_t)chisquarebest;
590 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
591
592 Int_t ntrackincurrentsel=0;
593 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
594
595 if(! usetrack[icsq]) continue;
596
597 ntrackincurrentsel++;
598 }
599
600 // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
601
602 // Pick up only those with C.L. >1%
603 if(confLevel>0.01 && ngoodsetsSel<200){
604 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
605 confLevelbestSel[ngoodsetsSel]=confLevel;
606 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
607 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
608 t0bestallSel += t0best/eT0best/eT0best;
609 sumWt0bestallSel += 1./eT0best/eT0best;
610 ngoodsetsSel++;
611 ngoodtrktrulyused+=ntracksinsetmyCut;
612 // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
613 }
614 else{
615 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
616 }
617 }
618 fTracksT0->Clear();
619 nsets++;
620
621 } // end for the current set
622
623 //Redo the computation of the best in order to esclude very bad samples
624 if(ngoodsetsSel > 1){
625 Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
626 Int_t nsamples=ngoodsetsSel;
627 ngoodsetsSel=0;
628 t0bestallSel=0;
629 sumWt0bestallSel=0;
630 for (Int_t itz=0; itz<nsamples;itz++) {
631 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
632 t0bestallSel += t0bestSel[itz];
633 sumWt0bestallSel += eT0bestSel[itz];
634 ngoodsetsSel++;
635 // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
636 }
637 else{
638 // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
639 }
640 }
641 }
642 if(ngoodsetsSel < 1){
643 sumWt0bestallSel = 0.0;
644 }
645 //--------------------------------End recomputation
646
647 nUsedTracks = ngoodtrkt0;
648 if(strstr(option,"all")){
649 if(sumAllweightspi>0.){
650 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
651 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
652 }
653
654 // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
655
656 if(sumWt0bestallSel>0){
657 t0bestallSel = t0bestallSel/sumWt0bestallSel;
658 eT0bestallSel = sqrt(1./sumWt0bestallSel);
659 // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
660 }// end of if(sumWt0bestallSel>0){
661
662 }
663
664 t0def=t0bestallSel;
665 deltat0def=eT0bestallSel;
666
667 fT0SigmaT0def[0]=t0def;
668 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
669 fT0SigmaT0def[2]=ngoodtrkt0;
670 fT0SigmaT0def[3]=ngoodtrktrulyused;
671 }
672 }
673
674 fGTracks->Clear();
675
676 if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
677
678 bench->Stop("t0computation");
679
680 fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
681 fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
682
683// bench->Print("t0computation");
684// 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]));
685
686 delete bench;
687 bench=NULL;
688
689 return fT0SigmaT0def;
690 }
691//__________________________________________________________________
692Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
693{
694 // Take the error extimate for the TOF time in the track reconstruction
695
696 static const Double_t kMasses[]={
697 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
698 };
699
700 Double_t mass=kMasses[index+2];
701
702 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
703
704 return sigma;
705}
706
707//__________________________________________________________________
708Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
709{
710
711 /* TPC refit */
712 if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
713 /* do not accept kink daughters */
714 if (track->GetKinkIndex(0)>0) return kFALSE;
715 /* N clusters TPC */
716 if (track->GetTPCclusters(0) < 50) return kFALSE;
717 /* chi2 TPC */
718 if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
719 /* sigma to vertex */
720 if (GetSigmaToVertex(track) > 4.) return kFALSE;
721
722 /* accept track */
723 return kTRUE;
724
725}
726
727//____________________________________________________________________
728Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
729{
730 // Calculates the number of sigma to the vertex.
731
732 Float_t b[2];
733 Float_t bRes[2];
734 Float_t bCov[3];
735 esdTrack->GetImpactParameters(b,bCov);
736
737 if (bCov[0]<=0 || bCov[2]<=0) {
738 bCov[0]=0; bCov[2]=0;
739 }
740 bRes[0] = TMath::Sqrt(bCov[0]);
741 bRes[1] = TMath::Sqrt(bCov[2]);
742
743 // -----------------------------------
744 // How to get to a n-sigma cut?
745 //
746 // The accumulated statistics from 0 to d is
747 //
748 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
749 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
750 //
751 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
752 // Can this be expressed in a different way?
753
754 if (bRes[0] == 0 || bRes[1] ==0)
755 return -1;
756
757 //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
758 Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
759
760 // work around precision problem
761 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
762 // 1e-15 corresponds to nsigma ~ 7.7
763 if (TMath::Exp(-d * d / 2) < 1e-15)
764 return 1000;
765
766 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
767 return nSigma;
768}
769//____________________________________________________________________
770
771Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
772 Bool_t status = kFALSE;
773
774 Double_t exptimes[5];
775 track->GetIntegratedTimes(exptimes);
776
777 Float_t dedx = track->GetTPCsignal();
778
779 Double_t ptpc[3];
780 track->GetInnerPxPyPz(ptpc);
781 Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
782
783 if(imass > 2 || imass < 0) return status;
784 Int_t i = imass+2;
785
786 AliPID::EParticleType type=AliPID::EParticleType(i);
787
788 Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
789 Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
790
791 if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
792 status = kTRUE;
793 }
794
795 return status;
796}
797
798//____________________________________________________________________
799Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const
800{
801 //
802 // Returns base^exponent
803 //
804
805 Float_t power=1.;
806
807 for (Int_t ii=exponent; ii>0; ii--)
808 power=power*base;
809
810 return power;
811
812}
813//____________________________________________________________________
814Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const
815{
816 //
817 // Returns base^exponent
818 //
819
820 Int_t power=1;
821
822 for (Int_t ii=exponent; ii>0; ii--)
823 power=power*base;
824
825 return power;
826
827}