Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GB01BOptrChangeCrossSection.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 //
28 //
32 
33 #include "G4ParticleDefinition.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4VProcess.hh"
36 
37 #include "Randomize.hh"
38 
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
44  G4String name)
45  : G4VBiasingOperator(name),
46  fSetup(true)
47 {
48  fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
49 
50  if ( fParticleToBias == 0 )
51  {
53  ed << "Particle `" << particleName << "' not found !" << G4endl;
54  G4Exception("GB01BOptrChangeCrossSection(...)",
55  "exGB01.01",
57  ed);
58  }
59 }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
64 {
65  for ( std::map< const G4BiasingProcessInterface*, G4BOptnChangeCrossSection* >::iterator
66  it = fChangeCrossSectionOperations.begin() ;
67  it != fChangeCrossSectionOperations.end() ;
68  it++ ) delete (*it).second;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
74 {
75  // --------------
76  // -- Setup stage:
77  // ---------------
78  // -- Start by collecting processes under biasing, create needed biasing
79  // -- operations and associate these operations to the processes:
80  if ( fSetup )
81  {
82  const G4ProcessManager* processManager = fParticleToBias->GetProcessManager();
83  const G4BiasingProcessSharedData* sharedData =
85  if ( sharedData ) // -- sharedData tested, as is can happen a user attaches an operator to a
86  // -- volume but without defined BiasingProcessInterface processes.
87  {
88  for ( size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++ )
89  {
90  const G4BiasingProcessInterface* wrapperProcess =
91  (sharedData->GetPhysicsBiasingProcessInterfaces())[i];
92  G4String operationName = "XSchange-" +
93  wrapperProcess->GetWrappedProcess()->GetProcessName();
94  fChangeCrossSectionOperations[wrapperProcess] =
95  new G4BOptnChangeCrossSection(operationName);
96  }
97  }
98  fSetup = false;
99  }
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
105 GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(const G4Track* track,
107  callingProcess)
108 {
109 
110  // -----------------------------------------------------
111  // -- Check if current particle type is the one to bias:
112  // -----------------------------------------------------
113  if ( track->GetDefinition() != fParticleToBias ) return 0;
114 
115 
116  // ---------------------------------------------------------------------
117  // -- select and setup the biasing operation for current callingProcess:
118  // ---------------------------------------------------------------------
119  // -- Check if the analog cross-section well defined : for example, the conversion
120  // -- process for a gamma below e+e- creation threshold has an DBL_MAX interaction
121  // -- length. Nothing is done in this case (ie, let analog process to deal with the case)
122  G4double analogInteractionLength =
123  callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
124  if ( analogInteractionLength > DBL_MAX/10. ) return 0;
125 
126  // -- Analog cross-section is well-defined:
127  G4double analogXS = 1./analogInteractionLength;
128 
129  // -- Choose a constant cross-section bias. But at this level, this factor can be made
130  // -- direction dependent, like in the exponential transform MCNP case, or it
131  // -- can be chosen differently, depending on the process, etc.
132  G4double XStransformation = 2.0 ;
133 
134  // -- fetch the operation associated to this callingProcess:
135  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
136  // -- get the operation that was proposed to the process in the previous step:
137  G4VBiasingOperation* previousOperation = callingProcess->GetPreviousOccurenceBiasingOperation();
138 
139  // -- now setup the operation to be returned to the process: this
140  // -- consists in setting the biased cross-section, and in asking
141  // -- the operation to sample its exponential interaction law.
142  // -- To do this, to first order, the two lines:
143  // operation->SetBiasedCrossSection( XStransformation * analogXS );
144  // operation->Sample();
145  // -- are correct and sufficient.
146  // -- But, to avoid having to shoot a random number each time, we sample
147  // -- only on the first time the operation is proposed, or if the interaction
148  // -- occured. If the interaction did not occur for the process in the previous,
149  // -- we update the number of interaction length instead of resampling.
150  if ( previousOperation == 0 )
151  {
152  operation->SetBiasedCrossSection( XStransformation * analogXS );
153  operation->Sample();
154  }
155  else
156  {
157  if ( previousOperation != operation )
158  {
159  // -- should not happen !
161  ed << " Logic problem in operation handling !" << G4endl;
162  G4Exception("GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(...)",
163  "exGB01.02",
164  JustWarning,
165  ed);
166  return 0;
167  }
168  if ( operation->GetInteractionOccured() )
169  {
170  operation->SetBiasedCrossSection( XStransformation * analogXS );
171  operation->Sample();
172  }
173  else
174  {
175  // -- update the 'interaction length' and underneath 'number of interaction lengths'
176  // -- for past step (this takes into accout the previous step cross-section value)
177  operation->UpdateForStep( callingProcess->GetPreviousStepSize() );
178  // -- update the cross-section value:
179  operation->SetBiasedCrossSection( XStransformation * analogXS );
180  // -- forces recomputation of the 'interaction length' taking into account above
181  // -- new cross-section value [tricky : to be improved]
182  operation->UpdateForStep( 0.0 );
183  }
184  }
185 
186  return operation;
187 
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
192 void GB01BOptrChangeCrossSection::
193 OperationApplied(const G4BiasingProcessInterface* callingProcess,
195  G4VBiasingOperation* occurenceOperationApplied,
196  G4double,
198  const G4VParticleChange* )
199 {
200  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
201  if ( operation == occurenceOperationApplied ) operation->SetInteractionOccured();
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ParticleDefinition * GetDefinition() const
const XML_Char * name
Definition: expat.h:151
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const std::vector< const G4BiasingProcessInterface * > & GetPhysicsBiasingProcessInterfaces() const
GB01BOptrChangeCrossSection(G4String particleToBias, G4String name="ChangeXS")
G4BiasingAppliedCase
G4VProcess * GetWrappedProcess() const
G4double GetCurrentInteractionLength() const
Definition: G4VProcess.hh:462
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ParticleTable * GetParticleTable()
G4ProcessManager * GetProcessManager() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Definition of the GB01BOptrChangeCrossSection class.
#define DBL_MAX
Definition: templates.hh:83
void UpdateForStep(G4double stepLength)
G4VBiasingOperation * GetPreviousOccurenceBiasingOperation() const
const G4BiasingProcessSharedData * GetSharedData() const