source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticBaseFS.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 16.9 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// 080801 Give a warning message for irregular mass value in data file by T. Koi
31// Introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
32// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33// 101111 Add Special treatment for Be9(n,2n)Be8(2a) case by T. Koi
34//
35#include "G4NeutronHPInelasticBaseFS.hh"
36#include "G4Nucleus.hh"
37#include "G4NucleiProperties.hh"
38#include "G4He3.hh"
39#include "G4Alpha.hh"
40#include "G4Electron.hh"
41#include "G4NeutronHPDataUsed.hh"
42
43#include "G4ParticleTable.hh"
44
45void G4NeutronHPInelasticBaseFS::InitGammas(G4double AR, G4double ZR)
46{
47 // char the[100] = {""};
48 // std::ostrstream ost(the, 100, std::ios::out);
49 // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
50 // G4String * aName = new G4String(the);
51 // std::ifstream from(*aName, std::ios::in);
52
53 std::ostringstream ost;
54 ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
55 G4String aName = ost.str();
56 std::ifstream from(aName, std::ios::in);
57
58 if(!from) return; // no data found for this isotope
59 // std::ifstream theGammaData(*aName, std::ios::in);
60 std::ifstream theGammaData(aName, std::ios::in);
61
62 G4double eps = 0.001;
63 theNuclearMassDifference =
64 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(AR+eps),static_cast<G4int>(ZR+eps)) -
65 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps));
66 theGammas.Init(theGammaData);
67 // delete aName;
68}
69
70void G4NeutronHPInelasticBaseFS::Init (G4double A, G4double Z, G4String & dirName, G4String & bit)
71{
72 gammaPath = "/Inelastic/Gammas/";
73 if(!getenv("G4NEUTRONHPDATA"))
74 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
75 G4String tBase = getenv("G4NEUTRONHPDATA");
76 gammaPath = tBase+gammaPath;
77 G4String tString = dirName;
78 G4bool dbool;
79 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), tString, bit, dbool);
80 G4String filename = aFile.GetName();
81 theBaseA = aFile.GetA();
82 theBaseZ = aFile.GetZ();
83 theNDLDataA = (int)aFile.GetA();
84 theNDLDataZ = aFile.GetZ();
85 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
86 {
87 if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
88 hasAnyData = false;
89 hasFSData = false;
90 hasXsec = false;
91 return;
92 }
93 theBaseA = A;
94 theBaseZ = G4int(Z+.5);
95 std::ifstream theData(filename, std::ios::in);
96 if(!(theData))
97 {
98 hasAnyData = false;
99 hasFSData = false;
100 hasXsec = false;
101 theData.close();
102 return; // no data for exactly this isotope and FS
103 }
104 // here we go
105 G4int infoType, dataType, dummy=INT_MAX;
106 hasFSData = false;
107 while (theData >> infoType)
108 {
109 theData >> dataType;
110 if(dummy==INT_MAX) theData >> dummy >> dummy;
111 if(dataType==3)
112 {
113 G4int total;
114 theData >> total;
115 theXsection->Init(theData, total, eV);
116 }
117 else if(dataType==4)
118 {
119 theAngularDistribution = new G4NeutronHPAngular;
120 theAngularDistribution->Init(theData);
121 hasFSData = true;
122 }
123 else if(dataType==5)
124 {
125 theEnergyDistribution = new G4NeutronHPEnergyDistribution;
126 theEnergyDistribution->Init(theData);
127 hasFSData = true;
128 }
129 else if(dataType==6)
130 {
131 theEnergyAngData = new G4NeutronHPEnAngCorrelation;
132 theEnergyAngData->Init(theData);
133 hasFSData = true;
134 }
135 else if(dataType==12)
136 {
137 theFinalStatePhotons = new G4NeutronHPPhotonDist;
138 theFinalStatePhotons->InitMean(theData);
139 hasFSData = true;
140 }
141 else if(dataType==13)
142 {
143 theFinalStatePhotons = new G4NeutronHPPhotonDist;
144 theFinalStatePhotons->InitPartials(theData);
145 hasFSData = true;
146 }
147 else if(dataType==14)
148 {
149 theFinalStatePhotons->InitAngular(theData);
150 hasFSData = true;
151 }
152 else if(dataType==15)
153 {
154 theFinalStatePhotons->InitEnergies(theData);
155 hasFSData = true;
156 }
157 else
158 {
159 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticBaseFS");
160 }
161 }
162 theData.close();
163}
164
165void G4NeutronHPInelasticBaseFS::BaseApply(const G4HadProjectile & theTrack,
166 G4ParticleDefinition ** theDefs,
167 G4int nDef)
168{
169
170// prepare neutron
171 theResult.Clear();
172 G4double eKinetic = theTrack.GetKineticEnergy();
173 const G4HadProjectile *incidentParticle = &theTrack;
174 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
175 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
176 theNeutron.SetKineticEnergy( eKinetic );
177
178// prepare target
179 G4double targetMass;
180 G4double eps = 0.0001;
181 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
182 G4Neutron::Neutron()->GetPDGMass();
183 if(theEnergyAngData!=0)
184 { targetMass = theEnergyAngData->GetTargetMass(); }
185 if(theAngularDistribution!=0)
186 { targetMass = theAngularDistribution->GetTargetMass(); }
187//080731a
188if ( targetMass == 0 ) G4cout << "080731a It looks like something wrong value in G4NDL, please update the latest version. If you use the latest, then please report this problem to Geant4 Hyper news." << G4endl;
189 G4Nucleus aNucleus;
190 G4ReactionProduct theTarget;
191 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
192 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
193
194// prepare energy in target rest frame
195 G4ReactionProduct boosted;
196 boosted.Lorentz(theNeutron, theTarget);
197 eKinetic = boosted.GetKineticEnergy();
198 G4double orgMomentum = boosted.GetMomentum().mag();
199
200// Take N-body phase-space distribution, if no other data present.
201 if(!HasFSData()) // adding the residual is trivial here @@@
202 {
203 G4NeutronHPNBodyPhaseSpace thePhaseSpaceDistribution;
204 G4double aPhaseMass=0;
205 G4int ii;
206 for(ii=0; ii<nDef; ii++)
207 {
208 aPhaseMass+=theDefs[ii]->GetPDGMass();
209 }
210 thePhaseSpaceDistribution.Init(aPhaseMass, nDef);
211 thePhaseSpaceDistribution.SetNeutron(&theNeutron);
212 thePhaseSpaceDistribution.SetTarget(&theTarget);
213 for(ii=0; ii<nDef; ii++)
214 {
215 G4double massCode = 1000.*std::abs(theDefs[ii]->GetPDGCharge());
216 massCode += theDefs[ii]->GetBaryonNumber();
217 G4double dummy = 0;
218 G4ReactionProduct * aSec = thePhaseSpaceDistribution.Sample(eKinetic, massCode, dummy);
219 aSec->Lorentz(*aSec, -1.*theTarget);
220 G4DynamicParticle * aPart = new G4DynamicParticle();
221 aPart->SetDefinition(aSec->GetDefinition());
222 aPart->SetMomentum(aSec->GetMomentum());
223 delete aSec;
224 theResult.AddSecondary(aPart);
225 }
226 theResult.SetStatusChange(stopAndKill);
227 return;
228 }
229
230// set target and neutron in the relevant exit channel
231 if(theAngularDistribution!=0)
232 {
233 theAngularDistribution->SetTarget(theTarget);
234 theAngularDistribution->SetNeutron(theNeutron);
235 }
236 else if(theEnergyAngData!=0)
237 {
238 theEnergyAngData->SetTarget(theTarget);
239 theEnergyAngData->SetNeutron(theNeutron);
240 }
241
242 G4ReactionProductVector * tmpHadrons = 0;
243 G4int ii, dummy;
244 unsigned int i;
245 if(theEnergyAngData != 0)
246 {
247 tmpHadrons = theEnergyAngData->Sample(eKinetic);
248 }
249 else if(theAngularDistribution!= 0)
250 {
251 G4bool * Done = new G4bool[nDef];
252 G4int i0;
253 for(i0=0; i0<nDef; i0++) Done[i0] = false;
254 if(tmpHadrons == 0)
255 {
256 tmpHadrons = new G4ReactionProductVector;
257 }
258 else
259 {
260 for(i=0; i<tmpHadrons->size(); i++)
261 {
262 for(ii=0; ii<nDef; ii++)
263 if(!Done[ii] && tmpHadrons->operator[](i)->GetDefinition() == theDefs[ii])
264 Done[ii] = true;
265 }
266 }
267 G4ReactionProduct * aHadron;
268 G4double localMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps)));
269 G4ThreeVector bufferedDirection(0,0,0);
270 for(i0=0; i0<nDef; i0++)
271 {
272 if(!Done[i0])
273 {
274 aHadron = new G4ReactionProduct;
275 if(theEnergyDistribution!=0)
276 {
277 aHadron->SetDefinition(theDefs[i0]);
278 aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
279 }
280 else if(nDef == 1)
281 {
282 aHadron->SetDefinition(theDefs[i0]);
283 aHadron->SetKineticEnergy(eKinetic);
284 }
285 else if(nDef == 2)
286 {
287 aHadron->SetDefinition(theDefs[i0]);
288 aHadron->SetKineticEnergy(50*MeV);
289 }
290 else
291 {
292 throw G4HadronicException(__FILE__, __LINE__, "No energy distribution to sample from in InelasticBaseFS::BaseApply");
293 }
294 theAngularDistribution->SampleAndUpdate(*aHadron);
295 if(theEnergyDistribution==0 && nDef == 2)
296 {
297 if(i0==0)
298 {
299 G4double m1 = theDefs[0]->GetPDGMass();
300 G4double m2 = theDefs[1]->GetPDGMass();
301 G4double mn = G4Neutron::Neutron()->GetPDGMass();
302 G4int z1 = static_cast<G4int>(theBaseZ+eps-theDefs[0]->GetPDGCharge()-theDefs[1]->GetPDGCharge());
303 G4int a1 = static_cast<G4int>(theBaseA+eps)-theDefs[0]->GetBaryonNumber()-theDefs[1]->GetBaryonNumber();
304 G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
305 G4double availableEnergy = eKinetic+mn+localMass-m1-m2-concreteMass;
306 // available kinetic energy in CMS (non relativistic)
307 G4double emin = availableEnergy+m1+m2 - std::sqrt((m1+m2)*(m1+m2)+orgMomentum*orgMomentum);
308 G4double p1=std::sqrt(2.*m2*emin);
309 bufferedDirection = p1*aHadron->GetMomentum().unit();
310 if(getenv("HTOKEN")) // @@@@@ verify the nucleon counting...
311 {
312 G4cout << "HTOKEN "<<z1<<" "<<theBaseZ<<" "<<a1<<" "<<theBaseA<<" "<<availableEnergy<<" "
313 << emin<<G4endl;
314 }
315 }
316 else
317 {
318 bufferedDirection = -bufferedDirection;
319 }
320 // boost from cms to lab
321 if(getenv("HTOKEN"))
322 {
323 G4cout << " HTOKEN "<<bufferedDirection.mag2()<<G4endl;
324 }
325 aHadron->SetTotalEnergy( std::sqrt(aHadron->GetMass()*aHadron->GetMass()
326 +bufferedDirection.mag2()) );
327 aHadron->SetMomentum(bufferedDirection);
328 aHadron->Lorentz(*aHadron, -1.*(theTarget+theNeutron));
329 if(getenv("HTOKEN"))
330 {
331 G4cout << " HTOKEN "<<aHadron->GetTotalEnergy()<<" "<<aHadron->GetMomentum()<<G4endl;
332 }
333 }
334 tmpHadrons->push_back(aHadron);
335 }
336 }
337 delete [] Done;
338 }
339 else
340 {
341 throw G4HadronicException(__FILE__, __LINE__, "No data to create the neutrons in NInelasticFS");
342 }
343
344 G4ReactionProductVector * thePhotons = 0;
345 if(theFinalStatePhotons!=0)
346 {
347 // the photon distributions are in the Nucleus rest frame.
348 G4ReactionProduct boosted;
349 boosted.Lorentz(theNeutron, theTarget);
350 G4double anEnergy = boosted.GetKineticEnergy();
351 thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
352 if(thePhotons!=0)
353 {
354 for(i=0; i<thePhotons->size(); i++)
355 {
356 // back to lab
357 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1.*theTarget);
358 }
359 }
360 }
361 else if(theEnergyAngData!=0)
362 {
363
364 G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
365 G4double anEnergy = boosted.GetKineticEnergy();
366 theGammaEnergy = anEnergy-theGammaEnergy;
367 theGammaEnergy += theNuclearMassDifference;
368 G4double eBindProducts = 0;
369 G4double eBindN = 0;
370 G4double eBindP = 0;
371 G4double eBindD = G4NucleiProperties::GetBindingEnergy(2,1);
372 G4double eBindT = G4NucleiProperties::GetBindingEnergy(3,1);
373 G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3,2);
374 G4double eBindA = G4NucleiProperties::GetBindingEnergy(4,2);
375 G4int ia=0;
376 for(i=0; i<tmpHadrons->size(); i++)
377 {
378 if(tmpHadrons->operator[](i)->GetDefinition() == G4Neutron::Neutron())
379 {
380 eBindProducts+=eBindN;
381 }
382 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Proton::Proton())
383 {
384 eBindProducts+=eBindP;
385 }
386 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Deuteron::Deuteron())
387 {
388 eBindProducts+=eBindD;
389 }
390 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Triton::Triton())
391 {
392 eBindProducts+=eBindT;
393 }
394 else if(tmpHadrons->operator[](i)->GetDefinition() == G4He3::He3())
395 {
396 eBindProducts+=eBindHe3;
397 }
398 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Alpha::Alpha())
399 {
400 eBindProducts+=eBindA;
401 ia++;
402 }
403 }
404
405 theGammaEnergy += eBindProducts;
406
407//101111
408//Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
409if ( (G4int)(theBaseZ+eps) == 4 && (G4int)(theBaseA+eps) == 9 )
410{
411 // This only valid for G4NDL3.13,,,
412 if ( std::abs( theNuclearMassDifference -
413 ( G4NucleiProperties::GetBindingEnergy( 8 , 4 ) -
414 G4NucleiProperties::GetBindingEnergy( 9 , 4 ) ) ) < 1*keV
415 && ia == 2 )
416 {
417 theGammaEnergy -= (2*eBindA);
418 }
419}
420
421 G4ReactionProductVector * theOtherPhotons = 0;
422 G4int iLevel;
423 while(theGammaEnergy>=theGammas.GetLevelEnergy(0))
424 {
425 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
426 {
427 if(theGammas.GetLevelEnergy(iLevel)<theGammaEnergy) break;
428 }
429 if(iLevel==0||iLevel==theGammas.GetNumberOfLevels()-1)
430 {
431 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
432 }
433 else
434 {
435 G4double random = G4UniformRand();
436 G4double eLow = theGammas.GetLevelEnergy(iLevel);
437 G4double eHigh = theGammas.GetLevelEnergy(iLevel+1);
438 if(random > (eHigh-eLow)/(theGammaEnergy-eLow)) iLevel++;
439 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
440 }
441 if(thePhotons==0) thePhotons = new G4ReactionProductVector;
442 if(theOtherPhotons != 0)
443 {
444 for(unsigned int ii=0; ii<theOtherPhotons->size(); ii++)
445 {
446 thePhotons->push_back(theOtherPhotons->operator[](ii));
447 }
448 delete theOtherPhotons;
449 }
450 theGammaEnergy -= theGammas.GetLevelEnergy(iLevel);
451 if(iLevel == -1) break;
452 }
453 }
454
455// fill the result
456 unsigned int nSecondaries = tmpHadrons->size();
457 unsigned int nPhotons = 0;
458 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
459 nSecondaries += nPhotons;
460 G4DynamicParticle * theSec;
461
462 for(i=0; i<nSecondaries-nPhotons; i++)
463 {
464 theSec = new G4DynamicParticle;
465 theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
466 theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
467 theResult.AddSecondary(theSec);
468 delete tmpHadrons->operator[](i);
469 }
470 if(thePhotons != 0)
471 {
472 for(i=0; i<nPhotons; i++)
473 {
474 theSec = new G4DynamicParticle;
475 theSec->SetDefinition(thePhotons->operator[](i)->GetDefinition());
476 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
477 theResult.AddSecondary(theSec);
478 delete thePhotons->operator[](i);
479 }
480 }
481
482// some garbage collection
483 delete thePhotons;
484 delete tmpHadrons;
485
486//080721
487 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
488 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
489 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
490 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
491 adjust_final_state ( init_4p_lab );
492
493// clean up the primary neutron
494 theResult.SetStatusChange(stopAndKill);
495}
Note: See TracBrowser for help on using the repository browser.