Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4NeutronElectronElModel.cc
이 파일의 문서화 페이지로 가기
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 // $Id: G4NeutronElectronElModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
27 //
28 // Geant4 Header : G4NeutronElectronElModel
29 //
30 // 16.5.17: V.Grichine
31 //
32 
34 #include "G4SystemOfUnits.hh"
35 #include "G4ParticleTable.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4IonTable.hh"
38 #include "Randomize.hh"
39 #include "G4Integrator.hh"
40 #include "G4Electron.hh"
41 #include "G4PhysicsTable.hh"
42 #include "G4PhysicsLogVector.hh"
43 #include "G4PhysicsFreeVector.hh"
44 
45 
46 using namespace std;
47 using namespace CLHEP;
48 
50  : G4HadronElastic(name)
51 {
52  // neutron magneton squared
53 
54  fM = neutron_mass_c2; // neutron mass
55  fM2 = fM*fM;
57  fme2 = fme*fme;
58  fMv2 = 0.7056*GeV*GeV;
59 
60  SetMinEnergy( 0.001*GeV );
61  SetMaxEnergy( 10.*TeV );
63 
65  // PDG2016: sin^2 theta Weinberg
66 
67  fEnergyBin = 200;
68  fMinEnergy = 1.*MeV;
69  fMaxEnergy = 10000.*GeV;
71 
72  fAngleBin = 500;
73  fAngleTable = 0;
74 
75  fCutEnergy = 0.; // default value
76 
77  Initialise();
78 }
79 
81 
83 {
84  if( fEnergyVector )
85  {
86  delete fEnergyVector;
87  fEnergyVector = 0;
88  }
89  if( fAngleTable )
90  {
92  delete fAngleTable;
93  fAngleTable = nullptr;
94  }
95 }
96 
98 
99 void G4NeutronElectronElModel::ModelDescription(std::ostream& outFile) const
100 {
101 
102  outFile << "G4NeutronElectronElModel is a neutrino-electron (neutral current) elastic scattering\n"
103  << "model which uses the standard model \n"
104  << "transfer parameterization. The model is fully relativistic\n";
105 
106 }
107 
109 
111  G4Nucleus & targetNucleus)
112 {
113  G4bool result = false;
114  G4String pName = aTrack.GetDefinition()->GetParticleName();
115  // G4double minEnergy = 0.;
116  G4double energy = aTrack.GetTotalEnergy();
117 
118  if( fCutEnergy > 0. ) // min detected recoil electron energy
119  {
120  // minEnergy = 0.5*(fCutEnergy+sqrt(fCutEnergy*(fCutEnergy+2.*electron_mass_c2)));
121  }
122  if( pName == "neutron" &&
123  energy >= fMinEnergy && energy <= fMaxEnergy )
124  {
125  result = true;
126  }
127  G4int Z = targetNucleus.GetZ_asInt();
128  Z *= 1;
129 
130  return result;
131 }
132 
134 
136 {
137  G4double result = 0., sum, Tkin, dt, t1, t2;
138  G4int iTkin, jTransfer;
140 
142 
143  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
144  {
145  Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
146  fAm = CalculateAm(Tkin);
147  dt = 1./fAngleBin;
148 
150 
151  sum = 0.;
152 
153  for( jTransfer = 0; jTransfer < fAngleBin; jTransfer++)
154  {
155  t1 = dt*jTransfer;
156  t2 = t1 + dt;
157 
158  result = integral.Legendre96( this, &G4NeutronElectronElModel::XscIntegrand, t1, t2 );
159 
160  sum += result;
161  // G4cout<<sum<<", ";
162  vectorT->PutValue(jTransfer, t1, sum);
163  }
164  // G4cout<<G4endl;
165  fAngleTable->insertAt(iTkin,vectorT);
166  }
167  return;
168 }
169 
171 //
172 // sample recoil electron energy in lab frame
173 
175 {
176  G4double result = 0., position;
177  G4int iTkin, iTransfer;
178 
179  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
180  {
181  if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
182  }
183  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
184  if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
185 
186  position = (*(*fAngleTable)(iTkin))(fAngleBin-1)*G4UniformRand();
187 
188  // G4cout<<"position = "<<position<<G4endl;
189 
190  for( iTransfer = 0; iTransfer < fAngleBin; iTransfer++)
191  {
192  if( position <= (*(*fAngleTable)(iTkin))(iTransfer) ) break;
193  }
194  if (iTransfer >= fAngleBin-1) iTransfer = fAngleBin-1;
195 
196  // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
197 
198  result = GetTransfer(iTkin, iTransfer, position);
199 
200  // G4cout<<"t = "<<t<<G4endl;
201 
202 
203  return result;
204 }
205 
207 
208 G4double
210 {
211  G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
212 
213  if( iTransfer == 0 || iTransfer == fAngleBin-1 )
214  {
215  randTransfer = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer);
216  // iTransfer++;
217  }
218  else
219  {
220  if ( iTransfer >= G4int((*fAngleTable)(iTkin)->GetVectorLength()) )
221  {
222  iTransfer = (*fAngleTable)(iTkin)->GetVectorLength() - 1;
223  }
224  y1 = (*(*fAngleTable)(iTkin))(iTransfer-1);
225  y2 = (*(*fAngleTable)(iTkin))(iTransfer);
226 
227  x1 = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
228  x2 = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer);
229 
230  delta = y2 - y1;
231  mean = y2 + y1;
232 
233  if ( x1 == x2 ) randTransfer = x2;
234  else
235  {
236  // if ( y1 == y2 )
237 
238  if ( delta < epsilon*mean )
239  {
240  randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
241  }
242  else
243  {
244  randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
245  }
246  }
247  }
248  return randTransfer;
249 }
250 
252 //
253 // Rosenbluth relation (ultra-relativistic!) in the neutron rest frame,
254 // x = sin^2(theta/2), theta is the electron scattering angle
255 // Magnetic form factor in the dipole approximation.
256 
258 {
259  G4double result = 1., q2, znq2, znf, znf2, znf4;
260 
261  znq2 = 1. + 2.*fee*x/fM;
262 
263  q2 = 4.*fee2*x/znq2;
264 
265  znf = 1 + q2/fMv2;
266  znf2 = znf*znf;
267  znf4 = znf2*znf2;
268 
269  result /= ( x + fAm )*znq2*znq2*znf4;
270 
271  result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x;
272 
273  return result;
274 }
275 
277 //
278 //
279 
281  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
282 {
284 
285  const G4HadProjectile* aParticle = &aTrack;
286  G4double Tkin = aParticle->GetKineticEnergy();
287  fAm = CalculateAm( Tkin);
288  // G4double En = aParticle->GetTotalEnergy();
289 
290  if( Tkin <= LowestEnergyLimit() )
291  {
294  return &theParticleChange;
295  }
296  // sample e-scattering angle and make final state in lab frame
297 
298  G4double sin2ht = SampleSin2HalfTheta( Tkin); // in n-rrest frame
299 
300  // G4cout<<"sin2ht = "<<sin2ht<<G4endl;
301 
302  G4double eTkin = fee; // fM;
303 
304  eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2*sin2ht;
305 
306  eTkin -= fme;
307 
308  // G4cout<<"eTkin = "<<eTkin<<G4endl;
309 
310  if( eTkin > fCutEnergy )
311  {
312  G4double ePlab = sqrt( eTkin*(eTkin + 2.*fme) );
313 
314  // G4cout<<"ePlab = "<<ePlab<<G4endl;
315 
316  G4double cost = 1. - 2*sin2ht;
317 
318  if( cost > 1. ) cost = 1.;
319  if( cost < -1. ) cost = -1.;
320 
321  G4double sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
323 
324  G4ThreeVector eP( sint*std::cos(phi), sint*std::sin(phi), cost );
325  eP *= ePlab;
326  G4LorentzVector lvt2( eP, eTkin + electron_mass_c2 ); // recoil e- in n-rest frame
327 
328  G4LorentzVector lvp1 = aParticle->Get4Momentum();
329  G4LorentzVector lvt1(0.,0.,0.,electron_mass_c2);
330  G4LorentzVector lvsum = lvp1+lvt1;
331 
332  G4ThreeVector bst = lvp1.boostVector();
333  lvt2.boost(bst);
334 
335  // G4cout<<"lvt2 = "<<lvt2<<G4endl;
336 
337  G4DynamicParticle * aSec = new G4DynamicParticle( theElectron, lvt2 );
339 
340  G4LorentzVector lvp2 = lvsum-lvt2;
341 
342  // G4cout<<"lvp2 = "<<lvp2<<G4endl;
343 
344  G4double Tkin2 = lvp2.e()-aParticle->GetDefinition()->GetPDGMass();
347  }
348  else if( eTkin > 0.0 )
349  {
351  Tkin -= eTkin;
352 
353  if( Tkin > 0. )
354  {
357  }
358  }
359  else
360  {
363  }
364  G4int Z = targetNucleus.GetZ_asInt();
365  Z *= 1;
366 
367  return &theParticleChange;
368 }
369 
370 //
371 //
Float_t x
Definition: compare.C:6
const XML_Char * name
Definition: expat.h:151
void SetMinEnergy(G4double anEnergy)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void SetMaxEnergy(const G4double anEnergy)
TTree * t1
Definition: plottest35.C:26
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double GetLowEdgeEnergy(size_t binNumber) const
Float_t y1[n_points_granero]
Definition: compare.C:5
virtual G4bool IsApplicable(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4ParticleDefinition * theElectron
void SetMomentumChange(const G4ThreeVector &aV)
G4double LowestEnergyLimit() const
Float_t x1[n_points_granero]
Definition: compare.C:5
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
const G4String & GetParticleName() const
G4PhysicsLogVector * fEnergyVector
G4double GetTransfer(G4int iTkin, G4int iTransfer, G4double position)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void insertAt(size_t, G4PhysicsVector *)
G4double GetTotalEnergy() const
G4double GetPDGMass() const
Float_t y2[n_points_geant4]
Definition: compare.C:26
static constexpr double neutron_mass_c2
Float_t Z
void SetEnergyChange(G4double anEnergy)
static constexpr double TeV
Definition: G4SIunits.hh:218
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define position
Definition: xmlparse.cc:622
G4double CalculateAm(G4double momentum)
static constexpr double electron_mass_c2
#define G4UniformRand()
Definition: Randomize.hh:53
double energy
Definition: plottest35.C:25
static constexpr double eV
Definition: G4SIunits.hh:215
static G4Electron * Electron()
Definition: G4Electron.cc:94
virtual void ModelDescription(std::ostream &) const
const G4LorentzVector & Get4Momentum() const
Hep3Vector unit() const
G4double G4ParticleHPJENDLHEData::G4double result
TTree * t2
Definition: plottest35.C:36
G4double GetKineticEnergy() const
double epsilon(double density, double temperature)
static constexpr double twopi
Definition: SystemOfUnits.h:55
const G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
G4NeutronElectronElModel(const G4String &name="nu-e-elastic")
G4double SampleSin2HalfTheta(G4double Tkin)
void SetLowestEnergyLimit(G4double value)
void clearAndDestroy()
Hep3Vector boostVector() const
Hep3Vector vect() const
void PutValue(size_t index, G4double energy, G4double dataValue)
void SetLocalEnergyDeposit(G4double aE)
Float_t x2[n_points_geant4]
Definition: compare.C:26
static constexpr double GeV
Definition: G4SIunits.hh:217
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
HepLorentzVector & boost(double, double, double)