Merge pull request cms-sw#126 from nucleosynthesis/nckw_add_groups_of…

…_nuisances

Add groups of nuisances

data/tutorials/groupnuisA.txt

@@ -0,0 +1,46 @@
+# -S 0 
+# --fastScan (MultiDimFit)
+# --profilingMode (ChannelCompatibilityCheck)
+
+
+imax 1  # number of channels
+jmax 1  # number of processes-1
+kmax *  # number of nuisance parameters
+--------------------------
+bin             hww       
+observation     18        
+-------------------------- 
+bin             hww hww   
+process         ggH bkgww
+process           0   1   
+rate             11   7   
+
+lumi lnN        1.05  -
+theo_ggH lnN	1.10  -
+theo_bkgww lnN	- 1.10
+exp_systww lnN		    1.05 1.05
+MCstatww lnN       1.01 -
+
+theo group = theo_bkgww theo_ggH
+theoB group = theo_bkgww
+theoS group = theo_ggH
+lumi group = lumi
+stat group = MCstatww
+
+#does not exist: error in modeltools
+#dummy group = foobar
+
+#empty error
+#testempty group
+
+#missing [+]=
+#testsyntax group bla bla bla
+#testsyntax group -= bla bla bla
+
+#add before define
+#add group += lumi
+
+equals group = lumi
+equals group += lumi MCstatww
+#redefinition error
+#equals group = theo_ggH

data/tutorials/groupnuisB.txt

@@ -0,0 +1,24 @@
+imax 1  # number of channels
+jmax 1  # number of processes-1
+kmax *  # number of nuisance parameters
+--------------------------
+bin             hgg
+observation     25
+-------------------------- 
+bin             hgg hgg
+process         ggH bkggg
+process           0   1
+rate             10   8
+
+lumi lnN        1.05  -
+theo_ggH lnN	1.10  -
+theo_bkggg lnN	- 1.10
+exp_systgg lnN		    1.05 1.05
+MCstatgg lnN       1.01 -
+
+theo group = theo_bkggg theo_ggH
+theoB group = theo_bkggg
+theoS group = theo_ggH
+lumi group = lumi
+stat group = MCstatgg
+

interface/Combine.h

@@ -3,6 +3,8 @@
 #include <TString.h>
 #include <boost/program_options.hpp>
 #include "RooArgSet.h"
+#include <boost/algorithm/string/split.hpp>
+#include <boost/algorithm/string/classification.hpp>
 
 class TDirectory;
 class TTree;
