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