fhCorr->Divide(fhGene, fhMeas, 1, 1, "B");
+ Int_t emptyBins = 0;
+ for (Int_t x=1; x<=fhCorr->GetNbinsX(); ++x)
+ for (Int_t y=1; y<=fhCorr->GetNbinsY(); ++y)
+ for (Int_t z=1; z<=fhCorr->GetNbinsZ(); ++z)
+ if (fhCorr->GetBinContent(x, y, z) == 0)
+ ++emptyBins;
+
+ if (emptyBins > 0)
+ printf("INFO: In %s we have %d empty bins (of %d) in the correction map\n", GetName(), emptyBins, fhCorr->GetNbinsX() * fhCorr->GetNbinsY() * fhCorr->GetNbinsZ());
}
//____________________________________________________________________
//
// loads the histograms from a file
//
-
- TFile* fin = TFile::Open(fileName);
-
+
+ TFile* fin = TFile::Open(fileName);
+
if(!fin) {
//Info("LoadHistograms",Form(" %s file does not exist",fileName));
return kFALSE;
class AliCorrectionMatrix : public TNamed
{
-public:
+protected: // do not create this baseclass
AliCorrectionMatrix();
AliCorrectionMatrix(const Char_t* name, const Char_t* title);
AliCorrectionMatrix(const AliCorrectionMatrix& c);
virtual ~AliCorrectionMatrix();
+public:
AliCorrectionMatrix& operator=(const AliCorrectionMatrix& corrMatrix);
virtual void Copy(TObject& c) const;
virtual Long64_t Merge(TCollection* list);
fhCorr->Sumw2();
}
+AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
+ Int_t nBinX, Float_t Xmin, Float_t Xmax,
+ Int_t nBinY, Float_t Ymin, Float_t Ymax,
+ Int_t nBinZ, const Float_t* zbins)
+ : AliCorrectionMatrix(name, title)
+{
+ // constructor with variable bin sizes
+
+ Float_t* binLimitsX = new Float_t[nBinX+1];
+ for (Int_t i=0; i<=nBinX; ++i)
+ binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
+
+ Float_t* binLimitsY = new Float_t[nBinY+1];
+ for (Int_t i=0; i<=nBinY; ++i)
+ binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
+
+ fhMeas = new TH3F(Form("meas_%s",name), Form("meas_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
+ fhGene = new TH3F(Form("gene_%s",name), Form("gene_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
+ fhCorr = new TH3F(Form("corr_%s",name), Form("corr_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
+
+ delete[] binLimitsX;
+ delete[] binLimitsY;
+
+ fhMeas->Sumw2();
+ fhGene->Sumw2();
+ fhCorr->Sumw2();
+}
+
//____________________________________________________________________
AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
{
AliCorrectionMatrix::SaveHistograms();
- AliPWG0Helper::CreateProjections(GetMeasuredHistogram());
- AliPWG0Helper::CreateProjections(GetGeneratedHistogram());
+ //AliPWG0Helper::CreateProjections(GetMeasuredHistogram());
+ //AliPWG0Helper::CreateProjections(GetGeneratedHistogram());
if (GetCorrectionHistogram())
- AliPWG0Helper::CreateProjections(GetCorrectionHistogram());
+ AliPWG0Helper::CreateDividedProjections(GetMeasuredHistogram(), GetGeneratedHistogram());
}
Int_t nBinY=10, Float_t Ymin=0., Float_t Ymax=10.,
Int_t nBinZ=10, Float_t Zmin=0., Float_t Zmax=10.);
+ AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
+ Int_t nBinX, Float_t Xmin, Float_t Xmax,
+ Int_t nBinY, Float_t Ymin, Float_t Ymax,
+ Int_t nBinZ, const Float_t* zbins);
+
virtual ~AliCorrectionMatrix3D();
TH3F* GetGeneratedHistogram();
proj->SetXTitle(hist->GetYaxis()->GetTitle());
proj->SetYTitle(hist->GetZaxis()->GetTitle());
}
+
+//____________________________________________________________________
+void AliPWG0Helper::CreateDividedProjections(TH3F* hist, TH3F* hist2, const char* axis)
+{
+ // create projections of the 3d hists divides them
+ // axis decides to which plane, if axis is 0 to all planes
+ // the histograms are not returned, just use them from memory or use this to create them in a file
+
+ if (axis == 0)
+ {
+ CreateDividedProjections(hist, hist2, "yx");
+ CreateDividedProjections(hist, hist2, "zx");
+ CreateDividedProjections(hist, hist2, "zy");
+
+ return;
+ }
+
+ TH1* proj = hist->Project3D(axis);
+ proj->SetXTitle(hist->GetXaxis()->GetTitle());
+ proj->SetYTitle(hist->GetYaxis()->GetTitle());
+
+ TH1* proj2 = hist2->Project3D(axis);
+ proj2->SetXTitle(hist2->GetXaxis()->GetTitle());
+ proj2->SetYTitle(hist2->GetYaxis()->GetTitle());
+
+ TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
+ division->Divide(proj2);
+}
+
static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries);
static void CreateProjections(TH3F* hist);
+ static void CreateDividedProjections(TH3F* hist, TH3F* hist2, const char* axis = 0);
protected:
ClassDef(AliPWG0Helper, 0)
AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() :
AliSelector(),
+ fdNdEtaAnalysisMBVtx(0),
+ fdNdEtaAnalysisMB(0),
+ fdNdEtaAnalysis(0),
fEsdTrackCuts(0),
fdNdEtaCorrection(0)
{
// Constructor. Initialization of pointers
//
- //AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
+ AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
}
AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
if (!fEsdTrackCuts && fInput)
fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
+ if (!fEsdTrackCuts && tree)
+ fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
+
if (!fEsdTrackCuts)
AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
+
+ if (!fdNdEtaCorrection && fInput)
+ fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fInput->FindObject("dndeta_correction"));
+
+ if (!fdNdEtaCorrection && fTree)
+ fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fTree->GetUserInfo()->FindObject("dndeta_correction"));
+
+ if (!fdNdEtaCorrection)
+ AliDebug(AliLog::kError, "ERROR: Could not read dndeta_correction from input list.");
+
+
+ fdNdEtaAnalysisMBVtx = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
+ fdNdEtaAnalysisMB = new dNdEtaAnalysis("dndeta_mb", "dndeta_mb");
fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
}
// read the user objects
AliSelector::Init(tree);
-
- if (!fEsdTrackCuts && fTree)
- fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fTree->GetUserInfo()->FindObject("AliESDtrackCuts"));
-
- if (!fEsdTrackCuts)
- AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info.");
-
- if (!fdNdEtaCorrection && fTree)
- fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fTree->GetUserInfo()->FindObject("dndeta_correction"));
}
Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
return kTRUE;
}
+ Float_t triggerCorr = fdNdEtaCorrection->GetTriggerCorrection(vtx[2], nGoodTracks);
+ if (triggerCorr <= 0)
+ {
+ AliDebug(AliLog::kError, Form("INFO: Skipping event because triggerCorr is <= 0 (%f)", triggerCorr));
+ delete list;
+ return kTRUE;
+ }
+
// loop over esd tracks
for (Int_t t=0; t<nGoodTracks; t++)
{
Float_t track2particleCorr = fdNdEtaCorrection->GetTrack2ParticleCorrection(vtx[2], eta, pt);
- Float_t weight = vertexRecoCorr * track2particleCorr;
+ Float_t weight = track2particleCorr * vertexRecoCorr * triggerCorr;
if (weight <= 0)
{
- AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f)", t, weight));
+ AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f) (vtx %f, eta %f, pt %f)", t, weight, vtx[2], eta, pt));
continue;
}
+ fdNdEtaAnalysisMBVtx->FillTrack(vtx[2], eta, pt, track2particleCorr);
+ fdNdEtaAnalysisMB->FillTrack(vtx[2], eta, pt, track2particleCorr * vertexRecoCorr);
fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, weight);
} // end of track loop
list = 0;
// for event count per vertex
- fdNdEtaAnalysis->FillEvent(vtx[2], vertexRecoCorr);
+ fdNdEtaAnalysisMBVtx->FillEvent(vtx[2], 1);
+ fdNdEtaAnalysisMB->FillEvent(vtx[2], vertexRecoCorr);
+ fdNdEtaAnalysis->FillEvent(vtx[2], vertexRecoCorr * triggerCorr);
return kTRUE;
}
if (fdNdEtaAnalysis)
fdNdEtaAnalysis->SaveHistograms();
+ if (fdNdEtaAnalysisMB)
+ fdNdEtaAnalysisMB->SaveHistograms();
+
+ if (fdNdEtaAnalysisMBVtx)
+ fdNdEtaAnalysisMBVtx->SaveHistograms();
+
if (fEsdTrackCuts)
fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
virtual void Terminate();
protected:
- dNdEtaAnalysis* fdNdEtaAnalysis; // contains the target histograms
+ dNdEtaAnalysis* fdNdEtaAnalysisMBVtx; // contains the histograms for the triggered events with vertex
+ dNdEtaAnalysis* fdNdEtaAnalysisMB; // contains the histograms corrected with vtx recon eff
+ dNdEtaAnalysis* fdNdEtaAnalysis; // contains the histograms corrected with vtx recon eff and trigger bias eff
+
AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
AlidNdEtaCorrection* fdNdEtaCorrection; // correction maps
// constructor
//
- fTrack2ParticleCorrection = new AliCorrectionMatrix3D("nTrackToNPart", "nTrackToNPart",80,-20,20,120,-6,6, 100, 0, 10);
+ Float_t binLimitsPt[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 10.0, 100.0};
- Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ fTrack2ParticleCorrection = new AliCorrectionMatrix3D("nTrackToNPart", "nTrackToNPart", 40, -20, 20, 60, -6, 6, 14, binLimitsPt);
+
+ Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
-
+
fVertexRecoCorrection = new AliCorrectionMatrix2D("vtxReco", "vtxReco",10,binLimitsVtx ,22,binLimitsN);
+ fTriggerCorrection = new AliCorrectionMatrix2D("trigger", "trigger",10,binLimitsVtx ,22,binLimitsN);
fTriggerBiasCorrection = new AliCorrectionMatrix2D("triggerBias", "triggerBias",120,-6,6,100, 0, 10);
fTrack2ParticleCorrection ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T}");
fVertexRecoCorrection ->SetAxisTitles("vtx z [cm]", "n particles/tracks/tracklets?");
+ fTriggerCorrection ->SetAxisTitles("vtx z [cm]", "n particles/tracks/tracklets?");
fTriggerBiasCorrection ->SetAxisTitles("#eta", "p_{T} [GeV/c]");
}
//
// divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
-
fTrack2ParticleCorrection->Divide();
+ TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
+ Int_t emptyBins = 0;
+ for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
+ for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
+ for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
+ if (hist->GetBinContent(x, y, z) == 0)
+ {
+ printf("Empty bin in fTrack2ParticleCorrection at vtx = %f, eta = %f, pt = %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
+ ++emptyBins;
+ }
+
+ printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
+
fVertexRecoCorrection->Divide();
+ fTriggerCorrection->Divide();
+ if (fNEvents == 0)
+ {
+ printf("ERROR: fNEvents is empty. Cannot scale histogram. Skipping processing of fTriggerBiasCorrection\n");
+ return;
+ }
fTriggerBiasCorrection->GetMeasuredHistogram()->Scale(Double_t(fNTriggeredEvents)/Double_t(fNEvents));
fTriggerBiasCorrection->Divide();
-
}
//____________________________________________________________________
fTrack2ParticleCorrection ->LoadHistograms(fileName, dir);
fVertexRecoCorrection ->LoadHistograms(fileName, dir);
+ fTriggerCorrection ->LoadHistograms(fileName, dir);
fTriggerBiasCorrection ->LoadHistograms(fileName, dir);
-
+
return kTRUE;
}
gDirectory->mkdir(fName.Data());
gDirectory->cd(fName.Data());
- fTrack2ParticleCorrection ->SaveHistograms();
- fVertexRecoCorrection ->SaveHistograms();
- fTriggerBiasCorrection ->SaveHistograms();
+ fTrack2ParticleCorrection->SaveHistograms();
+ fVertexRecoCorrection->SaveHistograms();
+ fTriggerCorrection->SaveHistograms();
+ fTriggerBiasCorrection->SaveHistograms();
gDirectory->cd("../");
}
fTrack2ParticleCorrection ->DrawHistograms();
fVertexRecoCorrection ->DrawHistograms();
+ fTriggerCorrection ->DrawHistograms();
fTriggerBiasCorrection ->DrawHistograms();
}
// calculates the fraction of particles measured (some are missed due to the pt cut off)
// uses the generated particle histogram from fTrack2ParticleCorrection
- TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
+ const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
// find eta borders, if eta is negative assume -0.8 ... 0.8
Int_t etaBegin = 0;
Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
+ ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
Int_t ptBin = ptProj->FindBin(ptCutOff);
Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
// - add documentation
// - add status: generate or use maps
// - add functionality to set the bin sizes
+// - update MERge function
//
#include <TNamed.h>
public:
AlidNdEtaCorrection(Char_t* name="dndeta_correction");
- // fVertexRecoCorrection
- void FillEvent(Float_t vtx, Float_t n) {fVertexRecoCorrection->FillGene(vtx, n);}
- void FillEventWithReconstructedVertex(Float_t vtx, Float_t n) {fVertexRecoCorrection->FillMeas(vtx, n);}
+ // fVertexRecoCorrection, fTriggerCorrection
+ void FillEvent(Float_t vtx, Float_t n) {fTriggerCorrection->FillGene(vtx, n);}
+ void FillEventWithTrigger(Float_t vtx, Float_t n) {fVertexRecoCorrection->FillGene(vtx, n); fTriggerCorrection->FillMeas(vtx, n);}
+ void FillEventWithTriggerWithReconstructedVertex(Float_t vtx, Float_t n) {fVertexRecoCorrection->FillMeas(vtx, n);}
// fTrack2ParticleCorrection
void FillParticle(Float_t vtx, Float_t eta, Float_t pt) {fTrack2ParticleCorrection->FillGene(vtx, eta, pt);}
Float_t GetVertexRecoCorrection(Float_t vtx, Float_t n) {return fVertexRecoCorrection->GetCorrection(vtx, n);}
+ Float_t GetTriggerCorrection(Float_t vtx, Float_t n) {return fTriggerCorrection->GetCorrection(vtx, n);}
+
Float_t GetTriggerBiasCorrection(Float_t eta, Float_t pt=0) {return fTriggerBiasCorrection->GetCorrection(eta, pt);}
Float_t GetMeasuredFraction(Float_t ptCutOff, Float_t eta = -1, Bool_t debug = kFALSE);
protected:
AliCorrectionMatrix3D* fTrack2ParticleCorrection; // handles the track-to-particle correction, function of vtx_z, eta, pt
AliCorrectionMatrix2D* fVertexRecoCorrection; // handles the vertex reconstruction efficiency, function of n_clustersITS and vtx_z
+ AliCorrectionMatrix2D* fTriggerCorrection; // handles the trigger efficiency efficiency, function of n_clustersITS and vtx_z
AliCorrectionMatrix2D* fTriggerBiasCorrection; // MB to desired sample
} // end of track loop
+ fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
if (eventTriggered)
{
- fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
+ fdNdEtaCorrection->FillEventWithTrigger(vtxMC[2], nGoodTracks);
if (vertexReconstructed)
- fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
+ fdNdEtaCorrection->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
}
return kTRUE;
// this macro combines the correction and the analysis and draws them
-void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd.root")
+void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd.root", const char* object = "dndeta")
{
gSystem->Load("libPWG0base");
//dNdEtaCorrection->RemoveEdges(2, 0, 2);
}
- fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+ fdNdEtaAnalysis = new dNdEtaAnalysis(object, object);
TFile* file = TFile::Open(filename);
if (!file)
}
fdNdEtaAnalysis->LoadHistograms();
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, (correct) ? 0.3 : -1);
fdNdEtaAnalysis->DrawHistograms();
}
//____________________________________________________________________
ClassImp(dNdEtaAnalysis)
+//____________________________________________________________________
+dNdEtaAnalysis::dNdEtaAnalysis() :
+ TNamed(),
+ fData(0),
+ fDataUncorrected(0),
+ fVtx(0),
+ fPtDist(0)
+{
+ // default constructor
+
+ for (Int_t i=0; i<kVertexBinning; ++i)
+ {
+ fdNdEta[i] = 0;
+ fdNdEtaPtCutOffCorrected[i] = 0;
+ }
+}
+
//____________________________________________________________________
dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
TNamed(name, title),
fData(0),
fDataUncorrected(0),
- fNEvents(0),
- fVtx(0)
+ fVtx(0),
+ fPtDist(0)
{
// constructor
- fData = new TH3F(Form("%s_analysis", name),"dNdEtaAnalysis",80,-20,20,120,-6,6,100, 0, 10);
+ fData = new TH3F(Form("%s_analysis", name),"Input data",80,-20,20,120,-6,6,100, 0, 10);
fData->SetXTitle("vtx z [cm]");
fData->SetYTitle("#eta");
fData->SetZTitle("p_{T}");
fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
+ fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle()));
fVtx = dynamic_cast<TH1D*> (fData->Project3D("x"));
+ fVtx->SetName(Form("%s_vtx", name));
+ fVtx->SetTitle("Vertex distribution");
+ fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
+ fdNdEta[0]->SetName(Form("%s_dNdEta", name));
+ fdNdEta[0]->SetTitle("dN/d#eta");
+ fdNdEta[0]->GetXaxis()->SetTitle(fData->GetYaxis()->GetTitle());
+ fdNdEta[0]->SetYTitle("dN/d#eta");
+
+ fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
+
for (Int_t i=1; i<kVertexBinning; ++i)
{
fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
- fdNdEta[i]->SetYTitle("dN/d#eta");
+ fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
}
+ fPtDist = dynamic_cast<TH1D*> (fData->Project3D("z"));
+ fPtDist->SetName(Form("%s_pt", name));
+ fPtDist->SetTitle("p_{T}");
+ fPtDist->GetXaxis()->SetTitle(fData->GetZaxis()->GetTitle());
+ fPtDist->SetYTitle("#frac{1}{p_{T}} #frac{dN}{d#eta dp_{T}}");
+
fData->Sumw2();
fVtx->Sumw2();
}
{
delete fdNdEta[i];
fdNdEta[i] = 0;
+ delete fdNdEtaPtCutOffCorrected[i];
+ fdNdEtaPtCutOffCorrected[i] = 0;
}
+
+ delete fPtDist;
+ fPtDist = 0;
}
//_____________________________________________________________________________
TNamed(c),
fData(0),
fDataUncorrected(0),
- fNEvents(0),
- fVtx(0)
+ fVtx(0),
+ fPtDist(0)
{
//
// dNdEtaAnalysis copy constructor
target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
for (Int_t i=0; i<kVertexBinning; ++i)
+ {
target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
+ target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[i]->Clone());
+ }
- target.fNEvents = fNEvents;
+ target.fPtDist = dynamic_cast<TH1D*> (fPtDist->Clone());
TNamed::Copy((TNamed &) c);
}
// fills an event into the histograms
fVtx->Fill(vtx, weight);
-
- fNEvents += weight;
}
//____________________________________________________________________
// In fData we have the track2particle and vertex reconstruction efficiency correction already applied
+
+ // create pt hist
+ {
+ // reset all ranges
+ fData->GetXaxis()->SetRange(0, 0);
+ fData->GetYaxis()->SetRange(0, 0);
+ fData->GetZaxis()->SetRange(0, 0);
+
+ // vtx cut
+ Int_t vertexBinBegin = fData->GetXaxis()->FindBin(-5);
+ Int_t vertexBinEnd = fData->GetXaxis()->FindBin(5);
+ fData->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
+ Float_t nEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd);
+
+ // eta cut
+ fData->GetYaxis()->SetRange(fData->GetYaxis()->FindBin(-0.8), fData->GetYaxis()->FindBin(0.8));
+ Float_t etaWidth = 1.6;
+
+ TH1D* ptHist = dynamic_cast<TH1D*> (fData->Project3D("ze"));
+
+ for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
+ {
+ Float_t binSize = fPtDist->GetBinWidth(i);
+ fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
+ fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
+ }
+
+ delete ptHist;
+ }
+
+ // reset all ranges
+ fData->GetXaxis()->SetRange(0, 0);
+ fData->GetYaxis()->SetRange(0, 0);
+ fData->GetZaxis()->SetRange(0, 0);
+
// integrate over pt (with pt cut)
- fData->GetZaxis()->SetRange(fData->GetZaxis()->FindBin(ptCut), fData->GetZaxis()->GetNbins());
- TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2"));
+ Int_t ptLowBin = 1;
+ if (ptCut > 0)
+ ptLowBin = fData->GetZaxis()->FindBin(ptCut);
+
+ fData->GetZaxis()->SetRange(ptLowBin, fData->GetZaxis()->GetNbins());
+ TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2e"));
+ vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
+ vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle());
+
if (vtxVsEta == 0)
{
printf("ERROR: pt integration failed\n");
}
}
- Float_t ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
- //ptCutOffCorrection = 1;
+ Float_t ptCutOffCorrection = 1;
+ if (correction)
+ ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
+
if (ptCutOffCorrection <= 0)
{
printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
continue;
}
- Float_t dndeta = sum / totalEvents / ptCutOffCorrection;
- Float_t error = TMath::Sqrt(sumError2) / totalEvents / ptCutOffCorrection;
+ Float_t dndeta = sum / totalEvents;
+ Float_t error = TMath::Sqrt(sumError2) / totalEvents;
dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
error = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
fdNdEta[vertexPos]->SetBinError(iEta, error);
+
+ dndeta /= ptCutOffCorrection;
+ error /= ptCutOffCorrection;
+
+ fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(iEta, dndeta);
+ fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(iEta, error);
}
}
}
gDirectory->mkdir(GetName());
gDirectory->cd(GetName());
- fData ->Write();
- AliPWG0Helper::CreateProjections(fData);
- fDataUncorrected->Write();
- AliPWG0Helper::CreateProjections(fDataUncorrected);
+ if (fData)
+ {
+ fData->Write();
+ AliPWG0Helper::CreateProjections(fData);
+ }
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
+
+ if (fDataUncorrected)
+ {
+ fDataUncorrected->Write();
+ AliPWG0Helper::CreateProjections(fDataUncorrected);
+ }
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fDataUncorrected is 0\n");
+
+ if (fData)
+ fVtx ->Write();
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fVtx is 0\n");
- fVtx ->Write();
for (Int_t i=0; i<kVertexBinning; ++i)
- fdNdEta[i] ->Write();
+ {
+ if (fdNdEta[i])
+ fdNdEta[i]->Write();
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
+
+ if (fdNdEtaPtCutOffCorrected[i])
+ fdNdEtaPtCutOffCorrected[i]->Write();
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
+ }
gDirectory->cd("../");
}
fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
for (Int_t i=0; i<kVertexBinning; ++i)
+ {
fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
+ fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
+ }
gDirectory->cd("../");
}
if (fData)
fData->Draw("COLZ");
- canvas->cd(2);
+ /*canvas->cd(2);
if (fDataUncorrected)
- fDataUncorrected->Draw("COLZ");
+ fDataUncorrected->Draw("COLZ");*/
- canvas->cd(3);
+ canvas->cd(2);
if (fVtx)
fVtx->Draw();
- canvas->cd(4);
+ canvas->cd(3);
+ if (fdNdEtaPtCutOffCorrected[0])
+ fdNdEtaPtCutOffCorrected[0]->Draw();
+
if (fdNdEta[0])
- fdNdEta[0]->Draw();
+ {
+ fdNdEta[0]->SetLineColor(kRed);
+ fdNdEta[0]->Draw("SAME");
+ }
+
+ canvas->cd(4);
+ if (fPtDist)
+ fPtDist->Draw();
// histograms for different vertices?
if (kVertexBinning > 0)
{
//canvas2->cd(i-1);
//printf("%d\n", i);
- if (fdNdEta[i])
+ if (fdNdEtaPtCutOffCorrected[i])
{
- fdNdEta[i]->SetLineColor(i+1);
- fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
- legend->AddEntry(fdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
+ fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
+ fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
+ legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
}
}
TObject* obj;
// sub collections
- const Int_t nCollections = kVertexBinning + 3;
+ const Int_t nCollections = 2 * kVertexBinning + 4; // 4 standalone hists, two arrays of size kVertexBinning
TList* collections[nCollections];
for (Int_t i=0; i<nCollections; ++i)
collections[i] = new TList;
collections[0]->Add(entry->fData);
collections[1]->Add(entry->fDataUncorrected);
collections[2]->Add(entry->fVtx);
+ collections[3]->Add(entry->fPtDist);
for (Int_t i=0; i<kVertexBinning; ++i)
- collections[3+i]->Add(entry->fdNdEta[i]);
+ {
+ collections[4+i]->Add(entry->fdNdEta[i]);
+ collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
+ }
++count;
}
fData->Merge(collections[0]);
fDataUncorrected->Merge(collections[1]);
fVtx->Merge(collections[2]);
+ fPtDist->Merge(collections[3]);
for (Int_t i=0; i<kVertexBinning; ++i)
- fdNdEta[i]->Merge(collections[3+i]);
+ {
+ fdNdEta[i]->Merge(collections[4+i]);
+ fdNdEta[i]->Merge(collections[4+kVertexBinning+i]);
+ }
for (Int_t i=0; i<nCollections; ++i)
delete collections[i];
// - add functionality to set the bin sizes
// - figure out correct way to treat the errors
// - add functionality to make dn/deta for different mult classes?
-// - implement destructor
#include <TNamed.h>
class dNdEtaAnalysis : public TNamed
{
public:
- enum { kVertexBinning = 1+4 }; // the first is for the whole vertex range, the others divide the vertex range
+ enum { kVertexBinning = 1+3 }; // the first is for the whole vertex range, the others divide the vertex range
+ dNdEtaAnalysis();
dNdEtaAnalysis(Char_t* name, Char_t* title);
virtual ~dNdEtaAnalysis();
TH3F* fData; // histogram Eta vs vtx (track count)
TH3F* fDataUncorrected; // uncorrected histograms Eta vs vtx (track count)
- Float_t fNEvents; // number of events (float because corrected by weight)
-
TH1D* fVtx; // vtx histogram (event count)
- TH1D* fdNdEta[kVertexBinning];// dndeta results for different vertex bins (0 = full range)
+
+ TH1D* fPtDist; // pt distribution
+
+ TH1D* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
+ TH1D* fdNdEtaPtCutOffCorrected[kVertexBinning]; // dndeta results for different vertex bins (0 = full range), pt cut off corrected
ClassDef(dNdEtaAnalysis, 0)
};
{
gSystem->Load("libPWG0base");
- dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
+ AlidNdEtaCorrection* dNdEtaMap = new AlidNdEtaCorrection();
dNdEtaMap->LoadCorrection("correction_map.root");
-
+
dNdEtaMap->DrawHistograms();
+
+ dNdEtaMap->GetMeasuredFraction(0.3, -1, kTRUE);
}
--- /dev/null
+/* $Id$ */
+
+void drawPlots()
+//void Track2Particle2D()
+{
+ TFile* file = TFile::Open("correction_map.root");
+
+ TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_nTrackToNPart_yx"));
+ TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_nTrackToNPart_zx"));
+ TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_nTrackToNPart_zy"));
+
+ Prepare2DPlot(corrYX);
+ Prepare2DPlot(corrZX);
+ Prepare2DPlot(corrZY);
+
+ SetRanges(corrYX);
+ SetRanges(corrZX);
+ SetRanges(corrZY);
+
+ TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
+ canvas->Divide(3, 1);
+
+ canvas->cd(1);
+ InitPadCOLZ();
+ corrYX->Draw("COLZ");
+
+ canvas->cd(2);
+ InitPadCOLZ();
+ corrZX->Draw("COLZ");
+
+ canvas->cd(3);
+ InitPadCOLZ();
+ corrZY->Draw("COLZ");
+}
+
+void Track2Particle3D()
+{
+ // get left margin proper
+
+ TFile* file = TFile::Open("correction_map.root");
+
+ TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
+ TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
+ TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
+
+ gene->SetTitle("Generated Particles");
+ meas->SetTitle("Measured Tracks");
+ corr->SetTitle("Correction Factor");
+
+ Prepare3DPlot(gene);
+ Prepare3DPlot(meas);
+ Prepare3DPlot(corr);
+
+ TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 1200, 400);
+ canvas->Divide(3, 1);
+
+ canvas->cd(1);
+ InitPad();
+ gene->Draw();
+
+ canvas->cd(2);
+ meas->Draw();
+
+ canvas->cd(3);
+ corr->Draw();
+}
+
+void SetRanges(TH1* hist)
+{
+ SetRanges(hist->GetXaxis());
+ SetRanges(hist->GetYaxis());
+ SetRanges(hist->GetZaxis());
+}
+
+void SetRanges(TAxis* axis)
+{
+ if (strcmp(axis->GetTitle(), "#eta") == 0)
+ axis->SetRangeUser(-0.8, 0.79999);
+ if (strcmp(axis->GetTitle(), "p_{T}") == 0)
+ axis->SetRangeUser(0, 9.9999);
+ if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
+ axis->SetRangeUser(-10, 9.9999);
+}
+
+void Prepare3DPlot(TH3* hist)
+{
+ hist->GetXaxis()->SetTitleOffset(1.5);
+ hist->GetYaxis()->SetTitleOffset(1.5);
+ hist->GetZaxis()->SetTitleOffset(1.5);
+}
+
+void Prepare2DPlot(TH2* hist)
+{
+ hist->SetStats(kFALSE);
+}
+
+void InitPad()
+{
+ gPad->Range(0, 0, 1, 1);
+ gPad->SetLeftMargin(0);
+ gPad->SetRightMargin(0.05);
+ gPad->SetTopMargin(0.13);
+ gPad->SetBottomMargin(0.1);
+
+ gPad->SetGridx();
+ gPad->SetGridy();
+}
+
+void InitPadCOLZ()
+{
+ gPad->Range(0, 0, 1, 1);
+ gPad->SetRightMargin(0.15);
+ gPad->SetLeftMargin(0.12);
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff