Skip to content

Commit

Permalink
Merge pull request #9 from musella/topic_diff_xs_sigmae
Browse files Browse the repository at this point in the history
Tweak sigmaE/E variations
  • Loading branch information
musella committed Aug 19, 2014
2 parents a011c98 + eb43339 commit 2eeac1a
Show file tree
Hide file tree
Showing 13 changed files with 138 additions and 55 deletions.
2 changes: 2 additions & 0 deletions AnalysisScripts/common/photonanalysis.dat
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ puTarget=aux/HCP2012_freezing_v2.json.69400.observed.pileup.root
energyCorrectionMethod=Bendavid
doEcorrectionSmear=1
doEcorrectionSyst=0
sigmaEoEUncEB=0.10
sigmaEoEUncEE=0.10

## energy scale adjustement and errors
scale_offset_file=energy_scale_offsets.dat
Expand Down
1 change: 1 addition & 0 deletions AnalysisScripts/common/systs_settings.dat
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
doSystematics=1

# do 1-sigma shifts and book only mH=125
systRange=1
Expand Down
5 changes: 5 additions & 0 deletions AnalysisScripts/diffanalysis/differentialnalysis.dat
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,11 @@ phoidMvaCutEE=0.10
# 3 categories decorrelated
sigmaMoMCategoryBoundaries=0.0183,0.0128,0.0079,0

doSystematics=0
sigmaEoEUncEB=0.05
sigmaEoEUncEE=0.10
doPhotonMvaIdSyst=0

massMin=100.
massMax=180.
nDataBins=160
Expand Down
1 change: 1 addition & 0 deletions AnalysisScripts/diffanalysis/unfolding.dat
Original file line number Diff line number Diff line change
Expand Up @@ -15,3 +15,4 @@ JetEtaForDiffAnalysis=2.5

doUnfoldHisto=1
doProcessSplitting=0

4 changes: 2 additions & 2 deletions EnergySmearer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -190,9 +190,9 @@ bool EnergySmearer::smearPhoton(PhotonReducedInfo & aPho, float & weight, int ru
// leave energy alone, bus change resolution (10% uncertainty on sigmaE/E scaling)
float newSigma;
if (fabs(aPho.caloPosition().Eta())<1.5){
newSigma = aPho.rawCorrEnergyErr()*(1.+syst_shift*0.1);
newSigma = aPho.rawCorrEnergyErr()*(1.+syst_shift*myParameters_.sigmaEoEUncEB);
} else {
newSigma = aPho.rawCorrEnergyErr()*(1.+syst_shift*0.1);
newSigma = aPho.rawCorrEnergyErr()*(1.+syst_shift*myParameters_.sigmaEoEUncEE);
}
aPho.setCorrEnergyErr(newSigma);
} else {
Expand Down
2 changes: 2 additions & 0 deletions EnergySmearer.h
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,8 @@ class EnergySmearer : public BaseSmearer
std::string efficiency_file;
// errors on correction will be a fraction of the correction itself
float corrRelErr;
float sigmaEoEUncEB;
float sigmaEoEUncEE;
};

EnergySmearer(EnergySmearer * orig, const std::vector<PhotonCategory> & presel=std::vector<PhotonCategory>());
Expand Down
2 changes: 1 addition & 1 deletion Macros/FinalResults/makeCombinePlots.py
Original file line number Diff line number Diff line change
Expand Up @@ -736,7 +736,7 @@ def plot1DNLL(returnErrors=False,xvar="", ext=""):
eplus2 = h2-m
eminus2 = m-l2

print "%15s : %4.3f +%4.3g -%4.3g" % ( ntitle+" "+ext, xmin, eplus , eminus )
print "%15s : %4.4f +%4.4g -%4.4g" % ( ntitle+" "+ext, xmin, eplus , eminus )

if returnErrors:
return [xmin,eplus,eminus,eplus2,eminus2]
Expand Down
4 changes: 3 additions & 1 deletion PhotonAnalysis/interface/PhotonAnalysis.h
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,9 @@ class PhotonAnalysis : public BaseAnalysis

string sigEoEtransformFile;
bool doSigEoEtransform;

float sigmaEoEUncEB;
float sigmaEoEUncEE;

int nEtaCategories, nR9Categories, nPtCategories, nPtOverMCategories, nVtxCategories;
float R9CatBoundary;
bool usePUjetveto;
Expand Down
4 changes: 4 additions & 0 deletions PhotonAnalysis/src/PhotonAnalysis.cc
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,8 @@ PhotonAnalysis::PhotonAnalysis() :

