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