diff --git a/data/sigetrans_v1.root b/data/sigetrans_v1.root new file mode 100644 index 0000000..6113f48 Binary files /dev/null and b/data/sigetrans_v1.root differ diff --git a/interface/SigETransform.h b/interface/SigETransform.h new file mode 100644 index 0000000..c3678ec --- /dev/null +++ b/interface/SigETransform.h @@ -0,0 +1,63 @@ +//-------------------------------------------------------------------------------------------------- +// $Id $ +// +// SigETransform +// +// Helper Class for applying transformation to SigmaE/E to decorrelated dependence on photon energy +// +// Authors: J.Bendavid +//-------------------------------------------------------------------------------------------------- + +#ifndef EGAMMATOOLS_SigETransform_H +#define EGAMMATOOLS_SigETransform_H + +#include "RooArgList.h" + +class RooWorkspace; +class RooRealVar; +class RooAbsPdf; +class RooAbsReal; +class HybridGBRForest; +class HybridGBRForestD; +class EcalClusterLazyTools; + +class SigETransform { + public: + SigETransform(); + ~SigETransform(); + + void Initialize(std::string regweights, int version); + Bool_t IsInitialized() const { return _isInitialized; } + + double sigEoverETranformed(double sigEoverE, double sceta, double energy); + + + protected: + + RooRealVar *_eerrvarEB; + RooRealVar *_eerrvarEE; + RooRealVar *_eerrvarinvEB; + RooRealVar *_eerrvarinvEE; + + RooRealVar *_scetavarEB; + RooRealVar *_scetavarEE; + RooRealVar *_scetavarinvEB; + RooRealVar *_scetavarinvEE; + + RooRealVar *_energyvarEB; + RooRealVar *_energyvarEE; + + RooAbsReal *_cdfEB; + RooAbsReal *_cdfEE; + RooAbsReal *_cdfinvEB; + RooAbsReal *_cdfinvEE; + + RooArgList _args; + + Bool_t _isInitialized; + + + }; + + +#endif diff --git a/src/SigETransform.cc b/src/SigETransform.cc new file mode 100644 index 0000000..d1f1efa --- /dev/null +++ b/src/SigETransform.cc @@ -0,0 +1,114 @@ +#include +#include "../interface/SigETransform.h" +#include "RooWorkspace.h" +#include "RooArgList.h" +#include "RooRealVar.h" +#include "RooAbsReal.h" +#include "RooAbsPdf.h" +#include "RooConstVar.h" +#include "TStreamerInfo.h" +#include "HiggsAnalysis/GBRLikelihood/interface/RooHybridBDTAutoPdf.h" +#include + +//-------------------------------------------------------------------------------------------------- +SigETransform::SigETransform() : +_eerrvarEB(0), +_eerrvarEE(0), +_eerrvarinvEB(0), +_eerrvarinvEE(0), +_scetavarEB(0), +_scetavarEE(0), +_scetavarinvEB(0), +_scetavarinvEE(0), +_energyvarEB(0), +_energyvarEE(0), +_cdfEB(0), +_cdfEE(0), +_cdfinvEB(0), +_cdfinvEE(0), +_isInitialized(kFALSE) +{ + // Constructor. +} + + +//-------------------------------------------------------------------------------------------------- +SigETransform::~SigETransform() +{ + +} + +//-------------------------------------------------------------------------------------------------- +void SigETransform::Initialize(std::string regweights, int version) { + + if (version==1) { + TFile *fin = TFile::Open(regweights.c_str(),"READ"); + RooWorkspace *win = static_cast(fin->Get("wssigetrans")); + + RooCondAddPdf *pdfEB = static_cast(win->pdf("eerrpdf_EB")); + RooCondAddPdf *pdfEE = static_cast(win->pdf("eerrpdf_EE")); + RooCondAddPdf *pdfinvEB = static_cast(win->pdf("eerrpdf_invEB")); + RooCondAddPdf *pdfinvEE = static_cast(win->pdf("eerrpdf_invEE")); + + _eerrvarEB = win->var("eerr_EB"); + _eerrvarEE = win->var("eerr_EE"); + _eerrvarinvEB = win->var("eerr_invEB"); + _eerrvarinvEE = win->var("eerr_invEE"); + + _scetavarEB = win->var("sceta_EB"); + _scetavarEE = win->var("sceta_EE"); + _scetavarinvEB = win->var("sceta_invEB"); + _scetavarinvEE = win->var("sceta_invEE"); + + _energyvarEB = win->var("energy_EB"); + _energyvarEE = win->var("energy_EE"); + + _cdfEB = pdfEB->createCDF(*_eerrvarEB); + _cdfEE = pdfEE->createCDF(*_eerrvarEE); + _cdfinvEB = pdfinvEB->createCDF(*_eerrvarinvEB); + _cdfinvEE = pdfinvEE->createCDF(*_eerrvarinvEE); + + _args.addOwned(*_cdfEB); + _args.addOwned(*_cdfEE); + _args.addOwned(*_cdfinvEB); + _args.addOwned(*_cdfinvEE); + + _isInitialized = true; + } + +} + + +//-------------------------------------------------------------------------------------------------- +double SigETransform::sigEoverETranformed(double sigEoverE, double sceta, double energy) { + + bool isbarrel = std::abs(sceta)<1.479; + + if (isbarrel) { + _eerrvarEB->setVal(sigEoverE); + _eerrvarinvEB->setVal(sigEoverE); + + _scetavarEB->setVal(sceta); + _scetavarinvEB->setVal(sceta); + + _energyvarEB->setVal(energy); + + double cdfval = _cdfEB->getVal(); + double sigEoverEtransformed = _cdfinvEB->findRoot(*_eerrvarinvEB,0.,1.,cdfval); + return sigEoverEtransformed; + } + else { + _eerrvarEE->setVal(sigEoverE); + _eerrvarinvEE->setVal(sigEoverE); + + _scetavarEE->setVal(sceta); + _scetavarinvEE->setVal(sceta); + + _energyvarEE->setVal(energy); + + double cdfval = _cdfEE->getVal(); + double sigEoverEtransformed = _cdfinvEE->findRoot(*_eerrvarinvEE,0.,1.,cdfval); + return sigEoverEtransformed; + } + +} \ No newline at end of file diff --git a/test/sigtrans.C b/test/sigtrans.C new file mode 100644 index 0000000..f56ca9c --- /dev/null +++ b/test/sigtrans.C @@ -0,0 +1,121 @@ +#include "RooRealVar.h" +#include "RooAbsPdf.h" +#include "RooExponential.h" +#include "RooGaussian.h" +#include "RooPlot.h" +#include "TCanvas.h" +#include "RooConstVar.h" +#include "RooDataSet.h" +#include "RooHybridBDTAutoPdf.h" +#include "RooFormulaVar.h" +#include "RooProdPdf.h" +#include "RooUniform.h" +#include "TRandom.h" +#include "TGraph.h" +#include "RooAddPdf.h" +#include "RooNDKeysPdf.h" +#include "RooExtendPdf.h" +#include "RooMinimizer.h" +#include "TFile.h" +#include "TNtuple.h" +#include "HybridGBRForest.h" +#include "RooProduct.h" +#include "RooChebychev.h" +#include "RooBernstein.h" +#include "RooPolynomial.h" +#include "RooGenericPdf.h" +//#include "HZZ2L2QRooPdfs.h" +#include "RooArgSet.h" +#include "RooArgList.h" +#include "RooCBShape.h" +//#include "RooCBShapeModified.h" +//#include "RooBernsteinFast.h" +#include "RooWorkspace.h" +#include "TH1D.h" +#include "TChain.h" +#include "TCut.h" +#include "TLine.h" +#include "TLegend.h" +#include "RooRandom.h" +#include "RooAddition.h" +#include "TSystem.h" +#include "RooCBFast.h" +#include "RooDoubleCBFast.h" +#include "RooGaussianFast.h" +//#include "RooBernsteinFast.h" +#include "TProfile.h" +#include "RooRevCBFast.h" +#include "RooBifurGauss.h" +#include "SigETransform.h" + +SigETransform *g_sigtrans; +double sigetransformed(double sige,double sceta,double energy) { + + return g_sigtrans->sigEoverETranformed(sige,sceta,energy); + +} + +void sigtrans() { + +TString dirname = "/scratch/bendavid/root/bare/sigtransplotsMay7/"; +gSystem->mkdir(dirname,true); +gSystem->cd(dirname); + + g_sigtrans = new SigETransform; + g_sigtrans->Initialize("/home/bendavid/CMSSW_6_1_2/src/HiggsAnalysis/GBRLikelihoodEGTools/data/sigetrans_v1.root",1); + + TChain *tree = new TChain("RunLumiSelectionMod/MCProcessSelectionMod/HLTModP/GoodPVFilterMod/PhotonMvaMod/JetPub/JetCorrectionMod/SeparatePileUpMod/ElectronIDMod/MuonIDMod/PhotonPairSelectorPresel/PhotonTreeWriterPresel/hPhotonTree"); + tree->Add("/home/mingyang/cms/hist/hgg-2013Final8TeV_train/merged/hgg-2013Final8TeV_s12-diphoj-m60-v7n_noskim.root"); + + //tree->Print(); + + TCut selcut = "(ph1.pt > (mass/3.0) && ph2.pt > (mass/4.0) && mass>100. && mass<180. && ph1.idmva>-0.2 && ph2.idmva>-0.2)"; + TCut selweight = "mcweight"; + TCut prescale100 = "(evt%10==0)"; + TCut ebcentral = "abs(ph1.sceta)<1. && abs(ph2.sceta)<1."; + + tree->SetAlias("ph1.sigeoetrans","sigetransformed(ph1.eerr/ph1.e,ph1.sceta,ph1.e)"); + tree->SetAlias("ph2.sigeoetrans","sigetransformed(ph2.eerr/ph2.e,ph2.sceta,ph2.e)"); + tree->SetAlias("sigmom","0.5*sqrt( (ph1.eerr/ph1.e)^2 + (ph1.esmearing/ph1.e)^2 + (ph2.eerr/ph2.e)^2 + (ph2.esmearing/ph2.e)^2)"); + tree->SetAlias("sigmomtrans","0.5*sqrt( (ph1.sigeoetrans)^2 + (ph1.esmearing/ph1.e)^2 + (ph2.sigeoetrans)^2 + (ph2.esmearing/ph2.e)^2)"); + + //tree->Draw("sigmom/(masserrsmeared/mass)",selcut*selweight*prescale100); + + TCanvas *c0 = new TCanvas; + tree->Draw("sigmom:sigmomtrans>>htmp0(100,0.,0.03,100,0.,0.03)",selcut*selweight*prescale100,"COLZ"); + c0->SaveAs("sigmvtrans.eps"); + + TCanvas *c1 = new TCanvas; + tree->Draw("sigmom:mass>>htmp1(100,100.,180.,100,0.,0.03)",selcut*selweight*prescale100,"COLZ"); + c1->SaveAs("sigmvmass.eps"); + + TCanvas *c2 = new TCanvas; + tree->Draw("sigmomtrans:mass>>htmp2(100,100.,180.,100,0.,0.03)",selcut*selweight*prescale100,"COLZ"); + c2->SaveAs("sigmtransvmass.eps"); + + TCanvas *c3 = new TCanvas; + tree->Draw("sigmom:TMath::Max(abs(ph1.sceta),abs(ph2.sceta))>>htmp3(100,0.,2.5,100,0.,0.03)",selcut*selweight*prescale100,"COLZ"); + c3->SaveAs("sigmvmaxeta.eps"); + + TCanvas *c4 = new TCanvas; + tree->Draw("sigmomtrans:TMath::Max(abs(ph1.sceta),abs(ph2.sceta))>>htmp4(100,0.,2.5,100,0.,0.03)",selcut*selweight*prescale100,"COLZ"); + c4->SaveAs("sigmtransvmaxeta.eps"); + + TCanvas *c5 = new TCanvas; + tree->Draw("sigmom:mass>>htmp5(100,100.,180.,100,0.,0.03)",selcut*selweight*prescale100*ebcentral,"COLZ"); + c5->SaveAs("sigmvmassEBcentral.eps"); + + TCanvas *c6 = new TCanvas; + tree->Draw("sigmomtrans:mass>>htmp6(100,100.,180.,100,0.,0.03)",selcut*selweight*prescale100*ebcentral,"COLZ"); + c6->SaveAs("sigmtransvmassEBcentral.eps"); + + //tree->Draw("sigmom>>(100,0.,0.07)",selcut*selweight*prescale100,"HIST"); + + + + //new TCanvas; + //tree->Draw("sigmomtrans",selcut*selweight*prescale100,"ESAME"); + + return; + +} \ No newline at end of file