Line data Source code
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 :
16 : /* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
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.
22 : // We select in the ESD some "primary" particles - or tracks in the following -
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
25 : // of the track when the TOF is reached,
26 : // the momentum and the time of flight
27 : // given by the TOF detector.
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
55 : //-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
56 : //-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table
57 : // for Power3)
58 : //////////////////////////////////////////////////////////////////////////////
59 :
60 : #include "AliESDtrack.h"
61 : #include "AliESDEvent.h"
62 : #include "AliTOFT0v1.h"
63 : #include "TBenchmark.h"
64 : #include "AliPID.h"
65 : #include "AliESDpid.h"
66 :
67 26 : ClassImp(AliTOFT0v1)
68 :
69 : //____________________________________________________________________________
70 : AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
71 8 : TObject(),
72 8 : fLowerMomBound(0.5),
73 8 : fUpperMomBound(3),
74 8 : fTimeCorr(0.),
75 8 : fEvent(0x0),
76 8 : fPIDesd(extPID),
77 24 : fTracks(new TObjArray(10)),
78 24 : fGTracks(new TObjArray(10)),
79 24 : fTracksT0(new TObjArray(10)),
80 8 : fOptFlag(kFALSE)
81 40 : {
82 : //
83 : // default constructor
84 : //
85 16 : if(AliPID::ParticleMass(0) == 0) new AliPID();
86 :
87 8 : if(!fPIDesd){
88 0 : fPIDesd = new AliESDpid();
89 0 : }
90 :
91 8 : Init(NULL);
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 256 : for (Int_t i=0; i<15; ++i) {
97 120 : fLookupPowerThree[i]=ToCalculatePower(3,i);
98 : }
99 16 : }
100 :
101 : //____________________________________________________________________________
102 : AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
103 0 : TObject(),
104 0 : fLowerMomBound(0.5),
105 0 : fUpperMomBound(3.0),
106 0 : fTimeCorr(0.),
107 0 : fEvent(event),
108 0 : fPIDesd(extPID),
109 0 : fTracks(new TObjArray(10)),
110 0 : fGTracks(new TObjArray(10)),
111 0 : fTracksT0(new TObjArray(10)),
112 0 : fOptFlag(kFALSE)
113 0 : {
114 : //
115 : // real constructor
116 : //
117 0 : if(AliPID::ParticleMass(0) == 0) new AliPID();
118 :
119 0 : if(!fPIDesd){
120 0 : fPIDesd = new AliESDpid();
121 0 : }
122 :
123 0 : Init(event);
124 : //initialise lookup table for power 3
125 0 : for (Int_t i=0; i<15; ++i) {
126 0 : fLookupPowerThree[i]=Int_t(TMath::Power(3,i));
127 : }
128 0 : }
129 : //____________________________________________________________________________
130 : AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
131 : {
132 : //
133 : // assign. operator
134 : //
135 :
136 0 : if (this == &tzero)
137 0 : return *this;
138 :
139 0 : fLowerMomBound=tzero.fLowerMomBound;
140 0 : fUpperMomBound=tzero.fUpperMomBound;
141 0 : fTimeCorr=tzero.fTimeCorr;
142 0 : fEvent=tzero.fEvent;
143 0 : fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
144 0 : fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
145 0 : fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
146 0 : fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
147 :
148 0 : fTracks=tzero.fTracks;
149 0 : fGTracks=tzero.fGTracks;
150 0 : fTracksT0=tzero.fTracksT0;
151 :
152 0 : for (Int_t ii=0; ii<tzero.fTracks->GetEntries(); ii++)
153 0 : fTracks->AddLast(tzero.fTracks->At(ii));
154 :
155 0 : for (Int_t ii=0; ii<tzero.fGTracks->GetEntries(); ii++)
156 0 : fGTracks->AddLast(tzero.fGTracks->At(ii));
157 :
158 0 : for (Int_t ii=0; ii<tzero.fTracksT0->GetEntries(); ii++)
159 0 : fTracksT0->AddLast(tzero.fTracksT0->At(ii));
160 :
161 0 : fOptFlag=tzero.fOptFlag;
162 :
163 0 : return *this;
164 0 : }
165 :
166 : //____________________________________________________________________________
167 : AliTOFT0v1::~AliTOFT0v1()
168 48 : {
169 : // dtor
170 8 : fEvent=NULL;
171 :
172 8 : if (fTracks) {
173 8 : fTracks->Clear();
174 16 : delete fTracks;
175 8 : fTracks=0x0;
176 8 : }
177 :
178 8 : if (fGTracks) {
179 8 : fGTracks->Clear();
180 16 : delete fGTracks;
181 8 : fGTracks=0x0;
182 8 : }
183 :
184 8 : if (fTracksT0) {
185 8 : fTracksT0->Clear();
186 16 : delete fTracksT0;
187 8 : fTracksT0=0x0;
188 8 : }
189 :
190 :
191 24 : }
192 : //____________________________________________________________________________
193 :
194 : void
195 : AliTOFT0v1::Init(AliESDEvent *event)
196 : {
197 :
198 : /*
199 : * init
200 : */
201 :
202 32 : fEvent = event;
203 16 : fT0SigmaT0def[0]=0.;
204 16 : fT0SigmaT0def[1]=0.6;
205 16 : fT0SigmaT0def[2]=0.;
206 16 : fT0SigmaT0def[3]=0.;
207 :
208 16 : }
209 : //____________________________________________________________________________
210 : Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
211 : {
212 176 : TBenchmark *bench=new TBenchmark();
213 88 : bench->Start("t0computation");
214 :
215 : // Caluclate the Event Time using the ESD TOF time
216 :
217 88 : fT0SigmaT0def[0]=0.;
218 88 : fT0SigmaT0def[1]=0.600;
219 88 : fT0SigmaT0def[2]=0.;
220 88 : fT0SigmaT0def[3]=0.;
221 :
222 : const Int_t nmaxtracksinsetMax=10;
223 : Int_t nmaxtracksinset=10;
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 88 : Float_t t0bestSel[300];
234 88 : Float_t eT0bestSel[300];
235 88 : Float_t chiSquarebestSel[300];
236 88 : 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 88 : Int_t ntrk=fEvent->GetNumberOfTracks();
249 :
250 : Int_t ngoodtrk=0;
251 88 : Int_t ngoodtrkt0 =0;
252 : Float_t meantime =0;
253 :
254 : // First Track loop, Selection of good tracks
255 :
256 88 : fTracks->Clear();
257 3520 : for (Int_t itrk=0; itrk<ntrk; itrk++) {
258 1672 : AliESDtrack *t=fEvent->GetTrack(itrk);
259 1672 : Double_t momOld=t->GetP();
260 1672 : Double_t mom=momOld-0.0036*momOld;
261 1672 : if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
262 2453 : if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
263 891 : Double_t time=t->GetTOFsignal();
264 :
265 891 : time*=1.E-3; // tof given in nanoseconds
266 1782 : if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
267 :
268 737 : if (!AcceptTrack(t)) continue;
269 :
270 693 : if(t->GetIntegratedLength() < 350)continue; //skip decays
271 1237 : if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
272 :
273 611 : meantime+=time;
274 611 : fTracks->AddLast(t);
275 611 : ngoodtrk++;
276 611 : }
277 :
278 176 : if(ngoodtrk > 1) meantime /= ngoodtrk;
279 :
280 88 : if(ngoodtrk>22) nmaxtracksinset = 6;
281 :
282 88 : fGTracks->Clear();
283 1398 : for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
284 611 : AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
285 : // Double_t time=t->GetTOFsignal();
286 : // if((time-meantime*1.E3)<50.E3){ // For pp and per
287 611 : fGTracks->AddLast(t);
288 611 : ngoodtrkt0++;
289 : // }
290 : }
291 :
292 88 : fTracks->Clear();
293 :
294 88 : Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
295 88 : Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
296 88 : if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
297 :
298 :
299 88 : if(ngoodtrkt0<2){
300 : t0def=-999;
301 : deltat0def=0.600;
302 0 : fT0SigmaT0def[0]=t0def;
303 0 : fT0SigmaT0def[1]=deltat0def;
304 0 : fT0SigmaT0def[2]=ngoodtrkt0;
305 0 : fT0SigmaT0def[3]=ngoodtrkt0;
306 : //goto finish;
307 0 : }
308 88 : if(ngoodtrkt0>=2){
309 : // Decide how many tracks in set
310 88 : Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
311 : Int_t nset=1;
312 :
313 88 : if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
314 :
315 : // Loop over selected sets
316 :
317 88 : if(nset>=1){
318 440 : for (Int_t i=0; i< nset; i++) {
319 : // printf("Set %i of %i\n",i+1,nset);
320 : Float_t t0best=999.;
321 : Float_t eT0best=999.;
322 : Float_t chisquarebest=99999.;
323 : Int_t npionbest=0;
324 :
325 88 : fTracksT0->Clear();
326 : Int_t ntracksinsetmy=0;
327 1398 : for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
328 611 : Int_t index = itrk+i*ntracksinset;
329 611 : if(index < fGTracks->GetEntries()){
330 611 : AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
331 611 : fTracksT0->AddLast(t);
332 611 : ntracksinsetmy++;
333 611 : }
334 : }
335 :
336 : // Analyse it
337 :
338 88 : Int_t assparticle[nmaxtracksinsetMax];
339 88 : Float_t exptof[nmaxtracksinsetMax][3];
340 88 : Float_t momErr[nmaxtracksinsetMax][3];
341 88 : Float_t timeofflight[nmaxtracksinsetMax];
342 88 : Float_t momentum[nmaxtracksinsetMax];
343 88 : Float_t timezero[nmaxtracksinsetMax];
344 88 : Float_t weightedtimezero[nmaxtracksinsetMax];
345 88 : Float_t beta[nmaxtracksinsetMax];
346 88 : Float_t texp[nmaxtracksinsetMax];
347 88 : Float_t dtexp[nmaxtracksinsetMax];
348 88 : Float_t sqMomError[nmaxtracksinsetMax];
349 88 : Float_t sqTrackError[nmaxtracksinsetMax];
350 88 : Float_t massarray[3]={0.13957,0.493677,0.9382723};
351 88 : Float_t tracktoflen[nmaxtracksinsetMax];
352 88 : Float_t besttimezero[nmaxtracksinsetMax];
353 88 : Float_t besttexp[nmaxtracksinsetMax];
354 88 : Float_t besttimeofflight[nmaxtracksinsetMax];
355 88 : Float_t bestmomentum[nmaxtracksinsetMax];
356 88 : Float_t bestchisquare[nmaxtracksinsetMax];
357 88 : Float_t bestweightedtimezero[nmaxtracksinsetMax];
358 88 : Float_t bestsqTrackError[nmaxtracksinsetMax];
359 88 : Int_t imass[nmaxtracksinsetMax];
360 :
361 1398 : for (Int_t j=0; j<ntracksinset; j++) {
362 611 : assparticle[j] = 3;
363 611 : timeofflight[j] = 0;
364 611 : momentum[j] = 0;
365 611 : timezero[j] = 0;
366 611 : weightedtimezero[j] = 0;
367 611 : beta[j] = 0;
368 611 : texp[j] = 0;
369 611 : dtexp[j] = 0;
370 611 : sqMomError[j] = 0;
371 611 : sqTrackError[j] = 0;
372 611 : tracktoflen[j] = 0;
373 611 : besttimezero[j] = 0;
374 611 : besttexp[j] = 0;
375 611 : besttimeofflight[j] = 0;
376 611 : bestmomentum[j] = 0;
377 611 : bestchisquare[j] = 0;
378 611 : bestweightedtimezero[j] = 0;
379 611 : bestsqTrackError[j] = 0;
380 611 : imass[j] = 1;
381 : }
382 :
383 1398 : for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
384 611 : AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
385 611 : Double_t momOld=t->GetP();
386 611 : Double_t mom=momOld-0.0036*momOld;
387 611 : Double_t time=t->GetTOFsignal();
388 :
389 611 : time*=1.E-3; // tof given in nanoseconds
390 611 : Double_t exptime[AliPID::kSPECIESC];
391 611 : t->GetIntegratedTimes(exptime,AliPID::kSPECIESC);
392 611 : Double_t toflen=t->GetIntegratedLength();
393 611 : toflen=toflen/100.; // toflen given in m
394 :
395 611 : timeofflight[j]=time;
396 611 : tracktoflen[j]=toflen;
397 611 : exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
398 611 : exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
399 611 : exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
400 611 : momentum[j]=mom;
401 611 : assparticle[j]=3;
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 611 : if (fOptFlag) {
406 0 : momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]);
407 0 : momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]);
408 0 : momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]);
409 0 : }
410 611 : } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
411 :
412 1398 : for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
413 1222 : beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
414 611 : +momentum[itz]*momentum[itz]);
415 611 : 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]));
416 611 : sqTrackError[itz]=sqMomError[itz]; //in ns
417 611 : timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
418 611 : weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
419 611 : sumAllweightspi+=1./sqTrackError[itz];
420 611 : meantzeropi+=weightedtimezero[itz];
421 : } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
422 :
423 :
424 : // Then, Combinatorial Algorithm
425 :
426 88 : if(ntracksinsetmy<2 )break;
427 :
428 1398 : for (Int_t j=0; j<ntracksinsetmy; j++) {
429 611 : imass[j] = 3;
430 : }
431 :
432 : Int_t ncombinatorial;
433 88 : if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
434 88 : else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
435 :
436 :
437 : // Loop on mass hypotheses
438 1488956 : for (Int_t k=0; k < ncombinatorial;k++) {
439 15314670 : for (Int_t j=0; j<ntracksinsetmy; j++) {
440 6912945 : imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
441 6912945 : texp[j]=exptof[j][imass[j]];
442 6912945 : if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
443 6912945 : else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
444 : }
445 :
446 : Float_t sumAllweights=0.;
447 : Float_t meantzero=0.;
448 : Float_t eMeanTzero=0.;
449 :
450 15314670 : for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
451 6912945 : sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
452 :
453 6912945 : timezero[itz]=texp[itz]-timeofflight[itz];// in ns
454 :
455 6912945 : weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
456 6912945 : sumAllweights+=1./sqTrackError[itz];
457 6912945 : meantzero+=weightedtimezero[itz];
458 :
459 : } // end loop for (Int_t itz=0; itz<15;itz++)
460 :
461 744390 : meantzero=meantzero/sumAllweights; // it is given in [ns]
462 744390 : eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
463 :
464 : // calculate chisquare
465 : Float_t chisquare=0.;
466 15314670 : for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
467 6912945 : chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
468 : } // end loop for (Int_t icsq=0; icsq<15;icsq++)
469 :
470 744390 : if(chisquare<=chisquarebest){
471 13580 : for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
472 :
473 5984 : bestsqTrackError[iqsq]=sqTrackError[iqsq];
474 5984 : besttimezero[iqsq]=timezero[iqsq];
475 5984 : bestmomentum[iqsq]=momentum[iqsq];
476 5984 : besttimeofflight[iqsq]=timeofflight[iqsq];
477 5984 : besttexp[iqsq]=texp[iqsq];
478 5984 : bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
479 5984 : bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
480 : }
481 :
482 : Int_t npion=0;
483 13580 : for (Int_t j=0; j<ntracksinsetmy; j++) {
484 5984 : assparticle[j]=imass[j];
485 10012 : if(imass[j] == 0) npion++;
486 : }
487 : npionbest=npion;
488 : chisquarebest=chisquare;
489 : t0best=meantzero;
490 : eT0best=eMeanTzero;
491 806 : } // close if(dummychisquare<=chisquare)
492 : }
493 :
494 88 : Double_t chi2cut[nmaxtracksinsetMax];
495 88 : chi2cut[0] = 0;
496 88 : chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
497 1046 : for (Int_t j=2; j<ntracksinset; j++) {
498 435 : chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
499 : }
500 :
501 88 : Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
502 :
503 : // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
504 :
505 : Bool_t kRedoT0 = kFALSE;
506 : ntracksinsetmyCut = ntracksinsetmy;
507 88 : Bool_t usetrack[nmaxtracksinsetMax];
508 1398 : for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
509 611 : usetrack[icsq] = kTRUE;
510 1213 : if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
511 : kRedoT0 = kTRUE;
512 96 : ntracksinsetmyCut--;
513 96 : usetrack[icsq] = kFALSE;
514 : // printf("tracks chi2 = %f\n",bestchisquare[icsq]);
515 96 : }
516 : } // end loop for (Int_t icsq=0; icsq<15;icsq++)
517 :
518 : // Loop on mass hypotheses Redo
519 88 : if(kRedoT0 && ntracksinsetmyCut > 1){
520 : // printf("Redo T0\n");
521 760020 : for (Int_t k=0; k < ncombinatorial;k++) {
522 7815042 : for (Int_t j=0; j<ntracksinsetmy; j++) {
523 3527550 : imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
524 3527550 : texp[j]=exptof[j][imass[j]];
525 3527550 : if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
526 3527550 : else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
527 : }
528 :
529 : Float_t sumAllweights=0.;
530 : Float_t meantzero=0.;
531 : Float_t eMeanTzero=0.;
532 :
533 7815042 : for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
534 3527550 : if(! usetrack[itz]) continue;
535 2162700 : sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
536 :
537 2162700 : timezero[itz]=texp[itz]-timeofflight[itz];// in ns
538 :
539 2162700 : weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
540 2162700 : sumAllweights+=1./sqTrackError[itz];
541 2162700 : meantzero+=weightedtimezero[itz];
542 :
543 2162700 : } // end loop for (Int_t itz=0; itz<15;itz++)
544 :
545 379971 : meantzero=meantzero/sumAllweights; // it is given in [ns]
546 379971 : eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
547 :
548 : // calculate chisquare
549 :
550 : Float_t chisquare=0.;
551 7815042 : for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
552 3527550 : if(! usetrack[icsq]) continue;
553 2162700 : chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
554 2162700 : } // end loop for (Int_t icsq=0; icsq<15;icsq++)
555 :
556 : Int_t npion=0;
557 7815042 : for (Int_t j=0; j<ntracksinsetmy; j++) {
558 3527550 : assparticle[j]=imass[j];
559 4703400 : if(imass[j] == 0) npion++;
560 : }
561 :
562 379971 : if(chisquare<=chisquarebest && npion>0){
563 18846 : for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
564 8412 : if(! usetrack[iqsq]) continue;
565 5052 : bestsqTrackError[iqsq]=sqTrackError[iqsq];
566 5052 : besttimezero[iqsq]=timezero[iqsq];
567 5052 : bestmomentum[iqsq]=momentum[iqsq];
568 5052 : besttimeofflight[iqsq]=timeofflight[iqsq];
569 5052 : besttexp[iqsq]=texp[iqsq];
570 5052 : bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
571 5052 : bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
572 5052 : }
573 :
574 : npionbest=npion;
575 : chisquarebest=chisquare;
576 : t0best=meantzero;
577 : eT0best=eMeanTzero;
578 1011 : } // close if(dummychisquare<=chisquare)
579 :
580 : }
581 39 : }
582 :
583 : // filling histos
584 : Float_t confLevel=999;
585 :
586 : // Sets with decent chisquares
587 : // printf("Chi2best of the set = %f \n",chisquarebest);
588 :
589 88 : if(chisquarebest<999.){
590 : Double_t dblechisquare=(Double_t)chisquarebest;
591 88 : confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
592 :
593 : Int_t ntrackincurrentsel=0;
594 1398 : for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
595 :
596 611 : if(! usetrack[icsq]) continue;
597 :
598 515 : ntrackincurrentsel++;
599 515 : }
600 :
601 : // printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
602 :
603 : // Pick up only those with C.L. >1%
604 88 : if(confLevel>0.01 && ngoodsetsSel<200){
605 88 : chiSquarebestSel[ngoodsetsSel]=chisquarebest;
606 88 : confLevelbestSel[ngoodsetsSel]=confLevel;
607 88 : t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
608 88 : eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
609 88 : t0bestallSel += t0best/eT0best/eT0best;
610 88 : sumWt0bestallSel += 1./eT0best/eT0best;
611 88 : ngoodsetsSel++;
612 88 : ngoodtrktrulyused+=ntracksinsetmyCut;
613 : // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
614 88 : }
615 : else{
616 : // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
617 : }
618 88 : }
619 88 : fTracksT0->Clear();
620 88 : nsets++;
621 :
622 176 : } // end for the current set
623 :
624 : //Redo the computation of the best in order to esclude very bad samples
625 88 : if(ngoodsetsSel > 1){
626 0 : Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
627 : Int_t nsamples=ngoodsetsSel;
628 : ngoodsetsSel=0;
629 : t0bestallSel=0;
630 : sumWt0bestallSel=0;
631 0 : for (Int_t itz=0; itz<nsamples;itz++) {
632 0 : if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
633 0 : t0bestallSel += t0bestSel[itz];
634 0 : sumWt0bestallSel += eT0bestSel[itz];
635 0 : ngoodsetsSel++;
636 : // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
637 0 : }
638 : else{
639 : // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
640 : }
641 : }
642 0 : }
643 88 : if(ngoodsetsSel < 1){
644 : sumWt0bestallSel = 0.0;
645 0 : }
646 : //--------------------------------End recomputation
647 :
648 : nUsedTracks = ngoodtrkt0;
649 88 : if(strstr(option,"all")){
650 88 : if(sumAllweightspi>0.){
651 : meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
652 88 : eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
653 88 : }
654 :
655 : // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
656 :
657 88 : if(sumWt0bestallSel>0){
658 88 : t0bestallSel = t0bestallSel/sumWt0bestallSel;
659 88 : eT0bestallSel = sqrt(1./sumWt0bestallSel);
660 : // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
661 88 : }// end of if(sumWt0bestallSel>0){
662 :
663 : }
664 :
665 88 : t0def=t0bestallSel;
666 88 : deltat0def=eT0bestallSel;
667 :
668 88 : fT0SigmaT0def[0]=t0def;
669 88 : fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
670 88 : fT0SigmaT0def[2]=ngoodtrkt0;
671 88 : fT0SigmaT0def[3]=ngoodtrktrulyused;
672 88 : }
673 88 : }
674 :
675 88 : fGTracks->Clear();
676 :
677 88 : if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
678 :
679 88 : bench->Stop("t0computation");
680 :
681 88 : fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
682 88 : 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 :
687 176 : delete bench;
688 : bench=NULL;
689 :
690 176 : return fT0SigmaT0def;
691 88 : }
692 : //__________________________________________________________________
693 : Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
694 : {
695 : // Take the error extimate for the TOF time in the track reconstruction
696 :
697 : static const Double_t kMasses[]={
698 : 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
699 : };
700 :
701 20880990 : Double_t mass=kMasses[index+2];
702 :
703 10440495 : Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
704 :
705 10440495 : return sigma;
706 : }
707 :
708 : //__________________________________________________________________
709 : Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
710 : {
711 :
712 : /* TPC refit */
713 1430 : if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
714 : /* do not accept kink daughters */
715 715 : if (track->GetKinkIndex(0)>0) return kFALSE;
716 : /* N clusters TPC */
717 715 : if (track->GetTPCclusters(0) < 50) return kFALSE;
718 : /* chi2 TPC */
719 715 : if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
720 : /* sigma to vertex */
721 737 : if (GetSigmaToVertex(track) > 4.) return kFALSE;
722 :
723 : /* accept track */
724 693 : return kTRUE;
725 :
726 715 : }
727 :
728 : //____________________________________________________________________
729 : Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
730 : {
731 : // Calculates the number of sigma to the vertex.
732 :
733 1430 : Float_t b[2];
734 : Float_t bRes[2];
735 715 : Float_t bCov[3];
736 715 : esdTrack->GetImpactParameters(b,bCov);
737 :
738 1430 : if (bCov[0]<=0 || bCov[2]<=0) {
739 0 : bCov[0]=0; bCov[2]=0;
740 0 : }
741 715 : bRes[0] = TMath::Sqrt(bCov[0]);
742 715 : 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 1430 : if (bRes[0] == 0 || bRes[1] ==0)
756 0 : return -1;
757 :
758 : //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
759 715 : Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
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 715 : if (TMath::Exp(-d * d / 2) < 1e-15)
765 22 : return 1000;
766 :
767 693 : Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
768 : return nSigma;
769 715 : }
770 : //____________________________________________________________________
771 :
772 : Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
773 : Bool_t status = kFALSE;
774 :
775 0 : Double_t exptimes[AliPID::kSPECIESC];
776 0 : track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC);
777 :
778 0 : Float_t dedx = track->GetTPCsignal();
779 :
780 0 : Double_t ptpc[3];
781 0 : track->GetInnerPxPyPz(ptpc);
782 0 : Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
783 :
784 0 : if(imass > 2 || imass < 0) return status;
785 0 : Int_t i = imass+2;
786 :
787 : AliPID::EParticleType type=AliPID::EParticleType(i);
788 :
789 0 : Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
790 0 : Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
791 :
792 0 : if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
793 : status = kTRUE;
794 0 : }
795 :
796 0 : return status;
797 0 : }
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 10010 : for (Int_t ii=exponent; ii>0; ii--)
809 2860 : power=power*base;
810 :
811 1430 : 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 3526 : for (Int_t ii=exponent; ii>0; ii--)
824 1451 : power=power*base;
825 :
826 208 : return power;
827 :
828 : }
|
