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