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