| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4AdjointSimManager
27 //
28 // Class description:
29 //
30 // This class represents the Manager of an adjoint/reverse MC simulation.
31 // An adjoint run is divided in a serie of alternative adjoint and forward
32 // tracking of adjoint and normal particles.
33 //
34 // Reverse tracking phase:
35 // -----------------------
36 // An adjoint particle of a given type (adjoint_e-, adjoint_gamma,...) is
37 // first generated on the so called adjoint source with a random energy (1/E
38 // distribution) and direction. The adjoint source is the external surface
39 // of a user defined volume or of a user defined sphere. The adjoint source
40 // should contain one or several sensitive volumes and should be small compared
41 // to the entire geometry. The user can set the min and max energy of the
42 // adjoint source. After its generation the adjoint primary particle is tracked
43 // bacward in the geometry till a user defined external surface (spherical or
44 // boundary of a volume) or is killed before if it reaches a user defined
45 // upper energy limit that represents the maximum energy of the external
46 // source. During the reverse tracking, reverse processes take place where
47 // the adjoint particle being tracked can be either scattered or transformed
48 // in another type of adjoint paticle. During the reverse tracking the
49 // G4SimulationManager replaces the user defined Primary, Run, ... actions, by
50 // its own actions.
51 //
52 // Forward tracking phase
53 // -----------------------
54 // When an adjoint particle reaches the external surface its weight,type,
55 // position, and directions are registered and a normal primary particle
56 // with a type equivalent to the last generated primary adjoint is generated
57 // with the same energy, position but opposite direction and is tracked
58 // normally in the sensitive region as in a fwd MC simulation. During this
59 // forward tracking phase the event, stacking, stepping, tracking actions
60 // defined by the user for its general fwd application are used. By this clear
61 // separation between adjoint and fwd tracking phases, the code of the user
62 // developed for a fwd simulation should be only slightly modified to adapt it
63 // for an adjoint simulation. Indeed the computation of the signal is done by
64 // the same actions or classes that the one used in the fwd simulation mode.
65 //
66 // Modification to bring in an existing G4 application to use the ReverseMC
67 // ------------------------------------------------------------------------
68 // In order to be able to use the ReverseMC method in his simulation, the
69 // user should modify its code as such:
70 // 1) Adapt its physics list to use ReverseProcesses for adjoint particles.
71 // An example of such physics list is provided in an extended example.
72 // 2) Create an instance of G4AdjointSimManager somewhere in the main code.
73 // 3) Modify the analysis part of the code to normalise the signal computed
74 // during the fwd phase to the weight of the last adjoint particle that
75 // reaches the external surface. This is done by using the following
76 // method of G4AdjointSimManager:
77 //
78 // G4int GetIDOfLastAdjParticleReachingExtSource()
79 // G4ThreeVector GetPositionAtEndOfLastAdjointTrack()
80 // G4ThreeVector GetDirectionAtEndOfLastAdjointTrack()
81 // G4double GetEkinAtEndOfLastAdjointTrack()
82 // G4double GetEkinNucAtEndOfLastAdjointTrack()
83 // G4double GetWeightAtEndOfLastAdjointTrack()
84 // G4double GetCosthAtEndOfLastAdjointTrack()
85 // G4String GetFwdParticleNameAtEndOfLastAdjointTrack()
86 // G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack()
87 // G4int GetFwdParticleIndexAtEndOfLastAdjointTrack().
88 //
89 // In order to have a code working for both forward and adjoint simulation
90 // mode, the extra code needed in user actions for the adjoint simulation
91 // mode can be separated from the code needed only for the normal forward
92 // simulation by using the following method:
93 // G4bool GetAdjointSimMode()
94 // that returns true if an adjoint simulation is running and false if not!
95 //
96 // Example of modification in the analysis part of the code
97 // --------------------------------------------------------
98 // Let's say that in the forward simulation a G4 application computes the
99 // energy deposited in a volume. The user wants to normalise its results for an
100 // external isotropic source of e- with differential spectrum given by f(E). A
101 // possible modification of the code where the deposited energy Edep during an
102 // event is registered would be the following:
103 //
104 // G4AdjointSimManager* theAdjSimManager = G4AdjointSimManager::GetInstance();
105 // if (theAdjSimManager->GetAdjointSimMode())
106 // {
107 // // code of the user that should be consider only for forward simulation
108 // G4double normalised_edep = 0.;
109 // if (theAdjSimManager->GetFwdParticleNameAtEndOfLastAdjointTrack()=="e-")
110 // {
111 // G4double ekin_prim =
112 // theAdjSimManager->GetEkinAtEndOfLastAdjointTrack();
113 // G4double weight_prim =
114 // theAdjSimManager->GetWeightAtEndOfLastAdjointTrack();
115 // normalised_edep = weight_prim*f(ekin_prim);
116 // }
117 // // then follow the code where normalised_edep is printed, or registered
118 // // or whatever ....
119 // }
120 // else
121 // {
122 // // code that should be considered only for forward simulation
123 // }
124 //
125 // Note that in this example a normalisation to only primary e- with only
126 // one spectrum f(E) is considered. The example code could be easily
127 // adapted for a normalisation to several spectra and several types of
128 // primary particles in the same simulation.
129
130 // --------------------------------------------------------------------
131 // Class Name: G4AdjointSimManager
132 // Author: L. Desorgher, 2007-2009
133 // Organisation: SpaceIT GmbH
134 // Contract: ESA contract 21435/08/NL/AT
135 // Customer: ESA/ESTEC
136 // --------------------------------------------------------------------
137 #ifndef G4AdjointSimManager_hh
138 #define G4AdjointSimManager_hh 1
139
140 #include "G4ThreeVector.hh"
141 #include "G4UserRunAction.hh"
142 #include "globals.hh"
143
144 #include <vector>
145
146 class G4UserEventAction;
147 class G4VUserPrimaryGeneratorAction;
148 class G4UserTrackingAction;
149 class G4UserSteppingAction;
150 class G4UserStackingAction;
151 class G4AdjointRunAction;
152 class G4AdjointPrimaryGeneratorAction;
153 class G4AdjointSteppingAction;
154 class G4AdjointEventAction;
155 class G4AdjointStackingAction;
156 class G4AdjointTrackingAction;
157 class G4ParticleDefinition;
158 class G4AdjointSimMessenger;
159 class G4PhysicsLogVector;
160 class G4Run;
161
162 class G4AdjointSimManager : public G4UserRunAction
163 {
164 public:
165 static G4AdjointSimManager* GetInstance();
166
167 void BeginOfRunAction(const G4Run* aRun) override;
168 void EndOfRunAction(const G4Run* aRun) override;
169 void RunAdjointSimulation(G4int nb_evt);
170
171 inline G4int GetNbEvtOfLastRun() { return nb_evt_of_last_run; }
172
173 void SetAdjointTrackingMode(G4bool aBool);
174 G4bool GetAdjointTrackingMode();
175 // true if an adjoint track is being processed
176
177 inline G4bool GetAdjointSimMode()
178 {
179 return adjoint_sim_mode;
180 } // true if an adjoint simulation is running
181
182 G4bool GetDidAdjParticleReachTheExtSource();
183 void RegisterAtEndOfAdjointTrack();
184 void RegisterAdjointPrimaryWeight(G4double aWeight);
185 void ResetDidOneAdjPartReachExtSourceDuringEvent();
186
187 inline G4int GetIDOfLastAdjParticleReachingExtSource()
188 {
189 return ID_of_last_particle_that_reach_the_ext_source;
190 }
191
192 G4ThreeVector GetPositionAtEndOfLastAdjointTrack(std::size_t i = 0);
193 G4ThreeVector GetDirectionAtEndOfLastAdjointTrack(std::size_t i = 0);
194 G4double GetEkinAtEndOfLastAdjointTrack(std::size_t i = 0);
195 G4double GetEkinNucAtEndOfLastAdjointTrack(std::size_t i = 0);
196 G4double GetWeightAtEndOfLastAdjointTrack(std::size_t i = 0);
197 G4double GetCosthAtEndOfLastAdjointTrack(std::size_t i = 0);
198 const G4String& GetFwdParticleNameAtEndOfLastAdjointTrack();
199 G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(std::size_t i = 0);
200 G4int GetFwdParticleIndexAtEndOfLastAdjointTrack(std::size_t i = 0);
201 std::size_t GetNbOfAdointTracksReachingTheExternalSurface();
202 void ClearEndOfAdjointTrackInfoVectors();
203 G4ParticleDefinition* GetLastGeneratedFwdPrimaryParticle();
204
205 std::vector<G4ParticleDefinition*>* GetListOfPrimaryFwdParticles();
206 std::size_t GetNbOfPrimaryFwdParticles();
207
208 G4bool DefineSphericalExtSource(G4double radius, G4ThreeVector pos);
209 G4bool DefineSphericalExtSourceWithCentreAtTheCentreOfAVolume(G4double radius,
210 const G4String& volume_name);
211 G4bool DefineExtSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name);
212 void SetExtSourceEmax(G4double Emax);
213
214 // Definition of adjoint source
215 //----------------------------
216 G4bool DefineSphericalAdjointSource(G4double radius, G4ThreeVector pos);
217 G4bool DefineSphericalAdjointSourceWithCentreAtTheCentreOfAVolume(G4double radius,
218 const G4String& volume_name);
219 G4bool DefineAdjointSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name);
220 void SetAdjointSourceEmin(G4double Emin);
221 void SetAdjointSourceEmax(G4double Emax);
222 inline G4double GetAdjointSourceArea() { return area_of_the_adjoint_source; }
223 void ConsiderParticleAsPrimary(const G4String& particle_name);
224 void NeglectParticleAsPrimary(const G4String& particle_name);
225 void SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon);
226 const G4String& GetPrimaryIonName();
227
228 inline void SetNormalisationMode(G4int n) { normalisation_mode = n; }
229 inline G4int GetNormalisationMode() { return normalisation_mode; }
230 inline G4double GetNumberNucleonsInIon() { return nb_nuc; }
231
232 // Definition of user actions for the adjoint tracking phase
233 //----------------------------
234 void SetAdjointEventAction(G4UserEventAction* anAction);
235 void SetAdjointSteppingAction(G4UserSteppingAction* anAction);
236 void SetAdjointStackingAction(G4UserStackingAction* anAction);
237 void SetAdjointRunAction(G4UserRunAction* anAction);
238
239 // Set methods for user run actions
240 //--------------------------------
241 inline void UseUserStackingActionInFwdTrackingPhase(G4bool aBool)
242 {
243 use_user_StackingAction = aBool;
244 }
245 inline void UseUserTrackingActionInFwdTrackingPhase(G4bool aBool)
246 {
247 use_user_TrackingAction = aBool;
248 }
249
250 // Set nb of primary fwd gamma
251 void SetNbOfPrimaryFwdGammasPerEvent(G4int);
252
253 // Set nb of adjoint primaries for reverse splitting
254 void SetNbAdjointPrimaryGammasPerEvent(G4int);
255 void SetNbAdjointPrimaryElectronsPerEvent(G4int);
256
257 // Convergence test
258 //-----------------------
259 /*
260 void RegisterSignalForConvergenceTest(G4double aSignal);
261 void DefineExponentialPrimarySpectrumForConvergenceTest(
262 G4ParticleDefinition* aPartDef, G4double E0);
263 void DefinePowerLawPrimarySpectrumForConvergenceTest(
264 G4ParticleDefinition* aPartDef, G4double alpha);
265 */
266
267 void SwitchToAdjointSimulationMode();
268 void BackToFwdSimulationMode();
269
270 private: // methods
271 static G4ThreadLocal G4AdjointSimManager* instance;
272
273 void SetRestOfAdjointActions();
274 void SetAdjointPrimaryRunAndStackingActions();
275 void SetAdjointActions();
276 void ResetRestOfUserActions();
277 void ResetUserPrimaryRunAndStackingActions();
278 void ResetUserActions();
279 void DefineUserActions();
280
281 // private constructor and destructor
282 G4AdjointSimManager();
283 ~G4AdjointSimManager() override;
284
285 private: // attributes
286 // Messenger
287 //----------
288 G4AdjointSimMessenger* theMessenger = nullptr;
289
290 // user defined actions for the normal fwd simulation.
291 // Taken from the G4RunManager
292 //-------------------------------------------------
293 G4bool user_action_already_defined = false;
294 G4UserRunAction* fUserRunAction = nullptr;
295 G4UserEventAction* fUserEventAction = nullptr;
296 G4VUserPrimaryGeneratorAction* fUserPrimaryGeneratorAction = nullptr;
297 G4UserTrackingAction* fUserTrackingAction = nullptr;
298 G4UserSteppingAction* fUserSteppingAction = nullptr;
299 G4UserStackingAction* fUserStackingAction = nullptr;
300 G4bool use_user_StackingAction = false; // only for fwd part of adjoint sim
301 G4bool use_user_TrackingAction = false;
302
303 // action for adjoint simulation
304 //-----------------------------
305 G4UserRunAction* theAdjointRunAction = nullptr;
306 G4UserEventAction* theAdjointEventAction = nullptr;
307 G4AdjointPrimaryGeneratorAction* theAdjointPrimaryGeneratorAction = nullptr;
308 G4AdjointTrackingAction* theAdjointTrackingAction = nullptr;
309 G4AdjointSteppingAction* theAdjointSteppingAction = nullptr;
310 G4AdjointStackingAction* theAdjointStackingAction = nullptr;
311
312 // adjoint mode
313 //-------------
314 G4bool adjoint_tracking_mode = false;
315 G4bool adjoint_sim_mode = false;
316
317 // adjoint particle information on the external surface
318 //-----------------------------
319 std::vector<G4ThreeVector> last_pos_vec;
320 std::vector<G4ThreeVector> last_direction_vec;
321 std::vector<G4double> last_ekin_vec;
322 std::vector<G4double> last_ekin_nuc_vec;
323 std::vector<G4double> last_cos_th_vec;
324 std::vector<G4double> last_weight_vec;
325 std::vector<G4int> last_fwd_part_PDGEncoding_vec;
326 std::vector<G4int> last_fwd_part_index_vec;
327 std::vector<G4int> ID_of_last_particle_that_reach_the_ext_source_vec;
328
329 G4ThreeVector last_pos;
330 G4ThreeVector last_direction;
331 G4double last_ekin = 0.0, last_ekin_nuc = 0.0;
332 // last_ekin_nuc=last_ekin/nuc, nuc is 1 if not a nucleus
333 G4double last_cos_th = 0.0;
334 G4String last_fwd_part_name;
335 G4int last_fwd_part_PDGEncoding = 0;
336 G4int last_fwd_part_index = 0;
337 G4double last_weight = 0.0;
338 G4int ID_of_last_particle_that_reach_the_ext_source = 0;
339
340 G4int nb_evt_of_last_run = 0;
341 G4int normalisation_mode = 3;
342
343 // Adjoint source
344 //--------------
345 G4double area_of_the_adjoint_source = 0.0;
346 G4double nb_nuc = 1.0;
347 G4double theAdjointPrimaryWeight = 0.0;
348
349 G4bool welcome_message = true;
350 };
351
352 #endif
353