Merge pull request cms-sw#90 from amarini/topic_tmva

Implement Tmva And Categorization

dat/configDY.dat

@@ -8,41 +8,15 @@ sub=tag|v1.4
 sub=amarini|/store/user/amarini/Nero
 sub=klute|/store/user/klute/Nero
 
+include=dat/config.dat
+
 ################
 ### DATASETS ###
 ################
 
-#### BKG ###
-Files=%(klute)s/%(tag)s/SingleMuon
-#/store/user/amarini/Nero/%(tag)s/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/WJetsToLNu_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8
-#addfiles=/store/user/amarini/Nero/v1.2/TTJets_TuneCUETP8M1_13TeV-madgraphMLM-pythia8
-##addfiles=/store/user/amarini/Nero/%(tag)s/TT_TuneCUETP8M1_13TeV-powheg-pythia8
-#### WW/WZ/ZZ###
-#addfiles=/store/user/amarini/Nero/%(tag)s/WW_TuneCUETP8M1_13TeV-pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/WZ_TuneCUETP8M1_13TeV-pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/ZZ_TuneCUETP8M1_13TeV-pythia8
-#### SIGNAL ###
-#### QCD
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_15to30_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_30to50_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_50to80_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_80to120_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_120to170_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_170to300_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_300to470_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_470to600_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_600to800_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_800to1000_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_1000to1400_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_1400to1800_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_1800to2400_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_2400to3200_TuneCUETP8M1_13TeV_pythia8
-#addfiles=/store/user/amarini/Nero/%(tag)s/QCD_Pt_3200toInf_TuneCUETP8M1_13TeV_pythia8
-#####
-##### DATA
-#addfiles=/store/user/amarini/Nero/%(tag)s/SingleMuon
-#addfiles=/store/user/amarini/Nero/%(tag)s/SingleElectron
+Files=
+include=dat/catalog/2015/dy.dat
+#addfiles=%(amarini)s/%(tag)s/SingleMuon
 
 #__________________________________________________________________
 #____________________ COMMON SETTINGS _____________________________
@@ -52,16 +26,6 @@ Output=DY.root
 ## Preselection:1,Preselection:2
 Analysis=JsonAnalysis,DYAnalysis
 
-MCDB=dat/mc_database.txt
-SFDB=aux/sf_database.txt
-Branches=dat/branches.txt
-###
-pileup=aux/pileup.root
-pileupRun=
-pileupLumi=
-#Lumi=2110 ## doing nothing
-Lumi=1
-###
 # additional JSON FILE for Data selection
 
 # Smearer
@@ -74,7 +38,6 @@ Correct=NONE
 #________________________________________________________________
 #____________________ CONFIGURATION _____________________________
 #
-config=MitPhiCorrector|fileName="aux/MetPhi.root"
 ##########  ANALYSIS  ################
 
 #############
@@ -91,8 +54,3 @@ addConfig=DYAnalysis|AddLabel('WW'),AddLabel('WZ'),AddLabel('ZZ'),AddLabel('TT')
 include=dat/json.dat
 
 
-#configure Analysis
-#config=AnalysisBase|Init()
-#config=Preselection|
-#config=MetPhiCorrector|fileName=""
-

dat/configTauNuTheory.dat

@@ -1,13 +1,12 @@
 include=dat/configTauNu.dat
 sub=tag|v1.4.lhe
-Files=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-200_13TeV_amcatnlo_pythia8
-#addfiles=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-250_13TeV_amcatnlo_pythia8
-addfiles=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-300_13TeV_amcatnlo_pythia8
-addfiles=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-350_13TeV_amcatnlo_pythia8
-addfiles=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-400_13TeV_amcatnlo_pythia8
-addfiles=%(amarini)s/%(tag)s/ChargedHiggs_HplusTB_HplusToTauNu_M-500_13TeV_amcatnlo_pythia8
-MCDB=dat/mc_db_scales.dat 
+Files=
+include=dat/catalog/2015/taunu.dat
+#MCDB=dat/mc_db_scales.dat 
 #Smear=NONE
 ## the %% is because the is a substitution before with the dictionaries
 Smear=NONE,@SmearScales("ScaleR"!1!0),@SmearScales("ScaleF"!0!1),@SmearScales("ScaleRF"!1!1),*[0!50]SmearPdfs(%%d)
 Output=ChHiggsTheory.root
+Final=true
+## be less peaky, given the n. of systs
+EntryPerSecond=10

dat/configTmvaTrainer.dat

@@ -2,11 +2,15 @@
 include=dat/config.dat
 ###########################################
 # --small set for test
