// 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 "AliESDtrack.h"
fPIDesd(extPID),
fTracks(new TObjArray(10)),
fGTracks(new TObjArray(10)),
- fTracksT0(new TObjArray(10))
+ fTracksT0(new TObjArray(10)),
+ fOptFlag(kFALSE)
{
//
// default constructor
}
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);
+ }
}
//____________________________________________________________________________
fPIDesd(extPID),
fTracks(new TObjArray(10)),
fGTracks(new TObjArray(10)),
- fTracksT0(new TObjArray(10))
+ fTracksT0(new TObjArray(10)),
+ fOptFlag(kFALSE)
{
//
// real constructor
}
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)
for (Int_t ii=0; ii<tzero.fTracksT0->GetEntries(); ii++)
fTracksT0->AddLast(tzero.fTracksT0->At(ii));
+
+ fOptFlag=tzero.fOptFlag;
return *this;
}
{
// dtor
fEvent=NULL;
+
+ if (fTracks) {
+ fTracks->Clear();
+ delete fTracks;
+ fTracks=0x0;
+ }
- fTracks->Delete();
- delete fTracks;
- fTracks=0x0;
+ if (fGTracks) {
+ fGTracks->Clear();
+ delete fGTracks;
+ fGTracks=0x0;
+ }
- fGTracks->Delete();
- delete fGTracks;
- fGTracks=0x0;
+ if (fTracksT0) {
+ fTracksT0->Clear();
+ delete fTracksT0;
+ fTracksT0=0x0;
+ }
- fTracksT0->Delete();
- delete fTracksT0;
- fTracksT0=0x0;
}
//____________________________________________________________________________
fT0SigmaT0def[3]=0.;
}
-
-//____________________________________________________________________________
-Double_t * AliTOFT0v1::DefineT0(Option_t *option)
-{
- // Caluclate the Event Time using the ESD TOF time
-
- TBenchmark *bench=new TBenchmark();
- bench->Start("t0computation");
-
- fT0SigmaT0def[0]=0.;
- fT0SigmaT0def[1]=0.600;
- fT0SigmaT0def[2]=0.;
- fT0SigmaT0def[3]=0.;
-
- 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 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;
-
- Int_t ntrk=fEvent->GetNumberOfTracks();
-
- //AliESDtrack **fTracks=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; itrk<ntrk; itrk++) {
- AliESDtrack *t=fEvent->GetTrack(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;
- //fTracks[ngoodtrk]=t;
- fTracks->AddLast(t);
- ngoodtrk++;
- }
-
- if(ngoodtrk>22) nmaxtracksinset = 6;
-
-// 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 **fGTracks=new AliESDtrack*[ngoodtrk];
-
- //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
- for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
- //AliESDtrack *t=fTracks[jtrk];
- AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
- Double_t time=t->GetTOFsignal();
-
- if((time-mintime*1.E3)<50.E3){ // For pp and per
- //fGTracks[ngoodtrkt0]=t;
- fGTracks->AddLast(t);
- ngoodtrkt0++;
- }
- }
-
- //delete[] fTracks;
- 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;
- fT0SigmaT0def[1]=deltat0def;
- fT0SigmaT0def[2]=ngoodtrkt0;
- fT0SigmaT0def[3]=ngoodtrkt0;
- //goto finish;
- }
- if(ngoodtrkt0>=2){
- // Decide how many tracks in set
- Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
- Int_t nset=1;
-
- if(ngoodtrkt0>nmaxtracksinsetCurrent) {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 **fTracksT0=new AliESDtrack*[ntracksinset];
- for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
- Int_t index = itrk+i*ntracksinset;
- //if(index < ngoodtrkt0){
- if(index < fGTracks->GetEntries()){
- //AliESDtrack *t=fGTracks[index];
- AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
- //fTracksT0[itrk]=t;
- fTracksT0->AddLast(t);
- ntracksinsetmy++;
- }
- }
-
- // Analyse it
-
- Int_t assparticle[nmaxtracksinsetMax];
- Float_t exptof[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[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; j<ntracksinset; j++) {
- assparticle[j] = 3;
- timeofflight[j] = 0;
- momentum[j] = 0;
- timezero[j] = 0;
- weightedtimezero[j] = 0;
- beta[j] = 0;
- texp[j] = 0;
- dtexp[j] = 0;
- sqMomError[j] = 0;
- sqTrackError[j] = 0;
- tracktoflen[j] = 0;
- besttimezero[j] = 0;
- besttexp[j] = 0;
- besttimeofflight[j] = 0;
- bestmomentum[j] = 0;
- bestchisquare[j] = 0;
- bestweightedtimezero[j] = 0;
- bestsqTrackError[j] = 0;
- imass[j] = 1;
- }
-
- //for (Int_t j=0; j<ntracksinsetmy; j++) {
- for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
- //AliESDtrack *t=fTracksT0[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 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<ntracksinsetmy; j++) {
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
- +momentum[itz]*momentum[itz]);
- 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]));
- sqTrackError[itz]=sqMomError[itz]; //in ns
- timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllweightspi+=1./sqTrackError[itz];
- meantzeropi+=weightedtimezero[itz];
- } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
-
-
- // Then, Combinatorial Algorithm
-
- if(ntracksinsetmy<2 )break;
-
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- imass[j] = 3;
- }
-
- //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
- Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
-
- // Loop on mass hypotheses
- for (Int_t k=0; k < ncombinatorial;k++) {
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
- texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
- }
-
- Float_t sumAllweights=0.;
- Float_t meantzero=0.;
- Float_t eMeanTzero=0.;
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
-
- timezero[itz]=texp[itz]-timeofflight[itz];// in ns
-
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllweights+=1./sqTrackError[itz];
- meantzero+=weightedtimezero[itz];
-
- } // end loop for (Int_t itz=0; itz<15;itz++)
-
- meantzero=meantzero/sumAllweights; // it is given in [ns]
- eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
-
- // calculate chisquare
-
- Float_t chisquare=0.;
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
- // else chisquare+=1000;
- } // end loop for (Int_t icsq=0; icsq<15;icsq++)
-
- if(chisquare<=chisquarebest){
- for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
-
- bestsqTrackError[iqsq]=sqTrackError[iqsq];
- besttimezero[iqsq]=timezero[iqsq];
- bestmomentum[iqsq]=momentum[iqsq];
- besttimeofflight[iqsq]=timeofflight[iqsq];
- besttexp[iqsq]=texp[iqsq];
- bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
- bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
- // else bestchisquare[iqsq]=1000;
- }
-
- Int_t npion=0;
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- assparticle[j]=imass[j];
- if(imass[j] == 0) npion++;
- }
- npionbest=npion;
- chisquarebest=chisquare;
- t0best=meantzero;
- eT0best=eMeanTzero;
- } // close if(dummychisquare<=chisquare)
-
- }
-
- Double_t chi2cut[nmaxtracksinsetMax];
- chi2cut[0] = 0;
- chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
- for (Int_t j=2; j<ntracksinset; j++) {
- chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
- }
-
- Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
-
-// 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[nmaxtracksinsetMax];
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- usetrack[icsq] = kTRUE;
- if((bestchisquare[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<ntracksinsetmy; j++) {
- //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
- texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
- }
-
- Float_t sumAllweights=0.;
- Float_t meantzero=0.;
- Float_t eMeanTzero=0.;
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- if(! usetrack[itz]) continue;
- sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
-
- timezero[itz]=texp[itz]-timeofflight[itz];// in ns
-
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllweights+=1./sqTrackError[itz];
- meantzero+=weightedtimezero[itz];
-
- } // end loop for (Int_t itz=0; itz<15;itz++)
-
- meantzero=meantzero/sumAllweights; // it is given in [ns]
- eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
-
- // calculate chisquare
-
- Float_t chisquare=0.;
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- if(! usetrack[icsq]) continue;
- chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
- // else chisquare+=1000;
- } // end loop for (Int_t icsq=0; icsq<15;icsq++)
-
- Int_t npion=0;
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- assparticle[j]=imass[j];
- if(imass[j] == 0) npion++;
- }
-
- if(chisquare<=chisquarebest && npion>0){
- for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
- if(! usetrack[iqsq]) continue;
- bestsqTrackError[iqsq]=sqTrackError[iqsq];
- besttimezero[iqsq]=timezero[iqsq];
- bestmomentum[iqsq]=momentum[iqsq];
- besttimeofflight[iqsq]=timeofflight[iqsq];
- besttexp[iqsq]=texp[iqsq];
- bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
- bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
- // else bestchisquare[iqsq]=1000;
- }
-
- npionbest=npion;
- chisquarebest=chisquare;
- t0best=meantzero;
- eT0best=eMeanTzero;
- } // close if(dummychisquare<=chisquare)
-
- }
- }
-
- // filling histos
- Float_t confLevel=999;
-
- // Sets with decent chisquares
-
- if(chisquarebest<999.){
- Double_t dblechisquare=(Double_t)chisquarebest;
- confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
-// cout << " Set Number " << nsets << endl;
-// cout << "Best Assignment, selection " << assparticle[0] <<
-// assparticle[1] << assparticle[2] <<
-// assparticle[3] << assparticle[4] <<
-// assparticle[5] << endl;
-// cout << " Chisquare of the set "<< chisquarebest <<endl;
-// cout << " C.L. of the set "<< confLevel <<endl;
-// cout << " T0 for this set (in ns) " << t0best << endl;
-
- for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
-
- if(! usetrack[icsq]) continue;
-
-// cout << "Track # " << icsq << " T0 offsets = "
-// << besttimezero[icsq]-t0best <<
-// " track error = " << bestsqTrackError[icsq]
-// << " Chisquare = " << bestchisquare[icsq]
-// << " Momentum = " << bestmomentum[icsq]
-// << " TOF = " << besttimeofflight[icsq]
-// << " TOF tracking = " << besttexp[icsq]
-// << " is used = " << usetrack[icsq] << endl;
- }
-
- // Pick up only those with C.L. >1%
- // 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;
- ngoodsetsSel++;
- ngoodtrktrulyused+=ntracksinsetmyCut;
-
- // printf("t0 of a set = %f +/- %f\n",t0best,eT0best);
- }
- else{
- // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
- }
- }
- //delete[] fTracksT0;
- 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; itz<nsamples;itz++) {
- if(TMath::Abs(t0bestSel[itz]/eT0bestSel[itz]-t0BestStep1)<1.0){
- t0bestallSel += t0bestSel[itz];
- sumWt0bestallSel += eT0bestSel[itz];
- ngoodsetsSel++;
-// printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
- }
- else{
-// printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
- }
- }
- }
- if(ngoodsetsSel < 1){
- sumWt0bestallSel = 0.0;
- }
- //--------------------------------End recomputation
-
- nUsedTracks = ngoodtrkt0;
- if(strstr(option,"all")){
- if(sumAllweightspi>0.){
- 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){
-
-// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
- }
-
- t0def=t0bestallSel;
- deltat0def=eT0bestallSel;
-
- fT0SigmaT0def[0]=t0def;
- fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
- fT0SigmaT0def[2]=ngoodtrkt0;
- fT0SigmaT0def[3]=ngoodtrktrulyused;
- }
- }
-
- //delete [] fGTracks;
- fGTracks->Clear();
-
- // if(strstr(option,"tim") || strstr(option,"all")){
- // cout << "AliTOFT0v1:" << endl ;
- //}
-
- 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("%f %f\n",fT0SigmaT0def[4],fT0SigmaT0def[5]);
-
-
- delete bench;
- bench=NULL;
-
- return fT0SigmaT0def;
- }
-//__________________________________________________________________
+//____________________________________________________________________________
Double_t * AliTOFT0v1::DefineT0(Option_t *option,Float_t pMinCut,Float_t pMaxCut)
{
TBenchmark *bench=new TBenchmark();
Int_t ntrk=fEvent->GetNumberOfTracks();
- //AliESDtrack **fTracks=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; itrk<ntrk; itrk++) {
AliESDtrack *t=fEvent->GetTrack(itrk);
Double_t momOld=t->GetP();
if(t->GetIntegratedLength() < 350)continue; //skip decays
if(t->GetP() > pMinCut && t->GetP() < pMaxCut) continue;
- if(time <= mintime) mintime=time;
- //fTracks[ngoodtrk]=t;
+
+ meantime+=time;
fTracks->AddLast(t);
ngoodtrk++;
}
-
+
+ if(ngoodtrk > 1) meantime /= ngoodtrk;
+
if(ngoodtrk>22) nmaxtracksinset = 6;
-
-// 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 **fGTracks=new AliESDtrack*[ngoodtrk];
-
- //for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
+ fGTracks->Clear();
for (Int_t jtrk=0; jtrk< fTracks->GetEntries(); jtrk++) {
- //AliESDtrack *t=fTracks[jtrk];
AliESDtrack *t=(AliESDtrack*)fTracks->At(jtrk);
- Double_t time=t->GetTOFsignal();
-
- if((time-mintime*1.E3)<50.E3){ // For pp and per
- //fGTracks[ngoodtrkt0]=t;
- fGTracks->AddLast(t);
- ngoodtrkt0++;
- }
+ // Double_t time=t->GetTOFsignal();
+ // if((time-meantime*1.E3)<50.E3){ // For pp and per
+ fGTracks->AddLast(t);
+ ngoodtrkt0++;
+ // }
}
- //delete [] fTracks;
fTracks->Clear();
-// cout << " N. of preselected tracks t0 : " << ngoodtrkt0 << endl;
-
Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
-// cout << " N. of sets : " << nseteq << endl;
if(ngoodtrkt0<2){
-// cout << "less than 2 tracks, skip event " << endl;
t0def=-999;
deltat0def=0.600;
fT0SigmaT0def[0]=t0def;
Float_t chisquarebest=99999.;
Int_t npionbest=0;
+ fTracksT0->Clear();
Int_t ntracksinsetmy=0;
- //AliESDtrack **fTracksT0=new AliESDtrack*[ntracksinset];
for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
Int_t index = itrk+i*ntracksinset;
- //if(index < ngoodtrkt0){
if(index < fGTracks->GetEntries()){
- //AliESDtrack *t=fGTracks[index];
AliESDtrack *t=(AliESDtrack*)fGTracks->At(index);
- //fTracksT0[itrk]=t;
fTracksT0->AddLast(t);
ntracksinsetmy++;
}
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];
imass[j] = 1;
}
- //for (Int_t j=0; j<ntracksinsetmy; j++) {
for (Int_t j=0; j<fTracksT0->GetEntries(); j++) {
- //AliESDtrack *t=fTracksT0[j];
AliESDtrack *t=(AliESDtrack*)fTracksT0->At(j);
Double_t momOld=t->GetP();
Double_t mom=momOld-0.0036*momOld;
exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
momentum[j]=mom;
assparticle[j]=3;
-
- } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
+
+ // in principle GetMomError only depends on two indices k=imass[j] and j itslef (see blow in the ncombinatorial loop)
+ // so it should be possible to make a lookup in order to speed up the code:
+ if (fOptFlag) {
+ momErr[j][0]=GetMomError(0, momentum[j], exptof[j][0]);
+ momErr[j][1]=GetMomError(1, momentum[j], exptof[j][1]);
+ momErr[j][2]=GetMomError(2, momentum[j], exptof[j][2]);
+ }
+ } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
for (Int_t j=0; j<ntracksinsetmy; j++) {
imass[j] = 3;
}
-
- //Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
- Int_t ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
-
+
+ Int_t ncombinatorial;
+ if (fOptFlag) ncombinatorial = fLookupPowerThree[ntracksinsetmy];
+ else ncombinatorial = ToCalculatePower(3,ntracksinsetmy);
+
+
// Loop on mass hypotheses
for (Int_t k=0; k < ncombinatorial;k++) {
for (Int_t j=0; j<ntracksinsetmy; j++) {
- //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j))/ToCalculatePower(3,ntracksinsetmy-j-1);
+ imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j])/fLookupPowerThree[ntracksinsetmy-j-1];
texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
+ if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
+ else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
}
Float_t sumAllweights=0.;
Float_t chisquare=0.;
for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
- // else chisquare+=1000;
} // end loop for (Int_t icsq=0; icsq<15;icsq++)
if(chisquare<=chisquarebest){
besttexp[iqsq]=texp[iqsq];
bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
- // else bestchisquare[iqsq]=1000;
}
Int_t npion=0;
Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
-// printf("tracks removed with a chi2 > %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;
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; j<ntracksinsetmy; j++) {
- //imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- imass[j] = (k % ToCalculatePower(3,ntracksinsetmy-j)) / ToCalculatePower(3,ntracksinsetmy-j-1);
+ imass[j] = (k % fLookupPowerThree[ntracksinsetmy-j]) / fLookupPowerThree[ntracksinsetmy-j-1];
texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- // if(! CheckTPCMatching((AliESDtrack*)fTracksT0->At(j),imass[j])) dtexp[j]*=100;
+ if (fOptFlag) dtexp[j]=momErr[j][imass[j]]; // see comments above in the initialisation of momErr
+ else dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
}
Float_t sumAllweights=0.;
for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
if(! usetrack[icsq]) continue;
chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(icsq),imass[icsq])) chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq]; // require TPC agreement
- // else chisquare+=1000;
} // end loop for (Int_t icsq=0; icsq<15;icsq++)
Int_t npion=0;
besttexp[iqsq]=texp[iqsq];
bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- // if(CheckTPCMatching((AliESDtrack*)fTracksT0->At(iqsq),imass[iqsq])) bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; // require TPC agreement
- // else bestchisquare[iqsq]=1000;
}
npionbest=npion;
if(chisquarebest<999.){
Double_t dblechisquare=(Double_t)chisquarebest;
confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
-// cout << " Set Number " << nsets << endl;
-// cout << "Best Assignment, selection " << assparticle[0] <<
-// assparticle[1] << assparticle[2] <<
-// assparticle[3] << assparticle[4] <<
-// assparticle[5] << endl;
-// cout << " Chisquare of the set "<< chisquarebest <<endl;
-// cout << " C.L. of the set "<< confLevel <<endl;
-// cout << " T0 for this set (in ns) " << t0best << endl;
Int_t ntrackincurrentsel=0;
for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
if(! usetrack[icsq]) continue;
-// cout << "Track # " << icsq << " T0 offsets = "
-// << besttimezero[icsq]-t0best <<
-// " track error = " << bestsqTrackError[icsq]
-// << " Chisquare = " << bestchisquare[icsq]
-// << " Momentum = " << bestmomentum[icsq]
-// << " TOF = " << besttimeofflight[icsq]
-// << " TOF tracking = " << besttexp[icsq]
-// << " is used = " << usetrack[icsq] << endl;
ntrackincurrentsel++;
}
// printf("%i) CL(Chi2) = %f < 0.01\n",ngoodsetsSel,confLevel);
// Pick up only those with C.L. >1%
- // if(confLevel>0.01 && ngoodsetsSel<200){
if(confLevel>0.01 && ngoodsetsSel<200){
chiSquarebestSel[ngoodsetsSel]=chisquarebest;
confLevelbestSel[ngoodsetsSel]=confLevel;
sumWt0bestallSel += 1./eT0best/eT0best;
ngoodsetsSel++;
ngoodtrktrulyused+=ntracksinsetmyCut;
-// printf("T0best = %f +/- %f (%i-%i) -- conflevel = %f\n",t0best,eT0best,ntrackincurrentsel,npionbest,confLevel);
+ // 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[] fTracksT0;
fTracksT0->Clear();
nsets++;
t0bestallSel += t0bestSel[itz];
sumWt0bestallSel += eT0bestSel[itz];
ngoodsetsSel++;
-// printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
+ // printf("not rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
}
else{
-// printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
+ // printf("rejected %f +/- %f\n",t0bestSel[itz]/eT0bestSel[itz],1./TMath::Sqrt(eT0bestSel[itz]));
}
}
}
// printf("Final) t0bestallSel = %f -- eT0bestallSel = %f\n",t0bestallSel,eT0bestallSel);
}// end of if(sumWt0bestallSel>0){
-// cout << "T0 all " << t0bestallSel << " +/- " << eT0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
}
t0def=t0bestallSel;
}
}
- //delete [] fGTracks;
fGTracks->Clear();
- // if(strstr(option,"tim") || strstr(option,"all")){
- // cout << "AliTOFT0v1:" << endl ;
- //}
-
if(fT0SigmaT0def[1] < 0.00001) fT0SigmaT0def[1] = 0.6;
bench->Stop("t0computation");
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff