ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPCaptureFS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPCaptureFS.cc
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 // 110430 Temporary solution in the case of being MF6 final state in Capture reaction (MT102)
35 //
36 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
37 //
38 #include "G4ParticleHPCaptureFS.hh"
39 #include "G4ParticleHPManager.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4Gamma.hh"
43 #include "G4ReactionProduct.hh"
44 #include "G4Nucleus.hh"
45 #include "G4PhotonEvaporation.hh"
46 #include "G4Fragment.hh"
47 #include "G4IonTable.hh"
48 #include "G4ParticleHPDataUsed.hh"
49 
51  {
52 
53  if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
54  theResult.Get()->Clear();
55 
56  G4int i;
57 
58 // prepare neutron
59  G4double eKinetic = theTrack.GetKineticEnergy();
60  const G4HadProjectile *incidentParticle = &theTrack;
61  G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition() ) );
62  theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
63  theNeutron.SetKineticEnergy( eKinetic );
64 
65  // Prepare target
67  G4Nucleus aNucleus;
68  G4double eps = 0.0001;
69  if (targetMass < 500*MeV) targetMass =
70  (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps) )) /
72  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
73  G4double temperature = theTrack.GetMaterial()->GetTemperature();
74  theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
76 
77  // Put neutron in nucleus rest system
78  theNeutron.Lorentz(theNeutron, theTarget);
79  eKinetic = theNeutron.GetKineticEnergy();
80 
81  // Sample the photons
82  G4ReactionProductVector * thePhotons = 0;
83  if ( HasFSData() && !G4ParticleHPManager::GetInstance()->GetUseOnlyPhotoEvaporation() )
84  {
85  //NDL has final state data
86  if ( hasExactMF6 ) {
87  theMF6FinalState.SetTarget(theTarget);
89  thePhotons = theMF6FinalState.Sample( eKinetic );
90  } else {
91  thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
92  }
93  if ( thePhotons == NULL ) {
94  throw G4HadronicException(__FILE__, __LINE__, "Final state data for photon is not properly allocated");
95  }
96  }
97  else
98  {
99  //NDL does not have final state data or forced to use PhotoEvaporation model
100  G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum();
101  G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
102  G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4);
103  G4PhotonEvaporation photonEvaporation;
104  // T. K. add
105  photonEvaporation.SetICM( TRUE );
106  G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
107  G4FragmentVector::iterator it;
108  thePhotons = new G4ReactionProductVector;
109  for(it=products->begin(); it!=products->end(); it++)
110  {
111  G4ReactionProduct * theOne = new G4ReactionProduct;
112  // T. K. add
113  if ( (*it)->GetParticleDefinition() != 0 )
114  theOne->SetDefinition( (*it)->GetParticleDefinition() );
115  else
116  theOne->SetDefinition( G4Gamma::Gamma() ); // this definiion will be over writen
117 
118  // T. K. comment out below line
119  //theOne->SetDefinition( G4Gamma::Gamma() );
120  G4IonTable* theTable = G4IonTable::GetIonTable();
121  if ( (*it)->GetMomentum().mag() > 10*MeV)
122  theOne->SetDefinition(theTable->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0) );
123 
124  if ( (*it)->GetExcitationEnergy() > 1.0e-2*eV) {
125  G4double ex = (*it)->GetExcitationEnergy();
126  G4ReactionProduct* aPhoton = new G4ReactionProduct;
127  aPhoton->SetDefinition( G4Gamma::Gamma() );
128  aPhoton->SetMomentum( (*it)->GetMomentum().vect().unit() * ex );
129  //aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum
130  thePhotons->push_back(aPhoton);
131  }
132 
133  theOne->SetMomentum( (*it)->GetMomentum().vect() * ( (*it)->GetMomentum().t() - (*it)->GetExcitationEnergy() ) / (*it)->GetMomentum().t() ) ;
134  thePhotons->push_back(theOne);
135  delete *it;
136  }
137  delete products;
138  }
139 
140  // Add them to the final state
141  G4int nPhotons = 0;
142  nPhotons=thePhotons->size();
143 
145  if ( DoNotAdjustFinalState() ) {
146 //Make at least one photon
147 //101203 TK
148  if ( nPhotons == 0 )
149  {
150  G4ReactionProduct* theOne = new G4ReactionProduct;
151  theOne->SetDefinition( G4Gamma::Gamma() );
152  // Bug #1745 DHW G4double theta = pi*G4UniformRand();
153  G4double costheta = 2.*G4UniformRand()-1.;
154  G4double theta = std::acos(costheta);
156  G4double sinth = std::sin(theta);
157  G4ThreeVector direction(sinth*std::cos(phi), sinth*std::sin(phi), costheta);
158  theOne->SetMomentum(direction);
159  thePhotons->push_back(theOne);
160  nPhotons++; // 0 -> 1
161  }
162 //One photon case: energy set to Q-value
163 //101203 TK
164  //if ( nPhotons == 1 )
165  if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
166  {
167  G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
168 
169  G4double Q = G4IonTable::GetIonTable()->GetIonMass(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0) + G4Neutron::Neutron()->GetPDGMass()
170  - G4IonTable::GetIonTable()->GetIonMass(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
171 
172  thePhotons->operator[](0)->SetMomentum( Q*direction );
173  }
174 //
175  }
176 
177  // back to lab system
178  for(i=0; i<nPhotons; i++)
179  {
180  thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1*theTarget);
181  }
182 
183  // Recoil, if only one gamma
184  //if (1==nPhotons)
185  if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
186  {
187  G4DynamicParticle * theOne = new G4DynamicParticle;
189  ->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
190  theOne->SetDefinition(aRecoil);
191  // Now energy;
192  // Can be done slightly better @
193  G4ThreeVector aMomentum = theTrack.Get4Momentum().vect()
194  +theTarget.GetMomentum()
195  -thePhotons->operator[](0)->GetMomentum();
196 
197  //TKDB 140520
198  //G4ThreeVector theMomUnit = aMomentum.unit();
199  //G4double aKinEnergy = theTrack.GetKineticEnergy()
200  // +theTarget.GetKineticEnergy(); // gammas come from Q-value
201  //G4double theResMass = aRecoil->GetPDGMass();
202  //G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
203  //G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
204  //G4ThreeVector theMomentum = theAbsMom*theMomUnit;
205  //theOne->SetMomentum(theMomentum);
206 
207  theOne->SetMomentum(aMomentum);
208  theResult.Get()->AddSecondary(theOne);
209  }
210 
211  // Now fill in the gammas.
212  for(i=0; i<nPhotons; i++)
213  {
214  // back to lab system
215  G4DynamicParticle * theOne = new G4DynamicParticle;
216  theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
217  theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
218  theResult.Get()->AddSecondary(theOne);
219  delete thePhotons->operator[](i);
220  }
221  delete thePhotons;
222 
223 //101203TK
224  G4bool residual = false;
226  ->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
227  for ( G4int j = 0 ; j != theResult.Get()->GetNumberOfSecondaries() ; j++ )
228  {
229  if ( theResult.Get()->GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
230  }
231 
232  if ( residual == false )
233  {
234  G4int nNonZero = 0;
235  G4LorentzVector p_photons(0,0,0,0);
236  for ( G4int j = 0 ; j != theResult.Get()->GetNumberOfSecondaries() ; j++ )
237  {
238  p_photons += theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum();
239  // To many 0 momentum photons -> Check PhotonDist
240  if ( theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0 ) nNonZero++;
241  }
242 
243  // Can we include kinetic energy here?
244  G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
245  - ( p_photons.e() + aRecoil->GetPDGMass() );
246 
247 //Add photons
248  if ( nPhotons - nNonZero > 0 )
249  {
250  //G4cout << "TKDB G4ParticleHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl;
251  std::vector<G4double> vRand;
252  vRand.push_back( 0.0 );
253  for ( G4int j = 0 ; j != nPhotons - nNonZero - 1 ; j++ )
254  {
255  vRand.push_back( G4UniformRand() );
256  }
257  vRand.push_back( 1.0 );
258  std::sort( vRand.begin(), vRand.end() );
259 
260  std::vector<G4double> vEPhoton;
261  for ( G4int j = 0 ; j < (G4int)vRand.size() - 1 ; j++ )
262  {
263  vEPhoton.push_back( deltaE * ( vRand[j+1] - vRand[j] ) );
264  }
265  std::sort( vEPhoton.begin(), vEPhoton.end() );
266 
267  for ( G4int j = 0 ; j < nPhotons - nNonZero - 1 ; j++ )
268  {
269  //Isotopic in LAB OK?
270  // Bug # 1745 DHW G4double theta = pi*G4UniformRand();
271  G4double costheta = 2.*G4UniformRand()-1.;
272  G4double theta = std::acos(costheta);
274  G4double sinth = std::sin(theta);
275  G4double en = vEPhoton[j];
276  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costheta);
277 
278  p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
279  G4DynamicParticle * theOne = new G4DynamicParticle;
280  theOne->SetDefinition( G4Gamma::Gamma() );
281  theOne->SetMomentum( tempVector );
282  theResult.Get()->AddSecondary(theOne);
283  }
284 
285 // Add last photon
286  G4DynamicParticle * theOne = new G4DynamicParticle;
287  theOne->SetDefinition( G4Gamma::Gamma() );
288 // For better momentum conservation
289  G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
290  p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
291  theOne->SetMomentum( lastPhoton );
292  theResult.Get()->AddSecondary(theOne);
293  }
294 
295 //Add residual
296  G4DynamicParticle * theOne = new G4DynamicParticle;
297  G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
298  - p_photons.vect();
299  theOne->SetDefinition(aRecoil);
300  theOne->SetMomentum( aMomentum );
301  theResult.Get()->AddSecondary(theOne);
302 
303  }
304 //101203TK END
305 
306 // clean up the primary neutron
308  return theResult.Get();
309  }
310 
311 #include <sstream>
313  {
314 
315  //TK110430 BEGIN
316  std::stringstream ss;
317  ss << static_cast<G4int>(Z);
318  G4String sZ;
319  ss >> sZ;
320  ss.clear();
321  ss << static_cast<G4int>(A);
322  G4String sA;
323  ss >> sA;
324 
325  ss.clear();
326  G4String sM;
327  if ( M > 0 )
328  {
329  ss << "m";
330  ss << M;
331  ss >> sM;
332  ss.clear();
333  }
334 
335  G4String element_name = theNames.GetName( static_cast<G4int>(Z)-1 );
336  G4String filenameMF6 = dirName+"/FSMF6/"+sZ+"_"+sA+sM+"_"+element_name;
337  //std::ifstream dummyIFS(filenameMF6, std::ios::in);
338  //if ( dummyIFS.good() == true ) hasExactMF6=true;
339  std::istringstream theData(std::ios::in);
340  G4ParticleHPManager::GetInstance()->GetDataStream(filenameMF6,theData);
341 
342  //TK110430 Only use MF6MT102 which has exactly same A and Z
343  //Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0
344  if ( theData.good() == true ) {
345  hasExactMF6=true;
346  theMF6FinalState.Init(theData);
347  //theData.close();
348  return;
349  }
350  //TK110430 END
351 
352 
353  G4String tString = "/FS";
354  G4bool dbool;
355  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
356 
357  G4String filename = aFile.GetName();
358  SetAZMs( A, Z, M, aFile );
359  //theBaseA = A;
360  //theBaseZ = G4int(Z+.5);
361  if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
362  {
363  hasAnyData = false;
364  hasFSData = false;
365  hasXsec = false;
366  return;
367  }
368  //std::ifstream theData(filename, std::ios::in);
369  //std::istringstream theData(std::ios::in);
370  theData.clear();
371  G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
373  if(hasFSData)
374  {
378  }
379  //theData.close();
380  }