ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointPrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointPrimaryGeneratorAction.cc
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 //
28 // Class Name: G4AdjointPrimaryGeneratorAction
29 // Author: L. Desorgher
30 // Organisation: SpaceIT GmbH
31 // Contract: ESA contract 21435/08/NL/AT
32 // Customer: ESA/ESTEC
34 
36 
37 #include "G4PhysicalConstants.hh"
38 #include "G4Event.hh"
39 #include "G4ParticleTable.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4AdjointSimManager.hh"
44 #include "G4Gamma.hh"
45 
47 //
49  : Emin(0.), Emax(0.), EminIon(0.), EmaxIon(0.),
50  index_particle(100000),
51  radius_spherical_source(0.), fwd_ion(0), adj_ion(0),
52  ion_name("not_defined")
53 {
55 
60 
66 
67 }
69 //
71 {
73 }
75 //
77 {
78 
79  G4int evt_id=anEvent->GetEventID();
80  size_t n=ListOfPrimaryAdjParticles.size();
81  index_particle=size_t(evt_id)-n*(size_t(evt_id)/n);
82 
83 
84  G4double E1=Emin;
85  G4double E2=Emax;
86  if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
87 
88  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
89  E1=EminIon;
90  E2=EmaxIon;
91  }
92  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
94  E1=EminIon*A;
95  E2=EmaxIon*A;
96  }
97  //Generate first the forwrad primaries
99  G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
100 
101  p=fwdPrimVertex->GetPrimary()->GetMomentum();
102  pos=fwdPrimVertex->GetPosition();
103  G4double pmag=p.mag();
105  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
106 
107  G4double weight_correction=1.;
108  //For gamma generate the particle along the backward ray
109  G4ThreeVector dir=-p/p.mag();
110 
111  /*if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma"){
112 
113  theAdjointPrimaryGenerator
114  ->ComputeAccumulatedDepthVectorAlongBackRay(pos,dir,ekin,ListOfPrimaryAdjParticles[index_particle]);
115 
116 
117  G4double distance = theAdjointPrimaryGenerator
118  ->SampleDistanceAlongBackRayAndComputeWeightCorrection(weight_correction);
119 
120  //pos=pos+dir*distance;
121  //fwdPrimVertex->SetPosition(pos[0],pos[1],pos[2]);
122  }
123  */
124  weight_correction=1.;
125 
128  fwdPrimVertex->SetWeight(weight);
129  for (int i=0;i<nb_fwd_gammas_per_event-1;i++){
130  G4PrimaryVertex* newFwdPrimVertex = new G4PrimaryVertex();
131  newFwdPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
132  newFwdPrimVertex->SetT0(0.);
133  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
134  p.x(),p.y(),p.z());
135  newFwdPrimVertex->SetPrimary(aPrimParticle);
136  newFwdPrimVertex->SetWeight(weight);
137  anEvent->AddPrimaryVertex(newFwdPrimVertex);
138  }
139  }
140 
141  //Now generate the adjoint primaries
142  G4PrimaryVertex* adjPrimVertex = new G4PrimaryVertex();
143  adjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
144  adjPrimVertex->SetT0(0.);
145  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
146  -p.x(),-p.y(),-p.z());
147 
148  adjPrimVertex->SetPrimary(aPrimParticle);
149  anEvent->AddPrimaryVertex(adjPrimVertex);
150 
151  //The factor pi is to normalise the weight to the directional flux
153  G4double adjoint_weight = weight_correction*ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
154  //if (ListOfPrimaryFwdParticles[index_particle] ==G4Gamma::Gamma()) adjoint_weight = adjoint_weight/3.;
155  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma") {
156  //The weight will be corrected at the end of the track if splitted tracks are used
157  adjoint_weight = adjoint_weight/nb_adj_primary_gammas_per_event;
158  for (int i=0;i<nb_adj_primary_gammas_per_event-1;i++){
159  G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
160  newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
161  newAdjPrimVertex->SetT0(0.);
162  aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
163  -p.x(),-p.y(),-p.z());
164  newAdjPrimVertex->SetPrimary(aPrimParticle);
165  newAdjPrimVertex->SetWeight(adjoint_weight);
166  anEvent->AddPrimaryVertex(newAdjPrimVertex);
167  }
168  }
169  else if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_electron") {
170  //The weight will be corrected at the end of the track if splitted tracks are used
171  adjoint_weight = adjoint_weight/nb_adj_primary_electrons_per_event;
172  for (int i=0;i<nb_adj_primary_electrons_per_event-1;i++){
173  G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
174  newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
175  newAdjPrimVertex->SetT0(0.);
176  aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
177  -p.x(),-p.y(),-p.z());
178  newAdjPrimVertex->SetPrimary(aPrimParticle);
179  newAdjPrimVertex->SetWeight(adjoint_weight);
180  anEvent->AddPrimaryVertex(newAdjPrimVertex);
181  }
182  }
183  adjPrimVertex->SetWeight(adjoint_weight);
184 
185  //Call some methods of G4AdjointSimManager
189 
190  /* if ( !last_generated_part_was_adjoint ) {
191 
192  index_particle++;
193  if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
194 
195 
196  G4double E1=Emin;
197  G4double E2=Emax;
198  if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
199 
200  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
201  E1=EminIon;
202  E2=EmaxIon;
203  }
204  if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
205  G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
206  E1=EminIon*A;
207  E2=EmaxIon*A;
208  }
209  theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
210  ListOfPrimaryAdjParticles[index_particle],
211  E1,E2);
212  G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
213 
214 
215  p=aPrimVertex->GetPrimary()->GetMomentum();
216  pos=aPrimVertex->GetPosition();
217  G4double pmag=p.mag();
218 
219  G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
220  G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
221 
222 
223  //The factor pi is to normalise the weight to the directional flux
224  G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
225  G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
226 
227  aPrimVertex->SetWeight(adjoint_weight);
228 
229  last_generated_part_was_adjoint =true;
230  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
231  G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
232  }
233  else {
234  //fwd particle equivalent to the last generated adjoint particle ios generated
235  G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
236  aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
237  aPrimVertex->SetT0(0.);
238  G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
239  -p.x(),-p.y(),-p.z());
240 
241  aPrimVertex->SetPrimary(aPrimParticle);
242  anEvent->AddPrimaryVertex(aPrimVertex);
243  last_generated_part_was_adjoint =false;
244  G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
245  */
246 
247 }
249 //
251 {
252  Emin=val;
253  EminIon=val;
254 }
256 //
258 {
259  Emax=val;
260  EmaxIon=val;
261 }
263 //
265 {
266  EminIon=val;
267 }
269 //
271 {
272  EmaxIon=val;
273 }
275 //
277 {
278  // We generate N numbers of primaries with a 1/E energy law distribution.
279  // We have therefore an energy distribution function
280  // f(E)=C/E (1)
281  // with C a constant that is such that
282  // N=Integral(f(E),E1,E2)=C.std::log(E2/E1) (2)
283  // Therefore from (2) we get
284  // C=N/ std::log(E2/E1) (3)
285  // and
286  // f(E)=N/ std::log(E2/E1)/E (4)
287  //For the adjoint simulation we need a energy distribution f'(E)=1..
288  //To get that we need therefore to apply a weight to the primary
289  // W=1/f(E)=E*std::log(E2/E1)/N
290  //
291  return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
292 
293 }
295 //
297 {
299  center_spherical_source = center_pos;
300  type_of_adjoint_source ="Spherical";
302 }
304 //
306 {
307  type_of_adjoint_source ="ExternalSurfaceOfAVolume";
309 }
311 //
313 {
314  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
315  PrimariesConsideredInAdjointSim[particle_name]=true;
316  }
318 }
320 //
322 {
323  if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
324  PrimariesConsideredInAdjointSim[particle_name]= false;
325  }
327 }
329 //
331 {
335  std::map<G4String, G4bool>::iterator iter;
336  for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
337  if(iter->second) {
338  G4String fwd_particle_name = iter->first;
339  if ( fwd_particle_name != "ion") {
340  G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
341  ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
342  ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
343  /*
344  if ( fwd_particle_name == "gamma") {
345  for (G4int i=0;i<2;i++){
346  ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
347  ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
348  }
349  }
350  */
351  }
352  else {
353  if (fwd_ion ){
355  G4String adj_ion_name=G4String("adj_") +ion_name;
358  }
359  else {
360  ListOfPrimaryFwdParticles.push_back(0);
361  ListOfPrimaryAdjParticles.push_back(0);
362 
363  }
364  }
365  }
366  }
367 }
369 //
371 {
372  fwd_ion = fwdIon;
373  adj_ion = adjointIon;
375 }
376