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