CMSBoostedTauSeedingAlgorithm.cc

Go to the documentation of this file.

// CMSBoostedTauSeedingAlgorithm Package
//
//----------------------------------------------------------------------
// This file is part of FastJet contrib.
//
// It is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.
//
// It is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this code. If not, see <http://www.gnu.org/licenses/>.

//CG: From CMSSW 8_0_X cut on eta subjets <3. Anyways taus are cut at 2.3 

//----------------------------------------------------------------------

#include "RecoJets/JetAlgorithms/interface/CMSBoostedTauSeedingAlgorithm.h"

FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

namespace contrib{


  // constructor
  CMSBoostedTauSeedingAlgorithm::CMSBoostedTauSeedingAlgorithm(double iminPt,                  
                                   double iminMassDrop, double imaxMassDrop,
                                   double iminY, double imaxY,
                                   double iminDeltaR, double imaxDeltaR,
                                   int maxDepth, 
                                   int verbosity) 
    : ptMin_(iminPt), 
      muMin_(iminMassDrop), muMax_(imaxMassDrop), 
      yMin_(iminY), yMax_(imaxY), 
      dRMin_(iminDeltaR), dRMax_(imaxDeltaR),
      maxDepth_(maxDepth),
      verbosity_(verbosity)
  {}

  // description
  std::string CMSBoostedTauSeedingAlgorithm::description() const {
    std::ostringstream oss;
    oss << "CMSBoostedTauSeedingAlgorithm algorithm";
    return oss.str();
  }


  void CMSBoostedTauSeedingAlgorithm::dumpSubJetStructure(const fastjet::PseudoJet& jet, int depth, int maxDepth, const std::string& depth_and_idx_string) const
  {
    if ( maxDepth != -1 && depth > maxDepth ) return;
    fastjet::PseudoJet subjet1, subjet2; 
    bool hasSubjets = jet.has_parents(subjet1, subjet2);
    if ( !hasSubjets ) return;
    std::string depth_and_idx_string_subjet1 = depth_and_idx_string;
    if ( depth_and_idx_string_subjet1.length() > 0 ) depth_and_idx_string_subjet1.append(".");
    depth_and_idx_string_subjet1.append("0");
    for ( int iSpace = 0; iSpace < depth; ++iSpace ) {
      std::cout << " ";
    }
    std::cout << " jetConstituent #" << depth_and_idx_string_subjet1 << " (depth = " << depth << "): Pt = " << subjet1.pt() << "," 
          << " eta = " << subjet1.eta() << ", phi = " << subjet1.phi() << ", mass = " << subjet1.m() 
          << " (constituents = " << subjet1.constituents().size() << ")" << std::endl;
    dumpSubJetStructure(subjet1, depth + 1, maxDepth, depth_and_idx_string_subjet1);
    std::string depth_and_idx_string_subjet2 = depth_and_idx_string;
    if ( depth_and_idx_string_subjet2.length() > 0 ) depth_and_idx_string_subjet2.append(".");
    depth_and_idx_string_subjet2.append("1");
    for ( int iSpace = 0; iSpace < depth; ++iSpace ) {
      std::cout << " ";
    }
    std::cout << " jetConstituent #" << depth_and_idx_string_subjet2 << " (depth = " << depth << "): Pt = " << subjet2.pt() << "," 
          << " eta = " << subjet2.eta() << ", phi = " << subjet2.phi() << ", mass = " << subjet2.m() 
          << " (constituents = " << subjet2.constituents().size() << ")" << std::endl;
    dumpSubJetStructure(subjet2, depth + 1, maxDepth, depth_and_idx_string_subjet2);
    for ( int iSpace = 0; iSpace < depth; ++iSpace ) {
      std::cout << " ";
    }
    double dR = subjet1.delta_R(subjet2);
    std::cout << " (mass-drop @ " << depth_and_idx_string << " = " << std::max(subjet1.m(), subjet2.m())/jet.m() << ", dR = " << dR << ")" << std::endl;
  }

