Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 102566 2017-02-09 08:35:24Z 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"
45 #include "G4Gamma.hh"
46 
48 //
50  : Emin(0.), Emax(0.), EminIon(0.), EmaxIon(0.),
51  index_particle(100000),
52  radius_spherical_source(0.), fwd_ion(0), adj_ion(0),
53  ion_name("not_defined")
54 {
55  theAdjointPrimaryGenerator= new G4AdjointPrimaryGenerator();
56 
57  PrimariesConsideredInAdjointSim[G4String("e-")]=false;
58  PrimariesConsideredInAdjointSim[G4String("gamma")]=false;
59  PrimariesConsideredInAdjointSim[G4String("proton")]=false;
60  PrimariesConsideredInAdjointSim[G4String("ion")]=false;
61 
62  ListOfPrimaryFwdParticles.clear();
63  ListOfPrimaryAdjParticles.clear();
64  nb_fwd_gammas_per_event = 1;
65  nb_adj_primary_gammas_per_event = 1;
66  nb_adj_primary_electrons_per_event = 1;
67 
68 }
70 //
72 {
73  delete theAdjointPrimaryGenerator;
74 }
76 //
78 {
79 
80  G4int evt_id=anEvent->GetEventID();
81  size_t n=ListOfPrimaryAdjParticles.size();
82  index_particle=size_t(evt_id)-n*(size_t(evt_id)/n);
83 
84 
85  G4double E1=Emin;
86  G4double E2=Emax;
87  if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
88 
89  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
90  E1=EminIon;
91  E2=EmaxIon;
92  }
93  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
94  G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
95  E1=EminIon*A;
96  E2=EmaxIon*A;
97  }
98  //Generate first the forwrad primaries
99  theAdjointPrimaryGenerator->GenerateFwdPrimaryVertex(anEvent,ListOfPrimaryFwdParticles[index_particle],E1,E2);
100  G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
101 
102  p=fwdPrimVertex->GetPrimary()->GetMomentum();
103  pos=fwdPrimVertex->GetPosition();
104  G4double pmag=p.mag();
105  G4double m0=ListOfPrimaryFwdParticles[index_particle]->GetPDGMass();
106  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
107 
108  G4double weight_correction=1.;
109  //For gamma generate the particle along the backward ray
110  G4ThreeVector dir=-p/p.mag();
111 
112  /*if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma"){
113 
114  theAdjointPrimaryGenerator
115  ->ComputeAccumulatedDepthVectorAlongBackRay(pos,dir,ekin,ListOfPrimaryAdjParticles[index_particle]);
116 
117 
118  G4double distance = theAdjointPrimaryGenerator
119  ->SampleDistanceAlongBackRayAndComputeWeightCorrection(weight_correction);
120 
121  //pos=pos+dir*distance;
122  //fwdPrimVertex->SetPosition(pos[0],pos[1],pos[2]);
123  }
124  */
125  weight_correction=1.;
126 
127  if (ListOfPrimaryFwdParticles[index_particle] ==G4Gamma::Gamma() && nb_fwd_gammas_per_event>1 ){
128  G4double weight = (1./nb_fwd_gammas_per_event);
129  fwdPrimVertex->SetWeight(weight);
130  for (int i=0;i<nb_fwd_gammas_per_event-1;i++){
131  G4PrimaryVertex* newFwdPrimVertex = new G4PrimaryVertex();
132  newFwdPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
133  newFwdPrimVertex->SetT0(0.);
134  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
135  p.x(),p.y(),p.z());
136  newFwdPrimVertex->SetPrimary(aPrimParticle);
137  newFwdPrimVertex->SetWeight(weight);
138  anEvent->AddPrimaryVertex(newFwdPrimVertex);
139  }
140  }
141 
142  //Now generate the adjoint primaries
143  G4PrimaryVertex* adjPrimVertex = new G4PrimaryVertex();
144  adjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
145  adjPrimVertex->SetT0(0.);
146  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
147  -p.x(),-p.y(),-p.z());
148 
149  adjPrimVertex->SetPrimary(aPrimParticle);
150  anEvent->AddPrimaryVertex(adjPrimVertex);
151 
152  //The factor pi is to normalise the weight to the directional flux
154  G4double adjoint_weight = weight_correction*ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
155  //if (ListOfPrimaryFwdParticles[index_particle] ==G4Gamma::Gamma()) adjoint_weight = adjoint_weight/3.;
156  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma") {
157  //The weight will be corrected at the end of the track if splitted tracks are used
158  adjoint_weight = adjoint_weight/nb_adj_primary_gammas_per_event;
159  for (int i=0;i<nb_adj_primary_gammas_per_event-1;i++){
160  G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
161  newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
162  newAdjPrimVertex->SetT0(0.);
163  aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
164  -p.x(),-p.y(),-p.z());
165  newAdjPrimVertex->SetPrimary(aPrimParticle);
166  newAdjPrimVertex->SetWeight(adjoint_weight);
167  anEvent->AddPrimaryVertex(newAdjPrimVertex);
168  }
169  }
170  else if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_electron") {
171  //The weight will be corrected at the end of the track if splitted tracks are used
172  adjoint_weight = adjoint_weight/nb_adj_primary_electrons_per_event;
173  for (int i=0;i<nb_adj_primary_electrons_per_event-1;i++){
174  G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
175  newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
176  newAdjPrimVertex->SetT0(0.);
177  aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
178  -p.x(),-p.y(),-p.z());
179  newAdjPrimVertex->SetPrimary(aPrimParticle);
180  newAdjPrimVertex->SetWeight(adjoint_weight);
181  anEvent->AddPrimaryVertex(newAdjPrimVertex);
182  }
183  }
184  adjPrimVertex->SetWeight(adjoint_weight);
185 
186  //Call some methods of G4AdjointSimManager
190 
191  /* if ( !last_generated_part_was_adjoint ) {
192 
193  index_particle++;
194  if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
195 
196 
197  G4double E1=Emin;
198  G4double E2=Emax;
199  if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
200 
201  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
202  E1=EminIon;
203  E2=EmaxIon;
204  }
205  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
206  G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
207  E1=EminIon*A;
208  E2=EmaxIon*A;
209  }
210  theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
211  ListOfPrimaryAdjParticles[index_particle],
212  E1,E2);
213  G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
214 
215 
216  p=aPrimVertex->GetPrimary()->GetMomentum();
217  pos=aPrimVertex->GetPosition();
218  G4double pmag=p.mag();
219 
220  G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
221  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
222 
223 
224  //The factor pi is to normalise the weight to the directional flux
225  G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
226  G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
227 
228  aPrimVertex->SetWeight(adjoint_weight);
229 
230  last_generated_part_was_adjoint =true;
231  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
232  G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
233  }
234  else {
235  //fwd particle equivalent to the last generated adjoint particle ios generated
236  G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
237  aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
238  aPrimVertex->SetT0(0.);
239  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
240  -p.x(),-p.y(),-p.z());
241 
242  aPrimVertex->SetPrimary(aPrimParticle);
243  anEvent->AddPrimaryVertex(aPrimVertex);
244  last_generated_part_was_adjoint =false;
245  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
246  */
247 
248 }
250 //
252 {
253  Emin=val;
254  EminIon=val;
255 }
257 //
259 {
260  Emax=val;
261  EmaxIon=val;
262 }
264 //
266 {
267  EminIon=val;
268 }
270 //
272 {
273  EmaxIon=val;
274 }
276 //
277 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
278 {
279  // We generate N numbers of primaries with a 1/E energy law distribution.
280  // We have therefore an energy distribution function
281  // f(E)=C/E (1)
282  // with C a constant that is such that
283  // N=Integral(f(E),E1,E2)=C.std::log(E2/E1) (2)
284  // Therefore from (2) we get
285  // C=N/ std::log(E2/E1) (3)
286  // and
287  // f(E)=N/ std::log(E2/E1)/E (4)
288  //For the adjoint simulation we need a energy distribution f'(E)=1..
289  //To get that we need therefore to apply a weight to the primary
290  // W=1/f(E)=E*std::log(E2/E1)/N
291  //
292  return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
293 
294 }
296 //
298 {
299  radius_spherical_source = radius;
300  center_spherical_source = center_pos;
301  type_of_adjoint_source ="Spherical";
302  theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
303 }
305 //
307 {
308  type_of_adjoint_source ="ExternalSurfaceOfAVolume";
309  theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
310 }
312 //
314 {
315  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
316  PrimariesConsideredInAdjointSim[particle_name]=true;
317  }
319 }
321 //
323 {
324  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
325  PrimariesConsideredInAdjointSim[particle_name]= false;
326  }
328 }
330 //
332 {
334  ListOfPrimaryFwdParticles.clear();
335  ListOfPrimaryAdjParticles.clear();
336  std::map<G4String, G4bool>::iterator iter;
337  for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
338  if(iter->second) {
339  G4String fwd_particle_name = iter->first;
340  if ( fwd_particle_name != "ion") {
341  G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
342  ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
343  ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
344  /*
345  if ( fwd_particle_name == "gamma") {
346  for (G4int i=0;i<2;i++){
347  ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
348  ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
349  }
350  }
351  */
352  }
353  else {
354  if (fwd_ion ){
355  ion_name=fwd_ion->GetParticleName();
356  G4String adj_ion_name=G4String("adj_") +ion_name;
357  ListOfPrimaryFwdParticles.push_back(fwd_ion);
358  ListOfPrimaryAdjParticles.push_back(adj_ion);
359  }
360  else {
361  ListOfPrimaryFwdParticles.push_back(0);
362  ListOfPrimaryAdjParticles.push_back(0);
363 
364  }
365  }
366  }
367  }
368 }
370 //
372 {
373  fwd_ion = fwdIon;
374  adj_ion = adjointIon;
376 }
377 
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:154
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:151
void SetPrimaryIon(G4ParticleDefinition *adjointIon, G4ParticleDefinition *fwdIon)
G4PrimaryParticle * GetPrimary(G4int i=0) const
void ConsiderParticleAsPrimary(const G4String &particle_name)
double A(double temperature)
void SetAdjointTrackingMode(G4bool aBool)
const G4int n
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition: G4Event.hh:167
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
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
static const G4double Emin
void SetPrimary(G4PrimaryParticle *pp)
void SetPosition(G4double x0, G4double y0, G4double z0)
static const G4double Emax
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
void SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector pos)
double mag() const
void ResetDidOneAdjPartReachExtSourceDuringEvent()