ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEnp.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LEnp.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 
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 
45 G4LEnp::G4LEnp():
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 
56 G4LEnp::~G4LEnp()
57 {
58  theParticleChange.Clear();
59 }
60 
61 G4HadFinalState*
62 G4LEnp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
63 {
64  theParticleChange.Clear();
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); /* Loop checking, 10.08.2015, A.Ribon */
122  G4double delab = elab[je2] - elab[je1];
123 
124  // Sample the angle
125 
126  G4double 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); /* Loop checking, 10.08.2015, A.Ribon */
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 
265  G4DynamicParticle* newP = new G4DynamicParticle;
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 
300  theParticleChange.SetMomentumChange(newP->GetMomentumDirection());
301  theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
302  delete newP;
303  theParticleChange.AddSecondary(targetParticle);
304 
305  return &theParticleChange;
306 }
307 
309 //
310 // sample momentum transfer using Lab. momentum
311 
312 G4double G4LEnp::SampleInvariantT(const G4ParticleDefinition* p,
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); /* Loop checking, 10.08.2015, A.Ribon */
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); /* Loop checking, 10.08.2015, A.Ribon */
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