---

DefaultReconstructedParticleMaker.cxx

Go to the documentation of this file.

// ================================================
// DefaultReconstructedParticleMaker class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast
//

#include "AtlfastAlgs/DefaultReconstructedParticleMaker.h"
#include "AtlfastAlgs/DefaultSmearer.h"
#include "AtlfastAlgs/ElectronSmearer.h"
#include "AtlfastAlgs/MuonSmearer.h"
#include "AtlfastAlgs/PhotonSmearer.h"
#include "AtlfastAlgs/GlobalEventData.h" 

#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/MCparticleCollection.h"

#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/HepMC_helper/IMCselector.h"
#include "AtlfastUtils/HepMC_helper/SelectType.h"
#include "AtlfastUtils/HepMC_helper/IsFinalState.h"
#include "AtlfastUtils/HepMC_helper/MCCuts.h"
#include "AtlfastUtils/HepMC_helper/NCutter.h"
#include "AtlfastUtils/FunctionObjects.h"

#include <cmath> 

// Generic algorithms
#include <algorithm>

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

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

namespace Atlfast {
  using std::abs;
  using std::vector;
  using std::sort;
  //--------------------------------
  // Constructor 
  //--------------------------------
  
  DefaultReconstructedParticleMaker::DefaultReconstructedParticleMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ){
    
    // Set the parameter defaults.
    m_particleType      = 11;
    m_mcPtMin           = 0.0;
    m_mcEtaMax          = 100.0;
    m_PtMin             = 5.0;
    m_EtaMax            = 2.5;
    m_doSmearing        = true;
    m_muSmearKey        = 1;
    m_MC_eventLocation  = "/Event/McEventCollection" ;
    m_outputLocation    = "/Event/AtlfastReconstructedParticle" ;
    
    // Declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    declareProperty( "ParticleType", m_particleType ) ;
    declareProperty( "mcMinimumPt", m_mcPtMin ) ;
    declareProperty( "mcMaximumEta", m_mcEtaMax ) ;
    declareProperty( "MinimumPt", m_PtMin ) ;
    declareProperty( "MaximumEta", m_EtaMax ) ;
    declareProperty( "DoSmearing",m_doSmearing ); 
    declareProperty( "MC_eventLocation", m_MC_eventLocation ) ;
    declareProperty( "OutputLocation", m_outputLocation ) ;
    declareProperty( "MuonSmearKey",m_muSmearKey ); 
    
  }
  
  //--------------------
  // Destructor
  //--------------------
  
  DefaultReconstructedParticleMaker::~DefaultReconstructedParticleMaker() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "Destructor Called" << endreq;
    
    if (m_smearer) {
      log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_smearer;
    }
    if (m_tesIO) {
      log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_tesIO;
    }
    if (m_ncutter) {
      log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_ncutter;
    }
  }
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode DefaultReconstructedParticleMaker::initialize()
  {
    MsgStream log( messageService(), name() ) ;
    
    // This method is called by Athena once only, after all 
    // properties have been set.
    
    // Set up NCutter to select the required HepMC::GenParticles
    typedef HepMC_helper::IMCselector Selector;

    Selector* typeSelector = new HepMC_helper::SelectType( m_particleType);
    Selector* kineSelector = new HepMC_helper::MCCuts(m_mcPtMin, m_mcEtaMax);
    Selector* fstaSelector = new HepMC_helper::IsFinalState();

    vector<HepMC_helper::IMCselector*> selectors;
    selectors.push_back(fstaSelector);
    selectors.push_back(typeSelector);
    selectors.push_back(kineSelector);
    
    m_ncutter = new  HepMC_helper::NCutter(selectors);
    
    delete typeSelector;
    delete fstaSelector;
    delete kineSelector;
    //=======================================================
    
    
    //    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();
    log << MSG::DEBUG << "got lumi " << endreq;
    int randSeed = ged->randSeed() ;
    if (m_doSmearing) this->getSmearer(lumi, randSeed, log);
    else m_smearer = NULL;
    
    HeaderPrinter hp("Atlfast ReconstructedParticle Maker:", log);
    hp.add("Particle Type           ", m_particleType);
    hp.add("Luminosity              ", lumi);        
    hp.add("Minimum four-vector Pt  ", m_mcPtMin);
    hp.add("Maximum four-vector Eta ", m_mcEtaMax);
    hp.add("Minimum particle    Pt  ", m_PtMin);
    hp.add("Maximum particle    Eta ", m_EtaMax);
    hp.add("Do Smearing             ", m_doSmearing);
    hp.add("Muon  Smearing Flag     ", m_muSmearKey);
    hp.add("Random Number Seed      ", randSeed);
    hp.add("MC location             ", m_MC_eventLocation);
    hp.add("Output Location         ", m_outputLocation);    
    hp.print();
    

      
    
        
    // We have to build an appropriate smearer object here.
    // The code which follows is to be regarded as a temporary fix to
    // get around the fact that Algorithms cannot be specialised by
    // constructor arguments.
    // When tools become avialable we will probably instead just pass 
    // an appropriate tool via job options
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode DefaultReconstructedParticleMaker::finalize()
  {
    
    MsgStream log( messageService(), name() ) ;
    
    log << MSG::INFO << "Finalizing" << endreq;  
    
    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  //  This algorithm creates smeared ReconstructedParticles.
  // 
  //  It takes HepMC::GenParticles from the TES as input. 
  //  Only HepMC::GenParticles satisfying the selection requirements
  //  are used.
  //
  //  A candidate ReconstructedParticle object is made by smearing the 
  //  four vector of the HepMC::GenParticle. 
  // 
  //  It then applies kinematic selection criteria to the properties of the 
  //  ReconstructedParticle. Only those which pass the cuts are kept.
  //
  //  Finally, all surviving ReconstructedParticles are added to the event 
  //  store.
  //
  //
  
  StatusCode DefaultReconstructedParticleMaker::execute( )
  {
    
    MsgStream log( messageService(), name() ) ;
    std::string mess;
    log << MSG::DEBUG << "Execute() " << endreq;
    

    //.............................
    // Obtain the truth HepMC::GenParticles which we want to process.
    // We fill a local collection with pointers to these Particles.
    // The private method which is used applies selection criteria
    // to the HepMC::GenParticles before placing them in the collection
    
    MCparticleCollection  my_MC_particles ;
    MCparticleCollectionCIter src ;
    
    TesIoStat stat = m_tesIO->getMC( my_MC_particles, m_ncutter ) ;
    mess = stat? "Retrieved MC from TES ":"Failed MC retrieve from TES";
    log << MSG::DEBUG << mess << endreq;
    
    // ......................
    // Make a container object which will store pointers to all  
    // ReconstructedParticles which are successfully created.  
    // Since it is going into the TES, it has to be a special Athena type 
    // of container. This is defined in an include file.
    
    ReconstructedParticleCollection* myReconstructedParticles 
      = new ReconstructedParticleCollection ;
    
    
    //.........................................................
    // From each mc truth Particle, create a ReconstructedParticle candidate.
    // If it satisfies requirements add to the output collection
    
    ReconstructedParticle* candidate ;
    
    for(src = my_MC_particles.begin() ; src != my_MC_particles.end() ; ++src){ 
        
      
      candidate = this->create( log, *src ); 
      
      if( this->isAcceptable( log, candidate ) ) {
        myReconstructedParticles->push_back( candidate ) ;     
      }else{
        delete candidate ;
      }
        
    }
    
    
    //.....................................
    // Sort into ascending order of pT
    
    if( myReconstructedParticles->size() > 0 ){
      sort( myReconstructedParticles->begin(),
            myReconstructedParticles->end(),
            SortAttribute::DescendingPT()
            ) ;
    }
    
    
    //......................................
    // Deal with resulting ReconstructedParticles (if any)
    
    stat =m_tesIO->store( myReconstructedParticles, m_outputLocation );
    mess = stat? "Store to TES success":"Store to TES failure";    
    log << MSG::DEBUG << mess << endreq;

    log << MSG::DEBUG << "Ending<-------- " << endreq;
    
    StatusCode sc=StatusCode::SUCCESS;  
    return sc ;
    
  }
  
  
  
  //----------------------------------
  // create()
  //----------------------------------
  
  // Takes a single MC Particle and creates a reconstructed
  // ReconstructedParticle according to the desired criteria
  //
  // Note that we must NEW these particles 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
  
  ReconstructedParticle* 
  DefaultReconstructedParticleMaker::create (
                                             MsgStream& log, 
                                             const HepMC::GenParticle* src ){
    
    ReconstructedParticle* candidate ;
    
    if (m_doSmearing) {
      
      // Ask our Smearer make a smeared 4-vector
      HepLorentzVector evec = m_smearer->smear( src->momentum() );
      
      log << MSG::DEBUG 
          << "\n\t"
          << "Unsmeared four-vector: " 
          << src->momentum()
          << "\n\t" 
          << "Smeared four-vector  : " 
          << evec 
          << endreq;
      
      // Note:It would really make sense to apply the kinematic cuts here. 
      // rather than in isAcceptable method.
      // However we leave this to next design phase as we also need to
      // consider where to put Isolation cuts.
      
      //Make a new ReconstructedParticle with smeared 4-vector
      //RMS Fix unsigned pdg_id
      candidate =  new ReconstructedParticle
        (src->pdg_id() /*m_particleType*/, evec, src ) ;
    }else{
      // Just use the original unsmeared four-momentum
      candidate = new ReconstructedParticle
        ( src->pdg_id() /*m_particleType*/, src->momentum(), src );
    }
    
    
    log << MSG::DEBUG 
        << "Created ReconstructedParticle: \n\t " 
        << candidate 
        << endreq;
    
    return candidate ;
    
  }
  
  
  //-------------------------------------------
  // 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 DefaultReconstructedParticleMaker::isAcceptable 
  (MsgStream& log, const ReconstructedParticle* candidate ){
    
    // Apply kinematic criteria to the candidate DefaultReconstructedParticle
    
    if( candidate->pT() < m_PtMin ) {        
      log << MSG::DEBUG 
          << "Candidate failed pt cut: pt was " 
          << candidate->pT() 
          << " cut was " 
          << m_PtMin 
          << endreq;
      return false ;
    }
    
    if( abs(candidate->eta()) > m_EtaMax ) {
      log << MSG::DEBUG 
          << "Candidate failed max eta cut: eta was " 
          << candidate->eta() 
          << " cut was " 
          << m_EtaMax 
          << endreq;
      return false ;
    }
    
    log << MSG::DEBUG 
        << "Candidate passed selection cuts "
        << endreq ;
    
    return true ;
  }
  
  //-------------------------------------------
  // getSmearer
  //------------------------------------------
  void DefaultReconstructedParticleMaker::getSmearer(int lumi, 
                                                     int randSeed, 
                                                     MsgStream& log){
    
    if (m_particleType == 11) 
      {
        log << MSG::DEBUG << "instantiating an ElectronSmearer" << endreq;
        m_smearer = new ElectronSmearer(randSeed, lumi);
      } 
    else if (m_particleType == 13) 
      {
        log << MSG::DEBUG << "instantiating a MuonSmearer" << endreq;
        m_smearer = new MuonSmearer(randSeed, lumi, m_muSmearKey);
      } 
    else if (m_particleType == 22) 
      {
        log << MSG::DEBUG << "instantiating a PhotonSmearer" << endreq;
        m_smearer = new PhotonSmearer(randSeed, lumi);
      } 
    else 
      {
        log << MSG::ERROR 
            << "No specific smearer known for particle type " 
            << m_particleType 
            << " using a default which is probably not what you need!" 
            << endreq;
        m_smearer = new DefaultSmearer(randSeed);
      }
  }
} // end namespace bracket

---