Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4LEnp.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 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4LEnp.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4LEnp.cc (Version 10.1.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26                                                    26 
 27 // G4 Low energy model: n-p scattering             27 // G4 Low energy model: n-p scattering
 28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch       28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch
 29                                                    29 
 30 // 11-OCT-2007 F.W. Jones: removed erroneous c     30 // 11-OCT-2007 F.W. Jones: removed erroneous code for identity
 31 //             exchange of particles.              31 //             exchange of particles.
 32 // FWJ 27-AUG-2010: extended to 5 GeV by Tony      32 // FWJ 27-AUG-2010: extended to 5 GeV by Tony Kwan TRIUMF
 33                                                    33 
 34 #include "G4LEnp.hh"                               34 #include "G4LEnp.hh"
 35 #include "G4PhysicalConstants.hh"                  35 #include "G4PhysicalConstants.hh"
 36 #include "G4SystemOfUnits.hh"                      36 #include "G4SystemOfUnits.hh"
 37 #include "Randomize.hh"                            37 #include "Randomize.hh"
 38 #include "G4ios.hh"                                38 #include "G4ios.hh"
 39                                                    39 
 40 // Initialization of static data arrays:           40 // Initialization of static data arrays:
 41 #include "G4LEnpData.hh"                           41 #include "G4LEnpData.hh"
 42 #include "Randomize.hh"                            42 #include "Randomize.hh"
 43                                                    43 
 44 #include "G4PhysicsModelCatalog.hh"            << 
 45                                                << 
 46                                                    44 
 47 G4LEnp::G4LEnp():                                  45 G4LEnp::G4LEnp():
 48  G4HadronElastic("G4LEnp")  // G4HadronicInter     46  G4HadronElastic("G4LEnp")  // G4HadronicInteraction("G4LEnp")
 49 {                                                  47 {
 50   secID = G4PhysicsModelCatalog::GetModelID( " << 
 51   //    theParticleChange.SetNumberOfSecondari     48   //    theParticleChange.SetNumberOfSecondaries(1);
 52                                                    49   
 53   //    SetMinEnergy(10.*MeV);                     50   //    SetMinEnergy(10.*MeV);
 54   //    SetMaxEnergy(1200.*MeV);                   51   //    SetMaxEnergy(1200.*MeV);
 55   SetMinEnergy(0.);                                52   SetMinEnergy(0.);
 56   SetMaxEnergy(5.*GeV);                            53   SetMaxEnergy(5.*GeV);
 57 }                                                  54 }
 58                                                    55 
 59 G4LEnp::~G4LEnp()                                  56 G4LEnp::~G4LEnp()
 60 {                                                  57 {
 61       theParticleChange.Clear();                   58       theParticleChange.Clear();
 62 }                                                  59 }
 63                                                    60 
 64 G4HadFinalState*                                   61 G4HadFinalState*
 65 G4LEnp::ApplyYourself(const G4HadProjectile& a     62 G4LEnp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 66 {                                                  63 {
 67     theParticleChange.Clear();                     64     theParticleChange.Clear();
 68     const G4HadProjectile* aParticle = &aTrack     65     const G4HadProjectile* aParticle = &aTrack;
 69                                                    66 
 70     G4double P = aParticle->GetTotalMomentum()     67     G4double P = aParticle->GetTotalMomentum();
 71     G4double Px = aParticle->Get4Momentum().x(     68     G4double Px = aParticle->Get4Momentum().x();
 72     G4double Py = aParticle->Get4Momentum().y(     69     G4double Py = aParticle->Get4Momentum().y();
 73     G4double Pz = aParticle->Get4Momentum().z(     70     G4double Pz = aParticle->Get4Momentum().z();
 74     G4double ek = aParticle->GetKineticEnergy(     71     G4double ek = aParticle->GetKineticEnergy();
 75     G4ThreeVector theInitial = aParticle->Get4     72     G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
 76                                                    73 
 77     if (verboseLevel > 1) {                        74     if (verboseLevel > 1) {
 78       G4double E = aParticle->GetTotalEnergy()     75       G4double E = aParticle->GetTotalEnergy();
 79       G4double E0 = aParticle->GetDefinition()     76       G4double E0 = aParticle->GetDefinition()->GetPDGMass();
 80       G4double Q = aParticle->GetDefinition()-     77       G4double Q = aParticle->GetDefinition()->GetPDGCharge();
 81       G4int A = targetNucleus.GetA_asInt();        78       G4int A = targetNucleus.GetA_asInt();
 82       G4int Z = targetNucleus.GetZ_asInt();        79       G4int Z = targetNucleus.GetZ_asInt();
 83       G4cout << "G4LEnp:ApplyYourself: inciden     80       G4cout << "G4LEnp:ApplyYourself: incident particle: "
 84              << aParticle->GetDefinition()->Ge     81              << aParticle->GetDefinition()->GetParticleName() << G4endl;
 85       G4cout << "P = " << P/GeV << " GeV/c"        82       G4cout << "P = " << P/GeV << " GeV/c"
 86              << ", Px = " << Px/GeV << " GeV/c     83              << ", Px = " << Px/GeV << " GeV/c"
 87              << ", Py = " << Py/GeV << " GeV/c     84              << ", Py = " << Py/GeV << " GeV/c"
 88              << ", Pz = " << Pz/GeV << " GeV/c     85              << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
 89       G4cout << "E = " << E/GeV << " GeV"          86       G4cout << "E = " << E/GeV << " GeV"
 90              << ", kinetic energy = " << ek/Ge     87              << ", kinetic energy = " << ek/GeV << " GeV"
 91              << ", mass = " << E0/GeV << " GeV     88              << ", mass = " << E0/GeV << " GeV"
 92              << ", charge = " << Q << G4endl;      89              << ", charge = " << Q << G4endl;
 93       G4cout << "G4LEnp:ApplyYourself: materia     90       G4cout << "G4LEnp:ApplyYourself: material:" << G4endl;
 94       G4cout << "A = " << A                        91       G4cout << "A = " << A
 95              << ", Z = " << Z                      92              << ", Z = " << Z
 96              << ", atomic mass "                   93              << ", atomic mass " 
 97              <<  G4Proton::Proton()->GetPDGMas     94              <<  G4Proton::Proton()->GetPDGMass()/GeV << "GeV" 
 98              << G4endl;                            95              << G4endl;
 99       //                                           96       //
100       // GHEISHA ADD operation to get total en     97       // GHEISHA ADD operation to get total energy, mass, charge
101       //                                           98       //
102       E += proton_mass_c2;                         99       E += proton_mass_c2;
103       G4double E02 = E*E - P*P;                   100       G4double E02 = E*E - P*P;
104       E0 = std::sqrt(std::abs(E02));              101       E0 = std::sqrt(std::abs(E02));
105       if (E02 < 0)E0 *= -1;                       102       if (E02 < 0)E0 *= -1;
106       Q += Z;                                     103       Q += Z;
107       G4cout << "G4LEnp:ApplyYourself: total:"    104       G4cout << "G4LEnp:ApplyYourself: total:" << G4endl;
108       G4cout << "E = " << E/GeV << " GeV"         105       G4cout << "E = " << E/GeV << " GeV"
109              << ", mass = " << E0/GeV << " GeV    106              << ", mass = " << E0/GeV << " GeV"
110              << ", charge = " << Q << G4endl;     107              << ", charge = " << Q << G4endl;
111     }                                             108     }
112                                                   109 
113     // Find energy bin                            110     // Find energy bin
114                                                   111 
115     G4int je1 = 0;                                112     G4int je1 = 0;
116     G4int je2 = NENERGY - 1;                      113     G4int je2 = NENERGY - 1;
117     ek = ek/GeV;                                  114     ek = ek/GeV;
118     do {                                          115     do {
119       G4int midBin = (je1 + je2)/2;               116       G4int midBin = (je1 + je2)/2;
120       if (ek < elab[midBin])                      117       if (ek < elab[midBin])
121         je2 = midBin;                             118         je2 = midBin;
122       else                                        119       else
123         je1 = midBin;                             120         je1 = midBin;
124     } while (je2 - je1 > 1);  /* Loop checking << 121     } while (je2 - je1 > 1); 
125     G4double delab = elab[je2] - elab[je1];       122     G4double delab = elab[je2] - elab[je1];
126                                                   123 
127     // Sample the angle                           124     // Sample the angle
128                                                   125 
129     G4double sample = G4UniformRand();         << 126     G4float sample = G4UniformRand();
130     G4int ke1 = 0;                                127     G4int ke1 = 0;
131     G4int ke2 = NANGLE - 1;                       128     G4int ke2 = NANGLE - 1;
132     G4double dsig = sig[je2][0] - sig[je1][0];    129     G4double dsig = sig[je2][0] - sig[je1][0];
133     G4double rc = dsig/delab;                     130     G4double rc = dsig/delab;
134     G4double b = sig[je1][0] - rc*elab[je1];      131     G4double b = sig[je1][0] - rc*elab[je1];
135     G4double sigint1 = rc*ek + b;                 132     G4double sigint1 = rc*ek + b;
136     G4double sigint2 = 0.;                        133     G4double sigint2 = 0.;
137                                                   134 
138     if (verboseLevel > 1) {                       135     if (verboseLevel > 1) {
139       G4cout << "sample=" << sample << G4endl     136       G4cout << "sample=" << sample << G4endl
140        << ke1 << " " << ke2 << " "                137        << ke1 << " " << ke2 << " " 
141        << sigint1 << " " << sigint2 << G4endl;    138        << sigint1 << " " << sigint2 << G4endl;
142     }                                             139     }
143     do {                                          140     do {
144       G4int midBin = (ke1 + ke2)/2;               141       G4int midBin = (ke1 + ke2)/2;
145       dsig = sig[je2][midBin] - sig[je1][midBi    142       dsig = sig[je2][midBin] - sig[je1][midBin];
146       rc = dsig/delab;                            143       rc = dsig/delab;
147       b = sig[je1][midBin] - rc*elab[je1];        144       b = sig[je1][midBin] - rc*elab[je1];
148       G4double sigint = rc*ek + b;                145       G4double sigint = rc*ek + b;
149       if (sample < sigint) {                      146       if (sample < sigint) {
150         ke2 = midBin;                             147         ke2 = midBin;
151         sigint2 = sigint;                         148         sigint2 = sigint;
152       }                                           149       }
153       else {                                      150       else {
154         ke1 = midBin;                             151         ke1 = midBin;
155         sigint1 = sigint;                         152         sigint1 = sigint;
156       }                                           153       }
157       if (verboseLevel > 1) {                     154       if (verboseLevel > 1) {
158   G4cout << ke1 << " " << ke2 << " "              155   G4cout << ke1 << " " << ke2 << " " 
159          << sigint1 << " " << sigint2 << G4end    156          << sigint1 << " " << sigint2 << G4endl;
160       }                                           157       }
161     } while (ke2 - ke1 > 1);  /* Loop checking << 158     } while (ke2 - ke1 > 1); 
162                                                   159 
163     dsig = sigint2 - sigint1;                     160     dsig = sigint2 - sigint1;
164     rc = 1./dsig;                                 161     rc = 1./dsig;
165     b = ke1 - rc*sigint1;                         162     b = ke1 - rc*sigint1;
166     G4double kint = rc*sample + b;                163     G4double kint = rc*sample + b;
167     G4double theta = (0.5 + kint)*pi/180.;        164     G4double theta = (0.5 + kint)*pi/180.;
168                                                   165 
169     if (verboseLevel > 1) {                       166     if (verboseLevel > 1) {
170       G4cout << "   energy bin " << je1 << " e    167       G4cout << "   energy bin " << je1 << " energy=" << elab[je1] << G4endl;
171       G4cout << "   angle bin " << kint << " a    168       G4cout << "   angle bin " << kint << " angle=" << theta/degree << G4endl;
172     }                                             169     }
173                                                   170 
174     // Get the target particle                    171     // Get the target particle
175                                                   172 
176     G4DynamicParticle* targetParticle = target    173     G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
177                                                   174 
178     G4double E1 = aParticle->GetTotalEnergy();    175     G4double E1 = aParticle->GetTotalEnergy();
179     G4double M1 = aParticle->GetDefinition()->    176     G4double M1 = aParticle->GetDefinition()->GetPDGMass();
180     G4double E2 = targetParticle->GetTotalEner    177     G4double E2 = targetParticle->GetTotalEnergy();
181     G4double M2 = targetParticle->GetDefinitio    178     G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
182     G4double totalEnergy = E1 + E2;               179     G4double totalEnergy = E1 + E2;
183     G4double pseudoMass = std::sqrt(totalEnerg    180     G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
184                                                   181 
185     // Transform into centre of mass system       182     // Transform into centre of mass system
186                                                   183 
187     G4double px = (M2/pseudoMass)*Px;             184     G4double px = (M2/pseudoMass)*Px;
188     G4double py = (M2/pseudoMass)*Py;             185     G4double py = (M2/pseudoMass)*Py;
189     G4double pz = (M2/pseudoMass)*Pz;             186     G4double pz = (M2/pseudoMass)*Pz;
190     G4double p = std::sqrt(px*px + py*py + pz*    187     G4double p = std::sqrt(px*px + py*py + pz*pz);
191                                                   188 
192     if (verboseLevel > 1) {                       189     if (verboseLevel > 1) {
193       G4cout << "  E1, M1 (GeV) " << E1/GeV <<    190       G4cout << "  E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
194       G4cout << "  E2, M2 (GeV) " << E2/GeV <<    191       G4cout << "  E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
195       G4cout << "  particle  1 momentum in CM     192       G4cout << "  particle  1 momentum in CM " << px/GeV << " " << py/GeV << " "
196            << pz/GeV << " " << p/GeV << G4endl    193            << pz/GeV << " " << p/GeV << G4endl;
197     }                                             194     }
198                                                   195 
199     // First scatter w.r.t. Z axis                196     // First scatter w.r.t. Z axis
200     G4double phi = G4UniformRand()*twopi;         197     G4double phi = G4UniformRand()*twopi;
201     G4double pxnew = p*std::sin(theta)*std::co    198     G4double pxnew = p*std::sin(theta)*std::cos(phi);
202     G4double pynew = p*std::sin(theta)*std::si    199     G4double pynew = p*std::sin(theta)*std::sin(phi);
203     G4double pznew = p*std::cos(theta);           200     G4double pznew = p*std::cos(theta);
204                                                   201 
205     // Rotate according to the direction of th    202     // Rotate according to the direction of the incident particle
206     if (px*px + py*py > 0) {                      203     if (px*px + py*py > 0) {
207       G4double cost, sint, ph, cosp, sinp;        204       G4double cost, sint, ph, cosp, sinp;
208       cost = pz/p;                                205       cost = pz/p;
209       sint = (std::sqrt(std::fabs((1-cost)*(1+    206       sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
210       py < 0 ? ph = 3*halfpi : ph = halfpi;       207       py < 0 ? ph = 3*halfpi : ph = halfpi;
211       if (std::abs(px) > 0.000001*GeV) ph = st    208       if (std::abs(px) > 0.000001*GeV) ph = std::atan2(py,px);
212       cosp = std::cos(ph);                        209       cosp = std::cos(ph);
213       sinp = std::sin(ph);                        210       sinp = std::sin(ph);
214       px = (cost*cosp*pxnew - sinp*pynew + sin    211       px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
215       py = (cost*sinp*pxnew + cosp*pynew + sin    212       py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
216       pz = (-sint*pxnew                  + cos    213       pz = (-sint*pxnew                  + cost*pznew);
217     }                                             214     }
218     else {                                        215     else {
219       px = pxnew;                                 216       px = pxnew;
220       py = pynew;                                 217       py = pynew;
221       pz = pznew;                                 218       pz = pznew;
222     }                                             219     }
223                                                   220 
224     if (verboseLevel > 1) {                       221     if (verboseLevel > 1) {
225       G4cout << "  AFTER SCATTER..." << G4endl    222       G4cout << "  AFTER SCATTER..." << G4endl;
226       G4cout << "  particle 1 momentum in CM "    223       G4cout << "  particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
227            << pz/GeV << " " << p/GeV << G4endl    224            << pz/GeV << " " << p/GeV << G4endl;
228     }                                             225     }
229                                                   226 
230     // Transform to lab system                    227     // Transform to lab system
231                                                   228 
232     G4double E1pM2 = E1 + M2;                     229     G4double E1pM2 = E1 + M2;
233     G4double betaCM  = P/E1pM2;                   230     G4double betaCM  = P/E1pM2;
234     G4double betaCMx = Px/E1pM2;                  231     G4double betaCMx = Px/E1pM2;
235     G4double betaCMy = Py/E1pM2;                  232     G4double betaCMy = Py/E1pM2;
236     G4double betaCMz = Pz/E1pM2;                  233     G4double betaCMz = Pz/E1pM2;
237     G4double gammaCM = E1pM2/std::sqrt(E1pM2*E    234     G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
238                                                   235 
239     if (verboseLevel > 1) {                       236     if (verboseLevel > 1) {
240       G4cout << "  betaCM " << betaCMx << " "     237       G4cout << "  betaCM " << betaCMx << " " << betaCMy << " "
241              << betaCMz << " " << betaCM << G4    238              << betaCMz << " " << betaCM << G4endl;
242       G4cout << "  gammaCM " << gammaCM << G4e    239       G4cout << "  gammaCM " << gammaCM << G4endl;
243     }                                             240     }
244                                                   241 
245     // Now following GLOREN...                    242     // Now following GLOREN...
246                                                   243 
247     G4double BETA[5], PA[5], PB[5];               244     G4double BETA[5], PA[5], PB[5];
248     BETA[1] = -betaCMx;                           245     BETA[1] = -betaCMx;
249     BETA[2] = -betaCMy;                           246     BETA[2] = -betaCMy;
250     BETA[3] = -betaCMz;                           247     BETA[3] = -betaCMz;
251     BETA[4] = gammaCM;                            248     BETA[4] = gammaCM;
252                                                   249 
253     //The incident particle...                    250     //The incident particle...
254                                                   251 
255     PA[1] = px;                                   252     PA[1] = px;
256     PA[2] = py;                                   253     PA[2] = py;
257     PA[3] = pz;                                   254     PA[3] = pz;
258     PA[4] = std::sqrt(M1*M1 + p*p);               255     PA[4] = std::sqrt(M1*M1 + p*p);
259                                                   256 
260     G4double BETPA  = BETA[1]*PA[1] + BETA[2]*    257     G4double BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
261     G4double BPGAM  = (BETPA * BETA[4]/(BETA[4    258     G4double BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
262                                                   259 
263     PB[1] = PA[1] + BPGAM  * BETA[1];             260     PB[1] = PA[1] + BPGAM  * BETA[1];
264     PB[2] = PA[2] + BPGAM  * BETA[2];             261     PB[2] = PA[2] + BPGAM  * BETA[2];
265     PB[3] = PA[3] + BPGAM  * BETA[3];             262     PB[3] = PA[3] + BPGAM  * BETA[3];
266     PB[4] = (PA[4] - BETPA) * BETA[4];            263     PB[4] = (PA[4] - BETPA) * BETA[4];
267                                                   264 
268     G4DynamicParticle* newP = new G4DynamicPar    265     G4DynamicParticle* newP = new G4DynamicParticle;
269     newP->SetDefinition(aParticle->GetDefiniti    266     newP->SetDefinition(aParticle->GetDefinition());
270     newP->SetMomentum(G4ThreeVector(PB[1], PB[    267     newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
271                                                   268 
272     //The target particle...                      269     //The target particle...
273                                                   270 
274     PA[1] = -px;                                  271     PA[1] = -px;
275     PA[2] = -py;                                  272     PA[2] = -py;
276     PA[3] = -pz;                                  273     PA[3] = -pz;
277     PA[4] = std::sqrt(M2*M2 + p*p);               274     PA[4] = std::sqrt(M2*M2 + p*p);
278                                                   275 
279     BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + B    276     BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
280     BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) -    277     BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
281                                                   278 
282     PB[1] = PA[1] + BPGAM  * BETA[1];             279     PB[1] = PA[1] + BPGAM  * BETA[1];
283     PB[2] = PA[2] + BPGAM  * BETA[2];             280     PB[2] = PA[2] + BPGAM  * BETA[2];
284     PB[3] = PA[3] + BPGAM  * BETA[3];             281     PB[3] = PA[3] + BPGAM  * BETA[3];
285     PB[4] = (PA[4] - BETPA) * BETA[4];            282     PB[4] = (PA[4] - BETPA) * BETA[4];
286                                                   283 
287     targetParticle->SetMomentum(G4ThreeVector(    284     targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
288                                                   285 
289     if (verboseLevel > 1) {                       286     if (verboseLevel > 1) {
290       G4cout << "  particle 1 momentum in LAB     287       G4cout << "  particle 1 momentum in LAB " 
291            << newP->GetMomentum()*(1./GeV)        288            << newP->GetMomentum()*(1./GeV) 
292            << " " << newP->GetTotalMomentum()/    289            << " " << newP->GetTotalMomentum()/GeV << G4endl;
293       G4cout << "  particle 2 momentum in LAB     290       G4cout << "  particle 2 momentum in LAB " 
294            << targetParticle->GetMomentum()*(1    291            << targetParticle->GetMomentum()*(1./GeV) 
295            << " " << targetParticle->GetTotalM    292            << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
296       G4cout << "  TOTAL momentum in LAB "        293       G4cout << "  TOTAL momentum in LAB " 
297            << (newP->GetMomentum()+targetParti    294            << (newP->GetMomentum()+targetParticle->GetMomentum())*(1./GeV) 
298            << " "                                 295            << " " 
299            << (newP->GetMomentum()+targetParti    296            << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
300            << G4endl;                             297            << G4endl;
301     }                                             298     }
302                                                   299 
303     theParticleChange.SetMomentumChange(newP->    300     theParticleChange.SetMomentumChange(newP->GetMomentumDirection());
304     theParticleChange.SetEnergyChange(newP->Ge    301     theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
305     delete newP;                                  302     delete newP;
306     theParticleChange.AddSecondary(targetParti << 303     theParticleChange.AddSecondary(targetParticle);    
307                                                   304 
308     return &theParticleChange;                    305     return &theParticleChange;
309 }                                                 306 }
310                                                   307 
311 //////////////////////////////////////////////    308 ////////////////////////////////////////////////////////////////////
312 //                                                309 //
313 // sample momentum transfer using Lab. momentu    310 // sample momentum transfer using Lab. momentum
314                                                   311 
315 G4double G4LEnp::SampleInvariantT(const G4Part    312 G4double G4LEnp::SampleInvariantT(const G4ParticleDefinition* p, 
316           G4double plab, G4int , G4int )          313           G4double plab, G4int , G4int )
317 {                                                 314 {
318   G4double nMass = p->GetPDGMass(); // 939.565    315   G4double nMass = p->GetPDGMass(); // 939.565346*MeV;
319   G4double ek = std::sqrt(plab*plab+nMass*nMas    316   G4double ek = std::sqrt(plab*plab+nMass*nMass) - nMass;
320                                                   317 
321     // Find energy bin                            318     // Find energy bin
322                                                   319 
323   G4int je1 = 0;                                  320   G4int je1 = 0;
324   G4int je2 = NENERGY - 1;                        321   G4int je2 = NENERGY - 1;
325   ek = ek/GeV;                                    322   ek = ek/GeV;
326                                                   323 
327   do                                              324   do 
328   {                                               325   {
329       G4int midBin = (je1 + je2)/2;               326       G4int midBin = (je1 + je2)/2;
330       if (ek < elab[midBin])                      327       if (ek < elab[midBin])
331         je2 = midBin;                             328         je2 = midBin;
332       else                                        329       else
333         je1 = midBin;                             330         je1 = midBin;
334   } while (je2 - je1 > 1);  /* Loop checking,  << 331   } while (je2 - je1 > 1); 
335                                                   332 
336   G4double delab = elab[je2] - elab[je1];         333   G4double delab = elab[je2] - elab[je1];
337                                                   334 
338     // Sample the angle                           335     // Sample the angle
339                                                   336 
340   G4double sample = G4UniformRand();              337   G4double sample = G4UniformRand();
341   G4int ke1 = 0;                                  338   G4int ke1 = 0;
342   G4int ke2 = NANGLE - 1;                         339   G4int ke2 = NANGLE - 1;
343   G4double dsig = sig[je2][0] - sig[je1][0];      340   G4double dsig = sig[je2][0] - sig[je1][0];
344   G4double rc = dsig/delab;                       341   G4double rc = dsig/delab;
345   G4double b = sig[je1][0] - rc*elab[je1];        342   G4double b = sig[je1][0] - rc*elab[je1];
346   G4double sigint1 = rc*ek + b;                   343   G4double sigint1 = rc*ek + b;
347   G4double sigint2 = 0.;                          344   G4double sigint2 = 0.;
348                                                   345 
349   do                                              346   do
350   {                                               347   {
351       G4int midBin = (ke1 + ke2)/2;               348       G4int midBin = (ke1 + ke2)/2;
352       dsig = sig[je2][midBin] - sig[je1][midBi    349       dsig = sig[je2][midBin] - sig[je1][midBin];
353       rc = dsig/delab;                            350       rc = dsig/delab;
354       b = sig[je1][midBin] - rc*elab[je1];        351       b = sig[je1][midBin] - rc*elab[je1];
355       G4double sigint = rc*ek + b;                352       G4double sigint = rc*ek + b;
356                                                   353 
357       if (sample < sigint)                        354       if (sample < sigint) 
358       {                                           355       {
359         ke2 = midBin;                             356         ke2 = midBin;
360         sigint2 = sigint;                         357         sigint2 = sigint;
361       }                                           358       }
362       else                                        359       else 
363       {                                           360       {
364         ke1 = midBin;                             361         ke1 = midBin;
365         sigint1 = sigint;                         362         sigint1 = sigint;
366       }                                           363       }
367   } while (ke2 - ke1 > 1);  /* Loop checking,  << 364   } while (ke2 - ke1 > 1); 
368                                                   365 
369   dsig = sigint2 - sigint1;                       366   dsig = sigint2 - sigint1;
370   rc = 1./dsig;                                   367   rc = 1./dsig;
371   b = ke1 - rc*sigint1;                           368   b = ke1 - rc*sigint1;
372                                                   369 
373   G4double kint = rc*sample + b;                  370   G4double kint = rc*sample + b;
374   G4double theta = (0.5 + kint)*pi/180.;          371   G4double theta = (0.5 + kint)*pi/180.;
375   G4double t = 0.5*plab*plab*(1-std::cos(theta    372   G4double t = 0.5*plab*plab*(1-std::cos(theta));
376                                                   373 
377   return t;                                       374   return t;
378 }                                                 375 }
379                                                   376 
380  // end of file                                   377  // end of file
381                                                   378