ReconstructedParticleHistogramMaker.cxx

Go to the documentation of this file.

// ================================================
// ReconstructedParticleHistogramMaker class Implementation
// ================================================
//
//
// Namespace Atlfast::
//
#include "AtlfastAlgs/ReconstructedParticleHistogramMaker.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/MCparticleCollection.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 "AtlfastAlgs/GlobalEventData.h"

// library includes
#include <cmath> 

// Gaudi includes
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/DataSvc.h"
#include "AIDA/IHistogram1D.h"
#include "GaudiKernel/IHistogramSvc.h"

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


//-------------------------------
// Gaudi magic lines.
// These are so that it can make instances of this algorithm
// You never need to know why or how !!!
//-------------------------------

static const AlgFactory<Atlfast::ReconstructedParticleHistogramMaker> s_factory ;
const IAlgFactory& ReconstructedParticleHistogramMakerFactory = s_factory ;


namespace Atlfast {
//--------------------------------
// Constructors and destructors
//-------------------------------

ReconstructedParticleHistogramMaker::ReconstructedParticleHistogramMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){




  // Default paths for entities in the TES
  m_inputLocation     = "/Event/Electrons";
  m_mcTruthLocation   = "/Event/McEventCollection";
  

  m_histStart         = 0;            // ID of first histogram
  m_particleType      = 11;           // Default particle type to examine
  m_histTitle         = "Electron ";  // Default prefixfor histogram titles

  declareProperty( "InputLocation", m_inputLocation ) ;
  declareProperty( "McTruthLocation", m_mcTruthLocation ); 
  declareProperty( "HistogramStart", m_histStart );
  declareProperty( "ParticleType", m_particleType );
  declareProperty( "HistogramTitle", m_histTitle );

}


ReconstructedParticleHistogramMaker::~ReconstructedParticleHistogramMaker(){
  MsgStream log(messageService(), name());
  log << MSG::INFO << "destructor" << endreq;

}


//---------------------------------
// initialise() 
//---------------------------------

StatusCode ReconstructedParticleHistogramMaker::initialize()
{
  // .. put any initialisation in here...
  MsgStream log(messageService(), name());
  log << MSG::DEBUG << "in initialize()" << endreq;
  
  m_nHist = m_histStart;
  
  book(m_h_multiplicity, string("Multiplicity"),20,0.0,20.0,-5.0, 5.0);
  log << MSG::DEBUG << "Made multiplicity histos" << endreq;  

  book(m_h_energy,       string("Energy"),100,0.0,500.0,-10.0, 10.0);
  log << MSG::DEBUG << "Made energy histos" << endreq;

  book(m_h_pt,           string("Pt"), 100,0.0,500.0,-10.0,10.0);
  log << MSG::DEBUG << "Made pt histos" << endreq;

  book(m_h_eta,          string("Eta"), 100,-3.0,3.0,-0.2,0.2);
  log << MSG::DEBUG << "Made eta histos" << endreq;

  book(m_h_phi,          string("Phi"),20,-M_PI,M_PI,-0.1,0.1);
  log << MSG::DEBUG << "Made phi histos" << endreq;

  book(m_h_theta,        string("Theta"),20,0.0,M_PI,-0.1,0.1);
  log << MSG::DEBUG << "Made thetahistos" << endreq;


  log << MSG::DEBUG << "Booked histograms OK" << endreq;

  GlobalEventData* ged = GlobalEventData::Instance();
  // load the location of the MC in StoreGate
  m_mcLocation       = ged -> mcLocation();


  //  m_tesIO = new TesIO(eventDataService());
  m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
  
  // 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;
  //=======================================================
  
  return StatusCode::SUCCESS ;
}



//---------------------------------
// finalise() 
//---------------------------------

StatusCode ReconstructedParticleHistogramMaker::finalize()
{
  // .. put any finalisation in here...
  MsgStream log(messageService(), name());
  log << MSG::INFO << "finalizing" << endreq;

  return StatusCode::SUCCESS ;
}




//----------------------------------------------
// execute() method called once per event
//----------------------------------------------

