DDG4/lcio/LCIOEventReader.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // @author  P.Kostka (main author)
0011 // @author  M.Frank  (code reshuffeling into new DDG4 scheme)
0012 //
0013 //====================================================================
0014
0015 // Framework include files
0016 #include "LCIOEventReader.h"
0017 #include <DD4hep/Printout.h>
0018 #include <DDG4/Geant4Primary.h>
0019 #include <DDG4/Geant4Context.h>
0020 #include <DDG4/Factories.h>
0021
0022 #include <G4ParticleTable.hh>
0023 #include <EVENT/MCParticle.h>
0024 #include <EVENT/LCCollection.h>
0025
0026 #include <G4Event.hh>
0027
0028 using namespace std;
0029 using namespace dd4hep;
0030 using namespace dd4hep::sim;
0031 typedef dd4hep::detail::ReferenceBitMask<int> PropertyMask;
0032
0033 // Neede for backwards compatibility:
0034 namespace dd4hep{namespace sim{typedef Geant4InputAction LCIOInputAction;}}
0035 DECLARE_GEANT4ACTION(LCIOInputAction)
0036
0037
0038 namespace {
0039   inline int GET_ENTRY(const map<EVENT::MCParticle*,int>& mcparts, EVENT::MCParticle* part)  {
0040     map<EVENT::MCParticle*,int>::const_iterator ip=mcparts.find(part);
0041     if ( ip == mcparts.end() )  {
0042       throw runtime_error("Unknown particle identifier look-up!");
0043     }
0044     return (*ip).second;
0045   }
0046 }
0047
0048
0049 /// Initializing constructor
0050 LCIOEventReader::LCIOEventReader(const string& nam)
0051   : Geant4EventReader(nam)
0052 {
0053 }
0054
0055 /// Default destructor
0056 LCIOEventReader::~LCIOEventReader()   {
0057 }
0058
0059
0060 /// Read an event and fill a vector of MCParticles.
0061 LCIOEventReader::EventReaderStatus
0062 LCIOEventReader::readParticles(int event_number,
0063                                Vertices& vertices,
0064                                vector<Particle*>& particles)
0065 {
0066   CollectionOwner primaries(nullptr, [](EVENT::LCCollection*){});
0067   map<EVENT::MCParticle*,int> mcparts;
0068   vector<EVENT::MCParticle*>  mcpcoll;
0069   EventReaderStatus ret = readParticleCollection(event_number, primaries);
0070
0071   if (ret != EVENT_READER_OK) {
0072     return ret;
0073   }
0074   int NHEP = primaries->getNumberOfElements();
0075   // check if there is at least one particle
0076   if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
0077
0078   mcpcoll.resize(NHEP,0);
0079   for(int i=0; i<NHEP; ++i ) {
0080     EVENT::MCParticle* p = dynamic_cast<EVENT::MCParticle*>(primaries->getElementAt(i));
0081     mcparts[p] = i;
0082     mcpcoll[i] = p;
0083   }
0084
0085   // build collection of MCParticles
0086   for(int i=0; i<NHEP; ++i )   {
0087     EVENT::MCParticle* mcp = mcpcoll[i];
0088     Geant4ParticleHandle p(new Particle(i));
0089     const double *mom   = mcp->getMomentum();
0090     const double *vsx   = mcp->getVertex();
0091     const double *vex   = mcp->getEndpoint();
0092     const float  *spin  = mcp->getSpin();
0093     const int    *color = mcp->getColorFlow();
0094     const int     pdg   = mcp->getPDG();
0095     p->pdgID        = pdg;
0096     p->charge       = int(mcp->getCharge()*3.0);
0097     p->psx          = mom[0]*CLHEP::GeV;
0098     p->psy          = mom[1]*CLHEP::GeV;
0099     p->psz          = mom[2]*CLHEP::GeV;
0100     p->time         = mcp->getTime()*CLHEP::ns;
0101     p->properTime   = mcp->getTime()*CLHEP::ns;
0102     p->vsx          = vsx[0]*CLHEP::mm;
0103     p->vsy          = vsx[1]*CLHEP::mm;
0104     p->vsz          = vsx[2]*CLHEP::mm;
0105     p->vex          = vex[0]*CLHEP::mm;
0106     p->vey          = vex[1]*CLHEP::mm;
0107     p->vez          = vex[2]*CLHEP::mm;
0108     p->process      = 0;
0109     p->spin[0]      = spin[0];
0110     p->spin[1]      = spin[1];
0111     p->spin[2]      = spin[2];
0112     p->colorFlow[0] = color[0];
0113     p->colorFlow[1] = color[1];
0114     p->mass         = mcp->getMass()*CLHEP::GeV;
0115     const EVENT::MCParticleVec &par = mcp->getParents(), &dau=mcp->getDaughters();
0116     for(int num=dau.size(),k=0; k<num; ++k)
0117       p->daughters.insert(GET_ENTRY(mcparts,dau[k]));
0118     for(int num=par.size(),k=0; k<num; ++k)
0119       p->parents.insert(GET_ENTRY(mcparts,par[k]));
0120
0121     PropertyMask status(p->status);
0122     int genStatus = mcp->getGeneratorStatus();
0123     // Copy raw generator status
0124     p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
0125     if(m_inputAction) {
0126       // in some tests we do not set up the inputAction
0127       m_inputAction->setGeneratorStatus(genStatus, status);
0128     }
0129
0130     //fg: we simply add all particles without parents as with their own vertex.
0131     //    This might include the incoming beam particles, e.g. in
0132     //    the case of lcio files written with Whizard2, which is slightly odd,
0133     //    however should be treated correctly in Geant4InputHandling.cpp.
0134     //    We no longer make an attempt to identify the incoming particles
0135     //    based on the generator status, as this varies widely with different
0136     //    generators.
0137
0138     if ( p->parents.size() == 0 )  {
0139
0140       Geant4Vertex* vtx = new Geant4Vertex ;
0141       vertices.emplace_back( vtx );
0142       vtx->x = p->vsx;
0143       vtx->y = p->vsy;
0144       vtx->z = p->vsz;
0145       vtx->time = p->time;
0146
0147       vtx->out.insert(p->id) ;
0148     }
0149
0150     if ( mcp->isCreatedInSimulation() )       status.set(G4PARTICLE_SIM_CREATED);
0151     if ( mcp->isBackscatter() )               status.set(G4PARTICLE_SIM_BACKSCATTER);
0152     if ( mcp->vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED);
0153     if ( mcp->isDecayedInTracker() )          status.set(G4PARTICLE_SIM_DECAY_TRACKER);
0154     if ( mcp->isDecayedInCalorimeter() )      status.set(G4PARTICLE_SIM_DECAY_CALO);
0155     if ( mcp->hasLeftDetector() )             status.set(G4PARTICLE_SIM_LEFT_DETECTOR);
0156     if ( mcp->isStopped() )                   status.set(G4PARTICLE_SIM_STOPPED);
0157     if ( mcp->isOverlay() )                   status.set(G4PARTICLE_SIM_OVERLAY);
0158     particles.emplace_back(p);
0159   }
0160   return EVENT_READER_OK;
0161 }
0162