ITS/UPGRADE/macros/FastVsSlowSim.C

   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