]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muondep/MuonTrackingEfficiency.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
remove 2 histograms and terminate method, adjust some histo ranges ranges, add more...
[u/mrichter/AliRoot.git] / PWG3 / muondep / MuonTrackingEfficiency.C
CommitLineData
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
23This 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
25Macro created by Lizardo Valencia Palomo
26
27lizardo.valencia.palomo@cern.ch
28
29*/
30
31
32
f785a99a 33TH1D* GraphToHist(Bool_t ToBePlotted, const TGraphAsymmErrors* Graph, Int_t Chamber);
1bc310ed 34std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors);
35std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chj, std::vector<Double_t> Chi, Int_t NumberOfDEs, Int_t SymmetricGroups);
36void 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);
37Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh);
38
39
f785a99a 40void 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
360std::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
375std::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
458Double_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
474void 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 588TH1D* 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}