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

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

fichier ajoutes
File size: 8.9 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 "G4AdjointPhotoElectricModel.hh"
27#include "G4AdjointCSManager.hh"
28
29
30#include "G4Integrator.hh"
31#include "G4TrackStatus.hh"
32#include "G4ParticleChange.hh"
33#include "G4AdjointElectron.hh"
34#include  "G4Gamma.hh"
35#include "G4AdjointGamma.hh"
36
37////////////////////////////////////////////////////////////////////////////////
38//
39G4AdjointPhotoElectricModel::G4AdjointPhotoElectricModel():
40 G4VEmAdjointModel("AdjointPEEffect")
41
42{ SetUseMatrix(false);
43  current_eEnergy =0.;
44  totAdjointCS=0.;
45}
46////////////////////////////////////////////////////////////////////////////////
47//
48G4AdjointPhotoElectricModel::~G4AdjointPhotoElectricModel()
49{;}
50
51////////////////////////////////////////////////////////////////////////////////
52//
53void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
54                                G4bool IsScatProjToProjCase,
55                                G4ParticleChange* fParticleChange)
56{ if (IsScatProjToProjCase) return ;
57
58  //Compute the totAdjointCS vectors if not already done for the current couple and electron energy
59  const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
60  const G4DynamicParticle* aDynPart =  aTrack.GetDynamicParticle() ;
61  G4double electronEnergy = aDynPart->GetKineticEnergy();
62  G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
63  totAdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
64                               
65 
66  //Sample gamma energy
67  //-------------
68        /////////////////////////////////////////////////////////////////////////////////
69//      Module:         G4ContinuousGainOfEnergy.hh
70//      Author:         L. Desorgher
71//      Date:           1 September 2007
72//      Organisation:   SpaceIT GmbH
73//      Customer:       ESA/ESTEC
74/////////////////////////////////////////////////////////////////////////////////
75//
76// CHANGE HISTORY
77// --------------
78//      ChangeHistory:
79//              1 September 2007 creation by L. Desorgher               
80//
81//-------------------------------------------------------------
82//      Documentation:
83//              Modell for the adjoint compton scattering
84//
85
86  //Sample element
87  //-------------
88   const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
89   const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
90   size_t nelm =  currentMaterial->GetNumberOfElements();
91   G4double rand_CS= totAdjointCS*G4UniformRand();
92   for (index_element=0; index_element<nelm-1; index_element++){
93        if (rand_CS<xsec[index_element]) break;
94   }
95       
96   //Sample shell and binding energy
97   //-------------
98   rand_CS= totAdjointCS*G4UniformRand()/theAtomNumDensityVector[index_element];
99   G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
100   G4int i  = 0; 
101   for (i=0; i<nShells-1; i++){
102        if (rand_CS<shell_prob[index_element][i]) break;
103   }
104   G4double gammaEnergy= electronEnergy+(*theElementVector)[index_element]->GetAtomicShell(i);
105       
106  //Sample cos theta
107  //Copy of the G4PEEffectModel cos theta sampling method ElecCosThetaDistribution.   
108  //This method cannot be used directly from G4PEEffectModel because it is a friend method. I should ask Vladimir to change that 
109  //------------------------------------------------------------------------------------------------   
110  //G4double cos_theta = theDirectPEEffectModel->ElecCosThetaDistribution(electronEnergy);
111       
112   G4double  cos_theta = 1.;
113   G4double gamma   = 1. + electronEnergy/electron_mass_c2;
114   if (gamma <= 5.) {
115        G4double beta  = std::sqrt(gamma*gamma-1.)/gamma;
116        G4double b     = 0.5*gamma*(gamma-1.)*(gamma-2);
117
118        G4double rndm,term,greject,grejsup;
119        if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*b);
120        else            grejsup = gamma*gamma*(1.+b+beta*b);
121
122        do {    rndm = 1.-2*G4UniformRand();
123                cos_theta = (rndm+beta)/(rndm*beta+1.);
124                term = 1.-beta*cos_theta;
125                greject = (1.-cos_theta*cos_theta)*(1.+b*term)/(term*term);
126        } while(greject < G4UniformRand()*grejsup);
127  }
128       
129  // direction of the adjoint gamma electron
130  //---------------------------------------
131 
132     
133  G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
134  G4double Phi     = twopi * G4UniformRand();
135  G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
136  G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz);
137  adjoint_gammaDirection.rotateUz(electronDirection);
138 
139 
140 
141  //Weight correction
142 //-----------------------                                         
143  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy);       
144 
145 
146 
147  //Create secondary and modify fParticleChange
148  //--------------------------------------------
149  G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
150                       G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
151  fParticleChange->ProposeTrackStatus(fStopAndKill);
152  fParticleChange->AddSecondary(anAdjointGamma);
153     
154       
155
156 
157
158}       
159
160////////////////////////////////////////////////////////////////////////////////
161//                     
162
163G4double G4AdjointPhotoElectricModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
164                                G4double electronEnergy,
165                                G4bool IsScatProjToProjCase)
166{ if (IsScatProjToProjCase) return 0.;
167  if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
168        totAdjointCS = 0.;
169        DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
170        const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
171        const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
172        size_t nelm =  currentMaterial->GetNumberOfElements();
173        for (index_element=0;index_element<nelm;index_element++){
174               
175                totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
176                xsec[index_element] = totAdjointCS;
177        } 
178  }     
179  return totAdjointCS;
180 
181 
182}
183////////////////////////////////////////////////////////////////////////////////
184//                     
185
186G4double G4AdjointPhotoElectricModel::AdjointCrossSectionPerAtom(const G4Element* anElement,G4double electronEnergy)
187{ 
188  G4int nShells = anElement->GetNbOfAtomicShells();
189  G4double Z= anElement->GetZ();
190  G4double N= anElement->GetN();
191  G4int i  = 0; 
192  G4double B0=anElement->GetAtomicShell(0);
193  G4double gammaEnergy = electronEnergy+B0;
194  G4double adjointCS = theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,N,0.,0.)*electronEnergy/gammaEnergy; 
195  shell_prob[index_element][0] = adjointCS;                                         
196  for (i=1;i<nShells;i++){
197        //G4cout<<i<<std::endl;
198        G4double Bi_= anElement->GetAtomicShell(i-1);
199        G4double Bi = anElement->GetAtomicShell(i);
200        //G4cout<<Bi_<<'\t'<<Bi<<std::endl;
201        if (electronEnergy <Bi_-Bi) {
202                gammaEnergy = electronEnergy+Bi;
203                adjointCS +=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,anElement->GetZ(),N,0.,0.)*electronEnergy/gammaEnergy;
204        }
205        shell_prob[index_element][i] = adjointCS;       
206 
207  }
208 
209  return adjointCS;
210 
211}                               
212////////////////////////////////////////////////////////////////////////////////
213//                     
214
215void G4AdjointPhotoElectricModel::DefineCurrentMaterialAndElectronEnergy(const G4MaterialCutsCouple* couple, G4double anEnergy)
216{ currentCouple   = const_cast<G4MaterialCutsCouple*> (couple);
217  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
218  currentCoupleIndex = couple->GetIndex();
219  currentMaterialIndex = currentMaterial->GetIndex();
220  current_eEnergy = anEnergy;   
221}
Note: See TracBrowser for help on using the repository browser.