]>
Commit | Line | Data |
---|---|---|
1bc310ed | 1 | #include <fstream> |
2 | #include <stdio.h> | |
3 | #include <stdlib.h> | |
4 | #include <TH1I.h> | |
5 | #include <TH2I.h> | |
6 | #include <TH1D.h> | |
7 | #include <TH2D.h> | |
8 | #include <TAxis.h> | |
9 | #include <TLegend.h> | |
10 | #include <TString.h> | |
11 | #include <Riostream.h> | |
12 | #include <TFile.h> | |
13 | #include <TList.h> | |
14 | #include <TCanvas.h> | |
15 | #include <TSystem.h> | |
16 | #include <TStyle.h> | |
17 | #include <TGraphAsymmErrors.h> | |
18 | #include <vector> | |
19 | #include <TGaxis.h> | |
20 | ||
21 | /* | |
22 | ||
23 | This macro is meant to be used with the output of the Efficiency Task. It computes the uncorrelated efficiencies per chammber, station and globally. It then looks for correlated dead arrays and corrects the previous calculated efficiency. | |
24 | ||
25 | Macro created by Lizardo Valencia Palomo | |
26 | ||
27 | lizardo.valencia.palomo@cern.ch | |
28 | ||
29 | */ | |
30 | ||
31 | ||
32 | ||
f785a99a | 33 | TH1D* GraphToHist(Bool_t ToBePlotted, const TGraphAsymmErrors* Graph, Int_t Chamber); |
1bc310ed | 34 | std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors); |
35 | std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chj, std::vector<Double_t> Chi, Int_t NumberOfDEs, Int_t SymmetricGroups); | |
36 | void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station); | |
37 | Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh); | |
38 | ||
39 | ||
f785a99a | 40 | void MuonTrackingEfficiency(TString fileName = "./AnalysisResults.root") |
1bc310ed | 41 | { |
42 | ||
43 | TFile *file = new TFile(fileName, "read"); | |
44 | TDirectoryFile *Directory = (TDirectoryFile*)((file->GetDirectory("MUON_Efficiency"))); | |
45 | ||
46 | std::vector<TH1D*> CorrectedHistos; | |
47 | std::vector<TH1D*> CorrectedChiHisto; | |
48 | std::vector<TH1D*> CorrectedChjHisto; | |
49 | ||
50 | TList *listTT = (TList *)Directory->Get("TotalTracksPerChamber;1"); | |
51 | TList *listTD = (TList *)Directory->Get("TracksDetectedPerChamber;1"); | |
f785a99a | 52 | TList *listSD = (TList *)Directory->Get("SingleDetectedPerChamber;1"); |
1bc310ed | 53 | |
54 | TH1D *EfficiencyHisto = 0x0; | |
55 | TH1D *TD = 0x0; | |
56 | TH1D *TT = 0x0; | |
57 | ||
58 | for (Int_t Chamber = 1; Chamber < 12; Chamber++) //Loop over all chambers + 1 | |
59 | { | |
60 | ||
61 | TH2F *TD2Dim = (TH2F *) listTD->At(Chamber - 1); | |
62 | TH2F *TT2Dim = (TH2F *) listTT->At(Chamber - 1); | |
63 | TH2F *SD2Dim = (TH2F *) listSD->At(Chamber - 1); | |
64 | ||
65 | Int_t nDEch = 0; | |
66 | nDEch = TD2Dim->GetNbinsX(); | |
67 | ||
68 | TH1D *Ncluster = 0x0; | |
69 | TH1D *Corrected = 0x0; | |
70 | TH1D *SD = 0x0; | |
71 | ||
72 | TH1D *EfficiencyHistoCh11 = 0x0; | |
73 | TH1D *TDch11 = 0x0; | |
74 | TH1D *TTch11 = 0x0; | |
75 | ||
76 | if (Chamber == 11) | |
77 | { | |
78 | ||
79 | TDch11 = new TH1D(Form("TDch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5); | |
80 | TTch11 = new TH1D(Form("TTch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5); | |
81 | ||
82 | } | |
83 | ||
84 | TD = new TH1D(Form("TD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5); | |
85 | TT = new TH1D(Form("TT %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5); | |
86 | SD = new TH1D(Form("SD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5); | |
87 | ||
88 | std::vector<Double_t> TDentries (nDEch, 0.0); | |
89 | std::vector<Double_t> TTentries (nDEch, 0.0); | |
90 | std::vector<Double_t> SDentries (nDEch, 0.0); | |
91 | ||
92 | for (Int_t ix = 1; ix <= nDEch; ix++) | |
93 | { | |
94 | for (Int_t iy = 1; iy <= TD2Dim->GetNbinsY(); iy++) | |
95 | { | |
96 | ||
97 | TDentries[ix-1] += TD2Dim->GetBinContent(ix,iy); | |
98 | TTentries[ix-1] += TT2Dim->GetBinContent(ix,iy); | |
99 | SDentries[ix-1] += SD2Dim->GetBinContent(ix,iy); | |
100 | ||
101 | } | |
102 | ||
103 | TD->SetBinContent(ix,TDentries[ix-1]); | |
104 | TD->SetBinError(ix,sqrt(TDentries[ix-1])); | |
105 | TT->SetBinContent(ix,TTentries[ix-1]); | |
106 | TT->SetBinError(ix,sqrt(TTentries[ix-1])); | |
107 | SD->SetBinContent(ix,SDentries[ix-1]); | |
108 | SD->SetBinError(ix,sqrt(SDentries[ix-1])); | |
109 | ||
110 | if (Chamber != 11) continue; | |
111 | ||
112 | TDch11->SetBinContent(ix,TDentries[ix-1]); | |
113 | TDch11->SetBinError(ix,sqrt(TDentries[ix-1])); | |
114 | TTch11->SetBinContent(ix,TTentries[ix-1]); | |
115 | TTch11->SetBinError(ix,sqrt(TTentries[ix-1])); | |
116 | ||
117 | } | |
118 | ||
119 | TGraphAsymmErrors *EfficiencyGraph = new TGraphAsymmErrors(TD, TT); | |
120 | ||
f785a99a | 121 | EfficiencyHisto = GraphToHist(kFALSE,EfficiencyGraph,Chamber); |
1bc310ed | 122 | |
123 | Ncluster = new TH1D(Form("Ncluster %d",Chamber), "N_{Cluster}", nDEch, - 0.5, nDEch - 0.5); | |
124 | Ncluster->Sumw2(); | |
125 | Ncluster->Add(TD,SD,1,1); | |
126 | ||
127 | Corrected = new TH1D(Form("Corrected %d",Chamber), Form("Corrected N_{Cluster} Ch %d",Chamber), nDEch, - 0.5, nDEch - 0.5); | |
128 | Corrected->Sumw2(); | |
129 | Corrected->Divide(Ncluster,EfficiencyHisto,1,1); | |
130 | ||
131 | CorrectedHistos.push_back(Corrected); | |
132 | ||
133 | if (Chamber != 11) continue; | |
134 | ||
135 | TGraphAsymmErrors *EfficiencyGraphCh11 = new TGraphAsymmErrors(TDch11, TTch11); | |
136 | ||
f785a99a | 137 | EfficiencyHistoCh11 = GraphToHist(kTRUE,EfficiencyGraphCh11,Chamber); |
1bc310ed | 138 | |
139 | TCanvas *CanvasEffPerCh = new TCanvas("Efficiency per chambers", "Efficiency per chambers"); | |
140 | gStyle->SetOptStat(0); | |
141 | CanvasEffPerCh->SetFillColor(0); | |
142 | CanvasEffPerCh->SetBorderSize(1); | |
143 | CanvasEffPerCh->SetFrameLineWidth(2); | |
144 | ||
145 | EfficiencyHistoCh11->SetTitle(kFALSE); | |
146 | EfficiencyHistoCh11->SetMarkerStyle(20); | |
147 | EfficiencyHistoCh11->SetMarkerColor(1); | |
148 | EfficiencyHistoCh11->SetMarkerSize(1.3); | |
149 | EfficiencyHistoCh11->GetXaxis()->SetTitle("Chamber"); | |
150 | EfficiencyHistoCh11->GetXaxis()->SetNdivisions(Chamber); | |
151 | EfficiencyHistoCh11->GetXaxis()->SetTitleSize(0.05); | |
152 | EfficiencyHistoCh11->GetXaxis()->SetTitleOffset(0.91); | |
153 | EfficiencyHistoCh11->GetXaxis()->SetLabelSize(0.05); | |
154 | EfficiencyHistoCh11->GetYaxis()->SetTitle("Efficiency"); | |
155 | EfficiencyHistoCh11->GetYaxis()->CenterTitle(kTRUE); | |
156 | EfficiencyHistoCh11->GetYaxis()->SetTitleSize(0.06); | |
157 | EfficiencyHistoCh11->GetYaxis()->SetTitleOffset(0.79); | |
158 | EfficiencyHistoCh11->GetYaxis()->SetLabelSize(0.05); | |
159 | EfficiencyHistoCh11->Draw("e"); | |
160 | ||
161 | } | |
162 | ||
163 | /*************************Correlated Dead Arrays check****************************/ | |
164 | ||
165 | ||
166 | Double_t MissingTracks [5] = {0.0};//Station one is in slot cero!!! | |
167 | ||
168 | for (Int_t CDAperSt = 0; CDAperSt < 5; CDAperSt++) | |
169 | { | |
170 | ||
171 | std::vector<Double_t> FilledVectorChi = FillVector(CorrectedHistos[CDAperSt*2], kFALSE); | |
172 | std::vector<Double_t> FilledVectorChj = FillVector(CorrectedHistos[(CDAperSt*2)+1], kFALSE); | |
173 | std::vector<Double_t> FilledVectorChiErr = FillVector(CorrectedHistos[CDAperSt*2], kTRUE); | |
174 | std::vector<Double_t> FilledVectorChjErr = FillVector(CorrectedHistos[(CDAperSt*2)+1], kTRUE); | |
175 | ||
176 | Int_t DEs = FilledVectorChi.size(); | |
177 | Int_t SimGroups = 0; | |
178 | Double_t ToFillMissingTracks = 0.0; | |
179 | ||
180 | if ( CDAperSt < 2 ) SimGroups = static_cast<Int_t>(DEs/2); | |
181 | else SimGroups = static_cast<Int_t>(DEs/2) + 1; | |
182 | ||
183 | std::vector<Double_t> CDAs = CorrelatedDeadArray(FilledVectorChi, FilledVectorChj, DEs, SimGroups); | |
184 | Drawing(FilledVectorChi, FilledVectorChj, FilledVectorChiErr, FilledVectorChjErr, DEs, SimGroups, CDAperSt); | |
185 | ||
186 | for (Int_t DEst = 0; DEst < DEs; DEst++) | |
187 | { | |
188 | if ( CDAs[DEst] != 0.0 ) | |
189 | { | |
190 | ||
191 | printf ("Looks like there is a CDA in DE %d of Station %d: N00 = %.2f \n", DEst, CDAperSt+1, CDAs[DEst]); | |
192 | ToFillMissingTracks += CDAs[DEst]; | |
193 | ||
194 | } | |
195 | } | |
196 | ||
197 | MissingTracks[CDAperSt] = ToFillMissingTracks; | |
198 | ||
199 | } | |
200 | ||
201 | ||
202 | ||
203 | /**********************************Efficiency from Efficiency Task************************************/ | |
204 | ||
205 | ||
206 | vector<Double_t> chambersEff; | |
207 | vector<Double_t> chambersEffErr; | |
208 | vector<Double_t> StEfficiency; | |
209 | vector<Double_t> StEfficiencyErr; | |
210 | ||
211 | chambersEff.push_back(0.0); | |
212 | chambersEffErr.push_back(0.0); | |
213 | StEfficiency.push_back(0.0); | |
214 | StEfficiencyErr.push_back(0.0); | |
215 | ||
216 | for (Int_t EfCh = 0; EfCh < 10; EfCh++) | |
217 | { | |
218 | ||
219 | chambersEff.push_back(EfficiencyHisto->GetBinContent(EfCh+1)); | |
220 | chambersEffErr.push_back(EfficiencyHisto->GetBinError(EfCh+1)); | |
221 | ||
222 | cout << "Efficiency chamber " << EfCh+1 << " = " << chambersEff[EfCh+1] << " +- " << chambersEffErr[EfCh+1] << endl; | |
223 | ||
224 | } | |
225 | ||
226 | for (Int_t EffSt = 0; EffSt < 4; EffSt++) | |
227 | { | |
228 | ||
229 | if (EffSt < 3) | |
230 | { | |
231 | ||
232 | StEfficiency.push_back(1.0 - (1.0-chambersEff[(2*EffSt)+1])*(1.0-chambersEff[(2*EffSt)+2])); | |
233 | StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1] + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2])); | |
234 | ||
235 | cout << "Efficiency station " << EffSt+1 <<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl; | |
236 | ||
237 | } | |
238 | ||
239 | else | |
240 | { | |
241 | ||
242 | StEfficiency.push_back(chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4] | |
243 | + (1.0 - chambersEff[(2*EffSt)+1]) * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4] | |
244 | + chambersEff[(2*EffSt)+1] * (1.0 - chambersEff[(2*EffSt)+2]) * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4] | |
245 | + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * (1.0 - chambersEff[(2*EffSt)+3]) * chambersEff[(2*EffSt)+4] | |
246 | + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * (1.0 - chambersEff[(2*EffSt)+4])); | |
247 | ||
248 | StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1] | |
249 | + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2] + chambersEffErr[(2*EffSt)+3]/chambersEff[(2*EffSt)+3] | |
250 | + chambersEffErr[(2*EffSt)+4]/chambersEff[(2*EffSt)+4])); | |
251 | ||
252 | cout << "Efficiency station 45"<<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl; | |
253 | //Notice efficiency from st 45 is in slot 4 of the vector!!! | |
254 | ||
255 | } | |
256 | ||
257 | } | |
258 | ||
259 | StEfficiency.push_back(StEfficiency[1]*StEfficiency[2]*StEfficiency[3]*StEfficiency[4]); | |
260 | StEfficiencyErr.push_back(StEfficiencyErr[1]+StEfficiencyErr[2]+StEfficiencyErr[3]+StEfficiencyErr[4]); | |
261 | ||
262 | cout << "Total efficiency = " <<StEfficiency[5]<< " +- " <<StEfficiencyErr[5]<< endl; | |
263 | ||
264 | ||
265 | ||
266 | /*************************************Correcting efficiencies********************************************/ | |
267 | ||
268 | ||
269 | TH1D *TTcorrected = new TH1D("TTcorrected", "TDcorrected", 10, -0.5, 9.5); | |
270 | ||
271 | for (Int_t St = 0; St < 5; St++) | |
272 | { | |
273 | ||
274 | Double_t NewBinContentChi = TT->GetBinContent((2*St)+1) + MissingTracks[St]; | |
275 | TTcorrected->SetBinContent((2*St)+1,NewBinContentChi); | |
276 | TTcorrected->SetBinContent((2*St)+2,TT->GetBinContent((2*St)+2)); | |
277 | ||
278 | } | |
279 | ||
280 | TGraphAsymmErrors *EfficiencyGraphCorrected = new TGraphAsymmErrors(TD,TTcorrected); | |
281 | ||
f785a99a | 282 | TH1D *EfficiencyHistoCorrected = GraphToHist(1,EfficiencyGraphCorrected,11);//In this case doesnt matter which value is used in the first slot. |
1bc310ed | 283 | |
284 | vector<Double_t> chambersEffCorrected; | |
285 | vector<Double_t> chambersEffErrCorrected; | |
286 | ||
287 | chambersEffCorrected.push_back(0.0); | |
288 | chambersEffErrCorrected.push_back(0.0); | |
289 | ||
290 | for (Int_t EfChCorr = 0; EfChCorr < 10; EfChCorr++) | |
291 | { | |
292 | ||
293 | chambersEffCorrected.push_back(EfficiencyHistoCorrected->GetBinContent(EfChCorr+1)); | |
294 | chambersEffErrCorrected.push_back(EfficiencyHistoCorrected->GetBinError(EfChCorr+1)); | |
295 | ||
296 | if ( (chambersEffCorrected[EfChCorr+1]/chambersEff[EfChCorr+1]) == 1 ) continue; | |
297 | ||
298 | cout << "Corrected efficiency chamber " << EfChCorr+1 << " = " << chambersEffCorrected[EfChCorr+1] << " +- " << chambersEffErrCorrected[EfChCorr+1] << endl; | |
299 | ||
300 | } | |
301 | ||
302 | vector<Double_t> StEfficiencyCorrected; | |
303 | vector<Double_t> StEfficiencyErrCorrected; | |
304 | ||
305 | StEfficiencyCorrected.push_back(0.0); | |
306 | StEfficiencyErrCorrected.push_back(0.0); | |
307 | ||
308 | //By default CDA in Station 4 is not corrected. | |
309 | ||
310 | for (Int_t Station = 0; Station < 4; Station++) | |
311 | { | |
312 | ||
313 | if (MissingTracks[Station] == 0.0) | |
314 | { | |
315 | ||
316 | StEfficiencyCorrected.push_back(StEfficiency[Station+1]); | |
317 | StEfficiencyErrCorrected.push_back(StEfficiencyErr[Station+1]); | |
318 | ||
319 | } | |
320 | ||
321 | else | |
322 | { | |
323 | ||
324 | if (Station < 3) | |
325 | { | |
326 | ||
327 | StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+1]/chambersEff[(2*Station)+1])); | |
328 | StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1])*(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1]) + (chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1])*(chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1]))); | |
329 | ||
330 | cout<<"Corrected efficiency station "<<Station+1<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl; | |
331 | ||
332 | } | |
333 | ||
334 | if (Station == 3) | |
335 | { | |
336 | ||
337 | StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+3]/chambersEff[(2*Station)+3])); | |
338 | StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3])*(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3]) + (chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3])*(chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3]))); | |
339 | ||
340 | cout<<"Corrected efficiency station 45"<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl; | |
341 | ||
342 | } | |
343 | ||
344 | ||
345 | } | |
346 | ||
347 | } | |
348 | ||
349 | StEfficiencyCorrected.push_back(StEfficiencyCorrected[1]*StEfficiencyCorrected[2]*StEfficiencyCorrected[3]*StEfficiencyCorrected[4]); | |
350 | StEfficiencyErrCorrected.push_back(StEfficiencyErrCorrected[1]+StEfficiencyErrCorrected[2]+StEfficiencyErrCorrected[3]+StEfficiencyErrCorrected[4]); | |
351 | ||
352 | cout<<"Corrected total efficiency"<<" = "<<StEfficiencyCorrected[5]<< " +- " <<StEfficiencyErrCorrected[5]<< endl; | |
353 | ||
354 | } | |
355 | ||
356 | ||
357 | /*************************************Functions************************************/ | |
358 | ||
359 | ||
360 | std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors) | |
361 | { | |
362 | ||
363 | Int_t NumOfDEs = Chamber->GetNbinsX(); | |
364 | std::vector<Double_t> ChVector (NumOfDEs, 0); | |
365 | for (Int_t i = 1; i <= NumOfDEs; i++) | |
366 | { | |
367 | if (Errors == kFALSE) ChVector.at(i-1) = Chamber->GetBinContent(i); | |
368 | if (Errors == kTRUE) ChVector.at(i-1) = Chamber->GetBinError(i); | |
369 | } | |
370 | ||
371 | return ChVector; | |
372 | ||
373 | } | |
374 | ||
375 | std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chi, std::vector<Double_t> Chj, Int_t NumberOfDEs, Int_t SymmetricGroups) | |
376 | { | |
377 | ||
378 | std::vector<Double_t> Average (SymmetricGroups, 0.0); | |
379 | std::vector<Double_t> CDA (NumberOfDEs, 0.0); | |
380 | ||
381 | for (Int_t Sim = 0; Sim < SymmetricGroups/2; Sim++) | |
382 | { | |
383 | if (Sim == 0) | |
384 | { | |
385 | ||
386 | if ( (Sim == 0) && (SymmetricGroups/2 == 1) ) | |
387 | { | |
388 | ||
389 | Average[Sim] = (Chi[Sim]+Chi[Sim+1]+Chi[Sim+2]+Chi[Sim+3])/4; | |
390 | Average[Sim+(SymmetricGroups/2)] = (Chj[Sim]+Chj[Sim+1]+Chj[Sim+2]+Chj[Sim+3])/4; | |
391 | ||
392 | for (Int_t i = 1; i < 5; i++) | |
393 | { | |
394 | ||
395 | if ( (Chi[i-1] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[i-1] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
396 | { | |
397 | CDA[i-1] = N00(Average[Sim], Chi[i-1], Sim, SymmetricGroups/2); | |
398 | } | |
399 | ||
400 | } | |
401 | ||
402 | } | |
403 | ||
404 | else | |
405 | { | |
406 | ||
407 | Average[Sim] = (Chi[Sim] + Chi[Sim+(NumberOfDEs/2)])/2; | |
408 | Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[Sim+(NumberOfDEs/2)])/2; | |
409 | ||
410 | ||
411 | if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
412 | { | |
413 | CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2); | |
414 | } | |
415 | if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
416 | { | |
417 | CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2); | |
418 | } | |
419 | ||
420 | } | |
421 | ||
422 | } | |
423 | ||
424 | else | |
425 | { | |
426 | ||
427 | Average[Sim] = (Chi[Sim] + Chi[(NumberOfDEs/2)-Sim] + Chi[Sim+(NumberOfDEs/2)] + Chi[NumberOfDEs-Sim])/4; | |
428 | Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[(NumberOfDEs/2)-Sim] + Chj[Sim+(NumberOfDEs/2)] + Chj[NumberOfDEs-Sim])/4; | |
429 | ||
430 | if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
431 | { | |
432 | CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2); | |
433 | } | |
434 | if ( (Chi[(NumberOfDEs/2)-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[(NumberOfDEs/2)-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
435 | { | |
436 | CDA[(NumberOfDEs/2)-Sim] = N00(Average[Sim], Chi[(NumberOfDEs/2)-Sim], Sim, SymmetricGroups/2); | |
437 | } | |
438 | if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
439 | { | |
440 | CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2); | |
441 | } | |
442 | ||
443 | if ( (Chi[NumberOfDEs-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[NumberOfDEs-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) ) | |
444 | { | |
445 | CDA[NumberOfDEs-Sim] = N00(Average[Sim], Chi[NumberOfDEs-Sim], Sim, SymmetricGroups/2); | |
446 | ||
447 | } | |
448 | ||
449 | } | |
450 | ||
451 | } | |
452 | ||
453 | return CDA; | |
454 | ||
455 | } | |
456 | ||
457 | ||
458 | Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh) | |
459 | { | |
460 | ||
461 | Double_t RecomputedAverage = 0.0; | |
462 | Double_t ToBeAdded = 0.0; | |
463 | ||
464 | if ( (SimGroup == 0) && (SimGroupsPerCh != 1) ) RecomputedAverage = (2*Average) - CDA; | |
465 | else RecomputedAverage = ((4*Average) - CDA)/3; | |
466 | ||
467 | ToBeAdded = RecomputedAverage - CDA; | |
468 | ||
469 | return ToBeAdded; | |
470 | ||
471 | } | |
472 | ||
473 | ||
474 | void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station) | |
475 | { | |
476 | ||
477 | TCanvas *Canvas = new TCanvas(Form("Station %d",Station+1), Form("Station %d",Station+1)); | |
478 | gStyle->SetOptStat(0); | |
479 | Canvas->SetFillColor(0); | |
480 | if ( (Station+1) < 3) Canvas->Divide(SymGroups/2); | |
481 | if ( (Station+1) == 3) Canvas->Divide(Station+1,Station); | |
482 | if ( (Station+1) > 3) Canvas->Divide(3,3); | |
483 | Canvas->SetBorderSize(1); | |
484 | Canvas->SetFrameLineWidth(2); | |
485 | ||
486 | for (Int_t Sim = 0; Sim < SymGroups/2; Sim++) | |
487 | { | |
488 | ||
489 | TH1D *Chamberi = new TH1D(Form("Chi SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5); | |
490 | TGaxis::SetMaxDigits(2); | |
491 | Chamberi->SetTitle(kFALSE); | |
492 | Chamberi->SetMarkerStyle(20); | |
493 | Chamberi->SetMarkerColor(1); | |
494 | Chamberi->SetMarkerSize(1.3); | |
495 | Chamberi->GetXaxis()->SetTitle("Detection Element"); | |
496 | Chamberi->GetXaxis()->SetNdivisions(Bins+1); | |
497 | Chamberi->GetXaxis()->SetTitleSize(0.05); | |
498 | Chamberi->GetXaxis()->SetTitleOffset(0.91); | |
499 | Chamberi->GetXaxis()->SetLabelSize(0.05); | |
500 | Chamberi->GetYaxis()->SetTitle("Corrected N_{Cluster}"); | |
501 | Chamberi->GetYaxis()->CenterTitle(kTRUE); | |
502 | Chamberi->GetYaxis()->SetTitleSize(0.06); | |
503 | Chamberi->GetYaxis()->SetTitleOffset(0.79); | |
504 | Chamberi->GetYaxis()->SetLabelSize(0.05); | |
505 | TH1D *Chamberj = new TH1D(Form("Chj SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5); | |
506 | Chamberj->SetTitle(kFALSE); | |
507 | Chamberj->SetMarkerStyle(20); | |
508 | Chamberj->SetMarkerColor(2); | |
509 | Chamberj->SetMarkerSize(1.3); | |
510 | ||
511 | if (Sim == 0) | |
512 | { | |
513 | ||
514 | if ( (Sim == 0) && (SymGroups/2 == 1) ) | |
515 | { | |
516 | ||
517 | for (Int_t i = 1; i < 5; i++) | |
518 | { | |
519 | ||
520 | Chamberi->SetBinContent(i, Chi[i-1]); | |
521 | Chamberi->SetBinError(i, ChiErr[i-1]); | |
522 | Chamberj->SetBinContent(i, Chj[i-1]); | |
523 | Chamberj->SetBinError(i, ChjErr[i-1]); | |
524 | ||
525 | } | |
526 | ||
527 | } | |
528 | ||
529 | else | |
530 | { | |
531 | ||
532 | Chamberi->SetBinContent(Sim+1, Chi[Sim]); | |
533 | Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]); | |
534 | Chamberi->SetBinError(Sim+1, ChiErr[Sim]); | |
535 | Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]); | |
536 | Chamberj->SetBinContent(Sim+1, Chj[Sim]); | |
537 | Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]); | |
538 | Chamberj->SetBinError(Sim+1, ChjErr[Sim]); | |
539 | Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]); | |
540 | ||
541 | } | |
542 | ||
543 | } | |
544 | ||
545 | else | |
546 | { | |
547 | ||
548 | Chamberi->SetBinContent(Sim+1, Chi[Sim]); | |
549 | Chamberi->SetBinContent((Bins/2)-(Sim-1), Chi[(Bins/2)-Sim]); | |
550 | Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]); | |
551 | Chamberi->SetBinContent(Bins-(Sim-1), Chi[Bins-Sim]); | |
552 | ||
553 | Chamberi->SetBinError(Sim+1, ChiErr[Sim]); | |
554 | Chamberi->SetBinError((Bins/2)-(Sim-1), ChiErr[(Bins/2)-Sim]); | |
555 | Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]); | |
556 | Chamberi->SetBinError(Bins-(Sim-1), ChiErr[Bins-Sim]); | |
557 | ||
558 | Chamberj->SetBinContent(Sim+1, Chj[Sim]); | |
559 | Chamberj->SetBinContent((Bins/2)-(Sim-1), Chj[(Bins/2)-Sim]); | |
560 | Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]); | |
561 | Chamberj->SetBinContent(Bins-(Sim-1), Chj[Bins-Sim]); | |
562 | ||
563 | Chamberj->SetBinError(Sim+1, ChjErr[Sim]); | |
564 | Chamberj->SetBinError((Bins/2)-(Sim-1), ChjErr[(Bins/2)-Sim]); | |
565 | Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]); | |
566 | Chamberj->SetBinError(Bins-(Sim-1), ChjErr[Bins-Sim]); | |
567 | ||
568 | } | |
569 | ||
570 | Canvas->cd(Sim+1); | |
571 | Chamberi->Draw("e"); | |
572 | Chamberj->Draw("esame"); | |
573 | ||
574 | TLegend *Legend = new TLegend (0.75, 0.8, 1.0, 1.0); | |
575 | Legend->SetBorderSize(0); | |
576 | Legend->SetFillColor(0); | |
577 | if ( (Station+1) < 3 ) Legend->SetTextSize(0.05); | |
578 | if ( (Station+1) > 2 ) Legend->SetTextSize(0.05); | |
579 | Legend->AddEntry(Chamberi,Form("Chamber %d",2*Station+1),"p"); | |
580 | Legend->AddEntry(Chamberj,Form("Chamber %d",2*Station+2),"p"); | |
581 | Legend->Draw(); | |
582 | ||
583 | } | |
584 | ||
585 | } | |
586 | ||
587 | ||
f785a99a | 588 | TH1D* GraphToHist(Bool_t ToBePlotted,const TGraphAsymmErrors* Graph,Int_t Chamber) |
1bc310ed | 589 | { |
590 | ||
f785a99a | 591 | Int_t nBins = 0; |
592 | if ( Chamber < 5 ) nBins = 4; | |
593 | if ( (Chamber == 5) || (Chamber == 6) ) nBins = 18; | |
594 | if ( (Chamber > 6) && (Chamber < 11) ) nBins = 26; | |
595 | if ( Chamber == 11 ) nBins = 10; | |
1bc310ed | 596 | TArrayF bins(nBins+1); |
597 | TArrayF y(nBins); | |
598 | TArrayF ey(nBins); | |
599 | Double_t dx = 1; | |
600 | Double_t xmin = 10000; | |
601 | Double_t xmax = -10000; | |
602 | for (Int_t i = 0; i < nBins; i++) { | |
603 | Double_t x = Graph->GetX()[i]; | |
604 | Double_t exl = Graph->GetEXlow()[i]; | |
605 | Double_t exh = Graph->GetEXhigh()[i]; | |
606 | xmin = TMath::Min(x-exl, xmin); | |
607 | xmax = TMath::Max(x+exh, xmax); | |
608 | bins.fArray[i] = x-exl; | |
609 | bins.fArray[i+1] = x+exh; | |
610 | Double_t dxi = exh+exl; | |
611 | if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1]; | |
612 | if (dx == 0) dx = dxi; | |
613 | else if (dxi != dx) dx = 0; | |
614 | ||
615 | y.fArray[i] = Graph->GetY()[i]; | |
616 | ey.fArray[i] = TMath::Max(Graph->GetEYlow()[i],Graph->GetEYhigh()[i]); | |
617 | ||
618 | } | |
619 | TString name(Graph->GetName()); | |
620 | TString title(Graph->GetTitle()); | |
621 | TH1D* h = 0; | |
622 | if (ToBePlotted == kTRUE) { | |
623 | h = new TH1D(name.Data(), title.Data(), nBins, bins[0]-dx/2 + 1.5, bins[nBins]-dx/2 + 1.5); | |
624 | } | |
625 | else { | |
f785a99a | 626 | if ( Chamber < 5 ) h = new TH1D(name.Data(), title.Data(), 4, -0.5, 3.5); |
627 | if ( (Chamber == 5) || (Chamber == 6) ) h = new TH1D(name.Data(), title.Data(), 18, -0.5, 17.5); | |
628 | if ( (Chamber > 6) && (Chamber < 11) ) h = new TH1D(name.Data(), title.Data(), 26, -0.5, 25.5); | |
629 | if ( Chamber == 11 ) h = new TH1D(name.Data(), title.Data(), 10, -0.5, 9.5); | |
1bc310ed | 630 | } |
631 | for (Int_t i = 1; i <= nBins; i++) { | |
632 | h->SetBinContent(i, y.fArray[i-1]); | |
633 | h->SetBinError(i, ey.fArray[i-1]); | |
634 | } | |
635 | h->SetMarkerStyle(Graph->GetMarkerStyle()); | |
636 | h->SetMarkerColor(Graph->GetMarkerColor()); | |
637 | h->SetMarkerSize(Graph->GetMarkerSize()); | |
638 | h->SetDirectory(0); | |
639 | ||
640 | return h; | |
641 | } |