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

Last change on this file since 1340 was 962, checked in by garnier, 15 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.