Warning, /geant4/examples/extended/biasing/ReverseMC01/README.md is written in an unsupported language. File is not indexed.
0001 \page ExampleReverseMC01 Example ReverseMC01
0002
0003 This example illustrates the use of Reverse Monte Carlo in Geant4.
0004
0005 ## Author
0006
0007 This example code and the adjoint classes in the G4 toolkit have been developed by L.Desorgher (SpaceIT GmbH)
0008 under the ESA contract 21435/08/NL/AT. For any (reasonable) question you may contact the author
0009 at the following email address : desorgher@spaceit.ch
0010
0011
0012
0013 ## Abstract
0014
0015 This is the README file for the first G4 example illustrating the use of the Reverse Monte Carlo (RMC) mode in a Geant4
0016 application. The Reverse Monte Carlo method is also known as the Adjoint Monte Carlo (AMC) method and
0017 in this document we will alternate both Reverse and Adjoint terms.
0018
0019 ## Other documentation
0020
0021 See also the section 3.7.3 Adjoint/Reverse Monte carlo in the
0022 Geant4 User guide for application developers.
0023
0024
0025 ## Definition of Reverse/Adjoint Monte Carlo
0026
0027 When the sensitive part of a detector is small compared to its entire size and to the size of the
0028 external extended primary particle source, a lot of computing time is spent during a normal Monte Carlo run
0029 in the simulation of particle showers that are not contributing to the detector signal.
0030 In such particular case the Reverse Monte Carlo (RMC) method, also known as the
0031 Adjoint Monte Carlo method, can be used.
0032 In this method particles are generated in or on the external surface of the sensitive volume
0033 of the instrument and then are tracked backward in the geometry till they reach the source surface,
0034 or exceed an energy threshold. During the reverse tracking reverse reactions are applied to the particles.
0035
0036
0037
0038 ## The Reverse Monte Carlo mode in Geant4 (since G4.9.3 release)
0039
0040 (See also the section 3.7.3 Adjoint/Reverse Monte carlo in the
0041 Geant4 User guide for application developers.)
0042
0043 Different G4Adjoint classes have been implemented into the Geant4
0044 toolkit to run an adjoint/reverse simulation in a Geant4 application.
0045 In this implementation an adjoint run is divided in a succession
0046 of alternative adjoint and forward tracking of adjoint and normal particles.
0047 One Geant4 event treats one of this tracking phase.
0048
0049
0050 ### Reverse tracking phase
0051
0052 Adjoint particles (adjoint_e-, adjoint_gamma,...) are generated one by one on the so called
0053 adjoint source with random position, energy (1/E distribution) and direction. The adjoint
0054 source is the external surface of a user defined volume or of a user defined sphere. The
0055 adjoint source should contain one or several sensitive volumes and should be small
0056 compared to the entire geometry. The user can set the minimum and maximum energy of the
0057 adjoint source. After its generation the adjoint primary particle is tracked backward in
0058 the geometry till a user defined external surface (spherical or boundary of a volume)
0059 or is killed before if it reaches a user defined upper energy limit that represents the
0060 maximum energy of the external source. During the reverse tracking, reverse processes take
0061 place where the adjoint particle being tracked can be either scattered or transformed in
0062 another type of adjoint particle. During the reverse tracking the
0063 G4AdjointSimulationManager replaces the user defined primary, run, stepping, ... actions,
0064 by its own actions.
0065
0066 ### Forward tracking phase:
0067
0068 When an adjoint particle reaches the external surface its weight, type, position,
0069 and direction are registered and a normal primary particle with a type equivalent
0070 to the last generated adjoint primary is generated with the same energy,
0071 position but opposite direction and is tracked in the forward direction
0072 in the sensitive region as in a forward MC simulation.
0073 During this forward tracking phase the event, stacking, stepping, tracking actions defined
0074 by the user for its general forward application are used.
0075 By this clear separation between adjoint and forward tracking phases, the code of the
0076 user developed for a forward simulation should be only slightly
0077 modified to adapt it for an adjoint simulation. Indeed the computation of the signal
0078 is done by the same user actions or analysis classes that the one used in the forward
0079 simulation mode. Before the G4.10.0 release the reverse and forward tracking mode
0080 took place in separated events. Since the G4.10.0 release,
0081 in order to preapre to the migration of the
0082 ReverseMC to the G4 Multiple Threading mode, the reverse and forward tracking
0083 phase of corresponding adjoint and forward primaries have been merged in the same
0084 event.
0085
0086
0087 ### Reverse Processes
0088
0089 During the reverse tracking phase reverse processes act on the adjoint particles.
0090 The Reverse processes that are available at the moment in Geant4 are the:
0091 - Reverse discrete Ionization for e-, proton and ions
0092 - Continuous gain of energy by ionization and bremsstrahlung for e- and by ionization for protons and ions
0093 - Reverse discrete e- bremsstrahlung
0094 - Reverse photoelectric effect
0095 - Reverse Compton scattering
0096 - Approximated multiple scattering (MS) (see section 5.3)
0097
0098
0099 It is important to note that the electromagnetic reverse processes are cut dependent
0100 as their equivalent forward processes. The implementation of the reverse processes is
0101 based on the forward processes
0102 implemented in the G4 standard electromagnetic package.
0103
0104
0105 ### Remark on Nb of adjoint particle types and Nb of G4 events considered in an adjoint simulation
0106
0107 The list of type of adjoint and forward particles that are generated on the adjoint source
0108 and considered in the simulation is a function of the adjoint processes declared in the
0109 physics list. For example if only the e- and gamma electromagnetic processes are considered
0110 , only adjoint e- and adjoint gamma will be considered as primaries. In this case an
0111 adjoint event will be divided in two G4 events. The first event will consist
0112 into the coupled reverse and forward tracking of an adjoint e- and its equivalent
0113 forward e-, while the second events will process the reverse and forward trackings
0114 of corresponsing adjoint and forward primary gamms. In this case a
0115 run of 100 adjoint events will consist into 200 Geant4 events. If the proton ionization is
0116 also considered adjoint and forward protons are also generated as primaries
0117 and 300 Geant4 events are processed for 100 adjoint events.
0118
0119 ### Modifications to bring in a existing G4 application to use the Reverse MC method
0120
0121 (for more details see also the section 3.7.3 Adjoint/Reverse Monte carlo in the
0122 Geant4 User guide for application developers.)
0123
0124 Due the clear separation between the reverse and forward tracking phase only few modifications are needed
0125 to an existing Geant4 application in order to adapt it for the use of the reverse simulation mode.
0126 Except in the physics list where all the reverse processes and their forward equivalent
0127 have to be declared, the principal code modifications are needed only in the analysis phase at the end
0128 of the forward tracking where computed signals have to be multiplied by the weight
0129 of the last reverse tracks and then normalized to different user defined spectra and angular distribution representing
0130 the external source.
0131 The weight of the adjoint tracks is computed by the G4Adjoint classes and the user needs
0132 only to multiply them to the primary differential, directional spectrum of its choice.
0133 The adjoint weight a the end of tracks can be also registered if needed in answer matrices.
0134
0135 More precisely, in order to be able to use the Reverse MC method in his simulation, the user should modify
0136 its code as such:
0137
0138 - Adapt its physics list to use Reverse Processes for adjoint particles. An example of such physics list is provided in an extended
0139 example.
0140 - Create an instance of G4AdjointSimManager somewhere in the main() code.
0141
0142 - Modify the analysis part of the code to normalize the signal computed during the forward phase to the weight of the last adjoint particle
0143 that reaches the external surface. This is done by using the following method of G4AdjointSimManager:
0144 - G4int GetIDOfLastAdjParticleReachingExtSource()
0145 - G4ThreeVector GetPositionAtEndOfLastAdjointTrack(){ return last_pos;}
0146 - G4ThreeVector GetDirectionAtEndOfLastAdjointTrack(){ return last_direction;}
0147 - G4double GetEkinAtEndOfLastAdjointTrack(){ return last_ekin;}
0148 - G4double GetEkinNucAtEndOfLastAdjointTrack(){ return last_ekin_nuc;}
0149 - G4double GetWeightAtEndOfLastAdjointTrack(){return last_weight;}
0150 - G4double GetCosthAtEndOfLastAdjointTrack(){return last_cos_th;}
0151 - G4String GetFwdParticleNameAtEndOfLastAdjointTrack(){return last_fwd_part_name;}
0152 - G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(){return last_fwd_part_PDGEncoding;}
0153 - G4int GetFwdParticleIndexAtEndOfLastAdjointTrack().
0154
0155 In order to have a code working for both forward and adjoint simulation mode, the extra code needed in user actions for the adjoint
0156 simulation mode can be separated to the code needed only for the normal forward simulation by using the following method:
0157 - G4bool GetAdjointSimMode() that return true if an adjoint simulation is running and false if not!
0158
0159
0160
0161 ## exampleRMC01
0162
0163 The example RMC01 illustrates how to modify a G4 application in order to use
0164 both forward and reverse MC modes in the same code.
0165
0166
0167 ### Geometry
0168
0169 The following simple geometry is considered:
0170 - sensitive Silicon cylinder at the center of an Aluminum spherical shielding with 10 cm Radius.
0171 - two 0.5mm thick Tantalum plates set horizontally above and below the Sensitive Cylinder
0172
0173 The free parameters of the geometry that can bes set by the user are:
0174 - the thickness of the Aluminum shielding
0175 - the height of the sensitive Si cylinder
0176 - the radius of the sensitive Si cylinder
0177
0178
0179
0180 ### Physics
0181
0182 The physical processes considered are:
0183 - Reverse and forward discrete Ionization for e- and proton
0184 - Continuous gain and loss of energy by ionization and bremsstrahlung for e- and by ionization for protons
0185 - Reverse and forward discrete e- bremsstrahlung
0186 - Reverse and forward photoelectric effect
0187 - Reverse and forward Compton scattering
0188 - Reverse and forward Multiple scattering
0189
0190 These processes are implemented in the class G4AdjointPhysicsList distributed with the example. The G4AdjointPhysicsMessenger allows the user
0191 to switch on/off some processes for testing purpose. By default all processes cited above are considered except the proton ionization that
0192 has to be specifically switch on in the macro file by the user.
0193
0194
0195
0196 ### Analysis and output of the code
0197
0198 The example computes the energy deposited in the sensitive Si cylinder and the current of e-, protons, and gamma
0199 entering this cylinder.
0200 The Hits are registered in the sensitive detector class RMC01SD that is a typical G4 sensitive detector class
0201 used in a forward simulation and is not modified at all
0202 for the adjoint simulation mode.
0203 The analysis of the registered hits during forward events is done by the RMCO1AnalysisManager.
0204 That is the class that illustrates how to adapt an analysis code of a fwd simulation in order to use it also for
0205 an adjoint simulation.
0206 In this class during a forward simulation the method EndOfEventForForwardSimulation is used at the end of an event
0207 while during an adjoint simulation at the end of fwd tracking event the method EndOfEventForAdjointSimulation is called.
0208 By looking at the source of RMCO1AnalysisManager and more particularly to its method EndOfEventForAdjointSimulation the user will
0209 learn how to adapt its G4 analysis code for an adjoint simulation.
0210
0211 The outputs of an adjoint simulation are:
0212
0213 - The total energy deposited and particle current entering the sensitive cylinder normalized
0214 automatically to a user defined primary spectrum(exponential or power law) .These results are stored in the files:
0215 - Adj_Edep_vs_EkinPrim.txt
0216 - Adj_ElectronCurrent.txt
0217 - Adj_GammaCurrent.txt
0218 - Adj_ProtonCurrent.txt
0219 - ConvergenceOfAdjointSimulationResults.txt:
0220 The total normalized edep and its relative error registered every 5000 adjoint events
0221
0222
0223 - The answer matrix of the energy deposited and particles current on the sensitive cylinder in function of primary energy of e-, gamma and
0224 protons. These results are stored in the files Adj********_Answer.txt
0225
0226
0227
0228 The outputs of a forward simulation are:
0229 - The mean energy deposited and particle current entering the sensitive cylinder per event.
0230 These results are stored in the files:
0231 - Fwd_Edep_vs_EkinPrim.txt
0232 - Fwd_ElectronCurrent.txt
0233 - Fwd_GammaCurrent.txt
0234 - Fwd_ProtonCurrent.txt
0235
0236
0237
0238 ### Run macrofiles
0239
0240 The following example run macro files are distributed with the code:
0241 - run_adjoint_simulation_electron.mac and run_adjoint_simulation_proton.mac for adjoint simulations
0242
0243 - run_forward_simulation_electron.mac and run_forward_simulation_proton.mac for forward simulations
0244
0245
0246 ### Comparison of adjoint and forward simulation results
0247
0248 It is the responsibility of the user to select in the macro file the same external spectrum
0249 for both the forward and adjoint simulations and to normalize the per event results of the forward simulation
0250 to the fluence considered in the adjoint simulation.
0251
0252 For the macro files that are provided with the examples it consists into multiplying the forward results by pi*100.
0253 This normalization factor is explained by the following:
0254
0255 - For the forward simulation the results are given per number of events. It corresponds
0256 to a normalization to a fluence of 1 particle emanating from the external source.
0257
0258 - In run_fwd_simulation.mac the source is set on a sphere of 10 cm radius (see /gps commands in
0259 macrofile).Therefore the omnidirectional fluence for the fwd simulation is 1./(pi*R^2) with R=10cm.
0260
0261 - The adjoint results are normalized to a fluence of 1/cm2.
0262 (See command /RMC01/analysis/SetExponentialSpectrumForAdjointSim in macrofile)
0263
0264 - In conclusion to compare the adjoint and forward results, the forward results should
0265 be multiplied by pi*R^2/cm2= pi*100.
0266
0267
0268
0269 ## Control of the adjoint simulation and the RMC01 code by G4 macro UI commands
0270
0271 Different G4 macro UI commands are provided to control the RMC01 example and the adjoint simulation.
0272 Some macro commands are provided within the geant4 toolkit and appears in a G4 application when the singleton
0273 class G4AdjointSimManager is called somewhere in the code, the other macro commands are
0274 declared in the code distributed within the example.
0275
0276
0277 ### G4UI commands in the directory /adjoint
0278
0279 The macro command directory /adjoint appears in a user application when the singleton
0280 class G4AdjointSimManager is called somewhere in the code.
0281 It allows to control the adjoint source, the external source and start an adjoint simulation.
0282
0283 The command to start an adjoint run is:
0284
0285 - /adjoint/start_run nb \n
0286 Start an adjoint simulation with a number of events given by nb. It is important to note that the total number of events in the sense of G4
0287 will be nb*2*nb_primary_considered (see 3.4.)
0288
0289
0290 The commands to control the adjoint source are:
0291
0292 - /adjoint/DefineSphericalAdjSource R X Y Z unit_length \n
0293 The adjoint source is set on a sphere with radius R and centered on position (X,Y,Z)
0294
0295 - /adjoint/DefineSphericalAdjSourceCenteredOnAVolume phys_vol_name R unit_length \n
0296 The external source is set on a sphere with radius R and with its center position located at the center of the
0297 the physical volume specified by the name phys_vol_name.
0298 - /adjoint/DefineAdjSourceOnExtSurfaceOfAVolume phys_vol_name \n
0299 The external surface is set as the external boundary of a the physical volume with name phys_vol_name
0300
0301 - /adjoint/SetAdjSourceEmin Emin energy_unit \n
0302 Set the minimum energy of the external source
0303
0304 - /adjoint/SetAdjSourceEmax Emax energy_unit \n
0305 Set the maximum energy of the external source
0306
0307 - /adjoint/ConsiderAsPrimary particle_name \n
0308 The type of particle specified by "particle_name" will be added in the list of primary adjoint particles.
0309 The list of candidates depends on the reverse physics processes considered in the simulation. At the most the
0310 potential candidates are (e-, gamma, proton , ion). For this example only e-, gamma, proton
0311 can be chosen. As the proton ionization is not considered by default, the default list of particles is
0312 [e-,gamma]. To have also the proton as candidate the proton ionization should
0313 be switch on (/adjoint_physics/UseProtonIonisation true).
0314
0315 - /adjoint/NeglectAsPrimary particle_name \n
0316
0317 The type of particle specified by "particle_name" will be removed from the list of primary adjoint particles.
0318 The list of candidates depends on the reverse physics processes considered in the simulation. At the most the
0319 potential candidates are (e-, gamma, proton , ion). For this example only e-, gamma, proton
0320 can be chosen. As the proton ionization is not considered by default, the default list of particles is
0321 [e-,gamma].To have also the proton as candidate the proton ionization should
0322 be switch on (/adjoint_physics/UseProtonIonisation true).
0323
0324
0325 The commands to control the external source are:
0326
0327 - /adjoint/DefineSphericalExtSource R X Y Z unit_length:\n
0328 The external source is set on a sphere with radius R and centered on position (X,Y,Z)
0329
0330 - /adjoint/DefineSphericalExtSourceCenteredOnAVolume phys_vol_name R unit_length\n
0331 The external source is set on a sphere with radius R and with its center position located at the center of the
0332 the physical volume specified by the name phys_vol_name.
0333
0334 - /adjoint/DefineExtSourceOnExtSurfaceOfAVolume phys_vol_name \n
0335 The external surface is set as the external boundary of a the physical volume with name phys_vol_name
0336
0337 - /adjoint/SetExtSourceEmax Emax energy_unit \n
0338 Set the maximum energy of the external source. An adjoint track will be stop when a an adjoint particle get an energy higher than this maximum energy.
0339
0340
0341
0342 ### G4UI commands in the directory /adjoint_physics
0343
0344 These commands allow to control the electromagnetic processes that will be considered in the simulation.
0345
0346 The processes that can be used are:
0347 - Reverse and forward e- continuous and discrete Ionization. Always switch on
0348 - Reverse and forward e- Bremsstrahlung. Switch on by default
0349 - Reverse and forward Compton scattering. Switch on by default
0350 - Reverse and forward photo electric effect. Switch on by default
0351 - Reverse and forward photo electric effect. Switch on by default
0352 - Reverse and forward multiple scattering. Switch on by default
0353 - Reverse and forward proton continuous and discrete Ionization. Switch off by default
0354 - Forward e-e+ pair production. Switch off by default.
0355
0356
0357 The commands that can be used to switch on of these processes are:
0358
0359 - /adjoint_physics/UseProtonIonisation true/false \n
0360 Switch on/off the reverse and forward proton ionization. Off by default.
0361
0362 - /adjoint_physics/UseBremsstrahlung true/false \n
0363 Switch on/off the reverse and forward e- bremsstrahlung. On by default.
0364
0365 - /adjoint_physics/UseCompton true/false \n
0366 Switch on/off the Compton scattering. On by default.
0367
0368
0369 - /adjoint_physics/UseMS true/false \n
0370 Switch on/off the multiple scattering. On by default.
0371
0372
0373 - /adjoint_physics/UseEgainElossFluctuation true/false \n
0374 Switch on/off the fluctuation in the continuous energy loss/gain. On by default. Only for test purpose.
0375
0376 - /adjoint_physics/UsePEEffect true/false \n
0377 Switch on/off the photo electric effect. On by default.
0378
0379
0380 - /adjoint_physics/UseGammaConversion true/false \n
0381 Switch on/off the forward e-e+ pair production from gamma. Off by default. When On all the e+
0382 electromagnetic physics is considered.
0383
0384
0385 The user can also fix the maximum energy Emax and minimum energy Emin of the adjoint physical processes used
0386 in the simulation. The adjoint process will be applied to particles within the energy range [Emin, Emax]
0387 and will produce adjoint secondary only in this energy range. It is recommended to fix Emin to the minimum
0388 energy of the adjoint source and fix Emax to the maximum energy of the external source.
0389 The commands controlling Emin and Emax are:
0390
0391 - /adjoint_physics/SetEminForAdjointModels Emin Energy_unit \n
0392 Set the minimum energy of the adjoint processes/models.
0393
0394 - /adjoint_physics/SetEmaxForAdjointModels Emin Energy_unit \n
0395 Set the maximum energy of the adjoint processes/models.
0396
0397
0398 ### G4UI commands in the directory /RMC01
0399
0400
0401 Commands/RMC01/geometry/ to control the geometry:
0402
0403 - /RMC01/geometry/SetSensitiveVolumeHeight H length_unit \n
0404 Set the height H of the Si sensitive cylinder.
0405
0406
0407 - /RMC01/geometry/SetSensitiveVolumeRadius R length_unit \n
0408 Set the radius R of the Si sensitive cylinder.
0409
0410 - /RMC01/geometry/SetShieldingThickness D length_unit \n
0411 Set the thickness D of the aluminum shielding.
0412
0413 Commands /RMC01/analysis/ to control the primary spectrum used for the normalization of the
0414 adjoint simulation results and fix the expected precision of the computed Edep:
0415
0416 - /RMC01/analysis/SetPowerLawPrimSpectrumForAdjointSim particle_name F F_unit alpha Emin Emax E_unit \n
0417 Set the primary spectrum to which the adjoint simulation results will be normalised to a power law
0418 spectrum E^(-alpha) of particle defined by particle_name, with an omnidirectional fluence F, and
0419 energy range [Emin,Emax]. The fluence unit candidates for F_unit are [1/cm2, 1/m2, cm-2, m-2].
0420
0421
0422 - /RMC01/analysis/SetExponentialSpectrumForAdjointSim particle_name F F_unit E0 Emin Emax E_unit \n
0423 Set the primary spectrum to which the adjoint simulation results will be normalised to an exponential
0424 spectrum exp(-E/E0) of particle defined by particle_name, with an omnidirectional fluence F, and
0425 energy range [Emin,Emax]. The fluence unit candidates for F_unit are [1/cm2, 1/m2, cm-2, m-2].
0426
0427
0428
0429 - /RMC01/analysis/SetExpectedPrecisionOfResults precision \n
0430 Set the expected precision in % for the computed energy deposited in the sensitive volume
0431 for both the forward and adjoint simulation case. When the relative statistical error
0432 of the computed energy deposited reach this precision the run is aborted and the results are registered.
0433 Otherwise the run continue till the nb of events specified by the user are processed. By default the precision is set
0434 to 0. meaning that the run will not be aborted in this case.
0435
0436
0437
0438
0439
0440
0441 ## Known issues
0442
0443 ### Rare too high weight in the adjoint simulation
0444
0445 In rare cases an adjoint track may get a much too high weight when reaching the external source.
0446 While this happen not often it may corrupt the simulation results significantly. The reason of this high weight is
0447 the joint use at low e- and gamma energy of both the photoelectric and bremsstrahlung processes.
0448 Unfortunately we still need some investigations to remove this problem at the level of physical processes.
0449 However this problem can be solved at the level of event action in the user code by adding a test on the adjoint
0450 weight. Such test has been implemented in the example RMC01.
0451 In this implementation an event is rejected when the relative error of the computed normalised edep
0452 increase during one event by more than 50% when the precision is already below 10%.
0453
0454
0455 ### Limitation of the reverse bremsstrahlung
0456
0457 The difference between the differential cross sections used in the adjoint and forward bremsstrahlung
0458 models is the source of a higher flux of >100 keV gamma in the reverse simulation compared to the forward simulation.
0459 The adjoint processes/models should make use of the direct differential cross section to sample
0460 the adjoint secondaries and compute the adjoint cross section.
0461 The differential cross section used in G4AdjointeBremstrahlungModel is obtained by the numerical derivation
0462 over the cut energy of the direct cross section provided by G4eBremsstrahlungModel.
0463 This would be a correct procedure if the distribution of secondary in G4eBremsstrahlungModel
0464 would match this differential cross section. Unfortunately it is not the case as independent parameterization are used
0465 in G4eBremsstrahlungModel for both the cross sections and the sample of secondary. (It means that in the forward case
0466 if one would integrate the effective differential cross section considered in the simulation we would not find back
0467 the used cross section).
0468 In the future we plan to correct this problem by using an extra weight correction factor after the occurrence of a reverse
0469 bremsstrahlung. This weight factor should be the ratio between the differential CS used in the adjoint simulation and the
0470 one effectively used in the forward processes. As it is impossible to have access to the forward differential CS
0471 in G4eBremsstrahlungModel we are investigating the feasibility to use the differential CS considered in
0472 G4Penelope models.
0473
0474
0475 ### Limitation of the reverse multiple scattering
0476
0477 For the reverse multiple scattering we are using the same models than for the forward case.
0478 This approximation makes that the discrepancy between the adjoint and forward
0479 simulation cases can get to a level of ~ 10-15% relative differences in the test cases that we have considered.
0480 In the future we plan to improve the adjoint multiple scattering models by forcing the computation of
0481 multiple scattering effect at the end of an adjoint step.