Skip to content

Commit

Permalink
one giant commit to solve the rebase conflict
Browse files Browse the repository at this point in the history
SanghyunKo authored and atolosadelgado committed Jan 30, 2025
1 parent 211e0de commit 8bd7462
Showing 10 changed files with 833 additions and 586 deletions.
60 changes: 30 additions & 30 deletions FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml
Original file line number Diff line number Diff line change
@@ -16,15 +16,15 @@

<define>
<constant name="world_side" value="6100*mm"/> <!-- Used in LumiCal and HOMAbsorber, != world_size -->
<constant name="CrossingAngle" value="0.030*rad"/>
<constant name="CrossingAngle" value="0.030*rad"/>

<constant name="GlobalTrackerReadoutID" type="string" value="system:5,side:-2,layer:3,module:16,sensor:6"/>

<constant name="SolenoidField" value="2*tesla"/>

<!-- Detector IDs -->
<constant name="DetID_NOTUSED" value=" 0"/>

<constant name="DetID_VXD_Barrel" value=" 1"/>
<constant name="DetID_VXD_Disks" value=" 2"/>

@@ -40,10 +40,10 @@
<constant name="DetID_HCAL_Barrel" value=" 10"/>
<constant name="DetID_HCAL_Endcap" value=" 11"/>
<constant name="DetID_HCAL_Ring" value=" 12"/>

<constant name="DetID_Yoke_Barrel" value=" 13"/>
<constant name="DetID_Yoke_Endcap" value=" 14"/>

<constant name="DetID_LumiCal" value=" 15"/>
<constant name="DetID_LumiCalInstrumentation" value=" 16"/>
<constant name="DetID_LumiCalCooling" value=" 17"/>
@@ -68,7 +68,7 @@
<constant name="SeparatedBeamPipe_z" value="1190.0*mm"/> <!-- was 1159.97*mm in FCCDetectors/> -->
<constant name="CentralBeamPipe_rmax" value="10.0*mm"/>
<constant name="ConeBeamPipe_Rmax" value="28.9*mm" />
<constant name="BeamPipeGoldWidth" value="0.005*mm" />
<constant name="BeamPipeGoldWidth" value="0.005*mm" />
<constant name="BeamPipeGoldTolerance" value="0.001*mm" /> <!-- dummy tolerance, some small non zero value -->
<constant name="BeamPipeConeHalfAngle" value="(ConeBeamPipe_Rmax + BeamPipeWidthFirstCone - CentralBeamPipe_rmax ) / (SeparatedBeamPipe_z - CentralBeamPipe_zmax)" />

@@ -82,7 +82,7 @@
<constant name="MiddleOfSRMask_z" value="2.1*m" />
<constant name="SynchRadMaskSize" value="5*mm" /> <!-- mask tip is at 10 mm from the beamline -->
<constant name="mask_epsilon" value="0.001*mm" />

<!-- Vertex detector. Changing the values here it not enough to resize the detector, contact expert (Armin Ilg) -->
<constant name="VertexClearanceTheta" value="0.1035*rad"/> <!-- Clearance of vertex detector in mrad from IP. !!! Too high currently, lumical acceptance is 107 mrad so consider 110 mrad as limit -> Need to adapt vertex disks!!! -->
<constant name="VTXIB_r_min_clearance" value="1.0*mm"/> <!-- Clearance of vertex detector in radius, used for definiton of vertex DD4hep_SubdetectorAssembly -->
@@ -103,7 +103,7 @@
<constant name="VTXD_r_min" value="34.5*mm"/> <!-- Start of VTX disks in r -->
<constant name="VTXD_r_max" value="315.0*mm"/> <!-- Start of VTX disks in r -->

