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