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 ]

Diff markup

Differences between /run/src/G4AdjointPrimaryGeneratorAction.cc (Version 11.3.0) and /run/src/G4AdjointPrimaryGeneratorAction.cc (Version 5.0)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // G4AdjointPrimaryGeneratorAction implementat    
 27 //                                                
 28 // -------------------------------------------    
 29 //   Class Name:   G4AdjointPrimaryGeneratorAc    
 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::G4AdjointPrim    
 49 {                                                 
 50   theAdjointPrimaryGenerator = new G4AdjointPr    
 51                                                   
 52   PrimariesConsideredInAdjointSim[G4String("e-    
 53   PrimariesConsideredInAdjointSim[G4String("ga    
 54   PrimariesConsideredInAdjointSim[G4String("pr    
 55   PrimariesConsideredInAdjointSim[G4String("io    
 56                                                   
 57   ListOfPrimaryFwdParticles.clear();              
 58   ListOfPrimaryAdjParticles.clear();              
 59 }                                                 
 60                                                   
 61 // -------------------------------------------    
 62 //                                                
 63 G4AdjointPrimaryGeneratorAction::~G4AdjointPri    
 64 {                                                 
 65   delete theAdjointPrimaryGenerator;              
 66 }                                                 
 67                                                   
 68 // -------------------------------------------    
 69 //                                                
 70 void G4AdjointPrimaryGeneratorAction::Generate    
 71 {                                                 
 72   G4int evt_id = anEvent->GetEventID();           
 73   std::size_t n = ListOfPrimaryAdjParticles.si    
 74   index_particle = std::size_t(evt_id) - n * (    
 75                                                   
 76   G4double E1 = Emin;                             
 77   G4double E2 = Emax;                             
 78   if (ListOfPrimaryAdjParticles[index_particle    
 79     UpdateListOfPrimaryParticles();  // ion ha    
 80                                                   
 81   if (ListOfPrimaryAdjParticles[index_particle    
 82     E1 = EminIon;                                 
 83     E2 = EmaxIon;                                 
 84   }                                               
 85   if (ListOfPrimaryAdjParticles[index_particle    
 86     G4int A = ListOfPrimaryAdjParticles[index_    
 87     E1 = EminIon * A;                             
 88     E2 = EmaxIon * A;                             
 89   }                                               
 90   // Generate first the forwrad primaries         
 91   theAdjointPrimaryGenerator->GenerateFwdPrima    
 92     anEvent, ListOfPrimaryFwdParticles[index_p    
 93   G4PrimaryVertex* fwdPrimVertex = anEvent->Ge    
 94                                                   
 95   p = fwdPrimVertex->GetPrimary()->GetMomentum    
 96   pos = fwdPrimVertex->GetPosition();             
 97   G4double pmag = p.mag();                        
 98   G4double m0 = ListOfPrimaryFwdParticles[inde    
 99   G4double ekin = std::sqrt(m0 * m0 + pmag * p    
100                                                   
101   G4double weight_correction = 1.;                
102   // For gamma generate the particle along the    
103   G4ThreeVector dir = -p / p.mag();               
104                                                   
105   weight_correction = 1.;                         
106                                                   
107   if (ListOfPrimaryFwdParticles[index_particle    
108   {                                               
109     G4double weight = (1. / nb_fwd_gammas_per_    
110     fwdPrimVertex->SetWeight(weight);             
111     for (G4int i = 0; i < nb_fwd_gammas_per_ev    
112       auto newFwdPrimVertex = new G4PrimaryVer    
113       newFwdPrimVertex->SetPosition(pos.x(), p    
114       newFwdPrimVertex->SetT0(0.);                
115       auto aPrimParticle =                        
116         new G4PrimaryParticle(ListOfPrimaryFwd    
117       newFwdPrimVertex->SetPrimary(aPrimPartic    
118       newFwdPrimVertex->SetWeight(weight);        
119       anEvent->AddPrimaryVertex(newFwdPrimVert    
120     }                                             
121   }                                               
122                                                   
123   // Now generate the adjoint primaries           
124   auto adjPrimVertex = new G4PrimaryVertex();     
125   adjPrimVertex->SetPosition(pos.x(), pos.y(),    
126   adjPrimVertex->SetT0(0.);                       
127   auto aPrimParticle =                            
128     new G4PrimaryParticle(ListOfPrimaryAdjPart    
129                                                   
130   adjPrimVertex->SetPrimary(aPrimParticle);       
131   anEvent->AddPrimaryVertex(adjPrimVertex);       
132                                                   
133   // The factor pi is to normalise the weight     
134   G4double adjoint_source_area = G4AdjointSimM    
135   G4double adjoint_weight =                       
136     weight_correction * ComputeEnergyDistWeigh    
137   if (ListOfPrimaryAdjParticles[index_particle    
138     // The weight will be corrected at the end    
139     // are used                                   
140     adjoint_weight = adjoint_weight / nb_adj_p    
141     for (G4int i = 0; i < nb_adj_primary_gamma    
142       auto newAdjPrimVertex = new G4PrimaryVer    
143       newAdjPrimVertex->SetPosition(pos.x(), p    
144       newAdjPrimVertex->SetT0(0.);                
145       aPrimParticle =                             
146         new G4PrimaryParticle(ListOfPrimaryAdj    
147       newAdjPrimVertex->SetPrimary(aPrimPartic    
148       newAdjPrimVertex->SetWeight(adjoint_weig    
149       anEvent->AddPrimaryVertex(newAdjPrimVert    
150     }                                             
151   }                                               
152   else if (ListOfPrimaryAdjParticles[index_par    
153     // The weight will be corrected at the end    
154     // are used                                   
155     adjoint_weight = adjoint_weight / nb_adj_p    
156     for (G4int i = 0; i < nb_adj_primary_elect    
157       auto newAdjPrimVertex = new G4PrimaryVer    
158       newAdjPrimVertex->SetPosition(pos.x(), p    
159       newAdjPrimVertex->SetT0(0.);                
160       aPrimParticle =                             
161         new G4PrimaryParticle(ListOfPrimaryAdj    
162       newAdjPrimVertex->SetPrimary(aPrimPartic    
163       newAdjPrimVertex->SetWeight(adjoint_weig    
164       anEvent->AddPrimaryVertex(newAdjPrimVert    
165     }                                             
166   }                                               
167   adjPrimVertex->SetWeight(adjoint_weight);       
168                                                   
169   // Call some methods of G4AdjointSimManager     
170   G4AdjointSimManager::GetInstance()->SetAdjoi    
171   G4AdjointSimManager::GetInstance()->ClearEnd    
172   G4AdjointSimManager::GetInstance()->ResetDid    
173 }                                                 
174                                                   
175 // -------------------------------------------    
176 //                                                
177 void G4AdjointPrimaryGeneratorAction::SetEmin(    
178 {                                                 
179   Emin = val;                                     
180   EminIon = val;                                  
181 }                                                 
182                                                   
183 // -------------------------------------------    
184 //                                                
185 void G4AdjointPrimaryGeneratorAction::SetEmax(    
186 {                                                 
187   Emax = val;                                     
188   EmaxIon = val;                                  
189 }                                                 
190                                                   
191 // -------------------------------------------    
192 //                                                
193 void G4AdjointPrimaryGeneratorAction::SetEminI    
194 {                                                 
195   EminIon = val;                                  
196 }                                                 
197                                                   
198 // -------------------------------------------    
199 //                                                
200 void G4AdjointPrimaryGeneratorAction::SetEmaxI    
201 {                                                 
202   EmaxIon = val;                                  
203 }                                                 
204                                                   
205 // -------------------------------------------    
206 //                                                
207 G4double G4AdjointPrimaryGeneratorAction::Comp    
208                                                   
209 {                                                 
210   // We generate N numbers of primaries with      
211   // We have therefore  an energy distribution    
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    
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 ener    
220   // To get that we need therefore to apply a     
221   //   W=1/f(E)=E*std::log(E2/E1)/N               
222   //                                              
223   return std::log(E2 / E1) * E / G4AdjointSimM    
224 }                                                 
225                                                   
226 // -------------------------------------------    
227 //                                                
228 void G4AdjointPrimaryGeneratorAction::SetSpher    
229                                                   
230 {                                                 
231   radius_spherical_source = radius;               
232   center_spherical_source = center_pos;           
233   type_of_adjoint_source = "Spherical";           
234   theAdjointPrimaryGenerator->SetSphericalAdjo    
235 }                                                 
236                                                   
237 // -------------------------------------------    
238 //                                                
239 void G4AdjointPrimaryGeneratorAction::SetAdjoi    
240   const G4String& volume_name)                    
241 {                                                 
242   type_of_adjoint_source = "ExternalSurfaceOfA    
243   theAdjointPrimaryGenerator->SetAdjointPrimar    
244 }                                                 
245                                                   
246 // -------------------------------------------    
247 //                                                
248 void G4AdjointPrimaryGeneratorAction::Consider    
249 {                                                 
250   if (PrimariesConsideredInAdjointSim.find(par    
251   {                                               
252     PrimariesConsideredInAdjointSim[particle_n    
253   }                                               
254   UpdateListOfPrimaryParticles();                 
255 }                                                 
256                                                   
257 // -------------------------------------------    
258 //                                                
259 void G4AdjointPrimaryGeneratorAction::NeglectP    
260 {                                                 
261   if (PrimariesConsideredInAdjointSim.find(par    
262   {                                               
263     PrimariesConsideredInAdjointSim[particle_n    
264   }                                               
265   UpdateListOfPrimaryParticles();                 
266 }                                                 
267                                                   
268 // -------------------------------------------    
269 //                                                
270 void G4AdjointPrimaryGeneratorAction::UpdateLi    
271 {                                                 
272   G4ParticleTable* theParticleTable = G4Partic    
273   ListOfPrimaryFwdParticles.clear();              
274   ListOfPrimaryAdjParticles.clear();              
275   for (const auto& iter : PrimariesConsideredI    
276     if (iter.second) {                            
277       G4String fwd_particle_name = iter.first;    
278       if (fwd_particle_name != "ion") {           
279         G4String adj_particle_name = G4String(    
280         ListOfPrimaryFwdParticles.push_back(th    
281         ListOfPrimaryAdjParticles.push_back(th    
282       }                                           
283       else {                                      
284         if (fwd_ion != nullptr) {                 
285           ion_name = fwd_ion->GetParticleName(    
286           G4String adj_ion_name = G4String("ad    
287           ListOfPrimaryFwdParticles.push_back(    
288           ListOfPrimaryAdjParticles.push_back(    
289         }                                         
290         else {                                    
291           ListOfPrimaryFwdParticles.push_back(    
292           ListOfPrimaryAdjParticles.push_back(    
293         }                                         
294       }                                           
295     }                                             
296   }                                               
297 }                                                 
298                                                   
299 // -------------------------------------------    
300 //                                                
301 void G4AdjointPrimaryGeneratorAction::SetPrima    
302                                                   
303 {                                                 
304   fwd_ion = fwdIon;                               
305   adj_ion = adjointIon;                           
306   UpdateListOfPrimaryParticles();                 
307 }                                                 
308