| 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 <<
27 #include "G4AdjointBremsstrahlungModel.hh" 26 #include "G4AdjointBremsstrahlungModel.hh"
28 <<
29 #include "G4AdjointCSManager.hh" 27 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointElectron.hh" << 28 #include "G4Integrator.hh"
31 #include "G4AdjointGamma.hh" <<
32 #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" 29 #include "G4TrackStatus.hh"
>> 30 #include "G4ParticleChange.hh"
>> 31 #include "G4AdjointElectron.hh"
>> 32 #include "G4Timer.hh"
40 33
41 ////////////////////////////////////////////// 34 ////////////////////////////////////////////////////////////////////////////////
42 G4AdjointBremsstrahlungModel::G4AdjointBremsst << 35 //
43 : G4VEmAdjointModel("AdjointeBremModel") << 36 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
44 { << 37 G4VEmAdjointModel("AdjointBremModel"),
45 fDirectModel = aModel; << 38 probsup(1.0),
46 Initialize(); << 39 MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
47 } << 40 LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
>> 41 theLPMflag(true)
48 42
>> 43 { isElectron= true;
>> 44 SetUseMatrix(true);
>> 45 SetUseMatrixPerElement(false);
>> 46 SetApplyCutInRange(true);
>> 47 SetIsIonisation(false);
>> 48 highKinEnergy= 100.*TeV;
>> 49 lowKinEnergy = 1.0*keV;
>> 50 theTimer =new G4Timer();
>> 51
>> 52 theTimer->Start();
>> 53 InitialiseParameters();
>> 54 theTimer->Stop();
>> 55 G4cout<<"Time elapsed in second for the initialidation of AdjointBrem "<<theTimer->GetRealElapsed()<<std::endl;
>> 56
>> 57 ModeldCS="MODEL1";
>> 58
>> 59 }
49 ////////////////////////////////////////////// 60 ////////////////////////////////////////////////////////////////////////////////
50 G4AdjointBremsstrahlungModel::G4AdjointBremsst << 61 //
51 : G4VEmAdjointModel("AdjointeBremModel") << 62 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
>> 63 {;}
>> 64 ////////////////////////////////////////////////////////////////////////////////
>> 65 //
>> 66 /*G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(
>> 67 const G4Material* aMaterial,
>> 68 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 69 G4double kinEnergyProd // kinetic energy of the secondary particle
>> 70 )
52 { 71 {
53 fDirectModel = new G4SeltzerBergerModel(); << 72
54 Initialize(); << 73 static const G4double
55 } << 74 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
56 << 75 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
>> 76 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
>> 77
>> 78 static const G4double
>> 79 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
>> 80 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
>> 81 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
>> 82
>> 83 static const G4double
>> 84 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
>> 85 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
>> 86 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
>> 87
>> 88 static const G4double
>> 89 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
>> 90 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
>> 91 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
>> 92
>> 93 static const G4double tlow = 1.*MeV;
>> 94
>> 95 G4double dCrossEprod=0.;
>> 96 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
>> 97 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
>> 98
>> 99
>> 100 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
>> 101
>> 102 G4double cross = 0.0;
>> 103
>> 104
>> 105
>> 106 G4double E1=kinEnergyProd;
>> 107 G4double E2=kinEnergyProd*1.000000001;
>> 108 G4double dE=(E2-E1);
>> 109
>> 110 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 111 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 112 G4double dum=0.;
>> 113
>> 114 for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
>> 115
>> 116 G4double fac=
>> 117
>> 118 cross += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
>> 119 kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
>> 120
>> 121
>> 122
>> 123 }
>> 124 dCrossEprod=(cross1-cross2)/dE; //first term
>> 125
>> 126 //Now come the correction
>> 127 //-----------------------
>> 128
>> 129 //First compute fsig for E1
>> 130 //-------------------------
>> 131
>> 132
>> 133 G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
>> 134 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
>> 135 *(aMaterial->GetElectronDensity());
>> 136
>> 137 G4double fsig = 0.;
>> 138 G4int nmax = 100;
>> 139 G4double vmin=std::log(E1);
>> 140 G4double vmax=std::log(kinEnergyProj) ;
>> 141 G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
>> 142 G4double u,fac,c,v,dv,y ;
>> 143 if(nn > 0) {
>> 144
>> 145 dv = (vmax-vmin)/nn ;
>> 146 v = vmin-dv ;
>> 147 for(G4int n=0; n<=nn; n++) {
>> 148
>> 149 v += dv;
>> 150 u = std::exp(v);
>> 151 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
>> 152 y = u/kinEnergyProj;
>> 153 fac *= (4.-4.*y+3.*y*y)/3.;
>> 154 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 155
>> 156 if ((n==0)||(n==nn)) c=0.5;
>> 157 else c=1. ;
>> 158
>> 159 fac *= c;
>> 160 fsig += fac;
>> 161 }
>> 162 y = E1/kinEnergyProj ;
>> 163 fsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 164
>> 165 }
>> 166 else {
>> 167 fsig = 1.;
>> 168 }
>> 169 if (fsig > 1.) fsig = 1.;
>> 170
>> 171 dCrossEprod*=fsig;
>> 172 //return dCrossEprod;
>> 173 //Now we compute dfsig
>> 174 //-------------------------
>> 175 G4double dfsig = 0.;
>> 176 nn=20;
>> 177 vmax=std::log(E2) ;
>> 178 dv = (vmax-vmin)/nn ;
>> 179 v = vmin-dv ;
>> 180 for(G4int n=0; n<=nn; n++) {
>> 181 v += dv;
>> 182 u = std::exp(v);
>> 183 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
>> 184 y = u/kinEnergyProj;
>> 185 fac *= (4.-4.*y+3.*y*y)/3.;
>> 186 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 187
>> 188 if ((n==0)||(n==nn)) c=0.5;
>> 189 else c=1. ;
>> 190
>> 191 fac *= c;
>> 192 dfsig += fac;
>> 193 }
>> 194 y = E1/kinEnergyProj;
>> 195 dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 196
>> 197 dCrossEprod+=dfsig*cross1/dE;
>> 198
>> 199
>> 200
>> 201
>> 202
>> 203 }
>> 204 return dCrossEprod;
>> 205
>> 206 }
>> 207 */
>> 208 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
>> 209 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 210 G4double kinEnergyProd // kinetic energy of the secondary particle
>> 211 )
>> 212 {if (ModeldCS=="MODEL2") return DiffCrossSectionPerVolumePrimToSecond2(aMaterial,
>> 213 kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 214 kinEnergyProd);
>> 215 if (ModeldCS=="MODEL3") return DiffCrossSectionPerVolumePrimToSecond3(aMaterial,
>> 216 kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 217 kinEnergyProd);
>> 218 return DiffCrossSectionPerVolumePrimToSecond1(aMaterial,
>> 219 kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 220 kinEnergyProd);
>> 221 }
57 ////////////////////////////////////////////// 222 ////////////////////////////////////////////////////////////////////////////////
58 void G4AdjointBremsstrahlungModel::Initialize( << 223 // the one used till now
>> 224 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond1(
>> 225 const G4Material* aMaterial,
>> 226 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 227 G4double kinEnergyProd // kinetic energy of the secondary particle
>> 228 )
59 { 229 {
60 SetUseMatrix(false); << 230 G4double dCrossEprod=0.;
61 SetUseMatrixPerElement(false); << 231 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
62 << 232 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
63 fEmModelManagerForFwdModels = new G4EmModelM << 233
64 fEmModelManagerForFwdModels->AddEmModel(1, f << 234
65 SetApplyCutInRange(true); << 235 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
>> 236
>> 237 G4double cross1 = 0.0;
>> 238 G4double cross2 = 0.0;
>> 239
>> 240
>> 241 G4double E1=kinEnergyProd;
>> 242 G4double E2=kinEnergyProd*1.01;
>> 243 G4double dE=(E2-E1);
>> 244
>> 245 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 246 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 247 G4double dum=0.;
>> 248
>> 249 for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
>> 250
>> 251 cross1 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
>> 252 kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
>> 253
>> 254 cross2 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
>> 255 kinEnergyProj, (*theElementVector)[i]->GetZ(), dum, E2);
>> 256
>> 257 }
>> 258 dCrossEprod=(cross1-cross2)/dE; //first term
>> 259
>> 260 //Now come the correction
>> 261 //-----------------------
>> 262
>> 263 //First compute fsig for E1
>> 264 //-------------------------
>> 265
>> 266
>> 267 G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
>> 268 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
>> 269 *(aMaterial->GetElectronDensity());
>> 270
>> 271 G4double fsig1 = 0.;
>> 272 G4int nmax = 100;
>> 273 G4double vmin=std::log(E1);
>> 274 G4double vmax=std::log(kinEnergyProj) ;
>> 275 G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
>> 276 G4double u,fac,c,v,dv,y ;
>> 277 if(nn > 0) {
>> 278
>> 279 dv = (vmax-vmin)/nn ;
>> 280 v = vmin-dv ;
>> 281 for(G4int n=0; n<=nn; n++) {
>> 282
>> 283 v += dv;
>> 284 u = std::exp(v);
>> 285 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
>> 286 y = u/kinEnergyProj;
>> 287 fac *= (4.-4.*y+3.*y*y)/3.;
>> 288 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 289
>> 290 if ((n==0)||(n==nn)) c=0.5;
>> 291 else c=1. ;
>> 292
>> 293 fac *= c;
>> 294 fsig1 += fac;
>> 295 }
>> 296 y = E1/kinEnergyProj ;
>> 297 fsig1 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 298
>> 299 }
>> 300 else {
>> 301 fsig1 = 1.;
>> 302 }
>> 303 if (fsig1 > 1.) fsig1 = 1.;
>> 304
>> 305 dCrossEprod*=fsig1;
>> 306
>> 307
>> 308 G4double fsig2 = 0.;
>> 309 vmin=std::log(E2);
>> 310 nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
>> 311 if(nn > 0) {
>> 312
>> 313 dv = (vmax-vmin)/nn ;
>> 314 v = vmin-dv ;
>> 315 for(G4int n=0; n<=nn; n++) {
>> 316
>> 317 v += dv;
>> 318 u = std::exp(v);
>> 319 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
>> 320 y = u/kinEnergyProj;
>> 321 fac *= (4.-4.*y+3.*y*y)/3.;
>> 322 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 323
>> 324 if ((n==0)||(n==nn)) c=0.5;
>> 325 else c=1. ;
>> 326
>> 327 fac *= c;
>> 328 fsig2 += fac;
>> 329 }
>> 330 y = E2/kinEnergyProj ;
>> 331 fsig2 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 332
>> 333 }
>> 334 else {
>> 335 fsig2 = 1.;
>> 336 }
>> 337 if (fsig2 > 1.) fsig2 = 1.;
>> 338
>> 339
>> 340 G4double dfsig=(fsig2-fsig1);
>> 341 dCrossEprod+=dfsig*cross1/dE;
>> 342
>> 343 dCrossEprod=(fsig1*cross1-fsig2*cross2)/dE;
>> 344
>> 345
>> 346
>> 347
>> 348
>> 349 /*if (fsig < 1.){
>> 350 //Now we compute dfsig
>> 351 //-------------------------
>> 352 G4double dfsig = 0.;
>> 353 nn=20;
>> 354 vmax=std::log(E2) ;
>> 355 dv = (vmax-vmin)/nn ;
>> 356 v = vmin-dv ;
>> 357 for(G4int n=0; n<=nn; n++) {
>> 358 v += dv;
>> 359 u = std::exp(v);
>> 360 fac = SupressionFunction(aMaterial, kinEnergyProj, u);
>> 361 y = u/kinEnergyProj;
>> 362 fac *= (4.-4.*y+3.*y*y)/3.;
>> 363 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 364
>> 365 if ((n==0)||(n==nn)) c=0.5;
>> 366 else c=1. ;
>> 367
>> 368 fac *= c;
>> 369 dfsig += fac;
>> 370 }
>> 371 y = E1/kinEnergyProj;
>> 372 dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 373 dCrossEprod+=dfsig*cross1/dE;
>> 374
>> 375 }
>> 376 */
>> 377
>> 378
>> 379
>> 380
>> 381
>> 382
>> 383
>> 384 }
>> 385 return dCrossEprod;
>> 386
>> 387 }
66 388
67 fElectron = G4Electron::Electron(); <<
68 fGamma = G4Gamma::Gamma(); <<
69 389
70 fAdjEquivDirectPrimPart = G4AdjointElectro << 390 ////////////////////////////////////////////////////////////////////////////////
71 fAdjEquivDirectSecondPart = G4AdjointGamma:: << 391 //
72 fDirectPrimaryPart = fElectron; << 392 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2(
73 fSecondPartSameType = false; << 393 const G4Material* aMaterial,
74 << 394 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
75 fCSManager = G4AdjointCSManager::GetAdjointC << 395 G4double kinEnergyProd // kinetic energy of the secondary particle
>> 396 )
>> 397 {
>> 398 G4double dCrossEprod=0.;
>> 399 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
>> 400 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
>> 401
>> 402
>> 403 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
>> 404
>> 405 G4double dEdX1 = 0.0;
>> 406 G4double dEdX2 = 0.0;
>> 407
>> 408
>> 409 G4double E1=kinEnergyProd;
>> 410 G4double E2=kinEnergyProd*1.001;
>> 411 G4double dE=(E2-E1);
>> 412 //G4double dum=0.;
>> 413
>> 414 dEdX1 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E1);
>> 415 dEdX2 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E2);
>> 416 dCrossEprod=(dEdX2-dEdX1)/dE/E1;
>> 417
>> 418
>> 419
>> 420
>> 421
>> 422
>> 423
>> 424 }
>> 425 return dCrossEprod;
>> 426
76 } 427 }
77 <<
78 ////////////////////////////////////////////// 428 ////////////////////////////////////////////////////////////////////////////////
79 G4AdjointBremsstrahlungModel::~G4AdjointBremss << 429 //
>> 430 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond3(
>> 431 const G4Material* aMaterial,
>> 432 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 433 G4double kinEnergyProd // kinetic energy of the secondary particle
>> 434 )
80 { 435 {
81 if(fEmModelManagerForFwdModels) << 436
82 delete fEmModelManagerForFwdModels; << 437 return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
>> 438 kinEnergyProj, // kinetic energy of the primary particle before the interaction
>> 439 kinEnergyProd);
>> 440
83 } 441 }
84 442
85 ////////////////////////////////////////////// 443 ////////////////////////////////////////////////////////////////////////////////
86 void G4AdjointBremsstrahlungModel::SampleSecon << 444 //
87 const G4Track& aTrack, G4bool isScatProjToPr << 445 G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material,
88 G4ParticleChange* fParticleChange) << 446 G4double kineticEnergy, G4double gammaEnergy)
89 { 447 {
90 if(!fUseMatrix) << 448 // supression due to the LPM effect+polarisation of the medium/
91 return RapidSampleSecondaries(aTrack, isSc << 449 // supression due to the polarisation alone
92 450
93 const G4DynamicParticle* theAdjointPrimary = <<
94 DefineCurrentMaterial(aTrack.GetMaterialCuts <<
95 451
96 G4double adjointPrimKinEnergy = theAdjoint << 452 G4double totEnergy = kineticEnergy+electron_mass_c2 ;
97 G4double adjointPrimTotalEnergy = theAdjoint << 453 G4double totEnergySquare = totEnergy*totEnergy ;
98 454
99 if(adjointPrimKinEnergy > GetHighEnergyLimit << 455 G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
100 { <<
101 return; <<
102 } <<
103 456
104 G4double projectileKinEnergy = << 457 G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
105 SampleAdjSecEnergyFromCSMatrix(adjointPrim <<
106 458
107 // Weight correction << 459 G4double electronDensity = material->GetElectronDensity();
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 460
119 // Angle of the gamma direction with the pro << 461 G4double sp = gammaEnergySquare/
120 // G4eBremsstrahlungModel << 462 (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
121 G4double u; <<
122 if(0.25 > G4UniformRand()) <<
123 u = -std::log(G4UniformRand() * G4UniformR <<
124 else <<
125 u = -std::log(G4UniformRand() * G4UniformR <<
126 <<
127 G4double theta = u * electron_mass_c2 / proj <<
128 G4double sint = std::sin(theta); <<
129 G4double cost = std::cos(theta); <<
130 <<
131 G4double phi = twopi * G4UniformRand(); <<
132 <<
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 463
149 projectileMomentum.rotateUz(theAdjointPrimar << 464 G4double supr = 1.0;
150 465
151 if(!isScatProjToProj) << 466 if (theLPMflag) {
152 { // kill the primary and add a secondary <<
153 fParticleChange->ProposeTrackStatus(fStopA <<
154 fParticleChange->AddSecondary( <<
155 new G4DynamicParticle(fAdjEquivDirectPri <<
156 } <<
157 else <<
158 { <<
159 fParticleChange->ProposeEnergy(projectileK <<
160 fParticleChange->ProposeMomentumDirection( <<
161 } <<
162 } <<
163 467
164 ////////////////////////////////////////////// << 468 G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
165 void G4AdjointBremsstrahlungModel::RapidSample <<
166 const G4Track& aTrack, G4bool isScatProjToPr <<
167 G4ParticleChange* fParticleChange) <<
168 { <<
169 const G4DynamicParticle* theAdjointPrimary = <<
170 DefineCurrentMaterial(aTrack.GetMaterialCuts <<
171 469
172 G4double adjointPrimKinEnergy = theAdjoint << 470 if (s2lpm < 1.) {
173 G4double adjointPrimTotalEnergy = theAdjoint <<
174 471
175 if(adjointPrimKinEnergy > GetHighEnergyLimit << 472 G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
176 { << 473 G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
177 return; << 474 G4double splim = LPMgEnergyLimit2/
178 } << 475 (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
>> 476 G4double w = 1.+1./splim ;
179 477
180 G4double projectileKinEnergy = 0.; << 478 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
181 G4double gammaEnergy = 0.; << 479 else w = s2lpm*(1.+1./sp);
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 480
210 // Weight correction: << 481 supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ;
211 // First w_corr is set to the ratio between << 482 supr /= sp;
212 // if this has to be done in the model. << 483 }
213 // For the case of forced interaction this w << 484
214 // the forced interaction. It is important << 485 }
215 // creation of the secondary << 486 return supr;
216 G4double w_corr = fOutsideWeightFactor; << 487 }
217 if(fInModelWeightCorr) <<
218 { <<
219 w_corr = fCSManager->GetPostStepWeightCorr <<
220 } <<
221 488
222 // Then another correction is needed due to << 489 ////////////////////////////////////////////////////////////////////////////////
223 // differential CS has been used rather than << 490 //
224 // direct model Here we consider the true di << 491 void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
225 // numerical differentiation over Tcut of th << 492 G4bool IsScatProjToProjCase,
226 // Migdal term. Basically any other differen << 493 G4ParticleChange* fParticleChange)
227 // (example Penelope). << 494 {
228 G4double diffCS = DiffCrossSectionPerVolumeP << 495
229 fCurrentMaterial, projectileKinEnergy, gam << 496 //G4cout<<"Adjoint Brem"<<std::endl;
230 w_corr *= diffCS / diffCSUsed; << 497 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
231 << 498
232 G4double new_weight = aTrack.GetWeight() * w << 499 size_t ind=0;
233 fParticleChange->SetParentWeightByProcess(fa << 500
234 fParticleChange->SetSecondaryWeightByProcess << 501 if (UseMatrixPerElement ) { //Select Material
235 fParticleChange->ProposeParentWeight(new_wei << 502 std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
236 << 503 if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
237 // Kinematic << 504 G4double rand_var= G4UniformRand();
238 G4double projectileM0 = fAdjEquivDi << 505 G4double SumCS=0.;
239 G4double projectileTotalEnergy = projectileM << 506 for (size_t i=0;i<CS_Vs_Element->size();i++){
240 G4double projectileP2 = << 507 SumCS+=(*CS_Vs_Element)[i];
241 projectileTotalEnergy * projectileTotalEne << 508 if (rand_var<=SumCS/lastCS){
242 G4double projectileP = std::sqrt(projectileP << 509 ind=i;
243 << 510 break;
244 // Use the angular model of the forward mode << 511 }
245 // Dummy dynamic particle to use the model << 512 }
246 G4DynamicParticle* aDynPart = << 513 }
247 new G4DynamicParticle(fElectron, G4ThreeVe << 514 else {
248 << 515 ind = currentMaterialIndex;
249 // Get the element from the direct model << 516 }
250 const G4Element* elm = fDirectModel->SelectR << 517
251 fCurrentCouple, fElectron, projectileKinEn << 518
252 G4int Z = elm->GetZasInt(); << 519 //Elastic inverse scattering modified compared to general G4VEmAdjointModel
253 G4double energy = aDynPart->GetTotalEnergy() << 520 //---------------------------
254 G4ThreeVector projectileMomentum = << 521 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
255 fDirectModel->GetAngularDistribution()->Sa << 522 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
256 fC << 523 //G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
257 G4double phi = projectileMomentum.getPhi(); << 524 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
258 << 525 return;
259 if(isScatProjToProj) << 526 }
260 { // the adjoint primary is the scattered e << 527
261 G4ThreeVector gammaMomentum = << 528 //Sample secondary energy
262 (projectileTotalEnergy - adjointPrimTota << 529 //-----------------------
263 G4ThreeVector(0., 0., 1.); << 530
264 G4ThreeVector dirProd = projectileMomentum << 531 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
265 G4double cost1 = std::cos(dirProd.a << 532 adjointPrimKinEnergy,
266 G4double sint1 = std::sqrt(1. - cos << 533 IsScatProjToProjCase);
267 projectileMomentum = << 534
268 G4ThreeVector(std::cos(phi) * sint1, std << 535
269 projectileP; << 536
270 } << 537
>> 538 //Weight correction
>> 539 //-----------------------
>> 540 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);
>> 541
>> 542
>> 543 //Kinematic
>> 544 //---------
>> 545
>> 546 G4double projectileM0 = electron_mass_c2;
>> 547 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
>> 548 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
>> 549 G4double projectileP = std::sqrt(projectileP2);
>> 550
>> 551
>> 552 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
>> 553 //------------------------------------------------
>> 554 G4double u;
>> 555 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
271 556
272 projectileMomentum.rotateUz(theAdjointPrimar << 557 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
>> 558 else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
273 559
274 if(!isScatProjToProj) << 560 G4double theta = u*electron_mass_c2/projectileTotalEnergy;
275 { // kill the primary and add a secondary <<
276 fParticleChange->ProposeTrackStatus(fStopA <<
277 fParticleChange->AddSecondary( <<
278 new G4DynamicParticle(fAdjEquivDirectPri <<
279 } <<
280 else <<
281 { <<
282 fParticleChange->ProposeEnergy(projectileK <<
283 fParticleChange->ProposeMomentumDirection( <<
284 } <<
285 } <<
286 561
287 ////////////////////////////////////////////// << 562 G4double sint = std::sin(theta);
288 G4double G4AdjointBremsstrahlungModel::DiffCro << 563 G4double cost = std::cos(theta);
289 const G4Material* aMaterial, << 564
290 G4double kinEnergyProj, // kin energy of pr << 565 G4double phi = twopi * G4UniformRand() ;
291 G4double kinEnergyProd // kinetic energy o << 566
292 ) << 567 G4ThreeVector projectileMomentum;
293 { << 568 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
294 if(!fIsDirectModelInitialised) << 569 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
295 { << 570 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
296 fEmModelManagerForFwdModels->Initialise(fE << 571 G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
297 fIsDirectModelInitialised = true; << 572 G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
>> 573 G4double sint1 = std::sqrt(1.-cost1*cost1);
>> 574 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
>> 575
298 } 576 }
299 return G4VEmAdjointModel::DiffCrossSectionPe << 577
300 aMaterial, kinEnergyProj, kinEnergyProd); << 578 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
301 } << 579
302 << 580
>> 581
>> 582 if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
>> 583 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 584 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
>> 585 //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
>> 586 }
>> 587 else {
>> 588 fParticleChange->ProposeEnergy(projectileKinEnergy);
>> 589 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
>> 590 //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
>> 591 }
>> 592 }
303 ////////////////////////////////////////////// 593 ////////////////////////////////////////////////////////////////////////////////
304 G4double G4AdjointBremsstrahlungModel::Adjoint << 594 //
305 const G4MaterialCutsCouple* aCouple, G4doubl << 595 void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel)
306 G4bool isScatProjToProj) << 596 {theDirectBremModel=aModel;
307 { << 597 DefineDirectEMModel(aModel);
308 static constexpr G4double maxEnergy = 100. * << 598 }
309 // 2.78.. == std::exp(1.) << 599 ////////////////////////////////////////////////////////////////////////////////
310 if(!fIsDirectModelInitialised) << 600 //
311 { << 601 void G4AdjointBremsstrahlungModel::InitialiseParameters()
312 fEmModelManagerForFwdModels->Initialise(fE << 602 {
313 fIsDirectModelInitialised = true; << 603 static const G4double
314 } << 604 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
315 if(fUseMatrix) << 605 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
316 return G4VEmAdjointModel::AdjointCrossSect << 606 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
317 << 607
318 DefineCurrentMaterial(aCouple); << 608 static const G4double
319 G4double Cross = 0.; << 609 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
320 // this gives the constant above << 610 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
321 fLastCZ = fDirectModel->CrossSectionPerVolum << 611 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
322 aCouple->GetMaterial(), fDirectPrimaryPart << 612
323 << 613 /* static const G4double
324 if(!isScatProjToProj) << 614 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
325 { << 615 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
326 G4double Emax_proj = GetSecondAdjEnergyMax << 616 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
327 G4double Emin_proj = GetSecondAdjEnergyMin << 617
328 if(Emax_proj > Emin_proj && primEnergy > f << 618 static const G4double
329 Cross = fCsBiasingFactor * fLastCZ * std << 619 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
330 } << 620 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
331 else << 621 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/
332 { << 622
333 G4double Emax_proj = GetSecondAdjEnergyMax << 623
334 G4double Emin_proj = << 624 const G4ElementTable* theElementTable = G4Element::GetElementTable();
335 GetSecondAdjEnergyMinForScatProjToProj(p << 625 FZ.clear();
336 if(Emax_proj > Emin_proj) << 626 ah1.clear();
337 Cross = fLastCZ * std::log((Emax_proj - << 627 ah2.clear();
338 Emax_proj / ( << 628 ah3.clear();
339 } << 629
340 return Cross; << 630 bh1.clear();
>> 631 bh2.clear();
>> 632 bh3.clear();
>> 633
>> 634 al0.clear();
>> 635 al1.clear();
>> 636 al2.clear();
>> 637
>> 638 bl0.clear();
>> 639 bl1.clear();
>> 640 bl2.clear();
>> 641 SigmaPerAtom.clear();
>> 642
>> 643 for (size_t j=0; j<theElementTable->size();j++){
>> 644
>> 645 G4Element* anElement=(*theElementTable)[j];
>> 646 G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
>> 647 FZ.push_back(lnZ* (4.- 0.55*lnZ));
>> 648 G4double ZZ = anElement->GetIonisation()->GetZZ3();
>> 649
>> 650 ah1.push_back(ah10 + ZZ* (ah11 + ZZ* ah12));
>> 651 ah2.push_back(ah20 + ZZ* (ah21 + ZZ* ah22));
>> 652 ah3.push_back(ah30 + ZZ* (ah31 + ZZ* ah32));
>> 653
>> 654 bh1.push_back(bh10 + ZZ* (bh11 + ZZ* bh12));
>> 655 bh2.push_back(bh20 + ZZ* (bh21 + ZZ* bh22));
>> 656 bh3.push_back(bh30 + ZZ* (bh31 + ZZ* bh32));
>> 657 /*SigmaPerAtom.push_back(theDirectEMModel->ComputeCrossSectionPerAtom(
>> 658 theDirectPrimaryPartDef,GetHighEnergyLimit()/2.,
>> 659 anElement->GetZ(),1.,GetLowEnergyLimit(),1.e20));*/
>> 660
>> 661
>> 662
>> 663 }
341 } 664 }
342 665