source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc@ 1199

Last change on this file since 1199 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 22.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// neutron_hp -- source file
27// J.P. Wellisch, Nov-1996
28// A prototype of the low energy neutron transport model.
29//
30// 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31// 070606 bug fix and migrate to enable to Partial cases by T. Koi
32// 080603 bug fix for Hadron Hyper News #932 by T. Koi
33// 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34// 080717 bug fix of calculation of residual momentum by T. Koi
35// 080801 protect negative avalable energy by T. Koi
36// introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38// 090514 Fix bug in IC electron emission case
39// Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40//
41#include "G4NeutronHPInelasticCompFS.hh"
42#include "G4Nucleus.hh"
43#include "G4NucleiProperties.hh"
44#include "G4He3.hh"
45#include "G4Alpha.hh"
46#include "G4Electron.hh"
47#include "G4NeutronHPDataUsed.hh"
48#include "G4ParticleTable.hh"
49
50void G4NeutronHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
51{
52 // char the[100] = {""};
53 // std::ostrstream ost(the, 100, std::ios::out);
54 // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
55 // G4String * aName = new G4String(the);
56 // std::ifstream from(*aName, std::ios::in);
57
58 std::ostringstream ost;
59 ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
60 G4String aName = ost.str();
61 std::ifstream from(aName, std::ios::in);
62
63 if(!from) return; // no data found for this isotope
64 // std::ifstream theGammaData(*aName, std::ios::in);
65 std::ifstream theGammaData(aName, std::ios::in);
66
67 theGammas.Init(theGammaData);
68 // delete aName;
69}
70
71void G4NeutronHPInelasticCompFS::Init (G4double A, G4double Z, G4String & dirName, G4String & aFSType)
72{
73
74 gammaPath = "/Inelastic/Gammas/";
75 if(!getenv("G4NEUTRONHPDATA"))
76 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
77 G4String tBase = getenv("G4NEUTRONHPDATA");
78 gammaPath = tBase+gammaPath;
79 G4String tString = dirName;
80 G4bool dbool;
81 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), tString, aFSType, dbool);
82 G4String filename = aFile.GetName();
83 theBaseA = aFile.GetA();
84 theBaseZ = aFile.GetZ();
85 theNDLDataA = (int)aFile.GetA();
86 theNDLDataZ = aFile.GetZ();
87 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
88 {
89 if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
90 hasAnyData = false;
91 hasFSData = false;
92 hasXsec = false;
93 return;
94 }
95 theBaseA = A;
96 theBaseZ = G4int(Z+.5);
97 std::ifstream theData(filename, std::ios::in);
98 if(!theData)
99 {
100 hasAnyData = false;
101 hasFSData = false;
102 hasXsec = false;
103 theData.close();
104 return;
105 }
106 // here we go
107 G4int infoType, dataType, dummy;
108 G4int sfType, it;
109 hasFSData = false;
110 while (theData >> infoType)
111 {
112 hasFSData = true;
113 theData >> dataType;
114 theData >> sfType >> dummy;
115 it = 50;
116 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
117 if(dataType==3)
118 {
119 theData >> dummy >> dummy;
120 theXsection[it] = new G4NeutronHPVector;
121 G4int total;
122 theData >> total;
123 theXsection[it]->Init(theData, total, eV);
124 //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
125 }
126 else if(dataType==4)
127 {
128 theAngularDistribution[it] = new G4NeutronHPAngular;
129 theAngularDistribution[it]->Init(theData);
130 }
131 else if(dataType==5)
132 {
133 theEnergyDistribution[it] = new G4NeutronHPEnergyDistribution;
134 theEnergyDistribution[it]->Init(theData);
135 }
136 else if(dataType==6)
137 {
138 theEnergyAngData[it] = new G4NeutronHPEnAngCorrelation;
139 theEnergyAngData[it]->Init(theData);
140 }
141 else if(dataType==12)
142 {
143 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
144 theFinalStatePhotons[it]->InitMean(theData);
145 }
146 else if(dataType==13)
147 {
148 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
149 theFinalStatePhotons[it]->InitPartials(theData);
150 }
151 else if(dataType==14)
152 {
153 theFinalStatePhotons[it]->InitAngular(theData);
154 }
155 else if(dataType==15)
156 {
157 theFinalStatePhotons[it]->InitEnergies(theData);
158 }
159 else
160 {
161 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS");
162 }
163 }
164 theData.close();
165}
166
167G4int G4NeutronHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
168{
169
170// it = 0 has without Photon
171 G4double running[50];
172 running[0] = 0;
173 unsigned int i;
174 for(i=0; i<50; i++)
175 {
176 if(i!=0) running[i]=running[i-1];
177 if(theXsection[i] != 0)
178 {
179 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
180 }
181 }
182 G4double random = G4UniformRand();
183 G4double sum = running[49];
184 G4int it = 50;
185 if(0!=sum)
186 {
187 G4int i0;
188 for(i0=0; i0<50; i0++)
189 {
190 it = i0;
191 if(random < running[i0]/sum) break;
192 }
193 }
194//debug: it = 1;
195 return it;
196}
197
198
199 //n,p,d,t,he3,a
200void G4NeutronHPInelasticCompFS::CompositeApply(const G4HadProjectile & theTrack, G4ParticleDefinition * aDefinition)
201{
202
203// prepare neutron
204 theResult.Clear();
205 G4double eKinetic = theTrack.GetKineticEnergy();
206 const G4HadProjectile *incidentParticle = &theTrack;
207 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
208 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
209 theNeutron.SetKineticEnergy( eKinetic );
210
211// prepare target
212 G4int i;
213 for(i=0; i<50; i++)
214 { if(theXsection[i] != 0) { break; } }
215
216 G4double targetMass=0;
217 G4double eps = 0.0001;
218 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
219 G4Neutron::Neutron()->GetPDGMass();
220// if(theEnergyAngData[i]!=0)
221// targetMass = theEnergyAngData[i]->GetTargetMass();
222// else if(theAngularDistribution[i]!=0)
223// targetMass = theAngularDistribution[i]->GetTargetMass();
224// else if(theFinalStatePhotons[50]!=0)
225// targetMass = theFinalStatePhotons[50]->GetTargetMass();
226 G4Nucleus aNucleus;
227 G4ReactionProduct theTarget;
228 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
229 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
230
231// prepare the residual mass
232 G4double residualMass=0;
233 G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
234 G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
235 residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
236 G4Neutron::Neutron()->GetPDGMass();
237
238// prepare energy in target rest frame
239 G4ReactionProduct boosted;
240 boosted.Lorentz(theNeutron, theTarget);
241 eKinetic = boosted.GetKineticEnergy();
242// G4double momentumInCMS = boosted.GetTotalMomentum();
243
244// select exit channel for composite FS class.
245 G4int it = SelectExitChannel( eKinetic );
246
247// set target and neutron in the relevant exit channel
248 InitDistributionInitialState(theNeutron, theTarget, it);
249
250 G4ReactionProductVector * thePhotons = 0;
251 G4ReactionProductVector * theParticles = 0;
252 G4ReactionProduct aHadron;
253 aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
254 G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
255 (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
256//080730c
257 if ( availableEnergy < 0 )
258 {
259 //G4cout << "080730c Adjust availavleEnergy " << G4endl;
260 availableEnergy = 0;
261 }
262 G4int nothingWasKnownOnHadron = 0;
263 G4int dummy;
264 G4double eGamm = 0;
265 G4int iLevel=it-1;
266
267// TK without photon has it = 0
268 if( 50 == it )
269 {
270
271// TK Excitation level is not determined
272 iLevel=-1;
273 aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
274 (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
275
276 aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
277 std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
278 aHadron.GetMass()*aHadron.GetMass()));
279
280/*
281 G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-aHadron.GetMass()*aHadron.GetMass() );
282 G4double p = 0.0;
283 if ( p2 > 0.0 )
284 {
285 p = std::sqrt( p );
286 }
287 aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
288*/
289
290 }
291 else
292 {
293 while( iLevel!=-1 && theGammas.GetLevel(iLevel)==0 ) { iLevel--; }
294 }
295
296
297 if ( theAngularDistribution[it] != 0 ) // MF4
298 {
299 if(theEnergyDistribution[it]!=0) // MF5
300 {
301 aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
302 G4double eSecN = aHadron.GetKineticEnergy();
303 eGamm = eKinetic-eSecN;
304 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
305 {
306 if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
307 }
308 G4double random = 2*G4UniformRand();
309 iLevel+=G4int(random);
310 if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
311 }
312 else
313 {
314 G4double eExcitation = 0;
315 if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
316 while (eKinetic-eExcitation < 0 && iLevel>0)
317 {
318 iLevel--;
319 eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
320 }
321
322 if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0)
323 {
324 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
325 }
326 if(eKinetic-eExcitation < 0) eExcitation = 0;
327 if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
328
329 }
330 theAngularDistribution[it]->SampleAndUpdate(aHadron);
331
332 if( theFinalStatePhotons[it] == 0 )
333 {
334// TK comment Most n,n* eneter to this
335 thePhotons = theGammas.GetDecayGammas(iLevel);
336 eGamm -= theGammas.GetLevelEnergy(iLevel);
337 if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
338 {
339 G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
340 theRestEnergy->SetDefinition(G4Gamma::Gamma());
341 theRestEnergy->SetKineticEnergy(eGamm);
342 G4double costh = 2.*G4UniformRand()-1.;
343 G4double phi = twopi*G4UniformRand();
344 theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
345 eGamm*std::sin(std::acos(costh))*std::sin(phi),
346 eGamm*costh);
347 if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
348 thePhotons->push_back(theRestEnergy);
349 }
350 }
351 }
352 else if(theEnergyAngData[it] != 0) // MF6
353 {
354 theParticles = theEnergyAngData[it]->Sample(eKinetic);
355 }
356 else
357 {
358 // @@@ what to do, if we have photon data, but no info on the hadron itself
359 nothingWasKnownOnHadron = 1;
360 }
361
362 //G4cout << "theFinalStatePhotons it " << it << G4endl;
363 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
364 //G4cout << "theFinalStatePhotons it " << it << G4endl;
365 //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
366 //G4cout << "thePhotons " << thePhotons << G4endl;
367
368 if ( theFinalStatePhotons[it] != 0 )
369 {
370 // the photon distributions are in the Nucleus rest frame.
371 // TK residual rest frame
372 G4ReactionProduct boosted;
373 boosted.Lorentz(theNeutron, theTarget);
374 G4double anEnergy = boosted.GetKineticEnergy();
375 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
376 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
377 G4double testEnergy = 0;
378 if(thePhotons!=0 && thePhotons->size()!=0)
379 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
380 if(theFinalStatePhotons[it]->NeedsCascade())
381 {
382 while(aBaseEnergy>0.01*keV)
383 {
384 // cascade down the levels
385 G4bool foundMatchingLevel = false;
386 G4int closest = 2;
387 G4double deltaEold = -1;
388 for(G4int i=1; i<it; i++)
389 {
390 if(theFinalStatePhotons[i]!=0)
391 {
392 testEnergy = theFinalStatePhotons[i]->GetLevelEnergy();
393 }
394 else
395 {
396 testEnergy = 0;
397 }
398 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
399 if(deltaE<0.1*keV)
400 {
401 G4ReactionProductVector * theNext =
402 theFinalStatePhotons[i]->GetPhotons(anEnergy);
403 thePhotons->push_back(theNext->operator[](0));
404 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
405 delete theNext;
406 foundMatchingLevel = true;
407 break; // ===>
408 }
409 if(theFinalStatePhotons[i]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
410 {
411 closest = i;
412 deltaEold = deltaE;
413 }
414 } // <=== the break goes here.
415 if(!foundMatchingLevel)
416 {
417 G4ReactionProductVector * theNext =
418 theFinalStatePhotons[closest]->GetPhotons(anEnergy);
419 thePhotons->push_back(theNext->operator[](0));
420 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
421 delete theNext;
422 }
423 }
424 }
425 }
426 unsigned int i0;
427 if(thePhotons!=0)
428 {
429 for(i0=0; i0<thePhotons->size(); i0++)
430 {
431 // back to lab
432 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
433 }
434 }
435 //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
436 if(nothingWasKnownOnHadron)
437 {
438 G4double totalPhotonEnergy = 0;
439 if(thePhotons!=0)
440 {
441 unsigned int nPhotons = thePhotons->size();
442 unsigned int i0;
443 for(i0=0; i0<nPhotons; i0++)
444 {
445 totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy();
446 }
447 }
448 availableEnergy -= totalPhotonEnergy;
449 residualMass += totalPhotonEnergy/G4Neutron::Neutron()->GetPDGMass();
450 aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
451 (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
452 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
453 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
454 G4double Phi = twopi*G4UniformRand();
455 G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta);
456 //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
457 // aHadron.GetMass()*aHadron.GetMass()));
458 G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass();
459
460 G4double p = 0.0;
461 if ( p2 > 0.0 )
462 p = std::sqrt ( p2 );
463
464 aHadron.SetMomentum( Vector*p );
465
466 }
467
468// fill the result
469// Beware - the recoil is not necessarily in the particles...
470// Can be calculated from momentum conservation?
471// The idea is that the particles ar emitted forst, and the gammas only once the
472// recoil is on the residual; assumption is that gammas do not contribute to
473// the recoil.
474// This needs more design @@@
475
476 G4int nSecondaries = 2; // the hadron and the recoil
477 G4bool needsSeparateRecoil = false;
478 G4int totalBaryonNumber = 0;
479 G4int totalCharge = 0;
480 G4ThreeVector totalMomentum(0);
481 if(theParticles != 0)
482 {
483 nSecondaries = theParticles->size();
484 G4ParticleDefinition * aDef;
485 unsigned int i0;
486 for(i0=0; i0<theParticles->size(); i0++)
487 {
488 aDef = theParticles->operator[](i0)->GetDefinition();
489 totalBaryonNumber+=aDef->GetBaryonNumber();
490 totalCharge+=G4int(aDef->GetPDGCharge()+eps);
491 totalMomentum += theParticles->operator[](i0)->GetMomentum();
492 }
493 if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()))
494 {
495 needsSeparateRecoil = true;
496 nSecondaries++;
497 residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()
498 -totalBaryonNumber);
499 residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge()
500 -totalCharge);
501 }
502 }
503
504 G4int nPhotons = 0;
505 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
506 nSecondaries += nPhotons;
507
508 G4DynamicParticle * theSec;
509
510 if( theParticles==0 )
511 {
512 theSec = new G4DynamicParticle;
513 theSec->SetDefinition(aHadron.GetDefinition());
514 theSec->SetMomentum(aHadron.GetMomentum());
515 theResult.AddSecondary(theSec);
516
517 aHadron.Lorentz(aHadron, theTarget);
518 G4ReactionProduct theResidual;
519 theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
520 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
521 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
522
523 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
524 //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
525 G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
526 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
527
528 theResidual.Lorentz(theResidual, -1.*theTarget);
529 G4ThreeVector totalPhotonMomentum(0,0,0);
530 if(thePhotons!=0)
531 {
532 for(i=0; i<nPhotons; i++)
533 {
534 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
535 }
536 }
537 theSec = new G4DynamicParticle;
538 theSec->SetDefinition(theResidual.GetDefinition());
539 theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
540 theResult.AddSecondary(theSec);
541 }
542 else
543 {
544 for(i0=0; i0<theParticles->size(); i0++)
545 {
546 theSec = new G4DynamicParticle;
547 theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
548 theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
549 theResult.AddSecondary(theSec);
550 delete theParticles->operator[](i0);
551 }
552 delete theParticles;
553 if(needsSeparateRecoil && residualZ!=0)
554 {
555 G4ReactionProduct theResidual;
556 theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
557 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
558 G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass();
559 resiualKineticEnergy += totalMomentum*totalMomentum;
560 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
561// cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
562 theResidual.SetKineticEnergy(resiualKineticEnergy);
563
564 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
565 //theResidual.SetMomentum(-1.*totalMomentum);
566 //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
567 //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
568//080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
569 theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
570
571 theSec = new G4DynamicParticle;
572 theSec->SetDefinition(theResidual.GetDefinition());
573 theSec->SetMomentum(theResidual.GetMomentum());
574 theResult.AddSecondary(theSec);
575 }
576 }
577 if(thePhotons!=0)
578 {
579 for(i=0; i<nPhotons; i++)
580 {
581 theSec = new G4DynamicParticle;
582 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
583 //theSec->SetDefinition(G4Gamma::Gamma());
584 theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
585 //But never cause real effect at least with G4NDL3.13 TK
586 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
587 theResult.AddSecondary(theSec);
588 delete thePhotons->operator[](i);
589 }
590// some garbage collection
591 delete thePhotons;
592 }
593
594//080721
595 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
596 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
597 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
598 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
599 adjust_final_state ( init_4p_lab );
600
601// clean up the primary neutron
602 theResult.SetStatusChange(stopAndKill);
603}
Note: See TracBrowser for help on using the repository browser.