Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plotpirej.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plotpirej.C
1 void plotpirej()
2 {
3 gStyle->SetOptStat(0);
4 
5 const int nbins = 5;
6 double bins[nbins+1];
7 for(int i=0; i<=nbins; i++) { bins[i] = 2. + 2.*double(i); }
8 //cout << "bin edges:" << endl;
9 //for(int i=0; i<=nbins; i++) { cout << bins[i] << endl; }
10 
11 int neopbins = 100;
12 char hname[999];
13 TH1D* heop_pi[nbins+1];
14 TH1D* heop_e[nbins+1];
15 for(int i=0; i<=nbins; i++) {
16  sprintf(hname,"heop_pi_%d",i);
17  heop_pi[i] = new TH1D(hname,hname,neopbins,0.,2.);
18  heop_pi[i]->Sumw2();
19  sprintf(hname,"heop_e_%d",i);
20  heop_e[i] = new TH1D(hname,hname,neopbins,0.,2.);
21  heop_e[i]->Sumw2();
22 }
23  TH1D* hhoe_e = new TH1D("hhoe_e","hhoe_e",150,0.,1.5);
24  TH1D* hhoe_pi = new TH1D("hhoe_pi","hhoe_pi",150,0.,1.5);
25  TH2D* hhhoe_e = new TH2D("hhhoe_e","hhhoe_e",100,0.,2.0,100,0.,2.0);
26  TH2D* hhhoe_pi = new TH2D("hhhoe_pi","hhhoe_pi",100,0.,2.0,100,0.,2.0);
27 
28  TH1D* hprob_e = new TH1D("hprob_e","hprob_e",100,0.,1.);
29  TH1D* hprob_pi = new TH1D("hprob_pi","hprob_pi",100,0.,1.);
30  TH1D* hchi2_e = new TH1D("hchi2_e","hchi2_e",100,0.,10.);
31  TH1D* hchi2_pi = new TH1D("hchi2_pi","hchi2_pi",100,0.,10.);
32 
33 TChain* ntp_track_pi = new TChain("ntppid");
34  ntp_track_pi->AddFile("pions0.root");
35  ntp_track_pi->AddFile("pions1.root");
36  ntp_track_pi->AddFile("pions2.root");
37  ntp_track_pi->AddFile("pions4.root");
38 
39 //string infname_piminus;
40 //ifstream ifs_piminus;
41 //ifs_piminus.open("piminuslist.txt");
42 //if ( ifs_piminus.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs_piminus.close(); return; }
43 //int tmpcount=0;
44 //while(!ifs_piminus.eof()) {
45 // ifs_piminus >> infname_piminus;
46 // ntp_track_pi->AddFile(infname_piminus.c_str());
47 // tmpcount++;
48 // //if(tmpcount>20) break;
49 //}
50 //ifs_piminus.close();
51 
52 TChain* ntp_track_e = new TChain("ntppid");
53  ntp_track_e->AddFile("electrons0.root");
54  ntp_track_e->AddFile("electrons1.root");
55  ntp_track_e->AddFile("electrons2.root");
56  ntp_track_e->AddFile("electrons3.root");
57 
58 //string infname_e;
59 //ifstream ifs_e;
60 //ifs_e.open("electronlist.txt");
61 //if ( ifs_e.fail() ) { cout << "FAIL TO READ INPUT FILE 2" << endl; ifs_e.close(); return; }
62 //tmpcount=0;
63 //while(!ifs_e.eof()) {
64 // ifs_e >> infname_e;
65 // ntp_track_e->AddFile(infname_e.c_str());
66 // tmpcount++;
67 // //if(tmpcount>20) break;
68 //}
69 //ifs_e.close();
70 
71 cout << "number of entries for pions = " << ntp_track_pi->GetEntries() << endl;
72 cout << "number of entries for electrons = " << ntp_track_e->GetEntries() << endl;
73 
74 float px,py,pz,pt,quality,nmvtx,ntpc,eta,chi2,ndf,mom;
75 float cemc_e,cemc_e3x3,cemc_dphi,cemc_deta,cemc_ecore,cemc_prob,cemc_chi2;
76 float dca2d,hcalin_e3x3;
77 
78 ntp_track_pi->SetBranchAddress("nmvtx",&nmvtx);
79 ntp_track_pi->SetBranchAddress("ntpc",&ntpc);
80 ntp_track_pi->SetBranchAddress("pt",&pt);
81 ntp_track_pi->SetBranchAddress("eta",&eta);
82 ntp_track_pi->SetBranchAddress("mom",&mom);
83 ntp_track_pi->SetBranchAddress("quality",&quality);
84 ntp_track_pi->SetBranchAddress("dca2d",&dca2d);
85 ntp_track_pi->SetBranchAddress("cemc_ecore",&cemc_ecore);
86 ntp_track_pi->SetBranchAddress("cemc_prob",&cemc_prob);
87 ntp_track_pi->SetBranchAddress("cemc_chi2",&cemc_chi2);
88 ntp_track_pi->SetBranchAddress("cemc_e3x3",&cemc_e3x3);
89 ntp_track_pi->SetBranchAddress("cemc_dphi",&cemc_dphi);
90 ntp_track_pi->SetBranchAddress("cemc_deta",&cemc_deta);
91 ntp_track_pi->SetBranchAddress("hcalin_e3x3",&hcalin_e3x3);
92 
93  for(int j=0; j<ntp_track_pi->GetEntries(); j++) {
94  ntp_track_pi->GetEvent(j);
95  if(quality>5) continue;
96  if(nmvtx<1) continue;
97  if(ntpc<20) continue;
98  if(mom<=0.) continue;
99  //if(fabs(cemc_dphi)>0.015) continue;
100  //if(fabs(cemc_deta)>0.015) continue;
101  if(cemc_ecore==0.) continue;
102  hhoe_pi->Fill(hcalin_e3x3/cemc_ecore);
103  float eop = cemc_ecore/mom;
104  hhhoe_pi->Fill(eop,hcalin_e3x3/cemc_ecore);
105  if(eop>=2.0) eop = 1.999;
106  heop_pi[nbins]->Fill(eop);
107  for(int j=0; j<nbins; j++) { if(pt>bins[j] && pt<bins[j+1]) { heop_pi[j]->Fill(eop); } }
108  hprob_pi->Fill(cemc_prob);
109  hchi2_pi->Fill(cemc_chi2);
110  }
111 
112 //----------------------------------------------------------------------------------------------
113 
114 ntp_track_e->SetBranchAddress("nmvtx",&nmvtx);
115 ntp_track_e->SetBranchAddress("ntpc",&ntpc);
116 ntp_track_e->SetBranchAddress("pt",&pt);
117 ntp_track_e->SetBranchAddress("eta",&eta);
118 ntp_track_e->SetBranchAddress("mom",&mom);
119 ntp_track_e->SetBranchAddress("quality",&quality);
120 ntp_track_e->SetBranchAddress("dca2d",&dca2d);
121 ntp_track_e->SetBranchAddress("cemc_ecore",&cemc_ecore);
122 ntp_track_e->SetBranchAddress("cemc_prob",&cemc_prob);
123 ntp_track_e->SetBranchAddress("cemc_chi2",&cemc_chi2);
124 ntp_track_e->SetBranchAddress("cemc_e3x3",&cemc_e3x3);
125 ntp_track_e->SetBranchAddress("cemc_dphi",&cemc_dphi);
126 ntp_track_e->SetBranchAddress("cemc_deta",&cemc_deta);
127 ntp_track_e->SetBranchAddress("hcalin_e3x3",&hcalin_e3x3);
128 
129  for(int j=0; j<ntp_track_e->GetEntries(); j++) {
130  ntp_track_e->GetEvent(j);
131  if(quality>5) continue;
132  if(nmvtx<1) continue;
133  if(ntpc<20) continue;
134  if(mom<=0.) continue;
135  //if(fabs(cemc_dphi)>0.015) continue;
136  //if(fabs(cemc_deta)>0.015) continue;
137  if(cemc_ecore==0.) continue;
138  hhoe_e->Fill(hcalin_e3x3/cemc_ecore);
139  float eop = cemc_ecore/mom;
140  hhhoe_e->Fill(eop,hcalin_e3x3/cemc_ecore);
141  if(eop>=2.0) eop = 1.999;
142  heop_e[nbins]->Fill(eop);
143  for(int j=0; j<nbins; j++) { if(pt>bins[j] && pt<bins[j+1]) { heop_e[j]->Fill(eop); } }
144  hprob_e->Fill(cemc_prob);
145  hchi2_e->Fill(cemc_chi2);
146  }
147 
148 //--------------------------------------------------------------------------
149 
150 TCanvas* choe = new TCanvas("choe","hcalin/cemc",50,50,600,600);
151 choe->SetLogy();
152 hhoe_e->SetLineColor(kRed);
153 hhoe_pi->Draw();
154 hhoe_e->Draw("same");
155 
156 TCanvas* choe2 = new TCanvas("choe2","hcalin/cemc vs e/p",100,100,600,600);
157 hhhoe_e->SetLineColor(kRed);
158 hhhoe_pi->Draw("box");
159 hhhoe_e->Draw("boxsame");
160 
161 TCanvas* cprob = new TCanvas("cprob","prob",150,150,600,500);
162 hprob_e->SetLineColor(kRed);
163 hprob_pi->Draw();
164 hprob_e->Draw("same");
165 
166 TCanvas* cchi2 = new TCanvas("cchi2","cemc chi2",200,200,600,500);
167 hchi2_e->SetLineColor(kRed);
168 hchi2_pi->Draw();
169 hchi2_e->Draw("same");
170 
171 TCanvas* ceop = new TCanvas("ceop","E/p",250,250,600,500);
172 ceop->SetLogy();
173 heop_pi[nbins]->Draw();
174 heop_e[nbins]->SetLineColor(kRed);
175 heop_e[nbins]->Draw("same");
176 
177 double vsume[nbins];
178 double vsumpi[nbins];
179 for(int i=0; i<nbins; i++) {vsume[i]=0.; vsumpi[i]=0.;}
180 double sume = 0.;
181 double sumpi = 0.;
182 int uplim = heop_e[nbins]->GetNbinsX();
183 for(int i=1; i<=uplim; i++) {
184  sume += heop_e[nbins]->GetBinContent(i);
185  sumpi += heop_pi[nbins]->GetBinContent(i);
186 }
187 cout << "total number of pions = " << sumpi << endl;
188 cout << "total number of electrons = " << sume << endl;
189 
190  for(int j=0; j<nbins; j++) {
191  for(int i=1; i<=uplim; i++) {
192  vsume[j] += heop_e[j]->GetBinContent(i);
193  vsumpi[j] += heop_pi[j]->GetBinContent(i);
194  }
195  cout << "bin " << j << " " << vsume[j] << " " << vsumpi[j] << endl;
196  }
197 
198 double eleff = 0.9;
199 
200 double vsumecut[nbins];
201 int vbincut[nbins];
202 for(int i = 0; i<nbins; i++) { vsumecut[i] = 0.; vbincut[i] = 0; }
203 double sumecut = 0.;
204 int bincut = 0;
205 for(int i=uplim; i>=1; i--) {
206  sumecut += heop_e[nbins]->GetBinContent(i);
207  if(sumecut>sume*eleff) {
208  cout << "90% cut: " << i << " " << heop_e[nbins]->GetBinLowEdge(i) << endl;
209  bincut = i;
210  break;
211  }
212 }
213  for(int j=0; j<nbins; j++) {
214  for(int i=uplim; i>=1; i--) {
215  vsumecut[j] += heop_e[j]->GetBinContent(i);
216  if(vsumecut[j]>vsume[j]*eleff) {
217  cout << "bin " << j << " 90% cut: " << i << " " << heop_e[j]->GetBinLowEdge(i) << endl;
218  vbincut[j] = i;
219  break;
220  }
221  }
222  }
223 
224 double vsumpicut[nbins];
225 for(int i=0; i<nbins; i++) {vsumpicut[i]=0.;}
226 double sumpicut = 0.;
227 for(int i=bincut; i<=heop_pi[nbins]->GetNbinsX(); i++) {
228  sumpicut += heop_pi[nbins]->GetBinContent(i);
229 }
230 cout << "number of misidentified pions = " << sumpicut << endl;
231 cout << "rejection factor = " << sumpi/sumpicut << " +- " << (sumpi/sumpicut)*(1./sqrt(sumpicut)) << endl;
232 
233  for(int j=0; j<nbins; j++) {
234  for(int i=vbincut[j]; i<=heop_pi[j]->GetNbinsX(); i++) {
235  vsumpicut[j] += heop_pi[j]->GetBinContent(i);
236  }
237  cout << "bin " << j << " number of misidentified pions = " << vsumpicut[j] << endl;
238  }
239 
240 double rejpi[nbins];
241 double errejpi[nbins];
242 double xx[nbins];
243 for(int i=0; i<nbins; i++) {
244  rejpi[i] = vsumpicut[i]/vsumpi[i];
245  errejpi[i] = rejpi[i]/sqrt(vsumpicut[i]);
246  xx[i] = (bins[i]+bins[i+1])/2.;
247 }
248 
249 TCanvas* c4 = new TCanvas("c4","pion rejection",120,120,600,500);
250 c4->SetLogy();
251 double lowlim = 0.001;
252 if(eleff==0.7) lowlim = 0.002;
253 TH2D* hh = new TH2D("hh","",10,1.,13.,10,lowlim,0.1);
254 hh->SetTitle(";Transverse momentum, [GeV] ; Inverse pion rejection factor");
255 hh->Draw();
256 
257 TGraphErrors* gr1 = new TGraphErrors(nbins,xx,rejpi,0,errejpi);
258 gr1->SetMarkerStyle(20);
259 gr1->SetMarkerColor(kBlack);
260 gr1->SetMarkerSize(1.5);
261 gr1->SetLineColor(kBlack);
262 gr1->SetLineWidth(2);
263 gr1->Draw("p");
264 
265 }
266