---

TrackMaker.cxx

Go to the documentation of this file.

// ================================================
// TrackMaker class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast::
//
#include "AtlfastAlgs/TrackMaker.h"
#include "AtlfastAlgs/TrackSmearer.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "AtlfastUtils/HepMC_helper/IMCselector.h"
#include "AtlfastUtils/HepMC_helper/IsCharged.h"
#include "AtlfastUtils/HepMC_helper/IsFinalState.h"
#include "AtlfastUtils/HepMC_helper/NCutter.h"
#include "AtlfastUtils/HepMC_helper/MCCuts.h"

#include "CLHEP/Matrix/Matrix.h"

#include <cmath> 
#include <algorithm>
#include <fstream>

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

// Generator includes
//#include "GeneratorObjects/McEventCollection.h"

//--------------------------------
// Constructors and destructors
//--------------------------------

using std::abs;
using ::HepMatrix;

namespace Atlfast {
  
  TrackMaker::TrackMaker ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ), m_smearer(NULL){
    
    // Setting the parameter defaults - these can be changed by user via
    // jobOptions
    m_mcPtMin      = 0.5;
    m_mcEtaMax     = 2.5;
    m_vlMax        = 40.0;
    m_vtMax        = 3.0;
    m_doSmearing   = true;
    m_bField       = 2.0;
    m_muConfig     = "000";
    
    // Default paths for entities in the TES
    m_MC_eventLocation  = "/Event/McEventCollection";
    m_outputLocation    = "/Event/AtlfastTracks";
    
    // Declare the paramemters to Athena so that
    // they can be over-written via the job options file
    declareProperty( "McPtMinimum",       m_mcPtMin ) ;
    declareProperty( "McEtaMaximum",      m_mcEtaMax ) ;
    declareProperty( "vlMaximum",         m_vlMax ) ;
    declareProperty( "vtMaximum",         m_vtMax ) ;
    declareProperty( "DoSmearing",        m_doSmearing ); 
    declareProperty( "Bfield",            m_bField );
    declareProperty( "MC_eventLocation",  m_MC_eventLocation ) ;
    declareProperty( "OutputLocation",    m_outputLocation ) ;
    declareProperty( "MuonConfiguration", m_muConfig );
  }
  
  //----------------------------
  // Destructor
  //----------------------------
  TrackMaker::~TrackMaker() {
    // These ought to be here, but appear to cause the job
    // to hang at termination time.
    // if( m_smearer    ) delete m_smearer ; 
    // if( m_mcSelector ) delete m_mcSelector ;
    delete m_ncutter;
    delete m_tesIO;
  } 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode TrackMaker::initialize(){
    MsgStream log( messageService(), name() ) ;
    
    // Make a particle properties table. This will become a Gaudi service.
    HepMC::IO_PDG_ParticleDataTable pdg_io("PDGTABLE");
    m_particleTable = *pdg_io.read_particle_data_table();
    m_particleTable.make_antiparticles_from_particles();
    
    //=======================================================
    // Set up NCutter to select the required HepMC::GenParticles
    HepMC_helper::IMCselector* charged = new 
      HepMC_helper::IsCharged() ;
    HepMC_helper::IMCselector* cuts = new 
      HepMC_helper::MCCuts( m_mcPtMin, m_mcEtaMax ) ;
    HepMC_helper::IMCselector* finalState = new 
      HepMC_helper::IsFinalState();
    vector<HepMC_helper::IMCselector*> selectors;
    selectors.push_back(finalState);
    selectors.push_back(charged);
    selectors.push_back(cuts);
    
    m_ncutter = new  HepMC_helper::NCutter(selectors);
    
    delete cuts;
    delete charged;
    delete finalState;
    //=======================================================

    m_tesIO = new TesIO();

    //get the Global Event Data from singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    int lumi             = ged->lumi();
    int randSeed         = ged->randSeed() ;
    m_smearer            = new TrackSmearer(m_muConfig, randSeed, log);
    log << MSG::DEBUG << "got lumi " << endreq;
    HeaderPrinter hp("Atlfast Track Maker:", log);
    hp.add("Luminosity              ", lumi);        
    hp.add("Minimum four-vector Pt  ", m_mcPtMin);
    hp.add("Maximum four-vector Eta ", m_mcEtaMax);
    hp.add("Maximum Transverse Vertex ", m_vtMax);
    hp.add("Maximum Longitudinal Vertex ", m_vlMax);
    hp.add("Do Smearing             ", m_doSmearing);
    hp.add("B Field                 ", m_bField);
    hp.add("Muon Configuration      ", m_muConfig);
    hp.add("Random Number Seed      ", randSeed);
    hp.add("MC location             ", m_MC_eventLocation);
    hp.add("Output Location         ", m_outputLocation);    
    hp.print();
    log << MSG::INFO << "Initialised successfully " << endreq ;    
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode TrackMaker::finalize(){
    
    MsgStream log( messageService(), name() ) ;
    
    log << MSG::INFO << "Finalised successfully " << endreq ;
    
    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  //  This algorithm creates smeared Tracks passing cuts. 
  //  It scans the list of HepMC::GenParticles from the Monte Carlo truth. 
  //  It creates a Track object by smearing the Trajectory track parameters
  //  corresponding to the HepMC::GenParticle. 
  // 
  //  It then applies kinematic criteria on the properties of the Track
  //  Those which pass the cuts are kept.
  //  Finally all successful ReconstructedParticles are added to 
  //the event store.
  //
  
  StatusCode TrackMaker::execute( ){

    //................................
    // make a message logging stream
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Execute() " << endreq;
    

    //.........................................................
    // We now need to access the HepMC event record(s) in the event store
    // and iterate over them to extract HepMC::GenParticles.
    // This method will only select the particles which are required
    // and add them to a local store (mcParticles)
    
    MCparticleCollection  mcParticles ;
    TesIoStat stat; 
    std::string message;
    
    stat = m_tesIO->getMC( mcParticles, m_ncutter );
    message = stat?  "Found MCparitcles in TES":"No MC particles found in TES";
    log<<MSG::DEBUG << message <<" "<<mcParticles.size()<<endreq;
    // ......................
    // Make a container object which will store pointers to all isolated 
    // Tracks which are successfully created.  Since it is going to go 
    // into the event store, it has to be a special Athena container. 
    // This is typedefined in an include file
    
    TrackCollection* tracks = new TrackCollection ;
    

    //.........................................................
    // From each mc Track create a reconstructed Track. 
    // If all requirements are satisfied add to the collection
    
    Track* candidate ;
    MCparticleCollectionCIter src ;
    
    for( src = mcParticles.begin() ; src != mcParticles.end() ; ++src ) { 
      
      log << MSG::DEBUG << *src << endreq;
      candidate = this->create( *src ); 
      log << MSG::DEBUG << "Created: " << candidate << endreq;
      
      if( this->isAcceptable( candidate ) ) {    
      log << MSG::DEBUG << "Passed cuts" << endreq;
      tracks->push_back( candidate ) ; 
      log << MSG::DEBUG << "pushed back track" <<endreq;
      }
      else {
        log << MSG::DEBUG << "Failed cuts" << endreq;
      delete candidate ;
      }
    }
    
    //......................................
    // Register any Tracks which have been successfilly created in the 
    // transient event store. If there are none then we delete the empty list.

    stat = m_tesIO->store(tracks, m_outputLocation);
    message = stat?  "Stored tracks in TES":"Error storing tracks in TES";
    log<<MSG::DEBUG << message <<" "<<tracks->size()<<endreq;
    
    return stat; 

  }
  
  
  
  //----------------------------------
  // create()
  //----------------------------------
  
  // Takes a single MC Particle and creates a Track
  // 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
  
  Track* TrackMaker::create ( const HepMC::GenParticle* src ){
    
    MsgStream log( messageService(), name() ) ;
    // Construct the Trajectory parameters corresponding to the source particle
    TrackTrajectory originalTrajectory = this->extractTrajectory( src ) ;

    if (m_doSmearing) {
      
      // Ask our Smearer make a smeared set of Trajectory parameters
      // If this needs to be different for different particles then
      // it can be done here. Best  method of dispatching to be decided later
      Track myTrack(originalTrajectory, src);
      Track smearedTrack = m_smearer->smear(myTrack);
      
      // Return a new track
      Track* newTrack =  new Track( smearedTrack );
      return newTrack;
    } 
    else 
      {
        // Just return the track with unsmeared four-momentum
        return new Track( originalTrajectory, src ) ;
      }
  }
  
  
  
  //-------------------------------------------
  // 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
  // however there are no post creation conditions applied to tracks
  // in this version 
  
  //  bool TrackMaker::isAcceptable ( Track* candidate )
  bool TrackMaker::isAcceptable ( Track* candidate ){
    HepPoint3D vertex = candidate->trajectory().startPoint();
    int pdg = abs(candidate->pdg_id() );
    if ( vertex.z() > m_vlMax ) return false;

    // test track for vl and vt
    if ( pdg == 11 || pdg == 13) {
      if ( vertex.perp() > m_vtMax ) return false;
    }
    return true ;
  }
  
  
  // Extract Trajectory parameters from HepMC::GenParticle
  TrackTrajectory 
  TrackMaker::extractTrajectory( const HepMC::GenParticle* particle ) {
    // look up vertex
    Hep3Vector vertex = (particle->production_vertex())->position().getV();
    // Get 4-momentum and explicitly convert to 3-momentum
    Hep3Vector threeMomentum = particle->momentum().getV() ;
    // Obtain charge of particle
    HepMC::ParticleData* data = m_particleTable[particle->pdg_id()] ;
    double charge = data->charge() ; 
    
    return TrackTrajectory( threeMomentum, vertex, charge, m_bField ) ;
    
  }
  
  
} // end namespace bracket

---