  std::pair<PseudoJet, PseudoJet> CMSBoostedTauSeedingAlgorithm::findSubjets(const PseudoJet& jet, int depth, bool& subjetsFound) const 
  {
    const float etaMax_ = 3.0; // cut on eta subjets <3. Anyway taus are cut at 2.3 
    if ( verbosity_ >= 2 ) {
      std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
      std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << ", mass = " << jet.m() << std::endl;
    }

    PseudoJet subjet1, subjet2;
    bool hasSubjets = jet.has_parents(subjet1, subjet2);
    if ( hasSubjets && (maxDepth_ == -1 || depth <= maxDepth_) ) {
      // make subjet1 the more massive jet
      if ( subjet1.m2() < subjet2.m2() ) {
    std::swap(subjet1, subjet2);
      }
      double dR = subjet1.delta_R(subjet2);
      double kT = subjet1.kt_distance(subjet2);
      double mu = ( jet.m() > 0. ) ? 
    sqrt(std::max(subjet1.m2(), subjet2.m2())/jet.m2()) : -1.;
      // check if subjets pass selection required for seeding boosted tau reconstruction
      if ( subjet1.pt() > ptMin_ && fabs(subjet1.eta()) < etaMax_  && subjet2.pt() > ptMin_ && fabs(subjet2.eta()) < etaMax_ && dR > dRMin_ && dR < dRMax_ && mu > muMin_ && mu < muMax_ && kT < (yMax_*jet.m2()) && kT > (yMin_*jet.m2()) ) {
    subjetsFound = true;
    return std::make_pair(subjet1, subjet2); 
      } else if ( subjet1.pt() > ptMin_ ) {
    return findSubjets(subjet1, depth + 1, subjetsFound);
      } else if ( subjet2.pt() > ptMin_ ) {
    return findSubjets(subjet2, depth + 1, subjetsFound);
      } 
    }
    subjetsFound = false;
    PseudoJet dummy_subjet1, dummy_subjet2;
    return std::make_pair(dummy_subjet1, dummy_subjet2); 
  }

  PseudoJet CMSBoostedTauSeedingAlgorithm::result(const PseudoJet& jet) const 
  {
    if ( verbosity_ >= 1 ) {
      std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
      std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << ", mass = " << jet.m() << std::endl;
    }

    if ( verbosity_ >= 2 ) {
      dumpSubJetStructure(jet, 0, maxDepth_, "");
    }

    bool subjetsFound = false;
    std::pair<PseudoJet, PseudoJet> subjets = findSubjets(jet, 0, subjetsFound);
    if ( subjetsFound ) {
      // fill structure for returning result
      PseudoJet subjet1 = subjets.first;
      PseudoJet subjet2 = subjets.second;
      if ( verbosity_ >= 1 ) {
    std::cout << "before recombination:" << std::endl;
    std::cout << " subjet #1: Pt = " << subjet1.pt() << ", eta = " << subjet1.eta() << ", phi = " << subjet1.phi() << ", mass = " << subjet1.m() << std::endl;
    std::cout << " subjet #2: Pt = " << subjet2.pt() << ", eta = " << subjet2.eta() << ", phi = " << subjet2.phi() << ", mass = " << subjet2.m() << std::endl;
      }

      const JetDefinition::Recombiner* rec = jet.associated_cluster_sequence()->jet_def().recombiner();
      PseudoJet result_local = join(subjet1, subjet2, *rec);
      if ( verbosity_ >= 1 ) {
    std::cout << "after recombination:" << std::endl;
    std::vector<fastjet::PseudoJet> subjets = result_local.pieces();
    int idx_subjet = 0;
    for ( std::vector<fastjet::PseudoJet>::const_iterator subjet = subjets.begin();
          subjet != subjets.end(); ++subjet ) {
      std::cout << " subjet #" << idx_subjet << ": Pt = " << subjet->pt() << ", eta = " << subjet->eta() << ", phi = " << subjet->phi() << ", mass = " << subjet->m() 
            << " (#constituents = " << subjet->constituents().size() << ")" << std::endl;
      std::vector<fastjet::PseudoJet> constituents = subjet->constituents();
      int idx_constituent = 0;
      for ( std::vector<fastjet::PseudoJet>::const_iterator constituent = constituents.begin();
        constituent != constituents.end(); ++constituent ) {
        if ( constituent->pt() < 1.e-3 ) continue; // CV: skip ghosts
        std::cout << "  constituent #" << idx_constituent << ": Pt = " << constituent->pt() << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << "," 
              << " mass = " << constituent->m() << std::endl;
        ++idx_constituent;
      }
      ++idx_subjet;
    }
      }

      CMSBoostedTauSeedingAlgorithmStructure* s = new CMSBoostedTauSeedingAlgorithmStructure(result_local);
      //s->_original_jet = jet;
      s->_mu = ( jet.m2() > 0. ) ? sqrt(std::max(subjet1.m2(), subjet2.m2())/jet.m2()) : 0.;
      s->_y  = ( jet.m2() > 0. ) ? subjet1.kt_distance(subjet2)/jet.m2() : 0.;
      s->_dR = subjet1.delta_R(subjet2);
      s->_pt = subjet2.pt();
            
      result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));

      return result_local;
    } else {
      // no subjets for seeding boosted tau reconstruction found, return an empty PseudoJet
      if ( verbosity_  >= 1 ) {
    std::cout << "No subjets found." << std::endl;
      }
      return PseudoJet();
    }
  }


} // namespace contrib
 
FASTJET_END_NAMESPACE