/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

HepMCtoG4.cc
Go to the documentation of this file.
00001 
00010 #include "HepMCtoG4.h"
00011 #include "G4DataHelpers/G4DhPrimaryVertexInformation.h"
00012 #include "G4DataHelpers/G4DhPrimaryParticleInformation.h"
00013 
00014 #include "HepMC/GenParticle.h"
00015 #include "HepMC/GenEvent.h"
00016 #include "HepMC/GenVertex.h"
00017 
00018 #include "GaudiAlg/GaudiCommon.h"  //djaffe
00019 #include "GaudiKernel/ParticleProperty.h" // particle mass
00020 
00021 
00022 #include "G4PrimaryParticle.hh"
00023 #include "G4PrimaryVertex.hh"
00024 #include "G4OpticalPhoton.hh"
00025 
00026 #include "CLHEP/Vector/LorentzVector.h"
00027 #include "CLHEP/Geometry/Normal3D.h"
00028 
00029 
00030 
00031 HepMCtoG4::HepMCtoG4(const std::string& type,
00032                      const std::string& name,
00033                      const IInterface* parent)
00034     : GaudiTool(type,name,parent)
00035     , m_pps(0) // djaffe
00036 {
00037     declareInterface<IHepMCtoG4>(this);
00038     
00039     declareProperty("ZeroTime",m_zeroTime = false,
00040                     "Set true and all times will be reset so the signal vertex is at t=0.");
00041 }
00042 
00043 HepMCtoG4::~HepMCtoG4()
00044 {
00045 }
00046 
00047 
00048 StatusCode HepMCtoG4::initialize()
00049 {
00050     m_pps = svc<IParticlePropertySvc>( "ParticlePropertySvc", true ) ; //djaffe
00051     return StatusCode::SUCCESS;
00052 }
00053 
00054 StatusCode HepMCtoG4::finalize()
00055 {
00056   return StatusCode::SUCCESS;
00057 }
00058 // djaffe: Set optical photon mass to be zero and get the mass of all other particles
00059 //         from the particle property service. If the mass isn't available from the
00060 //         service or the calculated mass doesn't agree with the mass from the service
00061 //         then issue a warning.
00062 //         No warning is issued for nuclei that have been created by GenDecay with zero energy,
00063 //         set Trac ticket #364 for details.
00064 G4PrimaryParticle*
00065 HepMCtoG4::makeG4particle(const HepMC::GenParticle& part, const HepMC::GenEvent* event, const HepMC::GenVertex* vertex)
00066 {
00067     HepMC::FourVector p = part.momentum();
00068 
00069     G4PrimaryParticle* g4part = 0;
00070     
00071     double m, cm ;
00072     if (part.pdg_id() == 20022 /*opticalphoton*/) {
00073         g4part = new G4PrimaryParticle(G4OpticalPhoton::Definition(),
00074                                        p.px(),p.py(),p.pz());
00075         m = 0. ;
00076     }
00077     else {
00078         g4part = new G4PrimaryParticle(part.pdg_id(),p.px(),p.py(),p.pz());
00079         m = cm = sqrt(fabs( p.t()*p.t() - (p.x()*p.x() + p.y()*p.y() + p.z()*p.z()) ));
00080 
00081         ParticleProperty*  ppp = m_pps->findByStdHepID( part.pdg_id() );
00082         if ( ppp  ) {
00083             m = ppp->mass();
00084             double r = 0.;
00085             if ( m > 0 ) { r = std::abs(cm-m)/m; } 
00086             if ( r > 0.001 || std::abs(cm-m) > 0.01*MeV ) {
00087                 warning() << "calculated mass = " << cm
00088                           << " differs from IParticlePropertySvc mass " << m
00089                           << " for pdg_id " << part.pdg_id()
00090                           << ". Using IParticlePropertySvc mass. " << endreq;
00091             }
00092         }
00093         else {
00094             if ( p.t() < 0.01*MeV && std::abs(part.pdg_id()) > 1000000000 ) {
00095                 debug() << "As expected, there is no IParticlePropertySvc mass for pdg_id "
00096                         << part.pdg_id() << " which is a zero-energy, recoil nucleus created by GenDecay. Use calculated mass " << m 
00097                         << " from particle energy " << p.t() << " momentum " << sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z())
00098                         << endreq;
00099             }
00100             else {
00101                 warning() << "cannot find IParticlePropertySvc mass for pdg_id "
00102                           << part.pdg_id() << " and status " << part.status() << " . Using calculated mass " << m 
00103                           << " from particle energy " << p.t() << " momentum " << sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z())
00104                           << endreq;
00105             }
00106         }
00107     }
00108     g4part->SetMass(m);
00109     HepMC::ThreeVector pol = part.polarization().normal3d();
00110     g4part->SetPolarization(pol.x(),pol.y(),pol.z());
00111 
00112     g4part->SetUserInformation(new G4DhPrimaryParticleInformation(event,vertex,&part));
00113     return g4part;        
00114 }
00115 
00116 StatusCode HepMCtoG4::convert(const HepMC::GenEvent& input,
00117                               std::vector<G4PrimaryVertex*>& output)
00118 {
00119 
00120     // Catch HepMC::GenEvent weighting
00121     bool hasEventWeight = false;
00122     double eventWeight = 1.0;
00123     // If HepMC event weights container has only one value, apply to
00124     // all vertices as described in HepMC v2 documentation.
00125     if( input.weights().size() == 1 ){ 
00126         hasEventWeight = true; 
00127         eventWeight *= input.weights().front();
00128     }
00129 
00130     // Loop over vertices in the event
00131     HepMC::GenEvent::vertex_const_iterator
00132         iVtx = input.vertices_begin(),
00133         doneVtx = input.vertices_end();
00134     for (/*nop*/; doneVtx != iVtx; ++iVtx) {
00135             
00136         G4PrimaryVertex* g4vtx = new G4PrimaryVertex;
00137         bool first_time = true;
00138         
00139         // Catch HepMC::GenVertex weighting
00140         bool hasVertexWeight = false;
00141         double vertexWeight = 1.0;
00142         // If HepMC::GenVertex weights container has only one value,
00143         // apply to current vertex as described in HepMC v2
00144         // documentation.
00145         if( (*iVtx)->weights().size() == 1 ){ 
00146             hasVertexWeight = true;
00147             vertexWeight *= (*iVtx)->weights().front();
00148         }
00149         if( hasEventWeight || hasVertexWeight ){
00150             g4vtx->SetWeight( eventWeight*vertexWeight );
00151         }
00152 
00153         // Loop over particles in the vertex
00154         HepMC::GenVertex::particles_out_const_iterator 
00155             iPart = (*iVtx)->particles_out_const_begin(),
00156             donePart = (*iVtx)->particles_out_const_end();
00157         for (/*nop*/; donePart != iPart; ++iPart) {
00158                 
00159             // Only keep particles that are important for tracking
00160             if ((*iPart)->status() != 1) continue;
00161 
00162             if (first_time) {
00163                 first_time = false;
00164                 HepMC::FourVector v = (*iPart)->production_vertex()->position();
00165                 g4vtx->SetPosition(v.x(),v.y(),v.z());
00166                 if (m_zeroTime) 
00167                     g4vtx->SetT0(0.0);
00168                 else
00169                     g4vtx->SetT0(v.t());
00170                 g4vtx->SetUserInformation(new G4DhPrimaryVertexInformation(&input,*iVtx));
00171             }
00172             g4vtx->SetPrimary(makeG4particle(*(*iPart),&input,*iVtx));
00173         }
00174 
00175         if (!g4vtx->GetNumberOfParticle()) {
00176             delete g4vtx;
00177             continue;
00178         }
00179 
00180         output.push_back(g4vtx);
00181     }
00182     return StatusCode::SUCCESS;
00183 }
00184 
00185 
00186 // Local Variables: **
00187 // c-basic-offset:4 **
00188 // End: **
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:05:42 for G4DataHelpers by doxygen 1.7.4