source: trunk/source/run/src/G4AdjointSimManager.cc @ 1331

Last change on this file since 1331 was 1197, checked in by garnier, 15 years ago

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