]>
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 | |
fc9b31a7 | 390 | Double_t exptime[AliPID::kSPECIESC]; |
391 | t->GetIntegratedTimes(exptime,AliPID::kSPECIESC); | |
5b4ed716 | 392 | Double_t toflen=t->GetIntegratedLength(); |
393 | toflen=toflen/100.; // toflen given in m | |
394 | ||
395 | timeofflight[j]=time; | |
396 | tracktoflen[j]=toflen; | |
397 | exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns | |
398 | exptof[j][1]=exptime[3]*1.E-3+fTimeCorr; | |
399 | exptof[j][2]=exptime[4]*1.E-3+fTimeCorr; | |
400 | momentum[j]=mom; | |
401 | assparticle[j]=3; | |
a558aa69 | 402 | |
403 | // in principle GetMomError only depends on two indices k=imass[j] and j itslef (see blow in the ncombinatorial loop) | |
404 | // so it should be possible to make a lookup in order to speed up the code: | |
405 | if (fOptFlag) { | |
406 | momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]); | |
407 | momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]); | |
408 | momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]); | |
409 | } | |
410 | } //end for (Int_t j=0; j<ntracksinsetmy; j++) { | |
5b4ed716 | 411 | |
412 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
413 | beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0] | |
414 | +momentum[itz]*momentum[itz]); | |
415 | sqMomError[itz]= ((1.-beta[itz]*beta[itz])*0.01)*((1.-beta[itz]*beta[itz])*0.01)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz])); | |
2a258f40 | 416 | sqTrackError[itz]=sqMomError[itz]; //in ns |
5b4ed716 | 417 | timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns |
418 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
419 | sumAllweightspi+=1./sqTrackError[itz]; | |
420 | meantzeropi+=weightedtimezero[itz]; | |
421 | } // end loop for (Int_t itz=0; itz< ntracksinset;itz++) | |
422 | ||
423 | ||
424 | // Then, Combinatorial Algorithm | |
425 | ||
426 | if(ntracksinsetmy<2 )break; | |
427 | ||
428 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
429 | imass[j] = 3; | |
430 | } | |
a558aa69 | 431 | |
432 | Int_t ncombinatorial; | |
433 | if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy]; | |
434 | else ncombinatorial = ToCalculatePower(3,ntracksinsetmy); | |
435 | ||
436 | ||
5b4ed716 | 437 | // Loop on mass hypotheses |
438 | for (Int_t k=0; k < ncombinatorial;k++) { | |
439 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
a558aa69 | 440 | imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1]; |
5b4ed716 | 441 | texp[j]=exptof[j][imass[j]]; |
a558aa69 | 442 | if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr |
443 | else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]); | |
5b4ed716 | 444 | } |
90beec49 | 445 | |
5b4ed716 | 446 | Float_t sumAllweights=0.; |
447 | Float_t meantzero=0.; | |
448 | Float_t eMeanTzero=0.; | |
449 | ||
450 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
2a258f40 | 451 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
5b4ed716 | 452 | |
453 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
454 | ||
455 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
456 | sumAllweights+=1./sqTrackError[itz]; | |
457 | meantzero+=weightedtimezero[itz]; | |
458 | ||
459 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
460 | ||
461 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
462 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] | |
463 | ||
464 | // calculate chisquare | |
5b4ed716 | 465 | Float_t chisquare=0.; |
466 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
1fde62a1 | 467 | chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; |
5b4ed716 | 468 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) |
469 | ||
470 | if(chisquare<=chisquarebest){ | |
471 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { | |
472 | ||
473 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
474 | besttimezero[iqsq]=timezero[iqsq]; | |
475 | bestmomentum[iqsq]=momentum[iqsq]; | |
476 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
477 | besttexp[iqsq]=texp[iqsq]; | |
478 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
1fde62a1 | 479 | bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; |
5b4ed716 | 480 | } |
481 | ||
482 | Int_t npion=0; | |
483 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
484 | assparticle[j]=imass[j]; | |
485 | if(imass[j] == 0) npion++; | |
486 | } | |
487 | npionbest=npion; | |
488 | chisquarebest=chisquare; | |
489 | t0best=meantzero; | |
490 | eT0best=eMeanTzero; | |
491 | } // close if(dummychisquare<=chisquare) | |
5b4ed716 | 492 | } |
493 | ||
2a258f40 | 494 | Double_t chi2cut[nmaxtracksinsetMax]; |
5b4ed716 | 495 | chi2cut[0] = 0; |
496 | chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01 | |
497 | for (Int_t j=2; j<ntracksinset; j++) { | |
498 | chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.); | |
499 | } | |
500 | ||
501 | Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy; | |
502 | ||
721ed687 | 503 | // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); |
5b4ed716 | 504 | |
505 | Bool_t kRedoT0 = kFALSE; | |
506 | ntracksinsetmyCut = ntracksinsetmy; | |
2a258f40 | 507 | Bool_t usetrack[nmaxtracksinsetMax]; |
5b4ed716 | 508 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { |
509 | usetrack[icsq] = kTRUE; | |
510 | if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){ | |
90beec49 | 511 | kRedoT0 = kTRUE; |
512 | ntracksinsetmyCut--; | |
513 | usetrack[icsq] = kFALSE; | |
721ed687 | 514 | // printf("tracks chi2 = %f\n",bestchisquare[icsq]); |
5b4ed716 | 515 | } |
516 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) | |
517 | ||
5b4ed716 | 518 | // Loop on mass hypotheses Redo |
519 | if(kRedoT0 && ntracksinsetmyCut > 1){ | |
520 | // printf("Redo T0\n"); | |
521 | for (Int_t k=0; k < ncombinatorial;k++) { | |
522 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
a558aa69 | 523 | imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1]; |
5b4ed716 | 524 | texp[j]=exptof[j][imass[j]]; |
a558aa69 | 525 | if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr |
526 | else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]); | |
5b4ed716 | 527 | } |
528 | ||
529 | Float_t sumAllweights=0.; | |
530 | Float_t meantzero=0.; | |
531 | Float_t eMeanTzero=0.; | |
532 | ||
533 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
534 | if(! usetrack[itz]) continue; | |
2a258f40 | 535 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
5b4ed716 | 536 | |
537 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
538 | ||
539 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
540 | sumAllweights+=1./sqTrackError[itz]; | |
541 | meantzero+=weightedtimezero[itz]; | |
542 | ||
543 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
544 | ||
545 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
546 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] | |
547 | ||
548 | // calculate chisquare | |
549 | ||
550 | Float_t chisquare=0.; | |
551 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
552 | if(! usetrack[icsq]) continue; | |
1fde62a1 | 553 | chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; |
5b4ed716 | 554 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) |
555 | ||
556 | Int_t npion=0; | |
557 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
558 | assparticle[j]=imass[j]; | |
559 | if(imass[j] == 0) npion++; | |
560 | } | |
561 | ||
90beec49 | 562 | if(chisquare<=chisquarebest && npion>0){ |
5b4ed716 | 563 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { |
564 | if(! usetrack[iqsq]) continue; | |
565 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
566 | besttimezero[iqsq]=timezero[iqsq]; | |
567 | bestmomentum[iqsq]=momentum[iqsq]; | |
568 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
569 | besttexp[iqsq]=texp[iqsq]; | |
570 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
1fde62a1 | 571 | bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; |
5b4ed716 | 572 | } |
573 | ||
574 | npionbest=npion; | |
575 | chisquarebest=chisquare; | |
576 | t0best=meantzero; | |
577 | eT0best=eMeanTzero; | |
578 | } // close if(dummychisquare<=chisquare) | |
579 | ||
580 | } | |
581 | } | |
582 | ||
583 | // filling histos | |
584 | Float_t confLevel=999; | |
585 | ||
586 | // Sets with decent chisquares | |
2a258f40 | 587 | // printf("Chi2best of the set = %f \n",chisquarebest); |
5b4ed716 | 588 | |
589 | if(chisquarebest<999.){ | |
590 | Double_t dblechisquare=(Double_t)chisquarebest; | |
591 | confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); | |
5b4ed716 | 592 | |
90beec49 | 593 | Int_t ntrackincurrentsel=0; |
5b4ed716 | 594 | for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){ |
595 | ||
596 | if(! usetrack[icsq]) continue; | |
597 | ||
90beec49 | 598 | ntrackincurrentsel++; |
5b4ed716 | 599 | } |
600 | ||
2a258f40 | 601 | // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel); |
602 | ||
5b4ed716 | 603 | // Pick up only those with C.L. >1% |
5b4ed716 | 604 | if(confLevel>0.01 && ngoodsetsSel<200){ |
605 | chiSquarebestSel[ngoodsetsSel]=chisquarebest; | |
606 | confLevelbestSel[ngoodsetsSel]=confLevel; | |
607 | t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best; | |
608 | eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best; | |
609 | t0bestallSel += t0best/eT0best/eT0best; | |
610 | sumWt0bestallSel += 1./eT0best/eT0best; | |
611 | ngoodsetsSel++; | |
2a258f40 | 612 | ngoodtrktrulyused+=ntracksinsetmyCut; |
721ed687 | 613 | // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel); |
5b4ed716 | 614 | } |
615 | else{ | |
616 | // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy); | |
617 | } | |
618 | } | |
62f44f07 | 619 | fTracksT0->Clear(); |
5b4ed716 | 620 | nsets++; |
621 | ||
622 | } // end for the current set | |
623 | ||
90beec49 | 624 | //Redo the computation of the best in order to esclude very bad samples |
625 | if(ngoodsetsSel > 1){ | |
626 | Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel; | |
627 | Int_t nsamples=ngoodsetsSel; | |
628 | ngoodsetsSel=0; | |
629 | t0bestallSel=0; | |
630 | sumWt0bestallSel=0; | |
631 | for (Int_t itz=0; itz<nsamples;itz++) { | |
632 | if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){ | |
633 | t0bestallSel += t0bestSel[itz]; | |
634 | sumWt0bestallSel += eT0bestSel[itz]; | |
635 | ngoodsetsSel++; | |
721ed687 | 636 | // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); |
90beec49 | 637 | } |
638 | else{ | |
721ed687 | 639 | // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); |
90beec49 | 640 | } |
641 | } | |
642 | } | |
643 | if(ngoodsetsSel < 1){ | |
644 | sumWt0bestallSel = 0.0; | |
645 | } | |
646 | //--------------------------------End recomputation | |
647 | ||
5b4ed716 | 648 | nUsedTracks = ngoodtrkt0; |
649 | if(strstr(option,"all")){ | |
650 | if(sumAllweightspi>0.){ | |
651 | meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns] | |
652 | eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns] | |
653 | } | |
654 | ||
2a258f40 | 655 | // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel); |
656 | ||
5b4ed716 | 657 | if(sumWt0bestallSel>0){ |
658 | t0bestallSel = t0bestallSel/sumWt0bestallSel; | |
659 | eT0bestallSel = sqrt(1./sumWt0bestallSel); | |
2a258f40 | 660 | // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel); |
5b4ed716 | 661 | }// end of if(sumWt0bestallSel>0){ |
662 | ||
5b4ed716 | 663 | } |
664 | ||
665 | t0def=t0bestallSel; | |
666 | deltat0def=eT0bestallSel; | |
5b4ed716 | 667 | |
668 | fT0SigmaT0def[0]=t0def; | |
2a258f40 | 669 | fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1)); |
5b4ed716 | 670 | fT0SigmaT0def[2]=ngoodtrkt0; |
671 | fT0SigmaT0def[3]=ngoodtrktrulyused; | |
672 | } | |
673 | } | |
62f44f07 | 674 | |
62f44f07 | 675 | fGTracks->Clear(); |
676 | ||
2a258f40 | 677 | if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6; |
678 | ||
679 | bench->Stop("t0computation"); | |
680 | ||
681 | fT0SigmaT0def[4]=bench->GetRealTime("t0computation"); | |
682 | fT0SigmaT0def[5]=bench->GetCpuTime("t0computation"); | |
683 | ||
684 | // bench->Print("t0computation"); | |
685 | // printf("(%4.1f < p < %4.1f GeV/c) T0-TOF =%9.1f +/- %5.1f ps (n_track = %i)\n\n",pMinCut,pMaxCut,-fT0SigmaT0def[0]*1000,fT0SigmaT0def[1]*1000,Int_t(fT0SigmaT0def[3])); | |
686 | ||
90beec49 | 687 | delete bench; |
688 | bench=NULL; | |
5b4ed716 | 689 | |
536031f2 | 690 | return fT0SigmaT0def; |
691 | } | |
692 | //__________________________________________________________________ | |
8f589502 | 693 | Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const |
536031f2 | 694 | { |
8f589502 | 695 | // Take the error extimate for the TOF time in the track reconstruction |
536031f2 | 696 | |
536031f2 | 697 | static const Double_t kMasses[]={ |
698 | 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 | |
699 | }; | |
700 | ||
701 | Double_t mass=kMasses[index+2]; | |
536031f2 | 702 | |
2a258f40 | 703 | Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass); |
536031f2 | 704 | |
705 | return sigma; | |
706 | } | |
14b2cbea | 707 | |
708 | //__________________________________________________________________ | |
709 | Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track) | |
710 | { | |
711 | ||
712 | /* TPC refit */ | |
713 | if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE; | |
714 | /* do not accept kink daughters */ | |
715 | if (track->GetKinkIndex(0)>0) return kFALSE; | |
716 | /* N clusters TPC */ | |
717 | if (track->GetTPCclusters(0) < 50) return kFALSE; | |
718 | /* chi2 TPC */ | |
719 | if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE; | |
720 | /* sigma to vertex */ | |
721 | if (GetSigmaToVertex(track) > 4.) return kFALSE; | |
722 | ||
723 | /* accept track */ | |
724 | return kTRUE; | |
725 | ||
726 | } | |
727 | ||
728 | //____________________________________________________________________ | |
8f589502 | 729 | Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const |
14b2cbea | 730 | { |
731 | // Calculates the number of sigma to the vertex. | |
732 | ||
733 | Float_t b[2]; | |
734 | Float_t bRes[2]; | |
735 | Float_t bCov[3]; | |
736 | esdTrack->GetImpactParameters(b,bCov); | |
737 | ||
738 | if (bCov[0]<=0 || bCov[2]<=0) { | |
739 | bCov[0]=0; bCov[2]=0; | |
740 | } | |
741 | bRes[0] = TMath::Sqrt(bCov[0]); | |
742 | bRes[1] = TMath::Sqrt(bCov[2]); | |
743 | ||
744 | // ----------------------------------- | |
745 | // How to get to a n-sigma cut? | |
746 | // | |
747 | // The accumulated statistics from 0 to d is | |
748 | // | |
749 | // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma) | |
750 | // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma) | |
751 | // | |
752 | // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2) | |
753 | // Can this be expressed in a different way? | |
754 | ||
755 | if (bRes[0] == 0 || bRes[1] ==0) | |
756 | return -1; | |
757 | ||
62f44f07 | 758 | //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); |
759 | Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2)); | |
14b2cbea | 760 | |
761 | // work around precision problem | |
762 | // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :( | |
763 | // 1e-15 corresponds to nsigma ~ 7.7 | |
764 | if (TMath::Exp(-d * d / 2) < 1e-15) | |
765 | return 1000; | |
766 | ||
767 | Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); | |
768 | return nSigma; | |
769 | } | |
90beec49 | 770 | //____________________________________________________________________ |
771 | ||
772 | Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{ | |
773 | Bool_t status = kFALSE; | |
774 | ||
d2dbb0fe | 775 | Double_t exptimes[AliPID::kSPECIESC]; |
fc9b31a7 | 776 | track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC); |
90beec49 | 777 | |
778 | Float_t dedx = track->GetTPCsignal(); | |
779 | ||
780 | Double_t ptpc[3]; | |
781 | track->GetInnerPxPyPz(ptpc); | |
782 | Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]); | |
783 | ||
784 | if(imass > 2 || imass < 0) return status; | |
785 | Int_t i = imass+2; | |
786 | ||
787 | AliPID::EParticleType type=AliPID::EParticleType(i); | |
788 | ||
789 | Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type); | |
790 | Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type); | |
791 | ||
792 | if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){ | |
793 | status = kTRUE; | |
794 | } | |
795 | ||
796 | return status; | |
797 | } | |
62f44f07 | 798 | |
799 | //____________________________________________________________________ | |
800 | Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const | |
801 | { | |
802 | // | |
803 | // Returns base^exponent | |
804 | // | |
805 | ||
806 | Float_t power=1.; | |
807 | ||
808 | for (Int_t ii=exponent; ii>0; ii--) | |
809 | power=power*base; | |
810 | ||
811 | return power; | |
812 | ||
813 | } | |
814 | //____________________________________________________________________ | |
815 | Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const | |
816 | { | |
817 | // | |
818 | // Returns base^exponent | |
819 | // | |
820 | ||
821 | Int_t power=1; | |
822 | ||
823 | for (Int_t ii=exponent; ii>0; ii--) | |
824 | power=power*base; | |
825 | ||
826 | return power; | |
827 | ||
828 | } |