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