File indexing completed on 2025-08-03 08:13:10
0001 #include <TH2D.h>
0002 #include <TFile.h>
0003 #include <TF1.h>
0004 #include <TGraph.h>
0005 #include <TCanvas.h>
0006 #include <TColor.h>
0007 #include <TROOT.h>
0008 #include <TStyle.h>
0009 #include <TLegend.h>
0010 #include <TLatex.h>
0011 #include <TLine.h>
0012 #include <TEfficiency.h>
0013 #include <TGraphAsymmErrors.h>
0014
0015 #include "sPhenixStyle.C"
0016
0017 void look_purity()
0018 {
0019 SetsPhenixStyle();
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 bool hijing = false;
0031
0032
0033 int n_maps_layer = 0;
0034
0035 int n_intt_layer = 0;
0036 int n_tpc_layer = 48;
0037
0038 double ptmax = 40.0;
0039 double hptmax = 40.0;
0040
0041
0042
0043
0044 double ptstep = 0.4;
0045
0046
0047
0048 bool pt_resolution_fit = true;
0049
0050
0051 TFile *fin = new TFile("root_files/purity_out.root");
0052
0053
0054
0055 if(!fin)
0056 {
0057 cout << "Failed to find input file" << endl;
0058 exit(1);
0059 }
0060
0061 TCanvas *c1 = new TCanvas("c1","c1",5,5,900,500);
0062 c1->Divide(3,1);
0063
0064
0065
0066 TH2D *hpt_dca2d = 0;
0067 fin->GetObject("hpt_dca2d",hpt_dca2d);
0068 c1->cd(1);
0069 gPad->SetLeftMargin(0.12);
0070 hpt_dca2d->GetYaxis()->SetTitleOffset(1.5);
0071 hpt_dca2d->SetMarkerStyle(20);
0072 hpt_dca2d->SetMarkerSize(0.3);
0073 hpt_dca2d->Draw();
0074
0075
0076 TH2D *hpt_dcaZ = 0;
0077 fin->GetObject("hpt_dcaZ",hpt_dcaZ);
0078 c1->cd(2);
0079 gPad->SetLeftMargin(0.12);
0080 hpt_dcaZ->SetMarkerStyle(20);
0081 hpt_dcaZ->SetMarkerSize(0.3);
0082 hpt_dcaZ->GetYaxis()->SetTitleOffset(1.5);
0083 hpt_dcaZ->Draw();
0084
0085 TH2D *hpt_compare = 0;
0086 fin->GetObject("hpt_compare",hpt_compare);
0087 c1->cd(3);
0088 gPad->SetLeftMargin(0.12);
0089 hpt_compare->GetYaxis()->SetTitleOffset(1.5);
0090 hpt_compare->SetMarkerStyle(20);
0091 hpt_compare->SetMarkerSize(0.3);
0092 hpt_compare->Draw();
0093
0094
0095
0096
0097
0098 TCanvas *c2 = new TCanvas("c2","c2",20,20,800,600);
0099
0100 int NPT = (int) (ptmax/ptstep);
0101
0102 double pT[200];
0103 double dca2d[200];
0104 for(int i = 0;i<NPT;i++)
0105 {
0106
0107
0108 double ptlo = (double) i * ptstep + 0.0;
0109 double pthi = (double) i * ptstep + ptstep;
0110
0111 int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0112 int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0113
0114 std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0115
0116 TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
0117 hpt_dca2d->ProjectionY("h",binlo,binhi);
0118 h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0119 h->GetXaxis()->SetTitle("#Delta dca2d (cm)");
0120 h->GetXaxis()->SetTitleOffset(1.0);
0121
0122 h->DrawCopy();
0123
0124 double mean = h->GetMean();
0125 double sigma = h->GetRMS();
0126 double yield = h->Integral();
0127
0128 pT[i] = (ptlo + pthi) / 2.0;
0129
0130 double low = -0.01, high=0.01;
0131 if(n_maps_layer == 0)
0132 {
0133 low = 3.0 * low;
0134 high = 3.0 * high;
0135 }
0136 if(pT[i] < 6.0)
0137 {
0138 low = 3.0*low;
0139 high = 3.0*high;
0140 }
0141
0142 TF1 *f = new TF1("f","gaus",low, high);
0143 f->SetParameter(1, yield/100.0);
0144 f->SetParameter(2, 0.0);
0145 f->SetParameter(3,0.002);
0146 h->Fit(f,"R");
0147
0148
0149
0150 dca2d[i] = f->GetParameter(2);
0151 cout << " pT " << pT[i] << " dca2d " << dca2d[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0152 }
0153
0154
0155
0156
0157 TCanvas *c3 = new TCanvas("c3","c3",100,100,800,600);
0158 TGraph *grdca2d = new TGraph(NPT,pT,dca2d);
0159 grdca2d->SetMarkerStyle(20);
0160 grdca2d->SetMarkerSize(1.2);
0161 grdca2d->SetMarkerColor(kRed);
0162 grdca2d->SetName("dca2d_resolution");
0163 grdca2d->SetTitle("dca2d resolution");
0164
0165 TH1D *hdummy = new TH1D("hdummy","#Delta dca2d vs p_{T}",100,0.0,hptmax);
0166 hdummy->SetMinimum(0);
0167 if(n_maps_layer == 0)
0168 hdummy->SetMaximum(0.015);
0169 else
0170 hdummy->SetMaximum(0.0051);
0171 hdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0172 hdummy->GetYaxis()->SetTitle("DCA(r#phi) resolution (cm)");
0173 hdummy->GetYaxis()->SetTitleOffset(1.5);
0174 gPad->SetLeftMargin(0.15);
0175 hdummy->Draw();
0176 grdca2d->Draw("p");
0177
0178 char dcalab[500];
0179 if(hijing)
0180
0181 sprintf(dcalab, "sPHENIX HIJING b=0-12 fm, 100 kHz pileup");
0182 else
0183 sprintf(dcalab, "sPHENIX 100 pions");
0184 TLegend *ldca = new TLegend(0.3, 0.75, 1.05, 0.90,dcalab,"NDC");
0185 ldca->SetBorderSize(0);
0186 ldca->SetFillColor(0);
0187 ldca->SetFillStyle(0);
0188 char lstr1[500];
0189
0190 sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0191 ldca->AddEntry(grdca2d, lstr1, "p");
0192 ldca->Draw();
0193
0194
0195
0196
0197
0198
0199
0200 TCanvas *c2a = new TCanvas("c2a","c2a",20,20,800,600);
0201
0202 double dcaZ[200];
0203 for(int i = 0;i<NPT;i++)
0204 {
0205
0206
0207 double ptlo = (double) i * ptstep + 0.0;
0208 double pthi = (double) i * ptstep + ptstep;
0209
0210 int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
0211 int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
0212
0213 std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0214 TH1D *h = new TH1D("h","DCA (Z) resolution (cm)",2000, -0.1, 0.1);
0215 hpt_dcaZ->ProjectionY("h",binlo,binhi);
0216 h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0217 h->GetYaxis()->SetTitle("DCA(Z) resolution (cm)");
0218 h->GetYaxis()->SetTitleOffset(1.0);
0219
0220 h->DrawCopy();
0221
0222 double mean = h->GetMean();
0223 double sigma = h->GetRMS();
0224 double yield = h->Integral();
0225
0226 pT[i] = (ptlo + pthi) / 2.0;
0227
0228 double low = -0.01, high=0.01;
0229 if(n_maps_layer == 0)
0230 {
0231 low = 3.0 * low;
0232 high = 3.0 * high;
0233 }
0234
0235 if(pT[i] < 6.0)
0236 {
0237 low = 3.0*low;
0238 high = 3.0*high;
0239 }
0240
0241 TF1 *f = new TF1("f","gaus",low, high);
0242 f->SetParameter(1, yield/100.0);
0243 f->SetParameter(2, 0.0);
0244 f->SetParameter(3,0.002);
0245 h->Fit(f,"R");
0246
0247 dcaZ[i] = f->GetParameter(2);
0248 cout << " pT " << pT[i] << " dcaZ " << dcaZ[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0249 }
0250
0251
0252
0253
0254 TCanvas *c3a = new TCanvas("c3a","c3a",100,100,800,600);
0255 TGraph *grdcaZ = new TGraph(NPT,pT,dcaZ);
0256 grdcaZ->SetMarkerStyle(20);
0257 grdcaZ->SetMarkerSize(1.2);
0258 grdcaZ->SetMarkerColor(kRed);
0259 grdcaZ->SetName("dcaZ_resolution");
0260 grdcaZ->SetTitle("dcaZ resolution");
0261
0262 TH1D *hdummya = new TH1D("hdummya","#Delta dcaZ vs p_{T}",100,0.0,hptmax);
0263 hdummya->SetMinimum(0);
0264
0265 if(n_maps_layer == 0)
0266 hdummya->SetMaximum(0.051);
0267 else
0268 hdummya->SetMaximum(0.0051);
0269 hdummya->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0270 hdummya->GetYaxis()->SetTitle("DCA(Z) resolution (cm)");
0271 hdummya->GetYaxis()->SetTitleOffset(1.5);
0272 gPad->SetLeftMargin(0.15);
0273 hdummya->Draw();
0274 grdcaZ->Draw("p");
0275
0276 TLegend *ldcaZ = new TLegend(0.3, 0.75, 1.05, 0.90, dcalab,"NDC");
0277 ldcaZ->SetBorderSize(0);
0278 ldcaZ->SetFillColor(0);
0279 ldcaZ->SetFillStyle(0);
0280
0281
0282 sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0283 ldcaZ->AddEntry(grdcaZ, lstr1, "p");
0284 ldcaZ->Draw();
0285
0286
0287
0288
0289
0290
0291 TCanvas *c4 = new TCanvas("c4","c4",60,60,800,600);
0292
0293 double dpT[200];
0294
0295 dpT[0] = 1.0;
0296 for(int i = 1;i<NPT;i++)
0297
0298 {
0299
0300
0301 double ptlo = (double) i * ptstep + 0.0;
0302 double pthi = (double) i * ptstep + ptstep;
0303
0304 int binlo = hpt_compare->GetXaxis()->FindBin(ptlo);
0305 int binhi = hpt_compare->GetXaxis()->FindBin(pthi);
0306
0307 TH1D *hpt = new TH1D("hpt","pT resolution ",200, 0, 2.0);
0308 hpt_compare->ProjectionY("hpt",binlo,binhi);
0309 hpt->GetXaxis()->SetTitle("#Delta p_{T}/p_{T}");
0310 hpt->GetXaxis()->SetTitleOffset(1.0);
0311 if(i>30) hpt->Rebin(4);
0312 hpt->DrawCopy();
0313
0314 std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << " integral " << hpt->Integral() << std::endl;
0315
0316 TF1 *f = new TF1("f","gaus",0.8,1.2);
0317 hpt->Fit(f,"R");
0318
0319 pT[i] = (ptlo + pthi) / 2.0;
0320
0321 dpT[i] = f->GetParameter(2);
0322 cout << " pT " << pT[i] << " dpT " << dpT[i] << " integral " << hpt->Integral() << std::endl;
0323 }
0324
0325
0326
0327
0328 TCanvas *c5 = new TCanvas("c5","c5",100,100,800,800);
0329 c5->SetLeftMargin(0.14);
0330 TGraph *grdpt = new TGraph(NPT,pT,dpT);
0331 grdpt->SetMarkerStyle(20);
0332 grdpt->SetMarkerSize(1.1);
0333 grdpt->SetName("pt_resolution");
0334 grdpt->SetTitle("pT resolution");
0335
0336 TH1D *hdummy2 = new TH1D("hdummy2","#Delta p_{T} vs p_{T}",100,0.0,hptmax);
0337 hdummy2->SetMinimum(0);
0338
0339
0340 hdummy2->SetMaximum(0.3);
0341 hdummy2->GetXaxis()->SetTitle("p_{T}");
0342 hdummy2->GetYaxis()->SetTitle("#Delta p_{T}/p_{T}");
0343 hdummy2->GetYaxis()->SetTitleOffset(1.5);
0344 hdummy2->Draw();
0345 grdpt->Draw("p");
0346
0347 TLegend *ldpt= new TLegend(0.25, 0.75, 0.95, 0.90, dcalab, "NDC");
0348 ldpt->SetBorderSize(0);
0349 ldpt->SetFillColor(0);
0350 ldpt->SetFillStyle(0);
0351
0352
0353 sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0354 ldpt->AddEntry(grdpt, lstr1, "p");
0355 ldpt->Draw();
0356
0357
0358
0359
0360 TF1 *fpt = new TF1("fpt","[0] + [1]*x", 2.0, hptmax);
0361 fpt->SetParameter(0,0.005);
0362
0363 fpt->SetParameter(1,0.0015);
0364 if(pt_resolution_fit)
0365 grdpt->Fit(fpt,"R");
0366
0367 char lab[1000];
0368
0369 sprintf(lab,"#frac{#Deltap_{T}}{p_{T}} = %.4f + %.6f #times p_{T}", fpt->GetParameter(0), fpt->GetParameter(1));
0370 TLatex *mres = new TLatex(0.45,0.25,lab);
0371
0372 mres->SetNDC();
0373 if(pt_resolution_fit)
0374 mres->Draw();
0375
0376
0377 double pT_sigmas = 6.0;
0378 double const_term = fpt->GetParameter(0);
0379 double linear_term = fpt->GetParameter(1);
0380
0381 TH1D *hpt_truth;
0382 fin->GetObject("hpt_truth",hpt_truth);
0383 if(!hpt_truth)
0384 {
0385 cout << "Failed to get hpt_truth, quit!" << endl;
0386 exit(1);
0387 }
0388
0389
0390
0391 TH1D *hpt_matched = new TH1D("hpt_matched","hpt_matched", 500, 0.0, hptmax);
0392 double eff_pt[200];
0393
0394 for(int i = 0;i<NPT;i++)
0395 {
0396
0397
0398 double ptlo = (double) i * ptstep + 0.0;
0399 double pthi = (double) i * ptstep + ptstep;
0400 double ptval = (ptlo+pthi)/2.0;
0401
0402 int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0403 int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0404
0405 TH1D *hpt1 = new TH1D("hpt1","pT resolution ",2000, 0.0, 2.0);
0406 hpt_compare->ProjectionY("hpt1",binlo,binhi);
0407
0408 double ptres = sqrt(pow(const_term,2) + pow(linear_term * ptval,2));
0409 double ptreslo = 1.0 - ptres * pT_sigmas;
0410 double ptreshi = 1.0 + ptres * pT_sigmas;
0411
0412 int momlo = hpt1->GetXaxis()->FindBin(ptreslo);
0413 int momhi = hpt1->GetXaxis()->FindBin(ptreshi);
0414
0415 hpt_matched->Fill(ptval, hpt1->Integral(momlo, momhi));
0416
0417 int tlo = hpt_truth->FindBin(ptlo);
0418 int thi = hpt_truth->FindBin(pthi);
0419 double truth_yield = hpt_truth->Integral(tlo, thi);;
0420 eff_pt[i] = hpt1->Integral(momlo, momhi) / truth_yield;
0421
0422 cout << " pT " << ptval << " ptres " << ptres << " ptreslo " << ptreslo << " ptreshi " << ptreshi << " momlo " << momlo << " momhi " << momhi << endl;
0423 cout << " truth_yield " << truth_yield << " yield " << hpt1->Integral(momlo,momhi) << " eff_pt " << eff_pt[i] << endl;
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446 cout << " ptval " << ptval
0447 << " ptreslo " << ptreslo
0448 << " ptreshi " << ptreshi
0449 << " total " << hpt1->Integral()
0450 << " momlo " << momlo
0451 << " momhi " << momhi
0452 << " cut " << hpt1->Integral(momlo,momhi)
0453 << " tlo " << tlo
0454 << " thi " << thi
0455 << " truth " << truth_yield
0456 << " eff_pt " << eff_pt[i]
0457 << endl;
0458
0459 delete hpt1;
0460
0461 }
0462
0463 cout << " create canvas c7" << endl;
0464 TCanvas *c7 = new TCanvas("c7","c7",60,60,800,800);
0465
0466 TH1F *hd = new TH1F("hd","hd",100, 0.0, hptmax);
0467 hd->SetMinimum(0.0);
0468 hd->SetMaximum(1.0);
0469 hd->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0470 hd->GetYaxis()->SetTitle("Single track efficiency");
0471 hd->GetYaxis()->SetTitleOffset(1.0);
0472 hd->Draw();
0473
0474 TGraph *gr_eff = new TGraph(NPT,pT,eff_pt);
0475 gr_eff->SetName("single_track_efficiency");
0476 gr_eff->SetMarkerStyle(20);
0477 gr_eff->SetMarkerSize(1);
0478 gr_eff->SetMarkerColor(kRed);
0479
0480 gr_eff->Draw("p");
0481
0482 TLegend *leff = new TLegend(0.40, 0.35, 0.99, 0.50, dcalab, "NDC");
0483 leff->SetBorderSize(0);
0484 leff->SetFillColor(0);
0485 leff->SetFillStyle(0);
0486 char lstr[500];
0487
0488 sprintf(lstr,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0489 leff->AddEntry(gr_eff, lstr, "p");
0490
0491 leff->Draw();
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574 TCanvas *cpq = new TCanvas("cpq","cpq",40,40,1200,600);
0575 TH1D * hquality;
0576 fin->GetObject("hquality",hquality);
0577 hquality->Draw();
0578
0579
0580
0581 TCanvas *c8 = new TCanvas("c8","c8",40,40,1200,600);
0582 TH1D * hnhits;
0583 fin->GetObject("hnhits",hnhits);
0584 hnhits->GetXaxis()->SetTitle("hits per track");
0585 hnhits->Draw();
0586 cout << "hnhits integral = " << hnhits->Integral() << endl;
0587
0588
0589
0590
0591
0592 TCanvas *cdcadist = new TCanvas("cdcadist","cdcadist",5,20,1200,800);
0593 cdcadist->Divide(3,1);
0594
0595 double dcaptbinlo[3] = {0.5, 1.0, 2.0};
0596 double dcaptbinhi[3] = {1.0, 2.0, 50.0};
0597 for(int i = 0;i<3;i++)
0598 {
0599 double ptlo = dcaptbinlo[i];
0600 double pthi = dcaptbinhi[i];
0601
0602 int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0603 int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0604
0605 std::cout << "DCA distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0606
0607 TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
0608 hpt_dca2d->ProjectionY("h",binlo,binhi);
0609 h->GetXaxis()->SetTitle("DCA (r#phi) (cm)");
0610 h->GetXaxis()->SetTitleOffset(1.0);
0611 h->SetMinimum(0.5);
0612
0613 cdcadist->cd(i+1);
0614 gPad->SetLogy(1);
0615 h->Sumw2();
0616 h->Rebin(4);
0617
0618 h->SetMarkerSize(0.8);
0619 h->GetXaxis()->SetNdivisions(505);
0620 h->DrawCopy("E");
0621
0622 double mean = h->GetMean();
0623 double sigma = h->GetRMS();
0624 double yield = h->Integral();
0625
0626 double low = -0.01, high=0.01;
0627 if(n_maps_layer == 0)
0628 {
0629 low = 3.0 * low;
0630 high = 3.0 * high;
0631 }
0632
0633 if(pthi < 6.0)
0634 {
0635 low = 3.0*low;
0636 high = 3.0*high;
0637 }
0638
0639 TF1 *f = new TF1("f","gaus",low, high);
0640 f->SetParameter(1, yield/100.0);
0641 f->SetParameter(2, 0.0);
0642 f->SetParameter(3,0.002);
0643 f->SetLineColor(kRed);
0644 h->Fit(f,"R");
0645
0646 char fitr[500];
0647 sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
0648 TLatex *l1 = new TLatex(0.2,0.971,fitr);
0649 l1->SetNDC();
0650 l1->SetTextSize(0.07);
0651 l1->Draw();
0652
0653 sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
0654 TLatex *l1a = new TLatex(0.2,0.892,fitr);
0655 l1a->SetNDC();
0656 l1a->SetTextSize(0.07);
0657 l1a->Draw();
0658
0659 cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0660 }
0661
0662
0663
0664
0665 TCanvas *cdcazdist = new TCanvas("cdcazdist","cdcazdist",5,20,1200,800);
0666 cdcazdist->Divide(3,1);
0667
0668 for(int i = 0;i<3;i++)
0669 {
0670 double ptlo = dcaptbinlo[i];
0671 double pthi = dcaptbinhi[i];
0672
0673 int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
0674 int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
0675
0676 std::cout << "DCA Z distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0677
0678 TH1D *h = new TH1D("h","dcaZ resolution",2000, -0.1, 0.1);
0679 hpt_dcaZ->ProjectionY("h",binlo,binhi);
0680 h->GetXaxis()->SetTitle("DCA (Z) (cm)");
0681 h->GetXaxis()->SetTitleOffset(1.0);
0682
0683 cdcazdist->cd(i+1);
0684 gPad->SetLogy(1);
0685 h->Sumw2();
0686 h->Rebin(4);
0687 h->SetMarkerSize(0.8);
0688 h->GetXaxis()->SetNdivisions(505);
0689 h->SetMinimum(0.5);
0690 h->DrawCopy("E");
0691
0692 double mean = h->GetMean();
0693 double sigma = h->GetRMS();
0694 double yield = h->Integral();
0695
0696 double low = -0.01, high=0.01;
0697 if(n_maps_layer == 0)
0698 {
0699 low = 3.0 * low;
0700 high = 3.0 * high;
0701 }
0702
0703 if(pthi < 6.0)
0704 {
0705 low = 3.0*low;
0706 high = 3.0*high;
0707 }
0708
0709 TF1 *f = new TF1("f","gaus",low, high);
0710 f->SetParameter(1, yield/100.0);
0711 f->SetParameter(2, 0.0);
0712 f->SetParameter(3,0.002);
0713 f->SetLineColor(kRed);
0714 h->Fit(f,"R");
0715
0716 char fitr[500];
0717 sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
0718 TLatex *l1 = new TLatex(0.2,0.971,fitr);
0719 l1->SetNDC();
0720 l1->SetTextSize(0.07);
0721 l1->Draw();
0722
0723 sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
0724 TLatex *l1a = new TLatex(0.2,0.892,fitr);
0725 l1a->SetNDC();
0726 l1a->SetTextSize(0.07);
0727 l1a->Draw();
0728
0729 cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0730 }
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809 TCanvas *ceta = new TCanvas("ceta","ceta",10,10,600,600);
0810
0811 TH1D *hgeta = 0;
0812 fin->GetObject("hgeta",hgeta);
0813 if(!hgeta)
0814 {
0815 cout << "Did not get hgeta" << endl;
0816 exit(1);
0817 }
0818 TH1D *hreta = 0;
0819 fin->GetObject("hreta",hreta);
0820
0821 hreta->GetXaxis()->SetTitle("#eta");
0822 hreta->GetYaxis()->SetNdivisions(505);
0823 hreta->SetLineColor(kRed);
0824 hreta->SetMinimum(0.0);
0825 hreta->Draw();
0826
0827 hgeta->Draw("same");
0828
0829
0830 TLegend *leta = new TLegend(0.3, 0.35, 1.0, 0.60, dcalab, "NDC");
0831 leta->SetBorderSize(0);
0832 leta->SetFillColor(0);
0833 leta->SetFillStyle(0);
0834 leta->AddEntry(hgeta,"Truth","l");
0835 leta->AddEntry(hreta,"Reconstructed","l");
0836 leta->Draw();
0837 ceta->Update();
0838
0839 TCanvas *cfake = new TCanvas("cfake","cfake",4,4,800,600);
0840 cfake->SetLeftMargin(0.15);
0841 TH2D *hpt_nfake = 0;
0842 fin->GetObject("hpt_nfake",hpt_nfake);
0843 if(!hpt_nfake)
0844 {
0845 cout << "Did not get hpt_nfake" << endl;
0846 exit(1);
0847 }
0848
0849 int n_noise = 4;
0850 double ynoise[10][200];
0851 double ynoise_ref[200];
0852
0853 double dptfake = 0.2;
0854 double fake_ptmax = 20.0;
0855 int nptfake = (int) (fake_ptmax/dptfake);
0856 double ptfake[200];
0857 for(int ipt=0;ipt<nptfake;ipt++)
0858 {
0859 double ptlo = dptfake * ipt;
0860 double pthi = dptfake*ipt + dptfake;
0861 ptfake[ipt] = (ptlo+pthi) / 2.0;
0862 int binlo = hpt_nfake->GetXaxis()->FindBin(ptlo);
0863 int binhi = hpt_nfake->GetXaxis()->FindBin(pthi);
0864
0865 TH1D *hnoise = new TH1D("hnoise","",73,-5,68);
0866 hpt_nfake->ProjectionY("hnoise",binlo,binhi);
0867 double ynoise_ref = hnoise->Integral();
0868 for(int j=0;j<n_noise;j++)
0869 {
0870 int bin = hnoise->FindBin(j);
0871 ynoise[j][ipt]= hnoise->GetBinContent(bin);
0872 ynoise[j][ipt] /= ynoise_ref;
0873 }
0874 hnoise->Delete();
0875 }
0876
0877
0878 TH1D *hfdummy = new TH1D("hfdummy","",100,0.0,fake_ptmax);
0879 hfdummy->SetMaximum(1.05);
0880 hfdummy->SetMinimum(0.001);
0881 hfdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0882 hfdummy->GetYaxis()->SetTitle("Fraction of tracks");
0883 gPad->SetLogy(1);
0884 hfdummy->Draw();
0885
0886 TLegend *lfake = new TLegend(0.35,0.55,0.87,0.89,"","NDC");
0887 lfake->SetBorderSize(1);
0888 TGraph *hnoise[10];
0889 int fake_col[5] = {kBlack,kRed,kBlue,kMagenta,kViolet};
0890 for(int j=0;j<n_noise;j++)
0891 {
0892 hnoise[j] = new TGraph(nptfake,ptfake,ynoise[j]);
0893 hnoise[j]->SetMarkerStyle(20);
0894 hnoise[j]->SetMarkerSize(1.0);
0895 hnoise[j]->SetMarkerColor(fake_col[j]);
0896 if(j == 0)
0897 hnoise[j]->Draw("p");
0898 else
0899 hnoise[j]->Draw("p same");
0900 char lab[500];
0901 sprintf(lab,"%i mismatched hits",j);
0902 lfake->AddEntry(hnoise[j],lab,"p");
0903 }
0904 lfake->Draw();
0905
0906 TH2D *hcorr_nfake_nmaps = 0;
0907 fin->GetObject("hcorr_nfake_nmaps",hcorr_nfake_nmaps);
0908 if(hcorr_nfake_nmaps)
0909 {
0910 TCanvas *cmm = new TCanvas("cmm","cmm",4,4,800,600);
0911 cmm->SetLeftMargin(0.15);
0912
0913 gPad->SetLogy(1);
0914
0915 TLegend *lmm = new TLegend(0.35,0.62,0.9,0.90,"","NDC");
0916 lmm->SetBorderSize(1);
0917 for(int imaps =0;imaps<n_maps_layer+1;imaps++)
0918 {
0919 double nmapslo = -0.5 + 1.0 * imaps;
0920 double nmapshi = nmapslo + 1.0;
0921
0922 int binlo = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapslo);
0923 int binhi = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapshi);
0924
0925 TH1D *h = new TH1D("h","h",40,-1,4);
0926 hcorr_nfake_nmaps->ProjectionX("h",binlo,binhi);
0927
0928 h->SetMarkerStyle(20);
0929 h->SetMarkerSize(1.5);
0930 h->SetMarkerColor(fake_col[imaps]);
0931 h->SetMinimum(10);
0932 h->SetMaximum(5.0e6);
0933 char lab[500];
0934 sprintf(lab,"%i missed maps layers",imaps);
0935 lmm->AddEntry(h,lab,"p");
0936
0937 if(imaps == 0)
0938 {
0939 h->GetXaxis()->SetTitle("Mismatched hits");
0940 h->GetYaxis()->SetTitle("Tracks");
0941 h->DrawCopy("p");
0942
0943 }
0944 else
0945 h->DrawCopy("same p");
0946 }
0947
0948 lmm->Draw();
0949 }
0950
0951 TCanvas *cvtx = new TCanvas("cvtx","cvtx",4,4,800,600);
0952 TH1D *hzevt = 0;
0953 fin->GetObject("hzevt",hzevt);
0954 if(!hzevt)
0955 {
0956 cout << "Did not get hzevt" << endl;
0957 exit(1);
0958 }
0959 hzevt->Draw();
0960
0961 TCanvas *cvtx_res = new TCanvas("cvtx_res","cvtx_res",4,4,1200,800);
0962 cvtx_res->Divide(2,2);
0963
0964 TH1D *hzvtx_res = 0;
0965 fin->GetObject("hzvtx_res",hzvtx_res);
0966 if(!hzvtx_res)
0967 {
0968 cout << "Did not get hzvtx_res" << endl;
0969 }
0970
0971 if(hzvtx_res)
0972 {
0973 cvtx_res->cd(1);
0974
0975 hzvtx_res->SetNdivisions(505);
0976 hzvtx_res->Draw();
0977
0978 TF1 *fvtxres = new TF1("fvtxres","gaus");
0979
0980 }
0981
0982 TH1D *hxvtx_res = 0;
0983 fin->GetObject("hxvtx_res",hxvtx_res);
0984 if(!hxvtx_res)
0985 {
0986 cout << "Did not get hxvtx_res" << endl;
0987 }
0988
0989 if(hxvtx_res)
0990 {
0991 cvtx_res->cd(2);
0992
0993 hxvtx_res->SetNdivisions(505);
0994 hxvtx_res->Draw();
0995
0996 TF1 *fvtxres = new TF1("fvtxres","gaus");
0997
0998 }
0999
1000 TH1D *hyvtx_res = 0;
1001 fin->GetObject("hyvtx_res",hyvtx_res);
1002 if(!hyvtx_res)
1003 {
1004 cout << "Did not get hyvtx_res" << endl;
1005 }
1006
1007 if(hyvtx_res)
1008 {
1009 cvtx_res->cd(3);
1010
1011 hyvtx_res->SetNdivisions(505);
1012 hyvtx_res->Draw();
1013
1014 TF1 *fvtxres = new TF1("fvtxres","gaus");
1015
1016 }
1017
1018 TH1D *hdcavtx_res = 0;
1019 fin->GetObject("hdcavtx_res",hdcavtx_res);
1020 if(!hdcavtx_res)
1021 {
1022 cout << "Did not get hdcavtx_res" << endl;
1023 }
1024
1025 if(hdcavtx_res)
1026 {
1027 cvtx_res->cd(4);
1028
1029 hdcavtx_res->SetNdivisions(505);
1030 hdcavtx_res->Draw();
1031
1032 TF1 *fvtxres = new TF1("fvtxres","gaus");
1033
1034 }
1035
1036
1037 TCanvas *g4ntr = new TCanvas("g4ntr","g4ntr",5,10,1200,800);
1038 g4ntr->Divide(3,1);
1039 g4ntr->cd(1);
1040 TH2D *hg4ntrack = 0;
1041 fin->GetObject("hg4ntrack",hg4ntrack);
1042 if(hg4ntrack)
1043 {
1044 hg4ntrack->GetXaxis()->SetTitle(" Tracks/event from truth");
1045 hg4ntrack->GetYaxis()->SetTitle(" Tracks/event reconstructed");
1046 hg4ntrack->Draw();
1047 }
1048 else
1049 cout << "Did not get hgn4track" << endl;
1050
1051 g4ntr->cd(2);
1052 TH1D *hg4ntr = new TH1D("hg4ntr","hg4ntr",200, 0.0, 2000.0);
1053 hg4ntrack->ProjectionX("hg4ntr");
1054 hg4ntr->GetXaxis()->SetTitle(" Tracks/event from truth");
1055 hg4ntr->GetYaxis()->SetTitle(" Yield");
1056 hg4ntr->Draw();
1057
1058 g4ntr->cd(3);
1059 double cent[20];
1060 double nhits[20];
1061 nhits[0] = 0;
1062 double ytot = hg4ntr->Integral();
1063 for(int i=1;i<20;i++)
1064 {
1065 nhits[i] = i*120.0;
1066 int binhi = hg4ntr->FindBin(nhits[i]);
1067 double y = hg4ntr->Integral(1,binhi);
1068 cent[i] = y / ytot;
1069 cout << "nhits = " << nhits[i] << " binhi = " << binhi << " y = " << y << " ytot = " << ytot << " cent = " << cent[i] << endl;
1070 }
1071 TGraph *gy = new TGraph(20,nhits,cent);
1072 gy->GetXaxis()->SetTitle(" Tracks/event from truth");
1073 gy->GetYaxis()->SetTitle(" integral fraction of events");
1074 gy->Draw();
1075
1076 TF1 *fgy = new TF1("fgy","[0]*x+[1]*x*x+[2]*x*x*x",0,2300);
1077 gy->Fit(fgy,"R");
1078
1079
1080
1081 TFile *fout = new TFile("root_files/look_purity_out.root","recreate");
1082 gr_eff->Write();
1083 grdca2d->Write();
1084 grdcaZ->Write();
1085 grdpt->Write();
1086
1087 }