]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/AliTOFT0v1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix Coverity defects
[u/mrichter/AliRoot.git] / TOF / AliTOFT0v1.cxx
CommitLineData
536031f2 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
6a4e212e 16/* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
536031f2 17
18//_________________________________________________________________________
19// This is a TTask that made the calculation of the Time zero using TOF.
20// Description: The algorithm used to calculate the time zero of interaction
21// using TOF detector is the following.
6a4e212e 22// We select in the ESD some "primary" particles - or tracks in the following -
536031f2 23// that strike the TOF detector (the larger part are pions, kaons or protons).
24// We choose a set of 10 selected tracks, for each track You have the length
6a4e212e 25// of the track when the TOF is reached,
26// the momentum and the time of flight
536031f2 27// given by the TOF detector.
536031f2 28// Let consider now only one set of 10 tracks (the algorithm is the same for all sets).
29// Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton,
30// we consider all the 3 at 10 possible cases.
31// For each track in each (mass) configuration
32// (a configuration can be e.g. pion/pion/kaon/proton/pion/proton/kaon/kaon/pion/pion)
33// we calculate the time zero (we know in fact the velocity of the track after
34// the assumption about its mass, the time of flight given by the TOF, and the
35// corresponding path travelled till the TOF detector). Then for each mass configuration we have
36// 10 time zero and we can calculate the ChiSquare for the current configuration using the
37// weighted mean over all 10 time zero.
38// We call the best assignment the mass configuration that gives the minimum value of the ChiSquare.
39// We plot the weighted mean over all 10 time zero for the best assignment,
40// the ChiSquare for the best assignment and the corresponding confidence level.
41// The strong assumption is the MC selection of primary particles. It will be introduced
42// in the future also some more realistic simulation about this point.
43// Use case:
44// root [0] AliTOFT0v1 * tzero = new AliTOFT0v1("galice.root")
45// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
46// root [1] tzero->ExecuteTask()
47// root [2] tzero->ExecuteTask("tim")
48// // available parameters:
49// tim - print benchmarking information
50// all - print usefull informations about the number of misidentified tracks
51// and a comparison about the true configuration (known from MC) and the best
52// assignment
53// Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions
54//-- Author: F. Pierella
6a4e212e 55//-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
536031f2 56//////////////////////////////////////////////////////////////////////////////
57
536031f2 58#include "AliESDtrack.h"
8f589502 59#include "AliESDEvent.h"
03bd764d 60#include "AliTOFT0v1.h"
2a258f40 61#include "TBenchmark.h"
62#include "AliPID.h"
63#include "AliESDpid.h"
536031f2 64
65ClassImp(AliTOFT0v1)
66
67//____________________________________________________________________________
2a258f40 68AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
5b4ed716 69 TObject(),
03bd764d 70 fLowerMomBound(0.5),
5b4ed716 71 fUpperMomBound(3),
8f589502 72 fTimeCorr(0.),
2a258f40 73 fEvent(0x0),
62f44f07 74 fPIDesd(extPID),
75 fTracks(new TObjArray(10)),
76 fGTracks(new TObjArray(10)),
77 fTracksT0(new TObjArray(10))
536031f2 78{
8f589502 79 //
80 // default constructor
81 //
2a258f40 82 if(AliPID::ParticleMass(0) == 0) new AliPID();
83
84 if(!fPIDesd){
85 fPIDesd = new AliESDpid();
86 }
536031f2 87
5b4ed716 88 Init(NULL);
89
8f589502 90}
91
8f589502 92//____________________________________________________________________________
2a258f40 93AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
5b4ed716 94 TObject(),
03bd764d 95 fLowerMomBound(0.5),
5b4ed716 96 fUpperMomBound(3.0),
8f589502 97 fTimeCorr(0.),
2a258f40 98 fEvent(event),
62f44f07 99 fPIDesd(extPID),
100 fTracks(new TObjArray(10)),
101 fGTracks(new TObjArray(10)),
102 fTracksT0(new TObjArray(10))
8f589502 103{
104 //
105 // real constructor
106 //
2a258f40 107 if(AliPID::ParticleMass(0) == 0) new AliPID();
536031f2 108
2a258f40 109 if(!fPIDesd){
110 fPIDesd = new AliESDpid();
111 }
5b4ed716 112
2a258f40 113 Init(event);
8f589502 114
536031f2 115}
8f589502 116//____________________________________________________________________________
117AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
118{
119 //
120 // assign. operator
121 //
122
123 if (this == &tzero)
124 return *this;
125
126 fLowerMomBound=tzero.fLowerMomBound;
127 fUpperMomBound=tzero.fUpperMomBound;
8f589502 128 fTimeCorr=tzero.fTimeCorr;
129 fEvent=tzero.fEvent;
8f589502 130 fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
131 fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
132 fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
133 fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
134
62f44f07 135 fTracks=tzero.fTracks;
136 fGTracks=tzero.fGTracks;
137 fTracksT0=tzero.fTracksT0;
138
139 for (Int_t ii=0; ii<tzero.fTracks->GetEntries(); ii++)
140 fTracks->AddLast(tzero.fTracks->At(ii));
141
142 for (Int_t ii=0; ii<tzero.fGTracks->GetEntries(); ii++)
143 fGTracks->AddLast(tzero.fGTracks->At(ii));
144
145 for (Int_t ii=0; ii<tzero.fTracksT0->GetEntries(); ii++)
146 fTracksT0->AddLast(tzero.fTracksT0->At(ii));
147
8f589502 148 return *this;
149}
5b4ed716 150
536031f2 151//____________________________________________________________________________
152AliTOFT0v1::~AliTOFT0v1()
153{
154 // dtor
8f589502 155 fEvent=NULL;
156
62f44f07 157 fTracks->Delete();
158 delete fTracks;
159 fTracks=0x0;
160
161 fGTracks->Delete();
162 delete fGTracks;
163 fGTracks=0x0;
164
165 fTracksT0->Delete();
166 delete fTracksT0;
167 fTracksT0=0x0;
168
536031f2 169}
5b4ed716 170//____________________________________________________________________________
171
172void
173AliTOFT0v1::Init(AliESDEvent *event)
174{
175
176 /*
177 * init
178 */
179
180 fEvent = event;
181 fT0SigmaT0def[0]=0.;
182 fT0SigmaT0def[1]=0.6;
183 fT0SigmaT0def[2]=0.;
184 fT0SigmaT0def[3]=0.;
185
186}
187
536031f2 188//____________________________________________________________________________
189Double_t * AliTOFT0v1::DefineT0(Option_t *option)
190{
8f589502 191 // Caluclate the Event Time using the ESD TOF time
192
2a258f40 193 TBenchmark *bench=new TBenchmark();
194 bench->Start("t0computation");
195
5b4ed716 196 fT0SigmaT0def[0]=0.;
197 fT0SigmaT0def[1]=0.600;
198 fT0SigmaT0def[2]=0.;
199 fT0SigmaT0def[3]=0.;
536031f2 200
2a258f40 201 const Int_t nmaxtracksinsetMax=10;
202 Int_t nmaxtracksinset=10;
6a4e212e 203// if(strstr(option,"all")){
204// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
205// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
206// }
536031f2 207
536031f2 208 Int_t nsets=0;
8f589502 209 Int_t nUsedTracks=0;
536031f2 210 Int_t ngoodsetsSel= 0;
211 Float_t t0bestSel[300];
8f589502 212 Float_t eT0bestSel[300];
213 Float_t chiSquarebestSel[300];
214 Float_t confLevelbestSel[300];
536031f2 215 Float_t t0bestallSel=0.;
8f589502 216 Float_t eT0bestallSel=0.;
536031f2 217 Float_t sumWt0bestallSel=0.;
8f589502 218 Float_t eMeanTzeroPi=0.;
536031f2 219 Float_t meantzeropi=0.;
220 Float_t sumAllweightspi=0.;
221 Double_t t0def=-999;
222 Double_t deltat0def=999;
223 Int_t ngoodtrktrulyused=0;
224 Int_t ntracksinsetmyCut = 0;
225
226 Int_t ntrk=fEvent->GetNumberOfTracks();
227
62f44f07 228 //AliESDtrack **fTracks=new AliESDtrack*[ntrk];
536031f2 229 Int_t ngoodtrk=0;
230 Int_t ngoodtrkt0 =0;
231 Float_t mintime =1E6;
232
233 // First Track loop, Selection of good tracks
234
235 for (Int_t itrk=0; itrk<ntrk; itrk++) {
236 AliESDtrack *t=fEvent->GetTrack(itrk);
237 Double_t momOld=t->GetP();
238 Double_t mom=momOld-0.0036*momOld;
239 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
03bd764d 240 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
536031f2 241 Double_t time=t->GetTOFsignal();
242
243 time*=1.E-3; // tof given in nanoseconds
244 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
2a258f40 245
14b2cbea 246 if (!AcceptTrack(t)) continue;
536031f2 247
2a258f40 248 if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
90beec49 249
536031f2 250 if(time <= mintime) mintime=time;
62f44f07 251 //fTracks[ngoodtrk]=t;
252 fTracks->AddLast(t);
536031f2 253 ngoodtrk++;
254 }
62f44f07 255
2a258f40 256 if(ngoodtrk>22) nmaxtracksinset = 6;
257
03bd764d 258// cout << " N. of ESD tracks : " << ntrk << endl;
259// cout << " N. of preselected tracks : " << ngoodtrk << endl;
260// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
536031f2 261
62f44f07 262 //AliESDtrack **fGTracks=new AliESDtrack*[ngoodtrk];
536031f2 263
62f44f07 264 //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
265 for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
266 //AliESDtrack *t=fTracks[jtrk];
267 AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
536031f2 268 Double_t time=t->GetTOFsignal();
269
270 if((time-mintime*1.E3)<50.E3){ // For pp and per
62f44f07 271 //fGTracks[ngoodtrkt0]=t;
272 fGTracks->AddLast(t);
536031f2 273 ngoodtrkt0++;
274 }
275 }
276
62f44f07 277 //delete[] fTracks;
278 fTracks->Clear();
536031f2 279
280 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
281 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
282 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
283
284 if(ngoodtrkt0<2){
6a4e212e 285// cout << "less than 2 tracks, skip event " << endl;
536031f2 286 t0def=-999;
287 deltat0def=0.600;
288 fT0SigmaT0def[0]=t0def;
289 fT0SigmaT0def[1]=deltat0def;
290 fT0SigmaT0def[2]=ngoodtrkt0;
291 fT0SigmaT0def[3]=ngoodtrkt0;
292 //goto finish;
293 }
294 if(ngoodtrkt0>=2){
295 // Decide how many tracks in set
03bd764d 296 Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
536031f2 297 Int_t nset=1;
298
299 if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
6a4e212e 300
536031f2 301 // Loop over selected sets
302
303 if(nset>=1){
304 for (Int_t i=0; i< nset; i++) {
305
306 Float_t t0best=999.;
8f589502 307 Float_t eT0best=999.;
536031f2 308 Float_t chisquarebest=99999.;
309 Int_t npionbest=0;
310
311 Int_t ntracksinsetmy=0;
62f44f07 312 //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset];
536031f2 313 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
314 Int_t index = itrk+i*ntracksinset;
62f44f07 315 //if(index < ngoodtrkt0){
316 if(index < fGTracks->GetEntries()){
317 //AliESDtrack *t=fGTracks[index];
318 AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
319 //fTracksT0[itrk]=t;
320 fTracksT0->AddLast(t);
536031f2 321 ntracksinsetmy++;
322 }
323 }
324
325 // Analyse it
326
2a258f40 327 Int_t assparticle[nmaxtracksinsetMax];
328 Float_t exptof[nmaxtracksinsetMax][3];
329 Float_t timeofflight[nmaxtracksinsetMax];
330 Float_t momentum[nmaxtracksinsetMax];
331 Float_t timezero[nmaxtracksinsetMax];
332 Float_t weightedtimezero[nmaxtracksinsetMax];
333 Float_t beta[nmaxtracksinsetMax];
334 Float_t texp[nmaxtracksinsetMax];
335 Float_t dtexp[nmaxtracksinsetMax];
336 Float_t sqMomError[nmaxtracksinsetMax];
337 Float_t sqTrackError[nmaxtracksinsetMax];
536031f2 338 Float_t massarray[3]={0.13957,0.493677,0.9382723};
2a258f40 339 Float_t tracktoflen[nmaxtracksinsetMax];
340 Float_t besttimezero[nmaxtracksinsetMax];
341 Float_t besttexp[nmaxtracksinsetMax];
342 Float_t besttimeofflight[nmaxtracksinsetMax];
343 Float_t bestmomentum[nmaxtracksinsetMax];
344 Float_t bestchisquare[nmaxtracksinsetMax];
345 Float_t bestweightedtimezero[nmaxtracksinsetMax];
346 Float_t bestsqTrackError[nmaxtracksinsetMax];
347 Int_t imass[nmaxtracksinsetMax];
536031f2 348
349 for (Int_t j=0; j<ntracksinset; j++) {
350 assparticle[j] = 3;
351 timeofflight[j] = 0;
352 momentum[j] = 0;
353 timezero[j] = 0;
354 weightedtimezero[j] = 0;
355 beta[j] = 0;
356 texp[j] = 0;
357 dtexp[j] = 0;
358 sqMomError[j] = 0;
359 sqTrackError[j] = 0;
360 tracktoflen[j] = 0;
361 besttimezero[j] = 0;
362 besttexp[j] = 0;
363 besttimeofflight[j] = 0;
364 bestmomentum[j] = 0;
365 bestchisquare[j] = 0;
366 bestweightedtimezero[j] = 0;
367 bestsqTrackError[j] = 0;
368 imass[j] = 1;
369 }
370
62f44f07 371 //for (Int_t j=0; j<ntracksinsetmy; j++) {
372 for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
373 //AliESDtrack *t=fTracksT0[j];
374 AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
536031f2 375 Double_t momOld=t->GetP();
376 Double_t mom=momOld-0.0036*momOld;
377 Double_t time=t->GetTOFsignal();
378
379 time*=1.E-3; // tof given in nanoseconds
380 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
381 Double_t toflen=t->GetIntegratedLength();
382 toflen=toflen/100.; // toflen given in m
383
8f589502 384 timeofflight[j]=time;
385 tracktoflen[j]=toflen;
536031f2 386 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
387 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
388 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
389 momentum[j]=mom;
390 assparticle[j]=3;
391
392 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
393
536031f2 394 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
395 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
396 +momentum[itz]*momentum[itz]);
397 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 398 sqTrackError[itz]=sqMomError[itz]; //in ns
536031f2 399 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
400 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
401 sumAllweightspi+=1./sqTrackError[itz];
402 meantzeropi+=weightedtimezero[itz];
403 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
404
405
406 // Then, Combinatorial Algorithm
407
408 if(ntracksinsetmy<2 )break;
409
410 for (Int_t j=0; j<ntracksinsetmy; j++) {
411 imass[j] = 3;
412 }
413
62f44f07 414 //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
415 Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
536031f2 416
417 // Loop on mass hypotheses
418 for (Int_t k=0; k < ncombinatorial;k++) {
419 for (Int_t j=0; j<ntracksinsetmy; j++) {
62f44f07 420 //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
421 imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
536031f2 422 texp[j]=exptof[j][imass[j]];
423 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
1fde62a1 424 // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
536031f2 425 }
90beec49 426
536031f2 427 Float_t sumAllweights=0.;
428 Float_t meantzero=0.;
8f589502 429 Float_t eMeanTzero=0.;
536031f2 430
431 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
2a258f40 432 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
536031f2 433
434 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
435
436 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
437 sumAllweights+=1./sqTrackError[itz];
438 meantzero+=weightedtimezero[itz];
439
440 } // end loop for (Int_t itz=0; itz<15;itz++)
441
442 meantzero=meantzero/sumAllweights; // it is given in [ns]
8f589502 443 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
536031f2 444
445 // calculate chisquare
446
447 Float_t chisquare=0.;
448 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1fde62a1 449 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
450 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
451 // else chisquare+=1000;
536031f2 452 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
453
454 if(chisquare<=chisquarebest){
455 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
456
457 bestsqTrackError[iqsq]=sqTrackError[iqsq];
458 besttimezero[iqsq]=timezero[iqsq];
459 bestmomentum[iqsq]=momentum[iqsq];
460 besttimeofflight[iqsq]=timeofflight[iqsq];
461 besttexp[iqsq]=texp[iqsq];
462 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 463 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
464 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
465 // else bestchisquare[iqsq]=1000;
536031f2 466 }
467
468 Int_t npion=0;
469 for (Int_t j=0; j<ntracksinsetmy; j++) {
470 assparticle[j]=imass[j];
471 if(imass[j] == 0) npion++;
472 }
473 npionbest=npion;
474 chisquarebest=chisquare;
475 t0best=meantzero;
8f589502 476 eT0best=eMeanTzero;
536031f2 477 } // close if(dummychisquare<=chisquare)
478
479 }
480
2a258f40 481 Double_t chi2cut[nmaxtracksinsetMax];
536031f2 482 chi2cut[0] = 0;
483 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
484 for (Int_t j=2; j<ntracksinset; j++) {
485 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
486 }
487
488 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
489
03bd764d 490// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
536031f2 491
492 Bool_t kRedoT0 = kFALSE;
493 ntracksinsetmyCut = ntracksinsetmy;
2a258f40 494 Bool_t usetrack[nmaxtracksinsetMax];
536031f2 495 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
496 usetrack[icsq] = kTRUE;
497 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
498 kRedoT0 = kTRUE;
499 ntracksinsetmyCut--;
500 usetrack[icsq] = kFALSE;
501 }
502 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
503
6a4e212e 504 // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
536031f2 505
506 // Loop on mass hypotheses Redo
507 if(kRedoT0 && ntracksinsetmyCut > 1){
6a4e212e 508 // printf("Redo T0\n");
536031f2 509 for (Int_t k=0; k < ncombinatorial;k++) {
510 for (Int_t j=0; j<ntracksinsetmy; j++) {
62f44f07 511 //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
512 imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
536031f2 513 texp[j]=exptof[j][imass[j]];
514 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
1fde62a1 515 // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
536031f2 516 }
517
518 Float_t sumAllweights=0.;
519 Float_t meantzero=0.;
8f589502 520 Float_t eMeanTzero=0.;
536031f2 521
522 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
523 if(! usetrack[itz]) continue;
2a258f40 524 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
536031f2 525
526 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
527
528 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
529 sumAllweights+=1./sqTrackError[itz];
530 meantzero+=weightedtimezero[itz];
531
532 } // end loop for (Int_t itz=0; itz<15;itz++)
533
534 meantzero=meantzero/sumAllweights; // it is given in [ns]
8f589502 535 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
536031f2 536
537 // calculate chisquare
538
539 Float_t chisquare=0.;
540 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
541 if(! usetrack[icsq]) continue;
1fde62a1 542 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
543 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
544 // else chisquare+=1000;
536031f2 545 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
546
547 Int_t npion=0;
548 for (Int_t j=0; j<ntracksinsetmy; j++) {
549 assparticle[j]=imass[j];
550 if(imass[j] == 0) npion++;
551 }
552
90beec49 553 if(chisquare<=chisquarebest && npion>0){
536031f2 554 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
555 if(! usetrack[iqsq]) continue;
556 bestsqTrackError[iqsq]=sqTrackError[iqsq];
557 besttimezero[iqsq]=timezero[iqsq];
558 bestmomentum[iqsq]=momentum[iqsq];
559 besttimeofflight[iqsq]=timeofflight[iqsq];
560 besttexp[iqsq]=texp[iqsq];
561 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 562 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
563 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
564 // else bestchisquare[iqsq]=1000;
536031f2 565 }
566
567 npionbest=npion;
568 chisquarebest=chisquare;
569 t0best=meantzero;
8f589502 570 eT0best=eMeanTzero;
536031f2 571 } // close if(dummychisquare<=chisquare)
572
573 }
574 }
6a4e212e 575
536031f2 576 // filling histos
577 Float_t confLevel=999;
578
579 // Sets with decent chisquares
580
581 if(chisquarebest<999.){
582 Double_t dblechisquare=(Double_t)chisquarebest;
583 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
03bd764d 584// cout << " Set Number " << nsets << endl;
585// cout << "Best Assignment, selection " << assparticle[0] <<
586// assparticle[1] << assparticle[2] <<
587// assparticle[3] << assparticle[4] <<
588// assparticle[5] << endl;
589// cout << " Chisquare of the set "<< chisquarebest <<endl;
590// cout << " C.L. of the set "<< confLevel <<endl;
591// cout << " T0 for this set (in ns) " << t0best << endl;
536031f2 592
593 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
594
595 if(! usetrack[icsq]) continue;
596
03bd764d 597// cout << "Track # " << icsq << " T0 offsets = "
598// << besttimezero[icsq]-t0best <<
599// " track error = " << bestsqTrackError[icsq]
600// << " Chisquare = " << bestchisquare[icsq]
601// << " Momentum = " << bestmomentum[icsq]
602// << " TOF = " << besttimeofflight[icsq]
603// << " TOF tracking = " << besttexp[icsq]
604// << " is used = " << usetrack[icsq] << endl;
536031f2 605 }
606
607 // Pick up only those with C.L. >1%
608 // if(confLevel>0.01 && ngoodsetsSel<200){
609 if(confLevel>0.01 && ngoodsetsSel<200){
8f589502 610 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
611 confLevelbestSel[ngoodsetsSel]=confLevel;
612 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
613 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
614 t0bestallSel += t0best/eT0best/eT0best;
615 sumWt0bestallSel += 1./eT0best/eT0best;
536031f2 616 ngoodsetsSel++;
617 ngoodtrktrulyused+=ntracksinsetmyCut;
2a258f40 618
619 // printf("t0 of a set = %f +/- %f\n",t0best,eT0best);
536031f2 620 }
621 else{
6a4e212e 622 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
536031f2 623 }
624 }
62f44f07 625 //delete[] fTracksT0;
626 fTracksT0->Clear();
536031f2 627 nsets++;
628
629 } // end for the current set
630
90beec49 631 //Redo the computation of the best in order to esclude very bad samples
632 if(ngoodsetsSel > 1){
633 Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
634 Int_t nsamples=ngoodsetsSel;
635 ngoodsetsSel=0;
636 t0bestallSel=0;
637 sumWt0bestallSel=0;
638 for (Int_t itz=0; itz<nsamples;itz++) {
639 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
640 t0bestallSel += t0bestSel[itz];
641 sumWt0bestallSel += eT0bestSel[itz];
642 ngoodsetsSel++;
643// printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
644 }
645 else{
646// printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
647 }
648 }
649 }
650 if(ngoodsetsSel < 1){
651 sumWt0bestallSel = 0.0;
652 }
653 //--------------------------------End recomputation
654
8f589502 655 nUsedTracks = ngoodtrkt0;
536031f2 656 if(strstr(option,"all")){
657 if(sumAllweightspi>0.){
658 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
8f589502 659 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
536031f2 660 }
661
662 if(sumWt0bestallSel>0){
663 t0bestallSel = t0bestallSel/sumWt0bestallSel;
8f589502 664 eT0bestallSel = sqrt(1./sumWt0bestallSel);
536031f2 665
666 }// end of if(sumWt0bestallSel>0){
667
6a4e212e 668// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
536031f2 669 }
670
671 t0def=t0bestallSel;
8f589502 672 deltat0def=eT0bestallSel;
536031f2 673
674 fT0SigmaT0def[0]=t0def;
2a258f40 675 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
6a4e212e 676 fT0SigmaT0def[2]=ngoodtrkt0;
536031f2 677 fT0SigmaT0def[3]=ngoodtrktrulyused;
678 }
679 }
62f44f07 680
681 //delete [] fGTracks;
682 fGTracks->Clear();
683
6a4e212e 684 // if(strstr(option,"tim") || strstr(option,"all")){
685 // cout << "AliTOFT0v1:" << endl ;
686 //}
536031f2 687
2a258f40 688 if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
689
690 bench->Stop("t0computation");
691
692 fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
693 fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
694
695// bench->Print("t0computation");
696// printf("%f %f\n",fT0SigmaT0def[4],fT0SigmaT0def[5]);
697
90beec49 698
699 delete bench;
700 bench=NULL;
5b4ed716 701
702 return fT0SigmaT0def;
703 }
704//__________________________________________________________________
705Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
706{
2a258f40 707 TBenchmark *bench=new TBenchmark();
708 bench->Start("t0computation");
709
5b4ed716 710 // Caluclate the Event Time using the ESD TOF time
711
712 fT0SigmaT0def[0]=0.;
713 fT0SigmaT0def[1]=0.600;
714 fT0SigmaT0def[2]=0.;
715 fT0SigmaT0def[3]=0.;
5b4ed716 716
2a258f40 717 const Int_t nmaxtracksinsetMax=10;
718 Int_t nmaxtracksinset=10;
5b4ed716 719// if(strstr(option,"all")){
720// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
721// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
722// }
723
724
725 Int_t nsets=0;
726 Int_t nUsedTracks=0;
727 Int_t ngoodsetsSel= 0;
728 Float_t t0bestSel[300];
729 Float_t eT0bestSel[300];
730 Float_t chiSquarebestSel[300];
731 Float_t confLevelbestSel[300];
732 Float_t t0bestallSel=0.;
733 Float_t eT0bestallSel=0.;
734 Float_t sumWt0bestallSel=0.;
735 Float_t eMeanTzeroPi=0.;
736 Float_t meantzeropi=0.;
737 Float_t sumAllweightspi=0.;
738 Double_t t0def=-999;
739 Double_t deltat0def=999;
740 Int_t ngoodtrktrulyused=0;
741 Int_t ntracksinsetmyCut = 0;
742
743 Int_t ntrk=fEvent->GetNumberOfTracks();
744
62f44f07 745 //AliESDtrack **fTracks=new AliESDtrack*[ntrk];
5b4ed716 746 Int_t ngoodtrk=0;
747 Int_t ngoodtrkt0 =0;
748 Float_t mintime =1E6;
749
750 // First Track loop, Selection of good tracks
751
752 for (Int_t itrk=0; itrk<ntrk; itrk++) {
753 AliESDtrack *t=fEvent->GetTrack(itrk);
754 Double_t momOld=t->GetP();
755 Double_t mom=momOld-0.0036*momOld;
756 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
757 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
758 Double_t time=t->GetTOFsignal();
759
760 time*=1.E-3; // tof given in nanoseconds
761 if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
2a258f40 762
5b4ed716 763 if (!AcceptTrack(t)) continue;
764
9c89b8bd 765 if(t->GetIntegratedLength() < 350)continue; //skip decays
5b4ed716 766 if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
767 if(time <= mintime) mintime=time;
62f44f07 768 //fTracks[ngoodtrk]=t;
769 fTracks->AddLast(t);
5b4ed716 770 ngoodtrk++;
771 }
772
2a258f40 773 if(ngoodtrk>22) nmaxtracksinset = 6;
774
5b4ed716 775
2a258f40 776// cout << " N. of ESD tracks : " << ntrk << endl;
777// cout << " N. of preselected tracks : " << ngoodtrk << endl;
778// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
5b4ed716 779
62f44f07 780 //AliESDtrack **fGTracks=new AliESDtrack*[ngoodtrk];
5b4ed716 781
62f44f07 782 //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
783 for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
784 //AliESDtrack *t=fTracks[jtrk];
785 AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
5b4ed716 786 Double_t time=t->GetTOFsignal();
787
788 if((time-mintime*1.E3)<50.E3){ // For pp and per
62f44f07 789 //fGTracks[ngoodtrkt0]=t;
790 fGTracks->AddLast(t);
5b4ed716 791 ngoodtrkt0++;
792 }
793 }
62f44f07 794
795 //delete [] fTracks;
796 fTracks->Clear();
797
2a258f40 798// cout << " N. of preselected tracks t0 : " << ngoodtrkt0 << endl;
5b4ed716 799
800 Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
801 Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
802 if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
803
2a258f40 804// cout << " N. of sets : " << nseteq << endl;
805
5b4ed716 806 if(ngoodtrkt0<2){
807// cout << "less than 2 tracks, skip event " << endl;
808 t0def=-999;
809 deltat0def=0.600;
810 fT0SigmaT0def[0]=t0def;
811 fT0SigmaT0def[1]=deltat0def;
812 fT0SigmaT0def[2]=ngoodtrkt0;
813 fT0SigmaT0def[3]=ngoodtrkt0;
814 //goto finish;
815 }
816 if(ngoodtrkt0>=2){
817 // Decide how many tracks in set
818 Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
819 Int_t nset=1;
820
821 if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
822
823 // Loop over selected sets
824
825 if(nset>=1){
826 for (Int_t i=0; i< nset; i++) {
2a258f40 827 // printf("Set %i of %i\n",i+1,nset);
5b4ed716 828 Float_t t0best=999.;
829 Float_t eT0best=999.;
830 Float_t chisquarebest=99999.;
831 Int_t npionbest=0;
832
833 Int_t ntracksinsetmy=0;
62f44f07 834 //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset];
5b4ed716 835 for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
836 Int_t index = itrk+i*ntracksinset;
62f44f07 837 //if(index < ngoodtrkt0){
838 if(index < fGTracks->GetEntries()){
839 //AliESDtrack *t=fGTracks[index];
840 AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
841 //fTracksT0[itrk]=t;
842 fTracksT0->AddLast(t);
5b4ed716 843 ntracksinsetmy++;
844 }
845 }
62f44f07 846
5b4ed716 847 // Analyse it
848
2a258f40 849 Int_t assparticle[nmaxtracksinsetMax];
850 Float_t exptof[nmaxtracksinsetMax][3];
851 Float_t timeofflight[nmaxtracksinsetMax];
852 Float_t momentum[nmaxtracksinsetMax];
853 Float_t timezero[nmaxtracksinsetMax];
854 Float_t weightedtimezero[nmaxtracksinsetMax];
855 Float_t beta[nmaxtracksinsetMax];
856 Float_t texp[nmaxtracksinsetMax];
857 Float_t dtexp[nmaxtracksinsetMax];
858 Float_t sqMomError[nmaxtracksinsetMax];
859 Float_t sqTrackError[nmaxtracksinsetMax];
5b4ed716 860 Float_t massarray[3]={0.13957,0.493677,0.9382723};
2a258f40 861 Float_t tracktoflen[nmaxtracksinsetMax];
862 Float_t besttimezero[nmaxtracksinsetMax];
863 Float_t besttexp[nmaxtracksinsetMax];
864 Float_t besttimeofflight[nmaxtracksinsetMax];
865 Float_t bestmomentum[nmaxtracksinsetMax];
866 Float_t bestchisquare[nmaxtracksinsetMax];
867 Float_t bestweightedtimezero[nmaxtracksinsetMax];
868 Float_t bestsqTrackError[nmaxtracksinsetMax];
869 Int_t imass[nmaxtracksinsetMax];
5b4ed716 870
871 for (Int_t j=0; j<ntracksinset; j++) {
872 assparticle[j] = 3;
873 timeofflight[j] = 0;
874 momentum[j] = 0;
875 timezero[j] = 0;
876 weightedtimezero[j] = 0;
877 beta[j] = 0;
878 texp[j] = 0;
879 dtexp[j] = 0;
880 sqMomError[j] = 0;
881 sqTrackError[j] = 0;
882 tracktoflen[j] = 0;
883 besttimezero[j] = 0;
884 besttexp[j] = 0;
885 besttimeofflight[j] = 0;
886 bestmomentum[j] = 0;
887 bestchisquare[j] = 0;
888 bestweightedtimezero[j] = 0;
889 bestsqTrackError[j] = 0;
890 imass[j] = 1;
891 }
892
62f44f07 893 //for (Int_t j=0; j<ntracksinsetmy; j++) {
894 for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
895 //AliESDtrack *t=fTracksT0[j];
896 AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
5b4ed716 897 Double_t momOld=t->GetP();
898 Double_t mom=momOld-0.0036*momOld;
899 Double_t time=t->GetTOFsignal();
900
901 time*=1.E-3; // tof given in nanoseconds
902 Double_t exptime[10]; t->GetIntegratedTimes(exptime);
903 Double_t toflen=t->GetIntegratedLength();
904 toflen=toflen/100.; // toflen given in m
905
906 timeofflight[j]=time;
907 tracktoflen[j]=toflen;
908 exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
909 exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
910 exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
911 momentum[j]=mom;
912 assparticle[j]=3;
913
914 } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
915
916 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
917 beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
918 +momentum[itz]*momentum[itz]);
919 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 920 sqTrackError[itz]=sqMomError[itz]; //in ns
5b4ed716 921 timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
922 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
923 sumAllweightspi+=1./sqTrackError[itz];
924 meantzeropi+=weightedtimezero[itz];
925 } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
926
927
928 // Then, Combinatorial Algorithm
929
930 if(ntracksinsetmy<2 )break;
931
932 for (Int_t j=0; j<ntracksinsetmy; j++) {
933 imass[j] = 3;
934 }
935
62f44f07 936 //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
937 Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
5b4ed716 938
939 // Loop on mass hypotheses
940 for (Int_t k=0; k < ncombinatorial;k++) {
941 for (Int_t j=0; j<ntracksinsetmy; j++) {
62f44f07 942 //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
943 imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
5b4ed716 944 texp[j]=exptof[j][imass[j]];
945 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
1fde62a1 946 // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
5b4ed716 947 }
90beec49 948
5b4ed716 949 Float_t sumAllweights=0.;
950 Float_t meantzero=0.;
951 Float_t eMeanTzero=0.;
952
953 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
2a258f40 954 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 955
956 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
957
958 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
959 sumAllweights+=1./sqTrackError[itz];
960 meantzero+=weightedtimezero[itz];
961
962 } // end loop for (Int_t itz=0; itz<15;itz++)
963
964 meantzero=meantzero/sumAllweights; // it is given in [ns]
965 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
966
967 // calculate chisquare
5b4ed716 968 Float_t chisquare=0.;
969 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1fde62a1 970 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
971 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
972 // else chisquare+=1000;
5b4ed716 973 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
974
975 if(chisquare<=chisquarebest){
976 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
977
978 bestsqTrackError[iqsq]=sqTrackError[iqsq];
979 besttimezero[iqsq]=timezero[iqsq];
980 bestmomentum[iqsq]=momentum[iqsq];
981 besttimeofflight[iqsq]=timeofflight[iqsq];
982 besttexp[iqsq]=texp[iqsq];
983 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 984 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
985 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
986 // else bestchisquare[iqsq]=1000;
5b4ed716 987 }
988
989 Int_t npion=0;
990 for (Int_t j=0; j<ntracksinsetmy; j++) {
991 assparticle[j]=imass[j];
992 if(imass[j] == 0) npion++;
993 }
994 npionbest=npion;
995 chisquarebest=chisquare;
996 t0best=meantzero;
997 eT0best=eMeanTzero;
998 } // close if(dummychisquare<=chisquare)
5b4ed716 999 }
1000
2a258f40 1001 Double_t chi2cut[nmaxtracksinsetMax];
5b4ed716 1002 chi2cut[0] = 0;
1003 chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
1004 for (Int_t j=2; j<ntracksinset; j++) {
1005 chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
1006 }
1007
1008 Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
1009
90beec49 1010// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
5b4ed716 1011
1012 Bool_t kRedoT0 = kFALSE;
1013 ntracksinsetmyCut = ntracksinsetmy;
2a258f40 1014 Bool_t usetrack[nmaxtracksinsetMax];
5b4ed716 1015 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1016 usetrack[icsq] = kTRUE;
1017 if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
90beec49 1018 kRedoT0 = kTRUE;
1019 ntracksinsetmyCut--;
1020 usetrack[icsq] = kFALSE;
5b4ed716 1021 }
1022 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1023
90beec49 1024// printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
2a258f40 1025
5b4ed716 1026 // Loop on mass hypotheses Redo
1027 if(kRedoT0 && ntracksinsetmyCut > 1){
1028 // printf("Redo T0\n");
1029 for (Int_t k=0; k < ncombinatorial;k++) {
1030 for (Int_t j=0; j<ntracksinsetmy; j++) {
62f44f07 1031 //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
1032 imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
5b4ed716 1033 texp[j]=exptof[j][imass[j]];
1034 dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
1fde62a1 1035 // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
5b4ed716 1036 }
1037
1038 Float_t sumAllweights=0.;
1039 Float_t meantzero=0.;
1040 Float_t eMeanTzero=0.;
1041
1042 for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
1043 if(! usetrack[itz]) continue;
2a258f40 1044 sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
5b4ed716 1045
1046 timezero[itz]=texp[itz]-timeofflight[itz];// in ns
1047
1048 weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
1049 sumAllweights+=1./sqTrackError[itz];
1050 meantzero+=weightedtimezero[itz];
1051
1052 } // end loop for (Int_t itz=0; itz<15;itz++)
1053
1054 meantzero=meantzero/sumAllweights; // it is given in [ns]
1055 eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
1056
1057 // calculate chisquare
1058
1059 Float_t chisquare=0.;
1060 for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1061 if(! usetrack[icsq]) continue;
1fde62a1 1062 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
1063 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
1064 // else chisquare+=1000;
5b4ed716 1065 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1066
1067 Int_t npion=0;
1068 for (Int_t j=0; j<ntracksinsetmy; j++) {
1069 assparticle[j]=imass[j];
1070 if(imass[j] == 0) npion++;
1071 }
1072
90beec49 1073 if(chisquare<=chisquarebest && npion>0){
5b4ed716 1074 for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
1075 if(! usetrack[iqsq]) continue;
1076 bestsqTrackError[iqsq]=sqTrackError[iqsq];
1077 besttimezero[iqsq]=timezero[iqsq];
1078 bestmomentum[iqsq]=momentum[iqsq];
1079 besttimeofflight[iqsq]=timeofflight[iqsq];
1080 besttexp[iqsq]=texp[iqsq];
1081 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1fde62a1 1082 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
1083 // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
1084 // else bestchisquare[iqsq]=1000;
5b4ed716 1085 }
1086
1087 npionbest=npion;
1088 chisquarebest=chisquare;
1089 t0best=meantzero;
1090 eT0best=eMeanTzero;
1091 } // close if(dummychisquare<=chisquare)
1092
1093 }
1094 }
1095
1096 // filling histos
1097 Float_t confLevel=999;
1098
1099 // Sets with decent chisquares
2a258f40 1100 // printf("Chi2best of the set = %f \n",chisquarebest);
5b4ed716 1101
1102 if(chisquarebest<999.){
1103 Double_t dblechisquare=(Double_t)chisquarebest;
1104 confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
1105// cout << " Set Number " << nsets << endl;
1106// cout << "Best Assignment, selection " << assparticle[0] <<
1107// assparticle[1] << assparticle[2] <<
1108// assparticle[3] << assparticle[4] <<
1109// assparticle[5] << endl;
1110// cout << " Chisquare of the set "<< chisquarebest <<endl;
1111// cout << " C.L. of the set "<< confLevel <<endl;
1112// cout << " T0 for this set (in ns) " << t0best << endl;
1113
90beec49 1114 Int_t ntrackincurrentsel=0;
5b4ed716 1115 for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1116
1117 if(! usetrack[icsq]) continue;
1118
1119// cout << "Track # " << icsq << " T0 offsets = "
1120// << besttimezero[icsq]-t0best <<
1121// " track error = " << bestsqTrackError[icsq]
1122// << " Chisquare = " << bestchisquare[icsq]
1123// << " Momentum = " << bestmomentum[icsq]
1124// << " TOF = " << besttimeofflight[icsq]
1125// << " TOF tracking = " << besttexp[icsq]
1126// << " is used = " << usetrack[icsq] << endl;
90beec49 1127 ntrackincurrentsel++;
5b4ed716 1128 }
1129
2a258f40 1130 // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
1131
5b4ed716 1132 // Pick up only those with C.L. >1%
1133 // if(confLevel>0.01 && ngoodsetsSel<200){
1134 if(confLevel>0.01 && ngoodsetsSel<200){
1135 chiSquarebestSel[ngoodsetsSel]=chisquarebest;
1136 confLevelbestSel[ngoodsetsSel]=confLevel;
1137 t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
1138 eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
1139 t0bestallSel += t0best/eT0best/eT0best;
1140 sumWt0bestallSel += 1./eT0best/eT0best;
1141 ngoodsetsSel++;
2a258f40 1142 ngoodtrktrulyused+=ntracksinsetmyCut;
90beec49 1143// printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
5b4ed716 1144 }
1145 else{
1146 // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
1147 }
1148 }
62f44f07 1149 //delete[] fTracksT0;
1150 fTracksT0->Clear();
5b4ed716 1151 nsets++;
1152
1153 } // end for the current set
1154
90beec49 1155 //Redo the computation of the best in order to esclude very bad samples
1156 if(ngoodsetsSel > 1){
1157 Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
1158 Int_t nsamples=ngoodsetsSel;
1159 ngoodsetsSel=0;
1160 t0bestallSel=0;
1161 sumWt0bestallSel=0;
1162 for (Int_t itz=0; itz<nsamples;itz++) {
1163 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
1164 t0bestallSel += t0bestSel[itz];
1165 sumWt0bestallSel += eT0bestSel[itz];
1166 ngoodsetsSel++;
1167// printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
1168 }
1169 else{
1170// printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
1171 }
1172 }
1173 }
1174 if(ngoodsetsSel < 1){
1175 sumWt0bestallSel = 0.0;
1176 }
1177 //--------------------------------End recomputation
1178
5b4ed716 1179 nUsedTracks = ngoodtrkt0;
1180 if(strstr(option,"all")){
1181 if(sumAllweightspi>0.){
1182 meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
1183 eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
1184 }
1185
2a258f40 1186 // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
1187
5b4ed716 1188 if(sumWt0bestallSel>0){
1189 t0bestallSel = t0bestallSel/sumWt0bestallSel;
1190 eT0bestallSel = sqrt(1./sumWt0bestallSel);
2a258f40 1191 // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
5b4ed716 1192 }// end of if(sumWt0bestallSel>0){
1193
1194// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
1195 }
1196
1197 t0def=t0bestallSel;
1198 deltat0def=eT0bestallSel;
5b4ed716 1199
1200 fT0SigmaT0def[0]=t0def;
2a258f40 1201 fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
5b4ed716 1202 fT0SigmaT0def[2]=ngoodtrkt0;
1203 fT0SigmaT0def[3]=ngoodtrktrulyused;
1204 }
1205 }
62f44f07 1206
1207 //delete [] fGTracks;
1208 fGTracks->Clear();
1209
5b4ed716 1210 // if(strstr(option,"tim") || strstr(option,"all")){
1211 // cout << "AliTOFT0v1:" << endl ;
1212 //}
1213
2a258f40 1214 if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
1215
1216 bench->Stop("t0computation");
1217
1218 fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
1219 fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
1220
1221// bench->Print("t0computation");
1222// 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]));
1223
90beec49 1224 delete bench;
1225 bench=NULL;
5b4ed716 1226
536031f2 1227 return fT0SigmaT0def;
1228 }
1229//__________________________________________________________________
8f589502 1230Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
536031f2 1231{
8f589502 1232 // Take the error extimate for the TOF time in the track reconstruction
536031f2 1233
536031f2 1234 static const Double_t kMasses[]={
1235 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
1236 };
1237
1238 Double_t mass=kMasses[index+2];
536031f2 1239
2a258f40 1240 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
536031f2 1241
1242 return sigma;
1243}
14b2cbea 1244
1245//__________________________________________________________________
1246Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
1247{
1248
1249 /* TPC refit */
1250 if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1251 /* do not accept kink daughters */
1252 if (track->GetKinkIndex(0)>0) return kFALSE;
1253 /* N clusters TPC */
1254 if (track->GetTPCclusters(0) < 50) return kFALSE;
1255 /* chi2 TPC */
1256 if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
1257 /* sigma to vertex */
1258 if (GetSigmaToVertex(track) > 4.) return kFALSE;
1259
1260 /* accept track */
1261 return kTRUE;
1262
1263}
1264
1265//____________________________________________________________________
8f589502 1266Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
14b2cbea 1267{
1268 // Calculates the number of sigma to the vertex.
1269
1270 Float_t b[2];
1271 Float_t bRes[2];
1272 Float_t bCov[3];
1273 esdTrack->GetImpactParameters(b,bCov);
1274
1275 if (bCov[0]<=0 || bCov[2]<=0) {
1276 bCov[0]=0; bCov[2]=0;
1277 }
1278 bRes[0] = TMath::Sqrt(bCov[0]);
1279 bRes[1] = TMath::Sqrt(bCov[2]);
1280
1281 // -----------------------------------
1282 // How to get to a n-sigma cut?
1283 //
1284 // The accumulated statistics from 0 to d is
1285 //
1286 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1287 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1288 //
1289 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1290 // Can this be expressed in a different way?
1291
1292 if (bRes[0] == 0 || bRes[1] ==0)
1293 return -1;
1294
62f44f07 1295 //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1296 Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
14b2cbea 1297
1298 // work around precision problem
1299 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1300 // 1e-15 corresponds to nsigma ~ 7.7
1301 if (TMath::Exp(-d * d / 2) < 1e-15)
1302 return 1000;
1303
1304 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1305 return nSigma;
1306}
90beec49 1307//____________________________________________________________________
1308
1309Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
1310 Bool_t status = kFALSE;
1311
1312 Double_t exptimes[5];
1313 track->GetIntegratedTimes(exptimes);
1314
1315 Float_t dedx = track->GetTPCsignal();
1316
1317 Double_t ptpc[3];
1318 track->GetInnerPxPyPz(ptpc);
1319 Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
1320
1321 if(imass > 2 || imass < 0) return status;
1322 Int_t i = imass+2;
1323
1324 AliPID::EParticleType type=AliPID::EParticleType(i);
1325
1326 Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
1327 Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
1328
1329 if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
1330 status = kTRUE;
1331 }
1332
1333 return status;
1334}
62f44f07 1335
1336//____________________________________________________________________
1337Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const
1338{
1339 //
1340 // Returns base^exponent
1341 //
1342
1343 Float_t power=1.;
1344
1345 for (Int_t ii=exponent; ii>0; ii--)
1346 power=power*base;
1347
1348 return power;
1349
1350}
1351//____________________________________________________________________
1352Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const
1353{
1354 //
1355 // Returns base^exponent
1356 //
1357
1358 Int_t power=1;
1359
1360 for (Int_t ii=exponent; ii>0; ii--)
1361 power=power*base;
1362
1363 return power;
1364
1365}