// ROOT includes
#include <Riostream.h>
#include "TMath.h"
-#include "TClonesArray.h"
+#include "TObjArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
Double_t langaufun(Double_t *x, Double_t *par);
void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
-void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gSigma);
+void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2);
TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3);
TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2 = 0x0, const char* fitting = "");
//------------------------------------------------------------------------------------
-void MUONRecoCheck (Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
- const char* ocdbPath = "local://$ALICE_ROOT/OCDB", Int_t absorberRegion = -1)
+void MUONRecoCheck(Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
+ const char* ocdbPath = "local://$ALICE_ROOT/OCDB", const Double_t pMin = 0., const Double_t pMax = 300.,
+ const Int_t pNBins = 30, Int_t absorberRegion = -1)
{
/// Associate the reconstructed tracks with the simulated ones and check the quality of the reconstruction
/// (tracking/trigger efficiency; momentum, slope,... resolutions at first cluster and at vertex; cluster resolution).
- /// You can limit the calculation of track resolution at vertex to the tracks crossing the absorber in a given region
+ /// - You can choose the momentum range and number of bins used to study the track resolution versus momentum.
+ /// - You can limit the calculation of track resolution at vertex to the tracks crossing the absorber in a given region
/// with the flag "absorberRegion": -1=all, 1=[2,3]deg, 2=[3,10]deg.
Double_t aAbsLimits[2];
aAbsLimits[1] = 90.;
}
+ cout<<endl<<"Momentum range for track resolution studies: "<<pNBins<<" bins in ["<<pMin<<","<<pMax<<"] GeV/c"<<endl;
+ if (pMin >= pMax || pNBins <= 0) {
+ cout<<"--> incorrect momentum range"<<endl<<endl;
+ return;
+ } else cout<<endl;
+
AliLog::SetClassDebugLevel("AliMCEvent",-1);
// ###################################### define histograms ###################################### //
histoFile->mkdir("momentumAtVertex","momentumAtVertex");
histoFile->cd("momentumAtVertex");
- const Int_t pNBins = 30;
- const Double_t pEdges[2] = {0., 300.};
+ const Double_t pEdges[2] = {pMin, pMax};
const Int_t deltaPAtVtxNBins = 250;
- const Double_t deltaPAtVtxEdges[2] = {-35., 15.};
+ Double_t deltaPAtVtxEdges[2];
+ if (pMax < 50.) { deltaPAtVtxEdges[0] = -20.; deltaPAtVtxEdges[1] = 5.; }
+ else { deltaPAtVtxEdges[0] = -35.; deltaPAtVtxEdges[1] = 15.; }
TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
histoFile->cd("momentumAtFirstCluster");
const Int_t deltaPAtFirstClNBins = 500;
- const Double_t deltaPAtFirstClEdges[2] = {-25., 25.};
+ Double_t deltaPAtFirstClEdges[2];
+ if (pMax < 25.) { deltaPAtFirstClEdges[0] = -5.; deltaPAtFirstClEdges[1] = 5.; }
+ else if (pMax < 50.) { deltaPAtFirstClEdges[0] = -10.; deltaPAtFirstClEdges[1] = 10.; }
+ else { deltaPAtFirstClEdges[0] = -25.; deltaPAtFirstClEdges[1] = 25.; }
TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
const Int_t deltaPDCANBins = 500;
const Double_t deltaPDCAEdges[2] = {0., 1000.};
- const Double_t deltaPMCSAngEdges[2] = {-0.5, 0.5};
+ const Double_t deltaPMCSAngEdges[2] = {-1.5, 1.5};
TH1F *hPDCA = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
TH2D *hPDCAVsMom_2_3_Deg = new TH2D("hPDCAVsMom_2_3_Deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
TH2D *hPDCAVsAngle = new TH2D("hPDCAVsAngle","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
TH2D *hPDCAVsMCAngle = new TH2D("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
+ TGraphAsymmErrors* gMeanPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
+ gMeanPDCAVsMom_2_3_Deg->SetName("gMeanPDCAVsMom_2_3_Deg");
+ gMeanPDCAVsMom_2_3_Deg->SetTitle("<p #times DCA> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
TGraphAsymmErrors* gSigmaPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
gSigmaPDCAVsMom_2_3_Deg->SetName("gSigmaPDCAVsMom_2_3_Deg");
gSigmaPDCAVsMom_2_3_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
+ TGraphAsymmErrors* gMeanPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
+ gMeanPDCAVsMom_3_10_Deg->SetName("gMeanPDCAVsMom_3_10_Deg");
+ gMeanPDCAVsMom_3_10_Deg->SetTitle("<p #times DCA> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
TGraphAsymmErrors* gSigmaPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
gSigmaPDCAVsMom_3_10_Deg->SetName("gSigmaPDCAVsMom_3_10_Deg");
gSigmaPDCAVsMom_3_10_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
AliMUONTrackParam *trackParam;
Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
+ Double_t dPhi;
Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
Double_t xDCA,yDCA,dca,pU;
Double_t aMCSMoy = 0., aMCS2Moy = 0., dMCSMoy = 0., dMCS2Moy = 0., adMCSMoy = 0.;
eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
+ dPhi = phi2-phi1;
+ if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
+ else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
+
AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First()));
pU = trackParamAtDCA.P();
AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
hResSlopeYVertex->Fill(slopey2-slopey1);
hPDCA->Fill(0.5*(p2+pU)*dca);
hResEtaVertex->Fill(eta2-eta1);
- hResPhiVertex->Fill(phi2-phi1);
+ hResPhiVertex->Fill(dPhi);
if (aMC >= aAbsLimits[0] && aMC <= aAbsLimits[1]) {
hResMomVertexVsMom->Fill(p1,p2-p1);
hResSlopeXVertexVsMom->Fill(p1,slopex2-slopex1);
hResSlopeYVertexVsMom->Fill(p1,slopey2-slopey1);
hResEtaVertexVsMom->Fill(p1,eta2-eta1);
- hResPhiVertexVsMom->Fill(p1,phi2-phi1);
+ hResPhiVertexVsMom->Fill(p1,dPhi);
}
hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
if (aAbs > 2. && aAbs < 3.) {
hResSlopeYVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopey2-slopey1);
hPDCAVsPosAbsEnd_0_2_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
hResEtaVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,eta2-eta1);
- hResPhiVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,phi2-phi1);
+ hResPhiVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,dPhi);
}
else if (aMC >= 2. && aMC < 3) {
hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
hResSlopeYVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopey2-slopey1);
hPDCAVsPosAbsEnd_2_3_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
hResEtaVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,eta2-eta1);
- hResPhiVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,phi2-phi1);
+ hResPhiVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,dPhi);
}
else if (aMC >= 3. && aMC < 10.) {
hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
hResSlopeYVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopey2-slopey1);
hPDCAVsPosAbsEnd_3_10_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
hResEtaVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,eta2-eta1);
- hResPhiVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,phi2-phi1);
+ hResPhiVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,dPhi);
}
hResMomVertexVsAngle->Fill(aAbs,p2-p1);
hResSlopeXVertexVsAngle->Fill(aAbs,slopex2-slopex1);
hResSlopeYVertexVsAngle->Fill(aAbs,slopey2-slopey1);
hPDCAVsAngle->Fill(aAbs,0.5*(p2+pU)*dca);
hResEtaVertexVsAngle->Fill(aAbs,eta2-eta1);
- hResPhiVertexVsAngle->Fill(aAbs,phi2-phi1);
+ hResPhiVertexVsAngle->Fill(aAbs,dPhi);
hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
hResSlopeXVertexVsMCAngle->Fill(aMC,slopex2-slopex1);
hResSlopeYVertexVsMCAngle->Fill(aMC,slopey2-slopey1);
hPDCAVsMCAngle->Fill(aMC,0.5*(p2+pU)*dca);
hResEtaVertexVsMCAngle->Fill(aMC,eta2-eta1);
- hResPhiVertexVsMCAngle->Fill(aMC,phi2-phi1);
+ hResPhiVertexVsMCAngle->Fill(aMC,dPhi);
trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
x1 = trackParam->GetNonBendingCoor();
Double_t sigmaErr = f2->GetParError(3);
Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
hResMomVertexVsMom->GetXaxis()->SetRange(i-rebinFactorX+1,i);
- Double_t p = hResMomVertexVsMom->GetMean();
+ Double_t p = (tmp->GetEntries() > 0) ? hResMomVertexVsMom->GetMean() : 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
hResMomVertexVsMom->GetXaxis()->SetRange();
Double_t pErr[2] = {p-hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1), hResMomVertexVsMom->GetBinLowEdge(i+1)-p};
gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
FitGausResVsMom(hResSlopeYFirstClusterVsMom, pNBins, 0., 2.e-4, "slopeY residuals at first cluster", gMeanResSlopeYFirstClusterVsMom, gSigmaResSlopeYFirstClusterVsMom);
// compute p*DCA resolution in the region [2,3] deg at absorber end
- FitPDCAVsMom(hPDCAVsMom_2_3_Deg, pNBins, "p*DCA (tracks in [2,3] deg.)", gSigmaPDCAVsMom_2_3_Deg);
+ FitPDCAVsMom(hPDCAVsMom_2_3_Deg, pNBins, "p*DCA (tracks in [2,3] deg.)", gMeanPDCAVsMom_2_3_Deg, gSigmaPDCAVsMom_2_3_Deg);
// compute p*DCA resolution in the region [3,10] deg at absorber end
- FitPDCAVsMom(hPDCAVsMom_3_10_Deg, pNBins, "p*DCA (tracks in [3,10] deg.)", gSigmaPDCAVsMom_3_10_Deg);
+ FitPDCAVsMom(hPDCAVsMom_3_10_Deg, pNBins, "p*DCA (tracks in [3,10] deg.)", gMeanPDCAVsMom_3_10_Deg, gSigmaPDCAVsMom_3_10_Deg);
// compute MCS angular dispersion in the region [2,3] deg at absorber end
FitGausResVsMom(hPMCSAngVsMom_2_3_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsMom_2_3_Deg, gSigmaPMCSAngVsMom_2_3_Deg);
cResSlopeYVsPos->Write();
histoFile->cd("DCA");
+ gMeanPDCAVsMom_2_3_Deg->Write();
gSigmaPDCAVsMom_2_3_Deg->Write();
+ gMeanPDCAVsMom_3_10_Deg->Write();
gSigmaPDCAVsMom_3_10_Deg->Write();
gMeanPMCSAngVsMom_2_3_Deg->Write();
gSigmaPMCSAngVsMom_2_3_Deg->Write();
tmp->Rebin(rebin);
tmp->Fit("fGaus","NQ");
h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
- Double_t p = h->GetMean();
+ Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
h->GetXaxis()->SetRange();
Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
}
//------------------------------------------------------------------------------------
-void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gSigma)
+void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
{
/// generic function to fit p*DCA distributions
static TF1* fPGaus = 0x0;
for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
- fPGaus->SetParameters(1.,-100.,100.);
- Int_t rebin = 50.*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1)));
+ fPGaus->SetParameters(1.,0.,80.);
+ Int_t rebin = 25.*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1)));
while (tmp->GetNbinsX()%rebin!=0) rebin--;
tmp->Rebin(rebin);
tmp->Fit("fPGaus","NQ");
h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
- Double_t p = h->GetMean();
+ Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
h->GetXaxis()->SetRange();
Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
+ gMean->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(1));
+ gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(1), fPGaus->GetParError(1));
gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
delete tmp;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff