Geant4_10
G4AdjointPrimaryGeneratorAction.cc
Go to the documentation of this file.
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 // $Id: G4AdjointPrimaryGeneratorAction.cc 76245 2013-11-08 11:14:32Z gcosmo $
27 //
29 // Class Name: G4AdjointPrimaryGeneratorAction
30 // Author: L. Desorgher
31 // Organisation: SpaceIT GmbH
32 // Contract: ESA contract 21435/08/NL/AT
33 // Customer: ESA/ESTEC
35 
37 
38 #include "G4PhysicalConstants.hh"
39 #include "G4Event.hh"
40 #include "G4ParticleTable.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4AdjointSimManager.hh"
46 //
48  : Emin(0.), Emax(0.), EminIon(0.), EmaxIon(0.),
49  index_particle(100000),
50  radius_spherical_source(0.), fwd_ion(0), adj_ion(0),
51  ion_name("not_defined")
52 {
53  theAdjointPrimaryGenerator= new G4AdjointPrimaryGenerator();
54 
55  PrimariesConsideredInAdjointSim[G4String("e-")]=false;
56  PrimariesConsideredInAdjointSim[G4String("gamma")]=false;
57  PrimariesConsideredInAdjointSim[G4String("proton")]=false;
58  PrimariesConsideredInAdjointSim[G4String("ion")]=false;
59 
60  ListOfPrimaryFwdParticles.clear();
61  ListOfPrimaryAdjParticles.clear();
62 }
64 //
66 {
67  delete theAdjointPrimaryGenerator;
68 }
70 //
72 { G4int evt_id=anEvent->GetEventID();
73  size_t n=ListOfPrimaryAdjParticles.size();
74  index_particle=size_t(evt_id)-n*(size_t(evt_id)/n);
75 
76 
77  G4double E1=Emin;
78  G4double E2=Emax;
79  if (!ListOfPrimaryAdjParticles[index_particle]) 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  theAdjointPrimaryGenerator->GenerateFwdPrimaryVertex(anEvent,ListOfPrimaryFwdParticles[index_particle],E1,E2);
91  G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
92 
93  p=fwdPrimVertex->GetPrimary()->GetMomentum();
94  pos=fwdPrimVertex->GetPosition();
95  G4double pmag=p.mag();
96 
97 
98  G4double m0=ListOfPrimaryFwdParticles[index_particle]->GetPDGMass();
99  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
100 
101  G4PrimaryVertex* adjPrimVertex = new G4PrimaryVertex();
102  adjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
103  adjPrimVertex->SetT0(0.);
104  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
105  -p.x(),-p.y(),-p.z());
106 
107  adjPrimVertex->SetPrimary(aPrimParticle);
108  anEvent->AddPrimaryVertex(adjPrimVertex);
109 
110  //The factor pi is to normalise the weight to the directional flux
112  G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
113  adjPrimVertex->SetWeight(adjoint_weight);
114 
115 
116 
118 
119 
120 
121  /* if ( !last_generated_part_was_adjoint ) {
122 
123  index_particle++;
124  if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
125 
126 
127  G4double E1=Emin;
128  G4double E2=Emax;
129  if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
130 
131  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
132  E1=EminIon;
133  E2=EmaxIon;
134  }
135  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
136  G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
137  E1=EminIon*A;
138  E2=EmaxIon*A;
139  }
140  theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
141  ListOfPrimaryAdjParticles[index_particle],
142  E1,E2);
143  G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
144 
145 
146  p=aPrimVertex->GetPrimary()->GetMomentum();
147  pos=aPrimVertex->GetPosition();
148  G4double pmag=p.mag();
149 
150  G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
151  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
152 
153 
154  //The factor pi is to normalise the weight to the directional flux
155  G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
156  G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
157 
158  aPrimVertex->SetWeight(adjoint_weight);
159 
160  last_generated_part_was_adjoint =true;
161  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
162  G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
163  }
164  else {
165  //fwd particle equivalent to the last generated adjoint particle ios generated
166  G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
167  aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
168  aPrimVertex->SetT0(0.);
169  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
170  -p.x(),-p.y(),-p.z());
171 
172  aPrimVertex->SetPrimary(aPrimParticle);
173  anEvent->AddPrimaryVertex(aPrimVertex);
174  last_generated_part_was_adjoint =false;
175  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
176  */
177 
178 }
180 //
182 {
183  Emin=val;
184  EminIon=val;
185 }
187 //
189 {
190  Emax=val;
191  EmaxIon=val;
192 }
194 //
196 {
197  EminIon=val;
198 }
200 //
202 {
203  EmaxIon=val;
204 }
206 //
207 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
208 {
209  // We generate N numbers of primaries with a 1/E energy law distribution.
210  // We have therefore an energy distribution function
211  // f(E)=C/E (1)
212  // with C a constant that is such that
213  // N=Integral(f(E),E1,E2)=C.std::log(E2/E1) (2)
214  // Therefore from (2) we get
215  // C=N/ std::log(E2/E1) (3)
216  // and
217  // f(E)=N/ std::log(E2/E1)/E (4)
218  //For the adjoint simulation we need a energy distribution f'(E)=1..
219  //To get that we need therefore to apply a weight to the primary
220  // W=1/f(E)=E*std::log(E2/E1)/N
221  //
222  return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
223 
224 }
226 //
228 {
229  radius_spherical_source = radius;
230  center_spherical_source = center_pos;
231  type_of_adjoint_source ="Spherical";
232  theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
233 }
235 //
237 {
238  type_of_adjoint_source ="ExternalSurfaceOfAVolume";
239  theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
240 }
242 //
244 {
245  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
246  PrimariesConsideredInAdjointSim[particle_name]=true;
247  }
249 }
251 //
253 {
254  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
255  PrimariesConsideredInAdjointSim[particle_name]= false;
256  }
258 }
260 //
262 {
264  ListOfPrimaryFwdParticles.clear();
265  ListOfPrimaryAdjParticles.clear();
266  std::map<G4String, G4bool>::iterator iter;
267  for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
268  if(iter->second) {
269  G4String fwd_particle_name = iter->first;
270  if ( fwd_particle_name != "ion") {
271  G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
272  ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
273  ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
274  }
275  else {
276  if (fwd_ion ){
277  ion_name=fwd_ion->GetParticleName();
278  G4String adj_ion_name=G4String("adj_") +ion_name;
279  ListOfPrimaryFwdParticles.push_back(fwd_ion);
280  ListOfPrimaryAdjParticles.push_back(adj_ion);
281  }
282  else {
283  ListOfPrimaryFwdParticles.push_back(0);
284  ListOfPrimaryAdjParticles.push_back(0);
285 
286  }
287  }
288  }
289  }
290 }
292 //
294 {
295  fwd_ion = fwdIon;
296  adj_ion = adjointIon;
298 }
299 
static G4AdjointSimManager * GetInstance()
void SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector pos)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4ThreeVector GetMomentum() const
G4ThreeVector GetPosition() const
double x() const
G4int first(char) const
void SetWeight(G4double w)
const char * p
Definition: xmltok.h:285
void AddPrimaryVertex(G4PrimaryVertex *aPrimaryVertex)
Definition: G4Event.hh:143
void SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String &volume_name)
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
double z() const
G4int GetEventID() const
Definition: G4Event.hh:140
void SetPrimaryIon(G4ParticleDefinition *adjointIon, G4ParticleDefinition *fwdIon)
G4PrimaryParticle * GetPrimary(G4int i=0) const
Char_t n[5]
void ConsiderParticleAsPrimary(const G4String &particle_name)
void SetAdjointTrackingMode(G4bool aBool)
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition: G4Event.hh:156
void GenerateFwdPrimaryVertex(G4Event *anEvt, G4ParticleDefinition *adj_part, G4double E1, G4double E2)
static G4ParticleTable * GetParticleTable()
void NeglectParticleAsPrimary(const G4String &particle_name)
void SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String &volume_name)
void SetT0(G4double t0)
double y() const
void SetPrimary(G4PrimaryParticle *pp)
void SetPosition(G4double x0, G4double y0, G4double z0)
double G4double
Definition: G4Types.hh:76
void SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector pos)
double mag() const