Geant4_10
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 //
29 
30 #include "G4ParticleDefinition.hh"
31 #include "G4ParticleTable.hh"
32 #include "G4VProcess.hh"
33 
34 #include "Randomize.hh"
35 
37  G4String name)
38  : G4VBiasingOperator(name),
39  fFirstProcess(0),
40  fLastProcess(0),
41  fSetup(true)
42 {
43  fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
44 
45  if ( fParticleToBias == 0 )
46  {
48  ed << "Particle `" << particleName << "' not found !" << G4endl;
49  G4Exception("GB01BOptrChangeCrossSection(...)",
50  "exGB01.01",
52  ed);
53  }
54 }
55 
57 {
58  for ( std::map< const G4BiasingProcessInterface*, G4BOptnChangeCrossSection* >::iterator
59  it = fChangeCrossSectionOperations.begin() ;
60  it != fChangeCrossSectionOperations.end() ;
61  it++ ) delete (*it).second;
62 }
63 
65 GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(const G4Track* track,
67  callingProcess)
68 {
69 
70  // -- Check if current particle type is the one to bias:
71  if ( track->GetDefinition() != fParticleToBias ) return 0;
72 
73  // -----------------------------------------------------------------------
74  // -- Setup stage ( Might consider moving this in a less called method. ):
75  // -----------------------------------------------------------------------
76  // -- Start by collecting processes under biasing, create needed biasing operations
77  // -- and assocation operations to these processes:
78  if ( fSetup )
79  {
80  if ( ( fFirstProcess == 0 ) &&
81  ( callingProcess->GetIsFirstPostStepGPILInterface() ) ) fFirstProcess = callingProcess;
82  if ( fLastProcess == 0 )
83  {
84  G4String operationName = "XSchange-" +
85  callingProcess->GetWrappedProcess()->GetProcessName();
86  fChangeCrossSectionOperations[callingProcess] =
87  new G4BOptnChangeCrossSection(operationName);
88  if ( callingProcess->GetIsLastPostStepGPILInterface() )
89  {
90  fLastProcess = callingProcess;
91  fSetup = false;
92  }
93  }
94  }
95 
96 
97 
98  // ---------------------------------------------------------------------
99  // -- select and setup the biasing operation for current callingProcess:
100  // ---------------------------------------------------------------------
101  // -- Check if the analog cross-section well defined : for example, the conversion
102  // -- process for a gamma below e+e- creation threshold has an DBL_MAX interaction
103  // -- length. Nothing is done in this case (ie, let analog process to deal with the case)
104  G4double analogInteractionLength =
105  callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
106  if ( analogInteractionLength > DBL_MAX/10. ) return 0;
107 
108  // -- Analog cross-section is well-defined:
109  G4double analogXS = 1./analogInteractionLength;
110 
111  // -- Choose a constant cross-section bias. But at this level, this factor can be made
112  // -- direction dependent, like in the exponential transform MCNP case, or it
113  // -- can be chosen differently, depending on the process, etc.
114  G4double XStransformation = 2.0 ;
115 
116  // -- fetch the operation associated to this callingProcess:
117  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
118  // -- get the operation that was proposed to the process in the previous step:
119  G4VBiasingOperation* previousOperation = callingProcess->GetPreviousOccurenceBiasingOperation();
120 
121  // -- now setup the operation to be returned to the process: this
122  // -- consists in setting the biased cross-section, and in asking
123  // -- the operation to sample its exponential interaction law.
124  // -- To do this, to first order, the two lines:
125  // operation->SetBiasedCrossSection( XStransformation * analogXS );
126  // operation->Sample();
127  // -- are correct and sufficient.
128  // -- But, to avoid having to shoot a random number each time, we sample
129  // -- only on the first time the operation is proposed, or if the interaction
130  // -- occured. If the interaction did not occur for the process in the previous,
131  // -- we update the number of interaction length instead of resampling.
132  if ( previousOperation == 0 )
133  {
134  operation->SetBiasedCrossSection( XStransformation * analogXS );
135  operation->Sample();
136  }
137  else
138  {
139  if ( previousOperation != operation )
140  {
141  // -- should not happen !
143  ed << " Logic problem in operation handling !" << G4endl;
144  G4Exception("GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(...)",
145  "exGB01.02",
146  JustWarning,
147  ed);
148  return 0;
149  }
150  if ( operation->GetInteractionOccured() )
151  {
152  operation->SetBiasedCrossSection( XStransformation * analogXS );
153  operation->Sample();
154  }
155  else
156  {
157  operation->UpdateForStep( callingProcess->GetPreviousStepSize() );
158  operation->SetBiasedCrossSection( XStransformation * analogXS );
159  }
160  }
161 
162  // operation->SetBiasedCrossSection( XStransformation * analogXS );
163  // operation->Sample();
164 
165  return operation;
166 
167 }
168 
169 
170 void GB01BOptrChangeCrossSection::
171 OperationApplied(const G4BiasingProcessInterface* callingProcess,
173  G4VBiasingOperation* occurenceOperationApplied,
174  G4double,
176  const G4VParticleChange* )
177 {
178  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
179  if ( operation == occurenceOperationApplied ) operation->SetInteractionOccured();
180 }
181 
G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const XML_Char * name
Definition: expat.h:151
GB01BOptrChangeCrossSection(G4String particleToBias, G4String name="ChangeXS")
G4BiasingAppliedCase
G4VProcess * GetWrappedProcess() const
G4bool GetIsFirstPostStepGPILInterface(G4bool physOnly=true) const
G4bool GetIsLastPostStepGPILInterface(G4bool physOnly=true) 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()
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
void UpdateForStep(G4double stepLength)
G4VBiasingOperation * GetPreviousOccurenceBiasingOperation() const