add tool for transformation to energy-decorrelated sigmaE/E variable

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

src/SigETransform.cc

@@ -0,0 +1,114 @@
+#include <TFile.h>
+#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 <cmath>
+
+//--------------------------------------------------------------------------------------------------
+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<RooWorkspace*>(fin->Get("wssigetrans"));
+
+      RooCondAddPdf *pdfEB = static_cast<RooCondAddPdf*>(win->pdf("eerrpdf_EB"));
+      RooCondAddPdf *pdfEE = static_cast<RooCondAddPdf*>(win->pdf("eerrpdf_EE"));
+      RooCondAddPdf *pdfinvEB = static_cast<RooCondAddPdf*>(win->pdf("eerrpdf_invEB"));
+      RooCondAddPdf *pdfinvEE = static_cast<RooCondAddPdf*>(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;
+  }
+
+}

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; 
+
+}