X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TOF%2FAliTOFT0v1.cxx;h=701e110839e50d3d3ed630f8ad13cd7de8f7c3b0;hb=73cf8fa924039c2b510eccc6fcc9a965d2c93724;hp=47c3f02cc966b4ca5f2d52bc872dc4136a88d349;hpb=10d6f841a0d30b3326455ae8486d38837537e7d1;p=u%2Fmrichter%2FAliRoot.git diff --git a/TOF/AliTOFT0v1.cxx b/TOF/AliTOFT0v1.cxx index 47c3f02cc96..701e110839e 100644 --- a/TOF/AliTOFT0v1.cxx +++ b/TOF/AliTOFT0v1.cxx @@ -13,28 +13,18 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id: AliTOFT0v1.cxx,v 1.8 2003/10/23 16:32:20 hristov Exp $ */ +/* $Id: AliTOFT0v1.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */ //_________________________________________________________________________ // This is a TTask that made the calculation of the Time zero using TOF. // Description: The algorithm used to calculate the time zero of interaction // using TOF detector is the following. -// We select in the MonteCarlo some primary particles - or tracks in the following - +// We select in the ESD some "primary" particles - or tracks in the following - // that strike the TOF detector (the larger part are pions, kaons or protons). // We choose a set of 10 selected tracks, for each track You have the length -// of the track when the TOF is reached (a standard TOF hit does not contain this -// additional information, this is the reason why we implemented a new time zero -// dedicated TOF hit class AliTOFhitT0; in order to store this type of hit You -// have to use the AliTOFv4T0 as TOF class in Your Config.C. In AliTOFv4T0 the -// StepManager was modified in order to fill the TOF hit branch with this type -// of hits; in fact the AliTOF::AddT0Hit is called rather that the usual AliTOF::AddHit), -// the momentum at generation (from TreeK) and the time of flight +// of the track when the TOF is reached, +// the momentum and the time of flight // given by the TOF detector. -// (Observe that the ctor of the AliTOF class, when the AliTOFv4T0 class is used, is called -// with the "tzero" option: it is in order create the fHits TClonesArray filled with -// AliTOFhitT0 objects, rather than with normal AliTOFhit) -// Then Momentum and time of flight for each track are smeared according to -// known experimental resolution (all sources of error have been token into account). // Let consider now only one set of 10 tracks (the algorithm is the same for all sets). // Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton, // we consider all the 3 at 10 possible cases. @@ -62,94 +52,192 @@ // assignment // Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions //-- Author: F. Pierella -//-- Mod By SA. +//-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni +//-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table +// for Power3) ////////////////////////////////////////////////////////////////////////////// -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "AliTOFT0v1.h" -#include "AliMC.h" #include "AliESDtrack.h" -#include "AliESD.h" -#include -#include "AliESDtrackCuts.h" -#include "AliTOFcalibHisto.h" +#include "AliESDEvent.h" +#include "AliTOFT0v1.h" +#include "TBenchmark.h" +#include "AliPID.h" +#include "AliESDpid.h" ClassImp(AliTOFT0v1) //____________________________________________________________________________ -AliTOFT0v1::AliTOFT0v1(AliESDEvent* event) +AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID): + TObject(), + fLowerMomBound(0.5), + fUpperMomBound(3), + fTimeCorr(0.), + fEvent(0x0), + fPIDesd(extPID), + fTracks(new TObjArray(10)), + fGTracks(new TObjArray(10)), + fTracksT0(new TObjArray(10)), + fOptFlag(kFALSE) { - fLowerMomBound=0.4; // [GeV/c] default value pp - fUpperMomBound=2.0 ; // [GeV/c] default value pp - fTimeResolution = 0.80e-10; // 80 ps by default - fTimeCorr = 0.0; // in ns by default - fLOffset=0.0; - fT0Offset=0.0; - fEvent=event; - - fT0SigmaT0def[0]=-999.; - fT0SigmaT0def[1]=999.; - fT0SigmaT0def[2]=-999.; - fT0SigmaT0def[3]=-999.; + // + // default constructor + // + if(AliPID::ParticleMass(0) == 0) new AliPID(); - fDeltaTfromMisallinement = 0.; // in ps + if(!fPIDesd){ + fPIDesd = new AliESDpid(); + } + + Init(NULL); + + //initialise lookup table for power 3 + // a set should only have 10 tracks a t maximum + // so up to 15 should be more than enough + for (Int_t i=0; i<15; ++i) { + fLookupPowerThree[i]=ToCalculatePower(3,i); + } } //____________________________________________________________________________ -AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero) +AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID): + TObject(), + fLowerMomBound(0.5), + fUpperMomBound(3.0), + fTimeCorr(0.), + fEvent(event), + fPIDesd(extPID), + fTracks(new TObjArray(10)), + fGTracks(new TObjArray(10)), + fTracksT0(new TObjArray(10)), + fOptFlag(kFALSE) +{ + // + // real constructor + // + if(AliPID::ParticleMass(0) == 0) new AliPID(); + + if(!fPIDesd){ + fPIDesd = new AliESDpid(); + } + + Init(event); + //initialise lookup table for power 3 + for (Int_t i=0; i<15; ++i) { + fLookupPowerThree[i]=Int_t(TMath::Power(3,i)); + } +} +//____________________________________________________________________________ +AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero) { - ( (AliTOFT0v1 &)tzero ).Copy(*this); + // + // assign. operator + // + + if (this == &tzero) + return *this; + + fLowerMomBound=tzero.fLowerMomBound; + fUpperMomBound=tzero.fUpperMomBound; + fTimeCorr=tzero.fTimeCorr; + fEvent=tzero.fEvent; + fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0]; + fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1]; + fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2]; + fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3]; + + fTracks=tzero.fTracks; + fGTracks=tzero.fGTracks; + fTracksT0=tzero.fTracksT0; + + for (Int_t ii=0; iiGetEntries(); ii++) + fTracks->AddLast(tzero.fTracks->At(ii)); + + for (Int_t ii=0; iiGetEntries(); ii++) + fGTracks->AddLast(tzero.fGTracks->At(ii)); + + for (Int_t ii=0; iiGetEntries(); ii++) + fTracksT0->AddLast(tzero.fTracksT0->At(ii)); + + fOptFlag=tzero.fOptFlag; + + return *this; } //____________________________________________________________________________ AliTOFT0v1::~AliTOFT0v1() { // dtor + fEvent=NULL; + if (fTracks) { + fTracks->Clear(); + delete fTracks; + fTracks=0x0; + } + + if (fGTracks) { + fGTracks->Clear(); + delete fGTracks; + fGTracks=0x0; + } + + if (fTracksT0) { + fTracksT0->Clear(); + delete fTracksT0; + fTracksT0=0x0; + } + + } -void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){ - fTimeResolution=timeresolution; +//____________________________________________________________________________ + +void +AliTOFT0v1::Init(AliESDEvent *event) +{ + + /* + * init + */ + + fEvent = event; + fT0SigmaT0def[0]=0.; + fT0SigmaT0def[1]=0.6; + fT0SigmaT0def[2]=0.; + fT0SigmaT0def[3]=0.; + } -//____________________________________________________________________________ -//____________________________________________________________________________ -Double_t * AliTOFT0v1::DefineT0(Option_t *option) +//____________________________________________________________________________ +Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut) { - Float_t timeresolutioninns=fTimeResolution*(1.e+9) * TMath::Sqrt(2.); // convert in [ns] + TBenchmark *bench=new TBenchmark(); + bench->Start("t0computation"); + + // Caluclate the Event Time using the ESD TOF time + + fT0SigmaT0def[0]=0.; + fT0SigmaT0def[1]=0.600; + fT0SigmaT0def[2]=0.; + fT0SigmaT0def[3]=0.; - const Int_t nmaxtracksinset=10; - if(strstr(option,"all")){ - cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl; - cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl; - } + const Int_t nmaxtracksinsetMax=10; + Int_t nmaxtracksinset=10; +// if(strstr(option,"all")){ +// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl; +// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl; +// } Int_t nsets=0; - Int_t NusedTracks=0; + Int_t nUsedTracks=0; Int_t ngoodsetsSel= 0; Float_t t0bestSel[300]; - Float_t Et0bestSel[300]; - Float_t ChisquarebestSel[300]; - Float_t ConfLevelbestSel[300]; + Float_t eT0bestSel[300]; + Float_t chiSquarebestSel[300]; + Float_t confLevelbestSel[300]; Float_t t0bestallSel=0.; - Float_t Et0bestallSel=0.; + Float_t eT0bestallSel=0.; Float_t sumWt0bestallSel=0.; - Float_t Emeantzeropi=0.; + Float_t eMeanTzeroPi=0.; Float_t meantzeropi=0.; Float_t sumAllweightspi=0.; Double_t t0def=-999; @@ -159,76 +247,56 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) Int_t ntrk=fEvent->GetNumberOfTracks(); - AliESDtrack **tracks=new AliESDtrack*[ntrk]; Int_t ngoodtrk=0; Int_t ngoodtrkt0 =0; - Float_t mintime =1E6; + Float_t meantime =0; // First Track loop, Selection of good tracks + fTracks->Clear(); for (Int_t itrk=0; itrkGetTrack(itrk); Double_t momOld=t->GetP(); Double_t mom=momOld-0.0036*momOld; if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue; + if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue; Double_t time=t->GetTOFsignal(); time*=1.E-3; // tof given in nanoseconds if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue; - -#if 0 - AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts(); - Bool_t tpcRefit = kTRUE; - Double_t nSigma = 4; - Bool_t sigmaToVertex = kTRUE; - esdTrackCuts->SetRequireSigmaToVertex(sigmaToVertex); - if (sigmaToVertex) { - esdTrackCuts->SetMaxNsigmaToVertex(nSigma); - } - else{ - esdTrackCuts->SetMaxDCAToVertexZ(3.0); - esdTrackCuts->SetMaxDCAToVertexXY(3.0); - } - esdTrackCuts->SetRequireTPCRefit(tpcRefit); - esdTrackCuts->SetAcceptKinkDaughters(kFALSE); - esdTrackCuts->SetMinNClustersTPC(50); - esdTrackCuts->SetMaxChi2PerClusterTPC(3.5); - Bool_t accepted; - accepted=esdTrackCuts->AcceptTrack(t); - if(!accepted) continue; -#endif - - if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays - if(time <= mintime) mintime=time; - tracks[ngoodtrk]=t; + + if (!AcceptTrack(t)) continue; + + if(t->GetIntegratedLength() < 350)continue; //skip decays + if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue; + + meantime+=time; + fTracks->AddLast(t); ngoodtrk++; } - - - cout << " T0 offset selected for this sample : " << fT0Offset << endl; - cout << " N. of ESD tracks : " << ntrk << endl; - cout << " N. of preselected tracks : " << ngoodtrk << endl; - cout << " Minimum tof time in set (in ns) : " << mintime << endl; - - AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk]; - - for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) { - AliESDtrack *t=tracks[jtrk]; - Double_t time=t->GetTOFsignal(); - if((time-mintime*1.E3)<50.E3){ // For pp and per - gtracks[ngoodtrkt0]=t; - ngoodtrkt0++; - } + if(ngoodtrk > 1) meantime /= ngoodtrk; + + if(ngoodtrk>22) nmaxtracksinset = 6; + + fGTracks->Clear(); + for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) { + AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk); + // Double_t time=t->GetTOFsignal(); + // if((time-meantime*1.E3)<50.E3){ // For pp and per + fGTracks->AddLast(t); + ngoodtrkt0++; + // } } - + + fTracks->Clear(); Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1; Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq; if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++; + if(ngoodtrkt0<2){ - cout << "less than 2 tracks, skip event " << endl; t0def=-999; deltat0def=0.600; fT0SigmaT0def[0]=t0def; @@ -239,59 +307,56 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) } if(ngoodtrkt0>=2){ // Decide how many tracks in set - Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent); + Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent); Int_t nset=1; if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} - - cout << " number of sets = " << nset << endl; - - if(strstr(option,"tim") || strstr(option,"all"))gBenchmark->Start("TOFT0v1"); - + // Loop over selected sets if(nset>=1){ for (Int_t i=0; i< nset; i++) { - + // printf("Set %i of %i\n",i+1,nset); Float_t t0best=999.; - Float_t Et0best=999.; + Float_t eT0best=999.; Float_t chisquarebest=99999.; Int_t npionbest=0; + fTracksT0->Clear(); Int_t ntracksinsetmy=0; - AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset]; for (Int_t itrk=0; itrkGetEntries()){ + AliESDtrack *t=(AliESDtrack*)fGTracks->At(index); + fTracksT0->AddLast(t); ntracksinsetmy++; } } - + // Analyse it - Int_t assparticle[nmaxtracksinset]; - Float_t exptof[nmaxtracksinset][3]; - Float_t timeofflight[nmaxtracksinset]; - Float_t momentum[nmaxtracksinset]; - Float_t timezero[nmaxtracksinset]; - Float_t weightedtimezero[nmaxtracksinset]; - Float_t beta[nmaxtracksinset]; - Float_t texp[nmaxtracksinset]; - Float_t dtexp[nmaxtracksinset]; - Float_t sqMomError[nmaxtracksinset]; - Float_t sqTrackError[nmaxtracksinset]; + Int_t assparticle[nmaxtracksinsetMax]; + Float_t exptof[nmaxtracksinsetMax][3]; + Float_t momErr[nmaxtracksinsetMax][3]; + Float_t timeofflight[nmaxtracksinsetMax]; + Float_t momentum[nmaxtracksinsetMax]; + Float_t timezero[nmaxtracksinsetMax]; + Float_t weightedtimezero[nmaxtracksinsetMax]; + Float_t beta[nmaxtracksinsetMax]; + Float_t texp[nmaxtracksinsetMax]; + Float_t dtexp[nmaxtracksinsetMax]; + Float_t sqMomError[nmaxtracksinsetMax]; + Float_t sqTrackError[nmaxtracksinsetMax]; Float_t massarray[3]={0.13957,0.493677,0.9382723}; - Float_t tracktoflen[nmaxtracksinset]; - Float_t besttimezero[nmaxtracksinset]; - Float_t besttexp[nmaxtracksinset]; - Float_t besttimeofflight[nmaxtracksinset]; - Float_t bestmomentum[nmaxtracksinset]; - Float_t bestchisquare[nmaxtracksinset]; - Float_t bestweightedtimezero[nmaxtracksinset]; - Float_t bestsqTrackError[nmaxtracksinset]; - Int_t imass[nmaxtracksinset]; + Float_t tracktoflen[nmaxtracksinsetMax]; + Float_t besttimezero[nmaxtracksinsetMax]; + Float_t besttexp[nmaxtracksinsetMax]; + Float_t besttimeofflight[nmaxtracksinsetMax]; + Float_t bestmomentum[nmaxtracksinsetMax]; + Float_t bestchisquare[nmaxtracksinsetMax]; + Float_t bestweightedtimezero[nmaxtracksinsetMax]; + Float_t bestsqTrackError[nmaxtracksinsetMax]; + Int_t imass[nmaxtracksinsetMax]; for (Int_t j=0; jGetEntries(); j++) { + AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j); Double_t momOld=t->GetP(); Double_t mom=momOld-0.0036*momOld; Double_t time=t->GetTOFsignal(); time*=1.E-3; // tof given in nanoseconds - Double_t exptime[10]; t->GetIntegratedTimes(exptime); + Double_t exptime[AliPID::kSPECIESC]; + t->GetIntegratedTimes(exptime,AliPID::kSPECIESC); Double_t toflen=t->GetIntegratedLength(); toflen=toflen/100.; // toflen given in m - timeofflight[j]=time+fT0Offset; - tracktoflen[j]=toflen+fLOffset; + timeofflight[j]=time; + tracktoflen[j]=toflen; exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns exptof[j][1]=exptime[3]*1.E-3+fTimeCorr; exptof[j][2]=exptime[4]*1.E-3+fTimeCorr; momentum[j]=mom; assparticle[j]=3; - - } //end for (Int_t j=0; j %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); + // printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); Bool_t kRedoT0 = kFALSE; ntracksinsetmyCut = ntracksinsetmy; - Bool_t usetrack[nmaxtracksinset]; + Bool_t usetrack[nmaxtracksinsetMax]; for (Int_t icsq=0; icsq chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){ - kRedoT0 = kTRUE; - ntracksinsetmyCut--; - usetrack[icsq] = kFALSE; + kRedoT0 = kTRUE; + ntracksinsetmyCut--; + usetrack[icsq] = kFALSE; + // printf("tracks chi2 = %f\n",bestchisquare[icsq]); } } // end loop for (Int_t icsq=0; icsq<15;icsq++) - printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut); - // Loop on mass hypotheses Redo if(kRedoT0 && ntracksinsetmyCut > 1){ - printf("Redo T0\n"); + // printf("Redo T0\n"); for (Int_t k=0; k < ncombinatorial;k++) { for (Int_t j=0; j0){ for(Int_t iqsq = 0; iqsq= 999){ - printf("How is it possible (chi2 = %f)? T0best = %f\n",chisquarebest,t0best); - - for(Int_t icsq=0; icsq1% - // if(confLevel>0.01 && ngoodsetsSel<200){ if(confLevel>0.01 && ngoodsetsSel<200){ - ChisquarebestSel[ngoodsetsSel]=chisquarebest; - ConfLevelbestSel[ngoodsetsSel]=confLevel; - t0bestSel[ngoodsetsSel]=t0best/Et0best/Et0best; - Et0bestSel[ngoodsetsSel]=1./Et0best/Et0best; - t0bestallSel += t0best/Et0best/Et0best; - sumWt0bestallSel += 1./Et0best/Et0best; + chiSquarebestSel[ngoodsetsSel]=chisquarebest; + confLevelbestSel[ngoodsetsSel]=confLevel; + t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best; + eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best; + t0bestallSel += t0best/eT0best/eT0best; + sumWt0bestallSel += 1./eT0best/eT0best; ngoodsetsSel++; - ngoodtrktrulyused+=ntracksinsetmyCut; + ngoodtrktrulyused+=ntracksinsetmyCut; + // printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel); } else{ - printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy); + // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy); } } - delete[] tracksT0; + fTracksT0->Clear(); nsets++; } // end for the current set - NusedTracks = ngoodtrkt0; + //Redo the computation of the best in order to esclude very bad samples + if(ngoodsetsSel > 1){ + Double_t t0BestStep1 = t0bestallSel/sumWt0bestallSel; + Int_t nsamples=ngoodsetsSel; + ngoodsetsSel=0; + t0bestallSel=0; + sumWt0bestallSel=0; + for (Int_t itz=0; itz0.){ meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns] - Emeantzeropi=sqrt(1./sumAllweightspi); // it is given in [ns] + eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns] } + // printf("t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,sumWt0bestallSel); + if(sumWt0bestallSel>0){ t0bestallSel = t0bestallSel/sumWt0bestallSel; - Et0bestallSel = sqrt(1./sumWt0bestallSel); - + eT0bestallSel = sqrt(1./sumWt0bestallSel); + // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel); }// end of if(sumWt0bestallSel>0){ - cout << "T0 all " << t0bestallSel << " +/- " << Et0bestallSel << "Number of tracks used: "<Stop("TOFT0v1"); - cout << "AliTOFT0v1:" << endl ; - cout << " took " << gBenchmark->GetCpuTime("TOFT0v1") << " seconds in order to calculate T0 " << endl ; - } - printf("T0 from TOF = %f ns\n",fT0SigmaT0def[0]); - + + fGTracks->Clear(); + + if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6; + + bench->Stop("t0computation"); + + fT0SigmaT0def[4]=bench->GetRealTime("t0computation"); + fT0SigmaT0def[5]=bench->GetCpuTime("t0computation"); + +// bench->Print("t0computation"); +// printf("(%4.1f < p < %4.1f GeV/c) T0-TOF =%9.1f +/- %5.1f ps (n_track = %i)\n\n",pMinCut,pMaxCut,-fT0SigmaT0def[0]*1000,fT0SigmaT0def[1]*1000,Int_t(fT0SigmaT0def[3])); + + delete bench; + bench=NULL; + return fT0SigmaT0def; } //__________________________________________________________________ -Double_t * AliTOFT0v1::DefineT0RawCorrection(Option_t *option) -{ - Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns] - - const Int_t nmaxtracksinset=10; - if(strstr(option,"all")){ - cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl; - cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl; - } - - Float_t stripmean = 0; - - Int_t nsets=0; - Int_t NusedTracks=0; - Int_t ngoodsetsSel= 0; - Float_t t0bestSel[300]; - Float_t Et0bestSel[300]; - Float_t ChisquarebestSel[300]; - Float_t ConfLevelbestSel[300]; - Float_t t0bestallSel=0.; - Float_t Et0bestallSel=0.; - Float_t sumWt0bestallSel=0.; - Float_t Emeantzeropi=0.; - Float_t meantzeropi=0.; - Float_t sumAllweightspi=0.; - Double_t t0def=-999; - Double_t deltat0def=999; - Int_t ngoodtrktrulyused=0; - Int_t ntracksinsetmyCut = 0; +Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const +{ + // Take the error extimate for the TOF time in the track reconstruction - Int_t ntrk=fEvent->GetNumberOfTracks(); - - AliESDtrack **tracks=new AliESDtrack*[ntrk]; - Int_t ngoodtrk=0; - Int_t ngoodtrkt0 =0; - Float_t mintime =1E6; - - // First Track loop, Selection of good tracks + static const Double_t kMasses[]={ + 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 + }; - for (Int_t itrk=0; itrkGetTrack(itrk); - Double_t momOld=t->GetP(); - Double_t mom=momOld-0.0036*momOld; - if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue; - Double_t tot = t->GetTOFsignalToT(); - Double_t time=t->GetTOFsignalRaw(); - Int_t chan = t->GetTOFCalChannel(); - Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0); - time -= corr*1000.; - - Int_t crate = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(crate == 63 || crate == 62){ - time += 9200; - } - - Int_t strip = fCalib->GetCalibMap(AliTOFcalibHisto::kSectorStrip,chan); + Double_t mass=kMasses[index+2]; - time*=1.E-3; // tof given in nanoseconds - if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue; - -#if 0 - AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts(); - Bool_t tpcRefit = kTRUE; - Double_t nSigma = 4; - Bool_t sigmaToVertex = kTRUE; - esdTrackCuts->SetRequireSigmaToVertex(sigmaToVertex); - if (sigmaToVertex) { - esdTrackCuts->SetMaxNsigmaToVertex(nSigma); - } - else{ - esdTrackCuts->SetMaxDCAToVertexZ(3.0); - esdTrackCuts->SetMaxDCAToVertexXY(3.0); - } - esdTrackCuts->SetRequireTPCRefit(tpcRefit); - esdTrackCuts->SetAcceptKinkDaughters(kFALSE); - esdTrackCuts->SetMinNClustersTPC(50); - esdTrackCuts->SetMaxChi2PerClusterTPC(3.5); - Bool_t accepted; - accepted=esdTrackCuts->AcceptTrack(t); - if(!accepted) continue; -#endif - - if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays - if(time <= mintime) mintime=time; - tracks[ngoodtrk]=t; - ngoodtrk++; - stripmean += strip; - } - if(ngoodtrk) stripmean /= ngoodtrk; - - cout << " T0 offset selected for this sample : " << fT0Offset << endl; - cout << " N. of ESD tracks : " << ntrk << endl; - cout << " N. of preselected tracks : " << ngoodtrk << endl; - cout << " Minimum tof time in set (in ns) : " << mintime << endl; - - AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk]; - - for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) { - AliESDtrack *t=tracks[jtrk]; - Double_t tot = t->GetTOFsignalToT(); - Double_t time=t->GetTOFsignalRaw(); - Int_t chan = t->GetTOFCalChannel(); - Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0); - time -= corr*1000.; - - Int_t create = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(create == 63 || create == 62){ - time += 9200; - } - - if((time-mintime*1.E3)<50.E3){ // For pp and per - gtracks[ngoodtrkt0]=t; - ngoodtrkt0++; - } - } + Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass); + + return sigma; +} + +//__________________________________________________________________ +Bool_t AliTOFT0v1::AcceptTrack(AliESDtrack *track) +{ + + /* TPC refit */ + if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE; + /* do not accept kink daughters */ + if (track->GetKinkIndex(0)>0) return kFALSE; + /* N clusters TPC */ + if (track->GetTPCclusters(0) < 50) return kFALSE; + /* chi2 TPC */ + if (track->GetTPCchi2()/Float_t(track->GetTPCclusters(0)) > 3.5) return kFALSE; + /* sigma to vertex */ + if (GetSigmaToVertex(track) > 4.) return kFALSE; - Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1; - Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq; - if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++; + /* accept track */ + return kTRUE; - if(ngoodtrkt0 > 1){ - Int_t nlastset = (ngoodtrkt0 % nmaxtracksinsetCurrent); +} - while(nlastset-nseteq+1 > 2 ){ - nmaxtracksinsetCurrent++; - nlastset -= nseteq-1; - } - if(nmaxtracksinsetCurrent > nmaxtracksinset) nmaxtracksinsetCurrent = nmaxtracksinset; - } +//____________________________________________________________________ +Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const +{ + // Calculates the number of sigma to the vertex. - if(ngoodtrkt0<2){ - cout << "less than 2 tracks, skip event " << endl; - t0def=-999; - deltat0def=0.600; - fT0SigmaT0def[0]=t0def; - fT0SigmaT0def[1]=deltat0def; - fT0SigmaT0def[2]=ngoodtrkt0; - fT0SigmaT0def[3]=ngoodtrkt0; - //goto finish; + Float_t b[2]; + Float_t bRes[2]; + Float_t bCov[3]; + esdTrack->GetImpactParameters(b,bCov); + + if (bCov[0]<=0 || bCov[2]<=0) { + bCov[0]=0; bCov[2]=0; } - if(ngoodtrkt0>=2){ - // Decide how many tracks in set - Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent); - Int_t nset=1; - if(ngoodtrkt0>nmaxtracksinset) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;} - - cout << " number of sets = " << nset << endl; + bRes[0] = TMath::Sqrt(bCov[0]); + bRes[1] = TMath::Sqrt(bCov[2]); + + // ----------------------------------- + // How to get to a n-sigma cut? + // + // The accumulated statistics from 0 to d is + // + // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma) + // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma) + // + // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2) + // Can this be expressed in a different way? + + if (bRes[0] == 0 || bRes[1] ==0) + return -1; + + //Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); + Float_t d = TMath::Sqrt(ToCalculatePower(b[0]/bRes[0],2) + ToCalculatePower(b[1]/bRes[1],2)); + + // work around precision problem + // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :( + // 1e-15 corresponds to nsigma ~ 7.7 + if (TMath::Exp(-d * d / 2) < 1e-15) + return 1000; + + Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); + return nSigma; +} +//____________________________________________________________________ + +Bool_t AliTOFT0v1::CheckTPCMatching(AliESDtrack *track,Int_t imass) const{ + Bool_t status = kFALSE; - if(strstr(option,"tim") || strstr(option,"all"))gBenchmark->Start("TOFT0v1"); + Double_t exptimes[AliPID::kSPECIESC]; + track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC); - // Loop over selected sets + Float_t dedx = track->GetTPCsignal(); - if(nset>=1){ - for (Int_t i=0; i< nset; i++) { - - Float_t t0best=999.; - Float_t Et0best=999.; - Float_t chisquarebest=99999.; - Int_t npionbest=0; - - Int_t ntracksinsetmy=0; - AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset]; - for (Int_t itrk=0; itrkGetP(); - Double_t mom=momOld-0.0036*momOld; - Double_t tot = t->GetTOFsignalToT(); - Double_t time=t->GetTOFsignalRaw(); - Int_t chan = t->GetTOFCalChannel(); - Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0); - time -= corr*1000.; - - Int_t create = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(create == 63 || create == 62){ - time += 9200; - } + Double_t ptpc[3]; + track->GetInnerPxPyPz(ptpc); + Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]); - time*=1.E-3; // tof given in nanoseconds - Double_t exptime[10]; t->GetIntegratedTimes(exptime); - Double_t toflen=t->GetIntegratedLength(); - toflen=toflen/100.; // toflen given in m - - timeofflight[j]=time+fT0Offset; - tracktoflen[j]=toflen+fLOffset; - exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns - exptof[j][1]=exptime[3]*1.E-3+fTimeCorr; - exptof[j][2]=exptime[4]*1.E-3+fTimeCorr; - momentum[j]=mom; - assparticle[j]=3; - - } //end for (Int_t j=0; j %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]); - - Bool_t kRedoT0 = kFALSE; - ntracksinsetmyCut = ntracksinsetmy; - Bool_t usetrack[nmaxtracksinset]; - for (Int_t icsq=0; icsq chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){ - kRedoT0 = kTRUE; - ntracksinsetmyCut--; - usetrack[icsq] = kFALSE; - } - } // end loop for (Int_t icsq=0; icsq<15;icsq++) - - printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut); - - // Loop on mass hypotheses Redo - if(kRedoT0 && ntracksinsetmyCut > 1){ - printf("Redo T0\n"); - for (Int_t k=0; k < ncombinatorial;k++) { - for (Int_t j=0; j= 999){ - printf("How is it possible (chi2 = %f)? T0best = %f\n",chisquarebest,t0best); - - for(Int_t icsq=0; icsq0.){ - meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns] - Emeantzeropi=sqrt(1./sumAllweightspi); // it is given in [ns] - } - - if(sumWt0bestallSel>0){ - t0bestallSel = t0bestallSel/sumWt0bestallSel; - Et0bestallSel = sqrt(1./sumWt0bestallSel); - }// end of if(sumWt0bestallSel>0){ + if(imass > 2 || imass < 0) return status; + Int_t i = imass+2; + + AliPID::EParticleType type=AliPID::EParticleType(i); + + Float_t dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,type); + Float_t resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),type); - cout << "T0 all " << t0bestallSel << " +/- " << Et0bestallSel << "Number of tracks used: "<Stop("TOFT0v1"); - cout << "AliTOFT0v1:" << endl ; - cout << " took " << gBenchmark->GetCpuTime("TOFT0v1") << " seconds in order to calculate T0 " << endl ; - } - printf("T0 from TOF = %f ns\n",fT0SigmaT0def[0]); - - return fT0SigmaT0def; - } -//__________________________________________________________________ -void AliTOFT0v1::SetTZeroFile(char * file ){ - cout << "Destination file : " << file << endl ; - fT0File=file; -} -//__________________________________________________________________ -void AliTOFT0v1::Print(Option_t* /*option*/)const -{ - cout << "------------------- "<< GetName() << " -------------" << endl ; - if(!fT0File.IsNull()) - cout << " Writing T0 Distribution to file " << (char*) fT0File.Data() << endl ; + + return status; } -//__________________________________________________________________ -Bool_t AliTOFT0v1::operator==( AliTOFT0v1 const &tzero )const +//____________________________________________________________________ +Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const { - // Equal operator. - // - if( (fTimeResolution==tzero.fTimeResolution)&&(fLowerMomBound==tzero.fLowerMomBound)&&(fUpperMomBound==tzero.fUpperMomBound)) - return kTRUE ; - else - return kFALSE ; -} + // + // Returns base^exponent + // + Float_t power=1.; -//__________________________________________________________________ -Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) + for (Int_t ii=exponent; ii>0; ii--) + power=power*base; + + return power; + +} +//____________________________________________________________________ +Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const { - static const Double_t kMasses[]={ - 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 - }; + // + // Returns base^exponent + // - Double_t mass=kMasses[index+2]; - Double_t dpp=0.01; //mean relative pt resolution; - Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass)); + Int_t power=1; - sigma =TMath::Sqrt(sigma*sigma); + for (Int_t ii=exponent; ii>0; ii--) + power=power*base; + + return power; - return sigma; }