]>
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 |
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 | |
65 | ClassImp(AliTOFT0v1) | |
66 | ||
67 | //____________________________________________________________________________ | |
2a258f40 | 68 | AliTOFT0v1::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 | 93 | AliTOFT0v1::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 | //____________________________________________________________________________ |
117 | AliTOFT0v1& 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 | //____________________________________________________________________________ |
152 | AliTOFT0v1::~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 | ||
172 | void | |
173 | AliTOFT0v1::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 | //____________________________________________________________________________ |
189 | Double_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]); | |
62f44f07 | 424 | //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100; |
425 | if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100; | |
536031f2 | 426 | } |
90beec49 | 427 | |
536031f2 | 428 | Float_t sumAllweights=0.; |
429 | Float_t meantzero=0.; | |
8f589502 | 430 | Float_t eMeanTzero=0.; |
536031f2 | 431 | |
432 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
2a258f40 | 433 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
536031f2 | 434 | |
435 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
436 | ||
437 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
438 | sumAllweights+=1./sqTrackError[itz]; | |
439 | meantzero+=weightedtimezero[itz]; | |
440 | ||
441 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
442 | ||
443 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
8f589502 | 444 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] |
536031f2 | 445 | |
446 | // calculate chisquare | |
447 | ||
448 | Float_t chisquare=0.; | |
449 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
62f44f07 | 450 | //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement |
451 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement | |
452 | else chisquare+=1000; | |
536031f2 | 453 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) |
454 | ||
455 | if(chisquare<=chisquarebest){ | |
456 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { | |
457 | ||
458 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
459 | besttimezero[iqsq]=timezero[iqsq]; | |
460 | bestmomentum[iqsq]=momentum[iqsq]; | |
461 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
462 | besttexp[iqsq]=texp[iqsq]; | |
463 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
62f44f07 | 464 | //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement |
465 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement | |
90beec49 | 466 | else bestchisquare[iqsq]=1000; |
536031f2 | 467 | } |
468 | ||
469 | Int_t npion=0; | |
470 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
471 | assparticle[j]=imass[j]; | |
472 | if(imass[j] == 0) npion++; | |
473 | } | |
474 | npionbest=npion; | |
475 | chisquarebest=chisquare; | |
476 | t0best=meantzero; | |
8f589502 | 477 | eT0best=eMeanTzero; |
536031f2 | 478 | } // close if(dummychisquare<=chisquare) |
479 | ||
480 | } | |
481 | ||
2a258f40 | 482 | Double_t chi2cut[nmaxtracksinsetMax]; |
536031f2 | 483 | chi2cut[0] = 0; |
484 | chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01 | |
485 | for (Int_t j=2; j<ntracksinset; j++) { | |
486 | chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.); | |
487 | } | |
488 | ||
489 | Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy; | |
490 | ||
03bd764d | 491 | // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); |
536031f2 | 492 | |
493 | Bool_t kRedoT0 = kFALSE; | |
494 | ntracksinsetmyCut = ntracksinsetmy; | |
2a258f40 | 495 | Bool_t usetrack[nmaxtracksinsetMax]; |
536031f2 | 496 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { |
497 | usetrack[icsq] = kTRUE; | |
498 | if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){ | |
499 | kRedoT0 = kTRUE; | |
500 | ntracksinsetmyCut--; | |
501 | usetrack[icsq] = kFALSE; | |
502 | } | |
503 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) | |
504 | ||
6a4e212e | 505 | // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut); |
536031f2 | 506 | |
507 | // Loop on mass hypotheses Redo | |
508 | if(kRedoT0 && ntracksinsetmyCut > 1){ | |
6a4e212e | 509 | // printf("Redo T0\n"); |
536031f2 | 510 | for (Int_t k=0; k < ncombinatorial;k++) { |
511 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
62f44f07 | 512 | //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1)); |
513 | imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1); | |
536031f2 | 514 | texp[j]=exptof[j][imass[j]]; |
515 | dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]); | |
62f44f07 | 516 | //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100; |
517 | if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100; | |
536031f2 | 518 | } |
519 | ||
520 | Float_t sumAllweights=0.; | |
521 | Float_t meantzero=0.; | |
8f589502 | 522 | Float_t eMeanTzero=0.; |
536031f2 | 523 | |
524 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
525 | if(! usetrack[itz]) continue; | |
2a258f40 | 526 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
536031f2 | 527 | |
528 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
529 | ||
530 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
531 | sumAllweights+=1./sqTrackError[itz]; | |
532 | meantzero+=weightedtimezero[itz]; | |
533 | ||
534 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
535 | ||
536 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
8f589502 | 537 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] |
536031f2 | 538 | |
539 | // calculate chisquare | |
540 | ||
541 | Float_t chisquare=0.; | |
542 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
543 | if(! usetrack[icsq]) continue; | |
62f44f07 | 544 | //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement |
545 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement | |
90beec49 | 546 | else chisquare+=1000; |
536031f2 | 547 | |
548 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) | |
549 | ||
550 | Int_t npion=0; | |
551 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
552 | assparticle[j]=imass[j]; | |
553 | if(imass[j] == 0) npion++; | |
554 | } | |
555 | ||
90beec49 | 556 | if(chisquare<=chisquarebest && npion>0){ |
536031f2 | 557 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { |
558 | if(! usetrack[iqsq]) continue; | |
559 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
560 | besttimezero[iqsq]=timezero[iqsq]; | |
561 | bestmomentum[iqsq]=momentum[iqsq]; | |
562 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
563 | besttexp[iqsq]=texp[iqsq]; | |
564 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
62f44f07 | 565 | //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement |
566 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement | |
90beec49 | 567 | else bestchisquare[iqsq]=1000; |
536031f2 | 568 | } |
569 | ||
570 | npionbest=npion; | |
571 | chisquarebest=chisquare; | |
572 | t0best=meantzero; | |
8f589502 | 573 | eT0best=eMeanTzero; |
536031f2 | 574 | } // close if(dummychisquare<=chisquare) |
575 | ||
576 | } | |
577 | } | |
6a4e212e | 578 | |
536031f2 | 579 | // filling histos |
580 | Float_t confLevel=999; | |
581 | ||
582 | // Sets with decent chisquares | |
583 | ||
584 | if(chisquarebest<999.){ | |
585 | Double_t dblechisquare=(Double_t)chisquarebest; | |
586 | confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); | |
03bd764d | 587 | // cout << " Set Number " << nsets << endl; |
588 | // cout << "Best Assignment, selection " << assparticle[0] << | |
589 | // assparticle[1] << assparticle[2] << | |
590 | // assparticle[3] << assparticle[4] << | |
591 | // assparticle[5] << endl; | |
592 | // cout << " Chisquare of the set "<< chisquarebest <<endl; | |
593 | // cout << " C.L. of the set "<< confLevel <<endl; | |
594 | // cout << " T0 for this set (in ns) " << t0best << endl; | |
536031f2 | 595 | |
596 | for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){ | |
597 | ||
598 | if(! usetrack[icsq]) continue; | |
599 | ||
03bd764d | 600 | // cout << "Track # " << icsq << " T0 offsets = " |
601 | // << besttimezero[icsq]-t0best << | |
602 | // " track error = " << bestsqTrackError[icsq] | |
603 | // << " Chisquare = " << bestchisquare[icsq] | |
604 | // << " Momentum = " << bestmomentum[icsq] | |
605 | // << " TOF = " << besttimeofflight[icsq] | |
606 | // << " TOF tracking = " << besttexp[icsq] | |
607 | // << " is used = " << usetrack[icsq] << endl; | |
536031f2 | 608 | } |
609 | ||
610 | // Pick up only those with C.L. >1% | |
611 | // if(confLevel>0.01 && ngoodsetsSel<200){ | |
612 | if(confLevel>0.01 && ngoodsetsSel<200){ | |
8f589502 | 613 | chiSquarebestSel[ngoodsetsSel]=chisquarebest; |
614 | confLevelbestSel[ngoodsetsSel]=confLevel; | |
615 | t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best; | |
616 | eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best; | |
617 | t0bestallSel += t0best/eT0best/eT0best; | |
618 | sumWt0bestallSel += 1./eT0best/eT0best; | |
536031f2 | 619 | ngoodsetsSel++; |
620 | ngoodtrktrulyused+=ntracksinsetmyCut; | |
2a258f40 | 621 | |
622 | // printf("t0 of a set = %f +/- %f\n",t0best,eT0best); | |
536031f2 | 623 | } |
624 | else{ | |
6a4e212e | 625 | // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy); |
536031f2 | 626 | } |
627 | } | |
62f44f07 | 628 | //delete[] fTracksT0; |
629 | fTracksT0->Clear(); | |
536031f2 | 630 | nsets++; |
631 | ||
632 | } // end for the current set | |
633 | ||
90beec49 | 634 | //Redo the computation of the best in order to esclude very bad samples |
635 | if(ngoodsetsSel > 1){ | |
636 | Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel; | |
637 | Int_t nsamples=ngoodsetsSel; | |
638 | ngoodsetsSel=0; | |
639 | t0bestallSel=0; | |
640 | sumWt0bestallSel=0; | |
641 | for (Int_t itz=0; itz<nsamples;itz++) { | |
642 | if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){ | |
643 | t0bestallSel += t0bestSel[itz]; | |
644 | sumWt0bestallSel += eT0bestSel[itz]; | |
645 | ngoodsetsSel++; | |
646 | // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); | |
647 | } | |
648 | else{ | |
649 | // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); | |
650 | } | |
651 | } | |
652 | } | |
653 | if(ngoodsetsSel < 1){ | |
654 | sumWt0bestallSel = 0.0; | |
655 | } | |
656 | //--------------------------------End recomputation | |
657 | ||
8f589502 | 658 | nUsedTracks = ngoodtrkt0; |
536031f2 | 659 | if(strstr(option,"all")){ |
660 | if(sumAllweightspi>0.){ | |
661 | meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns] | |
8f589502 | 662 | eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns] |
536031f2 | 663 | } |
664 | ||
665 | if(sumWt0bestallSel>0){ | |
666 | t0bestallSel = t0bestallSel/sumWt0bestallSel; | |
8f589502 | 667 | eT0bestallSel = sqrt(1./sumWt0bestallSel); |
536031f2 | 668 | |
669 | }// end of if(sumWt0bestallSel>0){ | |
670 | ||
6a4e212e | 671 | // cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl; |
536031f2 | 672 | } |
673 | ||
674 | t0def=t0bestallSel; | |
8f589502 | 675 | deltat0def=eT0bestallSel; |
536031f2 | 676 | |
677 | fT0SigmaT0def[0]=t0def; | |
2a258f40 | 678 | fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1)); |
6a4e212e | 679 | fT0SigmaT0def[2]=ngoodtrkt0; |
536031f2 | 680 | fT0SigmaT0def[3]=ngoodtrktrulyused; |
681 | } | |
682 | } | |
62f44f07 | 683 | |
684 | //delete [] fGTracks; | |
685 | fGTracks->Clear(); | |
686 | ||
6a4e212e | 687 | // if(strstr(option,"tim") || strstr(option,"all")){ |
688 | // cout << "AliTOFT0v1:" << endl ; | |
689 | //} | |
536031f2 | 690 | |
2a258f40 | 691 | if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6; |
692 | ||
693 | bench->Stop("t0computation"); | |
694 | ||
695 | fT0SigmaT0def[4]=bench->GetRealTime("t0computation"); | |
696 | fT0SigmaT0def[5]=bench->GetCpuTime("t0computation"); | |
697 | ||
698 | // bench->Print("t0computation"); | |
699 | // printf("%f %f\n",fT0SigmaT0def[4],fT0SigmaT0def[5]); | |
700 | ||
90beec49 | 701 | |
702 | delete bench; | |
703 | bench=NULL; | |
5b4ed716 | 704 | |
705 | return fT0SigmaT0def; | |
706 | } | |
707 | //__________________________________________________________________ | |
708 | Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut) | |
709 | { | |
2a258f40 | 710 | TBenchmark *bench=new TBenchmark(); |
711 | bench->Start("t0computation"); | |
712 | ||
5b4ed716 | 713 | // Caluclate the Event Time using the ESD TOF time |
714 | ||
715 | fT0SigmaT0def[0]=0.; | |
716 | fT0SigmaT0def[1]=0.600; | |
717 | fT0SigmaT0def[2]=0.; | |
718 | fT0SigmaT0def[3]=0.; | |
5b4ed716 | 719 | |
2a258f40 | 720 | const Int_t nmaxtracksinsetMax=10; |
721 | Int_t nmaxtracksinset=10; | |
5b4ed716 | 722 | // if(strstr(option,"all")){ |
723 | // cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl; | |
724 | // cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl; | |
725 | // } | |
726 | ||
727 | ||
728 | Int_t nsets=0; | |
729 | Int_t nUsedTracks=0; | |
730 | Int_t ngoodsetsSel= 0; | |
731 | Float_t t0bestSel[300]; | |
732 | Float_t eT0bestSel[300]; | |
733 | Float_t chiSquarebestSel[300]; | |
734 | Float_t confLevelbestSel[300]; | |
735 | Float_t t0bestallSel=0.; | |
736 | Float_t eT0bestallSel=0.; | |
737 | Float_t sumWt0bestallSel=0.; | |
738 | Float_t eMeanTzeroPi=0.; | |
739 | Float_t meantzeropi=0.; | |
740 | Float_t sumAllweightspi=0.; | |
741 | Double_t t0def=-999; | |
742 | Double_t deltat0def=999; | |
743 | Int_t ngoodtrktrulyused=0; | |
744 | Int_t ntracksinsetmyCut = 0; | |
745 | ||
746 | Int_t ntrk=fEvent->GetNumberOfTracks(); | |
747 | ||
62f44f07 | 748 | //AliESDtrack **fTracks=new AliESDtrack*[ntrk]; |
5b4ed716 | 749 | Int_t ngoodtrk=0; |
750 | Int_t ngoodtrkt0 =0; | |
751 | Float_t mintime =1E6; | |
752 | ||
753 | // First Track loop, Selection of good tracks | |
754 | ||
755 | for (Int_t itrk=0; itrk<ntrk; itrk++) { | |
756 | AliESDtrack *t=fEvent->GetTrack(itrk); | |
757 | Double_t momOld=t->GetP(); | |
758 | Double_t mom=momOld-0.0036*momOld; | |
759 | if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue; | |
760 | if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue; | |
761 | Double_t time=t->GetTOFsignal(); | |
762 | ||
763 | time*=1.E-3; // tof given in nanoseconds | |
764 | if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue; | |
2a258f40 | 765 | |
5b4ed716 | 766 | if (!AcceptTrack(t)) continue; |
767 | ||
9c89b8bd | 768 | if(t->GetIntegratedLength() < 350)continue; //skip decays |
5b4ed716 | 769 | if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue; |
770 | if(time <= mintime) mintime=time; | |
62f44f07 | 771 | //fTracks[ngoodtrk]=t; |
772 | fTracks->AddLast(t); | |
5b4ed716 | 773 | ngoodtrk++; |
774 | } | |
775 | ||
2a258f40 | 776 | if(ngoodtrk>22) nmaxtracksinset = 6; |
777 | ||
5b4ed716 | 778 | |
2a258f40 | 779 | // cout << " N. of ESD tracks : " << ntrk << endl; |
780 | // cout << " N. of preselected tracks : " << ngoodtrk << endl; | |
781 | // cout << " Minimum tof time in set (in ns) : " << mintime << endl; | |
5b4ed716 | 782 | |
62f44f07 | 783 | //AliESDtrack **fGTracks=new AliESDtrack*[ngoodtrk]; |
5b4ed716 | 784 | |
62f44f07 | 785 | //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) { |
786 | for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) { | |
787 | //AliESDtrack *t=fTracks[jtrk]; | |
788 | AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk); | |
5b4ed716 | 789 | Double_t time=t->GetTOFsignal(); |
790 | ||
791 | if((time-mintime*1.E3)<50.E3){ // For pp and per | |
62f44f07 | 792 | //fGTracks[ngoodtrkt0]=t; |
793 | fGTracks->AddLast(t); | |
5b4ed716 | 794 | ngoodtrkt0++; |
795 | } | |
796 | } | |
62f44f07 | 797 | |
798 | //delete [] fTracks; | |
799 | fTracks->Clear(); | |
800 | ||
2a258f40 | 801 | // cout << " N. of preselected tracks t0 : " << ngoodtrkt0 << endl; |
5b4ed716 | 802 | |
803 | Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1; | |
804 | Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq; | |
805 | if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++; | |
806 | ||
2a258f40 | 807 | // cout << " N. of sets : " << nseteq << endl; |
808 | ||
5b4ed716 | 809 | if(ngoodtrkt0<2){ |
810 | // cout << "less than 2 tracks, skip event " << endl; | |
811 | t0def=-999; | |
812 | deltat0def=0.600; | |
813 | fT0SigmaT0def[0]=t0def; | |
814 | fT0SigmaT0def[1]=deltat0def; | |
815 | fT0SigmaT0def[2]=ngoodtrkt0; | |
816 | fT0SigmaT0def[3]=ngoodtrkt0; | |
817 | //goto finish; | |
818 | } | |
819 | if(ngoodtrkt0>=2){ | |
820 | // Decide how many tracks in set | |
821 | Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent); | |
822 | Int_t nset=1; | |
823 | ||
824 | if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} | |
825 | ||
826 | // Loop over selected sets | |
827 | ||
828 | if(nset>=1){ | |
829 | for (Int_t i=0; i< nset; i++) { | |
2a258f40 | 830 | // printf("Set %i of %i\n",i+1,nset); |
5b4ed716 | 831 | Float_t t0best=999.; |
832 | Float_t eT0best=999.; | |
833 | Float_t chisquarebest=99999.; | |
834 | Int_t npionbest=0; | |
835 | ||
836 | Int_t ntracksinsetmy=0; | |
62f44f07 | 837 | //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset]; |
5b4ed716 | 838 | for (Int_t itrk=0; itrk<ntracksinset; itrk++) { |
839 | Int_t index = itrk+i*ntracksinset; | |
62f44f07 | 840 | //if(index < ngoodtrkt0){ |
841 | if(index < fGTracks->GetEntries()){ | |
842 | //AliESDtrack *t=fGTracks[index]; | |
843 | AliESDtrack *t=(AliESDtrack*)fGTracks->At(index); | |
844 | //fTracksT0[itrk]=t; | |
845 | fTracksT0->AddLast(t); | |
5b4ed716 | 846 | ntracksinsetmy++; |
847 | } | |
848 | } | |
62f44f07 | 849 | |
5b4ed716 | 850 | // Analyse it |
851 | ||
2a258f40 | 852 | Int_t assparticle[nmaxtracksinsetMax]; |
853 | Float_t exptof[nmaxtracksinsetMax][3]; | |
854 | Float_t timeofflight[nmaxtracksinsetMax]; | |
855 | Float_t momentum[nmaxtracksinsetMax]; | |
856 | Float_t timezero[nmaxtracksinsetMax]; | |
857 | Float_t weightedtimezero[nmaxtracksinsetMax]; | |
858 | Float_t beta[nmaxtracksinsetMax]; | |
859 | Float_t texp[nmaxtracksinsetMax]; | |
860 | Float_t dtexp[nmaxtracksinsetMax]; | |
861 | Float_t sqMomError[nmaxtracksinsetMax]; | |
862 | Float_t sqTrackError[nmaxtracksinsetMax]; | |
5b4ed716 | 863 | Float_t massarray[3]={0.13957,0.493677,0.9382723}; |
2a258f40 | 864 | Float_t tracktoflen[nmaxtracksinsetMax]; |
865 | Float_t besttimezero[nmaxtracksinsetMax]; | |
866 | Float_t besttexp[nmaxtracksinsetMax]; | |
867 | Float_t besttimeofflight[nmaxtracksinsetMax]; | |
868 | Float_t bestmomentum[nmaxtracksinsetMax]; | |
869 | Float_t bestchisquare[nmaxtracksinsetMax]; | |
870 | Float_t bestweightedtimezero[nmaxtracksinsetMax]; | |
871 | Float_t bestsqTrackError[nmaxtracksinsetMax]; | |
872 | Int_t imass[nmaxtracksinsetMax]; | |
5b4ed716 | 873 | |
874 | for (Int_t j=0; j<ntracksinset; j++) { | |
875 | assparticle[j] = 3; | |
876 | timeofflight[j] = 0; | |
877 | momentum[j] = 0; | |
878 | timezero[j] = 0; | |
879 | weightedtimezero[j] = 0; | |
880 | beta[j] = 0; | |
881 | texp[j] = 0; | |
882 | dtexp[j] = 0; | |
883 | sqMomError[j] = 0; | |
884 | sqTrackError[j] = 0; | |
885 | tracktoflen[j] = 0; | |
886 | besttimezero[j] = 0; | |
887 | besttexp[j] = 0; | |
888 | besttimeofflight[j] = 0; | |
889 | bestmomentum[j] = 0; | |
890 | bestchisquare[j] = 0; | |
891 | bestweightedtimezero[j] = 0; | |
892 | bestsqTrackError[j] = 0; | |
893 | imass[j] = 1; | |
894 | } | |
895 | ||
62f44f07 | 896 | //for (Int_t j=0; j<ntracksinsetmy; j++) { |
897 | for (Int_t j=0; j<fTracksT0->GetEntries(); j++) { | |
898 | //AliESDtrack *t=fTracksT0[j]; | |
899 | AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j); | |
5b4ed716 | 900 | Double_t momOld=t->GetP(); |
901 | Double_t mom=momOld-0.0036*momOld; | |
902 | Double_t time=t->GetTOFsignal(); | |
903 | ||
904 | time*=1.E-3; // tof given in nanoseconds | |
905 | Double_t exptime[10]; t->GetIntegratedTimes(exptime); | |
906 | Double_t toflen=t->GetIntegratedLength(); | |
907 | toflen=toflen/100.; // toflen given in m | |
908 | ||
909 | timeofflight[j]=time; | |
910 | tracktoflen[j]=toflen; | |
911 | exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns | |
912 | exptof[j][1]=exptime[3]*1.E-3+fTimeCorr; | |
913 | exptof[j][2]=exptime[4]*1.E-3+fTimeCorr; | |
914 | momentum[j]=mom; | |
915 | assparticle[j]=3; | |
916 | ||
917 | } //end for (Int_t j=0; j<ntracksinsetmy; j++) { | |
918 | ||
919 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
920 | beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0] | |
921 | +momentum[itz]*momentum[itz]); | |
922 | 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 | 923 | sqTrackError[itz]=sqMomError[itz]; //in ns |
5b4ed716 | 924 | timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns |
925 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
926 | sumAllweightspi+=1./sqTrackError[itz]; | |
927 | meantzeropi+=weightedtimezero[itz]; | |
928 | } // end loop for (Int_t itz=0; itz< ntracksinset;itz++) | |
929 | ||
930 | ||
931 | // Then, Combinatorial Algorithm | |
932 | ||
933 | if(ntracksinsetmy<2 )break; | |
934 | ||
935 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
936 | imass[j] = 3; | |
937 | } | |
938 | ||
62f44f07 | 939 | //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy)); |
940 | Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy); | |
5b4ed716 | 941 | |
942 | // Loop on mass hypotheses | |
943 | for (Int_t k=0; k < ncombinatorial;k++) { | |
944 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
62f44f07 | 945 | //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1)); |
946 | imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1); | |
5b4ed716 | 947 | texp[j]=exptof[j][imass[j]]; |
948 | dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]); | |
62f44f07 | 949 | //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100; |
950 | if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100; | |
5b4ed716 | 951 | } |
90beec49 | 952 | |
5b4ed716 | 953 | Float_t sumAllweights=0.; |
954 | Float_t meantzero=0.; | |
955 | Float_t eMeanTzero=0.; | |
956 | ||
957 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
2a258f40 | 958 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
5b4ed716 | 959 | |
960 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
961 | ||
962 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
963 | sumAllweights+=1./sqTrackError[itz]; | |
964 | meantzero+=weightedtimezero[itz]; | |
965 | ||
966 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
967 | ||
968 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
969 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] | |
970 | ||
971 | // calculate chisquare | |
5b4ed716 | 972 | Float_t chisquare=0.; |
973 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
62f44f07 | 974 | //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement |
975 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement | |
976 | else chisquare+=1000; | |
5b4ed716 | 977 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) |
978 | ||
979 | if(chisquare<=chisquarebest){ | |
980 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { | |
981 | ||
982 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
983 | besttimezero[iqsq]=timezero[iqsq]; | |
984 | bestmomentum[iqsq]=momentum[iqsq]; | |
985 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
986 | besttexp[iqsq]=texp[iqsq]; | |
987 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
62f44f07 | 988 | //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement |
989 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement | |
90beec49 | 990 | else bestchisquare[iqsq]=1000; |
5b4ed716 | 991 | } |
992 | ||
993 | Int_t npion=0; | |
994 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
995 | assparticle[j]=imass[j]; | |
996 | if(imass[j] == 0) npion++; | |
997 | } | |
998 | npionbest=npion; | |
999 | chisquarebest=chisquare; | |
1000 | t0best=meantzero; | |
1001 | eT0best=eMeanTzero; | |
1002 | } // close if(dummychisquare<=chisquare) | |
5b4ed716 | 1003 | } |
1004 | ||
2a258f40 | 1005 | Double_t chi2cut[nmaxtracksinsetMax]; |
5b4ed716 | 1006 | chi2cut[0] = 0; |
1007 | chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01 | |
1008 | for (Int_t j=2; j<ntracksinset; j++) { | |
1009 | chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.); | |
1010 | } | |
1011 | ||
1012 | Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy; | |
1013 | ||
90beec49 | 1014 | // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); |
5b4ed716 | 1015 | |
1016 | Bool_t kRedoT0 = kFALSE; | |
1017 | ntracksinsetmyCut = ntracksinsetmy; | |
2a258f40 | 1018 | Bool_t usetrack[nmaxtracksinsetMax]; |
5b4ed716 | 1019 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { |
1020 | usetrack[icsq] = kTRUE; | |
1021 | if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){ | |
90beec49 | 1022 | kRedoT0 = kTRUE; |
1023 | ntracksinsetmyCut--; | |
1024 | usetrack[icsq] = kFALSE; | |
5b4ed716 | 1025 | } |
1026 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) | |
1027 | ||
90beec49 | 1028 | // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut); |
2a258f40 | 1029 | |
5b4ed716 | 1030 | // Loop on mass hypotheses Redo |
1031 | if(kRedoT0 && ntracksinsetmyCut > 1){ | |
1032 | // printf("Redo T0\n"); | |
1033 | for (Int_t k=0; k < ncombinatorial;k++) { | |
1034 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
62f44f07 | 1035 | //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1)); |
1036 | imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1); | |
5b4ed716 | 1037 | texp[j]=exptof[j][imass[j]]; |
1038 | dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]); | |
62f44f07 | 1039 | //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100; |
1040 | if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100; | |
5b4ed716 | 1041 | } |
1042 | ||
1043 | Float_t sumAllweights=0.; | |
1044 | Float_t meantzero=0.; | |
1045 | Float_t eMeanTzero=0.; | |
1046 | ||
1047 | for (Int_t itz=0; itz<ntracksinsetmy;itz++) { | |
1048 | if(! usetrack[itz]) continue; | |
2a258f40 | 1049 | sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2 |
5b4ed716 | 1050 | |
1051 | timezero[itz]=texp[itz]-timeofflight[itz];// in ns | |
1052 | ||
1053 | weightedtimezero[itz]=timezero[itz]/sqTrackError[itz]; | |
1054 | sumAllweights+=1./sqTrackError[itz]; | |
1055 | meantzero+=weightedtimezero[itz]; | |
1056 | ||
1057 | } // end loop for (Int_t itz=0; itz<15;itz++) | |
1058 | ||
1059 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
1060 | eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns] | |
1061 | ||
1062 | // calculate chisquare | |
1063 | ||
1064 | Float_t chisquare=0.; | |
1065 | for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) { | |
1066 | if(! usetrack[icsq]) continue; | |
62f44f07 | 1067 | //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement |
1068 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement | |
90beec49 | 1069 | else chisquare+=1000; |
5b4ed716 | 1070 | } // end loop for (Int_t icsq=0; icsq<15;icsq++) |
1071 | ||
1072 | Int_t npion=0; | |
1073 | for (Int_t j=0; j<ntracksinsetmy; j++) { | |
1074 | assparticle[j]=imass[j]; | |
1075 | if(imass[j] == 0) npion++; | |
1076 | } | |
1077 | ||
90beec49 | 1078 | if(chisquare<=chisquarebest && npion>0){ |
5b4ed716 | 1079 | for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) { |
1080 | if(! usetrack[iqsq]) continue; | |
1081 | bestsqTrackError[iqsq]=sqTrackError[iqsq]; | |
1082 | besttimezero[iqsq]=timezero[iqsq]; | |
1083 | bestmomentum[iqsq]=momentum[iqsq]; | |
1084 | besttimeofflight[iqsq]=timeofflight[iqsq]; | |
1085 | besttexp[iqsq]=texp[iqsq]; | |
1086 | bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; | |
62f44f07 | 1087 | //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement |
1088 | if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement | |
90beec49 | 1089 | else bestchisquare[iqsq]=1000; |
5b4ed716 | 1090 | } |
1091 | ||
1092 | npionbest=npion; | |
1093 | chisquarebest=chisquare; | |
1094 | t0best=meantzero; | |
1095 | eT0best=eMeanTzero; | |
1096 | } // close if(dummychisquare<=chisquare) | |
1097 | ||
1098 | } | |
1099 | } | |
1100 | ||
1101 | // filling histos | |
1102 | Float_t confLevel=999; | |
1103 | ||
1104 | // Sets with decent chisquares | |
2a258f40 | 1105 | // printf("Chi2best of the set = %f \n",chisquarebest); |
5b4ed716 | 1106 | |
1107 | if(chisquarebest<999.){ | |
1108 | Double_t dblechisquare=(Double_t)chisquarebest; | |
1109 | confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1); | |
1110 | // cout << " Set Number " << nsets << endl; | |
1111 | // cout << "Best Assignment, selection " << assparticle[0] << | |
1112 | // assparticle[1] << assparticle[2] << | |
1113 | // assparticle[3] << assparticle[4] << | |
1114 | // assparticle[5] << endl; | |
1115 | // cout << " Chisquare of the set "<< chisquarebest <<endl; | |
1116 | // cout << " C.L. of the set "<< confLevel <<endl; | |
1117 | // cout << " T0 for this set (in ns) " << t0best << endl; | |
1118 | ||
90beec49 | 1119 | Int_t ntrackincurrentsel=0; |
5b4ed716 | 1120 | for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){ |
1121 | ||
1122 | if(! usetrack[icsq]) continue; | |
1123 | ||
1124 | // cout << "Track # " << icsq << " T0 offsets = " | |
1125 | // << besttimezero[icsq]-t0best << | |
1126 | // " track error = " << bestsqTrackError[icsq] | |
1127 | // << " Chisquare = " << bestchisquare[icsq] | |
1128 | // << " Momentum = " << bestmomentum[icsq] | |
1129 | // << " TOF = " << besttimeofflight[icsq] | |
1130 | // << " TOF tracking = " << besttexp[icsq] | |
1131 | // << " is used = " << usetrack[icsq] << endl; | |
90beec49 | 1132 | ntrackincurrentsel++; |
5b4ed716 | 1133 | } |
1134 | ||
2a258f40 | 1135 | // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel); |
1136 | ||
5b4ed716 | 1137 | // Pick up only those with C.L. >1% |
1138 | // if(confLevel>0.01 && ngoodsetsSel<200){ | |
1139 | if(confLevel>0.01 && ngoodsetsSel<200){ | |
1140 | chiSquarebestSel[ngoodsetsSel]=chisquarebest; | |
1141 | confLevelbestSel[ngoodsetsSel]=confLevel; | |
1142 | t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best; | |
1143 | eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best; | |
1144 | t0bestallSel += t0best/eT0best/eT0best; | |
1145 | sumWt0bestallSel += 1./eT0best/eT0best; | |
1146 | ngoodsetsSel++; | |
2a258f40 | 1147 | ngoodtrktrulyused+=ntracksinsetmyCut; |
90beec49 | 1148 | // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel); |
5b4ed716 | 1149 | } |
1150 | else{ | |
1151 | // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy); | |
1152 | } | |
1153 | } | |
62f44f07 | 1154 | //delete[] fTracksT0; |
1155 | fTracksT0->Clear(); | |
5b4ed716 | 1156 | nsets++; |
1157 | ||
1158 | } // end for the current set | |
1159 | ||
90beec49 | 1160 | //Redo the computation of the best in order to esclude very bad samples |
1161 | if(ngoodsetsSel > 1){ | |
1162 | Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel; | |
1163 | Int_t nsamples=ngoodsetsSel; | |
1164 | ngoodsetsSel=0; | |
1165 | t0bestallSel=0; | |
1166 | sumWt0bestallSel=0; | |
1167 | for (Int_t itz=0; itz<nsamples;itz++) { | |
1168 | if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){ | |
1169 | t0bestallSel += t0bestSel[itz]; | |
1170 | sumWt0bestallSel += eT0bestSel[itz]; | |
1171 | ngoodsetsSel++; | |
1172 | // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); | |
1173 | } | |
1174 | else{ | |
1175 | // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz])); | |
1176 | } | |
1177 | } | |
1178 | } | |
1179 | if(ngoodsetsSel < 1){ | |
1180 | sumWt0bestallSel = 0.0; | |
1181 | } | |
1182 | //--------------------------------End recomputation | |
1183 | ||
5b4ed716 | 1184 | nUsedTracks = ngoodtrkt0; |
1185 | if(strstr(option,"all")){ | |
1186 | if(sumAllweightspi>0.){ | |
1187 | meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns] | |
1188 | eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns] | |
1189 | } | |
1190 | ||
2a258f40 | 1191 | // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel); |
1192 | ||
5b4ed716 | 1193 | if(sumWt0bestallSel>0){ |
1194 | t0bestallSel = t0bestallSel/sumWt0bestallSel; | |
1195 | eT0bestallSel = sqrt(1./sumWt0bestallSel); | |
2a258f40 | 1196 | // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel); |
5b4ed716 | 1197 | }// end of if(sumWt0bestallSel>0){ |
1198 | ||
1199 | // cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl; | |
1200 | } | |
1201 | ||
1202 | t0def=t0bestallSel; | |
1203 | deltat0def=eT0bestallSel; | |
5b4ed716 | 1204 | |
1205 | fT0SigmaT0def[0]=t0def; | |
2a258f40 | 1206 | fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1)); |
5b4ed716 | 1207 | fT0SigmaT0def[2]=ngoodtrkt0; |
1208 | fT0SigmaT0def[3]=ngoodtrktrulyused; | |
1209 | } | |
1210 | } | |
62f44f07 | 1211 | |
1212 | //delete [] fGTracks; | |
1213 | fGTracks->Clear(); | |
1214 | ||
5b4ed716 | 1215 | // if(strstr(option,"tim") || strstr(option,"all")){ |
1216 | // cout << "AliTOFT0v1:" << endl ; | |
1217 | //} | |
1218 | ||
2a258f40 | 1219 | if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6; |
1220 | ||
1221 | bench->Stop("t0computation"); | |
1222 | ||
1223 | fT0SigmaT0def[4]=bench->GetRealTime("t0computation"); | |
1224 | fT0SigmaT0def[5]=bench->GetCpuTime("t0computation"); | |
1225 | ||
1226 | // bench->Print("t0computation"); | |
1227 | // 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])); | |
1228 | ||
90beec49 | 1229 | delete bench; |
1230 | bench=NULL; | |
5b4ed716 | 1231 | |
536031f2 | 1232 | return fT0SigmaT0def; |
1233 | } | |
1234 | //__________________________________________________________________ | |
8f589502 | 1235 | Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const |
536031f2 | 1236 | { |
8f589502 | 1237 | // Take the error extimate for the TOF time in the track reconstruction |
536031f2 | 1238 | |
536031f2 | 1239 | static const Double_t kMasses[]={ |
1240 | 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 | |
1241 | }; | |
1242 | ||
1243 | Double_t mass=kMasses[index+2]; | |
536031f2 | 1244 | |
2a258f40 | 1245 | Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass); |
536031f2 | 1246 | |
1247 | return sigma; | |
1248 | } | |
14b2cbea | 1249 | |
1250 | //__________________________________________________________________ | |
1251 | Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track) | |
1252 | { | |
1253 | ||
1254 | /* TPC refit */ | |
1255 | if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE; | |
1256 | /* do not accept kink daughters */ | |
1257 | if (track->GetKinkIndex(0)>0) return kFALSE; | |
1258 | /* N clusters TPC */ | |
1259 | if (track->GetTPCclusters(0) < 50) return kFALSE; | |
1260 | /* chi2 TPC */ | |
1261 | if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE; | |
1262 | /* sigma to vertex */ | |
1263 | if (GetSigmaToVertex(track) > 4.) return kFALSE; | |
1264 | ||
1265 | /* accept track */ | |
1266 | return kTRUE; | |
1267 | ||
1268 | } | |
1269 | ||
1270 | //____________________________________________________________________ | |
8f589502 | 1271 | Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const |
14b2cbea | 1272 | { |
1273 | // Calculates the number of sigma to the vertex. | |
1274 | ||
1275 | Float_t b[2]; | |
1276 | Float_t bRes[2]; | |
1277 | Float_t bCov[3]; | |
1278 | esdTrack->GetImpactParameters(b,bCov); | |
1279 | ||
1280 | if (bCov[0]<=0 || bCov[2]<=0) { | |
1281 | bCov[0]=0; bCov[2]=0; | |
1282 | } | |
1283 | bRes[0] = TMath::Sqrt(bCov[0]); | |
1284 | bRes[1] = TMath::Sqrt(bCov[2]); | |
1285 | ||
1286 | // ----------------------------------- | |
1287 | // How to get to a n-sigma cut? | |
1288 | // | |
1289 | // The accumulated statistics from 0 to d is | |
1290 | // | |
1291 | // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma) | |
1292 | // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma) | |
1293 | // | |
1294 | // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2) | |
1295 | // Can this be expressed in a different way? | |
1296 | ||
1297 | if (bRes[0] == 0 || bRes[1] ==0) | |
1298 | return -1; | |
1299 | ||
62f44f07 | 1300 | //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); |
1301 | Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2)); | |
14b2cbea | 1302 | |
1303 | // work around precision problem | |
1304 | // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :( | |
1305 | // 1e-15 corresponds to nsigma ~ 7.7 | |
1306 | if (TMath::Exp(-d * d / 2) < 1e-15) | |
1307 | return 1000; | |
1308 | ||
1309 | Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); | |
1310 | return nSigma; | |
1311 | } | |
90beec49 | 1312 | //____________________________________________________________________ |
1313 | ||
1314 | Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{ | |
1315 | Bool_t status = kFALSE; | |
1316 | ||
1317 | Double_t exptimes[5]; | |
1318 | track->GetIntegratedTimes(exptimes); | |
1319 | ||
1320 | Float_t dedx = track->GetTPCsignal(); | |
1321 | ||
1322 | Double_t ptpc[3]; | |
1323 | track->GetInnerPxPyPz(ptpc); | |
1324 | Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]); | |
1325 | ||
1326 | if(imass > 2 || imass < 0) return status; | |
1327 | Int_t i = imass+2; | |
1328 | ||
1329 | AliPID::EParticleType type=AliPID::EParticleType(i); | |
1330 | ||
1331 | Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type); | |
1332 | Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type); | |
1333 | ||
1334 | if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){ | |
1335 | status = kTRUE; | |
1336 | } | |
1337 | ||
1338 | return status; | |
1339 | } | |
62f44f07 | 1340 | |
1341 | //____________________________________________________________________ | |
1342 | Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const | |
1343 | { | |
1344 | // | |
1345 | // Returns base^exponent | |
1346 | // | |
1347 | ||
1348 | Float_t power=1.; | |
1349 | ||
1350 | for (Int_t ii=exponent; ii>0; ii--) | |
1351 | power=power*base; | |
1352 | ||
1353 | return power; | |
1354 | ||
1355 | } | |
1356 | //____________________________________________________________________ | |
1357 | Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const | |
1358 | { | |
1359 | // | |
1360 | // Returns base^exponent | |
1361 | // | |
1362 | ||
1363 | Int_t power=1; | |
1364 | ||
1365 | for (Int_t ii=exponent; ii>0; ii--) | |
1366 | power=power*base; | |
1367 | ||
1368 | return power; | |
1369 | ||
1370 | } |