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

Last change on this file since 847 was 819, checked in by garnier, 16 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.