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 {
  76   //
  77   // default constructor
  78   //
  79   if(AliPID::ParticleMass(0) == 0) new AliPID();
  80
  81   if(!fPIDesd){
  82     fPIDesd = new AliESDpid();
  83   }
  84
  85   Init(NULL);
  86
  87 }
  88
  89 //____________________________________________________________________________
  90 AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
  91   TObject(),
  92   fLowerMomBound(0.5),
  93   fUpperMomBound(3.0),
  94   fTimeCorr(0.),
  95   fEvent(event),
  96   fPIDesd(extPID)
  97 {
  98   //
  99   // real constructor
 100   //
 101   if(AliPID::ParticleMass(0) == 0) new AliPID();
 102
 103   if(!fPIDesd){
 104     fPIDesd = new AliESDpid();
 105   }
 106
 107   Init(event);
 108
 109 }
 110 //____________________________________________________________________________
 111 AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero)
 112 {
 113  //
 114   // assign. operator
 115   //
 116
 117   if (this == &tzero)
 118     return *this;
 119
 120   fLowerMomBound=tzero.fLowerMomBound;
 121   fUpperMomBound=tzero.fUpperMomBound;
 122   fTimeCorr=tzero.fTimeCorr;
 123   fEvent=tzero.fEvent;
 124   fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
 125   fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
 126   fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
 127   fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3];
 128
 129   return *this;
 130 }
 131
 132 //____________________________________________________________________________
 133 AliTOFT0v1::~AliTOFT0v1()
 134 {
 135   // dtor
 136   fEvent=NULL;
 137
 138 }
 139 //____________________________________________________________________________
 140
 141 void
 142 AliTOFT0v1::Init(AliESDEvent *event)
 143 {
 144
 145   /*
 146    * init
 147    */
 148
 149   fEvent = event;
 150   fT0SigmaT0def[0]=0.;
 151   fT0SigmaT0def[1]=0.6;
 152   fT0SigmaT0def[2]=0.;
 153   fT0SigmaT0def[3]=0.;
 154
 155 }
 156
 157 //____________________________________________________________________________
 158 Double_t * AliTOFT0v1::DefineT0(Option_t *option)
 159 {
 160   // Caluclate the Event Time using the ESD TOF time
 161
 162   TBenchmark *bench=new TBenchmark();
 163   bench->Start("t0computation");
 164
 165   fT0SigmaT0def[0]=0.;
 166   fT0SigmaT0def[1]=0.600;
 167   fT0SigmaT0def[2]=0.;
 168   fT0SigmaT0def[3]=0.;
 169
 170   const Int_t nmaxtracksinsetMax=10;
 171   Int_t nmaxtracksinset=10;
 172 //   if(strstr(option,"all")){
 173 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
 174 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
 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   if(ngoodtrk>22) nmaxtracksinset = 6;
 224
 225 //    cout << " N. of ESD tracks                    : " << ntrk << endl;
 226 //    cout << " N. of preselected tracks            : " << ngoodtrk << endl;
 227 //    cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
 228
 229   AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
 230
 231   for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
 232     AliESDtrack *t=tracks[jtrk];
 233     Double_t time=t->GetTOFsignal();
 234
 235     if((time-mintime*1.E3)<50.E3){ // For pp and per
 236       gtracks[ngoodtrkt0]=t;
 237       ngoodtrkt0++;
 238     }
 239   }
 240
 241
 242   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
 243   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
 244   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
 245
 246   if(ngoodtrkt0<2){
 247 //     cout << "less than 2 tracks, skip event " << endl;
 248     t0def=-999;
 249     deltat0def=0.600;
 250     fT0SigmaT0def[0]=t0def;
 251     fT0SigmaT0def[1]=deltat0def;
 252     fT0SigmaT0def[2]=ngoodtrkt0;
 253     fT0SigmaT0def[3]=ngoodtrkt0;
 254     //goto finish;
 255   }
 256   if(ngoodtrkt0>=2){
 257   // Decide how many tracks in set
 258     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
 259     Int_t nset=1;
 260
 261     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
 262
 263     // Loop over selected sets
 264
 265     if(nset>=1){
 266       for (Int_t i=0; i< nset; i++) {
 267
 268         Float_t t0best=999.;
 269         Float_t eT0best=999.;
 270         Float_t chisquarebest=99999.;
 271         Int_t npionbest=0;
 272
 273         Int_t ntracksinsetmy=0;
 274         AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
 275         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 276           Int_t index = itrk+i*ntracksinset;
 277           if(index < ngoodtrkt0){
 278             AliESDtrack *t=gtracks[index];
 279             tracksT0[itrk]=t;
 280             ntracksinsetmy++;
 281           }
 282         }
 283
 284         // Analyse it
 285
 286         Int_t   assparticle[nmaxtracksinsetMax];
 287         Float_t exptof[nmaxtracksinsetMax][3];
 288         Float_t timeofflight[nmaxtracksinsetMax];
 289         Float_t momentum[nmaxtracksinsetMax];
 290         Float_t timezero[nmaxtracksinsetMax];
 291         Float_t weightedtimezero[nmaxtracksinsetMax];
 292         Float_t beta[nmaxtracksinsetMax];
 293         Float_t texp[nmaxtracksinsetMax];
 294         Float_t dtexp[nmaxtracksinsetMax];
 295         Float_t sqMomError[nmaxtracksinsetMax];
 296         Float_t sqTrackError[nmaxtracksinsetMax];
 297         Float_t massarray[3]={0.13957,0.493677,0.9382723};
 298         Float_t tracktoflen[nmaxtracksinsetMax];
 299         Float_t besttimezero[nmaxtracksinsetMax];
 300         Float_t besttexp[nmaxtracksinsetMax];
 301         Float_t besttimeofflight[nmaxtracksinsetMax];
 302         Float_t bestmomentum[nmaxtracksinsetMax];
 303         Float_t bestchisquare[nmaxtracksinsetMax];
 304         Float_t bestweightedtimezero[nmaxtracksinsetMax];
 305         Float_t bestsqTrackError[nmaxtracksinsetMax];
 306         Int_t imass[nmaxtracksinsetMax];
 307
 308         for (Int_t j=0; j<ntracksinset; j++) {
 309           assparticle[j] = 3;
 310           timeofflight[j] = 0;
 311           momentum[j] = 0;
 312           timezero[j] = 0;
 313           weightedtimezero[j] = 0;
 314           beta[j] = 0;
 315           texp[j] = 0;
 316           dtexp[j] = 0;
 317           sqMomError[j] = 0;
 318           sqTrackError[j] = 0;
 319           tracktoflen[j] = 0;
 320           besttimezero[j] = 0;
 321           besttexp[j] = 0;
 322           besttimeofflight[j] = 0;
 323           bestmomentum[j] = 0;
 324           bestchisquare[j] = 0;
 325           bestweightedtimezero[j] = 0;
 326           bestsqTrackError[j] = 0;
 327           imass[j] = 1;
 328         }
 329
 330         for (Int_t j=0; j<ntracksinsetmy; j++) {
 331           AliESDtrack *t=tracksT0[j];
 332           Double_t momOld=t->GetP();
 333           Double_t mom=momOld-0.0036*momOld;
 334           Double_t time=t->GetTOFsignal();
 335
 336           time*=1.E-3; // tof given in nanoseconds
 337           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
 338           Double_t toflen=t->GetIntegratedLength();
 339           toflen=toflen/100.; // toflen given in m
 340
 341           timeofflight[j]=time;
 342           tracktoflen[j]=toflen;
 343           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
 344           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
 345           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
 346           momentum[j]=mom;
 347           assparticle[j]=3;
 348
 349         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
 350
 351         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 352           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
 353                                        +momentum[itz]*momentum[itz]);
 354           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]));
 355           sqTrackError[itz]=sqMomError[itz]; //in ns
 356           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
 357           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 358           sumAllweightspi+=1./sqTrackError[itz];
 359           meantzeropi+=weightedtimezero[itz];
 360         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
 361
 362
 363         // Then, Combinatorial Algorithm
 364
 365         if(ntracksinsetmy<2 )break;
 366
 367         for (Int_t j=0; j<ntracksinsetmy; j++) {
 368           imass[j] = 3;
 369         }
 370
 371         Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
 372
 373         // Loop on mass hypotheses
 374         for (Int_t k=0; k < ncombinatorial;k++) {
 375           for (Int_t j=0; j<ntracksinsetmy; j++) {
 376             imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 377             texp[j]=exptof[j][imass[j]];
 378             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 379           }
 380           Float_t sumAllweights=0.;
 381           Float_t meantzero=0.;
 382           Float_t eMeanTzero=0.;
 383
 384           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 385             sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 386
 387             timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 388
 389             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 390             sumAllweights+=1./sqTrackError[itz];
 391             meantzero+=weightedtimezero[itz];
 392
 393           } // end loop for (Int_t itz=0; itz<15;itz++)
 394
 395           meantzero=meantzero/sumAllweights; // it is given in [ns]
 396           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 397
 398           // calculate chisquare
 399
 400           Float_t chisquare=0.;
 401           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 402             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 403
 404           } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 405
 406           if(chisquare<=chisquarebest){
 407             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 408
 409               bestsqTrackError[iqsq]=sqTrackError[iqsq];
 410               besttimezero[iqsq]=timezero[iqsq];
 411               bestmomentum[iqsq]=momentum[iqsq];
 412               besttimeofflight[iqsq]=timeofflight[iqsq];
 413               besttexp[iqsq]=texp[iqsq];
 414               bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 415               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 416             }
 417
 418             Int_t npion=0;
 419             for (Int_t j=0; j<ntracksinsetmy; j++) {
 420               assparticle[j]=imass[j];
 421               if(imass[j] == 0) npion++;
 422             }
 423             npionbest=npion;
 424             chisquarebest=chisquare;
 425             t0best=meantzero;
 426             eT0best=eMeanTzero;
 427           } // close if(dummychisquare<=chisquare)
 428
 429         }
 430
 431         Double_t chi2cut[nmaxtracksinsetMax];
 432         chi2cut[0] = 0;
 433         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
 434         for (Int_t j=2; j<ntracksinset; j++) {
 435           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
 436         }
 437
 438         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
 439
 440 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
 441
 442         Bool_t kRedoT0 = kFALSE;
 443         ntracksinsetmyCut = ntracksinsetmy;
 444         Bool_t usetrack[nmaxtracksinsetMax];
 445         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 446           usetrack[icsq] = kTRUE;
 447           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
 448             kRedoT0 = kTRUE;
 449             ntracksinsetmyCut--;
 450             usetrack[icsq] = kFALSE;
 451           }
 452         } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 453
 454         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
 455
 456         // Loop on mass hypotheses Redo
 457         if(kRedoT0 && ntracksinsetmyCut > 1){
 458           //      printf("Redo T0\n");
 459           for (Int_t k=0; k < ncombinatorial;k++) {
 460             for (Int_t j=0; j<ntracksinsetmy; j++) {
 461               imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 462               texp[j]=exptof[j][imass[j]];
 463               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 464             }
 465
 466             Float_t sumAllweights=0.;
 467             Float_t meantzero=0.;
 468             Float_t eMeanTzero=0.;
 469
 470             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 471               if(! usetrack[itz]) continue;
 472               sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 473
 474               timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 475
 476               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 477               sumAllweights+=1./sqTrackError[itz];
 478               meantzero+=weightedtimezero[itz];
 479
 480             } // end loop for (Int_t itz=0; itz<15;itz++)
 481
 482             meantzero=meantzero/sumAllweights; // it is given in [ns]
 483             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 484
 485             // calculate chisquare
 486
 487             Float_t chisquare=0.;
 488             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 489               if(! usetrack[icsq]) continue;
 490               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 491
 492             } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 493
 494             Int_t npion=0;
 495             for (Int_t j=0; j<ntracksinsetmy; j++) {
 496               assparticle[j]=imass[j];
 497               if(imass[j] == 0) npion++;
 498             }
 499
 500             if(chisquare<=chisquarebest){
 501               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 502                 if(! usetrack[iqsq]) continue;
 503                 bestsqTrackError[iqsq]=sqTrackError[iqsq];
 504                 besttimezero[iqsq]=timezero[iqsq];
 505                 bestmomentum[iqsq]=momentum[iqsq];
 506                 besttimeofflight[iqsq]=timeofflight[iqsq];
 507                 besttexp[iqsq]=texp[iqsq];
 508                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 509                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 510               }
 511
 512               npionbest=npion;
 513               chisquarebest=chisquare;
 514               t0best=meantzero;
 515               eT0best=eMeanTzero;
 516             } // close if(dummychisquare<=chisquare)
 517
 518           }
 519         }
 520
 521         // filling histos
 522         Float_t confLevel=999;
 523
 524         // Sets with decent chisquares
 525
 526         if(chisquarebest<999.){
 527           Double_t dblechisquare=(Double_t)chisquarebest;
 528           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
 529 //        cout << " Set Number " << nsets << endl;
 530 //        cout << "Best Assignment, selection " << assparticle[0] <<
 531 //          assparticle[1] << assparticle[2] <<
 532 //          assparticle[3] << assparticle[4] <<
 533 //          assparticle[5] << endl;
 534 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
 535 //        cout << " C.L. of the set "<< confLevel <<endl;
 536 //        cout << " T0 for this set (in ns)  " << t0best << endl;
 537
 538           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
 539
 540             if(! usetrack[icsq]) continue;
 541
 542 //          cout << "Track # " << icsq  << " T0 offsets = "
 543 //               << besttimezero[icsq]-t0best <<
 544 //            " track error = "  << bestsqTrackError[icsq]
 545 //               << " Chisquare = " << bestchisquare[icsq]
 546 //               << " Momentum  = " << bestmomentum[icsq]
 547 //               << " TOF   = "     << besttimeofflight[icsq]
 548 //               << " TOF tracking  = " << besttexp[icsq]
 549 //               << " is used = " << usetrack[icsq] << endl;
 550           }
 551
 552           // Pick up only those with C.L. >1%
 553           //      if(confLevel>0.01 && ngoodsetsSel<200){
 554           if(confLevel>0.01 && ngoodsetsSel<200){
 555             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
 556             confLevelbestSel[ngoodsetsSel]=confLevel;
 557             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
 558             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
 559             t0bestallSel += t0best/eT0best/eT0best;
 560             sumWt0bestallSel += 1./eT0best/eT0best;
 561             ngoodsetsSel++;
 562             ngoodtrktrulyused+=ntracksinsetmyCut;
 563
 564             //      printf("t0 of a set = %f +/- %f\n",t0best,eT0best);
 565           }
 566           else{
 567             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
 568           }
 569         }
 570         delete[] tracksT0;
 571         nsets++;
 572
 573       } // end for the current set
 574
 575       nUsedTracks =  ngoodtrkt0;
 576       if(strstr(option,"all")){
 577         if(sumAllweightspi>0.){
 578           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
 579           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
 580         }
 581
 582         if(sumWt0bestallSel>0){
 583           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
 584           eT0bestallSel = sqrt(1./sumWt0bestallSel);
 585
 586         }// end of if(sumWt0bestallSel>0){
 587
 588 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
 589       }
 590
 591       t0def=t0bestallSel;
 592       deltat0def=eT0bestallSel;
 593
 594       fT0SigmaT0def[0]=t0def;
 595       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
 596       fT0SigmaT0def[2]=ngoodtrkt0;
 597       fT0SigmaT0def[3]=ngoodtrktrulyused;
 598     }
 599   }
 600
 601   //   if(strstr(option,"tim") || strstr(option,"all")){
 602   //     cout << "AliTOFT0v1:" << endl ;
 603   //}
 604
 605   if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
 606
 607   bench->Stop("t0computation");
 608
 609   fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
 610   fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
 611
 612 //   bench->Print("t0computation");
 613 //   printf("%f %f\n",fT0SigmaT0def[4],fT0SigmaT0def[5]);
 614
 615   bench->Reset();
 616
 617   return fT0SigmaT0def;
 618   }
 619 //__________________________________________________________________
 620 Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
 621 {
 622   TBenchmark *bench=new TBenchmark();
 623   bench->Start("t0computation");
 624
 625   // Caluclate the Event Time using the ESD TOF time
 626
 627   fT0SigmaT0def[0]=0.;
 628   fT0SigmaT0def[1]=0.600;
 629   fT0SigmaT0def[2]=0.;
 630   fT0SigmaT0def[3]=0.;
 631
 632   const Int_t nmaxtracksinsetMax=10;
 633   Int_t nmaxtracksinset=10;
 634 //   if(strstr(option,"all")){
 635 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
 636 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
 637 //   }
 638
 639
 640   Int_t nsets=0;
 641   Int_t nUsedTracks=0;
 642   Int_t ngoodsetsSel= 0;
 643   Float_t t0bestSel[300];
 644   Float_t eT0bestSel[300];
 645   Float_t chiSquarebestSel[300];
 646   Float_t confLevelbestSel[300];
 647   Float_t t0bestallSel=0.;
 648   Float_t eT0bestallSel=0.;
 649   Float_t sumWt0bestallSel=0.;
 650   Float_t eMeanTzeroPi=0.;
 651   Float_t meantzeropi=0.;
 652   Float_t sumAllweightspi=0.;
 653   Double_t t0def=-999;
 654   Double_t deltat0def=999;
 655   Int_t ngoodtrktrulyused=0;
 656   Int_t ntracksinsetmyCut = 0;
 657
 658   Int_t ntrk=fEvent->GetNumberOfTracks();
 659
 660   AliESDtrack **tracks=new AliESDtrack*[ntrk];
 661   Int_t ngoodtrk=0;
 662   Int_t ngoodtrkt0 =0;
 663   Float_t mintime =1E6;
 664
 665   // First Track loop, Selection of good tracks
 666
 667   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 668     AliESDtrack *t=fEvent->GetTrack(itrk);
 669     Double_t momOld=t->GetP();
 670     Double_t mom=momOld-0.0036*momOld;
 671     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 672     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 673     Double_t time=t->GetTOFsignal();
 674
 675     time*=1.E-3; // tof given in nanoseconds
 676     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
 677
 678     if (!AcceptTrack(t)) continue;
 679
 680     if(t->GetIntegratedLength() < 350)continue; //skip decays
 681     if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
 682     if(time <= mintime) mintime=time;
 683     tracks[ngoodtrk]=t;
 684     ngoodtrk++;
 685   }
 686
 687   if(ngoodtrk>22) nmaxtracksinset = 6;
 688
 689
 690 //     cout << " N. of ESD tracks                    : " << ntrk << endl;
 691 //     cout << " N. of preselected tracks            : " << ngoodtrk << endl;
 692 //     cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
 693
 694   AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
 695
 696   for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
 697     AliESDtrack *t=tracks[jtrk];
 698     Double_t time=t->GetTOFsignal();
 699
 700     if((time-mintime*1.E3)<50.E3){ // For pp and per
 701       gtracks[ngoodtrkt0]=t;
 702       ngoodtrkt0++;
 703     }
 704   }
 705
 706 //     cout << " N. of preselected tracks t0         : " << ngoodtrkt0 << endl;
 707
 708   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
 709   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
 710   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
 711
 712 //   cout << " N. of sets                        : " << nseteq << endl;
 713
 714   if(ngoodtrkt0<2){
 715 //     cout << "less than 2 tracks, skip event " << endl;
 716     t0def=-999;
 717     deltat0def=0.600;
 718     fT0SigmaT0def[0]=t0def;
 719     fT0SigmaT0def[1]=deltat0def;
 720     fT0SigmaT0def[2]=ngoodtrkt0;
 721     fT0SigmaT0def[3]=ngoodtrkt0;
 722     //goto finish;
 723   }
 724   if(ngoodtrkt0>=2){
 725   // Decide how many tracks in set
 726     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
 727     Int_t nset=1;
 728
 729     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
 730
 731     // Loop over selected sets
 732
 733     if(nset>=1){
 734       for (Int_t i=0; i< nset; i++) {
 735         //      printf("Set %i of %i\n",i+1,nset);
 736         Float_t t0best=999.;
 737         Float_t eT0best=999.;
 738         Float_t chisquarebest=99999.;
 739         Int_t npionbest=0;
 740
 741         Int_t ntracksinsetmy=0;
 742         AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
 743         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 744           Int_t index = itrk+i*ntracksinset;
 745           if(index < ngoodtrkt0){
 746             AliESDtrack *t=gtracks[index];
 747             tracksT0[itrk]=t;
 748             ntracksinsetmy++;
 749           }
 750         }
 751
 752         // Analyse it
 753
 754         Int_t   assparticle[nmaxtracksinsetMax];
 755         Float_t exptof[nmaxtracksinsetMax][3];
 756         Float_t timeofflight[nmaxtracksinsetMax];
 757         Float_t momentum[nmaxtracksinsetMax];
 758         Float_t timezero[nmaxtracksinsetMax];
 759         Float_t weightedtimezero[nmaxtracksinsetMax];
 760         Float_t beta[nmaxtracksinsetMax];
 761         Float_t texp[nmaxtracksinsetMax];
 762         Float_t dtexp[nmaxtracksinsetMax];
 763         Float_t sqMomError[nmaxtracksinsetMax];
 764         Float_t sqTrackError[nmaxtracksinsetMax];
 765         Float_t massarray[3]={0.13957,0.493677,0.9382723};
 766         Float_t tracktoflen[nmaxtracksinsetMax];
 767         Float_t besttimezero[nmaxtracksinsetMax];
 768         Float_t besttexp[nmaxtracksinsetMax];
 769         Float_t besttimeofflight[nmaxtracksinsetMax];
 770         Float_t bestmomentum[nmaxtracksinsetMax];
 771         Float_t bestchisquare[nmaxtracksinsetMax];
 772         Float_t bestweightedtimezero[nmaxtracksinsetMax];
 773         Float_t bestsqTrackError[nmaxtracksinsetMax];
 774         Int_t imass[nmaxtracksinsetMax];
 775
 776         for (Int_t j=0; j<ntracksinset; j++) {
 777           assparticle[j] = 3;
 778           timeofflight[j] = 0;
 779           momentum[j] = 0;
 780           timezero[j] = 0;
 781           weightedtimezero[j] = 0;
 782           beta[j] = 0;
 783           texp[j] = 0;
 784           dtexp[j] = 0;
 785           sqMomError[j] = 0;
 786           sqTrackError[j] = 0;
 787           tracktoflen[j] = 0;
 788           besttimezero[j] = 0;
 789           besttexp[j] = 0;
 790           besttimeofflight[j] = 0;
 791           bestmomentum[j] = 0;
 792           bestchisquare[j] = 0;
 793           bestweightedtimezero[j] = 0;
 794           bestsqTrackError[j] = 0;
 795           imass[j] = 1;
 796         }
 797
 798         for (Int_t j=0; j<ntracksinsetmy; j++) {
 799           AliESDtrack *t=tracksT0[j];
 800           Double_t momOld=t->GetP();
 801           Double_t mom=momOld-0.0036*momOld;
 802           Double_t time=t->GetTOFsignal();
 803
 804           time*=1.E-3; // tof given in nanoseconds
 805           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
 806           Double_t toflen=t->GetIntegratedLength();
 807           toflen=toflen/100.; // toflen given in m
 808
 809           timeofflight[j]=time;
 810           tracktoflen[j]=toflen;
 811           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
 812           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
 813           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
 814           momentum[j]=mom;
 815           assparticle[j]=3;
 816
 817         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
 818
 819         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 820           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
 821                                        +momentum[itz]*momentum[itz]);
 822           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]));
 823           sqTrackError[itz]=sqMomError[itz]; //in ns
 824           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
 825           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 826           sumAllweightspi+=1./sqTrackError[itz];
 827           meantzeropi+=weightedtimezero[itz];
 828         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
 829
 830
 831         // Then, Combinatorial Algorithm
 832
 833         if(ntracksinsetmy<2 )break;
 834
 835         for (Int_t j=0; j<ntracksinsetmy; j++) {
 836           imass[j] = 3;
 837         }
 838
 839         Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
 840
 841         // Loop on mass hypotheses
 842         for (Int_t k=0; k < ncombinatorial;k++) {
 843           for (Int_t j=0; j<ntracksinsetmy; j++) {
 844             imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 845             texp[j]=exptof[j][imass[j]];
 846             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 847           }
 848           Float_t sumAllweights=0.;
 849           Float_t meantzero=0.;
 850           Float_t eMeanTzero=0.;
 851
 852           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 853             sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 854
 855             timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 856
 857             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 858             sumAllweights+=1./sqTrackError[itz];
 859             meantzero+=weightedtimezero[itz];
 860
 861           } // end loop for (Int_t itz=0; itz<15;itz++)
 862
 863           meantzero=meantzero/sumAllweights; // it is given in [ns]
 864           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 865
 866           // calculate chisquare
 867
 868           Float_t chisquare=0.;
 869           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 870             chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 871           } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 872
 873           if(chisquare<=chisquarebest){
 874             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 875
 876               bestsqTrackError[iqsq]=sqTrackError[iqsq];
 877               besttimezero[iqsq]=timezero[iqsq];
 878               bestmomentum[iqsq]=momentum[iqsq];
 879               besttimeofflight[iqsq]=timeofflight[iqsq];
 880               besttexp[iqsq]=texp[iqsq];
 881               bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 882               bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 883             }
 884
 885             Int_t npion=0;
 886             for (Int_t j=0; j<ntracksinsetmy; j++) {
 887               assparticle[j]=imass[j];
 888               if(imass[j] == 0) npion++;
 889             }
 890             npionbest=npion;
 891             chisquarebest=chisquare;
 892             t0best=meantzero;
 893             eT0best=eMeanTzero;
 894           } // close if(dummychisquare<=chisquare)
 895         }
 896
 897         Double_t chi2cut[nmaxtracksinsetMax];
 898         chi2cut[0] = 0;
 899         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
 900         for (Int_t j=2; j<ntracksinset; j++) {
 901           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
 902         }
 903
 904         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
 905
 906         //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
 907
 908         Bool_t kRedoT0 = kFALSE;
 909         ntracksinsetmyCut = ntracksinsetmy;
 910         Bool_t usetrack[nmaxtracksinsetMax];
 911         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 912           usetrack[icsq] = kTRUE;
 913           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
 914             kRedoT0 = kTRUE;
 915             ntracksinsetmyCut--;
 916             usetrack[icsq] = kFALSE;
 917           }
 918         } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 919
 920         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
 921
 922         // Loop on mass hypotheses Redo
 923         if(kRedoT0 && ntracksinsetmyCut > 1){
 924           //      printf("Redo T0\n");
 925           for (Int_t k=0; k < ncombinatorial;k++) {
 926             for (Int_t j=0; j<ntracksinsetmy; j++) {
 927               imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 928               texp[j]=exptof[j][imass[j]];
 929               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 930             }
 931
 932             Float_t sumAllweights=0.;
 933             Float_t meantzero=0.;
 934             Float_t eMeanTzero=0.;
 935
 936             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 937               if(! usetrack[itz]) continue;
 938               sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 939
 940               timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 941
 942               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 943               sumAllweights+=1./sqTrackError[itz];
 944               meantzero+=weightedtimezero[itz];
 945
 946             } // end loop for (Int_t itz=0; itz<15;itz++)
 947
 948             meantzero=meantzero/sumAllweights; // it is given in [ns]
 949             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 950
 951             // calculate chisquare
 952
 953             Float_t chisquare=0.;
 954             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 955               if(! usetrack[icsq]) continue;
 956               chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 957
 958             } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 959
 960             Int_t npion=0;
 961             for (Int_t j=0; j<ntracksinsetmy; j++) {
 962               assparticle[j]=imass[j];
 963               if(imass[j] == 0) npion++;
 964             }
 965
 966             if(chisquare<=chisquarebest){
 967               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 968                 if(! usetrack[iqsq]) continue;
 969                 bestsqTrackError[iqsq]=sqTrackError[iqsq];
 970                 besttimezero[iqsq]=timezero[iqsq];
 971                 bestmomentum[iqsq]=momentum[iqsq];
 972                 besttimeofflight[iqsq]=timeofflight[iqsq];
 973                 besttexp[iqsq]=texp[iqsq];
 974                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 975                 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 976               }
 977
 978               npionbest=npion;
 979               chisquarebest=chisquare;
 980               t0best=meantzero;
 981               eT0best=eMeanTzero;
 982             } // close if(dummychisquare<=chisquare)
 983
 984           }
 985         }
 986
 987         // filling histos
 988         Float_t confLevel=999;
 989
 990         // Sets with decent chisquares
 991         //      printf("Chi2best of the set = %f \n",chisquarebest);
 992
 993         if(chisquarebest<999.){
 994           Double_t dblechisquare=(Double_t)chisquarebest;
 995           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
 996 //        cout << " Set Number " << nsets << endl;
 997 //        cout << "Best Assignment, selection " << assparticle[0] <<
 998 //          assparticle[1] << assparticle[2] <<
 999 //          assparticle[3] << assparticle[4] <<
1000 //          assparticle[5] << endl;
1001 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
1002 //        cout << " C.L. of the set "<< confLevel <<endl;
1003 //        cout << " T0 for this set (in ns)  " << t0best << endl;
1004
1005           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1006
1007             if(! usetrack[icsq]) continue;
1008
1009 //          cout << "Track # " << icsq  << " T0 offsets = "
1010 //               << besttimezero[icsq]-t0best <<
1011 //            " track error = "  << bestsqTrackError[icsq]
1012 //               << " Chisquare = " << bestchisquare[icsq]
1013 //               << " Momentum  = " << bestmomentum[icsq]
1014 //               << " TOF   = "     << besttimeofflight[icsq]
1015 //               << " TOF tracking  = " << besttexp[icsq]
1016 //               << " is used = " << usetrack[icsq] << endl;
1017           }
1018
1019           //      printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
1020
1021           // Pick up only those with C.L. >1%
1022           //      if(confLevel>0.01 && ngoodsetsSel<200){
1023           if(confLevel>0.01 && ngoodsetsSel<200){
1024             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
1025             confLevelbestSel[ngoodsetsSel]=confLevel;
1026             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
1027             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
1028             t0bestallSel += t0best/eT0best/eT0best;
1029             sumWt0bestallSel += 1./eT0best/eT0best;
1030             ngoodsetsSel++;
1031             ngoodtrktrulyused+=ntracksinsetmyCut;
1032             //      printf("T0 set = %f +/- %f\n",t0best,eT0best);
1033           }
1034           else{
1035             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
1036           }
1037         }
1038         delete[] tracksT0;
1039         nsets++;
1040
1041       } // end for the current set
1042
1043       nUsedTracks =  ngoodtrkt0;
1044       if(strstr(option,"all")){
1045         if(sumAllweightspi>0.){
1046           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
1047           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
1048         }
1049
1050         //      printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
1051
1052         if(sumWt0bestallSel>0){
1053           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
1054           eT0bestallSel = sqrt(1./sumWt0bestallSel);
1055           //      printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
1056         }// end of if(sumWt0bestallSel>0){
1057
1058 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
1059       }
1060
1061       t0def=t0bestallSel;
1062       deltat0def=eT0bestallSel;
1063
1064       fT0SigmaT0def[0]=t0def;
1065       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
1066       fT0SigmaT0def[2]=ngoodtrkt0;
1067       fT0SigmaT0def[3]=ngoodtrktrulyused;
1068     }
1069   }
1070
1071   //   if(strstr(option,"tim") || strstr(option,"all")){
1072   //     cout << "AliTOFT0v1:" << endl ;
1073   //}
1074
1075   if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
1076
1077   bench->Stop("t0computation");
1078
1079   fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
1080   fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
1081
1082 //   bench->Print("t0computation");
1083 //   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]));
1084
1085   bench->Reset();
1086
1087   return fT0SigmaT0def;
1088   }
1089 //__________________________________________________________________
1090 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
1091 {
1092   // Take the error extimate for the TOF time in the track reconstruction
1093
1094   static const Double_t kMasses[]={
1095     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
1096   };
1097
1098   Double_t mass=kMasses[index+2];
1099
1100   Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
1101
1102   return sigma;
1103 }
1104
1105 //__________________________________________________________________
1106 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
1107 {
1108
1109   /* TPC refit */
1110   if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1111   /* do not accept kink daughters */
1112   if (track->GetKinkIndex(0)>0) return kFALSE;
1113   /* N clusters TPC */
1114   if (track->GetTPCclusters(0) < 50) return kFALSE;
1115   /* chi2 TPC */
1116   if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
1117   /* sigma to vertex */
1118   if (GetSigmaToVertex(track) > 4.) return kFALSE;
1119
1120   /* accept track */
1121   return kTRUE;
1122
1123 }
1124
1125 //____________________________________________________________________
1126 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
1127 {
1128   // Calculates the number of sigma to the vertex.
1129
1130   Float_t b[2];
1131   Float_t bRes[2];
1132   Float_t bCov[3];
1133   esdTrack->GetImpactParameters(b,bCov);
1134
1135   if (bCov[0]<=0 || bCov[2]<=0) {
1136     bCov[0]=0; bCov[2]=0;
1137   }
1138   bRes[0] = TMath::Sqrt(bCov[0]);
1139   bRes[1] = TMath::Sqrt(bCov[2]);
1140
1141   // -----------------------------------
1142   // How to get to a n-sigma cut?
1143   //
1144   // The accumulated statistics from 0 to d is
1145   //
1146   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1147   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1148   //
1149   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1150   // Can this be expressed in a different way?
1151
1152   if (bRes[0] == 0 || bRes[1] ==0)
1153     return -1;
1154
1155   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1156
1157   // work around precision problem
1158   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1159   // 1e-15 corresponds to nsigma ~ 7.7
1160   if (TMath::Exp(-d * d / 2) < 1e-15)
1161     return 1000;
1162
1163   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1164   return nSigma;
1165 }