@@ -62,6 +64,7 @@ class Combine {
   std::string setPhysicsModelParameterRangeExpression_;
   std::string redefineSignalPOIs_;
   std::string freezeNuisances_;
+  std::string freezeNuisanceGroups_;
 
   // input-output related variables
   bool saveWorkspace_;

python/DatacardParser.py

@@ -55,6 +55,7 @@ def parseCard(file, options):
         raise RuntimeError, "You should pass as argument to parseCards a file object, stream or a list of lines, not a string"
     ret = Datacard()
     ret.discretes=[]
+    ret.groups={}
     #
     nbins      = -1; 
     nprocesses = -1; 
@@ -178,6 +179,31 @@ def parseCard(file, options):
             elif pdf=="discrete":
                 args = f[2:]
                 ret.discretes.append(lsyst)
+                continue
+            elif pdf=="group":
+                # This is not really a pdf type, but a way to be able to name groups of nuisances together
+                groupName = lsyst
+                groupNuisances = numbers
+
+                if not groupNuisances:
+                    raise RuntimeError, "Syntax error for group '%s': empty line after 'group'." % groupName
+
+                defToks = ('=','+=')
+                defTok = groupNuisances.pop(0)
+                if defTok not in defToks:
+                    raise RuntimeError, "Syntax error for group '%s': first thing after 'group' is not '[+]=' but '%s'." % (groupName,defTok)
+
+                if groupName not in ret.groups:
+                    if defTok=='=':
+                        ret.groups[groupName] = set(groupNuisances)
+                    else:
+                        raise RuntimeError, "Cannot append to group '%s' as it was not yet defined." % groupName                                                                                                    
+                else:
+                    if defTok=='+=' :
+                        ret.groups[groupName].update( set(groupNuisances) )
+                    else:
+                        raise RuntimeError, "Will not redefine group '%s'. It previously contained '%s' and you now wanted it to contain '%s'." % (groupName,ret.groups[groupName],groupNuisances)                        
+
                 continue
             else:
                 raise RuntimeError, "Unsupported pdf %s" % pdf

python/ModelTools.py

@@ -70,6 +70,14 @@ def setPhysics(self,physicsModel):
     def doModel(self):
         self.doObservables()
         self.physics.doParametersOfInterest()
+
+        # set a group attribute on POI variables
+        poiIter = self.out.set('POI').createIterator()
+        poi = poiIter.Next()
+        while poi:
+            self.out.var(poi.GetName()).setAttribute('group_POI',True)
+            poi = poiIter.Next()
+
         self.physics.preProcessNuisances(self.DC.systs)
         self.doNuisances()
         self.doExpectedEvents()
@@ -78,18 +86,32 @@ def doModel(self):
         self.physics.done()
         if self.options.bin:
             self.doModelConfigs()
-            if self.options.verbose > 2: self.out.Print("V")
-            if self.options.verbose > 3: 
-                print "Wrote GraphVizTree of model_s to ",self.options.out+".dot"
+            if self.options.verbose > 1: self.out.Print("tv")
+            if self.options.verbose > 2: 
                 self.out.pdf("model_s").graphVizTree(self.options.out+".dot", "\\n")
+                print "Wrote GraphVizTree of model_s to ",self.options.out+".dot"
     def doObservables(self):
         """create pdf_bin<X> and pdf_bin<X>_bonly for each bin"""
         raise RuntimeError, "Not implemented in ModelBuilder"
     def doNuisances(self):
         if len(self.DC.systs) == 0: return
         self.doComment(" ----- nuisances -----")
         globalobs = []
-	for cpar in self.DC.discretes: self.addDiscrete(cpar) 
+        # Prepare a dictionary of which group a certain nuisance belongs to
+        groupsFor = {}
+        existingNuisanceNames = tuple(syst[0] for syst in self.DC.systs)
+        for groupName,nuisanceNames in self.DC.groups.iteritems():
+            for nuisanceName in nuisanceNames:
+                if nuisanceName not in existingNuisanceNames:
+                    raise RuntimeError, 'Nuisance group "%(groupName)s" refers to nuisance "%(nuisanceName)s" but it does not exist. Perhaps you misspelled it.' % locals()
+                if nuisanceName in groupsFor:
+                    groupsFor[nuisanceName].append(groupName)
+                else:
+                    groupsFor[nuisanceName] = [ groupName ]
+
+        #print self.DC.groups
+        #print groupsFor
+        for cpar in self.DC.discretes: self.addDiscrete(cpar)
         for (n,nofloat,pdf,args,errline) in self.DC.systs: 
             if pdf == "lnN" or pdf.startswith("shape"):
                 r = "-4,4" if pdf == "shape" else "-7,7"
@@ -225,6 +247,14 @@ def doNuisances(self):
             else: raise RuntimeError, "Unsupported pdf %s" % pdf
             if nofloat: 
               self.out.var("%s" % n).setAttribute("globalConstrained",True)
+            # set an attribute related to the group(s) this nuisance belongs to
+            if n in groupsFor:
+                groupNames = groupsFor[n]
+                if self.options.verbose > 1:
+                    print 'Nuisance "%(n)s" is assigned to the following nuisance groups: %(groupNames)s' % locals()
+                for groupName in groupNames:
+                    self.out.var(n).setAttribute('group_'+groupName,True)
+            #self.out.var(n).Print('V')
         if self.options.bin:
             nuisPdfs = ROOT.RooArgList()
             nuisVars = ROOT.RooArgSet()
@@ -242,6 +272,13 @@ def doNuisances(self):
             self.doSet("nuisances", ",".join(["%s"    % n for (n,nf,p,a,e) in self.DC.systs]))
             self.doObj("nuisancePdf", "PROD", ",".join(["%s_Pdf" % n for (n,nf,p,a,e) in self.DC.systs]))
             self.doSet("globalObservables", ",".join(globalobs))
+
+        for groupName,nuisanceNames in self.DC.groups.iteritems():
+	    nuisanceargset = ROOT.RooArgSet()
+            for nuisanceName in nuisanceNames:
+		nuisanceargset.add(self.out.var(nuisanceName))
+	    self.out.defineSet("group_%s"%groupName,nuisanceargset)
+
     def doExpectedEvents(self):
         self.doComment(" --- Expected events in each bin, for each process ----")
         for b in self.DC.bins:

scripts/combineCards.py

@@ -145,7 +145,13 @@
             bout = label if singlebin else label+b
             if isVetoed(bout,options.channelVetos): continue
 	    if not isIncluded(bout,options.channelIncludes): continue
-            obsline += [str(DC.obs[b])]; 
+            obsline += [str(DC.obs[b])];
+    #get the groups - keep nuisances in a set so that they are never repetitions
+    for groupName,nuisanceNames in DC.groups.iteritems():
+        if groupName in groups:
+            groups[groupName].update(set(nuisanceNames))
+        else:
+            groups[groupName] = set(nuisanceNames)
 
 bins = []
 for (b,p,s) in keyline:
@@ -215,6 +221,8 @@
 
 for pname in flatParamNuisances.iterkeys(): 
     print "%-12s  flatParam" % pname
-
 for dname in discreteNuisances.iterkeys(): 
-    print "%-12s  discrete" % dname
+     print "%-12s  discrete" % dname
+for groupName,nuisanceNames in groups.iteritems():
+    nuisances = ' '.join(nuisanceNames)
+    print '%(groupName)s group = %(nuisances)s' % locals()

src/Combine.cc

@@ -103,6 +103,7 @@ Combine::Combine() :
       ("setPhysicsModelParameterRanges", po::value<string>(&setPhysicsModelParameterRangeExpression_)->default_value(""), "Set the range of relevant physics model parameters. Give a colon separated list of parameter ranges. Example: CV=0.0,2.0:CF=0.0,5.0")      
       ("redefineSignalPOIs", po::value<string>(&redefineSignalPOIs_)->default_value(""), "Redefines the POIs to be this comma-separated list of variables from the workspace.")      
       ("freezeNuisances", po::value<string>(&freezeNuisances_)->default_value(""), "Set as constant all these nuisance parameters.")      
+      ("freezeNuisanceGroups", po::value<string>(&freezeNuisanceGroups_)->default_value(""), "Set as constant all these groups of nuisance parameters.")      
       ;
     ioOptions_.add_options()
       ("saveWorkspace", "Save workspace to output root file")
@@ -424,6 +425,22 @@ void Combine::run(TString hlfFile, const std::string &dataset, double &limit, do
           nuisances = mc->GetNuisanceParameters();
       }
   }
+  if (freezeNuisanceGroups_ != "") {
+      std::vector<string> nuisanceGroups;
+      boost::algorithm::split(nuisanceGroups,freezeNuisanceGroups_,boost::algorithm::is_any_of(","));
+      for (std::vector<string>::iterator ng_it=nuisanceGroups.begin();ng_it!=nuisanceGroups.end();ng_it++){
+        RooArgSet toFreeze(*(w->set(Form("group_%s",(*ng_it).c_str()))));
+        if (verbose > 0) std::cout << "Freezing the following nuisance parameters: "; toFreeze.Print("");
+        utils::setAllConstant(toFreeze, true);
+        if (nuisances) {
+          RooArgSet newnuis(*nuisances);
+          newnuis.remove(toFreeze, /*silent=*/true, /*byname=*/true);      
+          mc->SetNuisanceParameters(newnuis);
+          mc_bonly->SetNuisanceParameters(newnuis);
+          nuisances = mc->GetNuisanceParameters();
+       }
+      }
+  }
 
   if (mc->GetPriorPdf() == 0) {
       if (prior_ == "flat") {