read_particles_from_HepMC.cc

Go to the documentation of this file.

#include "GeneratorInterface/TauolaInterface/interface/read_particles_from_HepMC.h"
using namespace TauSpinner;
/*******************************************************************************
  Get daughters of HepMC::GenParticle

  Recursively searches for final-state daughters of 'x'
*******************************************************************************/
inline std::vector<SimpleParticle> *getDaughters(HepMC::GenParticle *x)
{
  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;
}


/*******************************************************************************
  Find last self

  Recursively finds the last particle with the same PDG ID
  on the list of its decay products
*******************************************************************************/
inline HepMC::GenParticle *findLastSelf(HepMC::GenParticle *x)
{
  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;
}

inline bool isFirst(HepMC::GenParticle *x){
  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;
}

/*******************************************************************************
  Read HepMC::GenEvent.

  Read HepMC event from data file
  and return particles filtered out from the event.
  
  This routine is prepared for use with files generated by Pythia8.
  Fills:
  
  'X'              - Heavy particle (W+/-, H+/-, H, Z)
  'tau'            - first tau
  'tau2'           - second tau or nu_tau, if 'X' decays to one tau only
  'tau_daughters'  - daughters of 'tau'
  'tau2_daughters' - daughters of 'tau2' or empty list, if 'tau2' is nu_tau.
  
  Returns:
  0 - no more events to read               (finished processing the file)
  1 - no decay found in the event          (finished processing current event)
  2 - decay found and processed correctly.
      Event will continue to be processed
      with next function call. 
*******************************************************************************/
int readParticlesFromHepMC(const HepMC::GenEvent *event, SimpleParticle &X, SimpleParticle &tau, SimpleParticle &tau2, std::vector<SimpleParticle> &tau_daughters, std::vector<SimpleParticle> &tau2_daughters)
{
  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
          {
            std::cout<<"TauSpiner: three taus in one decay"<<std::endl;
            return 1;
          }
        }
        if (abs( (*it2)->pdg_id() )==16 )
        {
          if(!hTau2) hTau2 = *it2;
          else
          {
            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)
  {
    //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;
}