#101318: Patch for various problems in AliROOT
[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
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;
a558aa69 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++) {
5b4ed716 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]));
2a258f40 415 sqTrackError[itz]=sqMomError[itz]; //in ns
5b4ed716 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 }
a558aa69 430
431 Int_t ncombinatorial;
432 if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
433 else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
434
435
5b4ed716 436 // Loop on mass hypotheses
437 for (Int_t k=0; k < ncombinatorial;k++) {
438 for (Int_t j=0; j<ntracksinsetmy; j++) {
a558aa69 439 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
5b4ed716 440 texp[j]=exptof[j][imass[j]];
a558aa69 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]);
5b4ed716 443 }
90beec49 444
5b4ed716 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++) {
2a258f40 450 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 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
5b4ed716 464 Float_t chisquare=0.;
465 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1fde62a1 466 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
5b4ed716 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];
1fde62a1 478 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
5b4ed716 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)
5b4ed716 491 }
492
2a258f40 493 Double_t chi2cut[nmaxtracksinsetMax];
5b4ed716 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
721ed687 502 // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
5b4ed716 503
504 Bool_t kRedoT0 = kFALSE;
505 ntracksinsetmyCut = ntracksinsetmy;
2a258f40 506 Bool_t usetrack[nmaxtracksinsetMax];
5b4ed716 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)){
90beec49 510 kRedoT0 = kTRUE;
511 ntracksinsetmyCut--;
512 usetrack[icsq] = kFALSE;
721ed687 513 // printf("tracks chi2 = %f\n",bestchisquare[icsq]);
5b4ed716 514 }
515 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
516
5b4ed716 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++) {
a558aa69 522 imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
5b4ed716 523 texp[j]=exptof[j][imass[j]];
a558aa69 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]);
5b4ed716 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;
2a258f40 534 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 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;
1fde62a1 552 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
5b4ed716 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
90beec49 561 if(chisquare<=chisquarebest && npion>0){
5b4ed716 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];
1fde62a1 570 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
5b4ed716 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
2a258f40 586 // printf("Chi2best of the set = %f \n",chisquarebest);
5b4ed716 587
588 if(chisquarebest<999.){
589 Double_t dblechisquare=(Double_t)chisquarebest;
590 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
5b4ed716 591
90beec49 592 Int_t ntrackincurrentsel=0;
5b4ed716 593 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
594
595 if(! usetrack[icsq]) continue;
596
90beec49 597 ntrackincurrentsel++;
5b4ed716 598 }
599
2a258f40 600 // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
601
5b4ed716 602 // Pick up only those with C.L. >1%
5b4ed716 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++;
2a258f40 611 ngoodtrktrulyused+=ntracksinsetmyCut;
721ed687 612 // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
5b4ed716 613 }
614 else{
615 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
616 }
617 }
62f44f07 618 fTracksT0->Clear();
5b4ed716 619 nsets++;
620
621 } // end for the current set
622
90beec49 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++;
721ed687 635 // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
90beec49 636 }
637 else{
721ed687 638 // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
90beec49 639 }
640 }
641 }
642 if(ngoodsetsSel < 1){
643 sumWt0bestallSel = 0.0;
644 }
645 //--------------------------------End recomputation
646
5b4ed716 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
2a258f40 654 // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
655
5b4ed716 656 if(sumWt0bestallSel>0){
657 t0bestallSel = t0bestallSel/sumWt0bestallSel;
658 eT0bestallSel = sqrt(1./sumWt0bestallSel);
2a258f40 659 // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
5b4ed716 660 }// end of if(sumWt0bestallSel>0){
661
5b4ed716 662 }
663
664 t0def=t0bestallSel;
665 deltat0def=eT0bestallSel;
5b4ed716 666
667 fT0SigmaT0def[0]=t0def;
2a258f40 668 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
5b4ed716 669 fT0SigmaT0def[2]=ngoodtrkt0;
670 fT0SigmaT0def[3]=ngoodtrktrulyused;
671 }
672 }
62f44f07 673
62f44f07 674 fGTracks->Clear();
675
2a258f40 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
90beec49 686 delete bench;
687 bench=NULL;
5b4ed716 688
536031f2 689 return fT0SigmaT0def;
690 }
691//__________________________________________________________________
8f589502 692Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
536031f2 693{
8f589502 694 // Take the error extimate for the TOF time in the track reconstruction
536031f2 695
536031f2 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];
536031f2 701
2a258f40 702 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
536031f2 703
704 return sigma;
705}
14b2cbea 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//____________________________________________________________________
8f589502 728Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
14b2cbea 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
62f44f07 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));
14b2cbea 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}
90beec49 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}
62f44f07 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}