ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEpp.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LEpp.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-n or p-p scattering
28  // F.W. Jones, L.G. Greeniaus, H.P. Wellisch
29 
30 // FWJ 27-AUG-2010: extended Coulomb-suppressed data to 5 GeV
31 
32 #include "G4LEpp.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "Randomize.hh"
36 #include "G4ios.hh"
37 
38 // Initialization of static data arrays:
39 #include "G4LEppData.hh"
40 
41 G4LEpp::G4LEpp():G4HadronElastic("G4LEpp")
42 {
43  SetMinEnergy(0.);
44  SetMaxEnergy(5.*GeV);
45 }
46 
47 G4LEpp::~G4LEpp()
48 {}
49 
50 G4HadFinalState*
51 G4LEpp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
52 {
53  theParticleChange.Clear();
54  const G4HadProjectile* aParticle = &aTrack;
55 
56  G4double P = aParticle->GetTotalMomentum();
57  G4double Px = aParticle->Get4Momentum().x();
58  G4double Py = aParticle->Get4Momentum().y();
59  G4double Pz = aParticle->Get4Momentum().z();
60  G4double E = aParticle->GetTotalEnergy();
61  G4ThreeVector theInitial = aParticle->Get4Momentum().vect().unit();
62 
63  if (verboseLevel > 1) {
64  G4double ek = aParticle->GetKineticEnergy();
65  G4double E0 = aParticle->GetDefinition()->GetPDGMass();
66  G4double Q = aParticle->GetDefinition()->GetPDGCharge();
67  G4int A = targetNucleus.GetA_asInt();
68  G4int Z = targetNucleus.GetZ_asInt();
69  G4cout << "G4LEpp:ApplyYourself: incident particle: "
70  << aParticle->GetDefinition()->GetParticleName() << G4endl;
71  G4cout << "P = " << P/GeV << " GeV/c"
72  << ", Px = " << Px/GeV << " GeV/c"
73  << ", Py = " << Py/GeV << " GeV/c"
74  << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
75  G4cout << "E = " << E/GeV << " GeV"
76  << ", kinetic energy = " << ek/GeV << " GeV"
77  << ", mass = " << E0/GeV << " GeV"
78  << ", charge = " << Q << G4endl;
79  G4cout << "G4LEpp:ApplyYourself: material:" << G4endl;
80  G4cout << "A = " << A
81  << ", Z = " << Z
82  << ", atomic mass "
83  << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
84  << G4endl;
85  //
86  // GHEISHA ADD operation to get total energy, mass, charge
87  //
88  E += proton_mass_c2;
89  G4double E02 = E*E - P*P;
90  E0 = std::sqrt(std::fabs(E02));
91  if (E02 < 0)E0 *= -1;
92  Q += Z;
93  G4cout << "G4LEpp:ApplyYourself: total:" << G4endl;
94  G4cout << "E = " << E/GeV << " GeV"
95  << ", mass = " << E0/GeV << " GeV"
96  << ", charge = " << Q << G4endl;
97  }
98  G4double t = SampleInvariantT(aParticle->GetDefinition(), P, 0, 0);
99  G4double cost = 1.0 - 2*t/(P*P);
100  if(cost > 1.0) { cost = 1.0; }
101  if(cost <-1.0) { cost =-1.0; }
102  G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
103  G4double phi = twopi*G4UniformRand();
104  // Get the target particle
105  G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
106 
107  G4double E1 = aParticle->GetTotalEnergy();
108  G4double M1 = aParticle->GetDefinition()->GetPDGMass();
109  G4double E2 = targetParticle->GetTotalEnergy();
110  G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
111  G4double totalEnergy = E1 + E2;
112  G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
113 
114  // Transform into centre of mass system
115 
116  G4double px = (M2/pseudoMass)*Px;
117  G4double py = (M2/pseudoMass)*Py;
118  G4double pz = (M2/pseudoMass)*Pz;
119  G4double p = std::sqrt(px*px + py*py + pz*pz);
120 
121  if (verboseLevel > 1) {
122  G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
123  G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
124  G4cout << " particle 1 momentum in CM " << px/GeV
125  << " " << py/GeV << " "
126  << pz/GeV << " " << p/GeV << G4endl;
127  }
128 
129  // First scatter w.r.t. Z axis
130  G4double pxnew = p*sint*std::cos(phi);
131  G4double pynew = p*sint*std::sin(phi);
132  G4double pznew = p*cost;
133 
134  // Rotate according to the direction of the incident particle
135  if (px*px + py*py > 0) {
136  G4double ph, cosp, sinp;
137  cost = pz/p;
138  sint = (std::sqrt((1-cost)*(1+cost)) + std::sqrt(px*px+py*py)/p)/2;
139  py < 0 ? ph = 3*halfpi : ph = halfpi;
140  if (std::fabs(px) > 0.000001*GeV) ph = std::atan2(py,px);
141  cosp = std::cos(ph);
142  sinp = std::sin(ph);
143  px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
144  py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
145  pz = (-sint*pxnew + cost*pznew);
146  }
147  else {
148  px = pxnew;
149  py = pynew;
150  pz = pznew;
151  }
152 
153  if (verboseLevel > 1) {
154  G4cout << " AFTER SCATTER..." << G4endl;
155  G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
156  << pz/GeV << " " << p/GeV << G4endl;
157  }
158 
159  // Transform to lab system
160 
161  G4double E1pM2 = E1 + M2;
162  G4double betaCM = P/E1pM2;
163  G4double betaCMx = Px/E1pM2;
164  G4double betaCMy = Py/E1pM2;
165  G4double betaCMz = Pz/E1pM2;
166  G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
167 
168  if (verboseLevel > 1) {
169  G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
170  << betaCMz << " " << betaCM << G4endl;
171  G4cout << " gammaCM " << gammaCM << G4endl;
172  }
173 
174  // Now following GLOREN...
175 
176  G4double BETA[5], PA[5], PB[5];
177  BETA[1] = -betaCMx;
178  BETA[2] = -betaCMy;
179  BETA[3] = -betaCMz;
180  BETA[4] = gammaCM;
181 
182  //The incident particle...
183 
184  PA[1] = px;
185  PA[2] = py;
186  PA[3] = pz;
187  PA[4] = std::sqrt(M1*M1 + p*p);
188 
189  G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
190  G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
191 
192  PB[1] = PA[1] + BPGAM * BETA[1];
193  PB[2] = PA[2] + BPGAM * BETA[2];
194  PB[3] = PA[3] + BPGAM * BETA[3];
195  PB[4] = (PA[4] - BETPA) * BETA[4];
196 
197  G4DynamicParticle* newP = new G4DynamicParticle;
198  newP->SetDefinition(aParticle->GetDefinition());
199  newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
200 
201  //The target particle...
202 
203  PA[1] = -px;
204  PA[2] = -py;
205  PA[3] = -pz;
206  PA[4] = std::sqrt(M2*M2 + p*p);
207 
208  BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
209  BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
210 
211  PB[1] = PA[1] + BPGAM * BETA[1];
212  PB[2] = PA[2] + BPGAM * BETA[2];
213  PB[3] = PA[3] + BPGAM * BETA[3];
214  PB[4] = (PA[4] - BETPA) * BETA[4];
215 
216  targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
217 
218  if (verboseLevel > 1) {
219  G4cout << " particle 1 momentum in LAB "
220  << newP->GetMomentum()/GeV
221  << " " << newP->GetTotalMomentum()/GeV << G4endl;
222  G4cout << " particle 2 momentum in LAB "
223  << targetParticle->GetMomentum()/GeV
224  << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
225  G4cout << " TOTAL momentum in LAB "
226  << (newP->GetMomentum()+targetParticle->GetMomentum())/GeV
227  << " "
228  << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
229  << G4endl;
230  }
231 
232  theParticleChange.SetMomentumChange( newP->GetMomentumDirection());
233  theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
234  delete newP;
235 
236  // Recoil particle
237  theParticleChange.AddSecondary(targetParticle);
238  return &theParticleChange;
239 }
240 
242 //
243 // sample momentum transfer using Lab. momentum
244 
245 G4double G4LEpp::SampleInvariantT(const G4ParticleDefinition* p,
246  G4double plab, G4int , G4int )
247 {
248  G4double nMass = p->GetPDGMass(); // 939.565346*MeV;
249  G4double ek = std::sqrt(plab*plab+nMass*nMass) - nMass;
250 
251  // Find energy bin
252 
253  G4int je1 = 0;
254  G4int je2 = NENERGY - 1;
255  ek /= GeV;
256 
257  do
258  {
259  G4int midBin = (je1 + je2)/2;
260 
261  if (ek < elab[midBin]) je2 = midBin;
262  else je1 = midBin;
263  }
264  while (je2 - je1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
265 
266  G4double delab = elab[je2] - elab[je1];
267 
268  // Sample the angle
269 
270  G4double sample = G4UniformRand();
271  G4int ke1 = 0;
272  G4int ke2 = NANGLE - 1;
273  G4double dsig, b, rc;
274 
275  dsig = Sig[je2][0] - Sig[je1][0];
276  rc = dsig/delab;
277  b = Sig[je1][0] - rc*elab[je1];
278 
279  G4double sigint1 = rc*ek + b;
280  G4double sigint2 = 0.;
281 
282  do
283  {
284  G4int midBin = (ke1 + ke2)/2;
285  dsig = Sig[je2][midBin] - Sig[je1][midBin];
286  rc = dsig/delab;
287  b = Sig[je1][midBin] - rc*elab[je1];
288  G4double sigint = rc*ek + b;
289 
290  if (sample < sigint)
291  {
292  ke2 = midBin;
293  sigint2 = sigint;
294  }
295  else
296  {
297  ke1 = midBin;
298  sigint1 = sigint;
299  }
300  }
301  while (ke2 - ke1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
302 
303  dsig = sigint2 - sigint1;
304  rc = 1./dsig;
305  b = ke1 - rc*sigint1;
306 
307  G4double kint = rc*sample + b;
308  G4double theta = (0.5 + kint)*pi/180.;
309  G4double t = 0.5*plab*plab*(1 - std::cos(theta));
310 
311  return t;
312 }
313 // end of file