Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:57:52

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4AdjointSimManager
0027 //
0028 // Class description:
0029 //
0030 // This class represents the Manager of an adjoint/reverse MC simulation.
0031 // An adjoint run is divided in a serie of alternative adjoint and forward
0032 // tracking of adjoint and normal particles.
0033 //
0034 //   Reverse tracking phase:
0035 //   -----------------------
0036 // An adjoint particle of a given type  (adjoint_e-, adjoint_gamma,...) is
0037 // first generated on the so called adjoint source with a random energy (1/E
0038 // distribution) and direction. The adjoint source is the external surface
0039 // of a user defined volume or of a user defined sphere. The adjoint source
0040 // should contain one or several sensitive volumes and should be small compared
0041 // to the entire geometry. The user can set the min and max energy of the
0042 // adjoint source. After its generation the adjoint primary  particle is tracked
0043 // bacward in the geometry till a user  defined external surface (spherical or
0044 // boundary of a volume) or is killed before if it reaches a user defined
0045 // upper energy limit that represents the maximum energy of the external
0046 // source. During the reverse tracking, reverse processes take place where
0047 // the adjoint particle being tracked can be either scattered or transformed
0048 // in another type of adjoint paticle. During the reverse tracking the
0049 // G4SimulationManager replaces the user defined Primary, Run, ... actions, by
0050 // its own actions.
0051 //
0052 //   Forward tracking phase
0053 //   -----------------------
0054 // When an adjoint particle reaches the external surface its weight,type,
0055 // position, and directions are registered and a normal primary particle
0056 // with a type equivalent to the last generated primary adjoint is generated
0057 // with the same energy, position but opposite direction and is tracked
0058 // normally in the sensitive region as in a fwd MC simulation. During this
0059 // forward tracking phase the event, stacking, stepping, tracking actions
0060 // defined by the user for its general fwd application are used. By this clear
0061 // separation between adjoint and fwd tracking phases, the code of the user
0062 // developed for a fwd simulation should be only slightly modified to adapt it
0063 // for an adjoint simulation. Indeed  the computation of the signal is done by
0064 // the same actions or classes that the one used in the fwd simulation mode.
0065 //
0066 //   Modification to bring in an existing G4 application to use the ReverseMC
0067 //   ------------------------------------------------------------------------
0068 // In order to be able to use the ReverseMC method in his simulation, the
0069 // user should modify its code as such:
0070 //  1) Adapt its physics list to use ReverseProcesses for adjoint particles.
0071 //     An example of such physics list is provided in an extended example.
0072 //  2) Create an instance of G4AdjointSimManager somewhere in the main code.
0073 //  3) Modify the analysis part of the code to normalise the signal computed
0074 //     during the fwd phase to the weight of the last adjoint particle that
0075 //     reaches the external surface. This is done by using the following
0076 //     method of G4AdjointSimManager:
0077 //
0078 //     G4int GetIDOfLastAdjParticleReachingExtSource()
0079 //     G4ThreeVector GetPositionAtEndOfLastAdjointTrack()
0080 //     G4ThreeVector GetDirectionAtEndOfLastAdjointTrack()
0081 //     G4double GetEkinAtEndOfLastAdjointTrack()
0082 //     G4double GetEkinNucAtEndOfLastAdjointTrack()
0083 //     G4double GetWeightAtEndOfLastAdjointTrack()
0084 //     G4double GetCosthAtEndOfLastAdjointTrack()
0085 //     G4String GetFwdParticleNameAtEndOfLastAdjointTrack()
0086 //     G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack()
0087 //     G4int GetFwdParticleIndexAtEndOfLastAdjointTrack().
0088 //
0089 // In order to have a code working for both forward and adjoint simulation
0090 // mode, the extra code needed in user actions for the adjoint simulation
0091 // mode can be separated from the code needed only for the normal forward
0092 // simulation by using the following method:
0093 //   G4bool GetAdjointSimMode()
0094 // that returns true if an adjoint simulation is running and false if not!
0095 //
0096 //   Example of modification in the analysis part of the code
0097 //   --------------------------------------------------------
0098 // Let's say that in the forward simulation a G4 application computes the
0099 // energy deposited in a volume. The user wants to normalise its results for an
0100 // external isotropic source of e- with differential spectrum given by f(E). A
0101 // possible modification of the code where the deposited energy Edep during an
0102 // event is registered would be the following:
0103 //
0104 //   G4AdjointSimManager* theAdjSimManager = G4AdjointSimManager::GetInstance();
0105 //   if (theAdjSimManager->GetAdjointSimMode())
0106 //   {
0107 //     // code of the user that should be consider only for forward simulation
0108 //     G4double normalised_edep = 0.;
0109 //     if (theAdjSimManager->GetFwdParticleNameAtEndOfLastAdjointTrack()=="e-")
0110 //     {
0111 //       G4double ekin_prim =
0112 //         theAdjSimManager->GetEkinAtEndOfLastAdjointTrack();
0113 //       G4double weight_prim =
0114 //         theAdjSimManager->GetWeightAtEndOfLastAdjointTrack();
0115 //       normalised_edep = weight_prim*f(ekin_prim);
0116 //     }
0117 //     // then follow the code where normalised_edep is printed, or registered
0118 //     // or whatever ....
0119 //   }
0120 //   else
0121 //   {
0122 //     // code that should be considered only for forward simulation
0123 //   }
0124 //
0125 // Note that in this example a normalisation to only primary e- with only
0126 // one spectrum f(E) is considered. The example code could be easily
0127 // adapted for a normalisation to several spectra and several types of
0128 // primary particles in the same simulation.
0129 
0130 // --------------------------------------------------------------------
0131 //      Class Name:        G4AdjointSimManager
0132 //        Author:               L. Desorgher, 2007-2009
0133 //         Organisation:         SpaceIT GmbH
0134 //        Contract:        ESA contract 21435/08/NL/AT
0135 //         Customer:             ESA/ESTEC
0136 // --------------------------------------------------------------------
0137 #ifndef G4AdjointSimManager_hh
0138 #define G4AdjointSimManager_hh 1
0139 
0140 #include "G4ThreeVector.hh"
0141 #include "G4UserRunAction.hh"
0142 #include "globals.hh"
0143 
0144 #include <vector>
0145 
0146 class G4UserEventAction;
0147 class G4VUserPrimaryGeneratorAction;
0148 class G4UserTrackingAction;
0149 class G4UserSteppingAction;
0150 class G4UserStackingAction;
0151 class G4AdjointRunAction;
0152 class G4AdjointPrimaryGeneratorAction;
0153 class G4AdjointSteppingAction;
0154 class G4AdjointEventAction;
0155 class G4AdjointStackingAction;
0156 class G4AdjointTrackingAction;
0157 class G4ParticleDefinition;
0158 class G4AdjointSimMessenger;
0159 class G4PhysicsLogVector;
0160 class G4Run;
0161 
0162 class G4AdjointSimManager : public G4UserRunAction
0163 {
0164   public:
0165     static G4AdjointSimManager* GetInstance();
0166 
0167     void BeginOfRunAction(const G4Run* aRun) override;
0168     void EndOfRunAction(const G4Run* aRun) override;
0169     void RunAdjointSimulation(G4int nb_evt);
0170 
0171     inline G4int GetNbEvtOfLastRun() { return nb_evt_of_last_run; }
0172 
0173     void SetAdjointTrackingMode(G4bool aBool);
0174     G4bool GetAdjointTrackingMode();
0175     // true if an adjoint track is being processed
0176 
0177     inline G4bool GetAdjointSimMode()
0178     {
0179       return adjoint_sim_mode;
0180     }  // true if an adjoint simulation is running
0181 
0182     G4bool GetDidAdjParticleReachTheExtSource();
0183     void RegisterAtEndOfAdjointTrack();
0184     void RegisterAdjointPrimaryWeight(G4double aWeight);
0185     void ResetDidOneAdjPartReachExtSourceDuringEvent();
0186 
0187     inline G4int GetIDOfLastAdjParticleReachingExtSource()
0188     {
0189       return ID_of_last_particle_that_reach_the_ext_source;
0190     }
0191 
0192     G4ThreeVector GetPositionAtEndOfLastAdjointTrack(std::size_t i = 0);
0193     G4ThreeVector GetDirectionAtEndOfLastAdjointTrack(std::size_t i = 0);
0194     G4double GetEkinAtEndOfLastAdjointTrack(std::size_t i = 0);
0195     G4double GetEkinNucAtEndOfLastAdjointTrack(std::size_t i = 0);
0196     G4double GetWeightAtEndOfLastAdjointTrack(std::size_t i = 0);
0197     G4double GetCosthAtEndOfLastAdjointTrack(std::size_t i = 0);
0198     const G4String& GetFwdParticleNameAtEndOfLastAdjointTrack();
0199     G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(std::size_t i = 0);
0200     G4int GetFwdParticleIndexAtEndOfLastAdjointTrack(std::size_t i = 0);
0201     std::size_t GetNbOfAdointTracksReachingTheExternalSurface();
0202     void ClearEndOfAdjointTrackInfoVectors();
0203     G4ParticleDefinition* GetLastGeneratedFwdPrimaryParticle();
0204 
0205     std::vector<G4ParticleDefinition*>* GetListOfPrimaryFwdParticles();
0206     std::size_t GetNbOfPrimaryFwdParticles();
0207 
0208     G4bool DefineSphericalExtSource(G4double radius, G4ThreeVector pos);
0209     G4bool DefineSphericalExtSourceWithCentreAtTheCentreOfAVolume(G4double radius,
0210                                                                   const G4String& volume_name);
0211     G4bool DefineExtSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name);
0212     void SetExtSourceEmax(G4double Emax);
0213 
0214     // Definition of adjoint source
0215     //----------------------------
0216     G4bool DefineSphericalAdjointSource(G4double radius, G4ThreeVector pos);
0217     G4bool DefineSphericalAdjointSourceWithCentreAtTheCentreOfAVolume(G4double radius,
0218                                                                       const G4String& volume_name);
0219     G4bool DefineAdjointSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name);
0220     void SetAdjointSourceEmin(G4double Emin);
0221     void SetAdjointSourceEmax(G4double Emax);
0222     inline G4double GetAdjointSourceArea() { return area_of_the_adjoint_source; }
0223     void ConsiderParticleAsPrimary(const G4String& particle_name);
0224     void NeglectParticleAsPrimary(const G4String& particle_name);
0225     void SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon);
0226     const G4String& GetPrimaryIonName();
0227 
0228     inline void SetNormalisationMode(G4int n) { normalisation_mode = n; }
0229     inline G4int GetNormalisationMode() { return normalisation_mode; }
0230     inline G4double GetNumberNucleonsInIon() { return nb_nuc; }
0231 
0232     // Definition of user actions for the adjoint tracking phase
0233     //----------------------------
0234     void SetAdjointEventAction(G4UserEventAction* anAction);
0235     void SetAdjointSteppingAction(G4UserSteppingAction* anAction);
0236     void SetAdjointStackingAction(G4UserStackingAction* anAction);
0237     void SetAdjointRunAction(G4UserRunAction* anAction);
0238 
0239     // Set methods for user run actions
0240     //--------------------------------
0241     inline void UseUserStackingActionInFwdTrackingPhase(G4bool aBool)
0242     {
0243       use_user_StackingAction = aBool;
0244     }
0245     inline void UseUserTrackingActionInFwdTrackingPhase(G4bool aBool)
0246     {
0247       use_user_TrackingAction = aBool;
0248     }
0249 
0250     // Set nb of primary fwd gamma
0251     void SetNbOfPrimaryFwdGammasPerEvent(G4int);
0252 
0253     // Set nb of adjoint primaries for reverse splitting
0254     void SetNbAdjointPrimaryGammasPerEvent(G4int);
0255     void SetNbAdjointPrimaryElectronsPerEvent(G4int);
0256 
0257     // Convergence test
0258     //-----------------------
0259     /*
0260      void RegisterSignalForConvergenceTest(G4double aSignal);
0261      void DefineExponentialPrimarySpectrumForConvergenceTest(
0262           G4ParticleDefinition* aPartDef, G4double E0);
0263      void DefinePowerLawPrimarySpectrumForConvergenceTest(
0264           G4ParticleDefinition* aPartDef, G4double alpha);
0265     */
0266 
0267     void SwitchToAdjointSimulationMode();
0268     void BackToFwdSimulationMode();
0269 
0270   private:  // methods
0271     static G4ThreadLocal G4AdjointSimManager* instance;
0272 
0273     void SetRestOfAdjointActions();
0274     void SetAdjointPrimaryRunAndStackingActions();
0275     void SetAdjointActions();
0276     void ResetRestOfUserActions();
0277     void ResetUserPrimaryRunAndStackingActions();
0278     void ResetUserActions();
0279     void DefineUserActions();
0280 
0281     // private constructor and destructor
0282     G4AdjointSimManager();
0283     ~G4AdjointSimManager() override;
0284 
0285   private:  // attributes
0286     // Messenger
0287     //----------
0288     G4AdjointSimMessenger* theMessenger = nullptr;
0289 
0290     // user defined actions for the normal fwd simulation.
0291     // Taken from the G4RunManager
0292     //-------------------------------------------------
0293     G4bool user_action_already_defined = false;
0294     G4UserRunAction* fUserRunAction = nullptr;
0295     G4UserEventAction* fUserEventAction = nullptr;
0296     G4VUserPrimaryGeneratorAction* fUserPrimaryGeneratorAction = nullptr;
0297     G4UserTrackingAction* fUserTrackingAction = nullptr;
0298     G4UserSteppingAction* fUserSteppingAction = nullptr;
0299     G4UserStackingAction* fUserStackingAction = nullptr;
0300     G4bool use_user_StackingAction = false;  // only for fwd part of adjoint sim
0301     G4bool use_user_TrackingAction = false;
0302 
0303     // action for adjoint simulation
0304     //-----------------------------
0305     G4UserRunAction* theAdjointRunAction = nullptr;
0306     G4UserEventAction* theAdjointEventAction = nullptr;
0307     G4AdjointPrimaryGeneratorAction* theAdjointPrimaryGeneratorAction = nullptr;
0308     G4AdjointTrackingAction* theAdjointTrackingAction = nullptr;
0309     G4AdjointSteppingAction* theAdjointSteppingAction = nullptr;
0310     G4AdjointStackingAction* theAdjointStackingAction = nullptr;
0311 
0312     // adjoint mode
0313     //-------------
0314     G4bool adjoint_tracking_mode = false;
0315     G4bool adjoint_sim_mode = false;
0316 
0317     // adjoint particle information on the external surface
0318     //-----------------------------
0319     std::vector<G4ThreeVector> last_pos_vec;
0320     std::vector<G4ThreeVector> last_direction_vec;
0321     std::vector<G4double> last_ekin_vec;
0322     std::vector<G4double> last_ekin_nuc_vec;
0323     std::vector<G4double> last_cos_th_vec;
0324     std::vector<G4double> last_weight_vec;
0325     std::vector<G4int> last_fwd_part_PDGEncoding_vec;
0326     std::vector<G4int> last_fwd_part_index_vec;
0327     std::vector<G4int> ID_of_last_particle_that_reach_the_ext_source_vec;
0328 
0329     G4ThreeVector last_pos;
0330     G4ThreeVector last_direction;
0331     G4double last_ekin = 0.0, last_ekin_nuc = 0.0;
0332     // last_ekin_nuc=last_ekin/nuc, nuc is 1 if not a nucleus
0333     G4double last_cos_th = 0.0;
0334     G4String last_fwd_part_name;
0335     G4int last_fwd_part_PDGEncoding = 0;
0336     G4int last_fwd_part_index = 0;
0337     G4double last_weight = 0.0;
0338     G4int ID_of_last_particle_that_reach_the_ext_source = 0;
0339 
0340     G4int nb_evt_of_last_run = 0;
0341     G4int normalisation_mode = 3;
0342 
0343     // Adjoint source
0344     //--------------
0345     G4double area_of_the_adjoint_source = 0.0;
0346     G4double nb_nuc = 1.0;
0347     G4double theAdjointPrimaryWeight = 0.0;
0348 
0349     G4bool welcome_message = true;
0350 };
0351 
0352 #endif