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