FinalStateParticleHistogramMaker.cxx

Go to the documentation of this file.

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

#include <cmath> 
#include <algorithm>
#include <iomanip>
#include <float.h>

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


#include "CLHEP/Units/SystemOfUnits.h"

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

namespace Atlfast {
  using std::vector;
  using std::cout;
  using std::endl;
  using std::ios;
  
  // ================================================
  // FillHistograms methods
  // ================================================

  void FillHistograms::operator()(pair< string, HepMC_helper::MCselectorWrapper*> &pairPtr){

    // Called once per MCselectorWrapper, to make plots corresponding
    // to one particular set of selectors
    
    string identifier = pairPtr.first;
    HepMC_helper::MCselectorWrapper* wrapperPtr = pairPtr.second;

    MCparticleCollection mcParticles;
    const HepMC_helper::IMCselector* imcSelector = 
      wrapperPtr->asIMCselector();
    TesIoStat stat = m_tesIO->getMC( mcParticles, imcSelector) ;
    if (!stat){
      // Problem getting the GenParticles
      m_log << MSG::DEBUG << "Problem in getMC using " << identifier << " selectors" << endreq;
    }
    
    m_log << MSG::DEBUG << mcParticles.size() << " GenParticles selected with " 
          << identifier << " selectors" << endreq;

    Fill(mcParticles,identifier);

  }

  void FillHistograms::Fill(MCparticleCollection &mcParticles, string &identifier){
    
    // Fill histograms
    
    // Histograms have a prescribed order from 
    // FinalStateParticleHistogramMaker::BookHistograms
    
    // Histogram 0 - Multiplicity
    m_histogramContainer[identifier].at(0)->fill(float(mcParticles.size()),1.0);
    
    float sumE = 0., sumPx = 0., sumPy = 0., sumPz = 0.;
    
    MCparticleCollection::const_iterator genItr = mcParticles.begin();
    for (;genItr!=mcParticles.end();++genItr){

      // Histograms 1-4 - GenParticle energy, px, py and pz
      m_histogramContainer[identifier].at(1)->fill((*genItr)->momentum().e(),1.0);
      m_histogramContainer[identifier].at(2)->fill((*genItr)->momentum().px(),1.0);
      m_histogramContainer[identifier].at(3)->fill((*genItr)->momentum().py(),1.0);
      m_histogramContainer[identifier].at(4)->fill((*genItr)->momentum().pz(),1.0);

      // Add to sums
      sumE += (*genItr)->momentum().e();
      sumPx += (*genItr)->momentum().px();
      sumPy += (*genItr)->momentum().py();
      sumPz += (*genItr)->momentum().pz();

    }
    
    // Histograms 5-8 - E, px, py and pz sums of all selected GenParticles
    m_histogramContainer[identifier].at(5)->fill(sumE,1.0);
    m_histogramContainer[identifier].at(6)->fill(sumPx,1.0);
    m_histogramContainer[identifier].at(7)->fill(sumPy,1.0);
    m_histogramContainer[identifier].at(8)->fill(sumPz,1.0);
    
    return;

  }
  
  // ================================================
  // FinalStateParticleHistogramMaker methods
  // ================================================
  
