X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TOF%2FAliTOFT0v1.cxx;h=701e110839e50d3d3ed630f8ad13cd7de8f7c3b0;hb=b75653e3ede36cb78808d5b6c8553a1d7f9dea79;hp=11ccf42aa970f7763d15c0215024acb92efaf673;hpb=6a4e212e7b03d039a2cee463c3464e206b272728;p=u%2Fmrichter%2FAliRoot.git diff --git a/TOF/AliTOFT0v1.cxx b/TOF/AliTOFT0v1.cxx index 11ccf42aa97..701e110839e 100644 --- a/TOF/AliTOFT0v1.cxx +++ b/TOF/AliTOFT0v1.cxx @@ -53,79 +53,79 @@ // Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions //-- Author: F. Pierella //-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni +//-- AliTOFT0v1 contains code speed up provided by Jens Wiechula (look-up table +// for Power3) ////////////////////////////////////////////////////////////////////////////// -#include "Riostream.h" - -#include "AliTOFT0v1.h" #include "AliESDtrack.h" -#include "AliTOFcalibHisto.h" #include "AliESDEvent.h" +#include "AliTOFT0v1.h" +#include "TBenchmark.h" +#include "AliPID.h" +#include "AliESDpid.h" ClassImp(AliTOFT0v1) //____________________________________________________________________________ -AliTOFT0v1::AliTOFT0v1(): - fLowerMomBound(0.4), - fUpperMomBound(2.0), - fTimeResolution(0.80e-10), +AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID): + TObject(), + fLowerMomBound(0.5), + fUpperMomBound(3), fTimeCorr(0.), fEvent(0x0), - fCalib(0x0) + fPIDesd(extPID), + fTracks(new TObjArray(10)), + fGTracks(new TObjArray(10)), + fTracksT0(new TObjArray(10)), + fOptFlag(kFALSE) { // // default constructor // - - fT0SigmaT0def[0]=-999.; - fT0SigmaT0def[1]=999.; - fT0SigmaT0def[2]=-999.; - fT0SigmaT0def[3]=-999.; + if(AliPID::ParticleMass(0) == 0) new AliPID(); + + 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(AliESDEvent* event): - fLowerMomBound(0.4), - fUpperMomBound(2.0), - fTimeResolution(0.80e-10), +AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID): + TObject(), + fLowerMomBound(0.5), + fUpperMomBound(3.0), fTimeCorr(0.), fEvent(event), - fCalib(0x0) + fPIDesd(extPID), + fTracks(new TObjArray(10)), + fGTracks(new TObjArray(10)), + fTracksT0(new TObjArray(10)), + fOptFlag(kFALSE) { // // real constructor // - - fT0SigmaT0def[0]=-999.; - fT0SigmaT0def[1]= 999.; - fT0SigmaT0def[2]=-999.; - fT0SigmaT0def[3]=-999.; + if(AliPID::ParticleMass(0) == 0) new AliPID(); -} - -//____________________________________________________________________________ -AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero): - TObject(), - fLowerMomBound(tzero.fLowerMomBound), - fUpperMomBound(tzero.fUpperMomBound), - fTimeResolution(tzero.fTimeResolution), - fTimeCorr(tzero.fTimeCorr), - fEvent(tzero.fEvent), - fCalib(tzero.fCalib) -{ - // - // copy constructor - // - - fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0]; - fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1]; - fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2]; - fT0SigmaT0def[3]=tzero.fT0SigmaT0def[3]; + 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) { @@ -138,39 +138,89 @@ AliTOFT0v1& AliTOFT0v1::operator=(const AliTOFT0v1 &tzero) fLowerMomBound=tzero.fLowerMomBound; fUpperMomBound=tzero.fUpperMomBound; - fTimeResolution=tzero.fTimeResolution; fTimeCorr=tzero.fTimeCorr; fEvent=tzero.fEvent; - fCalib=tzero.fCalib; 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 - fCalib=NULL; 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){ - // Set the TOF time resolution - 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) { + TBenchmark *bench=new TBenchmark(); + bench->Start("t0computation"); + // Caluclate the Event Time using the ESD TOF time - Float_t timeresolutioninns=fTimeResolution*(1.e+9) * TMath::Sqrt(2.); // convert in [ns] + fT0SigmaT0def[0]=0.; + fT0SigmaT0def[1]=0.600; + fT0SigmaT0def[2]=0.; + fT0SigmaT0def[3]=0.; - const Int_t nmaxtracksinset=10; + 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; @@ -197,55 +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 (!AcceptTrack(t)) continue; - if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays - if(time <= mintime) mintime=time; - tracks[ngoodtrk]=t; + if(t->GetIntegratedLength() < 350)continue; //skip decays + if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue; + + meantime+=time; + fTracks->AddLast(t); ngoodtrk++; } - - -// 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; @@ -256,7 +307,7 @@ 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;} @@ -265,46 +316,47 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) 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 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 @@ -346,14 +399,21 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) 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"); for (Int_t k=0; k < ncombinatorial;k++) { for (Int_t j=0; j0){ for(Int_t iqsq = 0; iqsq1% - // if(confLevel>0.01 && ngoodsetsSel<200){ if(confLevel>0.01 && ngoodsetsSel<200){ chiSquarebestSel[ngoodsetsSel]=chisquarebest; confLevelbestSel[ngoodsetsSel]=confLevel; @@ -566,17 +609,42 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) 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); } } - delete[] tracksT0; + fTracksT0->Clear(); nsets++; } // end for the current set + //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.){ @@ -584,550 +652,43 @@ Double_t * AliTOFT0v1::DefineT0(Option_t *option) 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); - + // printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel); }// end of if(sumWt0bestallSel>0){ -// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<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 - - 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 = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(crate == 63 || crate == 62){ - time += 9200; - } - - Int_t strip = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kSectorStrip,chan); - - time*=1.E-3; // tof given in nanoseconds - if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue; - - if (!AcceptTrack(t)) continue; - - 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 << " 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 crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(crate == 63 || crate == 62){ - time += 9200; - } - - if((time-mintime*1.E3)<50.E3){ // For pp and per - gtracks[ngoodtrkt0]=t; - ngoodtrkt0++; - } - } - - Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1; - Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq; - if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++; - if(ngoodtrkt0 > 1){ - Int_t nlastset = (ngoodtrkt0 % nmaxtracksinsetCurrent); + fGTracks->Clear(); - while(nlastset-nseteq+1 > 2 ){ - nmaxtracksinsetCurrent++; - nlastset -= nseteq-1; - } - if(nmaxtracksinsetCurrent > nmaxtracksinset) nmaxtracksinsetCurrent = nmaxtracksinset; - } + if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6; - 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; - } - 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;} - - // Loop over selected sets - - 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 crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan); - if(crate == 63 || crate == 62){ - time += 9200; - } - - 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; - 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]); - - 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; jStop("t0computation"); - for (Int_t itz=0; itzGetRealTime("t0computation"); + fT0SigmaT0def[5]=bench->GetCpuTime("t0computation"); - // calculate chisquare - - Float_t chisquare=0.; - for (Int_t icsq=0; icsq= 999){ - printf("How is it possible (chi2 = %f)? T0best = %f\n",chisquarebest,t0best); - -// for(Int_t icsq=0; icsqPrint("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])); - 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){ - - } - - t0def=t0bestallSel; - deltat0def=eT0bestallSel; - if ((TMath::Abs(t0bestallSel)<0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){ - t0def=-999; deltat0def=0.600; - } - - fT0SigmaT0def[0]=t0def; - fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1))); - fT0SigmaT0def[2]=ngoodtrkt0; - fT0SigmaT0def[3]=ngoodtrktrulyused; - } - } - return fT0SigmaT0def; } - //__________________________________________________________________ Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const { @@ -1138,10 +699,8 @@ Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const }; 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)); - sigma =TMath::Sqrt(sigma*sigma); + Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,texp,mass); return sigma; } @@ -1196,7 +755,8 @@ Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const 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(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 :( @@ -1207,3 +767,62 @@ Float_t AliTOFT0v1::GetSigmaToVertex(AliESDtrack* esdTrack) const 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; + + Double_t exptimes[AliPID::kSPECIESC]; + track->GetIntegratedTimes(exptimes,AliPID::kSPECIESC); + + Float_t dedx = track->GetTPCsignal(); + + Double_t ptpc[3]; + track->GetInnerPxPyPz(ptpc); + Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]); + + 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); + + if(TMath::Abs(dedx - dedxExp) < 5 * resolutionTPC){ + status = kTRUE; + } + + return status; +} + +//____________________________________________________________________ +Float_t AliTOFT0v1::ToCalculatePower(Float_t base, Int_t exponent) const +{ + // + // Returns base^exponent + // + + Float_t power=1.; + + for (Int_t ii=exponent; ii>0; ii--) + power=power*base; + + return power; + +} +//____________________________________________________________________ +Int_t AliTOFT0v1::ToCalculatePower(Int_t base, Int_t exponent) const +{ + // + // Returns base^exponent + // + + Int_t power=1; + + for (Int_t ii=exponent; ii>0; ii--) + power=power*base; + + return power; + +}