read_particles_from_HepMC.cc File Reference

#include "GeneratorInterface/TauolaInterface/interface/read_particles_from_HepMC.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

Definition at line 44 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

Definition at line 8 of file read_particles_from_HepMC.cc.

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

Referenced by getDaughters(), 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 58 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 90 of file read_particles_from_HepMC.cc.

References funct::abs(), visDQMUpload::buf, gather_cfg::cout, findLastSelf(), GenParticle::GenParticle, getDaughters(), and cuy::isFirst.

                                                                      {
  if (event == nullptr)
    return 1;

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

  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 = nullptr;

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