Skip to content

Instantly share code, notes, and snippets.

@ram1123
Last active July 4, 2018 11:11
Show Gist options
  • Save ram1123/f6fce437307cce7b3fba2b25a3c9cc02 to your computer and use it in GitHub Desktop.
Save ram1123/f6fce437307cce7b3fba2b25a3c9cc02 to your computer and use it in GitHub Desktop.
Issue: Compare Integral of histogram and equivalent function

(RooFit) How to get integral of a RooRealVal

import ROOT as r

r.gSystem.Load("PDFs/PdfDiagonalizer_cc.so")
r.gSystem.Load("PDFs/Util_cxx.so")
r.gSystem.Load("PDFs/HWWLVJRooPdfs_cxx.so")

from ROOT import RooErfExpPdf

msgservice = r.RooMsgService.instance()
msgservice.setGlobalKillBelow(r.RooFit.FATAL)

fin = r.TFile.Open("root://cmseos.fnal.gov//store/user/rasharma/SecondStep/WWTree_2018_01_07_12h02/HaddedFiles/Hadds_for_BkgEstimation/WWTree_WJets.root")

treeIn = fin.Get("otree")

c1 = r.TCanvas("c1","c1",800,600)

h1 = r.TH1F("h1","",22,40,150)
h1.Sumw2()
h1.SetMinimum(0.)
h1.SetMaximum(850.)

#cutString = #"wSampleWeight*35867.06*genWeight*" \
	    # "pu_Weight*trig_eff_Weight*id_eff_Weight*btag0Wgt*" \
cutString =  "(l_pt2<0 && " \
	     "l_pt1>30 && abs(l_eta1)<2.5 && " \
	     "pfMET_Corr>50 && " \
	     "PuppiAK8_jet_tau2tau1<0.55 && " \
	     "ungroomed_PuppiAK8_jet_pt>200 && abs(ungroomed_PuppiAK8_jet_eta)<2.4 && " \
	     "mass_lvj_type0>600 && mass_lvj_type0<4000 && " \
	     "nBTagJet_loose==0 && " \
	     "BosonCentrality_type0>1.0 && " \
	     "ZeppenfeldWL_type0/vbf_maxpt_jj_Deta<0.3 && " \
	     "ZeppenfeldWH/vbf_maxpt_jj_Deta <0.3 && " \
	     "vbf_maxpt_jj_m>800 && " \
	     "vbf_maxpt_jj_Deta>4.0 && " \
	     "vbf_maxpt_j1_pt>30 && vbf_maxpt_j2_pt>30) && " \
	     "PuppiAK8_jet_mass_so>40 && PuppiAK8_jet_mass_so<150"


treeIn.Draw("PuppiAK8_jet_mass_so>>h1",cutString)

rrv_x		= r.RooRealVar("rrv_x","AK8 invariant mass",40,150);
rrv_c_ErfExp      = r.RooRealVar("rrv_c_ErfExp","rrv_c_ErfExp",-0.05,-0.1,-1e-4);
rrv_offset_ErfExp = r.RooRealVar("rrv_offset_ErfExp","rrv_offset_ErfExp",60.,30.,120);
rrv_width_ErfExp  = r.RooRealVar("rrv_width_ErfExp","rrv_width_ErfExp",30.,10, 60.);
model_pdf         = r.RooErfExpPdf("model_pdf","model_pdf",rrv_x,rrv_c_ErfExp,rrv_offset_ErfExp,rrv_width_ErfExp);


print "Entries from MC = ",h1.GetEntries(),h1.Integral()

dh =  r.RooDataHist("dh","plotOn test data with x",r.RooArgList(rrv_x),h1)

h1.Delete()
c1.Clear()


model_pdf.fitTo(dh)

rfresult =  model_pdf.createIntegral(r.RooArgSet(rrv_x))

print "*"*30
print "*"
print rfresult.getVal() 
print "*"
print "*"*30

frame = rrv_x.frame() 
dh.plotOn(frame)
dh.statOn(frame);  #//this will display hist stat on canvas
model_pdf.plotOn(frame)
#model_pdf.plotOn(frame,r.MarkerColor(2),r.MarkerSize(0.9),r.MarkerStyle(21))
model_pdf.paramOn(frame);
frame.Draw()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment