From 1afae8ff2009d0df7b55b38dec35c3c08058dabd Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Wed, 14 Jun 2006 17:40:55 +0000 Subject: [PATCH] a lot of work on the analysis --- PWG0/AliCorrectionMatrix.cxx | 15 +- PWG0/AliCorrectionMatrix.h | 3 +- PWG0/AliCorrectionMatrix3D.cxx | 34 +++- PWG0/AliCorrectionMatrix3D.h | 5 + PWG0/AliPWG0Helper.cxx | 29 +++ PWG0/AliPWG0Helper.h | 1 + PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx | 54 +++-- PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h | 5 +- PWG0/dNdEta/AlidNdEtaCorrection.cxx | 45 ++++- PWG0/dNdEta/AlidNdEtaCorrection.h | 11 +- PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx | 5 +- PWG0/dNdEta/CreatedNdEta.C | 6 +- PWG0/dNdEta/dNdEtaAnalysis.cxx | 202 ++++++++++++++++--- PWG0/dNdEta/dNdEtaAnalysis.h | 12 +- PWG0/dNdEta/drawCorrection.C | 6 +- PWG0/dNdEta/drawPlots.C | 114 +++++++++++ 16 files changed, 468 insertions(+), 79 deletions(-) create mode 100644 PWG0/dNdEta/drawPlots.C diff --git a/PWG0/AliCorrectionMatrix.cxx b/PWG0/AliCorrectionMatrix.cxx index 85070ddecb3..04686cef5b2 100644 --- a/PWG0/AliCorrectionMatrix.cxx +++ b/PWG0/AliCorrectionMatrix.cxx @@ -162,6 +162,15 @@ void AliCorrectionMatrix::Divide() 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()); } //____________________________________________________________________ @@ -170,9 +179,9 @@ Bool_t AliCorrectionMatrix::LoadHistograms(Char_t* fileName, Char_t* dir) // // 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; diff --git a/PWG0/AliCorrectionMatrix.h b/PWG0/AliCorrectionMatrix.h index 3507b586d70..e0e8ef6a08e 100644 --- a/PWG0/AliCorrectionMatrix.h +++ b/PWG0/AliCorrectionMatrix.h @@ -20,12 +20,13 @@ class TH1; 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); diff --git a/PWG0/AliCorrectionMatrix3D.cxx b/PWG0/AliCorrectionMatrix3D.cxx index b856ae6d42b..f69ea505acc 100644 --- a/PWG0/AliCorrectionMatrix3D.cxx +++ b/PWG0/AliCorrectionMatrix3D.cxx @@ -52,6 +52,34 @@ AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* t 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() { @@ -120,9 +148,9 @@ void AliCorrectionMatrix3D::SaveHistograms() AliCorrectionMatrix::SaveHistograms(); - AliPWG0Helper::CreateProjections(GetMeasuredHistogram()); - AliPWG0Helper::CreateProjections(GetGeneratedHistogram()); + //AliPWG0Helper::CreateProjections(GetMeasuredHistogram()); + //AliPWG0Helper::CreateProjections(GetGeneratedHistogram()); if (GetCorrectionHistogram()) - AliPWG0Helper::CreateProjections(GetCorrectionHistogram()); + AliPWG0Helper::CreateDividedProjections(GetMeasuredHistogram(), GetGeneratedHistogram()); } diff --git a/PWG0/AliCorrectionMatrix3D.h b/PWG0/AliCorrectionMatrix3D.h index f6619f6cb2b..7edadfec5b8 100644 --- a/PWG0/AliCorrectionMatrix3D.h +++ b/PWG0/AliCorrectionMatrix3D.h @@ -23,6 +23,11 @@ public: 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(); diff --git a/PWG0/AliPWG0Helper.cxx b/PWG0/AliPWG0Helper.cxx index 9f036a0e5b0..8f488a097f8 100644 --- a/PWG0/AliPWG0Helper.cxx +++ b/PWG0/AliPWG0Helper.cxx @@ -119,3 +119,32 @@ void AliPWG0Helper::CreateProjections(TH3F* hist) 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 (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName()))); + division->Divide(proj2); +} + diff --git a/PWG0/AliPWG0Helper.h b/PWG0/AliPWG0Helper.h index 3e79132e279..8b4516349b3 100644 --- a/PWG0/AliPWG0Helper.h +++ b/PWG0/AliPWG0Helper.h @@ -19,6 +19,7 @@ class AliPWG0Helper : public TObject 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) diff --git a/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx b/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx index e424e3fbda6..3ea1ea7501b 100644 --- a/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx +++ b/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx @@ -22,6 +22,9 @@ ClassImp(AlidNdEtaAnalysisESDSelector) AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() : AliSelector(), + fdNdEtaAnalysisMBVtx(0), + fdNdEtaAnalysisMB(0), + fdNdEtaAnalysis(0), fEsdTrackCuts(0), fdNdEtaCorrection(0) { @@ -29,7 +32,7 @@ AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() : // Constructor. Initialization of pointers // - //AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug); + AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug); } AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector() @@ -59,9 +62,25 @@ void AlidNdEtaAnalysisESDSelector::SlaveBegin(TTree* tree) if (!fEsdTrackCuts && fInput) fEsdTrackCuts = dynamic_cast (fInput->FindObject("AliESDtrackCuts")); + if (!fEsdTrackCuts && tree) + fEsdTrackCuts = dynamic_cast (tree->GetUserInfo()->FindObject("AliESDtrackCuts")); + if (!fEsdTrackCuts) AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list."); + + if (!fdNdEtaCorrection && fInput) + fdNdEtaCorrection = dynamic_cast (fInput->FindObject("dndeta_correction")); + + if (!fdNdEtaCorrection && fTree) + fdNdEtaCorrection = dynamic_cast (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"); } @@ -70,15 +89,6 @@ void AlidNdEtaAnalysisESDSelector::Init(TTree* tree) // read the user objects AliSelector::Init(tree); - - if (!fEsdTrackCuts && fTree) - fEsdTrackCuts = dynamic_cast (fTree->GetUserInfo()->FindObject("AliESDtrackCuts")); - - if (!fEsdTrackCuts) - AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info."); - - if (!fdNdEtaCorrection && fTree) - fdNdEtaCorrection = dynamic_cast (fTree->GetUserInfo()->FindObject("dndeta_correction")); } Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry) @@ -147,6 +157,14 @@ 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; tGetTrack2ParticleCorrection(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 @@ -181,7 +201,9 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry) 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; } @@ -225,6 +247,12 @@ void AlidNdEtaAnalysisESDSelector::Terminate() if (fdNdEtaAnalysis) fdNdEtaAnalysis->SaveHistograms(); + if (fdNdEtaAnalysisMB) + fdNdEtaAnalysisMB->SaveHistograms(); + + if (fdNdEtaAnalysisMBVtx) + fdNdEtaAnalysisMBVtx->SaveHistograms(); + if (fEsdTrackCuts) fEsdTrackCuts->SaveHistograms("esd_tracks_cuts"); diff --git a/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h b/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h index f894c661544..6dffb26f7bb 100644 --- a/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h +++ b/PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h @@ -21,7 +21,10 @@ class AlidNdEtaAnalysisESDSelector : public AliSelector { 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 diff --git a/PWG0/dNdEta/AlidNdEtaCorrection.cxx b/PWG0/dNdEta/AlidNdEtaCorrection.cxx index 6256bccfe0c..545ca2bae0f 100644 --- a/PWG0/dNdEta/AlidNdEtaCorrection.cxx +++ b/PWG0/dNdEta/AlidNdEtaCorrection.cxx @@ -18,18 +18,22 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(Char_t* name) // 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]"); } @@ -42,14 +46,31 @@ AlidNdEtaCorrection::Finish() { // // 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(); - } //____________________________________________________________________ @@ -107,8 +128,9 @@ AlidNdEtaCorrection::LoadHistograms(Char_t* fileName, Char_t* dir) { fTrack2ParticleCorrection ->LoadHistograms(fileName, dir); fVertexRecoCorrection ->LoadHistograms(fileName, dir); + fTriggerCorrection ->LoadHistograms(fileName, dir); fTriggerBiasCorrection ->LoadHistograms(fileName, dir); - + return kTRUE; } @@ -123,9 +145,10 @@ AlidNdEtaCorrection::SaveHistograms() { gDirectory->mkdir(fName.Data()); gDirectory->cd(fName.Data()); - fTrack2ParticleCorrection ->SaveHistograms(); - fVertexRecoCorrection ->SaveHistograms(); - fTriggerBiasCorrection ->SaveHistograms(); + fTrack2ParticleCorrection->SaveHistograms(); + fVertexRecoCorrection->SaveHistograms(); + fTriggerCorrection->SaveHistograms(); + fTriggerBiasCorrection->SaveHistograms(); gDirectory->cd("../"); } @@ -138,6 +161,7 @@ void AlidNdEtaCorrection::DrawHistograms() fTrack2ParticleCorrection ->DrawHistograms(); fVertexRecoCorrection ->DrawHistograms(); + fTriggerCorrection ->DrawHistograms(); fTriggerBiasCorrection ->DrawHistograms(); } @@ -148,7 +172,7 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, // 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; @@ -168,6 +192,7 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Int_t vertexEnd = generated->GetXaxis()->FindBin(10); TH1D* ptProj = dynamic_cast (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()); diff --git a/PWG0/dNdEta/AlidNdEtaCorrection.h b/PWG0/dNdEta/AlidNdEtaCorrection.h index 8c754b2612c..212c602fb81 100644 --- a/PWG0/dNdEta/AlidNdEtaCorrection.h +++ b/PWG0/dNdEta/AlidNdEtaCorrection.h @@ -15,6 +15,7 @@ // - add documentation // - add status: generate or use maps // - add functionality to set the bin sizes +// - update MERge function // #include @@ -27,9 +28,10 @@ class AlidNdEtaCorrection : public TNamed 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);} @@ -64,6 +66,8 @@ public: 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); @@ -71,6 +75,7 @@ public: 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 diff --git a/PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx b/PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx index 64b46cb8d72..9f37e7c0883 100644 --- a/PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx +++ b/PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx @@ -214,11 +214,12 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry) } // 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; diff --git a/PWG0/dNdEta/CreatedNdEta.C b/PWG0/dNdEta/CreatedNdEta.C index ce97337a786..19e207bb96f 100644 --- a/PWG0/dNdEta/CreatedNdEta.C +++ b/PWG0/dNdEta/CreatedNdEta.C @@ -1,6 +1,6 @@ // 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"); @@ -12,7 +12,7 @@ void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd //dNdEtaCorrection->RemoveEdges(2, 0, 2); } - fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta"); + fdNdEtaAnalysis = new dNdEtaAnalysis(object, object); TFile* file = TFile::Open(filename); if (!file) @@ -22,7 +22,7 @@ void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd } fdNdEtaAnalysis->LoadHistograms(); - fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3); + fdNdEtaAnalysis->Finish(dNdEtaCorrection, (correct) ? 0.3 : -1); fdNdEtaAnalysis->DrawHistograms(); } diff --git a/PWG0/dNdEta/dNdEtaAnalysis.cxx b/PWG0/dNdEta/dNdEtaAnalysis.cxx index fad1fa0fd0d..da5f3076639 100644 --- a/PWG0/dNdEta/dNdEtaAnalysis.cxx +++ b/PWG0/dNdEta/dNdEtaAnalysis.cxx @@ -19,31 +19,65 @@ //____________________________________________________________________ ClassImp(dNdEtaAnalysis) +//____________________________________________________________________ +dNdEtaAnalysis::dNdEtaAnalysis() : + TNamed(), + fData(0), + fDataUncorrected(0), + fVtx(0), + fPtDist(0) +{ + // default constructor + + for (Int_t i=0; iSetXTitle("vtx z [cm]"); fData->SetYTitle("#eta"); fData->SetZTitle("p_{T}"); fDataUncorrected = dynamic_cast (fData->Clone(Form("%s_analysis_uncorrected", name))); + fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle())); fVtx = dynamic_cast (fData->Project3D("x")); + fVtx->SetName(Form("%s_vtx", name)); + fVtx->SetTitle("Vertex distribution"); + fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle()); fdNdEta[0] = dynamic_cast (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 (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName()))); + for (Int_t i=1; i (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i))); - fdNdEta[i]->SetYTitle("dN/d#eta"); + fdNdEtaPtCutOffCorrected[i] = dynamic_cast (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i))); } + fPtDist = dynamic_cast (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(); } @@ -66,7 +100,12 @@ dNdEtaAnalysis::~dNdEtaAnalysis() { delete fdNdEta[i]; fdNdEta[i] = 0; + delete fdNdEtaPtCutOffCorrected[i]; + fdNdEtaPtCutOffCorrected[i] = 0; } + + delete fPtDist; + fPtDist = 0; } //_____________________________________________________________________________ @@ -74,8 +113,8 @@ dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) : TNamed(c), fData(0), fDataUncorrected(0), - fNEvents(0), - fVtx(0) + fVtx(0), + fPtDist(0) { // // dNdEtaAnalysis copy constructor @@ -109,9 +148,12 @@ void dNdEtaAnalysis::Copy(TObject &c) const target.fVtx = dynamic_cast (fVtx->Clone()); for (Int_t i=0; i (fdNdEta[i]->Clone()); + target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast (fdNdEtaPtCutOffCorrected[i]->Clone()); + } - target.fNEvents = fNEvents; + target.fPtDist = dynamic_cast (fPtDist->Clone()); TNamed::Copy((TNamed &) c); } @@ -131,8 +173,6 @@ void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight) // fills an event into the histograms fVtx->Fill(vtx, weight); - - fNEvents += weight; } //____________________________________________________________________ @@ -146,9 +186,51 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut) // 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 (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 (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 (fData->Project3D("yx2e")); + vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle()); + vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle()); + if (vtxVsEta == 0) { printf("ERROR: pt integration failed\n"); @@ -192,22 +274,30 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut) } } - 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); } } } @@ -220,14 +310,39 @@ void dNdEtaAnalysis::SaveHistograms() 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; iWrite(); + { + 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("../"); } @@ -244,7 +359,10 @@ void dNdEtaAnalysis::LoadHistograms() fVtx = dynamic_cast (gDirectory->Get(fVtx->GetName())); for (Int_t i=0; i (gDirectory->Get(fdNdEta[i]->GetName())); + fdNdEtaPtCutOffCorrected[i] = dynamic_cast (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName())); + } gDirectory->cd("../"); } @@ -261,17 +379,27 @@ void dNdEtaAnalysis::DrawHistograms() 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) @@ -288,11 +416,11 @@ void dNdEtaAnalysis::DrawHistograms() { //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)); } } @@ -316,7 +444,7 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list) 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; iAdd(entry->fData); collections[1]->Add(entry->fDataUncorrected); collections[2]->Add(entry->fVtx); + collections[3]->Add(entry->fPtDist); for (Int_t i=0; iAdd(entry->fdNdEta[i]); + { + collections[4+i]->Add(entry->fdNdEta[i]); + collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]); + } ++count; } @@ -341,8 +473,12 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list) fData->Merge(collections[0]); fDataUncorrected->Merge(collections[1]); fVtx->Merge(collections[2]); + fPtDist->Merge(collections[3]); for (Int_t i=0; iMerge(collections[3+i]); + { + fdNdEta[i]->Merge(collections[4+i]); + fdNdEta[i]->Merge(collections[4+kVertexBinning+i]); + } for (Int_t i=0; i @@ -28,8 +27,9 @@ class AlidNdEtaCorrection; 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(); @@ -57,10 +57,12 @@ protected: 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) }; diff --git a/PWG0/dNdEta/drawCorrection.C b/PWG0/dNdEta/drawCorrection.C index 4a819bc52d0..1ccf1f1fca5 100644 --- a/PWG0/dNdEta/drawCorrection.C +++ b/PWG0/dNdEta/drawCorrection.C @@ -4,8 +4,10 @@ void drawCorrection() { gSystem->Load("libPWG0base"); - dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection(); + AlidNdEtaCorrection* dNdEtaMap = new AlidNdEtaCorrection(); dNdEtaMap->LoadCorrection("correction_map.root"); - + dNdEtaMap->DrawHistograms(); + + dNdEtaMap->GetMeasuredFraction(0.3, -1, kTRUE); } diff --git a/PWG0/dNdEta/drawPlots.C b/PWG0/dNdEta/drawPlots.C new file mode 100644 index 00000000000..2065574430d --- /dev/null +++ b/PWG0/dNdEta/drawPlots.C @@ -0,0 +1,114 @@ +/* $Id$ */ + +void drawPlots() +//void Track2Particle2D() +{ + TFile* file = TFile::Open("correction_map.root"); + + TH2* corrYX = dynamic_cast (file->Get("dndeta_correction/corr_nTrackToNPart_yx")); + TH2* corrZX = dynamic_cast (file->Get("dndeta_correction/corr_nTrackToNPart_zx")); + TH2* corrZY = dynamic_cast (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 (file->Get("dndeta_correction/gene_nTrackToNPart")); + TH3* meas = dynamic_cast (file->Get("dndeta_correction/meas_nTrackToNPart")); + TH3* corr = dynamic_cast (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); +} -- 2.43.5