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 << 26 // 27 #include "G4AdjointBremsstrahlungModel.hh" 27 #include "G4AdjointBremsstrahlungModel.hh" 28 << 29 #include "G4AdjointCSManager.hh" 28 #include "G4AdjointCSManager.hh" >> 29 >> 30 #include "G4PhysicalConstants.hh" >> 31 #include "G4SystemOfUnits.hh" >> 32 >> 33 #include "G4Integrator.hh" >> 34 #include "G4TrackStatus.hh" >> 35 #include "G4ParticleChange.hh" 30 #include "G4AdjointElectron.hh" 36 #include "G4AdjointElectron.hh" 31 #include "G4AdjointGamma.hh" 37 #include "G4AdjointGamma.hh" 32 #include "G4Electron.hh" 38 #include "G4Electron.hh" 33 #include "G4EmModelManager.hh" << 39 #include "G4Timer.hh" 34 #include "G4Gamma.hh" << 35 #include "G4ParticleChange.hh" << 36 #include "G4PhysicalConstants.hh" << 37 #include "G4SeltzerBergerModel.hh" 40 #include "G4SeltzerBergerModel.hh" 38 #include "G4SystemOfUnits.hh" << 39 #include "G4TrackStatus.hh" << 40 41 41 ////////////////////////////////////////////// << 42 G4AdjointBremsstrahlungModel::G4AdjointBremsst << 43 : G4VEmAdjointModel("AdjointeBremModel") << 44 { << 45 fDirectModel = aModel; << 46 Initialize(); << 47 } << 48 42 49 ////////////////////////////////////////////// 43 //////////////////////////////////////////////////////////////////////////////// 50 G4AdjointBremsstrahlungModel::G4AdjointBremsst << 44 // 51 : G4VEmAdjointModel("AdjointeBremModel") << 45 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(G4VEmModel* aModel): 52 { << 46 G4VEmAdjointModel("AdjointeBremModel") 53 fDirectModel = new G4SeltzerBergerModel(); << 47 { 54 Initialize(); << 55 } << 56 << 57 ////////////////////////////////////////////// << 58 void G4AdjointBremsstrahlungModel::Initialize( << 59 { << 60 SetUseMatrix(false); 48 SetUseMatrix(false); 61 SetUseMatrixPerElement(false); 49 SetUseMatrixPerElement(false); >> 50 >> 51 theDirectStdBremModel = aModel; >> 52 theDirectEMModel=theDirectStdBremModel; >> 53 theEmModelManagerForFwdModels = new G4EmModelManager(); >> 54 isDirectModelInitialised = false; >> 55 G4VEmFluctuationModel* f=0; >> 56 G4Region* r=0; >> 57 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r); 62 58 63 fEmModelManagerForFwdModels = new G4EmModelM << 64 fEmModelManagerForFwdModels->AddEmModel(1, f << 65 SetApplyCutInRange(true); 59 SetApplyCutInRange(true); >> 60 highKinEnergy= 1.*GeV; >> 61 lowKinEnergy = 1.0*keV; 66 62 67 fElectron = G4Electron::Electron(); << 63 lastCZ =0.; 68 fGamma = G4Gamma::Gamma(); << 69 64 70 fAdjEquivDirectPrimPart = G4AdjointElectro << 65 71 fAdjEquivDirectSecondPart = G4AdjointGamma:: << 66 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); 72 fDirectPrimaryPart = fElectron; << 67 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); 73 fSecondPartSameType = false; << 68 theDirectPrimaryPartDef=G4Electron::Electron(); >> 69 second_part_of_same_type=false; 74 70 75 fCSManager = G4AdjointCSManager::GetAdjointC << 71 >> 72 CS_biasing_factor =1.; >> 73 >> 74 76 } 75 } >> 76 //////////////////////////////////////////////////////////////////////////////// >> 77 // >> 78 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(): >> 79 G4VEmAdjointModel("AdjointeBremModel") >> 80 { >> 81 SetUseMatrix(false); >> 82 SetUseMatrixPerElement(false); 77 83 >> 84 theDirectStdBremModel = new G4SeltzerBergerModel(); >> 85 theDirectEMModel=theDirectStdBremModel; >> 86 theEmModelManagerForFwdModels = new G4EmModelManager(); >> 87 isDirectModelInitialised = false; >> 88 G4VEmFluctuationModel* f=0; >> 89 G4Region* r=0; >> 90 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r); >> 91 // theDirectPenelopeBremModel =0; >> 92 SetApplyCutInRange(true); >> 93 highKinEnergy= 1.*GeV; >> 94 lowKinEnergy = 1.0*keV; >> 95 lastCZ =0.; >> 96 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); >> 97 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); >> 98 theDirectPrimaryPartDef=G4Electron::Electron(); >> 99 second_part_of_same_type=false; >> 100 } 78 ////////////////////////////////////////////// 101 //////////////////////////////////////////////////////////////////////////////// >> 102 // 79 G4AdjointBremsstrahlungModel::~G4AdjointBremss 103 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel() 80 { << 104 {if (theDirectStdBremModel) delete theDirectStdBremModel; 81 if(fEmModelManagerForFwdModels) << 105 if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels; 82 delete fEmModelManagerForFwdModels; << 83 } 106 } 84 107 85 ////////////////////////////////////////////// 108 //////////////////////////////////////////////////////////////////////////////// 86 void G4AdjointBremsstrahlungModel::SampleSecon << 109 // 87 const G4Track& aTrack, G4bool isScatProjToPr << 110 void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack, 88 G4ParticleChange* fParticleChange) << 111 G4bool IsScatProjToProjCase, >> 112 G4ParticleChange* fParticleChange) 89 { 113 { 90 if(!fUseMatrix) << 114 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 91 return RapidSampleSecondaries(aTrack, isSc << 92 115 93 const G4DynamicParticle* theAdjointPrimary = << 116 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 94 DefineCurrentMaterial(aTrack.GetMaterialCuts << 117 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 95 << 118 96 G4double adjointPrimKinEnergy = theAdjoint << 119 97 G4double adjointPrimTotalEnergy = theAdjoint << 120 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 98 << 121 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); 99 if(adjointPrimKinEnergy > GetHighEnergyLimit << 122 100 { << 123 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 101 return; << 124 return; 102 } << 125 } >> 126 >> 127 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, >> 128 IsScatProjToProjCase); >> 129 //Weight correction >> 130 //----------------------- >> 131 CorrectPostStepWeight(fParticleChange, >> 132 aTrack.GetWeight(), >> 133 adjointPrimKinEnergy, >> 134 projectileKinEnergy, >> 135 IsScatProjToProjCase); >> 136 >> 137 >> 138 //Kinematic >> 139 //--------- >> 140 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); >> 141 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; >> 142 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; >> 143 G4double projectileP = std::sqrt(projectileP2); >> 144 >> 145 >> 146 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel >> 147 //------------------------------------------------ >> 148 G4double u; >> 149 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; 103 150 104 G4double projectileKinEnergy = << 151 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; 105 SampleAdjSecEnergyFromCSMatrix(adjointPrim << 152 else u = - std::log(G4UniformRand()*G4UniformRand())/a2; 106 153 107 // Weight correction << 154 G4double theta = u*electron_mass_c2/projectileTotalEnergy; 108 CorrectPostStepWeight(fParticleChange, aTrac << 109 adjointPrimKinEnergy, << 110 isScatProjToProj); << 111 << 112 // Kinematic << 113 G4double projectileM0 = fAdjEquivDi << 114 G4double projectileTotalEnergy = projectileM << 115 G4double projectileP2 = << 116 projectileTotalEnergy * projectileTotalEne << 117 G4double projectileP = std::sqrt(projectileP << 118 155 119 // Angle of the gamma direction with the pro << 156 G4double sint = std::sin(theta); 120 // G4eBremsstrahlungModel << 157 G4double cost = std::cos(theta); 121 G4double u; << 158 122 if(0.25 > G4UniformRand()) << 159 G4double phi = twopi * G4UniformRand() ; 123 u = -std::log(G4UniformRand() * G4UniformR << 160 124 else << 161 G4ThreeVector projectileMomentum; 125 u = -std::log(G4UniformRand() * G4UniformR << 162 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame 126 << 163 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- 127 G4double theta = u * electron_mass_c2 / proj << 164 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); 128 G4double sint = std::sin(theta); << 165 G4ThreeVector dirProd=projectileMomentum-gammaMomentum; 129 G4double cost = std::cos(theta); << 166 G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); 130 << 167 G4double sint1 = std::sqrt(1.-cost1*cost1); 131 G4double phi = twopi * G4UniformRand(); << 168 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; 132 << 169 133 G4ThreeVector projectileMomentum = << 134 G4ThreeVector(std::cos(phi) * sint, std::s << 135 projectileP; // gamma frame << 136 if(isScatProjToProj) << 137 { // the adjoint primary is the scattered e << 138 G4ThreeVector gammaMomentum = << 139 (projectileTotalEnergy - adjointPrimTota << 140 G4ThreeVector(0., 0., 1.); << 141 G4ThreeVector dirProd = projectileMomentum << 142 G4double cost1 = std::cos(dirProd.a << 143 G4double sint1 = std::sqrt(1. - cos << 144 projectileMomentum = << 145 G4ThreeVector(std::cos(phi) * sint1, std << 146 projectileP; << 147 } 170 } 148 << 171 149 projectileMomentum.rotateUz(theAdjointPrimar 172 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); 150 << 173 151 if(!isScatProjToProj) << 174 152 { // kill the primary and add a secondary << 175 153 fParticleChange->ProposeTrackStatus(fStopA << 176 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 154 fParticleChange->AddSecondary( << 177 fParticleChange->ProposeTrackStatus(fStopAndKill); 155 new G4DynamicParticle(fAdjEquivDirectPri << 178 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 156 } << 179 } 157 else << 180 else { 158 { << 181 fParticleChange->ProposeEnergy(projectileKinEnergy); 159 fParticleChange->ProposeEnergy(projectileK << 182 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 160 fParticleChange->ProposeMomentumDirection( << 183 161 } << 184 } 162 } << 185 } 163 << 164 ////////////////////////////////////////////// 186 //////////////////////////////////////////////////////////////////////////////// 165 void G4AdjointBremsstrahlungModel::RapidSample << 187 // 166 const G4Track& aTrack, G4bool isScatProjToPr << 188 void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack, 167 G4ParticleChange* fParticleChange) << 189 G4bool IsScatProjToProjCase, 168 { << 190 G4ParticleChange* fParticleChange) 169 const G4DynamicParticle* theAdjointPrimary = << 191 { 170 DefineCurrentMaterial(aTrack.GetMaterialCuts << 192 >> 193 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); >> 194 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); >> 195 >> 196 >> 197 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); >> 198 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); >> 199 >> 200 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ >> 201 return; >> 202 } >> 203 >> 204 G4double projectileKinEnergy =0.; >> 205 G4double gammaEnergy=0.; >> 206 G4double diffCSUsed=0.; >> 207 if (!IsScatProjToProjCase){ >> 208 gammaEnergy=adjointPrimKinEnergy; >> 209 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); >> 210 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; >> 211 if (Emin>=Emax) return; >> 212 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand()); >> 213 diffCSUsed=CS_biasing_factor*lastCZ/projectileKinEnergy; >> 214 >> 215 } >> 216 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); >> 217 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); >> 218 if (Emin>=Emax) return; >> 219 G4double f1=(Emin-adjointPrimKinEnergy)/Emin; >> 220 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; >> 221 projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand())); >> 222 gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy; >> 223 diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy; >> 224 >> 225 } >> 226 >> 227 >> 228 >> 229 >> 230 //Weight correction >> 231 //----------------------- >> 232 //First w_corr is set to the ratio between adjoint total CS and fwd total CS >> 233 //if this has to be done in the model >> 234 //For the case of forced interaction this will be done in the PostStepDoIt of the >> 235 //forced interaction >> 236 //It is important to set the weight before the vreation of the secondary >> 237 // >> 238 G4double w_corr=additional_weight_correction_factor_for_post_step_outside_model; >> 239 if (correct_weight_for_post_step_in_model) { >> 240 w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); >> 241 } >> 242 //G4cout<<"Correction factor start in brem model "<<w_corr<<std::endl; >> 243 >> 244 >> 245 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model >> 246 //Here we consider the true diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term. >> 247 //Basically any other differential CS diffCS could be used here (example Penelope). >> 248 >> 249 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy); >> 250 /*G4cout<<"diffCS "<<diffCS <<std::endl; >> 251 G4cout<<"diffCS_Used "<<diffCSUsed <<std::endl;*/ >> 252 w_corr*=diffCS/diffCSUsed; >> 253 >> 254 >> 255 G4double new_weight = aTrack.GetWeight()*w_corr; >> 256 /*G4cout<<"New weight brem "<<new_weight<<std::endl; >> 257 G4cout<<"Weight correction brem "<<w_corr<<std::endl;*/ >> 258 fParticleChange->SetParentWeightByProcess(false); >> 259 fParticleChange->SetSecondaryWeightByProcess(false); >> 260 fParticleChange->ProposeParentWeight(new_weight); >> 261 >> 262 //Kinematic >> 263 //--------- >> 264 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); >> 265 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; >> 266 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; >> 267 G4double projectileP = std::sqrt(projectileP2); >> 268 >> 269 >> 270 //Use the angular model of the forward model to generate the gamma direction >> 271 //--------------------------------------------------------------------------- >> 272 //Dum dynamic particle to use the model >> 273 G4DynamicParticle * aDynPart = new G4DynamicParticle(G4Electron::Electron(),G4ThreeVector(0.,0.,1.)*projectileP); >> 274 >> 275 //Get the element from the direct model >> 276 const G4Element* elm = theDirectEMModel->SelectRandomAtom(currentCouple,G4Electron::Electron(), >> 277 projectileKinEnergy,currentTcutForDirectSecond); >> 278 G4int Z=elm->GetZasInt(); >> 279 G4double energy = aDynPart->GetTotalEnergy()-gammaEnergy; >> 280 G4ThreeVector projectileMomentum = >> 281 theDirectEMModel->GetAngularDistribution()->SampleDirection(aDynPart,energy,Z,currentMaterial)*projectileP; >> 282 G4double phi = projectileMomentum.getPhi(); 171 283 172 G4double adjointPrimKinEnergy = theAdjoint << 284 /* 173 G4double adjointPrimTotalEnergy = theAdjoint << 285 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel >> 286 //------------------------------------------------ >> 287 G4double u; >> 288 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; 174 289 175 if(adjointPrimKinEnergy > GetHighEnergyLimit << 290 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; 176 { << 291 else u = - std::log(G4UniformRand()*G4UniformRand())/a2; 177 return; << 178 } << 179 292 180 G4double projectileKinEnergy = 0.; << 293 G4double theta = u*electron_mass_c2/projectileTotalEnergy; 181 G4double gammaEnergy = 0.; << 182 G4double diffCSUsed = 0.; << 183 if(!isScatProjToProj) << 184 { << 185 gammaEnergy = adjointPrimKinEnergy; << 186 G4double Emax = GetSecondAdjEnergyMaxForPr << 187 G4double Emin = GetSecondAdjEnergyMinForPr << 188 if(Emin >= Emax) << 189 return; << 190 projectileKinEnergy = Emin * std::pow(Emax << 191 diffCSUsed = fCsBiasingFactor * f << 192 } << 193 else << 194 { << 195 G4double Emax = << 196 GetSecondAdjEnergyMaxForScatProjToProj(a << 197 G4double Emin = << 198 GetSecondAdjEnergyMinForScatProjToProj(a << 199 if(Emin >= Emax) << 200 return; << 201 G4double f1 = (Emin - adjointPrimKinEnergy << 202 G4double f2 = (Emax - adjointPrimKinEnergy << 203 projectileKinEnergy = << 204 adjointPrimKinEnergy / (1. - f1 * std::p << 205 gammaEnergy = projectileKinEnergy - adjoin << 206 diffCSUsed = << 207 fLastCZ * adjointPrimKinEnergy / project << 208 } << 209 294 210 // Weight correction: << 295 G4double sint = std::sin(theta); 211 // First w_corr is set to the ratio between << 296 G4double cost = std::cos(theta); 212 // if this has to be done in the model. << 297 213 // For the case of forced interaction this w << 298 G4double phi = twopi * G4UniformRand() ; 214 // the forced interaction. It is important << 299 G4ThreeVector projectileMomentum; 215 // creation of the secondary << 300 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame 216 G4double w_corr = fOutsideWeightFactor; << 301 */ 217 if(fInModelWeightCorr) << 302 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- 218 { << 303 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); 219 w_corr = fCSManager->GetPostStepWeightCorr << 304 G4ThreeVector dirProd=projectileMomentum-gammaMomentum; 220 } << 305 G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); 221 << 306 G4double sint1 = std::sqrt(1.-cost1*cost1); 222 // Then another correction is needed due to << 307 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; 223 // differential CS has been used rather than << 224 // direct model Here we consider the true di << 225 // numerical differentiation over Tcut of th << 226 // Migdal term. Basically any other differen << 227 // (example Penelope). << 228 G4double diffCS = DiffCrossSectionPerVolumeP << 229 fCurrentMaterial, projectileKinEnergy, gam << 230 w_corr *= diffCS / diffCSUsed; << 231 << 232 G4double new_weight = aTrack.GetWeight() * w << 233 fParticleChange->SetParentWeightByProcess(fa << 234 fParticleChange->SetSecondaryWeightByProcess << 235 fParticleChange->ProposeParentWeight(new_wei << 236 << 237 // Kinematic << 238 G4double projectileM0 = fAdjEquivDi << 239 G4double projectileTotalEnergy = projectileM << 240 G4double projectileP2 = << 241 projectileTotalEnergy * projectileTotalEne << 242 G4double projectileP = std::sqrt(projectileP << 243 << 244 // Use the angular model of the forward mode << 245 // Dummy dynamic particle to use the model << 246 G4DynamicParticle* aDynPart = << 247 new G4DynamicParticle(fElectron, G4ThreeVe << 248 << 249 // Get the element from the direct model << 250 const G4Element* elm = fDirectModel->SelectR << 251 fCurrentCouple, fElectron, projectileKinEn << 252 G4int Z = elm->GetZasInt(); << 253 G4double energy = aDynPart->GetTotalEnergy() << 254 G4ThreeVector projectileMomentum = << 255 fDirectModel->GetAngularDistribution()->Sa << 256 fC << 257 G4double phi = projectileMomentum.getPhi(); << 258 << 259 if(isScatProjToProj) << 260 { // the adjoint primary is the scattered e << 261 G4ThreeVector gammaMomentum = << 262 (projectileTotalEnergy - adjointPrimTota << 263 G4ThreeVector(0., 0., 1.); << 264 G4ThreeVector dirProd = projectileMomentum << 265 G4double cost1 = std::cos(dirProd.a << 266 G4double sint1 = std::sqrt(1. - cos << 267 projectileMomentum = << 268 G4ThreeVector(std::cos(phi) * sint1, std << 269 projectileP; << 270 } 308 } 271 309 272 projectileMomentum.rotateUz(theAdjointPrimar 310 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); 273 << 311 274 if(!isScatProjToProj) << 312 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 275 { // kill the primary and add a secondary << 313 fParticleChange->ProposeTrackStatus(fStopAndKill); 276 fParticleChange->ProposeTrackStatus(fStopA << 314 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 277 fParticleChange->AddSecondary( << 315 } 278 new G4DynamicParticle(fAdjEquivDirectPri << 316 else { 279 } << 317 fParticleChange->ProposeEnergy(projectileKinEnergy); 280 else << 318 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 281 { << 319 } 282 fParticleChange->ProposeEnergy(projectileK << 320 } 283 fParticleChange->ProposeMomentumDirection( << 321 //////////////////////////////////////////////////////////////////////////////// 284 } << 322 // 285 } << 323 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial, >> 324 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction >> 325 G4double kinEnergyProd // kinetic energy of the secondary particle >> 326 ) >> 327 {if (!isDirectModelInitialised) { >> 328 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); >> 329 isDirectModelInitialised =true; >> 330 } >> 331 /* >> 332 return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial, >> 333 kinEnergyProj, >> 334 kinEnergyProd); >> 335 */ >> 336 return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial, >> 337 kinEnergyProj, >> 338 kinEnergyProd); >> 339 } 286 340 287 ////////////////////////////////////////////// 341 //////////////////////////////////////////////////////////////////////////////// 288 G4double G4AdjointBremsstrahlungModel::DiffCro << 342 // 289 const G4Material* aMaterial, << 343 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1( 290 G4double kinEnergyProj, // kin energy of pr << 344 const G4Material* aMaterial, 291 G4double kinEnergyProd // kinetic energy o << 345 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction 292 ) << 346 G4double kinEnergyProd // kinetic energy of the secondary particle >> 347 ) 293 { 348 { 294 if(!fIsDirectModelInitialised) << 349 G4double dCrossEprod=0.; 295 { << 350 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); 296 fEmModelManagerForFwdModels->Initialise(fE << 351 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); 297 fIsDirectModelInitialised = true; << 352 298 } << 353 299 return G4VEmAdjointModel::DiffCrossSectionPe << 354 //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution 300 aMaterial, kinEnergyProj, kinEnergyProd); << 355 //This is what is applied in the discrete standard model before the rejection test that make a correction >> 356 //The application of the same rejection function is not possible here. >> 357 //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the >> 358 // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and >> 359 // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. >> 360 // In the future we plan to use the brem secondary spectra from the G4Penelope implementation >> 361 >> 362 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ >> 363 G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV); >> 364 dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV); >> 365 } >> 366 return dCrossEprod; >> 367 301 } 368 } 302 369 303 ////////////////////////////////////////////// 370 //////////////////////////////////////////////////////////////////////////////// 304 G4double G4AdjointBremsstrahlungModel::Adjoint << 371 // 305 const G4MaterialCutsCouple* aCouple, G4doubl << 372 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2( 306 G4bool isScatProjToProj) << 373 const G4Material* material, >> 374 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction >> 375 G4double kinEnergyProd // kinetic energy of the secondary particle >> 376 ) 307 { 377 { 308 static constexpr G4double maxEnergy = 100. * << 378 //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor 309 // 2.78.. == std::exp(1.) << 379 //used in the direct model 310 if(!fIsDirectModelInitialised) << 380 311 { << 381 G4double dCrossEprod=0.; 312 fEmModelManagerForFwdModels->Initialise(fE << 382 313 fIsDirectModelInitialised = true; << 383 const G4ElementVector* theElementVector = material->GetElementVector(); >> 384 const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector(); >> 385 G4double dum=0.; >> 386 G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001; >> 387 G4double dE=E2-E1; >> 388 for (size_t i=0; i<material->GetNumberOfElements(); i++) { >> 389 G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1); >> 390 G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2); >> 391 dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE; >> 392 } >> 393 >> 394 return dCrossEprod; >> 395 >> 396 } >> 397 //////////////////////////////////////////////////////////////////////////////// >> 398 // >> 399 G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple, >> 400 G4double primEnergy, >> 401 G4bool IsScatProjToProjCase) >> 402 { if (!isDirectModelInitialised) { >> 403 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); >> 404 isDirectModelInitialised =true; 314 } 405 } 315 if(fUseMatrix) << 406 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); 316 return G4VEmAdjointModel::AdjointCrossSect << 317 << 318 DefineCurrentMaterial(aCouple); 407 DefineCurrentMaterial(aCouple); 319 G4double Cross = 0.; << 408 G4double Cross=0.; 320 // this gives the constant above << 409 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above 321 fLastCZ = fDirectModel->CrossSectionPerVolum << 410 322 aCouple->GetMaterial(), fDirectPrimaryPart << 411 if (!IsScatProjToProjCase ){ 323 << 412 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); 324 if(!isScatProjToProj) << 413 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); 325 { << 414 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= CS_biasing_factor*lastCZ*std::log(Emax_proj/Emin_proj); 326 G4double Emax_proj = GetSecondAdjEnergyMax << 415 } 327 G4double Emin_proj = GetSecondAdjEnergyMin << 416 else { 328 if(Emax_proj > Emin_proj && primEnergy > f << 417 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); 329 Cross = fCsBiasingFactor * fLastCZ * std << 418 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond); 330 } << 419 if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy)); 331 else << 420 332 { << 421 } 333 G4double Emax_proj = GetSecondAdjEnergyMax << 422 return Cross; 334 G4double Emin_proj = << 423 } 335 GetSecondAdjEnergyMinForScatProjToProj(p << 424 336 if(Emax_proj > Emin_proj) << 425 G4double G4AdjointBremsstrahlungModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple, 337 Cross = fLastCZ * std::log((Emax_proj - << 426 G4double primEnergy, 338 Emax_proj / ( << 427 G4bool IsScatProjToProjCase) 339 } << 428 { 340 return Cross; << 429 return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); >> 430 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above >> 431 return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); >> 432 341 } 433 } >> 434 >> 435 >> 436 342 437