]>
Commit | Line | Data |
---|---|---|
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 | |
67 | ClassImp(AliTOFT0v1) | |
68 | ||
69 | //____________________________________________________________________________ | |
2a258f40 | 70 | AliTOFT0v1::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 | 102 | AliTOFT0v1::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 | } |
8f589502 | 129 | //____________________________________________________________________________ |
130 | AliTOFT0v1& 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 | |
536031f2 | 166 | //____________________________________________________________________________ |
167 | AliTOFT0v1::~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 | } |
5b4ed716 | 192 | //____________________________________________________________________________ |
193 | ||
194 | void | |
195 | AliTOFT0v1::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 | 210 | Double_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 | 692 | Float_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 | //__________________________________________________________________ | |
708 | Bool_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 | 728 | Float_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 | ||
771 | Bool_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 | //____________________________________________________________________ | |
799 | Float_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 | //____________________________________________________________________ | |
814 | Int_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 | } |