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