<constant name="VTX_r_min" value="VTXIB_r_min_layer-VTXIB_r_min_clearance"/>
<constant name="VTX_r_min" value="VTXIB_r_min_layer-VTXIB_r_min_clearance"/>
<constant name="VTX_r_max" value="VTXOB_r_max_layer+VTXOB_rmax_clearance"/>
<constant name="VTX_z_max" value="VTXD_z_max+4.5*cm"/>
<!-- End of VTX parameters -->
@@ -113,7 +113,7 @@
<constant name="DCH_inner_cyl_R_total" value=" 349.8 * mm" />
<constant name="DCH_outer_cyl_R_total" value=" 2015. * mm" />
<constant name="DCH_half_length_total" value=" 2250. * mm" />

<constant name="Solenoid_inner_radius" value="2100*mm"/>
<constant name="Solenoid_outer_radius" value="2400*mm"/>
<constant name="Solenoid_half_length" value="2380*mm"/>
@@ -122,7 +122,7 @@

<!-- Silicon wrapper. Changing the values of the disk parameters is not enough, please contact the expert (Armin Ilg) -->
<constant name="SiWrB_inner_radius" value="2040*mm"/>
<constant name="SiWrB_outer_radius" value="2080*mm"/>
<constant name="SiWrB_outer_radius" value="2080*mm"/>
<constant name="SiWrB_half_length" value="2400*mm"/>
<constant name="SiWrD_inner_radius" value="350.0*mm"/>
<constant name="SiWrD_outer_radius" value="2040.0*mm"/>
@@ -135,82 +135,82 @@
<constant name="EndPlateAbsorber_outer_radius" value="2030*mm"/>
<constant name="EndPlateAbsorber_z_min" value="2350*mm"/>
<constant name="EndPlateAbsorber_z_half_length" value="4.209/2.0*mm"/>

<constant name="YokeBarrel_inner_radius" value="4479*mm"/>
<constant name="YokeBarrel_outer_radius" value="6000*mm"/>
<constant name="YokeBarrel_half_length" value="3755*mm"/>
<constant name="YokeBarrel_symmetry" value="12"/>

<constant name="YokeEndcap_inner_radius" value="400*mm"/>
<constant name="YokeEndcap_outer_radius" value="6000*mm"/>
<constant name="YokeEndcap_min_z" value="3755*mm"/>
<constant name="YokeEndcap_max_z" value="5300*mm"/>
<constant name="YokeEndcap_outer_symmetry" value="12"/>
<constant name="YokeEndcap_inner_symmetry" value="0"/>

<constant name="CompSol_min_z" value="1230*mm"/>

<constant name="env_safety" value="0.1*mm"/>

<constant name="LumiCal_max_z" value="1186.5*mm" />
<constant name="LumiCal_min_z" value="1074*mm"/>

<constant name="LumiCal_dz" value="(LumiCal_max_z-LumiCal_min_z)/2.0"/>

<constant name="LumiCal_inner_radius" value="55.0*mm"/>
<constant name="LumiCal_outer_radius" value="112.0*mm- env_safety"/>

<constant name="LumiCal_Instr_thickness" value="20*mm"/>
<constant name="LumiCal_Instr_inner_radius" value="LumiCal_outer_radius"/>
<constant name="LumiCal_Instr_outer_radius" value="LumiCal_outer_radius+LumiCal_Instr_thickness - env_safety"/>

<constant name="LumiCal_Cool_thickness" value="9.75*mm"/>
<constant name="LumiCal_Cool_inner_radius" value="LumiCal_Instr_outer_radius"/>
<constant name="LumiCal_Cool_outer_radius" value="LumiCal_Instr_outer_radius+LumiCal_Cool_thickness"/>

