// ================================================
 // DefaultReconstructedParticleMaker class Implementation
 // ================================================
 //
 // THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
 //
 // Namespace Atlfast
 //
 
 #include "AtlfastAlgs/DefaultReconstructedParticleMaker.h"
 #include "AtlfastAlgs/ReconstructorFactory.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 "TruthHelper/GenIMCselector.h"
 #include "TruthHelper/IsGenType.h"
 #include "TruthHelper/IsGenStable.h" 
 #include "AtlfastUtils/HepMC_helper/MCCuts.h"
 #include "TruthHelper/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_reconstructor(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;
     m_schemaUsesParameters = false;
     m_detEffectsFile    = "";
     
     // 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 );
     declareProperty( "SchemaUsesParameters", m_schemaUsesParameters );
 
     // Option to apply efficiencies at particle creation stage
     declareProperty( "ApplyEfficiencies", m_applyEfficiencies );
 
   }
   
   //--------------------
   // Destructor
   //--------------------
   
   DefaultReconstructedParticleMaker::~DefaultReconstructedParticleMaker() {
 
     m_log << MSG::INFO << "Destructor Called" << endreq;
     
     if (m_reconstructor) {
       m_log << MSG::INFO << "Deleting reconstructor" << endreq;
       delete m_reconstructor;
     }
     if (m_acceptor) {
       m_log << MSG::INFO << "Deleting acceptor" << endreq;
       delete m_acceptor;
     }
     if (m_tesIO) {
       m_log << MSG::INFO << "Deleting tesIO" << endreq;
       delete m_tesIO;
     }
     if (m_ncutter) {
       m_log << MSG::INFO << "Deleting truth helpers" << 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 TruthHelper::GenIMCselector Selector;
 
     Selector* typeSelector = new TruthHelper::IsGenType( m_particleType);
     Selector* kineSelector = new HepMC_helper::MCCuts(m_mcPtMin, m_mcEtaMax);
     Selector* fstaSelector = new TruthHelper::IsGenStable();
 
     vector<TruthHelper::GenIMCselector*> selectors;
     selectors.push_back(fstaSelector);
     selectors.push_back(typeSelector);
     selectors.push_back(kineSelector);
     
      m_ncutter = new TruthHelper::NCutter(selectors);
     
     delete typeSelector;
     delete kineSelector;
     delete fstaSelector;
     //=======================================================
     
     
     
 
     //get the Global Event Data using singleton pattern
     GlobalEventData*   ged = GlobalEventData::Instance();
     int lumi =         ged->lumi();
     int randSeed =     ged->randSeed();
     m_detEffectsFile = ged->detEffectsFileName();
 
     m_tesIO = new TesIO(ged->mcLocation(), ged->justHardScatter());
     
     //=====================
     if (m_doSmearing){
       ReconstructorFactory rf(m_detEffectsFile, m_log);
       m_reconstructor = rf.makeReconstructor( m_particleType,
                                               randSeed,
                                               lumi,
                                               m_muSmearKey,
                                               m_muonResFile,
                                               m_smearParamArray,
                                               m_smearParamSchema );
     } else {
       m_reconstructor = 0;
     }
     
     //=====================
     if (m_applyEfficiencies) this->getAcceptor();
     else m_acceptor = 0;
     
     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_schemaUsesParameters){
        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;
        }
    }else{
         m_log << MSG::INFO << "Schema requires no input parameters"<< 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(src->momentum().px(),
                           src->momentum().py(),
                           src->momentum().pz(),
                           src->momentum().e());
     
     ReconstructedParticle tmpStackParticle(src->pdg_id(),evec,src);
     ReconstructedParticle *candidate;
 
     if(m_doSmearing){
       candidate = new ReconstructedParticle(m_reconstructor->reconstruct(tmpStackParticle));
       m_log << MSG::DEBUG 
             << "\n\t"
               << "Unsmeared four-vector: " 
             << src->momentum()
             << "\n\t" 
             << "Smeared four-vector  : " 
             << candidate->momentum()
             << endreq;
     }
     else{
       candidate = new ReconstructedParticle(tmpStackParticle);
     }
 
     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 = 0;
     }
 
   }
   
 } // end namespace bracket
 
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.