Merge pull request cms-sw#25 from amarini/topic_combine

Combine

Makefile

@@ -25,6 +25,27 @@ all:
 fast:
 	$(MAKE) libChargedHiggs.so
 
+# check if CMSSW is defined
+ifndef CMSSW_BASE
+$(info No CMSSSW !!!!)
+$(info I ll sleep 3s to let you acknowledge it)
+$(shell sleep 3s)
+else
+$(info CMSSW found: $(CMSSW_BASE) )
+COMBINELIBFILE = $(wildcard $(CMSSW_BASE)/lib/$(SCRAM_ARCH)/libHiggsAnalysisCombinedLimit.so)
+COMBINELIB = HiggsAnalysisCombinedLimit
+COMBINELIBDIR = $(CMSSW_BASE)/lib/$(SCRAM_ARCH)/
+endif
+
+# check if Combine is present and compiled 
+ifeq ("$(COMBINELIBFILE)", "") 
+$(info No Combine Package found)
+else
+$(info Using combine: $(COMBINELIB))
+CXXFLAGS += -L$(ROOFITSYS)/lib -lRooFit -lRooFitCore -I$(ROOFITSYS)/include
+CXXFLAGS += -D HAVE_COMBINE -L$(COMBINELIBDIR) -l$(COMBINELIB) -Wl,-rpath=$(COMBINELIBDIR)
+endif
+
 libChargedHiggs.so: $(OBJ) Dict | $(BINDIR)
 	$(GCC) $(CXXFLAGS) $(SOFLAGS) -o $(BINDIR)/$@ $(OBJ) $(BINDIR)/dict.o
 
@@ -35,7 +56,7 @@ $(OBJ) : $(BINDIR)/%.o : $(SRCDIR)/%.cpp interface/%.hpp | $(BINDIR)
 Dict: $(BINDIR)/dict.o
 
 $(BINDIR)/dict.o: $(SRC) | $(BINDIR)
-	cd $(BINDIR) && rootcint -v4 -f dict.cc -c -I../../ -I../ $(HPPLINKDEF)  ../interface/LinkDef.hpp 
+	cd $(BINDIR) && rootcint -v4 -f dict.cc -c -I../../ -I../  $(HPPLINKDEF)  ../interface/LinkDef.hpp 
 	cd $(BINDIR) && $(GCC) -c -o dict.o $(CXXFLAGS) -I../../ dict.cc
 
 $(BINDIR):

datacards/cms_datacard_chhiggs_taunu.txt

@@ -0,0 +1,31 @@
+# shapes process channel file histogram [ histogram_with_systematics ] 
+# PROCESS is replaced with the process name (or "data_obs" for the observed data)
+# $CHANNEL is replaced with the channel name
+# $SYSTEMATIC is replaced with the name of the systematic + (Up, Down)
+# $MASS is replaced with the higgs mass value which is passed as option in the command line used to run the limit tool 
+-------------------------------------
+imax *
+jmax *
+kmax *
+-------------------------------------
+shapes data_obs cat0  ChHiggs_interpolated.root  Mt_Data
+shapes Hplus 	cat0  ChHiggs_interpolated.root  Mt_HplusToTauNu-M$MASS_interpolated
+#shapes *  	cat0  ChHiggs_interpolated.root  Mt_$PROCESS
+shapes QCD  	cat0  ChHiggs_interpolated.root  Mt_QCD
+shapes TTJets  	cat0  ChHiggs_interpolated.root  Mt_TTJets
+shapes ZZ  	cat0  ChHiggs_interpolated.root  Mt_ZZ
+shapes WW  	cat0  ChHiggs_interpolated.root  Mt_WW
+shapes WZ  	cat0  ChHiggs_interpolated.root  Mt_WZ
+shapes DY  	cat0  ChHiggs_interpolated.root  Mt_DY
+-------------------------------------
+bin cat0
+observation -1
+-------------------------------------
+bin 	cat0	cat0	cat0	cat0 	cat0	cat0	cat0
+process Hplus	QCD	TTJets	WW	WZ	DY	ZZ
+process 0	1	2	3	4	5	6
+rate	-1	-1	-1	-1	-1	-1	-1
+-------------------------------------
+
+
+

interface/Fitter.hpp

