EvtHepMCEvent.cpp

Go to the documentation of this file.

 
/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
*                                                                      *
* This file is part of EvtGen.                                         *
*                                                                      *
* EvtGen is free software: you can redistribute it and/or modify       *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* EvtGen is distributed in the hope that it will be useful,            *
* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
* GNU General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *
* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
***********************************************************************/
 
#include "EvtGenBase/EvtHepMCEvent.hh"
 
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
 
EvtHepMCEvent::EvtHepMCEvent() :
    m_theEvent( nullptr ), m_translation( 0.0, 0.0, 0.0, 0.0 )
{
}
 
EvtHepMCEvent::~EvtHepMCEvent()
{
    this->deleteEvent();
}
 
void EvtHepMCEvent::deleteEvent()
{
    if ( m_theEvent != nullptr ) {
        m_theEvent->clear();
        delete m_theEvent;
        m_theEvent = nullptr;
    }
}
 
void EvtHepMCEvent::constructEvent( EvtParticle* baseParticle )
{
    EvtVector4R origin( 0.0, 0.0, 0.0, 0.0 );
    this->constructEvent( baseParticle, origin );
}
 
void EvtHepMCEvent::constructEvent( EvtParticle* baseParticle,
                                    EvtVector4R& translation )
{
    // This class does not take ownership of the base particle pointer.
    // Rather, it uses the base particle to construct the event.
 
    this->deleteEvent();
    if ( baseParticle == nullptr ) {
        return;
    }
 
    m_theEvent = new GenEvent( Units::GEV, Units::MM );
    m_translation = translation;
 
    // Use the recursive function addVertex to add a vertex with incoming/outgoing
    // particles. Adds a new vertex for any EvtParticles with decay daughters.
    // All particles are in the rest frame of the base particle ("lab frame").
 
    GenParticlePtr hepMCGenParticle =
        this->createGenParticle( baseParticle, EvtHepMCEvent::LAB );
 
    this->addVertex( baseParticle, hepMCGenParticle );
}
 
GenParticlePtr EvtHepMCEvent::createGenParticle( EvtParticle* theParticle,
                                                 int frameType )
{
    // Create an HepMC GenParticle, with the 4-momenta in the frame given by the frameType integer
    GenParticlePtr genParticle{ nullptr };
 
    if ( theParticle != nullptr ) {
        // Set the particle status integer to either stable or decayed
        int status( EvtHepMCEvent::STABLE );
        int nDaug = theParticle->getNDaug();
        if ( nDaug > 0 ) {
            status = EvtHepMCEvent::DECAYED;
        }
 
        // Get the 4-momentum (E, px, py, pz) for the EvtParticle.
        EvtVector4R p4( 0.0, 0.0, 0.0, 0.0 );
 
        if ( frameType == EvtHepMCEvent::RESTFRAME ) {
            p4 = theParticle->getP4Restframe();
        } else if ( frameType == EvtHepMCEvent::LAB ) {
            p4 = theParticle->getP4Lab();
        } else {
            p4 = theParticle->getP4();
        }
 
        // Convert this to the HepMC 4-momentum
        double E = p4.get( 0 );
        double px = p4.get( 1 );
        double py = p4.get( 2 );
        double pz = p4.get( 3 );
 
        FourVector hepMC_p4( px, py, pz, E );
 
        // Get the particle PDG integer id
        int PDGInt = EvtPDL::getStdHep( theParticle->getId() );
 
        genParticle = newGenParticlePtr( hepMC_p4, PDGInt, status );
    }
 
    return genParticle;
}
 
void EvtHepMCEvent::addVertex( EvtParticle* inEvtParticle,
                               GenParticlePtr inGenParticle )
{
    // This is a recursive function that adds GenVertices to the GenEvent for
    // the incoming EvtParticle and its daughters. We use two separate
    // pointers for the EvtParticle and GenParticle information: the former
    // to obtain the PDGId, 4-momenta, daughter and vertex positions, the latter to
    // set the incoming particle to the vertex. Note that the outgoing particle for
    // one vertex might be the incoming particle for another vertex - this needs to
    // be the same GenParticle pointer, hence the reason for using it as a 2nd argument
    // in this function.
 
    if ( m_theEvent == nullptr || inEvtParticle == nullptr ||
         inGenParticle == nullptr ) {
        return;
    }
 
    // Create the decay vertex
    FourVector vtxCoord = this->getVertexCoord( inEvtParticle );
    GenVertexPtr theVertex = newGenVertexPtr( vtxCoord );
 
    // Add the vertex to the event
    m_theEvent->add_vertex( theVertex );
 
    // Set the incoming particle
    theVertex->add_particle_in( inGenParticle );
 
    // Set the outgoing particles (decay products)
    int nDaug = inEvtParticle->getNDaug();
    int iDaug( 0 );
    // Loop over the daughters
    for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
        EvtParticle* evtDaughter = inEvtParticle->getDaug( iDaug );
        GenParticlePtr genDaughter =
            this->createGenParticle( evtDaughter, EvtHepMCEvent::LAB );
 
        if ( genDaughter != nullptr ) {
            // Add a new GenParticle (outgoing) particle daughter to the vertex
            theVertex->add_particle_out( genDaughter );
 
            // Find out if the daughter also has decay products.
            // If so, recursively run this function again.
            int nDaugProducts = evtDaughter->getNDaug();
 
            if ( nDaugProducts > 0 ) {
                // Recursively process daughter particles and add their vertices to the event
                this->addVertex( evtDaughter, genDaughter );
 
            }    // Have daughter products
 
        }    // hepMCDaughter != 0
 
    }    // Loop over daughters
}
 
FourVector EvtHepMCEvent::getVertexCoord( EvtParticle* theParticle )
{
    FourVector vertexCoord( 0.0, 0.0, 0.0, 0.0 );
 
    if ( theParticle != nullptr && theParticle->getNDaug() != 0 ) {
        // Get the position (t,x,y,z) of the EvtParticle, offset by the translation vector.
        // This position will be the point where the particle decays. So we ask
        // the position of the (1st) daughter particle.
        EvtParticle* daugParticle = theParticle->getDaug( 0 );
 
        if ( daugParticle != nullptr ) {
            EvtVector4R vtxPosition = daugParticle->get4Pos() + m_translation;
 
            // Create the HepMC 4 vector of the position (x,y,z,t)
            vertexCoord.setX( vtxPosition.get( 1 ) );
            vertexCoord.setY( vtxPosition.get( 2 ) );
            vertexCoord.setZ( vtxPosition.get( 3 ) );
            vertexCoord.setT( vtxPosition.get( 0 ) );
        }
    }
 
    return vertexCoord;
}