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 "AtlfastEvent/ChargeService.h"

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

#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Units/SystemOfUnits.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;

namespace Atlfast 
{

  TrackMaker::TrackMaker ( const std::string& name, ISvcLocator* pSvcLocator ) : 
    Algorithm( name, pSvcLocator ), 
    m_tesIO(0),
    m_ncutter(0),
    m_smearer(0),
    m_chargeService(0)
  {
    
    // Setting the parameter defaults - these can be changed by user via
    // jobOptions
    m_mcPtMin           = 0.5*GeV;
    m_mcEtaMax          = 2.5;
    m_vlMax             = 40.0*cm;
    m_vtMax             = 5.05*cm;
    m_doSmearing        = true;
    m_bField            = 2.0;
    m_muonParamFile     = "atlfastDatafiles/Atlfast_MuonResParam_CSC.dat";
    m_pionParamFile     = "atlfastDatafiles/Atlfast_PionResParam_DC1_NewUnits.dat";
    m_electronParamFile = "atlfastDatafiles/Atlfast_ElectronResParam_CSC.dat";
    m_bremParamFile     = "atlfastDatafiles/Atlfast_ElectronBremParam_CSC.dat";
    
    // Default paths for entities in the TES
    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( "OutputLocation",        m_outputLocation ) ;
    declareProperty( "MuonParamFile",         m_muonParamFile );
    declareProperty( "PionParamFile",         m_pionParamFile );
    declareProperty( "ElectronParamFile",     m_electronParamFile );
    declareProperty( "ElectronBremParamFile", m_bremParamFile );
  }
  
  //----------------------------
  // Destructor
  //----------------------------
  TrackMaker::~TrackMaker() 
  {

    if (m_ncutter) delete m_ncutter;
    if (m_tesIO) delete m_tesIO;
    if (m_smearer) delete m_smearer;
    if (m_chargeService) delete m_chargeService;
  } 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode TrackMaker::initialize()
  {
    MsgStream log( messageService(), name() ) ;
    
    //get the Global Event Data from singleton pattern
    GlobalEventData* ged         = GlobalEventData::Instance();
    int lumi                     = ged->lumi();
    int randSeed                 = ged->randSeed() ;
    m_mcLocation                 = ged->mcLocation();  
    std::vector<int> invisibles  = ged->invisibles();

    log << MSG::INFO << "Got invisibles " << invisibles << endreq ;

    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
    
    //=======================================================
    // Set up NCutter to select the required HepMC::GenParticles
    HepMC_helper::IMCselector* charged;
    try
    {
      charged = new HepMC_helper::IsCharged() ;
    }
    catch(std::string errMsg)
    {
      log << MSG::ERROR << "Error making IsCharged " << errMsg<< endreq ;
    }
    catch(...)
    {
      log << MSG::ERROR << "Unknown Error making IsCharged " << endreq ;
    }
    HepMC_helper::IMCselector* cuts            = new HepMC_helper::MCCuts( m_mcPtMin, m_mcEtaMax ) ;
    HepMC_helper::IMCselector* finalState      = new HepMC_helper::IsFinalState();
    HepMC_helper::IMCselector* fromHardScatter = new HepMC_helper::IsFromHardScatter();
    HepMC_helper::IMCselector* invisible       = new HepMC_helper::RejectType(invisibles);


    vector<HepMC_helper::IMCselector*> selectors;
    selectors.push_back(fromHardScatter);
    selectors.push_back(finalState);
    selectors.push_back(charged);
    selectors.push_back(cuts);
    selectors.push_back(invisible);
    
    m_ncutter = new HepMC_helper::NCutter(selectors);
    
    delete cuts;
    delete charged;
    delete finalState;
    delete fromHardScatter;
    delete invisible;


    try
    {
      m_chargeService = new ChargeService() ;
    }
    catch(std::string errMsg)
    {
      log << MSG::ERROR << "Error making ChargeService " << errMsg<< endreq ;
    }
    catch(...)
    {
      log << MSG::ERROR << "Unknown Error making ChargeService " << endreq ;
    }
    m_smearer = new TrackSmearer( randSeed, log, m_muonParamFile, 
                                                 m_pionParamFile,
                                                 m_electronParamFile,
                                                 m_bremParamFile );
    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 Parameters         ",     m_muonParamFile);
    hp.add("Pion Parameters         ",     m_pionParamFile);    
    hp.add("Electron Parameters     ",     m_electronParamFile);    
    hp.add("Electron Brem Parameters",     m_bremParamFile);    
    hp.add("Random Number Seed      ",     randSeed);
    hp.add("MC location             ",     m_mcLocation);
    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 MCparticles 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 )
  {
    
    // 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 )
  {
    HepPoint3D vertex = candidate->trajectory().startPoint();
    int pdg = abs( candidate->pdg_id() );
    if ( abs( 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
    double charge=m_chargeService->operator()(particle);
    
    return TrackTrajectory( threeMomentum, vertex, charge, m_bField ) ;
    
  }
  
} // end namespace bracket

---