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