PrimaryGeneratorAction.cc

Go to the documentation of this file.

// $Header: /cvs/lcd/slic/src/PrimaryGeneratorAction.cc,v 1.32 2008/09/10 00:04:34 jeremy Exp $

#include "PrimaryGeneratorAction.hh"

// lcdd
#include "StringUtil.hh"

// slic
#include "LcioManager.hh"
#include "SlicApplication.hh"
#include "LcioMcpManager.hh"
#include "LcioMcpFilter.hh"
#include "EventSourceManager.hh"

// lcio
#include "IMPL/MCParticleImpl.h"
#include "EVENT/LCCollection.h"

// geant4
#include "G4StateManager.hh"

namespace slic
{
    PrimaryGeneratorAction::PrimaryGeneratorAction()
        : Module("PrimaryGeneratorAction", false)
    {}

    PrimaryGeneratorAction::~PrimaryGeneratorAction()
    {}

    void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent)
    {
        EventSourceManager* mgr = EventSourceManager::instance();

        // print begin message
#ifdef SLIC_LOG
        printBeginEventMessage( anEvent );
#endif

        // needs to be called here for setup purposes rather than in EventAction
        mgr->beginEvent(anEvent);

        // delegate event generation to the manager
        mgr->GeneratePrimaryVertex(anEvent);

        applyLorentzTransformation(anEvent);    

        if ( mgr->isEOF() ) {
            SlicApplication::instance()->setAborting(true);
        }
    }

    void PrimaryGeneratorAction::printBeginEventMessage(G4Event* anEvent)
    {
        log() << LOG::okay << ">>>> BeginEvent <" << StringUtil::toString( anEvent->GetEventID() ) << ">" << LOG::done;
    }

    void PrimaryGeneratorAction::applyLorentzTransformation(G4Event *evt)
    {
        const G4double alpha = EventSourceManager::instance()->getLorentzTransformationAngle();
 
        if (alpha == 0) return; // nothing to do

//#ifdef SLIC_LOG
//        log() << LOG::always << "applying Lorentz Transformation angle <" << alpha << LOG::done;
//#endif

        // parameters of the Lorentz transformation matrix
        const G4double gamma = sqrt(1 + sqr(tan(alpha)));
        const G4double betagamma = tan(alpha);

        // scan through all vertices and all valid primary particles
        for (int iVertex = 0; iVertex < evt->GetNumberOfPrimaryVertex(); iVertex++) 
        {
            G4PrimaryVertex *vertex = evt->GetPrimaryVertex(iVertex);
            for (int iParticle = 0; iParticle < vertex->GetNumberOfParticle(); iParticle++) 
            {
                G4PrimaryParticle *particle = vertex->GetPrimary(iParticle);
                // does the particle have a valid particle definition attached?
                if (particle->GetG4code()) 
                {
                    // before the transformation
                    const G4double px = particle->GetPx();
                    const G4double py = particle->GetPy();
                    const G4double pz = particle->GetPz();

                    //FG: take generator mass ( not necessarily euqal to PDG mass ) 
                    const G4double m  = particle->GetMass() ;

                    // after the transformation (boost in x-direction)
                    const G4double pxPrime = betagamma * sqrt(sqr(px) + sqr(py) + sqr(pz) + sqr(m)) + gamma * px;

                    // py and pz remain the same, E changes implicitly with px
                    particle->SetMomentum(pxPrime, py, pz);
                }
            }
            // the position of the vertex is not transformed here
        }

        //FG: now we need to apply the same transformation to the MCParticles:
        EVENT::LCCollection* col = LcioMcpManager::instance()->getInitialMcpCollection();

        if( col !=0 ) 
        { // if particle gun is used no collection exists ...

            int nMCP = col->getNumberOfElements() ;
    
            for(int i=0; i < nMCP ; ++i ) 
            {
      
                IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>( col->getElementAt(i) ) ;
      
                const double* p = mcp->getMomentum() ;
      
                // before the transformation
                double pPrime[3] ;
      
                const G4double m  = mcp->getMass() ;
      
                // after the transformation (boost in x-direction)
                pPrime[0] = betagamma * sqrt(sqr(p[0] ) + sqr(p[1] ) + sqr( p[2]) + sqr(m) ) + gamma * p[0] ;
      
                pPrime[1] = p[1] ;
                pPrime[2] = p[2] ;
      
                // py and pz remain the same, E changes implicitly with px
                mcp->setMomentum( pPrime );
            }
        }
    }
}

---