---

JetMaker.cxx

Go to the documentation of this file.

// ================================================
// JetMaker class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast::
//
#include "AtlfastAlgs/JetMaker.h"
#include "AtlfastAlgs/ISmearer.h"
#include "AtlfastAlgs/JetSmearer.h"
#include "AtlfastAlgs/MissingMomentum.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "AtlfastEvent/Jet.h"
#include "AtlfastEvent/ICell.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastEvent/MCparticleCollection.h"

#include "AtlfastUtils/TesIO.h"
#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/HepMC_helper/SelectJetTag.h"
#include "AtlfastUtils/HepMC_helper/SelectTauTag.h"
#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastEvent/SimpleKinematic.h"

#include <cmath> 
#include <algorithm>
#include <numeric>

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "StoreGate/DataHandle.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"



// CLHEP,HepMC
#include "GeneratorObjects/McEventCollection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"


namespace Atlfast{

  // STL helper function object    
    class SmearCell{
    public:
      SmearCell(ISmearer* is): 
        m_smearer(is), m_doSmear((is== 0)? false:true) {}
      HepLorentzVector operator()(const HepLorentzVector& sumSoFar, 
                                  const ICell* c){
        
        HepLorentzVector sum(0., 0., 0., 0.);
        
        if(m_doSmear){
          sum = sumSoFar+(m_smearer->smear( c->momentum() ) );
        }else{
          sum = sumSoFar+(c->momentum());
        }
        
        return sum; 
      }
    private:
      ISmearer* m_smearer;
      bool      m_doSmear;
    };
      
    

//--------------------------------
// Constructors and destructors
//--------------------------------
  using std::abs;
  using std::sort;
  
