Geant4  10.01.p02
G4LEnp.cc
Go to the documentation of this file.
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 
27 // G4 Low energy model: n-p scattering
28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch
29 
30 // 11-OCT-2007 F.W. Jones: removed erroneous code for identity
31 // exchange of particles.
32 // FWJ 27-AUG-2010: extended to 5 GeV by Tony Kwan TRIUMF
33 
34 #include "G4LEnp.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 #include "G4ios.hh"
39 
40 // Initialization of static data arrays:
41 #include "G4LEnpData.hh"
42 #include "Randomize.hh"
43 
44 
46  G4HadronElastic("G4LEnp") // G4HadronicInteraction("G4LEnp")
47 {
48  // theParticleChange.SetNumberOfSecondaries(1);
49 
50  // SetMinEnergy(10.*MeV);
51  // SetMaxEnergy(1200.*MeV);
52  SetMinEnergy(0.);
53  SetMaxEnergy(5.*GeV);
54 }
55 
57 {
59 }
60 
62 G4LEnp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
63 {
65  const G4HadProjectile* aParticle = &aTrack;
66 
67  G4double P = aParticle->GetTotalMomentum();
68  G4double Px = aParticle->Get4Momentum().x();
69  G4double Py = aParticle->Get4Momentum().y();
70  G4double Pz = aParticle->Get4Momentum().z();
71  G4double ek = aParticle->GetKineticEnergy();
72  G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
73 
74  if (verboseLevel > 1) {
75  G4double E = aParticle->GetTotalEnergy();
76  G4double E0 = aParticle->GetDefinition()->GetPDGMass();
77  G4double Q = aParticle->GetDefinition()->GetPDGCharge();
78  G4int A = targetNucleus.GetA_asInt();
79  G4int Z = targetNucleus.GetZ_asInt();
80  G4cout << "G4LEnp:ApplyYourself: incident particle: "
81  << aParticle->GetDefinition()->GetParticleName() << G4endl;
82  G4cout << "P = " << P/GeV << " GeV/c"
83  << ", Px = " << Px/GeV << " GeV/c"
84  << ", Py = " << Py/GeV << " GeV/c"
85  << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
86  G4cout << "E = " << E/GeV << " GeV"
87  << ", kinetic energy = " << ek/GeV << " GeV"
88  << ", mass = " << E0/GeV << " GeV"
89  << ", charge = " << Q << G4endl;
90  G4cout << "G4LEnp:ApplyYourself: material:" << G4endl;
91  G4cout << "A = " << A
92  << ", Z = " << Z
93  << ", atomic mass "
94  << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
95  << G4endl;
96  //
97  // GHEISHA ADD operation to get total energy, mass, charge
98  //
99  E += proton_mass_c2;
100  G4double E02 = E*E - P*P;
101  E0 = std::sqrt(std::abs(E02));
102  if (E02 < 0)E0 *= -1;
103  Q += Z;
104  G4cout << "G4LEnp:ApplyYourself: total:" << G4endl;
105  G4cout << "E = " << E/GeV << " GeV"
106  << ", mass = " << E0/GeV << " GeV"
107  << ", charge = " << Q << G4endl;
108  }
109 
110  // Find energy bin
111 
112  G4int je1 = 0;
113  G4int je2 = NENERGY - 1;
114  ek = ek/GeV;
115  do {
116  G4int midBin = (je1 + je2)/2;
117  if (ek < elab[midBin])
118  je2 = midBin;
119  else
120  je1 = midBin;
121  } while (je2 - je1 > 1);
122  G4double delab = elab[je2] - elab[je1];
123 
124  // Sample the angle
125 
126  G4float sample = G4UniformRand();
127  G4int ke1 = 0;
128  G4int ke2 = NANGLE - 1;
129  G4double dsig = sig[je2][0] - sig[je1][0];
130  G4double rc = dsig/delab;
131  G4double b = sig[je1][0] - rc*elab[je1];
132  G4double sigint1 = rc*ek + b;
133  G4double sigint2 = 0.;
134 
135  if (verboseLevel > 1) {
136  G4cout << "sample=" << sample << G4endl
137  << ke1 << " " << ke2 << " "
138  << sigint1 << " " << sigint2 << G4endl;
139  }
140  do {
141  G4int midBin = (ke1 + ke2)/2;
142  dsig = sig[je2][midBin] - sig[je1][midBin];
143  rc = dsig/delab;
144  b = sig[je1][midBin] - rc*elab[je1];
145  G4double sigint = rc*ek + b;
146  if (sample < sigint) {
147  ke2 = midBin;
148  sigint2 = sigint;
149  }
150  else {
151  ke1 = midBin;
152  sigint1 = sigint;
153  }
154  if (verboseLevel > 1) {
155  G4cout << ke1 << " " << ke2 << " "
156  << sigint1 << " " << sigint2 << G4endl;
157  }
158  } while (ke2 - ke1 > 1);
159 
160  dsig = sigint2 - sigint1;
161  rc = 1./dsig;
162  b = ke1 - rc*sigint1;
163  G4double kint = rc*sample + b;
164  G4double theta = (0.5 + kint)*pi/180.;
165 
166  if (verboseLevel > 1) {
167  G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
168  G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
169  }
170 
171  // Get the target particle
172 
173  G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
174 
175  G4double E1 = aParticle->GetTotalEnergy();
176  G4double M1 = aParticle->GetDefinition()->GetPDGMass();
177  G4double E2 = targetParticle->GetTotalEnergy();
178  G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
179  G4double totalEnergy = E1 + E2;
180  G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
181 
182  // Transform into centre of mass system
183 
184  G4double px = (M2/pseudoMass)*Px;
185  G4double py = (M2/pseudoMass)*Py;
186  G4double pz = (M2/pseudoMass)*Pz;
187  G4double p = std::sqrt(px*px + py*py + pz*pz);
188 
189  if (verboseLevel > 1) {
190  G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
191  G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
192  G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
193  << pz/GeV << " " << p/GeV << G4endl;
194  }
195 
196  // First scatter w.r.t. Z axis
197  G4double phi = G4UniformRand()*twopi;
198  G4double pxnew = p*std::sin(theta)*std::cos(phi);
199  G4double pynew = p*std::sin(theta)*std::sin(phi);
200  G4double pznew = p*std::cos(theta);
201 
202  // Rotate according to the direction of the incident particle
203  if (px*px + py*py > 0) {
204  G4double cost, sint, ph, cosp, sinp;
205  cost = pz/p;
206  sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
207  py < 0 ? ph = 3*halfpi : ph = halfpi;
208  if (std::abs(px) > 0.000001*GeV) ph = std::atan2(py,px);
209  cosp = std::cos(ph);
210  sinp = std::sin(ph);
211  px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
212  py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
213  pz = (-sint*pxnew + cost*pznew);
214  }
215  else {
216  px = pxnew;
217  py = pynew;
218  pz = pznew;
219  }
220 
221  if (verboseLevel > 1) {
222  G4cout << " AFTER SCATTER..." << G4endl;
223  G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
224  << pz/GeV << " " << p/GeV << G4endl;
225  }
226 
227  // Transform to lab system
228 
229  G4double E1pM2 = E1 + M2;
230  G4double betaCM = P/E1pM2;
231  G4double betaCMx = Px/E1pM2;
232  G4double betaCMy = Py/E1pM2;
233  G4double betaCMz = Pz/E1pM2;
234  G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
235 
236  if (verboseLevel > 1) {
237  G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
238  << betaCMz << " " << betaCM << G4endl;
239  G4cout << " gammaCM " << gammaCM << G4endl;
240  }
241 
242  // Now following GLOREN...
243 
244  G4double BETA[5], PA[5], PB[5];
245  BETA[1] = -betaCMx;
246  BETA[2] = -betaCMy;
247  BETA[3] = -betaCMz;
248  BETA[4] = gammaCM;
249 
250  //The incident particle...
251 
252  PA[1] = px;
253  PA[2] = py;
254  PA[3] = pz;
255  PA[4] = std::sqrt(M1*M1 + p*p);
256 
257  G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
258  G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
259 
260  PB[1] = PA[1] + BPGAM * BETA[1];
261  PB[2] = PA[2] + BPGAM * BETA[2];
262  PB[3] = PA[3] + BPGAM * BETA[3];
263  PB[4] = (PA[4] - BETPA) * BETA[4];
264 
266  newP->SetDefinition(aParticle->GetDefinition());
267  newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
268 
269  //The target particle...
270 
271  PA[1] = -px;
272  PA[2] = -py;
273  PA[3] = -pz;
274  PA[4] = std::sqrt(M2*M2 + p*p);
275 
276  BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
277  BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
278 
279  PB[1] = PA[1] + BPGAM * BETA[1];
280  PB[2] = PA[2] + BPGAM * BETA[2];
281  PB[3] = PA[3] + BPGAM * BETA[3];
282  PB[4] = (PA[4] - BETPA) * BETA[4];
283 
284  targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
285 
286  if (verboseLevel > 1) {
287  G4cout << " particle 1 momentum in LAB "
288  << newP->GetMomentum()*(1./GeV)
289  << " " << newP->GetTotalMomentum()/GeV << G4endl;
290  G4cout << " particle 2 momentum in LAB "
291  << targetParticle->GetMomentum()*(1./GeV)
292  << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
293  G4cout << " TOTAL momentum in LAB "
294  << (newP->GetMomentum()+targetParticle->GetMomentum())*(1./GeV)
295  << " "
296  << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
297  << G4endl;
298  }
299 
302  delete newP;
303  theParticleChange.AddSecondary(targetParticle);
304 
305  return &theParticleChange;
306 }
307 
309 //
310 // sample momentum transfer using Lab. momentum
311 
313  G4double plab, G4int , G4int )
314 {
315  G4double nMass = p->GetPDGMass(); // 939.565346*MeV;
316  G4double ek = std::sqrt(plab*plab+nMass*nMass) - nMass;
317 
318  // Find energy bin
319 
320  G4int je1 = 0;
321  G4int je2 = NENERGY - 1;
322  ek = ek/GeV;
323 
324  do
325  {
326  G4int midBin = (je1 + je2)/2;
327  if (ek < elab[midBin])
328  je2 = midBin;
329  else
330  je1 = midBin;
331  } while (je2 - je1 > 1);
332 
333  G4double delab = elab[je2] - elab[je1];
334 
335  // Sample the angle
336 
337  G4double sample = G4UniformRand();
338  G4int ke1 = 0;
339  G4int ke2 = NANGLE - 1;
340  G4double dsig = sig[je2][0] - sig[je1][0];
341  G4double rc = dsig/delab;
342  G4double b = sig[je1][0] - rc*elab[je1];
343  G4double sigint1 = rc*ek + b;
344  G4double sigint2 = 0.;
345 
346  do
347  {
348  G4int midBin = (ke1 + ke2)/2;
349  dsig = sig[je2][midBin] - sig[je1][midBin];
350  rc = dsig/delab;
351  b = sig[je1][midBin] - rc*elab[je1];
352  G4double sigint = rc*ek + b;
353 
354  if (sample < sigint)
355  {
356  ke2 = midBin;
357  sigint2 = sigint;
358  }
359  else
360  {
361  ke1 = midBin;
362  sigint1 = sigint;
363  }
364  } while (ke2 - ke1 > 1);
365 
366  dsig = sigint2 - sigint1;
367  rc = 1./dsig;
368  b = ke1 - rc*sigint1;
369 
370  G4double kint = rc*sample + b;
371  G4double theta = (0.5 + kint)*pi/180.;
372  G4double t = 0.5*plab*plab*(1-std::cos(theta));
373 
374  return t;
375 }
376 
377  // end of file
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4LEnp.cc:312
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
void SetMomentum(const G4ThreeVector &momentum)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetTotalEnergy() const
~G4LEnp()
Definition: G4LEnp.cc:56
const G4double pi
float G4float
Definition: G4Types.hh:77
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:227
const G4String & GetParticleName() const
G4double GetTotalMomentum() const
void SetMinEnergy(G4double anEnergy)
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static const double GeV
Definition: G4SIunits.hh:196
static const G4float elab[NENERGY]
Definition: G4LEnp.hh:86
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
void SetEnergyChange(G4double anEnergy)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4LEnp.cc:62
G4double GetPDGMass() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
static const double degree
Definition: G4SIunits.hh:125
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double GetPDGCharge() const
void SetMomentumChange(const G4ThreeVector &aV)
static const G4float sig[NENERGY][NANGLE]
Definition: G4LEnp.hh:85
G4LEnp()
Definition: G4LEnp.cc:45
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
G4double GetTotalEnergy() const