<constant name="Lcal_services_rmax" value="LumiCal_outer_radius+30*mm"/>
<constant name="Lcal_offset_phi" value=" 0."/>
<!--preliminary LumiCal shielding-->
<!--back shielding-->
<constant name="LumiCal_Shield_inner_radius" value="LumiCal_inner_radius"/>
<constant name="LumiCal_Shield_outer_radius" value="LumiCal_outer_radius+LumiCal_Instr_thickness+LumiCal_Cool_thickness"/>
<constant name="LumiCal_shield_dz" value="1.75*mm"/>
<constant name="LumiCal_shield_dz" value="1.75*mm"/>
<!--nose-->
<constant name="LumiCal_NoseShield_inner_radius" value="LumiCal_inner_radius-5*mm"/>
<constant name="LumiCal_NoseShield_outer_radius" value="LumiCal_inner_radius+10*mm"/>
<constant name="LumiCal_nose_shield_dz" value="12*mm"/>

<constant name="BeamCal_inner_radius" value="32*mm"/>
<constant name="BeamCal_outer_radius" value="150*mm"/>
<constant name="BeamCal_min_z" value="3181*mm"/>
<constant name="BeamCal_max_z" value="3441*mm"/>
<constant name="BeamCal_dz" value="(BeamCal_max_z-BeamCal_min_z)/2.0"/>

<constant name="Kicker_inner_radius" value="4*mm"/>
<constant name="Kicker_outer_radius" value="25*mm"/>
<constant name="Kicker_min_z" value="3480*mm"/>
<constant name="Kicker_max_z" value="3780*mm"/>

<constant name="BPM_inner_radius" value="36*mm"/>
<constant name="BPM_outer_radius" value="55*mm"/>
<constant name="BPM_min_z" value="3790*mm"/>
<constant name="BPM_max_z" value="3880*mm"/>
<constant name="BPM_max_z" value="3880*mm"/>

<constant name="QD0_min_z" value="2000*mm"/>
<constant name="QD0_max_z" value="5200*mm"/>
<constant name="QD0Coil_outer_radius" value="30*mm"/>
<constant name="CollimatorInFrontOfQD0_dz" value="20*cm"/>
<constant name="CollimatorInFrontOfQD0_dz" value="20*cm"/>
<constant name="CollimatorInFrontOfQD0_radius" value="10*mm"/>
<constant name="CollimatorInFrontOfQD0_dr" value="16*mm"/>

<constant name="tracker_region_zmax" value="DCH_half_length_total"/>
<constant name="tracker_region_rmax" value="DCH_outer_cyl_R_total"/>

<!-- Pre-shower Parameters-->
<constant name = "psNumSides" value = "32"/> <!-- The number of sides of the pre-shower -->
<constant name = "psNumSides" value = "32"/> <!-- The number of sides of the pre-shower -->
<!-- Barrel -->
<constant name = "psBarrelFirstLayerRadius" value = "2420*mm"/> <!-- 1st Barrel microRWELL detector inner radius-> its the start point of thicknesses of the microRWELL material. In our case the shape is Polygon, so the radius is in the middle of the polygon side. -->
<constant name = "psBarrelLength" value = "4900*mm"/> <!--Barrel detector length, in the description of the detctor we always use the half-length -->
<constant name = "psBarrelLength" value = "4900*mm"/> <!--Barrel detector length, in the description of the detctor we always use the half-length -->
<!-- Endcap -->
<constant name = "psEndcapFirstLayerZOffset" value = "2400*mm"/> <!-- 1st Endcap microRWELL detector inner ZOffset-> its the start point of thicknesses of the microRWELL volume -->
<constant name = "psEndcapLayersInnerRadius" value = "390*mm"/> <!--Endcap detector inner radius, its the start point of thicknesses of the detector material ** it applies for both detector layers and yoke-->
@@ -232,18 +232,18 @@
<constant name="FiberDRCalo_endcap_nphi" value="144"/>

