PrimaryGeneratorAction.cc

Go to the documentation of this file.
00001 // $Header: /cvs/lcd/slic/src/PrimaryGeneratorAction.cc,v 1.32 2008/09/10 00:04:34 jeremy Exp $
00002 
00003 #include "PrimaryGeneratorAction.hh"
00004 
00005 // lcdd
00006 #include "StringUtil.hh"
00007 
00008 // slic
00009 #include "LcioManager.hh"
00010 #include "SlicApplication.hh"
00011 #include "LcioMcpManager.hh"
00012 #include "LcioMcpFilter.hh"
00013 #include "EventSourceManager.hh"
00014 
00015 // lcio
00016 #include "IMPL/MCParticleImpl.h"
00017 #include "EVENT/LCCollection.h"
00018 
00019 // geant4
00020 #include "G4StateManager.hh"
00021 
00022 namespace slic
00023 {
00024     PrimaryGeneratorAction::PrimaryGeneratorAction()
00025         : Module("PrimaryGeneratorAction", false)
00026     {}
00027 
00028     PrimaryGeneratorAction::~PrimaryGeneratorAction()
00029     {}
00030 
00031     void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent)
00032     {
00033         EventSourceManager* mgr = EventSourceManager::instance();
00034 
00035         // print begin message
00036 #ifdef SLIC_LOG
00037         printBeginEventMessage( anEvent );
00038 #endif
00039 
00040         // needs to be called here for setup purposes rather than in EventAction
00041         mgr->beginEvent(anEvent);
00042 
00043         // delegate event generation to the manager
00044         mgr->GeneratePrimaryVertex(anEvent);
00045 
00046         applyLorentzTransformation(anEvent);    
00047 
00048         if ( mgr->isEOF() ) {
00049             SlicApplication::instance()->setAborting(true);
00050         }
00051     }
00052 
00053     void PrimaryGeneratorAction::printBeginEventMessage(G4Event* anEvent)
00054     {
00055         log() << LOG::okay << ">>>> BeginEvent <" << StringUtil::toString( anEvent->GetEventID() ) << ">" << LOG::done;
00056     }
00057 
00058     void PrimaryGeneratorAction::applyLorentzTransformation(G4Event *evt)
00059     {
00060         const G4double alpha = EventSourceManager::instance()->getLorentzTransformationAngle();
00061  
00062         if (alpha == 0) return; // nothing to do
00063 
00064 //#ifdef SLIC_LOG
00065 //        log() << LOG::always << "applying Lorentz Transformation angle <" << alpha << LOG::done;
00066 //#endif
00067 
00068         // parameters of the Lorentz transformation matrix
00069         const G4double gamma = sqrt(1 + sqr(tan(alpha)));
00070         const G4double betagamma = tan(alpha);
00071 
00072         // scan through all vertices and all valid primary particles
00073         for (int iVertex = 0; iVertex < evt->GetNumberOfPrimaryVertex(); iVertex++) 
00074         {
00075             G4PrimaryVertex *vertex = evt->GetPrimaryVertex(iVertex);
00076             for (int iParticle = 0; iParticle < vertex->GetNumberOfParticle(); iParticle++) 
00077             {
00078                 G4PrimaryParticle *particle = vertex->GetPrimary(iParticle);
00079                 // does the particle have a valid particle definition attached?
00080                 if (particle->GetG4code()) 
00081                 {
00082                     // before the transformation
00083                     const G4double px = particle->GetPx();
00084                     const G4double py = particle->GetPy();
00085                     const G4double pz = particle->GetPz();
00086 
00087                     //FG: take generator mass ( not necessarily euqal to PDG mass ) 
00088                     const G4double m  = particle->GetMass() ;
00089 
00090                     // after the transformation (boost in x-direction)
00091                     const G4double pxPrime = betagamma * sqrt(sqr(px) + sqr(py) + sqr(pz) + sqr(m)) + gamma * px;
00092 
00093                     // py and pz remain the same, E changes implicitly with px
00094                     particle->SetMomentum(pxPrime, py, pz);
00095                 }
00096             }
00097             // the position of the vertex is not transformed here
00098         }
00099 
00100         //FG: now we need to apply the same transformation to the MCParticles:
00101         EVENT::LCCollection* col = LcioMcpManager::instance()->getInitialMcpCollection();
00102 
00103         if( col !=0 ) 
00104         { // if particle gun is used no collection exists ...
00105 
00106             int nMCP = col->getNumberOfElements() ;
00107     
00108             for(int i=0; i < nMCP ; ++i ) 
00109             {
00110       
00111                 IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>( col->getElementAt(i) ) ;
00112       
00113                 const double* p = mcp->getMomentum() ;
00114       
00115                 // before the transformation
00116                 double pPrime[3] ;
00117       
00118                 const G4double m  = mcp->getMass() ;
00119       
00120                 // after the transformation (boost in x-direction)
00121                 pPrime[0] = betagamma * sqrt(sqr(p[0] ) + sqr(p[1] ) + sqr( p[2]) + sqr(m) ) + gamma * p[0] ;
00122       
00123                 pPrime[1] = p[1] ;
00124                 pPrime[2] = p[2] ;
00125       
00126                 // py and pz remain the same, E changes implicitly with px
00127                 mcp->setMomentum( pPrime );
00128             }
00129         }
00130     }
00131 }

Generated on Mon Jun 7 17:45:21 2010 for Simulator for the Linear Collider by  doxygen 1.5.4