diff --git a/FullSimLight/README.md b/FullSimLight/README.md
index ddfa402b3b430b03be8d515720db5bcf592097cf..a25f8b0e1119c72cc278151cd0f903c4b649833b 100644
--- a/FullSimLight/README.md
+++ b/FullSimLight/README.md
@@ -1,7 +1,7 @@
FullSimLight Project
==========================
-FullSimLight project consists of different tools based on [Geant4](https://geant4.web.cern.ch) toolkit, that can be run on multiple geometries:
+FullSimLight project consists of different tools based on [Geant4](https://geant4.web.cern.ch) toolkit, that can be run on multiple geometries:
- fullSimLight: a Geant4 based light simulation (geometry, transport in magnetic field and basic physics scoring)
- gmclash: a tool that runs clash detection on your input geometry, producing a json file report
@@ -15,14 +15,14 @@ The supported geometry formats are SQLite (.db), GDML (.gdml) and dual-use plugi
## Dependencies:
-FullSimLight repository has been recently merged in a big monorepository under the main [GeoModel repository](https://gitlab.cern.ch/GeoModelDev/GeoModel).
+FullSimLight repository has been recently merged in a big monorepository under the main [GeoModel repository](https://gitlab.cern.ch/GeoModelDev/GeoModel).
FullSimLight project depends on [GeoModelCore](https://gitlab.cern.ch/GeoModelDev/GeoModel/GeoModelCore), [GeoModelIO](https://gitlab.cern.ch/GeoModelDev/GeoModel/GeoModelIO),
-[Geant4](https://geant4.web.cern.ch),
-[GeoModelG4](https://gitlab.cern.ch/GeoModelDev/GeoModel/GeoModelG4) and
+[Geant4](https://geant4.web.cern.ch),
+[GeoModelG4](https://gitlab.cern.ch/GeoModelDev/GeoModel/GeoModelG4) and
[nlohmann_json](https://github.com/nlohmann/json.git)
-GeoModelCore, GeoModelIO and GeoModelG4 will be built as part of the monorepository, and also nhlomann_json can be built within the monorepository. Before installing them you need to have Geant4 installed on you machine.
-Following are the installation instructions.
+GeoModelCore, GeoModelIO and GeoModelG4 will be built as part of the monorepository, and also nhlomann_json can be built within the monorepository. Before installing them you need to have Geant4 installed on you machine.
+Following are the installation instructions.
## Geant4:
@@ -56,7 +56,7 @@ mkdir build ; cd build
cmake -DCMAKE_INSTALL_PREFIX=../../install -DCMAKE_BUILD_TYPE=Release ../ -DGEANT4_INSTALL_DATA=ON -DGEANT4_USE_GDML=ON -DGEANT4_BUILD_MULTITHREADED=ON
make -j8 ; make install
```
-Before running ./fullSimLight and ./gmgeantino make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
+Before running ./fullSimLight and ./gmgeantino make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
```bash
source <path_to_geant4_install_dir>/bin/geant4.sh
```
@@ -95,7 +95,7 @@ Clone the new GeoModel monorepository at [GeoModel repo](https://gitlab.cern.ch/
git clone https://gitlab.cern.ch/GeoModelDev/GeoModel.git
cd GeoModel
mkdir build ; cd build
-cmake -DCMAKE_INSTALL_PREFIX=../../install -DCMAKE_BUILD_TYPE=Release ../ -DGEOMODEL_BUILD_FULLSIMLIGHT=1 ../
+cmake -DCMAKE_INSTALL_PREFIX=../../install -DCMAKE_BUILD_TYPE=Release ../ -DGEOMODEL_BUILD_FULLSIMLIGHT=1 ../
make -j8 ; make install
```
@@ -137,7 +137,7 @@ The detector can be built starting from a SQLite .db file, from a GDML file or f
## ATLAS Geometry Files:
-The .gdml and .SQLite files of ATLAS geometry tags ATLAS-R2-2016-01-00-01 are available at:
+The .gdml and .SQLite files of ATLAS geometry tags ATLAS-R2-2016-01-00-01 are available at:
```bash
wget https://gitlab.cern.ch/GeoModelATLAS/geometry-data/raw/master/geometry/geometry-ATLAS-R2-2016-01-00-01.gdml
wget https://gitlab.cern.ch/GeoModelATLAS/geometry-data/raw/master/geometry/geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
@@ -154,9 +154,9 @@ wget https://gitlab.cern.ch/GeoModelATLAS/geometry-data/raw/master/geometry/geom
The application can be built and used both with sequential and multithreaded
Geant4 builds. In case of multithreaded Geant4 toolkit, the applications will
-run in proper multithreaded mode. You can find the executables under the build/bin directory and/or under the install/bin dir.
+run in proper multithreaded mode. You can find the executables under the build/bin directory and/or under the install/bin dir.
-NB: Before running fullSimLight make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
+NB: Before running fullSimLight make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
```bash
source <path_to_geant4_install_dir>/bin/geant4.sh
```
@@ -180,60 +180,56 @@ export G4ENSDFSTATEDATA=$G4INSTALL/data/G4ENSDFSTATE2.2
Run the executable with the --help option to see the available options:
``` bash
--m : REQUIRED : the standard Geant4 macro file name
--g : REQUIRED : the Geometry file name
+-m : REQUIRED : the standard Geant4 macro file name
+-g : REQUIRED : the Geometry file name
-o : flag ==> run the geometry overlap check (default: FALSE)
--f : physics list name (default: FTFP_BERT)
+-f : physics list name (default: FTFP_BERT)
-P : generate events with Pythia [config. available: ttbar/higgs/minbias or use ascii input file]
--p : flag ==> run the application in performance mode i.e. no user actions
- : - ==> run the application in NON performance mode i.e. with user actions (default)
-```
+-p : flag ==> run the application in performance mode i.e. no user actions
+ : - ==> run the application in NON performance mode i.e. with user actions (default)
+```
FullSimLight uses by default the Geant4 particle gun as primary generator, but it supports also
input events from the Pythia generator (see the Primary generator section for more details)
-A minimal set of "observable" is collected during the simulation per-primary
-particle type: mean energy deposit, mean charged and neutral step lengths,
-mean number of steps made by charged and neutral particles, mean number of
-secondary e-, e+ and gamma particles. The result is reported at the end of
-each event for each primary particle that were transported in the given event.
-At the end of the simulation a final report is printed showing the run time,
-the primary generator and magnetic field settings used during the run, the
-total number of events and primary particles transported and the per-primary
-type simulation statistics of the above-mentioned quantities.
+A minimal set of "observable" is collected during the simulation: the mean energy
+deposit, mean charged and neutral track lengths, mean number of steps made by charged
+and neutral particles, mean number of secondary e-, e+ and gamma particles.
+These are reported at the end of the simulation together with the run time and
+total number of events transported.
The simulation can be executed in "performance" mode by providing the -p
input flag. No any user actions are created in this case beyond the
only one RunAction (only for the Master-thread in case of MT) that will
start and stop a timer at the beginning and the end of the simulation
-(initialization time won't be included). Therefore, there is no scoring
+(initialization time won't be included). Therefore, there is no any scoring
in this case.
## Examples
-During the installation a default macro file <macro.g4> will be installed in your share/FullSimLight directory.
+During the installation a default macro file <macro.g4> will be installed in your `share/FullSimLight` directory.
To execute the application using the <macro.g4> macro file, with the FTFP_BERT_ATL
Physics List, in performance mode and building the detector from the geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db file :
``` bash
./fullSimLight -m ../share/FullSimLight/macro.g4 -f FTFP_BERT_ATL -p -g geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
-```
+```
To execute the application using the <macro.g4> macro file, with the FTFP_BERT_ATL
Physics List, not in performance mode and building the detector from my geometry gdml file:
``` bash
./fullSimLight -m ../share/FullSimLight/macro.g4 -f FTFP_BERT_ATL -g mygeometry.gdml
-```
+```
To execute the application using the <macro.g4> macro file and building the detector with a geometry described in one of the [GeoModelPlugins repo](https://gitlab.cern.ch/atlas/GeoModelPlugins), i.e. *HGTDPlugin* :
``` bash
./fullSimLight -m ../share/FullSimLight/macro.g4 -g libHGTDPlugin.1.0.0.dylib
-```
+```
## Parameters settings via geant4 macro
-Fullsimlight and in general Geant4 based simulations, need a Geant4 macro to read some input parameters. The default macro used by fullSimLight is called 'macro.g4' and it should be found under the <install-path>/share/FullSimLight directory. The macro can be edited to change some parameters, i.e the verbosity, the number of threads, or to tune the simulation. The most relevant macro commands are explained in what follows.
+Fullsimlight and in general Geant4 based simulations, need a Geant4 macro to read some input parameters when executed in batch mode. The default macro used by `fullSimLight` is called 'macro.g4' and it should be found under the `<install-path>/share/FullSimLight` directory. The macro can be edited to change some parameters, i.e the verbosity, the number of threads, or to tune the simulation. The most relevant macro commands are explained in what follows.
## Magnetic field
-
+
A constant magnetic field can be set through the macro command:
-
+
``` bash
/mydet/setField <field-value> <unit>
```
@@ -257,26 +253,26 @@ Run the application with the --help flag to see the options:
-f : OPTIONAL: magnetic field filename [.data/.root] (default : use ATLAS magnetic field maps)
-r : FLAG: use root field map (default : false, use .data file)
--s : FLAG: set Solenoid Off
--t : FLAG: set Toroids Off
+-s : FLAG: set Solenoid Off
+-t : FLAG: set Toroids Off
```
By default the file tested is the bmagatlas_09_fullAsym20400.data. Use the -r option to test the ROOT file maps, add the -s flag to set the solenoid off and use the toroid_bfieldmap_0_20400_14m_version5.root:
``` bash
-./testMagneticField -r -s
+./testMagneticField -r -s
```
-Use the -t to set the Toroids off, and test the solenoid_bfieldmap_7730_0_14m_version5.root file.
+Use the -t to set the Toroids off, and test the solenoid_bfieldmap_7730_0_14m_version5.root file.
``` bash
-./testMagneticField -r -t
+./testMagneticField -r -t
```
## Primary Generator
The primary generator used by default is the Geant4 particle gun, but FullSimLight also supports the [Pythia generator](http://home.thep.lu.se/Pythia/)
-
+
## Particle gun
The particle gun used by default will generate primary particles
at the (0,0,0) position with the following options:
@@ -300,7 +296,7 @@ The primary particle energy can be set through the macro command:
By default, i.e. if it is not specified by the above command, the kinetic energy will be randomly selected for each individual primary particle from the [1 GeV, 100 GeV] uniformly.
### Primary particle direction:
-The primary particle momentum direction can be set through the macro command:
+The primary particle momentum direction can be set through the macro command:
``` bash
/mygun/direction <xdir-value> <ydir-value> <zdir-value>
@@ -312,34 +308,36 @@ The primary particle type can be set through the macro command:
``` bash
/mygun/particle <particle-name>
-```
+```
By default, i.e. if it is not specified by the above command, the type will be randomly selected from a pre-defined list for each individual primary particle uniformly. The current list of particles includes e-, e+ and gamma particles. It can be extended by adding more particles to the list in the MyPrimaryGeneratorAction class.
## Pythia generator
-
+
FullSimLight supports Pythia as primary particles generator. In order to use Pythia, the user should have it installed in their system and if Pythia is found FullSImLight will be compiled with the support on. There are three different default options available when using the -P or --pythia flag (i.e. ttbar, higgs and minbias):
``` bash
-P : generate events with Pythia [config. available: ttbar/higgs/minbias or use ascii input file]
- ```
- Alternatively the user can plug their own Pythia configuration file to simulate the desired events.
+ ```
+ Alternatively the user can plug their own Pythia configuration file to simulate the desired events.
For example, in order to simulate the default *ttbar* events, the command to be run is the following:
-
+
``` bash
-./fullSimLight -m ../share/FullSimLight/pythia.g4 -P ttbar -g geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
- ```
- The number of events that the user wants to simulate must be specified in the g4 macro file. A specific *pythia.g4* macro file can be found in the *share* directory, that should be used when simulating Pythia events and can be edited according to the user needs.
+./fullSimLight -m ../share/FullSimLight/pythia.g4 -P ttbar -g geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
+ ```
+ The number of events that the user wants to simulate must be specified in the g4 macro file. A specific *pythia.g4* macro file can be found in the *share* directory, that should be used when simulating Pythia events and can be edited according to the user needs.
## Physics List
- The Physics List can be specified as an input argument with the -f flag
- (e.g. -f FTFP_BERT). Notice that the name of the Geant4 built in Physics List
+ The Physics List can be specified as an input argument with the `-f` flag
+ (e.g. `-f FTFP_BERT`). Notice that the name of the Geant4 built in Physics List
must be in upper case, exactly as the corresponding header file. By default,
- i.e. if the Physics List name is not provided as an input argument, the FTFP_BERT
+ i.e. if the Physics List name is not provided as an input argument, the `FTFP_BERT`
Physics List will be used. FullSimLight has also the possibility to use a
- special custom Physics List MyGVPhysicsList (a Geant4 Physics
- List including only those particles, processes and models that are available
- and used in the corresponding GeantV application with exactly the same
- settings). In order to use it one need to specify -f GV as input parameter.
+ special custom Physics List `FTFP_BERT_ATL_WDCK`: a Geant4 Physics
+ List that corresponds to the `FTFP_BERT_ATL` reference physics list with the differences:
+ - the standard EM physics constructor is replaced by the local `StandardEmWithWoodcock` that can utilise *Woodcock*-tracking of gamma photons in a specified detector region
+ - the EM extra physics constructor is replaced with a local version that is identical to the one delivered by Geant4 with the only difference that its `G4GammaGeneralProcess` parts are replaced with the local `GammaGeneralProcess`
+ - the neutron tracking cut module is replaced by one that sets a 150 [ns] tracking cut
+ This latter is applied now in all cases of dedicated ATLAS physics lists (i.e. those with the `ATL` sub-string in their name) in order to use a similar setting that is in Athena.
# GeoModelClash: run and options
@@ -348,9 +346,9 @@ GeoModelClash (gmclash) allows to run geometry overlap checks on a geometry file
Run the executable with the --help option to see the available options:
``` bash
--g : MANDATORY: the Geometry file name [.db/.gdml/.dylib/.so]
+-g : MANDATORY: the Geometry file name [.db/.gdml/.dylib/.so]
-o : OPTIONAL : clashes report file name (default: gmclash_report.json)
-```
+```
The output json file format is the following:
@@ -366,7 +364,7 @@ The output json file format is the following:
"x": -1.736718203796568,
"y": -1263.348806272393,
"z": -166.75403155804725
-```
+```
where:
* *distance* is the minimum estimated distance of the overlap
* *typeOfClash* can be 0 for *withMother*, 1 for *withSister*, 2 for *fullyEncapsSister*, 3 for *invalidSolid*
@@ -377,22 +375,22 @@ where:
To run GeoModelClash one has to specify with the -g flag the geometry file (this is mandatory). By default gmclash writes out the clashes report in the *gmclash_report.json* file
``` bash
./gmclash -g geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
-```
+```
To execute a clash detection on a geometry described with the SQLite file *LArBarrel.db* and write out the clashes report in the *cr_LArBarrel.json* file :
``` bash
-./gmclash -g LArBarrel.db -o cr_LArBarrel.json
-```
+./gmclash -g LArBarrel.db -o cr_LArBarrel.json
+```
To execute a clash detection on a geometry described with one of the [GeoModelPlugins repo](https://gitlab.cern.ch/atlas/GeoModelPlugins), i.e. *HGTDPlugin* and write out the clashes report in the *cr_HGTD.json* file :
``` bash
-./gmclash -g libHGTDPlugin.1.0.0.dylib -o cr_HGTD.json
-```
+./gmclash -g libHGTDPlugin.1.0.0.dylib -o cr_HGTD.json
+```
# GeoModelGeantino: run and options
GeoModelGeantino (gmgeantino) is a Geant4 based application that allows you to produce geantino maps for the geometry specified as input. It supports .db/.gdml/.dylib/.so geometry formats and it writes out the geantino maps in a ROOT file. However, it does not depend on ROOT, cause it uses the G4AnalysisManager to create/fill/write 1D and 2D Profiles.
The 1D and 2D profiles are filled during the simulation, per Event or per Step. The creation of different profiles can be tuned with command line flags. In general X-Y, Z-R, eta-Phi RadiationLength/InteractionLength profiles can be created per DetectorVolume/Material/Element.
-gmgeantino uses a default geant4 macro, called *geantino.g4*, to take some input parameters. You should find this macro under your <install-path>/share/FullSimLight directory. You can edit this macro to change some parameters, like the verbosity, the n. of threads, the n. of primaries per event, the primary particle energy.. etc.
+gmgeantino uses a default geant4 macro, called *geantino.g4*, to take some input parameters. You should find this macro under your <install-path>/share/FullSimLight directory. You can edit this macro to change some parameters, like the verbosity, the n. of threads, the n. of primaries per event, the primary particle energy.. etc.
By default the primary particles shot by 'gmgeantino' are geantinos (this parameter should not be changed). By default the number of simulated geantinos is set to 1000. To increase the resolution of the maps, the n. of simulated events can be increased, editing the macro and changing the value at the following line:
@@ -400,7 +398,7 @@ By default the primary particles shot by 'gmgeantino' are geantinos (this parame
/run/beamOn <n. of events>
```
-NB: Before running gmgeantino make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
+NB: Before running gmgeantino make sure to source the *geant4.sh* file to set correctly all the Geant4 environment variables.
```bash
source <path_to_geant4_install_dir>/bin/geant4.sh
```
@@ -426,18 +424,18 @@ Run the executable with the --help option to see the available options:
``` bash
./gmgeantino --help
--g : [REQUIRED] the Geometry file name (supported extensions: .db/.gdml/.dylib/.so)
--m : [OPTIONAL] the standard Geant4 macro file name (default: 'geantino.g4')
--r : [OPTIONAL] r limit for geantino maps in mm (default: '12500')
--z : [OPTIONAL] z limit for geantino maps in mm (default: '23000')
--x : [OPTIONAL] x limit for geantino maps in mm (default: '12500')
--y : [OPTIONAL] y limit for geantino maps in mm (default: '12500')
--o : [OPTIONAL] output ROOT file name (supported extention: .root - default: 'geantinoMaps.root')
+-g : [REQUIRED] the Geometry file name (supported extensions: .db/.gdml/.dylib/.so)
+-m : [OPTIONAL] the standard Geant4 macro file name (default: 'geantino.g4')
+-r : [OPTIONAL] r limit for geantino maps in mm (default: '12500')
+-z : [OPTIONAL] z limit for geantino maps in mm (default: '23000')
+-x : [OPTIONAL] x limit for geantino maps in mm (default: '12500')
+-y : [OPTIONAL] y limit for geantino maps in mm (default: '12500')
+-o : [OPTIONAL] output ROOT file name (supported extention: .root - default: 'geantinoMaps.root')
-e : [FLAG] use this flag to create eta-phi radiation-interaction length 1D profile histograms (caveat: the process might run out of memory!)
-d : [FLAG] use this flag to create xy-rz radiation-interaction length 2D profile histograms for 'detectors' (caveat: the process might run out of memory!)
-a : [FLAG] use this flag to create xy-rz radiation-interaction length 2D profile histograms for 'materials' (caveat: the process might run out of memory!)
-l : [FLAG] use this flag to create xy-rz radiation-interaction length 2D profile histograms for 'elements' (caveat: the process might run out of memory!)
-```
+```
## Examples
@@ -445,13 +443,13 @@ Run the executable with the --help option to see the available options:
To run *gmgeantino* one has to specify with the -g flag the input geometry file (mandatory option). By default a 'geantinoMaps.root' file is created as output and it contains R-Z RadiationLenght and Interaction Lenght 2D profiles. To run gmgeantino on *LArBarrel.db* geometry, with the default *geantino.g4* macro file, and producing eta-phi maps and detector maps:
``` bash
-./gmgeantino -m ../share/FullSimLight/geantino.g4 -g LArBarrel.db -e -d
-```
+./gmgeantino -m ../share/FullSimLight/geantino.g4 -g LArBarrel.db -e -d
+```
To produce geantino maps of a geometry described by one of the [GeoModelPlugins repo](https://gitlab.cern.ch/atlas/GeoModelPlugins), i.e. *HGTDPlugin*, using a custom macro file *mymacro.g4*, activate detectors/materials and elements maps, and write out the geantino maps in the *geantinoMaps_HGTD.root* file :
``` bash
-./gmgeantino -m mymacro.g4 -g libHGTDPlugin.1.0.0.dylib -d -a -l -o geantinoMaps_HGTD.root
-```
+./gmgeantino -m mymacro.g4 -g libHGTDPlugin.1.0.0.dylib -d -a -l -o geantinoMaps_HGTD.root
+```
# GeoModelMassCalculator: run and options
@@ -462,12 +460,12 @@ GeoModelMassCalculator (gmmasscalc) is a command line tool that calculates the i
Run the executable with the --help option to see the available options:
``` bash
--g : [MANDATORY] the Geometry file name [.db/.gdml/.dylib/.so]
--p : [OPTIONAL] prefix of the Logical Volumes of interest (i.e. Pixel::)
--m : [OPTIONAL] material of interest (i.e. Aluminium)
+-g : [MANDATORY] the Geometry file name [.db/.gdml/.dylib/.so]
+-p : [OPTIONAL] prefix of the Logical Volumes of interest (i.e. Pixel::)
+-m : [OPTIONAL] material of interest (i.e. Aluminium)
-o : [OPTIONAL] mass report json file name (default: gmmasscalc_report.json)
-v : [OPTIONAL] verbose mode: 1 - print all the Logical Volume names and materials, 2 - print materials composition (default: off)
-```
+```
By default (if the optional flag are not used) *gmmasscalc*, takes the main *World Volume*, and calculates the inclusive and exclusive masses of the respective daughters, saving the calculated quantities in the output json file. At the end of the report, the total masses are reported for the whole *World Volume*. The output json file format is the following:
``` bash
@@ -478,7 +476,7 @@ By default (if the optional flag are not used) *gmmasscalc*, takes the main *Wor
"physicalVolumeName": "SCT",
"volumeCopyNo": 16969,
"volumeEntityType": "G4Tubs"
-```
+```
where:
- *exclusiveMass* is the mass of the considered volume only (from which the volumes occupied by the daughters volumes have been subtracted)
- *inclusiveMass* is the mass of the considered volume, comprehensive of the masses of the respective daughters (propagated in an iterative way to their daughter volumes).
@@ -495,19 +493,19 @@ At the end of the report, additional information about the whole geometry is rep
"logicalVolumeName": "newWorldLog",
"material": "Air",
"volumeEntityType": "World Volume"
-```
+```
where:
- *apparent weight in Air* by definition, the weight of a body as affected by the buoyancy of a fluid (such as air) in which it is immersed. It is calculated only on the total geometry, assuming that the World Volume is made of Air
- *exclusiveFilteredMass* is the sum of the exclusive masses of all the volumes with density>densityThreshold (0.02 g/cm3)
-*excludedFilteredMass* is the sum of the exclusive masses of all the volumes with density<densityThreshold (0.02 g/cm3)
-In addition to the default behaviour, *gmmasscalc* offers the possibility to apply 2 filters to the geometry, described in what follows.
+In addition to the default behaviour, *gmmasscalc* offers the possibility to apply 2 filters to the geometry, described in what follows.
The -p (--prefix) option allows to indicate the prefix of the volumes of interest. In this case *gmmasscalc* will loop over the geometry tree and calculate the mass of every Logical Volume that has the specified prefix in its name. For every Logical volume found, a different entry will be filled in the output report file.
-The -m (--material) option allows to specify to which material the user is interested. In this case *gmmasscalc* will loop over the geometry tree and calculate the mass of those volumes that are made of the material of interest.
+The -m (--material) option allows to specify to which material the user is interested. In this case *gmmasscalc* will loop over the geometry tree and calculate the mass of those volumes that are made of the material of interest.
-If both the -p and -m flags are used, *gmmasscalc* will combine the 2 filters and retrieve in the output file only the masses of the volumes containing the specified prefix and made of the desired material.
+If both the -p and -m flags are used, *gmmasscalc* will combine the 2 filters and retrieve in the output file only the masses of the volumes containing the specified prefix and made of the desired material.
At the end of the report, the total exclusive mass for the requested geometry is reported for all the volumes that satisfy the user request. The last item in the report will look like the following:
@@ -524,29 +522,29 @@ At the end of the report, the total exclusive mass for the requested geometry is
```
where:
-- *logicalVolumeName* and *material* are respectively the prefix used as filter for the Logical Volume names in the geometry tree navigation, and the material used as a second filter in the search.
+- *logicalVolumeName* and *material* are respectively the prefix used as filter for the Logical Volume names in the geometry tree navigation, and the material used as a second filter in the search.
-The -v flag is a verbosity flag. Two different levels of verbosity can be set: level 1 prints all the logical volume names and respective materials, level 2 adds also the material composition for each material in the geometry.
+The -v flag is a verbosity flag. Two different levels of verbosity can be set: level 1 prints all the logical volume names and respective materials, level 2 adds also the material composition for each material in the geometry.
## Examples
To run GeoModelMassCalculator one has to specify with the -g flag the geometry file (this is mandatory). By default *gmmasscalc* writes out the masses report in the *gmmasscalc_report.json* file:
``` bash
./gmmasscalc -g geometry-ATLAS-R2-2016-01-00-01_wSPECIALSHAPE.db
-```
+```
To calculate the mass of a geometry described with the SQLite file *LArBarrel.db* and write out the masses report in the *mass_LArBarrel.json* file :
``` bash
-./gmmasscalc -g LArBarrel.db -o mass_LArBarrel.json
-```
+./gmmasscalc -g LArBarrel.db -o mass_LArBarrel.json
+```
To calculate the mass of a geometry described with the SQLite file *SCT.db*, only for the logical volumes that have *BRLSensor* in their names, and only for those volumes that are made of *Silicon* and write out the masses report in the *mass_SCT.json* file:
``` bash
-./gmmasscalc -g SCT.db -p BRLSensor -m Silicon -o mass_SCT.json
-```
+./gmmasscalc -g SCT.db -p BRLSensor -m Silicon -o mass_SCT.json
+```
To calculate the mass of a geometry described with one of the [GeoModelPlugins repo](https://gitlab.cern.ch/atlas/GeoModelPlugins), i.e. *HGTDPlugin*, looking for volumes made of *Aluminium* and write out the masses report in the *mass_HGTD.json* file :
``` bash
-./gmmasscalc -g libHGTDPlugin.1.0.0.dylib -m Aluminium -o mass_HGTD.json
-```
+./gmmasscalc -g libHGTDPlugin.1.0.0.dylib -m Aluminium -o mass_HGTD.json
+```
# GeoModel2GDML: run and options
@@ -559,36 +557,36 @@ Run the executable with the --help option to see the available options:
``` bash
./gm2gdml --help
- **** Parameters:
--g : [MANDATORY] input geometry file name [.db/.gdml/.dylib/.so]
+ **** Parameters:
+-g : [MANDATORY] input geometry file name [.db/.gdml/.dylib/.so]
-o : [OPTIONAL] output GDML geometry file name (default: geometry.gdml)
-```
+```
## Examples
-To run GeoModel2GDML one has to specify with the -g flag the geometry file (this is mandatory). By default a 'geometry.gdml' output file is created, but the name of the file can be changed using the -o option.
+To run GeoModel2GDML one has to specify with the -g flag the geometry file (this is mandatory). By default a 'geometry.gdml' output file is created, but the name of the file can be changed using the -o option.
To run gm2gdml on *LArBarrel.db* geometry, and dump the output gdml file in the *LArBarrel.gdml* file:
``` bash
./gm2gdml -g LArBarrel.db -o LArBarrel.gdml
-```
+```
To dump your geometry described with a plugin, i.e. *HGTDPlugin*, and write out the geometry file in the *HGTDPlugin.gdml* file :
``` bash
./gm2gdml -g libHGTDPlugin.1.0.0.dylib -o HGTDPlugin.gdml
-```
+```
# Final Remark
To get the most performant Geant4 build, one needs to build the Geant4 toolkit
with the following CMake options (in addition to -DGEANT4_USE_GDML=ON):
-
+
``` bash
-DGEANT4_BUILD_MULTITHREADED=ON
-DGEANT4_BUILD_VERBOSE_CODE=OFF
-DGEANT4_BUILD_STORE_TRAJECTORY=OFF
-```
+```
Then set the number of threads in the macro.g4 macro to the maximum available
and execute the fullSimLight application with the -p flag.
diff --git a/FullSimLight/fullSimLight.cc b/FullSimLight/fullSimLight.cc
index 04dde09d648ef1075973edc2316e2fea967212e9..7ebd744ce6d88403aaa94e82c234df83f982a086 100644
--- a/FullSimLight/fullSimLight.cc
+++ b/FullSimLight/fullSimLight.cc
@@ -20,7 +20,9 @@
#include "Randomize.hh"
#include "MyDetectorConstruction.hh"
-#include "MyGVPhysicsList.hh"
+#include "StandardEmWithWoodcock.hh"
+#include "EmExtraPhysics.hh"
+#include "G4NeutronTrackingCut.hh"
#include "MyActionInitialization.hh"
#include "PythiaPrimaryGeneratorAction.hh"
@@ -78,30 +80,58 @@ int main(int argc, char** argv) {
G4RunManager* runManager = new G4RunManager;
#endif
- // set mandatory initialization classes
- // 1. Detector construction
- MyDetectorConstruction* detector = new MyDetectorConstruction;
-
- if (parRunOverlapCheck) detector->SetRunOverlapCheck(true);
-
- detector->SetGeometryFileName (geometryFileName);
- runManager->SetUserInitialization(detector);
-
- // 2. Physics list
+ // 1. Physics list
+ G4bool activateRegions = false;
+ G4VModularPhysicsList* physList = nullptr;
G4PhysListFactory factory;
if (factory.IsReferencePhysList(parPhysListName)) {
- G4VModularPhysicsList* physList = factory.GetReferencePhysList(parPhysListName);
- runManager->SetUserInitialization(physList);
- } else if (parPhysListName==G4String("GV")) {
- G4VUserPhysicsList* physList = new MyGVPhysicsList();
- runManager->SetUserInitialization(physList);
+ physList = factory.GetReferencePhysList(parPhysListName);
+ } else if (parPhysListName==G4String("FTFP_BERT_ATL_WDCK")) {
+ G4cout << "<<< Geant4 FTFP_BERT_ATL physics list with the local Woodcock settings " << G4endl;
+ physList = factory.GetReferencePhysList("FTFP_BERT_ATL");
+ // the local em-standard physics with Woodcock tracking for gamma
+ StandardEmWithWoodcock* em0AndWDCK = new StandardEmWithWoodcock;
+ // set the region name and low energy limit for Woodcock tracking
+ em0AndWDCK->SetRegionNameForWoodcockTracking("EMEC");
+ em0AndWDCK->SetLowEnergyLimitForWoodcockTracking(200.0*CLHEP::keV);
+ physList->ReplacePhysics(em0AndWDCK);
+ // the local version of the `G4EmExtraPhysics` that will use the local `GammaGeneralProcess`
+ G4VPhysicsConstructor* emExtra = new EmExtraPhysics;
+ physList->ReplacePhysics(emExtra);
+ //physList->RemovePhysics("G4GammaLeptoNuclearPhys");
+ // make sure that regions will also be added to the detector
+ activateRegions = true;
} else {
G4cerr << "ERROR: Physics List " << parPhysListName << " UNKNOWN!" << G4endl;
return -1;
}
+ // In cases of ATLAS physics lists, set the neutron tracking cut to be 150 [ns] as in Athena
+ if (parPhysListName.find("ATL") != std::string::npos) {
+ G4NeutronTrackingCut* neutronCut = new G4NeutronTrackingCut("neutronCutphysics", 1);
+ neutronCut->SetTimeLimit(150.0*CLHEP::ns);
+ physList->ReplacePhysics(neutronCut);
+ }
+
+ // register the final version of the physics list in the run manager
+ runManager->SetUserInitialization(physList);
+
+ // 2. Detector construction
+ MyDetectorConstruction* detector = new MyDetectorConstruction;
+
+ if (parRunOverlapCheck) detector->SetRunOverlapCheck(true);
+ if (activateRegions) detector->SetAddRegions(true);
+
+ detector->SetGeometryFileName (geometryFileName);
+ runManager->SetUserInitialization(detector);
// 3. User action
- runManager->SetUserInitialization(new MyActionInitialization(parIsPerformance));
+ MyActionInitialization* actInit = new MyActionInitialization(parIsPerformance);
+ // set the name of a region in which we are interested to see a very basic simulation
+ // stat e.g. "EMEC" (NOTE: only if the given region can be found and executed in
+ // non-perfomance mode)
+ const G4String nameSpecialScoringRegion = ""; //"EMEC"
+ actInit->SetSpecialScoringRegionName(nameSpecialScoringRegion);
+ runManager->SetUserInitialization(actInit);
// 4. Run the simulation in batch mode
G4UImanager* UI = G4UImanager::GetUIpointer();
diff --git a/FullSimLight/geantinoMaps.cc b/FullSimLight/geantinoMaps.cc
index 1f741ef99ffa359876b58ddd35c8c5a511be4f82..892cf4a92ea15fd5d6b4f7f677036226be82989e 100644
--- a/FullSimLight/geantinoMaps.cc
+++ b/FullSimLight/geantinoMaps.cc
@@ -20,7 +20,6 @@
#include "Randomize.hh"
#include "MyDetectorConstruction.hh"
-#include "MyGVPhysicsList.hh"
#include "MyActionInitialization.hh"
@@ -51,10 +50,10 @@ void Help();
int main(int argc, char** argv) {
-
+
// Get input arguments
GetInputArguments(argc, argv);
-
+
G4cout
<< " =============== Running geantinoMaps ================ " << G4endl
<< " Geant4 macro = " << parMacroFileName << G4endl
@@ -65,7 +64,7 @@ int main(int argc, char** argv) {
<< " Create Materials maps = " << parCreateMaterialsMaps << G4endl
<< " Create Elements maps = " << parCreateElementsMaps << G4endl
<< " ===================================================== " << G4endl;
-
+
// JFB: Check that the macro file exists and is readable:
if (access(parMacroFileName,R_OK)) {
G4cout << G4endl;
@@ -96,7 +95,7 @@ int main(int argc, char** argv) {
G4cout << "INFO: Exiting" <<G4endl;
return 1;
#endif
-
+
//choose the Random engine: set to MixMax explicitely (default form 10.4)
G4Random::setTheEngine(new CLHEP::MixMaxRng);
// set seed and print info
@@ -107,7 +106,7 @@ int main(int argc, char** argv) {
<< " Initial seed = " << G4Random::getTheSeed() << G4endl
<< " ===================================================== " << G4endl
<< G4endl;
-
+
// Construct the default run manager
#ifdef G4MULTITHREADED
G4MTRunManager* runManager = new G4MTRunManager;
@@ -117,29 +116,26 @@ int main(int argc, char** argv) {
#else
G4RunManager* runManager = new G4RunManager;
#endif
-
+
// set mandatory initialization classes
// 1. Detector construction
MyDetectorConstruction* detector = new MyDetectorConstruction;
-
+
if (parRunOverlapCheck) detector->SetRunOverlapCheck(true);
-
+
detector->SetGeometryFileName (parGeometryFileName);
runManager->SetUserInitialization(detector);
-
+
// 2. Physics list
G4PhysListFactory factory;
if (factory.IsReferencePhysList(parPhysListName)) {
G4VModularPhysicsList* physList = factory.GetReferencePhysList(parPhysListName);
runManager->SetUserInitialization(physList);
- } else if (parPhysListName==G4String("GV")) {
- G4VUserPhysicsList* physList = new MyGVPhysicsList();
- runManager->SetUserInitialization(physList);
} else {
G4cerr << "ERROR: Physics List " << parPhysListName << " UNKNOWN!" << G4endl;
return -1;
}
-
+
// 3. User action
MyActionInitialization* myAct = new MyActionInitialization(parIsPerformance, parCreateGeantinoMaps,parOutputFileName);
myAct->SetRlimit(parRlimit);
@@ -153,12 +149,12 @@ int main(int argc, char** argv) {
runManager->SetUserInitialization(myAct);
-
+
// 4. Run the simulation in batch mode
G4UImanager* UI = G4UImanager::GetUIpointer();
G4String command = "/control/execute ";
UI->ApplyCommand(command+parMacroFileName);
-
+
// Print out the final random number
G4cout << G4endl
<< " ================================================================= " << G4endl
diff --git a/FullSimLight/include/EmExtraPhysics.hh b/FullSimLight/include/EmExtraPhysics.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6dec17f5c45c07e39ad81eeb9ec3e0e654075325
--- /dev/null
+++ b/FullSimLight/include/EmExtraPhysics.hh
@@ -0,0 +1,292 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+//---------------------------------------------------------------------------
+//
+// ClassName: EmExtraPhysics
+//
+// Author: 2002 J.P. Wellisch
+//
+// Modified:
+// 10.11.2005 V.Ivanchenko edit to provide a standard
+// 19.06.2006 V.Ivanchenko add mu-nuclear process
+// 16.10.2012 A.Ribon: renamed G4EmExtraBertiniPhysics as EmExtraPhysics
+// 31.01.2018 V. Grichine: add neutrino-electron process and xsc
+//
+// NOTE: M. Novak 30.11.2021.
+// This is local version of the `G4EmExtraPhysics` constructor renamed to
+// `EmExtraPhysics` and all `G4GammaGeneralProcess` casts has been chaned
+// to `GammaGeneralProcess` in order ot be able to use with the local
+// `GammaGeneralProcess` when Woodcock tracking is utilised. This should
+// be removed (together with the `GammaGeneralProcess`) when moving to
+// Geant4.11.0 version of the `G4GammaGeneralProcess` as base class for
+// `WoodcockProcess`.
+//----------------------------------------------------------------------------
+//
+
+#ifndef EmExtraPhysics_h
+#define EmExtraPhysics_h 1
+
+#include "G4Version.hh"
+#include "G4VPhysicsConstructor.hh"
+#include "globals.hh"
+
+
+#if G4VERSION_NUMBER<1070
+
+ class G4PhotoNuclearProcess;
+ class G4CascadeInterface;
+
+ class EmExtraPhysics : public G4VPhysicsConstructor
+ {
+ public:
+
+ EmExtraPhysics(G4int ver = 1);
+
+ // obsolete
+ EmExtraPhysics(const G4String& name);
+
+ virtual ~EmExtraPhysics();
+
+ void ConstructParticle();
+ void ConstructProcess();
+
+ void Synch(G4bool val);
+ void SynchAll(G4bool val);
+ void GammaNuclear(G4bool val);
+ void LENDGammaNuclear(G4bool val);
+ void ElectroNuclear(G4bool val);
+ void MuonNuclear(G4bool val);
+ void GammaToMuMu(G4bool val);
+ void PositronToMuMu(G4bool val);
+ void PositronToHadrons(G4bool val);
+ void GammaToMuMuFactor(G4double val);
+ void PositronToMuMuFactor(G4double val);
+ void PositronToHadronsFactor(G4double val);
+
+ void NeutrinoActivated(G4bool val);
+ void NuETotXscActivated(G4bool val);
+ void SetNuEleCcBias(G4double bf);
+ void SetNuEleNcBias(G4double bf);
+ void SetNuNucleusBias(G4double bf);
+ void SetNuDetectorName(const G4String& dn);
+
+ private:
+
+ void ConstructGammaElectroNuclear();
+
+ void ConstructLENDGammaNuclear(G4CascadeInterface* cascade,
+ G4PhotoNuclearProcess* gnuc);
+
+ G4bool gnActivated;
+ G4bool eActivated;
+ G4bool gLENDActivated;
+ G4bool munActivated;
+ G4bool synActivated;
+ G4bool synActivatedForAll;
+ G4bool gmumuActivated;
+ G4bool pmumuActivated;
+ G4bool phadActivated;
+ G4bool fNuActivated;
+ G4bool fNuETotXscActivated;
+
+ G4double gmumuFactor;
+ G4double pmumuFactor;
+ G4double phadFactor;
+ G4double fNuEleCcBias;
+ G4double fNuEleNcBias;
+ G4double fNuNucleusBias;
+
+ G4String fNuDetectorName;
+
+ G4int verbose;
+ };
+
+#elif G4VERSION_NUMBER<1100
+
+ #include "G4VPhysicsConstructor.hh"
+ #include "globals.hh"
+
+ #include "G4EmMessenger.hh"
+
+ class G4CascadeInterface;
+ class G4PhotoNuclearProcess;
+
+ class EmExtraPhysics : public G4VPhysicsConstructor
+ {
+ public:
+
+ EmExtraPhysics(G4int ver = 1);
+
+ // obsolete
+ EmExtraPhysics(const G4String& name);
+
+ virtual ~EmExtraPhysics();
+
+ void ConstructParticle();
+ void ConstructProcess();
+
+ void Synch(G4bool val);
+ void SynchAll(G4bool val);
+ void GammaNuclear(G4bool val);
+ void LENDGammaNuclear(G4bool val);
+ void ElectroNuclear(G4bool val);
+ void MuonNuclear(G4bool val);
+ void GammaToMuMu(G4bool val);
+ void PositronToMuMu(G4bool val);
+ void PositronToHadrons(G4bool val);
+ void GammaToMuMuFactor(G4double val);
+ void PositronToMuMuFactor(G4double val);
+ void PositronToHadronsFactor(G4double val);
+ void GammaNuclearLEModelLimit(G4double val);
+
+ void NeutrinoActivated(G4bool val);
+ void NuETotXscActivated(G4bool val);
+ void SetUseGammaNuclearXS(G4bool val);
+ void SetNuEleCcBias(G4double bf);
+ void SetNuEleNcBias(G4double bf);
+ void SetNuNucleusBias(G4double bf);
+ void SetNuDetectorName(const G4String& dn);
+
+ private:
+
+ void ConstructGammaElectroNuclear();
+
+ void ConstructLENDGammaNuclear(G4CascadeInterface* cascade,
+ G4PhotoNuclearProcess* gnuc);
+
+ G4bool gnActivated;
+ G4bool eActivated;
+ G4bool gLENDActivated;
+ G4bool munActivated;
+ G4bool synActivated;
+ G4bool synActivatedForAll;
+ G4bool gmumuActivated;
+ G4bool pmumuActivated;
+ G4bool phadActivated;
+ G4bool fNuActivated;
+ G4bool fNuETotXscActivated;
+ G4bool fUseGammaNuclearXS;
+
+ G4double gmumuFactor;
+ G4double pmumuFactor;
+ G4double phadFactor;
+ G4double fNuEleCcBias;
+ G4double fNuEleNcBias;
+ G4double fNuNucleusBias;
+ G4double fGNLowEnergyLimit;
+
+ G4String fNuDetectorName;
+
+ G4EmMessenger* theMessenger;
+ G4int verbose;
+ };
+
+#else
+
+ #include "G4VPhysicsConstructor.hh"
+ #include "globals.hh"
+
+ #include "G4EmMessenger.hh"
+
+ class G4CascadeInterface;
+ class G4HadronInelasticProcess;
+
+ class EmExtraPhysics : public G4VPhysicsConstructor
+ {
+ public:
+
+ EmExtraPhysics(G4int ver = 1);
+
+ // obsolete
+ EmExtraPhysics(const G4String& name);
+
+ virtual ~EmExtraPhysics();
+
+ void ConstructParticle();
+ void ConstructProcess();
+
+ void Synch(G4bool val);
+ void SynchAll(G4bool val);
+ void GammaNuclear(G4bool val);
+ void LENDGammaNuclear(G4bool val);
+ void ElectroNuclear(G4bool val);
+ void MuonNuclear(G4bool val);
+ void GammaToMuMu(G4bool val);
+ void PositronToMuMu(G4bool val);
+ void PositronToHadrons(G4bool val);
+ void GammaToMuMuFactor(G4double val);
+ void PositronToMuMuFactor(G4double val);
+ void PositronToHadronsFactor(G4double val);
+ void GammaNuclearLEModelLimit(G4double val);
+
+ void NeutrinoActivated(G4bool val);
+ void NuETotXscActivated(G4bool val);
+ void SetUseGammaNuclearXS(G4bool val);
+ void SetNuEleCcBias(G4double bf);
+ void SetNuEleNcBias(G4double bf);
+ void SetNuNucleusBias(G4double bf);
+ void SetNuDetectorName(const G4String& dn);
+
+ private:
+
+ void ConstructGammaElectroNuclear();
+
+ void ConstructLENDGammaNuclear(G4CascadeInterface* cascade,
+ G4HadronInelasticProcess* gnuc);
+
+ G4bool gnActivated;
+ G4bool eActivated;
+ G4bool gLENDActivated;
+ G4bool munActivated;
+ G4bool synActivated;
+ G4bool synActivatedForAll;
+ G4bool gmumuActivated;
+ G4bool pmumuActivated;
+ G4bool phadActivated;
+ G4bool fNuActivated;
+ G4bool fNuETotXscActivated;
+ G4bool fUseGammaNuclearXS;
+
+ G4double gmumuFactor;
+ G4double pmumuFactor;
+ G4double phadFactor;
+ G4double fNuEleCcBias;
+ G4double fNuEleNcBias;
+ G4double fNuNucleusBias;
+ G4double fGNLowEnergyLimit;
+
+ G4String fNuDetectorName;
+
+ G4EmMessenger* theMessenger;
+ G4int verbose;
+ };
+
+
+#endif // G4VERSION_NUMBER
+
+
+#endif
diff --git a/FullSimLight/include/GammaGeneralProcess.hh b/FullSimLight/include/GammaGeneralProcess.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0b62030e951010c77e6ad18034f651588bef4eff
--- /dev/null
+++ b/FullSimLight/include/GammaGeneralProcess.hh
@@ -0,0 +1,242 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class header file
+//
+//
+// File name: GammaGeneralProcess
+//
+// Author: Vladimir Ivanchenko
+//
+// Creation date: 19.07.2018
+//
+// Modifications:
+//
+// Class Description:
+//
+// NOTE: M. Novak: 27.11.2021
+//
+// It is a local version of the `GammaGeneralProcess` that will be available
+// in Geant4-11.0. More exactly, a version that can be utilised on the top of
+// Geant4-11.0-beta as a base class for the `WoodcockProcess` (also local). When
+// using Geant4-11.0 or later, it is recommended to remove this local version
+// and derive the `WoodcockProcess` directly from the `GammaGeneralProcess`
+// delivered by the toolkit.
+
+// -------------------------------------------------------------------
+//
+
+#ifndef GammaGeneralProcess_h
+#define GammaGeneralProcess_h 1
+
+#include "G4VEmProcess.hh"
+#include "globals.hh"
+#include "G4EmDataHandler.hh"
+
+class G4Step;
+class G4Track;
+class G4ParticleDefinition;
+class G4VParticleChange;
+class G4GammaConversionToMuons;
+class G4HadronicProcess;
+class G4MaterialCutsCouple;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+class GammaGeneralProcess : public G4VEmProcess
+{
+public:
+
+ explicit GammaGeneralProcess(const G4String& pname="local-GammaGeneralProc");
+
+ ~GammaGeneralProcess() override;
+
+ G4bool IsApplicable(const G4ParticleDefinition&) override;
+
+ void AddEmProcess(G4VEmProcess*);
+
+ void AddMMProcess(G4GammaConversionToMuons*);
+
+ void AddHadProcess(G4HadronicProcess*);
+
+ void ProcessDescription(std::ostream& outFile) const override;
+
+protected:
+
+ void InitialiseProcess(const G4ParticleDefinition*) override;
+
+public:
+
+ // Initialise for build of tables
+ void PreparePhysicsTable(const G4ParticleDefinition&) override;
+
+ // Build physics table during initialisation
+ void BuildPhysicsTable(const G4ParticleDefinition&) override;
+
+ // Called before tracking of each new G4Track
+ void StartTracking(G4Track*) override;
+
+ // implementation of virtual method, specific for GammaGeneralProcess
+ G4double PostStepGetPhysicalInteractionLength(
+ const G4Track& track,
+ G4double previousStepSize,
+ G4ForceCondition* condition) override;
+
+ // implementation of virtual method, specific for GammaGeneralProcess
+ G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&) override;
+
+ // Store PhysicsTable in a file.
+ // Return false in case of failure at I/O
+ G4bool StorePhysicsTable(const G4ParticleDefinition*,
+ const G4String& directory,
+ G4bool ascii = false) override;
+
+ // Retrieve Physics from a file.
+ // (return true if the Physics Table can be build by using file)
+ // (return false if the process has no functionality or in case of failure)
+ // File name should is constructed as processName+particleName and the
+ // should be placed under the directory specifed by the argument.
+ G4bool RetrievePhysicsTable(const G4ParticleDefinition*,
+ const G4String& directory,
+ G4bool ascii) override;
+
+ // Temporary method
+ const G4String& GetSubProcessName() const;
+
+ // Temporary method
+ G4int GetSubProcessSubType() const;
+
+ G4VEmProcess* GetEmProcess(const G4String& name) override;
+
+ inline const G4VProcess* GetSelectedProcess() const;
+
+ // hide copy constructor and assignment operator
+ GammaGeneralProcess(GammaGeneralProcess &) = delete;
+ GammaGeneralProcess & operator=
+ (const GammaGeneralProcess &right) = delete;
+
+protected:
+
+ G4double GetMeanFreePath(const G4Track& track, G4double previousStepSize,
+ G4ForceCondition* condition) override;
+
+ inline G4double ComputeGeneralLambda(size_t idxe, size_t idxt);
+
+ inline G4double GetProbability(size_t idxt);
+
+ inline void SelectedProcess(const G4Step& step, G4VProcess* ptr);
+
+ inline void SelectEmProcess(const G4Step&, G4VEmProcess*);
+
+ void SelectHadProcess(const G4Track&, const G4Step&, G4HadronicProcess*);
+
+ // It returns the cross section per volume for energy/ material
+ G4double TotalCrossSectionPerVolume();
+
+private:
+
+ G4bool RetrieveTable(G4VEmProcess*, const G4String& directory,
+ G4bool ascii);
+
+protected:
+
+ G4HadronicProcess* theGammaNuclear = nullptr;
+ G4VProcess* selectedProc = nullptr;
+
+ G4double preStepLogE = 1.0;
+ G4double factor = 1.0;
+
+
+private:
+ static G4EmDataHandler* theHandler;
+ static const size_t nTables = 15;
+ static G4bool theT[nTables];
+ static G4String nameT[nTables];
+
+ G4VEmProcess* thePhotoElectric = nullptr;
+ G4VEmProcess* theCompton = nullptr;
+ G4VEmProcess* theConversionEE = nullptr;
+ G4VEmProcess* theRayleigh = nullptr;
+ G4GammaConversionToMuons* theConversionMM = nullptr;
+
+ G4double minPEEnergy;
+ G4double minEEEnergy;
+ G4double minMMEnergy;
+ G4double peLambda = 0.0;
+
+ size_t nLowE = 40;
+ size_t nHighE = 50;
+ size_t idxEnergy = 0;
+ G4bool splineFlag = false;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline G4double
+GammaGeneralProcess::ComputeGeneralLambda(size_t idxe, size_t idxt)
+{
+ idxEnergy = idxe;
+ return factor*theHandler->GetVector(idxt, basedCoupleIndex)
+ ->LogVectorValue(preStepKinEnergy, preStepLogE);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline G4double GammaGeneralProcess::GetProbability(size_t idxt)
+{
+ return theHandler->GetVector(idxt, basedCoupleIndex)
+ ->LogVectorValue(preStepKinEnergy, preStepLogE);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline void
+GammaGeneralProcess::SelectedProcess(const G4Step& step, G4VProcess* ptr)
+{
+ selectedProc = ptr;
+ step.GetPostStepPoint()->SetProcessDefinedStep(ptr);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline void GammaGeneralProcess::SelectEmProcess(const G4Step& step, G4VEmProcess* proc)
+{
+ proc->CurrentSetup(currentCouple, preStepKinEnergy);
+ SelectedProcess(step, proc);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline const G4VProcess* GammaGeneralProcess::GetSelectedProcess() const
+{
+ return selectedProc;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#endif
diff --git a/FullSimLight/include/MyActionInitialization.hh b/FullSimLight/include/MyActionInitialization.hh
index a4eaab0910c2db3787711a2b19474c7baf8874c7..bba4df966b1180a60e409bbe11cc31a021ca6de3 100644
--- a/FullSimLight/include/MyActionInitialization.hh
+++ b/FullSimLight/include/MyActionInitialization.hh
@@ -26,8 +26,11 @@ public:
void SetXlimit(G4double x) { fXlimit = x; }
void SetYlimit(G4double y) { fYlimit = y; }
+ void SetSpecialScoringRegionName(const G4String& rname) { fSpecialScoringRegionName = rname; }
+
+
private:
-
+
bool fIsPerformance;
bool fCreateGeantinoMaps;
bool fCreateEtaPhiMaps;
@@ -35,6 +38,7 @@ private:
bool fCreateMaterialsMaps;
bool fCreateElementsMaps;
G4String fGeantinoMapsFilename;
+ G4String fSpecialScoringRegionName;
G4double fRlimit;
G4double fZlimit;
G4double fXlimit;
diff --git a/FullSimLight/include/MyEventAction.hh b/FullSimLight/include/MyEventAction.hh
index f426fa8efe3f8a4454add5774dff9696ad9a1217..dea3d2f38cc91f16a483b7af804c80fbc4e33a7a 100644
--- a/FullSimLight/include/MyEventAction.hh
+++ b/FullSimLight/include/MyEventAction.hh
@@ -19,12 +19,17 @@ public:
void BeginOfEventAction(const G4Event* evt) override;
void EndOfEventAction(const G4Event* evt) override;
- void AddData(G4double edep, G4double length, G4bool ischarged);
- void AddSecondaryTrack(const G4Track* track);
+ void AddData(G4double edep, G4double length, G4bool ischarged, G4bool isspecial);
+ void AddSecondaryTrack(const G4Track* track, G4bool isspecial);
+
+ void SetIsSpecialScoringRegion(G4bool val) { fIsSpecialScoring = val; }
private:
+ G4bool fIsSpecialScoring;
+
MyEventData fEventData;
+ MyEventData fEventDataSpecialRegion;
};
diff --git a/FullSimLight/include/MyGVPhysicsList.hh b/FullSimLight/include/MyGVPhysicsList.hh
deleted file mode 100644
index 6e672ead0a67fbc14a6ea62a2be4340f44c87b11..0000000000000000000000000000000000000000
--- a/FullSimLight/include/MyGVPhysicsList.hh
+++ /dev/null
@@ -1,23 +0,0 @@
-
-#ifndef MyGVPhysicsList_h
-#define MyGVPhysicsList_h 1
-
-#include "G4VUserPhysicsList.hh"
-
-class MyGVPhysicsList: public G4VUserPhysicsList {
-public:
-
- MyGVPhysicsList();
- ~MyGVPhysicsList();
-
- virtual void ConstructParticle();
- virtual void ConstructProcess();
-// virtual void SetCuts();
-
-private:
-
- void BuildEMPhysics();
-
-};
-
-#endif
diff --git a/FullSimLight/include/MyRun.hh b/FullSimLight/include/MyRun.hh
index 0b7c2dbc086b6f0869bf6e2ecea79d0f4304ba32..cc7332c642351a824a5e1930d21b5cce1e1581e2 100644
--- a/FullSimLight/include/MyRun.hh
+++ b/FullSimLight/include/MyRun.hh
@@ -7,6 +7,7 @@
#include "MyRunData.hh"
class G4Track;
+class G4Region;
class MyEventData;
class MyRun : public G4Run {
@@ -18,15 +19,25 @@ public:
void Merge(const G4Run*) override;
- void FillPerEvent(const MyEventData&);
+ void SetSpecialScoringRegion(G4Region* reg) { fScoringRegion = reg; }
+
+ void FillPerEvent(const MyEventData&, G4bool);
void EndOfRun();
- const MyRunData& GetRunData() const { return fRunData; }
+ const MyRunData& GetRunData(G4bool isspecial) const { return isspecial ? fRunDataSpecialRegion : fRunData; }
private:
+ void AddData(const MyEventData& data, MyRunData& runData);
+ void PrintEndOfRunStat(MyRunData& runData, G4double norm);
+
+private:
+
+ G4Region* fScoringRegion;
+
MyRunData fRunData;
+ MyRunData fRunDataSpecialRegion;
};
diff --git a/FullSimLight/include/MyRunAction.hh b/FullSimLight/include/MyRunAction.hh
index dd35c06a3686138fcb850fbf57dadd2e7e6d506f..228c4778c54ed10abf75f9b50fbc98df5034ccea 100644
--- a/FullSimLight/include/MyRunAction.hh
+++ b/FullSimLight/include/MyRunAction.hh
@@ -10,6 +10,8 @@
class MyRun;
class G4Timer;
+class MySteppingAction;
+class MyTrackingAction;
class MyRunAction: public G4UserRunAction {
@@ -25,14 +27,23 @@ public:
void SetPerformanceFlag(G4bool val) { fIsPerformance = val; }
void SetPythiaConfig(G4String config) { fPythiaConfig = config; }
+
+ void SetSteppingAction(MySteppingAction* stpact) { fSteppingAction = stpact; }
+ void SetTrackingAction(MyTrackingAction* tract) { fTrackingAction = tract; }
+
+ void SetSpecialScoringRegionName(const G4String& rname) { fSpecialScoringRegionName = rname; }
+
private:
- G4bool fIsPerformance;
- G4bool fIsGeantino;
- MyRun* fRun;
- G4Timer* fTimer;
- G4String fGeantinoMapsFilename;
- G4String fPythiaConfig;
+ G4bool fIsPerformance;
+ G4bool fIsGeantino;
+ MyRun* fRun;
+ G4Timer* fTimer;
+ MySteppingAction* fSteppingAction;
+ MyTrackingAction* fTrackingAction;
+ G4String fGeantinoMapsFilename;
+ G4String fPythiaConfig;
+ G4String fSpecialScoringRegionName;
//TO DO: make private and add Get methods
diff --git a/FullSimLight/include/MySteppingAction.hh b/FullSimLight/include/MySteppingAction.hh
index 0af275c24e19c78f92dd2a405133e45afc37ed27..21bc46017539c01979608e78589e4f13754511a4 100644
--- a/FullSimLight/include/MySteppingAction.hh
+++ b/FullSimLight/include/MySteppingAction.hh
@@ -7,6 +7,7 @@
class MyEventAction;
class G4Step;
+class G4Region;
class MySteppingAction : public G4UserSteppingAction {
@@ -17,9 +18,12 @@ public:
void UserSteppingAction(const G4Step*) override;
+ void SetScoringRegion(G4Region* reg);
+
private:
MyEventAction* fEventAction;
+ G4Region* fScoringRegion;
};
diff --git a/FullSimLight/include/MyTrackingAction.hh b/FullSimLight/include/MyTrackingAction.hh
index 262d5558269989113ac60e444aeec8c3f906c92c..f471395d3e21e3e235029b7f3539ee27eb748af9 100644
--- a/FullSimLight/include/MyTrackingAction.hh
+++ b/FullSimLight/include/MyTrackingAction.hh
@@ -7,6 +7,7 @@
class G4Track;
class MyEventAction;
+class G4Region;
class MyTrackingAction : public G4UserTrackingAction {
@@ -15,11 +16,14 @@ public:
MyTrackingAction(MyEventAction*);
~MyTrackingAction() override;
+ void SetScoringRegion(G4Region* reg);
+
void PreUserTrackingAction(const G4Track*) override;
private:
- MyEventAction *fEventAction;
+ MyEventAction* fEventAction;
+ G4Region* fScoringRegion;
};
diff --git a/FullSimLight/include/StandardEmWithWoodcock.hh b/FullSimLight/include/StandardEmWithWoodcock.hh
new file mode 100644
index 0000000000000000000000000000000000000000..587fb53b8d02ec46066a583fa33c2209d28a6d9e
--- /dev/null
+++ b/FullSimLight/include/StandardEmWithWoodcock.hh
@@ -0,0 +1,95 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+//
+//---------------------------------------------------------------------------
+//
+// ClassName: StandardEmWithWoodcock
+//
+// Author: V.Ivanchenko 09.11.2005
+//
+// Modified:
+// 05.12.2005 V.Ivanchenko add controlled verbosity
+// 23.11.2006 V.Ivanchenko remove mscStepLimit option and improve cout
+// 17.11.2021 M. Novak: added Woodcock tracking in a given detector region
+//
+//
+// NOTE: M. Novak. 27.11.2021
+//
+// This EM physics constructor is almost identical to the EM standard (Em-opt0)
+// physics constructor with the possibility of utilising Woodcock-tracking in a
+// given detector region. The name of the region needs to be provided to the
+// `WoodcockProcess` at construction. The `WoodckockProcess` falls into the base
+// Gamma-General process (a local version of the `G4GammaGeneralProcess`) when:
+// - the given detector region cannot be found in the detector
+// - the gamma energy is below a given threshold (by default 0; can be set at
+// construction time)
+// As in the case of the base Gamma-General process, all the processes neeeds to
+// be assigned to the `G4Gamma` particle through the `WoodckockProcess` (and not
+// directly!).
+//
+//----------------------------------------------------------------------------
+//
+// This class provides construction of default EM standard physics
+//
+
+#ifndef StandardEmWithWoodcock_h
+#define StandardEmWithWoodcock_h 1
+
+#include "G4VPhysicsConstructor.hh"
+#include "G4EmParticleList.hh"
+#include "globals.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class StandardEmWithWoodcock : public G4VPhysicsConstructor
+{
+public:
+
+ explicit StandardEmWithWoodcock(G4int ver=1, const G4String& name="StandardEmWithWoodcock");
+
+ ~StandardEmWithWoodcock() override;
+
+ void ConstructParticle() override;
+ void ConstructProcess() override;
+
+
+ void SetRegionNameForWoodcockTracking(const G4String& rname) { fWDCKRegionName = rname; }
+ void SetLowEnergyLimitForWoodcockTracking(G4double val) { fWDCKLowEnergyThreshold = val; }
+
+private:
+
+ G4int verbose;
+ G4EmParticleList partList;
+
+ G4String fWDCKRegionName;
+ G4double fWDCKLowEnergyThreshold;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
diff --git a/FullSimLight/include/WoodcockProcess.hh b/FullSimLight/include/WoodcockProcess.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8bab51467ddee96ad6657510846a6dc759fdcd4b
--- /dev/null
+++ b/FullSimLight/include/WoodcockProcess.hh
@@ -0,0 +1,120 @@
+
+#ifndef WoodcockProcess_h
+#define WoodcockProcess_h 1
+
+#include "GammaGeneralProcess.hh"
+
+#include "G4AffineTransform.hh"
+
+#include <map>
+#include <vector>
+
+class G4Region;
+class G4VSolid;
+class G4Navigator;
+class G4MaterialCutsCouple;
+
+/**
+ * @file WoodcockProcess.hh
+ * @author Mihaly Novak
+ * @date 29.11.2021
+ *
+ * Geant4 process for utilising Woodcock tracking in a given detector region.
+ *
+ * It is assumed that Woodcock tracking is required in a single detector
+ * region and the process requires the name of this region. Then Woodcock
+ * tracking will be used whenever the Gamma particle is inside this detector
+ * region and its energy is above the low energy threshold (that is also
+ * an input argument of the process).
+ *
+ * Requires: (at least) Geant4-11.0 that contains an updated version of the base
+ * G4GammaGeneralProcess. In order to be able to utilise even before
+ * Geant4-11.0 is avaibale, this version of the `WoodcockProcess` is
+ * dervived from a local version of that `G4GammaGeneralProcess` with
+ * the name of `GammaGeneralProcess`. This `GammaGeneralProcess` is
+ * (almost) identical to that will be available in Geant4-11.0.
+ */
+
+class WoodcockProcess : public GammaGeneralProcess {
+
+public:
+
+ WoodcockProcess(const G4String& wdckRegionName, G4double wdckLowELimit=0.0);
+ ~WoodcockProcess() override;
+
+
+ void StartTracking(G4Track*) override;
+
+ void PreparePhysicsTable(const G4ParticleDefinition&) override;
+
+ void BuildPhysicsTable(const G4ParticleDefinition&) override;
+
+ G4double PostStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4ForceCondition* condition) override;
+
+ G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&) override;
+
+
+
+private:
+
+ // This will be called only once when the very first gamma makes the very first
+ // step inside the WDCK region and all the WDCK related initialisations, that
+ // requires a closed geometry (and a global point inside the WDCK region), will
+ // be done.
+ void InitWoodcockTracking(const G4ThreeVector& aPointInsideTheWDCKregion, const G4ThreeVector& aDirectionInto);
+
+ // After the above init, this method can be used to compute the Distance to the
+ // boundary of the (single) root volume of the WDCK region from the given (global)
+ // point along the given direction.
+ G4double ComputeStepToWDCKRegionBoundary(const G4ThreeVector& pGlobalpoint, const G4ThreeVector& pDirection);
+
+ // Computes the total (i.e. sum of individual processes) macroscopic cross section
+ // for the given material and photon energy (computed only if any of those has changed)
+ void ComputeTotalMacrCrossSection(const G4MaterialCutsCouple* coupe, const G4double energy, const G4double logenergy);
+
+ void FindWDCKRootLogicalVolume(const G4Track& track);
+
+private:
+
+ // WDCK tracking data for a single root logical volume of the WDCK region
+ struct WDCKDataPerRootLVolume {
+ G4AffineTransform fAffineTransform;
+ G4LogicalVolume* fLogicalVolume;
+ G4Navigator* fNavigator;
+ G4VSolid* fSolid;
+ G4bool fNeedsInit;
+ };
+
+ // name of the region where WDCK tracking is required: given at construction
+ // NOTE: if region is not found with the given name then no WDCK tracking is used
+ const G4String fWDCKRegionName;
+ // low energy limit below which normal tracking will be used even inside the WDCK region
+ G4double fWDCKLowEnergyThreshold;
+
+ // number of root logical volumes of the WDCK region
+ std::size_t fWDCKNumRootLVolume;
+ // a `map` of the WDCK root logical volume IDs to continous indices
+ std::map<G4int, G4int> fWDCKRootLVolIDToIndxMap;
+
+ // pointer to the WDCK region (single region is assumed)
+ G4Region* fWDCKRegion;
+ // the heaviest material-cuts of the (single) WDCK region
+ G4MaterialCutsCouple* fWDCKCouple;
+
+ // flag to indicate that a WDCK region step was limited by the WDCK region boundary
+ G4bool fOnWDCKRegionBoundary;
+ // flag to indicate that a step was performed inside the WDCK region
+ G4bool fInWDCKRegion;
+ // the current root logical volume index: under which the current tracking
+ // happens inside the WDCK region
+ G4int fWDCKRootLVolumeIndx;
+ // length of the step done in the WDCK region
+ G4double fWDCKStepLength;
+
+
+ // data for each of the root logical volume of the WDCK region
+ std::vector<WDCKDataPerRootLVolume> fWDCKDataPerRootLVolume;
+
+};
+
+#endif // WoodcockProcess_h
diff --git a/FullSimLight/src/EmExtraPhysics.cc b/FullSimLight/src/EmExtraPhysics.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a01daadd563e28e2b5b00a86c868327af47c1a65
--- /dev/null
+++ b/FullSimLight/src/EmExtraPhysics.cc
@@ -0,0 +1,1434 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+//---------------------------------------------------------------------------
+//
+// ClassName: EmExtraPhysics
+//
+// Author: 2002 J.P. Wellisch
+//
+// Modified:
+//
+// 10.11.2005 V.Ivanchenko edit to provide a standard
+// 19.06.2006 V.Ivanchenko add mu-nuclear process
+// 16.10.2012 A.Ribon: renamed G4EmExtraBertiniPhysics as EmExtraPhysics
+// 10.04.2014 A.Dotti: Add MT functionality for messenger
+// 24.04.2014 A.Ribon: switched on muon-nuclear by default
+// 29.01.2018 V.Grichine, adding neutrinos
+// 07.05.2019 V.Grichine, adding muon neutrino nucleus interactions
+//
+//
+///////////////////////////////////////////////////////////////
+
+#include "EmExtraPhysics.hh"
+
+#if G4VERSION_NUMBER<1070
+
+#include "G4SystemOfUnits.hh"
+
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+#include "G4AntiNeutrinoE.hh"
+#include "G4NeutrinoE.hh"
+#include "G4AntiNeutrinoMu.hh"
+#include "G4NeutrinoMu.hh"
+#include "G4AntiNeutrinoTau.hh"
+#include "G4NeutrinoTau.hh"
+
+#include "G4SynchrotronRadiation.hh"
+#include "G4MuonNuclearProcess.hh"
+#include "G4MuonVDNuclearModel.hh"
+#include "G4ElectroVDNuclearModel.hh"
+#include "G4TheoFSGenerator.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4QGSModel.hh"
+#include "G4GammaParticipants.hh"
+#include "G4QGSMFragmentation.hh"
+#include "G4ExcitedStringDecay.hh"
+#include "G4CascadeInterface.hh"
+
+#include "G4LENDorBERTModel.hh"
+#include "G4LENDCombinedCrossSection.hh"
+
+#include "G4GammaConversionToMuons.hh"
+#include "G4AnnihiToMuPair.hh"
+#include "G4eeToHadrons.hh"
+
+#include "G4PhotoNuclearProcess.hh"
+#include "G4ElectronNuclearProcess.hh"
+#include "G4PositronNuclearProcess.hh"
+
+#include "G4NeutrinoElectronProcess.hh"
+#include "G4NeutrinoElectronTotXsc.hh"
+#include "G4NeutrinoElectronCcModel.hh"
+#include "G4NeutrinoElectronNcModel.hh"
+#include "G4MuNeutrinoNucleusProcess.hh"
+#include "G4MuNeutrinoNucleusTotXsc.hh"
+#include "G4NuMuNucleusCcModel.hh"
+#include "G4NuMuNucleusNcModel.hh"
+#include "G4LossTableManager.hh"
+
+//#include "G4GammaGeneralProcess.hh"
+#include "GammaGeneralProcess.hh"
+
+#include "G4HadronicParameters.hh"
+#include "G4PhysicsListHelper.hh"
+#include "G4BuilderType.hh"
+
+
+//////////////////////////////////////
+
+EmExtraPhysics::EmExtraPhysics(G4int ver):
+ G4VPhysicsConstructor("GammaLeptoNuclearPhys"),
+ gnActivated (true),
+ eActivated (true),
+ gLENDActivated(false),
+ munActivated(true),
+ synActivated(false),
+ synActivatedForAll(false),
+ gmumuActivated(false),
+ pmumuActivated(false),
+ phadActivated (false),
+ fNuActivated (false),
+ fNuETotXscActivated (false),
+ gmumuFactor (1.0),
+ pmumuFactor (1.0),
+ phadFactor (1.0),
+ fNuEleCcBias(1.0),
+ fNuEleNcBias(1.0),
+ fNuNucleusBias(1.0),
+ fNuDetectorName("0"),
+ verbose(ver)
+{
+ SetPhysicsType(bEmExtra);
+ if(verbose > 1) G4cout << "### EmExtraPhysics" << G4endl;
+}
+
+EmExtraPhysics::EmExtraPhysics(const G4String&)
+ : EmExtraPhysics(1)
+{}
+
+EmExtraPhysics::~EmExtraPhysics()
+{}
+
+
+void EmExtraPhysics::Synch(G4bool val)
+{
+ synActivated = val;
+}
+
+void EmExtraPhysics::SynchAll(G4bool val)
+{
+ synActivatedForAll = val;
+ if(synActivatedForAll) { synActivated = true; }
+}
+
+void EmExtraPhysics::GammaNuclear(G4bool val)
+{
+ gnActivated = val;
+}
+
+void EmExtraPhysics::LENDGammaNuclear(G4bool val)
+{
+ gLENDActivated = val;
+}
+
+void EmExtraPhysics::ElectroNuclear(G4bool val)
+{
+ eActivated = val;
+}
+
+void EmExtraPhysics::MuonNuclear(G4bool val)
+{
+ munActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMu(G4bool val)
+{
+ gmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToMuMu(G4bool val)
+{
+ pmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToHadrons(G4bool val)
+{
+ phadActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMuFactor(G4double val)
+{
+ if(val > 0.0) gmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToMuMuFactor(G4double val)
+{
+ if(val > 0.0) pmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToHadronsFactor(G4double val)
+{
+ if(val > 0.0) phadFactor = val;
+}
+
+////////////////////////////////////////////////////
+
+void EmExtraPhysics::NeutrinoActivated(G4bool val)
+{
+ fNuActivated = val;
+}
+
+void EmExtraPhysics::NuETotXscActivated(G4bool val)
+{
+ fNuETotXscActivated = val;
+}
+
+void EmExtraPhysics::SetNuEleCcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleCcBias = bf;
+}
+
+void EmExtraPhysics::SetNuEleNcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleNcBias = bf;
+}
+
+void EmExtraPhysics::SetNuNucleusBias(G4double bf)
+{
+ if(bf > 0.0) fNuNucleusBias = bf;
+}
+
+void EmExtraPhysics::SetNuDetectorName(const G4String& dn)
+{
+ fNuDetectorName = dn;
+}
+
+/////////////////////////////////////////////////
+
+void EmExtraPhysics::ConstructParticle()
+{
+ G4Gamma::Gamma();
+ G4Electron::Electron();
+ G4Positron::Positron();
+ G4MuonPlus::MuonPlus();
+ G4MuonMinus::MuonMinus();
+
+ G4AntiNeutrinoE::AntiNeutrinoE();
+ G4NeutrinoE::NeutrinoE();
+ G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4NeutrinoMu::NeutrinoMu();
+ G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4NeutrinoTau::NeutrinoTau();
+}
+
+void EmExtraPhysics::ConstructProcess()
+{
+ G4ParticleDefinition* gamma = G4Gamma::Gamma();
+ G4ParticleDefinition* electron = G4Electron::Electron();
+ G4ParticleDefinition* positron = G4Positron::Positron();
+ G4ParticleDefinition* muonplus = G4MuonPlus::MuonPlus();
+ G4ParticleDefinition* muonminus = G4MuonMinus::MuonMinus();
+
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+
+ if(gnActivated) { ConstructGammaElectroNuclear(); }
+
+ if(munActivated) {
+ G4MuonNuclearProcess* muNucProcess = new G4MuonNuclearProcess();
+ G4MuonVDNuclearModel* muNucModel = new G4MuonVDNuclearModel();
+ muNucProcess->RegisterMe(muNucModel);
+ ph->RegisterProcess( muNucProcess, muonplus);
+ ph->RegisterProcess( muNucProcess, muonminus);
+ }
+ if(gmumuActivated) {
+ G4GammaConversionToMuons* theGammaToMuMu = new G4GammaConversionToMuons();
+ theGammaToMuMu->SetCrossSecFactor(gmumuFactor);
+ GammaGeneralProcess* sp =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(sp) {
+ sp->AddMMProcess(theGammaToMuMu);
+ } else {
+ ph->RegisterProcess(theGammaToMuMu, gamma);
+ }
+ }
+ if(pmumuActivated) {
+ G4AnnihiToMuPair* thePosiToMuMu = new G4AnnihiToMuPair();
+ thePosiToMuMu->SetCrossSecFactor(pmumuFactor);
+ ph->RegisterProcess(thePosiToMuMu, positron);
+ }
+ if(phadActivated) {
+ G4eeToHadrons* thePosiToHadrons = new G4eeToHadrons();
+ thePosiToHadrons->SetCrossSecFactor(phadFactor);
+ ph->RegisterProcess(thePosiToHadrons, positron);
+ }
+ if(synActivated) {
+ G4SynchrotronRadiation* theSynchRad = new G4SynchrotronRadiation();
+ ph->RegisterProcess( theSynchRad, electron);
+ ph->RegisterProcess( theSynchRad, positron);
+ if(synActivatedForAll) {
+ auto myParticleIterator=GetParticleIterator();
+ myParticleIterator->reset();
+ G4ParticleDefinition* particle = nullptr;
+
+ while( (*myParticleIterator)() ) {
+ particle = myParticleIterator->value();
+ if( particle->GetPDGStable() && particle->GetPDGCharge() != 0.0) {
+ if(verbose > 1) {
+ G4cout << "### G4SynchrotronRadiation for "
+ << particle->GetParticleName() << G4endl;
+ }
+ ph->RegisterProcess( theSynchRad, particle);
+ }
+ }
+ }
+ }
+ if( fNuActivated )
+ {
+ G4ParticleDefinition* anuelectron = G4AntiNeutrinoE::AntiNeutrinoE();
+ G4ParticleDefinition* nuelectron = G4NeutrinoE::NeutrinoE();
+ G4ParticleDefinition* anumuon = G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4ParticleDefinition* numuon = G4NeutrinoMu::NeutrinoMu();
+ G4ParticleDefinition* anutau = G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4ParticleDefinition* nutau = G4NeutrinoTau::NeutrinoTau();
+
+ G4NeutrinoElectronProcess* theNuEleProcess =
+ new G4NeutrinoElectronProcess(fNuDetectorName);
+ G4NeutrinoElectronTotXsc* theNuEleTotXsc = new G4NeutrinoElectronTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ G4double bftot = std::max(fNuEleCcBias,fNuEleNcBias);
+ theNuEleProcess->SetBiasingFactor(bftot);
+ }
+ else
+ {
+ theNuEleProcess->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ theNuEleTotXsc->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ }
+ theNuEleProcess->AddDataSet(theNuEleTotXsc);
+
+ G4NeutrinoElectronCcModel* ccModel = new G4NeutrinoElectronCcModel();
+ G4NeutrinoElectronNcModel* ncModel = new G4NeutrinoElectronNcModel();
+ theNuEleProcess->RegisterMe(ccModel);
+ theNuEleProcess->RegisterMe(ncModel);
+
+ ph->RegisterProcess(theNuEleProcess, anuelectron);
+ ph->RegisterProcess(theNuEleProcess, nuelectron);
+ ph->RegisterProcess(theNuEleProcess, anumuon);
+ ph->RegisterProcess(theNuEleProcess, numuon);
+ ph->RegisterProcess(theNuEleProcess, anutau);
+ ph->RegisterProcess(theNuEleProcess, nutau);
+
+ // nu_mu nucleus interactions
+ G4MuNeutrinoNucleusProcess* theNuMuNucleusProcess = new G4MuNeutrinoNucleusProcess(fNuDetectorName);
+ G4MuNeutrinoNucleusTotXsc* theNuMuNucleusTotXsc = new G4MuNeutrinoNucleusTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ theNuMuNucleusProcess->SetBiasingFactor(fNuNucleusBias);
+ }
+ theNuMuNucleusProcess->AddDataSet(theNuMuNucleusTotXsc);
+
+ G4NuMuNucleusCcModel* numunuclcc = new G4NuMuNucleusCcModel();
+ G4NuMuNucleusNcModel* numunuclnc = new G4NuMuNucleusNcModel();
+
+ theNuMuNucleusProcess->RegisterMe(numunuclcc);
+ theNuMuNucleusProcess->RegisterMe(numunuclnc);
+
+ ph->RegisterProcess(theNuMuNucleusProcess, anumuon);
+ ph->RegisterProcess(theNuMuNucleusProcess, numuon);
+ }
+}
+
+void EmExtraPhysics::ConstructGammaElectroNuclear()
+{
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+
+ G4PhotoNuclearProcess* gnuc = new G4PhotoNuclearProcess();
+
+ G4QGSModel< G4GammaParticipants >* theStringModel =
+ new G4QGSModel< G4GammaParticipants >;
+ G4QGSMFragmentation* theFrag = new G4QGSMFragmentation();
+ G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theFrag);
+ theStringModel->SetFragmentationModel(theStringDecay);
+
+ G4GeneratorPrecompoundInterface* theCascade =
+ new G4GeneratorPrecompoundInterface;
+
+ G4TheoFSGenerator* theModel = new G4TheoFSGenerator();
+ theModel->SetTransport(theCascade);
+ theModel->SetHighEnergyGenerator(theStringModel);
+
+ G4HadronicParameters* param = G4HadronicParameters::Instance();
+
+ G4CascadeInterface* cascade = new G4CascadeInterface;
+ cascade->SetMaxEnergy(param->GetMaxEnergyTransitionFTF_Cascade());
+ gnuc->RegisterMe(cascade);
+ theModel->SetMinEnergy(param->GetMinEnergyTransitionFTF_Cascade());
+ theModel->SetMaxEnergy(param->GetMaxEnergy());
+ gnuc->RegisterMe(theModel);
+
+ GammaGeneralProcess* sp =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(sp) {
+ sp->AddHadProcess(gnuc);
+ } else {
+ // LEND may be activated if the general process is not activated
+ ph->RegisterProcess(gnuc, G4Gamma::Gamma());
+ if(gLENDActivated) { ConstructLENDGammaNuclear(cascade, gnuc); }
+ }
+
+ if(eActivated) {
+ G4ElectronNuclearProcess* enuc = new G4ElectronNuclearProcess;
+ G4PositronNuclearProcess* pnuc = new G4PositronNuclearProcess;
+ G4ElectroVDNuclearModel* eModel = new G4ElectroVDNuclearModel;
+
+ enuc->RegisterMe(eModel);
+ ph->RegisterProcess(enuc, G4Electron::Electron());
+
+ pnuc->RegisterMe(eModel);
+ ph->RegisterProcess(enuc, G4Positron::Positron());
+ }
+}
+
+void EmExtraPhysics::ConstructLENDGammaNuclear(
+ G4CascadeInterface* cascade, G4PhotoNuclearProcess* gnuc)
+{
+ if (std::getenv("G4LENDDATA") == nullptr ) {
+ G4String message = "\n Skipping activation of Low Energy Nuclear Data (LEND) model for gamma nuclear interactions.\n The LEND model needs data files and they are available from ftp://gdo-nuclear.ucllnl.org/GND_after2013/GND_v1.3.tar.gz.\n Please set the environment variable G4LENDDATA to point to the directory named v1.3 extracted from the archive file.\n";
+ G4Exception( "EmExtraPhysics::ConstructLENDGammaNuclear()"
+ , "G4LENDBertiniGammaElectroNuclearBuilder001"
+ , JustWarning , message);
+ return;
+ }
+
+ cascade->SetMinEnergy(19.9*MeV);
+ G4LENDorBERTModel* theGammaReactionLowE =
+ new G4LENDorBERTModel( G4Gamma::Gamma() );
+ theGammaReactionLowE->DumpLENDTargetInfo(true);
+ G4LENDCombinedCrossSection* theGammaCrossSectionLowE =
+ new G4LENDCombinedCrossSection( G4Gamma::Gamma() );
+ theGammaReactionLowE->SetMaxEnergy(20*MeV);
+ gnuc->RegisterMe(theGammaReactionLowE);
+ gnuc->AddDataSet(theGammaCrossSectionLowE);
+}
+
+#elif G4VERSION_NUMBER<1100
+
+#include "G4SystemOfUnits.hh"
+
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+#include "G4AntiNeutrinoE.hh"
+#include "G4NeutrinoE.hh"
+#include "G4AntiNeutrinoMu.hh"
+#include "G4NeutrinoMu.hh"
+#include "G4AntiNeutrinoTau.hh"
+#include "G4NeutrinoTau.hh"
+
+#include "G4Proton.hh"
+#include "G4AntiProton.hh"
+#include "G4PionPlus.hh"
+#include "G4PionMinus.hh"
+#include "G4GenericIon.hh"
+
+#include "G4SynchrotronRadiation.hh"
+#include "G4MuonNuclearProcess.hh"
+#include "G4MuonVDNuclearModel.hh"
+#include "G4ElectroVDNuclearModel.hh"
+#include "G4TheoFSGenerator.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4QGSModel.hh"
+#include "G4GammaParticipants.hh"
+#include "G4QGSMFragmentation.hh"
+#include "G4ExcitedStringDecay.hh"
+#include "G4CascadeInterface.hh"
+#include "G4LowEGammaNuclearModel.hh"
+
+#include "G4LENDorBERTModel.hh"
+#include "G4LENDCombinedCrossSection.hh"
+
+#include "G4GammaConversionToMuons.hh"
+#include "G4AnnihiToMuPair.hh"
+#include "G4eeToHadrons.hh"
+
+#include "G4PhotoNuclearProcess.hh"
+#include "G4ElectronNuclearProcess.hh"
+#include "G4PositronNuclearProcess.hh"
+
+#include "G4NeutrinoElectronProcess.hh"
+#include "G4NeutrinoElectronTotXsc.hh"
+#include "G4NeutrinoElectronCcModel.hh"
+#include "G4NeutrinoElectronNcModel.hh"
+
+#include "G4MuNeutrinoNucleusProcess.hh"
+#include "G4ElNeutrinoNucleusProcess.hh"
+
+#include "G4MuNeutrinoNucleusTotXsc.hh"
+#include "G4ElNeutrinoNucleusTotXsc.hh"
+
+#include "G4NuMuNucleusCcModel.hh"
+#include "G4NuMuNucleusNcModel.hh"
+#include "G4ANuMuNucleusCcModel.hh"
+#include "G4ANuMuNucleusNcModel.hh"
+#include "G4NuElNucleusCcModel.hh"
+#include "G4NuElNucleusNcModel.hh"
+#include "G4ANuElNucleusCcModel.hh"
+#include "G4ANuElNucleusNcModel.hh"
+
+#include "GammaGeneralProcess.hh"
+#include "G4LossTableManager.hh"
+#include "G4GammaNuclearXS.hh"
+
+#include "G4HadronicParameters.hh"
+#include "G4PhysicsListHelper.hh"
+#include "G4BuilderType.hh"
+
+// factory
+#include "G4PhysicsConstructorFactory.hh"
+//
+G4_DECLARE_PHYSCONSTR_FACTORY(EmExtraPhysics);
+
+//////////////////////////////////////
+
+EmExtraPhysics::EmExtraPhysics(G4int ver):
+ G4VPhysicsConstructor("G4GammaLeptoNuclearPhys"),
+ gnActivated (true),
+ eActivated (true),
+ gLENDActivated(false),
+ munActivated(true),
+ synActivated(false),
+ synActivatedForAll(false),
+ gmumuActivated(false),
+ pmumuActivated(false),
+ phadActivated (false),
+ fNuActivated (false),
+ fNuETotXscActivated (false),
+ fUseGammaNuclearXS(false),
+ gmumuFactor (1.0),
+ pmumuFactor (1.0),
+ phadFactor (1.0),
+ fNuEleCcBias(1.0),
+ fNuEleNcBias(1.0),
+ fNuNucleusBias(1.0),
+ fGNLowEnergyLimit(200*CLHEP::MeV),
+ fNuDetectorName("0"),
+ verbose(ver)
+{
+ SetPhysicsType(bEmExtra);
+ if(verbose > 1) G4cout << "### EmExtraPhysics" << G4endl;
+}
+
+EmExtraPhysics::EmExtraPhysics(const G4String&)
+ : EmExtraPhysics(1)
+{}
+
+EmExtraPhysics::~EmExtraPhysics()
+{}
+
+void EmExtraPhysics::Synch(G4bool val)
+{
+ synActivated = val;
+}
+
+void EmExtraPhysics::SynchAll(G4bool val)
+{
+ synActivatedForAll = val;
+ if(synActivatedForAll) { synActivated = true; }
+}
+
+void EmExtraPhysics::GammaNuclear(G4bool val)
+{
+ gnActivated = val;
+}
+
+void EmExtraPhysics::LENDGammaNuclear(G4bool val)
+{
+ gLENDActivated = val;
+ // LEND cannot be used with low-energy model
+ if(val) { fGNLowEnergyLimit = 0.0; }
+}
+
+void EmExtraPhysics::ElectroNuclear(G4bool val)
+{
+ eActivated = val;
+}
+
+void EmExtraPhysics::MuonNuclear(G4bool val)
+{
+ munActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMu(G4bool val)
+{
+ gmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToMuMu(G4bool val)
+{
+ pmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToHadrons(G4bool val)
+{
+ phadActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMuFactor(G4double val)
+{
+ if(val > 0.0) gmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToMuMuFactor(G4double val)
+{
+ if(val > 0.0) pmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToHadronsFactor(G4double val)
+{
+ if(val > 0.0) phadFactor = val;
+}
+
+////////////////////////////////////////////////////
+
+void EmExtraPhysics::NeutrinoActivated(G4bool val)
+{
+ fNuActivated = val;
+}
+
+void EmExtraPhysics::NuETotXscActivated(G4bool val)
+{
+ fNuETotXscActivated = val;
+}
+
+void EmExtraPhysics::SetUseGammaNuclearXS(G4bool val)
+{
+ fUseGammaNuclearXS = val;
+}
+
+void EmExtraPhysics::SetNuEleCcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleCcBias = bf;
+}
+
+void EmExtraPhysics::SetNuEleNcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleNcBias = bf;
+}
+
+void EmExtraPhysics::SetNuNucleusBias(G4double bf)
+{
+ if(bf > 0.0) fNuNucleusBias = bf;
+}
+
+void EmExtraPhysics::GammaNuclearLEModelLimit(G4double val)
+{
+ if(val <= CLHEP::MeV) {
+ fGNLowEnergyLimit = 0.0;
+
+ // lowenergy model should not be applied at high energy
+ } else if(val <= CLHEP::GeV) {
+ fGNLowEnergyLimit = val;
+ gLENDActivated = false;
+ }
+}
+
+void EmExtraPhysics::SetNuDetectorName(const G4String& dn)
+{
+ fNuDetectorName = dn;
+}
+
+/////////////////////////////////////////////////
+
+void EmExtraPhysics::ConstructParticle()
+{
+ G4Gamma::Gamma();
+ G4Electron::Electron();
+ G4Positron::Positron();
+ G4MuonPlus::MuonPlus();
+ G4MuonMinus::MuonMinus();
+
+ G4AntiNeutrinoE::AntiNeutrinoE();
+ G4NeutrinoE::NeutrinoE();
+ G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4NeutrinoMu::NeutrinoMu();
+ G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4NeutrinoTau::NeutrinoTau();
+}
+
+void EmExtraPhysics::ConstructProcess()
+{
+ G4ParticleDefinition* gamma = G4Gamma::Gamma();
+ G4ParticleDefinition* electron = G4Electron::Electron();
+ G4ParticleDefinition* positron = G4Positron::Positron();
+ G4ParticleDefinition* muonplus = G4MuonPlus::MuonPlus();
+ G4ParticleDefinition* muonminus = G4MuonMinus::MuonMinus();
+
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+
+ if(gnActivated) { ConstructGammaElectroNuclear(); }
+
+ if(munActivated) {
+ G4MuonNuclearProcess* muNucProcess = new G4MuonNuclearProcess();
+ G4MuonVDNuclearModel* muNucModel = new G4MuonVDNuclearModel();
+ muNucProcess->RegisterMe(muNucModel);
+ ph->RegisterProcess( muNucProcess, muonplus);
+ ph->RegisterProcess( muNucProcess, muonminus);
+ }
+ if(gmumuActivated) {
+ G4GammaConversionToMuons* theGammaToMuMu = new G4GammaConversionToMuons();
+ theGammaToMuMu->SetCrossSecFactor(gmumuFactor);
+ GammaGeneralProcess* sp =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(sp) {
+ sp->AddMMProcess(theGammaToMuMu);
+ } else {
+ ph->RegisterProcess(theGammaToMuMu, gamma);
+ }
+ }
+ if(pmumuActivated) {
+ G4AnnihiToMuPair* thePosiToMuMu = new G4AnnihiToMuPair();
+ thePosiToMuMu->SetCrossSecFactor(pmumuFactor);
+ ph->RegisterProcess(thePosiToMuMu, positron);
+ }
+ if(phadActivated) {
+ G4eeToHadrons* thePosiToHadrons = new G4eeToHadrons();
+ thePosiToHadrons->SetCrossSecFactor(phadFactor);
+ ph->RegisterProcess(thePosiToHadrons, positron);
+ }
+ if(synActivated) {
+ G4SynchrotronRadiation* theSynchRad = new G4SynchrotronRadiation();
+ ph->RegisterProcess( theSynchRad, electron);
+ ph->RegisterProcess( theSynchRad, positron);
+ if(synActivatedForAll) {
+ ph->RegisterProcess( theSynchRad, muonplus);
+ ph->RegisterProcess( theSynchRad, muonminus);
+
+ ph->RegisterProcess( theSynchRad, G4Proton::Proton());
+ ph->RegisterProcess( theSynchRad, G4AntiProton::AntiProton());
+ ph->RegisterProcess( theSynchRad, G4PionPlus::PionPlus());
+ ph->RegisterProcess( theSynchRad, G4PionMinus::PionMinus());
+ ph->RegisterProcess( theSynchRad, G4GenericIon::GenericIon());
+ }
+ }
+ if( fNuActivated )
+ {
+ G4ParticleDefinition* anuelectron = G4AntiNeutrinoE::AntiNeutrinoE();
+ G4ParticleDefinition* nuelectron = G4NeutrinoE::NeutrinoE();
+ G4ParticleDefinition* anumuon = G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4ParticleDefinition* numuon = G4NeutrinoMu::NeutrinoMu();
+ G4ParticleDefinition* anutau = G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4ParticleDefinition* nutau = G4NeutrinoTau::NeutrinoTau();
+
+ G4NeutrinoElectronProcess* theNuEleProcess =
+ new G4NeutrinoElectronProcess(fNuDetectorName);
+ G4NeutrinoElectronTotXsc* theNuEleTotXsc = new G4NeutrinoElectronTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ G4double bftot = std::max(fNuEleCcBias,fNuEleNcBias);
+ theNuEleProcess->SetBiasingFactor(bftot);
+ }
+ else
+ {
+ theNuEleProcess->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ theNuEleTotXsc->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ }
+ theNuEleProcess->AddDataSet(theNuEleTotXsc);
+
+ G4NeutrinoElectronCcModel* ccModel = new G4NeutrinoElectronCcModel();
+ G4NeutrinoElectronNcModel* ncModel = new G4NeutrinoElectronNcModel();
+ theNuEleProcess->RegisterMe(ccModel);
+ theNuEleProcess->RegisterMe(ncModel);
+
+ ph->RegisterProcess(theNuEleProcess, anuelectron);
+ ph->RegisterProcess(theNuEleProcess, nuelectron);
+ ph->RegisterProcess(theNuEleProcess, anumuon);
+ ph->RegisterProcess(theNuEleProcess, numuon);
+ ph->RegisterProcess(theNuEleProcess, anutau);
+ ph->RegisterProcess(theNuEleProcess, nutau);
+
+ // nu_mu nucleus interactions
+
+ G4MuNeutrinoNucleusProcess* theNuMuNucleusProcess = new G4MuNeutrinoNucleusProcess(fNuDetectorName);
+ G4MuNeutrinoNucleusTotXsc* theNuMuNucleusTotXsc = new G4MuNeutrinoNucleusTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ theNuMuNucleusProcess->SetBiasingFactor(fNuNucleusBias);
+ }
+ theNuMuNucleusProcess->AddDataSet(theNuMuNucleusTotXsc);
+
+ G4NuMuNucleusCcModel* numunuclcc = new G4NuMuNucleusCcModel();
+ G4NuMuNucleusNcModel* numunuclnc = new G4NuMuNucleusNcModel();
+ G4ANuMuNucleusCcModel* anumunuclcc = new G4ANuMuNucleusCcModel();
+ G4ANuMuNucleusNcModel* anumunuclnc = new G4ANuMuNucleusNcModel();
+
+ theNuMuNucleusProcess->RegisterMe(numunuclcc);
+ theNuMuNucleusProcess->RegisterMe(numunuclnc);
+ theNuMuNucleusProcess->RegisterMe(anumunuclcc);
+ theNuMuNucleusProcess->RegisterMe(anumunuclnc);
+
+ ph->RegisterProcess(theNuMuNucleusProcess, anumuon);
+ ph->RegisterProcess(theNuMuNucleusProcess, numuon);
+
+ // nu_e nucleus interactions
+
+ G4ElNeutrinoNucleusProcess* theNuElNucleusProcess = new G4ElNeutrinoNucleusProcess(fNuDetectorName);
+ G4ElNeutrinoNucleusTotXsc* theNuElNucleusTotXsc = new G4ElNeutrinoNucleusTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ theNuElNucleusProcess->SetBiasingFactor(fNuNucleusBias);
+ }
+ theNuElNucleusProcess->AddDataSet(theNuElNucleusTotXsc);
+
+ G4NuElNucleusCcModel* nuelnuclcc = new G4NuElNucleusCcModel();
+ G4NuElNucleusNcModel* nuelnuclnc = new G4NuElNucleusNcModel();
+ G4ANuElNucleusCcModel* anuelnuclcc = new G4ANuElNucleusCcModel();
+ G4ANuElNucleusNcModel* anuelnuclnc = new G4ANuElNucleusNcModel();
+
+ theNuElNucleusProcess->RegisterMe(nuelnuclcc);
+ theNuElNucleusProcess->RegisterMe(nuelnuclnc);
+ theNuElNucleusProcess->RegisterMe(anuelnuclcc);
+ theNuElNucleusProcess->RegisterMe(anuelnuclnc);
+
+ ph->RegisterProcess(theNuElNucleusProcess, anuelectron);
+ ph->RegisterProcess(theNuElNucleusProcess, nuelectron);
+ }
+}
+
+void EmExtraPhysics::ConstructGammaElectroNuclear()
+{
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+
+ G4PhotoNuclearProcess* gnuc = new G4PhotoNuclearProcess();
+ if(fUseGammaNuclearXS) {
+ gnuc->AddDataSet(new G4GammaNuclearXS());
+ }
+
+ G4QGSModel< G4GammaParticipants >* theStringModel =
+ new G4QGSModel< G4GammaParticipants >;
+ G4QGSMFragmentation* theFrag = new G4QGSMFragmentation();
+ G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theFrag);
+ theStringModel->SetFragmentationModel(theStringDecay);
+
+ G4GeneratorPrecompoundInterface* theCascade =
+ new G4GeneratorPrecompoundInterface;
+
+ G4TheoFSGenerator* theModel = new G4TheoFSGenerator();
+ theModel->SetTransport(theCascade);
+ theModel->SetHighEnergyGenerator(theStringModel);
+
+ G4HadronicParameters* param = G4HadronicParameters::Instance();
+
+ G4CascadeInterface* cascade = new G4CascadeInterface;
+
+ // added low-energy model LEND disabled
+ if (fGNLowEnergyLimit > 0.0) {
+ G4LowEGammaNuclearModel* lemod = new G4LowEGammaNuclearModel();
+ lemod->SetMaxEnergy(fGNLowEnergyLimit);
+ gnuc->RegisterMe(lemod);
+ cascade->SetMinEnergy(fGNLowEnergyLimit - CLHEP::MeV);
+ }
+ cascade->SetMaxEnergy(param->GetMaxEnergyTransitionFTF_Cascade());
+ gnuc->RegisterMe(cascade);
+ theModel->SetMinEnergy(param->GetMinEnergyTransitionFTF_Cascade());
+ theModel->SetMaxEnergy(param->GetMaxEnergy());
+ gnuc->RegisterMe(theModel);
+
+ GammaGeneralProcess* gproc =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(gproc != nullptr) {
+ gproc->AddHadProcess(gnuc);
+ } else {
+ // LEND may be activated if the general process is not activated
+ ph->RegisterProcess(gnuc, G4Gamma::Gamma());
+ if(gLENDActivated) { ConstructLENDGammaNuclear(cascade, gnuc); }
+ }
+
+ if(eActivated) {
+ G4ElectronNuclearProcess* enuc = new G4ElectronNuclearProcess;
+ G4PositronNuclearProcess* pnuc = new G4PositronNuclearProcess;
+ G4ElectroVDNuclearModel* eModel = new G4ElectroVDNuclearModel;
+
+ GammaGeneralProcess* eproc =
+ (GammaGeneralProcess*)emManager->GetElectronGeneralProcess();
+ if(eproc != nullptr) {
+ eproc->AddHadProcess(enuc);
+ } else {
+ enuc->RegisterMe(eModel);
+ ph->RegisterProcess(enuc, G4Electron::Electron());
+ }
+
+ GammaGeneralProcess* pproc =
+ (GammaGeneralProcess*)emManager->GetPositronGeneralProcess();
+ if(pproc != nullptr) {
+ pproc->AddHadProcess(pnuc);
+ } else {
+ pnuc->RegisterMe(eModel);
+ ph->RegisterProcess(enuc, G4Positron::Positron());
+ }
+ }
+}
+
+void EmExtraPhysics::ConstructLENDGammaNuclear(
+ G4CascadeInterface* cascade, G4PhotoNuclearProcess* gnuc)
+{
+ if (std::getenv("G4LENDDATA") == nullptr ) {
+ G4String message = "\n Skipping activation of Low Energy Nuclear Data (LEND) model for gamma nuclear interactions.\n The LEND model needs data files and they are available from ftp://gdo-nuclear.ucllnl.org/GND_after2013/GND_v1.3.tar.gz.\n Please set the environment variable G4LENDDATA to point to the directory named v1.3 extracted from the archive file.\n";
+ G4Exception( "EmExtraPhysics::ConstructLENDGammaNuclear()"
+ , "G4LENDBertiniGammaElectroNuclearBuilder001"
+ , JustWarning , message);
+ return;
+ }
+
+ cascade->SetMinEnergy(19.9*MeV);
+ G4LENDorBERTModel* theGammaReactionLowE =
+ new G4LENDorBERTModel( G4Gamma::Gamma() );
+ theGammaReactionLowE->DumpLENDTargetInfo(true);
+ G4LENDCombinedCrossSection* theGammaCrossSectionLowE =
+ new G4LENDCombinedCrossSection( G4Gamma::Gamma() );
+ theGammaReactionLowE->SetMaxEnergy(20*MeV);
+ gnuc->RegisterMe(theGammaReactionLowE);
+ gnuc->AddDataSet(theGammaCrossSectionLowE);
+}
+
+#else
+
+#include "G4SystemOfUnits.hh"
+
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+#include "G4AntiNeutrinoE.hh"
+#include "G4NeutrinoE.hh"
+#include "G4AntiNeutrinoMu.hh"
+#include "G4NeutrinoMu.hh"
+#include "G4AntiNeutrinoTau.hh"
+#include "G4NeutrinoTau.hh"
+
+#include "G4Proton.hh"
+#include "G4AntiProton.hh"
+#include "G4PionPlus.hh"
+#include "G4PionMinus.hh"
+#include "G4GenericIon.hh"
+
+#include "G4SynchrotronRadiation.hh"
+#include "G4MuonNuclearProcess.hh"
+#include "G4MuonVDNuclearModel.hh"
+#include "G4ElectroVDNuclearModel.hh"
+#include "G4TheoFSGenerator.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4QGSModel.hh"
+#include "G4GammaParticipants.hh"
+#include "G4QGSMFragmentation.hh"
+#include "G4ExcitedStringDecay.hh"
+#include "G4CascadeInterface.hh"
+#include "G4LowEGammaNuclearModel.hh"
+
+#include "G4LENDorBERTModel.hh"
+#include "G4LENDCombinedCrossSection.hh"
+
+#include "G4GammaConversionToMuons.hh"
+#include "G4AnnihiToMuPair.hh"
+#include "G4eeToHadrons.hh"
+
+#include "G4HadronInelasticProcess.hh"
+#include "G4ElectronNuclearProcess.hh"
+#include "G4PositronNuclearProcess.hh"
+
+#include "G4NeutrinoElectronProcess.hh"
+#include "G4NeutrinoElectronTotXsc.hh"
+#include "G4NeutrinoElectronCcModel.hh"
+#include "G4NeutrinoElectronNcModel.hh"
+
+#include "G4MuNeutrinoNucleusProcess.hh"
+#include "G4ElNeutrinoNucleusProcess.hh"
+
+#include "G4MuNeutrinoNucleusTotXsc.hh"
+#include "G4ElNeutrinoNucleusTotXsc.hh"
+
+#include "G4NuMuNucleusCcModel.hh"
+#include "G4NuMuNucleusNcModel.hh"
+#include "G4ANuMuNucleusCcModel.hh"
+#include "G4ANuMuNucleusNcModel.hh"
+#include "G4NuElNucleusCcModel.hh"
+#include "G4NuElNucleusNcModel.hh"
+#include "G4ANuElNucleusCcModel.hh"
+#include "G4ANuElNucleusNcModel.hh"
+
+#include "GammaGeneralProcess.hh"
+#include "G4LossTableManager.hh"
+#include "G4PhotoNuclearCrossSection.hh"
+#include "G4GammaNuclearXS.hh"
+
+#include "G4HadronicParameters.hh"
+#include "G4PhysicsListHelper.hh"
+#include "G4BuilderType.hh"
+#include "G4CrossSectionDataSetRegistry.hh"
+
+// factory
+#include "G4PhysicsConstructorFactory.hh"
+//
+G4_DECLARE_PHYSCONSTR_FACTORY(EmExtraPhysics);
+
+//////////////////////////////////////
+
+EmExtraPhysics::EmExtraPhysics(G4int ver):
+ G4VPhysicsConstructor("G4GammaLeptoNuclearPhys"),
+ gnActivated (true),
+ eActivated (true),
+ gLENDActivated(false),
+ munActivated(true),
+ synActivated(false),
+ synActivatedForAll(false),
+ gmumuActivated(false),
+ pmumuActivated(false),
+ phadActivated (false),
+ fNuActivated (false),
+ fNuETotXscActivated (false),
+ fUseGammaNuclearXS(true),
+ gmumuFactor (1.0),
+ pmumuFactor (1.0),
+ phadFactor (1.0),
+ fNuEleCcBias(1.0),
+ fNuEleNcBias(1.0),
+ fNuNucleusBias(1.0),
+ fGNLowEnergyLimit(200*CLHEP::MeV),
+ fNuDetectorName("0"),
+ verbose(ver)
+{
+ SetPhysicsType(bEmExtra);
+ if(verbose > 1) G4cout << "### EmExtraPhysics" << G4endl;
+}
+
+EmExtraPhysics::EmExtraPhysics(const G4String&)
+ : EmExtraPhysics(1)
+{}
+
+EmExtraPhysics::~EmExtraPhysics()
+{}
+
+void EmExtraPhysics::Synch(G4bool val)
+{
+ synActivated = val;
+}
+
+void EmExtraPhysics::SynchAll(G4bool val)
+{
+ synActivatedForAll = val;
+ if(synActivatedForAll) { synActivated = true; }
+}
+
+void EmExtraPhysics::GammaNuclear(G4bool val)
+{
+ gnActivated = val;
+}
+
+void EmExtraPhysics::LENDGammaNuclear(G4bool val)
+{
+ gLENDActivated = val;
+ // LEND cannot be used with low-energy model
+ if(val) { fGNLowEnergyLimit = 0.0; }
+}
+
+void EmExtraPhysics::ElectroNuclear(G4bool val)
+{
+ eActivated = val;
+}
+
+void EmExtraPhysics::MuonNuclear(G4bool val)
+{
+ munActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMu(G4bool val)
+{
+ gmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToMuMu(G4bool val)
+{
+ pmumuActivated = val;
+}
+
+void EmExtraPhysics::PositronToHadrons(G4bool val)
+{
+ phadActivated = val;
+}
+
+void EmExtraPhysics::GammaToMuMuFactor(G4double val)
+{
+ if(val > 0.0) gmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToMuMuFactor(G4double val)
+{
+ if(val > 0.0) pmumuFactor = val;
+}
+
+void EmExtraPhysics::PositronToHadronsFactor(G4double val)
+{
+ if(val > 0.0) phadFactor = val;
+}
+
+////////////////////////////////////////////////////
+
+void EmExtraPhysics::NeutrinoActivated(G4bool val)
+{
+ fNuActivated = val;
+}
+
+void EmExtraPhysics::NuETotXscActivated(G4bool val)
+{
+ fNuETotXscActivated = val;
+}
+
+void EmExtraPhysics::SetUseGammaNuclearXS(G4bool val)
+{
+ fUseGammaNuclearXS = val;
+}
+
+void EmExtraPhysics::SetNuEleCcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleCcBias = bf;
+}
+
+void EmExtraPhysics::SetNuEleNcBias(G4double bf)
+{
+ if(bf > 0.0) fNuEleNcBias = bf;
+}
+
+void EmExtraPhysics::SetNuNucleusBias(G4double bf)
+{
+ if(bf > 0.0) fNuNucleusBias = bf;
+}
+
+void EmExtraPhysics::GammaNuclearLEModelLimit(G4double val)
+{
+ if(val <= CLHEP::MeV) {
+ fGNLowEnergyLimit = 0.0;
+
+ // lowenergy model should not be applied at high energy
+ } else if(val <= CLHEP::GeV) {
+ fGNLowEnergyLimit = val;
+ gLENDActivated = false;
+ }
+}
+
+void EmExtraPhysics::SetNuDetectorName(const G4String& dn)
+{
+ fNuDetectorName = dn;
+}
+
+/////////////////////////////////////////////////
+
+void EmExtraPhysics::ConstructParticle()
+{
+ G4Gamma::Gamma();
+ G4Electron::Electron();
+ G4Positron::Positron();
+ G4MuonPlus::MuonPlus();
+ G4MuonMinus::MuonMinus();
+
+ G4AntiNeutrinoE::AntiNeutrinoE();
+ G4NeutrinoE::NeutrinoE();
+ G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4NeutrinoMu::NeutrinoMu();
+ G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4NeutrinoTau::NeutrinoTau();
+}
+
+void EmExtraPhysics::ConstructProcess()
+{
+ G4ParticleDefinition* gamma = G4Gamma::Gamma();
+ G4ParticleDefinition* electron = G4Electron::Electron();
+ G4ParticleDefinition* positron = G4Positron::Positron();
+ G4ParticleDefinition* muonplus = G4MuonPlus::MuonPlus();
+ G4ParticleDefinition* muonminus = G4MuonMinus::MuonMinus();
+
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+
+ if(gnActivated) { ConstructGammaElectroNuclear(); }
+
+ if(munActivated) {
+ G4MuonNuclearProcess* muNucProcess = new G4MuonNuclearProcess();
+ G4MuonVDNuclearModel* muNucModel = new G4MuonVDNuclearModel();
+ muNucProcess->RegisterMe(muNucModel);
+ ph->RegisterProcess( muNucProcess, muonplus);
+ ph->RegisterProcess( muNucProcess, muonminus);
+ }
+ if(gmumuActivated) {
+ G4GammaConversionToMuons* theGammaToMuMu = new G4GammaConversionToMuons();
+ theGammaToMuMu->SetCrossSecFactor(gmumuFactor);
+ GammaGeneralProcess* sp =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(sp) {
+ sp->AddMMProcess(theGammaToMuMu);
+ } else {
+ ph->RegisterProcess(theGammaToMuMu, gamma);
+ }
+ }
+ if(pmumuActivated) {
+ G4AnnihiToMuPair* thePosiToMuMu = new G4AnnihiToMuPair();
+ thePosiToMuMu->SetCrossSecFactor(pmumuFactor);
+ ph->RegisterProcess(thePosiToMuMu, positron);
+ }
+ if(phadActivated) {
+ G4eeToHadrons* thePosiToHadrons = new G4eeToHadrons();
+ thePosiToHadrons->SetCrossSecFactor(phadFactor);
+ ph->RegisterProcess(thePosiToHadrons, positron);
+ }
+ if(synActivated) {
+ G4SynchrotronRadiation* theSynchRad = new G4SynchrotronRadiation();
+ ph->RegisterProcess( theSynchRad, electron);
+ ph->RegisterProcess( theSynchRad, positron);
+ if(synActivatedForAll) {
+ ph->RegisterProcess( theSynchRad, muonplus);
+ ph->RegisterProcess( theSynchRad, muonminus);
+
+ ph->RegisterProcess( theSynchRad, G4Proton::Proton());
+ ph->RegisterProcess( theSynchRad, G4AntiProton::AntiProton());
+ ph->RegisterProcess( theSynchRad, G4PionPlus::PionPlus());
+ ph->RegisterProcess( theSynchRad, G4PionMinus::PionMinus());
+ ph->RegisterProcess( theSynchRad, G4GenericIon::GenericIon());
+ }
+ }
+ if( fNuActivated )
+ {
+ G4ParticleDefinition* anuelectron = G4AntiNeutrinoE::AntiNeutrinoE();
+ G4ParticleDefinition* nuelectron = G4NeutrinoE::NeutrinoE();
+ G4ParticleDefinition* anumuon = G4AntiNeutrinoMu::AntiNeutrinoMu();
+ G4ParticleDefinition* numuon = G4NeutrinoMu::NeutrinoMu();
+ G4ParticleDefinition* anutau = G4AntiNeutrinoTau::AntiNeutrinoTau();
+ G4ParticleDefinition* nutau = G4NeutrinoTau::NeutrinoTau();
+
+ G4NeutrinoElectronProcess* theNuEleProcess =
+ new G4NeutrinoElectronProcess(fNuDetectorName);
+ G4NeutrinoElectronTotXsc* theNuEleTotXsc = new G4NeutrinoElectronTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ G4double bftot = std::max(fNuEleCcBias,fNuEleNcBias);
+ theNuEleProcess->SetBiasingFactor(bftot);
+ }
+ else
+ {
+ theNuEleProcess->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ theNuEleTotXsc->SetBiasingFactors(fNuEleCcBias,fNuEleNcBias);
+ }
+ theNuEleProcess->AddDataSet(theNuEleTotXsc);
+
+ G4NeutrinoElectronCcModel* ccModel = new G4NeutrinoElectronCcModel();
+ G4NeutrinoElectronNcModel* ncModel = new G4NeutrinoElectronNcModel();
+ theNuEleProcess->RegisterMe(ccModel);
+ theNuEleProcess->RegisterMe(ncModel);
+
+ ph->RegisterProcess(theNuEleProcess, anuelectron);
+ ph->RegisterProcess(theNuEleProcess, nuelectron);
+ ph->RegisterProcess(theNuEleProcess, anumuon);
+ ph->RegisterProcess(theNuEleProcess, numuon);
+ ph->RegisterProcess(theNuEleProcess, anutau);
+ ph->RegisterProcess(theNuEleProcess, nutau);
+
+ // nu_mu nucleus interactions
+
+ G4MuNeutrinoNucleusProcess* theNuMuNucleusProcess = new G4MuNeutrinoNucleusProcess(fNuDetectorName);
+ G4MuNeutrinoNucleusTotXsc* theNuMuNucleusTotXsc = new G4MuNeutrinoNucleusTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ theNuMuNucleusProcess->SetBiasingFactor(fNuNucleusBias);
+ }
+ theNuMuNucleusProcess->AddDataSet(theNuMuNucleusTotXsc);
+
+ G4NuMuNucleusCcModel* numunuclcc = new G4NuMuNucleusCcModel();
+ G4NuMuNucleusNcModel* numunuclnc = new G4NuMuNucleusNcModel();
+ G4ANuMuNucleusCcModel* anumunuclcc = new G4ANuMuNucleusCcModel();
+ G4ANuMuNucleusNcModel* anumunuclnc = new G4ANuMuNucleusNcModel();
+
+ theNuMuNucleusProcess->RegisterMe(numunuclcc);
+ theNuMuNucleusProcess->RegisterMe(numunuclnc);
+ theNuMuNucleusProcess->RegisterMe(anumunuclcc);
+ theNuMuNucleusProcess->RegisterMe(anumunuclnc);
+
+ ph->RegisterProcess(theNuMuNucleusProcess, anumuon);
+ ph->RegisterProcess(theNuMuNucleusProcess, numuon);
+
+ // nu_e nucleus interactions
+
+ G4ElNeutrinoNucleusProcess* theNuElNucleusProcess = new G4ElNeutrinoNucleusProcess(fNuDetectorName);
+ G4ElNeutrinoNucleusTotXsc* theNuElNucleusTotXsc = new G4ElNeutrinoNucleusTotXsc();
+
+ if(fNuETotXscActivated)
+ {
+ theNuElNucleusProcess->SetBiasingFactor(fNuNucleusBias);
+ }
+ theNuElNucleusProcess->AddDataSet(theNuElNucleusTotXsc);
+
+ G4NuElNucleusCcModel* nuelnuclcc = new G4NuElNucleusCcModel();
+ G4NuElNucleusNcModel* nuelnuclnc = new G4NuElNucleusNcModel();
+ G4ANuElNucleusCcModel* anuelnuclcc = new G4ANuElNucleusCcModel();
+ G4ANuElNucleusNcModel* anuelnuclnc = new G4ANuElNucleusNcModel();
+
+ theNuElNucleusProcess->RegisterMe(nuelnuclcc);
+ theNuElNucleusProcess->RegisterMe(nuelnuclnc);
+ theNuElNucleusProcess->RegisterMe(anuelnuclcc);
+ theNuElNucleusProcess->RegisterMe(anuelnuclnc);
+
+ ph->RegisterProcess(theNuElNucleusProcess, anuelectron);
+ ph->RegisterProcess(theNuElNucleusProcess, nuelectron);
+ }
+}
+
+void EmExtraPhysics::ConstructGammaElectroNuclear()
+{
+ G4LossTableManager* emManager = G4LossTableManager::Instance();
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+
+ G4HadronInelasticProcess* gnuc = new G4HadronInelasticProcess( "photonNuclear", G4Gamma::Gamma() );
+ auto xsreg = G4CrossSectionDataSetRegistry::Instance();
+ G4VCrossSectionDataSet* xs = nullptr;
+ if(fUseGammaNuclearXS) {
+ xs = xsreg->GetCrossSectionDataSet("GammaNuclearXS");
+ if(nullptr == xs) xs = new G4GammaNuclearXS();
+ } else {
+ xs = xsreg->GetCrossSectionDataSet("PhotoNuclearXS");
+ if(nullptr == xs) xs = new G4PhotoNuclearCrossSection();
+ }
+ gnuc->AddDataSet(xs);
+
+ G4QGSModel< G4GammaParticipants >* theStringModel =
+ new G4QGSModel< G4GammaParticipants >;
+ G4QGSMFragmentation* theFrag = new G4QGSMFragmentation();
+ G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theFrag);
+ theStringModel->SetFragmentationModel(theStringDecay);
+
+ G4GeneratorPrecompoundInterface* theCascade =
+ new G4GeneratorPrecompoundInterface;
+
+ G4TheoFSGenerator* theModel = new G4TheoFSGenerator();
+ theModel->SetTransport(theCascade);
+ theModel->SetHighEnergyGenerator(theStringModel);
+
+ G4HadronicParameters* param = G4HadronicParameters::Instance();
+
+ G4CascadeInterface* cascade = new G4CascadeInterface;
+
+ // added low-energy model LEND disabled
+ if (fGNLowEnergyLimit > 0.0) {
+ G4LowEGammaNuclearModel* lemod = new G4LowEGammaNuclearModel();
+ lemod->SetMaxEnergy(fGNLowEnergyLimit);
+ gnuc->RegisterMe(lemod);
+ cascade->SetMinEnergy(fGNLowEnergyLimit - CLHEP::MeV);
+ }
+ cascade->SetMaxEnergy(param->GetMaxEnergyTransitionFTF_Cascade());
+ gnuc->RegisterMe(cascade);
+ theModel->SetMinEnergy(param->GetMinEnergyTransitionFTF_Cascade());
+ theModel->SetMaxEnergy(param->GetMaxEnergy());
+ gnuc->RegisterMe(theModel);
+
+ GammaGeneralProcess* gproc =
+ (GammaGeneralProcess*)emManager->GetGammaGeneralProcess();
+ if(gproc != nullptr) {
+ gproc->AddHadProcess(gnuc);
+ } else {
+ // LEND may be activated if the general process is not activated
+ ph->RegisterProcess(gnuc, G4Gamma::Gamma());
+ if(gLENDActivated) { ConstructLENDGammaNuclear(cascade, gnuc); }
+ }
+
+ if(eActivated) {
+ G4ElectronNuclearProcess* enuc = new G4ElectronNuclearProcess;
+ G4PositronNuclearProcess* pnuc = new G4PositronNuclearProcess;
+ G4ElectroVDNuclearModel* eModel = new G4ElectroVDNuclearModel;
+
+ enuc->RegisterMe(eModel);
+ pnuc->RegisterMe(eModel);
+
+ GammaGeneralProcess* eproc =
+ (GammaGeneralProcess*)emManager->GetElectronGeneralProcess();
+ if(eproc != nullptr) {
+ eproc->AddHadProcess(enuc);
+ } else {
+ ph->RegisterProcess(enuc, G4Electron::Electron());
+ }
+
+ GammaGeneralProcess* pproc =
+ (GammaGeneralProcess*)emManager->GetPositronGeneralProcess();
+ if(pproc != nullptr) {
+ pproc->AddHadProcess(pnuc);
+ } else {
+ ph->RegisterProcess(pnuc, G4Positron::Positron());
+ }
+ }
+}
+
+void EmExtraPhysics::ConstructLENDGammaNuclear(
+ G4CascadeInterface* cascade, G4HadronInelasticProcess* gnuc)
+{
+ if (std::getenv("G4LENDDATA") == nullptr ) {
+ G4String message = "\n Skipping activation of Low Energy Nuclear Data (LEND) model for gamma nuclear interactions.\n The LEND model needs data files and they are available from ftp://gdo-nuclear.ucllnl.org/GND_after2013/GND_v1.3.tar.gz.\n Please set the environment variable G4LENDDATA to point to the directory named v1.3 extracted from the archive file.\n";
+ G4Exception( "EmExtraPhysics::ConstructLENDGammaNuclear()"
+ , "G4LENDBertiniGammaElectroNuclearBuilder001"
+ , JustWarning , message);
+ return;
+ }
+
+ cascade->SetMinEnergy(19.9*MeV);
+ G4LENDorBERTModel* theGammaReactionLowE =
+ new G4LENDorBERTModel( G4Gamma::Gamma() );
+ theGammaReactionLowE->DumpLENDTargetInfo(true);
+ G4LENDCombinedCrossSection* theGammaCrossSectionLowE =
+ new G4LENDCombinedCrossSection( G4Gamma::Gamma() );
+ theGammaReactionLowE->SetMaxEnergy(20*MeV);
+ gnuc->RegisterMe(theGammaReactionLowE);
+ gnuc->AddDataSet(theGammaCrossSectionLowE);
+}
+
+#endif
diff --git a/FullSimLight/src/GammaGeneralProcess.cc b/FullSimLight/src/GammaGeneralProcess.cc
new file mode 100644
index 0000000000000000000000000000000000000000..873043554c972cd0ab6e3550b95645ff4e204189
--- /dev/null
+++ b/FullSimLight/src/GammaGeneralProcess.cc
@@ -0,0 +1,796 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name: GammaGeneralProcess
+//
+// Author: Vladimir Ivanchenko
+//
+// Creation date: 19.07.2018
+//
+// Modifications:
+//
+// Class Description:
+//
+// NOTE: M. Novak: 27.11.2021
+//
+// It is a local version of the `GammaGeneralProcess` that will be available
+// in Geant4-11.0. More exactly, a version that can be utilised on the top
+// Geant4-11.0-beta as a base class for the `WoodcockProcess` (also local). When
+// using Geant4-11.0 or later, it is recommended to remove this local version
+// and derive the `WoodcockProcess` directly from the `GammaGeneralProcess`
+// delivered by the toolkit.
+
+// -------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#include "GammaGeneralProcess.hh"
+
+#include "G4VDiscreteProcess.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ProcessManager.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4LossTableBuilder.hh"
+#include "G4HadronicProcess.hh"
+#include "G4LossTableManager.hh"
+#include "G4Step.hh"
+#include "G4Track.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4DataVector.hh"
+#include "G4PhysicsTable.hh"
+#include "G4PhysicsLogVector.hh"
+#include "G4PhysicsLinearVector.hh"
+#include "G4VParticleChange.hh"
+#include "G4PhysicsTableHelper.hh"
+#include "G4EmParameters.hh"
+#include "G4EmProcessSubType.hh"
+#include "G4Material.hh"
+#include "G4MaterialCutsCouple.hh"
+#include "G4GammaConversionToMuons.hh"
+#include "G4Gamma.hh"
+
+#include "G4Version.hh"
+#include "G4Log.hh"
+#include <iostream>
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4EmDataHandler* GammaGeneralProcess::theHandler = nullptr;
+G4bool GammaGeneralProcess::theT[nTables] =
+ {true,false,true,true,true,false,true,true,true,
+ true,true,true,true,true,true};
+G4String GammaGeneralProcess::nameT[nTables] =
+ {"0","1","2","3","4","5","6","7","8",
+ "9","10","11","12","13","14"};
+
+GammaGeneralProcess::GammaGeneralProcess(const G4String& pname):
+ G4VEmProcess(pname, fElectromagnetic),
+ minPEEnergy(150*CLHEP::keV),
+ minEEEnergy(2*CLHEP::electron_mass_c2),
+ minMMEnergy(100*CLHEP::MeV)
+{
+ SetVerboseLevel(1);
+ SetParticle(G4Gamma::Gamma());
+ SetProcessSubType(fGammaGeneralProcess);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+GammaGeneralProcess::~GammaGeneralProcess()
+{
+ if(isTheMaster) {
+ delete theHandler;
+ theHandler = nullptr;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool GammaGeneralProcess::IsApplicable(const G4ParticleDefinition&)
+{
+ return true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::AddEmProcess(G4VEmProcess* ptr)
+{
+ if(nullptr == ptr) { return; }
+ G4int stype = ptr->GetProcessSubType();
+ if(stype == fRayleigh) { theRayleigh = ptr; }
+ else if(stype == fPhotoElectricEffect) { thePhotoElectric = ptr; }
+ else if(stype == fComptonScattering) { theCompton = ptr; }
+ else if(stype == fGammaConversion) { theConversionEE = ptr; }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::AddMMProcess(G4GammaConversionToMuons* ptr)
+{
+ theConversionMM = ptr;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::AddHadProcess(G4HadronicProcess* ptr)
+{
+ theGammaNuclear = ptr;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
+{
+ SetParticle(&part);
+ preStepLambda = 0.0;
+ idxEnergy = 0;
+ currentCouple = nullptr;
+ currentMaterial = nullptr;
+
+ G4EmParameters* param = G4EmParameters::Instance();
+ G4LossTableManager* man = G4LossTableManager::Instance();
+
+ isTheMaster = man->IsMaster();
+ if(isTheMaster) { SetVerboseLevel(param->Verbose()); }
+ else { SetVerboseLevel(param->WorkerVerbose()); }
+
+// G4LossTableBuilder* bld = man->GetTableBuilder();
+// baseMat = bld->GetBaseMaterialFlag();
+
+ if(1 < verboseLevel) {
+ G4cout << "local-GammaGeneralProcess::PreparePhysicsTable() for "
+ << GetProcessName()
+ << " and particle " << part.GetParticleName()
+ << " isMaster: " << isTheMaster << G4endl;
+ }
+
+ // 3 sub-processes must be always defined
+ if(thePhotoElectric == nullptr || theCompton == nullptr ||
+ theConversionEE == nullptr) {
+ G4ExceptionDescription ed;
+ ed << "### local-GeneralGammaProcess is initialized incorrectly"
+ << "\n Photoelectric: " << thePhotoElectric
+ << "\n Compton: " << theCompton
+ << "\n Conversion: " << theConversionEE;
+ G4Exception("G4GeneralGammaProcess","em0004",
+ FatalException, ed,"");
+ }
+
+ if(thePhotoElectric) { thePhotoElectric->PreparePhysicsTable(part); }
+ if(theCompton) { theCompton->PreparePhysicsTable(part); }
+ if(theConversionEE) { theConversionEE->PreparePhysicsTable(part); }
+ if(theRayleigh) { theRayleigh->PreparePhysicsTable(part); }
+ if(theGammaNuclear) { theGammaNuclear->PreparePhysicsTable(part); }
+ if(theConversionMM) { theConversionMM->PreparePhysicsTable(part); }
+
+ InitialiseProcess(&part);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::InitialiseProcess(const G4ParticleDefinition*)
+{
+ if(isTheMaster) {
+
+ G4EmParameters* param = G4EmParameters::Instance();
+ G4LossTableManager* man = G4LossTableManager::Instance();
+
+ // tables are created and its size is defined only once
+ if(nullptr == theHandler) {
+ theHandler = new G4EmDataHandler(nTables);
+ if(theRayleigh) { theT[1] = true; }
+ theHandler->SetMasterProcess(thePhotoElectric);
+ theHandler->SetMasterProcess(theCompton);
+ theHandler->SetMasterProcess(theConversionEE);
+ theHandler->SetMasterProcess(theRayleigh);
+ }
+ auto bld = man->GetTableBuilder();
+
+ const G4ProductionCutsTable* theCoupleTable=
+ G4ProductionCutsTable::GetProductionCutsTable();
+ size_t numOfCouples = theCoupleTable->GetTableSize();
+
+ G4double mine = param->MinKinEnergy();
+ G4double maxe = param->MaxKinEnergy();
+ G4int nd = param->NumberOfBinsPerDecade();
+ size_t nbin1 = std::max(5, nd*G4lrint(std::log10(minPEEnergy/mine)));
+ size_t nbin2 = std::max(5, nd*G4lrint(std::log10(maxe/minMMEnergy)));
+
+ G4PhysicsVector* vec = nullptr;
+#if G4VERSION_NUMBER>=1100
+ G4PhysicsLogVector aVector(mine,minPEEnergy,nbin1,splineFlag);
+ G4PhysicsLogVector bVector(minPEEnergy,minEEEnergy,nLowE,splineFlag);
+ G4PhysicsLogVector cVector(minEEEnergy,minMMEnergy,nHighE,splineFlag);
+ G4PhysicsLogVector dVector(minMMEnergy,maxe,nbin2,splineFlag);
+#else
+ G4PhysicsLogVector aVector(mine,minPEEnergy,nbin1);
+ G4PhysicsLogVector bVector(minPEEnergy,minEEEnergy,nLowE);
+ G4PhysicsLogVector cVector(minEEEnergy,minMMEnergy,nHighE);
+ G4PhysicsLogVector dVector(minMMEnergy,maxe,nbin2);
+ if(splineFlag) {
+ aVector.SetSpline(splineFlag);
+ bVector.SetSpline(splineFlag);
+ cVector.SetSpline(splineFlag);
+ dVector.SetSpline(splineFlag);
+ }
+#endif
+ for(size_t i=0; i<nTables; ++i) {
+ if(!theT[i]) { continue; }
+ //G4cout << "## PreparePhysTable " << i << "." << G4endl;
+ G4PhysicsTable* table = theHandler->MakeTable(i);
+ //G4cout << " make table " << table << G4endl;
+ for(size_t j=0; j<numOfCouples; ++j) {
+ vec = (*table)[j];
+ if (bld->GetFlag(j) && nullptr == vec) {
+ //G4cout <<" i= "<<i<<" j= "<< j <<" make new vector"<< G4endl;
+ if(i<=1) {
+ vec = new G4PhysicsVector(aVector);
+ } else if(i<=5) {
+ vec = new G4PhysicsVector(bVector);
+ } else if(i<=9) {
+ vec = new G4PhysicsVector(cVector);
+ } else {
+ vec = new G4PhysicsVector(dVector);
+ }
+ G4PhysicsTableHelper::SetPhysicsVector(table, j, vec);
+ }
+ }
+ }
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
+{
+ if(1 < verboseLevel) {
+ G4cout << "### G4VEmProcess::BuildPhysicsTable() for "
+ << GetProcessName()
+ << " and particle " << part.GetParticleName()
+ << G4endl;
+ }
+ if(!isTheMaster) {
+ thePhotoElectric->SetEmMasterProcess(theHandler->GetMasterProcess(0));
+// baseMat = theHandler->GetMasterProcess(0)->UseBaseMaterial();
+ }
+ thePhotoElectric->BuildPhysicsTable(part);
+
+ if(!isTheMaster) {
+ theCompton->SetEmMasterProcess(theHandler->GetMasterProcess(1));
+ }
+ theCompton->BuildPhysicsTable(part);
+
+ if(!isTheMaster) {
+ theConversionEE->SetEmMasterProcess(theHandler->GetMasterProcess(2));
+ }
+ theConversionEE->BuildPhysicsTable(part);
+
+ if(theRayleigh != nullptr) {
+ if(!isTheMaster) {
+ theRayleigh->SetEmMasterProcess(theHandler->GetMasterProcess(3));
+ }
+ theRayleigh->BuildPhysicsTable(part);
+ }
+ if(theGammaNuclear != nullptr) { theGammaNuclear->BuildPhysicsTable(part); }
+ if(theConversionMM != nullptr) { theConversionMM->BuildPhysicsTable(part); }
+
+ if(isTheMaster) {
+ const G4ProductionCutsTable* theCoupleTable=
+ G4ProductionCutsTable::GetProductionCutsTable();
+ size_t numOfCouples = theCoupleTable->GetTableSize();
+
+ G4LossTableBuilder* bld = G4LossTableManager::Instance()->GetTableBuilder();
+ const std::vector<G4PhysicsTable*>& tables = theHandler->GetTables();
+
+ G4CrossSectionDataStore* gn = (nullptr != theGammaNuclear)
+ ? theGammaNuclear->GetCrossSectionDataStore() : nullptr;
+ G4DynamicParticle* dynParticle =
+ new G4DynamicParticle(G4Gamma::Gamma(),G4ThreeVector(1,0,0),1.0);
+
+ G4double sigComp(0.), sigPE(0.), sigConv(0.), sigR(0.),
+ sigN(0.), sigM(0.), val(0.);
+
+ for(size_t i=0; i<numOfCouples; ++i) {
+
+ if (bld->GetFlag(i)) {
+
+// G4int idx = (!baseMat) ? i : DensityIndex(i);
+ G4int idx = (*theDensityIdx)[i];
+ const G4MaterialCutsCouple* couple =
+ theCoupleTable->GetMaterialCutsCouple(i);
+ const G4Material* material = couple->GetMaterial();
+
+ // energy interval 0
+ size_t nn = (*(tables[0]))[idx]->GetVectorLength();
+ if(1 < verboseLevel) {
+ G4cout << "======= Zone 0 ======= N= " << nn
+ << " for " << material->GetName() << G4endl;
+ }
+ for(size_t j=0; j<nn; ++j) {
+ G4double e = (*(tables[0]))[idx]->Energy(j);
+ G4double loge = G4Log(e);
+ sigComp = theCompton->GetLambda(e, couple, loge);
+ sigR = (nullptr != theRayleigh) ?
+ theRayleigh->GetLambda(e, couple, loge) : 0.0;
+ G4double sum = sigComp + sigR;
+ if(1 < verboseLevel) {
+ G4cout << j << ". E= " << e << " xs= " << sum
+ << " compt= " << sigComp << " Rayl= " << sigR << G4endl;
+ }
+ (*(tables[0]))[idx]->PutValue(j, sum);
+ if(theT[1]) {
+ val = sigR/sum;
+ (*(tables[1]))[idx]->PutValue(j, val);
+ }
+ }
+
+ // energy interval 1
+ nn = (*(tables[2]))[idx]->GetVectorLength();
+ if(1 < verboseLevel) {
+ G4cout << "======= Zone 1 ======= N= " << nn << G4endl;
+ }
+ for(size_t j=0; j<nn; ++j) {
+ G4double e = (*(tables[2]))[idx]->Energy(j);
+ G4double loge = G4Log(e);
+ sigComp = theCompton->GetLambda(e, couple, loge);
+ sigR = (nullptr != theRayleigh) ? theRayleigh->GetLambda(e, couple, loge) : 0.0;
+ sigPE = thePhotoElectric->GetLambda(e, couple, loge);
+ G4double sum = sigComp + sigR + sigPE;
+ if(1 < verboseLevel) {
+ G4cout << j << ". E= " << e << " xs= " << sum
+ << " compt= " << sigComp << " conv= " << sigConv
+ << " PE= " << sigPE << " Rayl= " << sigR
+ << " GN= " << sigN << G4endl;
+ }
+ (*(tables[2]))[idx]->PutValue(j, sum);
+
+ val = sigPE/sum;
+ (*(tables[3]))[idx]->PutValue(j, val);
+
+ val = (sigR > 0.0) ? (sigComp + sigPE)/sum : 1.0;
+ (*(tables[4]))[idx]->PutValue(j, val);
+ }
+
+ // energy interval 2
+ nn = (*(tables[6]))[idx]->GetVectorLength();
+ if(1 < verboseLevel) {
+ G4cout << "======= Zone 2 ======= N= " << nn << G4endl;
+ }
+ for(size_t j=0; j<nn; ++j) {
+ G4double e = (*(tables[6]))[idx]->Energy(j);
+ G4double loge = G4Log(e);
+ sigComp = theCompton->GetLambda(e, couple, loge);
+ sigConv = theConversionEE->GetLambda(e, couple, loge);
+ sigPE = thePhotoElectric->GetLambda(e, couple, loge);
+ sigN = 0.0;
+ if(nullptr != gn) {
+ dynParticle->SetKineticEnergy(e);
+ sigN = gn->ComputeCrossSection(dynParticle, material);
+ }
+ G4double sum = sigComp + sigConv + sigPE + sigN;
+ if(1 < verboseLevel) {
+ G4cout << j << ". E= " << e << " xs= " << sum
+ << " compt= " << sigComp << " conv= " << sigConv
+ << " PE= " << sigPE
+ << " GN= " << sigN << G4endl;
+ }
+ (*(tables[6]))[idx]->PutValue(j, sum);
+
+ val = sigConv/sum;
+ (*(tables[7]))[idx]->PutValue(j, val);
+
+ val = (sigConv + sigComp)/sum;
+ (*(tables[8]))[idx]->PutValue(j, val);
+
+ val = (sigN > 0.0) ? (sigConv + sigComp + sigPE)/sum : 1.0;
+ (*(tables[9]))[idx]->PutValue(j, val);
+ }
+
+ // energy interval 3
+ nn = (*(tables[10]))[idx]->GetVectorLength();
+ if(1 < verboseLevel) {
+ G4cout << "======= Zone 3 ======= N= " << nn
+ << " for " << material->GetName() << G4endl;
+ }
+ for(size_t j=0; j<nn; ++j) {
+ G4double e = (*(tables[10]))[idx]->Energy(j);
+ G4double loge = G4Log(e);
+ sigComp = theCompton->GetLambda(e, couple, loge);
+ sigConv = theConversionEE->GetLambda(e, couple, loge);
+ sigPE = thePhotoElectric->GetLambda(e, couple, loge);
+ sigN = 0.0;
+ if(nullptr != gn) {
+ dynParticle->SetKineticEnergy(e);
+ sigN = gn->ComputeCrossSection(dynParticle, material);
+ }
+ sigM = 0.0;
+ if(nullptr != theConversionMM) {
+ val = theConversionMM->ComputeMeanFreePath(e, material);
+ sigM = (val < DBL_MAX) ? 1./val : 0.0;
+ }
+ G4double sum = sigComp + sigConv + sigPE + sigN + sigM;
+ if(1 < verboseLevel) {
+ G4cout << j << ". E= " << e << " xs= " << sum
+ << " compt= " << sigComp << " conv= " << sigConv
+ << " PE= " << sigPE
+ << " GN= " << sigN << G4endl;
+ }
+ (*(tables[10]))[idx]->PutValue(j, sum);
+
+ val = (sigComp + sigPE + sigN + sigM)/sum;
+ (*(tables[11]))[idx]->PutValue(j, val);
+
+ val = (sigPE + sigN + sigM)/sum;
+ (*(tables[12]))[idx]->PutValue(j, val);
+
+ val = (sigN + sigM)/sum;
+ (*(tables[13]))[idx]->PutValue(j, val);
+
+ val = sigN/sum;
+ (*(tables[14]))[idx]->PutValue(j, val);
+ }
+ for(size_t k=0; k<nTables; ++k) {
+ if(splineFlag) {
+ (*(tables[k]))[idx]->FillSecondDerivatives();
+ }
+ }
+ }
+ }
+ delete dynParticle;
+ }
+
+ if(1 < verboseLevel) {
+ G4cout << "### G4VEmProcess::BuildPhysicsTable() done for "
+ << GetProcessName()
+ << " and particle " << part.GetParticleName()
+ << G4endl;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::StartTracking(G4Track*)
+{
+ theNumberOfInteractionLengthLeft = -1.0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double GammaGeneralProcess::PostStepGetPhysicalInteractionLength(
+ const G4Track& track,
+ G4double previousStepSize,
+ G4ForceCondition* condition)
+{
+ *condition = NotForced;
+ G4double x = DBL_MAX;
+
+ G4double energy = track.GetKineticEnergy();
+ const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
+
+ // compute mean free path
+ G4bool recompute = false;
+ if(couple != currentCouple) {
+ currentCouple = couple;
+ currentCoupleIndex = couple->GetIndex();
+ basedCoupleIndex = (*theDensityIdx)[currentCoupleIndex];
+ factor = (*theDensityFactor)[currentCoupleIndex];
+ currentMaterial = couple->GetMaterial();
+ recompute = true;
+ }
+ if(energy != preStepKinEnergy) {
+ preStepKinEnergy = energy;
+ preStepLogE = track.GetDynamicParticle()->GetLogKineticEnergy();
+ recompute = true;
+ }
+ if(recompute) {
+ preStepLambda = TotalCrossSectionPerVolume();
+
+ // zero cross section
+ if(preStepLambda <= 0.0) {
+ theNumberOfInteractionLengthLeft = -1.0;
+ currentInteractionLength = DBL_MAX;
+ }
+ }
+
+ // non-zero cross section
+ if(preStepLambda > 0.0) {
+
+ if (theNumberOfInteractionLengthLeft < 0.0) {
+
+ // beggining of tracking (or just after DoIt of this process)
+ theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
+ theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
+
+ } else if(currentInteractionLength < DBL_MAX) {
+
+ theNumberOfInteractionLengthLeft -=
+ previousStepSize/currentInteractionLength;
+ theNumberOfInteractionLengthLeft =
+ std::max(theNumberOfInteractionLengthLeft, 0.0);
+ }
+
+ // new mean free path and step limit for the next step
+ currentInteractionLength = 1.0/preStepLambda;
+ x = theNumberOfInteractionLengthLeft * currentInteractionLength;
+ }
+ /*
+ G4cout << "PostStepGetPhysicalInteractionLength: e= " << energy
+ << " idxe= " << idxEnergy << " xs= " << preStepLambda
+ << " x= " << x << G4endl;
+ */
+ return x;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double GammaGeneralProcess::TotalCrossSectionPerVolume()
+{
+ G4double cross = 0.0;
+
+// G4cout << "#Total: " << preStepKinEnergy << " " << minPEEnergy << " "
+// << minEEEnergy << " " << minMMEnergy<< G4endl;
+// G4cout << " idxE= " << idxEnergy
+// << " idxC= " << currentCoupleIndex << G4endl;
+
+ if(preStepKinEnergy < minPEEnergy) {
+ cross = ComputeGeneralLambda(0, 0);
+ peLambda = thePhotoElectric->GetLambda(preStepKinEnergy, currentCouple, preStepLogE);
+ cross += peLambda;
+
+ } else if(preStepKinEnergy < minEEEnergy) {
+ cross = ComputeGeneralLambda(1, 2);
+
+ } else if(preStepKinEnergy < minMMEnergy) {
+ cross = ComputeGeneralLambda(2, 6);
+ //G4cout << "XS4: " << cross << G4endl;
+
+ } else {
+ cross = ComputeGeneralLambda(3, 10);
+ //G4cout << "XS5: " << cross << G4endl;
+ }
+ /*
+ G4cout << "xs= " << cross << " idxE= " << idxEnergy
+ << " idxC= " << currentCoupleIndex
+ << " E= " << energy << G4endl;
+ */
+ return cross;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4VParticleChange* GammaGeneralProcess::PostStepDoIt(const G4Track& track,
+ const G4Step& step)
+{
+ // In all cases clear number of interaction lengths
+ theNumberOfInteractionLengthLeft = -1.0;
+ selectedProc = nullptr;
+ G4double q = G4UniformRand();
+ /*
+ G4cout << "PostStep: preStepLambda= " << preStepLambda
+ << " PE= " << peLambda << " q= " << q << " idxE= " << idxEnergy
+ << G4endl;
+ */
+ switch (idxEnergy) {
+ case 0:
+ if(preStepLambda*q <= peLambda) {
+ SelectEmProcess(step, thePhotoElectric);
+ } else {
+ if(theT[1] && preStepLambda*q < preStepLambda*GetProbability(1) + peLambda) {
+ SelectEmProcess(step, theRayleigh);
+ } else {
+ SelectEmProcess(step, theCompton);
+ }
+ }
+ break;
+
+ case 1:
+ if(q <= GetProbability(3)) {
+ SelectEmProcess(step, thePhotoElectric);
+ } else if(q <= GetProbability(4)) {
+ SelectEmProcess(step, theCompton);
+ } else if(theRayleigh) {
+ SelectEmProcess(step, theRayleigh);
+ }
+ break;
+
+ case 2:
+ if(q <= GetProbability(7)) {
+ SelectEmProcess(step, theConversionEE);
+ } else if(q <= GetProbability(8)) {
+ SelectEmProcess(step, theCompton);
+ } else if(q <= GetProbability(9)) {
+ SelectEmProcess(step, thePhotoElectric);
+ } else if(theGammaNuclear) {
+ SelectHadProcess(track, step, theGammaNuclear);
+ }
+ break;
+
+ case 3:
+ if(q + GetProbability(11) <= 1.0) {
+ SelectEmProcess(step, theConversionEE);
+ } else if(q + GetProbability(12) <= 1.0) {
+ SelectEmProcess(step, theCompton);
+ } else if(q + GetProbability(13) <= 1.0) {
+ SelectEmProcess(step, thePhotoElectric);
+ } else if(theGammaNuclear && q + GetProbability(14) <= 1.0) {
+ SelectHadProcess(track, step, theGammaNuclear);
+ } else if(theConversionMM) {
+ SelectedProcess(step, theConversionMM);
+ }
+ break;
+ }
+ // sample secondaries
+ if(selectedProc != nullptr) {
+ return selectedProc->PostStepDoIt(track, step);
+ }
+ // no interaction - exception case
+ fParticleChange.InitializeForPostStep(track);
+ return &fParticleChange;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::SelectHadProcess(const G4Track& track,
+ const G4Step& step, G4HadronicProcess* proc)
+{
+ SelectedProcess(step, proc);
+ proc->GetCrossSectionDataStore()->ComputeCrossSection(track.GetDynamicParticle(),
+ currentMaterial);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool GammaGeneralProcess::StorePhysicsTable(const G4ParticleDefinition* part,
+ const G4String& directory,
+ G4bool ascii)
+{
+ G4bool yes = true;
+ if(!isTheMaster) { return yes; }
+ if(!thePhotoElectric->StorePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(!theCompton->StorePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(!theConversionEE->StorePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(theRayleigh != nullptr &&
+ !theRayleigh->StorePhysicsTable(part, directory, ascii))
+ { yes = false; }
+
+ for(size_t i=0; i<nTables; ++i) {
+ if(theT[i]) {
+ G4String nam = (0==i || 2==i || 6==i || 10==i)
+ ? "LambdaGeneral" + nameT[i] : "ProbGeneral" + nameT[i];
+ G4String fnam = GetPhysicsTableFileName(part,directory,nam,ascii);
+ if(!theHandler->StorePhysicsTable(i, part, fnam, ascii)) { yes = false; }
+ }
+ }
+ return yes;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+// M . Novak: 27.11.2021
+// Deactivate `RetrievePhysicsTable` since it might require the splineFlag that
+// is not avaibale in case of earlier G4 versions (not in this way).
+G4bool
+GammaGeneralProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
+ const G4String& /*directory*/,
+ G4bool /*ascii*/)
+{
+ if(1 < verboseLevel) {
+ G4cout << "GammaGeneralProcess::RetrievePhysicsTable() for "
+ << part->GetParticleName() << " and process "
+ << GetProcessName() << " HAS BEN DEACTIVATED!!! "<< G4endl;
+ }
+ G4bool yes = true;
+/*
+ if(!thePhotoElectric->RetrievePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(!theCompton->RetrievePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(!theConversionEE->RetrievePhysicsTable(part, directory, ascii))
+ { yes = false; }
+ if(theRayleigh != nullptr &&
+ !theRayleigh->RetrievePhysicsTable(part, directory, ascii))
+ { yes = false; }
+
+ for(size_t i=0; i<nTables; ++i) {
+ if(theT[i]) {
+ G4String nam = (0==i || 2==i || 6==i || 10==i)
+ ? "LambdaGeneral" + nameT[i] : "ProbGeneral" + nameT[i];
+ G4String fnam = GetPhysicsTableFileName(part,directory,nam,ascii);
+ if(!theHandler->RetrievePhysicsTable(i, part, fnam, ascii, splineFlag))
+ { yes = false; }
+ }
+ }
+ if(yes) {
+ }
+*/
+ return yes;
+}
+
+//....Ooooo0ooooo ........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double GammaGeneralProcess::GetMeanFreePath(const G4Track& track,
+ G4double,
+ G4ForceCondition* condition)
+{
+ *condition = NotForced;
+ return MeanFreePath(track);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void GammaGeneralProcess::ProcessDescription(std::ostream& out) const
+{
+ thePhotoElectric->ProcessDescription(out);
+ theCompton->ProcessDescription(out);
+ theConversionEE->ProcessDescription(out);
+ if(theRayleigh) { theRayleigh->ProcessDescription(out); }
+ if(theGammaNuclear) { theGammaNuclear->ProcessDescription(out); }
+ if(theConversionMM) { theConversionMM->ProcessDescription(out); }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+const G4String& GammaGeneralProcess::GetSubProcessName() const
+{
+ return (selectedProc) ? selectedProc->GetProcessName()
+ : G4VProcess::GetProcessName();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4int GammaGeneralProcess::GetSubProcessSubType() const
+{
+ return (selectedProc) ? selectedProc->GetProcessSubType()
+ : fGammaGeneralProcess;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4VEmProcess* GammaGeneralProcess::GetEmProcess(const G4String& name)
+{
+ G4VEmProcess* proc = nullptr;
+ if(name == thePhotoElectric->GetProcessName()) {
+ proc = thePhotoElectric;
+ } else if(name == theCompton->GetProcessName()) {
+ proc = theCompton;
+ } else if(name == theConversionEE->GetProcessName()) {
+ proc = theConversionEE;
+ } else if(theRayleigh != nullptr && name == theRayleigh->GetProcessName()) {
+ proc = theRayleigh;
+ }
+ return proc;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
diff --git a/FullSimLight/src/MyActionInitialization.cc b/FullSimLight/src/MyActionInitialization.cc
index 3bcd4ed7680e6ca494230db0a914fac65fb84d98..1276b78216d4ee1168b37ea6134aef181143c8e9 100644
--- a/FullSimLight/src/MyActionInitialization.cc
+++ b/FullSimLight/src/MyActionInitialization.cc
@@ -22,7 +22,8 @@
//const G4AnalysisManager* MyActionInitialization::fMasterAnalysisManager = nullptr;
MyActionInitialization::MyActionInitialization(bool isperformance, bool createGeantinoMaps, G4String geantinoMapsFilename)
-: G4VUserActionInitialization(), fIsPerformance(isperformance), fCreateGeantinoMaps(createGeantinoMaps),fGeantinoMapsFilename(geantinoMapsFilename){}
+: G4VUserActionInitialization(), fIsPerformance(isperformance), fCreateGeantinoMaps(createGeantinoMaps),fGeantinoMapsFilename(geantinoMapsFilename),
+ fSpecialScoringRegionName("") {}
MyActionInitialization::~MyActionInitialization() {}
@@ -31,6 +32,7 @@ MyActionInitialization::~MyActionInitialization() {}
void MyActionInitialization::BuildForMaster() const {
MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps,fGeantinoMapsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
+ masterRunAct->SetSpecialScoringRegionName(fSpecialScoringRegionName);
#if USE_PYTHIA
if (use_pythia()) {
G4String str(get_pythia_config());
@@ -62,6 +64,7 @@ void MyActionInitialization::Build() const {
if (fIsPerformance) {
MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
+ masterRunAct->SetSpecialScoringRegionName(fSpecialScoringRegionName);
#if USE_PYTHIA
if (use_pythia()) {
G4String str(get_pythia_config());
@@ -75,12 +78,17 @@ void MyActionInitialization::Build() const {
if (!fIsPerformance) {
MyRunAction* runact = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
SetUserAction(runact);
+ runact->SetSpecialScoringRegionName(fSpecialScoringRegionName);
if(!fCreateGeantinoMaps){
- MyEventAction* evtact = new MyEventAction();
- SetUserAction(evtact);
- SetUserAction(new MyTrackingAction(evtact));
- SetUserAction(new MySteppingAction(evtact));
+ MyEventAction* evtAct = new MyEventAction();
+ MyTrackingAction* trAct = new MyTrackingAction(evtAct);
+ MySteppingAction* stpAct = new MySteppingAction(evtAct);
+ SetUserAction(evtAct);
+ SetUserAction(trAct);
+ SetUserAction(stpAct);
+ runact->SetTrackingAction(trAct);
+ runact->SetSteppingAction(stpAct);
}
#if G4VERSION_NUMBER>=1040
diff --git a/FullSimLight/src/MyEventAction.cc b/FullSimLight/src/MyEventAction.cc
index 275095007963d8c508f3a11ba23cec130b6ccc46..4ba76fbf24b6aa2fc3befe795197c06ba335f831 100644
--- a/FullSimLight/src/MyEventAction.cc
+++ b/FullSimLight/src/MyEventAction.cc
@@ -22,8 +22,11 @@
#include "MyPrimaryGeneratorAction.hh"
-MyEventAction::MyEventAction() : G4UserEventAction() {
+MyEventAction::MyEventAction() : G4UserEventAction(), fIsSpecialScoring(false) {
fEventData.Clear();
+ if (fIsSpecialScoring) {
+ fEventDataSpecialRegion.Clear();
+ }
}
@@ -32,18 +35,24 @@ MyEventAction::~MyEventAction() { }
void MyEventAction::BeginOfEventAction(const G4Event*) {
fEventData.Clear();
+ if (fIsSpecialScoring) {
+ fEventDataSpecialRegion.Clear();
+ }
}
void MyEventAction::EndOfEventAction(const G4Event*) {
//get the Run and add the data collected during the simulation of the event that has been completed
MyRun* run = static_cast<MyRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
- run->FillPerEvent(fEventData);
+ run->FillPerEvent(fEventData, false);
+ if (fIsSpecialScoring) {
+ run->FillPerEvent(fEventDataSpecialRegion, true);
+ }
}
// called from the SteppingAction
-void MyEventAction::AddData(G4double edep, G4double length, G4bool ischarged) {
+void MyEventAction::AddData(G4double edep, G4double length, G4bool ischarged, G4bool isspecial) {
fEventData.fEdep += edep;
if (ischarged) {
fEventData.fTrackLCh += length;
@@ -52,17 +61,40 @@ void MyEventAction::AddData(G4double edep, G4double length, G4bool ischarged) {
fEventData.fTrackLNe += length;
fEventData.fNeutralStep += 1.;
}
+ // return if this step was not done in the special region
+ if (!isspecial) {
+ return;
+ }
+ fEventDataSpecialRegion.fEdep += edep;
+ if (ischarged) {
+ fEventDataSpecialRegion.fTrackLCh += length;
+ fEventDataSpecialRegion.fChargedStep += 1.;
+ } else {
+ fEventDataSpecialRegion.fTrackLNe += length;
+ fEventDataSpecialRegion.fNeutralStep += 1.;
+ }
}
// called from the TrackingAction
-void MyEventAction::AddSecondaryTrack(const G4Track* track) {
- const G4ParticleDefinition* pdf = track->GetDefinition();
- if (pdf==G4Gamma::Gamma()) {
- fEventData.fNGamma += 1.;
- } else if (pdf==G4Electron::Electron()) {
- fEventData.fNElec += 1.;
- } else if (pdf==G4Positron::Positron()) {
- fEventData.fNPosit += 1.;
- }
+void MyEventAction::AddSecondaryTrack(const G4Track* track, G4bool isspecial) {
+ const G4ParticleDefinition* pdf = track->GetDefinition();
+ if (pdf==G4Gamma::Gamma()) {
+ fEventData.fNGamma += 1.;
+ } else if (pdf==G4Electron::Electron()) {
+ fEventData.fNElec += 1.;
+ } else if (pdf==G4Positron::Positron()) {
+ fEventData.fNPosit += 1.;
+ }
+ // return if this step was not done in the special region
+ if (!isspecial) {
+ return;
+ }
+ if (pdf==G4Gamma::Gamma()) {
+ fEventDataSpecialRegion.fNGamma += 1.;
+ } else if (pdf==G4Electron::Electron()) {
+ fEventDataSpecialRegion.fNElec += 1.;
+ } else if (pdf==G4Positron::Positron()) {
+ fEventDataSpecialRegion.fNPosit += 1.;
+ }
}
diff --git a/FullSimLight/src/MyGVPhysicsList.cc b/FullSimLight/src/MyGVPhysicsList.cc
deleted file mode 100644
index ae8db553fe90fe4dc88e240dd2dd111698c1d318..0000000000000000000000000000000000000000
--- a/FullSimLight/src/MyGVPhysicsList.cc
+++ /dev/null
@@ -1,120 +0,0 @@
-
-#include "MyGVPhysicsList.hh"
-#include "G4Version.hh"
-#include "G4UnitsTable.hh"
-#include "G4SystemOfUnits.hh"
-
-#include "G4ParticleDefinition.hh"
-#include "G4Electron.hh"
-#include "G4Positron.hh"
-#include "G4Gamma.hh"
-
-#include "G4LossTableManager.hh"
-#include "G4ProcessManager.hh"
-#include "G4PhysicsListHelper.hh"
-
-#include "G4ComptonScattering.hh"
-#include "G4GammaConversion.hh"
-#include "G4PhotoElectricEffect.hh"
-#include "G4LivermorePhotoElectricModel.hh"
-//#include "G4RayleighScattering.hh"
-
-#include "G4eMultipleScattering.hh"
-#include "G4GoudsmitSaundersonMscModel.hh"
-#include "G4eIonisation.hh"
-#include "G4eBremsstrahlung.hh"
-#include "G4eplusAnnihilation.hh"
-
-#include "G4EmParameters.hh"
-#include "G4MscStepLimitType.hh"
-
-
-MyGVPhysicsList::MyGVPhysicsList() : G4VUserPhysicsList() {
- SetDefaultCutValue(1.0);
- SetVerboseLevel(0);
-}
-
-
-MyGVPhysicsList::~MyGVPhysicsList() {}
-
-
-void MyGVPhysicsList::ConstructParticle() {
- G4Electron::ElectronDefinition();
- G4Positron::PositronDefinition();
- G4Gamma::GammaDefinition();
-}
-
-
-void MyGVPhysicsList::ConstructProcess() {
- // Transportation
- AddTransportation();
- // EM physics
- BuildEMPhysics();
-}
-
-
-void MyGVPhysicsList::BuildEMPhysics() {
- G4EmParameters* param = G4EmParameters::Instance();
-#if G4VERSION_NUMBER>=1040
- param->SetDefaults();
-#endif
- param->SetVerbose(1);
- // inactivate energy loss fluctuations
- param->SetLossFluctuations(false);
- // inactivate to use cuts as final range
- param->SetUseCutAsFinalRange(false);
- //
- // MSC options:
- param->SetMscStepLimitType(fUseSafety);
- param->SetMscSkin(3);
- param->SetMscRangeFactor(0.06);
- G4LossTableManager::Instance();
- //
- // Add standard EM physics processes to e-/e+ and gamma that GeantV has
- G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
-#if G4VERSION_NUMBER<1040
- auto aParticleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
-#else
- auto aParticleIterator = GetParticleIterator();
-#endif
- aParticleIterator->reset();
- while((*aParticleIterator)()) {
- G4ParticleDefinition* particle = aParticleIterator->value();
- G4String particleName = particle->GetParticleName();
- if (particleName=="gamma") {
-// ph->RegisterProcess(new G4PhotoElectricEffect, particle);
- ph->RegisterProcess(new G4ComptonScattering(), particle);
- ph->RegisterProcess(new G4GammaConversion(), particle);
- G4double LivermoreLowEnergyLimit = 1.*eV;
- G4double LivermoreHighEnergyLimit = 1.*TeV;
- G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
- G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel = new G4LivermorePhotoElectricModel();
- theLivermorePhotoElectricModel->SetLowEnergyLimit(LivermoreLowEnergyLimit);
- theLivermorePhotoElectricModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
- thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
- ph->RegisterProcess(thePhotoElectricEffect, particle);
- } else if (particleName =="e-") {
-// ph->RegisterProcess(new G4eMultipleScattering(), particle);
- G4eMultipleScattering* msc = new G4eMultipleScattering;
- G4GoudsmitSaundersonMscModel* msc1 = new G4GoudsmitSaundersonMscModel();
- msc->AddEmModel(0, msc1);
- ph->RegisterProcess(msc,particle);
- //
- G4eIonisation* eIoni = new G4eIonisation();
- ph->RegisterProcess(eIoni, particle);
- ph->RegisterProcess(new G4eBremsstrahlung(), particle);
- } else if (particleName=="e+") {
-// ph->RegisterProcess(new G4eMultipleScattering(), particle);
- G4eMultipleScattering* msc = new G4eMultipleScattering;
- G4GoudsmitSaundersonMscModel* msc1 = new G4GoudsmitSaundersonMscModel();
- msc->AddEmModel(0, msc1);
- ph->RegisterProcess(msc,particle);
- //
- G4eIonisation* eIoni = new G4eIonisation();
- ph->RegisterProcess(eIoni, particle);
- ph->RegisterProcess(new G4eBremsstrahlung(), particle);
- //
- ph->RegisterProcess(new G4eplusAnnihilation(), particle);
- }
- }
-}
diff --git a/FullSimLight/src/MyRun.cc b/FullSimLight/src/MyRun.cc
index 15a6d33e3765c5a38bb4dee08383a27230da90c9..2b68006013411ad8ffa5f9b2913d4c6b2b5d3655 100644
--- a/FullSimLight/src/MyRun.cc
+++ b/FullSimLight/src/MyRun.cc
@@ -3,33 +3,44 @@
#include "MyEventData.hh"
#include "G4SystemOfUnits.hh"
+#include "G4Region.hh"
#include <iomanip>
-MyRun::MyRun() : G4Run() {
+MyRun::MyRun() : G4Run(), fScoringRegion(nullptr) {
fRunData.Clear();
+ fRunDataSpecialRegion.Clear();
}
MyRun::~MyRun() {}
-void MyRun::FillPerEvent(const MyEventData& data) {
- fRunData.fEdep += data.fEdep; fRunData.fEdep2 += data.fEdep*data.fEdep;
- fRunData.fTrackLCh += data.fTrackLCh; fRunData.fTrackLCh2 += data.fTrackLCh*data.fTrackLCh;
- fRunData.fTrackLNe += data.fTrackLNe; fRunData.fTrackLNe2 += data.fTrackLNe*data.fTrackLNe;
- fRunData.fChargedStep += data.fChargedStep; fRunData.fChargedStep2 += data.fChargedStep*data.fChargedStep;
- fRunData.fNeutralStep += data.fNeutralStep; fRunData.fNeutralStep2 += data.fNeutralStep*data.fNeutralStep;
- fRunData.fNGamma += data.fNGamma; fRunData.fNGamma2 += data.fNGamma*data.fNGamma;
- fRunData.fNElec += data.fNElec; fRunData.fNElec2 += data.fNElec*data.fNElec;
- fRunData.fNPosit += data.fNPosit; fRunData.fNPosit2 += data.fNPosit*data.fNPosit;
+void MyRun::FillPerEvent(const MyEventData& data, G4bool isspecial) {
+ if (!isspecial) {
+ AddData(data, fRunData);
+ } else {
+ AddData(data, fRunDataSpecialRegion);
+ }
}
+void MyRun::AddData(const MyEventData& data, MyRunData& runData) {
+ runData.fEdep += data.fEdep; runData.fEdep2 += data.fEdep*data.fEdep;
+ runData.fTrackLCh += data.fTrackLCh; runData.fTrackLCh2 += data.fTrackLCh*data.fTrackLCh;
+ runData.fTrackLNe += data.fTrackLNe; runData.fTrackLNe2 += data.fTrackLNe*data.fTrackLNe;
+ runData.fChargedStep += data.fChargedStep; runData.fChargedStep2 += data.fChargedStep*data.fChargedStep;
+ runData.fNeutralStep += data.fNeutralStep; runData.fNeutralStep2 += data.fNeutralStep*data.fNeutralStep;
+ runData.fNGamma += data.fNGamma; runData.fNGamma2 += data.fNGamma*data.fNGamma;
+ runData.fNElec += data.fNElec; runData.fNElec2 += data.fNElec*data.fNElec;
+ runData.fNPosit += data.fNPosit; runData.fNPosit2 += data.fNPosit*data.fNPosit;
+}
+
void MyRun::Merge(const G4Run* run) {
const MyRun* localRun = static_cast<const MyRun*>(run);
if (localRun) {
- fRunData += localRun->GetRunData();
+ fRunData += localRun->GetRunData(false);
+ fRunDataSpecialRegion += localRun->GetRunData(true);
}
G4Run::Merge(run);
}
@@ -53,22 +64,37 @@ void MyRun::EndOfRun() {
}
//compute and print statistic
//
- const G4double meanEdep = fRunData.fEdep*norm;
- const G4double rmsEdep = std::sqrt(std::abs(fRunData.fEdep2*norm-meanEdep*meanEdep));
- const G4double meanLCh = fRunData.fTrackLCh*norm;
- const G4double rmsLCh = std::sqrt(std::abs(fRunData.fTrackLCh2*norm-meanLCh*meanLCh));
- const G4double meanLNe = fRunData.fTrackLNe*norm;
- const G4double rmsLNe = std::sqrt(std::abs(fRunData.fTrackLNe2*norm-meanLNe*meanLNe));
- const G4double meanStpCh = fRunData.fChargedStep*norm;
- const G4double rmsStpCh = std::sqrt(std::abs(fRunData.fChargedStep2*norm-meanStpCh*meanStpCh));
- const G4double meanStpNe = fRunData.fNeutralStep*norm;
- const G4double rmsStpNe = std::sqrt(std::abs(fRunData.fNeutralStep2*norm-meanStpNe*meanStpNe));
- const G4double meanNGam = fRunData.fNGamma*norm;
- const G4double rmsNGam = std::sqrt(std::abs(fRunData.fNGamma2*norm-meanNGam*meanNGam));
- const G4double meanNElec = fRunData.fNElec*norm;
- const G4double rmsNElec = std::sqrt(std::abs(fRunData.fNElec2*norm-meanNElec*meanNElec));
- const G4double meanNPos = fRunData.fNPosit*norm;
- const G4double rmsNPos = std::sqrt(std::abs(fRunData.fNPosit2*norm-meanNPos*meanNPos));
+ PrintEndOfRunStat(fRunData, norm);
+ if (fScoringRegion != nullptr) {
+ G4cout<< "\n --- In the Special Detector Region: " << fScoringRegion->GetName() << " \n" << G4endl;
+ PrintEndOfRunStat(fRunDataSpecialRegion, norm);
+ }
+ G4cout<< " ......................................................................................... \n" << G4endl;
+
+ G4cout<< " \n ======================================================================================== \n" << G4endl;
+
+ G4cout.setf(mode,std::ios::floatfield);
+ G4cout.precision(prec);
+}
+
+void MyRun::PrintEndOfRunStat(MyRunData& runData, G4double norm) {
+ //compute and print statistic
+ const G4double meanEdep = runData.fEdep*norm;
+ const G4double rmsEdep = std::sqrt(std::abs(runData.fEdep2*norm-meanEdep*meanEdep));
+ const G4double meanLCh = runData.fTrackLCh*norm;
+ const G4double rmsLCh = std::sqrt(std::abs(runData.fTrackLCh2*norm-meanLCh*meanLCh));
+ const G4double meanLNe = runData.fTrackLNe*norm;
+ const G4double rmsLNe = std::sqrt(std::abs(runData.fTrackLNe2*norm-meanLNe*meanLNe));
+ const G4double meanStpCh = runData.fChargedStep*norm;
+ const G4double rmsStpCh = std::sqrt(std::abs(runData.fChargedStep2*norm-meanStpCh*meanStpCh));
+ const G4double meanStpNe = runData.fNeutralStep*norm;
+ const G4double rmsStpNe = std::sqrt(std::abs(runData.fNeutralStep2*norm-meanStpNe*meanStpNe));
+ const G4double meanNGam = runData.fNGamma*norm;
+ const G4double rmsNGam = std::sqrt(std::abs(runData.fNGamma2*norm-meanNGam*meanNGam));
+ const G4double meanNElec = runData.fNElec*norm;
+ const G4double rmsNElec = std::sqrt(std::abs(runData.fNElec2*norm-meanNElec*meanNElec));
+ const G4double meanNPos = runData.fNPosit*norm;
+ const G4double rmsNPos = std::sqrt(std::abs(runData.fNPosit2*norm-meanNPos*meanNPos));
G4cout<< " Mean energy deposit per event = " << meanEdep/GeV << " +- " << rmsEdep/GeV << " [GeV]" <<G4endl;
G4cout<< G4endl;
@@ -83,9 +109,4 @@ void MyRun::EndOfRun() {
<< " Electrons = " << meanNElec << " +- " << rmsNElec << G4endl
<< " Positrons = " << meanNPos << " +- " << rmsNPos << G4endl;
G4cout<< " ......................................................................................... \n" << G4endl;
-
- G4cout<< " \n ======================================================================================== \n" << G4endl;
-
- G4cout.setf(mode,std::ios::floatfield);
- G4cout.precision(prec);
}
diff --git a/FullSimLight/src/MyRunAction.cc b/FullSimLight/src/MyRunAction.cc
index 6b624953fc28236828794dbef5418d68a954d5d9..38f5ccc9de0cc065f2e78dcfcb5c121a30cce63b 100644
--- a/FullSimLight/src/MyRunAction.cc
+++ b/FullSimLight/src/MyRunAction.cc
@@ -11,6 +11,8 @@
#include "globals.hh"
#include "MyRun.hh"
+#include "MySteppingAction.hh"
+#include "MyTrackingAction.hh"
#include "G4ProductionCutsTable.hh"
@@ -20,9 +22,9 @@
G4AnalysisManager* MyRunAction::fMasterAnalysisManager = nullptr;
MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename)
-: G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino),
- fRun(nullptr), fTimer(nullptr), fGeantinoMapsFilename(geantinoMapFilename),
- fPythiaConfig("") { }
+: G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino), fRun(nullptr), fTimer(nullptr),
+ fSteppingAction(nullptr), fTrackingAction(nullptr), fGeantinoMapsFilename(geantinoMapFilename),
+ fPythiaConfig(""), fSpecialScoringRegionName("") { }
MyRunAction::~MyRunAction() {
if(fIsGeantino)
@@ -103,6 +105,23 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
}
#endif
+
+ // set the special scoring region in the stepping and tracking actions
+ if (fRun && fSpecialScoringRegionName!="") {
+ std::vector<G4Region*>* theRegionVector = G4RegionStore::GetInstance();
+ for (std::size_t ir=0, nr=theRegionVector->size(); ir<nr; ++ir) {
+ G4Region* reg = (*theRegionVector)[ir];
+ if (reg->GetName() == fSpecialScoringRegionName) {
+ if (isMaster) {
+ G4cout << " --- Special scoring will be done in Region = " << fSpecialScoringRegionName << G4endl;
+ }
+ if (fSteppingAction) { fSteppingAction->SetScoringRegion(reg); }
+ if (fTrackingAction) { fTrackingAction->SetScoringRegion(reg); }
+ fRun->SetSpecialScoringRegion(reg);
+ }
+ }
+ }
+
if (isMaster) {
//G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
diff --git a/FullSimLight/src/MySteppingAction.cc b/FullSimLight/src/MySteppingAction.cc
index e33e45022bbd9d45932662dfe9c16ba7a9426f09..35b8ae8520b4a693bc55baaa308a770ff45fa26e 100644
--- a/FullSimLight/src/MySteppingAction.cc
+++ b/FullSimLight/src/MySteppingAction.cc
@@ -4,17 +4,26 @@
#include "MyEventAction.hh"
#include "G4Step.hh"
#include "G4ParticleDefinition.hh"
+#include "G4Region.hh"
-MySteppingAction::MySteppingAction(MyEventAction* evtact) : G4UserSteppingAction(), fEventAction(evtact) {}
+MySteppingAction::MySteppingAction(MyEventAction* evtact)
+: G4UserSteppingAction(), fEventAction(evtact), fScoringRegion(nullptr) {}
MySteppingAction::~MySteppingAction() {}
+void MySteppingAction::SetScoringRegion(G4Region* reg) {
+ fScoringRegion = reg;
+ fEventAction->SetIsSpecialScoringRegion(reg!=nullptr);
+}
+
void MySteppingAction::UserSteppingAction(const G4Step* theStep) {
- G4double edep = theStep->GetTotalEnergyDeposit();
- G4double stepl = theStep->GetStepLength();
- G4bool isCharged = (theStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.);
- fEventAction->AddData(edep, stepl, isCharged);
+ // scoring in the spcial given region if any
+ const G4bool isSpecialScoring = (fScoringRegion!=nullptr && fScoringRegion==theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion());
+ const G4double edep = theStep->GetTotalEnergyDeposit();
+ const G4double stepl = theStep->GetStepLength();
+ const G4bool isCharged = (theStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.);
+ fEventAction->AddData(edep, stepl, isCharged, isSpecialScoring);
}
diff --git a/FullSimLight/src/MyTrackingAction.cc b/FullSimLight/src/MyTrackingAction.cc
index 9e3c591ac9354281061086a61afff2d9dd1f6e62..73dd6ec70c7ad51099d80ac0c369b8ccee4d6458 100644
--- a/FullSimLight/src/MyTrackingAction.cc
+++ b/FullSimLight/src/MyTrackingAction.cc
@@ -2,18 +2,26 @@
#include "MyTrackingAction.hh"
#include "G4Track.hh"
+#include "G4Region.hh"
#include "MyEventAction.hh"
-MyTrackingAction::MyTrackingAction(MyEventAction* evtact) : G4UserTrackingAction(), fEventAction(evtact) {}
+MyTrackingAction::MyTrackingAction(MyEventAction* evtact)
+: G4UserTrackingAction(), fEventAction(evtact), fScoringRegion(nullptr) {}
MyTrackingAction::~MyTrackingAction() {}
+void MyTrackingAction::SetScoringRegion(G4Region* reg) {
+ fScoringRegion = reg;
+ fEventAction->SetIsSpecialScoringRegion(reg!=nullptr);
+}
+
void MyTrackingAction::PreUserTrackingAction(const G4Track* theTrack) {
// count secondaries
+ const G4bool isSpecialScoring = (fScoringRegion!=nullptr && fScoringRegion==theTrack->GetVolume()->GetLogicalVolume()->GetRegion());
if (theTrack->GetParentID()!=0) {
- fEventAction->AddSecondaryTrack(theTrack);
+ fEventAction->AddSecondaryTrack(theTrack, isSpecialScoring);
}
}
diff --git a/FullSimLight/src/StandardEmWithWoodcock.cc b/FullSimLight/src/StandardEmWithWoodcock.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7ee881482f69faa24a2c9f664b315fe520b0c586
--- /dev/null
+++ b/FullSimLight/src/StandardEmWithWoodcock.cc
@@ -0,0 +1,377 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+//---------------------------------------------------------------------------
+//
+// ClassName: StandardEmWithWoodcock
+//
+// Author: V.Ivanchenko 09.11.2005
+//
+// Modified:
+//
+// NOTE: M. Novak. 27.11.2021
+//
+// This EM physics constructor is almost identical to the EM standard (Em-opt0)
+// physics constructor with the possibility of utilising Woodcock-tracking in a
+// given detector region. The name of the region needs to be provided to the
+// `WoodcockProcess` at construction. The `WoodckockProcess` falls into the base
+// Gamma-General process (a local version of the `G4GammaGeneralProcess`) when:
+// - the given detector region cannot be found in the detector
+// - the gamma energy is below a given threshold (by default 0; can be set at
+// construction time)
+// As in the case of the base Gamma-General process, all the processes neeeds to
+// be assigned to the `G4Gamma` particle through the `WoodckockProcess` (and not
+// directly!).
+//----------------------------------------------------------------------------
+//
+
+
+#include "StandardEmWithWoodcock.hh"
+
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4EmParameters.hh"
+#include "G4LossTableManager.hh"
+
+#include "G4ComptonScattering.hh"
+#include "G4GammaConversion.hh"
+#include "G4PhotoElectricEffect.hh"
+#include "G4RayleighScattering.hh"
+#include "G4LivermorePhotoElectricModel.hh"
+
+#include "G4eMultipleScattering.hh"
+#include "G4MuMultipleScattering.hh"
+#include "G4hMultipleScattering.hh"
+#include "G4CoulombScattering.hh"
+#include "G4eCoulombScatteringModel.hh"
+#include "G4WentzelVIModel.hh"
+#include "G4UrbanMscModel.hh"
+
+#include "G4MuBremsstrahlungModel.hh"
+#include "G4MuPairProductionModel.hh"
+#include "G4hBremsstrahlungModel.hh"
+#include "G4hPairProductionModel.hh"
+
+#include "G4eIonisation.hh"
+#include "G4eBremsstrahlung.hh"
+#include "G4eplusAnnihilation.hh"
+#include "G4UAtomicDeexcitation.hh"
+
+#include "G4MuIonisation.hh"
+#include "G4MuBremsstrahlung.hh"
+#include "G4MuPairProduction.hh"
+#include "G4hBremsstrahlung.hh"
+#include "G4hPairProduction.hh"
+
+#include "G4hIonisation.hh"
+#include "G4ionIonisation.hh"
+#include "G4alphaIonisation.hh"
+
+#include "G4ParticleTable.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+#include "G4PionPlus.hh"
+#include "G4PionMinus.hh"
+#include "G4KaonPlus.hh"
+#include "G4KaonMinus.hh"
+#include "G4Proton.hh"
+#include "G4AntiProton.hh"
+#include "G4Deuteron.hh"
+#include "G4Triton.hh"
+#include "G4He3.hh"
+#include "G4Alpha.hh"
+#include "G4GenericIon.hh"
+
+#include "G4PhysicsListHelper.hh"
+#include "G4BuilderType.hh"
+#include "G4EmModelActivator.hh"
+
+// use the local version of the G4GammaGeneralProcess till Geant4-11.0
+#include "GammaGeneralProcess.hh"
+#include "WoodcockProcess.hh"
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+StandardEmWithWoodcock::StandardEmWithWoodcock(G4int ver, const G4String& name)
+ : G4VPhysicsConstructor(name), verbose(ver), fWDCKRegionName(""), fWDCKLowEnergyThreshold(0.0)
+{
+ SetVerboseLevel(ver);
+ G4EmParameters* param = G4EmParameters::Instance();
+ param->SetDefaults();
+ param->SetVerbose(ver);
+ // activate gamma-general process in order to add all (even non EM) proesses
+ param->SetGeneralProcessActive(true);
+ SetPhysicsType(bElectromagnetic);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+StandardEmWithWoodcock::~StandardEmWithWoodcock()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void StandardEmWithWoodcock::ConstructParticle()
+{
+ // gamma
+ G4Gamma::Gamma();
+ // leptons
+ G4Electron::Electron();
+ G4Positron::Positron();
+ G4MuonPlus::MuonPlus();
+ G4MuonMinus::MuonMinus();
+ // mesons
+ G4PionPlus::PionPlus();
+ G4PionMinus::PionMinus();
+ G4KaonPlus::KaonPlus();
+ G4KaonMinus::KaonMinus();
+ // barions
+ G4Proton::Proton();
+ G4AntiProton::AntiProton();
+ // ions
+ G4Deuteron::Deuteron();
+ G4Triton::Triton();
+ G4He3::He3();
+ G4Alpha::Alpha();
+ G4GenericIon::GenericIonDefinition();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//
+// NOTE: M. Novak
+// Using old style EM process construction that good for G4.10.6p03 as well
+void StandardEmWithWoodcock::ConstructProcess()
+{
+ if(verbose > 1) {
+ G4cout << "### " << GetPhysicsName() << " Construct Processes " << G4endl;
+ }
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+ G4LossTableManager* man = G4LossTableManager::Instance();
+
+ // muon & hadron bremsstrahlung and pair production
+ G4MuBremsstrahlung* mub = new G4MuBremsstrahlung();
+ G4MuPairProduction* mup = new G4MuPairProduction();
+ G4hBremsstrahlung* pib = new G4hBremsstrahlung();
+ G4hPairProduction* pip = new G4hPairProduction();
+ G4hBremsstrahlung* kb = new G4hBremsstrahlung();
+ G4hPairProduction* kp = new G4hPairProduction();
+ G4hBremsstrahlung* pb = new G4hBremsstrahlung();
+ G4hPairProduction* pp = new G4hPairProduction();
+
+ // muon & hadron multiple scattering
+ G4MuMultipleScattering* mumsc = new G4MuMultipleScattering();
+ mumsc->SetEmModel(new G4WentzelVIModel());
+ G4CoulombScattering* muss = new G4CoulombScattering();
+
+ G4hMultipleScattering* pimsc = new G4hMultipleScattering();
+ pimsc->SetEmModel(new G4WentzelVIModel());
+ G4CoulombScattering* piss = new G4CoulombScattering();
+
+ G4hMultipleScattering* kmsc = new G4hMultipleScattering();
+ kmsc->SetEmModel(new G4WentzelVIModel());
+ G4CoulombScattering* kss = new G4CoulombScattering();
+
+ G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
+
+ // high energy limit for e+- scattering models
+ G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
+
+ // Add standard EM Processes
+ G4ParticleTable* table = G4ParticleTable::GetParticleTable();
+ for(const auto& particleName : partList.PartNames()) {
+ G4ParticleDefinition* particle = table->FindParticle(particleName);
+ if (!particle) { continue; }
+ if (particleName == "gamma") {
+ //
+ // Create a Woodcock process, add all EM Gamma processes (in the same way as in
+ // case of the G4GammaGeneralProcess). Note that further processes (extra and
+ // non-EM like gamma-nuclear) might be added automatically in later physics
+ // constructors just as in case of the Gamma-Generel process. This is done
+ // automatically when the general-process is activated, e.g. by invoking
+ // `G4EmParameters::Instance()->SetGeneralProcessActive(true)` or through the
+ // `/process/em/UseGeneralProcess true` UI command.
+ //
+ //
+ // Woodcock tracking will be required inside the `theWDCKregionName`
+ // and only above `theWDCKLowEnergyThreshold` gamma energies.
+ // NOTE: the `WoodcockProcess` is equivalent to the `G4GammaGeneralProcess`
+ // when no detector region can be found with the given name.
+ WoodcockProcess* theWDCKProcess = new WoodcockProcess(fWDCKRegionName, fWDCKLowEnergyThreshold);
+ //
+ G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
+ pee->SetEmModel(new G4LivermorePhotoElectricModel());
+ theWDCKProcess->AddEmProcess(pee);
+ theWDCKProcess->AddEmProcess(new G4ComptonScattering);
+ theWDCKProcess->AddEmProcess(new G4GammaConversion);
+ theWDCKProcess->AddEmProcess(new G4RayleighScattering);
+ //
+ G4LossTableManager::Instance()->SetGammaGeneralProcess(theWDCKProcess);
+ ph->RegisterProcess(theWDCKProcess, particle);
+
+ } else if (particleName == "e-") {
+
+ G4eMultipleScattering* msc = new G4eMultipleScattering;
+ G4UrbanMscModel* msc1 = new G4UrbanMscModel();
+ G4WentzelVIModel* msc2 = new G4WentzelVIModel();
+ msc1->SetHighEnergyLimit(highEnergyLimit);
+ msc2->SetLowEnergyLimit(highEnergyLimit);
+ msc->SetEmModel(msc1);
+ msc->SetEmModel(msc2);
+
+ G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
+ G4CoulombScattering* ss = new G4CoulombScattering();
+ ss->SetEmModel(ssm);
+ ss->SetMinKinEnergy(highEnergyLimit);
+ ssm->SetLowEnergyLimit(highEnergyLimit);
+ ssm->SetActivationLowEnergyLimit(highEnergyLimit);
+
+ ph->RegisterProcess(msc, particle);
+ ph->RegisterProcess(new G4eIonisation(), particle);
+ ph->RegisterProcess(new G4eBremsstrahlung(), particle);
+ ph->RegisterProcess(ss, particle);
+
+ } else if (particleName == "e+") {
+
+ G4eMultipleScattering* msc = new G4eMultipleScattering;
+ G4UrbanMscModel* msc1 = new G4UrbanMscModel();
+ G4WentzelVIModel* msc2 = new G4WentzelVIModel();
+ msc1->SetHighEnergyLimit(highEnergyLimit);
+ msc2->SetLowEnergyLimit(highEnergyLimit);
+ msc->SetEmModel(msc1);
+ msc->SetEmModel(msc2);
+
+ G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
+ G4CoulombScattering* ss = new G4CoulombScattering();
+ ss->SetEmModel(ssm);
+ ss->SetMinKinEnergy(highEnergyLimit);
+ ssm->SetLowEnergyLimit(highEnergyLimit);
+ ssm->SetActivationLowEnergyLimit(highEnergyLimit);
+
+ ph->RegisterProcess(msc, particle);
+ ph->RegisterProcess(new G4eIonisation(), particle);
+ ph->RegisterProcess(new G4eBremsstrahlung(), particle);
+ ph->RegisterProcess(new G4eplusAnnihilation(), particle);
+ ph->RegisterProcess(ss, particle);
+
+ } else if (particleName == "mu+" ||
+ particleName == "mu-" ) {
+
+ ph->RegisterProcess(mumsc, particle);
+ ph->RegisterProcess(new G4MuIonisation(), particle);
+ ph->RegisterProcess(mub, particle);
+ ph->RegisterProcess(mup, particle);
+ ph->RegisterProcess(muss, particle);
+
+ } else if (particleName == "alpha" ||
+ particleName == "He3") {
+
+ ph->RegisterProcess(new G4hMultipleScattering(), particle);
+ ph->RegisterProcess(new G4ionIonisation(), particle);
+
+ } else if (particleName == "GenericIon") {
+
+ ph->RegisterProcess(hmsc, particle);
+ ph->RegisterProcess(new G4ionIonisation(), particle);
+
+ } else if (particleName == "pi+" ||
+ particleName == "pi-" ) {
+
+ ph->RegisterProcess(pimsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(pib, particle);
+ ph->RegisterProcess(pip, particle);
+ ph->RegisterProcess(piss, particle);
+
+ } else if (particleName == "kaon+" ||
+ particleName == "kaon-" ) {
+
+ ph->RegisterProcess(kmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(kb, particle);
+ ph->RegisterProcess(kp, particle);
+ ph->RegisterProcess(kss, particle);
+
+ } else if (particleName == "proton" ||
+ particleName == "anti_proton") {
+
+ G4hMultipleScattering* pmsc = new G4hMultipleScattering();
+ pmsc->SetEmModel(new G4WentzelVIModel());
+
+ ph->RegisterProcess(pmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(pb, particle);
+ ph->RegisterProcess(pp, particle);
+ ph->RegisterProcess(new G4CoulombScattering(), particle);
+
+ } else if (particleName == "B+" ||
+ particleName == "B-" ||
+ particleName == "D+" ||
+ particleName == "D-" ||
+ particleName == "Ds+" ||
+ particleName == "Ds-" ||
+ particleName == "anti_He3" ||
+ particleName == "anti_alpha" ||
+ particleName == "anti_deuteron" ||
+ particleName == "anti_lambda_c+" ||
+ particleName == "anti_omega-" ||
+ particleName == "anti_sigma_c+" ||
+ particleName == "anti_sigma_c++" ||
+ particleName == "anti_sigma+" ||
+ particleName == "anti_sigma-" ||
+ particleName == "anti_triton" ||
+ particleName == "anti_xi_c+" ||
+ particleName == "anti_xi-" ||
+ particleName == "deuteron" ||
+ particleName == "lambda_c+" ||
+ particleName == "omega-" ||
+ particleName == "sigma_c+" ||
+ particleName == "sigma_c++" ||
+ particleName == "sigma+" ||
+ particleName == "sigma-" ||
+ particleName == "tau+" ||
+ particleName == "tau-" ||
+ particleName == "triton" ||
+ particleName == "xi_c+" ||
+ particleName == "xi-" ) {
+ //
+ ph->RegisterProcess(hmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ }
+ }
+
+ // Deexcitation
+ //
+ man->SetAtomDeexcitation(new G4UAtomicDeexcitation());
+
+ G4EmModelActivator mact(GetPhysicsName());
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/FullSimLight/src/WoodcockProcess.cc b/FullSimLight/src/WoodcockProcess.cc
new file mode 100644
index 0000000000000000000000000000000000000000..00b41834f1f78ec57d3a1d0b609dc1eb8dcfc9b2
--- /dev/null
+++ b/FullSimLight/src/WoodcockProcess.cc
@@ -0,0 +1,370 @@
+#include "WoodcockProcess.hh"
+
+#include "G4Region.hh"
+#include "G4VSolid.hh"
+#include "G4RegionStore.hh"
+
+
+#include "G4TransportationManager.hh"
+#include "G4Navigator.hh"
+#include "G4MaterialCutsCouple.hh"
+
+
+WoodcockProcess::WoodcockProcess(const G4String& wdckRegionName, G4double wdckLowELimit)
+: GammaGeneralProcess("Woodcock-GG"),
+ fWDCKRegionName(wdckRegionName),
+ fWDCKLowEnergyThreshold(wdckLowELimit),
+ fWDCKNumRootLVolume(0),
+ fWDCKRegion(nullptr),
+ fWDCKCouple(nullptr),
+ fOnWDCKRegionBoundary(false),
+ fInWDCKRegion(false),
+ fWDCKRootLVolumeIndx(-1),
+ fWDCKStepLength(0.0) {
+ preStepKinEnergy = -1.0;
+ preStepLambda = 0.0;
+ currentCouple = nullptr;
+}
+
+
+WoodcockProcess::~WoodcockProcess() {
+ // delete the local navigator used to locate the post-step point in WDCK tracking
+ for (std::size_t iv=0; iv<fWDCKNumRootLVolume; ++iv) {
+ if (fWDCKDataPerRootLVolume[iv].fNavigator) {
+ delete fWDCKDataPerRootLVolume[iv].fNavigator;
+ }
+ }
+}
+
+
+void WoodcockProcess::StartTracking(G4Track* track) {
+ // set initial values of some flags
+ fOnWDCKRegionBoundary = false;
+ fWDCKRootLVolumeIndx = -1;
+ // call the base class method
+ GammaGeneralProcess::StartTracking(track);
+}
+
+
+void WoodcockProcess::PreparePhysicsTable(const G4ParticleDefinition& part) {
+ GammaGeneralProcess::PreparePhysicsTable(part);
+}
+
+// Init the WDCK tracking if the detector region, specified by the user giving
+// its name, was found. Only partial initialisation is done here: finding the
+// heaviest material of the region and setting the size of containers that stores
+// the individual root logical volume related information. The other part of the
+// init is done when the fist gamma track enters to the given volumes (this is
+// because the geometry must be closed for getting further require information).
+void WoodcockProcess::BuildPhysicsTable(const G4ParticleDefinition& part) {
+ // try to find the WDCK region in the store by its name given by the user
+ fWDCKRegion = G4RegionStore::GetInstance()->GetRegion(fWDCKRegionName);
+ if (fWDCKRegion != nullptr) {
+ fWDCKNumRootLVolume = fWDCKRegion->GetNumberOfRootVolumes();
+ G4AffineTransform aTransform;
+ fWDCKDataPerRootLVolume.resize(fWDCKNumRootLVolume, {aTransform, nullptr, nullptr, nullptr, true});
+ // build the root logical volume ID --> vector index (in the fWDCKDataPerRootLVolume) map
+ std::vector<G4LogicalVolume*>::const_iterator itrLV = fWDCKRegion->GetRootLogicalVolumeIterator();
+ for (std::size_t ilv = 0; ilv<fWDCKNumRootLVolume; ++ilv) {
+ fWDCKRootLVolIDToIndxMap[(*itrLV)->GetInstanceID()] = ilv;
+ fWDCKDataPerRootLVolume[ilv].fLogicalVolume = (*itrLV);
+ ++itrLV;
+ }
+ // loop over the materials of this region and get the couple of the heaviest
+ std::vector<G4Material*>::const_iterator itr = fWDCKRegion->GetMaterialIterator();
+ const G4int nmat = fWDCKRegion->GetNumberOfMaterials();
+ G4double maxDensity =0.0;
+ for (G4int im=0; im<nmat; ++im) {
+ const G4double density = (*itr)->GetDensity();
+ if (density>maxDensity) {
+ maxDensity = density;
+ fWDCKCouple = fWDCKRegion->FindCouple((*itr));
+ }
+ ++itr;
+ }
+ }
+ // call the underlying Gamma-general process `BuildPhysicsTable`
+ GammaGeneralProcess::BuildPhysicsTable(part);
+}
+
+
+G4double
+WoodcockProcess::PostStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4ForceCondition* condition) {
+ // check if this step will be performed within the Woodcock region
+ *condition = G4ForceCondition::NotForced;
+ G4double stepLength = DBL_MAX;
+ fInWDCKRegion = false;
+ // NOTE: this is in order to get out from the Woodcock Region
+ //
+ // Since we make a dirty hack of the navigation/transportation,
+ // the GPIL and DoIt methods of transportation are not called properly
+ // and boundary crossing is not recognised even if we `put` the track on
+ // the boundary. So we just detect when it happenes based on this flag
+ // set in the previous step limit:
+ // - setting the `fOnWDCKRegionBoundary` flag to be true when reaching
+ // the WDCK top volume boundary (see well below)
+ // - shorten a tiny bit of the overall step length inside the WDCK top
+ // volume in order to arrive in the vicinity of the WDCK top volume
+ // bundary (at the end of the pevious step)
+ // - also note, that the DoIt method detects this case and there is no
+ // interaction in at the end of the previous step
+ // So now we are very close to the WDCK top volume boundary and just
+ // giving a step length that ensures that Transportation wins that makes
+ // sure the appropriate exit from the WDCK region.
+ // Therefore, at the very end of the WDCK tracking of a photon in the
+ // WDCK root volume, this top volume is left in 2 "transportation" steps:
+ // one fake that just brings close (< 1.0E-4) to the boundary and a real
+ // one that will make an appropriate boundary crossing step.
+ if (fOnWDCKRegionBoundary) {
+ fOnWDCKRegionBoundary = false;
+ // we are already reached the WDCK region boundary at the end of the previous
+ // step: just give a step length > 1.0E-4 and let the transportation win this
+ // time normaly (this will set everything properly)
+ fWDCKStepLength = 10.0;
+ fWDCKRootLVolumeIndx = -1;
+ return fWDCKStepLength;
+ }
+
+ // Check if WDCK tracking was required ,i.e. if the WDCK region was found at init.,
+ // and the current step is performed inside this region with a high enough energy.
+ if (fWDCKRegion !=nullptr && fWDCKRegion == track.GetVolume()->GetLogicalVolume()->GetRegion() && track.GetKineticEnergy()>fWDCKLowEnergyThreshold) {
+ fInWDCKRegion = true;
+ // Need to find the root logical volume of this branch: only in the very
+ // first step in this branch, i.e. when entering or just first step after
+ // created in this branch.) After this, fWDCKRootLVolumeIndx is set for sure.
+ if (fWDCKNumRootLVolume > 1) {
+ FindWDCKRootLogicalVolume(track);
+ } else {
+ fWDCKRootLVolumeIndx = 0;
+ }
+ // Compute the distance to WDCK region root volume boundary
+ // - get direction and starting position
+ const G4ThreeVector& r0 = track.GetPosition();
+ const G4ThreeVector& v0 = track.GetMomentumDirection();
+ // This method is called only once when the very firts gamma that makes its
+ // fist step inside this branch (i.e. below this root logical volume) of the
+ // WDKC region: initialize the local navigation and transformation for this
+ // root logical volume of the WDCK region.
+ if (fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fNeedsInit) {
+ InitWoodcockTracking(r0, v0);
+ fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fNeedsInit = false;
+ }
+ G4double distToBoundary = ComputeStepToWDCKRegionBoundary(r0, v0);
+ // Compute the total macroscopic cross section in the WDCK material
+ // (result will be in `preStepLambda`)
+ const G4double ekin = track.GetKineticEnergy();
+ const G4double lekin = track.GetDynamicParticle()->GetLogKineticEnergy();
+ ComputeTotalMacrCrossSection(fWDCKCouple, ekin, lekin);
+ const G4double invWDCKMXsec = preStepLambda > 0.0 ? 1.0/preStepLambda : DBL_MAX;
+ // Set initial values for the step length and the flag that indicates ending
+ // the WDCK root volume boundary on boundary
+ fWDCKStepLength = 0.0;
+ fOnWDCKRegionBoundary = false;
+ stepLength = 0.0;
+ G4bool doStop = false;
+ G4Navigator* localNav = fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fNavigator;
+ const G4Material* wdckMaterial = fWDCKCouple->GetMaterial();
+ // Untill we either interact or reach the boundary of the WDCK region:
+ while (true) {
+ // Compute the step length till the next interaction in the WDCK material
+ const G4double pstep = invWDCKMXsec < DBL_MAX ? -G4Log( G4UniformRand() )*invWDCKMXsec : DBL_MAX;
+ // Take the minimum of this and the distance to the WDCK root volume boundary
+ // Check if this step will end up on the WDCK region boundary
+ if (distToBoundary < pstep) {
+ stepLength += distToBoundary;
+ fOnWDCKRegionBoundary = true;
+ doStop = true;
+ } else {
+ // We will move the particle by a step length of `pstep` so we reduce the
+ // distance to boundary with teh same length.
+ stepLength += pstep;
+ distToBoundary -= pstep;
+ // Locate the actual post step point in order to get the real material.
+ // Use relative search (true) since we already located a point at the
+ // local navigator init inside the WDCK region and ignore direction (true)
+ const G4VPhysicalVolume* pVol = localNav->LocateGlobalPointAndSetup(r0+stepLength*v0, nullptr, true, true);
+ const G4LogicalVolume* lVol = pVol->GetLogicalVolume();
+ const G4Material* postStepMat = lVol->GetMaterial();
+ // Check if the real material of the post-step point is the WDCK one
+ if (wdckMaterial != postStepMat) {
+ // The post step point is not in the WDCK material: need to check if we interact
+ // Compute the total macroscopic cross section for the real, non-wdck material
+ // (the result will be in `preStepLambda`)
+ const G4MaterialCutsCouple* couple = lVol->GetMaterialCutsCouple();
+ ComputeTotalMacrCrossSection(couple, ekin, lekin);
+ // P(interact) = preStepLambda/wdckMXsec note: preStepLambda <= wdckMXsec
+ doStop = (preStepLambda*invWDCKMXsec > G4UniformRand()); // interact
+ } else {
+ // Interact since the post step point is in the wdck material which cross section was used
+ // Make sure that all members of the underlying gamma-general set properly so re-compute
+ // the partial macroscopic cross sections (if needed)
+ ComputeTotalMacrCrossSection(fWDCKCouple, ekin, lekin);
+ doStop = true;
+ }
+ }
+ if (doStop) {
+ // Reached the end, i.e. either interaction happens at the current post step point
+ // or it is locating at the WDCK reagion boundary (fOnWDCKRegionBoundary=true).
+ // Update the track position and return with zero step lenght (transportation/navigation)
+ // will correctly re-locate the post-step point of the track while we get back to the DoIt()
+ // BUT we do it by our self in order to avoid Warnings!
+ // NOTE: the material-cuts couple and the material that will be used at
+ // the post-step DoIt in case of interaction are those that are in
+ // `currentCouple` and `currentMaterial`, etc. and not taken from
+ // the track. Therefore, as long as we are inside the given WDCK
+ // root volume, we do not care much if Transporation relocates or
+ // not the track. If the gamma leaves the given WDCK root volume,
+ // we make sure that a normal transportation step will transverse
+ // the volume boundary.
+ G4Track& theTrack = const_cast<G4Track&>(track);
+ // shorten a bit this step length if we would reach the boundary:
+ // - WoodcockProcess limits the step but we see in the DoIt that actually its boundary
+ // - the next step will be real Transportation since Woodcock will propose a longer step
+ // - this will bring us properly to the next boundary in the next step
+ double sl = stepLength;
+ if (fOnWDCKRegionBoundary) {
+ sl = sl > 1.0E-4 ? sl-1.0E-4 : sl*0.999;
+ }
+ fWDCKStepLength = sl;
+ theTrack.SetPosition(track.GetPosition()+sl*track.GetMomentumDirection());
+ // relocate the moved G4Track in order to avoid Navigator warnings
+ G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->LocateGlobalPointAndUpdateTouchableHandle(track.GetPosition(), track.GetMomentumDirection(), const_cast<G4TouchableHandle&>(track.GetTouchableHandle()));
+// theTrack.GetStep()->GetPostStepPoint()->SetTouchableHandle(track.GetTouchableHandle());
+ //
+ return 0.0;
+ }
+ // If the WDCK region boundary has not been reached in this step then delta interaction
+ // happend so just keepKeep moving the post-step point toward to the WDCK region boundary otherwise by
+ // sampling and performing yet an other step
+ };
+ } else {
+ // This step will be performed OUTSIDE the WDK region: compute the step exactly as the
+ // the Gamma-general process PostStepGetPhysicalInteractionLength
+ return GammaGeneralProcess::PostStepGetPhysicalInteractionLength(track, previousStepSize, condition);
+ }
+}
+
+
+// Computes the path length to the WDCK root volume boundary.
+// NOTE: assumes that there is a single WDCK region with a single root volume.
+G4double
+WoodcockProcess::ComputeStepToWDCKRegionBoundary(const G4ThreeVector& pGlobalPoint, const G4ThreeVector& pDirection) {
+ const G4AffineTransform& transform = fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fAffineTransform;
+ const G4ThreeVector localPoint = transform.TransformPoint(pGlobalPoint);
+ const G4ThreeVector localDirection = transform.TransformAxis(pDirection);
+ return fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fSolid->DistanceToOut(localPoint, localDirection );
+}
+
+
+G4VParticleChange*
+WoodcockProcess::PostStepDoIt(const G4Track& track, const G4Step& step) {
+ // Check if the step was done inside the WDCK region with DWCK tracking conditions
+ if (fInWDCKRegion) {
+ // Set the step length in the track in order to give correct value for the user
+ const_cast<G4Track&>(track).SetStepLength( fWDCKStepLength );
+ const_cast<G4Step&>(step).SetStepLength( fWDCKStepLength );
+ // Do nothing else if the step ended on the WDCK region boundary
+ if (fOnWDCKRegionBoundary) {
+ // boundary limited: no interaction
+ // set the number of interaction legth left to restart tracking correctly
+ theNumberOfInteractionLengthLeft = -1.0;
+ fParticleChange.InitializeForPostStep(track);
+// std::cout << " I should set Transportation as limiter becasue actually its the WoodcockProcess itself" << std::endl;
+ return &fParticleChange;
+ }
+ // NOTE: make the user SteppingAction of THIS APPLICATION to see the material
+ // in which the intercation happend! This is becuase of the WDCK!!!
+ step.GetPreStepPoint()->SetMaterialCutsCouple(currentCouple);
+ step.GetPreStepPoint()->SetMaterial(const_cast<G4Material*>(currentCouple->GetMaterial()));
+// step.GetPostStepPoint()->SetMaterialCutsCouple(currentCouple);
+// step.GetPostStepPoint()->SetMaterial(const_cast<G4Material*>(currentCouple->GetMaterial()));
+ return GammaGeneralProcess::PostStepDoIt(track, step);
+ } else {
+ // Normal, i.e. non-WDCK tracking otherwise
+ return GammaGeneralProcess::PostStepDoIt(track, step);
+ }
+ return &fParticleChange;
+}
+
+
+//
+// NOTE: - we need to get the touchable of the Woodcock region root volume in order
+// to obtain the transformations needed at run time to compute the distance
+// to the root volume boundary (i.e. calorimeter boundary)
+// - we obtain it by locating a point that we know that for sure is inside
+// the Woodckock region and we move up till we are in the Woodcock region
+// - the problem is that:
+// 1. we need to know a position that is inside
+// 2. we need to do this after the geometry is closed (i.e. cannot be done
+// at normal process init) only here at start tracking
+//
+// This is why we do these initialisations (only once) when the very first gamma
+// makes its very first step inside the given branch of the WDCK region, i.e. under
+// the given root logical volume of the WDCK region: --> we know a global point
+// that is inside this branch the geometry is already closed.
+//
+void
+WoodcockProcess::InitWoodcockTracking(const G4ThreeVector& aPointInsideTheWDCKregion, const G4ThreeVector& aDirectionInto) {
+ // create the local navigator for this branch of the WDCK region and locate
+ // a point inside this branch then obtain the root volume transformation of this branch
+ G4TouchableHandle theWDCKRegionTouchable;
+ G4Navigator* nav = new G4Navigator();
+ fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fNavigator = nav;
+ nav->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
+ nav->LocateGlobalPointAndUpdateTouchableHandle( aPointInsideTheWDCKregion, aDirectionInto, theWDCKRegionTouchable, false);
+ // Keep moving up till we reach the root logical volume (of the given branch) of the WDCK region
+ const G4LogicalVolume* wdckRootLVol = fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fLogicalVolume;
+ while (wdckRootLVol != theWDCKRegionTouchable->GetVolume()->GetLogicalVolume()) {
+ theWDCKRegionTouchable->MoveUpHistory();
+ }
+ // Obtain the transformation and the solid of the root volume of the WDCK region
+ // that are required for the run-time computation of the distance to out/boundary.
+ const G4NavigationHistory* navHistory = theWDCKRegionTouchable->GetHistory();
+ fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fAffineTransform = navHistory->GetTopTransform();
+ // Get the solid of the WDCK root volume
+ const G4VPhysicalVolume* motherPhysical = navHistory->GetTopVolume();
+ const G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
+ fWDCKDataPerRootLVolume[fWDCKRootLVolumeIndx].fSolid = motherLogical->GetSolid();
+}
+
+
+void WoodcockProcess::ComputeTotalMacrCrossSection(const G4MaterialCutsCouple* couple, const G4double energy, const G4double logenergy) {
+ // compute macroscopic cross section only if the material and/or energy has changed
+ G4bool isRecompute = false;
+ if (couple != currentCouple) {
+ currentCouple = couple;
+ currentCoupleIndex = couple->GetIndex();
+ currentMaterial = couple->GetMaterial();
+ basedCoupleIndex = (*theDensityIdx)[currentCoupleIndex];
+ factor = (*theDensityFactor)[currentCoupleIndex];
+ isRecompute = true;
+ }
+ if (energy != preStepKinEnergy) {
+ preStepKinEnergy = energy;
+ preStepLogE = logenergy;
+ isRecompute = true;
+ }
+ if (isRecompute) {
+ preStepLambda = std::max(0.0, TotalCrossSectionPerVolume());
+ }
+}
+
+
+// This method is supposed to be called only from inside the WDCK region and it
+// will set the index of the root logical volume inside which the current step
+// will be done.
+void WoodcockProcess::FindWDCKRootLogicalVolume(const G4Track& track) {
+ // no need to reset since we have not left the branch yet, i.e. the same root LV
+ if (fWDCKRootLVolumeIndx > -1) { return; }
+ // get the navigation history and keep moving up till a root logical volume of
+ // the WDCK region is found, then set the current root logical volume index
+ const G4VTouchable* tchb = track.GetTouchable();
+ std::map<G4int, G4int>::iterator itr;
+ G4int depth = -1;
+ do {
+ ++depth;
+ G4int id = tchb->GetVolume(depth)->GetLogicalVolume()->GetInstanceID();
+ itr = fWDCKRootLVolIDToIndxMap.find(id);
+ } while ( itr == fWDCKRootLVolIDToIndxMap.end() );
+ fWDCKRootLVolumeIndx = itr->second;
+}