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