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

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

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