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