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