source: trunk/source/run/src/G4AdjointPrimaryGeneratorAction.cc@ 1261

Last change on this file since 1261 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 10.3 KB
Line 
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,v 1.2 2009/11/18 18:02:06 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29/////////////////////////////////////////////////////////////////////////////
30// Class Name: G4AdjointCrossSurfChecker
31// Author: L. Desorgher
32// Organisation: SpaceIT GmbH
33// Contract: ESA contract 21435/08/NL/AT
34// Customer: ESA/ESTEC
35/////////////////////////////////////////////////////////////////////////////
36
37#include "G4AdjointPrimaryGeneratorAction.hh"
38#include "G4Event.hh"
39#include "G4ParticleTable.hh"
40#include "G4ParticleDefinition.hh"
41#include "G4ParticleTable.hh"
42#include "G4AdjointSimManager.hh"
43#include "G4AdjointPrimaryGenerator.hh"
44/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
45//
46G4AdjointPrimaryGeneratorAction::G4AdjointPrimaryGeneratorAction()
47{
48 theAdjointPrimaryGenerator= new G4AdjointPrimaryGenerator();
49
50 PrimariesConsideredInAdjointSim[G4String("e-")]=false;
51 PrimariesConsideredInAdjointSim[G4String("gamma")]=false;
52 PrimariesConsideredInAdjointSim[G4String("proton")]=false;
53 PrimariesConsideredInAdjointSim[G4String("ion")]=false;
54
55 ListOfPrimaryFwdParticles.clear();
56 ListOfPrimaryAdjParticles.clear();
57
58 last_generated_part_was_adjoint=false;
59 index_particle=100000;
60
61 ion_name="not_defined";
62 fwd_ion = 0;
63 adj_ion = 0;
64}
65/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
66//
67G4AdjointPrimaryGeneratorAction::~G4AdjointPrimaryGeneratorAction()
68{
69 delete theAdjointPrimaryGenerator;
70}
71/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
72//
73void G4AdjointPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
74{
75 if ( !last_generated_part_was_adjoint ) {
76
77 index_particle++;
78 if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
79
80
81 G4double E1=Emin;
82 G4double E2=Emax;
83 if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
84
85 if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
86 E1=EminIon;
87 E2=EmaxIon;
88 }
89 if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
90 G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
91 E1=EminIon*A;
92 E2=EmaxIon*A;
93 }
94 theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
95 ListOfPrimaryAdjParticles[index_particle],
96 E1,E2);
97 G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
98
99
100 p=aPrimVertex->GetPrimary()->GetMomentum();
101 pos=aPrimVertex->GetPosition();
102 G4double pmag=p.mag();
103
104 G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
105 G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
106
107
108 //The factor pi is to normalise the weight to the directional flux
109 G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
110 G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
111
112 aPrimVertex->SetWeight(adjoint_weight);
113
114 last_generated_part_was_adjoint =true;
115 G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
116 G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
117 }
118 else {
119 //fwd particle equivalent to the last generated adjoint particle ios generated
120 G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
121 aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
122 aPrimVertex->SetT0(0.);
123 G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
124 -p.x(),-p.y(),-p.z());
125
126 aPrimVertex->SetPrimary(aPrimParticle);
127 anEvent->AddPrimaryVertex(aPrimVertex);
128 last_generated_part_was_adjoint =false;
129 G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
130 }
131}
132/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
133//
134void G4AdjointPrimaryGeneratorAction::SetEmin(G4double val)
135{
136 Emin=val;
137 EminIon=val;
138}
139/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
140//
141void G4AdjointPrimaryGeneratorAction::SetEmax(G4double val)
142{
143 Emax=val;
144 EmaxIon=val;
145}
146/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
147//
148void G4AdjointPrimaryGeneratorAction::SetEminIon(G4double val)
149{
150 EminIon=val;
151}
152/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
153//
154void G4AdjointPrimaryGeneratorAction::SetEmaxIon(G4double val)
155{
156 EmaxIon=val;
157}
158/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
159//
160G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
161{
162 // We generate N numbers of primaries with a 1/E energy law distribution.
163 // We have therefore an energy distribution function
164 // f(E)=C/E (1)
165 // with C a constant that is such that
166 // N=Integral(f(E),E1,E2)=C.std::log(E2/E1) (2)
167 // Therefore from (2) we get
168 // C=N/ std::log(E2/E1) (3)
169 // and
170 // f(E)=N/ std::log(E2/E1)/E (4)
171 //For the adjoint simulation we need a energy distribution f'(E)=1..
172 //To get that we need therefore to apply a weight to the primary
173 // W=1/f(E)=E*std::log(E2/E1)/N
174 //
175 return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
176
177}
178/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
179//
180void G4AdjointPrimaryGeneratorAction::SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector center_pos)
181{
182 radius_spherical_source = radius;
183 center_spherical_source = center_pos;
184 type_of_adjoint_source ="Spherical";
185 theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
186}
187/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
188//
189void G4AdjointPrimaryGeneratorAction::SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String& volume_name)
190{
191 type_of_adjoint_source ="ExternalSurfaceOfAVolume";
192 theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
193}
194/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
195//
196void G4AdjointPrimaryGeneratorAction::ConsiderParticleAsPrimary(const G4String& particle_name)
197{
198 if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
199 PrimariesConsideredInAdjointSim[particle_name]=true;
200 }
201 UpdateListOfPrimaryParticles();
202}
203/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
204//
205void G4AdjointPrimaryGeneratorAction::NeglectParticleAsPrimary(const G4String& particle_name)
206{
207 if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
208 PrimariesConsideredInAdjointSim[particle_name]= false;
209 }
210 UpdateListOfPrimaryParticles();
211}
212/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
213//
214void G4AdjointPrimaryGeneratorAction::UpdateListOfPrimaryParticles()
215{
216 G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
217 ListOfPrimaryFwdParticles.clear();
218 ListOfPrimaryAdjParticles.clear();
219 std::map<G4String, G4bool>::iterator iter;
220 for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
221 if(iter->second) {
222 G4String fwd_particle_name = iter->first;
223 if ( fwd_particle_name != "ion") {
224 G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
225 ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
226 ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
227 }
228 else {
229 if (fwd_ion ){
230 ion_name=fwd_ion->GetParticleName();
231 G4String adj_ion_name=G4String("adj_") +ion_name;
232 ListOfPrimaryFwdParticles.push_back(fwd_ion);
233 ListOfPrimaryAdjParticles.push_back(adj_ion);
234 }
235 else {
236 ListOfPrimaryFwdParticles.push_back(0);
237 ListOfPrimaryAdjParticles.push_back(0);
238
239 }
240 }
241 }
242 }
243}
244/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
245//
246void G4AdjointPrimaryGeneratorAction::SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon)
247{
248 fwd_ion = fwdIon;
249 adj_ion = adjointIon;
250 UpdateListOfPrimaryParticles();
251}
252
Note: See TracBrowser for help on using the repository browser.