@@ -0,0 +1,124 @@
+#ifndef FITTER_H
+#define FITTER_H
+
+// --- STD ---
+#include <string>
+#include <vector>
+#include <algorithm>
+using namespace std;
+// --- ROOT ---
+#include "TFile.h"
+#include "TDirectory.h"
+#include "TH1D.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include "TROOT.h" 
+// --- ROOFIT ---
+#include "RooGlobalFunc.h"
+#include "RooPlot.h"
+#include "RooWorkspace.h"
+#include "RooGaussian.h"
+#include "RooVoigtian.h"
+#include "RooProduct.h"
+#include "RooAddition.h"
+#include "RooConstVar.h"
+#include "RooFormulaVar.h"
+#include "RooDataHist.h"
+#include "RooRealVar.h"
+#include "RooAddPdf.h"
+// do not use namespace for dictionaries
+//using namespace RooFit;
+class RooSpline1D;
+
+
+
+class Fitter{
+    /* Original Author: Andrea C. Marini 
+     * Date: 8th June 2015
+     * This class provide an interface to create 
+     * smooth workspace for combine
+     */
+
+    // private members end with _
+    //
+    RooWorkspace *w_;
+    RooRealVar *mh_; // truth
+    RooRealVar *x_; // observable
+
+    map< string, RooDataHist*> hist_;
+
+    string datasetMask_ ;
+    string xsecMask_;
+    string massMask_;
+
+    bool writeDatasets_ ; // write the RooDatasets into the ws
+
+    vector<float> startMean_;
+    vector<float> startSigma_;
+    vector<float> startFraction_;
+    vector<float> startBern_;
+
+    // save fit parameters vs mass category 
+    map<string, float> fitParameters_;
+    // interpolates fit parameters and norm vs mass
+    map<string, RooSpline1D*> splines_;
+    // save RooRealVars used for fit
+    map<string, RooRealVar*> vars_;
+
+    // -- RooAbsReal* getMeanWithSyst(string name, RooAbsReal*mean);
+    void info();
+
+    // ------------ INIT -------
+    void initGaus();
+    void initBern();
+
+
+    // ------------SIGNAL MODEL AND FIT -----------
+    // add Bern - Gaus To the fit Model
+    void addBernFitModel(RooArgList *pdfs,RooArgList *coeffs,bool isLast);
+    void addGausFitModel(RooArgList *pdfs,RooArgList *coeffs,bool isLast);
+
+    // -- save Coeff of the fit in fitParameters_
+    void saveCoefficientsGaus( int cat,string mass,bool isLast);
+    void saveCoefficientsBern( int cat,string mass,bool isLast);
+
+    void interpolateBern(int cat,bool isLast);
+    void interpolateGaus(int cat,bool isLast);
+
+    // --- FINAL MODEL ---
+
+    public:
+
+    // --- objects that can be set	
+    vector<float> mIn; // input masses
+    vector<string> inputMasks; // input FileName Mask. Must contain a replacement for float. One for each cat
+    string outputFileName; // 
+    string inputFileName;
+    string plotDir;
+    bool plot; // make plot fit
+    bool verbose;
+
+    int nGaussians;
+    int nBernstein;
+
+    float xmin;
+    float xmax;
+
+    // --- objects that can be called
+    Fitter();
+    void init();	
+    void fitSignal();
+    void write();
+    void finalModel();
+
+};
+
+// Local Variables:
+// mode:c++
+// indent-tabs-mode:nil
+// tab-width:4
+// c-basic-offset:4
+// End:
+// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 
+
+#endif

interface/LinkDef.hpp

@@ -13,6 +13,9 @@
 #pragma link C++ class Event+ ;
 #pragma link C++ class Looper+;
 #pragma link C++ class SmearBase+;
+
+#pragma link C++ class Fitter+;
+
 //#pragma link C++ class sigint_exception+;
 
 #endif

python/ParseDat.py

@@ -156,6 +156,7 @@ def PrintUsage():
 	print 'pileup = file'
 	print 'pileupRun  = ""/ -1 / 1,10,100,100'
 	print 'pileupLumi = ""/ -1 /  1.2,2.3,4.5 # partial luminosity '
+	print 'branches = brancfile'
 	print 'Analysis = AnalysisBase,Analysis2 ..'
 	print 'Smears = @SmearBase,JER,JES'
 	print 'config = AnalysisBase|a=1,b=2,c(3)'

python/fitter.py

