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