Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/run/src/G4AdjointPrimaryGeneratorAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4AdjointPrimaryGeneratorAction implementation
 27 //
 28 // --------------------------------------------------------------------
 29 //   Class Name:   G4AdjointPrimaryGeneratorAction
 30 //   Author:       L. Desorgher, 2007-2009
 31 //   Organisation: SpaceIT GmbH
 32 //   Contract:     ESA contract 21435/08/NL/AT
 33 //   Customer:     ESA/ESTEC
 34 // --------------------------------------------------------------------
 35 
 36 #include "G4AdjointPrimaryGeneratorAction.hh"
 37 
 38 #include "G4AdjointPrimaryGenerator.hh"
 39 #include "G4AdjointSimManager.hh"
 40 #include "G4Event.hh"
 41 #include "G4Gamma.hh"
 42 #include "G4ParticleDefinition.hh"
 43 #include "G4ParticleTable.hh"
 44 #include "G4PhysicalConstants.hh"
 45 
 46 // --------------------------------------------------------------------
 47 //
 48 G4AdjointPrimaryGeneratorAction::G4AdjointPrimaryGeneratorAction()
 49 {
 50   theAdjointPrimaryGenerator = new G4AdjointPrimaryGenerator();
 51 
 52   PrimariesConsideredInAdjointSim[G4String("e-")] = false;
 53   PrimariesConsideredInAdjointSim[G4String("gamma")] = false;
 54   PrimariesConsideredInAdjointSim[G4String("proton")] = false;
 55   PrimariesConsideredInAdjointSim[G4String("ion")] = false;
 56 
 57   ListOfPrimaryFwdParticles.clear();
 58   ListOfPrimaryAdjParticles.clear();
 59 }
 60 
 61 // --------------------------------------------------------------------
 62 //
 63 G4AdjointPrimaryGeneratorAction::~G4AdjointPrimaryGeneratorAction()
 64 {
 65   delete theAdjointPrimaryGenerator;
 66 }
 67 
 68 // --------------------------------------------------------------------
 69 //
 70 void G4AdjointPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 71 {
 72   G4int evt_id = anEvent->GetEventID();
 73   std::size_t n = ListOfPrimaryAdjParticles.size();
 74   index_particle = std::size_t(evt_id) - n * (std::size_t(evt_id) / n);
 75 
 76   G4double E1 = Emin;
 77   G4double E2 = Emax;
 78   if (ListOfPrimaryAdjParticles[index_particle] == nullptr)
 79     UpdateListOfPrimaryParticles();  // ion has not been created yet
 80 
 81   if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
 82     E1 = EminIon;
 83     E2 = EmaxIon;
 84   }
 85   if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
 86     G4int A = ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
 87     E1 = EminIon * A;
 88     E2 = EmaxIon * A;
 89   }
 90   // Generate first the forwrad primaries
 91   theAdjointPrimaryGenerator->GenerateFwdPrimaryVertex(
 92     anEvent, ListOfPrimaryFwdParticles[index_particle], E1, E2);
 93   G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
 94 
 95   p = fwdPrimVertex->GetPrimary()->GetMomentum();
 96   pos = fwdPrimVertex->GetPosition();
 97   G4double pmag = p.mag();
 98   G4double m0 = ListOfPrimaryFwdParticles[index_particle]->GetPDGMass();
 99   G4double ekin = std::sqrt(m0 * m0 + pmag * pmag) - m0;
