aglevine · aglevine · Sep 30, 2015 · Sep 30, 2015 · Sep 30, 2015
diff --git a/DataFormats/interface/PATFinalState.h b/DataFormats/interface/PATFinalState.h
@@ -166,6 +166,8 @@ class PATFinalState : public pat::PATObject<reco::LeafCandidate> {
     double ht(const std::string& sysTags) const;
     double ht() const;
 
+    double jetHt(const std::string& jetCuts) const;
+
     /// Compute the pZeta variable using the ith and jth legs as
     /// the "visible" objects.
     double pZeta(int i=0, int j=1) const;
@@ -194,9 +196,17 @@ class PATFinalState : public pat::PATObject<reco::LeafCandidate> {
     // and the same thing for a general Lorentz vector. Doesn't check sign.
     double zCompatibility(const PATFinalState::LorentzVector& p4) const;
 
-    // Try this method to see if I can get the parser to work when I reorder my candidates
-    double zCompatibilityFSR(int i, int j, std::string fsrLabel) const;
-
+    // closest Z value
+    double closestZ(int i, const std::string& filter, std::vector<const reco::Candidate*> legs) const;
+    double closestZElectron(int i, const std::string& filter) const;
+    double closestZMuon(int i, const std::string& filter) const;
+    double closestZTau(int i, const std::string& filter) const;
+
+    // smallest invariant mass
+    double smallestMll(int i, const std::string& filter, std::vector<const reco::Candidate*> legs) const;
+    double smallestMee(int i, const std::string& filter) const;
+    double smallestMmm(int i, const std::string& filter) const;
+    double smallestMtt(int i, const std::string& filter) const;
     /// Get the VBF selection variables.  The jet cuts are applied to the veto
     /// jets using dR of 0.3 away from the members.
     VBFVariables vbfVariables(const std::string& jetCuts) const;
@@ -255,22 +265,6 @@ class PATFinalState : public pat::PATObject<reco::LeafCandidate> {
 //    double massUsingSuperCluster(int electronIndex, int j,
 //        int x=-1, int y=-1, int z=-1) const;
 
-    /// Build a subcandidate
-    PATFinalStateProxy subcand(int i, int j,
-        int x=-1, int y=-1, int z=-1) const;
-
-    /// Build a subcandidate w/ fsr
-    PATFinalStateProxy subcandfsr( int i, int j, const std::string& fsrLabel="" ) const;
-
-    /// quad candidate p4 w/ fsr
-    LorentzVector p4fsr(const std::string& fsrLabel="") const;
-
-    /// Returns the index of this lepton's Z partner in 4l ordering
-    const inline size_t get4LPartner(size_t i) const
-    {
-      return i + (i%2 ? -1 : 1);
-    }
-
     /// Build a subcand using a tag string
     PATFinalStateProxy subcand(const std::string& tags) const;
 
@@ -279,9 +273,6 @@ class PATFinalState : public pat::PATObject<reco::LeafCandidate> {
         const std::string& tags, const std::string& extras,
         const std::string& filter="") const;
 
-    // Get the FSR candidate that moves the invariant mass of the lepton pair closest to nominal Z mass
-    const reco::CandidatePtr bestFSROfZ(int i, int j, const std::string& fsrLabel) const;
-
     // Abstract interface to derived classes.  Can return null stuff.
     virtual const reco::Candidate* daughterUnsafe(size_t i) const = 0;
     virtual const reco::CandidatePtr daughterPtrUnsafe(size_t i) const = 0;

diff --git a/DataFormats/interface/PATQuadFinalStateT.h b/DataFormats/interface/PATQuadFinalStateT.h
@@ -1,7 +1,14 @@
 #ifndef FinalStateAnalysis_DataFormats_PATQuadFinalStateT_h
 #define FinalStateAnalysis_DataFormats_PATQuadFinalStateT_h
 
++class PATFinalStateProxy;
++class PATMultiCandFinalState;
++
 #include "FinalStateAnalysis/DataFormats/interface/PATFinalState.h"
+#include "FinalStateAnalysis/DataFormats/interface/PATFinalStateProxy.h"
+#include "FinalStateAnalysis/DataFormats/interface/PATFinalStateEventFwd.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
 
 template<class T1, class T2, class T3, class T4>
 class PATQuadFinalStateT : public PATFinalState {
@@ -87,6 +94,201 @@ class PATQuadFinalStateT : public PATFinalState {
         << ") is null!" << std::endl;
     }
 
+    /// Build a subcandidate w/ fsr
+    PATFinalStateProxy subcandfsr( int i, int j, const std::string& fsrLabel="" ) const
+    {
+      std::vector<reco::CandidatePtr> output;
+      output.push_back( daughterPtr(i) );
+      output.push_back( daughterPtr(j) );
+
+      const reco::CandidatePtr fsrPho = bestFSROfZ(i, j, fsrLabel);
+      if(fsrPho.isNonnull() && fsrPho.isAvailable())
+	output.push_back(fsrPho);
+      return PATFinalStateProxy(new PATMultiCandFinalState(output, evt()));
+    }
+
+    /// quad candidate p4 w/ fsr
+    LorentzVector p4fsr(const std::string& fsrLabel="") const
+    {
+      // start with the 4-momentum of the first Z and add to it
+      PATFinalState::LorentzVector p = subcandfsr(0, 1, fsrLabel)->p4();
+
+      for(unsigned int i = 2; i+1 < numberOfDaughters(); i += 2)
+	{
+	  p += subcandfsr(i, i+1, fsrLabel)->p4();
+	}
+
+      return p;
+    }
+
+    /// Returns the index of this lepton's Z partner in 4l ordering
+    const inline size_t get4LPartner(size_t i) const
+    {
+      return i + (i%2 ? -1 : 1);
+    }
+
+    // Get the FSR candidate that moves the invariant mass of the lepton pair closest to nominal Z mass
+    const reco::CandidatePtr bestFSROfZ(int i, int j, const std::string& fsrLabel) const
+    {
+      bool iIsElectron = false;
+      bool jIsElectron = false;
+
+      edm::Ptr<pat::Electron> ei = daughterAsElectron(i);
+      edm::Ptr<pat::Muon> mi = daughterAsMuon(i);
+      if(daughter(i)->isElectron() && ei.isNonnull() && ei.isAvailable())
+	{
+	  iIsElectron = true;
+	}
+      else 
+	{
+	  if(!(daughter(i)->isMuon() && mi.isNonnull() && mi.isAvailable()))
+	    return reco::CandidatePtr();
+	}
+
+      edm::Ptr<pat::Electron> ej = daughterAsElectron(j);
+      edm::Ptr<pat::Muon> mj = daughterAsMuon(j);
+      if(daughter(j)->isElectron() && ej.isNonnull() && ej.isAvailable())
+	{
+	  jIsElectron = true;
+	}
+      else 
+	{
+	  if(!(daughter(j)->isMuon() && mj.isNonnull() && mj.isAvailable()))
+	    return reco::CandidatePtr();//edm::Ptr<reco::Candidate>(new const reco::Candidate());
+	}
+
+      int leptonOfBest; // index of daughter that has best FSR cand 
+      int bestFSR = -1; // index of best photon as userCand
+      double bestFSRPt = -1; // Pt of best photon
+      double bestFSRDeltaR = 1000; // delta R between best FSR and nearest lepton
+      float noFSRDist = zCompatibility(i,j); // must be better Z candidate than no FSR
+      if(noFSRDist == -1000) // Same sign leptons
+	return reco::CandidatePtr();
+      PATFinalState::LorentzVector p4NoFSR = daughter(i)->p4() + daughter(j)->p4();
+
+      if((iIsElectron?ei->hasUserInt("n"+fsrLabel):mi->hasUserInt("n"+fsrLabel)))
+	{
+	  for(int ind = 0; ind < (iIsElectron?ei->userInt("n"+fsrLabel):mi->userInt("n"+fsrLabel)); ++ind)
+	    {
+	      PATFinalState::LorentzVector fsrCandP4 = daughterUserCandP4(i, fsrLabel+std::to_string(ind));
+	      PATFinalState::LorentzVector zCandP4 = p4NoFSR + fsrCandP4;
+	      if(zCandP4.mass() < 4. || zCandP4.mass() > 100.) // overall mass cut
+		continue;
+	      if(zCompatibility(zCandP4) > noFSRDist) // Must bring us closer to on-shell Z
+		continue;
+	      // If any FSR candidate has pt > 4, pick the highest pt candidate. 
+	      // Otherwise, pick the one with the smallest deltaR to its lepton.
+	      if(bestFSRPt > 4. || fsrCandP4.pt() > 4.)
+		{
+		  if(fsrCandP4.pt() < bestFSRPt)
+		    continue;
+		}
+	      else if(reco::deltaR(daughter(i)->p4(), fsrCandP4) > bestFSRDeltaR)
+		continue;
+
+	      // This one looks like the best for now
+	      leptonOfBest = i;
+	      bestFSR = ind;
+	      bestFSRPt = fsrCandP4.pt();
+	      bestFSRDeltaR = reco::deltaR(daughter(i)->p4(), fsrCandP4);
+	    }
+	}
+      if((jIsElectron?ej->hasUserInt("n"+fsrLabel):mj->hasUserInt("n"+fsrLabel)))
+	{
+	  //      std::cout << "Found an embedded candidate in event " << evt()->evtId().event() << std::endl;
+	  for(int ind = 0; ind < (jIsElectron?ej->userInt("n"+fsrLabel):mj->userInt("n"+fsrLabel)); ++ind)
+	    {
+	      PATFinalState::LorentzVector fsrCandP4 = daughterUserCandP4(j, fsrLabel+std::to_string(ind));
+	      PATFinalState::LorentzVector zCandP4 = p4NoFSR + fsrCandP4;
+	      if(zCandP4.mass() < 4. || zCandP4.mass() > 100.) // overall mass cut
+		continue;
+	      if(zCompatibility(zCandP4) > noFSRDist) // Must bring us closer to on-shell Z
+		continue;
+	      // If any FSR candidate has pt > 4, pick the highest pt candidate. 
+	      // Otherwise, pick the one with the smallest deltaR to its lepton.
+	      if(bestFSRPt > 4. || fsrCandP4.pt() > 4.)
+		{
+		  if(fsrCandP4.pt() < bestFSRPt)
+		    continue;
+		}
+	      else if(reco::deltaR(daughter(j)->p4(), fsrCandP4) > bestFSRDeltaR)
+		continue;
+
+	      // This one looks like the best for now
+	      leptonOfBest = j;
+	      bestFSR = ind;
+	      bestFSRPt = fsrCandP4.pt();
+	      bestFSRDeltaR = reco::deltaR(daughter(j)->p4(), fsrCandP4);
+	    }
+	}
+      if(bestFSR != -1)
+	{
+	  //  std::cout << "Accepted FSR cand, event " << evt()->evtId().event() << std::endl;
+	  return daughterUserCand(leptonOfBest, fsrLabel+std::to_string(bestFSR));
+	}
+      //  std::cout << "Rejected FSR cand, event " << evt()->evtId().event() << std::endl;
+      return reco::CandidatePtr();//edm::Ptr<reco::Candidate>(new const reco::Candidate());
+    }
+
+    double zCompatibilityFSR(int i, int j, const std::string fsrLabel) const
+    {
+      PATFinalStateProxy z = subcandfsr(i, j, fsrLabel);
+      return zCompatibility(z);
+    }
+
+
+    float allFSRIsoContribution(const size_t i, const std::string& label,
+				const float dRMin=0.01, const float dRMax=0.4) const
+    {
+      float isoContrib = 0;
+
+      for(size_t l1 = 0; l1 < numberOfDaughters(); l1 += 2)
+	{
+	  reco::CandidatePtr fsr = bestFSROfZ(l1, l1+1, label);
+	  if(!(fsr.isNonnull() && fsr.isAvailable()))
+	    continue;
+
+	  float dR = reco::deltaR(daughter(i)->p4(), fsr->p4());
+
+	  if(dR > dRMin && dR < dRMax)
+	    isoContrib += fsr->pt();
+	}
+
+      return isoContrib;
+    }
+
+
+    /// Build a subcandidate w/ fsr of the main Zs FSR photons
+    PATFinalStateProxy subcandPrimaryFSR( size_t i, size_t j, const std::string& fsrLabel="" ) const
+    {
+      if(j == get4LPartner(i))
+	return subcandfsr(i,j,fsrLabel);
+
+      std::vector<reco::CandidatePtr> output;
+      output.push_back( daughterPtr(i) );
+      output.push_back( daughterPtr(j) );
+
+      size_t i2 = get4LPartner(i);
+      size_t j2 = get4LPartner(j);
+
+      const reco::CandidatePtr fsr1 = bestFSROfZ(i, i2, fsrLabel);
+      const reco::CandidatePtr fsr2 = bestFSROfZ(j, j2, fsrLabel);
+      if(fsr1.isNonnull() && fsr1.isAvailable())
+	{
+	  // photon is always matched with closer lepton
+	  if(reco::deltaR(daughter(i)->p4(), fsr1->p4()) < reco::deltaR(daughter(i2)->p4(), fsr1->p4()))
+	    output.push_back(fsr1);
+	}
+      if(fsr2.isNonnull() && fsr2.isAvailable())
+	{
+	  if(reco::deltaR(daughter(j)->p4(), fsr2->p4()) < reco::deltaR(daughter(j2)->p4(), fsr2->p4()))
+	    output.push_back(fsr2);
+	}
+      return PATFinalStateProxy(new PATMultiCandFinalState(output, evt()));
+    }
+
+
+
   private:
     edm::Ptr<T1> p1_;
     edm::Ptr<T2> p2_;