<!-- Muon System Parameters-->
<constant name = "numberOfSides" value = "8"/> <!-- The number of sides of the muon system e.g (Octagon, Hexagon, ...)-->
<constant name = "numberOfSides" value = "8"/> <!-- The number of sides of the muon system e.g (Octagon, Hexagon, ...)-->
<!-- Barrel -->
<constant name = "BarrelFirstLayerRadius" value = "4530*mm"/> <!-- 1st Barrel microRWELL detector inner radius-> its the start point of thicknesses of the microRWELL material. In our case the shape is Polygon, so the radius is in the middle of the polygon side. -->
<constant name = "BarrelLength" value = "9060*mm"/> <!--Barrel detector length, in the description of the detctor we always use the half-length -->
<constant name = "BarrelLength" value = "9060*mm"/> <!--Barrel detector length, in the description of the detctor we always use the half-length -->
<!-- Endcap -->
<constant name = "EndcapFirstLayerZOffset" value = "4530*mm"/> <!-- 1st Endcap microRWELL detector inner ZOffset-> its the start point of thicknesses of the microRWELL volume -->
<constant name = "EndcapLayersInnerRadius" value = "700*mm"/> <!--Endcap detector inner radius, its the start point of thicknesses of the detector material ** it applies for both detector layers and yoke-->
<constant name = "EndcapLayersOuterRadius" value = "5350*mm"/> <!--Endcap detector outer radius, its the end point of thicknesses of the detector material ** it applies for both detector layers and yoke-->
<!-- End of Muon system Parameters-->
</define>


<limits>
<limitset name="cal_limits">
<limit name="step_length_max" particles="*" value="5.0" unit="mm" />
Original file line number Diff line number Diff line change
@@ -13,7 +13,7 @@
The compact format for the dual-readout calorimeter (for FCCee IDEA)
</comment>
</info>

<!-- For optics -->
<properties>
<matrix name="RI_Air" coldim="2" values="
@@ -589,7 +589,7 @@
<sensitive type="DRcaloSiPMSD"/>
<sipmDim height="0.3*mm" material="PolyvinylChloride" vis="DRCGenericVis">
<sipmGlass material="DR_PyrexGlass" vis="DRCGlassVis"/>
<sipmWafer height="0.3*mm" material="Silicon" vis="DRCWaferVis" sensitive="true"/> <!-- Original height : 0.01 mm, change to 0.3 mm, same as sipmDim (to reduce memory)-->
<sipmWafer height="0.3*mm" material="Silicon" vis="DRCWaferVis" sensitive="false"/> <!-- Original height : 0.01 mm, change to 0.3 mm, same as sipmDim (to reduce memory)-->
</sipmDim>
<structure>
<dim distance="1.5*mm" dx="1.0*mm"/>
@@ -691,7 +691,7 @@
<readout name="DRcaloSiPMreadout">
<segmentation type="GridDRcalo_k4geo"/>
<!-- Mendatory to use the first 32 bits for tower infos & the last 32 bits for fiber/SiPM infos -->
<id>system:5,assembly:1,eta:-8,phi:9,x:32:-11,y:-9,c:1,module:2</id>
<id>system:5,assembly:1,eta:-8,phi:9,x:32:-12,y:-12,c:1,module:2</id>
</readout>
</readouts>

8 changes: 4 additions & 4 deletions FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
Original file line number Diff line number Diff line change
@@ -11,7 +11,7 @@
url="no"
status="development"
version="o1_v03">
<comment>
<comment>
Version o1_v03 of the IDEA detector
</comment>
</info>
@@ -37,9 +37,9 @@
<!-- shape based model of the beam pipe -->
<include ref="../../../MDI/compact/MDI_o1_v00/Beampipe_o4_v05.xml" />
<include ref="../../../MDI/compact/MDI_o1_v00/BeamInstrumentation_o1_v01.xml" />

<!-- engineered CAD model of the beam pipe -->
<!-- In order to use the CAD beampipe, build k4geo with the following CMake option:
<!-- In order to use the CAD beampipe, build k4geo with the following CMake option:
cmake -D INSTALL_BEAMPIPE_STL_FILES=ON which will download the files needed -->
<!-- <include ref="../../../MDI/compact/MDI_o1_v01/Beampipe_CADimport_o1_v02.xml" /> -->
<!-- <include ref="../../../MDI/compact/MDI_o1_v01/BeamInstrumentation_o1_v01.xml"/> -->
@@ -65,7 +65,7 @@

