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   fTracks->Delete();
 158   delete fTracks;
 159   fTracks=0x0;
 160
 161   fGTracks->Delete();
 162   delete fGTracks;
 163   fGTracks=0x0;
 164
 165   fTracksT0->Delete();
 166   delete fTracksT0;
 167   fTracksT0=0x0;
 168
 169 }
 170 //____________________________________________________________________________
 171
 172 void
 173 AliTOFT0v1::Init(AliESDEvent *event)
 174 {
 175
 176   /*
 177    * init
 178    */
 179
 180   fEvent = event;
 181   fT0SigmaT0def[0]=0.;
 182   fT0SigmaT0def[1]=0.6;
 183   fT0SigmaT0def[2]=0.;
 184   fT0SigmaT0def[3]=0.;
 185
 186 }
 187
 188 //____________________________________________________________________________
 189 Double_t * AliTOFT0v1::DefineT0(Option_t *option)
 190 {
 191   // Caluclate the Event Time using the ESD TOF time
 192
 193   TBenchmark *bench=new TBenchmark();
 194   bench->Start("t0computation");
 195
 196   fT0SigmaT0def[0]=0.;
 197   fT0SigmaT0def[1]=0.600;
 198   fT0SigmaT0def[2]=0.;
 199   fT0SigmaT0def[3]=0.;
 200
 201   const Int_t nmaxtracksinsetMax=10;
 202   Int_t nmaxtracksinset=10;
 203 //   if(strstr(option,"all")){
 204 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
 205 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
 206 //   }
 207
 208   Int_t nsets=0;
 209   Int_t nUsedTracks=0;
 210   Int_t ngoodsetsSel= 0;
 211   Float_t t0bestSel[300];
 212   Float_t eT0bestSel[300];
 213   Float_t chiSquarebestSel[300];
 214   Float_t confLevelbestSel[300];
 215   Float_t t0bestallSel=0.;
 216   Float_t eT0bestallSel=0.;
 217   Float_t sumWt0bestallSel=0.;
 218   Float_t eMeanTzeroPi=0.;
 219   Float_t meantzeropi=0.;
 220   Float_t sumAllweightspi=0.;
 221   Double_t t0def=-999;
 222   Double_t deltat0def=999;
 223   Int_t ngoodtrktrulyused=0;
 224   Int_t ntracksinsetmyCut = 0;
 225
 226   Int_t ntrk=fEvent->GetNumberOfTracks();
 227
 228   //AliESDtrack **fTracks=new AliESDtrack*[ntrk];
 229   Int_t ngoodtrk=0;
 230   Int_t ngoodtrkt0 =0;
 231   Float_t mintime =1E6;
 232
 233   // First Track loop, Selection of good tracks
 234
 235   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 236     AliESDtrack *t=fEvent->GetTrack(itrk);
 237     Double_t momOld=t->GetP();
 238     Double_t mom=momOld-0.0036*momOld;
 239     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 240     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 241     Double_t time=t->GetTOFsignal();
 242
 243     time*=1.E-3; // tof given in nanoseconds
 244     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
 245
 246     if (!AcceptTrack(t)) continue;
 247
 248     if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
 249
 250     if(time <= mintime) mintime=time;
 251     //fTracks[ngoodtrk]=t;
 252     fTracks->AddLast(t);
 253     ngoodtrk++;
 254   }
 255
 256   if(ngoodtrk>22) nmaxtracksinset = 6;
 257
 258 //    cout << " N. of ESD tracks                    : " << ntrk << endl;
 259 //    cout << " N. of preselected tracks            : " << ngoodtrk << endl;
 260 //    cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
 261
 262   //AliESDtrack **fGTracks=new AliESDtrack*[ngoodtrk];
 263
 264   //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
 265   for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
 266     //AliESDtrack *t=fTracks[jtrk];
 267     AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
 268     Double_t time=t->GetTOFsignal();
 269
 270     if((time-mintime*1.E3)<50.E3){ // For pp and per
 271       //fGTracks[ngoodtrkt0]=t;
 272       fGTracks->AddLast(t);
 273       ngoodtrkt0++;
 274     }
 275   }
 276
 277   //delete[] fTracks;
 278   fTracks->Clear();
 279
 280   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
 281   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
 282   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
 283
 284   if(ngoodtrkt0<2){
 285 //     cout << "less than 2 tracks, skip event " << endl;
 286     t0def=-999;
 287     deltat0def=0.600;
 288     fT0SigmaT0def[0]=t0def;
 289     fT0SigmaT0def[1]=deltat0def;
 290     fT0SigmaT0def[2]=ngoodtrkt0;
 291     fT0SigmaT0def[3]=ngoodtrkt0;
 292     //goto finish;
 293   }
 294   if(ngoodtrkt0>=2){
 295   // Decide how many tracks in set
 296     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
 297     Int_t nset=1;
 298
 299     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
 300
 301     // Loop over selected sets
 302
 303     if(nset>=1){
 304       for (Int_t i=0; i< nset; i++) {
 305
 306         Float_t t0best=999.;
 307         Float_t eT0best=999.;
 308         Float_t chisquarebest=99999.;
 309         Int_t npionbest=0;
 310
 311         Int_t ntracksinsetmy=0;
 312         //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset];
 313         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 314           Int_t index = itrk+i*ntracksinset;
 315           //if(index < ngoodtrkt0){
 316           if(index < fGTracks->GetEntries()){
 317             //AliESDtrack *t=fGTracks[index];
 318             AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
 319             //fTracksT0[itrk]=t;
 320             fTracksT0->AddLast(t);
 321             ntracksinsetmy++;
 322           }
 323         }
 324
 325         // Analyse it
 326
 327         Int_t   assparticle[nmaxtracksinsetMax];
 328         Float_t exptof[nmaxtracksinsetMax][3];
 329         Float_t timeofflight[nmaxtracksinsetMax];
 330         Float_t momentum[nmaxtracksinsetMax];
 331         Float_t timezero[nmaxtracksinsetMax];
 332         Float_t weightedtimezero[nmaxtracksinsetMax];
 333         Float_t beta[nmaxtracksinsetMax];
 334         Float_t texp[nmaxtracksinsetMax];
 335         Float_t dtexp[nmaxtracksinsetMax];
 336         Float_t sqMomError[nmaxtracksinsetMax];
 337         Float_t sqTrackError[nmaxtracksinsetMax];
 338         Float_t massarray[3]={0.13957,0.493677,0.9382723};
 339         Float_t tracktoflen[nmaxtracksinsetMax];
 340         Float_t besttimezero[nmaxtracksinsetMax];
 341         Float_t besttexp[nmaxtracksinsetMax];
 342         Float_t besttimeofflight[nmaxtracksinsetMax];
 343         Float_t bestmomentum[nmaxtracksinsetMax];
 344         Float_t bestchisquare[nmaxtracksinsetMax];
 345         Float_t bestweightedtimezero[nmaxtracksinsetMax];
 346         Float_t bestsqTrackError[nmaxtracksinsetMax];
 347         Int_t imass[nmaxtracksinsetMax];
 348
 349         for (Int_t j=0; j<ntracksinset; j++) {
 350           assparticle[j] = 3;
 351           timeofflight[j] = 0;
 352           momentum[j] = 0;
 353           timezero[j] = 0;
 354           weightedtimezero[j] = 0;
 355           beta[j] = 0;
 356           texp[j] = 0;
 357           dtexp[j] = 0;
 358           sqMomError[j] = 0;
 359           sqTrackError[j] = 0;
 360           tracktoflen[j] = 0;
 361           besttimezero[j] = 0;
 362           besttexp[j] = 0;
 363           besttimeofflight[j] = 0;
 364           bestmomentum[j] = 0;
 365           bestchisquare[j] = 0;
 366           bestweightedtimezero[j] = 0;
 367           bestsqTrackError[j] = 0;
 368           imass[j] = 1;
 369         }
 370
 371         //for (Int_t j=0; j<ntracksinsetmy; j++) {
 372         for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
 373           //AliESDtrack *t=fTracksT0[j];
 374           AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
 375           Double_t momOld=t->GetP();
 376           Double_t mom=momOld-0.0036*momOld;
 377           Double_t time=t->GetTOFsignal();
 378
 379           time*=1.E-3; // tof given in nanoseconds
 380           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
 381           Double_t toflen=t->GetIntegratedLength();
 382           toflen=toflen/100.; // toflen given in m
 383
 384           timeofflight[j]=time;
 385           tracktoflen[j]=toflen;
 386           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
 387           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
 388           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
 389           momentum[j]=mom;
 390           assparticle[j]=3;
 391
 392         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
 393
 394         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 395           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
 396                                        +momentum[itz]*momentum[itz]);
 397           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]));
 398           sqTrackError[itz]=sqMomError[itz]; //in ns
 399           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
 400           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 401           sumAllweightspi+=1./sqTrackError[itz];
 402           meantzeropi+=weightedtimezero[itz];
 403         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
 404
 405
 406         // Then, Combinatorial Algorithm
 407
 408         if(ntracksinsetmy<2 )break;
 409
 410         for (Int_t j=0; j<ntracksinsetmy; j++) {
 411           imass[j] = 3;
 412         }
 413
 414         //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
 415         Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
 416
 417         // Loop on mass hypotheses
 418         for (Int_t k=0; k < ncombinatorial;k++) {
 419           for (Int_t j=0; j<ntracksinsetmy; j++) {
 420             //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 421             imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
 422             texp[j]=exptof[j][imass[j]];
 423             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 424             //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100;
 425             if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
 426           }
 427
 428           Float_t sumAllweights=0.;
 429           Float_t meantzero=0.;
 430           Float_t eMeanTzero=0.;
 431
 432           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 433             sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 434
 435             timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 436
 437             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 438             sumAllweights+=1./sqTrackError[itz];
 439             meantzero+=weightedtimezero[itz];
 440
 441           } // end loop for (Int_t itz=0; itz<15;itz++)
 442
 443           meantzero=meantzero/sumAllweights; // it is given in [ns]
 444           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 445
 446           // calculate chisquare
 447
 448           Float_t chisquare=0.;
 449           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 450             //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 451             if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 452             else  chisquare+=1000;
 453           } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 454
 455           if(chisquare<=chisquarebest){
 456             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 457
 458               bestsqTrackError[iqsq]=sqTrackError[iqsq];
 459               besttimezero[iqsq]=timezero[iqsq];
 460               bestmomentum[iqsq]=momentum[iqsq];
 461               besttimeofflight[iqsq]=timeofflight[iqsq];
 462               besttexp[iqsq]=texp[iqsq];
 463               bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 464               //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 465               if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 466               else  bestchisquare[iqsq]=1000;
 467             }
 468
 469             Int_t npion=0;
 470             for (Int_t j=0; j<ntracksinsetmy; j++) {
 471               assparticle[j]=imass[j];
 472               if(imass[j] == 0) npion++;
 473             }
 474             npionbest=npion;
 475             chisquarebest=chisquare;
 476             t0best=meantzero;
 477             eT0best=eMeanTzero;
 478           } // close if(dummychisquare<=chisquare)
 479
 480         }
 481
 482         Double_t chi2cut[nmaxtracksinsetMax];
 483         chi2cut[0] = 0;
 484         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
 485         for (Int_t j=2; j<ntracksinset; j++) {
 486           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
 487         }
 488
 489         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
 490
 491 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
 492
 493         Bool_t kRedoT0 = kFALSE;
 494         ntracksinsetmyCut = ntracksinsetmy;
 495         Bool_t usetrack[nmaxtracksinsetMax];
 496         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 497           usetrack[icsq] = kTRUE;
 498           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
 499             kRedoT0 = kTRUE;
 500             ntracksinsetmyCut--;
 501             usetrack[icsq] = kFALSE;
 502           }
 503         } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 504
 505         //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
 506
 507         // Loop on mass hypotheses Redo
 508         if(kRedoT0 && ntracksinsetmyCut > 1){
 509           //      printf("Redo T0\n");
 510           for (Int_t k=0; k < ncombinatorial;k++) {
 511             for (Int_t j=0; j<ntracksinsetmy; j++) {
 512               //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 513               imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
 514               texp[j]=exptof[j][imass[j]];
 515               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 516               //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100;
 517               if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
 518             }
 519
 520             Float_t sumAllweights=0.;
 521             Float_t meantzero=0.;
 522             Float_t eMeanTzero=0.;
 523
 524             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 525               if(! usetrack[itz]) continue;
 526               sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 527
 528               timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 529
 530               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 531               sumAllweights+=1./sqTrackError[itz];
 532               meantzero+=weightedtimezero[itz];
 533
 534             } // end loop for (Int_t itz=0; itz<15;itz++)
 535
 536             meantzero=meantzero/sumAllweights; // it is given in [ns]
 537             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 538
 539             // calculate chisquare
 540
 541             Float_t chisquare=0.;
 542             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 543               if(! usetrack[icsq]) continue;
 544               //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 545               if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 546               else  chisquare+=1000;
 547
 548             } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 549
 550             Int_t npion=0;
 551             for (Int_t j=0; j<ntracksinsetmy; j++) {
 552               assparticle[j]=imass[j];
 553               if(imass[j] == 0) npion++;
 554             }
 555
 556             if(chisquare<=chisquarebest && npion>0){
 557               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 558                 if(! usetrack[iqsq]) continue;
 559                 bestsqTrackError[iqsq]=sqTrackError[iqsq];
 560                 besttimezero[iqsq]=timezero[iqsq];
 561                 bestmomentum[iqsq]=momentum[iqsq];
 562                 besttimeofflight[iqsq]=timeofflight[iqsq];
 563                 besttexp[iqsq]=texp[iqsq];
 564                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 565                 //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 566                 if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 567                 else  bestchisquare[iqsq]=1000;
 568               }
 569
 570               npionbest=npion;
 571               chisquarebest=chisquare;
 572               t0best=meantzero;
 573               eT0best=eMeanTzero;
 574             } // close if(dummychisquare<=chisquare)
 575
 576           }
 577         }
 578
 579         // filling histos
 580         Float_t confLevel=999;
 581
 582         // Sets with decent chisquares
 583
 584         if(chisquarebest<999.){
 585           Double_t dblechisquare=(Double_t)chisquarebest;
 586           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
 587 //        cout << " Set Number " << nsets << endl;
 588 //        cout << "Best Assignment, selection " << assparticle[0] <<
 589 //          assparticle[1] << assparticle[2] <<
 590 //          assparticle[3] << assparticle[4] <<
 591 //          assparticle[5] << endl;
 592 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
 593 //        cout << " C.L. of the set "<< confLevel <<endl;
 594 //        cout << " T0 for this set (in ns)  " << t0best << endl;
 595
 596           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
 597
 598             if(! usetrack[icsq]) continue;
 599
 600 //          cout << "Track # " << icsq  << " T0 offsets = "
 601 //               << besttimezero[icsq]-t0best <<
 602 //            " track error = "  << bestsqTrackError[icsq]
 603 //               << " Chisquare = " << bestchisquare[icsq]
 604 //               << " Momentum  = " << bestmomentum[icsq]
 605 //               << " TOF   = "     << besttimeofflight[icsq]
 606 //               << " TOF tracking  = " << besttexp[icsq]
 607 //               << " is used = " << usetrack[icsq] << endl;
 608           }
 609
 610           // Pick up only those with C.L. >1%
 611           //      if(confLevel>0.01 && ngoodsetsSel<200){
 612           if(confLevel>0.01 && ngoodsetsSel<200){
 613             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
 614             confLevelbestSel[ngoodsetsSel]=confLevel;
 615             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
 616             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
 617             t0bestallSel += t0best/eT0best/eT0best;
 618             sumWt0bestallSel += 1./eT0best/eT0best;
 619             ngoodsetsSel++;
 620             ngoodtrktrulyused+=ntracksinsetmyCut;
 621
 622             //      printf("t0 of a set = %f +/- %f\n",t0best,eT0best);
 623           }
 624           else{
 625             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
 626           }
 627         }
 628         //delete[] fTracksT0;
 629         fTracksT0->Clear();
 630         nsets++;
 631
 632       } // end for the current set
 633
 634       //Redo the computation of the best in order to esclude very bad samples
 635         if(ngoodsetsSel > 1){
 636             Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
 637             Int_t nsamples=ngoodsetsSel;
 638             ngoodsetsSel=0;
 639             t0bestallSel=0;
 640             sumWt0bestallSel=0;
 641             for (Int_t itz=0; itz<nsamples;itz++) {
 642                 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
 643                     t0bestallSel += t0bestSel[itz];
 644                     sumWt0bestallSel += eT0bestSel[itz];
 645                     ngoodsetsSel++;
 646 //            printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
 647                 }
 648                 else{
 649 //            printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
 650                 }
 651             }
 652         }
 653         if(ngoodsetsSel < 1){
 654             sumWt0bestallSel = 0.0;
 655         }
 656       //--------------------------------End recomputation
 657
 658       nUsedTracks =  ngoodtrkt0;
 659       if(strstr(option,"all")){
 660         if(sumAllweightspi>0.){
 661           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
 662           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
 663         }
 664
 665         if(sumWt0bestallSel>0){
 666           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
 667           eT0bestallSel = sqrt(1./sumWt0bestallSel);
 668
 669         }// end of if(sumWt0bestallSel>0){
 670
 671 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
 672       }
 673
 674       t0def=t0bestallSel;
 675       deltat0def=eT0bestallSel;
 676
 677       fT0SigmaT0def[0]=t0def;
 678       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
 679       fT0SigmaT0def[2]=ngoodtrkt0;
 680       fT0SigmaT0def[3]=ngoodtrktrulyused;
 681     }
 682   }
 683
 684   //delete [] fGTracks;
 685   fGTracks->Clear();
 686
 687   //   if(strstr(option,"tim") || strstr(option,"all")){
 688   //     cout << "AliTOFT0v1:" << endl ;
 689   //}
 690
 691   if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
 692
 693   bench->Stop("t0computation");
 694
 695   fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
 696   fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
 697
 698 //   bench->Print("t0computation");
 699 //   printf("%f %f\n",fT0SigmaT0def[4],fT0SigmaT0def[5]);
 700
 701
 702   delete bench;
 703   bench=NULL;
 704
 705   return fT0SigmaT0def;
 706   }
 707 //__________________________________________________________________
 708 Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
 709 {
 710   TBenchmark *bench=new TBenchmark();
 711   bench->Start("t0computation");
 712
 713   // Caluclate the Event Time using the ESD TOF time
 714
 715   fT0SigmaT0def[0]=0.;
 716   fT0SigmaT0def[1]=0.600;
 717   fT0SigmaT0def[2]=0.;
 718   fT0SigmaT0def[3]=0.;
 719
 720   const Int_t nmaxtracksinsetMax=10;
 721   Int_t nmaxtracksinset=10;
 722 //   if(strstr(option,"all")){
 723 //     cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
 724 //     cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
 725 //   }
 726
 727
 728   Int_t nsets=0;
 729   Int_t nUsedTracks=0;
 730   Int_t ngoodsetsSel= 0;
 731   Float_t t0bestSel[300];
 732   Float_t eT0bestSel[300];
 733   Float_t chiSquarebestSel[300];
 734   Float_t confLevelbestSel[300];
 735   Float_t t0bestallSel=0.;
 736   Float_t eT0bestallSel=0.;
 737   Float_t sumWt0bestallSel=0.;
 738   Float_t eMeanTzeroPi=0.;
 739   Float_t meantzeropi=0.;
 740   Float_t sumAllweightspi=0.;
 741   Double_t t0def=-999;
 742   Double_t deltat0def=999;
 743   Int_t ngoodtrktrulyused=0;
 744   Int_t ntracksinsetmyCut = 0;
 745
 746   Int_t ntrk=fEvent->GetNumberOfTracks();
 747
 748   //AliESDtrack **fTracks=new AliESDtrack*[ntrk];
 749   Int_t ngoodtrk=0;
 750   Int_t ngoodtrkt0 =0;
 751   Float_t mintime =1E6;
 752
 753   // First Track loop, Selection of good tracks
 754
 755   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 756     AliESDtrack *t=fEvent->GetTrack(itrk);
 757     Double_t momOld=t->GetP();
 758     Double_t mom=momOld-0.0036*momOld;
 759     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 760     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 761     Double_t time=t->GetTOFsignal();
 762
 763     time*=1.E-3; // tof given in nanoseconds
 764     if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
 765
 766     if (!AcceptTrack(t)) continue;
 767
 768     if(t->GetIntegratedLength() < 350)continue; //skip decays
 769     if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
 770     if(time <= mintime) mintime=time;
 771     //fTracks[ngoodtrk]=t;
 772     fTracks->AddLast(t);
 773     ngoodtrk++;
 774   }
 775
 776   if(ngoodtrk>22) nmaxtracksinset = 6;
 777
 778
 779 //     cout << " N. of ESD tracks                    : " << ntrk << endl;
 780 //     cout << " N. of preselected tracks            : " << ngoodtrk << endl;
 781 //     cout << " Minimum tof time in set (in ns)                 : " << mintime << endl;
 782
 783   //AliESDtrack **fGTracks=new AliESDtrack*[ngoodtrk];
 784
 785   //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
 786   for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
 787     //AliESDtrack *t=fTracks[jtrk];
 788     AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
 789     Double_t time=t->GetTOFsignal();
 790
 791     if((time-mintime*1.E3)<50.E3){ // For pp and per
 792       //fGTracks[ngoodtrkt0]=t;
 793       fGTracks->AddLast(t);
 794       ngoodtrkt0++;
 795     }
 796   }
 797
 798   //delete [] fTracks;
 799   fTracks->Clear();
 800
 801 //     cout << " N. of preselected tracks t0         : " << ngoodtrkt0 << endl;
 802
 803   Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
 804   Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
 805   if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
 806
 807 //   cout << " N. of sets                        : " << nseteq << endl;
 808
 809   if(ngoodtrkt0<2){
 810 //     cout << "less than 2 tracks, skip event " << endl;
 811     t0def=-999;
 812     deltat0def=0.600;
 813     fT0SigmaT0def[0]=t0def;
 814     fT0SigmaT0def[1]=deltat0def;
 815     fT0SigmaT0def[2]=ngoodtrkt0;
 816     fT0SigmaT0def[3]=ngoodtrkt0;
 817     //goto finish;
 818   }
 819   if(ngoodtrkt0>=2){
 820   // Decide how many tracks in set
 821     Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
 822     Int_t nset=1;
 823
 824     if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
 825
 826     // Loop over selected sets
 827
 828     if(nset>=1){
 829       for (Int_t i=0; i< nset; i++) {
 830         //      printf("Set %i of %i\n",i+1,nset);
 831         Float_t t0best=999.;
 832         Float_t eT0best=999.;
 833         Float_t chisquarebest=99999.;
 834         Int_t npionbest=0;
 835
 836         Int_t ntracksinsetmy=0;
 837         //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset];
 838         for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 839           Int_t index = itrk+i*ntracksinset;
 840           //if(index < ngoodtrkt0){
 841           if(index < fGTracks->GetEntries()){
 842             //AliESDtrack *t=fGTracks[index];
 843             AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
 844             //fTracksT0[itrk]=t;
 845             fTracksT0->AddLast(t);
 846             ntracksinsetmy++;
 847           }
 848         }
 849
 850         // Analyse it
 851
 852         Int_t   assparticle[nmaxtracksinsetMax];
 853         Float_t exptof[nmaxtracksinsetMax][3];
 854         Float_t timeofflight[nmaxtracksinsetMax];
 855         Float_t momentum[nmaxtracksinsetMax];
 856         Float_t timezero[nmaxtracksinsetMax];
 857         Float_t weightedtimezero[nmaxtracksinsetMax];
 858         Float_t beta[nmaxtracksinsetMax];
 859         Float_t texp[nmaxtracksinsetMax];
 860         Float_t dtexp[nmaxtracksinsetMax];
 861         Float_t sqMomError[nmaxtracksinsetMax];
 862         Float_t sqTrackError[nmaxtracksinsetMax];
 863         Float_t massarray[3]={0.13957,0.493677,0.9382723};
 864         Float_t tracktoflen[nmaxtracksinsetMax];
 865         Float_t besttimezero[nmaxtracksinsetMax];
 866         Float_t besttexp[nmaxtracksinsetMax];
 867         Float_t besttimeofflight[nmaxtracksinsetMax];
 868         Float_t bestmomentum[nmaxtracksinsetMax];
 869         Float_t bestchisquare[nmaxtracksinsetMax];
 870         Float_t bestweightedtimezero[nmaxtracksinsetMax];
 871         Float_t bestsqTrackError[nmaxtracksinsetMax];
 872         Int_t imass[nmaxtracksinsetMax];
 873
 874         for (Int_t j=0; j<ntracksinset; j++) {
 875           assparticle[j] = 3;
 876           timeofflight[j] = 0;
 877           momentum[j] = 0;
 878           timezero[j] = 0;
 879           weightedtimezero[j] = 0;
 880           beta[j] = 0;
 881           texp[j] = 0;
 882           dtexp[j] = 0;
 883           sqMomError[j] = 0;
 884           sqTrackError[j] = 0;
 885           tracktoflen[j] = 0;
 886           besttimezero[j] = 0;
 887           besttexp[j] = 0;
 888           besttimeofflight[j] = 0;
 889           bestmomentum[j] = 0;
 890           bestchisquare[j] = 0;
 891           bestweightedtimezero[j] = 0;
 892           bestsqTrackError[j] = 0;
 893           imass[j] = 1;
 894         }
 895
 896         //for (Int_t j=0; j<ntracksinsetmy; j++) {
 897         for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
 898           //AliESDtrack *t=fTracksT0[j];
 899           AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
 900           Double_t momOld=t->GetP();
 901           Double_t mom=momOld-0.0036*momOld;
 902           Double_t time=t->GetTOFsignal();
 903
 904           time*=1.E-3; // tof given in nanoseconds
 905           Double_t exptime[10]; t->GetIntegratedTimes(exptime);
 906           Double_t toflen=t->GetIntegratedLength();
 907           toflen=toflen/100.; // toflen given in m
 908
 909           timeofflight[j]=time;
 910           tracktoflen[j]=toflen;
 911           exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
 912           exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
 913           exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
 914           momentum[j]=mom;
 915           assparticle[j]=3;
 916
 917         } //end  for (Int_t j=0; j<ntracksinsetmy; j++) {
 918
 919         for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 920           beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
 921                                        +momentum[itz]*momentum[itz]);
 922           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]));
 923           sqTrackError[itz]=sqMomError[itz]; //in ns
 924           timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
 925           weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 926           sumAllweightspi+=1./sqTrackError[itz];
 927           meantzeropi+=weightedtimezero[itz];
 928         } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
 929
 930
 931         // Then, Combinatorial Algorithm
 932
 933         if(ntracksinsetmy<2 )break;
 934
 935         for (Int_t j=0; j<ntracksinsetmy; j++) {
 936           imass[j] = 3;
 937         }
 938
 939         //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
 940         Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
 941
 942         // Loop on mass hypotheses
 943         for (Int_t k=0; k < ncombinatorial;k++) {
 944           for (Int_t j=0; j<ntracksinsetmy; j++) {
 945             //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
 946             imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
 947             texp[j]=exptof[j][imass[j]];
 948             dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
 949             //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100;
 950             if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
 951           }
 952
 953           Float_t sumAllweights=0.;
 954           Float_t meantzero=0.;
 955           Float_t eMeanTzero=0.;
 956
 957           for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
 958             sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
 959
 960             timezero[itz]=texp[itz]-timeofflight[itz];// in ns
 961
 962             weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 963             sumAllweights+=1./sqTrackError[itz];
 964             meantzero+=weightedtimezero[itz];
 965
 966           } // end loop for (Int_t itz=0; itz<15;itz++)
 967
 968           meantzero=meantzero/sumAllweights; // it is given in [ns]
 969           eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
 970
 971           // calculate chisquare
 972           Float_t chisquare=0.;
 973           for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
 974             //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 975             if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
 976             else  chisquare+=1000;
 977           } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 978
 979           if(chisquare<=chisquarebest){
 980             for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
 981
 982               bestsqTrackError[iqsq]=sqTrackError[iqsq];
 983               besttimezero[iqsq]=timezero[iqsq];
 984               bestmomentum[iqsq]=momentum[iqsq];
 985               besttimeofflight[iqsq]=timeofflight[iqsq];
 986               besttexp[iqsq]=texp[iqsq];
 987               bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 988               //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 989               if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
 990               else  bestchisquare[iqsq]=1000;
 991             }
 992
 993             Int_t npion=0;
 994             for (Int_t j=0; j<ntracksinsetmy; j++) {
 995               assparticle[j]=imass[j];
 996               if(imass[j] == 0) npion++;
 997             }
 998             npionbest=npion;
 999             chisquarebest=chisquare;
