StRoot: StSpinPool/StWalgoB2009/macros/plCalcEffic.C Source File

00001 
00002 ofstream spreadsheet("./effic.csv");
00003 
00004 void plCalcEffic(int charge = 2){
00005   system("mkdir -p plots/");
00006 
00007   if(charge == 2){
00008     plCalcEfficX(0);
00009     plCalcEfficX(1);
00010   }
00011   else plCalcEfficX(charge);
00012   return;
00013 }
00014 
00015 void plCalcEfficX(int charge = 0, char* iPath="/star/u/stevens4/wAnalysis/efficXsec/outEmb/gainUp2/"){ 
00016 
00017   char* sign;
00018   if(charge==0) sign="+";
00019   else if(charge==1) sign="-";
00020   
00021   char* core0; //these have no vertex reweighting
00022   if(charge == 1) core0="Wminus";
00023   else if(charge == 0) core0="Wplus";
00024 
00025   ofstream latex(Form("./%s.txt",core0));
00026 
00027   //gStyle->SetOptFit(1);
00028   //gStyle->SetOptStat(0);
00029   gStyle->SetPalette(1);
00030 
00031   //load file
00032   TString fullInpName=iPath;  fullInpName+=core0; 
00033   fullInpName+=".wana.hist.root";
00034   fd=new TFile(fullInpName);
00035   if(!fd->IsOpen()) {
00036     printf("ERROR: input histo file not found, quit\n",fullInpName.Data());
00037     return;
00038   } else {
00039     printf("Opened: %s\n",fullInpName.Data());
00040   }
00041   
00042   //use rebin 1 for effic # calc, but rebin N for plot
00043   int etrebin=4;
00044   int etaRebin=50;
00045   int zRebin=4;
00046   int phiModuleRebin=6;
00047   int phiRebin=1;
00048   int zdcRebin=5;
00049 
00050   //get total efficiency and error
00051   float etplot[5]={1.0,1.0,1.0,1.0,1.0};
00052   float etplotErr[5]={1.0,1.0,1.0,1.0,0.0};
00053   float zdcplot[5]={2.0,2.0,2.0,2.0,2.0};
00054   float zdcplotErr[5]={1.0,1.0,1.0,1.0,0.0};
00055 
00056   //trigger efficiency
00057   TH1F *hTrigET=doEfficiency(fd,"MCeleETall","MCeleETtrig",Form("W%s trigger efficiency",sign),etrebin,etplot[0],etplotErr[0]);
00058   TH1F *hTrigEta=doEfficiency(fd,"MCeleEtaAll","MCeleEtaTrig",Form("W%s trigger efficiency",sign),etaRebin,0);
00059   TH1F *hTrigDetEta=doEfficiency(fd,"MCeleDetEtaAll","MCeleDetEtaTrig",Form("W%s trigger efficiency",sign),etaRebin,0);
00060   TH1F *hTrigZvert=doEfficiency(fd,"MCeleZvertAll","MCeleZvertTrig",Form("W%s trigger efficiency",sign),zRebin,0);
00061   TH1F *hTrigPhi=doEfficiency(fd,"MCelePhiModulePairAll","MCelePhiModulePairTrig",Form("W%s trigger efficiency",sign),phiModuleRebin,0); 
00062   TH1F *hTrigZDC=doEfficiency(fd,"MCeleZdcAll","MCeleZdcTrig",Form("W%s trigger efficiency",sign),zdcRebin,zdcplot[0],zdcplotErr[0]);
00063   //vertex efficiency
00064   TH1F *hVertET=doEfficiency(fd,"MCeleETtrig","MCeleETvert",Form("W%s vertex efficiency",sign),etrebin,etplot[1],etplotErr[1]);
00065   TH1F *hVertEta=doEfficiency(fd,"MCeleEtaTrig","MCeleEtaVert",Form("W%s vertex efficiency",sign),etaRebin,0);
00066   TH1F *hVertZvert=doEfficiency(fd,"MCeleZvertTrig","MCeleZvertVert",Form("W%s vertex efficiency",sign),zRebin,0);
00067   TH1F *hVertPhi=doEfficiency(fd,"MCelePhiTrig","MCelePhiVert",Form("W%s vertex efficiency",sign),phiRebin,0);
00068   TH1F *hVertZDC=doEfficiency(fd,"MCeleZdcTrig","MCeleZdcVert",Form("W%s vertex efficiency",sign),zdcRebin,zdcplot[1],zdcplotErr[1]);
00069   //tracking efficiency
00070   TH1F *hTrackET=doEfficiency(fd,"MCeleETvert","MCeleETTrack",Form("W%s tracking efficiency",sign),etrebin,etplot[2],etplotErr[2]);
00071   TH1F *hTrackEta=doEfficiency(fd,"MCeleEtaVert","MCeleEtaTrack",Form("W%s tracking efficiency",sign),etaRebin,0);
00072   TH1F *hTrackZvert=doEfficiency(fd,"MCeleZvertVert","MCeleZvertTrack",Form("W%s tracking efficiency",sign),zRebin,0);
00073   TH1F *hTrackPhi=doEfficiency(fd,"MCelePhiVert","MCelePhiTrack",Form("W%s tracking efficiency",sign),phiRebin,0);
00074   TH1F *hTrackZDC=doEfficiency(fd,"MCeleZdcVert","MCeleZdcTrack",Form("W%s tracking efficiency",sign),zdcRebin,zdcplot[2],zdcplotErr[2]);
00075   //algo efficiency
00076   TH1F *hRecoET=doEfficiency(fd,"MCeleETTrack","MCeleETreco",Form("W%s algo efficiency",sign),etrebin,etplot[3],etplotErr[3]);
00077   TH1F *hRecoEta=doEfficiency(fd,"MCeleEtaTrack","MCeleEtaReco",Form("W%s algo efficiency",sign),etaRebin,0);
00078   TH1F *hRecoZvert=doEfficiency(fd,"MCeleZvertTrack","MCeleZvertReco",Form("W%s algo efficiency",sign),zRebin,0);
00079   TH1F *hRecoPhi=doEfficiency(fd,"MCelePhiTrack","MCelePhiReco",Form("W%s algo efficiency",sign),phiRebin,0);
00080   TH1F *hRecoZDC=doEfficiency(fd,"MCeleZdcTrack","MCeleZdcReco",Form("W%s algo efficiency",sign),zdcRebin,zdcplot[3],zdcplotErr[3]);
00081 
00082   //total efficiency for overall stat uncertainty
00083   TH1F *hTotET=doEfficiency(fd,"MCeleETall","MCeleETreco",Form("W%s total efficiency",sign),etrebin,etplot[4],etplotErr[4]);
00084   TH1F *hTotZDC=doEfficiency(fd,"MCeleZdcAll","MCeleZdcReco",Form("W%s total efficiency",sign),zdcRebin,zdcplot[4],zdcplotErr[4]);
00085   float totEffic=etplot[0]*etplot[1]*etplot[2]*etplot[3];
00086   //cout<<totEffic<<endl;
00087   cout<<"Total Efficiency = "<<etplot[4]<<" $\\pm$ "<<etplotErr[4]<<endl;  
00088   //spreadsheet<<"Track Effic = "<<etplot[2]<<" +/- "<<etplotErr[2]<<endl;
00089 
00090   spreadsheet<<"W"<<sign<<" Total Efficiency"<<endl;
00091   for(int j=0; j<6; j++){
00092     spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hTotET->GetBinContent(7+j))<<","<<Form("%.4f",hTotET->GetBinError(7+j))<<endl;
00093   }
00094 
00095   //set axis range for ET plot
00096   int etAxisMin=0; int etAxisMax=70;
00097 
00098   //***************  Algo Efficiency plots ***************//
00099   cA=new TCanvas(Form("W%s algo effic",sign),"algo effic",800,600);
00100   cA->Divide(2,2);
00101   cA->cd(1);
00102   hRecoET->Draw(); 
00103   hRecoET->SetAxisRange(etAxisMin,etAxisMax);
00104   cA->cd(2);
00105   hRecoEta->Draw();
00106   cA->cd(3);
00107   hRecoZvert->Draw();
00108   cA->cd(4);
00109   hRecoZDC->Draw();
00110   //hRecoPhi->Draw();
00111   cA->Print(Form("plots/W%salgoEffic.eps",sign));
00112   cA->Print(Form("plots/W%salgoEffic.png",sign));
00113 
00114   //algo efficiency numbers
00115   spreadsheet<<"W"<<sign<<" Algo Efficiency"<<endl;
00116   latex<<"W"<<sign<<" Algo Efficiency"<<endl;
00117   for(int j=0; j<6; j++){
00118     //semi-latex format
00119     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;
00120     
00121     //for spreadsheet
00122     spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hRecoET->GetBinContent(7+j))<<","<<Form("%.4f",hRecoET->GetBinError(7+j))<<endl;
00123   }
00124 
00125   //***************  Tracking Efficiency plots ***********//
00126   cT=new TCanvas(Form("W%s track effic",sign),"track effic",800,600);
00127   cT->Divide(2,2);
00128   cT->cd(1);
00129   hTrackET->Draw(); 
00130   hTrackET->SetAxisRange(etAxisMin,etAxisMax);
00131   cT->cd(2);
00132   hTrackEta->Draw();
00133   cT->cd(3);
00134   //hTrackZvert->Draw();
00135   hTrackZDC->SetMinimum(0.6); 
00136   hTrackZDC->Draw();
00137   cT->cd(4);
00138   hTrackPhi->Draw();
00139   cT->Print(Form("plots/W%strackEffic.eps",sign));
00140   cT->Print(Form("plots/W%strackEffic.png",sign));
00141 
00142   //***************  Vertex Efficiency plots *************//
00143   //Draw trigger efficiency 
00144   cV=new TCanvas(Form("W%s vertex effic",sign),"vertex effic",800,600);
00145   cV->Divide(2,2);
00146   cV->cd(1);
00147   hVertET->Draw(); 
00148   hVertET->SetAxisRange(etAxisMin,etAxisMax);
00149   cV->cd(2);
00150   hVertEta->Draw();
00151   cV->cd(3);
00152   hVertZvert->Draw();
00153   cV->cd(4);
00154   hVertZDC->SetMinimum(0.6);
00155   hVertZDC->Draw();
00156   cV->Print(Form("plots/W%svertEffic.eps",sign));
00157   cV->Print(Form("plots/W%svertEffic.png",sign));
00158 
00159   //***************  Trigger Efficiency plots ************//
00160   //draw trigger efficiency to show why W- isn't const w/ ET
00161   c2D=new TCanvas(Form("W%s trigger effic ET bins",sign),"trigger effic ET bins",700,500);
00162   j1=(TH2F*)fd->Get("MCeleEta_ptPreTrig"); //j1->Rebin2D(2,5);
00163   TH1D* h25[10];
00164   c2D->Divide(2,2);
00165   spreadsheet<<"W"<<sign<<" Trigger Efficiency"<<endl;
00166   latex<<"W"<<sign<<" Trigger Efficiency"<<endl;
00167   for(int j=0; j<6; j++){
00168     h25[j] = j1->ProjectionX(Form("pt%d_%d",25+j*4,29+j*4),26+j*4,29+j*4);
00169     h25[j]->SetTitle(Form("Lepton detector #eta (from Geant): PT=[%d,%d]",25+j*4,29+j*4));
00170     h25[j]->Rebin(2); 
00171     
00172     //count frac. of events with |eta|<1 and effic for each bin
00173     float accepted = h25[j]->Integral(26,125);
00174     float total = h25[j]->Integral();
00175     //cout<<"ET=["<<25+j*4<<","<<29+j*4<<"] efficiency = "<<hTrigET->GetBinContent(7+j)<<" ; fraction with |eta| < 1 : "<<accepted<<" "<<total<<endl;
00176 
00177     //semi-latex format
00178     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;
00179     spreadsheet<<25+j*4<<","<<29+j*4<<","<<Form("%.4f",hTrigET->GetBinContent(7+j))<<","<<Form("%.4f",hTrigET->GetBinError(7+j))<<endl;
00180 
00181     //add lines for detector cutoff at eta=1
00182     h25[j]->Rebin(5);
00183     Lx=h25[j]->GetListOfFunctions();
00184     int max=h25[j]->GetMaximum();
00185     ln1=new TLine(-1,0,-1,max);
00186     ln2=new TLine(1,0,1,max);
00187     ln1->SetLineWidth(2); ln2->SetLineWidth(2);
00188     ln1->SetLineColor(kRed); ln1->SetLineStyle(2);
00189     ln2->SetLineColor(kRed); ln2->SetLineStyle(2);
00190     Lx->Add(ln1); Lx->Add(ln2);
00191     if(j<4){ //only plot first four ET bins
00192       c2D->cd(j+1);
00193       h25[j]->Draw("h");
00194     }
00195   }
00196   c2D->Print(Form("plots/W%strigEfficNonConst.eps",sign));
00197   c2D->Print(Form("plots/W%strigEfficNonConst.png",sign));
00198   
00199   //Draw trigger efficiency 
00200   c=new TCanvas(Form("W%s trigger effic",sign),"trigger effic",800,600);
00201   c->Divide(2,2);
00202   c->cd(1);
00203   hTrigET->Draw(); 
00204   hTrigET->SetAxisRange(etAxisMin,etAxisMax);
00205   c->cd(2);
00206   hTrigEta->Draw();
00207   c->cd(3);
00208   hTrigDetEta->Draw();
00209   c->cd(4);
00210   hTrigPhi->Draw();
00211   c->Print(Form("plots/W%strigEffic.eps",sign));
00212   c->Print(Form("plots/W%strigEffic.png",sign));
00213   
00214 
00215   //other efficiency related histos
00216   TH2F * g0=(TH2F * )fd->Get("muBdist1");  assert(g0);
00217   TH2F * g1=(TH2F * )fd->Get("muBclET24R_ET");  assert(g1);
00218   TH2F * g2=(TH2F * )fd->Get("muBclEjetE2D_ET");  assert(g2);
00219   TH2F * g3=(TH2F * )fd->Get("musPtBalance_clust");  assert(g3);
00220   cAlgo2=new TCanvas(Form("W%s algo effic 2 ",sign),"algo effic 2",800,600);
00221   cAlgo2->Divide(2,2);
00222   cAlgo2->cd(1);
00223   g0->Draw("colz");
00224   g0->GetXaxis()->SetRange(0,70);
00225   cAlgo2->cd(2);
00226   g1->Draw("colz");
00227   g1->GetXaxis()->SetRange(0,70);
00228   cAlgo2->cd(3);
00229   g2->Draw("colz");
00230   g2->GetXaxis()->SetRange(0,70);
00231   cAlgo2->cd(4);
00232   g3->Draw("colz");
00233   g3->GetXaxis()->SetRange(0,70);
00234   cAlgo2->Print(Form("plots/W%salgoEffic2.eps",sign));
00235   cAlgo2->Print(Form("plots/W%salgoEffic2.png",sign));
00236 
00237   //cout<<g2->GetYaxis()->FindBin(0.88)<<endl;
00238   //cout<<g2->Integral(26,30,0,1000)<<" "<<g2->Integral(26,30,1,g2->GetYaxis()->FindBin(0.88))<<endl;
00239   //cout<<g2->Integral(26,1000,1,1000)<<" "<<g2->Integral(26,1000,1,g2->GetYaxis()->FindBin(0.88))<<endl;
00240 
00241   return;
00242 }
00243 
00244 //calculate efficiencies
00245 TH1F* doEfficiency(TFile *fd,char* name0, char* name1, char *tit, int reb=1, float &etplot, float &etplotErr=0){
00246   
00247   TH1F * h0=(TH1F * )fd->Get(name0);  assert(h0);
00248   TH1F * h1=(TH1F * )fd->Get(name1);  assert(h1);
00249 
00250   ha=(TH1F*) h0->Clone(); assert(ha);
00251   hb=(TH1F*) h1->Clone(); assert(hb);
00252 
00253   ha->Rebin(reb);
00254   hb->Rebin(reb);
00255 
00256   hc=(TH1F*) hb->Clone(); assert(hc);
00257   hc->SetTitle(tit);
00258   hc->Reset();
00259 
00260   float num=0; float den=0;
00261   float matchErr2=0; float thrownErr2=0;
00262   int nb=hb->GetNbinsX();
00263   int startBin,endBin; 
00264   if(etplot==1.0) {
00265     startBin=hb->GetXaxis()->FindBin(25.);
00266     //endBin=hb->GetXaxis()->FindBin(48.0);
00267     //startBin=hb->GetXaxis()->FindBin(50.);
00268     endBin=hb->GetXaxis()->FindBin(100.0);
00269   }
00270   if(etplot==2.0) {
00271     endBin=hb->GetXaxis()->FindBin(200000.);
00272     startBin=hb->GetXaxis()->FindBin(0.);
00273   }
00274     
00275   //cout<<startBin<<endl;
00276   //cout<<hb->GetNbinsX()<<" "<<ha->GetNbinsX()<<endl;
00277   assert(nb==ha->GetNbinsX());
00278   for(int i=1; i<=nb; i++) {
00279     float n0=ha->GetBinContent(i);//n0: # of thrown tracks
00280     float n1=hb->GetBinContent(i);//n1: # of matched tracks
00281     float e0=ha->GetBinError(i);//e0: error on # of thrown tracks
00282     float e1=hb->GetBinError(i);//e1: error on # of matched tracks
00283     if(n0==0) continue; // no thown tracks found
00284     float eff=(float) n1/n0;
00285     
00286     //old method with un-weighted histos
00287     //float x=n1*(n0-n1);
00288     //if(n1==0 || n1==n0) x=n0;
00289     //float errEff=sqrt(x/n0/n0/n0); 
00290     
00291     float errEff=0;
00292     //method for using weighted histos (must use Sumw2)
00293     if(n1==n0) {//problem in effic=1 limit, statErr=0
00294       //cout<<name1<<" bin="<<i<<" problem with error"<<endl;
00295       errEff=1./n0; //just something big to see its low stats
00296     } 
00297     else errEff=sqrt(e1*e1*(n0-n1)*(n0-n1)+(e0-e1)*(e0-e1)*n1*n1)/n0/n0;
00298        
00299     //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);
00300     hc->SetBinContent(i,eff);
00301     hc->SetBinError(i,errEff);
00302     if(i>=startBin && i<=endBin){//define ET range for efficTot
00303       num+=n1;
00304       den+=n0;
00305       matchErr2+=e1*e1;
00306       thrownErr2+=e0*e0;
00307     }
00308   }
00309   
00310   if(etplot==1.0){
00311     //old method with un-weighted histos
00312     //cout<<name1<<"="<<(float) num/den<<"+/-"<<sqrt((float)num*(den-num)/den/den/den)<<endl;
00313     
00314     //method for using weighted histos (must use Sumw2)
00315     if(etplotErr>0.0) cout<<name1<<"="<< num/den <<" $\\pm$ "<< sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den <<endl;
00316     etplot=num/den;
00317     etplotErr=sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den;
00318   }
00319   if(etplot==2.0){
00320     //method for using weighted histos (must use Sumw2)
00321     //cout<<name1<<"="<< num/den <<" $\\pm$ "<< sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den <<endl;
00322 
00323     etplot=num/den;
00324     etplotErr=sqrt(matchErr2*(den-num)*(den-num) + (thrownErr2-matchErr2)*num*num)/den/den;
00325   }
00326   hc->SetMinimum(0.0);
00327   hc->SetMaximum(1.1);
00328 
00329   return hc;
00330 }