3 #include "dNdEtaAnalysis.h"
10 #include <TCollection.h>
11 #include <TIterator.h>
14 //____________________________________________________________________
15 ClassImp(dNdEtaAnalysis)
17 //____________________________________________________________________
18 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
21 hEtaVsVtx = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
22 hEtaVsVtx->SetXTitle("vtx z [cm]");
23 hEtaVsVtx->SetYTitle("#eta");
25 hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
26 hVtx = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
27 for (Int_t i=0; i<kVertexBinning; ++i)
29 hdNdEta[i] = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
30 hdNdEta[i]->SetYTitle("dN/d#eta");
37 //____________________________________________________________________
39 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
40 hEtaVsVtx->Fill(vtx, eta, weight);
41 hEtaVsVtxUncorrected->Fill(vtx,eta);
44 //____________________________________________________________________
46 dNdEtaAnalysis::FillEvent(Float_t vtx) {
50 //____________________________________________________________________
52 dNdEtaAnalysis::Finish() {
54 // first normalize with n events (per vtx)
55 for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) {
56 Float_t nEvents = hVtx->GetBinContent(i_vtx);
57 Float_t nEventsError = hVtx->GetBinError(i_vtx);
59 if (nEvents==0) continue;
61 for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
62 Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents;
63 if (value==0) continue;
64 Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents;
65 error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/
66 hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) +
67 TMath::Power(nEventsError/nEvents,2));
68 hEtaVsVtx->SetBinContent(i_vtx, i_eta, value);
69 hEtaVsVtx->SetBinError(i_vtx, i_eta, error);
73 // then take the wieghted average for each eta
74 // is this the right way to do it???
75 for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
77 Float_t sumError2 = 0;
78 Int_t nMeasurements = 0;
83 // do we have several histograms for different vertex positions?
84 Int_t vertexBinWidth = hVtx->GetNbinsX() / kVertexBinning;
85 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
87 Int_t vertexBinBegin = vertexBinWidth * vertexPos;
88 Int_t vertexBinEnd = vertexBinBegin + vertexBinWidth;
90 for (Int_t i_vtx=vertexBinBegin; i_vtx<=vertexBinEnd; i_vtx++) {
91 if (hVtx->GetBinContent(i_vtx)==0) continue;
92 if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue;
94 Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
95 sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w;
98 sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta);
99 sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
105 if (nMeasurements!=0) {
106 dndeta = sum/Float_t(nMeasurements);
107 error = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
109 dndeta = sum_xw/sum_w;
110 error = 1/TMath::Sqrt(sum_w);
112 dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(i_eta);
113 error = error/hdNdEta[vertexPos]->GetBinWidth(i_eta);
115 hdNdEta[vertexPos]->SetBinContent(i_eta, dndeta);
116 hdNdEta[vertexPos]->SetBinError(i_eta, error);
123 //____________________________________________________________________
125 dNdEtaAnalysis::SaveHistograms() {
127 gDirectory->mkdir(GetName());
128 gDirectory->cd(GetName());
132 for (Int_t i=0; i<kVertexBinning; ++i)
133 hdNdEta[i] ->Write();
135 gDirectory->cd("../");
138 //____________________________________________________________________
139 void dNdEtaAnalysis::DrawHistograms()
141 TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
142 canvas->Divide(2, 2);
146 hEtaVsVtx->Draw("COLZ");
149 if (hEtaVsVtxUncorrected)
150 hEtaVsVtxUncorrected->Draw("COLZ");
161 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
163 // Merges a list of dNdEtaAnalysis objects with this one.
164 // This is needed for PROOF.
165 // Returns the number of merged objects (including this)
173 TIterator* iter = list->MakeIterator();
177 const Int_t nCollections = kVertexBinning + 3;
178 TList* collections[nCollections];
179 for (Int_t i=0; i<nCollections; ++i)
180 collections[i] = new TList;
183 while ((obj = iter->Next()))
185 dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
189 collections[0]->Add(entry->hEtaVsVtx);
190 collections[1]->Add(entry->hEtaVsVtxUncorrected);
191 collections[2]->Add(entry->hVtx);
193 for (Int_t i=0; i<kVertexBinning; ++i)
194 collections[3+i]->Add(entry->hdNdEta[i]);
199 hEtaVsVtx->Merge(collections[0]);
200 hEtaVsVtxUncorrected->Merge(collections[1]);
201 hVtx->Merge(collections[2]);
202 for (Int_t i=0; i<kVertexBinning; ++i)
203 hdNdEta[i]->Merge(collections[3+i]);
205 for (Int_t i=0; i<nCollections; ++i)
206 delete collections[i];