sigEoEtransformFile="";
doSigEoEtransform = false;
sigmaEoEUncEB = 0.1;
sigmaEoEUncEE = 0.1;
}

// ----------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -1285,6 +1287,8 @@ void PhotonAnalysis::Init(LoopAll& l)
eSmearPars.byRun = false;
eSmearPars.n_categories = tmp_smear_smearing_cat.size();
eSmearPars.photon_categories = tmp_smear_smearing_cat;
eSmearPars.sigmaEoEUncEB = sigmaEoEUncEB;
eSmearPars.sigmaEoEUncEE = sigmaEoEUncEE;

eSmearPars.smearing_sigma = tmp_smear_smearing[0].scale_offset;
eSmearPars.smearing_sigma_error = tmp_smear_smearing[0].scale_offset_error;
Expand Down
6 changes: 4 additions & 2 deletions SimultaneousSignalFitting/interface/FinalModelConstruction.h
Original file line number Diff line number Diff line change
Expand Up @@ -146,13 +146,15 @@ class FinalModelConstruction {
std::map<std::string,std::vector<std::pair<int,float> > > globalScalesOpts;
std::map<std::string,std::vector<std::pair<int,float> > > globalScalesCorrOpts;
std::vector<std::string> systematicsList;

std::vector<float> systematicsCorr;
std::vector<int> systematicsIdx;

std::vector<std::string> photonCats;
std::map<std::string,std::map<int,std::map<std::string,double> > > meanSysts;
std::map<std::string,std::map<int,std::map<std::string,double> > > sigmaSysts;
std::map<std::string,std::map<int,std::map<std::string,double> > > rateSysts;

std::map<string,RooRealVar*> photonSystematics;
std::map<string,RooAbsReal*> photonSystematics;
std::map<string,RooConstVar*> photonSystematicConsts;

// utility funcs
Expand Down
5 changes: 3 additions & 2 deletions SimultaneousSignalFitting/scripts/checkNuisances.py
Original file line number Diff line number Diff line change
Expand Up @@ -158,15 +158,16 @@ def main(options,args):
v = ws.var(isyst)
if v:
roonuis[ isyst ] = v

print roonuis

perproc = None
if dat != "":
perproc = readDatFile(dat)
pprint( perproc )

for icat in range(ncat):
pdf = getPdf(icat, sqrts, ws)
pdf.Print()
pdf.Print("V")
vars = {}
expect = {}
diff = {}
Expand Down
100 changes: 81 additions & 19 deletions SimultaneousSignalFitting/src/FinalModelConstruction.cc
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#include <fstream>
#include <sstream>

#include "TVectorT.h"
#include "TMatrixTSym.h"
#include "TCanvas.h"
#include "TF1.h"
#include "RooPlot.h"
Expand Down Expand Up @@ -101,15 +103,28 @@ FinalModelConstruction::~FinalModelConstruction(){}

void FinalModelConstruction::addToSystematicsList(vector<string>::iterator begin, vector<string>::iterator end){
for (vector<string>::iterator it=begin; it!=end; it++){
if (find(systematicsList.begin(),systematicsList.end(),*it)!=systematicsList.end()) {
string name = *it;
float corr = 0.;
int index = -1;
if( it->find(":") != string::npos ) {
vector<string> toks;
split(toks,*it,boost::is_any_of(":"));
name = toks[0];
corr = boost::lexical_cast<float>(toks[1]);
index = boost::lexical_cast<int>(toks[2]);
/// *it = Form("%dTeV%s",sqrts_,name.c_str());
*it = name;
}
if (find(systematicsList.begin(),systematicsList.end(),name)!=systematicsList.end()) {
cout << "ERROR - duplicate systematic names! " << *it << " already found in systematics list." << endl;
exit(1);
}
else {
systematicsList.push_back(*it);
systematicsList.push_back(name);
systematicsCorr.push_back(corr);
systematicsIdx.push_back(index);
}
}

}

void FinalModelConstruction::addToSystematicsList(vector<string> systs){
Expand Down Expand Up @@ -209,7 +224,7 @@ void FinalModelConstruction::printSignalSystematics(){
}
// nuisance parameters
cout << "Implementing the following floating nuisance parameters" << endl;
for (map<string,RooRealVar*>::iterator sys=photonSystematics.begin(); sys!=photonSystematics.end(); sys++){
for (map<string,RooAbsReal*>::iterator sys=photonSystematics.begin(); sys!=photonSystematics.end(); sys++){
cout << "\t " << Form("%-50s",sys->first.c_str()) << " -- "; sys->second->Print();
}
// const parameters
Expand Down Expand Up @@ -241,42 +256,42 @@ void FinalModelConstruction::loadSignalSystematics(string filename){
if (starts_with(line,"photonCatScales=")){
line = line.substr(line.find("=")+1,string::npos);
if (line.empty()) continue;
split_append(photonCatScales,line,boost::is_any_of(","));
vector<string>::iterator beg = split_append(photonCatScales,line,boost::is_any_of(","));
addToSystematicsList(beg,photonCatScales.end());
if (verbosity_){
cout << "PhotonCatScales: ";
if (verbosity_) printVec(photonCatScales);
}
addToSystematicsList(photonCatScales);
}
else if (starts_with(line,"photonCatScalesCorr=")){
line = line.substr(line.find("=")+1,string::npos);
if (line.empty()) continue;
vector<string>::iterator beg = split_append(photonCatScalesCorr,line,boost::is_any_of(","));
addToSystematicsList(beg,photonCatScalesCorr.end());
if (verbosity_){
cout << "PhotonCatScalesCorr: ";
if (verbosity_) printVec(photonCatScalesCorr);
}
addToSystematicsList(beg,photonCatScalesCorr.end());
}
else if (starts_with(line,"photonCatSmears=")){
line = line.substr(line.find("=")+1,string::npos);
if (line.empty()) continue;
split_append(photonCatSmears,line,boost::is_any_of(","));
vector<string>::iterator beg = split_append(photonCatSmears,line,boost::is_any_of(","));
addToSystematicsList(beg,photonCatSmears.end());
if (verbosity_){
cout << "PhotonCatSmears: ";
if (verbosity_) printVec(photonCatSmears);
}
addToSystematicsList(photonCatSmears);
}
else if (starts_with(line,"photonCatSmearsCorr=")){
line = line.substr(line.find("=")+1,string::npos);
if (line.empty()) continue;
split_append(photonCatSmearsCorr,line,boost::is_any_of(","));
vector<string>::iterator beg = split_append(photonCatSmearsCorr,line,boost::is_any_of(","));
addToSystematicsList(beg,photonCatSmearsCorr.end());
if (verbosity_){
cout << "PhotonCatSmearsCorr: ";
if (verbosity_) printVec(photonCatSmearsCorr);
}
addToSystematicsList(photonCatSmearsCorr);
}
else if (starts_with(line,"globalScales=")){
line = line.substr(line.find("=")+1,string::npos);
Expand Down Expand Up @@ -373,10 +388,53 @@ void FinalModelConstruction::loadSignalSystematics(string filename){
datfile.close();

// now make the actual systematics
for (vector<string>::iterator it=systematicsList.begin(); it!=systematicsList.end(); it++){
RooRealVar *var = new RooRealVar(Form("CMS_hgg_nuisance_%s",it->c_str()),Form("CMS_hgg_nuisance_%s",it->c_str()),0.,-5.,5.);
var->setConstant(true);
photonSystematics.insert(make_pair(var->GetName(),var));
// for (vector<string>::iterator it=systematicsList.begin(); it!=systematicsList.end(); it++){
for (size_t is=0; is<systematicsList.size(); ++is) {
std::string & name = systematicsList[is];
float & corr = systematicsCorr[is];
int & idx = systematicsIdx[is];
if( corr != 0. ) {
TString nuname = name;
nuname = nuname.ReplaceAll(Form("%dTeV",sqrts_),"");
RooRealVar *var1 = new RooRealVar(Form("CMS_hgg_nuisance_eigen1_%s",nuname.Data()),Form("CMS_hgg_nuisance_eigen1_%s",nuname.Data()),0.,-5.,5.);
var1->setConstant(true);
photonSystematics.insert(make_pair(var1->GetName(),var1));
RooRealVar *var2 = new RooRealVar(Form("CMS_hgg_nuisance_eigen2_%s",nuname.Data()),Form("CMS_hgg_nuisance_eigen2_%s",nuname.Data()),0.,-5.,5.);
var2->setConstant(true);
photonSystematics.insert(make_pair(var2->GetName(),var2));

TMatrixTSym<double> varmat(2);
varmat[0][0] = 1.;//sigmaX*sigmaX;
varmat[0][1] = corr;//*sigmaX*sigmaY;
varmat[1][0] = corr;//*sigmaX*sigmaY;
varmat[1][1] = 1.;//sigmaY*sigmaY;

//// // Print Matrix
//// std::cout << "Covariance Matrix" << std::endl;
//// varmat.Print();
//// std::cout << "-----------------" << std::endl;
////
//// // 2 Get the eigenvalues/eigenvectors
TVectorT<double> eval;
TMatrixT<double> evec = varmat.EigenVectors(eval);
//// std::cout << "Matrix of Eigenvectors" << std::endl;
//// evec.Print();
//// std::cout << "----------------------" << std::endl;
double cs = evec[idx][0]*sqrt(eval[0]) , ss = evec[idx][1]*sqrt(eval[1]);
//// double tot = sqrt( cs*cs + ss*ss );
//// std::cout << "size of eigenvec " << tot << std::endl;
/// cs /= tot; ss /= tot;

RooFormulaVar *var = new RooFormulaVar(Form("CMS_hgg_nuisance_%s",name.c_str()),
Form("CMS_hgg_nuisance_%s",name.c_str()),
Form("%1.3g * @0 + %1.3g * @1",cs,ss), RooArgList(*var1,*var2) );
photonSystematics.insert(make_pair(var->GetName(),var));
// var->Print("V");
} else {
RooRealVar *var = new RooRealVar(Form("CMS_hgg_nuisance_%s",name.c_str()),Form("CMS_hgg_nuisance_%s",name.c_str()),0.,-5.,5.);
var->setConstant(true);
photonSystematics.insert(make_pair(var->GetName(),var));
}
}
}

Expand Down Expand Up @@ -520,7 +578,7 @@ RooAbsReal* FinalModelConstruction::getMeanWithPhotonSyst(RooAbsReal *dm, string
string syst = systematicsList[i];
int formPlace = dependents->getSize();
if (isGlobalSyst(syst)) {
RooRealVar *nuisScale = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
RooAbsReal *nuisScale = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
formula += Form("+@%d",formPlace);
// should check special extras
float additionalFactor = getRequiredAddtionalGlobalScaleFactor(syst);
Expand All @@ -536,7 +594,11 @@ RooAbsReal* FinalModelConstruction::getMeanWithPhotonSyst(RooAbsReal *dm, string
if (isPerCatSyst(syst)) {
if (photonSystematicConsts.find(Form("const_%s_cat%d_%dTeV_mean_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())) != photonSystematicConsts.end() ) {
RooConstVar *constVar = photonSystematicConsts[Form("const_%s_cat%d_%dTeV_mean_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())];
RooRealVar *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
RooAbsReal *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
if( verbosity_ ) {
std::cout << "Systematic " << syst << std::endl;
nuisVar->Print("V");
}
if ( fabs(constVar->getVal())>=5.e-5) {
hasEffect = true;
formula += Form("+@%d*@%d",formPlace,formPlace+1);
Expand Down Expand Up @@ -569,7 +631,7 @@ RooAbsReal* FinalModelConstruction::getSigmaWithPhotonSyst(RooAbsReal *sig_fit,
if (isPerCatSyst(syst)) {
if (photonSystematicConsts.find(Form("const_%s_cat%d_%dTeV_sigma_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())) != photonSystematicConsts.end() ) {
RooConstVar *constVar = photonSystematicConsts[Form("const_%s_cat%d_%dTeV_sigma_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())];
RooRealVar *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
RooAbsReal *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
if (constVar->getVal()>=1.e-4) {
hasEffect = true;
if( quadraticSigmaSum_ ) {
Expand Down Expand Up @@ -603,7 +665,7 @@ RooAbsReal* FinalModelConstruction::getRateWithPhotonSyst(string name){
if (isPerCatSyst(syst)) {
if (photonSystematicConsts.find(Form("const_%s_cat%d_%dTeV_rate_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())) != photonSystematicConsts.end() ) {
RooConstVar *constVar = photonSystematicConsts[Form("const_%s_cat%d_%dTeV_rate_%s",proc_.c_str(),cat_,sqrts_,syst.c_str())];
RooRealVar *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
RooAbsReal *nuisVar = photonSystematics[Form("CMS_hgg_nuisance_%s",syst.c_str())];
if (constVar->getVal()>=5.e-4) {
hasEffect = true;
formula += Form("+@%d*@%d",formPlace,formPlace+1);
Expand Down
Loading

0 comments on commit 2eeac1a

Please sign in to comment.