]>
Commit | Line | Data |
---|---|---|
75ec0f41 | 1 | #include "dNdEtaAnalysis.h" |
2 | ||
ceb5d1b5 | 3 | #include <TFile.h> |
4 | #include <TH2F.h> | |
5 | #include <TH1D.h> | |
6 | #include <TMath.h> | |
7 | #include <TCanvas.h> | |
8 | ||
75ec0f41 | 9 | //____________________________________________________________________ |
b7f4a1fd | 10 | ClassImp(dNdEtaAnalysis) |
75ec0f41 | 11 | |
12 | //____________________________________________________________________ | |
13 | dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name) { | |
14 | ||
15 | fName = TString(name); | |
16 | ||
ceb5d1b5 | 17 | hEtaVsVtx = new TH2F("eta_vs_vtx","",80,-20,20,120,-6,6); |
75ec0f41 | 18 | hEtaVsVtx->SetXTitle("vtx z [cm]"); |
19 | hEtaVsVtx->SetYTitle("#eta"); | |
4dd2ad81 | 20 | |
21 | hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone("eta_vs_vtx_uncorrected")); | |
22 | hVtx = hEtaVsVtx->ProjectionX("vtx"); | |
23 | hdNdEta = hEtaVsVtx->ProjectionY("dNdEta"); | |
24 | ||
75ec0f41 | 25 | hdNdEta->SetYTitle("dN/d#eta"); |
26 | ||
27 | hEtaVsVtx->Sumw2(); | |
28 | hVtx->Sumw2(); | |
75ec0f41 | 29 | } |
30 | ||
31 | //____________________________________________________________________ | |
32 | void | |
33 | dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) { | |
34 | hEtaVsVtx->Fill(vtx, eta, weight); | |
ceb5d1b5 | 35 | hEtaVsVtxUncorrected->Fill(vtx,eta); |
75ec0f41 | 36 | } |
37 | ||
38 | //____________________________________________________________________ | |
39 | void | |
40 | dNdEtaAnalysis::FillEvent(Float_t vtx) { | |
41 | hVtx->Fill(vtx); | |
42 | } | |
43 | ||
44 | //____________________________________________________________________ | |
45 | void | |
46 | dNdEtaAnalysis::Finish() { | |
47 | ||
48 | // first normalize with n events (per vtx) | |
49 | for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) { | |
50 | Float_t nEvents = hVtx->GetBinContent(i_vtx); | |
51 | Float_t nEventsError = hVtx->GetBinError(i_vtx); | |
52 | ||
53 | if (nEvents==0) continue; | |
54 | ||
55 | for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) { | |
56 | Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents; | |
57 | if (value==0) continue; | |
58 | Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents; | |
59 | error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/ | |
60 | hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) + | |
61 | TMath::Power(nEventsError/nEvents,2)); | |
62 | hEtaVsVtx->SetBinContent(i_vtx, i_eta, value); | |
63 | hEtaVsVtx->SetBinError(i_vtx, i_eta, error); | |
64 | } | |
65 | } | |
66 | ||
67 | // then take the wieghted average for each eta | |
68 | // is this the right way to do it??? | |
69 | for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) { | |
70 | Float_t sum = 0; | |
71 | Float_t sumError2 = 0; | |
72 | Int_t nMeasurements = 0; | |
73 | ||
74 | Float_t sum_xw = 0; | |
75 | Float_t sum_w = 0; | |
76 | ||
77 | for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) { | |
78 | if (hVtx->GetBinContent(i_vtx)==0) continue; | |
79 | if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue; | |
80 | ||
81 | Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2); | |
82 | sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w; | |
83 | sum_w = sum_w + w; | |
84 | ||
85 | sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta); | |
86 | sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2); | |
87 | nMeasurements++; | |
88 | } | |
89 | Float_t dndeta = 0; | |
90 | Float_t error = 0; | |
91 | ||
92 | if (nMeasurements!=0) { | |
93 | dndeta = sum/Float_t(nMeasurements); | |
94 | error = TMath::Sqrt(sumError2)/Float_t(nMeasurements); | |
95 | ||
96 | dndeta = sum_xw/sum_w; | |
97 | error = 1/TMath::Sqrt(sum_w); | |
98 | ||
99 | dndeta = dndeta/hdNdEta->GetBinWidth(i_eta); | |
100 | error = error/hdNdEta->GetBinWidth(i_eta); | |
101 | ||
ceb5d1b5 | 102 | hdNdEta->SetBinContent(i_eta, dndeta); |
103 | hdNdEta->SetBinError(i_eta, error); | |
75ec0f41 | 104 | } |
105 | ||
106 | } | |
107 | } | |
108 | ||
109 | ||
110 | ||
111 | //____________________________________________________________________ | |
112 | void | |
113 | dNdEtaAnalysis::SaveHistograms() { | |
114 | ||
115 | gDirectory->mkdir(fName.Data()); | |
116 | gDirectory->cd(fName.Data()); | |
117 | ||
118 | hEtaVsVtx ->Write(); | |
119 | hVtx ->Write(); | |
120 | hdNdEta ->Write(); | |
121 | ||
122 | gDirectory->cd("../"); | |
123 | } | |
124 | ||
ceb5d1b5 | 125 | //____________________________________________________________________ |
126 | void dNdEtaAnalysis::DrawHistograms() | |
127 | { | |
128 | TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800); | |
129 | canvas->Divide(2, 2); | |
130 | ||
131 | canvas->cd(1); | |
132 | if (hEtaVsVtx) | |
133 | hEtaVsVtx->Draw("COLZ"); | |
134 | ||
135 | canvas->cd(2); | |
136 | if (hEtaVsVtxUncorrected) | |
137 | hEtaVsVtxUncorrected->Draw("COLZ"); | |
138 | ||
139 | canvas->cd(3); | |
140 | if (hVtx) | |
141 | hVtx->Draw(); | |
142 | ||
143 | canvas->cd(4); | |
144 | if (hdNdEta) | |
145 | hdNdEta->Draw(); | |
146 | } |