HepMCSplitter Class Reference

#include <MyAna/HepMCSplitter/src/HepMCSplitter.cc>

Inheritance diagram for HepMCSplitter:

Public Member Functions
	HepMCSplitter (const edm::ParameterSet &)
	~HepMCSplitter ()
Public Member Functions inherited from edm::EDProducer
	EDProducer ()
ModuleDescription const &	moduleDescription () const
virtual	~EDProducer ()
Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
	ProducerBase ()
void	registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...
virtual	~ProducerBase ()
Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const
	EDConsumerBase ()
ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const
void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const
void	labelsForToken (EDGetToken iToken, Labels &oLabels) const
void	modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const
void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const
bool	registeredToConsumeMany (TypeID const &, BranchType) const
void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
virtual	~EDConsumerBase ()

Private Member Functions
virtual void	beginJob ()
int	cntStableParticles (const HepMC::GenEvent *evtUE)
virtual void	endJob ()
virtual void	produce (edm::Event &, const edm::EventSetup &)

Private Attributes
bool	_doUe
bool	_doZtautau
edm::InputTag	_input
std::vector< std::string >	_todo

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType
Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList
Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels
Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()
static void	fillDescriptions (ConfigurationDescriptions &descriptions)
static void	prevalidate (ConfigurationDescriptions &descriptions)
Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)
EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)
template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)
ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...
template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()
void	consumesMany (const TypeToGet &id)
template<BranchType B>
void	consumesMany (const TypeToGet &id)
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)
template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 74 of file HepMCSplitter.cc.

Constructor & Destructor Documentation

HepMCSplitter::HepMCSplitter ( const edm::ParameterSet &  iConfig )

explicit

Definition at line 106 of file HepMCSplitter.cc.

References _doUe, _doZtautau, _todo, Exception, and spr::find().

                                                          :
      _input(iConfig.getParameter<edm::InputTag>("input") ),
      _todo(iConfig.getParameter< std::vector<std::string> >("todo"))
{

  // Oscar producer uses string, not input tag :(
  
  
  _doUe =   std::find(_todo.begin(), _todo.end(), "UE")!=_todo.end();
  _doZtautau =   std::find(_todo.begin(), _todo.end(), "Ztautau")!=_todo.end();
  
  
  if (_todo.size()==1 ) {
    produces<edm::HepMCProduct>();
    
    if (!_doUe && !_doZtautau )
      throw cms::Exception("") << "todo wrong - select from Ztautau and UE \n";
  } else if (_todo.size()==2 ) {
    produces<edm::HepMCProduct>("Ztautau");
    produces<edm::HepMCProduct>("UE");
    if (!_doUe || !_doZtautau )
      throw cms::Exception("") << "todo wrong - select from Ztautau and UE \n";
  } else {
    throw cms::Exception("") << "todo wrong\n";
  }
}

HepMCSplitter::~HepMCSplitter

(

)

Definition at line 133 of file HepMCSplitter.cc.

