]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
major updates to plotting macros
authorjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Oct 2009 07:01:08 +0000 (07:01 +0000)
committerjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Oct 2009 07:01:08 +0000 (07:01 +0000)
PWG4/macros/electrons/makeCombinedData.C [new file with mode: 0644]
PWG4/macros/electrons/makeSysErr.C [new file with mode: 0644]
PWG4/macros/electrons/plotMCRates.C
PWG4/macros/electrons/plotNPERates.C

diff --git a/PWG4/macros/electrons/makeCombinedData.C b/PWG4/macros/electrons/makeCombinedData.C
new file mode 100644 (file)
index 0000000..4ece6c3
--- /dev/null
@@ -0,0 +1,367 @@
+/////////////////////////////////////////
+// A set of plotting utilities for the
+// electron chapter of the PPR
+//
+// J.L.Klay (Cal Poly)
+// 28-Oct-2009
+////////////////////////////////////////
+
+TH1F *alltte,    *alltrk,    *allemc;
+TH1F *sumtte,    *sumtrk,    *sumemc;
+TH1F *sumHFemc,  *sumHFmc;
+TH1F *btte,      *btrk,      *bptemc;
+TH1F *ctte,      *ctrk,      *cptemc;
+TH1F *cbtte,     *cbtrk,     *cbemc;
+TH1F *convtte,   *convtrk,   *convemc;
+TH1F *daltte,    *daltrk,    *dalemc;
+TH1F *wztte,     *wzpttrk,   *wzemc;
+TH1F *othtte,    *othtrk,    *othemc;
+TH1F *htte,      *htrk,      *hemc;
+
+TH1F *allmc;
+TH1F *sigemc, *bkgemc, *wallemc, *hijemc;
+TH1F* belemc, *celemc, *candbmc;
+TH1F *convmc, *dalmc, *wzmc, *othermc;
+TH1F* mchad;
+
+void makeData(char* hijfname = "data/scaled25Oct09/histosLHC08d6.root",
+              char* jjfname = "data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root",
+              char* bfname = "data/scaled25Oct09/histosscaledLHC09b4AODc.root",
+              char* wfname = "data/scaled25Oct09/histosWboson.root") {
+
+
+  //For HIJING need to divide by the number of events, which we
+  //can get from the file and do when we perform scaling
+  double hijscale = 0.05*(1.E6)*0.5*7700; //0-5% * seconds*lumi*PbPb
+                                         //x-section                      
+  //For bjet and jet-jet events
+  double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance
+  double bscale = pyscale; //Do we need to scale by Branching ratio
+                          //for forced                      
+  //semi-leptonic decays?                                             
+  double wscale = pyscale;
+
+  TFile* hijfile = new TFile(hijfname);
+  if(!hijfile) { printf("NO HIJING FILE\n"); return; }
+  TList* hijlist = (TList*)hijfile->Get("histos");
+  TH2F* hijtte = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleTTE");
+  TH2F* hijemc2d = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleEMCAL");
+  TH2F* hijtrk = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleTPCTRD");
+  TH1F* hijmult = (TH1F*)histos->FindObject("AnaElectron_hRefMult");
+  Int_t nEvt = hijmult->GetEntries();
+  if(nEvt == 0) { printf("NO HIJING EVENTS\n"); return; }
+  hijtte->Scale(hijscale/nEvt);
+  hijemc2d->Scale(hijscale/nEvt);
+  hijtrk->Scale(hijscale/nEvt);
+  TH2F* hijmcele = (TH2F*)histos->FindObject("AnaElectron_hPtMCElectron");
+  TH1F* hijmchad = (TH1F*)histos->FindObject("AnaElectron_hPtMCHadron");
+  hijmcele->Scale(hijscale/nEvt);
+  hijmchad->Scale(hijscale/nEvt);
+
+  TFile* jjfile = new TFile(jjfname);
+  if(!jjfile) { printf("NO JET-JET FILE\n"); return; }
+  TH2F* jjtte = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* jjemc = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleEMCALScaled");
+  TH2F* jjtrk = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleTPCTRDScaled");
+  TH1F* jjmult = (TH1F*)jjfile->Get("AnaElectron_hRefMultScaled");
+  Int_t nEvtJJ = jjmult->GetEntries();
+  jjtte->Scale(pyscale);
+  jjemc->Scale(pyscale);
+  jjtrk->Scale(pyscale);
+  TH2F* jjmcele = (TH2F*)jjfile->Get("AnaElectron_hPtMCElectronScaled");
+  TH1F* jjmchad = (TH1F*)jjfile->Get("AnaElectron_hPtMCHadronScaled");
+  jjmcele->Scale(pyscale);
+  jjmchad->Scale(pyscale);
+
+  TFile* bfile = new TFile(bfname);
+  if(!bfile) { printf("NO B-JET FILE\n"); return; }
+  TH2F* bjtte = (TH2F*)bfile->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* bjemc = (TH2F*)bfile->Get("AnaElectron_hPtNPEleEMCALScaled");
+  TH2F* bjtrk = (TH2F*)bfile->Get("AnaElectron_hPtNPEleTPCTRDScaled");
+  TH1F* bjmult = (TH1F*)bfile->Get("AnaElectron_hRefMultScaled");
+  Int_t nEvtB = bjmult->GetEntries();
+  bjtte->Scale(bscale);
+  bjemc->Scale(bscale);
+  bjtrk->Scale(bscale);
+  TH2F* bjmcele = (TH2F*)bfile->Get("AnaElectron_hPtMCElectronScaled");
+  TH1F* bjmchad = (TH1F*)bfile->Get("AnaElectron_hPtMCHadronScaled");
+  bjmcele->Scale(bscale);
+  bjmchad->Scale(bscale);
+
+  TFile* wfile = new TFile(wfname);
+  if(!wfile) { printf("NO W-BOSON FILE\n"); return; }
+  TH2F* wjtte = (TH2F*)wfile->Get("AnaElectron_hPtNPEleTTE");
+  TH2F* wjemc = (TH2F*)wfile->Get("AnaElectron_hPtNPEleEMCAL");
+  TH2F* wjtrk = (TH2F*)wfile->Get("AnaElectron_hPtNPEleTPCTRD");
+  TH1F* wjmult = (TH1F*)wfile->Get("AnaElectron_hRefMult");
+  Int_t nEvtW = wjmult->GetEntries();
+  wjtte->Scale(wscale); //already scaled by evts
+  wjemc->Scale(wscale); //already scaled by evts
+  wjtrk->Scale(wscale); //already scaled by evts
+  TH2F* wjmcele = (TH2F*)wfile->Get("AnaElectron_hPtMCElectron");
+  TH1F* wjmchad = (TH1F*)wfile->Get("AnaElectron_hPtMCHadron");
+  wjmcele->Scale(wscale);
+  wjmchad->Scale(wscale);
+
+  printf("Event statistics: %d (HIJING)  %d (JET-JET)  %d (B-JET)  %d (W-Boson)\n",nEvt,nEvtJJ,nEvtB,nEvtW);
+
+  TH2F* combined = (TH2F*)hijmcele->Clone();
+  combined->Add(jjmcele);
+  combined->Add(bjmcele);  
+  combined->Add(wjmcele);
+  combined->SetTitle("MC electrons in Pb+Pb in EMCAL acceptance");
+  combined->SetName("CombinedMCEle");
+  combined->SetXTitle("p_T (GeV/c)");
+
+  mchad = (TH1F*)hijmchad->Clone();
+  mchad->Add(jjmchad);
+  mchad->Add(bjmchad);
+  mchad->Add(wjmchad);
+  mchad->SetTitle("MC hadrons in Pb+Pb in EMCAL acceptance");
+  mchad->SetName("CombinedMCHad");
+  mchad->SetXTitle("p_T (GeV/c)");
+
+  allmc = (TH1F*)combined->ProjectionX("allmc",1,1);
+  belemc = (TH1F*)combined->ProjectionX("bmc",2,2);
+  celemc = (TH1F*)combined->ProjectionX("cmc",3,3);
+  candbmc = (TH1F*)combined->ProjectionX("candbmc",4,4);
+  convmc = (TH1F*)combined->ProjectionX("convmc",5,5);
+  dalmc = (TH1F*)combined->ProjectionX("dalmc",6,6);
+  wzmc = (TH1F*)combined->ProjectionX("wzmc",7,7);
+  othermc = (TH1F*)combined->ProjectionX("othermc",8,8);
+
+  //For comparing contributions                        
+  wallemc = (TH1F*)wjmcele->ProjectionX("wallemc",7,7);
+  sigemc = (TH1F*)bjmcele->ProjectionX("sigemc",1,1);
+  bkgemc = (TH1F*)jjmcele->ProjectionX("bkgemc",1,1);
+  hijemc = (TH1F*)hijmcele->ProjectionX("hijemc",1,1);
+
+  double myscale = 1.; //we already scaled them       
+  ScaleAndConfigure(allmc,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(belemc,myscale,kRed,kFALSE);
+  ScaleAndConfigure(sigemc,myscale,kRed,kFALSE);
+  ScaleAndConfigure(celemc,myscale,kBlue,kFALSE);
+  ScaleAndConfigure(candbmc,myscale,kViolet,kFALSE);
+  ScaleAndConfigure(convmc,myscale,kOrange-3,kFALSE);
+  ScaleAndConfigure(bkgemc,myscale,kOrange-3,kFALSE);
+  ScaleAndConfigure(dalmc,myscale,kGreen-3,kFALSE);
+  ScaleAndConfigure(wzmc,myscale,kOrange-7,kFALSE);
+  ScaleAndConfigure(wallemc,myscale,kOrange-7,kFALSE);
+  ScaleAndConfigure(mchad,myscale,kGreen+2,kFALSE);
+  ScaleAndConfigure(hijemc,myscale,kGreen+2,kFALSE);
+
+  TH2F* combTTE = (TH2F*)hijtte->Clone();
+  combTTE->Add(jjtte);
+  combTTE->Add(bjtte);
+  combTTE->Add(wjtte);
+  combTTE->SetTitle("Identified non-phot. electrons (TPC+TRD+EMCAL)");
+  combTTE->SetName("CombinedEleTTE");
+  combTTE->SetXTitle("p_T (GeV/c)");
+
+  alltte = (TH1F*)combTTE->ProjectionX("alltte",1,1);
+  btte = (TH1F*)combTTE->ProjectionX("btte",2,2);
+  ctte = (TH1F*)combTTE->ProjectionX("ctte",3,3);
+  cbtte = (TH1F*)combTTE->ProjectionX("cbtte",4,4);
+  convtte = (TH1F*)combTTE->ProjectionX("convtte",5,5);
+  daltte = (TH1F*)combTTE->ProjectionX("daltte",6,6);
+  wztte = (TH1F*)combTTE->ProjectionX("wztte",7,7);
+  othtte = (TH1F*)combTTE->ProjectionX("othtte",8,8);
+  sumtte = (TH1F*)btte->Clone(); sumtte->SetName("sumtte");
+  sumtte->Add(ctte); sumtte->Add(cbtte); sumtte->Add(convtte);
+  sumtte->Add(daltte); sumtte->Add(wztte); //sumtte->Add(othtte);  
+  htte = (TH1F*)combTTE->ProjectionX("htte",9,9);
+
+  ScaleAndConfigure(alltte,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(sumtte,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(btte,myscale,kRed,kFALSE);
+  ScaleAndConfigure(ctte,myscale,kBlue,kFALSE);
+  ScaleAndConfigure(cbtte,myscale,kViolet,kFALSE);
+  ScaleAndConfigure(convtte,myscale,kOrange-3,kFALSE);
+  ScaleAndConfigure(daltte,myscale,kGreen-3,kFALSE);
+  ScaleAndConfigure(htte,myscale,kGreen+2,kFALSE);
+  ScaleAndConfigure(wztte,myscale,kOrange-7,kFALSE);
+
+  TH2F* combEMC = (TH2F*)hijemc2d->Clone();
+  combEMC->Add(jjemc);
+  combEMC->Add(bjemc);
+  combEMC->Add(wjemc);
+  combEMC->SetTitle("Identified non-phot. electrons (EMCAL)");
+  combEMC->SetName("CombinedEleEMC");
+  combEMC->SetXTitle("p_T (GeV/c)");
+
+  allemc = (TH1F*)combEMC->ProjectionX("allemc",1,1);
+  bemc = (TH1F*)combEMC->ProjectionX("bemc",2,2);
+  cemc = (TH1F*)combEMC->ProjectionX("cemc",3,3);
+  cbemc = (TH1F*)combEMC->ProjectionX("cbemc",4,4);
+  convemc = (TH1F*)combEMC->ProjectionX("convemc",5,5);
+  dalemc = (TH1F*)combEMC->ProjectionX("dalemc",6,6);
+  wzemc = (TH1F*)combEMC->ProjectionX("wzemc",7,7);
+  othemc = (TH1F*)combEMC->ProjectionX("othemc",8,8);
+  hemc = (TH1F*)combEMC->ProjectionX("hemc",9,9);
+
+  sumemc = (TH1F*)bemc->Clone(); sumemc->SetName("sumemc");
+  sumemc->Add(cemc); sumemc->Add(cbemc); sumemc->Add(convemc);
+  sumemc->Add(dalemc); //sumemc->Add(othemc);
+  sumemc->Add(wzemc);
+  sumHFemc = (TH1F*)bemc->Clone(); sumHFemc->SetName("sumHFemc");
+  sumHFemc->Add(cemc); sumHFemc->Add(cbemc); sumHFemc->Add(wzemc);
+
+  ScaleAndConfigure(allemc,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(sumemc,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(sumHFemc,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(bemc,myscale,kRed,kFALSE);
+  ScaleAndConfigure(cemc,myscale,kBlue,kFALSE);
+  ScaleAndConfigure(cbemc,myscale,kViolet,kFALSE);
+  ScaleAndConfigure(convemc,myscale,kOrange-3,kFALSE);
+  ScaleAndConfigure(dalemc,myscale,kGreen-3,kFALSE);
+  ScaleAndConfigure(hemc,myscale,kGreen+2,kFALSE);
+  ScaleAndConfigure(wzemc,myscale,kOrange-7,kFALSE);
+
+  TH2F* combTRK = (TH2F*)hijtrk->Clone();
+  combTRK->Add(jjtrk);
+  combTRK->Add(bjtrk);
+  combTRK->Add(wjtrk);
+  combTRK->SetTitle("Identified non-phot. electrons (TPC+TRD)");
+  combTRK->SetName("CombinedEleTRK");
+  combTRK->SetXTitle("p_T (GeV/c)");
+
+  alltrk = (TH1F*)combTRK->ProjectionX("alltrk",1,1);
+  btrk = (TH1F*)combTRK->ProjectionX("btrk",2,2);
+  ctrk = (TH1F*)combTRK->ProjectionX("ctrk",3,3);
+  cbtrk = (TH1F*)combTRK->ProjectionX("cbtrk",4,4);
+  convtrk = (TH1F*)combTRK->ProjectionX("convtrk",5,5);
+  daltrk = (TH1F*)combTRK->ProjectionX("daltrk",6,6);
+  wztrk = (TH1F*)combTRK->ProjectionX("wztrk",7,7);
+  othtrk = (TH1F*)combTRK->ProjectionX("othtrk",8,8);
+  sumtrk = (TH1F*)btrk->Clone(); sumtrk->SetName("sumtrk");
+  sumtrk->Add(ctrk); sumtrk->Add(cbtrk); sumtrk->Add(convtrk);
+  sumtrk->Add(daltrk); //sumtrk->Add(othtrk);  sumtrk->Add(wztrk);
+  htrk = (TH1F*)combTRK->ProjectionX("htrk",9,9);
+
+  for(Int_t i = 1; i <= alltrk->GetNbinsX(); i++) {
+    Double_t myall = alltrk->GetBinContent(i);
+    Double_t partsum = btrk->GetBinContent(i) +
+      ctrk->GetBinContent(i) + cbtrk->GetBinContent(i) +
+      convtrk->GetBinContent(i) + daltrk->GetBinContent(i) +
+      wztrk->GetBinContent(i) + othtrk->GetBinContent(i);
+    Double_t mysum = partsum + htrk->GetBinContent(i);
+    printf("<%d> Compare bins All: %d, Sum: %d   Had: %d partSum: %d\n",i,myall,mysum,htrk->GetBinContent(i),partsum);
+  }
+
+  ScaleAndConfigure(alltrk,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(sumtrk,myscale,kBlack,kFALSE);
+  ScaleAndConfigure(btrk,myscale,kRed,kFALSE);
+  ScaleAndConfigure(ctrk,myscale,kBlue,kFALSE);
+  ScaleAndConfigure(cbtrk,myscale,kViolet,kFALSE);
+  ScaleAndConfigure(convtrk,myscale,kOrange-3,kFALSE);
+  ScaleAndConfigure(daltrk,myscale,kGreen-3,kFALSE);
+  ScaleAndConfigure(htrk,myscale,kGreen+2,kFALSE);
+  ScaleAndConfigure(wztrk,myscale,kOrange-7,kFALSE);
+
+  TFile* fout = new TFile("CombinedCocktailHistograms.root","RECREATE");
+  fout->cd();
+  combined->Write();
+  combTTE->Write();
+  combTRK->Write();
+  combEMC->Write();
+  alltte->Write();
+  alltrk->Write();
+  allemc->Write();
+  sumtte->Write();
+  sumtrk->Write();
+  sumemc->Write();
+  btte->Write();
+  btrk->Write();
+  bemc->Write();
+  ctte->Write();
+  ctrk->Write();
+  cemc->Write();
+  cbtte->Write();
+  cbtrk->Write();
+  cbemc->Write();
+  convtte->Write();
+  convtrk->Write();
+  convemc->Write();
+  daltte->Write();
+  daltrk->Write();
+  dalemc->Write();
+  othtte->Write();
+  othtrk->Write();
+  othemc->Write();
+  wztte->Write();
+  wztrk->Write();
+  wzemc->Write();
+  htte->Write();
+  htrk->Write();
+  hemc->Write();
+  allmc->Write();
+  belemc->Write();
+  celemc->Write();
+  candbmc->Write();
+  convmc->Write();
+  dalmc->Write();
+  wzmc->Write();
+  othermc->Write();
+  mchad->Write();
+  sigemc->Write();
+  bkgemc->Write();
+  wallemc->Write();
+  hijemc->Write();
+  fout->Close();
+
+}
+
+
+void ScaleAndConfigure(TH1F* hist,Double_t scale, Int_t color,Bool_t keepErr)
+{
+  hist->Scale(scale);
+  hist->SetLineColor(color);
+  hist->SetMarkerColor(color);
+  hist->SetLineWidth(2);
+  if(keepErr == kFALSE) {
+    //remove the error bars - useful for MC rates
+    for(Int_t i = 1; i <= hist->GetNbinsX(); i++) {
+      if(hist->GetBinContent(i) > 0.) {
+        if(hist->GetBinError(i)/hist->GetBinContent(i) > 0.5) {
+          Double_t avg = 0.;
+          if(i > 1 && i < hist->GetNbinsX())
+            avg = (hist->GetBinContent(i-1) + hist->GetBinContent(i+1))/2.;
+          hist->SetBinContent(i,avg);
+        }
+      }
+      hist->SetBinError(i,0.);
+    }
+  } else {
+    //Set the error bars to statistics of the bin content
+    for(Int_t i = 1; i <= hist->GetNbinsX(); i++) {
+      if(hist->GetBinContent(i) > 0.) {
+        if(hist->GetBinError(i)/hist->GetBinContent(i) > 0.5) {
+          Double_t avg = 0;
+          if(i > 1 && i < hist->GetNbinsX())
+            avg = (hist->GetBinContent(i-1) + hist->GetBinContent(i+1))/2.;
+          hist->SetBinContent(i,avg);
+        }
+      }
+      hist->SetBinError(i,TMath::Sqrt(hist->GetBinContent(i)));
+    }
+  }
+}
+
+TH1F* GetPtCutHisto(TH1F* input)
+{
+  //Given a rate histogram vs pt, return the histogram with yield
+  //above a given pTcut                                          
+
+  TH1F* result = (TH1F*)input->Clone();
+  char name[100];
+  sprintf(name,"%s_ptCut",result->GetName());
+  result->SetNameTitle(name,name);
+  for(Int_t i = 1; i <= result->GetNbinsX(); i++) {
+    Double_t val = input->Integral(i,result->GetNbinsX());
+    result->SetBinContent(i,val);
+    result->SetBinError(i,0.);
+  }
+
+  return result;
+
+}
diff --git a/PWG4/macros/electrons/makeSysErr.C b/PWG4/macros/electrons/makeSysErr.C
new file mode 100644 (file)
index 0000000..c6bc26d
--- /dev/null
@@ -0,0 +1,114 @@
+
+
+
+void makeSysErr() {
+
+  TFile* f0 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root");
+  TFile* f1 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-1.root");
+  TFile* f2 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-2.root");
+  TFile* f3 = new TFile("data/scaled25Oct09/TOTALhistosscaled-LHC09b2-3.root");
+
+  double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance
+
+  TH2F* tte0 = (TH2F*)f0->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* emc0 = (TH2F*)f0->Get("AnaElectron_hPtNPEleEMCALScaled");
+
+  TH2F* tte1 = (TH2F*)f1->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* emc1 = (TH2F*)f1->Get("AnaElectron_hPtNPEleEMCALScaled");
+
+  TH2F* tte2 = (TH2F*)f2->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* emc2 = (TH2F*)f2->Get("AnaElectron_hPtNPEleEMCALScaled");
+
+  TH2F* tte3 = (TH2F*)f3->Get("AnaElectron_hPtNPEleTTEScaled");
+  TH2F* emc3 = (TH2F*)f3->Get("AnaElectron_hPtNPEleEMCALScaled");
+
+  tte0->Scale(pyscale);
+  tte1->Scale(pyscale);
+  tte2->Scale(pyscale);
+  tte3->Scale(pyscale);
+  emc0->Scale(pyscale);
+  emc1->Scale(pyscale);
+  emc2->Scale(pyscale);
+  emc3->Scale(pyscale);
+
+  alltte0 = (TH1F*)tte0->ProjectionX("alltte0",1,1);
+  allemc0 = (TH1F*)emc0->ProjectionX("allemc0",1,1);
+  alltte0->Rebin(5);
+  allemc0->Rebin(5);
+
+  alltte1 = (TH1F*)tte1->ProjectionX("alltte1",1,1);
+  allemc1 = (TH1F*)emc1->ProjectionX("allemc1",1,1);
+  alltte1->Rebin(5);
+  allemc1->Rebin(5);
+
+  alltte2 = (TH1F*)tte2->ProjectionX("alltte2",1,1);
+  allemc2 = (TH1F*)emc2->ProjectionX("allemc2",1,1);
+  alltte2->Rebin(5);
+  allemc2->Rebin(5);
+
+  alltte3 = (TH1F*)tte3->ProjectionX("alltte3",1,1);
+  allemc3 = (TH1F*)emc3->ProjectionX("allemc3",1,1);
+  alltte3->Rebin(5);
+  allemc3->Rebin(5);
+  
+  Double_t tsum = 0.;
+  Double_t esum = 0.;
+  for(Int_t i = 1; i <= alltte0->GetNbinsX(); i++) {
+    Double_t t0 = alltte0->GetBinContent(i);
+    Double_t t1 = alltte1->GetBinContent(i);
+    Double_t t2 = alltte2->GetBinContent(i);
+    Double_t t3 = alltte3->GetBinContent(i);
+    Double_t e0 = allemc0->GetBinContent(i);
+    Double_t e1 = allemc1->GetBinContent(i);
+    Double_t e2 = allemc2->GetBinContent(i);
+    Double_t e3 = allemc3->GetBinContent(i);
+
+    Double_t td01 = (t0 - t1)/(t0+t1+t2+t3);
+    Double_t td02 = (t0 - t2)/(t0+t1+t2+t3);
+    Double_t td03 = (t0 - t3)/(t0+t1+t2+t3);
+
+    Double_t tave = (TMath::Abs(td01) + TMath::Abs(td02) + TMath::Abs(td03))/2.;
+    tsum += tave;
+
+    Double_t ed01 = (e0 - e1)/(e0+e1+e2+e3);
+    Double_t ed02 = (e0 - e2)/(e0+e1+e2+e3);
+    Double_t ed03 = (e0 - e3)/(e0+e1+e2+e3);
+
+    Double_t eave = (TMath::Abs(ed01) + TMath::Abs(ed02) + TMath::Abs(ed03))/2.;
+    esum += eave;
+  }
+  Double_t tfinal = tsum/alltte0->GetNbinsX();
+  Double_t efinal = esum/allemc0->GetNbinsX();
+  printf("tfinal %f, efinal %f\n",tfinal,efinal);
+  
+  TCanvas *c1 = new TCanvas();
+  gPad->SetLogy();
+  alltte0->SetLineWidth(2);
+  alltte0->Draw();
+
+  TCanvas *c2 = new TCanvas();
+  gPad->SetLogy();
+  allemc0->SetLineWidth(2);
+  allemc0->Draw();
+
+  TGraphErrors* terr = new TGraphErrors();
+  terr->SetName("tteErr");
+  TGraphErrors* eerr = new TGraphErrors();
+  eerr->SetName("emcErr");
+  for(Int_t i = 1; i <= alltte0->GetNbinsX(); i++) {
+    terr->SetPoint(i-1,alltte0->GetBinCenter(i),alltte0->GetBinContent(i));
+    terr->SetPointError(i-1,0.,(tfinal)*alltte0->GetBinContent(i));
+
+    eerr->SetPoint(i-1,allemc0->GetBinCenter(i),allemc0->GetBinContent(i));
+    eerr->SetPointError(i-1,0.,(efinal)*allemc0->GetBinContent(i));
+  }
+  c1->cd();
+  terr->SetFillColor(kRed-8);
+  terr->Draw("3same");
+  alltte0->Draw("same");
+
+  c2->cd();
+  eerr->SetFillColor(kRed-8);
+  eerr->Draw("3same");
+  allemc0->Draw("same");
+}
index f236a985e6611080920a1a895ba99ecb4be100a5..0c3ee751d7ef4ad887841e05cea1a6c5818ae819 100644 (file)
 //\r
 /////////////////////////////////////////////////\r
 \r
-TH1F* all;\r
-TH1F* bele;\r
-TH1F* cele;\r
-TH1F* candb;\r
-TH1F* conv;\r
-TH1F* dalitz;\r
-TH1F* wz;\r
-TH1F* other;\r
-TH1F* mchad;\r
-TH1F* sige;\r
-TH1F* bkge;\r
-TH1F* walle;\r
-TH1F* hije;\r
-\r
 TLegend* leg;\r
 \r
-void plotMCRates(char* hijfname = "data/scaled25Oct09/histosLHC08d6.root",\r
-                char* jjfname = "data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root",\r
-                char* bfname = "data/scaled25Oct09/histosscaledLHC09b4AODc.root",\r
-                char* wfname = "data/scaled25Oct09/histosWBoson.root") {\r
-\r
-  //For HIJING need to divide by the number of events, which we\r
-  //can get from the file and do when we perform scaling\r
-  double hijscale = 0.05*(1.E6)*0.5*7700; //0-5% * seconds*lumi*PbPb x-section\r
-  //For bjet and jet-jet events\r
-  double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance\r
-  double bscale = pyscale; //Do we need the Branching ratio for forced\r
-                               //semi-leptonic decays?\r
-  double wscale = pyscale; \r
-  \r
-  TFile* hijfile = new TFile(hijfname);\r
-  if(!hijfile) { printf("NO HIJING FILE\n"); return; }\r
-  TList* hijlist = (TList*)hijfile->Get("histos");\r
-  TH2F* hijmcele = (TH2F*)histos->FindObject("AnaElectron_hPtMCElectron");\r
-  TH1F* hijmchad = (TH1F*)histos->FindObject("AnaElectron_hPtMCHadron");\r
-  TH1F* hijmult = (TH1F*)histos->FindObject("AnaElectron_hRefMult");\r
-  Int_t nEvt = hijmult->GetEntries();\r
-  if(nEvt == 0) { printf("NO HIJING EVENTS\n"); return; }\r
-  hijmcele->Scale(hijscale/nEvt);\r
-  hijmchad->Scale(hijscale/nEvt);\r
-\r
-  TFile* jjfile = new TFile(jjfname);\r
-  if(!jjfile) { printf("NO JET-JET FILE\n"); return; }\r
-  TH2F* jjmcele = (TH2F*)jjfile->Get("AnaElectron_hPtMCElectronScaled");\r
-  TH1F* jjmchad = (TH1F*)jjfile->Get("AnaElectron_hPtMCHadronScaled");\r
-  TH1F* jjmult = (TH1F*)jjfile->Get("AnaElectron_hRefMultScaled");\r
-  Int_t nEvtJJ = jjmult->GetEntries();\r
-  jjmcele->Scale(pyscale);\r
-  jjmchad->Scale(pyscale);\r
-\r
-  TFile* bfile = new TFile(bfname);\r
-  if(!bfile) { printf("NO B-JET FILE\n"); return; }\r
-  TH2F* bmcele = (TH2F*)bfile->Get("AnaElectron_hPtMCElectronScaled");\r
-  TH1F* bmchad = (TH1F*)bfile->Get("AnaElectron_hPtMCHadronScaled");\r
-  TH1F* bmult = (TH1F*)bfile->Get("AnaElectron_hRefMultScaled");\r
-  Int_t nEvtB = bmult->GetEntries();\r
-  bmcele->Scale(bscale);\r
-  bmchad->Scale(bscale);\r
-\r
-  TFile* wfile = new TFile(wfname);\r
-  if(!wfile) { printf("NO W-BOSON FILE\n"); return; }\r
-  TH2F* wmcele = (TH2F*)wfile->Get("AnaElectron_hPtMCElectron");\r
-  TH1F* wmchad = (TH1F*)wfile->Get("AnaElectron_hPtMCHadron");\r
-  TH1F* wmult = (TH1F*)wfile->Get("AnaElectron_hRefMult");\r
-  Int_t nEvtW = wmult->GetEntries();\r
-  wmcele->Scale(wscale);\r
-  wmchad->Scale(wscale);\r
-\r
-  printf("Event statistics: %d (HIJING)  %d (JET-JET)  %d (B-JET)  %d (W-Boson)\n",nEvt,nEvtJJ,nEvtB,nEvtW);\r
-\r
-  TH2F* combined = (TH2F*)hijmcele->Clone();\r
-  combined->Add(jjmcele);\r
-  combined->Add(bmcele);\r
-  combined->SetTitle("MC electrons in Pb+Pb in EMCAL acceptance");\r
-  combined->SetName("CombinedMCEle");\r
-  combined->SetXTitle("p_T (GeV/c)");\r
-\r
-  mchad = (TH1F*)hijmchad->Clone();\r
-  mchad->Add(jjmchad);\r
-  mchad->Add(bmchad);\r
-  mchad->Add(wmchad);\r
-  mchad->SetTitle("MC hadrons in Pb+Pb in EMCAL acceptance");\r
-  mchad->SetName("CombinedMCHad");\r
-  mchad->SetXTitle("p_T (GeV/c)");\r
-\r
-  all = (TH1F*)combined->ProjectionX("all",1,1);\r
-  bele = (TH1F*)combined->ProjectionX("b",2,2);\r
-  cele = (TH1F*)combined->ProjectionX("c",3,3);\r
-  candb = (TH1F*)combined->ProjectionX("candb",4,4);\r
-  conv = (TH1F*)combined->ProjectionX("conv",5,5);\r
-  dalitz = (TH1F*)combined->ProjectionX("dalitz",6,6);\r
-  other = (TH1F*)combined->ProjectionX("other",8,8);\r
-\r
-  wz = (TH1F*)wmcele->ProjectionX("wz",7,7);\r
-  \r
-  all->Add(wz); //because it had to be done separately\r
-\r
-  //For comparing contributions\r
-  walle = (TH1F*)wmcele->ProjectionX("walle",7,7);\r
-  sige = (TH1F*)bmcele->ProjectionX("sige",1,1);\r
-  bkge = (TH1F*)jjmcele->ProjectionX("bkge",1,1);\r
-  hije = (TH1F*)hijmcele->ProjectionX("hije",1,1);\r
+void plotMCRates() {\r
 \r
-  double myscale = 1.; //we already scaled them\r
-  ScaleAndConfigure(all,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(bele,myscale,kRed,kFALSE);\r
-  ScaleAndConfigure(sige,myscale,kRed,kFALSE);\r
-  ScaleAndConfigure(cele,myscale,kBlue,kFALSE);\r
-  ScaleAndConfigure(candb,myscale,kViolet,kFALSE);\r
-  ScaleAndConfigure(conv,myscale,kOrange-3,kFALSE);\r
-  ScaleAndConfigure(bkge,myscale,kOrange-3,kFALSE);\r
-  ScaleAndConfigure(dalitz,myscale,kGreen-3,kFALSE);\r
-  ScaleAndConfigure(wz,myscale,kOrange-7,kFALSE);\r
-  ScaleAndConfigure(walle,myscale,kOrange-7,kFALSE);\r
-  ScaleAndConfigure(mchad,myscale,kGreen+2,kFALSE);\r
-  ScaleAndConfigure(hije,myscale,kGreen+2,kFALSE);\r
+  gROOT->LoadMacro("makeCombinedData.C");\r
+  makeData("data/scaled25Oct09/histosLHC08d6.root",\r
+           "data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root",\r
+           "data/scaled25Oct09/histosscaledLHC09b4AODc.root",\r
+           "data/scaled25Oct09/histosWboson.root");\r
 \r
   gStyle->SetOptStat(0);\r
   //drawXSRates();\r
@@ -131,43 +24,6 @@ void plotMCRates(char* hijfname = "data/scaled25Oct09/histosLHC08d6.root",
   drawHadEleRatios();\r
   drawSigBkg();\r
 \r
-  TFile* mcout = new TFile("MCElectrons.root","RECREATE");\r
-  all->Write();\r
-  bele->Write();\r
-  cele->Write();\r
-  candb->Write();\r
-  conv->Write();\r
-  dalitz->Write();\r
-  wz->Write();\r
-  other->Write();\r
-  mchad->Write();\r
-  sige->Write();\r
-  bkge->Write();\r
-  walle->Write();\r
-  hije->Write();\r
-  mcout->Close();\r
-\r
-}\r
-\r
-void ScaleAndConfigure(TH1F* hist,Double_t scale, Int_t color,Bool_t keepErr=kFALSE)\r
-{\r
-  hist->Scale(scale);\r
-  hist->SetLineColor(color);\r
-  hist->SetLineWidth(2);\r
-  if(keepErr == kFALSE) {\r
-    //remove the error bars - useful for MC rates\r
-    for(Int_t i = 1; i <= hist->GetNbinsX(); i++) {\r
-      if(hist->GetBinContent(i) > 0.) {\r
-        if(hist->GetBinError(i)/hist->GetBinContent(i) > 0.5) {\r
-          Double_t avg = 0.;\r
-          if(i > 1 && i < hist->GetNbinsX())\r
-            avg = (hist->GetBinContent(i-1) + hist->GetBinContent(i+1))/2.;\r
-          hist->SetBinContent(i,avg);\r
-        }\r
-      }\r
-      hist->SetBinError(i,0.);\r
-    }\r
-  }\r
 }\r
 \r
 void drawAnnualYields() {\r
@@ -175,28 +31,28 @@ void drawAnnualYields() {
   TCanvas* crates = new TCanvas();\r
   crates->cd();\r
   gPad->SetLogy();\r
-  all->SetXTitle("p_{T} (GeV/c)");\r
-  all->SetTitle("MC electrons in Pb+Pb, 5.5 TeV");\r
-  all->SetYTitle("Annual yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
-  all->GetYaxis()->SetRangeUser(1,2.E6);\r
-  all->GetXaxis()->SetRangeUser(10.,50.);\r
-  all->Draw();\r
-  bele->Draw("same");  \r
-  cele->Draw("same");  \r
-  candb->Draw("same");  \r
-  conv->Draw("same");  \r
-  dalitz->Draw("same");  \r
-  wz->Draw("same");\r
+  allmc->SetXTitle("p_{T} (GeV/c)");\r
+  allmc->SetTitle("MC electrons in Pb+Pb, 5.5 TeV");\r
+  allmc->SetYTitle("Annual yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
+  allmc->GetYaxis()->SetRangeUser(1,2.E6);\r
+  allmc->GetXaxis()->SetRangeUser(10.,50.);\r
+  allmc->Draw();\r
+  belemc->Draw("same");  \r
+  celemc->Draw("same");  \r
+  candbmc->Draw("same");  \r
+  convmc->Draw("same");  \r
+  dalmc->Draw("same");  \r
+  wzmc->Draw("same");\r
 \r
   leg = new TLegend(0.6,0.6,0.9,0.9);\r
   leg->SetTextSize(leg->GetTextSize()*1.2);\r
-  leg->AddEntry(all,"All MC electrons","l");\r
-  leg->AddEntry(bele,"Bottom e","l");\r
-  leg->AddEntry(cele,"Charm e","l");\r
-  leg->AddEntry(candb,"B-->C e","l");\r
-  leg->AddEntry(dalitz,"Dalitz e","l");\r
-  leg->AddEntry(conv,"Conversion e","l");\r
-  leg->AddEntry(wz,"W Boson e","l");\r
+  leg->AddEntry(allmc,"All MC electrons","l");\r
+  leg->AddEntry(belemc,"Bottom e","l");\r
+  leg->AddEntry(celemc,"Charm e","l");\r
+  leg->AddEntry(candbmc,"B-->C e","l");\r
+  leg->AddEntry(dalmc,"Dalitz e","l");\r
+  leg->AddEntry(convmc,"Conversion e","l");\r
+  leg->AddEntry(wzmc,"W Boson e","l");\r
   leg->Draw();\r
   crates->Print("MCRates_all.pdf");\r
 \r
@@ -207,24 +63,24 @@ void drawSigBkg() {
   TCanvas* csigbkg = new TCanvas();\r
   csigbkg->cd();\r
   gPad->SetLogy();\r
-  all->SetXTitle("p_{T} (GeV/c)");\r
-  all->SetTitle("MC electrons in Pb+Pb, 5.5 TeV");\r
-  all->SetYTitle("Annual Yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
-  all->GetYaxis()->SetRangeUser(1.,6.E8);\r
-  all->GetXaxis()->SetRangeUser(0.,50.);\r
-  all->Draw();\r
-  sige->Draw("same");  \r
-  bkge->Draw("same");  \r
-  hije->Draw("same");  \r
-  walle->Draw("same");\r
+  allmc->SetXTitle("p_{T} (GeV/c)");\r
+  allmc->SetTitle("MC electrons in Pb+Pb, 5.5 TeV");\r
+  allmc->SetYTitle("Annual Yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
+  allmc->GetYaxis()->SetRangeUser(1.,6.E8);\r
+  allmc->GetXaxis()->SetRangeUser(0.,50.);\r
+  allmc->Draw();\r
+  sigemc->Draw("same");  \r
+  bkgemc->Draw("same");  \r
+  hijemc->Draw("same");  \r
+  wallemc->Draw("same");\r
 \r
   TLegend* leg1 = new TLegend(0.6,0.6,0.9,0.9);\r
   leg1->SetTextSize(leg->GetTextSize()*1.2);\r
-  leg1->AddEntry(all,"All MC electrons","l");\r
-  leg1->AddEntry(sige,"B-Jet Events","l");\r
-  leg1->AddEntry(hije,"Pb+Pb Underlying Event","l");\r
-  leg1->AddEntry(bkge,"Jet-Jet Events","l");\r
-  leg1->AddEntry(walle,"W-decay Events","l");\r
+  leg1->AddEntry(allmc,"All MC electrons","l");\r
+  leg1->AddEntry(sigemc,"B-Jet Events","l");\r
+  leg1->AddEntry(hijemc,"Pb+Pb Underlying Event","l");\r
+  leg1->AddEntry(bkgemc,"Jet-Jet Events","l");\r
+  leg1->AddEntry(wallemc,"W-decay Events","l");\r
   leg1->Draw();\r
   csigbkg->Print("MCRates_byEventSource.pdf");\r
 \r
@@ -235,13 +91,13 @@ void drawPtCutRates() {
   TCanvas* cptcut = new TCanvas();\r
   cptcut->cd();\r
   gPad->SetLogy();\r
-  TH1F* alleptcut = GetPtCutHisto(all);\r
-  TH1F* beleptcut = GetPtCutHisto(bele);\r
-  TH1F* celeptcut = GetPtCutHisto(cele);\r
-  TH1F* cbeleptcut = GetPtCutHisto(candb);\r
-  TH1F* dalitzptcut = GetPtCutHisto(dalitz);\r
-  TH1F* convptcut = GetPtCutHisto(conv);\r
-  TH1F* wzptcut = GetPtCutHisto(wz);\r
+  TH1F* alleptcut = GetPtCutHisto(allmc);\r
+  TH1F* beleptcut = GetPtCutHisto(belemc);\r
+  TH1F* celeptcut = GetPtCutHisto(celemc);\r
+  TH1F* cbeleptcut = GetPtCutHisto(candbmc);\r
+  TH1F* dalitzptcut = GetPtCutHisto(dalmc);\r
+  TH1F* convptcut = GetPtCutHisto(convmc);\r
+  TH1F* wzptcut = GetPtCutHisto(wzmc);\r
   alleptcut->GetXaxis()->SetRangeUser(10,50);\r
   alleptcut->GetYaxis()->SetRangeUser(10,2.e6);\r
   alleptcut->SetXTitle("p_{T}^{cut} (GeV/c)");\r
@@ -265,14 +121,14 @@ void drawHadEleRatios() {
   ceh->cd();\r
   gPad->SetLogy();\r
   gStyle->SetOptStat(0);\r
-  TH1F* allratio = (TH1F*)all->Clone();\r
-  TH1F* behratio = (TH1F*)bele->Clone();\r
+  TH1F* allratio = (TH1F*)allmc->Clone();\r
+  TH1F* behratio = (TH1F*)belemc->Clone();\r
   allratio->SetTitle("MC hadrons and electrons in Pb+Pb, 5.5 TeV");\r
   allratio->SetXTitle("p_{T} (GeV/c)");\r
   allratio->SetYTitle("Hadrons/Electrons");\r
-  for(Int_t i = 1; i < all->GetNbinsX(); i++) {\r
-    Double_t vale = all->GetBinContent(i);\r
-    Double_t valb = bele->GetBinContent(i);\r
+  for(Int_t i = 1; i < allmc->GetNbinsX(); i++) {\r
+    Double_t vale = allmc->GetBinContent(i);\r
+    Double_t valb = belemc->GetBinContent(i);\r
     Double_t valh = mchad->GetBinContent(i);\r
     //printf("pT %.2f, Hadron %.1f, Electron %.1f, B-electron %.1f\n",all->GetBinCenter(i),valh,vale,valb);\r
     if(vale>0) allratio->SetBinContent(i,valh/vale);\r
@@ -303,22 +159,3 @@ void drawHadEleRatios() {
   ceh->Print("MCRates_heratio.pdf");\r
 }\r
 \r
-TH1F* GetPtCutHisto(TH1F* input) \r
-{\r
-  //Given a rate histogram vs pt, return the histogram with yield\r
-  //above a given pTcut\r
-\r
-  TH1F* result = (TH1F*)input->Clone();\r
-  char name[100];\r
-  sprintf(name,"%s_ptCut",result->GetName());\r
-  result->SetNameTitle(name,name);\r
-  for(Int_t i = 1; i <= result->GetNbinsX(); i++) {\r
-    Double_t val = input->Integral(i,result->GetNbinsX());\r
-    result->SetBinContent(i,val);\r
-    result->SetBinError(i,0.);\r
-  }\r
-\r
-  return result;\r
-\r
-}\r
-\r
index 1ea1ccfc08ec9411867a6c98685e6793ba1dcd7d..6c480e1dd97e5785c08af9ebab72b1cb8b5b9bb5 100644 (file)
 //\r
 /////////////////////////////////////////////////\r
 \r
-TH1F *alltte,      *alltrk,      *allemc;\r
-TH1F *sumtte,      *sumtrk,      *sumemc;\r
-TH1F *sumHFemc,    *sumHFmc;\r
-TH1F *bpttte,      *bpttrk,      *bptemc;\r
-TH1F *cpttte,      *cpttrk,      *cptemc;\r
-TH1F *candbpttte,  *candbpttrk,  *candbptemc;\r
-TH1F *convpttte,   *convpttrk,   *convptemc;\r
-TH1F *dalitzpttte, *dalitzpttrk, *dalitzptemc;\r
-TH1F *wzpttte,     *wzpttrk,     *wzptemc;\r
-TH1F *otherpttte,  *otherpttrk,  *otherptemc;\r
-TH1F *misidpttte,  *misidpttrk,  *misidptemc;\r
-\r
 TLegend* leg;\r
 \r
-void plotNPERates(char* which = "TTE", \r
-                 char* hijfname = "data/scaled25Oct09/histosLHC08d6.root",\r
-                 char* jjfname = "data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root",\r
-                 char* bfname = "data/scaled25Oct09/histosscaledLHC09b4AODc.root",\r
-                 char* wfname = "data/scaled25Oct09/histosWboson.root") {\r
-\r
-  //For HIJING need to divide by the number of events, which we\r
-  //can get from the file and do when we perform scaling\r
-  double hijscale = 0.05*(1.E6)*0.5*7700; //0-5% * seconds*lumi*PbPb x-section\r
-  //For bjet and jet-jet events\r
-  double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance\r
-  double bscale = pyscale; //Do we need to scale by Branching ratio for forced\r
-                               //semi-leptonic decays?\r
-  double wscale = pyscale;\r
-  \r
-  TFile* hijfile = new TFile(hijfname);\r
-  if(!hijfile) { printf("NO HIJING FILE\n"); return; }\r
-  TList* hijlist = (TList*)hijfile->Get("histos");\r
-  TH2F* hijtte = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleTTE");\r
-  TH2F* hijemc = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleEMCAL");\r
-  TH2F* hijtrk = (TH2F*)histos->FindObject("AnaElectron_hPtNPEleTPCTRD");\r
-  TH1F* hijmult = (TH1F*)histos->FindObject("AnaElectron_hRefMult");\r
-  Int_t nEvt = hijmult->GetEntries();\r
-  if(nEvt == 0) { printf("NO HIJING EVENTS\n"); return; }\r
-  hijtte->Scale(hijscale/nEvt);\r
-  hijemc->Scale(hijscale/nEvt);\r
-  hijtrk->Scale(hijscale/nEvt);\r
-\r
-  TFile* jjfile = new TFile(jjfname);\r
-  if(!jjfile) { printf("NO JET-JET FILE\n"); return; }\r
-  TH2F* jjtte = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleTTEScaled");\r
-  TH2F* jjemc = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleEMCALScaled");\r
-  TH2F* jjtrk = (TH2F*)jjfile->Get("AnaElectron_hPtNPEleTPCTRDScaled");\r
-  TH1F* jjmult = (TH1F*)jjfile->Get("AnaElectron_hRefMultScaled");\r
-  Int_t nEvtJJ = jjmult->GetEntries();\r
-  jjtte->Scale(pyscale);\r
-  jjemc->Scale(pyscale);\r
-  jjtrk->Scale(pyscale);\r
-\r
-  TFile* bfile = new TFile(bfname);\r
-  if(!bfile) { printf("NO B-JET FILE\n"); return; }\r
-  TH2F* btte = (TH2F*)bfile->Get("AnaElectron_hPtNPEleTTEScaled");\r
-  TH2F* bemc = (TH2F*)bfile->Get("AnaElectron_hPtNPEleEMCALScaled");\r
-  TH2F* btrk = (TH2F*)bfile->Get("AnaElectron_hPtNPEleTPCTRDScaled");\r
-  TH1F* bmult = (TH1F*)bfile->Get("AnaElectron_hRefMultScaled");\r
-  Int_t nEvtB = bmult->GetEntries();\r
-  btte->Scale(bscale);\r
-  bemc->Scale(bscale);\r
-  btrk->Scale(bscale);\r
-\r
-  TFile* wfile = new TFile(wfname);\r
-  if(!wfile) { printf("NO W-BOSON FILE\n"); return; }\r
-  TH2F* wtte = (TH2F*)wfile->Get("AnaElectron_hPtNPEleTTE");\r
-  TH2F* wemc = (TH2F*)wfile->Get("AnaElectron_hPtNPEleEMCAL");\r
-  TH2F* wtrk = (TH2F*)wfile->Get("AnaElectron_hPtNPEleTPCTRD");\r
-  TH1F* wmult = (TH1F*)wfile->Get("AnaElectron_hRefMult");\r
-  Int_t nEvtW = wmult->GetEntries();\r
-  wtte->Scale(wscale); //already scaled by evts\r
-  wemc->Scale(wscale); //already scaled by evts\r
-  wtrk->Scale(wscale); //already scaled by evts\r
-\r
-  printf("Event statistics: %d (HIJING)  %d (JET-JET)  %d (B-JET)  %d (W-Boson)\n",nEvt,nEvtJJ,nEvtB,nEvtW);\r
-\r
-  TH2F* combTTE = (TH2F*)hijtte->Clone();\r
-  combTTE->Add(jjtte);\r
-  combTTE->Add(btte);\r
-  combTTE->SetTitle("Identified non-phot. electrons (TPC+TRD+EMCAL)");\r
-  combTTE->SetName("CombinedEleTTE");\r
-  combTTE->SetXTitle("p_T (GeV/c)");\r
-\r
-  alltte = (TH1F*)combTTE->ProjectionX("alltte",1,1);\r
-  bpttte = (TH1F*)combTTE->ProjectionX("btte",2,2);\r
-  cpttte = (TH1F*)combTTE->ProjectionX("ctte",3,3);\r
-  candbpttte = (TH1F*)combTTE->ProjectionX("candbtte",4,4);\r
-  convpttte = (TH1F*)combTTE->ProjectionX("convtte",5,5);\r
-  dalitzpttte = (TH1F*)combTTE->ProjectionX("dalitztte",6,6);\r
-  otherpttte = (TH1F*)combTTE->ProjectionX("othertte",8,8);\r
-  sumtte = (TH1F*)bpttte->Clone(); sumtte->SetName("sumtte");\r
-  sumtte->Add(cpttte); sumtte->Add(candbpttte); sumtte->Add(convpttte);\r
-  sumtte->Add(dalitzpttte); //sumtte->Add(otherpttte);\r
-  misidpttte = (TH1F*)combTTE->ProjectionX("misidtte",9,9);\r
-\r
-  wzpttte = (TH1F*)wtte->ProjectionX("wz",7,7);\r
-  alltte->Add(wzpttte);\r
-  sumtte->Add(wzpttte);\r
-\r
-  double myscale = 1.; //we already scaled them\r
-  ScaleAndConfigure(alltte,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(sumtte,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(bpttte,myscale,kRed,kFALSE);\r
-  ScaleAndConfigure(cpttte,myscale,kBlue,kFALSE);\r
-  ScaleAndConfigure(candbpttte,myscale,kViolet,kFALSE);\r
-  ScaleAndConfigure(convpttte,myscale,kOrange-3,kFALSE);\r
-  ScaleAndConfigure(dalitzpttte,myscale,kGreen-3,kFALSE);\r
-  ScaleAndConfigure(misidpttte,myscale,kGreen+2,kFALSE);\r
-  ScaleAndConfigure(wzpttte,myscale,kOrange-7,kFALSE);\r
-\r
-  TH2F* combEMC = (TH2F*)hijemc->Clone();\r
-  combEMC->Add(jjemc);\r
-  combEMC->Add(bemc);\r
-  combEMC->SetTitle("Identified non-phot. electrons (EMCAL)");\r
-  combEMC->SetName("CombinedEleEMC");\r
-  combEMC->SetXTitle("p_T (GeV/c)");\r
+void plotNPERates(const char* which = "TTE") {\r
 \r
-  allemc = (TH1F*)combEMC->ProjectionX("allemc",1,1);\r
-  bptemc = (TH1F*)combEMC->ProjectionX("bemc",2,2);\r
-  cptemc = (TH1F*)combEMC->ProjectionX("cemc",3,3);\r
-  candbptemc = (TH1F*)combEMC->ProjectionX("candbemc",4,4);\r
-  convptemc = (TH1F*)combEMC->ProjectionX("convemc",5,5);\r
-  dalitzptemc = (TH1F*)combEMC->ProjectionX("dalitzemc",6,6);\r
-  otherptemc = (TH1F*)combEMC->ProjectionX("otheremc",8,8);\r
-  misidptemc = (TH1F*)combEMC->ProjectionX("misidemc",9,9);\r
-\r
-  wzptemc = (TH1F*)wemc->ProjectionX("wz",7,7);\r
-  allemc->Add(wzptemc);\r
-  sumemc = (TH1F*)bptemc->Clone(); sumemc->SetName("sumemc");\r
-  sumemc->Add(cptemc); sumemc->Add(candbptemc); sumemc->Add(convptemc);\r
-  sumemc->Add(dalitzptemc); //sumemc->Add(otherptemc);\r
-  sumemc->Add(wzptemc);\r
-  sumHFemc = (TH1F*)bptemc->Clone(); sumHFemc->SetName("sumHFemc");\r
-  sumHFemc->Add(cptemc); sumHFemc->Add(candbptemc); sumHFemc->Add(wzptemc);\r
-\r
-  double myscale = 1.; //we already scaled them\r
-  ScaleAndConfigure(allemc,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(sumemc,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(sumHFemc,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(bptemc,myscale,kRed,kFALSE);\r
-  ScaleAndConfigure(cptemc,myscale,kBlue,kFALSE);\r
-  ScaleAndConfigure(candbptemc,myscale,kViolet,kFALSE);\r
-  ScaleAndConfigure(convptemc,myscale,kOrange-3,kFALSE);\r
-  ScaleAndConfigure(dalitzptemc,myscale,kGreen-3,kFALSE);\r
-  ScaleAndConfigure(misidptemc,myscale,kGreen+2,kFALSE);\r
-  ScaleAndConfigure(wzptemc,myscale,kOrange-7,kFALSE);\r
-\r
-  TH2F* combTRK = (TH2F*)hijtrk->Clone();\r
-  combTRK->Add(jjtrk);\r
-  combTRK->Add(btrk);\r
-  combTRK->SetTitle("Identified non-phot. electrons (TPC+TRD)");\r
-  combTRK->SetName("CombinedEleTRK");\r
-  combTRK->SetXTitle("p_T (GeV/c)");\r
-\r
-  alltrk = (TH1F*)combTRK->ProjectionX("alltrk",1,1);\r
-  bpttrk = (TH1F*)combTRK->ProjectionX("btrk",2,2);\r
-  cpttrk = (TH1F*)combTRK->ProjectionX("ctrk",3,3);\r
-  candbpttrk = (TH1F*)combTRK->ProjectionX("candbtrk",4,4);\r
-  convpttrk = (TH1F*)combTRK->ProjectionX("convtrk",5,5);\r
-  dalitzpttrk = (TH1F*)combTRK->ProjectionX("dalitztrk",6,6);\r
-  otherpttrk = (TH1F*)combTRK->ProjectionX("othertrk",8,8);\r
-  sumtrk = (TH1F*)bpttrk->Clone(); sumtrk->SetName("sumtrk");\r
-  sumtrk->Add(cpttrk); sumtrk->Add(candbpttrk); sumtrk->Add(convpttrk);\r
-  sumtrk->Add(dalitzpttrk); //sumtrk->Add(otherpttrk);\r
-  misidpttrk = (TH1F*)combTRK->ProjectionX("misidtrk",9,9);\r
-\r
-  wzpttrk = (TH1F*)wtrk->ProjectionX("wztrk",7,7);\r
-  alltrk->Add(wzpttrk);\r
-  sumtrk->Add(wzpttrk);\r
-\r
-  double myscale = 1.; //we already scaled them\r
-  ScaleAndConfigure(alltrk,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(sumtrk,myscale,kBlack,kFALSE);\r
-  ScaleAndConfigure(bpttrk,myscale,kRed,kFALSE);\r
-  ScaleAndConfigure(cpttrk,myscale,kBlue,kFALSE);\r
-  ScaleAndConfigure(candbpttrk,myscale,kViolet,kFALSE);\r
-  ScaleAndConfigure(convpttrk,myscale,kOrange-3,kFALSE);\r
-  ScaleAndConfigure(dalitzpttrk,myscale,kGreen-3,kFALSE);\r
-  ScaleAndConfigure(misidpttrk,myscale,kGreen+2,kFALSE);\r
-  ScaleAndConfigure(wzpttrk,myscale,kOrange-7,kFALSE);\r
+  gROOT->LoadMacro("makeCombinedData.C");\r
+  makeData("data/scaled25Oct09/histosLHC08d6.root",\r
+           "data/scaled25Oct09/TOTALhistosscaled-LHC09b2-0.root",\r
+           "data/scaled25Oct09/histosscaledLHC09b4AODc.root",\r
+           "data/scaled25Oct09/histosWboson.root");\r
 \r
   //define common legend\r
   leg = new TLegend(0.5,0.6,0.9,0.9);\r
   leg->SetTextSize(leg->GetTextSize()*1.2);\r
   //leg->AddEntry(alltte,"All N-P e candidates","l");\r
   leg->AddEntry(sumtte,"All N-P electrons","l");\r
-  leg->AddEntry(bpttte,"Bottom e","l");\r
-  leg->AddEntry(cpttte,"Charm e","l");\r
-  leg->AddEntry(candbpttte,"B-->C e","l");\r
-  leg->AddEntry(dalitzpttte,"Dalitz e","l");\r
-  leg->AddEntry(convpttte,"Conversion e","l");\r
-  leg->AddEntry(wzpttte,"W Boson e","l");\r
-  //leg->AddEntry(misidpttte,"Mis-identified hadrons","l");\r
+  leg->AddEntry(btte,"Bottom e","l");\r
+  leg->AddEntry(ctte,"Charm e","l");\r
+  leg->AddEntry(cbtte,"B-->C e","l");\r
+  leg->AddEntry(daltte,"Dalitz e","l");\r
+  leg->AddEntry(convtte,"Conversion e","l");\r
+  leg->AddEntry(wztte,"W Boson e","l");\r
 \r
   gStyle->SetOptStat(0);\r
   drawAnnualYields(which);\r
-  drawComparePID();\r
   drawPtCutRates(which);\r
   drawCompareTruth();\r
 \r
 }\r
 \r
-void ScaleAndConfigure(TH1F* hist,Double_t scale, Int_t color,Bool_t keepErr)\r
-{\r
-  hist->Scale(scale);\r
-  hist->SetLineColor(color);\r
-  hist->SetMarkerColor(color);\r
-  hist->SetLineWidth(2);\r
-  if(keepErr == kFALSE) {\r
-    //remove the error bars - useful for MC rates\r
-    for(Int_t i = 1; i <= hist->GetNbinsX(); i++) {\r
-      if(hist->GetBinContent(i) > 0.) {\r
-       if(hist->GetBinError(i)/hist->GetBinContent(i) > 0.5) {\r
-         Double_t avg = 0.;\r
-         if(i > 1 && i < hist->GetNbinsX()) \r
-           avg = (hist->GetBinContent(i-1) + hist->GetBinContent(i+1))/2.;\r
-         hist->SetBinContent(i,avg);\r
-       }\r
-      }\r
-      hist->SetBinError(i,0.);\r
-    }\r
-  } else {\r
-    //Set the error bars to statistics of the bin content\r
-    for(Int_t i = 1; i <= hist->GetNbinsX(); i++) {\r
-      if(hist->GetBinContent(i) > 0.) {\r
-       if(hist->GetBinError(i)/hist->GetBinContent(i) > 0.5) {\r
-         Double_t avg = 0;\r
-         if(i > 1 && i < hist->GetNbinsX()) \r
-           avg = (hist->GetBinContent(i-1) + hist->GetBinContent(i+1))/2.;\r
-         hist->SetBinContent(i,avg);\r
-       }\r
-      }\r
-      hist->SetBinError(i,TMath::Sqrt(hist->GetBinContent(i)));\r
-    }\r
-  }\r
-  \r
-}\r
-\r
 void drawAnnualYields(char* which = "EMC") {\r
 \r
   TCanvas* crates = new TCanvas();\r
@@ -264,13 +54,13 @@ void drawAnnualYields(char* which = "EMC") {
     sumemc->SetTitle("Annual yield of non-phot. electrons (EMCAL pid)");\r
     sumemc->Rebin(2); sumemc->Scale(0.5);\r
 \r
-    bptemc->Rebin(2); bptemc->Scale(0.5);\r
-    cptemc->Rebin(2); cptemc->Scale(0.5);\r
-    candbptemc->Rebin(2); candbptemc->Scale(0.5);\r
-    convptemc->Rebin(2); convptemc->Scale(0.5);\r
-    dalitzptemc->Rebin(2); dalitzptemc->Scale(0.5);\r
-    wzptemc->Rebin(2); wzptemc->Scale(0.5);\r
-    misidptemc->Rebin(2); misidptemc->Scale(0.5);\r
+    bemc->Rebin(2); bemc->Scale(0.5);\r
+    cemc->Rebin(2); cemc->Scale(0.5);\r
+    cbemc->Rebin(2); cbemc->Scale(0.5);\r
+    convemc->Rebin(2); convemc->Scale(0.5);\r
+    dalemc->Rebin(2); dalemc->Scale(0.5);\r
+    wzemc->Rebin(2); wzemc->Scale(0.5);\r
+    hemc->Rebin(2); hemc->Scale(0.5);\r
     /*\r
     allemc->GetYaxis()->SetRangeUser(1.,2.e6);\r
     allemc->GetXaxis()->SetRangeUser(10.,49.);\r
@@ -279,13 +69,13 @@ void drawAnnualYields(char* which = "EMC") {
     sumemc->GetYaxis()->SetRangeUser(1.,2.e6);\r
     sumemc->GetXaxis()->SetRangeUser(10.,49.);\r
     sumemc->Draw();\r
-    bptemc->Draw("same");\r
-    cptemc->Draw("same");\r
-    candbptemc->Draw("same");\r
-    convptemc->Draw("same");\r
-    dalitzptemc->Draw("same");\r
-    wzptemc->Draw("same");\r
-    //misidptemc->Draw("same");\r
+    bemc->Draw("same");\r
+    cemc->Draw("same");\r
+    cbemc->Draw("same");\r
+    convemc->Draw("same");\r
+    dalemc->Draw("same");\r
+    wzemc->Draw("same");\r
+    //hemc->Draw("same");\r
     leg->Draw();\r
     crates->Print("NPERates_EMC_all.pdf");\r
   }\r
@@ -299,13 +89,13 @@ void drawAnnualYields(char* which = "EMC") {
     sumtte->SetYTitle("Annual yield");\r
     sumtte->SetTitle("Annual yield of non-phot. electrons (Tracking+EMCAL pid)");\r
     sumtte->Rebin(2); sumtte->Scale(0.5);\r
-    bpttte->Rebin(2); bpttte->Scale(0.5);\r
-    cpttte->Rebin(2); cpttte->Scale(0.5);\r
-    candbpttte->Rebin(2); candbpttte->Scale(0.5);\r
-    convpttte->Rebin(2); convpttte->Scale(0.5);\r
-    dalitzpttte->Rebin(2); dalitzpttte->Scale(0.5);\r
-    wzpttte->Rebin(2); wzpttte->Scale(0.5);\r
-    misidpttte->Rebin(2); misidpttte->Scale(0.5);\r
+    btte->Rebin(2); btte->Scale(0.5);\r
+    ctte->Rebin(2); ctte->Scale(0.5);\r
+    cbtte->Rebin(2); cbtte->Scale(0.5);\r
+    convtte->Rebin(2); convtte->Scale(0.5);\r
+    daltte->Rebin(2); daltte->Scale(0.5);\r
+    wztte->Rebin(2); wztte->Scale(0.5);\r
+    htte->Rebin(2); htte->Scale(0.5);\r
     /*    alltte->GetYaxis()->SetRangeUser(1.,2.e6);\r
     alltte->GetXaxis()->SetRangeUser(10.,49.);\r
     alltte->Draw();\r
@@ -313,95 +103,17 @@ void drawAnnualYields(char* which = "EMC") {
     sumtte->GetYaxis()->SetRangeUser(1.,2.e6);\r
     sumtte->GetXaxis()->SetRangeUser(10.,49.);\r
     sumtte->Draw();\r
-    bpttte->Draw("same");\r
-    cpttte->Draw("same");\r
-    candbpttte->Draw("same");\r
-    convpttte->Draw("same");\r
-    dalitzpttte->Draw("same");\r
-    wzpttte->Draw("same");\r
-    //misidpttte->Draw("same");\r
+    btte->Draw("same");\r
+    ctte->Draw("same");\r
+    cbtte->Draw("same");\r
+    convtte->Draw("same");\r
+    daltte->Draw("same");\r
+    wztte->Draw("same");\r
+    //htte->Draw("same");\r
     leg->Draw();\r
     crates->Print("NPERates_TTE_all.pdf");\r
   }\r
 \r
-\r
-}\r
-\r
-void drawComparePID() {\r
-\r
-  TCanvas* crates = new TCanvas("crates","",0,0,1600,600);\r
-  crates->Divide(2,1);\r
-  crates->cd(1);\r
-  gPad->SetLogy();\r
-  alltrk->SetXTitle("p_{T} (GeV/c)");\r
-  alltrk->SetYTitle("Annual yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
-  alltrk->SetTitle("PID comparison: Tracking only vs. EMCAL only");\r
-  alltrk->Smooth(2);\r
-  alltrk->Rebin(2); alltrk->Scale(0.5);\r
-  alltrk->GetYaxis()->SetRangeUser(10.,alltrk->GetMaximum()*2.);\r
-  alltrk->GetXaxis()->SetRangeUser(10.,49.);\r
-  alltrk->Draw();\r
-  misidpttrk->Smooth(2);\r
-  misidpttrk->Rebin(2); misidpttrk->Scale(0.5);\r
-  misidpttrk->Draw("same");\r
-  TH1F* tempallemc = (TH1F*)allemc->Clone();\r
-  tempallemc->SetNameTitle("tempallemc","tempallemc");\r
-  tempallemc->SetLineColor(kBlue);\r
-  tempallemc->Rebin(2); tempallemc->Scale(0.5);\r
-  tempallemc->Draw("same");\r
-  \r
-  TH1F* tempmisidemc = (TH1F*)misidptemc->Clone();\r
-  tempmisidemc->SetNameTitle("tempmisidemc","tempmisidemc");\r
-  tempmisidemc->SetLineColor(kOrange-3);\r
-  tempmisidemc->Rebin(2); tempmisidemc->Scale(0.5);\r
-  tempmisidemc->Draw("same");\r
-\r
-  TLegend* leg2 = new TLegend(0.35,0.6,0.9,0.9);\r
-  leg2->SetTextSize(leg->GetTextSize()*1.2);\r
-  leg2->AddEntry(alltrk,"Electron Candidates (Tracking PID only)","l");\r
-  leg2->AddEntry(misidpttrk,"Hadron Contamination (Tracking PID only)","l");\r
-  leg2->AddEntry(tempallemc,"Electron Candidates (EMCAL PID)","l");\r
-  leg2->AddEntry(tempmisidemc,"Hadron Contamination (EMCAL PID)","l");\r
-  leg2->Draw();\r
-\r
-  crates->cd(2);\r
-  gPad->SetLogy();\r
-  TH1F* subtrk = (TH1F*)alltrk->Clone();\r
-  for(Int_t i = 1; i <= alltrk->GetNbinsX(); i++) {\r
-    Double_t diff = alltrk->GetBinContent(i) - misidpttrk->GetBinContent(i);\r
-    Double_t unc = diff/TMath::Sqrt(diff+2.*misidpttrk->GetBinContent(i));\r
-    //    printf("Cand %d, Contam %d, diff %d, unc %d\n",alltrk->GetBinContent(i),misidpttrk->GetBinContent(i), diff, unc);\r
-    subtrk->SetBinContent(i,diff);\r
-    subtrk->SetBinError(i,unc);\r
-  }\r
-  subtrk->SetYTitle("Annual yield in EMCAL dN/dp_{T} (GeV/c)^{-1}");\r
-  subtrk->SetLineColor(kRed);\r
-  subtrk->Draw();\r
-  //  sumtrk->Rebin(2); sumtrk->Scale(0.5);\r
-  //sumtrk->Draw("same");\r
-\r
-  TH1F* subemc = (TH1F*)tempallemc->Clone();\r
-  for(Int_t i = 1; i <= tempallemc->GetNbinsX(); i++) {\r
-    Double_t diff = tempallemc->GetBinContent(i) - tempmisidemc->GetBinContent(i);\r
-    Double_t unc = 0.;\r
-    if(diff > 0.) \r
-      unc = diff/TMath::Sqrt(diff+2.*tempmisidemc->GetBinContent(i));\r
-    //printf("Cand %d, Contam %d, diff %d, unc %d\n",tempallemc->GetBinContent(i),tempmisidemc->GetBinContent(i), diff, unc);\r
-    subemc->SetBinContent(i,diff);\r
-    subemc->SetBinError(i,unc);\r
-  }\r
-  subemc->SetLineColor(kCyan);\r
-  subemc->Draw("same");\r
-  TLegend *legx = new TLegend(0.3,0.7,0.9,0.9);\r
-  legx->AddEntry(subtrk,"Signal (candidates - contam) (Tracking PID only)","l");\r
-  legx->AddEntry(subemc,"Signal (candidates - contam) (EMCAL PID)","l");\r
-  legx->Draw();\r
-\r
-  TLatex* latex = new TLatex(0.5,0.6,"Unc = #frac{S}{#sqrt{S+2B}}");\r
-  latex->SetNDC();\r
-  latex->Draw();\r
-\r
-  crates->Print("NPERates_PIDCompare_all.pdf");\r
 }\r
 \r
 void drawPtCutRates(char* which = "EMC") {\r
@@ -412,13 +124,13 @@ void drawPtCutRates(char* which = "EMC") {
   if(strcmp(which,"EMC")==0) {\r
     //    TH1F* alleptcut = GetPtCutHisto(allemc);\r
     TH1F* alleptcut = GetPtCutHisto(sumemc);\r
-    TH1F* beleptcut = GetPtCutHisto(bptemc);\r
-    TH1F* celeptcut = GetPtCutHisto(cptemc);\r
-    TH1F* cbeleptcut = GetPtCutHisto(candbptemc);\r
-    TH1F* dalitzptcut = GetPtCutHisto(dalitzptemc);\r
-    TH1F* convptcut = GetPtCutHisto(convptemc);\r
-    TH1F* wzptcut = GetPtCutHisto(wzptemc);\r
-    TH1F* misidptcut = GetPtCutHisto(misidptemc);\r
+    TH1F* beleptcut = GetPtCutHisto(bemc);\r
+    TH1F* celeptcut = GetPtCutHisto(cemc);\r
+    TH1F* cbeleptcut = GetPtCutHisto(cbemc);\r
+    TH1F* dalitzptcut = GetPtCutHisto(dalemc);\r
+    TH1F* convptcut = GetPtCutHisto(convemc);\r
+    TH1F* wzptcut = GetPtCutHisto(wzemc);\r
+    TH1F* misidptcut = GetPtCutHisto(hemc);\r
     alleptcut->SetXTitle("p_{T}^{cut} (GeV/c)");\r
     alleptcut->SetYTitle("Annual Yield in EMCAL for p_{T}>p_{T}^{cut}");\r
     alleptcut->SetTitle("Pb+Pb, 5.5 TeV reconstructed N-P electrons (EMCAL pid)");\r
@@ -438,13 +150,13 @@ void drawPtCutRates(char* which = "EMC") {
   if(strcmp(which,"TTE")==0) {\r
     //    TH1F* alleptcut = GetPtCutHisto(alltte);\r
     TH1F* alleptcut = GetPtCutHisto(sumtte);\r
-    TH1F* beleptcut = GetPtCutHisto(bpttte);\r
-    TH1F* celeptcut = GetPtCutHisto(cpttte);\r
-    TH1F* cbeleptcut = GetPtCutHisto(candbpttte);\r
-    TH1F* dalitzptcut = GetPtCutHisto(dalitzpttte);\r
-    TH1F* convptcut = GetPtCutHisto(convpttte);\r
-    TH1F* wzptcut = GetPtCutHisto(wzpttte);\r
-    TH1F* misidptcut = GetPtCutHisto(misidpttte);\r
+    TH1F* beleptcut = GetPtCutHisto(btte);\r
+    TH1F* celeptcut = GetPtCutHisto(ctte);\r
+    TH1F* cbeleptcut = GetPtCutHisto(cbtte);\r
+    TH1F* dalitzptcut = GetPtCutHisto(daltte);\r
+    TH1F* convptcut = GetPtCutHisto(convtte);\r
+    TH1F* wzptcut = GetPtCutHisto(wztte);\r
+    TH1F* misidptcut = GetPtCutHisto(htte);\r
     alleptcut->SetXTitle("p_{T}^{cut} (GeV/c)");\r
     alleptcut->SetYTitle("Annual Yield in EMCAL for p_{T}>p_{T}^{cut}");\r
     alleptcut->SetTitle("Pb+Pb, 5.5 TeV reconstructed N-P electrons (Tracking+EMCAL pid)");\r
@@ -466,15 +178,13 @@ void drawPtCutRates(char* which = "EMC") {
 }\r
 \r
 TH1F* drawCompareTruth() {\r
-  TFile* f = new TFile("MCElectrons.root");\r
-  if(!f) { printf("No MC electron file.  exiting\n"); return;}\r
-  TH1F* mctruth = (TH1F*)f->Get("b");\r
-  TH1F* temp = (TH1F*)f->Get("c");\r
-  mctruth->Add(temp);\r
-  temp = (TH1F*)f->Get("candb");\r
-  mctruth->Add(temp);\r
-  temp = (TH1F*)f->Get("wz");\r
-  mctruth->Add(temp);\r
+\r
+  TH1F* mctruth = (TH1F*)belemc->Clone();\r
+  mctruth->SetName("mctruth");\r
+  mctruth->Add(celemc);\r
+  mctruth->Add(candbmc);\r
+  mctruth->Add(wzmc);\r
+  mctruth->Rebin(2); mctruth->Scale(0.5);\r
 \r
   TFile* effic = new TFile("elec_eff.root");\r
   TH1F* heff = (TH1F*)effic->Get("h111");\r
@@ -488,6 +198,18 @@ TH1F* drawCompareTruth() {
     if(eff > 0.) corr = hcorr->GetBinContent(i)/eff;\r
     hcorr->SetBinContent(i,corr);\r
   }\r
+  hcorr->Rebin(2); hcorr->Scale(0.5);\r
+  sumHFemc->Rebin(2); sumHFemc->Scale(0.5);\r
+\r
+  Double_t efinal = 0.258;\r
+  TGraphErrors* eerr = new TGraphErrors();\r
+  eerr->SetName("emcErr");\r
+  for(Int_t i = 1; i <= hcorr->GetNbinsX(); i++) {\r
+    eerr->SetPoint(i-1,hcorr->GetBinCenter(i),hcorr->GetBinContent(i));\r
+    eerr->SetPointError(i-1,0.,efinal*hcorr->GetBinContent(i));\r
+  }\r
+  eerr->SetFillColor(kRed-8);\r
+\r
   TCanvas* ctruth = new TCanvas();\r
   ctruth->cd();\r
   gPad->SetLogy();\r
@@ -497,6 +219,7 @@ TH1F* drawCompareTruth() {
   mctruth->SetXTitle("p_{T} (GeV/c)");\r
   mctruth->GetXaxis()->SetRangeUser(10.,49.);\r
   mctruth->Draw();\r
+  eerr->Draw("3same");\r
   hcorr->SetMarkerColor(kRed);\r
   hcorr->SetLineColor(kRed);\r
   hcorr->Draw("same");\r
@@ -507,28 +230,9 @@ TH1F* drawCompareTruth() {
   legy->AddEntry(mctruth,"MC HF+W electrons","l");\r
   legy->AddEntry(sumHFemc,"Rec (EMCAL) HF+W electrons","l");\r
   legy->AddEntry(hcorr,"Eff. corrected reco HF+W electrons","l");\r
+  legy->AddEntry(eerr,"Systematic uncertainty","f");\r
   legy->Draw();\r
 \r
   ctruth->Print("NPERates_TruthComparison.pdf");\r
 \r
 }\r
-\r
-TH1F* GetPtCutHisto(TH1F* input) \r
-{\r
-  //Given a rate histogram vs pt, return the histogram with yield\r
-  //above a given pTcut\r
-\r
-  TH1F* result = (TH1F*)input->Clone();\r
-  char name[100];\r
-  sprintf(name,"%s_ptCut",result->GetName());\r
-  result->SetNameTitle(name,name);\r
-  for(Int_t i = 1; i <= result->GetNbinsX(); i++) {\r
-    Double_t val = input->Integral(i,result->GetNbinsX());\r
-    result->SetBinContent(i,val);\r
-    result->SetBinError(i,0.);\r
-  }\r
-\r
-  return result;\r
-\r
-}\r
-\r