read_particles_from_HepMC.cc File Reference

#include "GeneratorInterface/ExternalDecays/interface/read_particles_from_HepMC.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

Go to the source code of this file.

Functions
HepMC::GenParticle *	findLastSelf (HepMC::GenParticle *x)
std::vector< SimpleParticle > *	getDaughters (HepMC::GenParticle *x)
bool	isFirst (HepMC::GenParticle *x)
int	readParticlesFromHepMC (const HepMC::GenEvent *event, SimpleParticle &X, SimpleParticle &tau, SimpleParticle &tau2, std::vector< SimpleParticle > &tau_daughters, std::vector< SimpleParticle > &tau2_daughters)

Function Documentation

HepMC::GenParticle* findLastSelf ( HepMC::GenParticle *  x )

inline

Find last self

Recursively finds the last particle with the same PDG ID on the list of its decay products

Definition at line 48 of file read_particles_from_HepMC.cc.

References findLastSelf(), AlCaHLTBitMon_ParallelJobs::p, and x.

Referenced by findLastSelf(), and readParticlesFromHepMC().

{
  if(!x->end_vertex()) return x;

  for(HepMC::GenVertex::particle_iterator p = x->end_vertex()->particles_begin(HepMC::children); p!=x->end_vertex()->particles_end(HepMC::children); ++p)
  {
    if( (*p)->pdg_id()==x->pdg_id() ) return findLastSelf( *p );
  }

  return x;
}

std::vector<SimpleParticle>* getDaughters ( HepMC::GenParticle *  x )

inline

Get daughters of HepMC::GenParticle

Recursively searches for final-state daughters of 'x'

Definition at line 10 of file read_particles_from_HepMC.cc.

References configurableAnalysis::GenParticle, getDaughters(), AlCaHLTBitMon_ParallelJobs::p, and createTree::pp.

Referenced by getDaughters(), CandMatcherBase< C1, C2 >::initMaps(), and readParticlesFromHepMC().

{
  std::vector<SimpleParticle> *daughters = new std::vector<SimpleParticle>();
  if(!x->end_vertex()) return daughters;

  // Check decay products of 'x'
  for(HepMC::GenVertex::particles_out_const_iterator p = x->end_vertex()->particles_out_const_begin(); p!=x->end_vertex()->particles_out_const_end(); ++p)
  {
    HepMC::GenParticle *pp = *p;
    HepMC::FourVector   mm = pp->momentum();

    // If the daughter of 'x' has its end vertex - recursively read
    // all of its daughters.
    if( pp->end_vertex() && pp->pdg_id()!=111)
    {
      std::vector<SimpleParticle> *sub_daughters = getDaughters(pp);
      daughters->insert(daughters->end(),sub_daughters->begin(),sub_daughters->end());
      
      delete sub_daughters;
    }
    // Otherwise - add this particle to the list of daughters.
    else if(pp->pdg_id()!=x->pdg_id())
    {
      SimpleParticle tp( mm.px(), mm.py(), mm.pz(), mm.e(), pp->pdg_id() );
      daughters->push_back(tp);
    }
  }

  return daughters;
}

bool isFirst ( HepMC::GenParticle *  x )

inline

Definition at line 60 of file read_particles_from_HepMC.cc.

References AlCaHLTBitMon_ParallelJobs::p, and parents.

                                        {
  for(HepMC::GenVertex::particle_iterator p = x->production_vertex()->particles_begin(HepMC::parents); p!=x->production_vertex()->particles_end(HepMC::parents); ++p){
    if(x->pdg_id()==(*p)->pdg_id()) return false; 
  }
  return true;
}

int readParticlesFromHepMC	(	const HepMC::GenEvent *	event,
		SimpleParticle &	X,
		SimpleParticle &	tau,
		SimpleParticle &	tau2,
		std::vector< SimpleParticle > &	tau_daughters,
		std::vector< SimpleParticle > &	tau2_daughters
	)

Definition at line 89 of file read_particles_from_HepMC.cc.

