]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/macros/FastVsSlowSim.C
Update master to aliroot
[u/mrichter/AliRoot.git] / ITS / UPGRADE / macros / FastVsSlowSim.C
CommitLineData
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
67Int_t atLeastcorr = 3;
68
69
70void FastVsSlowSim(Bool_t extract=1) {
71
72 if (extract)
73 ExtractOutputHistos(0,1);
74 else
75 plotMerged();
76}
77
78TGraph *gr[5];
79Int_t colors[5]={1,2,3,4,6};
80Int_t width =2;
81
82
83// new ideal Pixel properties?
84
85Double_t etaCut = 0.9;
86Double_t X0 = 0.003;
87Double_t resRPhi = 0.0004;
88Double_t resZ = 0.0004;
89
90void 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
247void 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
384void 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
690void 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
701void 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
714void CountTrackableMCs(TH1F *hAllMC=0, Bool_t onlyPrims=0,Bool_t onlyPion=0);
715void CountPrimaries(TH1F *hMultCount);
716
717
718void 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
1072void 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
1182void 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
1215void 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