| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file GB01/src/GB01BOptrChangeCrossSection <<
27 /// \brief Implementation of the GB01BOptrChan <<
28 // <<
29 #include "GB01BOptrChangeCrossSection.hh" 26 #include "GB01BOptrChangeCrossSection.hh"
30 <<
31 #include "G4BOptnChangeCrossSection.hh" <<
32 #include "G4BiasingProcessInterface.hh" 27 #include "G4BiasingProcessInterface.hh"
33 #include "G4InteractionLawPhysical.hh" << 28 #include "G4BOptnChangeCrossSection.hh"
>> 29
34 #include "G4ParticleDefinition.hh" 30 #include "G4ParticleDefinition.hh"
35 #include "G4ParticleTable.hh" 31 #include "G4ParticleTable.hh"
36 #include "G4VProcess.hh" 32 #include "G4VProcess.hh"
37 #include "Randomize.hh" <<
38 33
39 //....oooOO0OOooo........oooOO0OOooo........oo << 34 #include "Randomize.hh"
40 35
41 GB01BOptrChangeCrossSection::GB01BOptrChangeCr << 36 GB01BOptrChangeCrossSection::GB01BOptrChangeCrossSection(G4String particleName,
42 : G4VBiasingOperator(name), fSetup(true) << 37 G4String name)
>> 38 : G4VBiasingOperator(name),
>> 39 fFirstProcess(0),
>> 40 fLastProcess(0),
>> 41 fSetup(true)
43 { 42 {
44 fParticleToBias = G4ParticleTable::GetPartic 43 fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
45 << 44
46 if (fParticleToBias == 0) { << 45 if ( fParticleToBias == 0 )
47 G4ExceptionDescription ed; << 46 {
48 ed << "Particle `" << particleName << "' n << 47 G4ExceptionDescription ed;
49 G4Exception("GB01BOptrChangeCrossSection(. << 48 ed << "Particle `" << particleName << "' not found !" << G4endl;
50 } << 49 G4Exception("GB01BOptrChangeCrossSection(...)",
>> 50 "exGB01.01",
>> 51 JustWarning,
>> 52 ed);
>> 53 }
51 } 54 }
52 55
53 //....oooOO0OOooo........oooOO0OOooo........oo <<
54 <<
55 GB01BOptrChangeCrossSection::~GB01BOptrChangeC 56 GB01BOptrChangeCrossSection::~GB01BOptrChangeCrossSection()
56 { 57 {
57 for (std::map<const G4BiasingProcessInterfac << 58 for ( std::map< const G4BiasingProcessInterface*, G4BOptnChangeCrossSection* >::iterator
58 fChangeCrossSectionOperations.begin() << 59 it = fChangeCrossSectionOperations.begin() ;
59 it != fChangeCrossSectionOperations.end << 60 it != fChangeCrossSectionOperations.end() ;
60 delete (*it).second; << 61 it++ ) delete (*it).second;
61 } 62 }
62 63
63 //....oooOO0OOooo........oooOO0OOooo........oo << 64 G4VBiasingOperation*
64 << 65 GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(const G4Track* track,
65 void GB01BOptrChangeCrossSection::StartRun() << 66 const G4BiasingProcessInterface*
>> 67 callingProcess)
66 { 68 {
67 // -------------- << 69
68 // -- Setup stage: << 70 // -- Check if current particle type is the one to bias:
69 // --------------- << 71 if ( track->GetDefinition() != fParticleToBias ) return 0;
70 // -- Start by collecting processes under bi << 72
71 // -- operations and associate these operati << 73 // -----------------------------------------------------------------------
72 if (fSetup) { << 74 // -- Setup stage ( Might consider moving this in a less called method. ):
73 const G4ProcessManager* processManager = f << 75 // -----------------------------------------------------------------------
74 const G4BiasingProcessSharedData* sharedDa << 76 // -- Start by collecting processes under biasing, create needed biasing operations
75 G4BiasingProcessInterface::GetSharedData << 77 // -- and assocation operations to these processes:
76 if (sharedData) // -- sharedData tested, << 78 if ( fSetup )
77 // -- volume but without <<
78 { 79 {
79 for (size_t i = 0; i < (sharedData->GetP << 80 if ( ( fFirstProcess == 0 ) &&
80 const G4BiasingProcessInterface* wrapp << 81 ( callingProcess->GetIsFirstPostStepGPILInterface() ) ) fFirstProcess = callingProcess;
81 (sharedData->GetPhysicsBiasingProces << 82 if ( fLastProcess == 0 )
82 G4String operationName = << 83 {
83 "XSchange-" + wrapperProcess->GetWra << 84 G4String operationName = "XSchange-" +
84 fChangeCrossSectionOperations[wrapperP << 85 callingProcess->GetWrappedProcess()->GetProcessName();
85 new G4BOptnChangeCrossSection(operat << 86 fChangeCrossSectionOperations[callingProcess] =
86 } << 87 new G4BOptnChangeCrossSection(operationName);
>> 88 if ( callingProcess->GetIsLastPostStepGPILInterface() )
>> 89 {
>> 90 fLastProcess = callingProcess;
>> 91 fSetup = false;
>> 92 }
>> 93 }
87 } 94 }
88 fSetup = false; <<
89 } <<
90 } <<
91 95
92 //....oooOO0OOooo........oooOO0OOooo........oo <<
93 96
94 G4VBiasingOperation* GB01BOptrChangeCrossSecti <<
95 const G4Track* track, const G4BiasingProcess <<
96 { <<
97 // ----------------------------------------- <<
98 // -- Check if current particle type is the <<
99 // ----------------------------------------- <<
100 if (track->GetDefinition() != fParticleToBia <<
101 97
102 // ----------------------------------------- 98 // ---------------------------------------------------------------------
103 // -- select and setup the biasing operation 99 // -- select and setup the biasing operation for current callingProcess:
104 // ----------------------------------------- 100 // ---------------------------------------------------------------------
105 // -- Check if the analog cross-section well 101 // -- Check if the analog cross-section well defined : for example, the conversion
106 // -- process for a gamma below e+e- creatio 102 // -- process for a gamma below e+e- creation threshold has an DBL_MAX interaction
107 // -- length. Nothing is done in this case ( 103 // -- length. Nothing is done in this case (ie, let analog process to deal with the case)
108 G4double analogInteractionLength = << 104 G4double analogInteractionLength =
109 callingProcess->GetWrappedProcess()->GetCu 105 callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
110 if (analogInteractionLength > DBL_MAX / 10.) << 106 if ( analogInteractionLength > DBL_MAX/10. ) return 0;
111 107
112 // -- Analog cross-section is well-defined: 108 // -- Analog cross-section is well-defined:
113 G4double analogXS = 1. / analogInteractionLe << 109 G4double analogXS = 1./analogInteractionLength;
114 110
115 // -- Choose a constant cross-section bias. 111 // -- Choose a constant cross-section bias. But at this level, this factor can be made
116 // -- direction dependent, like in the expon 112 // -- direction dependent, like in the exponential transform MCNP case, or it
117 // -- can be chosen differently, depending o 113 // -- can be chosen differently, depending on the process, etc.
118 G4double XStransformation = 2.0; << 114 G4double XStransformation = 2.0 ;
119 << 115
120 // -- fetch the operation associated to this 116 // -- fetch the operation associated to this callingProcess:
121 G4BOptnChangeCrossSection* operation = fChan << 117 G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
122 // -- get the operation that was proposed to 118 // -- get the operation that was proposed to the process in the previous step:
123 G4VBiasingOperation* previousOperation = cal 119 G4VBiasingOperation* previousOperation = callingProcess->GetPreviousOccurenceBiasingOperation();
124 << 120
125 // -- now setup the operation to be returned 121 // -- now setup the operation to be returned to the process: this
126 // -- consists in setting the biased cross-s 122 // -- consists in setting the biased cross-section, and in asking
127 // -- the operation to sample its exponentia 123 // -- the operation to sample its exponential interaction law.
128 // -- To do this, to first order, the two li 124 // -- To do this, to first order, the two lines:
129 // operation->SetBiasedCrossSection( 125 // operation->SetBiasedCrossSection( XStransformation * analogXS );
130 // operation->Sample(); 126 // operation->Sample();
131 // -- are correct and sufficient. 127 // -- are correct and sufficient.
132 // -- But, to avoid having to shoot a random 128 // -- But, to avoid having to shoot a random number each time, we sample
133 // -- only on the first time the operation i 129 // -- only on the first time the operation is proposed, or if the interaction
134 // -- occured. If the interaction did not oc 130 // -- occured. If the interaction did not occur for the process in the previous,
135 // -- we update the number of interaction le 131 // -- we update the number of interaction length instead of resampling.
136 if (previousOperation == 0) { << 132 if ( previousOperation == 0 )
137 operation->SetBiasedCrossSection(XStransfo << 133 {
138 operation->Sample(); << 134 operation->SetBiasedCrossSection( XStransformation * analogXS );
139 } <<
140 else { <<
141 if (previousOperation != operation) { <<
142 // -- should not happen ! <<
143 G4ExceptionDescription ed; <<
144 ed << " Logic problem in operation handl <<
145 G4Exception("GB01BOptrChangeCrossSection <<
146 JustWarning, ed); <<
147 return 0; <<
148 } <<
149 if (operation->GetInteractionOccured()) { <<
150 operation->SetBiasedCrossSection(XStrans <<
151 operation->Sample(); 135 operation->Sample();
152 } 136 }
153 else { << 137 else
154 // -- update the 'interaction length' an << 138 {
155 // -- for past step (this takes into ac << 139 if ( previousOperation != operation )
156 operation->UpdateForStep(callingProcess- << 140 {
157 // -- update the cross-section value: << 141 // -- should not happen !
158 operation->SetBiasedCrossSection(XStrans << 142 G4ExceptionDescription ed;
159 // -- forces recomputation of the 'inter << 143 ed << " Logic problem in operation handling !" << G4endl;
160 // -- new cross-section value [tricky : << 144 G4Exception("GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(...)",
161 operation->UpdateForStep(0.0); << 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 }
162 } 160 }
163 } << 161
>> 162 // operation->SetBiasedCrossSection( XStransformation * analogXS );
>> 163 // operation->Sample();
164 164
165 return operation; 165 return operation;
>> 166
166 } 167 }
167 168
168 //....oooOO0OOooo........oooOO0OOooo........oo <<
169 169
170 void GB01BOptrChangeCrossSection::OperationApp << 170 void GB01BOptrChangeCrossSection::
171 << 171 OperationApplied(const G4BiasingProcessInterface* callingProcess,
172 << 172 G4BiasingAppliedCase,
173 << 173 G4VBiasingOperation* occurenceOperationApplied,
174 << 174 G4double,
>> 175 G4VBiasingOperation*,
>> 176 const G4VParticleChange* )
175 { 177 {
176 G4BOptnChangeCrossSection* operation = fChan 178 G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
177 if (operation == occurenceOperationApplied) << 179 if ( operation == occurenceOperationApplied ) operation->SetInteractionOccured();
178 } 180 }
179 181
180 //....oooOO0OOooo........oooOO0OOooo........oo <<
181 182