StatusCode ReconstructedParticleHistogramMaker::execute( )
{

  //................................
  // make a message logging stream

  MsgStream log( messageService(), name() ) ;
  log << MSG::DEBUG << "Executing " << endreq;

  // Access the reconstructed particles in the TES
  std::vector<ReconstructedParticle*> myParticles;
  // Access the collection of "McEvent" records in the TES
  std::vector<McEvent*> mcEvents;


  //Get the MC truth list
    MCparticleCollection  myMC_particles ;
    MCparticleCollectionCIter src ;
    
    TesIoStat stat = m_tesIO->getMC( myMC_particles, m_ncutter ) ;
    mess = stat? "Retrieved MC from TES ":"Failed MC retrieve from TES";
    log << MSG::DEBUG << mess << endreq;
 

   //check to see whether any reconstructed particles are found at the requested location
  if(!m_tesIO->copy(myParticles, m_inputLocation)){
    log << MSG::DEBUG << "no particles found in the TES. This is probably normal" << endreq; 
    // fill any histograms that have meaning in the absence of rec particles
    fill(m_h_multiplicity,0.0,float(myMC_particles.size()) );
    //log << MSG::DEBUG 
    //<< "filled multiplicity 0.0 and " 
    //<< float(myMC_particles.size()) 
    //<< endreq;
  } else {
    // fill all histograms
    //log << MSG::DEBUG << "Fill multiplicity" << endreq;
    fill(m_h_multiplicity,
         double(myParticles.size()),
         double(myMC_particles.size()) );

    //iterate over the electrons

    std::vector<ReconstructedParticle*>::const_iterator iter = myParticles.begin();
    for (; iter != myParticles.end(); ++iter) {
      const HepMC::GenParticle* truth = (*iter)->truth();
      HepLorentzVector momentum = (*iter)->momentum();
      HepLorentzVector trueMomentum = truth->momentum();
      
      //log << MSG::DEBUG << "Fill per-electron histos" << endreq;
      fill(m_h_pt,    (*iter)->pT(),                trueMomentum.perp());
      fill(m_h_eta,   (*iter)->eta(),           trueMomentum.pseudoRapidity());
      fill(m_h_energy,((*iter)->momentum()).e(),    trueMomentum.e());
      fill(m_h_phi,   (*iter)->phi(),               trueMomentum.phi());
      fill(m_h_theta, ((*iter)->momentum()).theta(),trueMomentum.theta());
    }
  }
  log << MSG::DEBUG << "Done the filling OK" << endreq;
  return StatusCode::SUCCESS ;
}



//--------------------------------------------------------------
//
void ReconstructedParticleHistogramMaker::book
      (
       std::vector<IHistogram1D*>& start, 
       const std::string title, 
       const int nbins, 
       const double xmin, 
       const double xmax,
       const double xminDiff,
       const double xmaxDiff
      ) 
{
  
  // book the histograms in turn, incrementing the counter before using it
  start.push_back(histoSvc()->book(
                                   "/stat/simple1d/",
                                   ++m_nHist,
                                   m_histTitle+" Rec "+title,
                                   nbins,
                                   xmin,
                                   xmax
                                   )
                  );

  start.push_back(histoSvc()->book(
                                   "/stat/simple1d/",
                                   ++m_nHist,
                                   m_histTitle+" Tru "+title,
                                   nbins,
                                   xmin,
                                   xmax
                                   )
                  );

  start.push_back(histoSvc()->book(
                                   "/stat/simple1d/",
                                   ++m_nHist,
                                   m_histTitle+" Dif "+title,
                                   nbins*2,
                                   xminDiff,
                                   xmaxDiff
                                   )
                  );
  start.push_back(histoSvc()->book(
                                   "/stat/simple1d/",
                                   ++m_nHist,
                                   m_histTitle+" Res "+title,
                                   nbins*2,
                                   -2.0,
                                   2.0
                                   )
                  );

  return;
}


//-------------------------------------------------------------------

void ReconstructedParticleHistogramMaker::fill(std::vector<IHistogram1D*>& start, const double rec, const double tru) {
 
  MsgStream log( messageService(), name() ) ;

  if (!start[0] || !start[1] || !start[2] || !start[3]) {
    log << MSG::WARNING << "trying to fill to a null histogram pointer" << endreq;
    return;
   }
  start[0]->fill(rec,1.0);
  start[1]->fill(tru,1.0);
  start[2]->fill(tru-rec,1.0);
  if (tru != 0.0) {
    start[3]->fill((tru-rec)/tru,1.0);
  }
  return;
}

} // end of namespace bracket

---