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

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

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