Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointSimManager.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$
27 //
29 // Class Name: G4AdjointCrossSurfChecker
30 // Author: L. Desorgher
31 // Organisation: SpaceIT GmbH
32 // Contract: ESA contract 21435/08/NL/AT
33 // Customer: ESA/ESTEC
35 
36 #include "G4AdjointSimManager.hh"
37 #include "G4Run.hh"
38 #include "G4RunManager.hh"
39 
40 #include "G4UserEventAction.hh"
41 #include "G4VUserPrimaryGeneratorAction.hh"
42 #include "G4UserTrackingAction.hh"
43 #include "G4UserSteppingAction.hh"
44 #include "G4UserStackingAction.hh"
45 #include "G4UserRunAction.hh"
46 
47 #include "G4AdjointPrimaryGeneratorAction.hh"
48 #include "G4AdjointSteppingAction.hh"
49 #include "G4AdjointStackingAction.hh"
50 
51 #include "G4AdjointSimMessenger.hh"
52 
53 #include "G4AdjointCrossSurfChecker.hh"
54 
55 #include "G4ParticleTable.hh"
56 #include "G4PhysicsLogVector.hh"
57 
59 //
60 G4AdjointSimManager* G4AdjointSimManager::instance = 0;
61 
63 //
64 G4AdjointSimManager::G4AdjointSimManager()
65 {
66  //Create adjoint actions;
67  //----------------------
68 
69  theAdjointRunAction = 0;
70  theAdjointPrimaryGeneratorAction = new G4AdjointPrimaryGeneratorAction();
71  theAdjointSteppingAction = new G4AdjointSteppingAction();
72  theAdjointEventAction = 0;
73  theAdjointTrackingAction = 0;
74  theAdjointStackingAction = new G4AdjointStackingAction();
75 
76  //Create messenger
77  //----------------
78  theMessenger = new G4AdjointSimMessenger(this);
79 
80  user_action_already_defined=false;
81  use_user_StackingAction = false;
82 
83  fUserTrackingAction= 0;
84  fUserEventAction= 0;
85  fUserSteppingAction= 0;
86  fUserPrimaryGeneratorAction= 0;
87  fUserRunAction= 0;
88  fUserStackingAction= 0;
89 
90  adjoint_sim_mode = false;
91 
92  normalisation_mode=3;
93 
94  nb_nuc=1.;
95 
96  welcome_message =true;
97 
98  /*electron_last_weight_vector = new G4PhysicsLogVector(1.e-20,1.e20,400);
99  proton_last_weight_vector = new G4PhysicsLogVector(1.e-20,1.e20,400);
100  gamma_last_weight_vector = new G4PhysicsLogVector(1.e-20,1.e20,400);*/
101 }
103 //
104 G4AdjointSimManager::~G4AdjointSimManager()
105 {
106  if (theAdjointRunAction) delete theAdjointRunAction;
107  if (theAdjointPrimaryGeneratorAction) delete theAdjointPrimaryGeneratorAction;
108  if (theAdjointSteppingAction) delete theAdjointSteppingAction;
109  if (theAdjointEventAction) delete theAdjointEventAction;
110  if (theAdjointTrackingAction) delete theAdjointTrackingAction;
111  if (theAdjointStackingAction) delete theAdjointStackingAction;
112  if (theMessenger) delete theMessenger;
113 }
115 //
116 G4AdjointSimManager* G4AdjointSimManager::GetInstance()
117 {
118  if (instance == 0) instance = new G4AdjointSimManager;
119  return instance;
120 }
122 //
123 void G4AdjointSimManager::RunAdjointSimulation(G4int nb_evt)
124 {
125  if (welcome_message) {
126  G4cout<<"****************************************************************"<<std::endl;
127  G4cout<<"*** Geant4 Reverse/Adjoint Monte Carlo mode ***"<<std::endl;
128  G4cout<<"*** Author: L.Desorgher ***"<<std::endl;
129  G4cout<<"*** Company: SpaceIT GmbH, Bern, Switzerland ***"<<std::endl;
130  G4cout<<"*** Sponsored by: ESA/ESTEC contract contract 21435/08/NL/AT ***"<<std::endl;
131  G4cout<<"****************************************************************"<<std::endl;
132  welcome_message=false;
133  }
134 
135  //Replace the user defined actions by the adjoint actions
136  //---------------------------------------------------------
137  SetAdjointPrimaryRunAndStackingActions();
138  SetRestOfAdjointActions();
139 
140  //Update the list of primaries
141  //-----------------------------
142  theAdjointPrimaryGeneratorAction->UpdateListOfPrimaryParticles();
143 
144  adjoint_sim_mode=true;
145 
146  ID_of_last_particle_that_reach_the_ext_source=0;
147 
148  //Make the run
149  //------------
150 
151  nb_evt_of_last_run =nb_evt;
152  G4RunManager::GetRunManager()->BeamOn(theAdjointPrimaryGeneratorAction->GetNbOfAdjointPrimaryTypes()*2*nb_evt);
153 
154  //Restore the user defined actions
155  //--------------------------------
156  ResetRestOfUserActions();
157  ResetUserPrimaryRunAndStackingActions();
158  adjoint_sim_mode=false;
159 
160  /*
161  //Register the weight vector
162  //--------------------------
163  std::ofstream FileOutputElectronWeight("ElectronWeight.txt", std::ios::out);
164  FileOutputElectronWeight<<std::setiosflags(std::ios::scientific);
165  FileOutputElectronWeight<<std::setprecision(6);
166  G4bool aBool = electron_last_weight_vector->Store(FileOutputElectronWeight, true);
167  FileOutputElectronWeight.close();
168 
169  std::ofstream FileOutputProtonWeight("ProtonWeight.txt", std::ios::out);
170  FileOutputProtonWeight<<std::setiosflags(std::ios::scientific);
171  FileOutputProtonWeight<<std::setprecision(6);
172  aBool = proton_last_weight_vector->Store(FileOutputProtonWeight, true);
173  FileOutputProtonWeight.close();
174 
175  std::ofstream FileOutputGammaWeight("GammaWeight.txt", std::ios::out);
176  FileOutputGammaWeight<<std::setiosflags(std::ios::scientific);
177  FileOutputGammaWeight<<std::setprecision(6);
178  aBool = gamma_last_weight_vector->Store(FileOutputGammaWeight, true);
179  FileOutputGammaWeight.close();
180  */
181 }
183 //
184 void G4AdjointSimManager::SetRestOfAdjointActions()
185 {
186  G4RunManager* theRunManager = G4RunManager::GetRunManager();
187 
188  if (!user_action_already_defined) DefineUserActions();
189 
190  //Replace the user action by the adjoint actions
191  //-------------------------------------------------
192 
193  theRunManager->SetUserAction(theAdjointEventAction);
194  theRunManager->SetUserAction(theAdjointSteppingAction);
195  theRunManager->SetUserAction(theAdjointTrackingAction);
196 }
198 //
199 void G4AdjointSimManager::SetAdjointPrimaryRunAndStackingActions()
200 {
201  G4RunManager* theRunManager = G4RunManager::GetRunManager();
202 
203  if (!user_action_already_defined) DefineUserActions();
204 
205  //Replace the user action by the adjoint actions
206  //-------------------------------------------------
207 
208  theRunManager->SetUserAction(theAdjointRunAction);
209  theRunManager->SetUserAction(theAdjointPrimaryGeneratorAction);
210  theRunManager->SetUserAction(theAdjointStackingAction);
211  if (use_user_StackingAction) theAdjointStackingAction->SetUserFwdStackingAction(fUserStackingAction);
212  else theAdjointStackingAction->SetUserFwdStackingAction(0);
213 }
215 //
216 void G4AdjointSimManager::ResetRestOfUserActions()
217 {
218  G4RunManager* theRunManager = G4RunManager::GetRunManager();
219 
220  //Restore the user defined actions
221  //-------------------------------
222 
223  theRunManager->SetUserAction(fUserEventAction);
224  theRunManager->SetUserAction(fUserSteppingAction);
225  theRunManager->SetUserAction(fUserTrackingAction);
226 }
227 
229 //
230 void G4AdjointSimManager::ResetUserPrimaryRunAndStackingActions()
231 {
232  G4RunManager* theRunManager = G4RunManager::GetRunManager();
233  //Restore the user defined actions
234  //-------------------------------
235  theRunManager->SetUserAction(fUserRunAction);
236  theRunManager->SetUserAction(fUserPrimaryGeneratorAction);
237  theRunManager->SetUserAction(fUserStackingAction);
238 }
240 //
241 void G4AdjointSimManager::DefineUserActions()
242 {
243  G4RunManager* theRunManager = G4RunManager::GetRunManager();
244  fUserTrackingAction= const_cast<G4UserTrackingAction* >( theRunManager->GetUserTrackingAction() );
245  fUserEventAction= const_cast<G4UserEventAction* >( theRunManager->GetUserEventAction() );
246  fUserSteppingAction= const_cast<G4UserSteppingAction* >( theRunManager->GetUserSteppingAction() );
247  fUserPrimaryGeneratorAction= const_cast<G4VUserPrimaryGeneratorAction* >( theRunManager->GetUserPrimaryGeneratorAction() );
248  fUserRunAction= const_cast<G4UserRunAction*>( theRunManager->GetUserRunAction() );
249  fUserStackingAction= const_cast<G4UserStackingAction* >( theRunManager->GetUserStackingAction() );
250  user_action_already_defined=true;
251 }
253 //
254 void G4AdjointSimManager::SetAdjointTrackingMode(G4bool aBool)
255 {
256  adjoint_tracking_mode = aBool;
257 
258  if (adjoint_tracking_mode) {
259  SetRestOfAdjointActions();
260  theAdjointStackingAction->SetAdjointMode(true);
261  theAdjointStackingAction->SetKillTracks(false);
262 
263  }
264  else {
265 
266  ResetRestOfUserActions();
267  theAdjointStackingAction->SetAdjointMode(false);
268  if (GetDidAdjParticleReachTheExtSource()){
269  theAdjointStackingAction->SetKillTracks(false);
270  RegisterAtEndOfAdjointTrack();
271  }
272  else theAdjointStackingAction->SetKillTracks(true);
273  }
274 }
276 //
277 G4bool G4AdjointSimManager::GetDidAdjParticleReachTheExtSource()
278 {
279  return theAdjointSteppingAction->GetDidAdjParticleReachTheExtSource();
280 }
282 //
283 std::vector<G4ParticleDefinition*> G4AdjointSimManager::GetListOfPrimaryFwdParticles()
284 {
285  return theAdjointPrimaryGeneratorAction->GetListOfPrimaryFwdParticles();
286 }
288 //
289 void G4AdjointSimManager::RegisterAtEndOfAdjointTrack()
290 {
291  last_pos = theAdjointSteppingAction->GetLastPosition();
292  last_direction = theAdjointSteppingAction->GetLastMomentum();
293  last_direction /=last_direction.mag();
294  last_cos_th = last_direction.z();
295  G4ParticleDefinition* aPartDef= theAdjointSteppingAction->GetLastPartDef();
296 
297  last_fwd_part_name= aPartDef->GetParticleName();
298 
299  last_fwd_part_name.remove(0,4);
300 
301  last_fwd_part_PDGEncoding=G4ParticleTable::GetParticleTable()->FindParticle(last_fwd_part_name)->GetPDGEncoding();
302 
303  std::vector<G4ParticleDefinition*> aList = theAdjointPrimaryGeneratorAction->GetListOfPrimaryFwdParticles();
304  last_fwd_part_index=-1;
305  size_t i=0;
306  while(i<aList.size() && last_fwd_part_index<0) {
307  if (aList[i]->GetParticleName() == last_fwd_part_name) last_fwd_part_index=i;
308  i++;
309  }
310 
311  last_ekin = theAdjointSteppingAction->GetLastEkin();
312  last_ekin_nuc = last_ekin;
313  if (aPartDef->GetParticleType() == "adjoint_nucleus") {
314  nb_nuc=double(aPartDef->GetBaryonNumber());
315  last_ekin_nuc /=nb_nuc;
316  }
317 
318  last_weight = theAdjointSteppingAction->GetLastWeight();
319 
320  /* G4PhysicsLogVector* theWeightVector=0;
321  if (last_fwd_part_name =="e-") theWeightVector=electron_last_weight_vector;
322  else if (last_fwd_part_name =="gamma") theWeightVector=gamma_last_weight_vector;
323  else if (last_fwd_part_name =="proton") theWeightVector=proton_last_weight_vector;
324 
325  if (theWeightVector){
326 
327  size_t ind = size_t(std::log10(last_weight/theAdjointPrimaryWeight)*10. + 200);
328  G4double low_val =theWeightVector->GetLowEdgeEnergy(ind);
329  G4bool aBool = true;
330  G4double bin_weight = theWeightVector->GetValue(low_val, aBool)+1.;
331  theWeightVector->PutValue(ind, bin_weight);
332  }
333  */
334  /*if ((last_weight/theAdjointPrimaryWeight)>1.) last_weight*=1000. ;
335  else if ( (last_weight/theAdjointPrimaryWeight)>0.1) last_weight*=100. ;
336  else if ( (last_weight/theAdjointPrimaryWeight)>0.01) last_weight*=10. ;*/
337 
338 
339  //G4cout <<"Last Weight "<<last_weight<<'\t'<<theAdjointPrimaryWeight<<'\t'<<last_weight/theAdjointPrimaryWeight<<std::endl;
340  /*if (last_weight/theAdjointPrimaryWeight >10.) {
341  G4cout<<"Warning a weight increase by a factor : "<<last_weight/theAdjointPrimaryWeight<<std::endl;
342  }
343  */
344 
345  ID_of_last_particle_that_reach_the_ext_source++;
346 }
348 //
349 G4bool G4AdjointSimManager::DefineSphericalExtSource(G4double radius, G4ThreeVector pos)
350 {
351  G4double area;
352  return G4AdjointCrossSurfChecker::GetInstance()->AddaSphericalSurface("ExternalSource", radius, pos, area);
353 }
355 //
356 G4bool G4AdjointSimManager::DefineSphericalExtSourceWithCentreAtTheCentreOfAVolume(G4double radius, const G4String& volume_name)
357 {
358  G4double area;
359  G4ThreeVector center;
360  return G4AdjointCrossSurfChecker::GetInstance()->AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume( "ExternalSource", radius, volume_name,center, area);
361 }
363 //
364 G4bool G4AdjointSimManager::DefineExtSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name)
365 {
366  G4double area;
367  return G4AdjointCrossSurfChecker::GetInstance()->AddanExtSurfaceOfAvolume( "ExternalSource", volume_name,area);
368 }
370 //
371 void G4AdjointSimManager::SetExtSourceEmax(G4double Emax)
372 {
373  theAdjointSteppingAction->SetExtSourceEMax(Emax);
374 }
376 //
377 G4bool G4AdjointSimManager::DefineSphericalAdjointSource(G4double radius, G4ThreeVector pos)
378 {
379  G4double area;
380  G4bool aBool = G4AdjointCrossSurfChecker::GetInstance()->AddaSphericalSurface("AdjointSource", radius, pos, area);
381  theAdjointPrimaryGeneratorAction->SetSphericalAdjointPrimarySource(radius, pos);
382  area_of_the_adjoint_source=area;
383  return aBool;
384 }
386 //
387 G4bool G4AdjointSimManager::DefineSphericalAdjointSourceWithCentreAtTheCentreOfAVolume(G4double radius, const G4String& volume_name)
388 {
389  G4double area;
390  G4ThreeVector center;
391  G4bool aBool = G4AdjointCrossSurfChecker::GetInstance()->AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume( "AdjointSource", radius, volume_name,center, area);
392  theAdjointPrimaryGeneratorAction->SetSphericalAdjointPrimarySource(radius, center);
393  area_of_the_adjoint_source=area;
394  return aBool;
395 }
397 //
398 G4bool G4AdjointSimManager::DefineAdjointSourceOnTheExtSurfaceOfAVolume(const G4String& volume_name)
399 {
400  G4double area;
401  G4bool aBool = G4AdjointCrossSurfChecker::GetInstance()->AddanExtSurfaceOfAvolume( "AdjointSource", volume_name,area);
402  area_of_the_adjoint_source=area;
403  if (aBool) {
404  theAdjointPrimaryGeneratorAction->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
405  }
406  return aBool;
407 }
409 //
410 void G4AdjointSimManager::SetAdjointSourceEmin(G4double Emin)
411 {
412  theAdjointPrimaryGeneratorAction->SetEmin(Emin);
413 }
415 //
416 void G4AdjointSimManager::SetAdjointSourceEmax(G4double Emax)
417 {
418  theAdjointPrimaryGeneratorAction->SetEmax(Emax);
419 }
421 //
422 void G4AdjointSimManager::ConsiderParticleAsPrimary(const G4String& particle_name)
423 {
424  theAdjointPrimaryGeneratorAction->ConsiderParticleAsPrimary(particle_name);
425 }
427 //
428 void G4AdjointSimManager::NeglectParticleAsPrimary(const G4String& particle_name)
429 {
430  theAdjointPrimaryGeneratorAction->NeglectParticleAsPrimary(particle_name);
431 }
433 //
434 /*void G4AdjointSimManager::SetPrimaryIon(G4int Z, G4int A)
435 {
436  theAdjointPrimaryGeneratorAction->SetPrimaryIon(Z, A);
437 }
438 */
440 //
441 void G4AdjointSimManager::SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon)
442 {
443  theAdjointPrimaryGeneratorAction->SetPrimaryIon(adjointIon, fwdIon);
444 }
446 //
447 const G4String& G4AdjointSimManager::GetPrimaryIonName()
448 {
449  return theAdjointPrimaryGeneratorAction->GetPrimaryIonName();
450 }
452 //
453 void G4AdjointSimManager::RegisterAdjointPrimaryWeight(G4double aWeight)
454 {
455  theAdjointPrimaryWeight = aWeight;
456  theAdjointSteppingAction->SetPrimWeight(aWeight);
457 }
458 
460 //
461 void G4AdjointSimManager::SetAdjointEventAction(G4UserEventAction* anAction)
462 {
463  theAdjointEventAction = anAction;
464 }
466 //
467 void G4AdjointSimManager::SetAdjointSteppingAction(G4UserSteppingAction* anAction)
468 {
469  theAdjointSteppingAction->SetUserAdjointSteppingAction(anAction);
470 }
472 //
473 void G4AdjointSimManager::SetAdjointStackingAction(G4UserStackingAction* anAction)
474 {
475  theAdjointStackingAction->SetUserAdjointStackingAction(anAction);
476 }
478 //
479 void G4AdjointSimManager::SetAdjointTrackingAction(G4UserTrackingAction* anAction)
480 {
481  theAdjointTrackingAction=anAction;
482 }
484 //
485 void G4AdjointSimManager::SetAdjointRunAction(G4UserRunAction* anAction)
486 {
487  theAdjointRunAction=anAction;
488 }
490 //
Generated on Sat May 25 2013 14:34:11 for Geant4 by   doxygen 1.8.4