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

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

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