diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml b/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml
index 7a72d28eb..386c212c4 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml
@@ -16,15 +16,15 @@
-
+
-
+
-
+
@@ -40,10 +40,10 @@
-
+
-
+
@@ -68,7 +68,7 @@
-
+
@@ -82,7 +82,7 @@
-
+
@@ -103,7 +103,7 @@
-
+
@@ -113,7 +113,7 @@
-
+
@@ -122,7 +122,7 @@
-
+
@@ -135,71 +135,71 @@
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
@@ -207,10 +207,10 @@
-
+
-
+
@@ -232,18 +232,18 @@
-
+
-
+
-
-
+
+
diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml b/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml
index 8961ffac4..b97641e84 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml
@@ -13,7 +13,7 @@
The compact format for the dual-readout calorimeter (for FCCee IDEA)
-
+
-
+
@@ -691,7 +691,7 @@
- system:5,assembly:1,eta:-8,phi:9,x:32:-11,y:-9,c:1,module:2
+ system:5,assembly:1,eta:-8,phi:9,x:32:-12,y:-12,c:1,module:2
diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml b/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
index 2cdcd42f1..5fe33a2c8 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
@@ -11,7 +11,7 @@
url="no"
status="development"
version="o1_v03">
-
+
Version o1_v03 of the IDEA detector
@@ -37,9 +37,9 @@
-
+
-
@@ -65,7 +65,7 @@
-
+
diff --git a/detector/calorimeter/dual-readout/src/DRconstructor.cpp b/detector/calorimeter/dual-readout/src/DRconstructor.cpp
index c01afa77f..094a9e2d4 100644
--- a/detector/calorimeter/dual-readout/src/DRconstructor.cpp
+++ b/detector/calorimeter/dual-readout/src/DRconstructor.cpp
@@ -87,7 +87,7 @@ void ddDRcalo::DRconstructor::implementTowers(xml_comp_t& x_theta, dd4hep::DDSeg
dd4hep::Volume towerVol( "tower", tower, fDescription->material(x_theta.materialStr()) );
towerVol.setVisAttributes(*fDescription, x_theta.visStr());
-
+
implementFibers(x_theta, towerVol, tower, param);
xml_comp_t x_wafer ( fX_sipmDim.child( _Unicode(sipmWafer) ) );
@@ -123,7 +123,10 @@ void ddDRcalo::DRconstructor::placeAssembly(dd4hep::DDSegmentation::DRparamBase_
int towerId32 = fSegmentation->getFirst32bits(towerId64);
dd4hep::Position towerPos = param->GetTowerPos(nPhi) + dd4hep::Position(0, 0, -(fX_worldTube.height()/2.));
- AssemblyBoxVol.placeVolume( towerVol, towerId32, dd4hep::Transform3D( param->GetRotationZYX(nPhi), towerPos ) );
+ dd4hep::PlacedVolume towerPhys = AssemblyBoxVol.placeVolume( towerVol, towerId32, dd4hep::Transform3D( param->GetRotationZYX(nPhi), towerPos ) );
+ towerPhys.addPhysVolID("eta", towerNoLR);
+ towerPhys.addPhysVolID("phi", nPhi);
+ towerPhys.addPhysVolID("module", 2);
// Remove sipmLayer
dd4hep::Position sipmPos = param->GetSipmLayerPos(nPhi) + dd4hep::Position(0, 0, -(fX_worldTube.height()/2.));
@@ -250,6 +253,8 @@ void ddDRcalo::DRconstructor::implementFiber(dd4hep::Volume& towerVol, dd4hep::P
if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreC.visStr());
cladVol.placeVolume( coreVol );
+ // we use the region for the sensitive elements for
+ // manipulating optical photons (DRCaloFastSimModel)
coreVol.setRegion(*fDescription, fX_det.regionStr());
cladVol.setRegion(*fDescription, fX_det.regionStr());
} else { // s fiber
@@ -261,6 +266,8 @@ void ddDRcalo::DRconstructor::implementFiber(dd4hep::Volume& towerVol, dd4hep::P
if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreS.visStr());
cladVol.placeVolume( coreVol );
+ // we use the region for the sensitive elements for
+ // manipulating optical photons (DRCaloFastSimModel)
coreVol.setRegion(*fDescription, fX_det.regionStr());
cladVol.setRegion(*fDescription, fX_det.regionStr());
}
diff --git a/example/SteeringFile_IDEA_o1_v03.py b/example/SteeringFile_IDEA_o1_v03.py
index 102085aa5..0796bb359 100644
--- a/example/SteeringFile_IDEA_o1_v03.py
+++ b/example/SteeringFile_IDEA_o1_v03.py
@@ -16,11 +16,11 @@
## Macro file to execute for runType 'run' or 'vis'
SIM.macroFile = ""
## number of events to simulate, used in batch mode
-SIM.numberOfEvents = 1
+SIM.numberOfEvents = 10
## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported
SIM.outputFile = "testIDEA_o1_v03.root"
## Physics list to use in simulation
-SIM.physicsList = None
+SIM.physicsList = "FTFP_BERT"
## Verbosity use integers from 1(most) to 7(least) verbose
## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS
SIM.printLevel = 3
@@ -112,7 +112,7 @@
##
## SIM.action.mapActions['tpc'] = "TPCSDAction"
##
-SIM.action.mapActions = {"DRcalo": "DRCaloSDAction"}
+SIM.action.mapActions["DRcalo"] = "DRCaloSDAction"
## set the default run action
SIM.action.run = []
@@ -177,7 +177,8 @@
## default filter for calorimeter sensitive detectors;
## this is applied if no other filter is used for a calorimeter
##
-SIM.filter.calo = "edep0"
+# note: do not turn on the calo filter, otherwise all optical photons will be killed!
+SIM.filter.calo = ""
## list of filter objects: map between name and parameter dictionary
SIM.filter.filters = {
@@ -229,6 +230,8 @@
## Print information about Sensitives
SIM.geometry.enablePrintSensitives = False
+## configure regex SD
+SIM.geometry.regexSensitiveDetector["DRcalo"] = {"Match": ["(core|clad)"], "OutputLevel": 3}
################################################################################
## Configuration for the GuineaPig InputFiles
@@ -246,7 +249,7 @@
################################################################################
## direction of the particle gun, 3 vector
-SIM.gun.direction = (1.0, 1.0, 1.0)
+SIM.gun.direction = (1.0, 0.1, 0.1)
## choose the distribution of the random direction for theta
##
@@ -264,7 +267,7 @@
## Total energy (including mass) for the particle gun.
##
## If not None, it will overwrite the setting of momentumMin and momentumMax
-SIM.gun.energy = None
+SIM.gun.energy = 10.0 * GeV
## Maximal pseudorapidity for random distibution (overrides thetaMin)
SIM.gun.etaMax = None
@@ -280,12 +283,12 @@
SIM.gun.isotrop = False
## Maximal momentum when using distribution (default = 0.0)
-SIM.gun.momentumMax = 10000.0
+# SIM.gun.momentumMax = 10000.0
## Minimal momentum when using distribution (default = 0.0)
-SIM.gun.momentumMin = 0.0
+# SIM.gun.momentumMin = 0.0
SIM.gun.multiplicity = 1
-SIM.gun.particle = "mu-"
+SIM.gun.particle = "e-"
## Maximal azimuthal angle for random distribution
SIM.gun.phiMax = None
@@ -552,7 +555,7 @@ def Geant4Output2EDM4hep_DRC_plugin(dd4hepSimulation):
## Set printlevel to DEBUG to see a printout of all range cuts,
## but this only works if range cut is not "None"
##
-SIM.physics.rangecut = 0.7
+SIM.physics.rangecut = None
## Set of PDG IDs that will not be passed from the input record to Geant4.
##
@@ -613,17 +616,11 @@ def setupOpticalPhysics(kernel):
cerenkov.enableUI()
seq.adopt(cerenkov)
- scint = PhysicsList(kernel, "Geant4ScintillationPhysics/ScintillationPhys")
- scint.VerboseLevel = 1
- scint.TrackSecondariesFirst = True
- scint.BoundaryInvokeSD = True
- scint.enableUI()
- seq.adopt(scint)
-
opt = PhysicsList(kernel, "Geant4OpticalPhotonPhysics/OpticalGammaPhys")
opt.addParticleConstructor("G4OpticalPhoton")
opt.VerboseLevel = 1
- opt.BoundaryInvokeSD = True
+ # set BoundaryInvokeSD to true when using DRC wafer as the SD
+ # opt.BoundaryInvokeSD = True
opt.enableUI()
seq.adopt(opt)
@@ -656,7 +653,8 @@ def setupDRCFastSim(kernel):
phys.dump()
-SIM.physics.setupUserPhysics(setupDRCFastSim)
+# turn-on fastsim if the skipScint option of the SD is set to false
+# SIM.physics.setupUserPhysics(setupDRCFastSim)
################################################################################
## Properties for the random number generator
diff --git a/plugins/DRCaloFastSimModel.cpp b/plugins/DRCaloFastSimModel.cpp
index 5caa21456..c29789386 100644
--- a/plugins/DRCaloFastSimModel.cpp
+++ b/plugins/DRCaloFastSimModel.cpp
@@ -1,15 +1,7 @@
// Framework include files
#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
+#include "G4FastStep.hh"
// C/C++ include files
#include "DRCaloFastSimModel.h"
@@ -17,246 +9,9 @@
/// Namespace for the AIDA detector description toolkit
namespace dd4hep
{
-
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace sim
{
-
- class DRCFiberModel
- {
- public:
- // G4FastSimHitMaker hitMaker { };
- G4OpBoundaryProcess *pOpBoundaryProc{nullptr};
- G4OpAbsorption *pOpAbsorption{nullptr};
- G4OpWLS *pOpWLS{nullptr};
- G4bool fProcAssigned{false};
-
- FastFiberData mDataPrevious{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))};
- FastFiberData mDataCurrent{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))};
-
- G4int fSafety{2};
- G4double mNtransport{0.};
- G4double mTransportUnit{0.};
- G4ThreeVector mFiberPos{G4ThreeVector(0)};
- G4ThreeVector mFiberAxis{G4ThreeVector(0)};
- G4bool fKill{false};
- G4bool fTransported{false};
- G4bool fSwitch{true};
- G4int fVerbose{0};
-
- G4bool checkTotalInternalReflection(const G4Track *track)
- {
- if (!fProcAssigned)
- setPostStepProc(track); // locate OpBoundaryProcess only once
-
- if (track->GetTrackStatus() == fStopButAlive || track->GetTrackStatus() == fStopAndKill)
- return false;
-
- // accumulate step length
- mDataCurrent.AddStepLengthInterval(track->GetStepLength());
-
- G4int theStatus = pOpBoundaryProc->GetStatus();
-
- if (fVerbose > 1)
- {
- G4cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID();
- G4cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus;
- G4cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << G4endl;
- }
-
- // skip exceptional iteration with FresnelReflection
- if (theStatus == G4OpBoundaryProcessStatus::FresnelReflection)
- mDataCurrent.SetOpBoundaryStatus(theStatus);
-
- // some cases have a status StepTooSmall when the reflection happens between the boundary of cladding & outer volume (outside->cladding) since the outer volume is not a G4Region
- if (theStatus == G4OpBoundaryProcessStatus::TotalInternalReflection || theStatus == G4OpBoundaryProcessStatus::StepTooSmall)
- {
- if (theStatus != G4OpBoundaryProcessStatus::TotalInternalReflection)
- { // skip StepTooSmall if the track already has TotalInternalReflection
- if (mDataCurrent.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection)
- return false;
- if (mDataPrevious.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection)
- return false;
- }
-
- G4int trackID = track->GetTrackID();
- G4double kineticEnergy = track->GetKineticEnergy();
- G4double globalTime = track->GetGlobalTime();
- G4double pathLength = track->GetStepLength();
- G4ThreeVector globalPosition = track->GetPosition();
- G4ThreeVector momentumDirection = track->GetMomentumDirection();
- G4ThreeVector polarization = track->GetPolarization();
-
- auto candidate = FastFiberData(trackID, kineticEnergy, globalTime, pathLength, globalPosition, momentumDirection, polarization, theStatus);
- if (pOpAbsorption != nullptr)
- candidate.SetAbsorptionNILL(pOpAbsorption->GetNumberOfInteractionLengthLeft());
- if (pOpWLS != nullptr)
- candidate.SetWLSNILL(pOpWLS->GetNumberOfInteractionLengthLeft());
-
- G4bool repetitive = false;
- if (candidate.checkRepetitive(mDataCurrent, false) && mDataCurrent.checkRepetitive(mDataPrevious))
- repetitive = true;
-
- mDataPrevious = mDataCurrent;
- mDataCurrent = candidate;
-
- return repetitive;
- }
-
- return false;
- }
-
- G4bool checkAbsorption(const G4double prevNILL, const G4double currentNILL)
- {
- if (prevNILL < 0. || currentNILL < 0.)
- return false; // the number of interaction length left has to be reset
- if (prevNILL == currentNILL)
- return false; // no absorption
- if (prevNILL == DBL_MAX || currentNILL == DBL_MAX)
- return false; // NILL is re-initialized
-
- G4double deltaNILL = prevNILL - currentNILL;
-
- if (currentNILL - deltaNILL * (mNtransport + fSafety) < 0.)
- return true; // absorbed before reaching fiber end
-
- return false;
- }
-
- G4bool checkNILL()
- {
- if (!fTransported)
- return true; // do nothing if the track is not already transported
-
- G4double wlsNILL = DBL_MAX;
- G4double absorptionNILL = DBL_MAX;
-
- if (pOpWLS != nullptr)
- {
- wlsNILL = pOpWLS->GetNumberOfInteractionLengthLeft();
- if (mDataPrevious.GetWLSNILL() == DBL_MAX || mDataCurrent.GetWLSNILL() == DBL_MAX)
- return true; // NILL is re-initialized
- }
-
- if (pOpAbsorption != nullptr)
- {
- absorptionNILL = pOpAbsorption->GetNumberOfInteractionLengthLeft();
- if (mDataPrevious.GetAbsorptionNILL() == DBL_MAX || mDataCurrent.GetAbsorptionNILL() == DBL_MAX)
- return true; // NILL is re-initialized
- }
-
- if (wlsNILL < 0. || absorptionNILL < 0.)
- return true; // let GEANT4 to reset them
-
- G4double deltaWlsNILL = mDataPrevious.GetWLSNILL() - mDataCurrent.GetWLSNILL();
- G4double deltaAbsorptionNILL = mDataPrevious.GetAbsorptionNILL() - mDataCurrent.GetAbsorptionNILL();
-
- G4double finalWlsNILL = wlsNILL - deltaWlsNILL * fSafety;
- G4double finalAbsorptionNILL = absorptionNILL - deltaAbsorptionNILL * fSafety;
-
- // prevent double counting of the probability of getting absorbed (which already estimated before transportation)
- // reset NILL again
- if (finalWlsNILL < 0. || finalAbsorptionNILL < 0.)
- return false;
-
- return true;
- }
-
- void setPostStepProc(const G4Track *track)
- {
- G4ProcessManager *pm = track->GetDefinition()->GetProcessManager();
- auto postStepProcessVector = pm->GetPostStepProcessVector();
-
- for (unsigned int np = 0; np < postStepProcessVector->entries(); np++)
- {
- auto theProcess = (*postStepProcessVector)[np];
-
- auto theType = theProcess->GetProcessType();
-
- if (theType != fOptical)
- continue;
-
- if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpBoundary)
- pOpBoundaryProc = dynamic_cast(theProcess);
- else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpAbsorption)
- pOpAbsorption = dynamic_cast(theProcess);
- else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpWLS)
- pOpWLS = dynamic_cast(theProcess);
- }
-
- fProcAssigned = true;
-
- return;
- }
-
- void reset()
- {
- mNtransport = 0.;
- mTransportUnit = 0.;
- mFiberPos = G4ThreeVector(0);
- mFiberAxis = G4ThreeVector(0);
- fKill = false;
- fTransported = false;
- mDataCurrent.reset();
- mDataPrevious.reset();
- }
-
- void print()
- {
- if (fVerbose > 1)
- {
- G4cout << G4endl;
-
- G4cout << "mDataPrevious.trackID = " << mDataPrevious.trackID;
- G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataPrevious.GetOpBoundaryStatus();
- G4cout << " | .mStepLengthInterval = " << mDataPrevious.GetStepLengthInterval() << G4endl;
-
- if (fVerbose > 2)
- {
- G4cout << " | globalPosition = (" << std::setw(9) << mDataPrevious.globalPosition.x();
- G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.y();
- G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.z() << ")" << G4endl;
-
- G4cout << " | momentumDirection = (" << std::setw(9) << mDataPrevious.momentumDirection.x();
- G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.y();
- G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.z() << ")" << G4endl;
-
- G4cout << " | polarization = (" << std::setw(9) << mDataPrevious.polarization.x();
- G4cout << "," << std::setw(9) << mDataPrevious.polarization.y();
- G4cout << "," << std::setw(9) << mDataPrevious.polarization.z() << ")" << G4endl;
-
- G4cout << " | globalTime = " << std::setw(9) << mDataPrevious.globalTime << G4endl;
- G4cout << " | WLSNILL = " << std::setw(9) << mDataPrevious.GetWLSNILL() << G4endl;
- G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataPrevious.GetAbsorptionNILL() << G4endl;
- }
-
- G4cout << "mDataCurrent.trackID = " << mDataCurrent.trackID;
- G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataCurrent.GetOpBoundaryStatus() << G4endl;
-
- if (fVerbose > 2)
- {
- G4cout << " | globalPosition = (" << std::setw(9) << mDataCurrent.globalPosition.x();
- G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.y();
- G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.z() << ")" << G4endl;
-
- G4cout << " | momentumDirection = (" << std::setw(9) << mDataCurrent.momentumDirection.x();
- G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.y();
- G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.z() << ")" << G4endl;
-
- G4cout << " | polarization = (" << std::setw(9) << mDataCurrent.polarization.x();
- G4cout << "," << std::setw(9) << mDataCurrent.polarization.y();
- G4cout << "," << std::setw(9) << mDataCurrent.polarization.z() << ")" << G4endl;
-
- G4cout << " | globalTime = " << std::setw(9) << mDataCurrent.globalTime << G4endl;
- G4cout << " | WLSNILL = " << std::setw(9) << mDataCurrent.GetWLSNILL() << G4endl;
- G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataCurrent.GetAbsorptionNILL() << G4endl;
- }
-
- G4cout << G4endl;
- }
- }
- };
-
template <>
void Geant4FSShowerModel::initialize()
{
@@ -305,8 +60,7 @@ namespace dd4hep
}
template <>
- bool Geant4FSShowerModel::check_trigger(const G4FastTrack &fasttrack)
- {
+ bool Geant4FSShowerModel::check_trigger(const G4FastTrack &fasttrack) {
if (!locals.fSwitch)
return false; // turn on/off the model
@@ -316,55 +70,12 @@ namespace dd4hep
if (locals.mDataCurrent.trackID != track->GetTrackID())
locals.reset();
- // make sure that the track does not get absorbed after transportation, as number of interaction length left is reset when doing transportation
+ // make sure that the track does not get absorbed after transportation
+ // as number of interaction length left is reset when doing transportation
if (!locals.checkNILL())
return true; // track is already transported but did not pass NILL check, attempt to reset NILL
- if (locals.fTransported)
- { // track is already transported and did pass NILL check, nothing to do
- if (locals.mFiberAxis.dot(track->GetMomentumDirection()) * locals.mFiberAxis.dot(locals.mDataCurrent.momentumDirection) < 0) // different propagation direction (e.g. mirror)
- locals.reset();
-
- return false;
- }
-
- if (!locals.checkTotalInternalReflection(track))
- return false; // nothing to do if the track has no repetitive total internal reflection
-
- auto theTouchable = track->GetTouchableHandle();
- auto solid = theTouchable->GetSolid();
-
- if (solid->GetEntityType() != "G4Tubs")
- return false; // only works for G4Tubs at the moment
-
- if (locals.fVerbose > 0)
- locals.print(); // at this point, the track should have passed all prerequisites before entering computationally heavy operations
-
- G4Tubs *tubs = static_cast(solid);
- G4double fiberLen = 2. * tubs->GetZHalfLength();
-
- locals.mFiberPos = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(G4ThreeVector(0., 0., 0.));
- locals.mFiberAxis = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformAxis(G4ThreeVector(0., 0., 1.));
-
- auto delta = locals.mDataCurrent.globalPosition - locals.mDataPrevious.globalPosition;
- locals.mTransportUnit = delta.dot(locals.mFiberAxis);
-
- // estimate the number of expected total internal reflections before reaching fiber end
- auto fiberEnd = (locals.mTransportUnit > 0.) ? locals.mFiberPos + locals.mFiberAxis * fiberLen / 2. : locals.mFiberPos - locals.mFiberAxis * fiberLen / 2.;
- auto toEnd = fiberEnd - track->GetPosition();
- G4double toEndAxis = toEnd.dot(locals.mFiberAxis);
- G4double maxTransport = std::floor(toEndAxis / locals.mTransportUnit);
- locals.mNtransport = maxTransport - locals.fSafety;
-
- if (locals.mNtransport < 1.)
- return false; // require at least n = fSafety of total internal reflections at the end
-
- if (locals.checkAbsorption(locals.mDataPrevious.GetWLSNILL(), locals.mDataCurrent.GetWLSNILL()))
- return false; // do nothing if WLS happens before reaching fiber end
- if (locals.checkAbsorption(locals.mDataPrevious.GetAbsorptionNILL(), locals.mDataCurrent.GetAbsorptionNILL()))
- locals.fKill = true; // absorbed before reaching fiber end
-
- return true;
+ return locals.check_trigger(track);
}
typedef Geant4FSShowerModel Geant4DRCFiberModel;
diff --git a/plugins/DRCaloFastSimModel.h b/plugins/DRCaloFastSimModel.h
index a8f29cb14..3b9273ce2 100644
--- a/plugins/DRCaloFastSimModel.h
+++ b/plugins/DRCaloFastSimModel.h
@@ -1,17 +1,43 @@
-#include "G4OpBoundaryProcess.hh"
+#ifndef DRCaloFastSimModel_h
+#define DRCaloFastSimModel_h
+#include "G4OpBoundaryProcess.hh"
+#include "G4GenericMessenger.hh"
+#include "G4OpAbsorption.hh"
+#include "G4OpWLS.hh"
+#include "G4Material.hh"
+#include
+#include
+#include
+#include
+#include
+#include
-struct FastFiberData
-{
+struct FastFiberData {
public:
- FastFiberData(G4int, G4double, G4double, G4double, G4ThreeVector, G4ThreeVector, G4ThreeVector, G4int status = G4OpBoundaryProcessStatus::Undefined);
- FastFiberData& operator=(const FastFiberData& theData) = default;
-
- FastFiberData(const FastFiberData& theData) = default;
- ~FastFiberData() = default;
-
+ FastFiberData(G4int id, G4double en, G4double globTime, G4double path,
+ G4ThreeVector pos, G4ThreeVector mom, G4ThreeVector pol,
+ G4int status = G4OpBoundaryProcessStatus::Undefined) {
+ trackID = id;
+ kineticEnergy = en;
+ globalTime = globTime;
+ pathLength = path;
+ globalPosition = pos;
+ momentumDirection = mom;
+ polarization = pol;
+ mOpBoundaryStatus = status;
+ mOpAbsorptionNumIntLenLeft = DBL_MAX;
+ mOpWLSNumIntLenLeft = DBL_MAX;
+ mStepLengthInterval = 0.;
+ }
+ ~FastFiberData()=default;
- void reset();
+ void reset() {
+ this->mOpBoundaryStatus = G4OpBoundaryProcessStatus::Undefined;
+ this->mOpAbsorptionNumIntLenLeft = DBL_MAX;
+ this->mOpWLSNumIntLenLeft = DBL_MAX;
+ this->mStepLengthInterval = 0.;
+ }
G4double GetAbsorptionNILL() { return mOpAbsorptionNumIntLenLeft; }
void SetAbsorptionNILL(G4double in) { mOpAbsorptionNumIntLenLeft = in; }
@@ -25,7 +51,16 @@ struct FastFiberData
G4double GetStepLengthInterval() { return mStepLengthInterval; }
void AddStepLengthInterval(G4double in) { mStepLengthInterval += in; }
- G4bool checkRepetitive(const FastFiberData, G4bool checkInterval = true);
+ G4bool checkRepetitive(const FastFiberData theData, G4bool checkInterval = true) {
+ if (this->trackID != theData.trackID)
+ return false;
+ if (this->mOpBoundaryStatus != theData.mOpBoundaryStatus)
+ return false;
+ if (checkInterval && std::abs(this->mStepLengthInterval - theData.mStepLengthInterval) > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
+ return false;
+
+ return true;
+ }
G4int trackID;
G4double kineticEnergy;
@@ -42,37 +77,282 @@ struct FastFiberData
G4double mStepLengthInterval;
};
-FastFiberData::FastFiberData(G4int id, G4double en, G4double globTime, G4double path, G4ThreeVector pos, G4ThreeVector mom, G4ThreeVector pol, G4int status)
-{
- trackID = id;
- kineticEnergy = en;
- globalTime = globTime;
- pathLength = path;
- globalPosition = pos;
- momentumDirection = mom;
- polarization = pol;
- mOpBoundaryStatus = status;
- mOpAbsorptionNumIntLenLeft = DBL_MAX;
- mOpWLSNumIntLenLeft = DBL_MAX;
- mStepLengthInterval = 0.;
-}
-
-G4bool FastFiberData::checkRepetitive(const FastFiberData theData, G4bool checkInterval)
-{
- if (this->trackID != theData.trackID)
- return false;
- if (this->mOpBoundaryStatus != theData.mOpBoundaryStatus)
+class DRCFiberModel {
+public:
+ DRCFiberModel()=default;
+ ~DRCFiberModel()=default;
+
+ G4OpBoundaryProcess* pOpBoundaryProc = nullptr;
+ G4OpAbsorption* pOpAbsorption = nullptr;
+ G4OpWLS* pOpWLS = nullptr;
+ G4bool fProcAssigned = false;
+
+ FastFiberData mDataPrevious = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0));
+ FastFiberData mDataCurrent = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0));
+
+ G4int fSafety = 1;
+ G4double mNtransport = 0.;
+ G4double mTransportUnit = 0.;
+ G4ThreeVector mFiberPos = G4ThreeVector(0);
+ G4ThreeVector mFiberAxis = G4ThreeVector(0);
+ G4bool fKill = false;
+ G4bool fTransported = false;
+ G4bool fSwitch = true;
+ G4int fVerbose = 0;
+
+ G4bool checkTotalInternalReflection(const G4Track *track) {
+ if (!fProcAssigned)
+ setPostStepProc(track); // locate OpBoundaryProcess only once
+
+ G4int theStatus = pOpBoundaryProc->GetStatus();
+
+ if (fVerbose > 1) {
+ std::cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID();
+ std::cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus;
+ std::cout << " | Track status = " << track->GetTrackStatus();
+ std::cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << std::endl;
+ }
+
+ if (track->GetTrackStatus() == fStopButAlive || track->GetTrackStatus() == fStopAndKill)
+ return false;
+
+ // accumulate step length
+ mDataCurrent.AddStepLengthInterval(track->GetStepLength());
+
+ // skip exceptional iteration with FresnelReflection
+ if (theStatus == G4OpBoundaryProcessStatus::FresnelReflection)
+ mDataCurrent.SetOpBoundaryStatus(theStatus);
+
+ // some cases have a status StepTooSmall when the reflection happens
+ // between the boundary of cladding & outer volume (outside->cladding)
+ // since the outer volume is not a G4Region
+ if (theStatus == G4OpBoundaryProcessStatus::TotalInternalReflection || theStatus == G4OpBoundaryProcessStatus::StepTooSmall) {
+ if (theStatus != G4OpBoundaryProcessStatus::TotalInternalReflection) {
+ // skip StepTooSmall if the track already has TotalInternalReflection
+ if (mDataCurrent.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection)
+ return false;
+ if (mDataPrevious.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection)
+ return false;
+ }
+
+ G4int trackID = track->GetTrackID();
+ G4double kineticEnergy = track->GetKineticEnergy();
+ G4double globalTime = track->GetGlobalTime();
+ G4double pathLength = track->GetStepLength();
+ G4ThreeVector globalPosition = track->GetPosition();
+ G4ThreeVector momentumDirection = track->GetMomentumDirection();
+ G4ThreeVector polarization = track->GetPolarization();
+
+ auto candidate = FastFiberData(trackID, kineticEnergy, globalTime, pathLength,
+ globalPosition, momentumDirection, polarization, theStatus);
+
+ if (pOpAbsorption != nullptr)
+ candidate.SetAbsorptionNILL(pOpAbsorption->GetNumberOfInteractionLengthLeft());
+ if (pOpWLS != nullptr)
+ candidate.SetWLSNILL(pOpWLS->GetNumberOfInteractionLengthLeft());
+
+ G4bool repetitive = false;
+
+ if (candidate.checkRepetitive(mDataCurrent, false) && mDataCurrent.checkRepetitive(mDataPrevious))
+ repetitive = true;
+
+ mDataPrevious = mDataCurrent;
+ mDataCurrent = candidate;
+
+ return repetitive;
+ }
+
return false;
- if (checkInterval && std::abs(this->mStepLengthInterval - theData.mStepLengthInterval) > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
+ } // checkTotalInternalReflection
+
+ G4bool checkAbsorption(const G4double prevNILL, const G4double currentNILL) {
+ if (prevNILL < 0. || currentNILL < 0.)
+ return false; // the number of interaction length left has to be reset
+ if (prevNILL == currentNILL)
+ return false; // no absorption
+ if (prevNILL == DBL_MAX || currentNILL == DBL_MAX)
+ return false; // NILL is re-initialized
+
+ G4double deltaNILL = prevNILL - currentNILL;
+
+ if (currentNILL - deltaNILL * (mNtransport + fSafety) < 0.)
+ return true; // absorbed before reaching fiber end
+
return false;
+ } // checkAbsorption
+
+ G4bool checkNILL() {
+ if (!fTransported)
+ return true; // do nothing if the track is not already transported
+
+ G4double wlsNILL = DBL_MAX;
+ G4double absorptionNILL = DBL_MAX;
+
+ if (pOpWLS != nullptr) {
+ wlsNILL = pOpWLS->GetNumberOfInteractionLengthLeft();
+ if (mDataPrevious.GetWLSNILL() == DBL_MAX || mDataCurrent.GetWLSNILL() == DBL_MAX)
+ return true; // NILL is re-initialized
+ }
+
+ if (pOpAbsorption != nullptr) {
+ absorptionNILL = pOpAbsorption->GetNumberOfInteractionLengthLeft();
+ if (mDataPrevious.GetAbsorptionNILL() == DBL_MAX || mDataCurrent.GetAbsorptionNILL() == DBL_MAX)
+ return true; // NILL is re-initialized
+ }
+
+ if (wlsNILL < 0. || absorptionNILL < 0.)
+ return true; // let GEANT4 to reset them
+
+ G4double deltaWlsNILL = mDataPrevious.GetWLSNILL() - mDataCurrent.GetWLSNILL();
+ G4double deltaAbsorptionNILL = mDataPrevious.GetAbsorptionNILL() - mDataCurrent.GetAbsorptionNILL();
+
+ G4double finalWlsNILL = wlsNILL - deltaWlsNILL * fSafety;
+ G4double finalAbsorptionNILL = absorptionNILL - deltaAbsorptionNILL * fSafety;
+
+ // prevent double counting of the probability of getting absorbed
+ // (which already estimated before transportation)
+ // reset NILL again
+ if (finalWlsNILL < 0. || finalAbsorptionNILL < 0.)
+ return false;
+
+ return true;
+ } // checkNILL
+
+ bool check_trigger(const G4Track* track) {
+ if (fTransported) {
+ // different propagation direction (e.g. mirror)
+ if (mFiberAxis.dot(track->GetMomentumDirection()) * mFiberAxis.dot(mDataCurrent.momentumDirection) < 0)
+ reset();
+
+ // track is already transported and did pass NILL check, nothing to do
+ return false;
+ }
+
+ if (!checkTotalInternalReflection(track))
+ return false; // nothing to do if the track has no repetitive total internal reflection
+
+ auto theTouchable = track->GetTouchableHandle();
+ auto solid = theTouchable->GetSolid();
+
+ if (fVerbose > 0)
+ print(); // at this point, the track should have passed all prerequisites before entering computationally heavy operations
+
+ if (solid->GetEntityType() != "G4Tubs")
+ return false; // only works for G4Tubs at the moment
+
+ G4Tubs* tubs = static_cast(solid);
+ G4double fiberLen = 2. * tubs->GetZHalfLength();
+
+ mFiberPos = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(G4ThreeVector(0., 0., 0.));
+ mFiberAxis = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformAxis(G4ThreeVector(0., 0., 1.));
+
+ auto delta = mDataCurrent.globalPosition - mDataPrevious.globalPosition;
+ mTransportUnit = delta.dot(mFiberAxis);
+
+ // estimate the number of expected total internal reflections before reaching fiber end
+ auto fiberEnd = (mTransportUnit > 0.) ? mFiberPos + mFiberAxis * fiberLen / 2. : mFiberPos - mFiberAxis * fiberLen / 2.;
+ auto toEnd = fiberEnd - track->GetPosition();
+ G4double toEndAxis = toEnd.dot(mFiberAxis);
+ G4double maxTransport = std::floor(toEndAxis / mTransportUnit);
+ mNtransport = maxTransport - fSafety;
+
+ if (mNtransport < 1.)
+ return false; // require at least n = fSafety of total internal reflections at the end
+
+ if (checkAbsorption(mDataPrevious.GetWLSNILL(), mDataCurrent.GetWLSNILL()))
+ return false; // do nothing if WLS happens before reaching fiber end
+ if (checkAbsorption(mDataPrevious.GetAbsorptionNILL(), mDataCurrent.GetAbsorptionNILL()))
+ fKill = true; // absorbed before reaching fiber end
+
+ return true;
+ }
+
+ void setPostStepProc(const G4Track *track) {
+ G4ProcessManager *pm = track->GetDefinition()->GetProcessManager();
+ auto postStepProcessVector = pm->GetPostStepProcessVector();
+
+ for (unsigned int np = 0; np < postStepProcessVector->entries(); np++) {
+ auto theProcess = (*postStepProcessVector)[np];
+
+ auto theType = theProcess->GetProcessType();
+
+ if (theType != fOptical)
+ continue;
+
+ if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpBoundary)
+ pOpBoundaryProc = dynamic_cast(theProcess);
+ else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpAbsorption)
+ pOpAbsorption = dynamic_cast(theProcess);
+ else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpWLS)
+ pOpWLS = dynamic_cast(theProcess);
+ }
+
+ fProcAssigned = true;
+
+ return;
+ } // setPostStepProc
+
+ void reset() {
+ mNtransport = 0.;
+ mTransportUnit = 0.;
+ mFiberPos = G4ThreeVector(0);
+ mFiberAxis = G4ThreeVector(0);
+ fKill = false;
+ fTransported = false;
+ mDataCurrent.reset();
+ mDataPrevious.reset();
+ }
+
+ void print() {
+ if (fVerbose > 1) {
+ G4cout << G4endl;
+
+ G4cout << "mDataPrevious.trackID = " << mDataPrevious.trackID;
+ G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataPrevious.GetOpBoundaryStatus();
+ G4cout << " | .mStepLengthInterval = " << mDataPrevious.GetStepLengthInterval() << G4endl;
+
+ if (fVerbose > 2) {
+ G4cout << " | globalPosition = (" << std::setw(9) << mDataPrevious.globalPosition.x();
+ G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.y();
+ G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.z() << ")" << G4endl;
+
+ G4cout << " | momentumDirection = (" << std::setw(9) << mDataPrevious.momentumDirection.x();
+ G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.y();
+ G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.z() << ")" << G4endl;
+
+ G4cout << " | polarization = (" << std::setw(9) << mDataPrevious.polarization.x();
+ G4cout << "," << std::setw(9) << mDataPrevious.polarization.y();
+ G4cout << "," << std::setw(9) << mDataPrevious.polarization.z() << ")" << G4endl;
+
+ G4cout << " | globalTime = " << std::setw(9) << mDataPrevious.globalTime << G4endl;
+ G4cout << " | WLSNILL = " << std::setw(9) << mDataPrevious.GetWLSNILL() << G4endl;
+ G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataPrevious.GetAbsorptionNILL() << G4endl;
+ }
+
+ G4cout << "mDataCurrent.trackID = " << mDataCurrent.trackID;
+ G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataCurrent.GetOpBoundaryStatus() << G4endl;
+
+ if (fVerbose > 2) {
+ G4cout << " | globalPosition = (" << std::setw(9) << mDataCurrent.globalPosition.x();
+ G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.y();
+ G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.z() << ")" << G4endl;
+
+ G4cout << " | momentumDirection = (" << std::setw(9) << mDataCurrent.momentumDirection.x();
+ G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.y();
+ G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.z() << ")" << G4endl;
+
+ G4cout << " | polarization = (" << std::setw(9) << mDataCurrent.polarization.x();
+ G4cout << "," << std::setw(9) << mDataCurrent.polarization.y();
+ G4cout << "," << std::setw(9) << mDataCurrent.polarization.z() << ")" << G4endl;
+
+ G4cout << " | globalTime = " << std::setw(9) << mDataCurrent.globalTime << G4endl;
+ G4cout << " | WLSNILL = " << std::setw(9) << mDataCurrent.GetWLSNILL() << G4endl;
+ G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataCurrent.GetAbsorptionNILL() << G4endl;
+ }
- return true;
-}
+ G4cout << G4endl;
+ }
+ } // print
+}; // end of the class
-void FastFiberData::reset()
-{
- this->mOpBoundaryStatus = G4OpBoundaryProcessStatus::Undefined;
- this->mOpAbsorptionNumIntLenLeft = DBL_MAX;
- this->mOpWLSNumIntLenLeft = DBL_MAX;
- this->mStepLengthInterval = 0.;
-}
+#endif
diff --git a/plugins/FiberDRCaloSDAction.cpp b/plugins/FiberDRCaloSDAction.cpp
index 1d946ff21..37dc55d9c 100644
--- a/plugins/FiberDRCaloSDAction.cpp
+++ b/plugins/FiberDRCaloSDAction.cpp
@@ -1,6 +1,4 @@
// DD4hep Framework include files
-#include "DD4hep/Version.h"
-#include "DD4hep/Objects.h"
#include "DD4hep/Segmentations.h"
#include "DDG4/Geant4Random.h"
#include "DDG4/Geant4SensDetAction.inl"
@@ -9,8 +7,22 @@
#include "G4OpticalPhoton.hh"
// k4geo Framework include files
-
#include "FiberDRCaloSDAction.h"
+#include "DRCaloFastSimModel.h"
+
+// Geant4 include files
+#include "G4SystemOfUnits.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4Step.hh"
+#include "G4TouchableHistory.hh"
+#include "G4ThreeVector.hh"
+#include "G4HCofThisEvent.hh"
+#include "G4SDManager.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTypes.hh"
+#include "G4OpticalPhoton.hh"
+#include "G4OpBoundaryProcess.hh"
+#include "G4OpProcessSubType.hh"
#if DD4HEP_VERSION_GE(1, 21)
#define GEANT4_CONST_STEP const
@@ -29,15 +41,15 @@ namespace dd4hep
/*
* Geant4SensitiveAction sensitive detector for the Dual-readout calorimeter
*
- * \author Sungwon Kim
+ * \author Sungwon Kim, Sanghyun Ko
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
- struct DRCData
- {
- Geant4HitCollection *fHitCollection;
- G4int fHCID;
+ struct DRCData {
+ public:
+ bool skipScint = true;
+ DRCFiberModel fastfiber;
G4int fWavBin;
G4int fTimeBin;
@@ -47,63 +59,136 @@ namespace dd4hep
G4float fTimeEnd;
G4float fWavlenStep;
G4float fTimeStep;
- static const int fArraySize = 24;
- double fGraph_X[fArraySize] ={1.37760 * dd4hep::eV,
- 1.45864 * dd4hep::eV,
- 1.54980 * dd4hep::eV,
- 1.65312 * dd4hep::eV,
- 1.71013 * dd4hep::eV,
- 1.77120 * dd4hep::eV,
- 1.83680 * dd4hep::eV,
- 1.90745 * dd4hep::eV,
- 1.98375 * dd4hep::eV,
- 2.06640 * dd4hep::eV,
- 2.10143 * dd4hep::eV,
- 2.13766 * dd4hep::eV,
- 2.17516 * dd4hep::eV,
- 2.21400 * dd4hep::eV,
- 2.25426 * dd4hep::eV,
- 2.29600 * dd4hep::eV,
- 2.33932 * dd4hep::eV,
- 2.38431 * dd4hep::eV,
- 2.43106 * dd4hep::eV,
- 2.47968 * dd4hep::eV,
- 2.53029 * dd4hep::eV,
- 2.58300 * dd4hep::eV,
- 2.63796 * dd4hep::eV,
- 2.69531 * dd4hep::eV};
- double fGraph_Y[fArraySize] ={0.903 * 0.903,
- 0.903 * 0.903,
- 0.903 * 0.903,
- 0.903 * 0.903,
- 0.903 * 0.903,
- 0.903 * 0.903,
- 0.902 * 0.902,
- 0.901 * 0.901,
- 0.898 * 0.898,
- 0.895 * 0.895,
- 0.893 * 0.893,
- 0.891 * 0.891,
- 0.888 * 0.888,
- 0.883 * 0.883,
- 0.87 * 0.87,
- 0.838 * 0.838,
- 0.76 * 0.76,
- 0.62 * 0.62,
- 0.488 * 0.488,
- 0.345 * 0.345,
- 0.207 * 0.207,
- 0.083 * 0.083,
- 0.018 * 0.018,
- 0.};
-
- G4double wavToE(G4double wav) { return CLHEP::h_Planck * CLHEP::c_light / wav; }
-
- int findWavBin(G4double en)
- {
+
+ private:
+ // from 900 nm to 460 nm
+ const std::vector fGraph_X = {
+ 1.37760 * CLHEP::eV,
+ 1.45864 * CLHEP::eV,
+ 1.54980 * CLHEP::eV,
+ 1.65312 * CLHEP::eV,
+ 1.71013 * CLHEP::eV,
+ 1.77120 * CLHEP::eV,
+ 1.83680 * CLHEP::eV,
+ 1.90745 * CLHEP::eV,
+ 1.98375 * CLHEP::eV,
+ 2.06640 * CLHEP::eV,
+
+ 2.10143 * CLHEP::eV,
+ 2.13766 * CLHEP::eV,
+ 2.17516 * CLHEP::eV,
+ 2.21400 * CLHEP::eV,
+ 2.25426 * CLHEP::eV,
+ 2.29600 * CLHEP::eV,
+ 2.33932 * CLHEP::eV,
+ 2.38431 * CLHEP::eV,
+ 2.43106 * CLHEP::eV,
+ 2.47968 * CLHEP::eV,
+
+ 2.53029 * CLHEP::eV,
+ 2.58300 * CLHEP::eV,
+ 2.63796 * CLHEP::eV,
+ 2.69531 * CLHEP::eV,
+ 2.75520 * CLHEP::eV,
+ 2.81782 * CLHEP::eV,
+ 2.88335 * CLHEP::eV,
+ 2.95200 * CLHEP::eV,
+ 3.09960 * CLHEP::eV,
+ 3.54241 * CLHEP::eV,
+
+ 4.13281 * CLHEP::eV
+ };
+ // filter efficiency of the Kodak Wratten No.9 filter
+ const std::vector fKodakEff = {
+ 0.903,
+ 0.903,
+ 0.903,
+ 0.903,
+ 0.903,
+ 0.903,
+ 0.902,
+ 0.901,
+ 0.898,
+ 0.895,
+
+ 0.893,
+ 0.891,
+ 0.888,
+ 0.883,
+ 0.87,
+ 0.838,
+ 0.76,
+ 0.62,
+ 0.488,
+ 0.345,
+
+ 0.207,
+ 0.083,
+ 0.018,
+ 0.,
+ 0.,
+ 0.,
+ 0.,
+ 0.,
+ 0.,
+ 0.,
+
+ 0.
+ };
+
+ // SiPM efficiency Hamamatsu S14160-1310PS
+ // TODO migrate this part to the digitization step!
+ // Note: Ideally, this should be part of the digitization step.
+ // (not the simulation step)
+ // But, to do this, we need to store the distribution
+ // of the optical photon wavelength.
+ // While we can develop another feature to enable this,
+ // let's emulate the SiPM efficiency in the simulation step for now.
+ // We just need a working code and without this,
+ // the number of Cherenkov photon will be order of magnitude higher
+ const std::vector fSipmEff = {
+ 0.02,
+ 0.025,
+ 0.045,
+ 0.06,
+ 0.0675,
+ 0.075,
+ 0.0925,
+ 0.11,
+ 0.125,
+ 0.14,
+
+ 0.146,
+ 0.152,
+ 0.158,
+ 0.164,
+ 0.17,
+ 0.173,
+ 0.176,
+ 0.178,
+ 0.179,
+ 0.18,
+
+ 0.181,
+ 0.182,
+ 0.183,
+ 0.184,
+ 0.18,
+ 0.173,
+ 0.166,
+ 0.158,
+ 0.15,
+ 0.12,
+
+ 0.05
+ };
+
+ public:
+ G4double wavToE(G4double wav) const { return CLHEP::h_Planck * CLHEP::c_light / wav; }
+
+ int findWavBin(G4double en) const {
int i = 0;
- for (; i < fWavBin + 1; i++)
- {
+ for (; i < fWavBin + 1; i++) {
if (en < wavToE((fWavlenStart - static_cast(i) * fWavlenStep) * CLHEP::nm))
break;
}
@@ -111,11 +196,9 @@ namespace dd4hep
return fWavBin + 1 - i;
}
- int findTimeBin(G4double stepTime)
- {
+ int findTimeBin(G4double stepTime) const {
int i = 0;
- for (; i < fTimeBin + 1; i++)
- {
+ for (; i < fTimeBin + 1; i++) {
if (stepTime < ((fTimeStart + static_cast(i) * fTimeStep) * CLHEP::ns))
break;
}
@@ -123,37 +206,31 @@ namespace dd4hep
return i;
}
- // Linear interpolation function for calculating the efficiency of yellow filter, used in rejectedByYellowFilter
- double getFilterEff(G4double G4energy) {
- double energy = (G4energy * (dd4hep::MeV / CLHEP::MeV)); // Convert G4 MeV to dd4hep MeV
- if (energy <= fGraph_X[0]) // If the photon energy <= 1.37760 eV, than return maximum filter efficiency
- return fGraph_Y[0];
-
- for(int idx = 0; idx < fArraySize; ++idx) {
- if (energy <= fGraph_X[idx]) {
- double x1 = fGraph_X[idx - 1];
- double x2 = fGraph_X[idx];
- double y1 = fGraph_Y[idx - 1];
- double y2 = fGraph_Y[idx];
-
- return (y1 + ((y2 - y1) / (x2 - x1))*(energy - x1)); // return linear interpolated filter efficiency
+ // Linear interpolation function for calculating the efficiency of yellow filter
+ // used in rejectedByYellowFilter
+ double getFilterEff(const std::vector& yarray, const G4double G4energy) const {
+ // If the photon energy <= 1.37760 eV, than return maximum filter efficiency
+ if (G4energy <= fGraph_X.at(0))
+ return yarray.at(0);
+
+ for(unsigned idx = 1; idx < yarray.size(); ++idx) {
+ if (G4energy <= fGraph_X.at(idx)) {
+ double x1 = fGraph_X.at(idx-1);
+ double x2 = fGraph_X.at(idx);
+ double y1 = yarray.at(idx-1);
+ double y2 = yarray.at(idx);
+
+ // return linear interpolated filter efficiency
+ return (y1 + ((y2 - y1) / (x2 - x1))*(G4energy - x1));
}
}
- // This should not happen
- std::cout << "Error in Yellow filter efficiency with photon energy : " << energy << " MeV" << std::endl;
- std::cout << "Cannot find corresponding filter efficiency" << std::endl;
return 0.;
}
// If true, then the photon is rejected by yellow filter
- bool rejectedByYellowFilter(G4double G4energy, double rndVal)
- {
- double energy = (G4energy * (dd4hep::MeV / CLHEP::MeV)); // Convert G4 MeV to dd4hep MeV
- if ( energy >= fGraph_X[fArraySize-1] ) // Photon w/ E > 2.69531 eV rejected
- return true;
-
- double FilterEff = getFilterEff(G4energy); // Get efficiency of filter using photon's energy
+ bool rejectedByYellowFilter(G4double G4energy, double rndVal) const {
+ const double FilterEff = getFilterEff(fKodakEff,G4energy); // Get efficiency of filter using photon's energy
// filter efficiency == probability of photon accepted by filter
// == Probability of random value (from uniform distribution with range 0 ~ 1) smaller than filter efficiency
@@ -161,97 +238,299 @@ namespace dd4hep
return (rndVal > FilterEff);
}
+ // check sipm efficiency
+ // TODO migrate this to the digitization step
+ bool rejectedBySiPM(double G4energy, double rndVal) const {
+ return (rndVal > getFilterEff(fSipmEff,G4energy));
+ }
+
DRCData()
- : fHitCollection(0), fHCID(-1), fWavBin(120), fTimeBin(650),
- fWavlenStart(900.), fWavlenEnd(300.), fTimeStart(5.), fTimeEnd(70.)
+ : fWavBin(120), fTimeBin(650), fWavlenStart(900.),
+ fWavlenEnd(300.), fTimeStart(5.), fTimeEnd(70.)
{
fWavlenStep = (fWavlenStart - fWavlenEnd) / (float)fWavBin;
fTimeStep = (fTimeEnd - fTimeStart) / (float)fTimeBin;
}
- };
+ }; // struct DRCData
- /// Define collections created by this sensitive action object
template <>
- void Geant4SensitiveAction::defineCollections()
+ Geant4SensitiveAction::Geant4SensitiveAction(Geant4Context* ctxt,
+ const std::string& nam,
+ DetElement det,
+ Detector& desc)
+ : Geant4Sensitive(ctxt,nam,det,desc), m_collectionName(), m_collectionID(0)
{
+ declareProperty("skipScint", m_userData.skipScint = true);
+ declareProperty("ReadoutName", m_readoutName);
+ declareProperty("CollectionName", m_collectionName);
+ initialize();
+ InstanceCount::increment(this);
+
+ m_userData.fastfiber.fSafety = 1;
+ m_userData.fastfiber.fVerbose = 0;
+ }
+
+ /// Define collections created by this sensitive action object
+ template <>
+ void Geant4SensitiveAction::defineCollections() {
std::string readout_name = m_sensitive.readout().name();
m_collectionID = defineCollection(m_sensitive.readout().name());
std::cout << "defineCollection Geant4DRCalorimeter readout_name : " << readout_name << std::endl;
std::cout << "defineCollection Geant4DRCalorimeter m_collectionID : " << m_collectionID << std::endl;
+
+ if (m_userData.skipScint) {
+ defineCollection(std::string(m_sensitive.readout().name())+"_scint");
+ std::cout << "defineCollection Geant4Calorimeter readout_name : " << readout_name + "_scint" << std::endl;
+ std::cout << "defineCollection Geant4Calorimeter m_collectionID : " << m_collectionID + 1 << std::endl;
+ }
}
/// Method for generating hit(s) using the information of G4Step object.
template <>
- G4bool Geant4SensitiveAction::process(G4Step GEANT4_CONST_STEP *step, G4TouchableHistory*)
- {
- if (step->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
- return false;
-
- typedef Geant4DRCalorimeter::Hit Hit;
-
- Geant4StepHandler h(step);
- Geant4HitCollection *coll = collection(m_collectionID);
- HitContribution contrib = Hit::extractContribution(step);
-
- auto theTouchable = step->GetPostStepPoint()->GetTouchable();
- dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager();
- dd4hep::VolumeID volID = volMgr.volumeID(theTouchable);
- G4ThreeVector global = step->GetPostStepPoint()->GetPosition();
- G4ThreeVector local = theTouchable->GetHistory()->GetTopTransform().TransformPoint(global);
- // MoveUpHistory(GetHistoryDepth - 2) -> tower touchable, local position
- dd4hep::Position loc(local.x() * dd4hep::millimeter / CLHEP::millimeter, local.y() * dd4hep::millimeter / CLHEP::millimeter, local.z() * dd4hep::millimeter / CLHEP::millimeter);
- dd4hep::Position glob(global.x() * dd4hep::millimeter / CLHEP::millimeter, global.y() * dd4hep::millimeter / CLHEP::millimeter, global.z() * dd4hep::millimeter / CLHEP::millimeter);
-
- auto cID = m_segmentation->cellID(loc, glob, volID); // This returns cID corresponding to SiPM Wafer
- Hit *hit = coll->find(CellIDCompare(cID));
-
- G4double hitTime = step->GetPostStepPoint()->GetGlobalTime();
- G4double energy = step->GetTrack()->GetTotalEnergy();
-
- // Apply yellow filter here
- // Get random number from (0~1) uniform distribution
- // If the photon is from Scint. process, calculate the filter eff. based on its energy
- // If (random number) > (Eff), the photon is rejected by yellow filter (= do not make hit (or count photon) for that photon)
- dd4hep::DDSegmentation::VolumeID ceren = static_cast(m_segmentation->decoder()->get(cID, "c"));
- bool IsCeren = static_cast(ceren);
- if (!(IsCeren))
- {
- Geant4Event& evt = context()->event();
+ G4bool Geant4SensitiveAction::process(G4Step GEANT4_CONST_STEP *step, G4TouchableHistory *) {
+ // optical photons traveling through the cladding is not interesting
+ // but consume lots of computation power
+ // let's kill optical photons inside the cladding whose status is not StepTooSmall
+ if ( step->GetPreStepPoint() &&
+ step->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetNoDaughters()==1 &&
+ step->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition() ) {
+ // 1e-9 is the default tolerance
+ // for the warnings, see https://geant4-forum.web.cern.ch/t/error-occurs-when-an-optical-photon-hits-the-edge-of-a-cubic-scintillator/8748
+ // we're not interested in the optical photon position
+ // but only the number of photons (and timing)
+ // so we assume it's fine to ignore the warnings
+ if (step->GetStepLength() > 1e-8 * CLHEP::mm)
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+ }
+ // the fiber itself is the SD
+ // we skip the scintillation process and only account for the Cherenkov process
+ // remember to turn off the scintillation process!
+ if (m_userData.skipScint) {
+ // we need to move the touchable to the tower to retrieve volID
+ auto* touchable = const_cast(step->GetPostStepPoint()->GetTouchable());
+ const auto* logicalVol = touchable->GetVolume()->GetLogicalVolume();
+
+ // we're not interested in the world or assembly volume
+ if (touchable->GetHistoryDepth() < 2)
+ return false;
+
+ // we're only interested in the fiber core
+ if ( logicalVol->GetNoDaughters()!=0 )
+ return false;
+
+ // now let's make the touchable points to the tower
+ // world -> assembly -> tower
+ touchable->MoveUpHistory( touchable->GetHistoryDepth()-2 );
+
+ // now the touchable is the tower
+ // get local position and volumeID
+ // dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager();
+ // dd4hep::VolumeID volID = volMgr.volumeID(touchable);
+ // note the above method started to throw warnings since the PR DD4hep#1390
+ // so we switch to the copy number
+
+ // trick: we retrieve the tower's volID by using the fact that
+ // in the DRconstructor.cpp, the first 32bits of the volID
+ // contain the information relevant to the tower
+ // and this 32bit integer becomes the tower's copy number
+ auto copyNo = touchable->GetCopyNumber();
+
+ // ideally the below operation should be GridDRcalo_k4geo::convertFirst32to64
+ // but this is trivial, so let's avoid doing dynamic cast
+ dd4hep::VolumeID volID = static_cast(copyNo);
+ G4ThreeVector global = step->GetPostStepPoint()->GetPosition();
+ G4ThreeVector local = touchable->GetHistory()->GetTopTransform().TransformPoint(global);
+
+ // convert G4 position to dd4hep position
+ dd4hep::Position loc(local.x() * dd4hep::millimeter / CLHEP::millimeter, local.y() * dd4hep::millimeter / CLHEP::millimeter, local.z() * dd4hep::millimeter / CLHEP::millimeter);
+ dd4hep::Position glob(global.x() * dd4hep::millimeter / CLHEP::millimeter, global.y() * dd4hep::millimeter / CLHEP::millimeter, global.z() * dd4hep::millimeter / CLHEP::millimeter);
+
+ // retrieve cellID and distinguish Cherenkov & scintillation fibers
+ auto cID = m_segmentation->cellID(loc, glob, volID);
+ auto ceren = static_cast(m_segmentation->decoder()->get(cID, "c"));
+ bool IsCeren = static_cast(ceren);
+
+ Geant4HitCollection* coll = collection(m_collectionID);
+
+ if (IsCeren) {
+ // Cherenkov fiber
+ // skip anything else than optical photons
+ if (step->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
+ return false;
+
+ const auto* track = step->GetTrack();
+
+ // reset when moving to the next track
+ if (m_userData.fastfiber.mDataCurrent.trackID != track->GetTrackID())
+ m_userData.fastfiber.reset();
+
+ // need repetitive total internal reflection
+ if (!m_userData.fastfiber.check_trigger(track))
+ return false;
+
+ // absorption
+ if (m_userData.fastfiber.fKill) {
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+
+ return false;
+ }
+
+ // backward transportation
+ if (m_userData.fastfiber.mTransportUnit < 0.) {
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+
+ return false;
+ }
+
+ // for timing measurement
+ double timeUnit = m_userData.fastfiber.mDataCurrent.globalTime - m_userData.fastfiber.mDataPrevious.globalTime;
+ double timeShift = timeUnit * m_userData.fastfiber.mNtransport;
+
+ // check wheter the optical photon will be rejected by the SiPM
+ // TODO migrate this to the digitization step
+ Geant4Event& evt = context()->event();
+ Geant4Random& rnd = evt.random(); // Get random generator
+ double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1]
+ G4double energy = step->GetTrack()->GetTotalEnergy();
+
+ if (m_userData.rejectedBySiPM(energy, rndVal)) {
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+
+ return false;
+ }
+
+ // default hit (optical photon count)
+ Geant4DRCalorimeter::Hit* drHit = coll->find(CellIDCompare(cID));
+
+ if (!drHit) {
+ drHit = new Geant4DRCalorimeter::Hit(m_userData.fWavlenStep, m_userData.fTimeStep);
+ drHit->cellID = cID;
+ drHit->SetSiPMnum(cID);
+ drHit->SetTimeStart(m_userData.fTimeStart);
+ drHit->SetTimeEnd(m_userData.fTimeEnd);
+ drHit->SetWavlenMax(m_userData.fWavlenStart);
+ drHit->SetWavlenMin(m_userData.fWavlenEnd);
+ coll->add(cID, drHit);
+ }
+
+ // everything should be in the G4 unit
+ // (approximate) timing at the end of the fiber
+ G4double hitTime = step->GetPostStepPoint()->GetGlobalTime() + timeShift;
+
+ drHit->photonCount();
+ int wavBin = m_userData.findWavBin(energy);
+ drHit->CountWavlenSpectrum(wavBin);
+ int timeBin = m_userData.findTimeBin(hitTime);
+ drHit->CountTimeStruct(timeBin);
+
+ drHit->position = glob;
+
+ // finally kill optical photon
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+
+ return true;
+ } else {
+ // scintillation calo hit
+
+ // assume nPhoton_scint >> nPhoton_cherenkov
+ // kill optical photons (from Cherenkov process)
+ if (step->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
+ step->GetTrack()->SetTrackStatus(fStopAndKill);
+
+ return false;
+ }
+
+ // copy-paste of the dd4hep scintillation calorimeter SD
+ Geant4HitCollection* coll_scint = collection(m_collectionID+1);
+ Geant4Calorimeter::Hit* caloHit = coll_scint->find(CellIDCompare(cID));
+ HitContribution contrib = Geant4Calorimeter::Hit::extractContribution(step,true);
+
+ if ( !caloHit ) {
+ caloHit = new Geant4Calorimeter::Hit(glob);
+ caloHit->cellID = cID;
+ coll_scint->add(cID, caloHit);
+ }
+
+ caloHit->truth.emplace_back(contrib);
+ caloHit->energyDeposit += contrib.deposit;
+
+ Geant4StepHandler h(step);
+ mark(h.track);
+
+ return true;
+ } // !IsCeren
+ } else {
+ // no skipping optical photon propagation
+ // SiPM wafers are SD in this case
+ if (step->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
+ return false;
+
+ typedef Geant4DRCalorimeter::Hit Hit;
+
+ Geant4HitCollection *coll = collection(m_collectionID);
+
+ auto* touchable = step->GetPostStepPoint()->GetTouchable();
+ dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager();
+ dd4hep::VolumeID volID = volMgr.volumeID(touchable);
+ G4ThreeVector global = step->GetPostStepPoint()->GetPosition();
+ G4ThreeVector local = touchable->GetHistory()->GetTopTransform().TransformPoint(global);
+
+ dd4hep::Position loc(local.x() * dd4hep::millimeter / CLHEP::millimeter, local.y() * dd4hep::millimeter / CLHEP::millimeter, local.z() * dd4hep::millimeter / CLHEP::millimeter);
+ dd4hep::Position glob(global.x() * dd4hep::millimeter / CLHEP::millimeter, global.y() * dd4hep::millimeter / CLHEP::millimeter, global.z() * dd4hep::millimeter / CLHEP::millimeter);
+
+ auto cID = m_segmentation->cellID(loc, glob, volID); // This returns cID corresponding to SiPM Wafer
+ Hit *hit = coll->find(CellIDCompare(cID));
+
+ G4double hitTime = step->GetPostStepPoint()->GetGlobalTime();
+ G4double energy = step->GetTrack()->GetTotalEnergy();
+
+ // Apply yellow filter here
+ // Get random number from (0~1) uniform distribution
+ // If the photon is from Scint. process, calculate the filter eff. based on its energy
+ // If (random number) > (Eff), the photon is rejected by yellow filter (= do not make hit (or count photon) for that photon)
+ dd4hep::DDSegmentation::VolumeID ceren = static_cast(m_segmentation->decoder()->get(cID, "c"));
+ bool IsCeren = static_cast(ceren);
+
+ Geant4Event& evt = context()->event();
Geant4Random& rnd = evt.random(); // Get random generator
double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1]
- if (m_userData.rejectedByYellowFilter(energy, rndVal))
- return true;
- }
- if (!hit)
- {
- hit = new Geant4DRCalorimeter::Hit(m_userData.fWavlenStep, m_userData.fTimeStep);
- hit->cellID = cID;
- hit->SetSiPMnum(cID);
- hit->SetTimeStart(m_userData.fTimeStart);
- hit->SetTimeEnd(m_userData.fTimeEnd);
- hit->SetWavlenMax(m_userData.fWavlenStart);
- hit->SetWavlenMin(m_userData.fWavlenEnd);
- coll->add(hit);
- }
+ if (!IsCeren) {
+ if (m_userData.rejectedByYellowFilter(energy, rndVal))
+ return false;
+ }
- hit->photonCount();
- int wavBin = m_userData.findWavBin(energy);
- hit->CountWavlenSpectrum(wavBin);
- int timeBin = m_userData.findTimeBin(hitTime);
- hit->CountTimeStruct(timeBin);
+ // check wheter the optical photon will be rejected by the SiPM
+ // TODO migrate this to the digitization step
+ if (m_userData.rejectedBySiPM(energy, rndVal))
+ return false;
+
+ if (!hit) {
+ hit = new Geant4DRCalorimeter::Hit(m_userData.fWavlenStep, m_userData.fTimeStep);
+ hit->cellID = cID;
+ hit->SetSiPMnum(cID);
+ hit->SetTimeStart(m_userData.fTimeStart);
+ hit->SetTimeEnd(m_userData.fTimeEnd);
+ hit->SetWavlenMax(m_userData.fWavlenStart);
+ hit->SetWavlenMin(m_userData.fWavlenEnd);
+ coll->add(cID,hit);
+ }
- hit->position = glob;
- hit->energyDeposit += contrib.deposit;
- hit->truth.emplace_back(contrib);
+ hit->photonCount();
+ int wavBin = m_userData.findWavBin(energy);
+ hit->CountWavlenSpectrum(wavBin);
+ int timeBin = m_userData.findTimeBin(hitTime);
+ hit->CountTimeStruct(timeBin);
- return true;
- }
+ hit->position = glob;
+
+ return true;
+ } // !skipScint
+ } // Geant4SensitiveAction::process
typedef Geant4SensitiveAction DRCaloSDAction;
}
}
#include "DDG4/Factories.h"
-
DECLARE_GEANT4SENSITIVE(DRCaloSDAction)
diff --git a/plugins/FiberDRCaloSDAction.h b/plugins/FiberDRCaloSDAction.h
index 11c4e04ad..0815a2558 100644
--- a/plugins/FiberDRCaloSDAction.h
+++ b/plugins/FiberDRCaloSDAction.h
@@ -149,4 +149,5 @@ namespace dd4hep
}
}
+
#endif
diff --git a/plugins/Geant4Output2EDM4hep_DRC.cpp b/plugins/Geant4Output2EDM4hep_DRC.cpp
index fad2fde9f..864b9b0b4 100644
--- a/plugins/Geant4Output2EDM4hep_DRC.cpp
+++ b/plugins/Geant4Output2EDM4hep_DRC.cpp
@@ -75,7 +75,7 @@ namespace dd4hep {
int m_eventNo { 0 };
int m_eventNumberOffset { 0 };
bool m_filesByRun { false };
-
+
/// Data conversion interface for MC particles to EDM4hep format
void saveParticles(Geant4ParticleMap* particles);
/// Store the metadata frame with e.g. the cellID encoding strings
@@ -111,7 +111,7 @@ namespace dd4hep {
}
}
};
-
+
template <> void EventParameters::extractParameters(podio::Frame& frame) {
for(auto const& p: this->intParameters()) {
printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
@@ -161,7 +161,7 @@ namespace dd4hep {
#endif // DD4HEP_DDG4_Geant4Output2EDM4hep_DRC_H
//==========================================================================
-// AIDA Detector description implementation
+// AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
@@ -590,42 +590,13 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/,
}
}
//-------------------------------------------------------------------
- }
- else if( typeid( Geant4DRCalorimeter::Hit ) == coll->type().type() ){
- Geant4Sensitive* sd = coll->sensitive();
- int hit_creation_mode = sd->hitCreationMode();
+ } else if( typeid( Geant4DRCalorimeter::Hit ) == coll->type().type() ) {
// Create the hit container even if there are no entries!
- auto& hits = m_calorimeterHits[colName];
auto& DRhits = m_drcaloHits[colName];
auto& DRwaves = m_drcaloWaves[colName];
- for(unsigned i=0 ; i < nhits ; ++i){
- auto sch = hits.first->create();
+ for(unsigned i=0; i < nhits; ++i) {
const Geant4DRCalorimeter::Hit* hit = coll->hit(i);
- const auto& pos = hit->position;
- sch.setCellID( hit->cellID );
- sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
- sch.setEnergy( hit->energyDeposit/CLHEP::GeV );
-
- // now add the individual step contributions
- for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){
-
- auto sCaloHitCont = hits.second->create();
- sch.addToContributions( sCaloHitCont );
-
- const Geant4HitData::Contribution& c = *ci;
- int trackID = pm->particleID(c.trackID);
- auto mcp = m_particles.at(trackID);
- sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV );
- sCaloHitCont.setTime( c.time/CLHEP::ns );
- sCaloHitCont.setParticle( mcp );
-
- if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE ) {
- edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm);
- sCaloHitCont.setPDG( c.pdgID );
- sCaloHitCont.setStepPosition( p );
- }
- }
// For DRC raw calo hit & raw time series
auto rawCaloHits = DRhits.first->create();
@@ -636,7 +607,7 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/,
float timeStart = hit->GetTimeStart();
float timeEnd = hit->GetTimeEnd();
auto& timemap = hit->GetTimeStruct();
-
+
rawTimeStruct.setInterval(samplingT);
rawTimeStruct.setTime(timeStart);
rawTimeStruct.setCharge( static_cast(hit->GetPhotonCount()) );
@@ -691,4 +662,4 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/,
} else {
error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name());
}
-}
\ No newline at end of file
+}