From 1f3f34eb5a702bfcbc6f6824305071262e5573e4 Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Fri, 21 Nov 2014 16:49:22 +0000
Subject: [PATCH] Fix problem with tracking Geantinos through Geant4 --
reported by Nikiforos
---
DDG4/examples/CLICSidSimuMarkus.py | 14 +++++---
DDG4/include/DDG4/Geant4Particle.h | 12 +++++--
DDG4/lcio/Geant4Output2LCIO.cpp | 4 +--
DDG4/lcio/LCIOEventReader.cpp | 8 +++--
DDG4/src/Geant4InputHandling.cpp | 14 ++++++--
DDG4/src/Geant4IsotropeGenerator.cpp | 2 +-
DDG4/src/Geant4Particle.cpp | 53 +++++++++++++++++-----------
DDG4/src/Geant4ParticleGun.cpp | 4 +--
DDG4/src/Geant4ParticleHandler.cpp | 10 +++---
9 files changed, 77 insertions(+), 44 deletions(-)
diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py
index 653a2fb67..b5dea3f23 100644
--- a/DDG4/examples/CLICSidSimuMarkus.py
+++ b/DDG4/examples/CLICSidSimuMarkus.py
@@ -41,11 +41,12 @@ def run():
generator_output_level = Output.DEBUG
# Configure I/O
- evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
- evt_lcio.OutputLevel = generator_output_level
- evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
+ ##evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
+ ##evt_lcio.OutputLevel = generator_output_level
+ ##evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
gen = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit")
+ gen.OutputLevel = generator_output_level
kernel.generatorAction().adopt(gen)
#VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
@@ -65,7 +66,7 @@ def run():
#gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt"
#gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/sherpa-2.1.1_zjets.hepmc2g"
gen.Input = "LCIOFileReader|/afs/cern.ch/user/n/nikiforo/public/Markus/muons.slcio"
- gen.Input = "LCIOFileReader|/afs/cern.ch/user/n/nikiforo/public/Markus/geantinos.slcio"
+ #gen.Input = "LCIOFileReader|/afs/cern.ch/user/n/nikiforo/public/Markus/geantinos.slcio"
gen.OutputLevel = generator_output_level
gen.MomentumScale = 1.0
@@ -132,6 +133,11 @@ def run():
# Now build the physics list:
phys = simple.setupPhysics('QGSP_BERT')
+ ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics')
+ ph.addParticleConstructor('G4Geantino')
+ ph.addParticleConstructor('G4BosonConstructor')
+ ph.enableUI()
+ phys.adopt(ph)
phys.dump()
kernel.configure()
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 252cd5f8d..7511a3d8c 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -118,7 +118,7 @@ namespace DD4hep {
/// User data extension if required
std::auto_ptr<ParticleExtension> extension;
const G4VProcess *process; //! not persistent
- const G4ParticleDefinition *definition; //! not persistent
+ //const G4ParticleDefinition *definition; //! not persistent
/// Default constructor
Geant4Particle();
/// Constructor with ID initialization
@@ -173,7 +173,13 @@ namespace DD4hep {
/// Scalar particle energy
double energy() const;
/// Scalar particle momentum
- double momentum() const { return sqrt(momentum2()); }
+ double momentum() const { return sqrt(momentum2()); }
+ /// Geant4 charge of the particle
+ double charge() const { return double(particle->charge); }
+ /// Geant4 mass of the particle
+ double mass() const { return particle->mass; }
+ /// Geant4 time of the particle
+ double time() const { return particle->time; }
/// Access to the Geant4 particle name
std::string particleName() const;
/// Access to the Geant4 particle type
@@ -188,6 +194,8 @@ namespace DD4hep {
ThreeVector startVertex() const;
/// Access patricle momentum, energy as 4 vector
ThreeVector endVertex() const;
+ /// Access the Geant4 particle definition object (expensive!)
+ const G4ParticleDefinition *definition() const;
/// Various output formats:
diff --git a/DDG4/lcio/Geant4Output2LCIO.cpp b/DDG4/lcio/Geant4Output2LCIO.cpp
index 4f3544639..a2093c5fe 100644
--- a/DDG4/lcio/Geant4Output2LCIO.cpp
+++ b/DDG4/lcio/Geant4Output2LCIO.cpp
@@ -177,9 +177,9 @@ lcio::LCCollectionVec* Geant4Output2LCIO::saveParticles(Geant4ParticleMap* parti
// First create the particles
for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end();++i, ++cnt) {
int id = (*i).first;
- const Geant4Particle* p = (*i).second;
+ const Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->status);
- const G4ParticleDefinition* def = p->definition;
+ const G4ParticleDefinition* def = p.definition();
MCParticleImpl* q = new lcio::MCParticleImpl();
q->setPDG(p->pdgID);
diff --git a/DDG4/lcio/LCIOEventReader.cpp b/DDG4/lcio/LCIOEventReader.cpp
index 952b0593c..0272241b6 100644
--- a/DDG4/lcio/LCIOEventReader.cpp
+++ b/DDG4/lcio/LCIOEventReader.cpp
@@ -91,9 +91,11 @@ LCIOEventReader::readParticles(int event_number, vector<Particle*>& particles)
const double *vex = mcp->getEndpoint();
const float *spin = mcp->getSpin();
const int *color = mcp->getColorFlow();
- G4ParticleDefinition* def = tab->FindParticle(mcp->getPDG());
+ int pdg = mcp->getPDG();
+ //if ( 0 == pdg ) pdg = 211;
+ const G4ParticleDefinition* def = tab->FindParticle(pdg);
PropertyMask status(p->status);
- p->pdgID = mcp->getPDG();
+ p->pdgID = pdg;//mcp->getPDG();
p->charge = int(mcp->getCharge()*3.0);
p->psx = mom[0]*GeV;
p->psy = mom[1]*GeV;
@@ -106,7 +108,7 @@ LCIOEventReader::readParticles(int event_number, vector<Particle*>& particles)
p->vex = vex[0]*mm;
p->vey = vex[1]*mm;
p->vez = vex[2]*mm;
- p->definition = def;
+ //p->definition = def;
p->process = 0;
p->spin[0] = spin[0];
p->spin[1] = spin[1];
diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
index 53fa7e617..887609a7b 100644
--- a/DDG4/src/Geant4InputHandling.cpp
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -168,13 +168,13 @@ int DD4hep::Simulation::boostInteraction(const Geant4Action* /* caller */,
// Now move begin and end-vertex of all primary and generator particles accordingly
for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
- Geant4Particle* p = (*ip).second;
+ Geant4ParticleHandle p = (*ip).second;
double t = gamma * p->time + betagamma * p->vsx / c_light;
double x = gamma * p->vsx + betagamma * c_light * p->time;
double y = p->vsx;
double z = p->vsz;
- double m = p->definition->GetPDGMass();
+ double m = p.definition()->GetPDGMass();
double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
double px = betagamma * std::sqrt(e2) + gamma * p->psx;
double py = p->psy;
@@ -224,7 +224,15 @@ int DD4hep::Simulation::smearInteraction(const Geant4Action* /* caller */,
}
static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p) {
- G4PrimaryParticle* g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
+ G4PrimaryParticle* g4 = 0;
+ if ( 0 != p->pdgID ) {
+ g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
+ }
+ else {
+ const G4ParticleDefinition* def = p.definition();
+ g4 = new G4PrimaryParticle(def, p->psx, p->psy, p->psz);
+ g4->SetCharge(p.charge());
+ }
g4->SetMass(p->mass);
return g4;
}
diff --git a/DDG4/src/Geant4IsotropeGenerator.cpp b/DDG4/src/Geant4IsotropeGenerator.cpp
index dc4ef90f7..5f9c2e985 100644
--- a/DDG4/src/Geant4IsotropeGenerator.cpp
+++ b/DDG4/src/Geant4IsotropeGenerator.cpp
@@ -76,7 +76,7 @@ void Geant4IsotropeGenerator::operator()(G4Event*) {
p->status |= G4PARTICLE_GEN_STABLE;
p->mask = m_mask;
p->pdgID = m_particle->GetPDGEncoding();
- p->definition = m_particle;
+ //p->definition = m_particle;
p->psx = x1*momentum;
p->psy = x2*momentum;
p->psz = x3*momentum;
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index 996dedda3..4d9fe1abd 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -17,6 +17,9 @@
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4VProcess.hh"
+#include "G4ChargedGeantino.hh"
+#include "G4Geantino.hh"
+
#include <iostream>
using namespace DD4hep;
@@ -29,16 +32,16 @@ ParticleExtension::~ParticleExtension() {
/// Copy constructor
Geant4Particle::Geant4Particle(const Geant4Particle& c)
- : ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), mask(c.mask),
- steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
- status(c.status), charge(0),
- vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
- vex(c.vex), vey(c.vey), vez(c.vez),
- psx(c.psx), psy(c.psy), psz(c.psz),
- pex(c.pex), pey(c.pey), pez(c.pez),
- mass(c.mass), time(c.time), properTime(c.properTime),
- parents(c.parents), daughters(c.daughters), extension(),
- process(c.process), definition(c.definition)
+: ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), mask(c.mask),
+ steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
+ status(c.status), charge(0),
+ vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
+ vex(c.vex), vey(c.vey), vez(c.vez),
+ psx(c.psx), psy(c.psy), psz(c.psz),
+ pex(c.pex), pey(c.pey), pez(c.pez),
+ mass(c.mass), time(c.time), properTime(c.properTime),
+ parents(c.parents), daughters(c.daughters), extension(),
+ process(c.process)//, definition(c.definition)
{
InstanceCount::increment(this);
spin[0] = c.spin[0];
@@ -58,7 +61,7 @@ Geant4Particle::Geant4Particle()
psx(0.0), psy(0.0), psz(0.0),
pex(0.0), pey(0.0), pez(0.0),
mass(0.0), time(0.0), properTime(0.0),
- daughters(), extension(), process(0), definition(0)
+ daughters(), extension(), process(0)//, definition(0)
{
InstanceCount::increment(this);
spin[0] = spin[1] = spin[2] = 0;
@@ -75,7 +78,7 @@ Geant4Particle::Geant4Particle(int part_id)
psx(0.0), psy(0.0), psz(0.0),
pex(0.0), pey(0.0), pez(0.0),
mass(0.0), time(0.0), properTime(0.0),
- daughters(), extension(), process(0), definition(0)
+ daughters(), extension(), process(0)//, definition(0)
{
InstanceCount::increment(this);
spin[0] = spin[1] = spin[2] = 0;
@@ -123,7 +126,7 @@ Geant4Particle& Geant4Particle::get_data(Geant4Particle& c) {
time = c.time;
properTime = c.properTime;
process = c.process;
- definition = c.definition;
+ //definition = c.definition;
daughters = c.daughters;
parents = c.parents;
extension = c.extension;
@@ -138,13 +141,23 @@ void Geant4Particle::removeDaughter(int id_daughter) {
if ( j != daughters.end() ) daughters.erase(j);
}
-/// Access to the Geant4 particle name
-std::string Geant4ParticleHandle::particleName() const {
- if ( particle->definition ) return particle->definition->GetParticleName();
+/// Access the Geant4 particle definition object (expensive!)
+const G4ParticleDefinition* Geant4ParticleHandle::definition() const {
G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
+ if ( 0 == def && 0 == particle->pdgID ) {
+ if ( fabs(particle->charge) < 0.001 )
+ return G4Geantino::Definition();
+ return G4ChargedGeantino::Definition();
+ }
+ return def;
+}
+
+/// Access to the Geant4 particle name
+std::string Geant4ParticleHandle::particleName() const {
+ const G4ParticleDefinition* def = definition();
if ( def ) {
- particle->definition = def;
+ //particle->definition = def;
return def->GetParticleName();
}
#if 0
@@ -159,11 +172,9 @@ std::string Geant4ParticleHandle::particleName() const {
/// Access to the Geant4 particle type
std::string Geant4ParticleHandle::particleType() const {
- if ( particle->definition ) return particle->definition->GetParticleType();
- G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
- G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
+ const G4ParticleDefinition* def = definition();
if ( def ) {
- particle->definition = def;
+ //particle->definition = def;
return def->GetParticleType();
}
#if 0
diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp
index d2ec80adf..9a2f4893a 100644
--- a/DDG4/src/Geant4ParticleGun.cpp
+++ b/DDG4/src/Geant4ParticleGun.cpp
@@ -106,8 +106,8 @@ void Geant4ParticleGun::operator()(G4Event* event) {
p->vex = m_position.X();
p->vey = m_position.Y();
p->vez = m_position.Z();
- p->definition = m_particle;
- p->process = 0;
+ //p->definition = m_particle;
+ //p->process = 0;
p->spin[0] = 0;
p->spin[1] = 0;
p->spin[2] = 0;
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 09742b6c0..f0a1ce3aa 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -229,7 +229,6 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
m_currTrack.colorFlow[1] = prim_part->colorFlow[1];
m_currTrack.parents = prim_part->parents;
m_currTrack.daughters = prim_part->daughters;
- m_currTrack.definition = prim_part->definition;
m_currTrack.pdgID = prim_part->pdgID;
m_currTrack.mass = prim_part->mass;
}
@@ -245,7 +244,6 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
m_currTrack.colorFlow[1] = 0;
m_currTrack.parents.clear();
m_currTrack.daughters.clear();
- m_currTrack.definition = h.trackDef();
m_currTrack.pdgID = h.trackDef()->GetPDGEncoding();
m_currTrack.mass = h.trackDef()->GetPDGMass();
++m_globalParticleID;
@@ -440,14 +438,14 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
}
if ( ipar != m_particleMap.end() ) {
equivalents[(*i).first] = (*ipar).second->id; // requires (1) !
- Particle* p = (*ipar).second;
- const G4ParticleDefinition* def = p->definition;
+ Geant4ParticleHandle p = (*ipar).second;
+ const G4ParticleDefinition* def = p.definition();
int pdg = int(fabs(def->GetPDGEncoding())+0.1);
- if ( pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
+ if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
error("+++ ERROR: Geant4 particle for track:%d last known is:%d -- is gluon or quark!",(*i).second,g4_equiv);
}
pdg = int(fabs(p->pdgID)+0.1);
- if ( pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
+ if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
error("+++ ERROR(2): Geant4 particle for track:%d last known is:%d -- is gluon or quark!",(*i).second,g4_equiv);
}
}
--
GitLab