From 5af556491f181d7d35d88e5a851ddced198f7bd6 Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Wed, 31 May 2006 17:31:20 +0000 Subject: [PATCH] o) new way of calculating the final dndeta o) debugged the dndeta histograms which show only a selection of the vertex range --- PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx | 24 ++++++- PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h | 7 +- PWG0/dNdEta/dNdEtaAnalysis.cxx | 72 +++++++++------------ PWG0/dNdEta/dNdEtaAnalysis.h | 2 +- PWG0/dNdEta/drawVertexRecEff.C | 57 ++++++++++++++++ 5 files changed, 116 insertions(+), 46 deletions(-) create mode 100644 PWG0/dNdEta/drawVertexRecEff.C diff --git a/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx b/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx index b486000878b..5b511787501 100644 --- a/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx +++ b/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -22,7 +23,9 @@ ClassImp(AlidNdEtaAnalysisMCSelector) AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() : AlidNdEtaAnalysisSelector(), - fVertex(0) + fVertex(0), + fPartEta(0), + fEvents(0) { // // Constructor. Initialization of pointers @@ -42,7 +45,9 @@ void AlidNdEtaAnalysisMCSelector::Init(TTree *tree) tree->SetBranchStatus("ESD", 0); - fVertex = new TH3F("vertex", "vertex", 50, -50, 50, 50, -50, 50, 50, -50, 50); + fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50); + fPartEta = new TH1F("dndeta_check", "dndeta_check", 120, -6, 6); + fPartEta->Sumw2(); } Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry) @@ -102,9 +107,13 @@ Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry) fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), 1); fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz()); + + fPartEta->Fill(particle->Eta()); } fdNdEtaAnalysis->FillEvent(vtxMC[2]); + ++fEvents; + return kTRUE; } @@ -112,6 +121,15 @@ void AlidNdEtaAnalysisMCSelector::Terminate() { AlidNdEtaAnalysisSelector::Terminate(); - new TCanvas; + fPartEta->Scale(1.0/fEvents); + fPartEta->Scale(1.0/fPartEta->GetBinWidth(1)); + + TCanvas* canvas = new TCanvas("control", "control", 900, 450); + canvas->Divide(2, 1); + + canvas->cd(1); fVertex->Draw(); + + canvas->cd(2); + fPartEta->Draw(); } diff --git a/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h b/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h index 4f54eb8d77a..3aae19ec791 100644 --- a/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h +++ b/PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h @@ -6,6 +6,7 @@ #include "AlidNdEtaAnalysisSelector.h" class TH3F; +class TH1F; class AlidNdEtaAnalysisMCSelector : public AlidNdEtaAnalysisSelector { public: @@ -19,9 +20,11 @@ class AlidNdEtaAnalysisMCSelector : public AlidNdEtaAnalysisSelector { protected: private: - TH3F* fVertex; + TH3F* fVertex; //! vertex of counted particles + TH1F* fPartEta; //! counted particles as function of eta + Int_t fEvents; //! number of processed events - ClassDef(AlidNdEtaAnalysisMCSelector, 0); + ClassDef(AlidNdEtaAnalysisMCSelector, 0); }; #endif diff --git a/PWG0/dNdEta/dNdEtaAnalysis.cxx b/PWG0/dNdEta/dNdEtaAnalysis.cxx index ecb6f7bf67d..38d878de605 100644 --- a/PWG0/dNdEta/dNdEtaAnalysis.cxx +++ b/PWG0/dNdEta/dNdEtaAnalysis.cxx @@ -83,7 +83,7 @@ void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction) } // normalize with n events (per vtx) - for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) { +/* for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) { Float_t nEvents = hVtx->GetBinContent(iVtx); Float_t nEventsError = hVtx->GetBinError(iVtx); @@ -111,60 +111,52 @@ void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction) hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value); hEtaVsVtxCheck->SetBinError(iVtx, iEta, error); } - } + }*/ // then take the wieghted average for each eta // is this the right way to do it??? - for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) { - Float_t sum = 0; - Float_t sumError2 = 0; - Int_t nMeasurements = 0; - - Float_t sumXw = 0; - Float_t sumW = 0; - + for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) + { // do we have several histograms for different vertex positions? Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1); for (Int_t vertexPos=0; vertexPosGetNbinsX(); + Int_t vertexBinBegin = 1; + Int_t vertexBinEnd = hVtx->GetNbinsX() + 1; // the first histogram is always for the whole vertex range if (vertexPos > 0) { - vertexBinBegin = vertexBinWidth * (vertexPos-1); + vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1); vertexBinEnd = vertexBinBegin + vertexBinWidth; } - for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) { - if (hVtx->GetBinContent(iVtx)==0) continue; - if (hEtaVsVtxCheck->GetBinContent(iVtx, iEta)==0) continue; - - Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2); - sumXw = sumXw + hEtaVsVtxCheck->GetBinContent(iVtx, iEta)*w; - sumW = sumW + w; - - sum = sum + hEtaVsVtxCheck->GetBinContent(iVtx, iEta); - sumError2 = sumError2 + TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta),2); - nMeasurements++; + Float_t totalEvents = hVtx->Integral(vertexBinBegin, vertexBinEnd - 1); + if (totalEvents == 0) + { + printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd); + continue; } - Float_t dndeta = 0; - Float_t error = 0; - if (nMeasurements!=0) { - dndeta = sum/Float_t(nMeasurements); - error = TMath::Sqrt(sumError2)/Float_t(nMeasurements); + Float_t sum = 0; + Float_t sumError2 = 0; + for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++) + { + if (hEtaVsVtx->GetBinContent(iVtx, iEta) != 0) + { + sum = sum + hEtaVsVtx->GetBinContent(iVtx, iEta); + sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2); + } + } - dndeta = sumXw/sumW; - error = 1/TMath::Sqrt(sumW); + Float_t dndeta = sum / totalEvents; + Float_t error = TMath::Sqrt(sumError2) / totalEvents; - dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta); - error = error/hdNdEta[vertexPos]->GetBinWidth(iEta); + dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta); + error = error/hdNdEta[vertexPos]->GetBinWidth(iEta); - hdNdEta[vertexPos]->SetBinContent(iEta, dndeta); - hdNdEta[vertexPos]->SetBinError(iEta, error); - } + hdNdEta[vertexPos]->SetBinContent(iEta, dndeta); + hdNdEta[vertexPos]->SetBinError(iEta, error); } } } @@ -227,15 +219,15 @@ void dNdEtaAnalysis::DrawHistograms() TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9); - for (Int_t i=1; icd(i-1); //printf("%d\n", i); if (hdNdEta[i]) { - hdNdEta[i]->SetLineColor(i); - hdNdEta[i]->Draw((i == 1) ? "" : "SAME"); - legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1)); + hdNdEta[i]->SetLineColor(i+1); + hdNdEta[i]->Draw((i == 0) ? "" : "SAME"); + legend->AddEntry(hdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1)); } } diff --git a/PWG0/dNdEta/dNdEtaAnalysis.h b/PWG0/dNdEta/dNdEtaAnalysis.h index 754f7978206..ef90a1358fa 100644 --- a/PWG0/dNdEta/dNdEtaAnalysis.h +++ b/PWG0/dNdEta/dNdEtaAnalysis.h @@ -28,7 +28,7 @@ class dNdEtaCorrection; class dNdEtaAnalysis : public TNamed { public: - enum { kVertexBinning = 1+6 }; // the first is for the whole vertex range, the others divide the vertex range + enum { kVertexBinning = 1+4 }; // the first is for the whole vertex range, the others divide the vertex range dNdEtaAnalysis(Char_t* name, Char_t* title); diff --git a/PWG0/dNdEta/drawVertexRecEff.C b/PWG0/dNdEta/drawVertexRecEff.C new file mode 100644 index 00000000000..0b5ae24a715 --- /dev/null +++ b/PWG0/dNdEta/drawVertexRecEff.C @@ -0,0 +1,57 @@ +/* $Id$ */ + +// draws the result of AlidNdEtaVertexRecEffSelector + +void drawVertexRecEff() +{ + TFile* fout = TFile::Open("vertexRecEff.root"); + + TH1F* dNGen = dynamic_cast fout->Get("dNGen"); + TH1F* dNRec = dynamic_cast fout->Get("dNRec"); + + TH1F* vtxGen = dynamic_cast fout->Get("VtxGen"); + TH1F* vtxRec = dynamic_cast fout->Get("VtxRec"); + + // calculate ratios + TH1F* dNRatiodN = dynamic_cast (dNRec->Clone("dNRatiodN")); + dNRatiodN->SetTitle("Ratio"); + dNRatiodN->Divide(dNGen); + + //vtxGen->Rebin(4); + //vtxRec->Rebin(4); + + // calculate correction ratio number + Float_t sumGen = 0; + Float_t sumRec = 0; + for (Int_t i=1; i<=dNGen->GetNbinsX(); ++i) + { + sumGen += dNGen->GetBinCenter(i) * dNGen->GetBinContent(i); + sumRec += dNRec->GetBinCenter(i) * dNRec->GetBinContent(i); + } + Float_t ratio = sumRec / dNRec->Integral(1, dNGen->GetNbinsX()) * dNGen->Integral(1, dNGen->GetNbinsX()) / sumGen; + + cout << "Ratio: " << ratio << endl; + + TH1F* dNRatioVtx = dynamic_cast (vtxRec->Clone("dNRatioVtx")); + dNRatioVtx->SetTitle("Ratio"); + dNRatioVtx->Divide(vtxGen); + + TCanvas* canvas = new TCanvas("dN", "dN", 1000, 1000); + canvas->Divide(2, 2); + + canvas->cd(1); + dNGen->Draw(); + dNRec->SetLineColor(kRed); + dNRec->Draw("SAME"); + + canvas->cd(2); + dNRatiodN->Draw(); + + canvas->cd(3); + vtxGen->Draw(); + vtxRec->SetLineColor(kRed); + vtxRec->Draw("SAME"); + + canvas->cd(4); + dNRatioVtx->Draw(); +} -- 2.43.0