]>
Commit | Line | Data |
---|---|---|
dc740de4 | 1 | /* $Id$ */ |
2 | ||
75ec0f41 | 3 | #include "dNdEtaAnalysis.h" |
4 | ||
ceb5d1b5 | 5 | #include <TFile.h> |
6 | #include <TH2F.h> | |
7 | #include <TH1D.h> | |
8 | #include <TMath.h> | |
9 | #include <TCanvas.h> | |
7029240a | 10 | #include <TCollection.h> |
11 | #include <TIterator.h> | |
12 | #include <TList.h> | |
fcf2fb36 | 13 | #include <TLegend.h> |
14 | ||
15 | #include "dNdEtaCorrection.h" | |
ceb5d1b5 | 16 | |
75ec0f41 | 17 | //____________________________________________________________________ |
b7f4a1fd | 18 | ClassImp(dNdEtaAnalysis) |
75ec0f41 | 19 | |
20 | //____________________________________________________________________ | |
7029240a | 21 | dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) : |
22 | TNamed(name, title) | |
23 | { | |
24 | hEtaVsVtx = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6); | |
75ec0f41 | 25 | hEtaVsVtx->SetXTitle("vtx z [cm]"); |
26 | hEtaVsVtx->SetYTitle("#eta"); | |
4dd2ad81 | 27 | |
7029240a | 28 | hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name))); |
fcf2fb36 | 29 | hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name))); |
7029240a | 30 | hVtx = hEtaVsVtx->ProjectionX(Form("%s_vtx", name)); |
31 | for (Int_t i=0; i<kVertexBinning; ++i) | |
32 | { | |
33 | hdNdEta[i] = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i)); | |
34 | hdNdEta[i]->SetYTitle("dN/d#eta"); | |
35 | } | |
75ec0f41 | 36 | |
37 | hEtaVsVtx->Sumw2(); | |
38 | hVtx->Sumw2(); | |
75ec0f41 | 39 | } |
40 | ||
41 | //____________________________________________________________________ | |
42 | void | |
fcf2fb36 | 43 | dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) { |
ceb5d1b5 | 44 | hEtaVsVtxUncorrected->Fill(vtx,eta); |
fcf2fb36 | 45 | hEtaVsVtxCheck->Fill(vtx, eta, c); |
75ec0f41 | 46 | } |
47 | ||
48 | //____________________________________________________________________ | |
49 | void | |
50 | dNdEtaAnalysis::FillEvent(Float_t vtx) { | |
51 | hVtx->Fill(vtx); | |
52 | } | |
53 | ||
54 | //____________________________________________________________________ | |
fcf2fb36 | 55 | void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction) |
56 | { | |
57 | // correct with correction values if available | |
58 | // TODO what do we do with the error? | |
59 | if (!correction) | |
60 | printf("INFO: No correction applied\n"); | |
61 | ||
62 | // this can be replaced by TH2F::Divide if we agree that the binning will be always the same | |
63 | for (Int_t iVtx=0; iVtx<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++) | |
64 | { | |
65 | for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++) | |
66 | { | |
67 | Float_t correctionValue = 1; | |
68 | if (correction) | |
69 | correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta)); | |
70 | ||
71 | Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta); | |
72 | Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta); | |
73 | ||
74 | Float_t correctedValue = value * correctionValue; | |
75 | Float_t correctedError = error * correctionValue; | |
76 | ||
77 | if (correctedValue != 0) | |
78 | { | |
79 | hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue); | |
80 | hEtaVsVtx->SetBinError(iVtx, iEta, correctedError); | |
81 | } | |
82 | } | |
83 | } | |
84 | ||
85 | // normalize with n events (per vtx) | |
86 | for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) { | |
87 | Float_t nEvents = hVtx->GetBinContent(iVtx); | |
88 | Float_t nEventsError = hVtx->GetBinError(iVtx); | |
75ec0f41 | 89 | |
75ec0f41 | 90 | if (nEvents==0) continue; |
fcf2fb36 | 91 | |
92 | for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) { | |
93 | Float_t value = hEtaVsVtx->GetBinContent(iVtx, iEta) / nEvents; | |
94 | if (value==0) continue; | |
95 | Float_t error = hEtaVsVtx->GetBinError(iVtx, iEta)/nEvents; | |
96 | error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta)/ | |
97 | hEtaVsVtx->GetBinContent(iVtx, iEta),2) + | |
98 | TMath::Power(nEventsError/nEvents,2)); | |
99 | hEtaVsVtx->SetBinContent(iVtx, iEta, value); | |
100 | hEtaVsVtx->SetBinError(iVtx, iEta, error); | |
101 | } | |
102 | ||
103 | //debug | |
104 | for (Int_t iEta=0; iEta<=hEtaVsVtxCheck->GetNbinsY(); iEta++) { | |
105 | Float_t value = hEtaVsVtxCheck->GetBinContent(iVtx, iEta) / nEvents; | |
75ec0f41 | 106 | if (value==0) continue; |
fcf2fb36 | 107 | Float_t error = hEtaVsVtxCheck->GetBinError(iVtx, iEta)/nEvents; |
108 | error = TMath::Sqrt(TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta)/ | |
109 | hEtaVsVtxCheck->GetBinContent(iVtx, iEta),2) + | |
75ec0f41 | 110 | TMath::Power(nEventsError/nEvents,2)); |
fcf2fb36 | 111 | hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value); |
112 | hEtaVsVtxCheck->SetBinError(iVtx, iEta, error); | |
75ec0f41 | 113 | } |
114 | } | |
115 | ||
116 | // then take the wieghted average for each eta | |
117 | // is this the right way to do it??? | |
fcf2fb36 | 118 | for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) { |
75ec0f41 | 119 | Float_t sum = 0; |
120 | Float_t sumError2 = 0; | |
121 | Int_t nMeasurements = 0; | |
122 | ||
fcf2fb36 | 123 | Float_t sumXw = 0; |
124 | Float_t sumW = 0; | |
75ec0f41 | 125 | |
7029240a | 126 | // do we have several histograms for different vertex positions? |
fcf2fb36 | 127 | Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1); |
7029240a | 128 | for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos) |
129 | { | |
fcf2fb36 | 130 | Int_t vertexBinBegin = 0; |
131 | Int_t vertexBinEnd = hVtx->GetNbinsX(); | |
132 | ||
133 | // the first histogram is always for the whole vertex range | |
134 | if (vertexPos > 0) | |
135 | { | |
136 | vertexBinBegin = vertexBinWidth * (vertexPos-1); | |
137 | vertexBinEnd = vertexBinBegin + vertexBinWidth; | |
138 | } | |
7029240a | 139 | |
fcf2fb36 | 140 | for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) { |
141 | if (hVtx->GetBinContent(iVtx)==0) continue; | |
a8d3b360 | 142 | if (hEtaVsVtxCheck->GetBinContent(iVtx, iEta)==0) continue; |
7029240a | 143 | |
fcf2fb36 | 144 | Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2); |
a8d3b360 | 145 | sumXw = sumXw + hEtaVsVtxCheck->GetBinContent(iVtx, iEta)*w; |
fcf2fb36 | 146 | sumW = sumW + w; |
7029240a | 147 | |
a8d3b360 | 148 | sum = sum + hEtaVsVtxCheck->GetBinContent(iVtx, iEta); |
149 | sumError2 = sumError2 + TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta),2); | |
7029240a | 150 | nMeasurements++; |
151 | } | |
152 | Float_t dndeta = 0; | |
153 | Float_t error = 0; | |
154 | ||
155 | if (nMeasurements!=0) { | |
156 | dndeta = sum/Float_t(nMeasurements); | |
157 | error = TMath::Sqrt(sumError2)/Float_t(nMeasurements); | |
158 | ||
fcf2fb36 | 159 | dndeta = sumXw/sumW; |
160 | error = 1/TMath::Sqrt(sumW); | |
7029240a | 161 | |
fcf2fb36 | 162 | dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta); |
163 | error = error/hdNdEta[vertexPos]->GetBinWidth(iEta); | |
7029240a | 164 | |
fcf2fb36 | 165 | hdNdEta[vertexPos]->SetBinContent(iEta, dndeta); |
166 | hdNdEta[vertexPos]->SetBinError(iEta, error); | |
7029240a | 167 | } |
75ec0f41 | 168 | } |
75ec0f41 | 169 | } |
170 | } | |
171 | ||
75ec0f41 | 172 | //____________________________________________________________________ |
173 | void | |
174 | dNdEtaAnalysis::SaveHistograms() { | |
175 | ||
7029240a | 176 | gDirectory->mkdir(GetName()); |
177 | gDirectory->cd(GetName()); | |
75ec0f41 | 178 | |
179 | hEtaVsVtx ->Write(); | |
fcf2fb36 | 180 | hEtaVsVtxUncorrected->Write(); |
75ec0f41 | 181 | hVtx ->Write(); |
7029240a | 182 | for (Int_t i=0; i<kVertexBinning; ++i) |
183 | hdNdEta[i] ->Write(); | |
75ec0f41 | 184 | |
185 | gDirectory->cd("../"); | |
186 | } | |
187 | ||
ceb5d1b5 | 188 | //____________________________________________________________________ |
189 | void dNdEtaAnalysis::DrawHistograms() | |
190 | { | |
a8d3b360 | 191 | TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800); |
192 | canvas->Divide(3, 2); | |
ceb5d1b5 | 193 | |
194 | canvas->cd(1); | |
195 | if (hEtaVsVtx) | |
196 | hEtaVsVtx->Draw("COLZ"); | |
197 | ||
198 | canvas->cd(2); | |
fcf2fb36 | 199 | if (hEtaVsVtxCheck) |
200 | hEtaVsVtxCheck->Draw("COLZ"); | |
201 | ||
202 | canvas->cd(3); | |
ceb5d1b5 | 203 | if (hEtaVsVtxUncorrected) |
204 | hEtaVsVtxUncorrected->Draw("COLZ"); | |
205 | ||
a8d3b360 | 206 | canvas->cd(4); |
207 | TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone"); | |
208 | clone->Divide(hEtaVsVtx); | |
209 | clone->Draw("COLZ"); | |
210 | ||
211 | canvas->cd(5); | |
ceb5d1b5 | 212 | if (hVtx) |
a8d3b360 | 213 | hVtx->Draw(); |
ceb5d1b5 | 214 | |
a8d3b360 | 215 | canvas->cd(6); |
7029240a | 216 | if (hdNdEta[0]) |
217 | hdNdEta[0]->Draw(); | |
fcf2fb36 | 218 | |
219 | // histograms for different vertices? | |
220 | if (kVertexBinning > 0) | |
221 | { | |
222 | // doesnt work, but i dont get it, giving up... | |
223 | /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450); | |
224 | //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2); | |
225 | //printf("%d\n", yPads); | |
226 | //canvas2->Divide(2, yPads); | |
227 | ||
228 | TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9); | |
229 | ||
230 | for (Int_t i=1; i<kVertexBinning; ++i) | |
231 | { | |
232 | //canvas2->cd(i-1); | |
233 | //printf("%d\n", i); | |
234 | if (hdNdEta[i]) | |
235 | { | |
236 | hdNdEta[i]->SetLineColor(i); | |
237 | hdNdEta[i]->Draw((i == 1) ? "" : "SAME"); | |
238 | legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1)); | |
239 | } | |
240 | } | |
241 | ||
242 | legend->Draw(); | |
243 | } | |
7029240a | 244 | } |
245 | ||
246 | Long64_t dNdEtaAnalysis::Merge(TCollection* list) | |
247 | { | |
248 | // Merges a list of dNdEtaAnalysis objects with this one. | |
249 | // This is needed for PROOF. | |
250 | // Returns the number of merged objects (including this) | |
251 | ||
252 | if (!list) | |
253 | return 0; | |
254 | ||
255 | if (list->IsEmpty()) | |
256 | return 1; | |
257 | ||
258 | TIterator* iter = list->MakeIterator(); | |
259 | TObject* obj; | |
260 | ||
261 | // sub collections | |
262 | const Int_t nCollections = kVertexBinning + 3; | |
263 | TList* collections[nCollections]; | |
264 | for (Int_t i=0; i<nCollections; ++i) | |
265 | collections[i] = new TList; | |
266 | ||
267 | Int_t count = 0; | |
268 | while ((obj = iter->Next())) | |
269 | { | |
270 | dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj); | |
271 | if (entry == 0) | |
272 | continue; | |
273 | ||
274 | collections[0]->Add(entry->hEtaVsVtx); | |
275 | collections[1]->Add(entry->hEtaVsVtxUncorrected); | |
276 | collections[2]->Add(entry->hVtx); | |
277 | ||
278 | for (Int_t i=0; i<kVertexBinning; ++i) | |
279 | collections[3+i]->Add(entry->hdNdEta[i]); | |
280 | ||
281 | ++count; | |
282 | } | |
283 | ||
284 | hEtaVsVtx->Merge(collections[0]); | |
285 | hEtaVsVtxUncorrected->Merge(collections[1]); | |
286 | hVtx->Merge(collections[2]); | |
287 | for (Int_t i=0; i<kVertexBinning; ++i) | |
288 | hdNdEta[i]->Merge(collections[3+i]); | |
289 | ||
290 | for (Int_t i=0; i<nCollections; ++i) | |
291 | delete collections[i]; | |
292 | ||
293 | return count+1; | |
ceb5d1b5 | 294 | } |