| Geant4 Cross Reference |
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 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 080612 bug fix contribution from Benoit Pir 31 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
32 // 110505 protection for object is created but 32 // 110505 protection for object is created but not initialized
33 // 110510 delete above protection with more co 33 // 110510 delete above protection with more coordinated work to other classes
34 // 34 //
35 // P. Arce, June-2014 Conversion neutron_hp to 35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 // 36 //
37 #include "G4ParticleHPAngular.hh" 37 #include "G4ParticleHPAngular.hh"
38 38
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 41
42 void G4ParticleHPAngular::Init(std::istream& a 42 void G4ParticleHPAngular::Init(std::istream& aDataFile)
43 { 43 {
44 // G4cout << "here we are entering the Angu 44 // G4cout << "here we are entering the Angular Init"<<G4endl;
45 aDataFile >> theAngularDistributionType >> t 45 aDataFile >> theAngularDistributionType >> targetMass;
46 aDataFile >> frameFlag; 46 aDataFile >> frameFlag;
47 if (theAngularDistributionType == 0) { 47 if (theAngularDistributionType == 0) {
48 theIsoFlag = true; 48 theIsoFlag = true;
49 } 49 }
50 else if (theAngularDistributionType == 1) { 50 else if (theAngularDistributionType == 1) {
51 theIsoFlag = false; 51 theIsoFlag = false;
52 G4int nEnergy; 52 G4int nEnergy;
53 aDataFile >> nEnergy; 53 aDataFile >> nEnergy;
54 theCoefficients = new G4ParticleHPLegendre 54 theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
55 theCoefficients->InitInterpolation(aDataFi 55 theCoefficients->InitInterpolation(aDataFile);
56 G4double temp, energy; 56 G4double temp, energy;
57 G4int tempdep, nLegendre; 57 G4int tempdep, nLegendre;
58 G4int i, ii; 58 G4int i, ii;
59 for (i = 0; i < nEnergy; i++) { 59 for (i = 0; i < nEnergy; i++) {
60 aDataFile >> temp >> energy >> tempdep > 60 aDataFile >> temp >> energy >> tempdep >> nLegendre;
61 energy *= eV; 61 energy *= eV;
62 theCoefficients->Init(i, energy, nLegend 62 theCoefficients->Init(i, energy, nLegendre);
63 theCoefficients->SetTemperature(i, temp) 63 theCoefficients->SetTemperature(i, temp);
64 G4double coeff = 0; 64 G4double coeff = 0;
65 for (ii = 0; ii < nLegendre; ii++) { 65 for (ii = 0; ii < nLegendre; ii++) {
66 aDataFile >> coeff; 66 aDataFile >> coeff;
67 theCoefficients->SetCoeff(i, ii + 1, c 67 theCoefficients->SetCoeff(i, ii + 1, coeff);
68 } 68 }
69 } 69 }
70 } 70 }
71 else if (theAngularDistributionType == 2) { 71 else if (theAngularDistributionType == 2) {
72 theIsoFlag = false; 72 theIsoFlag = false;
73 G4int nEnergy; 73 G4int nEnergy;
74 aDataFile >> nEnergy; 74 aDataFile >> nEnergy;
75 theProbArray = new G4ParticleHPPartial(nEn 75 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
76 theProbArray->InitInterpolation(aDataFile) 76 theProbArray->InitInterpolation(aDataFile);
77 G4double temp, energy; 77 G4double temp, energy;
78 G4int tempdep; 78 G4int tempdep;
79 for (G4int i = 0; i < nEnergy; i++) { 79 for (G4int i = 0; i < nEnergy; i++) {
80 aDataFile >> temp >> energy >> tempdep; 80 aDataFile >> temp >> energy >> tempdep;
81 energy *= eV; 81 energy *= eV;
82 theProbArray->SetT(i, temp); 82 theProbArray->SetT(i, temp);
83 theProbArray->SetX(i, energy); 83 theProbArray->SetX(i, energy);
84 theProbArray->InitData(i, aDataFile); 84 theProbArray->InitData(i, aDataFile);
85 } 85 }
86 } 86 }
87 else { 87 else {
88 theIsoFlag = false; 88 theIsoFlag = false;
89 G4cout << "unknown distribution found for 89 G4cout << "unknown distribution found for Angular: " << theAngularDistributionType << G4endl;
90 throw G4HadronicException(__FILE__, __LINE 90 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
91 } 91 }
92 } 92 }
93 93
94 void G4ParticleHPAngular::SampleAndUpdate(G4Re 94 void G4ParticleHPAngular::SampleAndUpdate(G4ReactionProduct& aHadron)
95 { 95 {
96 //****************************************** 96 //********************************************************************
97 // EMendoza -> sampling can be isotropic in 97 // EMendoza -> sampling can be isotropic in LAB or in CMS
98 /* 98 /*
99 if(theIsoFlag) 99 if(theIsoFlag)
100 { 100 {
101 // G4cout << "Angular result "<<aHadron.GetTo 101 // G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
102 // @@@ add code for isotropic emission in CMS. 102 // @@@ add code for isotropic emission in CMS.
103 G4double costheta = 2.*G4UniformRand()-1 103 G4double costheta = 2.*G4UniformRand()-1;
104 G4double theta = std::acos(costheta); 104 G4double theta = std::acos(costheta);
105 G4double phi = twopi*G4UniformRand(); 105 G4double phi = twopi*G4UniformRand();
106 G4double sinth = std::sin(theta); 106 G4double sinth = std::sin(theta);
107 G4double en = aHadron.GetTotalMomentum() 107 G4double en = aHadron.GetTotalMomentum();
108 G4ThreeVector temp(en*sinth*std::cos(phi 108 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
109 aHadron.SetMomentum( temp ); 109 aHadron.SetMomentum( temp );
110 aHadron.Lorentz(aHadron, -1.*theTarget); 110 aHadron.Lorentz(aHadron, -1.*theTarget);
111 } 111 }
112 else 112 else
113 { 113 {
114 */ 114 */
115 //****************************************** 115 //********************************************************************
116 if (frameFlag == 1) // LAB 116 if (frameFlag == 1) // LAB
117 { 117 {
118 G4double en = aHadron.GetTotalMomentum(); 118 G4double en = aHadron.GetTotalMomentum();
119 G4ReactionProduct boosted; 119 G4ReactionProduct boosted;
120 boosted.Lorentz(*fCache.Get().theProjectil 120 boosted.Lorentz(*fCache.Get().theProjectileRP, *fCache.Get().theTarget);
121 G4double kineticEnergy = boosted.GetKineti 121 G4double kineticEnergy = boosted.GetKineticEnergy();
122 G4double cosTh = 0.0; 122 G4double cosTh = 0.0;
123 //**************************************** 123 //********************************************************************
124 // EMendoza --> sampling can be also isotr 124 // EMendoza --> sampling can be also isotropic
125 /* 125 /*
126 if(theAngularDistributionType == 1) cosTh 126 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
127 if(theAngularDistributionType == 2) cosTh 127 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
128 */ 128 */
129 //**************************************** 129 //********************************************************************
130 if (theIsoFlag) { 130 if (theIsoFlag) {
131 cosTh = 2. * G4UniformRand() - 1; 131 cosTh = 2. * G4UniformRand() - 1;
132 } 132 }
133 else if (theAngularDistributionType == 1) 133 else if (theAngularDistributionType == 1) {
134 cosTh = theCoefficients->SampleMax(kinet 134 cosTh = theCoefficients->SampleMax(kineticEnergy);
135 } 135 }
136 else if (theAngularDistributionType == 2) 136 else if (theAngularDistributionType == 2) {
137 cosTh = theProbArray->Sample(kineticEner 137 cosTh = theProbArray->Sample(kineticEnergy);
138 } 138 }
139 else { 139 else {
140 G4cout << "unknown distribution found fo 140 G4cout << "unknown distribution found for Angular: " << theAngularDistributionType << G4endl;
141 throw G4HadronicException(__FILE__, __LI 141 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
142 } 142 }
143 //**************************************** 143 //********************************************************************
144 G4double theta = std::acos(cosTh); 144 G4double theta = std::acos(cosTh);
145 G4double phi = twopi * G4UniformRand(); 145 G4double phi = twopi * G4UniformRand();
146 G4double sinth = std::sin(theta); 146 G4double sinth = std::sin(theta);
147 G4ThreeVector temp(en * sinth * std::cos(p 147 G4ThreeVector temp(en * sinth * std::cos(phi), en * sinth * std::sin(phi),
148 en * std::cos(theta)); 148 en * std::cos(theta));
149 aHadron.SetMomentum(temp); 149 aHadron.SetMomentum(temp);
150 } 150 }
151 else if (frameFlag == 2) // costh in CMS 151 else if (frameFlag == 2) // costh in CMS
152 { 152 {
153 G4ReactionProduct boostedN; 153 G4ReactionProduct boostedN;
154 boostedN.Lorentz(*fCache.Get().theProjecti 154 boostedN.Lorentz(*fCache.Get().theProjectileRP, *fCache.Get().theTarget);
155 G4double kineticEnergy = boostedN.GetKinet 155 G4double kineticEnergy = boostedN.GetKineticEnergy();
156 156
157 G4double cosTh = 0.0; 157 G4double cosTh = 0.0;
158 //**************************************** 158 //********************************************************************
159 // EMendoza --> sampling can be also isotr 159 // EMendoza --> sampling can be also isotropic
160 /* 160 /*
161 if(theAngularDistributionType == 1) cosTh 161 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
162 if(theAngularDistributionType == 2) cosTh 162 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
163 */ 163 */
164 //**************************************** 164 //********************************************************************
165 if (theIsoFlag) { 165 if (theIsoFlag) {
166 cosTh = 2. * G4UniformRand() - 1; 166 cosTh = 2. * G4UniformRand() - 1;
167 } 167 }
168 else if (theAngularDistributionType == 1) 168 else if (theAngularDistributionType == 1) {
169 cosTh = theCoefficients->SampleMax(kinet 169 cosTh = theCoefficients->SampleMax(kineticEnergy);
170 } 170 }
171 else if (theAngularDistributionType == 2) 171 else if (theAngularDistributionType == 2) {
172 cosTh = theProbArray->Sample(kineticEner 172 cosTh = theProbArray->Sample(kineticEnergy);
173 } 173 }
174 else { 174 else {
175 G4cout << "unknown distribution found fo 175 G4cout << "unknown distribution found for Angular: " << theAngularDistributionType << G4endl;
176 throw G4HadronicException(__FILE__, __LI 176 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
177 } 177 }
178 //**************************************** 178 //********************************************************************
179 // 080612TK bug fix contribution from Beno 179 // 080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
180 /* 180 /*
181 if(theAngularDistributionType == 1) // 181 if(theAngularDistributionType == 1) // LAB
182 { 182 {
183 G4double en = aHadron.GetTotalMoment 183 G4double en = aHadron.GetTotalMomentum();
184 G4ReactionProduct boosted; 184 G4ReactionProduct boosted;
185 boosted.Lorentz(theProjectile, theTa 185 boosted.Lorentz(theProjectile, theTarget);
186 G4double kineticEnergy = boosted.Get 186 G4double kineticEnergy = boosted.GetKineticEnergy();
187 G4double cosTh = theCoefficients->Sa 187 G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
188 G4double theta = std::acos(cosTh); 188 G4double theta = std::acos(cosTh);
189 G4double phi = twopi*G4UniformRand() 189 G4double phi = twopi*G4UniformRand();
190 G4double sinth = std::sin(theta); 190 G4double sinth = std::sin(theta);
191 G4ThreeVector temp(en*sinth*std::cos 191 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
192 aHadron.SetMomentum( temp ); 192 aHadron.SetMomentum( temp );
193 } 193 }
194 else if(theAngularDistributionType == 194 else if(theAngularDistributionType == 2) // costh in CMS {
195 } 195 }
196 */ 196 */
197 197
198 // G4ReactionProduct boostedN; 198 // G4ReactionProduct boostedN;
199 // boostedN.Lorentz(theProjectile, th 199 // boostedN.Lorentz(theProjectile, theTarget);
200 // G4double kineticEnergy = boostedN. 200 // G4double kineticEnergy = boostedN.GetKineticEnergy();
201 // G4double cosTh = theProbArray->Sam 201 // G4double cosTh = theProbArray->Sample(kineticEnergy);
202 202
203 G4double theta = std::acos(cosTh); 203 G4double theta = std::acos(cosTh);
204 G4double phi = twopi * G4UniformRand(); 204 G4double phi = twopi * G4UniformRand();
205 G4double sinth = std::sin(theta); 205 G4double sinth = std::sin(theta);
206 206
207 G4ThreeVector temp(sinth * std::cos(phi), 207 G4ThreeVector temp(sinth * std::cos(phi), sinth * std::sin(phi), std::cos(theta)); // CMS
208 208
209 // 080612TK bug fix contribution from Beno 209 // 080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
210 /* 210 /*
211 G4double en = aHadron.GetTotalEnergy 211 G4double en = aHadron.GetTotalEnergy(); // Target rest
212 212
213 // get trafo from Target rest frame 213 // get trafo from Target rest frame to CMS
214 G4ReactionProduct boostedT; 214 G4ReactionProduct boostedT;
215 boostedT.Lorentz(theTarget, theTarge 215 boostedT.Lorentz(theTarget, theTarget);
216 216
217 G4ThreeVector the3IncidentParticle = 217 G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
218 G4double nEnergy = boostedN.GetTotal 218 G4double nEnergy = boostedN.GetTotalEnergy();
219 G4ThreeVector the3Target = boostedT. 219 G4ThreeVector the3Target = boostedT.GetMomentum();
220 G4double tEnergy = boostedT.GetTotal 220 G4double tEnergy = boostedT.GetTotalEnergy();
221 G4double totE = nEnergy+tEnergy; 221 G4double totE = nEnergy+tEnergy;
222 G4ThreeVector the3trafo = -the3Targe 222 G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
223 G4ReactionProduct trafo; // for tran 223 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
224 trafo.SetMomentum(the3trafo); 224 trafo.SetMomentum(the3trafo);
225 G4double cmsMom = std::sqrt(the3traf 225 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
226 G4double sqrts = std::sqrt((totE-cms 226 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
227 trafo.SetMass(sqrts); 227 trafo.SetMass(sqrts);
228 trafo.SetTotalEnergy(totE); 228 trafo.SetTotalEnergy(totE);
229 229
230 G4double gamma = trafo.GetTotalEnerg 230 G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
231 G4double cosalpha = temp*trafo.GetMo 231 G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
232 G4double fac = cosalpha*trafo.GetTot 232 G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
233 fac*=gamma; 233 fac*=gamma;
234 234
235 G4double mom; 235 G4double mom;
236 // For G4FPE_DEBUG ON 236 // For G4FPE_DEBUG ON
237 G4double mom2 = ( en*fac*en*fac - 237 G4double mom2 = ( en*fac*en*fac -
238 (fac*fac - gamma*gamma) 238 (fac*fac - gamma*gamma)*
239 (en*en - gamma*gamma*aH 239 (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
240 ); 240 );
241 if ( mom2 > 0.0 ) 241 if ( mom2 > 0.0 )
242 mom = std::sqrt( mom2 ); 242 mom = std::sqrt( mom2 );
243 else 243 else
244 mom = 0.0; 244 mom = 0.0;
245 245
246 mom = -en*fac - mom; 246 mom = -en*fac - mom;
247 mom /= (fac*fac-gamma*gamma); 247 mom /= (fac*fac-gamma*gamma);
248 temp = mom*temp; 248 temp = mom*temp;
249 249
250 aHadron.SetMomentum( temp ); // now 250 aHadron.SetMomentum( temp ); // now all in CMS
251 aHadron.SetTotalEnergy( std::sqrt( m 251 aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
252 aHadron.Lorentz(aHadron, trafo); // 252 aHadron.Lorentz(aHadron, trafo); // now in target rest frame
253 */ 253 */
254 // Determination of the hadron kinetic ene 254 // Determination of the hadron kinetic energy in CMS
255 // aHadron.GetKineticEnergy() is actually 255 // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
256 // kineticEnergy is incident neutron kinet 256 // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
257 G4double QValue = aHadron.GetKineticEnergy 257 G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
258 G4double A1 = fCache.Get().theTarget->GetM 258 G4double A1 = fCache.Get().theTarget->GetMass() / boostedN.GetMass();
259 G4double A1prim = aHadron.GetMass() / boos 259 G4double A1prim = aHadron.GetMass() / boostedN.GetMass();
260 G4double kinE = 260 G4double kinE =
261 (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1) 261 (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1) * (A1 * kineticEnergy + (1 + A1) * QValue);
262 G4double totalE = kinE + aHadron.GetMass() 262 G4double totalE = kinE + aHadron.GetMass();
263 G4double mom2 = totalE * totalE - aHadron. 263 G4double mom2 = totalE * totalE - aHadron.GetMass() * aHadron.GetMass();
264 G4double mom; 264 G4double mom;
265 if (mom2 > 0.0) 265 if (mom2 > 0.0)
266 mom = std::sqrt(mom2); 266 mom = std::sqrt(mom2);
267 else 267 else
268 mom = 0.0; 268 mom = 0.0;
269 269
270 aHadron.SetMomentum(mom * temp); // Set m 270 aHadron.SetMomentum(mom * temp); // Set momentum in CMS
271 aHadron.SetKineticEnergy(kinE); // Set ki 271 aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
272 272
273 // get trafo from Target rest frame to CMS 273 // get trafo from Target rest frame to CMS
274 G4ReactionProduct boostedT; 274 G4ReactionProduct boostedT;
275 boostedT.Lorentz(*fCache.Get().theTarget, 275 boostedT.Lorentz(*fCache.Get().theTarget, *fCache.Get().theTarget);
276 276
277 G4ThreeVector the3IncidentParticle = boost 277 G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
278 G4double nEnergy = boostedN.GetTotalEnergy 278 G4double nEnergy = boostedN.GetTotalEnergy();
279 G4ThreeVector the3Target = boostedT.GetMom 279 G4ThreeVector the3Target = boostedT.GetMomentum();
280 G4double tEnergy = boostedT.GetTotalEnergy 280 G4double tEnergy = boostedT.GetTotalEnergy();
281 G4double totE = nEnergy + tEnergy; 281 G4double totE = nEnergy + tEnergy;
282 G4ThreeVector the3trafo = -the3Target - th 282 G4ThreeVector the3trafo = -the3Target - the3IncidentParticle;
283 G4ReactionProduct trafo; // for transform 283 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
284 trafo.SetMomentum(the3trafo); 284 trafo.SetMomentum(the3trafo);
285 G4double cmsMom = std::sqrt(the3trafo * th 285 G4double cmsMom = std::sqrt(the3trafo * the3trafo);
286 G4double sqrts = std::sqrt((totE - cmsMom) 286 G4double sqrts = std::sqrt((totE - cmsMom) * (totE + cmsMom));
287 trafo.SetMass(sqrts); 287 trafo.SetMass(sqrts);
288 trafo.SetTotalEnergy(totE); 288 trafo.SetTotalEnergy(totE);
289 289
290 aHadron.Lorentz(aHadron, trafo); 290 aHadron.Lorentz(aHadron, trafo);
291 } 291 }
292 else { 292 else {
293 throw G4HadronicException(__FILE__, __LINE 293 throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
294 } 294 }
295 aHadron.Lorentz(aHadron, -1. * (*fCache.Get( 295 aHadron.Lorentz(aHadron, -1. * (*fCache.Get().theTarget));
296 // G4cout << aHadron.GetMomentum()<<" "; 296 // G4cout << aHadron.GetMomentum()<<" ";
297 // G4cout << aHadron.GetTotalMomentum()<<G4 297 // G4cout << aHadron.GetTotalMomentum()<<G4endl;
298 } 298 }
299 299