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

Last change on this file since 1248 was 1197, checked in by garnier, 16 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.