#include "ConvertMCTruthToG4.h"

Inheritance diagram for g4b::ConvertMCTruthToG4:

Public Member Functions
	ConvertMCTruthToG4 ()
	Standard constructor and destructor. More...

virtual	~ConvertMCTruthToG4 ()

void	Reset ()
	Get ready to convert a new set of MCTruth objects. More...

void	Append (art::Ptr< simb::MCTruth > &mct)

void	Append (const simb::MCTruth *mct)

virtual void	GeneratePrimaries (G4Event *)

Private Attributes
std::vector< const simb::MCTruth * >	fConvertList
	List of MCTruth objects to convert for this spill. More...

std::map< G4int, G4int >	fUnknownPDG
	map of unknown PDG codes to instances More...

Static Private Attributes
static G4ParticleTable *	fParticleTable = 0
	Geant4's table of particle definitions. More...

Detailed Description

Definition at line 36 of file ConvertMCTruthToG4.h.

Constructor & Destructor Documentation

g4b::ConvertMCTruthToG4::ConvertMCTruthToG4 ( )

Standard constructor and destructor.

Definition at line 39 of file ConvertMCTruthToG4.cxx.

40 {

41 }

g4b::ConvertMCTruthToG4::~ConvertMCTruthToG4 ( )

virtual

Definition at line 44 of file ConvertMCTruthToG4.cxx.

References fUnknownPDG.

   {
     // Print out a list of "unknown" PDG codes we saw in the input.
     if ( ! fUnknownPDG.empty() ){
       std::ostringstream badtxt;
       std::map<G4int, G4int>::iterator i = fUnknownPDG.begin();
       for ( ; i != fUnknownPDG.end(); ++i ){
         int pdg = (*i).first;
         badtxt << "\n   Unknown PDG code = " << pdg
                 << ", seen " << (*i).second << " times.";
 
         const int genieLo = 2000000001;
         const int genieHi = 2000000202;
         if ( pdg >= genieLo && pdg <= genieHi ) {
           badtxt << "  (GENIE specific)";
         }
       }
 
       mf::LogWarning("ConvertPrimaryToGeant4")
         << "The following unknown PDG codes were present in the "
         << "simb::MCTruth input.\n"
         << "They were not processed by Geant4."
         << badtxt.str();
     }
 
   }

Member Function Documentation

void g4b::ConvertMCTruthToG4::Append ( art::Ptr< simb::MCTruth > & mct )

Add a new MCTruth object to the list of primary particles to be appended to the Geant4 event.

Definition at line 78 of file ConvertMCTruthToG4.cxx.

References art::Ptr< T >::get().

Referenced by g4b::G4Helper::G4Run().

   {
     this->Append( mct.get() );
   }

void g4b::ConvertMCTruthToG4::Append ( const simb::MCTruth * mct )

Add a new MCTruth object to the list of primary particles to be appended to the Geant4 event.

Definition at line 84 of file ConvertMCTruthToG4.cxx.

References fConvertList.

   {
     fConvertList.push_back( mct );
   }

void g4b::ConvertMCTruthToG4::GeneratePrimaries ( G4Event * event )

virtual

Required of any class that inherits from G4VUserPrimaryGeneratorAction; append primary particles to a G4Event object. This method is invoked by Geant4, and is not directly called by the user application.

Definition at line 90 of file ConvertMCTruthToG4.cxx.

