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

Last change on this file since 1168 was 966, checked in by garnier, 15 years ago

fichier ajoutes
File size: 21.4 KB
Line 
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.