  JetMaker::JetMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){
    
    // Setting the parameter defaults.
    m_minPT            = 10;         // Min jet trans momentum
    m_maxEta           = 5;          //Max jet eta 
    m_doSmearing       = true;       //smearing is done
    m_rconeb           = 0.400;
    m_rconef           = 0.400;
    
    m_bPtMin           = 5.0;        //min pt for bquark to tag jet
    m_bMaxDeltaR       = 0.2;        //max dR for bquark to tag jet
    
    m_cPtMin           = 5.0;        //min pt for cquark to tag jet
    m_cMaxDeltaR       = 0.2;
    
    m_tauPtMin         = 10.0;
    m_tauMaxDeltaR     = 0.3;
    
    m_etaTagMax        = 2.5;

    m_tauJetPtRatio    = 0.9;
    
    
    // Default paths for entities in the TES
    m_inputLocation             = "/Event/AtlfastClusters";
    m_outputLocation            = "/Event/AtlfastJets";
    m_muonLocation              = "/Event/AtlfastNonIsolatedMuons";
    m_unusedCellLocation        = "/Event/AtlfastUnusedCells";
    m_missingMomentumLocation   = "/Event/AtlfastMissingMomentum";
    m_mcTruthLocation           = "/Event/McEventCollection";
    
    // This is how you declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    declareProperty( "MinimumPT",       m_minPT ) ;
    declareProperty( "MaximumEta",      m_maxEta ) ;
    declareProperty( "DoSmearing",      m_doSmearing ); 
    declareProperty( "RconeB",          m_rconeb );
    declareProperty( "RconeF",          m_rconef );
    declareProperty( "bPtMin",          m_bPtMin );
    declareProperty( "bMaxDeltaR",      m_bMaxDeltaR );
    declareProperty( "cPtMin",          m_cPtMin );
    declareProperty( "cMaxDeltaR",      m_cMaxDeltaR );
    declareProperty( "tauPtMin",        m_tauPtMin );
    declareProperty( "tauMaxDeltaR",    m_tauMaxDeltaR );
    declareProperty( "etaTagMax",       m_etaTagMax );
    declareProperty( "tauJetPtRatio",   m_tauJetPtRatio);
    
    declareProperty( "InputLocation",           m_inputLocation ) ;
    declareProperty( "OutputLocation",          m_outputLocation ) ;
    declareProperty( "UnusedCellLocation",      m_unusedCellLocation ) ;
    declareProperty( "MuonLocation",            m_muonLocation );
    declareProperty( "MissingMomentumLocation", m_missingMomentumLocation );
    declareProperty( "McTruthLocation",         m_mcTruthLocation );
  }
  
  JetMaker::~JetMaker() {
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "Destructor Called" << endreq;
    
    if (m_smearer) {
      log << MSG::INFO << "Deleting jet Smearer" << endreq;
      delete m_smearer;
    }
    if (m_cellSmearer) {
      log << MSG::INFO << "Deleting Cell Smearer" << endreq;
      delete m_cellSmearer;
    }
    if (m_tesIO) {
      delete m_tesIO;
    }
  } 
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode JetMaker::initialize(){
    MsgStream log( messageService(), name() ) ;
    
    log << MSG::DEBUG << "instantiating a JetSmearer" << endreq;
    
    //    m_tesIO = new TesIO(eventDataService());  
    m_tesIO = new TesIO();  
    
    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    int lumi                  = ged->lumi();
    int randSeed              = ged->randSeed() ;
    m_barrelForwardEta = ged->barrelForwardEta();
    
    if(m_doSmearing){
      m_smearer                 = new JetSmearer(randSeed, 
                                                 lumi,
                                                 m_rconeb, 
                                                 m_rconef, 
                                                 m_barrelForwardEta);    
      m_cellSmearer             = new JetSmearer(randSeed, 
                                                 lumi,
                                                 0., 
                                                 0., 
                                                 m_barrelForwardEta);    
    } else {
      m_smearer = NULL;
      m_cellSmearer = NULL;
    }

    HeaderPrinter hp("Atlfast JetMaker:", log);
    
    hp.add("TES Locations:              ");
    hp.add(" Muons from                 ", m_muonLocation); 
    hp.add(" Clusters from              ", m_inputLocation); 
    hp.add(" GeneratorParticles from    ", m_mcTruthLocation);
    hp.add(" Jets to                    ", m_outputLocation);
    hp.add(" unused cell,cluster sum to ", m_missingMomentumLocation);
    hp.add("Smearing on                 ", m_doSmearing);
    hp.add("Cluster min Pt              ", m_minPT);
    hp.add("Cluster max Eta             ", m_maxEta);           
    hp.add("Luminosity                  ", lumi);
    hp.add("Initial random seed         ", randSeed);      
    hp.add("End cap Cone size           ", m_rconef);
    hp.add("Barrel  Cone size           ", m_rconeb);
    hp.add("Eta of start of end cap     ", m_barrelForwardEta);
    hp.add("Parameters for tags         ", "       ");
    hp.add(" Max eta                    ", m_etaTagMax);
    hp.add(" b-quark min pt             ", m_bPtMin);
    hp.add(" b-quark max DeltaR         ", m_bMaxDeltaR);
    hp.add(" c-quark min pt             ", m_cPtMin);
    hp.add(" c-quark max DeltaR         ", m_cMaxDeltaR);
    hp.add(" tau lep min pt             ", m_tauPtMin);
    hp.add(" tau lep max DeltaR         ", m_tauMaxDeltaR);
    hp.add(" tau/Jet Pt Ratio           ", m_tauJetPtRatio);
    hp.print();
    


    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode JetMaker::finalize(){

    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing" << endreq;    
    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  //  This algorithm creates smeared Jets passing cuts. 
  //  It scans the list of Clusters from the Monte Carlo truth. 
  //  If a Cluster of the right flavour is found it creates 
  //  a candidate Jet object by smearing the four vector
  //  of the Cluster and adding non-isolated muons. 
  // 
  //  It then applies kinematic criteria on the properties of the Jet
  //  Those which pass the cuts are kept.
  //  Finally all successful Jets are added to the event store.
  //
  //  
  //
  
  StatusCode JetMaker::execute( ){
    
    MsgStream log( messageService(), name() ) ;    

    std::string mess;
    // make a message logging stream
    

      

    //.............................
    // Make some locals stores which be used to keep pointers to
    // all of the entities from the event store which are needed by
    // this algorithm. These are all local and are typedefed in this class
    
    localClusterCollection      my_Clusters ;
    MuonCollection         my_Muons ;
    //................................
    
    //    if( ! m_tesIO->copy( my_Clusters, m_inputLocation ) ) {
    if( ! m_tesIO->copy<ClusterCollection> ( my_Clusters, m_inputLocation ) ) {
      log << MSG::INFO << "No Clusters in TES " << endreq;
    } else {
      log << MSG::DEBUG << my_Clusters.size()<<" Clusters from TES " << endreq;
    }
    // get non-isolated muons for jet construction
    
    if( ! m_tesIO->copy<ReconstructedParticleCollection>
        ( my_Muons, m_muonLocation ) ) {
      log << MSG::DEBUG << "No muons in TES " << endreq;
    }  

    // ......................
    // Make a container object which will store pointers to all isolated 
    // Jetss which are successfully created.  Since it is going to go 
    // into the event store, it has to be a special Gaudi type of container. 
    // This is typedefined in the entity include file
    // As it is going into the TES it needs to be created on the heap with new.
    
    JetCollection* myJets = new JetCollection ;
    
        //.........................................................
    // From each Cluster create a reconstructed Jet. 
    // If all requirements are satisfied add to the collection
    
    Jet* candidate ;
    localClusterIterator src ;
    HepLorentzVector missingMomentum(0.,0.,0.,0.);
    const int BQUARK(ParticleCodes::BQUARK);
    const int CQUARK(ParticleCodes::CQUARK);
    //Construct an object to select b, c and tau quarks with kine conditions
    const HepMC_helper::SelectJetTag bSelector(BQUARK, m_bPtMin, m_etaTagMax);
    const HepMC_helper::SelectJetTag cSelector(CQUARK, m_cPtMin, m_etaTagMax);
    const HepMC_helper::SelectTauTag tauSelector(m_tauPtMin, m_etaTagMax);

    MCparticleCollection  mcParticles_b ;
    MCparticleCollection  mcParticles_c ;
    MCparticleCollection  mcParticles_tau ;
    //get quarks used in tagging from TES using tesio
    m_tesIO->getMC( mcParticles_b, &bSelector ) ;
    m_tesIO->getMC( mcParticles_c, &cSelector ) ;
    m_tesIO->getMC( mcParticles_tau, &tauSelector ) ;
 
    for( src = my_Clusters.begin() ; src != my_Clusters.end() ; ++src ) 
      { 
        IAssociationManager* iAssocMan = *src;
        ReconstructedParticle rp;
        if (iAssocMan->associations(rp).size() == 0)
          {
            // Make a new jet
            candidate = this->create( *src );

            //Make a temp HepLorentzVector from candidates 4-mom
            //before muon is added in
            HepLorentzVector temp4Vec = candidate->momentum();

            // add isolated muons
            this->addMuons(candidate, my_Muons);

            
            // add to missing momentum vector if not used in a jet.
            // remember, these are smeared momenta if smearing is enabled
            if( this->isAcceptable( candidate ) ) {    

            // if the jet is acceptable flavor tag the jet
              this->tag(candidate,mcParticles_b,mcParticles_c,mcParticles_tau);

              log << MSG::DEBUG << "Pushing back a Jet" << endreq;
              myJets->push_back( candidate ) ; 
              log<<MSG::DEBUG<<"jet "<<*candidate<<endreq;
            }else {
              log << MSG::DEBUG << "Jet fails cuts" << endreq;
              //add in temp4Vec to avoid muon double counting
              missingMomentum += temp4Vec;
              delete candidate ;
            }
          }
      }
    
    log<<MSG::DEBUG<<"Missing momentum: clusters "<<missingMomentum<<endreq;
    missingMomentum += addUnusedCells(log);
    log<<MSG::DEBUG<<"Missing momentum: cells+clstrs "
       <<missingMomentum<<endreq;
    //......................................
    // Register any Jets which have been successfilly created in the 
    // transient event store. 
    TesIoStat stat;
    log<<MSG::DEBUG<<"Storing "<< myJets->size() <<" Jets"<<endreq;
    stat = m_tesIO->store( myJets, m_outputLocation ) ;
    if(!stat){
      log<<MSG::ERROR<<"Could not store jets"<<endreq;
      return stat;
    }


    MissingMomentum* mm = new MissingMomentum(missingMomentum);
    stat = m_tesIO->store(mm,m_missingMomentumLocation ) ;
    std::string message = stat? "Stored missing mom":"Failed to missing mom";
    log<<MSG::DEBUG<<message<<endreq;
  

    return StatusCode::SUCCESS;
  }
  
  
  
  //----------------------------------
  // create()
  //----------------------------------
  
  // Takes a single Cluster and creates a Jet
  // according to the desired criteria
  //
  // Note that we must NEW these jets so they can go to the TES
  // and if we decide that we do not want them, we must DELETE them
  // Once they are put to the TES, they are no longer our responsibility to delete
  
  Jet* JetMaker::create ( Cluster* src )
  {
    Jet* jet;
    MsgStream log( messageService(), name() ) ;
    HepLorentzVector vec = src->momentum();
    
    if (m_doSmearing) 
      {
        // Ask our Smearer make a smeared 4-vector
        vec = m_smearer->smear(vec);      
      }
    
    jet = new Jet(vec, *src);
    log << MSG::DEBUG << "Unsmeared Cluster: " << *src << endreq;
    log << MSG::DEBUG << "Smeared Jet      : " << jet << endreq ;
    
    // return Jet
    return jet;
    
    
  }
  //-------------------------------------------
  // addmuons
  //-------------------------------------------
  
  // adds non-isolated muons to the jets
  
  void JetMaker::addMuons(Jet* jet, MuonCollection& muons )
  {
    
    if (muons.size() == 0)  return;
    
    MsgStream log( messageService(), name() ) ;
    
    MuonIterator first = muons.begin() ;
    MuonIterator last  = muons.end() ;
    
    // Sort by distance from the jet
    sort( first, last, SortAttribute::DeltaR(jet) );
    
    // check that closest muon is within cone  
    // if it is, add muon to jet
    
    float muonJetRdist = m_kinehelp.deltaR( jet, *first ) ;
    
    log << MSG::DEBUG << "distance in R of muon to Jet " << muonJetRdist << endreq;
    
    if (
        (
         abs((*first)->eta()) <= m_barrelForwardEta && 
         muonJetRdist < m_rconeb
         ) 
        ||
        (
         abs((*first)->eta()) >  m_barrelForwardEta && 
         muonJetRdist < m_rconef
         )
        ) 
      {
        
        log << MSG::DEBUG 
            << "Adding muon to jet with momentum "
            << (*first)->momentum() 
            << endreq;
        
        log << MSG::DEBUG << "Old mom: " << *jet << endreq;
        
        jet->setMomentum(jet->momentum() + (*first)->momentum()); 
        
        log << MSG::DEBUG << "New mom: " << *jet << endreq ;
        
        // create association
        IAssociationManager* ia = jet;
        ia->associate(*first);
        
        // remove muon from list
        muons.erase(first);
        
        
      }
    else
      {        
        log << MSG::DEBUG << "Muon not in R cone" << endreq;
      }
    
  }
  
  
  //-------------------------------------------
  // isAcceptable
  //------------------------------------------
  
  // Decides whether a candidate is acceptable to this Maker, i.e. whether
  // it wants to proceed and add it to its output product collection
  
  bool JetMaker::isAcceptable ( Jet* candidate)
  {
    
    MsgStream log( messageService(), name() ) ;
    
    // Apply kinematic criteria to the candidate jet
    
    if( candidate->pT() < m_minPT ) {        
      log << MSG::DEBUG << "Candidate failed pt cut: pt " << candidate->pT() 
          << " min was " << m_minPT << endreq;
      return false ;
    }
    
    if( abs(candidate->eta()) > m_maxEta ) {
      log << MSG::DEBUG << "Candidate failed max eta cut: eta " 
          << candidate->eta() 
          << " max was " << m_maxEta << endreq;
      return false ;
    }
    
    return true ;
  }
  
  //-------------------------------------------
  // tag
  //------------------------------------------
  //
  
  void JetMaker::tag(Jet* jet, 
                     MCparticleCollection& mcParticles_b, 
                     MCparticleCollection& mcParticles_c, 
                     MCparticleCollection& mcParticles_tau)const {
 
    //See if jet can be b-tagged
    //Construct an object to select IKinematics within DeltaR of the current jet
    PartitionCondition::BelowThresholdDeltaR bTagSelector(jet, m_bMaxDeltaR); 
    int btag = getJetTags(mcParticles_b,bTagSelector);
    if(btag != 0 )
      {
        jet->setBTag();
        return;
      }
    //else see if jet can be c-tagged
    //Construct an object to select IKinematics within DeltaR of the current jet
    PartitionCondition::BelowThresholdDeltaR cTagSelector(jet, m_cMaxDeltaR); 
    int ctag = getJetTags( mcParticles_c, cTagSelector);
    if(ctag != 0 )
      { 
        jet->setCTag();
        return;
      }

    //finally see if jet can be tau-tagged
    //Construct an object to select IKinematics within DeltaR of the current jet
    PartitionCondition::BelowThresholdDeltaR tauTagSelector(jet, m_tauMaxDeltaR); 
    int tauTag=getTauTags( mcParticles_tau, tauTagSelector, jet->pT());
    if(tauTag !=0 ) jet->setTauTag(tauTag);
    return; 
  }
  //-------------------------------------------
  // getJetTags
  //------------------------------------------
  int JetMaker::getJetTags(MCparticleCollection&  mcParticles, 
                           const PartitionCondition::BelowThresholdDeltaR& tagSelector)const{
    
    MsgStream log( messageService(), name() ) ;

    MCparticleCollectionCIter src ;
        for( src = mcParticles.begin(); src != mcParticles.end(); ++src ){
          if(tagSelector(SimpleKinematic(*src))) return (*src)->pdg_id();
        }
        return 0;
  }
  //-------------------------------------------
  // getTauTags
  //------------------------------------------
  int JetMaker::getTauTags(
                           MCparticleCollection&  mcParticles,
                           const PartitionCondition::BelowThresholdDeltaR& 
                           tagSelector,
                           const double jetPt)const{
    
    MsgStream log( messageService(), name() ) ;
    MCparticleCollectionCIter src ;
    for( src = mcParticles.begin() ; src != mcParticles.end() ; ++src )
      {
        double tauPt = (*src)->momentum().perp();
        double neuPt=0.;
        // Find Tau children
        HepMC::GenVertex*  endVertex = (*src)->end_vertex();
        if(!endVertex) return 0;

        HepMC::GenVertex::particle_iterator firstChild = 
          endVertex->particles_begin(HepMC::children);
        
        HepMC::GenVertex::particle_iterator endChild = 
          endVertex->particles_end(HepMC::children);
        
        HepMC::GenVertex::particle_iterator thisChild = firstChild; 
        
        // Find Pt of Taus tau-neutrino
        for(; thisChild!=endChild; ++thisChild){
          if(  abs((*thisChild)->pdg_id()) == 16){
            neuPt = (*thisChild)->momentum().perp();
          }
        }
          
          if(tagSelector(SimpleKinematic(*src))      &&
             ((tauPt-neuPt)/jetPt) >= m_tauJetPtRatio) return (*src)->pdg_id();
        
        }
      return 0;
  }   

  HepLorentzVector JetMaker::addUnusedCells(MsgStream& log) const{

    HepLorentzVector temp(0., 0., 0., 0.);
    

    DataHandle<ITwoCptCellCollection> unusedCells;
    if( ! m_tesIO->getDH( unusedCells, m_unusedCellLocation) ){
      log << MSG::INFO << "No unused cells  in TES " << endreq;
      return temp;
    }  
    log << MSG::DEBUG << "No unused cells  in TES " << endreq;
    
    return std::accumulate(unusedCells->begin(), 
                           unusedCells->end(), 
                           temp, 
                           SmearCell(m_cellSmearer)
                           );

  }
    
} // end namespace bracket

---