#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 (const std::string &iProcessName, std::vector< const char * > &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.

394 {

395 }

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

private

Definition at line 402 of file HepMCSplitter.cc.

References run_regression::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.

399 {

400 }

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);
   }
 
   
 }