Merge pull request #20 from gpetruc/TTHAnalysis-SUSYScan-CMG_PAT_V5_1…

…8_from-CMSSW_5_3_14

TTHAnalysis for Susy: analyzer to log into the tree the masses of the Susy particles

CMGTools/TTHAnalysis/cfg/run_ttHLep8TeV_newNtpl_cfg.py

@@ -58,6 +58,10 @@
     PDFWeights = [ pdf for pdf,num in PDFWeights ]
     )
 
+susyScanAna = cfg.Analyzer(
+    'susyParameterScanAnalyzer',
+    )
+
 # Lepton Analyzer
 ttHLepAna = cfg.Analyzer(
     'ttHLepAnalyzerBase',
@@ -179,6 +183,7 @@
     triggerAna,
     pileUpAna,
     ttHGenAna,
+    susyScanAna,
     ttHVertexAna,
     ttHLepAna,
     ttHLepMCAna,

CMGTools/TTHAnalysis/python/analyzers/susyParameterScanAnalyzer.py

@@ -0,0 +1,99 @@
+from CMGTools.RootTools.fwlite.Analyzer import Analyzer
+#from CMGTools.RootTools.fwlite.Event import Event
+#from CMGTools.RootTools.statistics.Counter import Counter, Counters
+from CMGTools.RootTools.fwlite.AutoHandle import AutoHandle
+from math import floor
+import re
+
+class susyParameterScanAnalyzer( Analyzer ):
+    """Get information for susy parameter scans    """
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        super(susyParameterScanAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+        self.susyParticles = {
+            100001 : 'Squark',
+            (1000000 + 21) : 'Gluino',
+            (1000000 + 39) : 'Gravitino',
+            (1000000 + 5) : 'Sbottom',
+            (2000000 + 5) : 'Sbottom2',
+            (1000000 + 6) : 'Stop',
+            (2000000 + 6) : 'Stop2',
+            (1000000 + 22) : 'Neutralino',
+            (1000000 + 23) : 'Neutralino2',
+            (1000000 + 25) : 'Neutralino3',
+            (1000000 + 35) : 'Neutralino4',
+            (1000000 + 24) : 'Chargino',
+            (1000000 + 37) : 'Chargino2',
+        }
+    #---------------------------------------------
+    # DECLARATION OF HANDLES OF GEN LEVEL OBJECTS 
+    #---------------------------------------------
+
+
+    def declareHandles(self):
+        super(susyParameterScanAnalyzer, self).declareHandles()
+
+        #mc information
+        self.mchandles['genParticles'] = AutoHandle( 'genParticlesPruned',
+                                                     'std::vector<reco::GenParticle>' )
+        self.mchandles['lhe'] = AutoHandle( 'source', 'LHEEventProduct' )
+
+    def beginLoop(self):
+        super(susyParameterScanAnalyzer,self).beginLoop()
+
+
+    def findSusyMasses(self,event):
+        masses = {}
+        for p in event.genParticles:
+            if p.status() != 3: continue
+            id = abs(p.pdgId())
+            if (id / 1000000) % 10 in [1,2]:
+                particle = None
+                if id % 100 in [1,2,3,4]:
+                    particle = "Squark"
+                elif id in self.susyParticles:
+                    particle = self.susyParticles[id]
+                if particle != None:
+                    if particle not in masses: masses[particle] = []
+                    masses[particle].append(p.mass())
+        for p,ms in masses.iteritems():
+            avgmass = floor(sum(ms)/len(ms)+0.5)
+            setattr(event, "genSusyM"+p, avgmass)
+
+    def readLHE(self,event):
+        lheprod = self.mchandles['lhe'].product()
+        scanline = re.compile(r"#\s*model\s+([A-Za-z0-9]+)_((\d+\.?\d*)(_\d+\.?\d*)*)\s+(\d+\.?\d*)\s*")
+        for i in xrange(lheprod.comments_size()):
+            comment = lheprod.getComment(i)
+            m = re.match(scanline, comment) 
+            if m:
+                event.susyModel = m.group(1)
+                masses = [float(x) for x in m.group(2).split("_")]
+                if len(masses) >= 1: event.genSusyMScan1 = masses[0]
+                if len(masses) >= 2: event.genSusyMScan2 = masses[1]
+                if len(masses) >= 3: event.genSusyMScan3 = masses[2]
+                if len(masses) >= 4: event.genSusyMScan4 = masses[3]
+            elif "model" in comment:
+                if not hasattr(self,"warned_already"):
+                    print "ERROR: I can't understand the model: ",comment
+                    self.warned_already = True
+
+    def process(self, iEvent, event):
+        # if not MC, nothing to do
+        if not self.cfg_comp.isMC: 
+            return True
+
+        self.readCollections( iEvent )
+
+        # create parameters
+        event.susyModel = None
+        for id,X in self.susyParticles.iteritems():
+            setattr(event, "genSusyM"+X, -99.0)
+        event.genSusyMScan1 = 0.0
+        event.genSusyMScan2 = 0.0
+        event.genSusyMScan3 = 0.0
+        event.genSusyMScan4 = 0.0
+
+        # do MC level analysis
+        self.readLHE(event)
+        self.findSusyMasses(event)
+        return True

CMGTools/TTHAnalysis/python/analyzers/ttHLepTreeProducerTTH.py

@@ -52,6 +52,24 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
             NTupleVariable("LepEff_2lep", lambda ev : ev.LepEff_2lep, mcOnly=True, help="Lepton preselection SF (2 lep)"),
             NTupleVariable("LepEff_3lep", lambda ev : ev.LepEff_3lep, mcOnly=True, help="Lepton preselection SF (3 lep)"),
             NTupleVariable("LepEff_4lep", lambda ev : ev.LepEff_4lep, mcOnly=True, help="Lepton preselection SF (4 lep)"),
+            ##------------------------------------------------
+            NTupleVariable("GenSusyMScan1", lambda ev : ev.genSusyMScan1, int, mcOnly=True, help="Susy mass 1 in scan"),
+            NTupleVariable("GenSusyMScan2", lambda ev : ev.genSusyMScan2, int, mcOnly=True, help="Susy mass 2 in scan"),
+            NTupleVariable("GenSusyMScan3", lambda ev : ev.genSusyMScan3, int, mcOnly=True, help="Susy mass 3 in scan"),
+            NTupleVariable("GenSusyMScan4", lambda ev : ev.genSusyMScan4, int, mcOnly=True, help="Susy mass 4 in scan"),
+            NTupleVariable("GenSusyMGluino", lambda ev : ev.genSusyMGluino, int, mcOnly=True, help="Susy Gluino mass"),
+            NTupleVariable("GenSusyMGravitino", lambda ev : ev.genSusyMGravitino, int, mcOnly=True, help="Susy Gravitino mass"),
+            NTupleVariable("GenSusyMStop", lambda ev : ev.genSusyMStop, int, mcOnly=True, help="Susy Stop mass"),
+            NTupleVariable("GenSusyMSbottom", lambda ev : ev.genSusyMSbottom, int, mcOnly=True, help="Susy Sbottom mass"),
+            NTupleVariable("GenSusyMStop2", lambda ev : ev.genSusyMStop2, int, mcOnly=True, help="Susy Stop2 mass"),
+            NTupleVariable("GenSusyMSbottom2", lambda ev : ev.genSusyMSbottom2, int, mcOnly=True, help="Susy Sbottom2 mass"),
+            NTupleVariable("GenSusyMSquark", lambda ev : ev.genSusyMSquark, int, mcOnly=True, help="Susy Squark mass"),
+            NTupleVariable("GenSusyMNeutralino", lambda ev : ev.genSusyMNeutralino, int, mcOnly=True, help="Susy Neutralino mass"),
+            NTupleVariable("GenSusyMNeutralino2", lambda ev : ev.genSusyMNeutralino2, int, mcOnly=True, help="Susy Neutralino2 mass"),
+            NTupleVariable("GenSusyMNeutralino3", lambda ev : ev.genSusyMNeutralino3, int, mcOnly=True, help="Susy Neutralino3 mass"),
+            NTupleVariable("GenSusyMNeutralino4", lambda ev : ev.genSusyMNeutralino4, int, mcOnly=True, help="Susy Neutralino4 mass"),
+            NTupleVariable("GenSusyMChargino", lambda ev : ev.genSusyMChargino, int, mcOnly=True, help="Susy Chargino mass"),
+            NTupleVariable("GenSusyMChargino2", lambda ev : ev.genSusyMChargino2, int, mcOnly=True, help="Susy Chargino2 mass"),
         ]
 
         self.globalObjects = {

SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h

@@ -29,6 +29,11 @@ class LHEEventProduct {
 	size_type comments_size() const { return comments_.size(); }
 	comments_const_iterator comments_begin() const { return comments_.begin(); }
 	comments_const_iterator comments_end() const { return comments_.end(); }
+
+  const char* getComment(unsigned i) const {
+    if(comments_.size()<1 || i>=comments_.size()) return "";
+    else return (const char*) comments_[i].c_str();
+  }
 
 	class const_iterator {
 	    public: