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

Last change on this file since 1347 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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-04-beta-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.