diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index 0a266b801be9532b1d51ac35498af790a48f86c8..a9d739a68e1fa8bb16f8211ea46c0c0ba711e25e 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -2,6 +2,9 @@ #include "G4Event.hh" #include "G4THitsCollection.hh" +#include "G4EventManager.hh" +#include "G4TrackingManager.hh" +#include "G4SteppingManager.hh" #include "DD4hep/Detector.h" #include "DD4hep/Plugins.h" @@ -204,6 +207,11 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { } edm_trk_hit.setMCParticle(mcCol->at(pritrkid-1)); + + if (pritrkid != trackID) { + // If the track is a secondary, then the primary track id and current track id is different + edm_trk_hit.setProducedBySecondary(true); + } } dd4hep::sim::Geant4CalorimeterHit* cal_hit = dynamic_cast<dd4hep::sim::Geant4CalorimeterHit*>(h); @@ -311,6 +319,90 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) { stop_mom.z()/CLHEP::GeV); primary_particle.setMomentumAtEndpoint(mom_endpoint); + // =================================================================== + // Update the flags of the primary particles + // =================================================================== + + // processes + static G4VProcess* proc_decay = nullptr; + if (!proc_decay) { + G4ProcessManager* pm + = track->GetDefinition()->GetProcessManager(); + G4int nprocesses = pm->GetProcessListLength(); + G4ProcessVector* pv = pm->GetProcessList(); + + for(G4int i=0; i<nprocesses; ++i){ + if((*pv)[i]->GetProcessName()=="Decay"){ + proc_decay = (*pv)[i]; + } + } + + } + + // flags + bool is_decay = false; + + G4TrackingManager* tm = G4EventManager::GetEventManager() + -> GetTrackingManager(); + G4TrackVector* secondaries = tm->GimmeSecondaries(); + if(secondaries) { + + + size_t nSeco = secondaries->size(); + for (size_t i = 0; i < nSeco; ++i) { + G4Track* sectrk = (*secondaries)[i]; + G4ParticleDefinition* secparticle = sectrk->GetDefinition(); + const G4VProcess* creatorProcess = sectrk->GetCreatorProcess(); + + // select the necessary processes + if (creatorProcess==proc_decay) { + info() << "Creator Process is Decay for secondary particle: " + << " idx: " << i + << " trkid: " << sectrk->GetTrackID() // not valid until track + << " particle: " << secparticle->GetParticleName() + << " pdg: " << secparticle->GetPDGEncoding() + << " at position: " << sectrk->GetPosition() // + << " time: " << sectrk->GetGlobalTime() + << " momentum: " << sectrk->GetMomentum() // + << endmsg; + is_decay = true; + + // create secondaries in MC particles + // todo: convert the const collection to non-const + auto mcCol = const_cast<edm4hep::MCParticleCollection*>(m_mcParCol.get()); + edm4hep::MCParticle mcp = mcCol->create(); + mcp.setPDG(secparticle->GetPDGEncoding()); + mcp.setGeneratorStatus(0); // not created by Generator + mcp.setCreatedInSimulation(1); + mcp.setCharge(secparticle->GetPDGCharge()); + mcp.setTime(sectrk->GetGlobalTime()/CLHEP::ns); // todo + mcp.setMass(secparticle->GetPDGMass()); + + const G4ThreeVector& sec_init_pos = sectrk->GetPosition(); + double x=sec_init_pos.x()/CLHEP::mm; + double y=sec_init_pos.y()/CLHEP::mm; + double z=sec_init_pos.z()/CLHEP::mm; + + const G4ThreeVector& sec_init_mom = sectrk->GetMomentum(); + double px=sec_init_mom.x()/CLHEP::GeV; + double py=sec_init_mom.y()/CLHEP::GeV; + double pz=sec_init_mom.z()/CLHEP::GeV; + mcp.setVertex(edm4hep::Vector3d(x,y,z)); // todo + mcp.setEndpoint(edm4hep::Vector3d(x,y,z)); // todo + mcp.setMomentum(edm4hep::Vector3f(px,py,pz)); // todo + mcp.setMomentumAtEndpoint(edm4hep::Vector3f(px,py,pz)); //todo + + mcp.addToParents(primary_particle); + primary_particle.addToDaughters(mcp); + } + } + } + + // now update + if (is_decay) { + primary_particle.setDecayedInTracker(is_decay); + primary_particle.setDecayedInCalorimeter(is_decay); + } } else { // TODO: select other interested tracks