some coverity fixes
[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
a558aa69 56//-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table
57// for Power3)
536031f2 58//////////////////////////////////////////////////////////////////////////////
59
536031f2 60#include "AliESDtrack.h"
8f589502 61#include "AliESDEvent.h"
03bd764d 62#include "AliTOFT0v1.h"
2a258f40 63#include "TBenchmark.h"
64#include "AliPID.h"
65#include "AliESDpid.h"
536031f2 66
67ClassImp(AliTOFT0v1)
68
69//____________________________________________________________________________
2a258f40 70AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
5b4ed716 71 TObject(),
03bd764d 72 fLowerMomBound(0.5),
5b4ed716 73 fUpperMomBound(3),
8f589502 74 fTimeCorr(0.),
2a258f40 75 fEvent(0x0),
62f44f07 76 fPIDesd(extPID),
77 fTracks(new TObjArray(10)),
78 fGTracks(new TObjArray(10)),
a558aa69 79 fTracksT0(new TObjArray(10)),
80 fOptFlag(kFALSE)
536031f2 81{
8f589502 82 //
83 // default constructor
84 //
2a258f40 85 if(AliPID::ParticleMass(0) == 0) new AliPID();
86
87 if(!fPIDesd){
88 fPIDesd = new AliESDpid();
89 }
536031f2 90
5b4ed716 91 Init(NULL);
a558aa69 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 }
8f589502 99}
100
8f589502 101//____________________________________________________________________________
2a258f40 102AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
5b4ed716 103 TObject(),
03bd764d 104 fLowerMomBound(0.5),
5b4ed716 105 fUpperMomBound(3.0),
8f589502 106 fTimeCorr(0.),
2a258f40 107 fEvent(event),
62f44f07 108 fPIDesd(extPID),
109 fTracks(new TObjArray(10)),
110 fGTracks(new TObjArray(10)),
a558aa69 111 fTracksT0(new TObjArray(10)),
112 fOptFlag(kFALSE)
8f589502 113{
114 //
115 // real constructor
116 //
2a258f40 117 if(AliPID::ParticleMass(0) == 0) new AliPID();
536031f2 118
2a258f40 119 if(!fPIDesd){
120 fPIDesd = new AliESDpid();
121 }
5b4ed716 122
2a258f40 123 Init(event);
a558aa69 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 }
536031f2 128}
536031f2 129//____________________________________________________________________________
8f589502 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;
8f589502 141 fTimeCorr=tzero.fTimeCorr;
142 fEvent=tzero.fEvent;
8f589502 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
62f44f07 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));
a558aa69 160
161 fOptFlag=tzero.fOptFlag;
62f44f07 162
8f589502 163 return *this;
164}
5b4ed716 165
8f589502 166//____________________________________________________________________________
536031f2 167AliTOFT0v1::~AliTOFT0v1()
168{
169 // dtor
8f589502 170 fEvent=NULL;
e26296dc 171
172 if (fTracks) {
173 fTracks->Clear();
174 delete fTracks;
175 fTracks=0x0;
176 }
8f589502 177
e26296dc 178 if (fGTracks) {
179 fGTracks->Clear();
180 delete fGTracks;
181 fGTracks=0x0;
182 }
62f44f07 183
e26296dc 184 if (fTracksT0) {
185 fTracksT0->Clear();
186 delete fTracksT0;
187 fTracksT0=0x0;
188 }
62f44f07 189
a558aa69 190
536031f2 191}
8f589502 192//____________________________________________________________________________
5b4ed716 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}
721ed687 209//____________________________________________________________________________
5b4ed716 210Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
211{
2a258f40 212 TBenchmark *bench=new TBenchmark();
213 bench->Start("t0computation");
214
5b4ed716 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.;
5b4ed716 221
2a258f40 222 const Int_t nmaxtracksinsetMax=10;
223 Int_t nmaxtracksinset=10;
5b4ed716 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
5b4ed716 250 Int_t ngoodtrk=0;
251 Int_t ngoodtrkt0 =0;
721ed687 252 Float_t meantime =0;
5b4ed716 253
254 // First Track loop, Selection of good tracks
255
995c3398 256 fTracks->Clear();
5b4ed716 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;
2a258f40 267
5b4ed716 268 if (!AcceptTrack(t)) continue;
269
9c89b8bd 270 if(t->GetIntegratedLength() < 350)continue; //skip decays
5b4ed716 271 if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
721ed687 272
273 meantime+=time;
62f44f07 274 fTracks->AddLast(t);
5b4ed716 275 ngoodtrk++;
276 }
2a258f40 277
721ed687 278 if(ngoodtrk > 1) meantime /= ngoodtrk;
279
280 if(ngoodtrk>22) nmaxtracksinset = 6;
995c3398 281
282 fGTracks->Clear();
62f44f07 283 for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
62f44f07 284 AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
721ed687 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 // }
5b4ed716 290 }
62f44f07 291
62f44f07 292 fTracks->Clear();
293
5b4ed716 294 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
295 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
296 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
297
2a258f40 298
5b4ed716 299 if(ngoodtrkt0<2){
5b4ed716 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++) {
2a258f40 319 // printf("Set %i of %i\n",i+1,nset);
5b4ed716 320 Float_t t0best=999.;
321 Float_t eT0best=999.;
322 Float_t chisquarebest=99999.;
323 Int_t npionbest=0;
324
995c3398 325 fTracksT0->Clear();
5b4ed716 326 Int_t ntracksinsetmy=0;
5b4ed716 327 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
328 Int_t index = itrk+i*ntracksinset;
62f44f07 329 if(index < fGTracks->GetEntries()){
62f44f07 330 AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
62f44f07 331 fTracksT0->AddLast(t);
5b4ed716 332 ntracksinsetmy++;
333 }
334 }
62f44f07 335
5b4ed716 336 // Analyse it
337
2a258f40 338 Int_t assparticle[nmaxtracksinsetMax];
339 Float_t exptof[nmaxtracksinsetMax][3];
a558aa69 340 Float_t momErr[nmaxtracksinsetMax][3];
2a258f40 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];
5b4ed716 350 Float_t massarray[3]={0.13957,0.493677,0.9382723};
2a258f40 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];
5b4ed716 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
62f44f07 383 for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
62f44f07 384 AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
5b4ed716 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
fc9b31a7 390 Double_t exptime[AliPID::kSPECIESC];
391 t->GetIntegratedTimes(exptime,AliPID::kSPECIESC);
5b4ed716 392 Double_t toflen=t->GetIntegratedLength();
393 toflen=toflen/100.; // toflen given in m
394
395 timeofflight[j]=time;
396 tracktoflen[j]=toflen;
397 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
398 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
399 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
400 momentum[j]=mom;
401 assparticle[j]=3;
a558aa69 402
403 // in principle GetMomError only depends on two indices k=imass[j] and j itslef (see blow in the ncombinatorial loop)
404 // so it should be possible to make a lookup in order to speed up the code:
405 if (fOptFlag) {
406 momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]);
407 momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]);
408 momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]);
409 }
410 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
5b4ed716 411
412 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
413 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
414 +momentum[itz]*momentum[itz]);
415 sqMomError[itz]= ((1.-beta[itz]*beta[itz])*0.01)*((1.-beta[itz]*beta[itz])*0.01)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz]));
2a258f40 416 sqTrackError[itz]=sqMomError[itz]; //in ns
5b4ed716 417 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
418 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
419 sumAllweightspi+=1./sqTrackError[itz];
420 meantzeropi+=weightedtimezero[itz];
421 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
422
423
424 // Then, Combinatorial Algorithm
425
426 if(ntracksinsetmy<2 )break;
427
428 for (Int_t j=0; j<ntracksinsetmy; j++) {
429 imass[j] = 3;
430 }
a558aa69 431
432 Int_t ncombinatorial;
433 if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
434 else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
435
436
5b4ed716 437 // Loop on mass hypotheses
438 for (Int_t k=0; k < ncombinatorial;k++) {
439 for (Int_t j=0; j<ntracksinsetmy; j++) {
a558aa69 440 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
5b4ed716 441 texp[j]=exptof[j][imass[j]];
a558aa69 442 if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
443 else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
5b4ed716 444 }
90beec49 445
5b4ed716 446 Float_t sumAllweights=0.;
447 Float_t meantzero=0.;
448 Float_t eMeanTzero=0.;
449
450 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
2a258f40 451 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 452
453 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
454
455 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
456 sumAllweights+=1./sqTrackError[itz];
457 meantzero+=weightedtimezero[itz];
458
459 } // end loop for (Int_t itz=0; itz<15;itz++)
460
461 meantzero=meantzero/sumAllweights; // it is given in [ns]
462 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
463
464 // calculate chisquare
5b4ed716 465 Float_t chisquare=0.;
466 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1fde62a1 467 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
5b4ed716 468 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
469
470 if(chisquare<=chisquarebest){
471 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
472
473 bestsqTrackError[iqsq]=sqTrackError[iqsq];
474 besttimezero[iqsq]=timezero[iqsq];
475 bestmomentum[iqsq]=momentum[iqsq];
476 besttimeofflight[iqsq]=timeofflight[iqsq];
477 besttexp[iqsq]=texp[iqsq];
478 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 479 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
5b4ed716 480 }
481
482 Int_t npion=0;
483 for (Int_t j=0; j<ntracksinsetmy; j++) {
484 assparticle[j]=imass[j];
485 if(imass[j] == 0) npion++;
486 }
487 npionbest=npion;
488 chisquarebest=chisquare;
489 t0best=meantzero;
490 eT0best=eMeanTzero;
491 } // close if(dummychisquare<=chisquare)
5b4ed716 492 }
493
2a258f40 494 Double_t chi2cut[nmaxtracksinsetMax];
5b4ed716 495 chi2cut[0] = 0;
496 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
497 for (Int_t j=2; j<ntracksinset; j++) {
498 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
499 }
500
501 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
502
721ed687 503 // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
5b4ed716 504
505 Bool_t kRedoT0 = kFALSE;
506 ntracksinsetmyCut = ntracksinsetmy;
2a258f40 507 Bool_t usetrack[nmaxtracksinsetMax];
5b4ed716 508 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
509 usetrack[icsq] = kTRUE;
510 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
90beec49 511 kRedoT0 = kTRUE;
512 ntracksinsetmyCut--;
513 usetrack[icsq] = kFALSE;
721ed687 514 // printf("tracks chi2 = %f\n",bestchisquare[icsq]);
5b4ed716 515 }
516 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
517
5b4ed716 518 // Loop on mass hypotheses Redo
519 if(kRedoT0 && ntracksinsetmyCut > 1){
520 // printf("Redo T0\n");
521 for (Int_t k=0; k < ncombinatorial;k++) {
522 for (Int_t j=0; j<ntracksinsetmy; j++) {
a558aa69 523 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
5b4ed716 524 texp[j]=exptof[j][imass[j]];
a558aa69 525 if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
526 else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
5b4ed716 527 }
528
529 Float_t sumAllweights=0.;
530 Float_t meantzero=0.;
531 Float_t eMeanTzero=0.;
532
533 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
534 if(! usetrack[itz]) continue;
2a258f40 535 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 536
537 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
538
539 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
540 sumAllweights+=1./sqTrackError[itz];
541 meantzero+=weightedtimezero[itz];
542
543 } // end loop for (Int_t itz=0; itz<15;itz++)
544
545 meantzero=meantzero/sumAllweights; // it is given in [ns]
546 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
547
548 // calculate chisquare
549
550 Float_t chisquare=0.;
551 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
552 if(! usetrack[icsq]) continue;
1fde62a1 553 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
5b4ed716 554 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
555
556 Int_t npion=0;
557 for (Int_t j=0; j<ntracksinsetmy; j++) {
558 assparticle[j]=imass[j];
559 if(imass[j] == 0) npion++;
560 }
561
90beec49 562 if(chisquare<=chisquarebest && npion>0){
5b4ed716 563 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
564 if(! usetrack[iqsq]) continue;
565 bestsqTrackError[iqsq]=sqTrackError[iqsq];
566 besttimezero[iqsq]=timezero[iqsq];
567 bestmomentum[iqsq]=momentum[iqsq];
568 besttimeofflight[iqsq]=timeofflight[iqsq];
569 besttexp[iqsq]=texp[iqsq];
570 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 571 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
5b4ed716 572 }
573
574 npionbest=npion;
575 chisquarebest=chisquare;
576 t0best=meantzero;
577 eT0best=eMeanTzero;
578 } // close if(dummychisquare<=chisquare)
579
580 }
581 }
582
583 // filling histos
584 Float_t confLevel=999;
585
586 // Sets with decent chisquares
2a258f40 587 // printf("Chi2best of the set = %f \n",chisquarebest);
5b4ed716 588
589 if(chisquarebest<999.){
590 Double_t dblechisquare=(Double_t)chisquarebest;
591 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
5b4ed716 592
90beec49 593 Int_t ntrackincurrentsel=0;
5b4ed716 594 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
595
596 if(! usetrack[icsq]) continue;
597
90beec49 598 ntrackincurrentsel++;
5b4ed716 599 }
600
2a258f40 601 // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
602
5b4ed716 603 // Pick up only those with C.L. >1%
5b4ed716 604 if(confLevel>0.01 && ngoodsetsSel<200){
605 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
606 confLevelbestSel[ngoodsetsSel]=confLevel;
607 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
608 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
609 t0bestallSel += t0best/eT0best/eT0best;
610 sumWt0bestallSel += 1./eT0best/eT0best;
611 ngoodsetsSel++;
2a258f40 612 ngoodtrktrulyused+=ntracksinsetmyCut;
721ed687 613 // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
5b4ed716 614 }
615 else{
616 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
617 }
618 }
62f44f07 619 fTracksT0->Clear();
5b4ed716 620 nsets++;
621
622 } // end for the current set
623
90beec49 624 //Redo the computation of the best in order to esclude very bad samples
625 if(ngoodsetsSel > 1){
626 Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
627 Int_t nsamples=ngoodsetsSel;
628 ngoodsetsSel=0;
629 t0bestallSel=0;
630 sumWt0bestallSel=0;
631 for (Int_t itz=0; itz<nsamples;itz++) {
632 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
633 t0bestallSel += t0bestSel[itz];
634 sumWt0bestallSel += eT0bestSel[itz];
635 ngoodsetsSel++;
721ed687 636 // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
90beec49 637 }
638 else{
721ed687 639 // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
90beec49 640 }
641 }
642 }
643 if(ngoodsetsSel < 1){
644 sumWt0bestallSel = 0.0;
645 }
646 //--------------------------------End recomputation
647
5b4ed716 648 nUsedTracks = ngoodtrkt0;
649 if(strstr(option,"all")){
650 if(sumAllweightspi>0.){
651 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
652 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
653 }
654
2a258f40 655 // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
656
5b4ed716 657 if(sumWt0bestallSel>0){
658 t0bestallSel = t0bestallSel/sumWt0bestallSel;
659 eT0bestallSel = sqrt(1./sumWt0bestallSel);
2a258f40 660 // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
5b4ed716 661 }// end of if(sumWt0bestallSel>0){
662
5b4ed716 663 }
664
665 t0def=t0bestallSel;
666 deltat0def=eT0bestallSel;
5b4ed716 667
668 fT0SigmaT0def[0]=t0def;
2a258f40 669 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
5b4ed716 670 fT0SigmaT0def[2]=ngoodtrkt0;
671 fT0SigmaT0def[3]=ngoodtrktrulyused;
672 }
673 }
62f44f07 674
62f44f07 675 fGTracks->Clear();
676
2a258f40 677 if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
678
679 bench->Stop("t0computation");
680
681 fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
682 fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
683
684// bench->Print("t0computation");
685// 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]));
686
90beec49 687 delete bench;
688 bench=NULL;
5b4ed716 689
536031f2 690 return fT0SigmaT0def;
691 }
692//__________________________________________________________________
8f589502 693Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
536031f2 694{
8f589502 695 // Take the error extimate for the TOF time in the track reconstruction
536031f2 696
536031f2 697 static const Double_t kMasses[]={
698 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
699 };
700
701 Double_t mass=kMasses[index+2];
536031f2 702
2a258f40 703 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
536031f2 704
705 return sigma;
706}
14b2cbea 707
708//__________________________________________________________________
709Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
710{
711
712 /* TPC refit */
713 if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
714 /* do not accept kink daughters */
715 if (track->GetKinkIndex(0)>0) return kFALSE;
716 /* N clusters TPC */
717 if (track->GetTPCclusters(0) < 50) return kFALSE;
718 /* chi2 TPC */
719 if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
720 /* sigma to vertex */
721 if (GetSigmaToVertex(track) > 4.) return kFALSE;
722
723 /* accept track */
724 return kTRUE;
725
726}
727
728//____________________________________________________________________
8f589502 729Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
14b2cbea 730{
731 // Calculates the number of sigma to the vertex.
732
733 Float_t b[2];
734 Float_t bRes[2];
735 Float_t bCov[3];
736 esdTrack->GetImpactParameters(b,bCov);
737
738 if (bCov[0]<=0 || bCov[2]<=0) {
739 bCov[0]=0; bCov[2]=0;
740 }
741 bRes[0] = TMath::Sqrt(bCov[0]);
742 bRes[1] = TMath::Sqrt(bCov[2]);
743
744 // -----------------------------------
745 // How to get to a n-sigma cut?
746 //
747 // The accumulated statistics from 0 to d is
748 //
749 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
750 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
751 //
752 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
753 // Can this be expressed in a different way?
754
755 if (bRes[0] == 0 || bRes[1] ==0)
756 return -1;
757
62f44f07 758 //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
759 Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
14b2cbea 760
761 // work around precision problem
762 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
763 // 1e-15 corresponds to nsigma ~ 7.7
764 if (TMath::Exp(-d * d / 2) < 1e-15)
765 return 1000;
766
767 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
768 return nSigma;
769}
90beec49 770//____________________________________________________________________
771
772Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
773 Bool_t status = kFALSE;
774
d2dbb0fe 775 Double_t exptimes[AliPID::kSPECIESC];
fc9b31a7 776 track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC);
90beec49 777
778 Float_t dedx = track->GetTPCsignal();
779
780 Double_t ptpc[3];
781 track->GetInnerPxPyPz(ptpc);
782 Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
783
784 if(imass > 2 || imass < 0) return status;
785 Int_t i = imass+2;
786
787 AliPID::EParticleType type=AliPID::EParticleType(i);
788
789 Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
790 Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
791
792 if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
793 status = kTRUE;
794 }
795
796 return status;
797}
62f44f07 798
799//____________________________________________________________________
800Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const
801{
802 //
803 // Returns base^exponent
804 //
805
806 Float_t power=1.;
807
808 for (Int_t ii=exponent; ii>0; ii--)
809 power=power*base;
810
811 return power;
812
813}
814//____________________________________________________________________
815Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const
816{
817 //
818 // Returns base^exponent
819 //
820
821 Int_t power=1;
822
823 for (Int_t ii=exponent; ii>0; ii--)
824 power=power*base;
825
826 return power;
827
828}