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   EVENT::LCCollection*        primaries = 0;
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 ) return ret;
0072
0073   int NHEP = primaries->getNumberOfElements();
0074   // check if there is at least one particle
0075   if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
0076
0077   mcpcoll.resize(NHEP,0);
0078   for(int i=0; i<NHEP; ++i ) {
0079     EVENT::MCParticle* p = dynamic_cast<EVENT::MCParticle*>(primaries->getElementAt(i));
0080     mcparts[p] = i;
0081     mcpcoll[i] = p;
0082   }
0083
0084   // build collection of MCParticles
0085   for(int i=0; i<NHEP; ++i )   {
0086     EVENT::MCParticle* mcp = mcpcoll[i];
0087     Geant4ParticleHandle p(new Particle(i));
0088     const double *mom   = mcp->getMomentum();
0089     const double *vsx   = mcp->getVertex();
0090     const double *vex   = mcp->getEndpoint();
0091     const float  *spin  = mcp->getSpin();
0092     const int    *color = mcp->getColorFlow();
0093     const int     pdg   = mcp->getPDG();
0094     p->pdgID        = pdg;
0095     p->charge       = int(mcp->getCharge()*3.0);
0096     p->psx          = mom[0]*CLHEP::GeV;
0097     p->psy          = mom[1]*CLHEP::GeV;
0098     p->psz          = mom[2]*CLHEP::GeV;
0099     p->time         = mcp->getTime()*CLHEP::ns;
0100     p->properTime   = mcp->getTime()*CLHEP::ns;
0101     p->vsx          = vsx[0]*CLHEP::mm;
0102     p->vsy          = vsx[1]*CLHEP::mm;
0103     p->vsz          = vsx[2]*CLHEP::mm;
0104     p->vex          = vex[0]*CLHEP::mm;
0105     p->vey          = vex[1]*CLHEP::mm;
0106     p->vez          = vex[2]*CLHEP::mm;
0107     p->process      = 0;
0108     p->spin[0]      = spin[0];
0109     p->spin[1]      = spin[1];
0110     p->spin[2]      = spin[2];
0111     p->colorFlow[0] = color[0];
0112     p->colorFlow[1] = color[1];
0113     p->mass         = mcp->getMass()*CLHEP::GeV;
0114     const EVENT::MCParticleVec &par = mcp->getParents(), &dau=mcp->getDaughters();
0115     for(int num=dau.size(),k=0; k<num; ++k)
0116       p->daughters.insert(GET_ENTRY(mcparts,dau[k]));
0117     for(int num=par.size(),k=0; k<num; ++k)
0118       p->parents.insert(GET_ENTRY(mcparts,par[k]));
0119
0120     PropertyMask status(p->status);
0121     int genStatus = mcp->getGeneratorStatus();
0122     // Copy raw generator status
0123     p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
0124     m_inputAction->setGeneratorStatus(genStatus, status);
0125
0126     //fg: we simply add all particles without parents as with their own vertex.
0127     //    This might include the incoming beam particles, e.g. in
0128     //    the case of lcio files written with Whizard2, which is slightly odd,
0129     //    however should be treated correctly in Geant4InputHandling.cpp.
0130     //    We no longer make an attempt to identify the incoming particles
0131     //    based on the generator status, as this varies widely with different
0132     //    generators.
0133
0134     if ( p->parents.size() == 0 )  {
0135
0136       Geant4Vertex* vtx = new Geant4Vertex ;
0137       vertices.emplace_back( vtx );
0138       vtx->x = p->vsx;
0139       vtx->y = p->vsy;
0140       vtx->z = p->vsz;
0141       vtx->time = p->time;
0142
0143       vtx->out.insert(p->id) ;
0144     }
0145
0146     if ( mcp->isCreatedInSimulation() )       status.set(G4PARTICLE_SIM_CREATED);
0147     if ( mcp->isBackscatter() )               status.set(G4PARTICLE_SIM_BACKSCATTER);
0148     if ( mcp->vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED);
0149     if ( mcp->isDecayedInTracker() )          status.set(G4PARTICLE_SIM_DECAY_TRACKER);
0150     if ( mcp->isDecayedInCalorimeter() )      status.set(G4PARTICLE_SIM_DECAY_CALO);
0151     if ( mcp->hasLeftDetector() )             status.set(G4PARTICLE_SIM_LEFT_DETECTOR);
0152     if ( mcp->isStopped() )                   status.set(G4PARTICLE_SIM_STOPPED);
0153     if ( mcp->isOverlay() )                   status.set(G4PARTICLE_SIM_OVERLAY);
0154     particles.emplace_back(p);
0155   }
0156   return EVENT_READER_OK;
0157 }
0158