diff --git a/roofit/histfactory/inc/HistoToWorkspaceFactory.h b/roofit/histfactory/inc/HistoToWorkspaceFactory.h
index 6dac11d64b7aeac1acf2f7429e016aff8ca0c1b1..515cab92af409d91cd66ccc4fd74e3d81ec4c2e7 100644
--- a/roofit/histfactory/inc/HistoToWorkspaceFactory.h
+++ b/roofit/histfactory/inc/HistoToWorkspaceFactory.h
@@ -59,7 +59,7 @@ namespace HistFactory{
void Customize(RooWorkspace* proto, const char* pdfNameChar, map<string,string> renameMap);
- void EditSyst(RooWorkspace* proto, const char* pdfNameChar, map<string,double> gammaSyst, map<string,double> uniformSyst);
+ void EditSyst(RooWorkspace* proto, const char* pdfNameChar, map<string,double> gammaSyst, map<string,double> uniformSyst, map<string,double> logNormSyst);
void FormatFrameForLikelihood(RooPlot* frame, string XTitle=string("#sigma / #sigma_{SM}"), string YTitle=string("-log likelihood"));
diff --git a/roofit/histfactory/src/HistoToWorkspaceFactory.cxx b/roofit/histfactory/src/HistoToWorkspaceFactory.cxx
index 160d008fa056f33d20edd130104deb4220737235..1d59d437fefcf18b9966a2603d8b2b9a77e663a4 100644
--- a/roofit/histfactory/src/HistoToWorkspaceFactory.cxx
+++ b/roofit/histfactory/src/HistoToWorkspaceFactory.cxx
@@ -419,7 +419,7 @@ namespace HistFactory{
}
//_____________________________________________________________
- void HistoToWorkspaceFactory::EditSyst(RooWorkspace* proto, const char* pdfNameChar, map<string,double> gammaSyst, map<string,double> uniformSyst) {
+ void HistoToWorkspaceFactory::EditSyst(RooWorkspace* proto, const char* pdfNameChar, map<string,double> gammaSyst, map<string,double> uniformSyst,map<string,double> logNormSyst) {
// cout << "in edit, gammamap.size = " << gammaSyst.size() << ", unimap.size = " << uniformSyst.size() << endl;
string pdfName(pdfNameChar);
@@ -459,8 +459,7 @@ namespace HistFactory{
it->first.c_str())) ;
/*
- // this has some problems because N in poisson is rounded to nearest integer
-
+ // this has some problems because N in poisson is rounded to nearest integer
proto->factory(Form("Poisson::beta_%sConstraint(y_%s[%f],prod::taub_%s(taus_%s[%f],beta_%s[1,0,5]))",
it->first.c_str(),
it->first.c_str(),
@@ -475,13 +474,6 @@ namespace HistFactory{
// combined->factory(Form("expr::alphaOfBeta_%s('(beta_%s-1)/%f',beta_%s)",it->first.c_str(),it->first.c_str(),scale,it->first.c_str()));
proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1./scale,1./scale));
- // clean up constraints
- // cout << "got here, about to remove" << endl;
- // temp.remove(*proto->var(Form("alpha_%s",it->first.c_str())));
- // temp.add(*proto->var(Form("beta_%s",it->first.c_str())));
- // cout << "KC CHECK 2" << endl;
- // temp.Print();
-
// set beta const status to be same as alpha
if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
@@ -534,13 +526,73 @@ namespace HistFactory{
proto->factory(Form("Uniform::beta_%sConstraint(beta_%s)",it->first.c_str(),it->first.c_str()));
proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{-1,1})",it->first.c_str(),it->first.c_str()));
- // clean up constraints
- cout << "got here, about to remove" << endl;
- // temp.remove(*proto->var(Form("alpha_%s",it->first.c_str())));
- // temp.add(*proto->var(Form("beta_%s",it->first.c_str())));
- // cout << "KC CHECK 2" << endl;
- // temp.Print();
+ // set beta const status to be same as alpha
+ if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
+ proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
+ else
+ proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
+ // set alpha const status to true
+ // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
+
+ // replace alphas with alphaOfBeta and replace constraints
+ cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
+ editList+=preceed + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
+ preceed=",";
+ cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
+ editList+=preceed + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
+
+ if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
+ cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
+ else
+ cout << "NOT THERE" << endl;
+
+ // EDIT seems to die if the list of edits is too long. So chunck them up.
+ if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
+ edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
+ lastPdf+="_"; // append an underscore for the edit
+ editList=""; // reset edit list
+ preceed="";
+ cout << edit<< endl;
+ proto->factory( edit.c_str() );
+ RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
+ if(!newOne)
+ cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
+
+ }
+ }
+
+ /////////////////////////////////////////
+ ////////////////////////////////////
+
+ // add lognormal terms and their constraints
+ for(it=logNormSyst.begin(); it!=logNormSyst.end(); ++it) {
+ cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
+ if(! proto->var(("alpha_"+it->first).c_str())){
+ cout << "systematic not there" << endl;
+ nskipped++;
+ continue;
+ }
+ numReplacements++;
+
+ double relativeUncertainty = it->second;
+ double kappa = 1+relativeUncertainty;
+ // when transforming beta -> alpha, need alpha=1 to be +1sigma value.
+ // the P(beta>kappa*\hat(beta)) = 16%
+ // and \hat(beta) is 1, thus
+ double scale = relativeUncertainty;
+ //double scale = kappa;
+
+ // this is the LogNormal
+ proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
+ proto->factory(Form("kappa_%s[%f]",it->first.c_str(),kappa));
+ proto->factory(Form("Lognormal::beta_%sConstraint(beta_%s,one[1],kappa_%s)",
+ it->first.c_str(),
+ it->first.c_str(),
+ it->first.c_str())) ;
+ proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1./scale,1./scale));
+ // proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1.,1./scale));
+
// set beta const status to be same as alpha
if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
@@ -576,6 +628,9 @@ namespace HistFactory{
}
}
+ /////////////////////////////////////////
+ ////////////////////////////////////
+
// commit last bunch of edits
edit="EDIT::newSimPdf("+lastPdf+","+editList+")";
cout << edit<< endl;
diff --git a/roofit/histfactory/src/MakeModelAndMeasurements.cxx b/roofit/histfactory/src/MakeModelAndMeasurements.cxx
index 6c1da06b5497f87b399ee195b6fa65a0607e784b..79e530b3834b8658903e93af42799209e44ff38e 100644
--- a/roofit/histfactory/src/MakeModelAndMeasurements.cxx
+++ b/roofit/histfactory/src/MakeModelAndMeasurements.cxx
@@ -159,6 +159,7 @@ void topDriver(string input ){
vector<string> systToFix;
map<string,double> gammaSyst;
map<string,double> uniformSyst;
+ map<string,double> logNormSyst;
bool exportOnly = false;
// TListIter attribIt = node->GetAttributes();
@@ -233,6 +234,9 @@ void topDriver(string input ){
if (type=="Uniform" && rel!=0) {
for (vector<string>::const_iterator it=syst.begin(); it!=syst.end(); it++) uniformSyst[(*it).c_str()] = rel;
}
+ if (type=="LogNormal" && rel!=0) {
+ for (vector<string>::const_iterator it=syst.begin(); it!=syst.end(); it++) logNormSyst[(*it).c_str()] = rel;
+ }
}
mnode = mnode->GetNextNode();
}
@@ -288,9 +292,9 @@ void topDriver(string input ){
// Gamma/Uniform Constraints:
- // turn some Gaussian constraints into Gamma/Uniform constraints, rename model newSimPdf
- if(gammaSyst.size()>0 || uniformSyst.size()>0) {
- factory.EditSyst(ws,("model_"+oneChannel[0].channel).c_str(),gammaSyst,uniformSyst);
+ // turn some Gaussian constraints into Gamma/Uniform/LogNorm constraints, rename model newSimPdf
+ if(gammaSyst.size()>0 || uniformSyst.size()>0 || logNormSyst.size()>0) {
+ factory.EditSyst(ws,("model_"+oneChannel[0].channel).c_str(),gammaSyst,uniformSyst,logNormSyst);
proto_config->SetPdf(*ws->pdf("newSimPdf"));
}
@@ -317,9 +321,9 @@ void topDriver(string input ){
if(mode.find("comb")!=string::npos){
RooWorkspace* ws=factory.MakeCombinedModel(ch_names,chs);
// Gamma/Uniform Constraints:
- // turn some Gaussian constraints into Gamma/Uniform constraints, rename model newSimPdf
- if(gammaSyst.size()>0 || uniformSyst.size()>0)
- factory.EditSyst(ws,"simPdf",gammaSyst,uniformSyst);
+ // turn some Gaussian constraints into Gamma/Uniform/logNormal constraints, rename model newSimPdf
+ if(gammaSyst.size()>0 || uniformSyst.size()>0 || logNormSyst.size()>0)
+ factory.EditSyst(ws,"simPdf",gammaSyst,uniformSyst,logNormSyst);
//
// set parameter of interest according to the configuration
//
@@ -334,7 +338,7 @@ void topDriver(string input ){
ws->Print();
// Set new PDF if there are gamma/uniform constraint terms
- if(gammaSyst.size()>0 || uniformSyst.size()>0)
+ if(gammaSyst.size()>0 || uniformSyst.size()>0 || logNormSyst.size()>0)
combined_config->SetPdf(*ws->pdf("newSimPdf"));
RooAbsData* simData = ws->data("simData");
diff --git a/tutorials/histfactory/example.xml b/tutorials/histfactory/example.xml
index 600b0e796c22b1157bf18a6cda3e6ce8021309bd..28a5d2b810f0286135748e0f0d12d87c44a33cc7 100644
--- a/tutorials/histfactory/example.xml
+++ b/tutorials/histfactory/example.xml
@@ -18,10 +18,22 @@
<Input>./config/example_channel.xml</Input>
+ <Measurement Name="GaussExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
+ <POI>SigXsecOverSM</POI>
+ <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
+ <!-- don't need <ConstraintTerm> default is Gaussian-->
+ </Measurement>
+
<Measurement Name="GammaExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
<POI>SigXsecOverSM</POI>
<ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
- <ConstraintTerm Type="Gamma" RelativeUncertainty="1.4">syst2</ConstraintTerm>
+ <ConstraintTerm Type="Gamma" RelativeUncertainty=".3">syst2</ConstraintTerm>
+ </Measurement>
+
+ <Measurement Name="LogNormExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
+ <POI>SigXsecOverSM</POI>
+ <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
+ <ConstraintTerm Type="LogNormal" RelativeUncertainty=".3">syst2</ConstraintTerm>
</Measurement>
<Measurement Name="ConstExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" ExportOnly="True">
@@ -29,4 +41,5 @@
<ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
</Measurement>
+
</Combination>