  FinalStateParticleHistogramMaker::FinalStateParticleHistogramMaker
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){
    
    
    // Default paths for entities in the TES
    m_inputLocation    = "/Event/McEventCollection";
    
    declareProperty( "InputLocation", m_inputLocation ) ;

    SetHistogramLimits();

  }
  
  
  FinalStateParticleHistogramMaker::~FinalStateParticleHistogramMaker() {} 
  
  void FinalStateParticleHistogramMaker::PutIntoContainer(std::string &identifier, 
                                                          HepMC_helper::MCselectorWrapper* wrapper){
    std::pair<string,HepMC_helper::MCselectorWrapper*> thePair;
    thePair.first = identifier;
    thePair.second = wrapper;
    m_selectorContainer.push_back(thePair);
    return;
  }
  
  //---------------------------------
  // initialise() 
  //---------------------------------

  StatusCode FinalStateParticleHistogramMaker::initialize()
  {

    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Initialising" << endreq;
    
    string identifier;
    HepMC_helper::MCselectorWrapper* wrapper = 0;
    
    //set up MC selector.
    GlobalEventData* ged = GlobalEventData::Instance();

    AllSelectorSetups(identifier,wrapper);
    Z0SelectorSetups(identifier,wrapper);
    VisibleToCalSelectorSetups(identifier,wrapper,ged);
    InvisibleToAtlasSelectorSetups(identifier,wrapper,ged);
    IsFinalStateSelectorSetups(identifier,wrapper);
    BSelectorSetups(identifier,wrapper);
    ElectronSelectorSetups(identifier,wrapper);
    PhotonSelectorSetups(identifier,wrapper);
    MuonSelectorSetups(identifier,wrapper);

    m_mcLocation       = ged -> mcLocation();
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
    
    return StatusCode::SUCCESS ;

  }
  
  void FinalStateParticleHistogramMaker::SetHistogramLimits(){
    
    declareProperty("MultiplicityBins", m_multBins = 100);
    declareProperty("MultiplicityLow",  m_multLow = 0.);
    declareProperty("MultiplicityHigh", m_multHigh = 1000.);

    declareProperty("EnergyBins",       m_eBins = 100);
    declareProperty("EnergyLow",        m_eLow = 0.);
    declareProperty("EnergyHigh",       m_eHigh = 50.*GeV);
    declareProperty("PxBins",           m_pxBins = 100);
    declareProperty("PxLow",            m_pxLow = -20.*GeV);
    declareProperty("PxHigh",           m_pxHigh = 20.*GeV);
    declareProperty("PyBins",           m_pyBins = 100);
    declareProperty("PyLow",            m_pyLow = -20.*GeV);
    declareProperty("PyHigh",           m_pyHigh = 20.*GeV);
    declareProperty("PzBins",           m_pzBins = 100);
    declareProperty("PzLow",            m_pzLow = -20.*GeV);
    declareProperty("PzHigh",           m_pzHigh = 20.*GeV);

    declareProperty("EnergySumBins",    m_eSumBins = 100);
    declareProperty("EnergySumLow",     m_eSumLow = 13000.*GeV);
    declareProperty("EnergySumHigh",    m_eSumHigh = 15000.*GeV);
    declareProperty("PxSumBins",        m_pxSumBins = 100);
    declareProperty("PxSumLow",         m_pxSumLow = -2.*GeV);
    declareProperty("PxSumHigh",        m_pxSumHigh = 2.*GeV);
    declareProperty("PySumBins",        m_pySumBins = 100);
    declareProperty("PySumLow",         m_pySumLow = -2.*GeV);
    declareProperty("PySumHigh",        m_pySumHigh = 2.*GeV);
    declareProperty("PzSumBins",        m_pzSumBins = 100);
    declareProperty("PzSumLow",         m_pzSumLow = -2.*GeV);
    declareProperty("PzSumHigh",        m_pzSumHigh = 2.*GeV);

    return;

  }

  void FinalStateParticleHistogramMaker::AllSelectorSetups(string &identifier,
                                                           HepMC_helper::MCselectorWrapper* wrapper){
    identifier = "all";    
    wrapper = new HepMC_helper::MCselectorWrapper( new HepMC_helper::All() ) ;
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::Z0SelectorSetups(string &identifier,
                                                          HepMC_helper::MCselectorWrapper* wrapper){
    
    identifier = "z0selector"; 
    wrapper = new HepMC_helper::MCselectorWrapper( new HepMC_helper::SelectZ0(0.0, FLT_MAX));
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::VisibleToCalSelectorSetups(string &identifier,
                                                                    HepMC_helper::MCselectorWrapper* wrapper,
                                                                    GlobalEventData* ged){
 
    // Same selector as used in CellMaker

    identifier = "visibleToCal";
    wrapper = new HepMC_helper::MCselectorWrapper(ged -> visibleToCal());
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::InvisibleToAtlasSelectorSetups(string &identifier,
                                                                        HepMC_helper::MCselectorWrapper* wrapper,
                                                                        GlobalEventData* ged){
    
    // Same selector as used in EventHeaderMaker::escapedMomentum

    identifier = "invisibleToAtlas";
    wrapper = new HepMC_helper::MCselectorWrapper(ged -> visibleToAtlas());
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::IsFinalStateSelectorSetups(string &identifier,
                                                                    HepMC_helper::MCselectorWrapper* wrapper){
    identifier = "isFinalState";
    wrapper = new HepMC_helper::MCselectorWrapper(new HepMC_helper::IsFinalState() );
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;
    
  }

  void FinalStateParticleHistogramMaker::BSelectorSetups(string &identifier,
                                                         HepMC_helper::MCselectorWrapper* wrapper){
    
    identifier = "bSelector";    
    wrapper = 
      new HepMC_helper::MCselectorWrapper( new HepMC_helper::SelectJetTag (ParticleCodes::BQUARK, 
                                                                           5.0*GeV, 
                                                                           2.5));
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::ElectronSelectorSetups(string &identifier,
                                                                HepMC_helper::MCselectorWrapper* wrapper){
 
    // Same selectors as used in DefaultReconstructedParticleMaker/ElectronMaker

    identifier = "electronSelector";
    DefaultReconstructedParticleSelectorSetups(wrapper,11,0.0*GeV,100.0); // Current defaults
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::PhotonSelectorSetups(string &identifier,
                                                              HepMC_helper::MCselectorWrapper* wrapper){
    
    // Same selectors as used in DefaultReconstructedParticleMaker/PhotonMaker

    identifier = "photonSelector";
    DefaultReconstructedParticleSelectorSetups(wrapper,22,0.0*GeV,100.0); // Current defaults
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;
    
  }
  
  void FinalStateParticleHistogramMaker::MuonSelectorSetups(string &identifier,
                                                            HepMC_helper::MCselectorWrapper* wrapper){

    // Same selectors as used in DefaultReconstructedParticleMaker/MuonMaker

    identifier = "muonSelector"; 
    DefaultReconstructedParticleSelectorSetups(wrapper,13,0.5*GeV,100.0); // Current defaults
    BookHistograms(identifier);
    PutIntoContainer(identifier,wrapper);
    return;

  }

  void FinalStateParticleHistogramMaker::DefaultReconstructedParticleSelectorSetups
  (HepMC_helper::MCselectorWrapper* &wrapper, int particleType, double ptMin, double etaMax){
    
    // Common code for all DefaultReconstructedParticles

    HepMC_helper::IMCselector* typeSelector = 
      new HepMC_helper::SelectType(particleType);
    HepMC_helper::IMCselector* kineSelector = 
      new HepMC_helper::MCCuts(ptMin, etaMax);
    HepMC_helper::IMCselector* fstaSelector = 
      new HepMC_helper::IsFinalState();
    
    vector<HepMC_helper::IMCselector*> selectors;
    selectors.push_back(fstaSelector);
    selectors.push_back(typeSelector);
    selectors.push_back(kineSelector);   

    HepMC_helper::NCutter* ncutter = new HepMC_helper::NCutter(selectors);
    wrapper = new HepMC_helper::MCselectorWrapper(ncutter);    

    delete typeSelector;
    delete kineSelector;
    delete fstaSelector;
    return;
    
  }
  
  void FinalStateParticleHistogramMaker::BookHistograms(string &id){
    
    int nHistAll = 0;
    vector<IHistogram1D*> &hists = m_histogramContainer[id];
    string path = "/stat/finalstate_"+id+"/";
    string of = " of GenParticles from '"+id+"' selectors";

    hists.push_back(histoSvc()->book(path,++nHistAll,"Multiplicity"+of,m_multBins,m_multLow,m_multHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"Energy"+of,m_eBins,m_eLow,m_eHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"px"+of,m_pxBins,m_pxLow,m_pxHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"py"+of,m_pyBins,m_pyLow,m_pyHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"pz"+of,m_pzBins,m_pzLow,m_pzHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"Energy sum"+of,m_eSumBins,m_eSumLow,m_eSumHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"px sum"+of,m_pxSumBins,m_pxSumLow,m_pxSumHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"py sum"+of,m_pySumBins,m_pySumLow,m_pySumHigh));
    hists.push_back(histoSvc()->book(path,++nHistAll,"pz sum"+of,m_pzSumBins,m_pzSumLow,m_pzSumHigh));

    return;

  }

  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode FinalStateParticleHistogramMaker::finalize()
  {
    // .. put any finalisation in here...
    
    return StatusCode::SUCCESS ;
  }

  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  
  StatusCode FinalStateParticleHistogramMaker::execute( )
  {
    //................................
    // make a message logging stream
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Executing" << endreq;
    
    
    // apply a function to each element in m_selectorContainer
    //SB: no side effect in the constructor of this functor, and no
    //    useful public method to retrieve a particular "sum" variable from
    //    the STL algorithm, so it is not useful to retrieve the output of
    //    this for_each loop...
    /*FillHistograms mp =*/  
    std::for_each( m_selectorContainer.begin(),
                   m_selectorContainer.end(),
                   FillHistograms(m_tesIO,log,m_histogramContainer) );
    
    return StatusCode::SUCCESS;
  }
  
}  // end of namespace bracket

---