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/ISmearerFactory.h"
#include "AtlfastAlgs/GlobalEventData.h" 

//#include "AtlfastAlgs/CorrectionFactor.h"

#include "AtlfastAlgs/MuonAcceptor.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/IsFromHardScatter.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"

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

namespace Atlfast {
  using std::abs;
  using std::vector;
  using std::sort;
  //--------------------------------
  // Constructor 
  //--------------------------------
  
  DefaultReconstructedParticleMaker::DefaultReconstructedParticleMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ),
      m_ncutter(0),
      m_smearer(0),
      m_acceptor(0),
      m_lnkReconstructedParticle(0),
      m_tesIO(0),
      m_log( messageService(), name )
  {
    
    // Set the parameter defaults.
    m_particleType      = 11;
    m_mcPtMin           = 0.0*GeV;
    m_mcEtaMax          = 100.0;
    m_PtMin             = 5.0*GeV;
    m_EtaMax            = 2.5;
    m_doSmearing        = true;
    m_muSmearKey        = 1;
    m_outputLocation    = "/Event/AtlfastReconstructedParticle" ;
    m_muonResFile       = "atlfastDatafiles/MuonResolutionTable.xml" ;
    m_smearParamSchema  = 1;
    m_applyEfficiencies = false;
    
    // Declare the parameters 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( "OutputLocation",   m_outputLocation ) ;
    declareProperty( "MuonResFile",      m_muonResFile ) ;
    declareProperty( "MuonSmearKey",     m_muSmearKey ); 
    
    //cct: declare the smearing parameters and schema.
    // It is possible these JOs will not be found
    // (say, if someone uses AtlfastKStandardOptions.py).
    // In this case, the smearer would overwrite its constructor
    // defaults with -9999 and an STL vector of size 0, which will mess up
    // the smearing. Therefore if nothing is found, we stop athena running.
    //
    // I had originally thought it would be enough to just omit calling
    // the setSmearParamArray and/or setSmearParamSchema methods if no JOs
    // are found, but I decided it is not enough. This is because if, in the future,
    // the defaults change (in AtlfastStandardOptions, for example) then someone
    // runs using AtlfastKStandardOptions, the code picks up the defaults from the
    // constructor in the smearer. Unless someone has maintained the default values
    // in the smearer .h file, the defaults used will be the old, wrong ones.
    // in short, it is hard to maintain!
    //
    // The only assumption we make by having smearParamArray and smearParamSchema
    // in the jobOptions and failing if they are not fould is that if the defaults change,
    // they should be changed in ALL jobOptions.
    // Currently, this is:
    //  share/Atlfast_CBNT.py
    //  share/AtlfastIKinematicDumperStandardOptions.py
    //  share/AtlfastKStandardOptions.py
    //  share/AtlfastStandardOptions.py
    //
    //NB: the reason for all this is that the smearer is NOT an algorithm,
    //so can not have Gaudi Properties.
    //DefaultReconstructedParticleMaker is an algorithm so can have properties.
    //Defaults must be set before properties are declared.
    //Different particle smearers have different defaults, so we can not set these
    //to any sensible values before the properties are loaded. We can only set defaults
    //inside the smearer and these are then overridden by the 'set' method to potentially
    //give junk physics...
    declareProperty( "SmearParamArray", m_smearParamArray );
    declareProperty( "SmearParamSchema",m_smearParamSchema );

    // Option to apply efficiencies at particle creation stage
    declareProperty( "ApplyEfficiencies", m_applyEfficiencies );

  }
  
  //--------------------
  // Destructor
  //--------------------
  
  DefaultReconstructedParticleMaker::~DefaultReconstructedParticleMaker() {

    m_log << MSG::INFO << "Destructor Called" << endreq;
    
    if (m_smearer) {
      m_log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_smearer;
    }
    if (m_acceptor) {
      m_log << MSG::INFO << "Deleting acceptor" << endreq;
      delete m_acceptor;
    }
    if (m_tesIO) {
      m_log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_tesIO;
    }
    if (m_ncutter) {
      m_log << MSG::INFO << "Deleting smearer" << endreq;
      delete m_ncutter;
    }

    if (m_lnkReconstructedParticle) delete m_lnkReconstructedParticle;

  }
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode DefaultReconstructedParticleMaker::initialize()
  {
    
    // 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();
    Selector* fhscSelector = new HepMC_helper::IsFromHardScatter();

    vector<HepMC_helper::IMCselector*> selectors;
    selectors.push_back(fhscSelector);
    selectors.push_back(fstaSelector);
    selectors.push_back(typeSelector);
    selectors.push_back(kineSelector);
    
    m_ncutter = new HepMC_helper::NCutter(selectors);
    
    delete typeSelector;
    delete kineSelector;
    delete fstaSelector;
    delete fhscSelector;
    //=======================================================
    
    
    

    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    int lumi =       ged->lumi();
    int randSeed =   ged->randSeed();

    m_tesIO = new TesIO(ged->mcLocation(), ged->justHardScatter());
    
    //=====================
    if (m_doSmearing){
      m_smearer = ISmearerFactory::create( m_particleType,
                                           randSeed,
                                           lumi,
                                           m_muSmearKey,
                                           m_muonResFile,
                                           m_log,
                                           m_smearParamArray,
                                           m_smearParamSchema );
    } else {
      m_smearer = NULL;
    }
    
    //=====================
    if (m_applyEfficiencies) this->getAcceptor();
    else m_acceptor = NULL;
    
    HeaderPrinter hp("Atlfast ReconstructedParticle Maker:", m_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             ", ged->mcLocation());
    hp.add("Output Location         ", m_outputLocation);
    hp.add("Muon Resolution File    ", m_muonResFile);
    hp.print();
    
    //stop execution if no smearparamschema or smearParamArray are found.
    if(m_smearParamArray.size() == 0){
      m_log << MSG::WARNING <<"No smearing parameters found in ATLFAST jobOptions..."<< endreq;
    }else{
      m_log << MSG::INFO  << "Smearing Array m_smearParamArray set to: "<< m_smearParamArray << endreq;
    }
    if(m_smearParamSchema == -9999){
      m_log << MSG::FATAL <<"No smearing schema found in ATLFAST jobOptions..."<< endreq;
      return StatusCode::FAILURE;
    }else{
      m_log << MSG::INFO  << "Smearing Schema m_smearParamSchema set to: "<< m_smearParamSchema << endreq;
    }
    
    return StatusCode::SUCCESS;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode DefaultReconstructedParticleMaker::finalize()
  {
    
    m_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( )
  {
    
    std::string mess;
    m_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";
    m_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 ;

    // Booleans for overall accept decision and kinematic and 
    // pre-defined efficiency accept decisions separately
    bool accept = true, accept_kinematic = true, accept_efficiency = true;
    
    for(src = my_MC_particles.begin() ; src != my_MC_particles.end() ; ++src){ 
      
      candidate = this->create( *src ); 
      
      // Accept decision based on pre-defined efficiency function
      if (m_applyEfficiencies){ 
        bool instance_check = m_acceptor ? true : false;
        m_log << MSG::DEBUG << "Does m_acceptor point to anything sensible?: " << instance_check << endreq;
        accept_efficiency = ( m_acceptor && m_acceptor->accept( *candidate, m_log ) ) ? true : false;
        m_log << MSG::DEBUG << "m_acceptor && m_acceptor->accept( *candidate ) = " 
            << accept_efficiency << endreq;
      }
      // Accept decision based on kinematic cuts
      accept_kinematic = this->isAcceptable( candidate );
      
      accept = m_applyEfficiencies ? 
        accept_efficiency && accept_kinematic : 
        accept_kinematic;
      
      if ( accept ) {
        m_log << MSG::DEBUG << "Accepted" << endreq;
        myReconstructedParticles->push_back( candidate );     
      }else{
        m_log << MSG::DEBUG << "Rejected" << endreq;
        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";    
    m_log << MSG::DEBUG << mess << endreq;

    m_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 ( const HepMC::GenParticle* src ){
    
    HepLorentzVector evec;
    
    if (m_doSmearing) {
      
      // Ask our Smearer make a smeared 4-vector
      evec = m_smearer->smear( *src );
      
      m_log << MSG::DEBUG 
          << "\n\t"
          << "Unsmeared four-vector: " 
          << src->momentum()
          << "\n\t" 
          << "Smeared four-vector  : " 
          << evec 
          << endreq;
      
    }else{
      
      evec = src->momentum();
      
    }
    
    ReconstructedParticle* candidate = 
      new ReconstructedParticle( src->pdg_id(), evec, src );
    
    m_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 
  (const ReconstructedParticle* candidate ){
    
    // Apply kinematic criteria to the candidate DefaultReconstructedParticle
    
    if( candidate->pT() < m_PtMin ) {        
      m_log << MSG::DEBUG 
          << "Candidate failed pt cut: pt was " 
          << candidate->pT() 
          << " cut was " 
          << m_PtMin 
          << endreq;
      return false ;
    }
    
    if( abs(candidate->eta()) > m_EtaMax ) {
      m_log << MSG::DEBUG 
          << "Candidate failed max eta cut: eta was " 
          << candidate->eta() 
          << " cut was " 
          << m_EtaMax 
          << endreq;
      return false ;
    }
    
    m_log << MSG::DEBUG 
        << "Candidate passed selection cuts "
        << endreq ;
    
    return true ;
  }
  
  //-------------------------------------------
  // getAcceptor
  //------------------------------------------
  void DefaultReconstructedParticleMaker::getAcceptor(){

    m_log << MSG::DEBUG << "Getting an acceptor" << endreq;

    if (m_particleType == 13){
      m_acceptor = new MuonAcceptor(m_muSmearKey,m_log);
    }else{
      m_log << MSG::ERROR << "No Acceptor exists for this particle type!" << endreq;
      m_acceptor = NULL;
    }

  }
  
} // end namespace bracket

---