TOF/AliTOFT0v1.cxx

   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 //////////////////////////////////////////////////////////////////////////////
  57
  58 #include "AliESDtrack.h"
  59 #include "AliESDEvent.h"
  60 #include "AliTOFT0v1.h"
  61
  62 ClassImp(AliTOFT0v1)
  63
  64 //____________________________________________________________________________
  65 AliTOFT0v1::AliTOFT0v1():
  66   fLowerMomBound(0.5),
  67   fUpperMomBound(1.5),
  68   fTimeResolution(0.80e-10),
  69   fTimeCorr(0.),
  70   fEvent(0x0)
  71 //   fCalib(0x0)
  72 {
  73   //
  74   // default constructor
  75   //
  76
  77   fT0SigmaT0def[0]=-999.;
  78   fT0SigmaT0def[1]=999.;
  79   fT0SigmaT0def[2]=-999.;
  80   fT0SigmaT0def[3]=-999.;
  81
  82 }
  83
  84
  85 //____________________________________________________________________________
  86 AliTOFT0v1::AliTOFT0v1(AliESDEvent* event):
  87   fLowerMomBound(0.5),
  88   fUpperMomBound(1.5),
  89   fTimeResolution(0.80e-10),
  90   fTimeCorr(0.),
  91   fEvent(event)
  92 //   fCalib(0x0)
  93 {
  94   //
  95   // real constructor
  96   //
  97
  98   fT0SigmaT0def[0]=-999.;
  99   fT0SigmaT0def[1]= 999.;
 100   fT0SigmaT0def[2]=-999.;
 101   fT0SigmaT0def[3]=-999.;
 102
 103 }
 104
 105 //____________________________________________________________________________
 106 AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero):
 107   TObject(),
 108   fLowerMomBound(tzero.fLowerMomBound),
 109   fUpperMomBound(tzero.fUpperMomBound),
 110   fTimeResolution(tzero.fTimeResolution),
 111   fTimeCorr(tzero.fTimeCorr),
 112   fEvent(tzero.fEvent)
 113 //   fCalib(tzero.fCalib)
 114 {
 115   //
 116   // copy constructor
 117   //
 118
 119   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
 120   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
 121   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
 122   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
 123
 124 }
 125
 126 //____________________________________________________________________________
 127 AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
 128 {
 129  //
 130   // assign. operator
 131   //
 132
 133   if (this == &tzero)
 134     return *this;
 135
 136   fLowerMomBound=tzero.fLowerMomBound;
 137   fUpperMomBound=tzero.fUpperMomBound;
 138   fTimeResolution=tzero.fTimeResolution;
 139   fTimeCorr=tzero.fTimeCorr;
 140   fEvent=tzero.fEvent;
 141 //   fCalib=tzero.fCalib;
 142   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
 143   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
 144   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
 145   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
 146
 147   return *this;
 148 }
 149 //____________________________________________________________________________
 150 AliTOFT0v1::~AliTOFT0v1()
 151 {
 152   // dtor
 153 //   fCalib=NULL;
 154   fEvent=NULL;
 155
 156 }
 157 //____________________________________________________________________________
 158 void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
 159   // Set the TOF time resolution
 160   fTimeResolution=timeresolution;
 161 }
 162 //____________________________________________________________________________
 163 //____________________________________________________________________________
 164 Double_t * AliTOFT0v1::DefineT0(Option_t *option)
 165 {
 166   // Caluclate the Event Time using the ESD TOF time
 167
 168  Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
 169
 170   const Int_t nmaxtracksinset=10;
 171 //   if(strstr(option,"all")){
 172 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
 173 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
 174 //   }
 175
 176
 177   Int_t nsets=0;
 178   Int_t nUsedTracks=0;
 179   Int_t ngoodsetsSel= 0;
 180   Float_t t0bestSel[300];
 181   Float_t eT0bestSel[300];
 182   Float_t chiSquarebestSel[300];
 183   Float_t confLevelbestSel[300];
 184   Float_t t0bestallSel=0.;
 185   Float_t eT0bestallSel=0.;
 186   Float_t sumWt0bestallSel=0.;
 187   Float_t eMeanTzeroPi=0.;
 188   Float_t meantzeropi=0.;
 189   Float_t sumAllweightspi=0.;
 190   Double_t t0def=-999;
 191   Double_t deltat0def=999;
 192   Int_t ngoodtrktrulyused=0;
 193   Int_t ntracksinsetmyCut = 0;
 194
 195   Int_t ntrk=fEvent->GetNumberOfTracks();
 196
 197   AliESDtrack **tracks=new AliESDtrack*[ntrk];
 198   Int_t ngoodtrk=0;
 199   Int_t ngoodtrkt0 =0;
 200   Float_t mintime =1E6;
 201
 202   // First Track loop, Selection of good tracks
 203
 204   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 205     AliESDtrack *t=fEvent->GetTrack(itrk);
 206     Double_t momOld=t->GetP();
 207     Double_t mom=momOld-0.0036*momOld;
 208     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 209     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 210     Double_t time=t->GetTOFsignal();
 211
 212     time*=1.E-3; // tof given in nanoseconds
 213     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
 214
 215     if (!AcceptTrack(t)) continue;
 216
 217     if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
 218     if(time <= mintime) mintime=time;
 219     tracks[ngoodtrk]=t;
 220     ngoodtrk++;
 221   }
 222
 223
 224 //    cout << " N. of ESD tracks                    : " << ntrk << endl;
 225 //    cout << " N. of preselected tracks            : " << ngoodtrk << endl;
 226 //    cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
 227
 228   AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
 229
 230   for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
 231     AliESDtrack *t=tracks[jtrk];
 232     Double_t time=t->GetTOFsignal();
 233
 234     if((time-mintime*1.E3)<50.E3){ // For pp and per
 235       gtracks[ngoodtrkt0]=t;
 236       ngoodtrkt0++;
 237     }
 238   }
 239
 240
 241   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
 242   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
 243   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
 244
 245   if(ngoodtrkt0<2){
 246 //     cout << "less than 2 tracks, skip event " << endl;
 247     t0def=-999;
 248     deltat0def=0.600;
 249     fT0SigmaT0def[0]=t0def;
 250     fT0SigmaT0def[1]=deltat0def;
 251     fT0SigmaT0def[2]=ngoodtrkt0;
 252     fT0SigmaT0def[3]=ngoodtrkt0;
 253     //goto finish;
 254   }
 255   if(ngoodtrkt0>=2){
 256   // Decide how many tracks in set
 257     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
 258     Int_t nset=1;
 259
 260     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
 261
 262     // Loop over selected sets
 263
 264     if(nset>=1){
 265       for (Int_t i=0; i< nset; i++) {
 266
 267         Float_t t0best=999.;
 268         Float_t eT0best=999.;
 269         Float_t chisquarebest=99999.;
 270         Int_t npionbest=0;
 271
 272         Int_t ntracksinsetmy=0;
 273         AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
 274         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 275           Int_t index = itrk+i*ntracksinset;
 276           if(index < ngoodtrkt0){
 277             AliESDtrack *t=gtracks[index];
 278             tracksT0[itrk]=t;
 279             ntracksinsetmy++;
 280           }
 281         }
 282
 283         // Analyse it
 284
 285         Int_t   assparticle[nmaxtracksinset];
 286         Float_t exptof[nmaxtracksinset][3];
 287         Float_t timeofflight[nmaxtracksinset];
 288         Float_t momentum[nmaxtracksinset];
 289         Float_t timezero[nmaxtracksinset];
 290         Float_t weightedtimezero[nmaxtracksinset];
 291         Float_t beta[nmaxtracksinset];
 292         Float_t texp[nmaxtracksinset];
 293         Float_t dtexp[nmaxtracksinset];
 294         Float_t sqMomError[nmaxtracksinset];
 295         Float_t sqTrackError[nmaxtracksinset];
 296         Float_t massarray[3]={0.13957,0.493677,0.9382723};
 297         Float_t tracktoflen[nmaxtracksinset];
 298         Float_t besttimezero[nmaxtracksinset];
 299         Float_t besttexp[nmaxtracksinset];
 300         Float_t besttimeofflight[nmaxtracksinset];
 301         Float_t bestmomentum[nmaxtracksinset];
 302         Float_t bestchisquare[nmaxtracksinset];
 303         Float_t bestweightedtimezero[nmaxtracksinset];
 304         Float_t bestsqTrackError[nmaxtracksinset];
 305         Int_t imass[nmaxtracksinset];
 306
 307         for (Int_t j=0; j<ntracksinset; j++) {
 308           assparticle[j] = 3;
 309           timeofflight[j] = 0;
 310           momentum[j] = 0;
 311           timezero[j] = 0;
 312           weightedtimezero[j] = 0;
 313           beta[j] = 0;
 314           texp[j] = 0;
 315           dtexp[j] = 0;
 316           sqMomError[j] = 0;
 317           sqTrackError[j] = 0;
 318           tracktoflen[j] = 0;
 319           besttimezero[j] = 0;
 320           besttexp[j] = 0;
 321           besttimeofflight[j] = 0;
 322           bestmomentum[j] = 0;
 323           bestchisquare[j] = 0;
 324           bestweightedtimezero[j] = 0;
 325           bestsqTrackError[j] = 0;
 326           imass[j] = 1;
 327         }
 328
 329         for (Int_t j=0; j<ntracksinsetmy; j++) {
 330           AliESDtrack *t=tracksT0[j];
 331           Double_t momOld=t->GetP();
 332           Double_t mom=momOld-0.0036*momOld;
 333           Double_t time=t->GetTOFsignal();
 334
 335           time*=1.E-3; // tof given in nanoseconds
 336           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
 337           Double_t toflen=t->GetIntegratedLength();
 338           toflen=toflen/100.; // toflen given in m
 339
 340           timeofflight[j]=time;
 341           tracktoflen[j]=toflen;
 342           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
 343           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
 344           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
 345           momentum[j]=mom;
 346           assparticle[j]=3;
 347
 348         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
 349
 350         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 351           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
 352                                        +momentum[itz]*momentum[itz]);
 353           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]));
 354           sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
 355           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
 356           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 357           sumAllweightspi+=1./sqTrackError[itz];
 358           meantzeropi+=weightedtimezero[itz];
 359         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
 360
 361
 362         // Then, Combinatorial Algorithm
 363
 364         if(ntracksinsetmy<2 )break;
 365
 366         for (Int_t j=0; j<ntracksinsetmy; j++) {
 367           imass[j] = 3;
 368         }
 369
 370         Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
 371
 372         // Loop on mass hypotheses
 373         for (Int_t k=0; k < ncombinatorial;k++) {
 374           for (Int_t j=0; j<ntracksinsetmy; j++) {
 375             imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 376             texp[j]=exptof[j][imass[j]];
 377             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 378           }
 379           Float_t sumAllweights=0.;
 380           Float_t meantzero=0.;
 381           Float_t eMeanTzero=0.;
 382
 383           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 384             sqTrackError[itz]=
 385               (timeresolutioninns*
 386                timeresolutioninns
 387                +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
 388
 389             timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 390
 391             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 392             sumAllweights+=1./sqTrackError[itz];
 393             meantzero+=weightedtimezero[itz];
 394
 395           } // end loop for (Int_t itz=0; itz<15;itz++)
 396
 397           meantzero=meantzero/sumAllweights; // it is given in [ns]
 398           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 399
 400           // calculate chisquare
 401
 402           Float_t chisquare=0.;
 403           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 404             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 405
 406           } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 407
 408           if(chisquare<=chisquarebest){
 409             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 410
 411               bestsqTrackError[iqsq]=sqTrackError[iqsq];
 412               besttimezero[iqsq]=timezero[iqsq];
 413               bestmomentum[iqsq]=momentum[iqsq];
 414               besttimeofflight[iqsq]=timeofflight[iqsq];
 415               besttexp[iqsq]=texp[iqsq];
 416               bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 417               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 418             }
 419
 420             Int_t npion=0;
 421             for (Int_t j=0; j<ntracksinsetmy; j++) {
 422               assparticle[j]=imass[j];
 423               if(imass[j] == 0) npion++;
 424             }
 425             npionbest=npion;
 426             chisquarebest=chisquare;
 427             t0best=meantzero;
 428             eT0best=eMeanTzero;
 429           } // close if(dummychisquare<=chisquare)
 430
 431         }
 432
 433         Double_t chi2cut[nmaxtracksinset];
 434         chi2cut[0] = 0;
 435         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
 436         for (Int_t j=2; j<ntracksinset; j++) {
 437           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
 438         }
 439
 440         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
 441
 442 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
 443
 444         Bool_t kRedoT0 = kFALSE;
 445         ntracksinsetmyCut = ntracksinsetmy;
 446         Bool_t usetrack[nmaxtracksinset];
 447         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 448           usetrack[icsq] = kTRUE;
 449           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
 450             kRedoT0 = kTRUE;
 451             ntracksinsetmyCut--;
 452             usetrack[icsq] = kFALSE;
 453           }
 454         } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 455
 456         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
 457
 458         // Loop on mass hypotheses Redo
 459         if(kRedoT0 && ntracksinsetmyCut > 1){
 460           //      printf("Redo T0\n");
 461           for (Int_t k=0; k < ncombinatorial;k++) {
 462             for (Int_t j=0; j<ntracksinsetmy; j++) {
 463               imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 464               texp[j]=exptof[j][imass[j]];
 465               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 466             }
 467
 468             Float_t sumAllweights=0.;
 469             Float_t meantzero=0.;
 470             Float_t eMeanTzero=0.;
 471
 472             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 473               if(! usetrack[itz]) continue;
 474               sqTrackError[itz]=
 475                 (timeresolutioninns*
 476                  timeresolutioninns
 477                  +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
 478
 479               timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 480
 481               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 482               sumAllweights+=1./sqTrackError[itz];
 483               meantzero+=weightedtimezero[itz];
 484
 485             } // end loop for (Int_t itz=0; itz<15;itz++)
 486
 487             meantzero=meantzero/sumAllweights; // it is given in [ns]
 488             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 489
 490             // calculate chisquare
 491
 492             Float_t chisquare=0.;
 493             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 494               if(! usetrack[icsq]) continue;
 495               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 496
 497             } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 498
 499             Int_t npion=0;
 500             for (Int_t j=0; j<ntracksinsetmy; j++) {
 501               assparticle[j]=imass[j];
 502               if(imass[j] == 0) npion++;
 503             }
 504
 505             if(chisquare<=chisquarebest){
 506               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 507                 if(! usetrack[iqsq]) continue;
 508                 bestsqTrackError[iqsq]=sqTrackError[iqsq];
 509                 besttimezero[iqsq]=timezero[iqsq];
 510                 bestmomentum[iqsq]=momentum[iqsq];
 511                 besttimeofflight[iqsq]=timeofflight[iqsq];
 512                 besttexp[iqsq]=texp[iqsq];
 513                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 514                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 515               }
 516
 517               npionbest=npion;
 518               chisquarebest=chisquare;
 519               t0best=meantzero;
 520               eT0best=eMeanTzero;
 521             } // close if(dummychisquare<=chisquare)
 522
 523           }
 524         }
 525
 526         // filling histos
 527         Float_t confLevel=999;
 528
 529         // Sets with decent chisquares
 530
 531         if(chisquarebest<999.){
 532           Double_t dblechisquare=(Double_t)chisquarebest;
 533           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
 534 //        cout << " Set Number " << nsets << endl;
 535 //        cout << "Best Assignment, selection " << assparticle[0] <<
 536 //          assparticle[1] << assparticle[2] <<
 537 //          assparticle[3] << assparticle[4] <<
 538 //          assparticle[5] << endl;
 539 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
 540 //        cout << " C.L. of the set "<< confLevel <<endl;
 541 //        cout << " T0 for this set (in ns)  " << t0best << endl;
 542
 543           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
 544
 545             if(! usetrack[icsq]) continue;
 546
 547 //          cout << "Track # " << icsq  << " T0 offsets = "
 548 //               << besttimezero[icsq]-t0best <<
 549 //            " track error = "  << bestsqTrackError[icsq]
 550 //               << " Chisquare = " << bestchisquare[icsq]
 551 //               << " Momentum  = " << bestmomentum[icsq]
 552 //               << " TOF   = "     << besttimeofflight[icsq]
 553 //               << " TOF tracking  = " << besttexp[icsq]
 554 //               << " is used = " << usetrack[icsq] << endl;
 555           }
 556
 557           // Pick up only those with C.L. >1%
 558           //      if(confLevel>0.01 && ngoodsetsSel<200){
 559           if(confLevel>0.01 && ngoodsetsSel<200){
 560             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
 561             confLevelbestSel[ngoodsetsSel]=confLevel;
 562             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
 563             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
 564             t0bestallSel += t0best/eT0best/eT0best;
 565             sumWt0bestallSel += 1./eT0best/eT0best;
 566             ngoodsetsSel++;
 567             ngoodtrktrulyused+=ntracksinsetmyCut;
 568           }
 569           else{
 570             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
 571           }
 572         }
 573         delete[] tracksT0;
 574         nsets++;
 575
 576       } // end for the current set
 577
 578       nUsedTracks =  ngoodtrkt0;
 579       if(strstr(option,"all")){
 580         if(sumAllweightspi>0.){
 581           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
 582           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
 583         }
 584
 585         if(sumWt0bestallSel>0){
 586           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
 587           eT0bestallSel = sqrt(1./sumWt0bestallSel);
 588
 589         }// end of if(sumWt0bestallSel>0){
 590
 591 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
 592       }
 593
 594       t0def=t0bestallSel;
 595       deltat0def=eT0bestallSel;
 596       if ((TMath::Abs(t0bestallSel) < 0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
 597         t0def=-999; deltat0def=0.600;
 598       }
 599
 600       fT0SigmaT0def[0]=t0def;
 601       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
 602       fT0SigmaT0def[2]=ngoodtrkt0;
 603       fT0SigmaT0def[3]=ngoodtrktrulyused;
 604     }
 605   }
 606
 607   //   if(strstr(option,"tim") || strstr(option,"all")){
 608   //     cout << "AliTOFT0v1:" << endl ;
 609   //}
 610
 611   return fT0SigmaT0def;
 612   }
 613 //__________________________________________________________________
 614 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
 615 {
 616   // Take the error extimate for the TOF time in the track reconstruction
 617
 618   static const Double_t kMasses[]={
 619     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
 620   };
 621
 622   Double_t mass=kMasses[index+2];
 623   Double_t dpp=0.01;      //mean relative pt resolution;
 624   if(mom > 1) dpp = 0.01*mom;
 625   Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
 626
 627   sigma =TMath::Sqrt(sigma*sigma);
 628
 629   return sigma;
 630 }
 631
 632 //__________________________________________________________________
 633 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
 634 {
 635
 636   /* TPC refit */
 637   if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
 638   /* do not accept kink daughters */
 639   if (track->GetKinkIndex(0)>0) return kFALSE;
 640   /* N clusters TPC */
 641   if (track->GetTPCclusters(0) < 50) return kFALSE;
 642   /* chi2 TPC */
 643   if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
 644   /* sigma to vertex */
 645   if (GetSigmaToVertex(track) > 4.) return kFALSE;
 646
 647   /* accept track */
 648   return kTRUE;
 649
 650 }
 651
 652 //____________________________________________________________________
 653 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
 654 {
 655   // Calculates the number of sigma to the vertex.
 656
 657   Float_t b[2];
 658   Float_t bRes[2];
 659   Float_t bCov[3];
 660   esdTrack->GetImpactParameters(b,bCov);
 661
 662   if (bCov[0]<=0 || bCov[2]<=0) {
 663     bCov[0]=0; bCov[2]=0;
 664   }
 665   bRes[0] = TMath::Sqrt(bCov[0]);
 666   bRes[1] = TMath::Sqrt(bCov[2]);
 667
 668   // -----------------------------------
 669   // How to get to a n-sigma cut?
 670   //
 671   // The accumulated statistics from 0 to d is
 672   //
 673   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 674   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 675   //
 676   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
 677   // Can this be expressed in a different way?
 678
 679   if (bRes[0] == 0 || bRes[1] ==0)
 680     return -1;
 681
 682   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 683
 684   // work around precision problem
 685   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 686   // 1e-15 corresponds to nsigma ~ 7.7
 687   if (TMath::Exp(-d * d / 2) < 1e-15)
 688     return 1000;
 689
 690   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 691   return nSigma;
 692 }