File indexing completed on 2025-01-30 09:17:23
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #ifndef DD4HEP_DDG4_GEANT4OUTPUT2LCIO_H
0015 #define DD4HEP_DDG4_GEANT4OUTPUT2LCIO_H
0016
0017
0018 #include <DD4hep/VolumeManager.h>
0019 #include <DDG4/Geant4OutputAction.h>
0020
0021 #include <DDG4/EventParameters.h>
0022 #include <DDG4/RunParameters.h>
0023
0024 #include <G4Threading.hh>
0025 #include <G4AutoLock.hh>
0026
0027 #include <DD4hep/Detector.h>
0028 #include <G4Version.hh>
0029
0030
0031 #include <lcio.h>
0032 #include <IO/LCWriter.h>
0033 #include <IMPL/LCEventImpl.h>
0034 #include <IMPL/LCCollectionVec.h>
0035 #include <EVENT/LCParameters.h>
0036
0037 using namespace lcio ;
0038
0039
0040 namespace dd4hep {
0041
0042 class ComponentCast;
0043
0044
0045 namespace sim {
0046 template <class T=lcio::LCEventImpl> void EventParameters::extractParameters(T& event){
0047 auto& lcparameters = event.parameters();
0048 event.setRunNumber(this->runNumber());
0049 event.setEventNumber(this->eventNumber());
0050 for(auto const& ival: this->intParameters()) {
0051 lcparameters.setValues(ival.first, ival.second);
0052 }
0053 for(auto const& ival: this->fltParameters()) {
0054 lcparameters.setValues(ival.first, ival.second);
0055 }
0056 for(auto const& ival: this->strParameters()) {
0057 lcparameters.setValues(ival.first, ival.second);
0058 }
0059 #if LCIO_VERSION_GE(2, 17)
0060 for(auto const& ival: this->dblParameters()) {
0061 lcparameters.setValues(ival.first, ival.second);
0062 }
0063 #endif
0064 }
0065
0066 template <class T=lcio::LCRunHeaderImpl> void RunParameters::extractParameters(T& runHeader){
0067 auto& lcparameters = runHeader.parameters();
0068 for(auto const& ival: this->intParameters()) {
0069 lcparameters.setValues(ival.first, ival.second);
0070 }
0071 for(auto const& ival: this->fltParameters()) {
0072 lcparameters.setValues(ival.first, ival.second);
0073 }
0074 for(auto const& ival: this->strParameters()) {
0075 lcparameters.setValues(ival.first, ival.second);
0076 }
0077 #if LCIO_VERSION_GE(2, 17)
0078 for(auto const& ival: this->dblParameters()) {
0079 lcparameters.setValues(ival.first, ival.second);
0080 }
0081 #endif
0082 }
0083
0084 class Geant4ParticleMap;
0085
0086
0087
0088
0089
0090
0091
0092
0093 class Geant4Output2LCIO : public Geant4OutputAction {
0094 protected:
0095 lcio::LCWriter* m_file;
0096 int m_runNo;
0097 int m_runNumberOffset;
0098 int m_eventNumberOffset;
0099 std::map< std::string, std::string > m_runHeader;
0100 std::map< std::string, std::string > m_eventParametersInt;
0101 std::map< std::string, std::string > m_eventParametersFloat;
0102 std::map< std::string, std::string > m_eventParametersString;
0103
0104
0105 lcio::LCCollectionVec* saveParticles(Geant4ParticleMap* particles);
0106 public:
0107
0108 Geant4Output2LCIO(Geant4Context* ctxt, const std::string& nam);
0109
0110 virtual ~Geant4Output2LCIO();
0111
0112 virtual void beginRun(const G4Run* run);
0113
0114 virtual void endRun(const G4Run* run);
0115
0116
0117 virtual void saveRun(const G4Run* run);
0118
0119 virtual void saveEvent( OutputContext<G4Event>& ctxt);
0120
0121 virtual void saveCollection( OutputContext<G4Event>& ctxt, G4VHitsCollection* collection);
0122
0123 virtual void commit( OutputContext<G4Event>& ctxt);
0124
0125
0126 virtual void begin(const G4Event* event);
0127 protected:
0128
0129 template <typename T>
0130 void saveEventParameters(lcio::LCEventImpl* event, const std::map<std::string, std::string >& parameters);
0131 };
0132
0133
0134 template <typename T>
0135 inline void Geant4Output2LCIO::saveEventParameters(lcio::LCEventImpl* event, const std::map<std::string, std::string >& parameters) {
0136 for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
0137 T parameter;
0138 std::istringstream iss(iter->second);
0139 if ( (iss >> parameter).fail() ) {
0140 printout(FATAL,"saveEventParameters","+++ Event parameter %s: FAILED to convert to type :%s",iter->first.c_str(),typeid(T).name());
0141 continue;
0142 }
0143 event->parameters().setValue(iter->first,parameter);
0144 }
0145 }
0146
0147
0148 template <>
0149 inline void Geant4Output2LCIO::saveEventParameters<std::string>(lcio::LCEventImpl* event, const std::map<std::string, std::string >& parameters) {
0150 for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
0151 event->parameters().setValue(iter->first,iter->second);
0152 }
0153 }
0154
0155 }
0156 }
0157 #endif
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173 #include <DD4hep/InstanceCount.h>
0174 #include <DD4hep/Detector.h>
0175 #include <DDG4/Geant4HitCollection.h>
0176 #include <DDG4/Geant4DataConversion.h>
0177 #include <DDG4/Geant4Context.h>
0178 #include <DDG4/Geant4Particle.h>
0179 #include <DDG4/Geant4Data.h>
0180 #include <DDG4/Geant4Action.h>
0181
0182
0183 #include <G4ParticleDefinition.hh>
0184 #include <G4VProcess.hh>
0185 #include <G4Event.hh>
0186 #include <G4Run.hh>
0187
0188
0189 #include <IMPL/LCEventImpl.h>
0190 #include <IMPL/LCRunHeaderImpl.h>
0191 #include <IMPL/LCCollectionVec.h>
0192 #include <IMPL/ClusterImpl.h>
0193 #include <IMPL/SimTrackerHitImpl.h>
0194 #include <IMPL/SimCalorimeterHitImpl.h>
0195 #include <IMPL/MCParticleImpl.h>
0196 #include <UTIL/ILDConf.h>
0197
0198 using namespace dd4hep::sim;
0199 using namespace dd4hep;
0200 using namespace std;
0201 namespace {
0202 G4Mutex action_mutex=G4MUTEX_INITIALIZER;
0203 }
0204
0205 #include <DDG4/Factories.h>
0206 DECLARE_GEANT4ACTION(Geant4Output2LCIO)
0207
0208
0209 Geant4Output2LCIO::Geant4Output2LCIO(Geant4Context* ctxt, const string& nam)
0210 : Geant4OutputAction(ctxt,nam), m_file(0), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
0211 {
0212 declareProperty("RunHeader", m_runHeader);
0213 declareProperty("EventParametersInt", m_eventParametersInt);
0214 declareProperty("EventParametersFloat", m_eventParametersFloat);
0215 declareProperty("EventParametersString", m_eventParametersString);
0216 declareProperty("RunNumberOffset", m_runNumberOffset);
0217 declareProperty("EventNumberOffset", m_eventNumberOffset);
0218 InstanceCount::increment(this);
0219 }
0220
0221
0222 Geant4Output2LCIO::~Geant4Output2LCIO() {
0223 G4AutoLock protection_lock(&action_mutex);
0224 if ( m_file ) {
0225 m_file->close();
0226 detail::deletePtr(m_file);
0227 }
0228 InstanceCount::decrement(this);
0229 }
0230
0231
0232 void Geant4Output2LCIO::beginRun(const G4Run* run) {
0233 if ( 0 == m_file && !m_output.empty() ) {
0234 G4AutoLock protection_lock(&action_mutex);
0235 m_file = lcio::LCFactory::getInstance()->createLCWriter();
0236 m_file->open(m_output,lcio::LCIO::WRITE_NEW);
0237 }
0238
0239 saveRun(run);
0240 }
0241
0242
0243 void Geant4Output2LCIO::endRun(const G4Run* ) {
0244
0245 }
0246
0247
0248 void Geant4Output2LCIO::commit( OutputContext<G4Event>& ) {
0249 lcio::LCEventImpl* e = context()->event().extension<lcio::LCEventImpl>();
0250 if ( m_file ) {
0251 G4AutoLock protection_lock(&action_mutex);
0252 m_file->writeEvent(e);
0253 return;
0254 }
0255 except("+++ Failed to write output file. [Stream is not open]");
0256 }
0257
0258
0259 void Geant4Output2LCIO::saveRun(const G4Run* run) {
0260 G4AutoLock protection_lock(&action_mutex);
0261
0262 lcio::LCRunHeaderImpl* rh = new lcio::LCRunHeaderImpl;
0263 for (std::map< std::string, std::string >::iterator it = m_runHeader.begin(); it != m_runHeader.end(); ++it) {
0264 rh->parameters().setValue( it->first, it->second );
0265 }
0266 m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
0267 rh->parameters().setValue("GEANT4Version", G4Version);
0268 rh->parameters().setValue("DD4HEPVersion", versionString());
0269 rh->setRunNumber(m_runNo);
0270 rh->setDetectorName(context()->detectorDescription().header().name());
0271 auto* parameters = context()->run().extension<RunParameters>(false);
0272 if (parameters) {
0273 parameters->extractParameters(*rh);
0274 }
0275 m_file->writeRunHeader(rh);
0276 }
0277
0278 void Geant4Output2LCIO::begin(const G4Event* ) {
0279 lcio::LCEventImpl* e = new lcio::LCEventImpl;
0280
0281
0282 context()->event().addExtension<lcio::LCEventImpl>( e );
0283 }
0284
0285
0286 lcio::LCCollectionVec* Geant4Output2LCIO::saveParticles(Geant4ParticleMap* particles) {
0287 typedef detail::ReferenceBitMask<const int> PropertyMask;
0288 typedef Geant4ParticleMap::ParticleMap ParticleMap;
0289 const ParticleMap& pm = particles->particleMap;
0290 size_t nparts = pm.size();
0291 lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::MCPARTICLE);
0292 lc_coll->reserve(nparts);
0293 if ( nparts > 0 ) {
0294 size_t cnt = 0;
0295 map<int,int> p_ids;
0296 vector<const Geant4Particle*> p_part(pm.size(),0);
0297 vector<MCParticleImpl*> p_lcio(pm.size(),0);
0298
0299 for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end();++i, ++cnt) {
0300 int id = (*i).first;
0301 const Geant4ParticleHandle p = (*i).second;
0302 PropertyMask mask(p->status);
0303
0304 const G4ParticleDefinition* def = p.definition();
0305 MCParticleImpl* q = new lcio::MCParticleImpl();
0306 q->setPDG(p->pdgID);
0307
0308 float ps_fa[3] = {float(p->psx/CLHEP::GeV),float(p->psy/CLHEP::GeV),float(p->psz/CLHEP::GeV)};
0309 q->setMomentum( ps_fa );
0310
0311 #if LCIO_VERSION_GE(2,7)
0312 float pe_fa[3] = {float(p->pex/CLHEP::GeV),float(p->pey/CLHEP::GeV),float(p->pez/CLHEP::GeV)};
0313 q->setMomentumAtEndpoint( pe_fa );
0314 #endif
0315 double vs_fa[3] = { p->vsx/CLHEP::mm, p->vsy/CLHEP::mm, p->vsz/CLHEP::mm } ;
0316 q->setVertex( vs_fa );
0317
0318 double ve_fa[3] = { p->vex/CLHEP::mm, p->vey/CLHEP::mm, p->vez/CLHEP::mm } ;
0319 q->setEndpoint( ve_fa );
0320
0321 q->setTime(p->time/CLHEP::ns);
0322 q->setMass(p->mass/CLHEP::GeV);
0323 q->setCharge(def ? def->GetPDGCharge() : 0);
0324
0325
0326 q->setGeneratorStatus(0);
0327 if( p->genStatus ) {
0328 q->setGeneratorStatus( p->genStatus ) ;
0329 } else {
0330
0331 if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) q->setGeneratorStatus(1);
0332 else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) q->setGeneratorStatus(2);
0333 else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) q->setGeneratorStatus(3);
0334 else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) q->setGeneratorStatus(4);
0335 else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) q->setGeneratorStatus(9);
0336 }
0337
0338
0339
0340 q->setCreatedInSimulation( mask.isSet(G4PARTICLE_SIM_CREATED) );
0341 q->setBackscatter( mask.isSet(G4PARTICLE_SIM_BACKSCATTER) );
0342 q->setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) );
0343 q->setDecayedInTracker( mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) );
0344 q->setDecayedInCalorimeter( mask.isSet(G4PARTICLE_SIM_DECAY_CALO) );
0345 q->setHasLeftDetector( mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) );
0346 q->setStopped( mask.isSet(G4PARTICLE_SIM_STOPPED) );
0347 q->setOverlay( false );
0348
0349
0350 if( q->isCreatedInSimulation() )
0351 q->setGeneratorStatus( 0 ) ;
0352
0353 q->setSpin(p->spin);
0354 q->setColorFlow(p->colorFlow);
0355
0356 lc_coll->addElement(q);
0357 p_ids[id] = cnt;
0358 p_part[cnt] = p;
0359 p_lcio[cnt] = q;
0360 }
0361
0362
0363 for(size_t i=0, n=p_ids.size(); i<n; ++i) {
0364 map<int,int>::iterator k;
0365 const Geant4Particle* p = p_part[i];
0366 MCParticleImpl* q = p_lcio[i];
0367 const Geant4Particle::Particles& dau = p->daughters;
0368 for( Geant4Particle::Particles::const_iterator j=dau.begin(); j != dau.end(); ++j ) {
0369 int idau = *j;
0370 if ( (k=p_ids.find(idau)) == p_ids.end() ) {
0371 printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
0372 continue;
0373 }
0374 int iqdau = (*k).second;
0375 MCParticleImpl* qdau = p_lcio[iqdau];
0376 qdau->addParent(q);
0377 }
0378 const Geant4Particle::Particles& par = p->parents;
0379 for( Geant4Particle::Particles::const_iterator j=par.begin(); j != par.end(); ++j ) {
0380 int ipar = *j;
0381 if ( ipar >= 0 ) {
0382 if( (k=p_ids.find(ipar)) == p_ids.end() ) {
0383 printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
0384 continue;
0385 }
0386 int iqpar = (*k).second;
0387 MCParticleImpl* qpar = p_lcio[iqpar];
0388 q->addParent(qpar);
0389 }
0390 }
0391 }
0392 }
0393 return lc_coll;
0394 }
0395
0396
0397 void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) {
0398 lcio::LCEventImpl* e = context()->event().extension<lcio::LCEventImpl>();
0399 EventParameters* parameters = context()->event().extension<EventParameters>(false);
0400 int runNumber(0), eventNumber(0);
0401 const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0);
0402 const int runNumberOffset(m_runNumberOffset > 0 ? m_runNumberOffset : 0);
0403 double eventWeight{0};
0404
0405 if (parameters) {
0406 runNumber = parameters->runNumber() + runNumberOffset;
0407 eventNumber = parameters->eventNumber() + eventNumberOffset;
0408 parameters->extractParameters(*e);
0409 #if LCIO_VERSION_GE(2, 17)
0410 eventWeight = e->getParameters().getDoubleVal("EventWeights");
0411 #endif
0412 } else {
0413 runNumber = m_runNo + runNumberOffset;
0414 eventNumber = ctxt.context->GetEventID() + eventNumberOffset;
0415 }
0416 print("+++ Saving LCIO event %d run %d ....", eventNumber, runNumber);
0417 e->setRunNumber(runNumber);
0418 e->setEventNumber(eventNumber);
0419 e->setWeight(eventWeight);
0420 e->setDetectorName(context()->detectorDescription().header().name());
0421 saveEventParameters<int>(e, m_eventParametersInt);
0422 saveEventParameters<float>(e, m_eventParametersFloat);
0423 saveEventParameters<std::string>(e, m_eventParametersString);
0424 lcio::LCEventImpl* evt = context()->event().extension<lcio::LCEventImpl>();
0425 Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>(false);
0426 if ( part_map ) {
0427 print("+++ Saving %d LCIO particles....",int(part_map->particleMap.size()));
0428 if ( part_map->particleMap.size() > 0 ) {
0429 lcio::LCCollectionVec* col = saveParticles(part_map);
0430 evt->addCollection(col,lcio::LCIO::MCPARTICLE);
0431 }
0432 }
0433 }
0434
0435
0436 void Geant4Output2LCIO::saveCollection(OutputContext<G4Event>& , G4VHitsCollection* collection) {
0437 size_t nhits = collection->GetSize();
0438 std::string hc_nam = collection->GetName();
0439 print("+++ Saving LCIO collection %s with %d entries....",hc_nam.c_str(),int(nhits));
0440 typedef pair<const Geant4Context*,G4VHitsCollection*> _Args;
0441 typedef Geant4Conversion<lcio::LCCollectionVec,_Args> _C;
0442 const _C& cnv = _C::converter(typeid(Geant4HitCollection));
0443 cnv(_Args(context(),collection));
0444 }
0445
