]>
Commit | Line | Data |
---|---|---|
d4945a3b | 1 | #ifndef __CINT__ |
2 | #include <TTree.h> | |
3 | #include <TFile.h> | |
4 | #include <TClonesArray.h> | |
5 | #include <TStyle.h> | |
6 | #include <TH1F.h> | |
7 | #include <TH2F.h> | |
8 | #include <TProfile.h> | |
9 | #include <TGraph.h> | |
10 | #include <TCanvas.h> | |
11 | #include <TLegend.h> | |
12 | #include "AliCDBManager.h" | |
13 | #include "AliCDBEntry.h" | |
14 | #include "AliRawReader.h" | |
15 | #include <AliTPCRawStreamV3.h> | |
16 | #include "AliFMDRawReader.h" | |
17 | #include "AliFMDParameters.h" | |
18 | #include "AliFMDDigit.h" | |
19 | #include "AliTPCAltroMapping.h" | |
20 | #include "AliTriggerConfiguration.h" | |
21 | #include "TAlienCollection.h" | |
22 | #include "AliTriggerClass.h" | |
23 | #include "TGrid.h" | |
24 | #include "TPRegexp.h" | |
25 | #include "AliVZERORawStream.h" | |
26 | ||
27 | #include "AliITSsegmentationUpgrade.h" | |
28 | #include "AliRunLoader.h" | |
29 | #include "AliStack.h" | |
30 | #include "AliITSRecPointU.h" | |
31 | #include "AliITSDigitUpgrade.h" | |
32 | #include "TParticle.h" | |
33 | #include "AliESDtrack.h" | |
34 | #include "AliESDEvent.h" | |
35 | #include "AliTrackReference.h" | |
36 | #include <TROOT.h> | |
37 | #include <TStyle.h> | |
38 | #include "TGraphErrors.h" | |
39 | #include <TF1.h> | |
40 | #include <TGeoManager.h> | |
41 | ||
42 | #endif | |
43 | #include <iostream> | |
44 | #include <iomanip> | |
45 | #include <cstdio> | |
46 | ||
47 | /* | |
48 | .L ~/ITSupgrade/CDRpics/FastVsSlowSim.C | |
49 | ExtractOutputHistos(0,1); | |
50 | ||
51 | ||
52 | .L ~/ITSupgrade/BuildDetector/DetectorK.cxx++ | |
53 | .L FastVsSlowSim.C | |
54 | ||
55 | FastVsSlowSimRes(); | |
56 | FastVsSlowSimPtRes(); | |
57 | FastVsSlowSimEff(0,1); | |
58 | FastVsSlowSimEff(1,1); | |
59 | FastVsSlowSimEff(2,1); | |
60 | FastVsSlowSimEff(3,1); | |
61 | ||
62 | FastVsSlowSimEff(0); | |
63 | ||
64 | */ | |
65 | ||
66 | ||
67 | Int_t atLeastcorr = 3; | |
68 | ||
69 | ||
70 | void FastVsSlowSim(Bool_t extract=1) { | |
71 | ||
72 | if (extract) | |
73 | ExtractOutputHistos(0,1); | |
74 | else | |
75 | plotMerged(); | |
76 | } | |
77 | ||
78 | TGraph *gr[5]; | |
79 | Int_t colors[5]={1,2,3,4,6}; | |
80 | Int_t width =2; | |
81 | ||
82 | ||
83 | // new ideal Pixel properties? | |
84 | ||
85 | Double_t etaCut = 0.9; | |
86 | Double_t X0 = 0.003; | |
87 | Double_t resRPhi = 0.0004; | |
88 | Double_t resZ = 0.0004; | |
89 | ||
90 | void FastVsSlowSimRes() { | |
91 | ||
92 | Int_t plusTPC =0; | |
93 | ||
94 | gROOT->LoadMacro("~/fig_template.C"); // figure style | |
95 | myOptions(0); | |
96 | gROOT->ForceStyle(); | |
97 | ||
98 | TCanvas *myCan = new TCanvas("myCan"); | |
99 | myCan->Draw(); | |
100 | myCan->cd(); | |
101 | ||
102 | TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); | |
103 | myPadSetUp(myPad,0.15,0.04,0.04,0.15); | |
104 | myPad->Draw(); myPad->cd(); | |
105 | myPad->SetGridx(); myPad->SetGridy(); myPad->SetLogx(); | |
106 | ||
107 | // TLegend *leg = new TLegend(0.7,160,20,290,"","brCDN"); | |
108 | TLegend *leg = new TLegend(0.44,160,1.7,290,"","brCDN"); | |
109 | ||
110 | leg->SetFillColor(0); | |
111 | ||
112 | ||
113 | // Current ITS +++++++++++++++++++++++++++++++++++++++++ | |
114 | ||
115 | DetectorK its("ALICE","ITS"); | |
116 | its.MakeAliceCurrent(0,plusTPC); | |
117 | its.SetMaxRadiusOfSlowDetectors(0.1); | |
118 | its.SolveViaBilloir(0); | |
119 | Int_t color=1; Int_t linewidth=2; | |
120 | ||
121 | TGraph *c[6]; | |
122 | TGraph *d[6]; | |
123 | ||
124 | Int_t pi =0; | |
125 | d[pi] = its.GetGraphPointingResolution(1,color,linewidth); | |
126 | d[pi]->SetLineStyle(2); | |
127 | // d[pi]->GetYaxis()->SetTitle("Pointing resolution #sigma [#mum]"); | |
128 | // d[pi]->SetTitle("Pointing resolution .vs. Pt"); | |
129 | // d[pi]->Draw("AC"); | |
130 | ||
131 | c[pi] = its.GetGraphPointingResolution(0,color,linewidth); | |
132 | c[pi]->SetMinimum(-1); | |
133 | c[pi]->Draw("AC"); | |
134 | ||
135 | leg->AddEntry(c[pi],"FastTool: Current ITS","l"); | |
136 | // leg->AddEntry(d[pi],"in z - Current ITS","l"); | |
137 | ||
138 | ||
139 | ||
140 | ||
141 | // Current ITS +++++++++++++++++++++++++++++++++++++++++ | |
142 | ||
143 | Int_t color=3; Int_t linewidth=2; | |
144 | Int_t pi =2; | |
145 | ||
146 | DetectorK its("ALICE","ITS"); | |
147 | its.MakeAliceCurrent(0,plusTPC); | |
148 | ||
149 | its.SetRadius("bpipe",2.0); | |
150 | its.AddLayer("spd0", 2.2,1,1,1); | |
151 | ||
152 | its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); | |
153 | its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); | |
154 | its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); | |
155 | ||
156 | its.SetMaxRadiusOfSlowDetectors(0.1); | |
157 | its.SolveViaBilloir(0); | |
158 | ||
159 | d[pi] = its.GetGraphPointingResolution(1,color,linewidth); | |
160 | d[pi]->SetLineStyle(2); | |
161 | // d[pi]->Draw("C"); | |
162 | ||
163 | c[pi] = its.GetGraphPointingResolution(0,color,linewidth); | |
164 | c[pi]->Draw("C"); | |
165 | ||
166 | leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); | |
167 | // leg->AddEntry(d[pi],"in z - \"New SPDs\"","l"); | |
168 | ||
169 | ||
170 | ||
171 | // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ | |
172 | ||
173 | color=2; Int_t linewidth=2; | |
174 | Int_t pi =1; | |
175 | ||
176 | ||
177 | // for a 0.8,0.2 weight configuration | |
178 | ||
179 | DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); | |
180 | ||
181 | itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe | |
182 | itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation | |
183 | ||
184 | itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); | |
185 | itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); | |
186 | itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); | |
187 | itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); | |
188 | itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); | |
189 | itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); | |
190 | itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); | |
191 | ||
192 | if(plusTPC) itsU->AddTPC(0.1,0.1); | |
193 | itsU->SetMaxRadiusOfSlowDetectors(0.1); | |
194 | itsU->SolveViaBilloir(0); | |
195 | itsU->PrintLayout(); | |
196 | ||
197 | ||
198 | d[pi] = itsU->GetGraphPointingResolution(1,color,linewidth); | |
199 | d[pi]->SetLineStyle(2); | |
200 | // d[pi]->Draw("C"); | |
201 | ||
202 | c[pi] = itsU->GetGraphPointingResolution(0,color,linewidth); | |
203 | c[pi]->SetMaximum(150); | |
204 | c[pi]->Draw("C"); | |
205 | ||
206 | leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); | |
207 | // leg->AddEntry(d[pi],"in z - \"All New\" ","l"); | |
208 | ||
209 | ||
210 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
211 | ||
212 | ||
213 | ||
214 | TFile f1("root/FastVsSlow_CurrentITS-PbPb-fran.root"); | |
215 | TFile f2("root/FastVsSlow_NewSPDs-PbPb-fran.root"); | |
216 | TFile f3("root/FastVsSlow_AllNew-PbPb-fran.root"); | |
217 | TGraphErrors *dca1 = (TGraphErrors*)f1.Get("dca"); | |
218 | TGraphErrors *dca2 = (TGraphErrors*)f2.Get("dca"); | |
219 | TGraphErrors *dca3 = (TGraphErrors*)f3.Get("dca"); | |
220 | ||
221 | dca1->SetMarkerStyle(21); dca1->SetMarkerColor(1); | |
222 | dca2->SetMarkerStyle(21); dca2->SetMarkerColor(3); | |
223 | dca3->SetMarkerStyle(21); dca3->SetMarkerColor(2); | |
224 | ||
225 | leg->AddEntry(dca1,"FullMC: Current ITS","PE"); | |
226 | leg->AddEntry(dca2,"FullMC: \"New SPDs\"","PE"); | |
227 | leg->AddEntry(dca3,"FullMC: \"All New\" ","PE"); | |
228 | ||
229 | dca1->Draw("APE"); dca1->SetMinimum(-1); dca1->SetMaximum(300); | |
230 | dca2->Draw("PE"); | |
231 | dca3->Draw("PE"); | |
232 | c[0]->Draw("C"); | |
233 | c[1]->Draw("C"); | |
234 | c[2]->Draw("C"); | |
235 | ||
236 | leg->Draw(); | |
237 | ||
238 | myCan->SaveAs(Form("FastVsSlowSim-Res-%d.pdf",plusTPC)); | |
239 | myCan->SaveAs(Form("FastVsSlowSim-Res-%d.eps",plusTPC)); | |
240 | ||
241 | ||
242 | } | |
243 | ||
244 | // ============================================================================================== | |
245 | // ============================================================================================== | |
246 | ||
247 | void FastVsSlowSimPtRes() { | |
248 | ||
249 | Int_t plusTPC =0; | |
250 | ||
251 | gROOT->LoadMacro("~/fig_template.C"); // figure style | |
252 | myOptions(0); | |
253 | gROOT->ForceStyle(); | |
254 | ||
255 | TCanvas *myCan = new TCanvas("myCan"); | |
256 | myCan->Draw(); | |
257 | myCan->cd(); | |
258 | ||
259 | TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); | |
260 | myPadSetUp(myPad,0.15,0.04,0.04,0.15); | |
261 | myPad->Draw(); myPad->cd(); | |
262 | myPad->SetGridx(); myPad->SetGridy(); myPad->SetLogx(); myPad->SetLogy(); | |
263 | ||
264 | ||
265 | TLegend *leg = new TLegend(0.44,0.13,0.1.7,0.9,"","brCDN"); | |
266 | leg->SetFillColor(0); | |
267 | ||
268 | TGraph *c[6]; | |
269 | ||
270 | // Current ITS +++++++++++++++++++++++++++++++++++++++++ | |
271 | Int_t color=1; Int_t linewidth=2; | |
272 | ||
273 | Int_t pi =0; | |
274 | ||
275 | DetectorK its("ALICE","ITS"); | |
276 | its.MakeAliceCurrent(0,plusTPC); | |
277 | its.SetMaxRadiusOfSlowDetectors(0.1); | |
278 | its.SolveViaBilloir(0); | |
279 | Int_t color=1; Int_t linewidth=2; | |
280 | ||
281 | c[pi] = its.GetGraphMomentumResolution(color,linewidth); | |
282 | c[pi]->Draw("AC"); | |
283 | ||
284 | leg->AddEntry(c[pi],"FastTool: Current ITS","l"); | |
285 | ||
286 | ||
287 | // Current ITS +++++++++++++++++++++++++++++++++++++++++ | |
288 | ||
289 | Int_t color=3; Int_t linewidth=2; | |
290 | Int_t pi =2; | |
291 | ||
292 | DetectorK its("ALICE","ITS"); | |
293 | its.MakeAliceCurrent(0,plusTPC); | |
294 | ||
295 | its.SetRadius("bpipe",2.0); | |
296 | its.AddLayer("spd0", 2.2,1,1,1); | |
297 | ||
298 | its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); | |
299 | its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); | |
300 | its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); | |
301 | ||
302 | its.SetMaxRadiusOfSlowDetectors(0.1); | |
303 | its.SolveViaBilloir(0); | |
304 | ||
305 | c[pi] = its.GetGraphMomentumResolution(color,linewidth); | |
306 | c[pi]->Draw("C"); | |
307 | ||
308 | leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); | |
309 | ||
310 | ||
311 | // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ | |
312 | ||
313 | color=2; Int_t linewidth=2; | |
314 | Int_t pi =1; | |
315 | ||
316 | ||
317 | // for a 0.8,0.2 weight configuration | |
318 | ||
319 | DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); | |
320 | ||
321 | itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe | |
322 | itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation | |
323 | ||
324 | itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); | |
325 | itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); | |
326 | itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); | |
327 | itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); | |
328 | itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); | |
329 | itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); | |
330 | itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); | |
331 | ||
332 | if(plusTPC) itsU->AddTPC(0.1,0.1); | |
333 | itsU->SetMaxRadiusOfSlowDetectors(0.1); | |
334 | itsU->SolveViaBilloir(0); | |
335 | itsU->PrintLayout(); | |
336 | ||
337 | ||
338 | c[pi] = itsU->GetGraphMomentumResolution(color,linewidth); | |
339 | c[pi]->Draw("C"); | |
340 | ||
341 | leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); | |
342 | ||
343 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
344 | ||
345 | ||
346 | ||
347 | TFile f1("root/FastVsSlow_CurrentITS-PbPb-fran.root"); | |
348 | TFile f2("root/FastVsSlow_NewSPDs-PbPb-fran.root"); | |
349 | TFile f3("root/FastVsSlow_AllNew-PbPb-fran.root"); | |
350 | ||
351 | TGraphErrors *gpt1 = (TGraphErrors*)f1.Get("dPt"); | |
352 | TGraphErrors *gpt2 = (TGraphErrors*)f2.Get("dPt"); | |
353 | TGraphErrors *gpt3 = (TGraphErrors*)f3.Get("dPt"); | |
354 | ||
355 | gpt1->SetMarkerStyle(21); gpt1->SetMarkerColor(1); | |
356 | gpt2->SetMarkerStyle(21); gpt2->SetMarkerColor(3); | |
357 | gpt3->SetMarkerStyle(21); gpt3->SetMarkerColor(2); | |
358 | ||
359 | leg->AddEntry(gpt1,"FullMC: Current ITS","PE"); | |
360 | leg->AddEntry(gpt2,"FullMC: \"New SPDs\"","PE"); | |
361 | leg->AddEntry(gpt3,"FullMC: \"All New\" ","PE"); | |
362 | ||
363 | gpt1->Draw("APE"); gpt1->SetMinimum(0.1); gpt1->SetMaximum(20); | |
364 | gpt2->Draw("PE"); | |
365 | gpt3->Draw("PE"); | |
366 | c[0]->Draw("C"); | |
367 | c[1]->Draw("C"); | |
368 | c[2]->Draw("C"); | |
369 | ||
370 | leg->Draw(); | |
371 | ||
372 | myCan->SaveAs(Form("FastVsSlowSim-PtRes-%d.pdf",plusTPC)); | |
373 | myCan->SaveAs(Form("FastVsSlowSim-PtRes-%d.eps",plusTPC)); | |
374 | ||
375 | ||
376 | ||
377 | } | |
378 | ||
379 | ||
380 | ||
381 | // ============================================================================================== | |
382 | // ============================================================================================== | |
383 | ||
384 | void FastVsSlowSimEff(Int_t id=0,Int_t PbPb=0) { | |
385 | ||
386 | Int_t mult = 2400; // 2800 // deducted from "Frackable" | |
387 | if (PbPb) mult=2800; | |
388 | ||
389 | ||
390 | Int_t plusTPC =0; | |
391 | ||
392 | gROOT->LoadMacro("~/fig_template.C"); // figure style | |
393 | myOptions(0); | |
394 | gROOT->ForceStyle(); | |
395 | ||
396 | TCanvas *myCan = new TCanvas("myCan"); | |
397 | myCan->Draw(); | |
398 | myCan->cd(); | |
399 | ||
400 | TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); | |
401 | myPadSetUp(myPad,0.15,0.04,0.04,0.15); | |
402 | myPad->Draw(); myPad->cd(); | |
403 | myPad->SetGridx(); myPad->SetGridy();// myPad->SetLogx(); | |
404 | ||
405 | ||
406 | TLegend *leg = new TLegend(0.9,30,1.7,70,"","brCDN"); | |
407 | leg->SetFillColor(0); | |
408 | ||
409 | TGraph *c[6]; | |
410 | if (id!=2) { | |
411 | ||
412 | ||
413 | // Current ITS +++++++++++++++++++++++++++++++++++++++++ | |
414 | Int_t color=1; Int_t linewidth=2; | |
415 | ||
416 | Int_t pi =0; | |
417 | ||
418 | DetectorK its("ALICE","ITS"); | |
419 | its.MakeAliceCurrent(0,plusTPC); | |
420 | its.SetMaxRadiusOfSlowDetectors(0.01); | |
421 | its.SetAtLeastCorr(atLeastcorr); | |
422 | if (PbPb) its.SetdNdEtaCent(mult); | |
423 | its.SolveViaBilloir(0); | |
424 | Int_t color=1; Int_t linewidth=2; | |
425 | ||
426 | if (id==0) | |
427 | c[pi] = its.GetGraphRecoEfficiency(0,color,linewidth); | |
428 | else if (id==1) | |
429 | c[pi] = its.GetGraphRecoPurity(0,color,linewidth); | |
430 | else | |
431 | c[pi] = its.GetGraphRecoFakes(0,color,linewidth); | |
432 | ||
433 | c[pi]->Draw("AC"); | |
434 | ||
435 | leg->AddEntry(c[pi],"FastTool: Current ITS","l"); | |
436 | ||
437 | ||
438 | // NEW SPD +++++++++++++++++++++++++++++++++++++++++ | |
439 | ||
440 | Int_t color=3; Int_t linewidth=2; | |
441 | Int_t pi =2; | |
442 | ||
443 | DetectorK its("ALICE","ITS"); | |
444 | its.MakeAliceCurrent(0,plusTPC); | |
445 | its.SetAtLeastCorr(atLeastcorr); | |
446 | if (PbPb) its.SetdNdEtaCent(mult); | |
447 | its.SetRadius("bpipe",2.0); | |
448 | its.AddLayer("spd0", 2.2,1,1,1); | |
449 | ||
450 | its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); | |
451 | its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); | |
452 | its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); | |
453 | ||
454 | its.SetMaxRadiusOfSlowDetectors(0.1); | |
455 | its.SolveViaBilloir(0); | |
456 | ||
457 | if (id==0) | |
458 | c[pi] = its.GetGraphRecoEfficiency(0,color,linewidth); | |
459 | else if (id==1) | |
460 | c[pi] = its.GetGraphRecoPurity(0,color,linewidth); | |
461 | else | |
462 | c[pi] = its.GetGraphRecoFakes(0,color,linewidth); | |
463 | ||
464 | c[pi]->Draw("C"); | |
465 | ||
466 | leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); | |
467 | ||
468 | ||
469 | // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ | |
470 | ||
471 | color=4; Int_t linewidth=2; | |
472 | Int_t pi =1; | |
473 | ||
474 | ||
475 | // for a 0.8,0.2 weight configuration | |
476 | ||
477 | DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); | |
478 | itsU->SetAtLeastCorr(atLeastcorr); | |
479 | itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe | |
480 | itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation | |
481 | ||
482 | itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); | |
483 | itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); | |
484 | itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); | |
485 | itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); | |
486 | itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); | |
487 | itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); | |
488 | // itsU->AddLayer("ddd6", 42.6 , X0, resRPhi, resZ); | |
489 | itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); | |
490 | // itsU->AddLayer("ddd8", 43.4 , X0, resRPhi, resZ); | |
491 | ||
492 | ||
493 | if (PbPb) itsU->SetdNdEtaCent(mult); | |
494 | if(plusTPC) itsU->AddTPC(0.1,0.1); | |
495 | itsU->SetMaxRadiusOfSlowDetectors(0.1); | |
496 | itsU->SolveViaBilloir(0); | |
497 | itsU->PrintLayout(); | |
498 | ||
499 | if (id==0) | |
500 | c[pi] = itsU->GetGraphRecoEfficiency(0,color,linewidth); | |
501 | else if (id==1) | |
502 | c[pi] = itsU->GetGraphRecoPurity(0,color,linewidth); | |
503 | else | |
504 | c[pi] = itsU->GetGraphRecoFakes(0,color,linewidth); | |
505 | c[pi]->Draw("C"); | |
506 | ||
507 | leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); | |
508 | ||
509 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
510 | ||
511 | // ALL NEW - double outer layer +++++++++++++++++++++++++++++++++++++ | |
512 | ||
513 | color=2; Int_t linewidth=2; | |
514 | Int_t pi =3; | |
515 | ||
516 | ||
517 | // for a 0.8,0.2 weight configuration | |
518 | ||
519 | DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); | |
520 | itsU->SetAtLeastCorr(atLeastcorr); | |
521 | itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe | |
522 | itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation | |
523 | ||
524 | itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); | |
525 | itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); | |
526 | itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); | |
527 | itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); | |
528 | itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); | |
529 | itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); | |
530 | itsU->AddLayer("ddd8", 40.0 , X0, resRPhi, resZ); | |
531 | itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); | |
532 | itsU->AddLayer("ddd9", 43.4 , X0, resRPhi, resZ); | |
533 | ||
534 | ||
535 | if (PbPb) itsU->SetdNdEtaCent(mult); | |
536 | if(plusTPC) itsU->AddTPC(0.1,0.1); | |
537 | itsU->SetMaxRadiusOfSlowDetectors(0.1); | |
538 | itsU->SolveViaBilloir(0); | |
539 | itsU->PrintLayout(); | |
540 | ||
541 | if (id==0) | |
542 | c[pi] = itsU->GetGraphRecoEfficiency(0,color,linewidth); | |
543 | else if (id==1) | |
544 | c[pi] = itsU->GetGraphRecoPurity(0,color,linewidth); | |
545 | else | |
546 | c[pi] = itsU->GetGraphRecoFakes(0,color,linewidth); | |
547 | c[pi]->Draw("C"); | |
548 | ||
549 | leg->AddEntry(c[pi],"FastTool: \"All New\" (2x double layer)","l"); | |
550 | ||
551 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
552 | ||
553 | ||
554 | } | |
555 | ||
556 | char h[100]; | |
557 | if (PbPb==0) | |
558 | sprintf(h,"-fran"); | |
559 | else | |
560 | sprintf(h,"-Anna"); | |
561 | ||
562 | TFile f1(Form("root/FastVsSlow_CurrentITS-PbPb%s.root",h)); | |
563 | TFile f2(Form("root/FastVsSlow_NewSPDs-PbPb%s.root",h)); | |
564 | TFile f3(Form("root/FastVsSlow_AllNew-PbPb%s.root",h)); | |
565 | TFile f4(Form("root/FastVsSlow_AllNew-9-PbPb%s.root",h)); | |
566 | ||
567 | // TFile f1(Form("root/FastVsSlow_CurrentITS%s-fran.root",h)); | |
568 | // TFile f2(Form("root/FastVsSlow_NewSPDs%s-fran.root",h)); | |
569 | // TFile f3(Form("root/FastVsSlow_AllNew%s-fran.root",h)); | |
570 | ||
571 | TH1F *eff1 = 0; | |
572 | TH1F *eff2 = 0; | |
573 | TH1F *eff3 = 0; | |
574 | TH1F *eff4 = 0; | |
575 | if (id==0) { | |
576 | eff1 = (TH1F*)f1.Get("efficiency"); | |
577 | eff2 = (TH1F*)f2.Get("efficiency"); | |
578 | eff3 = (TH1F*)f3.Get("efficiency"); | |
579 | eff4 = (TH1F*)f4.Get("efficiency"); | |
580 | eff1->GetYaxis()->SetTitle("efficiency (%)"); | |
581 | } else if (id==1) { | |
582 | eff1 = (TH1F*)f1.Get("purity"); | |
583 | eff2 = (TH1F*)f2.Get("purity"); | |
584 | eff3 = (TH1F*)f3.Get("purity"); | |
585 | eff4 = (TH1F*)f4.Get("purity"); | |
586 | eff1->GetYaxis()->SetTitle("purity (%)"); | |
587 | } else if (id==2) { | |
588 | eff1 = (TH1F*)f1.Get("annaEff"); | |
589 | eff2 = (TH1F*)f2.Get("annaEff"); | |
590 | eff3 = (TH1F*)f3.Get("annaEff"); | |
591 | eff4 = (TH1F*)f4.Get("annaEff"); | |
592 | eff1->GetYaxis()->SetTitle("Overall efficiency (%)"); | |
593 | } else if (id==3) { | |
594 | eff1 = (TH1F*)f1.Get("fake"); | |
595 | eff2 = (TH1F*)f2.Get("fake"); | |
596 | eff3 = (TH1F*)f3.Get("fake"); | |
597 | eff4 = (TH1F*)f4.Get("fake"); | |
598 | eff1->GetYaxis()->SetTitle("Fake ratio (%)"); | |
599 | } | |
600 | ||
601 | eff1->SetMarkerStyle(21); eff1->SetMarkerColor(1); | |
602 | eff2->SetMarkerStyle(21); eff2->SetMarkerColor(3); | |
603 | eff3->SetMarkerStyle(21); eff3->SetMarkerColor(4); | |
604 | eff4->SetMarkerStyle(21); eff4->SetMarkerColor(2); | |
605 | ||
606 | leg->AddEntry(eff1,"FullMC: Current ITS","PE"); | |
607 | leg->AddEntry(eff2,"FullMC: \"New SPDs\"","PE"); | |
608 | leg->AddEntry(eff3,"FullMC: \"All New\" ","PE"); | |
609 | leg->AddEntry(eff4,"FullMC: \"All New\" (2x double layer)","PE"); | |
610 | ||
611 | eff1->SetMinimum(0.4); eff1->SetMaximum(100); | |
612 | eff1->DrawCopy("E"); | |
613 | eff2->DrawCopy("sameE"); | |
614 | eff4->DrawCopy("sameE"); | |
615 | eff3->DrawCopy("sameE"); | |
616 | if (id!=2) { | |
617 | c[0]->Draw("C"); | |
618 | c[1]->Draw("C"); | |
619 | c[2]->Draw("C"); | |
620 | c[3]->Draw("C"); | |
621 | } | |
622 | eff2->DrawCopy("sameE"); | |
623 | eff4->DrawCopy("sameE"); | |
624 | eff3->DrawCopy("sameE"); | |
625 | ||
626 | ||
627 | leg->Draw(); | |
628 | ||
629 | ||
630 | ||
631 | ||
632 | TPaveText *pt = 0; | |
633 | if (id!=3) | |
634 | pt = new TPaveText(0.4,0.1,1.76,30); | |
635 | else | |
636 | pt = new TPaveText(0.4,70,1.76,100); | |
637 | ||
638 | pt->SetBorderSize(1); // no shadow | |
639 | pt->SetTextFont(12); | |
640 | TText *t1 = pt->AddText("FastTool settings: "); t1->SetTextFont(32); // bold | |
641 | ||
642 | pt->AddText(Form(" Tracked particles: Pions; Average rapidity: 0.45; dN_{ch}/d#eta = %d ",mult)); | |
643 | ||
644 | // pt->AddText("\"New SPDs\": layer radii: r = {2.2,4.8,9.1} cm"); | |
645 | // pt->AddText("\"All New\: layer radii: r = {2.2,3.8,6.8,...} cm"); | |
646 | // pt->AddText(Form(" New layer prop.: X/X_{0}=%1.1lf%%; #sigma_{r#phi,z}=%1.0lf#mum",X0*100,resZ*1e4)); | |
647 | ||
648 | TText *t2 = pt->AddText("FullMC settings: "); t2->SetTextFont(32); // bold | |
649 | if (PbPb==0) { | |
650 | pt->AddText(" Generator: AliGenHIJINGpara (parametrized PbPb event)"); | |
651 | pt->AddText(" dN_{ch.pr.}/d#eta = 2070"); | |
652 | pt->AddText(" Track selection: Pions, |#eta|<0.9"); | |
653 | } else { | |
654 | pt->AddText(" Generator: AliGenHijing (modified); #sqrt{s_{NN}} = 5.5 TeV"); | |
655 | pt->AddText(" dN_{ch.pr.}/d#eta = 2410; Impactparameter range: b#in(0,5) #rightarrow central PbPb "); | |
656 | pt->AddText(" Track selection: Pions, |#eta|<0.9"); | |
657 | } | |
658 | // pt->SetLabel("Settings"); | |
659 | pt->SetTextAlign(12); | |
660 | pt->SetFillColor(0); | |
661 | pt->Draw(); | |
662 | ||
663 | ||
664 | ||
665 | ||
666 | ||
667 | if (PbPb==0) { | |
668 | myCan->SaveAs(Form("FastVsSlowSim-Eff-%d.pdf",id)); | |
669 | myCan->SaveAs(Form("FastVsSlowSim-Eff-%d.eps",id)); | |
670 | }else{ | |
671 | myCan->SaveAs(Form("FastVsSlowSim-Eff-PbPb-%d.pdf",id)); | |
672 | myCan->SaveAs(Form("FastVsSlowSim-Eff-PbPb-%d.eps",id)); | |
673 | } | |
674 | ||
675 | ||
676 | ||
677 | ||
678 | } | |
679 | ||
680 | ||
681 | ||
682 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
683 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
684 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
685 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
686 | // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
687 | // Extraction | |
688 | ||
689 | ||
690 | void GetDetectorRadii(TArrayD *rmin,TArrayD *rmax) { | |
691 | // Loads geometry of the ITS upgrade | |
692 | AliITSsegmentationUpgrade *seg=new AliITSsegmentationUpgrade(); | |
693 | Int_t nlayers = seg->GetNLayers(); | |
694 | rmin->Set(nlayers); rmax->Set(nlayers); | |
695 | for (Int_t i=0; i<nlayers; i++) { | |
696 | rmin->AddAt(seg->GetRadius(i),i); | |
697 | rmax->AddAt(seg->GetRadius(i)+seg->GetThickness(i),i); | |
698 | } | |
699 | } | |
700 | ||
701 | void PrintDetectorGeometry() { | |
702 | gSystem->Load("libITSUpgradeBase"); | |
703 | gSystem->Load("libITSUpgradeSim"); | |
704 | ||
705 | TArrayD rmin(0); | |
706 | TArrayD rmax(0); | |
707 | GetDetectorRadii(&rmin,&rmax); | |
708 | for (Int_t i=0; i<rmin.GetSize(); i++) { | |
709 | cout<<i<<": (rmin,rmax)=("<<rmin.At(i)<<","<<rmax.At(i)<<")cm"<<endl; | |
710 | } | |
711 | } | |
712 | ||
713 | ||
714 | void CountTrackableMCs(TH1F *hAllMC=0, Bool_t onlyPrims=0,Bool_t onlyPion=0); | |
715 | void CountPrimaries(TH1F *hMultCount); | |
716 | ||
717 | ||
718 | void ExtractOutputHistos(Bool_t onlyPrims=0,Bool_t onlyPion=0,Int_t plotFlag=0) { | |
719 | ||
720 | // gROOT->SetStyle("Plain"); | |
721 | gStyle->SetPalette(1); | |
722 | ||
723 | const Int_t nbins=20; | |
724 | Double_t ptmin=0.06;//04; | |
725 | Double_t ptmax=2.0;//GeV | |
726 | Double_t logxmin = TMath::Log10(ptmin); | |
727 | Double_t logxmax = TMath::Log10(ptmax); | |
728 | Double_t binwidth = (logxmax-logxmin)/(nbins+1); | |
729 | enum {nb=nbins+1}; | |
730 | Double_t xbins[nb]; | |
731 | xbins[0] = ptmin; | |
732 | for (Int_t i=1;i<=nbins;i++) { | |
733 | xbins[i] = ptmin + TMath::Power(10,logxmin+(i)*binwidth); | |
734 | // cout<<xbins[i]<<endl; | |
735 | } | |
736 | // TH1F *h = new TH1F("h","hist with log x axis",nbins,xbins); | |
737 | ||
738 | TH1F *hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.); | |
739 | hMultCount->GetXaxis()->SetTitle("eta"); | |
740 | hMultCount->GetYaxis()->SetTitle("N/d#eta"); | |
741 | ||
742 | TH1F *hAllMC = new TH1F("allMC","All Tracks MC primaries",nbins,xbins); | |
743 | TH1F *hAllFound = new TH1F("allFound","All Tracks found",nbins,xbins); | |
744 | TH1F *hImperfect = new TH1F("imperfect","Imperfect tracks",nbins,xbins); | |
745 | TH1F *hPerfect = new TH1F("perfect","Perfect tracks",nbins,xbins); | |
746 | TH1F *hEff = new TH1F("efficiency","Efficiency (Perfect tracks in \"ALL MC\")",nbins,xbins); | |
747 | TH1F *hFake = new TH1F("fake","Fake tracks (Inperfect tracks in \"ALL MC\")",nbins,xbins); | |
748 | TH1F *hPurity = new TH1F("purity","Purity (Perfect tracks in \"All Found\")",nbins,xbins); | |
749 | TH1F *hAnna = new TH1F("annaEff","AnnalisaEff ",nbins,xbins); | |
750 | TH1F *hNoMCTrack = new TH1F("noMCtrack","noMCtrack ",nbins,xbins); | |
751 | ||
752 | TH1F *hEta = new TH1F("","",50,-2,2); | |
753 | // TH1F *hEtaMC = new TH1F("","",50,-2,2); | |
754 | ||
755 | TH2D *h2Ddca = new TH2D("dca2D","DCAvsPt2D",nbins,xbins,50,-0.05,0.05); | |
756 | TH2D *h2Dpt = new TH2D("dPt2D","dPtdvsPt2D",nbins,xbins,50,-25,25); | |
757 | ||
758 | // open run loader and load gAlice, kinematics and header | |
759 | AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); | |
760 | if (!runLoader) { | |
761 | Error("Check kine", "getting run loader from file %s failed", | |
762 | "galice.root"); | |
763 | return; | |
764 | } | |
765 | runLoader->LoadgAlice(); | |
766 | gAlice = runLoader->GetAliRun(); | |
767 | if (!gAlice) { | |
768 | Error("Check kine", "no galice object found"); | |
769 | return; | |
770 | } | |
771 | runLoader->LoadHeader(); | |
772 | runLoader->LoadKinematics(); | |
773 | ||
774 | TFile* esdFile = TFile::Open("AliESDs.root"); | |
775 | if (!esdFile || !esdFile->IsOpen()) { | |
776 | Error("CheckESD", "opening ESD file %s failed", "AliESDs.root"); | |
777 | return; | |
778 | } | |
779 | AliESDEvent *esd = new AliESDEvent(); | |
780 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
781 | if (!tree) { | |
782 | Error("CheckESD", "no ESD tree found"); | |
783 | return; | |
784 | } | |
785 | esd->ReadFromTree(tree); | |
786 | ||
787 | Int_t nTrackTotalMC = 0; | |
788 | Int_t nTrackFound = 0; | |
789 | Int_t nTrackImperfect = 0; | |
790 | Int_t nTrackPerfect = 0; | |
791 | Int_t nNoMCTrack = 0; | |
792 | ||
793 | ||
794 | for(Int_t iEv =0; iEv<tree->GetEntries(); iEv++){ | |
795 | tree->GetEvent(iEv); | |
796 | runLoader->GetEvent(iEv); | |
797 | ||
798 | printf("+++ event %i (of %lld) +++++++++++++++++++++++ # ESDtracks: %d \n",iEv,tree->GetEntries()-1,esd->GetNumberOfTracks()); | |
799 | Int_t nESDtracks = esd->GetNumberOfTracks(); | |
800 | for (Int_t iTrack = 0; iTrack < nESDtracks; iTrack++) { | |
801 | AliESDtrack* track = esd->GetTrack(iTrack); | |
802 | if (!(iTrack%1000)) printf("event %i: ESD track count %d (of %d)\n",iEv,iTrack,nESDtracks); | |
803 | ||
804 | Int_t label = track->GetLabel(); | |
805 | ||
806 | Int_t idx[12]; | |
807 | // Int_t ncl = track->GetITSclusters(idx); | |
808 | ||
809 | if(label<0) { | |
810 | // cout<< " ESD track label " << label; | |
811 | // cout<<" ---> imperfect track (label "<<label<<"<0) !! -> track Pt: "<< track->Pt() << endl; | |
812 | } | |
813 | ||
814 | AliStack* stack = runLoader->Stack(); | |
815 | // nTrackTotalMC += stack->GetNprimary(); | |
816 | ||
817 | ||
818 | TParticle* particle = stack->Particle(TMath::Abs(label)); | |
819 | Double_t pt = track->Pt(); | |
820 | ||
821 | if(particle) { | |
822 | ||
823 | if (TMath::Abs(particle->Eta())>etaCut) continue; | |
824 | ||
825 | Double_t ptMC = particle->Pt(); | |
826 | ||
827 | // Efficiencies | |
828 | if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue; | |
829 | ||
830 | if ( (!onlyPrims) || stack->IsPhysicalPrimary(TMath::Abs(label))) { | |
831 | // cout<<" # clusters "<<ncl<<endl; | |
832 | ||
833 | nTrackFound++; | |
834 | hAllFound->Fill(ptMC); | |
835 | hEta->Fill(track->Eta()); | |
836 | ||
837 | if (label<0) { | |
838 | nTrackImperfect++; | |
839 | hImperfect->Fill(ptMC); | |
840 | } else { | |
841 | nTrackPerfect++; | |
842 | hPerfect->Fill(ptMC); | |
843 | } | |
844 | ||
845 | } | |
846 | ||
847 | ||
848 | // following only for "true tracks, pions | |
849 | ||
850 | if(particle->Pt() < 0.001)continue; | |
851 | if (TMath::Abs(particle->GetPdgCode())!=211) continue; | |
852 | if (label>0) { | |
853 | ||
854 | // Impact parameters for Pions only | |
855 | Double_t dca = track->GetD(0,0,0.5); | |
856 | h2Ddca->Fill(ptMC,dca); | |
857 | ||
858 | // Pt resolution for Pions only | |
859 | Double_t dPt = (pt-ptMC)/ptMC*100; | |
860 | h2Dpt->Fill(ptMC,dPt); | |
861 | } | |
862 | ||
863 | } else { | |
864 | nNoMCTrackFound++; | |
865 | hNoMCTrack->Fill(pt); | |
866 | cout<<" according MC particle not found"<<endl; | |
867 | } | |
868 | ||
869 | } //entries track esd | |
870 | ||
871 | ||
872 | }//entries tree | |
873 | runLoader->UnloadHeader(); | |
874 | runLoader->UnloadKinematics(); | |
875 | delete runLoader; | |
876 | ||
877 | ||
878 | // Count trackable MC tracks | |
879 | CountTrackableMCs(hAllMC, onlyPrims, onlyPion); | |
880 | ||
881 | ||
882 | // Count trackable MC tracks | |
883 | CountPrimaries(hMultCount); | |
884 | ||
885 | ||
886 | ||
887 | ||
888 | // Get Errors right | |
889 | hMultCount->Sumw2(); | |
890 | hAllMC->Sumw2(); | |
891 | hAllFound->Sumw2(); | |
892 | hPerfect->Sumw2(); | |
893 | hImperfect->Sumw2(); | |
894 | h2Dpt->Sumw2(); | |
895 | h2Ddca->Sumw2(); | |
896 | ||
897 | // -- Global efficienies | |
898 | ||
899 | nTrackTotalMC = hAllMC->GetEntries(); | |
900 | Double_t eff = ((Double_t)nTrackPerfect)/nTrackTotalMC; | |
901 | printf("-> Total number of events: %lld -> MCtracks %d -> nPerfect %d -> Eff: %3.2lf \n", | |
902 | tree->GetEntries(),nTrackTotalMC,nTrackPerfect,eff); | |
903 | ||
904 | Double_t purity = ((Double_t)nTrackPerfect)/nTrackFound; | |
905 | printf("-> Total number of events: %lld -> FoundTracks %d -> nPerfect %d -> Purity: %3.2lf \n", | |
906 | tree->GetEntries(),nTrackFound,nTrackPerfect,purity); | |
907 | ||
908 | // Efficiencies - and normalize to 100% | |
909 | ||
910 | TF1 f1("f1","100+x*0",0.,1.e3); | |
911 | ||
912 | hPurity->Divide(hPerfect,hAllFound,1,1,"b"); | |
913 | hPurity->Multiply(&f1); | |
914 | hPurity->SetMarkerColor(kGreen); | |
915 | hPurity->SetMarkerStyle(21); | |
916 | hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
917 | hPurity->SetStats(0); | |
918 | ||
919 | hPurity->GetYaxis()->SetRangeUser(0,100); | |
920 | hPurity->SetTitle("Efficiency & Purity"); | |
921 | ||
922 | hEff->Divide(hPerfect,hAllMC,1,1,"b"); | |
923 | hEff->Multiply(&f1); | |
924 | hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
925 | hEff->SetMarkerColor(kBlue); | |
926 | hEff->SetMarkerStyle(21); | |
927 | hEff->SetStats(0); | |
928 | ||
929 | hFake->Divide(hImperfect,hAllMC,1,1,"b"); | |
930 | hFake->Multiply(&f1); | |
931 | hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
932 | hFake->SetMarkerColor(kRed); | |
933 | hFake->SetMarkerStyle(21); | |
934 | hFake->SetStats(0); | |
935 | ||
936 | ||
937 | hAnna->Divide(hAllFound,hAllMC,1,1,"b"); | |
938 | hAnna->Multiply(&f1); | |
939 | hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
940 | hAnna->SetMarkerColor(kBlack); | |
941 | hAnna->SetMarkerStyle(21); | |
942 | hAnna->SetStats(0); | |
943 | ||
944 | TCanvas *c1 = new TCanvas("c1","NoMCTrackFound");//,200,10,900,900); | |
945 | TVirtualPad *pad = c1->cd(); | |
946 | pad->SetGridx(); pad->SetGridy(); | |
947 | hNoMCTrack->Draw(); | |
948 | ||
949 | TCanvas *c2 = new TCanvas("c2","Eff&Purity");//,200,10,900,900); | |
950 | TVirtualPad *pad = c2->cd(); | |
951 | pad->SetGridx(); pad->SetGridy(); | |
952 | // pad->SetLogx(); | |
953 | ||
954 | hPurity->Draw("E"); | |
955 | hEff->Draw("Same E"); | |
956 | hFake->Draw("Same E"); | |
957 | hAnna->Draw("Same E"); | |
958 | ||
959 | TLegend *leg = new TLegend(0.1,0.8,0.6,0.9);leg->SetFillColor(0); | |
960 | leg->AddEntry(hPurity,"Purity (\"Perfect tracks\" within \"Found Tracks\")","PE"); | |
961 | leg->AddEntry(hEff,"Efficiency (\"Perfect tracks\" within \"MC findable Tracks\")","PE"); | |
962 | leg->AddEntry(hFake,"Fake (\"Inperfect tracks\" within \"MC findable Tracks\")","PE"); | |
963 | leg->AddEntry(hAnna,"AnnaLisa - Efficiency (\"Found tracks\" within \"MC findable Tracks\")","PE"); | |
964 | leg->Draw(); | |
965 | ||
966 | ||
967 | if (plotFlag==1){ | |
968 | hAllMC->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
969 | hAllMC->Draw(); // MC pt distribution | |
970 | hAllFound->SetLineColor(2); | |
971 | hAllFound->Draw("same"); // MC pt distribution | |
972 | } | |
973 | ||
974 | ||
975 | /* | |
976 | ||
977 | .L ~/ITSupgrade/BuildDetector/DetectorK.cxx+ | |
978 | ||
979 | // All NEW | |
980 | DetectorK its("ALICE","ITS"); | |
981 | its.MakeAliceAllNew(0); | |
982 | its.SetMaxRadiusOfSlowDetectors(0.01); | |
983 | its.SolveViaBilloir(0); | |
984 | TGraph *c = its.GetGraphRecoEfficiency(0,3,2); | |
985 | c->Draw("C"); | |
986 | ||
987 | ||
988 | // Current | |
989 | DetectorK its("ALICE","ITS"); | |
990 | its.MakeAliceCurrent(0,0); | |
991 | its.SetMaxRadiusOfSlowDetectors(0.01); | |
992 | its.SolveViaBilloir(0); | |
993 | TGraph *c = its.GetGraphRecoEfficiency(0,4,2); | |
994 | c->Draw("C"); | |
995 | ||
996 | */ | |
997 | ||
998 | TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900); | |
999 | c3->Divide(2,1); c3->cd(1); | |
1000 | // Impact parameter | |
1001 | ||
1002 | // Impact parameter resolution --------------- | |
1003 | h2Ddca->Draw("colz"); | |
1004 | h2Ddca->FitSlicesY() ; | |
1005 | TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same"); | |
1006 | TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw(); | |
1007 | TGraphErrors *d0 = new TGraphErrors(); | |
1008 | for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) { | |
1009 | d0->SetPoint( ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns | |
1010 | d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns | |
1011 | } | |
1012 | d0->SetMarkerStyle(21); | |
1013 | d0->SetMaximum(200); d0->SetMinimum(0); | |
1014 | d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1015 | d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)"); | |
1016 | d0->SetName("dca"); d0->SetTitle("DCAvsPt"); | |
1017 | ||
1018 | c3->cd(1); h2Ddca->Draw("surf2"); | |
1019 | c3->cd(2); d0->Draw("APE"); | |
1020 | ||
1021 | // PT RESOLUTION ------------ | |
1022 | TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900); | |
1023 | c4->Divide(2,1); c4->cd(1); | |
1024 | // Impact parameter | |
1025 | h2Dpt->Draw("colz"); | |
1026 | h2Dpt->FitSlicesY() ; | |
1027 | TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same"); | |
1028 | TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw(""); | |
1029 | TGraphErrors *gPt = new TGraphErrors(); | |
1030 | for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) { | |
1031 | gPt->SetPoint( ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); | |
1032 | gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); | |
1033 | } | |
1034 | gPt->SetMarkerStyle(21); | |
1035 | gPt->SetMaximum(20); gPt->SetMinimum(0); | |
1036 | gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1037 | gPt->GetYaxis()->SetTitle("relative momentum resolution (%)"); | |
1038 | gPt->SetName("dPt"); gPt->SetTitle("DPTvsPt"); | |
1039 | ||
1040 | c4->cd(1); h2Dpt->Draw("surf2"); | |
1041 | c4->cd(2); gPt->Draw("APE"); | |
1042 | ||
1043 | ||
1044 | // EXPORT -------- | |
1045 | ||
1046 | TFile f("histos.root","RECREATE"); | |
1047 | ||
1048 | hMultCount->Write(); | |
1049 | hAllMC->Write(); | |
1050 | hAllFound->Write(); | |
1051 | hImperfect->Write(); | |
1052 | hPerfect->Write(); | |
1053 | hNoMCTrack->Write(); | |
1054 | ||
1055 | hPurity->Write(); | |
1056 | hEff->Write(); | |
1057 | hFake->Write(); | |
1058 | hAnna->Write(); | |
1059 | ||
1060 | h2Ddca->Write(); | |
1061 | d0->Write(); | |
1062 | ||
1063 | h2Dpt->Write(); | |
1064 | gPt->Write(); | |
1065 | ||
1066 | f.Close(); | |
1067 | ||
1068 | return; | |
1069 | ||
1070 | } | |
1071 | ||
1072 | void CountTrackableMCs(TH1F *hAllMC, Bool_t onlyPrims,Bool_t onlyPion) { | |
1073 | ||
1074 | gSystem->Load("libITSUpgradeBase"); | |
1075 | gSystem->Load("libITSUpgradeSim"); | |
1076 | ||
1077 | // open run loader and load gAlice, kinematics and header | |
1078 | AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); | |
1079 | if (!runLoader) { | |
1080 | Error("Check kine", "getting run loader from file %s failed", | |
1081 | "galice.root"); | |
1082 | return; | |
1083 | } | |
1084 | runLoader->LoadHeader(); | |
1085 | runLoader->LoadKinematics(); | |
1086 | runLoader->LoadTrackRefs(); | |
1087 | ||
1088 | AliLoader *dl = runLoader->GetDetectorLoader("ITS"); | |
1089 | ||
1090 | //Trackf | |
1091 | TTree *trackRefTree = 0x0; | |
1092 | TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000); | |
1093 | ||
1094 | // TH1F *hRef = new TH1F("","",100,0,100); | |
1095 | TH1F *hR = new TH1F("","",100,0,100); | |
1096 | if (hAllMC==0) hAllMC = new TH1F("","",100,0.1,2); | |
1097 | Float_t ptmin = hAllMC->GetBinCenter(1)-hAllMC->GetBinWidth(1)/2; | |
1098 | Float_t ptmax = hAllMC->GetBinCenter(hAllMC->GetNbinsX())+hAllMC->GetBinWidth(hAllMC->GetNbinsX())/2; | |
1099 | // Int_t nAllMC = 0; | |
1100 | ||
1101 | // Detector geometry | |
1102 | TArrayD rmin(0); TArrayD rmax(0); | |
1103 | GetDetectorRadii(&rmin,&rmax); | |
1104 | TArrayI nLaySigs(rmin.GetSize()); | |
1105 | ||
1106 | printf("Counting trackable MC tracks ...\n"); | |
1107 | ||
1108 | for(Int_t iEv =0; iEv<runLoader->GetNumberOfEvents(); iEv++){ | |
1109 | Int_t nTrackableTracks = 0; | |
1110 | runLoader->GetEvent(iEv); | |
1111 | AliStack* stack = runLoader->Stack(); | |
1112 | printf("+++ event %i (of %d) +++++++++++++++++++++++ # total MCtracks: %d \n",iEv,runLoader->GetNumberOfEvents()-1,stack->GetNtrack()); | |
1113 | ||
1114 | trackRefTree=runLoader->TreeTR(); | |
1115 | TBranch *br = trackRefTree->GetBranch("TrackReferences"); | |
1116 | if(!br) { | |
1117 | printf("no TR branch available , exiting \n"); | |
1118 | return; | |
1119 | } | |
1120 | br->SetAddress(&trackRef); | |
1121 | ||
1122 | // init the trackRef tree | |
1123 | trackRefTree=runLoader->TreeTR(); | |
1124 | trackRefTree->SetBranchAddress("TrackReferences",&trackRef); | |
1125 | ||
1126 | // Count trackable MC tracks | |
1127 | for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++) { | |
1128 | ||
1129 | TParticle* particle = stack->Particle(iMC); | |
1130 | if (TMath::Abs(particle->Eta())>etaCut) continue; | |
1131 | if (onlyPrims && !stack->IsPhysicalPrimary(iMC)) continue; | |
1132 | if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue; | |
1133 | ||
1134 | ||
1135 | Bool_t isTrackable = 0; | |
1136 | nLaySigs.Reset(0); | |
1137 | ||
1138 | trackRefTree->GetEntry(stack->TreeKEntry(iMC)); | |
1139 | Int_t nref=trackRef->GetEntriesFast(); | |
1140 | for(Int_t iref =0; iref<nref; iref++){ | |
1141 | AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref); | |
1142 | if(!trR) continue; | |
1143 | if(trR->DetectorId()!=AliTrackReference::kITS) continue; | |
1144 | Float_t radPos = trR->R(); | |
1145 | hR->Fill(radPos); | |
1146 | for (Int_t il=0; il<rmin.GetSize();il++) { | |
1147 | if (radPos>=rmin.At(il)-0.1 && radPos<=rmax.At(il)+0.1) { | |
1148 | // cout<<" in Layer "<<il<<" "<<radPos; | |
1149 | nLaySigs.AddAt(1.,il); | |
1150 | // cout<<" "<<nLaySigs.At(il)<<endl; | |
1151 | } | |
1152 | } | |
1153 | } | |
1154 | ||
1155 | if (nLaySigs.GetSum()>=3) { | |
1156 | isTrackable =1; | |
1157 | // cout<<nLaySigs.GetSum()<<endl; | |
1158 | } | |
1159 | ||
1160 | if (isTrackable) { | |
1161 | Double_t ptMC = particle->Pt(); | |
1162 | // Double_t etaMC = particle->Eta(); | |
1163 | // if (ptMC>ptmin&&ptMC<ptmax) {nTrackableTracks++;hAllMC->Fill(ptMC);} | |
1164 | if (ptMC>ptmin) {nTrackableTracks++;hAllMC->Fill(ptMC);} | |
1165 | ||
1166 | } | |
1167 | ||
1168 | ||
1169 | } // entries tracks MC | |
1170 | printf(" -> trackable MC tracks: %d (%d)\n",nTrackableTracks,hAllMC->GetEntries()); | |
1171 | }//entries Events | |
1172 | ||
1173 | ||
1174 | hR->DrawCopy(); | |
1175 | hAllMC->DrawCopy(); | |
1176 | runLoader->UnloadHeader(); | |
1177 | runLoader->UnloadKinematics(); | |
1178 | delete runLoader; | |
1179 | ||
1180 | } | |
1181 | ||
1182 | void CountPrimaries(TH1F *hMultCount) { | |
1183 | ||
1184 | if (hMultCount==0) hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.); | |
1185 | ||
1186 | AliRunLoader *rl = AliRunLoader::Open("galice.root"); | |
1187 | rl->SetKineFileName("Kinematics.root"); | |
1188 | rl->LoadHeader(); | |
1189 | rl->LoadKinematics(); | |
1190 | Int_t nEvents = rl->GetNumberOfEvents(); | |
1191 | cout<< "N events "<<nEvents<<endl; | |
1192 | for(Int_t iEv=0; iEv<nEvents; iEv++){ | |
1193 | rl->GetEvent(iEv); | |
1194 | AliStack *s = rl->Stack(); | |
1195 | for(Int_t iP=0; iP<s->GetNtrack(); iP++ ){ | |
1196 | TParticle *p = s->Particle(iP); | |
1197 | if (!(s->IsPhysicalPrimary(iP))) continue; | |
1198 | Float_t eta = p->Eta(); | |
1199 | if (p->Pt()>0.06) { | |
1200 | hMultCount->Fill(eta); | |
1201 | } | |
1202 | } | |
1203 | } | |
1204 | ||
1205 | hMultCount->DrawCopy(); | |
1206 | rl->UnloadHeader(); | |
1207 | rl->UnloadKinematics(); | |
1208 | delete rl; | |
1209 | ||
1210 | ||
1211 | ||
1212 | } | |
1213 | ||
1214 | ||
1215 | void plotMerged(Bool_t onlyPlot=0) { | |
1216 | ||
1217 | gStyle->SetPalette(1); | |
1218 | ||
1219 | TFile f("histoSum.root","UPDATE"); | |
1220 | ||
1221 | TH1F* hAllMC = f.Get("allMC"); | |
1222 | TH1F* hAllFound= f.Get("allFound"); | |
1223 | TH1F* hImperfect= f.Get("imperfect"); | |
1224 | TH1F* hPerfect= f.Get("perfect"); | |
1225 | TH1F* hNoMCTrack= f.Get("noMCtrack"); | |
1226 | ||
1227 | ||
1228 | // have to be recalculated | |
1229 | TH1F* hPurity = f.Get("purity"); | |
1230 | TH1F* hEff= f.Get("efficiency"); | |
1231 | TH1F* hFake= f.Get("fake"); | |
1232 | TH1F* hAnna= f.Get("annaEff"); | |
1233 | ||
1234 | TH2D* h2Ddca= f.Get("dca2D"); | |
1235 | TGraphErrors *d0= f.Get("dca"); | |
1236 | ||
1237 | TH2D* h2Dpt= f.Get("dPt2D"); | |
1238 | TGraphErrors *gPt= f.Get("dPt"); | |
1239 | ||
1240 | ||
1241 | if (!onlyPlot) { | |
1242 | /* // Get Errors right | |
1243 | hAllMC->Sumw2(); | |
1244 | hAllFound->Sumw2(); | |
1245 | hPerfect->Sumw2(); | |
1246 | hImperfect->Sumw2(); | |
1247 | h2Dpt->Sumw2(); | |
1248 | h2Ddca->Sumw2(); | |
1249 | */ | |
1250 | ||
1251 | // Efficiencies - and normalize to 100% | |
1252 | ||
1253 | TF1 f1("f1","100+x*0",0.,1.e3); | |
1254 | ||
1255 | hPurity->Divide(hPerfect,hAllFound,1,1,"b"); | |
1256 | hPurity->Multiply(&f1); | |
1257 | hPurity->SetMarkerColor(kGreen); | |
1258 | hPurity->SetMarkerStyle(21); | |
1259 | hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1260 | hPurity->SetStats(0); | |
1261 | ||
1262 | hPurity->GetYaxis()->SetRangeUser(0,100); | |
1263 | hPurity->SetTitle("Efficiency & Purity"); | |
1264 | ||
1265 | hEff->Divide(hPerfect,hAllMC,1,1,"b"); | |
1266 | hEff->Multiply(&f1); | |
1267 | hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1268 | hEff->SetMarkerColor(kBlue); | |
1269 | hEff->SetMarkerStyle(21); | |
1270 | hEff->SetStats(0); | |
1271 | ||
1272 | hFake->Divide(hImperfect,hAllMC,1,1,"b"); | |
1273 | hFake->Multiply(&f1); | |
1274 | hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1275 | hFake->SetMarkerColor(kRed); | |
1276 | hFake->SetMarkerStyle(21); | |
1277 | hFake->SetStats(0); | |
1278 | ||
1279 | hAnna->Divide(hAllFound,hAllMC,1,1,"b"); | |
1280 | hAnna->Multiply(&f1); | |
1281 | hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1282 | hAnna->SetMarkerColor(kBlack); | |
1283 | hAnna->SetMarkerStyle(21); | |
1284 | hAnna->SetStats(0); | |
1285 | ||
1286 | ||
1287 | // Impact parameter resolution --------------- | |
1288 | TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900); | |
1289 | c3->Divide(2,1); c3->cd(1); | |
1290 | h2Ddca->DrawCopy("colz"); | |
1291 | h2Ddca->FitSlicesY() ; | |
1292 | TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same"); | |
1293 | TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw(); | |
1294 | TGraphErrors *d0 = new TGraphErrors(); | |
1295 | for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) { | |
1296 | d0->SetPoint( ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns | |
1297 | d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns | |
1298 | } | |
1299 | d0->SetMarkerStyle(21); | |
1300 | d0->SetMaximum(200); d0->SetMinimum(0); | |
1301 | d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1302 | d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)"); | |
1303 | d0->SetName("dca"); d0->SetTitle("DCAvsPt"); | |
1304 | // c3->cd(1); h2Ddca->Draw("surf2"); | |
1305 | c3->cd(2); d0->Draw("APE"); | |
1306 | ||
1307 | // PT RESOLUTION ------------ | |
1308 | TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900); | |
1309 | c4->Divide(2,1); c4->cd(1); | |
1310 | h2Dpt->DrawCopy("colz"); | |
1311 | h2Dpt->FitSlicesY() ; | |
1312 | TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same"); | |
1313 | TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw(""); | |
1314 | TGraphErrors *gPt = new TGraphErrors(); | |
1315 | for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) { | |
1316 | gPt->SetPoint( ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); | |
1317 | gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); | |
1318 | } | |
1319 | gPt->SetMarkerStyle(21); | |
1320 | gPt->SetMaximum(20); gPt->SetMinimum(0); | |
1321 | gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); | |
1322 | gPt->GetYaxis()->SetTitle("relative momentum resolution (%)"); | |
1323 | gPt->SetName("dPt"); gPt->SetTitle("DPTvsPt"); | |
1324 | // c4->cd(1); h2Dpt->Draw("surf2"); | |
1325 | c4->cd(2); gPt->Draw("APE"); | |
1326 | ||
1327 | ||
1328 | // overwrite with normalized graphs | |
1329 | hPurity->Write(); | |
1330 | hEff->Write(); | |
1331 | hFake->Write(); | |
1332 | hAnna->Write(); | |
1333 | h2Ddca->Write(); | |
1334 | d0->Write(); | |
1335 | h2Dpt->Write(); | |
1336 | gPt->Write(); | |
1337 | ||
1338 | } | |
1339 | ||
1340 | // Plots | |
1341 | ||
1342 | TCanvas *c2 = new TCanvas("c2","Eff&Purity");//,200,10,900,900); | |
1343 | TVirtualPad *pad = c2->cd(); | |
1344 | pad->SetGridx(); pad->SetGridy(); | |
1345 | // pad->SetLogx(); | |
1346 | ||
1347 | TLegend *leg = new TLegend(0.1,0.8,0.6,0.9);leg->SetFillColor(0); | |
1348 | leg->AddEntry(hPurity,"Purity (\"Perfect tracks\" within \"Found Tracks\")","PE"); | |
1349 | leg->AddEntry(hEff,"Efficiency (\"Perfect tracks\" within \"MC findable Tracks\")","PE"); | |
1350 | leg->AddEntry(hFake,"Fake (\"Inperfect tracks\" within \"MC findable Tracks\")","PE"); | |
1351 | leg->AddEntry(hAnna,"AnnaLisa - Efficiency (\"Found tracks\" within \"MC findable Tracks\")","PE"); | |
1352 | ||
1353 | ||
1354 | hPurity->DrawCopy("E"); | |
1355 | hEff->DrawCopy("Same E"); | |
1356 | hFake->DrawCopy("Same E"); | |
1357 | hAnna->DrawCopy("Same E"); | |
1358 | leg->Draw(); | |
1359 | ||
1360 | c2->SaveAs("EffPlot.png"); | |
1361 | ||
1362 | f.Close(); | |
1363 | ||
1364 | ||
1365 | ||
1366 | } | |
1367 |