1000             t0best=meantzero;
1001             eT0best=eMeanTzero;
1002           } // close if(dummychisquare<=chisquare)
1003         }
1004
1005         Double_t chi2cut[nmaxtracksinsetMax];
1006         chi2cut[0] = 0;
1007         chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
1008         for (Int_t j=2; j<ntracksinset; j++) {
1009           chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
1010         }
1011
1012         Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
1013
1014 //      printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
1015
1016         Bool_t kRedoT0 = kFALSE;
1017         ntracksinsetmyCut = ntracksinsetmy;
1018         Bool_t usetrack[nmaxtracksinsetMax];
1019         for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1020           usetrack[icsq] = kTRUE;
1021           if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
1022               kRedoT0 = kTRUE;
1023               ntracksinsetmyCut--;
1024               usetrack[icsq] = kFALSE;
1025           }
1026         } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1027
1028 //      printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
1029
1030         // Loop on mass hypotheses Redo
1031         if(kRedoT0 && ntracksinsetmyCut > 1){
1032           //      printf("Redo T0\n");
1033           for (Int_t k=0; k < ncombinatorial;k++) {
1034             for (Int_t j=0; j<ntracksinsetmy; j++) {
1035               //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
1036               imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
1037               texp[j]=exptof[j][imass[j]];
1038               dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
1039               //if(! CheckTPCMatching(fTracksT0[j],imass[j])) dtexp[j]*=100;
1040               if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
1041             }
1042
1043             Float_t sumAllweights=0.;
1044             Float_t meantzero=0.;
1045             Float_t eMeanTzero=0.;
1046
1047             for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
1048               if(! usetrack[itz]) continue;
1049               sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
1050
1051               timezero[itz]=texp[itz]-timeofflight[itz];// in ns
1052
1053               weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
1054               sumAllweights+=1./sqTrackError[itz];
1055               meantzero+=weightedtimezero[itz];
1056
1057             } // end loop for (Int_t itz=0; itz<15;itz++)
1058
1059             meantzero=meantzero/sumAllweights; // it is given in [ns]
1060             eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
1061
1062             // calculate chisquare
1063
1064             Float_t chisquare=0.;
1065             for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
1066               if(! usetrack[icsq]) continue;
1067               //if(CheckTPCMatching(fTracksT0[icsq],imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
1068               if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
1069               else  chisquare+=1000;
1070             } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1071
1072             Int_t npion=0;
1073             for (Int_t j=0; j<ntracksinsetmy; j++) {
1074               assparticle[j]=imass[j];
1075               if(imass[j] == 0) npion++;
1076             }
1077
1078             if(chisquare<=chisquarebest && npion>0){
1079               for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
1080                 if(! usetrack[iqsq]) continue;
1081                 bestsqTrackError[iqsq]=sqTrackError[iqsq];
1082                 besttimezero[iqsq]=timezero[iqsq];
1083                 bestmomentum[iqsq]=momentum[iqsq];
1084                 besttimeofflight[iqsq]=timeofflight[iqsq];
1085                 besttexp[iqsq]=texp[iqsq];
1086                 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1087                 //if(CheckTPCMatching(fTracksT0[iqsq],imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
1088                 if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
1089                 else  bestchisquare[iqsq]=1000;
1090               }
1091
1092               npionbest=npion;
1093               chisquarebest=chisquare;
1094               t0best=meantzero;
1095               eT0best=eMeanTzero;
1096             } // close if(dummychisquare<=chisquare)
1097
1098           }
1099         }
1100
1101         // filling histos
1102         Float_t confLevel=999;
1103
1104         // Sets with decent chisquares
1105         //      printf("Chi2best of the set = %f \n",chisquarebest);
1106
1107         if(chisquarebest<999.){
1108           Double_t dblechisquare=(Double_t)chisquarebest;
1109           confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
1110 //        cout << " Set Number " << nsets << endl;
1111 //        cout << "Best Assignment, selection " << assparticle[0] <<
1112 //          assparticle[1] << assparticle[2] <<
1113 //          assparticle[3] << assparticle[4] <<
1114 //          assparticle[5] << endl;
1115 //        cout << " Chisquare of the set "<< chisquarebest <<endl;
1116 //        cout << " C.L. of the set "<< confLevel <<endl;
1117 //        cout << " T0 for this set (in ns)  " << t0best << endl;
1118
1119           Int_t ntrackincurrentsel=0;
1120           for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
1121
1122             if(! usetrack[icsq]) continue;
1123
1124 //          cout << "Track # " << icsq  << " T0 offsets = "
1125 //               << besttimezero[icsq]-t0best <<
1126 //            " track error = "  << bestsqTrackError[icsq]
1127 //               << " Chisquare = " << bestchisquare[icsq]
1128 //               << " Momentum  = " << bestmomentum[icsq]
1129 //               << " TOF   = "     << besttimeofflight[icsq]
1130 //               << " TOF tracking  = " << besttexp[icsq]
1131 //               << " is used = " << usetrack[icsq] << endl;
1132             ntrackincurrentsel++;
1133           }
1134
1135           //      printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
1136
1137           // Pick up only those with C.L. >1%
1138           //      if(confLevel>0.01 && ngoodsetsSel<200){
1139           if(confLevel>0.01 && ngoodsetsSel<200){
1140             chiSquarebestSel[ngoodsetsSel]=chisquarebest;
1141             confLevelbestSel[ngoodsetsSel]=confLevel;
1142             t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
1143             eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
1144             t0bestallSel += t0best/eT0best/eT0best;
1145             sumWt0bestallSel += 1./eT0best/eT0best;
1146             ngoodsetsSel++;
1147             ngoodtrktrulyused+=ntracksinsetmyCut;
1148 //          printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
1149           }
1150           else{
1151             //      printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
1152           }
1153         }
1154         //delete[] fTracksT0;
1155         fTracksT0->Clear();
1156         nsets++;
1157
1158       } // end for the current set
1159
1160       //Redo the computation of the best in order to esclude very bad samples
1161         if(ngoodsetsSel > 1){
1162             Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel;
1163             Int_t nsamples=ngoodsetsSel;
1164             ngoodsetsSel=0;
1165             t0bestallSel=0;
1166             sumWt0bestallSel=0;
1167             for (Int_t itz=0; itz<nsamples;itz++) {
1168                 if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
1169                     t0bestallSel += t0bestSel[itz];
1170                     sumWt0bestallSel += eT0bestSel[itz];
1171                     ngoodsetsSel++;
1172 //            printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
1173                 }
1174                 else{
1175 //            printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
1176                 }
1177             }
1178         }
1179         if(ngoodsetsSel < 1){
1180             sumWt0bestallSel = 0.0;
1181         }
1182       //--------------------------------End recomputation
1183
1184       nUsedTracks =  ngoodtrkt0;
1185       if(strstr(option,"all")){
1186         if(sumAllweightspi>0.){
1187           meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
1188           eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
1189         }
1190
1191         //      printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel);
1192
1193         if(sumWt0bestallSel>0){
1194           t0bestallSel  = t0bestallSel/sumWt0bestallSel;
1195           eT0bestallSel = sqrt(1./sumWt0bestallSel);
1196           //      printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
1197         }// end of if(sumWt0bestallSel>0){
1198
1199 //      cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
1200       }
1201
1202       t0def=t0bestallSel;
1203       deltat0def=eT0bestallSel;
1204
1205       fT0SigmaT0def[0]=t0def;
1206       fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
1207       fT0SigmaT0def[2]=ngoodtrkt0;
1208       fT0SigmaT0def[3]=ngoodtrktrulyused;
1209     }
1210   }
1211
1212   //delete [] fGTracks;
1213   fGTracks->Clear();
1214
1215   //   if(strstr(option,"tim") || strstr(option,"all")){
1216   //     cout << "AliTOFT0v1:" << endl ;
1217   //}
1218
1219   if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
1220
1221   bench->Stop("t0computation");
1222
1223   fT0SigmaT0def[4]=bench->GetRealTime("t0computation");
1224   fT0SigmaT0def[5]=bench->GetCpuTime("t0computation");
1225
1226 //   bench->Print("t0computation");
1227 //   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]));
1228
1229   delete bench;
1230   bench=NULL;
1231
1232   return fT0SigmaT0def;
1233   }
1234 //__________________________________________________________________
1235 Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
1236 {
1237   // Take the error extimate for the TOF time in the track reconstruction
1238
1239   static const Double_t kMasses[]={
1240     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
1241   };
1242
1243   Double_t mass=kMasses[index+2];
1244
1245   Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass);
1246
1247   return sigma;
1248 }
1249
1250 //__________________________________________________________________
1251 Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track)
1252 {
1253
1254   /* TPC refit */
1255   if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1256   /* do not accept kink daughters */
1257   if (track->GetKinkIndex(0)>0) return kFALSE;
1258   /* N clusters TPC */
1259   if (track->GetTPCclusters(0) < 50) return kFALSE;
1260   /* chi2 TPC */
1261   if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE;
1262   /* sigma to vertex */
1263   if (GetSigmaToVertex(track) > 4.) return kFALSE;
1264
1265   /* accept track */
1266   return kTRUE;
1267
1268 }
1269
1270 //____________________________________________________________________
1271 Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const
1272 {
1273   // Calculates the number of sigma to the vertex.
1274
1275   Float_t b[2];
1276   Float_t bRes[2];
1277   Float_t bCov[3];
1278   esdTrack->GetImpactParameters(b,bCov);
1279
1280   if (bCov[0]<=0 || bCov[2]<=0) {
1281     bCov[0]=0; bCov[2]=0;
1282   }
1283   bRes[0] = TMath::Sqrt(bCov[0]);
1284   bRes[1] = TMath::Sqrt(bCov[2]);
1285
1286   // -----------------------------------
1287   // How to get to a n-sigma cut?
1288   //
1289   // The accumulated statistics from 0 to d is
1290   //
1291   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1292   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1293   //
1294   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1295   // Can this be expressed in a different way?
1296
1297   if (bRes[0] == 0 || bRes[1] ==0)
1298     return -1;
1299
1300   //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1301   Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2));
1302
1303   // work around precision problem
1304   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1305   // 1e-15 corresponds to nsigma ~ 7.7
1306   if (TMath::Exp(-d * d / 2) < 1e-15)
1307     return 1000;
1308
1309   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1310   return nSigma;
1311 }
1312 //____________________________________________________________________
1313
1314 Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{
1315     Bool_t status = kFALSE;
1316
1317     Double_t exptimes[5];
1318     track->GetIntegratedTimes(exptimes);
1319
1320     Float_t dedx = track->GetTPCsignal();
1321
1322     Double_t ptpc[3];
1323     track->GetInnerPxPyPz(ptpc);
1324     Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
1325
1326     if(imass > 2 || imass < 0) return status;
1327     Int_t i = imass+2;
1328
1329     AliPID::EParticleType type=AliPID::EParticleType(i);
1330
1331     Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type);
1332     Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
1333
1334     if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){
1335         status = kTRUE;
1336     }
1337
1338     return status;
1339 }
1340
1341 //____________________________________________________________________
1342 Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const
1343 {
1344   //
1345   // Returns base^exponent
1346   //
1347
1348   Float_t power=1.;
1349
1350   for (Int_t ii=exponent; ii>0; ii--)
1351       power=power*base;
1352
1353   return power;
1354
1355 }
1356 //____________________________________________________________________
1357 Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const
1358 {
1359   //
1360   // Returns base^exponent
1361   //
1362
1363   Int_t power=1;
1364
1365   for (Int_t ii=exponent; ii>0; ii--)
1366       power=power*base;
1367
1368   return power;
1369
1370 }