<!-- Import fiber-based dual-readout calorimeter -->
<!-- (uncomment the following line to effectively include it) -->
<!-- <include ref="FiberDualReadoutCalo_o1_v01.xml"/> -->
<include ref="FiberDualReadoutCalo_o1_v01.xml"/>

<!-- Import muon system -->
<include ref="MuonSystem_o1_v01.xml"/>
11 changes: 9 additions & 2 deletions detector/calorimeter/dual-readout/src/DRconstructor.cpp
Original file line number Diff line number Diff line change
@@ -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());
}
36 changes: 17 additions & 19 deletions example/SteeringFile_IDEA_o1_v03.py
Original file line number Diff line number Diff line change
@@ -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
299 changes: 5 additions & 294 deletions plugins/DRCaloFastSimModel.cpp
Original file line number Diff line number Diff line change
@@ -1,262 +1,17 @@
// Framework include files
#include <DDG4/Geant4FastSimShowerModel.inl.h>

#include <G4VFastSimulationModel.hh>
#include <G4ParticleDefinition.hh>
#include <G4ProcessManager.hh>
#include <G4OpProcessSubType.hh>
#include <G4GeometryTolerance.hh>
#include <G4Tubs.hh>
#include <G4OpBoundaryProcess.hh>
#include <G4OpAbsorption.hh>
#include <G4OpWLS.hh>
#include "G4FastStep.hh"

// C/C++ include files
#include "DRCaloFastSimModel.h"

/// 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<G4OpBoundaryProcess *>(theProcess);
else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpAbsorption)
pOpAbsorption = dynamic_cast<G4OpAbsorption *>(theProcess);
else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpWLS)
pOpWLS = dynamic_cast<G4OpWLS *>(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<DRCFiberModel>::initialize()
{
@@ -305,8 +60,7 @@ namespace dd4hep
}

template <>
bool Geant4FSShowerModel<DRCFiberModel>::check_trigger(const G4FastTrack &fasttrack)
{
bool Geant4FSShowerModel<DRCFiberModel>::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<G4Tubs *>(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<DRCFiberModel> Geant4DRCFiberModel;
362 changes: 321 additions & 41 deletions plugins/DRCaloFastSimModel.h
Original file line number Diff line number Diff line change
@@ -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 <G4ParticleDefinition.hh>
#include <G4ParticleTypes.hh>
#include <G4ProcessManager.hh>
#include <G4OpProcessSubType.hh>
#include <G4GeometryTolerance.hh>
#include <G4Tubs.hh>

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<G4Tubs*>(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<G4OpBoundaryProcess *>(theProcess);
else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpAbsorption)
pOpAbsorption = dynamic_cast<G4OpAbsorption *>(theProcess);
else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpWLS)
pOpWLS = dynamic_cast<G4OpWLS *>(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
593 changes: 436 additions & 157 deletions plugins/FiberDRCaloSDAction.cpp

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions plugins/FiberDRCaloSDAction.h
Original file line number Diff line number Diff line change
@@ -149,4 +149,5 @@ namespace dd4hep

}
}

#endif
43 changes: 7 additions & 36 deletions plugins/Geant4Output2EDM4hep_DRC.cpp
Original file line number Diff line number Diff line change
@@ -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<G4Event>& /*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<G4Event>& /*ctxt*/,
float timeStart = hit->GetTimeStart();
float timeEnd = hit->GetTimeEnd();
auto& timemap = hit->GetTimeStruct();

rawTimeStruct.setInterval(samplingT);
rawTimeStruct.setTime(timeStart);
rawTimeStruct.setCharge( static_cast<float>(hit->GetPhotonCount()) );
@@ -691,4 +662,4 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext<G4Event>& /*ctxt*/,
} else {
error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name());
}
}
}

0 comments on commit 8bd7462

Please sign in to comment.