From 4514a018c22e1a4d75173bf6f519212d179d1b1a Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Tue, 16 Dec 2014 17:41:16 +0000
Subject: [PATCH] Change TrackerCombine SD to use energy weighted positions
---
DDDetectors/compact/SiD_Markus.xml | 8 +-
DDG4/examples/CLICSidSimuMarkus.py | 14 +-
DDG4/include/DDG4/DDG4Dict.h | 4 +-
DDG4/include/DDG4/Geant4Data.h | 4 +-
DDG4/include/DDG4/Geant4Kernel.h | 11 +-
DDG4/include/DDG4/Geant4SensDetAction.inl | 11 +-
DDG4/include/DDG4/Geant4StepHandler.h | 19 ++
DDG4/include/DDG4/Geant4TouchableHandler.h | 7 +-
DDG4/plugins/Geant4MaterialScanner.cpp | 216 ++++++++++++++
DDG4/plugins/Geant4Processes.cpp | 2 +-
DDG4/plugins/Geant4SDActions.cpp | 320 ++++++++++-----------
DDG4/python/DDG4.py | 12 +-
DDG4/src/Geant4Converter.cpp | 8 +-
DDG4/src/Geant4Data.cpp | 4 +-
DDG4/src/Geant4Kernel.cpp | 9 +-
DDG4/src/Geant4TouchableHandler.cpp | 11 +
16 files changed, 467 insertions(+), 193 deletions(-)
create mode 100644 DDG4/plugins/Geant4MaterialScanner.cpp
diff --git a/DDDetectors/compact/SiD_Markus.xml b/DDDetectors/compact/SiD_Markus.xml
index 291df9cd8..ec0f9b6e0 100644
--- a/DDDetectors/compact/SiD_Markus.xml
+++ b/DDDetectors/compact/SiD_Markus.xml
@@ -228,16 +228,16 @@
<!-- include ref="SiD/SiD_TrackerBarrel.xml" -->
<comment>Tracking detectors</comment>
- <include ref="${SiD_dir}/SiD_Tracker.xml"/>
+ <include ref="${SiD_dir}/SiD_Tracker.xml"/>
+
+ <include ref="SiD/SiD_Ecal.xml"/>
+<!--
<plugins>
<plugin name="DD4hep_SiTrackerBarrelSurfacePlugin">
<argument value="SiTrackerBarrel"/>
</plugin>
</plugins>
-<!--
-
- <include ref="SiD/SiD_Ecal.xml"/>
<plugins>
<plugin name="TestSurfaces"><argument value="EcalEndcap"/></plugin>
<plugin name="DD4hep_PolyhedraEndcapCalorimeterSurfacePlugin"><argument value="EcalEndcap"/></plugin>
diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py
index 9ca55c861..a929f349d 100644
--- a/DDG4/examples/CLICSidSimuMarkus.py
+++ b/DDG4/examples/CLICSidSimuMarkus.py
@@ -48,7 +48,7 @@ def run():
#VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
"""
Generation of primary particles from LCIO input files
- """
+
# First particle file reader
gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO1");
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep"
@@ -67,8 +67,11 @@ def run():
gen.Mask = 1
simple.buildInputStage([gen],output_level=generator_output_level)
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- #gen = simple.setupGun("Gun",particle='pi+',energy=20*GeV,position=(0.15*mm,0.12*mm,50*cm))
- #gen.OutputLevel = generator_output_level
+ """
+ gen = simple.setupGun("Gun",particle='geantino',energy=20*GeV,position=(0*mm,0*mm,0*cm),multiplicity=3)
+ gen.isotrop = False
+ gen.direction = (1,0,0)
+ gen.OutputLevel = generator_output_level
# And handle the simulation particles.
part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler")
@@ -112,8 +115,12 @@ def run():
seq,act = simple.setupCalorimeter('MuonEndcap')
seq,act = simple.setupCalorimeter('LumiCal')
seq,act = simple.setupCalorimeter('BeamCal')
+
"""
+ scan = DDG4.SteppingAction(kernel,'Geant4MaterialScanner/MaterialScan')
+ kernel.steppingAction().adopt(scan)
+
# Now build the physics list:
phys = simple.setupPhysics('QGSP_BERT')
ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics')
@@ -128,6 +135,7 @@ def run():
#DDG4.setPrintLevel(Output.DEBUG)
kernel.run()
+ print 'End of run. Terminating .......'
kernel.terminate()
if __name__ == "__main__":
diff --git a/DDG4/include/DDG4/DDG4Dict.h b/DDG4/include/DDG4/DDG4Dict.h
index 1ad0e8a53..67fed961a 100644
--- a/DDG4/include/DDG4/DDG4Dict.h
+++ b/DDG4/include/DDG4/DDG4Dict.h
@@ -101,7 +101,7 @@ namespace DD4hep {
/// Default destructor
inline Geant4HitData::~Geant4HitData() { }
/// Extract the MC contribution for a given hit from the step information
- inline Geant4HitData::Contribution Geant4HitData::extractContribution(G4Step*) { return Contribution(); }
+ inline Geant4HitData::Contribution Geant4HitData::extractContribution(const G4Step*) { return Contribution(); }
/// Default constructor
inline Geant4Tracker::Hit::Hit() { }
/// Initializing constructor
@@ -113,7 +113,7 @@ namespace DD4hep {
/// Clear hit content
inline Geant4Tracker::Hit& Geant4Tracker::Hit::clear() { return *this; }
/// Store Geant4 point and step information into tracker hit structure.
- inline Geant4Tracker::Hit& Geant4Tracker::Hit::storePoint(G4Step*, G4StepPoint*) { return *this;}
+ inline Geant4Tracker::Hit& Geant4Tracker::Hit::storePoint(const G4Step*, const G4StepPoint*) { return *this;}
/// Default constructor
inline Geant4Calorimeter::Hit::Hit() { }
/// Initializing constructor
diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h
index f4157072b..0abc1b2da 100644
--- a/DDG4/include/DDG4/Geant4Data.h
+++ b/DDG4/include/DDG4/Geant4Data.h
@@ -178,7 +178,7 @@ namespace DD4hep {
/// Default destructor
virtual ~Geant4HitData();
/// Extract the MC contribution for a given hit from the step information
- static Contribution extractContribution(G4Step* step);
+ static Contribution extractContribution(const G4Step* step);
};
/// Helper class to define structures used by the generic DDG4 tracker sensitive detector
@@ -221,7 +221,7 @@ namespace DD4hep {
/// Clear hit content
Hit& clear();
/// Store Geant4 point and step information into tracker hit structure.
- Hit& storePoint(G4Step* step, G4StepPoint* point);
+ Hit& storePoint(const G4Step* step, const G4StepPoint* point);
};
};
diff --git a/DDG4/include/DDG4/Geant4Kernel.h b/DDG4/include/DDG4/Geant4Kernel.h
index dd220a19d..4cc10552a 100644
--- a/DDG4/include/DDG4/Geant4Kernel.h
+++ b/DDG4/include/DDG4/Geant4Kernel.h
@@ -1,4 +1,4 @@
-// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+// $Id: $
//====================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------
@@ -289,11 +289,18 @@ namespace DD4hep {
}
/// Construct detector geometry using lcdd plugin
void loadGeometry(const std::string& compact_file);
- /// Run the simulation
+
+ /// Run the simulation: Configure Geant4
void configure();
+ /// Run the simulation: Initialize Geant4
void initialize();
+ /// Run the simulation: Simulate the number of events given by the property "NumEvents"
void run();
+ /// Run the simulation: Simulate the number of events "num_events" and modify the property "NumEvents"
+ void runEvents(int num_events);
+ /// Run the simulation: Terminate Geant4
void terminate();
+ /// Load XML file
void loadXML(const char* fname);
};
/// Declare property
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.inl b/DDG4/include/DDG4/Geant4SensDetAction.inl
index e56a8aeeb..16a8798c3 100644
--- a/DDG4/include/DDG4/Geant4SensDetAction.inl
+++ b/DDG4/include/DDG4/Geant4SensDetAction.inl
@@ -1,4 +1,4 @@
-// $Id: Geant4Field.cpp 513 2013-04-05 14:31:53Z gaede $
+// $Id: $
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
@@ -12,6 +12,7 @@
#include "DDG4/Geant4SensDetAction.h"
#endif
#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Data.h"
namespace DD4hep {
namespace Simulation {
@@ -63,5 +64,13 @@ namespace DD4hep {
Geant4Sensitive::clear(hce);
}
+ // Forward declarations
+ typedef Geant4HitData::Contribution HitContribution;
}
}
+
+
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4TouchableHandler.h"
+#include "DDG4/Geant4StepHandler.h"
+#include "DDG4/Geant4VolumeManager.h"
diff --git a/DDG4/include/DDG4/Geant4StepHandler.h b/DDG4/include/DDG4/Geant4StepHandler.h
index b17d69c79..85ad64586 100644
--- a/DDG4/include/DDG4/Geant4StepHandler.h
+++ b/DDG4/include/DDG4/Geant4StepHandler.h
@@ -84,6 +84,25 @@ namespace DD4hep {
const G4ThreeVector& postPosG4() const {
return post->GetPosition();
}
+ /// Average position vector of the step.
+ Position avgPosition() const {
+ const G4ThreeVector& p1 = pre->GetPosition();
+ const G4ThreeVector& p2 = post->GetPosition();
+ G4ThreeVector r = 0.5*(p2+p1);
+ return Position(r.x(),r.y(),r.z());
+ }
+ /// Direction calculated from the post- and pre-position ofthe step
+ Position direction() const {
+ const G4ThreeVector& p1 = pre->GetPosition();
+ const G4ThreeVector& p2 = post->GetPosition();
+ G4ThreeVector r = (p2-p1);
+ double len = r.r();
+ if ( len > 1e-15 ) {
+ r /= len;
+ return Position(r.x(),r.y(),r.z());
+ }
+ return Position();
+ }
Momentum preMom() const {
const G4ThreeVector& p = pre->GetMomentum();
return Momentum(p.x(), p.y(), p.z());
diff --git a/DDG4/include/DDG4/Geant4TouchableHandler.h b/DDG4/include/DDG4/Geant4TouchableHandler.h
index 9e1ade4d4..d67d12862 100644
--- a/DDG4/include/DDG4/Geant4TouchableHandler.h
+++ b/DDG4/include/DDG4/Geant4TouchableHandler.h
@@ -40,20 +40,23 @@ namespace DD4hep {
public:
typedef std::vector<const G4VPhysicalVolume*> Geant4PlacementPath;
-
+ /// Data member of the helper objects
const G4VTouchable* touchable;
/// Default constructor.
Geant4TouchableHandler(const G4VTouchable* t) : touchable(t) {}
/// Default constructor. Takes the step's pre-touchable
Geant4TouchableHandler(const G4Step* step);
+ /// Default constructor. Takes the step's pre-touchable
+ Geant4TouchableHandler(const G4Step* step, bool use_post_step_point);
/// Helper: Generate placement path from touchable object
Geant4PlacementPath placementPath(bool exception=false) const;
/// Helper: Access the placement path of a Geant4 touchable object as a string
std::string path() const;
-
+ /// Touchable history depth
+ int depth() const;
};
} // End namespace Simulation
diff --git a/DDG4/plugins/Geant4MaterialScanner.cpp b/DDG4/plugins/Geant4MaterialScanner.cpp
new file mode 100644
index 000000000..379bf106b
--- /dev/null
+++ b/DDG4/plugins/Geant4MaterialScanner.cpp
@@ -0,0 +1,216 @@
+// $Id:$
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+// Framework include files
+#include "DD4hep/Objects.h"
+#include "DDG4/Geant4SteppingAction.h"
+
+// Forward declarations
+class G4LogicalVolume;
+
+/// Namespace for the AIDA detector description toolkit
+namespace DD4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace Simulation {
+
+ /// Class to perform directional material scans using Geantinos.
+ /**
+ *
+ * \author M.Frank
+ * \version 1.0
+ * \ingroup DD4HEP_SIMULATION
+ */
+ class Geant4MaterialScanner : public Geant4SteppingAction {
+ protected:
+ typedef Geometry::Position Position;
+ /// Structure to hold the information of one simulation step.
+ class StepInfo {
+ public:
+ /// Pre-step and Post-step position
+ Position pre, post;
+ /// Reference to the logical volue
+ const G4LogicalVolume* volume;
+
+ /// Initializing constructor
+ StepInfo(const Position& pre, const Position& post, const G4LogicalVolume* volume);
+ /// Copy constructor
+ StepInfo(const StepInfo& c);
+ /// Default destructor
+ ~StepInfo() {}
+ /// Assignment operator
+ StepInfo& operator=(const StepInfo& c);
+ };
+ typedef std::vector<StepInfo*> Steps;
+
+ double m_sumX0;
+ double m_sumLambda;
+ double m_sumPath;
+ Steps m_steps;
+
+ public:
+ /// Standard constructor
+ Geant4MaterialScanner(Geant4Context* context, const std::string& name);
+ /// Default destructor
+ virtual ~Geant4MaterialScanner();
+ /// User stepping callback
+ virtual void operator()(const G4Step* step, G4SteppingManager* mgr);
+ /// Begin-of-tracking callback
+ virtual void begin(const G4Track* track);
+ /// End-of-tracking callback
+ virtual void end(const G4Track* track);
+ /// Registered callback on Begin-event
+ void beginEvent(const G4Event* event);
+ };
+ }
+}
+
+// $Id:$
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#include "DD4hep/InstanceCount.h"
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4TouchableHandler.h"
+#include "DDG4/Geant4StepHandler.h"
+#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4TrackingAction.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "G4LogicalVolume.hh"
+#include "G4Material.hh"
+
+using namespace std;
+using namespace CLHEP;
+
+using namespace DD4hep::Simulation;
+#include "DDG4/Factories.h"
+DECLARE_GEANT4ACTION(Geant4MaterialScanner)
+
+/// Initializing constructor
+Geant4MaterialScanner::StepInfo::StepInfo(const Position& prePos, const Position& postPos, const G4LogicalVolume* vol)
+: pre(prePos), post(postPos), volume(vol)
+{
+}
+
+/// Copy constructor
+Geant4MaterialScanner::StepInfo::StepInfo(const StepInfo& c)
+: pre(c.pre), post(c.post), volume(c.volume)
+{
+}
+
+/// Assignment operator
+Geant4MaterialScanner::StepInfo& Geant4MaterialScanner::StepInfo::operator=(const StepInfo& c) {
+ pre = c.pre;
+ post = c.post;
+ volume = c.volume;
+ return *this;
+}
+
+/// Standard constructor
+Geant4MaterialScanner::Geant4MaterialScanner(Geant4Context* context, const string& name)
+: Geant4SteppingAction(context,name)
+{
+ m_needsControl = true;
+ eventAction().callAtBegin(this,&Geant4MaterialScanner::beginEvent);
+ trackingAction().callAtEnd(this,&Geant4MaterialScanner::end);
+ trackingAction().callAtBegin(this,&Geant4MaterialScanner::begin);
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4MaterialScanner::~Geant4MaterialScanner() {
+ InstanceCount::decrement(this);
+}
+
+/// User stepping callback
+void Geant4MaterialScanner::operator()(const G4Step* step, G4SteppingManager*) {
+ Geant4StepHandler h(step);
+#if 0
+ Geant4TouchableHandler pre_handler(step);
+ string prePath = pre_handler.path();
+ Geant4TouchableHandler post_handler(step);
+ string postPath = post_handler.path();
+#endif
+ G4LogicalVolume* logVol = h.logvol(h.pre);
+ m_steps.push_back(new StepInfo(h.prePos(), h.postPos(), logVol));
+}
+
+/// Registered callback on Begin-event
+void Geant4MaterialScanner::beginEvent(const G4Event* /* event */) {
+ for_each(m_steps.begin(),m_steps.end(),DestroyObject<StepInfo*>());
+ m_steps.clear();
+ m_sumX0 = 0;
+ m_sumLambda = 0;
+ m_sumPath = 0;
+}
+
+/// Begin-of-tracking callback
+void Geant4MaterialScanner::begin(const G4Track* track) {
+ printP2("Starting tracking action for track ID=%d",track->GetTrackID());
+ for_each(m_steps.begin(),m_steps.end(),DestroyObject<StepInfo*>());
+ m_steps.clear();
+ m_sumX0 = 0;
+ m_sumLambda = 0;
+ m_sumPath = 0;
+}
+
+/// End-of-tracking callback
+void Geant4MaterialScanner::end(const G4Track* track) {
+ if ( !m_steps.empty() ) {
+ const char* line = " +--------------------------------------------------------------------------------------------------------------------------------------------------\n";
+ const char* fmt1 = " | %5d %-20s %3.0f %8.3f %8.4f %11.4f %11.4f %10.3f %8.2f %11.6f %11.6f (%7.2f,%7.2f,%7.2f)\n";
+ const char* fmt2 = " | %5d %-20s %3.0f %8.3f %8.4f %11.6g %11.6g %10.3f %8.2f %11.6f %11.6f (%7.2f,%7.2f,%7.2f)\n";
+ const Position& pre = m_steps[0]->pre;
+ const Position& post = m_steps[m_steps.size()-1]->post;
+
+ ::printf("%s + Material scan between: x_0 = (%7.2f,%7.2f,%7.2f) [cm] and x_1 = (%7.2f,%7.2f,%7.2f) [cm] TrackID:%d: \n%s",
+ line,pre.X()/cm,pre.Y()/cm,pre.Z()/cm,post.X()/cm,post.Y()/cm,post.Z()/cm,track->GetTrackID(),line);
+ ::printf(" | \\ %-11s Atomic Radiation Interaction Path Integrated Integrated Material\n","Material");
+ ::printf(" | Num. \\ %-11s Number/Z Mass/A Density Length Length Thickness Length X0 Lambda Endpoint \n","Name");
+ ::printf(" | Layer \\ %-11s [g/mole] [g/cm3] [cm] [cm] [cm] [cm] [cm] [cm] ( cm, cm, cm)\n","");
+ ::printf("%s",line);
+ int count = 1;
+ for(Steps::const_iterator i=m_steps.begin(); i!=m_steps.end(); ++i, ++count) {
+ const G4LogicalVolume* logVol = (*i)->volume;
+ G4Material* m = logVol->GetMaterial();
+ const Position& prePos = (*i)->pre;
+ const Position& postPos = (*i)->post;
+ Position direction = postPos - prePos;
+ double length = direction.R()/cm;
+ double intLen = m->GetNuclearInterLength()/cm;
+ double radLen = m->GetRadlen()/cm;
+ double density = m->GetDensity()/(gram/cm3);
+ double nLambda = length / intLen;
+ double nx0 = length / radLen;
+ double Aeff = 0.0;
+ double Zeff = 0.0;
+ const char* fmt = radLen >= 1e5 ? fmt2 : fmt1;
+ const double* fractions = m->GetFractionVector();
+ for(size_t i=0; i<m->GetNumberOfElements(); ++i) {
+ Zeff += fractions[i]*(m->GetElement(i)->GetZ());
+ Aeff += fractions[i]*(m->GetElement(i)->GetA())/gram;
+ }
+ m_sumX0 += nx0;
+ m_sumLambda += nLambda;
+ m_sumPath += length;
+ ::printf(fmt,count,m->GetName().c_str(),
+ Zeff, Aeff, density, radLen, intLen, length,
+ m_sumPath,m_sumX0,m_sumLambda,
+ postPos.X()/cm,postPos.Y()/cm,postPos.Z()/cm);
+ //cout << *m << endl;
+ }
+ for_each(m_steps.begin(),m_steps.end(),DestroyObject<StepInfo*>());
+ m_steps.clear();
+ }
+}
diff --git a/DDG4/plugins/Geant4Processes.cpp b/DDG4/plugins/Geant4Processes.cpp
index e09de97a7..c0b408e7a 100644
--- a/DDG4/plugins/Geant4Processes.cpp
+++ b/DDG4/plugins/Geant4Processes.cpp
@@ -1,4 +1,4 @@
-// $Id: Factories.h 797 2013-10-03 19:20:32Z markus.frank@cern.ch $
+// $Id:$
//====================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------
diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp
index 6fdc1b946..54a0704c5 100644
--- a/DDG4/plugins/Geant4SDActions.cpp
+++ b/DDG4/plugins/Geant4SDActions.cpp
@@ -1,4 +1,4 @@
-// $Id: Geant4Field.cpp 513 2013-04-05 14:31:53Z gaede $
+// $Id:$
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
@@ -8,28 +8,16 @@
//====================================================================
// Framework include files
#include "DDG4/Geant4SensDetAction.inl"
-#include "DDG4/Geant4Data.h"
-#include "DD4hep/Printout.h"
-
-using namespace std;
-
-#include "DDG4/Geant4TouchableHandler.h"
-#include "DDG4/Geant4StepHandler.h"
#include "DDG4/Geant4EventAction.h"
-#include "DDG4/Geant4VolumeManager.h"
-#include "DDG4/Geant4Mapping.h"
#include "G4OpticalPhoton.hh"
#include "G4VProcess.hh"
-
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
- typedef Geant4HitData::Contribution HitContribution;
-
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Geant4SensitiveAction<Geant4Tracker>
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -49,34 +37,34 @@ namespace DD4hep {
Position position = mean_direction(prePos,postPos);
double hit_len = direction.R();
- if ( step->GetTotalEnergyDeposit() < 1e-5 ) return true;
-
- if (hit_len > 0) {
- double new_len = mean_length(h.preMom(),h.postMom())/hit_len;
- direction *= new_len/hit_len;
+ if ( step->GetTotalEnergyDeposit() < std::numeric_limits<double>::epsilon() )
+ return true;
+ else if (hit_len > 0) {
+ double new_len = mean_length(h.preMom(),h.postMom())/hit_len;
+ direction *= new_len/hit_len;
}
print("Geant4Tracker","%s> Add hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f",
- c_name(),step->GetTotalEnergyDeposit(),position.X(),position.Y(),position.Z());
+ c_name(),step->GetTotalEnergyDeposit(),position.X(),position.Y(),position.Z());
Hit* hit = new Hit(h.trkID(), h.trkPdgID(), h.deposit(), h.track->GetGlobalTime());
if ( hit ) {
- HitContribution contrib = Hit::extractContribution(step);
- hit->cellID = cellID(step);
- hit->energyDeposit = contrib.deposit;
- hit->position = position;
- hit->momentum = direction;
- hit->length = hit_len;
- collection(m_collectionID)->add(hit);
- mark(h.track);
- if ( 0 == hit->cellID ) {
- hit->cellID = volumeID( step ) ;
- except("+++ Invalid CELL ID for hit!");
- }
- print("Geant4Tracker","%s> Hit with deposit:%f Pos:%f %f %f ID=%016X",
- c_name(),step->GetTotalEnergyDeposit(),position.X(),position.Y(),position.Z(),
- (void*)hit->cellID);
- Geant4TouchableHandler handler(step);
- print("Geant4Tracker","%s> Geant4 path:%s",c_name(),handler.path().c_str());
- return true;
+ HitContribution contrib = Hit::extractContribution(step);
+ hit->cellID = cellID(step);
+ hit->energyDeposit = contrib.deposit;
+ hit->position = position;
+ hit->momentum = direction;
+ hit->length = hit_len;
+ collection(m_collectionID)->add(hit);
+ mark(h.track);
+ if ( 0 == hit->cellID ) {
+ hit->cellID = volumeID( step ) ;
+ except("+++ Invalid CELL ID for hit!");
+ }
+ print("Geant4Tracker","%s> Hit with deposit:%f Pos:%f %f %f ID=%016X",
+ c_name(),step->GetTotalEnergyDeposit(),position.X(),position.Y(),position.Z(),
+ (void*)hit->cellID);
+ Geant4TouchableHandler handler(step);
+ print("Geant4Tracker","%s> Geant4 path:%s",c_name(),handler.path().c_str());
+ return true;
}
except("new() failed: Cannot allocate hit object");
return false;
@@ -100,22 +88,22 @@ namespace DD4hep {
long long int cell = cellID(step);
Hit* hit = coll->find<Hit>(CellIDCompare<Hit>(cell));
- if ( h.totalEnergy() < numeric_limits<double>::epsilon() ) {
- return true;
+ if ( h.totalEnergy() < std::numeric_limits<double>::epsilon() ) {
+ return true;
}
else if ( !hit ) {
- Geant4TouchableHandler handler(step);
- DDSegmentation::Vector3D pos = m_segmentation.position(cell);
- Position global = h.localToGlobal(pos);
- hit = new Hit(global);
- hit->cellID = cell;
- coll->add(hit);
- printM2("Geant4Calorimeter","%s> CREATE hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f %s",
- c_name(),contrib.deposit,pos.X,pos.Y,pos.Z,handler.path().c_str());
- if ( 0 == hit->cellID ) { // for debugging only!
- hit->cellID = cellID(step);
- except("+++ Invalid CELL ID for hit!");
- }
+ Geant4TouchableHandler handler(step);
+ DDSegmentation::Vector3D pos = m_segmentation.position(cell);
+ Position global = h.localToGlobal(pos);
+ hit = new Hit(global);
+ hit->cellID = cell;
+ coll->add(hit);
+ printM2("Geant4Calorimeter","%s> CREATE hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f %s",
+ c_name(),contrib.deposit,pos.X,pos.Y,pos.Z,handler.path().c_str());
+ if ( 0 == hit->cellID ) { // for debugging only!
+ hit->cellID = cellID(step);
+ except("+++ Invalid CELL ID for hit!");
+ }
}
hit->truth.push_back(contrib);
hit->energyDeposit += contrib.deposit;
@@ -145,33 +133,33 @@ namespace DD4hep {
G4Track * track = step->GetTrack();
// check that particle is optical photon:
if( track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition() ) {
- return false;
+ return false;
}
else if ( track->GetCreatorProcess()->G4VProcess::GetProcessName() != "Cerenkov") {
- track->SetTrackStatus(fStopAndKill);
- return false;
+ track->SetTrackStatus(fStopAndKill);
+ return false;
}
else {
- typedef Geant4Calorimeter::Hit Hit;
- StepHandler h(step);
- HitCollection* coll = collection(m_collectionID);
- HitContribution contrib = Hit::extractContribution(step);
- Position pos = h.prePos();
- Hit* hit = coll->find<Hit>(PositionCompare<Hit,Position>(pos));
- if ( !hit ) {
- hit = new Hit(pos);
- hit->cellID = volumeID(step);
- coll->add(hit);
- if ( 0 == hit->cellID ) {
- hit->cellID = volumeID(step);
- except("+++ Invalid CELL ID for hit!");
- }
- }
- hit->energyDeposit += contrib.deposit;
- hit->truth.push_back(contrib);
- track->SetTrackStatus(fStopAndKill); // don't step photon any further
- mark(h.track);
- return true;
+ typedef Geant4Calorimeter::Hit Hit;
+ StepHandler h(step);
+ HitCollection* coll = collection(m_collectionID);
+ HitContribution contrib = Hit::extractContribution(step);
+ Position pos = h.prePos();
+ Hit* hit = coll->find<Hit>(PositionCompare<Hit,Position>(pos));
+ if ( !hit ) {
+ hit = new Hit(pos);
+ hit->cellID = volumeID(step);
+ coll->add(hit);
+ if ( 0 == hit->cellID ) {
+ hit->cellID = volumeID(step);
+ except("+++ Invalid CELL ID for hit!");
+ }
+ }
+ hit->energyDeposit += contrib.deposit;
+ hit->truth.push_back(contrib);
+ track->SetTrackStatus(fStopAndKill); // don't step photon any further
+ mark(h.track);
+ return true;
}
}
typedef Geant4SensitiveAction<Geant4OpticalCalorimeter> Geant4OpticalCalorimeterAction;
@@ -189,128 +177,138 @@ namespace DD4hep {
struct TrackerCombine {
typedef Geant4HitCollection HitCollection;
Geant4Tracker::Hit pre, post;
+ Position mean_pos;
Geant4Sensitive* sensitive;
+ double mean_time;
double e_cut;
int current;
int combined;
long long int cell;
-
+
TrackerCombine() : pre(), post(), sensitive(0), e_cut(0.0), current(-1), combined(0), cell(0) {
}
/// Start a new hit
void start(G4Step* step, G4StepPoint* point) {
- pre.storePoint(step,point);
- pre.truth.deposit = 0.0;
- current = pre.truth.trackID;
- sensitive->mark(step->GetTrack());
- post = pre;
- combined = 0;
- cell = 0;
+ pre.storePoint(step,point);
+ pre.truth.deposit = 0.0;
+ current = pre.truth.trackID;
+ sensitive->mark(step->GetTrack());
+ mean_pos.SetXYZ(0,0,0);
+ mean_time = 0;
+ post = pre;
+ combined = 0;
+ cell = 0;
}
/// Update energy and track information during hit info accumulation
void update(G4Step* step) {
- post.storePoint(step,step->GetPostStepPoint());
- pre.truth.deposit += post.truth.deposit;
- if ( 0 == cell ) {
- cell = sensitive->cellID(step);
- if ( 0 == cell ) {
- cell = sensitive->volumeID(step) ;
- sensitive->except("+++ Invalid CELL ID for hit!");
- }
- }
- ++combined;
+ post.storePoint(step,step->GetPostStepPoint());
+ pre.truth.deposit += post.truth.deposit;
+ mean_pos.SetX(mean_pos.x()+post.position.x()*post.truth.deposit);
+ mean_pos.SetY(mean_pos.y()+post.position.y()*post.truth.deposit);
+ mean_pos.SetZ(mean_pos.z()+post.position.z()*post.truth.deposit);
+ mean_time += post.truth.time*post.truth.deposit;
+ if ( 0 == cell ) {
+ cell = sensitive->cellID(step);
+ if ( 0 == cell ) {
+ cell = sensitive->volumeID(step) ;
+ sensitive->except("+++ Invalid CELL ID for hit!");
+ }
+ }
+ ++combined;
}
/// Clear collected information and restart for new hit
void clear() {
- post.clear();
- pre.clear();
- current = -1;
- combined = 0;
- cell = 0;
+ mean_pos.SetXYZ(0,0,0);
+ mean_time = 0;
+ post.clear();
+ pre.clear();
+ current = -1;
+ combined = 0;
+ cell = 0;
}
/// Helper function to decide if the hit has to be extracted and saved in the collection
bool mustSaveTrack(const G4Track* tr) const {
- return current > 0 && current != tr->GetTrackID();
+ return current > 0 && current != tr->GetTrackID();
}
/// Extract hit information and add the created hit to the collection
void extractHit(Geant4HitCollection* collection) {
- if ( current == -1 ) {
- return;
- }
- double deposit = pre.truth.deposit, time = pre.truth.time;
- Position pos = 0.5 * (pre.position + post.position);
- Momentum mom = 0.5 * (pre.momentum + post.momentum);
- double path_len = (post.position - pre.position).R();
- Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID,
- pre.truth.pdgID,
- deposit,time);
- hit->position = pos;
- hit->momentum = mom;
- hit->length = path_len;
- hit->cellID = cell;
- collection->add(hit);
- sensitive->printM2("+++ TrackID:%6d [%s] CREATE hit combination with %2d deposit(s):"
- " %e MeV Pos:%8.2f %8.2f %8.2f",
- pre.truth.trackID,sensitive->c_name(),combined,pre.truth.deposit/MeV,
- pos.X()/mm,pos.Y()/mm,pos.Z()/mm);
- clear();
+ if ( current == -1 ) {
+ return;
+ }
+ double deposit = pre.truth.deposit, time = mean_time / deposit;
+ Position pos = mean_pos / deposit;
+ Momentum mom = 0.5 * (pre.momentum + post.momentum);
+ double path_len = (post.position - pre.position).R();
+ Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID,
+ pre.truth.pdgID,
+ deposit,time);
+ hit->position = pos;
+ hit->momentum = mom;
+ hit->length = path_len;
+ hit->cellID = cell;
+ collection->add(hit);
+ sensitive->printM2("+++ TrackID:%6d [%s] CREATE hit combination with %2d deposit(s):"
+ " %e MeV Pos:%8.2f %8.2f %8.2f",
+ pre.truth.trackID,sensitive->c_name(),combined,pre.truth.deposit/MeV,
+ pos.X()/mm,pos.Y()/mm,pos.Z()/mm);
+ clear();
}
/// Method for generating hit(s) using the information of G4Step object.
G4bool process(G4Step* step, G4TouchableHistory* ) {
- Geant4StepHandler h(step);
- void *prePV = h.volume(h.pre), *postPV = h.volume(h.post);
-
- Geant4HitCollection* coll = sensitive->collection(0);
- /// If we are handling a new track, then store the content of the previous one.
- if ( mustSaveTrack(h.track) ) {
- extractHit(coll);
- }
- /// There must be something in.
- if ( h.deposit()/keV <= 0 ) {
- return false;
- }
- /// Initialize the deposits of the next hit.
- if ( current < 0 ) {
- start(step, h.pre);
- }
- /// ....update .....
- update(step);
-
- if ( prePV != postPV ) {
- void* postSD = h.sd(h.post);
- extractHit(coll);
- if ( 0 != postSD ) {
- void* preSD = h.sd(h.pre);
- if ( preSD == postSD ) {
- start(step,h.post);
- }
- }
- }
- else if ( h.track->GetTrackStatus() == fStopAndKill ) {
- extractHit(coll);
- }
- return true;
+ Geant4StepHandler h(step);
+ void *prePV = h.volume(h.pre), *postPV = h.volume(h.post);
+
+ Geant4HitCollection* coll = sensitive->collection(0);
+ /// If we are handling a new track, then store the content of the previous one.
+ if ( mustSaveTrack(h.track) ) {
+ extractHit(coll);
+ }
+ /// There must be something in.
+ if ( h.deposit()/keV <= 0 ) {
+ return false;
+ }
+ /// Initialize the deposits of the next hit.
+ if ( current < 0 ) {
+ start(step, h.pre);
+ }
+ /// ....update .....
+ update(step);
+
+ if ( prePV != postPV ) {
+ void* postSD = h.sd(h.post);
+ extractHit(coll);
+ if ( 0 != postSD ) {
+ void* preSD = h.sd(h.pre);
+ if ( preSD == postSD ) {
+ start(step,h.post);
+ }
+ }
+ }
+ else if ( h.track->GetTrackStatus() == fStopAndKill ) {
+ extractHit(coll);
+ }
+ return true;
}
/// Post-event action callback
void endEvent(const G4Event* /* event */) {
- // We need to add the possibly last added hit to the collection here.
- // otherwise the last hit would be assigned to the next event and the
- // MC truth would be screwed.
- //
- // Alternatively the 'update' method would become rather CPU consuming,
- // beacuse the extract action would have to be recalculated over and over.
- if ( current > 0 ) {
- Geant4HitCollection* coll = sensitive->collection(0);
- extractHit(coll);
- }
+ // We need to add the possibly last added hit to the collection here.
+ // otherwise the last hit would be assigned to the next event and the
+ // MC truth would be screwed.
+ //
+ // Alternatively the 'update' method would become rather CPU consuming,
+ // beacuse the extract action would have to be recalculated over and over.
+ if ( current > 0 ) {
+ Geant4HitCollection* coll = sensitive->collection(0);
+ extractHit(coll);
+ }
}
};
diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py
index 99d178d94..79c8821f5 100644
--- a/DDG4/python/DDG4.py
+++ b/DDG4/python/DDG4.py
@@ -57,7 +57,7 @@ def importConstants(lcdd,namespace=None,debug=False):
todo = {}
strings = {}
for c in lcdd.constants():
- if c.second.dataType != 'string':
+ if c.second.dataType == 'string':
strings[c.first] = c.second.GetTitle()
else:
todo[c.first] = c.second.GetTitle().replace('(int)','')
@@ -240,11 +240,9 @@ class Simple:
if kernel is None:
self.kernel = Kernel()
self.lcdd = self.kernel.lcdd()
- self.calo = calo
- self.tracker = tracker
self.sensitive_types = {}
- self.sensitive_types['tracker'] = self.tracker
- self.sensitive_types['calorimeter'] = self.calo
+ self.sensitive_types['tracker'] = tracker
+ self.sensitive_types['calorimeter'] = calo
def execute(self):
self.kernel.configure()
@@ -276,13 +274,13 @@ class Simple:
def setupCalorimeter(self,name,type=None):
sd = self.lcdd.sensitiveDetector(name)
sd.setType('calorimeter')
- if type is None: type = self.calo
+ if type is None: type = self.sensitive_types['calorimeter']
return self.setupDetector(name,type)
def setupTracker(self,name,type=None):
sd = self.lcdd.sensitiveDetector(name)
sd.setType('tracker')
- if type is None: type = self.tracker
+ if type is None: type = self.sensitive_types['tracker']
return self.setupDetector(name,type)
def setupPhysics(self,name):
diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp
index 1e9f1f72f..2b81f0ffe 100644
--- a/DDG4/src/Geant4Converter.cpp
+++ b/DDG4/src/Geant4Converter.cpp
@@ -330,18 +330,22 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
}
if (m->IsMixture()) {
double A_total = 0.0;
+ double W_total = 0.0;
TGeoMixture* mix = (TGeoMixture*) m;
int nElements = mix->GetNelements();
mat = new G4Material(name, density, nElements, state, m->GetTemperature(), m->GetPressure());
- for (int i = 0; i < nElements; ++i)
+ for (int i = 0; i < nElements; ++i) {
A_total += (mix->GetAmixt())[i];
+ W_total += (mix->GetWmixt())[i];
+ }
for (int i = 0; i < nElements; ++i) {
TGeoElement* e = mix->GetElement(i);
G4Element* g4e = (G4Element*) handleElement(e->GetName(), Atom(e));
if (!g4e) {
printout(ERROR, "Material", "Missing component %s for material %s.", e->GetName(), mix->GetName());
}
- mat->AddElement(g4e, (mix->GetAmixt())[i] / A_total);
+ //mat->AddElement(g4e, (mix->GetAmixt())[i] / A_total);
+ mat->AddElement(g4e, (mix->GetWmixt())[i] / W_total);
}
}
else {
diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp
index 240f39951..3c81189d0 100644
--- a/DDG4/src/Geant4Data.cpp
+++ b/DDG4/src/Geant4Data.cpp
@@ -60,7 +60,7 @@ Geant4HitData::~Geant4HitData() {
}
/// Extract the MC contribution for a given hit from the step information
-Geant4HitData::Contribution Geant4HitData::extractContribution(G4Step* step) {
+Geant4HitData::Contribution Geant4HitData::extractContribution(const G4Step* step) {
Geant4StepHandler h(step);
double deposit =
(h.trackDef() == G4OpticalPhoton::OpticalPhotonDefinition()) ? h.trkEnergy() : h.totalEnergy();
@@ -109,7 +109,7 @@ Geant4Tracker::Hit& Geant4Tracker::Hit::clear() {
}
/// Store Geant4 point and step information into tracker hit structure.
-Geant4Tracker::Hit& Geant4Tracker::Hit::storePoint(G4Step* step, G4StepPoint* pnt) {
+Geant4Tracker::Hit& Geant4Tracker::Hit::storePoint(const G4Step* step, const G4StepPoint* pnt) {
G4Track* trk = step->GetTrack();
G4ThreeVector pos = pnt->GetPosition();
G4ThreeVector mom = pnt->GetMomentum();
diff --git a/DDG4/src/Geant4Kernel.cpp b/DDG4/src/Geant4Kernel.cpp
index fe4b9470e..df65d48c2 100644
--- a/DDG4/src/Geant4Kernel.cpp
+++ b/DDG4/src/Geant4Kernel.cpp
@@ -179,22 +179,23 @@ void Geant4Kernel::loadGeometry(const std::string& compact_file) {
void Geant4Kernel::loadXML(const char* fname) {
const char* args[] = { fname, 0 };
m_lcdd.apply("DD4hepXMLLoader", 1, (char**) args);
- //return *this;
}
void Geant4Kernel::configure() {
Geant4Exec::configure(*this);
- //return *this;
}
void Geant4Kernel::initialize() {
Geant4Exec::initialize(*this);
- //return *this;
}
void Geant4Kernel::run() {
Geant4Exec::run(*this);
- //return *this;
+}
+
+void Geant4Kernel::runEvents(int num_events) {
+ m_numEvent = num_events;
+ Geant4Exec::run(*this);
}
void Geant4Kernel::terminate() {
diff --git a/DDG4/src/Geant4TouchableHandler.cpp b/DDG4/src/Geant4TouchableHandler.cpp
index 7a4f423a7..4793a9264 100644
--- a/DDG4/src/Geant4TouchableHandler.cpp
+++ b/DDG4/src/Geant4TouchableHandler.cpp
@@ -19,11 +19,22 @@
using namespace DD4hep::Simulation;
+/// Default constructor. Takes the step's pre-touchable
+Geant4TouchableHandler::Geant4TouchableHandler(const G4Step* step, bool use_post_step_point) {
+ const G4StepPoint* p = use_post_step_point ? step->GetPostStepPoint() : step->GetPreStepPoint();
+ touchable = p->GetTouchable();
+}
+
/// Default constructor. Takes the step's pre-touchable
Geant4TouchableHandler::Geant4TouchableHandler(const G4Step* step) {
touchable = step->GetPreStepPoint()->GetTouchable();
}
+/// Touchable history depth
+int Geant4TouchableHandler::depth() const {
+ return touchable->GetHistoryDepth();
+}
+
/// Helper: Generate placement path from touchable object
Geant4TouchableHandler::Geant4PlacementPath Geant4TouchableHandler::placementPath(bool exception) const {
Geant4PlacementPath path;
--
GitLab