fExternalPIDFlag(kFALSE),
fTOFcalib(NULL),
fNoTOFT0(0),
+ fNmomBins(0),
fTimeResolution(100),
fT0sigma(1000),
fHmapChannel(0),
fKmask(0),
- fT0width(150.)
+ fT0width(150.),
+ fT0spreadExt(-1.),
+ fT0fillExt(0)
{
// ctr
fCalculated[0] = 0;
fCalculated[2] = 0;
fCalculated[3] = 0;
- fT0TOF = new AliTOFT0v1();
if(AliPID::ParticleMass(0) == 0) new AliPID();
fPIDesd = new AliESDpid();
- fPtCutMin[0] = 0.3;
- fPtCutMin[1] = 0.5;
- fPtCutMin[2] = 0.6;
- fPtCutMin[3] = 0.7;
- fPtCutMin[4] = 0.8;
- fPtCutMin[5] = 0.9;
- fPtCutMin[6] = 1;
- fPtCutMin[7] = 1.2;
- fPtCutMin[8] = 1.5;
- fPtCutMin[9] = 2;
-
- fPtCutMax[0] = 0.5;
- fPtCutMax[1] = 0.6;
- fPtCutMax[2] = 0.7;
- fPtCutMax[3] = 0.8;
- fPtCutMax[4] = 0.9;
- fPtCutMax[5] = 1;
- fPtCutMax[6] = 1.2;
- fPtCutMax[7] = 1.5;
- fPtCutMax[8] = 2;
- fPtCutMax[9] = 3;
-
- /* init arrays */
- for (Int_t i = 0; i < 10; i++) {
- fT0pt[i] = 0.;
- fT0ptSigma[i] = 0.;
- }
+ fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
+ SetTOFResponse();
+
+ fT0TOF = new AliTOFT0v1(fPIDesd);
+
}
//____________________________________________________________________________
AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
fExternalPIDFlag(kTRUE),
fTOFcalib(tofCalib),
fNoTOFT0(0),
+ fNmomBins(0),
fTimeResolution(100),
fT0sigma(1000),
fHmapChannel(0),
fKmask(0),
- fT0width(150.)
+ fT0width(150.),
+ fT0spreadExt(-1.),
+ fT0fillExt(0)
{
// ctr
fCalculated[0] = 0;
fCalculated[2] = 0;
fCalculated[3] = 0;
- fT0TOF = new AliTOFT0v1();
if(AliPID::ParticleMass(0) == 0) new AliPID();
if(!fPIDesd){
fExternalPIDFlag = kFALSE;
}
- fPtCutMin[0] = 0.3;
- fPtCutMin[1] = 0.5;
- fPtCutMin[2] = 0.6;
- fPtCutMin[3] = 0.7;
- fPtCutMin[4] = 0.8;
- fPtCutMin[5] = 0.9;
- fPtCutMin[6] = 1;
- fPtCutMin[7] = 1.2;
- fPtCutMin[8] = 1.5;
- fPtCutMin[9] = 2;
-
- fPtCutMax[0] = 0.5;
- fPtCutMax[1] = 0.6;
- fPtCutMax[2] = 0.7;
- fPtCutMax[3] = 0.8;
- fPtCutMax[4] = 0.9;
- fPtCutMax[5] = 1;
- fPtCutMax[6] = 1.2;
- fPtCutMax[7] = 1.5;
- fPtCutMax[8] = 2;
- fPtCutMax[9] = 3;
-
- /* init arrays */
- for (Int_t i = 0; i < 10; i++) {
- fT0pt[i] = 0.;
- fT0ptSigma[i] = 0.;
- }
-}
-
-/* copy-constructor and operator= suppressed
-
-
-//____________________________________________________________________________
-AliTOFT0maker::AliTOFT0maker(const AliTOFT0maker & t) :
- TObject(t),
- fT0TOF(t.fT0TOF),
- fPIDESD(t.fPIDESD),
- fNoTOFT0(t.fNoTOFT0),
- fTimeResolution(t.fTimeResolution),
- fT0sigma(t.fT0sigma),
- fHmapChannel(t.fHmapChannel),
- fKmask(t.fKmask)
- {
- // copy ctr
-}
-
-//____________________________________________________________________________
-AliTOFT0maker& AliTOFT0maker::operator=(const AliTOFT0maker &t)
-{
- //
- // assign. operator
- //
-
- if (this == &t)
- return *this;
+ fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
+ SetTOFResponse();
- TObject::operator=(t);
- fTimeResolution = t.fTimeResolution;
- fT0sigma = t.fT0sigma;
+ fT0TOF = new AliTOFT0v1(fPIDesd);
- return *this;
}
-*/
-
//____________________________________________________________________________
AliTOFT0maker::~AliTOFT0maker()
{
// Remake TOF PID probabilities
//
- Double_t t0tof[4];
+ Double_t t0tof[6];
if(fKmask) ApplyMask(esd);
+ Double_t t0fill = 0.;
+
+ fPIDesd->GetTOFResponse().ResetT0info();
+
/* get T0 spread from TOFcalib if available otherwise use default value */
if (fTOFcalib && esd) {
AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
if (runParams && runParams->GetTimestamp(0) != 0) {
Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
- SetT0FillWidth(t0spread);
+ if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
+ else{
+ SetT0FillWidth(t0spread);
+ t0fill = fT0fillExt;
+ }
+ }
+ else{
+ SetT0FillWidth(fT0spreadExt);
+ t0fill = fT0fillExt;
}
}
fT0TOF->Init(esd);
AliTOFT0v1* t0maker= fT0TOF;
- t0maker->SetTimeResolution(fTimeResolution*1e-12);
t0maker->DefineT0("all",1.5,3.0);
t0tof[0] = t0maker->GetResult(0);
t0tof[1] = t0maker->GetResult(1);
t0tof[2] = t0maker->GetResult(2);
t0tof[3] = t0maker->GetResult(3);
+ t0tof[4] = t0maker->GetResult(4);
+ t0tof[5] = t0maker->GetResult(5);
Float_t lT0Current=0.;
fT0sigma=1000;
// Int_t nrun = esd->GetRunNumber();
- Double_t t0fill = 0.;
t0time += t0fill;
fCalculated[6]=sigmaFill; // sigma t0 fill
fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
+ //statistics
+ fCalculated[8] = t0tof[4]; // real time in s
+ fCalculated[9] = t0tof[5]; // cpu time in s
+
if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500 && fCalculated[1] < fTimeResolution*1.2){
fT0sigma=fCalculated[1];
lT0Current=fCalculated[0];
}
// T0 pt bin
- for(Int_t i=0;i<10;i++){
- t0maker->DefineT0("all",fPtCutMin[i],fPtCutMax[i]);
- t0tof[0] = t0maker->GetResult(0);
- t0tof[1] = t0maker->GetResult(1);
- t0tof[2] = t0maker->GetResult(2);
- t0tof[3] = t0maker->GetResult(3);
- fT0pt[i] =-1000*t0tof[0]; // best t0
- fT0ptSigma[i] =1000*t0tof[1]; // sigma best t0
-
- if(fT0ptSigma[i] < sigmaFill && fT0ptSigma[i] < fTimeResolution * 1.2 && TMath::Abs(fT0pt[i] - t0fill) < 500){
- // Ok T0
+ if(fCalculated[7] < 100){
+ for(Int_t i=0;i<fNmomBins;i++){
+ t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
+ t0tof[0] = t0maker->GetResult(0);
+ t0tof[1] = t0maker->GetResult(1);
+ t0tof[2] = t0maker->GetResult(2);
+ t0tof[3] = t0maker->GetResult(3);
+
+
+ Float_t t0bin =-1000*t0tof[0]; // best t0
+ Float_t t0binRes =1000*t0tof[1]; // sigma best t0
+
+ if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < 500){
+ // Ok T0
+ if(t0sigma < 1000){
+ Double_t w1 = 1./t0sigma/t0sigma;
+ Double_t w2 = 1./t0binRes/t0binRes;
+
+ Double_t wtot = w1+w2;
+
+ t0bin = (w1*t0time + w2*t0bin) / wtot;
+ t0binRes = TMath::Sqrt(1./wtot);
+ }
+ }
+ else{
+ t0bin = t0fill;
+ t0binRes = sigmaFill;
+ if(t0sigma < 1000){
+ t0bin = t0time;
+ t0binRes= t0sigma;
+ }
+ }
+ fPIDesd->GetTOFResponse().SetT0bin(i,t0bin);
+ fPIDesd->GetTOFResponse().SetT0binRes(i,t0binRes);
}
- else{
- fT0pt[i] = t0fill;
- fT0ptSigma[i] = sigmaFill;
+ }
+ else{
+ for(Int_t i=0;i<fNmomBins;i++){
+ fPIDesd->GetTOFResponse().SetT0bin(i,lT0Current);
+ fPIDesd->GetTOFResponse().SetT0binRes(i,fT0sigma);
}
}
- //----
- SetTOFResponse();
return fCalculated;
}
//____________________________________________________________________________
Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
- Int_t i=0;
- while(p > fPtCutMin[i] && i < 10) i++;
- if(i > 0) i--;
+ Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
- fT0cur[0] = fT0pt[i];
- fT0cur[1] = fT0ptSigma[i];
+ fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
+ fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
+
return fT0cur;
}
//____________________________________________________________________________
void AliTOFT0maker::SetTOFResponse(){
- fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
+ fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
}
//____________________________________________________________________________
Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
- Double_t *sigmaT0 = GetT0p(mom);
- fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(sigmaT0[1]*sigmaT0[1] + fTimeResolution*fTimeResolution));
Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
- fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
return sigma;
}
#include "AliESDtrack.h"
#include "AliESDEvent.h"
#include "AliTOFT0v1.h"
+#include "TBenchmark.h"
+#include "AliPID.h"
+#include "AliESDpid.h"
ClassImp(AliTOFT0v1)
//____________________________________________________________________________
-AliTOFT0v1::AliTOFT0v1():
+AliTOFT0v1::AliTOFT0v1(AliESDpid *extPID):
TObject(),
fLowerMomBound(0.5),
fUpperMomBound(3),
- fTimeResolution(0.80e-10),
fTimeCorr(0.),
- fEvent(0x0)
-// fCalib(0x0)
+ fEvent(0x0),
+ fPIDesd(extPID)
{
//
// default constructor
//
+ if(AliPID::ParticleMass(0) == 0) new AliPID();
+
+ if(!fPIDesd){
+ fPIDesd = new AliESDpid();
+ }
Init(NULL);
}
-
//____________________________________________________________________________
-AliTOFT0v1::AliTOFT0v1(AliESDEvent* event):
+AliTOFT0v1::AliTOFT0v1(AliESDEvent* event,AliESDpid *extPID):
TObject(),
fLowerMomBound(0.5),
fUpperMomBound(3.0),
- fTimeResolution(0.80e-10),
fTimeCorr(0.),
- fEvent(event)
-// fCalib(0x0)
+ fEvent(event),
+ fPIDesd(extPID)
{
//
// real constructor
//
-
- Init(event);
-
-}
+ if(AliPID::ParticleMass(0) == 0) new AliPID();
-/* copy-constructor and operator= suppresed
+ if(!fPIDesd){
+ fPIDesd = new AliESDpid();
+ }
-//____________________________________________________________________________
-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];
+ Init(event);
}
-
//____________________________________________________________________________
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];
return *this;
}
-*/
//____________________________________________________________________________
AliTOFT0v1::~AliTOFT0v1()
{
// dtor
-// fCalib=NULL;
fEvent=NULL;
}
}
-//____________________________________________________________________________
-void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
- // Set the TOF time resolution
- fTimeResolution=timeresolution;
-}
-//____________________________________________________________________________
//____________________________________________________________________________
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.;
-
- Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
- 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;
time*=1.E-3; // tof given in nanoseconds
if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
-
+
if (!AcceptTrack(t)) continue;
- if(t->GetIntegratedLength() < 350)continue; //skip decays
+ if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
if(time <= mintime) mintime=time;
tracks[ngoodtrk]=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;
// 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 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; j<ntracksinset; j++) {
assparticle[j] = 3;
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]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
+ 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];
Float_t eMeanTzero=0.;
for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
+ sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
timezero[itz]=texp[itz]-timeofflight[itz];// in ns
}
- Double_t chi2cut[nmaxtracksinset];
+ 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++) {
Bool_t kRedoT0 = kFALSE;
ntracksinsetmyCut = ntracksinsetmy;
- Bool_t usetrack[nmaxtracksinset];
+ 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)){
for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
if(! usetrack[itz]) continue;
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
+ sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
timezero[itz]=texp[itz]-timeofflight[itz];// in ns
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);
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[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
fT0SigmaT0def[2]=ngoodtrkt0;
fT0SigmaT0def[3]=ngoodtrktrulyused;
}
// cout << "AliTOFT0v1:" << endl ;
//}
- if(fT0SigmaT0def[1] < 0.01) fT0SigmaT0def[1] = 0.6;
+ 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]);
+
+ bench->Reset();
return fT0SigmaT0def;
}
//__________________________________________________________________
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
fT0SigmaT0def[0]=0.;
fT0SigmaT0def[1]=0.600;
fT0SigmaT0def[2]=0.;
fT0SigmaT0def[3]=0.;
-
- Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
- 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;
time*=1.E-3; // tof given in nanoseconds
if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
-
+
if (!AcceptTrack(t)) continue;
if(t->GetIntegratedLength() < 350)continue; //skip decays
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;
+// 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];
}
}
+// 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;
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.;
// 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 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; j<ntracksinset; j++) {
assparticle[j] = 3;
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]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
+ 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];
Float_t eMeanTzero=0.;
for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
+ sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
timezero[itz]=texp[itz]-timeofflight[itz];// in ns
Float_t chisquare=0.;
for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
-
} // end loop for (Int_t icsq=0; icsq<15;icsq++)
if(chisquare<=chisquarebest){
t0best=meantzero;
eT0best=eMeanTzero;
} // close if(dummychisquare<=chisquare)
-
}
- Double_t chi2cut[nmaxtracksinset];
+ 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++) {
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;
- Bool_t usetrack[nmaxtracksinset];
+ 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)){
} // 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 itz=0; itz<ntracksinsetmy;itz++) {
if(! usetrack[itz]) continue;
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
+ sqTrackError[itz]=dtexp[itz]*dtexp[itz]*1E-6; //in ns2
timezero[itz]=texp[itz]-timeofflight[itz];// in ns
Float_t confLevel=999;
// Sets with decent chisquares
+ // printf("Chi2best of the set = %f \n",chisquarebest);
if(chisquarebest<999.){
Double_t dblechisquare=(Double_t)chisquarebest;
// << " is used = " << usetrack[icsq] << endl;
}
+ // 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){
t0bestallSel += t0best/eT0best/eT0best;
sumWt0bestallSel += 1./eT0best/eT0best;
ngoodsetsSel++;
- ngoodtrktrulyused+=ntracksinsetmyCut;
+ ngoodtrktrulyused+=ntracksinsetmyCut;
+ // printf("T0 set = %f +/- %f\n",t0best,eT0best);
}
else{
// printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
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: "<<ngoodtrktrulyused<<endl;
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[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
fT0SigmaT0def[2]=ngoodtrkt0;
fT0SigmaT0def[3]=ngoodtrktrulyused;
}
// cout << "AliTOFT0v1:" << endl ;
//}
- if(fT0SigmaT0def[1] < 0.01) fT0SigmaT0def[1] = 0.6;
+ 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]));
+
+ bench->Reset();
return fT0SigmaT0def;
}
};
Double_t mass=kMasses[index+2];
- Double_t dpp=0.02; //mean relative pt resolution;
- // if(mom > 1) dpp = 0.02*mom;
- 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;
}