FinalStateParticleDumper.cxx

Go to the documentation of this file.

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

#include <cmath> 
#include <algorithm>
#include <iomanip>
#include <float.h>
// Gaudi includes
#include "GaudiKernel/DataSvc.h"


// This algorithm includes
#include "AtlfastEvent/Cell.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/SimpleKinematic.h"
#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "AtlfastUtils/HepMC_helper/HepMC_helper.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "CLHEP/Units/SystemOfUnits.h"

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

namespace Atlfast {
  using std::vector;
  using std::cout;
  using std::endl;
  using std::ios;
    
  FinalStateParticleDumper::FinalStateParticleDumper
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){
    
    
    // Default paths for entities in the TES
    m_inputLocation  = "/Event/McEventCollection";
    m_selectorName   = "All";
    
    declareProperty( "InputLocation", m_inputLocation ) ;
    declareProperty( "SelectorName",  m_selectorName ) ;
    
  }
  
  
  FinalStateParticleDumper::~FinalStateParticleDumper() {} 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------

  StatusCode FinalStateParticleDumper::initialize()
  {

    MsgStream log( messageService(), name() ) ;
    log<<MSG::DEBUG<<"Initialising "<<endreq;
    
    //set up MC selector.
    log<< "m_mcSelector  set to "<<m_selectorName<<endreq;
    GlobalEventData* ged = GlobalEventData::Instance();
    if(m_selectorName == "All") {      
      m_mcSelector = 
        new HepMC_helper::MCselectorWrapper( new HepMC_helper::All() ) ;
    }else if(m_selectorName == "bSelector") {
      m_mcSelector =  
        new HepMC_helper::MCselectorWrapper( new 
                                             HepMC_helper::SelectJetTag 
                                             (ParticleCodes::BQUARK, 
                                              5.0*GeV, 
                                              2.5))
        ;
    }else if(m_selectorName == "ElectronSelector"){ 
      HepMC_helper::IMCselector* pidSelector =  
        new HepMC_helper::SelectPid(0.0, FLT_MAX, ParticleCodes::ELECTRON);
      HepMC_helper::IMCselector* finalStateSelector = 
        new HepMC_helper::IsFinalState();

      vector<HepMC_helper::IMCselector*> selectors;
      selectors.push_back(finalStateSelector);
      selectors.push_back(pidSelector);

    
      HepMC_helper::NCutter* ncutter  = new  HepMC_helper::NCutter(selectors);
      m_mcSelector = new HepMC_helper::MCselectorWrapper(ncutter) ;
    
      delete pidSelector;
      delete finalStateSelector;

    }else if(m_selectorName == "PhotonSelector"){ 
      HepMC_helper::IMCselector* pidSelector =  
        new HepMC_helper::SelectPid(0.0, FLT_MAX, ParticleCodes::PHOTON);
      HepMC_helper::IMCselector* finalStateSelector = 
        new HepMC_helper::IsFinalState();

      vector<HepMC_helper::IMCselector*> selectors;
      selectors.push_back(finalStateSelector);
      selectors.push_back(pidSelector);

    
      HepMC_helper::NCutter* ncutter  = new  HepMC_helper::NCutter(selectors);
      m_mcSelector = new HepMC_helper::MCselectorWrapper(ncutter) ;
    
      delete pidSelector;
      delete finalStateSelector;

    }else if(m_selectorName == "MuonSelector"){ 
      HepMC_helper::IMCselector* pidSelector =  
        new HepMC_helper::SelectPid(0.0, FLT_MAX, ParticleCodes::MUON);
      HepMC_helper::IMCselector* finalStateSelector = 
        new HepMC_helper::IsFinalState();

      vector<HepMC_helper::IMCselector*> selectors;
      selectors.push_back(finalStateSelector);
      selectors.push_back(pidSelector);

    
      HepMC_helper::NCutter* ncutter  = new  HepMC_helper::NCutter(selectors);
      m_mcSelector = new HepMC_helper::MCselectorWrapper(ncutter) ;
    
      delete pidSelector;
      delete finalStateSelector;

    }else if(m_selectorName == "Z0selector"){ 
      m_mcSelector =  
        new HepMC_helper::MCselectorWrapper( new 
                                           HepMC_helper::SelectZ0(0.0, FLT_MAX)
                                           )
     ;

    }else if(m_selectorName == "Unseen") {

      std::vector<int> requiredTypes;
      requiredTypes.push_back(66); //ERW's favorite LSP number.
      requiredTypes.push_back(ParticleCodes::NU_E);
      requiredTypes.push_back(ParticleCodes::NU_MU);
      requiredTypes.push_back(ParticleCodes::NU_TAU);

      HepMC_helper::IMCselector* unseen  = 
        new HepMC_helper::Unseen(requiredTypes);

      HepMC_helper::IMCselector* finalStateSelector = 
        new HepMC_helper::IsFinalState();

      vector<HepMC_helper::IMCselector*> selectors;
      selectors.push_back(finalStateSelector);
      selectors.push_back(unseen);

    
      HepMC_helper::NCutter* ncutter  = new  HepMC_helper::NCutter(selectors);
      m_mcSelector = new HepMC_helper::MCselectorWrapper(ncutter) ;
    
      delete unseen;
      delete finalStateSelector;

    }else if(m_selectorName ==  "VisibleToCal"){
      m_mcSelector = 
        new HepMC_helper::MCselectorWrapper(ged -> visibleToCal());
    }else{
      m_mcSelector = 
        new HepMC_helper::MCselectorWrapper(new HepMC_helper::IsFinalState() ) ;
    }
    m_mcLocation       = ged -> mcLocation();
    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode FinalStateParticleDumper::finalize()
  {
    // .. put any finalisation in here...
    
    return StatusCode::SUCCESS ;
  }

  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  
  StatusCode FinalStateParticleDumper::execute( )
  {
    //................................
    // make a message logging stream
    
    MsgStream log( messageService(), name() ) ;
    log<<MSG::DEBUG<<"Executing "<<endreq;
    std::string mess;
    
    // read MC particles from TES
    MCparticleCollection  mcParticles ;
    const HepMC_helper::IMCselector* imcSelector = 
      m_mcSelector->asIMCselector();
    TesIoStat stat = m_tesIO->getMC( mcParticles, imcSelector) ;
    mess = stat? "Retrieved MC from TES ":"Failed MC retrieve from TES";
    log << MSG::DEBUG << mess << endreq;
    
    std::cout<<"Dumping  particles "<<std::endl;
    cout<<"------------------- "<<std::endl;

    std::sort(mcParticles.begin(), 
              mcParticles.end(), 
              HepMC_helper::AscendingEta());

    cout<<"Number particles found = "<<mcParticles.size()<<std::endl;

    /*
    for(;ip!=mcParticles.end();++ip){
      string space = " ";
      this->DumpParticle(space, *ip);
    }
    */

    string space = " ";
    log<<MSG::INFO<<"Flat Dump"<<endreq<<endreq;
    this->DumpFlat(space, mcParticles.begin(), mcParticles.end() );
    //    log<<MSG::INFO<<"Tree Dump"<<endreq<<endreq;
    //    this->DumpTree(space, mcParticles.begin(), mcParticles.end() );
    //this->DumpTwoGenerations(space, mcParticles.begin(), mcParticles.end() );
    log<<MSG::INFO<<"End of execute "<<endreq;    
    return StatusCode::SUCCESS;
  }
  void FinalStateParticleDumper::DumpParticle(std::string s, 
                                              const HepMC::GenParticle* p){
    cout<<setiosflags(ios::fixed);
    cout<<std::setprecision(3);
    cout<<s
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().phi()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().pseudoRapidity()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().perp()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().e()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->pdg_id()<<" "
      <<"phi, eta, pt, e"
      <<std::endl;
    cout<<s
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().px()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().py()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().pz()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->momentum().e()<<" "
      <<std::setw(15)<<std::setprecision(3)<<p->pdg_id()<<" "
      <<"  x,   y,  z, e"
      <<std::endl;
  }
  void FinalStateParticleDumper::DumpFlat(std::string space,
                                          MCparticleCollection::const_iterator ib,
                                          MCparticleCollection::const_iterator ie){
    MCparticleCollection::const_iterator ip=ib;
    for(; ip!=ie; ++ip){
      this->DumpParticle(space, *ip);
    }
  }
  void 
  FinalStateParticleDumper::DumpTwoGenerations(
                                               std::string space,
                                               MCparticleCollection::const_iterator ib,
                                               MCparticleCollection::const_iterator ie){
    std::string space2 = space +" ";
    MCparticleCollection::const_iterator ip=ib;
    HepMC::GenVertex::particle_iterator ip2; 
    cout<<"-------first"<<std::endl;
    cout << "******dopo first *******"<<std::endl;
for(; ip!=ie; ++ip){
  cout << "****particle at ip *****" << endl;
  this->DumpParticle(space, *ip);

      if((*ip)->end_vertex()){
        HepMC::GenVertex::particle_iterator firstChild = 
          (*ip)->end_vertex()->particles_begin(HepMC::children);
        
        HepMC::GenVertex::particle_iterator endChild = 
          (*ip)->end_vertex()->particles_end(HepMC::children);

        cout<<"----------second"<<std::endl;
        for(ip2=firstChild; ip2!=endChild; ++ip2){
          this->DumpParticle(space2, *ip2);
        }       
        cout<<"-------first"<<std::endl;
      }
}
  }
  void FinalStateParticleDumper::DumpTree(std::string space,
                                          MCparticleCollection::const_iterator ib,
                                          MCparticleCollection::const_iterator ie){
    MCparticleCollection::const_iterator ip=ib;
    space = space+" ";
    cout<<"-----------"<<endl;
    for(; ip!=ie; ++ip){
      this->DumpParticle(space, *ip);

      if((*ip)->end_vertex()){
        HepMC::GenVertex::particle_iterator firstChild = 
          (*ip)->end_vertex()->particles_begin(HepMC::children);
        
        HepMC::GenVertex::particle_iterator endChild = 
          (*ip)->end_vertex()->particles_end(HepMC::children);
        /*
        HepMC::GenVertex::particle_iterator secondChild = firstChild;
        ++secondChild;
        
        if(secondChild!=endChild) this->DumpTree(space,secondChild,endChild);
        */
        this->DumpTree(space, firstChild, endChild);    
      }
    }
  }
  void FinalStateParticleDumper::DumpTree(std::string space,
                                          HepMC::GenVertex::particle_iterator ib,
                                          HepMC::GenVertex::particle_iterator ie){
    HepMC::GenVertex::particle_iterator ip=ib;
    space = space+" ";
    cout<<"-----------"<<endl;
    for(; ip!=ie; ++ip){
      this->DumpParticle(space, *ip);
      if((*ip)->end_vertex()){
        HepMC::GenVertex::particle_iterator firstChild = 
          (*ip)->end_vertex()->particles_begin(HepMC::children);
        
        HepMC::GenVertex::particle_iterator endChild = 
          (*ip)->end_vertex()->particles_end(HepMC::children);

        /*
        HepMC::GenVertex::particle_iterator secondChild = firstChild;
        ++secondChild;

        if(secondChild != endChild) this->DumpTree(space, secondChild, endChild);
        */
        this->DumpTree(space, firstChild, endChild);
      }
    }
  }   

}  // end of namespace bracket

---