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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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