References fConvertList, fParticleTable, fUnknownPDG, MF_LOG_DEBUG, simb::MCParticle::Momentum(), simb::MCParticle::PdgCode(), simb::MCParticle::Polarization(), g4b::PrimaryParticleInformation::SetMCTruth(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), x, y, Z, and z.

   {
     // Take the contents of MCTruth objects and use them to
     // initialize the G4Event.
   
     // A G4Event organizes its particles in terms of "vertices" and
     // "particles", like HepMC.  Unfortunately, ROOT doesn't use
     // HepMC, so the MCTruth objects aren't organized that way.
     // For most of the work that we'll ever do, there'll be only one
     // vertex in the event.  However, just in case there are multiple
     // vertices (e.g., overlays, double vertex studies) I want the
     // code to function properly.
   
     // So create a map of particle positions and associated
     // G4PrimaryVertex*.  Note that the map must use CLHEP's
     // LorentzVector, and not ROOT's, since ROOT does not define an
     // operator< for its physics vectors.
     std::map< CLHEP::HepLorentzVector, G4PrimaryVertex* >                  vertexMap;
     std::map< CLHEP::HepLorentzVector, G4PrimaryVertex* >::const_iterator  vi; 
 
     // For each MCTruth (probably only one, but you never know):
     // index keeps track of which MCTruth object you are using
     size_t index = 0;
     for( auto const* mct : fConvertList){
     
       // For each simb::MCParticle in the MCTruth:
       for ( int p = 0; p != mct->NParticles(); ++p ){
         // Implementation note: the following statement copies the
         // MCParticle from the MCTruth, instead of getting a
         // const reference.    
         simb::MCParticle const& particle = mct->GetParticle(p);
 
         // status code == 1 means "track this particle."  Any
         // other status code should be ignored by the Monte Carlo.
         if ( particle.StatusCode() != 1 ) continue;
       
         // Get the Particle Data Group code for the particle.
         G4int pdgCode = particle.PdgCode();
       
         // Get the vertex.  Note that LArSoft/ROOT uses cm, but
         // Geant4/CLHEP uses mm.
         G4double x = particle.Vx() * CLHEP::cm;
         G4double y = particle.Vy() * CLHEP::cm;
         G4double z = particle.Vz() * CLHEP::cm;
         G4double t = particle.T()  * CLHEP::ns;
       
         // Create a CLHEP four-vector from the particle's vertex.
         CLHEP::HepLorentzVector fourpos(x,y,z,t);
       
         // Is this vertex already in our map?
         G4PrimaryVertex* vertex = 0;
         std::map< CLHEP::HepLorentzVector, G4PrimaryVertex* >::const_iterator result = vertexMap.find( fourpos );
         if ( result == vertexMap.end() ){
           // No, it's not, so create a new vertex and add it to the
           // map.
           vertex = new G4PrimaryVertex(x, y, z, t);
           vertexMap[ fourpos ] = vertex;
 
           // Add the vertex to the G4Event.
           event->AddPrimaryVertex( vertex );
         }
         else{
           // Yes, it is, so use the existing vertex.
           vertex = (*result).second;
         }
       
         // Get additional particle information.
         TLorentzVector momentum = particle.Momentum(); // (px,py,pz,E)
         TVector3 polarization = particle.Polarization();
       
         // Get the particle table if necessary.  (Note: we're
         // doing this "late" because I'm not sure at what point
         // the G4 particle table is initialized in the loading process.
         if ( fParticleTable == 0 ){
           fParticleTable = G4ParticleTable::GetParticleTable();
         }
 
         // Get Geant4's definition of the particle.
         G4ParticleDefinition* particleDefinition;
       
         if(pdgCode==0){
           particleDefinition = fParticleTable->FindParticle("opticalphoton");
         }
         else
           particleDefinition = fParticleTable->FindParticle(pdgCode);
 
         if ( pdgCode > 1000000000) { // If the particle is a nucleus
           MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% Nuclear PDG code = " << pdgCode
                                               << " (x,y,z,t)=(" << x
                                               << "," << y
                                               << "," << z
                                               << "," << t << ")"
                                               << " P=" << momentum.P()
                                               << ", E=" << momentum.E();
           // If the particle table doesn't have a definition yet, ask the ion
           // table for one. This will create a new ion definition as needed.
           if (!particleDefinition) {
             int Z = (pdgCode % 10000000) / 10000; // atomic number
             int A = (pdgCode % 10000) / 10; // mass number
             particleDefinition = fParticleTable->GetIonTable()->GetIon(Z, A, 0.);
           }
         }
 
         // What if the PDG code is unknown?  This has been a known
         // issue with GENIE.
         if ( particleDefinition == 0 ){
           MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% Code not found = " << pdgCode;
           fUnknownPDG[ pdgCode ] += 1;
           continue;
         }
       
         // Create a Geant4 particle to add to the vertex.
         G4PrimaryParticle* g4particle = new G4PrimaryParticle( particleDefinition,
                                                                momentum.Px() * CLHEP::GeV,
                                                                momentum.Py() * CLHEP::GeV,
                                                                momentum.Pz() * CLHEP::GeV);
 
         // Add more particle information the Geant4 particle.
         G4double charge = particleDefinition->GetPDGCharge();
         g4particle->SetCharge( charge );
         g4particle->SetPolarization( polarization.x(),
                                      polarization.y(),
                                      polarization.z() );
       
         // Add the particle to the vertex.
         vertex->SetPrimary( g4particle );
 
         // Create a PrimaryParticleInformation object, and save
         // the MCTruth pointer in it.  This will allow the
         // ParticleActionList class to access MCTruth
         // information during Geant4's tracking.
         PrimaryParticleInformation* primaryParticleInfo = new PrimaryParticleInformation;
         primaryParticleInfo->SetMCTruth( mct, index, p );
           
         // Save the PrimaryParticleInformation in the
         // G4PrimaryParticle for access during tracking.
         g4particle->SetUserInformation( primaryParticleInfo );
 
         MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% primary PDG=" << pdgCode
                                             << ", (x,y,z,t)=(" << x
                                             << "," << y
                                             << "," << z
                                             << "," << t << ")"
                                             << " P=" << momentum.P()
                                             << ", E=" << momentum.E();
 
       } // for each particle in MCTruth
       ++index;
     } // for each MCTruth
   }// GeneratePrimaries