{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

Member Function Documentation

void HepMCSplitter::beginJob ( void  )

privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 393 of file HepMCSplitter.cc.

{
}

int HepMCSplitter::cntStableParticles ( const HepMC::GenEvent *  evtUE )

private

Definition at line 402 of file HepMCSplitter.cc.

References runTheMatrix::ret.

Referenced by produce().

                                                              {
  
  int ret = 0;
  HepMC::GenEvent::particle_const_iterator parIt, parItE;
  parIt = evt->particles_begin();
  parItE = evt->particles_end();
  for ( ;parIt != parItE; ++parIt){  
    if ( (*parIt)->status()==1) ++ ret;
  
  }

  return ret;
}

void HepMCSplitter::endJob ( void  )

privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 399 of file HepMCSplitter.cc.

                      {
}

void HepMCSplitter::produce ( edm::Event &  iEvent, const edm::EventSetup &  iSetup  )

privatevirtual

Implements edm::EDProducer.

Definition at line 148 of file HepMCSplitter.cc.

References _doUe, _doZtautau, _input, funct::abs(), cntStableParticles(), gather_cfg::cout, Exception, GenParticle::GenParticle, edm::Event::getByLabel(), edm::Event::put(), and dt_dqm_sourceclient_common_cff::reco.

{
  using namespace edm;
  using namespace reco;

  // Get the product 

  Handle<edm::HepMCProduct> prodIn;
  iEvent.getByLabel(_input, prodIn);

  

  //  find Z decay and prod vertex
  int ZprodVtx, ZdecayVtx = 0;
  bool vtxFound = false;
  HepMC::GenEvent::particle_const_iterator parIt, parItE;
  parIt = prodIn->GetEvent()->particles_begin();
  parItE = prodIn->GetEvent()->particles_end();
  for ( ;parIt != parItE; ++parIt){  
      int pdg = abs( (*parIt)->pdg_id() ) ;
      if (pdg == 23){
    ZdecayVtx = (*parIt)->end_vertex()->barcode();
    ZprodVtx = (*parIt)->production_vertex()->barcode();
        vtxFound=true;
        break;
      }
  }

  if ( !vtxFound) throw cms::Exception("Z vtx not found");
//   std::cout << " Found Z - prod vtx " << ZprodVtx 
//       << " decay vtx " << ZdecayVtx 
//       << std::endl;

  
  
  // iterate over all vertices, mark those comming from Z
  std::set<int> barcodesZ;
  barcodesZ.insert(ZdecayVtx);
  
  HepMC::GenEvent::vertex_const_iterator itVtx, itVtxE;
  itVtx = prodIn->GetEvent()->vertices_begin();
  itVtxE = prodIn->GetEvent()->vertices_end();

  for(;itVtx!=itVtxE;++itVtx){
    if ( barcodesZ.find( (*itVtx)->barcode() )!=barcodesZ.end()){
      
      // iterate over out particles. Copy their decay vertices.to Z vertices
      HepMC::GenVertex::particles_out_const_iterator parIt, parItE;
      parIt=(*itVtx)->particles_out_const_begin();
      parItE=(*itVtx)->particles_out_const_end();
      for ( ; parIt!=parItE; ++parIt){
        if ((*parIt)->end_vertex()){
          barcodesZ.insert((*parIt)->end_vertex()->barcode());
          //std::cout << " Inserted " << (*parIt)->end_vertex()->barcode() << std::endl;
        }
      }
      
    }
  }
  
  barcodesZ.insert(ZprodVtx);
  barcodesZ.insert(ZdecayVtx);
  
  //std::cout << " Z vertices identified "  << barcodesZ.size() << std::endl;

  // prepare final products
  //     std::auto_ptr<edm::HepMCProduct> prodTauTau(new edm::HepMCProduct());  

  std::auto_ptr<HepMC::GenEvent> evtUE(new HepMC::GenEvent());
  std::auto_ptr<HepMC::GenEvent> evtTauTau(new HepMC::GenEvent());
  
  // Copy the vertices
  itVtx = prodIn->GetEvent()->vertices_begin();
  itVtxE = prodIn->GetEvent()->vertices_end();
  int newBarcode = 1;
  
  std::map<int, HepMC::GenVertex*> barcodeMapUE;
  std::map<int, HepMC::GenVertex*> barcodeMapTauTau;
  for(;itVtx!=itVtxE;++itVtx){
    bool isZ = barcodesZ.find( (*itVtx)->barcode() )!=barcodesZ.end();
    
    bool leadToZ = false;
    if ( (*parIt)->production_vertex()  ) {
      leadToZ = std::abs( (*itVtx)->barcode()) < std::abs(ZprodVtx);
    }
    
    
    // special case - Z production vertex. Put in both
    if (isZ || (*itVtx)->barcode()==ZprodVtx || leadToZ ){
      HepMC::GenVertex* newvertex = new HepMC::GenVertex((*itVtx)->position(), 
                                               (*itVtx)->id(), 
                                               (*itVtx)->weights() );
      
      newvertex->suggest_barcode(  -(++newBarcode) );
      //std::cout << "Adding vtx " << (*itVtx)->barcode() << std::endl;
      barcodeMapTauTau[(*itVtx)->barcode()] = newvertex;
      evtTauTau->add_vertex(newvertex);
    } 
    
    if (!isZ || (*itVtx)->barcode()==ZprodVtx){
      HepMC::GenVertex* newvertex = new HepMC::GenVertex((*itVtx)->position(), 
                                             (*itVtx)->id(), 
                                               (*itVtx)->weights() );
      
      newvertex->suggest_barcode(  -(++newBarcode) );
      barcodeMapUE[ (*itVtx)->barcode()] = newvertex;
      evtUE->add_vertex(newvertex);
    }
  }
  
  
  
  //std::cout << " Vertices copied " << std::endl;
  

  
  // iterate over particles again. Copy and assign
  parIt = prodIn->GetEvent()->particles_begin();
  parItE = prodIn->GetEvent()->particles_end();
  for ( ;parIt != parItE; ++parIt){  
    
    // check if prod vertex is Z vertex.
    bool fromZprodVtx = false;
    if ( (*parIt)->production_vertex() ){
      fromZprodVtx = ( (*parIt)->production_vertex()->barcode() == ZprodVtx);
    }
    
    // check if from Z
    bool fromZ = false;
    if ( (*parIt)->production_vertex() ){
      fromZ = barcodesZ.find( (*parIt)->production_vertex()->barcode() )!=barcodesZ.end();
    }
    
    
    // incoming particle (proton) or interaction that lead to Z boson. Copy to both events
    bool leadToZ = false;
    //if ( (*parIt)->production_vertex()  && (*parIt)->end_vertex()) {
    if ( (*parIt)->end_vertex()) {
      //leadToZ = std::abs((*parIt)->production_vertex()->barcode()) < std::abs(ZprodVtx);
      leadToZ = std::abs((*parIt)->end_vertex()->barcode()) < std::abs(ZdecayVtx);
    }
    
    //if ( !(*parIt)->production_vertex() || leadToZ ) {
    if ( leadToZ ) {
    
      HepMC::GenParticle* newparticle = new HepMC::GenParticle(*(*parIt));
      newparticle->suggest_barcode(newBarcode++);
      if ( (*parIt)->end_vertex()){
        //XX Add xcheck if this barcode is in map
        barcodeMapUE[ (*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
      } else {
        //std::cout << "XXX particle without vertices!" << std::endl;
      }
      
      newparticle = new HepMC::GenParticle(**parIt);
      newparticle->suggest_barcode(newBarcode++);
      if ( (*parIt)->end_vertex()){
        //XX Add xcheck if this barcode is in map
        //std::cout << "Trying: " << (*parIt)->end_vertex()->barcode() << std::endl;
        barcodeMapTauTau[ (*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
      } else {
        //std::cout << "XXX particle without vertices!" << std::endl;
      }
    
    } else if (fromZprodVtx)     { // Special case - particle comes from Z prod vertex
      bool isZ = (*parIt)->pdg_id();
      HepMC::GenParticle* newparticle = new HepMC::GenParticle(**parIt);
      newparticle->suggest_barcode(newBarcode++);
      if (isZ){
        //XX Add xcheck if this barcode is in map
        if ( (*parIt)->end_vertex())  
          barcodeMapTauTau[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
        if ((*parIt)->production_vertex()) 
          barcodeMapTauTau[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
      } else {
        if ((*parIt)->end_vertex())  
          barcodeMapUE[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
        if ((*parIt)->production_vertex()) 
          barcodeMapUE[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
      }
    } else {
      // Generic case
      HepMC::GenParticle* newparticle = new HepMC::GenParticle(**parIt);
      newparticle->suggest_barcode(newBarcode++);
      if (fromZ){
          //XX Add xcheck if this barcode is in map
        if ((*parIt)->end_vertex())  
          barcodeMapTauTau[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
        if ((*parIt)->production_vertex()) 
          barcodeMapTauTau[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
      } else {
        if ((*parIt)->end_vertex())  
  //         std::cout << "trying " << (*parIt)->end_vertex()->barcode() << std::endl;
          barcodeMapUE[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
        if ((*parIt)->production_vertex()) 
          barcodeMapUE[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
      }
    }
  } // end particles iteration
  
//   std::cout << "XXXXXXXXXXXXXXXXXX" << std::endl;
  
  int orgCnt = cntStableParticles(prodIn->GetEvent());
  int tautauCnt = cntStableParticles(evtTauTau.get());
  int UEcnt = cntStableParticles(evtUE.get());
  if (orgCnt != tautauCnt + UEcnt){
    std::cout << " XXX Warning - wrong num particles - ORG: " 
      <<  orgCnt
      << " tautau " << tautauCnt
      << " UE "     << UEcnt
      << std::endl;
  }
  
//   std::cout << "=====================================\n";
//   std::cout << "EV Tau \n";
//   //evtTauTau->print();
//   std::cout << "=====================================\n";
//   std::cout << "EV org \n";
//   //prodIn->GetEvent()->print();
//   
//   std::cout << "=====================================\n";
//   std::cout << "EV EU \n";
  //   evtUE->print();
  
  if (_doUe){
    std::auto_ptr<edm::HepMCProduct> prodUE(new edm::HepMCProduct());  
    prodUE->addHepMCData( evtUE.release() );
    if (_doZtautau) iEvent.put( prodUE, "UE");
    else iEvent.put( prodUE);
  }

  if (_doZtautau) {
    std::auto_ptr<edm::HepMCProduct> prodTauTau(new edm::HepMCProduct());  
//     std::cout << "===================================== " << iEvent.id() << std::endl;
//     evtTauTau->print();
    prodTauTau->addHepMCData( evtTauTau.release() );
    if (_doUe ) iEvent.put( prodTauTau, "Ztautau");
    else iEvent.put( prodTauTau);
  }

  
}

Member Data Documentation

bool HepMCSplitter::_doUe

private

Definition at line 90 of file HepMCSplitter.cc.

Referenced by HepMCSplitter(), and produce().

bool HepMCSplitter::_doZtautau

private

Definition at line 91 of file HepMCSplitter.cc.

Referenced by HepMCSplitter(), and produce().

edm::InputTag HepMCSplitter::_input

private

Definition at line 88 of file HepMCSplitter.cc.

Referenced by produce().

std::vector<std::string> HepMCSplitter::_todo

private

Definition at line 89 of file HepMCSplitter.cc.

Referenced by HepMCSplitter().