source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPCaptureFS.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 12.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// 12-April-06 Enable IC electron emissions T. Koi
31// 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION flag
32// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33// 101203 Bugzilla/Geant4 Problem 1155 Lack of residual in some case
34//
35#include "G4NeutronHPCaptureFS.hh"
36#include "G4Gamma.hh"
37#include "G4ReactionProduct.hh"
38#include "G4Nucleus.hh"
39#include "G4PhotonEvaporation.hh"
40#include "G4Fragment.hh"
41#include "G4ParticleTable.hh"
42#include "G4NeutronHPDataUsed.hh"
43
44  G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself(const G4HadProjectile & theTrack)
45  {
46
47    G4int i;
48    theResult.Clear();
49// prepare neutron
50    G4double eKinetic = theTrack.GetKineticEnergy();
51    const G4HadProjectile *incidentParticle = &theTrack;
52    G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
53    theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
54    theNeutron.SetKineticEnergy( eKinetic );
55
56// prepare target
57    G4ReactionProduct theTarget; 
58    G4Nucleus aNucleus;
59    G4double eps = 0.0001;
60    if(targetMass<500*MeV)
61      targetMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theBaseA+eps) , static_cast<G4int>(theBaseZ+eps) )) /
62                     G4Neutron::Neutron()->GetPDGMass();
63    G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
64    G4double temperature = theTrack.GetMaterial()->GetTemperature();
65    theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
66
67// go to nucleus rest system
68    theNeutron.Lorentz(theNeutron, -1*theTarget);
69    eKinetic = theNeutron.GetKineticEnergy();
70
71// dice the photons
72
73    G4ReactionProductVector * thePhotons = 0;
74    if ( HasFSData() && !getenv ( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) ) 
75    { 
76      thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
77    }
78    else
79    {
80      G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum();
81      G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
82      G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4);
83      G4PhotonEvaporation photonEvaporation;
84      // T. K. add
85      photonEvaporation.SetICM( TRUE );
86      G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
87      G4FragmentVector::iterator i;
88      thePhotons = new G4ReactionProductVector;
89      for(i=products->begin(); i!=products->end(); i++)
90      {
91        G4ReactionProduct * theOne = new G4ReactionProduct;
92        // T. K. add
93        if ( (*i)->GetParticleDefinition() != 0 ) 
94           theOne->SetDefinition( (*i)->GetParticleDefinition() );
95        else
96           theOne->SetDefinition( G4Gamma::Gamma() ); // this definiion will be over writen
97       
98        // T. K. comment out below line
99        //theOne->SetDefinition( G4Gamma::Gamma() );
100        G4ParticleTable* theTable = G4ParticleTable::GetParticleTable();
101        if((*i)->GetMomentum().mag() > 10*MeV) 
102                 theOne->SetDefinition( 
103                 theTable->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)) );
104        theOne->SetMomentum( (*i)->GetMomentum().vect() ) ;
105        theOne->SetTotalEnergy( (*i)->GetMomentum().t() );
106        thePhotons->push_back(theOne);
107        delete *i;
108      } 
109      delete products;
110    }
111
112
113
114// add them to the final state
115
116    G4int nPhotons = 0;
117    if(thePhotons!=0) nPhotons=thePhotons->size();
118    G4int nParticles = nPhotons;
119    if(1==nPhotons) nParticles = 2;
120
121
122//Make at least one photon 
123//101203 TK
124    if ( nPhotons == 0 )
125    {
126       G4ReactionProduct * theOne = new G4ReactionProduct;
127       theOne->SetDefinition( G4Gamma::Gamma() ); 
128       G4double theta = pi*G4UniformRand();
129       G4double phi = twopi*G4UniformRand();
130       G4double sinth = std::sin(theta);
131       G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) );
132       theOne->SetMomentum( direction ) ;
133       thePhotons->push_back(theOne);
134       nPhotons++; // 0 -> 1
135    }
136//One photon case: energy set to Q-value
137//101203 TK
138    if ( nPhotons == 1 )
139    {
140       G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
141       G4double Q = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass()
142         - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass();
143       thePhotons->operator[](0)->SetMomentum( Q*direction );
144    } 
145//
146
147    // back to lab system
148    for(i=0; i<nPhotons; i++)
149    {
150      thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), theTarget);
151    }
152   
153    // Recoil, if only one gamma
154    if (1==nPhotons)
155    {
156       G4DynamicParticle * theOne = new G4DynamicParticle;
157       G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
158                                        ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
159       theOne->SetDefinition(aRecoil);
160       // Now energy;
161       // Can be done slightly better @
162       G4ThreeVector aMomentum =  theTrack.Get4Momentum().vect()
163                                 +theTarget.GetMomentum()
164                                 -thePhotons->operator[](0)->GetMomentum();
165
166       G4ThreeVector theMomUnit = aMomentum.unit();
167       G4double aKinEnergy =  theTrack.GetKineticEnergy()
168                             +theTarget.GetKineticEnergy(); // gammas come from Q-value
169       G4double theResMass = aRecoil->GetPDGMass();
170       G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
171       G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
172       G4ThreeVector theMomentum = theAbsMom*theMomUnit;
173       theOne->SetMomentum(theMomentum);
174       theResult.AddSecondary(theOne);
175    }
176
177    // Now fill in the gammas.
178    for(i=0; i<nPhotons; i++)
179    {
180      // back to lab system
181      G4DynamicParticle * theOne = new G4DynamicParticle;
182      theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
183      theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
184      theResult.AddSecondary(theOne);
185      delete thePhotons->operator[](i);
186    }
187    delete thePhotons; 
188
189//101203TK
190    G4bool residual = false;
191    G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
192                                   ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
193    for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
194    {
195       if ( theResult.GetSecondary(i)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
196    }
197
198    if ( residual == false )
199    {
200       G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
201                                        ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
202       G4int nNonZero = 0;
203       G4LorentzVector p_photons(0,0,0,0);
204       for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
205       {
206          p_photons += theResult.GetSecondary(i)->GetParticle()->Get4Momentum();
207          // To many 0 momentum photons -> Check PhotonDist
208          if ( theResult.GetSecondary(i)->GetParticle()->Get4Momentum() > 0 ) nNonZero++;
209       }
210
211       // Can we include kinetic energy here?
212       G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
213                       - ( p_photons.e() + aRecoil->GetPDGMass() );
214
215//Add photons
216       if ( nPhotons - nNonZero > 0 ) 
217       {
218              //G4cout << "TKDB G4NeutronHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl;
219          std::vector<G4double> vRand;
220          vRand.push_back( 0.0 );
221          for ( G4int i = 0 ; i != nPhotons - nNonZero - 1 ; i++ )
222          { 
223             vRand.push_back( G4UniformRand() );
224          }
225          vRand.push_back( 1.0 );
226          std::sort( vRand.begin(), vRand.end() );
227
228          std::vector<G4double> vEPhoton;
229          for ( G4int i = 0 ; i < (G4int)vRand.size() - 1 ; i++ )
230          {
231             vEPhoton.push_back( deltaE * ( vRand[i+1] - vRand[i] ) );
232          }
233          std::sort( vEPhoton.begin(), vEPhoton.end() );
234
235          for ( G4int i = 0 ; i < nPhotons - nNonZero - 1 ; i++ )
236          {
237             //Isotopic in LAB OK?
238             G4double theta = pi*G4UniformRand();
239             G4double phi = twopi*G4UniformRand();
240             G4double sinth = std::sin(theta);
241             G4double en = vEPhoton[i];
242             G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
243             
244             p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
245             G4DynamicParticle * theOne = new G4DynamicParticle;
246             theOne->SetDefinition( G4Gamma::Gamma() );
247             theOne->SetMomentum( tempVector );
248             theResult.AddSecondary(theOne);
249          }
250
251//        Add last photon
252          G4DynamicParticle * theOne = new G4DynamicParticle;
253          theOne->SetDefinition( G4Gamma::Gamma() );
254//        For better momentum conservation
255          G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
256          p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
257          theOne->SetMomentum( lastPhoton );
258          theResult.AddSecondary(theOne);
259       }
260
261//Add residual
262       G4DynamicParticle * theOne = new G4DynamicParticle;
263       G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
264                               - p_photons.vect();
265       theOne->SetDefinition(aRecoil);
266       theOne->SetMomentum( aMomentum );
267       theResult.AddSecondary(theOne);
268
269    }
270//101203TK END
271
272// clean up the primary neutron
273    theResult.SetStatusChange(stopAndKill);
274    return &theResult;
275  }
276
277  void G4NeutronHPCaptureFS::Init (G4double A, G4double Z, G4String & dirName, G4String & )
278  {
279    G4String tString = "/FS/";
280    G4bool dbool;
281    G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), dirName, tString, dbool);
282    G4String filename = aFile.GetName();
283    theBaseA = A;
284    theBaseZ = G4int(Z+.5);
285    if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
286    {
287      hasAnyData = false;
288      hasFSData = false; 
289      hasXsec = false;
290      return;
291    }
292    std::ifstream theData(filename, std::ios::in);
293   
294    hasFSData = theFinalStatePhotons.InitMean(theData); 
295    if(hasFSData)
296    {
297      targetMass = theFinalStatePhotons.GetTargetMass();
298      theFinalStatePhotons.InitAngular(theData); 
299      theFinalStatePhotons.InitEnergies(theData); 
300    }
301    theData.close();
302  }
Note: See TracBrowser for help on using the repository browser.