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

Last change on this file since 1353 was 1347, checked in by garnier, 14 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.