Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
trkeff.C
Go to the documentation of this file.
1 void trkeff(int pid)
2 {
3  gROOT->Reset();
4 
5  Double_t Red[5] = { 0.5,1.00,,0,0, 0 };
6  Double_t Green[5] = { 0.1, 0.3, 1., 0., 0.};
7  Double_t Blue[5] = { 0.1, 0.35, 0.65,0.8,1.0 };
8  Double_t Length[5] = { 0.00, 0.3, 0.5, 0.85, 1.00 };
9  Int_t nb=50;
10  TColor::CreateGradientColorTable(5,Length,Red,Green,Blue,nb);
11 
12  // If we are not already in the TRACKING directory, cd there
13  TDirectory *TRACKING = (TDirectory*)gROOT->FindObject("TRACKING");
14  if(TRACKING)TRACKING->cd();
15 
16  // Get pointer to trkeff tree
17  TTree *trkeff = (TTree*)gROOT->FindObject("trkeff");
18 
19  // Create histograms to calculate efficiency
20  TH2D *np = new TH2D("np","Numerator", 70,0,140,70,-0.05,6.95);
21  TH2D *dp = new TH2D("dp","Denominator",70,0,140,70,-0.05,6.95);
22  TH2D *npcut = new TH2D("npcut","numerator",70,0,140,70,-0.05,6.95);
23 
24  trkeff->Project("np", "F.pthrown.Mag():180./3.14159*F.pthrown.Theta()", "F.trktb.Nfdc>0 || F.trktb.Ncdc>0");
25  trkeff->Project("dp", "F.pthrown.Mag():180/3.14159*F.pthrown.Theta()", "");
26  trkeff->Project("npcut", "F.pthrown.Mag():180./3.14159*F.pthrown.Theta()", "(F.trktb.Nfdc>0 || F.trktb.Ncdc>0) && TMath::Prob(F.trktb.trk_chisq,F.trktb.trk_Ndof)>0.01");
27 
28  TH2D *effp = np->Clone("effp");
29 
30  TCanvas *c1 = new TCanvas("c1");
31  c1->SetTickx();
32  c1->SetTicky();
33  effp->Divide(dp);
34  effp->SetMarkerStyle(2);
35  effp->SetStats(0);
36  //gStyle->SetPalette(1,0);
37  effp->SetContour(nb);
38  if (pid==9)
39  effp->SetTitle("Pion reconstruction efficiency");
40  else{
41  effp->SetTitle("Proton reconstruction efficiency");
42  effp->SetAxisRange(0,88,"X");
43  effp->SetAxisRange(0,2,"Y");
44  }
45  effp->SetXTitle("#theta [degrees]");
46  effp->SetYTitle("p [GeV/c]");
47 
48  // Find the region where the efficiency drops through 50% and fit
49  // the resulting ridge with a polynomial.
50  TH2D *ridge = new TH2D("ridge","ridge", 70,0,140,70,-0.05,6.95);
51  for (unsigned int i=2;i<65;i++){
52  for (unsigned int j=2;j<70;j++){
53  double frac=effp->GetBinContent(i,j);
54  if (frac>0.0){
55  double min_weight=exp(-(0.8-0.5)*(0.8-0.5)/1);
56  double weight=exp(-(frac-0.5)*(frac-0.5)/1);
57  if (weight>min_weight){
58  ridge->SetBinContent(i,j,1000.*(weight-min_weight));
59 
60  }
61  }
62  }
63  }
64  ridge->Fit("pol3");
65  effp->Draw("colz");
66  double par[4];
67  TF1 *f1=ridge->GetFunction("pol3");
68  f1->GetParameters(par);
69  f1->SetRange(0.95,130.5);
70  f1->SetLineColor(5);
71  f1->Draw("same");
72 
73  double pmax=6.95;
74  if (pid==14) pmax=2.05;
75  double p1=par[0]+par[1]+par[2]+par[3];
76  double p2=par[0]+par[1]*130.+par[2]*130.*130.+par[3]*130.*130.*130.;
77  TLine *line1 = new TLine(1,p1,1,pmax);
78  line1->SetLineColor(5);
79  line1->SetLineWidth(2);
80  line1->Draw();
81  if (pid==9){
82  TLine *line2 = new TLine(130,p2,130,pmax);
83  line2->SetLineColor(5);
84  line2->SetLineWidth(2);
85  line2->Draw();
86  }
87  if (pid==9) c1->SaveAs("pion_efficiency.png");
88  else c1->SaveAs("proton_efficiency.png");
89 
90 
91  TH2D *effpcut = npcut->Clone("effpcut");
92 
93  TCanvas *c2 = new TCanvas("c2");
94  c2->SetTickx();
95  c2->SetTicky();
96  effpcut->Divide(dp);
97  effpcut->SetMarkerStyle(2);
98  effpcut->SetStats(0);
99 
100 
101  effpcut->SetContour(nb);
102  if (pid==9)
103  effpcut->SetTitle("Pion reconstruction efficiency, P(#chi^2)>0.01");
104  else{
105  effpcut->SetTitle("Proton reconstruction efficiency, P(#chi^2)>0.01");
106  effpcut->SetAxisRange(0,88,"X");
107  effpcut->SetAxisRange(0,2,"Y");
108  }
109  effpcut->SetXTitle("#theta [degrees]");
110  effpcut->SetYTitle("p [GeV/c]");
111  effpcut->Draw("colz");
112  f1->Draw("same");
113  line1->Draw();
114  if (pid==9) line2->Draw();
115 
116  if (pid==9) c2->SaveAs("pion_efficiency_with_chi2_cut.png");
117  else c2->SaveAs("proton_efficiency_with_chi2_cut.png");
118 
119  if (pid==9){
120  ofstream ofile("pion.txt");
121  ofile << "<html>" <<endl;
122  ofile<<"<h1>Single track reconstruction</h1>" <<endl;
123  ofile<<"<p> Yellow lines indicate recommended fiducial cuts </p>"
124  <<endl;
125  ofile<<"<h2>Pion events</h6>" <<endl;
126  }
127  else{
128  ofstream ofile("proton.txt");
129  ofile<<"<h2>Proton events</h6>" <<endl;
130  }
131  double eff = np->GetEntries()/dp->GetEntries();
132  ofile <<"<p>Overall efficiency = "<<eff<<"<br/>"<<endl;
133  double effcut = npcut->GetEntries()/dp->GetEntries();
134  ofile <<" Efficiency (Prob>0.01) = "<<effcut<<"<br/>" << endl;
135  ofile << "Momentum cut: p=" << par[0] << "&ensp;+&ensp;"
136  << par[1] << "&thinsp;&theta;&ensp;+&ensp;"
137  << par[2] << "&thinsp;&theta;<sup>2</sup>&ensp;+&ensp;" << par[3]
138  <<"&thinsp;&theta;<sup>3</sup></p>" <<endl;
139  if (pid==9){
140  ofile <<"<img src=\"pion_efficiency.png\">" <<endl;
141  ofile <<"<img src=\"pion_efficiency_with_chi2_cut.png\">" <<endl;
142  }
143  else{
144  ofile <<"<img src=\"proton_efficiency.png\">" <<endl;
145  ofile <<"<img src=\"proton_efficiency_with_chi2_cut.png\">" <<endl;
146  ofile << "</html>" <<endl;
147  }
148 
149  ofile.close();
150 }
trkeff
void trkeff(int pid)
Definition: trkeff.C:1
p1
double p1
Definition: JEventProcessor_PS_E_calib.cc:39
c1
Double_t c1[2][NMODULES]
Definition: tw_corr.C:68
f1
TFile f1(fnam)
c2
Double_t c2[2][NMODULES]
Definition: tw_corr.C:69
dp
int * dp
Definition: hddm_t.c:75
p2
double p2
Definition: JEventProcessor_PS_E_calib.cc:40