test_hepmc.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------
// 
// ----------------------------------------------------------------------

#include <iostream>
#include <sstream>
#include <vector>
#include <cmath>

#include "HepMC/GenEvent.h"
#include "HepMC/SimpleVector.h"

#include "test_hepmc.h"

using namespace std;

// ----------------------------------------------------------------------
bool testHepMC( const HepMC::GenEvent* genEvt )
{
  double totalPx = 0;
  double totalPy = 0;
  double totalPz = 0;
  double totalE  = 0;
  vector< int > negEnPart;
  vector< int > tacheons;
  vector< int > orphans;
  vector< int > unstNoEnd;
  vector< int > unDecPi0;
  
  const int k_pdg = 15;
  const double k_energy_diff = 1000.;
  const double k_cm_energy = 7000000.;
  //const double k_cm_energy = 14000000.;

  for ( HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin();
         pitr != genEvt->particles_end(); ++pitr ) {

      // Check for NaNs
    if ( std::isnan((*pitr)->momentum().px()) || std::isinf((*pitr)->momentum().px()) ||
      std::isnan((*pitr)->momentum().py()) || std::isinf((*pitr)->momentum().py()) ||
      std::isnan((*pitr)->momentum().pz()) || std::isinf((*pitr)->momentum().pz()) ||
      std::isnan((*pitr)->momentum().e())  || std::isinf((*pitr)->momentum().e()) ) {

      cerr << "NaN (Not A Number) found in the event record" << endl;
      return false;
    }

      // Check for undecayed pi0
    if ( (*pitr)->pdg_id() == 111 && !(*pitr)->end_vertex() ) {
      unDecPi0.push_back( (*pitr)->barcode());
    }


    // Check for unstables with no end vertex, such as undecayed gluons
    if ( (*pitr)->status() == 2 && !(*pitr)->end_vertex() ) {
      // NB. PYTHIA 6.4xx with MSTP(95)=2 produces zero momentum gluons
      bool notdumglu = true;

      /// @todo BPK: "temporarily" disabled: WHY?
      notdumglu = false;
      //if ((*pitr)->pdg_id() == 21 && (*pitr)->momentum().mag() == 0.) {   // new version of HepMC uses 'rho' insted of 'mag'
      if ((*pitr)->pdg_id() == 21 && (*pitr)->momentum().rho() == 0.) {
        cout << "Dummy gluon found! BARCODE=" << (*pitr)->barcode() << endl;
        notdumglu = false;
      }
      if (notdumglu) {
        unstNoEnd.push_back((*pitr)->barcode());
      }
    }


    // Check for production vertices of stable and unstable particles (documentary particles are not included in the check)
    if ( (*pitr)->status() == 1 || (*pitr)->status() == 2) {
      // Sequential resonance decays in PYTHIA externals (e.g. AcerMC tt~) don't have bosons hooked up to anywhere so skip
      // Also, a pythia bug has zero parent strings attached.
      // if (!(*pitr)->production_vertex()) orphans.push_back((*pitr)->barcode());

      // "Temporarily" disable check (WHY?)
      bool tmpdisable = true;
      if (!(*pitr)->production_vertex()) {
        if (!tmpdisable) {
          if (!(abs((*pitr)->pdg_id()) > 20 && abs((*pitr)->pdg_id()) < 27) &&
              !(abs((*pitr)->pdg_id()) < 10)) {
            if (abs((*pitr)->pdg_id()) != 2101 &&
                abs((*pitr)->pdg_id()) != 2103 && abs((*pitr)->pdg_id()) != 2203 &&
                abs((*pitr)->pdg_id()) != 2201) {
              orphans.push_back((*pitr)->barcode());
              cout << "Orphan found! BARCODE=" << (*pitr)->barcode() << ", ID= " << (*pitr)->pdg_id() << endl;
            }
          }
        }
      } else if (!((*pitr)->production_vertex()->particles_in_size()) &&  !tmpdisable) {
    // HERWIG record inserts initial particles with production vertex but no particle attached to it, check for px, py > 0
        if ((*pitr)->momentum().px() > 1e-7 ||  (*pitr)->momentum().py() > 1e-7) {
          orphans.push_back((*pitr)->barcode());
        }
      }
    }


    if ( (*pitr)->status() == 1 && !(*pitr)->end_vertex() ) {
      totalPx += (*pitr)->momentum().px();
      totalPy += (*pitr)->momentum().py();
      totalPz += (*pitr)->momentum().pz();
      totalE  += (*pitr)->momentum().e();
      if ( (*pitr)->momentum().e() < 0. ) negEnPart.push_back((*pitr)->barcode());
        //      double mass2 = (*pitr)->momentum().e()*(*pitr)->momentum().e() -
        //        (*pitr)->momentum().px()*(*pitr)->momentum().px() -
        //        (*pitr)->momentum().py()*(*pitr)->momentum().py() -
        //        (*pitr)->momentum().pz()*(*pitr)->momentum().pz();
        //      if ( mass2 < 0. ) tacheons.push_back((*pitr)->barcode());
      double aener = fabs((*pitr)->momentum().e());
      if ( aener < fabs((*pitr)->momentum().px()) ||
           aener < fabs((*pitr)->momentum().py()) ||
           aener < fabs((*pitr)->momentum().pz()) ) {
        tacheons.push_back((*pitr)->barcode());
      }
    }

    int p_id = (*pitr)->pdg_id();
    int tau_child = 0;
    if (abs(p_id) == k_pdg && (*pitr)->status() < 3) {
      HepMC::GenVertex *vtx = (*pitr)->end_vertex();
      if (vtx) {
        double p_energy = 0;
        for ( HepMC::GenVertex::particle_iterator desc = (vtx)->particles_begin(HepMC::descendants);
              desc != (vtx)->particles_end(HepMC::descendants); ++desc) {
          if (abs((*desc)->pdg_id() )== k_pdg) tau_child =1;
          if ((*desc)->status()==1) p_energy += (*desc)->momentum().e();
        }
        if (fabs( p_energy - (*pitr)->momentum().e()) > k_energy_diff && !tau_child) {
          cout << "EnergySum (decay products) - "
               << "Energy (original particle) > " << k_energy_diff << " MeV, "
               << "Event #" << genEvt->event_number() << ", "
               << "Barcode of the original particle = " << (*pitr)->barcode() << endl;
        }
      } else {
        cout << "UNDECAYED PARTICLE WITH PDG_ID = " << k_pdg << endl;
      }
    }
  }


    // Energy Balance
  double loste = fabs(totalE - k_cm_energy);
  if ( loste > k_energy_diff ) {
    cout << "ENERGY BALANCE FAILED : E-difference = " << loste << " MeV" << endl;
    return false;
  }

    // Momentum Balance
  if ( fabs(totalPx) > k_energy_diff || fabs(totalPy) > k_energy_diff || fabs(totalPz) > k_energy_diff ) {
    cout << "MOMENTUM BALANCE FAILED : SumPx = " << totalPx << " SumPy = " <<  totalPy << " SumPz = " <<  totalPz << " MeV"<< endl;
    return false;
  }

    // Negative energy particles
  if (negEnPart.size() > 0) {
    stringstream ss;
    ss << "NEGATIVE ENERGY PARTICLES FOUND : BARCODES =";
    for (vector<int>::const_iterator b = negEnPart.begin(); b != negEnPart.end(); ++b)
      ss << " " << *b;
    cout << ss.str() << endl;
    return false;
  }

    // Tacheons
  if (tacheons.size() > 0) {
    stringstream ss;
    ss << "PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND : BARCODES =";
    for (vector<int>::const_iterator b = tacheons.begin(); b != tacheons.end(); ++b)
      ss << " " << *b;
    cout << ss.str() << endl;
    return false;
  }

    // Stable or unstable particles with no parent (no documentaries)
  if (orphans.size() > 0) {
    stringstream ss;
    ss << "Particle with no parent found: BARCODES =";
    for (vector<int>::const_iterator b = orphans.begin(); b != orphans.end(); ++b)
      ss << " " << *b;
    cout << ss.str() << endl;
    return false;
  }

    // Unstable particles with no decay vertex
  if (unstNoEnd.size() > 0) {
    stringstream ss;
    ss << "Unstable particle with no decay vertex found: BARCODES =";
    for (vector<int>::const_iterator b = unstNoEnd.begin(); b != unstNoEnd.end(); ++b)
      ss << " " << *b;
    cout << ss.str() << endl;
    return false;
  }

    // Undecayed pi0
  if (unDecPi0.size() > 0) {
    stringstream ss;
    ss << "pi0 with no decay vertex found: BARCODES =";
    for (vector<int>::const_iterator b = unDecPi0.begin(); b != unDecPi0.end(); ++b)
      ss << " " << *b;
    cout << ss.str() << endl;;
    return false;
  }

  return true; 
}

---