source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc@ 1036

Last change on this file since 1036 was 966, checked in by garnier, 17 years ago

fichier ajoutes
File size: 21.4 KB
RevLine 
[966]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26#include "G4AdjointBremsstrahlungModel.hh"
27#include "G4AdjointCSManager.hh"
28#include "G4Integrator.hh"
29#include "G4TrackStatus.hh"
30#include "G4ParticleChange.hh"
31#include "G4AdjointElectron.hh"
32#include "G4Timer.hh"
33
34////////////////////////////////////////////////////////////////////////////////
35//
36G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
37 G4VEmAdjointModel("AdjointBremModel"),
38 probsup(1.0),
39 MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
40 LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
41 theLPMflag(true)
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}
60////////////////////////////////////////////////////////////////////////////////
61//
62G4AdjointBremsstrahlungModel::~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 )
71{
72
73 static const G4double
74 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
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*/
208G4double 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}
222////////////////////////////////////////////////////////////////////////////////
223// the one used till now
224G4double 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 )
229{
230 G4double dCrossEprod=0.;
231 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
232 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
233
234
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}
388
389
390////////////////////////////////////////////////////////////////////////////////
391//
392G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2(
393 const G4Material* aMaterial,
394 G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
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
427}
428////////////////////////////////////////////////////////////////////////////////
429//
430G4double 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 )
435{
436
437 return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
438 kinEnergyProj, // kinetic energy of the primary particle before the interaction
439 kinEnergyProd);
440
441}
442
443////////////////////////////////////////////////////////////////////////////////
444//
445G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material,
446 G4double kineticEnergy, G4double gammaEnergy)
447{
448 // supression due to the LPM effect+polarisation of the medium/
449 // supression due to the polarisation alone
450
451
452 G4double totEnergy = kineticEnergy+electron_mass_c2 ;
453 G4double totEnergySquare = totEnergy*totEnergy ;
454
455 G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
456
457 G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
458
459 G4double electronDensity = material->GetElectronDensity();
460
461 G4double sp = gammaEnergySquare/
462 (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
463
464 G4double supr = 1.0;
465
466 if (theLPMflag) {
467
468 G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
469
470 if (s2lpm < 1.) {
471
472 G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
473 G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
474 G4double splim = LPMgEnergyLimit2/
475 (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
476 G4double w = 1.+1./splim ;
477
478 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
479 else w = s2lpm*(1.+1./sp);
480
481 supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ;
482 supr /= sp;
483 }
484
485 }
486 return supr;
487}
488
489////////////////////////////////////////////////////////////////////////////////
490//
491void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
492 G4bool IsScatProjToProjCase,
493 G4ParticleChange* fParticleChange)
494{
495
496 //G4cout<<"Adjoint Brem"<<std::endl;
497 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
498
499 size_t ind=0;
500
501 if (UseMatrixPerElement ) { //Select Material
502 std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
503 if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
504 G4double rand_var= G4UniformRand();
505 G4double SumCS=0.;
506 for (size_t i=0;i<CS_Vs_Element->size();i++){
507 SumCS+=(*CS_Vs_Element)[i];
508 if (rand_var<=SumCS/lastCS){
509 ind=i;
510 break;
511 }
512 }
513 }
514 else {
515 ind = currentMaterialIndex;
516 }
517
518
519 //Elastic inverse scattering modified compared to general G4VEmAdjointModel
520 //---------------------------
521 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
522 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
523 //G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
524 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
525 return;
526 }
527
528 //Sample secondary energy
529 //-----------------------
530
531 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
532 adjointPrimKinEnergy,
533 IsScatProjToProjCase);
534
535
536
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. ;
556
557 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
558 else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
559
560 G4double theta = u*electron_mass_c2/projectileTotalEnergy;
561
562 G4double sint = std::sin(theta);
563 G4double cost = std::cos(theta);
564
565 G4double phi = twopi * G4UniformRand() ;
566
567 G4ThreeVector projectileMomentum;
568 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
569 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
570 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
571 G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
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
576 }
577
578 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
579
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}
593////////////////////////////////////////////////////////////////////////////////
594//
595void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel)
596{theDirectBremModel=aModel;
597 DefineDirectEMModel(aModel);
598}
599////////////////////////////////////////////////////////////////////////////////
600//
601void G4AdjointBremsstrahlungModel::InitialiseParameters()
602{
603 static const G4double
604 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
605 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
606 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
607
608 static const G4double
609 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
610 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
611 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
612
613 /* static const G4double
614 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
615 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
616 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
617
618 static const G4double
619 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
620 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
621 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/
622
623
624 const G4ElementTable* theElementTable = G4Element::GetElementTable();
625 FZ.clear();
626 ah1.clear();
627 ah2.clear();
628 ah3.clear();
629
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 }
664}
Note: See TracBrowser for help on using the repository browser.