Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4LEHadronProtonElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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