@@ -0,0 +1,59 @@
+import os,sys,re,time
+from glob import glob
+from subprocess import call
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-v","--verbose",dest="verbose",action='store_true',help="Verbose",default=False)
+
+(opts,args) = parser.parse_args()
+
+#  find the base directory
+
+if opts.verbose: print "-> Looking for basepath"
+basepath = ""
+mypath = os.path.abspath(os.getcwd())
+while mypath != "" and mypath != "/":
+	if "ChargedHiggs" in os.path.basename(mypath):
+		basepath = os.path.abspath(mypath)
+	mypath = os.path.dirname(mypath)
+################ IMPORT ROOT  ################
+if opts.verbose: print "-> Importing root",
+sys.argv=[]
+import ROOT as r
+r.gROOT.SetBatch()
+if opts.verbose: print "DONE"
+
+################ LOAD LIBRARY ###############
+if opts.verbose: print "-> Load Bare library"
+r.gSystem.Load( "../NeroProducer/Core/bin/libBare.so")
+if opts.verbose: print "-> Load ChargedHiggs library"
+r.gSystem.Load( "./bin/libChargedHiggs.so")
+if opts.verbose: print "DONE",
+
+################ CREATING FITTER ##########
+from ROOT import Fitter
+fitter = Fitter()
+
+################ CONFIGURING FITTER ##########
+
+call( "mkdir -p plot/sigfit", shell=True)
+fitter.plotDir = "plot/sigfit/"
+fitter.nGaussians = 1
+fitter.nBernstein = 3
+fitter.xmin = 150
+fitter.mIn.push_back(200)
+fitter.mIn.push_back(250)
+fitter.mIn.push_back(500)
+fitter.mIn.push_back(900)
+fitter.outputFileName = "sigfit.root"
+fitter.inputFileName = "ChHiggs.root"
+
+################ INIT FITTER ##########
+if opts.verbose: print "-> Init"
+fitter.init()
+if opts.verbose: print "-> Fit"
+fitter.fitSignal()
+if opts.verbose: print "-> FinalModel"
+fitter.finalModel()
+if opts.verbose: print "-> Write"
+fitter.write()

script/combine.py

