Merge pull request #162 from susy2015/nanoAODUpdates

Nano aod updates

TopTagger/interface/Constituent.h

@@ -9,16 +9,18 @@
 #include <map>
 #include <string>
 
-enum ConstituentType
-{
-    NOTYPE, AK4JET, AK6JET, AK8JET, CA8JET, AK8SUBJET, RESOLVEDTOPCAND
-};
-
 /**
  *Serves as a container class for jet objects including the 4-vector for each jet along with supporting information such as the b-tag discriminator and jet type
  */
 class Constituent
 {
+public:
+    enum ConstituentType
+    {
+        NOTYPE, AK4JET, AK6JET, AK8JET, CA8JET, AK8SUBJET, RESOLVEDTOPCAND
+    };
+
+
 private:
     TLorentzVector p_;
     ConstituentType type_;
@@ -46,15 +48,15 @@ class Constituent
     /** Construct an AK4 jet from TLorentzVector, b-tag discriminator, and quark-gluonlikelihood */
     Constituent(const TLorentzVector& p, const double& bTagDisc, const double& qgLikelihood);
     /** Construct generic constituent jet from TLorentzVector and type */
-    Constituent(const TLorentzVector& p, const ConstituentType& type);
+    Constituent(const TLorentzVector& p, const Constituent::ConstituentType& type);
     /** Construct an AK8 jet from TLorentzVector, Nsubjettiness inputs tau1/2/3, softdrop mass, softdrop subjets vector and wMassCorrection (if applicable). */
     Constituent(const TLorentzVector& p, const double& tau1, const double& tau2, const double& tau3, const double& softDropMass, const std::vector<Constituent>& subjets, const double& wMassCorr);
 
     void setPBtag(const TLorentzVector& p, const double& bTagDisc, const double& qgLikelihood);
     void setP(const TLorentzVector& p);
     void setBTag(const double&  bTagDisc);
     void setQGLikelihood(const double& qgLikelihood);
-    void setType(const ConstituentType type);
+    void setType(const Constituent::ConstituentType type);
     void setTau1(const double& tau1);
     void setTau2(const double& tau2);
     void setTau3(const double& tau3);

TopTagger/interface/TopTaggerUtilities.h

@@ -296,7 +296,7 @@ namespace ttUtility
                 {	        
                     for(unsigned int iSJ = 0; iSJ < (*vecSubjetsLVec_)[iJet].size(); ++iSJ)
                     {
-                        subjets.emplace_back((*vecSubjetsLVec_)[iJet][iSJ], AK8SUBJET);
+                        subjets.emplace_back((*vecSubjetsLVec_)[iJet][iSJ], Constituent::AK8SUBJET);
                         if(vecSubjetsBtag_)  subjets.back().setBTag((*vecSubjetsBtag_)[iJet][iSJ]);
                         if(vecSubjetsMult_)  subjets.back().setExtraVar("mult",  (*vecSubjetsMult_) [iJet][iSJ]);
                         if(vecSubjetsPtD_)   subjets.back().setExtraVar("ptD",   (*vecSubjetsPtD_)  [iJet][iSJ]);
@@ -312,7 +312,7 @@ namespace ttUtility
                         double myDR = ROOT::Math::VectorUtil::DeltaR((*jetsLVec_)[iJet], (*subjetsLVec_)[iSJ]);
                         if (myDR < 0.8)
                         {
-                            subjets.emplace_back((*subjetsLVec_)[iSJ], AK8SUBJET);
+                            subjets.emplace_back((*subjetsLVec_)[iSJ], Constituent::AK8SUBJET);
                             if(subjetsBtag_)  subjets.back().setBTag((*subjetsBtag_)[iSJ]);
                             if(subjetsMult_)  subjets.back().setExtraVar("mult", (*subjetsMult_)[iSJ]);
                             if(subjetsPtD_)   subjets.back().setExtraVar("ptD", (*subjetsPtD_)[iSJ]);
@@ -471,7 +471,7 @@ namespace ttUtility
             //Fill the constituent with the necessary information 
             for(unsigned int iTop = 0; iTop < topCandLVec_->size(); ++iTop)
             {
-                constituents.emplace_back((*topCandLVec_)[iTop], RESOLVEDTOPCAND);
+                constituents.emplace_back((*topCandLVec_)[iTop], Constituent::RESOLVEDTOPCAND);
                 auto& constituent = constituents.back();
                 constituent.setTopDisc((*topCandDisc_)[iTop]);
                 constituent.addJetRefIndex((*topCandJ1_)[iTop]);
@@ -517,11 +517,6 @@ namespace ttUtility
     ///Tool to calcualte MT2 from tagger results
     double calculateMT2(const TopTaggerResults& ttr, const TLorentzVector& metLVec);
 
-    //MVA helper functions
-    /* \fn std::map<std::string, double> ttUtility::createMVAInputs(const TopObject& topCand, const double csvThresh) */
-    ///Function to calculate MVA input variables from a TopObject. This function is depricated in favor of the new MVAInputCalculator based approach.
-    std::map<std::string, double> createMVAInputs(const TopObject& topCand, const double csvThresh);
-
     //New MVA variable helper class
     /**
      *Base class for MVA input variable calculator 

TopTagger/python/TopTagger.py

@@ -1,17 +1,23 @@
 import TopTaggerInterface as tti
 
 class Top:
-    def __init__(self, pt, eta, phi, mass, disc, type):
+    def __init__(self, pt, eta, phi, mass, disc, type, j1Idx = -999, j2Idx = -999, j3Idx = -999):
         self.pt = pt
         self.eta = eta
         self.phi = phi
         self.mass = mass
         self.disc = disc
         self.type = type
+        self.j1Idx = j1Idx
+        self.j2Idx = j2Idx
+        self.j3Idx = j3Idx
 
     def __str__(self):
         return "Top properties:   pt:  %7.3f,    eta:  %7.3f,    phi:  %7.3f,    mass:  %7.3f,    disc:  %7.3f,    type:  %3i"%(self.pt, self.eta, self.phi, self.mass, self.disc, self.type)
 
+    def __repr__(self):
+        return "Top(%f, %f, %f, %f, %f, %i, %i, %i, %i)"%(self.pt, self.eta, self.phi, self.mass, self.disc, self.type, self.j1Idx, self.j2Idx, self.j3Idx)
+
 class TopTaggerResult:
     def __init__(self, results):
         self.floatVals = results[0]
@@ -21,7 +27,7 @@ def __len__(self):
         return self.floatVals.shape[0]
 
     def __iter__(self):
-        for variables in zip(self.ptCol(), self.etaCol(), self.phiCol(), self.massCol(), self.discCol(), self.typeCol()):
+        for variables in zip(self.ptCol(), self.etaCol(), self.phiCol(), self.massCol(), self.discCol(), self.typeCol(), self.j1IdxCol(), self.j2IdxCol(), self.j3IdxCol()):
             yield Top(*variables)
 
     def ptCol(self):
@@ -42,6 +48,15 @@ def discCol(self):
     def typeCol(self):
         return self.intVals[:, 0]
 
+    def j1IdxCol(self):
+        return self.intVals[:, 1]
+
+    def j2IdxCol(self):
+        return self.intVals[:, 2]
+
+    def j3IdxCol(self):
+        return self.intVals[:, 3]
+
 class TopTagger:
 
     def __init__(self, cfgFile, workingDir = ""):
@@ -75,12 +90,11 @@ def run(self, *args, **kwargs):
         results = tti.getResults(self.tt)
         return TopTaggerResult(results)
 
-    def runFromNanoAOD(self, event):
+    def runFromNanoAOD(self, event, isFirstEvent = False):
         #This is a hack for the nanoAOD postprocessor to force it to read all necessary variables before passing them to C because each new branch accessed causes all branches to be reallocated 
         nHackLoop = 1
-        if self.firstEvent:
+        if isFirstEvent:
             nHackLoop = 2
-            self.firstEvent = False
 
         for i in xrange(nHackLoop):
 

TopTagger/src/TTMConstituentReqs.cpp

@@ -39,7 +39,7 @@ void TTMConstituentReqs::getParameters(const cfg::CfgDocument* cfgDoc, const std
 bool TTMConstituentReqs::passAK8WReqs(const Constituent& constituent) const
 {
     //check that it is an AK8 jet
-    if(constituent.getType() != AK8JET) return false;
+    if(constituent.getType() != Constituent::AK8JET) return false;
 
     //check that tau1 and 2 are valid
     if(constituent.getTau1() <= 0 || constituent.getTau2() <= 0) return false;
@@ -55,7 +55,7 @@ bool TTMConstituentReqs::passAK8WReqs(const Constituent& constituent) const
 bool TTMConstituentReqs::passAK4WReqs(const Constituent& constituent, const Constituent& constituentAK8) const
 {
     //basic AK4 jet requirements 
-    bool basicReqs = constituent.getType() == AK4JET && constituent.p().Pt() > minAK4WPt_;
+    bool basicReqs = constituent.getType() == Constituent::AK4JET && constituent.p().Pt() > minAK4WPt_;
 
     //check that the AK4 jet does not overlap with the selected AK8 subjets
     //to keep the algorithm from needing to do unnecessary calculations on all matches
@@ -68,7 +68,7 @@ bool TTMConstituentReqs::passAK4WReqs(const Constituent& constituent, const Cons
 bool TTMConstituentReqs::passAK8TopReqs(const Constituent& constituent) const
 {
     //check that it is an AK8 jet
-    if(constituent.getType() != AK8JET) return false;
+    if(constituent.getType() != Constituent::AK8JET) return false;
 
     //check that tau2 and 3 are valid
     if(constituent.getTau2() <= 0 || constituent.getTau3() <= 0) return false;
@@ -84,7 +84,7 @@ bool TTMConstituentReqs::passAK8TopReqs(const Constituent& constituent) const
 bool TTMConstituentReqs::passDeepAK8WReqs(const Constituent& constituent) const
 {
     //check that it is an AK8 jet
-    if(constituent.getType() != AK8JET) return false;
+    if(constituent.getType() != Constituent::AK8JET) return false;
 
     return constituent.p().Pt() > minAK8WPt_ &&
            constituent.getSoftDropMass() > minAK8WMass_  && 
@@ -95,7 +95,7 @@ bool TTMConstituentReqs::passDeepAK8WReqs(const Constituent& constituent) const
 bool TTMConstituentReqs::passDeepAK8TopReqs(const Constituent& constituent) const
 {
     //check that it is an AK8 jet
-    if(constituent.getType() != AK8JET) return false;
+    if(constituent.getType() != Constituent::AK8JET) return false;
 
     return constituent.p().Pt() > minAK8TopPt_ &&
            constituent.getSoftDropMass() > minAK8TopMass_  && 
@@ -105,5 +105,5 @@ bool TTMConstituentReqs::passDeepAK8TopReqs(const Constituent& constituent) cons
 
 bool TTMConstituentReqs::passAK4ResolvedReqs(const Constituent& constituent, const double minPt) const
 {
-    return constituent.getType() == AK4JET && constituent.p().Pt() > minPt;
+    return constituent.getType() == Constituent::AK4JET && constituent.p().Pt() > minPt;
 }

TopTagger/src/TTMFilterBase.cpp

@@ -15,20 +15,20 @@ bool TTMFilterBase::constituentsAreUsed(const std::vector<const Constituent*>& c
             //First return true if constituent is found (this covers all AK4 and most AK8 jets)
             return true;
         }
-        else if(constituent->getType() == AK8JET)
+        else if(constituent->getType() == Constituent::AK8JET)
         {
             //If the constituent is AK8 we will also check its subjets are not used
-            if(constituent->getSubjets().size() >= 1)
+            if(constituent->getSubjets().size() <= 1)
             {
-                // If this jet has only one subjet, use matching to the overall AK8 jet instaed 
+                // If this jet has only one subjet, use matching to the overall AK8 jet instead 
                 for(const auto& usedConstituent : usedConsts)
                 {
                     if(ROOT::Math::VectorUtil::DeltaR(constituent->p(), usedConstituent->p()) < dRMaxAK8)
                     {
                         //we found a match
                         return true;
                     }
-                }                
+                }
             }
             else
             {
@@ -59,9 +59,9 @@ void TTMFilterBase::markConstituentsUsed(const std::vector<const Constituent *>&
         usedConstituents.insert(constituent);
 
         //If the constituent is an AK8JET, then add AK4JETs matching its subjets as well 
-        if(constituent->getType() == AK8JET)
+        if(constituent->getType() == Constituent::AK8JET)
         {
-            if(constituent->getSubjets().size() >= 1)
+            if(constituent->getSubjets().size() <= 1)
             {
                 //If there is one or fewer subjets, instead match to the overall AK8 jet
                 for(const auto& matchConst : allConstituents) 

TopTagger/src/TTMHEPRequirements.cpp

@@ -74,18 +74,18 @@ void TTMHEPRequirements::run(TopTaggerResults& ttResults)
         }
         else if(doDijet_ && jets.size() == 2) //dijets
         {
-            double m23  = (jets[0]->getType() == AK8JET)?(jets[0]->getSoftDropMass() * jets[0]->getWMassCorr()):(jets[1]->getSoftDropMass() * jets[1]->getWMassCorr());
+            double m23  = (jets[0]->getType() == Constituent::AK8JET)?(jets[0]->getSoftDropMass() * jets[0]->getWMassCorr()):(jets[1]->getSoftDropMass() * jets[1]->getWMassCorr());
             //small hack for legacy tagger
-            if(jets[0]->getType() == AK4JET && jets[1]->getType() == AK4JET) m23 = jets[0]->p().M();
+            if(jets[0]->getType() == Constituent::AK4JET && jets[1]->getType() == Constituent::AK4JET) m23 = jets[0]->p().M();
 
             double m123 = topCand.p().M();
-            if(jets[0]->getType() == AK8JET)
+            if(jets[0]->getType() == Constituent::AK8JET)
             {
                 TLorentzVector psudoVec;
                 psudoVec.SetPtEtaPhiM(jets[0]->p().Pt(), jets[0]->p().Eta(), jets[0]->p().Phi(), jets[0]->getSoftDropMass() * jets[0]->getWMassCorr());
                 m123 = (psudoVec + jets[1]->p()).M();
             }
-            else if(jets[1]->getType() == AK8JET)
+            else if(jets[1]->getType() == Constituent::AK8JET)
             {
                 TLorentzVector psudoVec;
                 psudoVec.SetPtEtaPhiM(jets[1]->p().Pt(), jets[1]->p().Eta(), jets[1]->p().Phi(), jets[1]->getSoftDropMass() * jets[1]->getWMassCorr());
@@ -98,7 +98,7 @@ void TTMHEPRequirements::run(TopTaggerResults& ttResults)
         else if(doMonojet_ && jets.size() == 1) //monojets
         {
             //Monojets just get passed if they are active and AK4
-            if(jets[0]->getType() == AK4JET)
+            if(jets[0]->getType() == Constituent::AK4JET)
             {
                 passHEPRequirments = true;
             }

TopTagger/src/TTMLazyClusterAlgo.cpp

@@ -32,7 +32,7 @@ void TTMLazyClusterAlgo::run(TopTaggerResults& ttResults)
 
     for(unsigned int i = 0; i < constituents.size(); ++i)
     {
-        if(constituents[i].p().Pt() < minJetPt_ || constituents[i].getType() != AK4JET) continue;
+        if(constituents[i].p().Pt() < minJetPt_ || constituents[i].getType() != Constituent::AK4JET) continue;
 
         //singlet tops
         if(doMonojet_)
@@ -54,7 +54,7 @@ void TTMLazyClusterAlgo::run(TopTaggerResults& ttResults)
                 for(unsigned int j = 0; j < constituents.size(); ++j)
                 {
                     if(i == j) continue;
-                    if(constituents[j].p().Pt() < minJetPt_ || constituents[j].getType() != AK4JET) continue;
+                    if(constituents[j].p().Pt() < minJetPt_ || constituents[j].getType() != Constituent::AK4JET) continue;
 
                     TopObject topCand({&constituents[i], &constituents[j]});
 
@@ -71,11 +71,11 @@ void TTMLazyClusterAlgo::run(TopTaggerResults& ttResults)
         {
             for(unsigned int j = 0; j < i; ++j)
             {
-                if(constituents[j].p().Pt() < minJetPt_ || constituents[j].getType() != AK4JET) continue;
+                if(constituents[j].p().Pt() < minJetPt_ || constituents[j].getType() != Constituent::AK4JET) continue;
 
                 for(unsigned int k = 0; k < j; ++k)
                 {
-                    if(constituents[k].p().Pt() < minJetPt_ || constituents[k].getType() != AK4JET) continue;
+                    if(constituents[k].p().Pt() < minJetPt_ || constituents[k].getType() != Constituent::AK4JET) continue;
 
                     TopObject topCand({&constituents[k], &constituents[j], &constituents[i]});
 

TopTagger/src/TTMNanoAODClusterAlgo.cpp

@@ -59,7 +59,7 @@ void TTMNanoAODClusterAlgo::run(TopTaggerResults& ttResults)
         if(doTrijet_)
         {
             //We assume all clustering was done before, we only need to check the type
-            if(constituents[i].getType() == RESOLVEDTOPCAND)
+            if(constituents[i].getType() == Constituent::RESOLVEDTOPCAND)
             {
                 //Fill the resolved top 
 
@@ -75,7 +75,7 @@ void TTMNanoAODClusterAlgo::run(TopTaggerResults& ttResults)
                     std::vector<const Constituent*> resolvedTopConstituents;
                     for(const auto& constituent : constituents)
                     {
-                        if(constituent.getType() == AK4JET)
+                        if(constituent.getType() == Constituent::AK4JET)
                         {
                             int constIndex = constituent.getIndex();
                             if(jetIndex1 == constIndex || jetIndex2 == constIndex || jetIndex3 == constIndex)
@@ -85,6 +85,9 @@ void TTMNanoAODClusterAlgo::run(TopTaggerResults& ttResults)
                         }
                     }
 
+                    //Make sure that this resolved top has exactly three constituents 
+                    if(resolvedTopConstituents.size() != 3) continue;
+
                     TopObject topCand(resolvedTopConstituents, TopObject::RESOLVED_TOP);
                     topCand.setDiscriminator(constituents[i].getTopDisc());
 

TopTagger/src/TTMOverlapResolution.cpp

@@ -39,7 +39,7 @@ void TTMOverlapResolution::getParameters(const cfg::CfgDocument* cfgDoc, const s
                 m1 = t1->p().M();
                 break;
             case 2:
-                if(constVec1[0]->getType() == AK8JET)
+                if(constVec1[0]->getType() == Constituent::AK8JET)
                 {
                     TLorentzVector psudoVec;
                     psudoVec.SetPtEtaPhiM(constVec1[0]->p().Pt(), constVec1[0]->p().Eta(), constVec1[0]->p().Phi(), constVec1[0]->getSoftDropMass() * constVec1[0]->getWMassCorr());
@@ -64,7 +64,7 @@ void TTMOverlapResolution::getParameters(const cfg::CfgDocument* cfgDoc, const s
                 m2 = t2->p().M();
                 break;
             case 2:
-                if(constVec2[0]->getType() == AK8JET)
+                if(constVec2[0]->getType() == Constituent::AK8JET)
                 {
                     TLorentzVector psudoVec;
                     psudoVec.SetPtEtaPhiM(constVec2[0]->p().Pt(), constVec2[0]->p().Eta(), constVec2[0]->p().Phi(), constVec2[0]->getSoftDropMass() * constVec2[0]->getWMassCorr());

TopTagger/src/TTMRemainingSystem.cpp

@@ -47,13 +47,13 @@ void TTMRemainingSystem::run(TopTaggerResults& ttResults)
         if(usedJets.count(&jet)) continue;
 
         //Take first jet not in a top (the constituents are pt ordered)
-        if(seed == nullptr && jet.getType() == AK4JET)
+        if(seed == nullptr && jet.getType() == Constituent::AK4JET)
         {
             seed = &jet;
         }
 
         //switch to a W-jet if one is found
-        if(allowW_ && seed->getType() != AK8JET && passAK8WReqs(jet))
+        if(allowW_ && seed->getType() != Constituent::AK8JET && passAK8WReqs(jet))
         {
             seed = &jet;
         }
@@ -94,7 +94,7 @@ void TTMRemainingSystem::run(TopTaggerResults& ttResults)
             {
                 //if seed is AK4 jet, prefer AK8W over another AK4 jet
                 //else if seed is AK8, prefer AK4 matching jet
-                if(allowW_ && seed->getType() == AK4JET && passAK8WReqs(jet))
+                if(allowW_ && seed->getType() == Constituent::AK4JET && passAK8WReqs(jet))
                 {
                     if(!ak8wFound || dR < mindR)
                     {
@@ -105,7 +105,7 @@ void TTMRemainingSystem::run(TopTaggerResults& ttResults)
                 }
 
                 //if no AK8 jet is found use AK4 jet
-                if(!ak8wFound && jet.getType() == AK4JET && dR < mindR)
+                if(!ak8wFound && jet.getType() == Constituent::AK4JET && dR < mindR)
                 {
                     mindR = dR;
                     secondJet = &jet;

TopTagger/src/TopTaggerInterface.c

@@ -142,7 +142,7 @@ static int TopTaggerInterface_makeAK4Const(
             }
 
             //filter out jet if it is matched to a lepton, or if it has pt < 20 GeV
-            filterVec[iJet] = !isLep && jetsLV[iJet].Pt() > 20.0;
+            filterVec[iJet] = !isLep && jetsLV[iJet].Pt() >= 19.9; //better to underclean just a bit because of nanoAOD rounding
         }
     }
 
@@ -261,7 +261,7 @@ static int TopTaggerInterface_makeAK8Const(
     auto& jetWDisc = tempFloatBuffers.back();
 
     //Create the AK8 constituent helper
-    ak8ConstInputs.reset(new ttUtility::ConstAK8Inputs<Float_t, ttPython::Py_buffer_wrapper<Float_t>, ttPython::Py_buffer_wrapper<TLorentzVector>>(jetsLV, jetTopDisc, jetWDisc, jetSDMass, subjetsLV));
+    ak8ConstInputs.reset(new ttUtility::ConstAK8Inputs<Float_t, ttPython::Py_buffer_wrapper<Float_t>, ttPython::Py_buffer_wrapper<TLorentzVector>>(jetsLV, jetTopDisc, jetWDisc, jetSDMass, vecSubjetsLV));
 
     return 0;
 }