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