@@ -0,0 +1,119 @@
+#!env python
+
+import os,sys
+from glob import glob
+from subprocess import call, check_output
+import re
+
+def drange(start, stop, step):
+        ''' Return a floating range list. Start and stop are included in the list if possible.'''
+        r = start
+        while r <= stop:
+                yield r
+                r += step
+
+from optparse import OptionParser, OptionGroup
+
+usage = "usage: %prog [options]"
+parser=OptionParser(usage=usage)
+parser.add_option("-i","--input" ,dest='input',type='string',help="Input Datacard [Default=%default]",default="")
+parser.add_option("-d","--dir" ,dest='dir',type='string',help="Directory where to write the configuration [Default=%default]",default="submit")
+#parser.add_option("-n","--njobs" ,dest='njobs',type='int',help="Number of Job to submit [Default=%default]",default=50)
+parser.add_option("","--begin" ,dest='begin',type='float',help="Begin Mass [Default=%default]",default=200)
+parser.add_option("","--end" ,dest='end',type='float',help="End Mass [Default=%default]",default=900)
+parser.add_option("","--step" ,dest='step',type='float',help="Step [Default=%default]",default=100)
+parser.add_option("-q","--queue" ,dest='queue',type='string',help="Queue [Default=%default]",default="1nh")
+parser.add_option("","--dryrun" ,dest='dryrun',action='store_true',help="Do not Submit [Default=%default]",default=False)
+parser.add_option("-e","--expected" ,dest='exp',action='store_true',help="Run Expected [Default=%default]",default=True)
+parser.add_option("-u","--unblind" ,dest='exp',action='store_false',help="Unblind Results")
+
+(opts,args)=parser.parse_args()
+
+EOS='/afs/cern.ch/project/eos/installation/0.3.84-aquamarine/bin/eos.select'
+
+print "inserting in path cwd"
+sys.path.insert(0,os.getcwd())
+print "inserting in path cwd/python"
+sys.path.insert(0,os.getcwd()+'/python')
+
+print "inserting combine in path"
+if not os.path.exists(os.path.expandvars('$CMSSW_BASE/bin/$SCRAM_ARCH/combine')):
+        sys.exit('ERROR - CombinedLimit package must be installed')
+if not os.path.exists(os.path.expandvars('$CMSSW_BASE/bin/$SCRAM_ARCH/text2workspace.py')):
+        sys.exit('ERROR - CombinedLimit package must be installed')
+if not os.path.exists(os.path.expandvars('$CMSSW_BASE/bin/$SCRAM_ARCH/combineCards.py')):
+        sys.exit('ERROR - CombinedLimit package must be installed')
+
+
+call("[ -d %s ] && rmdir %s"%(opts.dir,opts.dir),shell=True)
+call("mkdir -p %s"%opts.dir,shell=True)
+
+#for iJob in range(0,opts.njobs):
+cmdFile = open("%s/cmdFile.sh"%opts.dir,"w")
+iJob = -1
+
+for mass in drange(opts.begin,opts.end,opts.step):
+	iJob += 1
+	basedir=opts.dir
+	if basedir[0] != '/': basedir=os.environ['PWD'] + "/" + opts.dir
+
+	sh=open("%s/sub%d.sh"%(opts.dir,iJob),"w")
+
+	call(["chmod","u+x","%s/sub%d.sh"%(opts.dir,iJob)])
+
+	## Open File and write Preamble
+	sh.write('#!/bin/bash\n')
+	sh.write('[ "$WORKDIR" == "" ] && export WORKDIR="/tmp/%s/" \n'%(os.environ['USER']))
+	sh.write('cd %s\n'%(os.getcwd() ) )
+	sh.write('LD_LIBRARY_PATH=%s:$LD_LIBRARY_PATH\n'%os.getcwd())
+
+	print " WORK-AROUND COMBINE: FIXME!!! "
+	sh.write('cd ~amarini/work/ProductionJanuary2014/CMSSW_6_1_1_CategoryFull/src\n')
+	sh.write('eval `scramv1 runtime -sh`\n') # cmsenv
+
+	## Touch control files
+	touch = "touch " + basedir + "/sub%d.pend"%iJob
+	call(touch,shell=True)
+	cmd = "rm " + basedir + "/sub%d.run 2>/dev/null >/dev/null"%iJob
+	call(cmd,shell=True)
+	cmd = "rm " + basedir + "/sub%d.done 2>/dev/null >/dev/null"%iJob
+	call(cmd,shell=True)
+	cmd = "rm " + basedir + "/sub%d.fail 2>/dev/null >/dev/null"%iJob
+	call(cmd,shell=True)
+
+	sh.write('date > %s/sub%d.run\n'%(basedir,iJob))
+	sh.write('rm %s/sub%d.done\n'%(basedir,iJob))
+	sh.write('rm %s/sub%d.pend\n'%(basedir,iJob))
+	sh.write('rm %s/sub%d.fail\n'%(basedir,iJob))
+
+	sh.write('cd $WORKDIR\n')
+
+	## Construct Datacard
+	datacard= re.sub( ".txt","_M%.0f.root"%mass,opts.input.split("/")[-1])
+	cmd = "text2workspace.py -m " + str(mass) + " -o " + opts.dir + "/" + datacard +  " " + opts.input
+	print "-> Running command :'"+cmd+"'"
+	call(cmd,shell=True)
+
+	## copy datacard in workdir
+	sh.write('cp -v '+ basedir + "/"+  datacard + " ./ \n" )
+	##Write combine line
+	#combine -M Asymptotic -m 200 -S 0 --run=expected --expectSignal=1 --expectSignalMass=200  cms_datacard_chhiggs_taunu.txt
+	combine = "combine -M Asymptotic -m "+ str(mass) + " -S 0 "
+	combine += "  --cminDefaultMinimizerType=Minuit2 "
+	if opts.exp : combine += " --run=expected --expectSignal=1 --expectSignalMass="+str(mass) + " "
+	combine += " --rMax=1000000 "
+	combine += datacard
+	sh.write( combine + "\n" )
+	sh.write('EXITCODE=$?\n')
+	sh.write('[ $EXITCODE == 0 ] && touch %s/sub%d.done\n'%(basedir,iJob))
+	sh.write('[ $EXITCODE != 0 ] && echo $EXITCODE > %s/sub%d.fail\n'%(basedir,iJob))
+
+	sh.write('cp higgs*root %s/\n'%basedir)
+	sh.write('rm %s/sub%d.run\n'%(basedir,iJob))
+
+	cmdline = "bsub -q " + opts.queue + " -o %s/log%d.txt"%(basedir,iJob) + " -J " + "%s/Job_%d"%(opts.dir,iJob) + " %s/sub%d.sh"%(basedir,iJob)
+	print cmdline
+	cmdFile.write(cmdline+"\n")
+
+	if not opts.dryrun: 
+		call(cmdline,shell=True)