100 
101   G4double weight_correction = 1.;
102   // For gamma generate the particle along the backward ray
103   G4ThreeVector dir = -p / p.mag();
104 
105   weight_correction = 1.;
106 
107   if (ListOfPrimaryFwdParticles[index_particle] == G4Gamma::Gamma() && nb_fwd_gammas_per_event > 1)
108   {
109     G4double weight = (1. / nb_fwd_gammas_per_event);
110     fwdPrimVertex->SetWeight(weight);
111     for (G4int i = 0; i < nb_fwd_gammas_per_event - 1; ++i) {
112       auto newFwdPrimVertex = new G4PrimaryVertex();
113       newFwdPrimVertex->SetPosition(pos.x(), pos.y(), pos.z());
114       newFwdPrimVertex->SetT0(0.);
115       auto aPrimParticle =
116         new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle], p.x(), p.y(), p.z());
117       newFwdPrimVertex->SetPrimary(aPrimParticle);
118       newFwdPrimVertex->SetWeight(weight);
119       anEvent->AddPrimaryVertex(newFwdPrimVertex);
120     }
121   }
122 
123   // Now generate the adjoint primaries
124   auto adjPrimVertex = new G4PrimaryVertex();
125   adjPrimVertex->SetPosition(pos.x(), pos.y(), pos.z());
126   adjPrimVertex->SetT0(0.);
127   auto aPrimParticle =
128     new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle], -p.x(), -p.y(), -p.z());
129 
130   adjPrimVertex->SetPrimary(aPrimParticle);
131   anEvent->AddPrimaryVertex(adjPrimVertex);
132 
133   // The factor pi is to normalise the weight to the directional flux
134   G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
135   G4double adjoint_weight =
136     weight_correction * ComputeEnergyDistWeight(ekin, E1, E2) * adjoint_source_area * pi;
137   if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma") {
138     // The weight will be corrected at the end of the track if splitted tracks
139     // are used
140     adjoint_weight = adjoint_weight / nb_adj_primary_gammas_per_event;
141     for (G4int i = 0; i < nb_adj_primary_gammas_per_event - 1; ++i) {
142       auto newAdjPrimVertex = new G4PrimaryVertex();
143       newAdjPrimVertex->SetPosition(pos.x(), pos.y(), pos.z());
144       newAdjPrimVertex->SetT0(0.);
145       aPrimParticle =
146         new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle], -p.x(), -p.y(), -p.z());
147       newAdjPrimVertex->SetPrimary(aPrimParticle);
148       newAdjPrimVertex->SetWeight(adjoint_weight);
149       anEvent->AddPrimaryVertex(newAdjPrimVertex);
150     }
151   }
152   else if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_electron") {
153     // The weight will be corrected at the end of the track if splitted tracks
154     // are used
155     adjoint_weight = adjoint_weight / nb_adj_primary_electrons_per_event;
156     for (G4int i = 0; i < nb_adj_primary_electrons_per_event - 1; ++i) {
157       auto newAdjPrimVertex = new G4PrimaryVertex();
158       newAdjPrimVertex->SetPosition(pos.x(), pos.y(), pos.z());
159       newAdjPrimVertex->SetT0(0.);
160       aPrimParticle =
161         new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle], -p.x(), -p.y(), -p.z());
162       newAdjPrimVertex->SetPrimary(aPrimParticle);
163       newAdjPrimVertex->SetWeight(adjoint_weight);
164       anEvent->AddPrimaryVertex(newAdjPrimVertex);
165     }
166   }
167   adjPrimVertex->SetWeight(adjoint_weight);
168 
169   // Call some methods of G4AdjointSimManager
170   G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
171   G4AdjointSimManager::GetInstance()->ClearEndOfAdjointTrackInfoVectors();
172   G4AdjointSimManager::GetInstance()->ResetDidOneAdjPartReachExtSourceDuringEvent();
173 }
174 
175 // --------------------------------------------------------------------
176 //
177 void G4AdjointPrimaryGeneratorAction::SetEmin(G4double val)
178 {
179   Emin = val;
180   EminIon = val;
181 }
182 
183 // --------------------------------------------------------------------
184 //
185 void G4AdjointPrimaryGeneratorAction::SetEmax(G4double val)
186 {
187   Emax = val;
188   EmaxIon = val;
189 }
190 
191 // --------------------------------------------------------------------
192 //
193 void G4AdjointPrimaryGeneratorAction::SetEminIon(G4double val)
194 {
195   EminIon = val;
196 }
197 
198 // --------------------------------------------------------------------
199 //
200 void G4AdjointPrimaryGeneratorAction::SetEmaxIon(G4double val)
201 {
202   EmaxIon = val;
203 }
204 
205 // --------------------------------------------------------------------
206 //
207 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E, G4double E1,
208                                                                   G4double E2)
209 {
210   // We generate N numbers of primaries with  a 1/E energy law distribution.
211   // We have therefore  an energy distribution function
212   //   f(E)=C/E  (1)
213   // with   C a constant that is such that
214   //   N=Integral(f(E),E1,E2)=C.std::log(E2/E1)  (2)
215   // Therefore from (2) we get
216   //   C=N/ std::log(E2/E1) (3)
217   // and
218   //   f(E)=N/ std::log(E2/E1)/E (4)
219   // For the adjoint simulation we need a energy distribution  f'(E)=1..
220   // To get that we need therefore to apply a weight to the primary
221   //   W=1/f(E)=E*std::log(E2/E1)/N
222   //
223   return std::log(E2 / E1) * E / G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
224 }
225 
226 // --------------------------------------------------------------------
227 //
228 void G4AdjointPrimaryGeneratorAction::SetSphericalAdjointPrimarySource(G4double radius,
229                                                                        G4ThreeVector center_pos)
230 {
231   radius_spherical_source = radius;
232   center_spherical_source = center_pos;
233   type_of_adjoint_source = "Spherical";
234   theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius, center_pos);
235 }
236 
237 // --------------------------------------------------------------------
238 //
239 void G4AdjointPrimaryGeneratorAction::SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(
240   const G4String& volume_name)
241 {
242   type_of_adjoint_source = "ExternalSurfaceOfAVolume";
243   theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
244 }
245 
246 // --------------------------------------------------------------------
247 //
248 void G4AdjointPrimaryGeneratorAction::ConsiderParticleAsPrimary(const G4String& particle_name)
249 {
250   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end())
251   {
252     PrimariesConsideredInAdjointSim[particle_name] = true;
253   }
254   UpdateListOfPrimaryParticles();
255 }
256 
257 // --------------------------------------------------------------------
258 //
259 void G4AdjointPrimaryGeneratorAction::NeglectParticleAsPrimary(const G4String& particle_name)
260 {
261   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end())
262   {
263     PrimariesConsideredInAdjointSim[particle_name] = false;
264   }
265   UpdateListOfPrimaryParticles();
266 }
267 
268 // --------------------------------------------------------------------
269 //
270 void G4AdjointPrimaryGeneratorAction::UpdateListOfPrimaryParticles()
271 {
272   G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
273   ListOfPrimaryFwdParticles.clear();
274   ListOfPrimaryAdjParticles.clear();
275   for (const auto& iter : PrimariesConsideredInAdjointSim) {
276     if (iter.second) {
277       G4String fwd_particle_name = iter.first;
278       if (fwd_particle_name != "ion") {
279         G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
280         ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
281         ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
282       }
283       else {
284         if (fwd_ion != nullptr) {
285           ion_name = fwd_ion->GetParticleName();
286           G4String adj_ion_name = G4String("adj_") + ion_name;
287           ListOfPrimaryFwdParticles.push_back(fwd_ion);
288           ListOfPrimaryAdjParticles.push_back(adj_ion);
289         }
290         else {
291           ListOfPrimaryFwdParticles.push_back(nullptr);
292           ListOfPrimaryAdjParticles.push_back(nullptr);
293         }
294       }
295     }
296   }
297 }
298 
299 // --------------------------------------------------------------------
300 //
301 void G4AdjointPrimaryGeneratorAction::SetPrimaryIon(G4ParticleDefinition* adjointIon,
302                                                     G4ParticleDefinition* fwdIon)
303 {
304   fwd_ion = fwdIon;
305   adj_ion = adjointIon;
306   UpdateListOfPrimaryParticles();
307 }
308