-addfiles=/store/user/amarini/Nero/%(tag)s/Tau
+#addfiles=/store/user/amarini/Nero/%(tag)s/Tau
+Files=
+include=dat/catalog/2015/taunu.dat
+include=dat/catalog/2015/tt.dat
+include=dat/catalog/2015/wnjets.dat 
 
 Output=TmvaOutput.root
 Analysis=TmvaTrainer
 config=TmvaTrainer|AddLabel('HplusToTauNu_M-180_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-200_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-220_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-250_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-300_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-400_13TeV_amcatnlo'),AddLabel('HplusToTauNu_M-500_13TeV_amcatnlo')
 addConfig=TmvaTrainer|AddLabel("DY"),AddLabel("QCD"),AddLabel("TTJets"),AddLabel("WJets"),AddLabel("WW"),AddLabel("WZ"),AddLabel("ZZ")
 ## just save the tree for training
-addConfig=TmvaTrainer|train=0
+#addConfig=TmvaTrainer|train=0

interface/Event.hpp

@@ -99,17 +99,19 @@ class Event{
     //
     inline Met GetMet( ) const { return met_;} // should be const,  -- FIXED? but noCorrPt is not set correctly without &
 
-    inline float Rho() { return rho_; }
-    inline int Npv() { return npv_ ;} 
-    inline float Ht()   { float ht=0 ; for(auto j : jets_ ) if( j->IsJet()  ) ht+= j->Pt() ; return ht;}
-
-    inline int Njets(){int n=0; for(auto j : jets_) if(j->IsJet()) n++; return n;}
-    inline int NjetsInvIso(){int n=0; for(auto j : jets_) if(j->IsJetInvIso()) n++; return n;}
-    inline int NcentralJets(){int n=0; for(auto j : jets_) if(j->IsCentralJet()) n++; return n;}
-    inline int NforwardJets(){int n=0; for(auto j : jets_) if(j->IsForwardJet()) n++; return n;}
-    inline int Bjets(){int n=0; for(auto j : jets_) if(j->IsBJet()) n++; return n;}
-    inline int Ljets(){int n=0; for(auto j : jets_) if(not j->IsBJet()) n++; return n;}
-    inline int BjetsInvIso(){int n=0; for(auto j : jets_) if(j->IsBJetInvIso()) n++; return n;}
+    inline float Rho() const { return rho_; }
+    inline int Npv() const { return npv_ ;} 
+    inline float Ht() const   { float ht=0 ; for(auto j : jets_ ) if( j->IsJet()  ) ht+= j->Pt() ; return ht;}
+    inline float HtInvIso() const   { float ht=0 ; for(auto j : jets_ ) if( j->IsJetInvIso()  ) ht+= j->Pt() ; return ht;}
+
+    inline int Njets()const {int n=0; for(auto j : jets_) if(j->IsJet()) n++; return n;}
+    inline int NjetsInvIso()const {int n=0; for(auto j : jets_) if(j->IsJetInvIso()) n++; return n;}
+    inline int NcentralJets() const {int n=0; for(auto j : jets_) if(j->IsCentralJet()) n++; return n;}
+    inline int NcentralJetsInvIso() const {int n=0; for(auto j : jets_) if(j->IsCentralJetInvIso()) n++; return n;}
+    inline int NforwardJets()const {int n=0; for(auto j : jets_) if(j->IsForwardJet()) n++; return n;}
+    inline int Bjets()const {int n=0; for(auto j : jets_) if(j->IsBJet()) n++; return n;}
+    inline int Ljets()const {int n=0; for(auto j : jets_) if(not j->IsBJet()) n++; return n;}
+    inline int BjetsInvIso()const {int n=0; for(auto j : jets_) if(j->IsBJetInvIso()) n++; return n;}
     inline int Ntaus(){int n=0; for(auto t : taus_) if(t->IsTau()) n++; return n;}
     inline int Nleps(){int n=0; for(auto t : leps_) if(t->IsLep()) n++; return n;}
     inline int NGenPar(){return genparticles_.size();}

interface/GenParticle.hpp

@@ -19,6 +19,7 @@ class GenParticle : virtual public Object
     unsigned GetFlag() const { return flags_;}
     void SetFlags(unsigned x) { flags_ = x;}
     bool IsPromptFinalState() const { return flags_ & BareMonteCarlo::PromptFinalState ;}
+    bool IsDressed() const { return flags_ & BareMonteCarlo::Dressed;}
 
     int GetParentPdgId() const { return mopdgid_;}
     int GetParentIdx() const {return moidx_;}

interface/Jet.hpp

@@ -108,6 +108,12 @@ class Jet : virtual public Object, virtual public SmearableComplex
         return 1;
     }
 
