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