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