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 8.2.p1)


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