+    inline int IsCentralJetInvIso() const {
+        if ( not IsJetInvIso() ) return 0;
+        if ( fabs(Eta()) >= etacutcentral_ ) return 0;
+        return 1;
+    }
+
     inline int IsForwardJet() const {
         if ( not IsJet() ) return 0;
         if ( fabs(Eta()) < etacutcentral_ ) return 0;

interface/TmvaTrainer.hpp

@@ -1,8 +1,6 @@
 #ifndef TMVA_TRAINER_H
 #define TMVA_TRAINER_H
 
-#include "TMVA/Factory.h"
-#include "TMVA/Tools.h"
 #include "TTree.h"
 
 #include "interface/AnalysisBase.hpp"
@@ -12,22 +10,17 @@
 class TmvaTrainer : virtual public AnalysisBase
 {
     private:
-        TMVA::Factory *factory_;
-        TFile *fout_;
         vector<Double_t> variables_;
 
         DataStore varsValues_;
     public:
-        TmvaTrainer();
-        ~TmvaTrainer();
+        TmvaTrainer():AnalysisBase(){}
+        ~TmvaTrainer(){}
         virtual int analyze(Event*,string systname);
         int analyzeInvIso(Event*e,string systname);
         virtual void Init();
-        virtual void End();
         virtual const string name(){return "TmvaTrainer";}
         virtual void AddVariable( string, char, double=0.,double=0.);
-
-        bool train;
 };
 
 #endif

interface/Weight.hpp

@@ -121,13 +121,13 @@ class Weight : virtual public SmearableBase {
         if (syst == MC::none and systPdf <0)
         return mcWeight_* mcXsec_ * lumi_ * sf_ * pu_.GetPUWeight(mcName_,puInt_,runNum_)/ nEvents_; 
         else if (systPdf<0) 
-        {
+        { // SCALES
         if (not scales_) {Log(__FUNCTION__,"ERROR","Scales reweighting uncorrect!!!");PrintInfo();}
         return scalesWeights_[syst] * scalesNeventReweight_[syst] * mcXsec_ * lumi_ * sf_ * pu_.GetPUWeight(mcName_,puInt_,runNum_)/ nEvents_; 
         }
         else { // pdfs
         if (not pdfs_) {Log(__FUNCTION__,"ERROR","Pdfs reweighting uncorrect!!!"); PrintInfo();}
-        return pdfsWeights_[syst] * pdfsNeventReweight_[syst] * mcXsec_ * lumi_ * sf_ * pu_.GetPUWeight(mcName_,puInt_,runNum_)/ nEvents_; 
+        return pdfsWeights_[systPdf] * pdfsNeventReweight_[systPdf] * mcXsec_ * lumi_ * sf_ * pu_.GetPUWeight(mcName_,puInt_,runNum_)/ nEvents_; 
         }
     }
     double weight(){

python/Loop.py

@@ -81,6 +81,9 @@ def applyJson(obj,fname):
 if opts.verbose:
 	PrintDat(cfg)
 
+if 'EntryPerSecond' in cfg:
+	loop.SetEntryPerSecond( cfg['EntryPerSecond'])
+
 ################ Loader ####
 if opts.verbose: print "-> Init Loader:",cfg['Loader']
 loop.InitLoader(cfg['Loader']);

python/ParseDat.py

@@ -115,6 +115,9 @@ def ParseDat(name):
 		####### FLOAT ##########
 		if key == 'Lumi':
 			config[key]=FloatKey(value)
+		####### INT ##########
+		if key == 'EntryPerSecond':
+			config[key]=IntKey(value)
 		####### V FLOAT ##########
 		if key == 'pileupLumi':
 			config[key]=vFloatKey(value)

python/mcDatabase.py

@@ -24,6 +24,7 @@
 	#if '/eos/user' in opts.eos: 
 	cmd += " --childcount "
 	cmd+=opts.eos
+	cmd += " || true" ##FIXME
 	list = check_output(cmd,shell=True);
 	for line in list.split('\n'):
 		if line =='': continue

script/buildworkspace.py

@@ -493,7 +493,7 @@ def ImportPdfStatUncFromTH1(tfile, name,syst="TTSTAT", target="", add=[]): ## w
 ## print "     FIX MH POINTS       "
 ## print "#########################"
 for syst in systSig:
- for sigMH in [ 200,250,300,350,400,500,750,1000,2000,3000]:
+ for sigMH in [ 200,250,300,350,400,500,750,1000,3000]:
  #for sigMH in [ 200,350,400,500]:
    for cat in range(0,opts.ncat):
 	sigStr="HplusToTauNu_M-"+str(sigMH)+"_13TeV_amcatnlo"
@@ -690,7 +690,7 @@ def writeNormSystList(name="lumi",valueL=["1.027","1.026"], regexpL=["TT","ST"])
 				regexp=regexp[1:]
 			if (match and not invert) or (not match and invert): 
 				idx=i
-				break
+				continue
 
 		if (idx>=0):
 		   datacard.write("\t"+valueL[idx])
@@ -716,7 +716,8 @@ def writeSyst(syst="JES"):
 	systOrig=syst[:]
 	syst=WorkspaceSubstitution(syst)
 	datacard.write(syst+" shape")
-	for proc in mcAll:
+	for cat in range(0,opts.ncat):
+	   for proc in mcAll:
 		if 'STAT' in systOrig:
 			if proc in re.sub('STAT','',systOrig):
 	        		datacard.write("\t1")

src/AnalysisChargedHiggsQCDPurity.cpp

@@ -91,16 +91,22 @@ void ChargedHiggsQCDPurity::Init()
 
         }
         // I don't need to split it by pt
-            Book( dir + "Mt"+"_"+ l  , ("Mt "+ l).c_str(),8000,0.,8000); // same binning in TauNu
-            Book( dir + "MtIsoInv"+"_"+ l  , ("MtIsoInv "+ l).c_str(),8000,0.,8000.);
-            Book( none + "EtMissIsoInv"+"_"+ l  , ("EtMissIsoInv "+ l).c_str(),1000,0.,1000.);
-            Book( none + "EtMiss"+"_"+ l  , ("EtMiss "+ l).c_str(),1000,0.,1000.); // copy of the Tau Nu ? 
+        Book( dir + "Mt"+"_"+ l  , ("Mt "+ l).c_str(),8000,0.,8000); // same binning in TauNu
+        Book( dir + "MtIsoInv"+"_"+ l  , ("MtIsoInv "+ l).c_str(),8000,0.,8000.);
+        Book( none + "EtMissIsoInv"+"_"+ l  , ("EtMissIsoInv "+ l).c_str(),1000,0.,1000.);
+        Book( none + "EtMiss"+"_"+ l  , ("EtMiss "+ l).c_str(),1000,0.,1000.); // copy of the Tau Nu ? 
 
         Book(    none+"RbbMin_"+l,"RbbMin "+l+";R_{bb}^{min}",100,0,2*TMath::Pi());
         Book(    none+"RCollMin_"+l,"RCollMin "+l+";R_{coll}^{min}",100,0,2*TMath::Pi());
 
         Book(    none+"RbbMinIsoInv_"+l,"RbbMin IsoInv "+l+";R_{bb}^{min}",100,0,2*TMath::Pi());
         Book(    none+"RCollMinIsoInv_"+l,"RCollMin IsoInv "+l+";R_{coll}^{min}",100,0,2*TMath::Pi());
+        // study categorization
+        for(int i=0;i<=1;++i)
+        {
+            Book( dir + "Mt_"+Form("cat%d",i)+"_"+ l  , ("Mt "+ l).c_str(),8000,0.,8000); // same binning in TauNu
+            Book( dir + "MtIsoInv_"+Form("cat%d",i)+"_"+ l  , ("MtIsoInv "+ l).c_str(),8000,0.,8000.);
+        }
     }
 
 }
@@ -314,6 +320,10 @@ int ChargedHiggsQCDPurity::analyze(Event*e,string systname)
             if (t->GetNProng()==1) hist+="_1p";
             else hist+="_3p";
             Fill( dir + hist +"_"+label,systname, e->GetMet().Pt(), e->weight() );
+
+            if( e->Bjets()> 1)hist = "Mt_cat0";
+            else hist = "Mt_cat1";
+            Fill(dir+hist+"_"+label,systname, e->Mt() ,e->weight(false));
     }
 
     // ---------------------- INF TAU SF 
@@ -395,6 +405,10 @@ int ChargedHiggsQCDPurity::analyze(Event*e,string systname)
             hist = HistName(pt,false, true,"Upar");
             Fill( dir+hist +"_"+label,systname, Uperp(e,tInv), e->weight(false) );
 
+            if( e->Bjets()> 1)hist = "MtIsoInv_cat0";
+            else hist = "MtIsoInv_cat1";
+            Fill(dir+hist+"_"+label,systname, e->Mt(Event::MtTauInv) ,e->weight(false));
+
             hist = HistName(pt, false, true) ;
             if (flavor == 15) hist += "_T";
             else if (flavor == 21 ) hist+="_G";