References funct::abs(), findLastSelf(), configurableAnalysis::GenParticle, getDaughters(), cuy::isFirst, and NULL.

{
  if(event==NULL) return 1;

  // Exctract particles from event
  HepMC::GenParticle *hX=NULL, *hTau=NULL, *hTau2=NULL;

  for(HepMC::GenEvent::particle_const_iterator it = event->particles_begin(); it!=event->particles_end(); ++it)
  {
    int pdgid = (*it)->pdg_id();
    if(
        (
          abs(pdgid)==37 ||
              pdgid ==25 ||
          abs(pdgid)==24 ||
              pdgid ==23          
        ) &&
        (*it)->end_vertex() &&
        (*it)->end_vertex()->particles_out_size()>=2
      )
    {
      hX = *it;
      HepMC::GenParticle *LhX=hX;
      if(!isFirst(hX)) continue;
      findLastSelf(LhX);
      hTau = hTau2 = NULL;

      for(HepMC::GenVertex::particle_iterator it2 = LhX->end_vertex()->particles_begin(HepMC::children); it2!=LhX->end_vertex()->particles_end(HepMC::children); ++it2)
      {
        if (abs( (*it2)->pdg_id() )==15  && (*it2)->end_vertex() ){
          if(!hTau)       hTau  = *it2;
          else if(!hTau2) hTau2 = *it2;
          else
          {
        edm::LogInfo("TauSpinnerCMSSW::read_particles_from_HepMC") << "TauSpiner: three taus in one decay" << std::endl;
        //std::cout<<"TauSpiner: three taus in one decay"<<std::endl;
            return 1;
          }
        }
        if (abs( (*it2)->pdg_id() )==16 )
        {
          if(!hTau2) hTau2 = *it2;
          else
          {
        edm::LogInfo("TauSpinnerCMSSW::read_particles_from_HepMC") << "TauSpiner: two neutrinos or two taus and neutrino in one decay" 
                                       << std::endl;
            //std::cout<<"TauSpiner: two neutrinos or two taus and neutrino in one decay"<<std::endl; 
            return 1;
          }
        }
      }
      if(hTau && hTau2) break;

    }
  }

  if(!hX) return 1;

  if(!hTau || !hTau2)
  {
    //edm::LogInfo("TauSpinnerCMSSW::read_particles_from_HepMC") << "TauSpiner: boson found but no proper tau pair or tau + neutrino found."<<std::endl;
    //std::cout<<"TauSpiner: boson found but no proper tau pair or tau + neutrino found."<<std::endl;
    return 1;
  }

  // Check for self-decaying taus
  hTau  = findLastSelf(hTau);
  hTau2 = findLastSelf(hTau2);

  // Fill SimpleParticles from HepMC particles
  X.setPx(hX->momentum().px());
  X.setPy(hX->momentum().py());
  X.setPz(hX->momentum().pz());
  X.setE (hX->momentum().e ());
  X.setPdgid(hX->pdg_id());

  tau.setPx(hTau->momentum().px());
  tau.setPy(hTau->momentum().py());
  tau.setPz(hTau->momentum().pz());
  tau.setE (hTau->momentum().e ());
  tau.setPdgid(hTau->pdg_id());

  tau2.setPx(hTau2->momentum().px());
  tau2.setPy(hTau2->momentum().py());
  tau2.setPz(hTau2->momentum().pz());
  tau2.setE (hTau2->momentum().e ());
  tau2.setPdgid(hTau2->pdg_id());

  // Create list of tau daughters
  std::vector<SimpleParticle> *buf = getDaughters(hTau);
  tau_daughters.clear();
  tau_daughters.insert(tau_daughters.end(),buf->begin(),buf->end());
  
  delete buf;

  // Second particle can be 2nd tau. In that case - read its daughters.
  // Otherwise it is nu_tau~
  buf = getDaughters(hTau2);
  tau2_daughters.clear();
  tau2_daughters.insert(tau2_daughters.end(),buf->begin(),buf->end());
  
  delete buf;

  return 0;
}