StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plCalcEffic.C
1 
2 ofstream spreadsheet("./effic.csv");
3 
4 void plCalcEffic(int charge = 2){
5  system("mkdir -p plots/");
6 
7  if(charge == 2){
8  plCalcEfficX(0);
9  plCalcEfficX(1);
10  }
11  else plCalcEfficX(charge);
12  return;
13 }
14 
15 void plCalcEfficX(int charge = 0, char* iPath="/star/u/stevens4/wAnalysis/efficXsec/outEmb/gainUp2/"){
16 
17  char* sign;
18  if(charge==0) sign="+";
19  else if(charge==1) sign="-";
20 
21  char* core0; //these have no vertex reweighting
22  if(charge == 1) core0="Wminus";
23  else if(charge == 0) core0="Wplus";
24 
25  ofstream latex(Form("./%s.txt",core0));
26 
27  //gStyle->SetOptFit(1);
28  //gStyle->SetOptStat(0);
29  gStyle->SetPalette(1);
30 
31  //load file
32  TString fullInpName=iPath; fullInpName+=core0;
33  fullInpName+=".wana.hist.root";
34  fd=new TFile(fullInpName);
35  if(!fd->IsOpen()) {
36  printf("ERROR: input histo file not found, quit\n",fullInpName.Data());
37  return;
38  } else {
39  printf("Opened: %s\n",fullInpName.Data());
40  }
41 
42  //use rebin 1 for effic # calc, but rebin N for plot
43  int etrebin=4;
44  int etaRebin=50;
45  int zRebin=4;
46  int phiModuleRebin=6;
47  int phiRebin=1;
48  int zdcRebin=5;
49 
50  //get total efficiency and error
51  float etplot[5]={1.0,1.0,1.0,1.0,1.0};
52  float etplotErr[5]={1.0,1.0,1.0,1.0,0.0};
53  float zdcplot[5]={2.0,2.0,2.0,2.0,2.0};
54  float zdcplotErr[5]={1.0,1.0,1.0,1.0,0.0};
55 
56  //trigger efficiency
57  TH1F *hTrigET=doEfficiency(fd,"MCeleETall","MCeleETtrig",Form("W%s trigger efficiency",sign),etrebin,etplot[0],etplotErr[0]);
58  TH1F *hTrigEta=doEfficiency(fd,"MCeleEtaAll","MCeleEtaTrig",Form("W%s trigger efficiency",sign),etaRebin,0);
59  TH1F *hTrigDetEta=doEfficiency(fd,"MCeleDetEtaAll","MCeleDetEtaTrig",Form("W%s trigger efficiency",sign),etaRebin,0);
60  TH1F *hTrigZvert=doEfficiency(fd,"MCeleZvertAll","MCeleZvertTrig",Form("W%s trigger efficiency",sign),zRebin,0);
61  TH1F *hTrigPhi=doEfficiency(fd,"MCelePhiModulePairAll","MCelePhiModulePairTrig",Form("W%s trigger efficiency",sign),phiModuleRebin,0);
62  TH1F *hTrigZDC=doEfficiency(fd,"MCeleZdcAll","MCeleZdcTrig",Form("W%s trigger efficiency",sign),zdcRebin,zdcplot[0],zdcplotErr[0]);
63  //vertex efficiency
64  TH1F *hVertET=doEfficiency(fd,"MCeleETtrig","MCeleETvert",Form("W%s vertex efficiency",sign),etrebin,etplot[1],etplotErr[1]);
65  TH1F *hVertEta=doEfficiency(fd,"MCeleEtaTrig","MCeleEtaVert",Form("W%s vertex efficiency",sign),etaRebin,0);
66  TH1F *hVertZvert=doEfficiency(fd,"MCeleZvertTrig","MCeleZvertVert",Form("W%s vertex efficiency",sign),zRebin,0);
67  TH1F *hVertPhi=doEfficiency(fd,"MCelePhiTrig","MCelePhiVert",Form("W%s vertex efficiency",sign),phiRebin,0);
68  TH1F *hVertZDC=doEfficiency(fd,"MCeleZdcTrig","MCeleZdcVert",Form("W%s vertex efficiency",sign),zdcRebin,zdcplot[1],zdcplotErr[1]);
69  //tracking efficiency
70  TH1F *hTrackET=doEfficiency(fd,"MCeleETvert","MCeleETTrack",Form("W%s tracking efficiency",sign),etrebin,etplot[2],etplotErr[2]);
71  TH1F *hTrackEta=doEfficiency(fd,"MCeleEtaVert","MCeleEtaTrack",Form("W%s tracking efficiency",sign),etaRebin,0);
72  TH1F *hTrackZvert=doEfficiency(fd,"MCeleZvertVert","MCeleZvertTrack",Form("W%s tracking efficiency",sign),zRebin,0);
73  TH1F *hTrackPhi=doEfficiency(fd,"MCelePhiVert","MCelePhiTrack",Form("W%s tracking efficiency",sign),phiRebin,0);
74  TH1F *hTrackZDC=doEfficiency(fd,"MCeleZdcVert","MCeleZdcTrack",Form("W%s tracking efficiency",sign),zdcRebin,zdcplot[2],zdcplotErr[2]);
75  //algo efficiency
76  TH1F *hRecoET=doEfficiency(fd,"MCeleETTrack","MCeleETreco",Form("W%s algo efficiency",sign),etrebin,etplot[3],etplotErr[3]);
77  TH1F *hRecoEta=doEfficiency(fd,"MCeleEtaTrack","MCeleEtaReco",Form("W%s algo efficiency",sign),etaRebin,0);
78  TH1F *hRecoZvert=doEfficiency(fd,"MCeleZvertTrack","MCeleZvertReco",Form("W%s algo efficiency",sign),zRebin,0);
79  TH1F *hRecoPhi=doEfficiency(fd,"MCelePhiTrack","MCelePhiReco",Form("W%s algo efficiency",sign),phiRebin,0);
80  TH1F *hRecoZDC=doEfficiency(fd,"MCeleZdcTrack","MCeleZdcReco",Form("W%s algo efficiency",sign),zdcRebin,zdcplot[3],zdcplotErr[3]);
81 
82  //total efficiency for overall stat uncertainty
83  TH1F *hTotET=doEfficiency(fd,"MCeleETall","MCeleETreco",Form("W%s total efficiency",sign),etrebin,etplot[4],etplotErr[4]);
84  TH1F *hTotZDC=doEfficiency(fd,"MCeleZdcAll","MCeleZdcReco",Form("W%s total efficiency",sign),zdcRebin,zdcplot[4],zdcplotErr[4]);
85  float totEffic=etplot[0]*etplot[1]*etplot[2]*etplot[3];
86  //cout<<totEffic<<endl;
87  cout<<"Total Efficiency = "<<etplot[4]<<" $\\pm$ "<<etplotErr[4]<<endl;
88  //spreadsheet<<"Track Effic = "<<etplot[2]<<" +/- "<<etplotErr[2]<<endl;
89 
90  spreadsheet<<"W"<<sign<<" Total Efficiency"<<endl;
91  for(int j=0; j<6; j++){
92  spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hTotET->GetBinContent(7+j))<<","<<Form("%.4f",hTotET->GetBinError(7+j))<<endl;
93  }
94 
95  //set axis range for ET plot
96  int etAxisMin=0; int etAxisMax=70;
97 
98  //*************** Algo Efficiency plots ***************//
99  cA=new TCanvas(Form("W%s algo effic",sign),"algo effic",800,600);
100  cA->Divide(2,2);
101  cA->cd(1);
102  hRecoET->Draw();
103  hRecoET->SetAxisRange(etAxisMin,etAxisMax);
104  cA->cd(2);
105  hRecoEta->Draw();
106  cA->cd(3);
107  hRecoZvert->Draw();
108  cA->cd(4);
109  hRecoZDC->Draw();
110  //hRecoPhi->Draw();
111  cA->Print(Form("plots/W%salgoEffic.eps",sign));
112  cA->Print(Form("plots/W%salgoEffic.png",sign));
113 
114  //algo efficiency numbers
115  spreadsheet<<"W"<<sign<<" Algo Efficiency"<<endl;
116  latex<<"W"<<sign<<" Algo Efficiency"<<endl;
117  for(int j=0; j<6; j++){
118  //semi-latex format
119  if(j<6) latex<<25+j*4<<"$<E_T<$"<<29+j*4<<" & "<<Form("%.4f",hRecoET->GetBinContent(7+j))<<" $\\pm$ "<<Form("%.4f",hRecoET->GetBinError(7+j))<<" $\\pm$ & \\\\"<<endl;
120 
121  //for spreadsheet
122  spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hRecoET->GetBinContent(7+j))<<","<<Form("%.4f",hRecoET->GetBinError(7+j))<<endl;
123  }
124 
125  //*************** Tracking Efficiency plots ***********//
126  cT=new TCanvas(Form("W%s track effic",sign),"track effic",800,600);
127  cT->Divide(2,2);
128  cT->cd(1);
129  hTrackET->Draw();
130  hTrackET->SetAxisRange(etAxisMin,etAxisMax);
131  cT->cd(2);
132  hTrackEta->Draw();
133  cT->cd(3);
134  //hTrackZvert->Draw();
135  hTrackZDC->SetMinimum(0.6);
136  hTrackZDC->Draw();
137  cT->cd(4);
138  hTrackPhi->Draw();
139  cT->Print(Form("plots/W%strackEffic.eps",sign));
140  cT->Print(Form("plots/W%strackEffic.png",sign));
141 
142  //*************** Vertex Efficiency plots *************//
143  //Draw trigger efficiency
144  cV=new TCanvas(Form("W%s vertex effic",sign),"vertex effic",800,600);
145  cV->Divide(2,2);
146  cV->cd(1);
147  hVertET->Draw();
148  hVertET->SetAxisRange(etAxisMin,etAxisMax);
149  cV->cd(2);
150  hVertEta->Draw();
151  cV->cd(3);
152  hVertZvert->Draw();
153  cV->cd(4);
154  hVertZDC->SetMinimum(0.6);
155  hVertZDC->Draw();
156  cV->Print(Form("plots/W%svertEffic.eps",sign));
157  cV->Print(Form("plots/W%svertEffic.png",sign));
158 
159  //*************** Trigger Efficiency plots ************//
160  //draw trigger efficiency to show why W- isn't const w/ ET
161  c2D=new TCanvas(Form("W%s trigger effic ET bins",sign),"trigger effic ET bins",700,500);
162  j1=(TH2F*)fd->Get("MCeleEta_ptPreTrig"); //j1->Rebin2D(2,5);
163  TH1D* h25[10];
164  c2D->Divide(2,2);
165  spreadsheet<<"W"<<sign<<" Trigger Efficiency"<<endl;
166  latex<<"W"<<sign<<" Trigger Efficiency"<<endl;
167  for(int j=0; j<6; j++){
168  h25[j] = j1->ProjectionX(Form("pt%d_%d",25+j*4,29+j*4),26+j*4,29+j*4);
169  h25[j]->SetTitle(Form("Lepton detector #eta (from Geant): PT=[%d,%d]",25+j*4,29+j*4));
170  h25[j]->Rebin(2);
171 
172  //count frac. of events with |eta|<1 and effic for each bin
173  float accepted = h25[j]->Integral(26,125);
174  float total = h25[j]->Integral();
175  //cout<<"ET=["<<25+j*4<<","<<29+j*4<<"] efficiency = "<<hTrigET->GetBinContent(7+j)<<" ; fraction with |eta| < 1 : "<<accepted<<" "<<total<<endl;
176 
177  //semi-latex format
178  if(j<6) latex<<25+j*4<<"$<E_T<$"<<29+j*4<<" & "<<Form("%.4f",hTrigET->GetBinContent(7+j))<<" $\\pm$ "<<Form("%.4f",hTrigET->GetBinError(7+j))<<" $\\pm$ & & "<<Form("%.4f",accepted/total)<<" & \\\\"<<endl;
179  spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hTrigET->GetBinContent(7+j))<<","<<Form("%.4f",hTrigET->GetBinError(7+j))<<endl;
180 
181  //add lines for detector cutoff at eta=1
182  h25[j]->Rebin(5);
183  Lx=h25[j]->GetListOfFunctions();
184  int max=h25[j]->GetMaximum();
185  ln1=new TLine(-1,0,-1,max);
186  ln2=new TLine(1,0,1,max);
187  ln1->SetLineWidth(2); ln2->SetLineWidth(2);
188  ln1->SetLineColor(kRed); ln1->SetLineStyle(2);
189  ln2->SetLineColor(kRed); ln2->SetLineStyle(2);
190  Lx->Add(ln1); Lx->Add(ln2);
191  if(j<4){ //only plot first four ET bins
192  c2D->cd(j+1);
193  h25[j]->Draw("h");
194  }
195  }
196  c2D->Print(Form("plots/W%strigEfficNonConst.eps",sign));
197  c2D->Print(Form("plots/W%strigEfficNonConst.png",sign));
198 
199  //Draw trigger efficiency
200  c=new TCanvas(Form("W%s trigger effic",sign),"trigger effic",800,600);
201  c->Divide(2,2);
202  c->cd(1);
203  hTrigET->Draw();
204  hTrigET->SetAxisRange(etAxisMin,etAxisMax);
205  c->cd(2);
206  hTrigEta->Draw();
207  c->cd(3);
208  hTrigDetEta->Draw();
209  c->cd(4);
210  hTrigPhi->Draw();
211  c->Print(Form("plots/W%strigEffic.eps",sign));
212  c->Print(Form("plots/W%strigEffic.png",sign));
213 
214 
215  //other efficiency related histos
216  TH2F * g0=(TH2F * )fd->Get("muBdist1"); assert(g0);
217  TH2F * g1=(TH2F * )fd->Get("muBclET24R_ET"); assert(g1);
218  TH2F * g2=(TH2F * )fd->Get("muBclEjetE2D_ET"); assert(g2);
219  TH2F * g3=(TH2F * )fd->Get("musPtBalance_clust"); assert(g3);
220  cAlgo2=new TCanvas(Form("W%s algo effic 2 ",sign),"algo effic 2",800,600);
221  cAlgo2->Divide(2,2);
222  cAlgo2->cd(1);
223  g0->Draw("colz");
224  g0->GetXaxis()->SetRange(0,70);
225  cAlgo2->cd(2);
226  g1->Draw("colz");
227  g1->GetXaxis()->SetRange(0,70);
228  cAlgo2->cd(3);
229  g2->Draw("colz");
230  g2->GetXaxis()->SetRange(0,70);
231  cAlgo2->cd(4);
232  g3->Draw("colz");
233  g3->GetXaxis()->SetRange(0,70);
234  cAlgo2->Print(Form("plots/W%salgoEffic2.eps",sign));
235  cAlgo2->Print(Form("plots/W%salgoEffic2.png",sign));
236 
237  //cout<<g2->GetYaxis()->FindBin(0.88)<<endl;
238  //cout<<g2->Integral(26,30,0,1000)<<" "<<g2->Integral(26,30,1,g2->GetYaxis()->FindBin(0.88))<<endl;
239  //cout<<g2->Integral(26,1000,1,1000)<<" "<<g2->Integral(26,1000,1,g2->GetYaxis()->FindBin(0.88))<<endl;
240 
241  return;
242 }
243 
244 //calculate efficiencies
245 TH1F* doEfficiency(TFile *fd,char* name0, char* name1, char *tit, int reb=1, float &etplot, float &etplotErr=0){
246 
247  TH1F * h0=(TH1F * )fd->Get(name0); assert(h0);
248  TH1F * h1=(TH1F * )fd->Get(name1); assert(h1);
249 
250  ha=(TH1F*) h0->Clone(); assert(ha);
251  hb=(TH1F*) h1->Clone(); assert(hb);
252 
253  ha->Rebin(reb);
254  hb->Rebin(reb);
255 
256  hc=(TH1F*) hb->Clone(); assert(hc);
257  hc->SetTitle(tit);
258  hc->Reset();
259 
260  float num=0; float den=0;
261  float matchErr2=0; float thrownErr2=0;
262  int nb=hb->GetNbinsX();
263  int startBin,endBin;
264  if(etplot==1.0) {
265  startBin=hb->GetXaxis()->FindBin(25.);
266  //endBin=hb->GetXaxis()->FindBin(48.0);
267  //startBin=hb->GetXaxis()->FindBin(50.);
268  endBin=hb->GetXaxis()->FindBin(100.0);
269  }
270  if(etplot==2.0) {
271  endBin=hb->GetXaxis()->FindBin(200000.);
272  startBin=hb->GetXaxis()->FindBin(0.);
273  }
274 
275  //cout<<startBin<<endl;
276  //cout<<hb->GetNbinsX()<<" "<<ha->GetNbinsX()<<endl;
277  assert(nb==ha->GetNbinsX());
278  for(int i=1; i<=nb; i++) {
279  float n0=ha->GetBinContent(i);//n0: # of thrown tracks
280  float n1=hb->GetBinContent(i);//n1: # of matched tracks
281  float e0=ha->GetBinError(i);//e0: error on # of thrown tracks
282  float e1=hb->GetBinError(i);//e1: error on # of matched tracks
283  if(n0==0) continue; // no thown tracks found
284  float eff=(float) n1/n0;
285 
286  //old method with un-weighted histos
287  //float x=n1*(n0-n1);
288  //if(n1==0 || n1==n0) x=n0;
289  //float errEff=sqrt(x/n0/n0/n0);
290 
291  float errEff=0;
292  //method for using weighted histos (must use Sumw2)
293  if(n1==n0) {//problem in effic=1 limit, statErr=0
294  //cout<<name1<<" bin="<<i<<" problem with error"<<endl;
295  errEff=1./n0; //just something big to see its low stats
296  }
297  else errEff=sqrt(e1*e1*(n0-n1)*(n0-n1)+(e0-e1)*(e0-e1)*n1*n1)/n0/n0;
298 
299  //printf("name=%s: bin#= %d, n0= %f, n1= %f, e0=%f, e1=%f, eff= %f, errEff=+/- %f\n",name0,i,n0,n1,e0,e1,eff,errEff);
300  hc->SetBinContent(i,eff);
301  hc->SetBinError(i,errEff);
302  if(i>=startBin && i<=endBin){//define ET range for efficTot
303  num+=n1;
304  den+=n0;
305  matchErr2+=e1*e1;
306  thrownErr2+=e0*e0;
307  }
308  }
309 
310  if(etplot==1.0){
311  //old method with un-weighted histos
312  //cout<<name1<<"="<<(float) num/den<<"+/-"<<sqrt((float)num*(den-num)/den/den/den)<<endl;
313 
314  //method for using weighted histos (must use Sumw2)
315  if(etplotErr>0.0) cout<<name1<<"="<< num/den <<" $\\pm$ "<< sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den <<endl;
316  etplot=num/den;
317  etplotErr=sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den;
318  }
319  if(etplot==2.0){
320  //method for using weighted histos (must use Sumw2)
321  //cout<<name1<<"="<< num/den <<" $\\pm$ "<< sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den <<endl;
322 
323  etplot=num/den;
324  etplotErr=sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den;
325  }
326  hc->SetMinimum(0.0);
327  hc->SetMaximum(1.1);
328 
329  return hc;
330 }
Generated by   doxygen 1.8.5