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