LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Public Attributes | Private Attributes | List of all members
ArParticleHPCaptureFS Class Reference

#include "ArParticleHPCaptureFS.hh"

Inheritance diagram for ArParticleHPCaptureFS:

Public Member Functions

 ArParticleHPCaptureFS ()
 
void Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *) override
 
G4HadFinalState * ApplyYourself (const G4HadProjectile &theTrack) override
 
G4ParticleHPFinalState * New () override
 

Public Attributes

bool useArCapGamma = true
 

Private Attributes

G4double targetMass {0}
 
G4ParticleHPPhotonDist theFinalStatePhotons
 
ArCaptureGammas theFinalgammas
 
G4ParticleHPEnAngCorrelation theMF6FinalState
 
G4bool hasExactMF6 {false}
 
G4ParticleHPNames theNames
 

Detailed Description

Definition at line 45 of file ArParticleHPCaptureFS.hh.

Constructor & Destructor Documentation

ArParticleHPCaptureFS::ArParticleHPCaptureFS ( )
inline

Definition at line 49 of file ArParticleHPCaptureFS.hh.

References ApplyYourself(), Init(), and Z.

Referenced by New().

50  {
51  hasXsec = false;
52  }

Member Function Documentation

G4HadFinalState * ArParticleHPCaptureFS::ApplyYourself ( const G4HadProjectile &  theTrack)
override

Definition at line 53 of file ArParticleHPCaptureFS.cc.

References ArCaptureGammas::GetGammas(), hasExactMF6, targetMass, theFinalgammas, theFinalStatePhotons, theMF6FinalState, theTarget, and useArCapGamma.

Referenced by ArParticleHPCaptureFS().

54 {
55 
56  if (theResult.Get() == NULL)
57  theResult.Put(new G4HadFinalState);
58  theResult.Get()->Clear();
59 
60  G4int i;
61 
62  // prepare neutron
63  G4double eKinetic = theTrack.GetKineticEnergy();
64  const G4HadProjectile* incidentParticle = &theTrack;
65  G4ReactionProduct theNeutron(
66  const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition()));
67  theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect());
68  theNeutron.SetKineticEnergy(eKinetic);
69 
70  // Prepare target
71  G4ReactionProduct theTarget;
72  G4Nucleus aNucleus;
73  G4double eps = 0.0001;
74  if (targetMass < 500 * MeV)
75  targetMass = (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA + eps),
76  static_cast<G4int>(theBaseZ + eps))) /
77  G4Neutron::Neutron()->GetPDGMass();
78  G4ThreeVector neutronVelocity =
79  1. / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
80  G4double temperature = theTrack.GetMaterial()->GetTemperature();
81  theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
82  theTarget.SetDefinitionAndUpdateE(
83  G4IonTable::GetIonTable()->GetIon(G4int(theBaseZ), G4int(theBaseA), 0.0));
84 
85  // Put neutron in nucleus rest system
86  theNeutron.Lorentz(theNeutron, theTarget);
87  eKinetic = theNeutron.GetKineticEnergy();
88 
89  // Sample the photons
90  G4ReactionProductVector* thePhotons = 0;
91  if (useArCapGamma) {
92  thePhotons = theFinalgammas.GetGammas();
93  } else {
94  if (HasFSData() && !G4ParticleHPManager::GetInstance()->GetUseOnlyPhotoEvaporation()) {
95  // NDL has final state data
96  if (hasExactMF6) {
97  theMF6FinalState.SetTarget(theTarget);
98  theMF6FinalState.SetProjectileRP(theNeutron);
99  thePhotons = theMF6FinalState.Sample(eKinetic);
100  } else {
101  thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
102  }
103  if (thePhotons == NULL) {
104  throw G4HadronicException(
105  __FILE__, __LINE__, "Final state data for photon is not properly allocated");
106  }
107  } else {
108  // NDL does not have final state data or forced to use PhotoEvaporation model
109  G4ThreeVector aCMSMomentum = theNeutron.GetMomentum() + theTarget.GetMomentum();
110  G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
111  G4Fragment nucleus(static_cast<G4int>(theBaseA + 1), static_cast<G4int>(theBaseZ), p4);
112  G4PhotonEvaporation photonEvaporation;
113  // T. K. add
114  photonEvaporation.SetICM(TRUE);
115  G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
116  G4FragmentVector::iterator it;
117  thePhotons = new G4ReactionProductVector;
118  for (it = products->begin(); it != products->end(); it++) {
119  G4ReactionProduct* theOne = new G4ReactionProduct;
120  // T. K. add
121  if ((*it)->GetParticleDefinition() != 0)
122  theOne->SetDefinition((*it)->GetParticleDefinition());
123  else
124  theOne->SetDefinition(G4Gamma::Gamma()); // this definiion will be over writen
125 
126  // T. K. comment out below line
127  // theOne->SetDefinition( G4Gamma::Gamma() );
128  G4IonTable* theTable = G4IonTable::GetIonTable();
129  if ((*it)->GetMomentum().mag() > 10 * MeV)
130  theOne->SetDefinition(
131  theTable->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0));
132 
133  if ((*it)->GetExcitationEnergy() > 1.0e-2 * eV) {
134  G4double ex = (*it)->GetExcitationEnergy();
135  G4ReactionProduct* aPhoton = new G4ReactionProduct;
136  aPhoton->SetDefinition(G4Gamma::Gamma());
137  aPhoton->SetMomentum((*it)->GetMomentum().vect().unit() * ex);
138  // aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum
139  thePhotons->push_back(aPhoton);
140  }
141 
142  theOne->SetMomentum((*it)->GetMomentum().vect() *
143  ((*it)->GetMomentum().t() - (*it)->GetExcitationEnergy()) /
144  (*it)->GetMomentum().t());
145  thePhotons->push_back(theOne);
146  delete *it;
147  }
148  delete products;
149  }
150  }
151 
152  // Add them to the final state
153  G4int nPhotons = 0;
154  nPhotons = thePhotons->size();
155 
157  if (DoNotAdjustFinalState()) {
158  // Make at least one photon
159  // 101203 TK
160  if (nPhotons == 0) {
161  G4ReactionProduct* theOne = new G4ReactionProduct;
162  theOne->SetDefinition(G4Gamma::Gamma());
163  // Bug #1745 DHW G4double theta = pi*G4UniformRand();
164  G4double costheta = 2. * G4UniformRand() - 1.;
165  G4double theta = std::acos(costheta);
166  G4double phi = twopi * G4UniformRand();
167  G4double sinth = std::sin(theta);
168  G4ThreeVector direction(sinth * std::cos(phi), sinth * std::sin(phi), costheta);
169  theOne->SetMomentum(direction);
170  thePhotons->push_back(theOne);
171  nPhotons++; // 0 -> 1
172  }
173  // One photon case: energy set to Q-value
174  // 101203 TK
175  // if ( nPhotons == 1 )
176  if (nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0) {
177  G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
178 
179  G4double Q = G4IonTable::GetIonTable()->GetIonMass(
180  static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0) +
181  G4Neutron::Neutron()->GetPDGMass() -
182  G4IonTable::GetIonTable()->GetIonMass(
183  static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);
184 
185  thePhotons->operator[](0)->SetMomentum(Q * direction);
186  }
187  //
188  }
189 
190  // back to lab system
191  for (i = 0; i < nPhotons; i++) {
192  thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1 * theTarget);
193  }
194 
195  // Recoil, if only one gamma
196  // if (1==nPhotons)
197  if (nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0) {
198  G4DynamicParticle* theOne = new G4DynamicParticle;
199  G4ParticleDefinition* aRecoil = G4IonTable::GetIonTable()->GetIon(
200  static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);
201  theOne->SetDefinition(aRecoil);
202  // Now energy;
203  // Can be done slightly better @
204  G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum() -
205  thePhotons->operator[](0)->GetMomentum();
206 
207  // TKDB 140520
208  // G4ThreeVector theMomUnit = aMomentum.unit();
209  // G4double aKinEnergy = theTrack.GetKineticEnergy()
210  // +theTarget.GetKineticEnergy(); // gammas come from Q-value
211  // G4double theResMass = aRecoil->GetPDGMass();
212  // G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
213  // G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
214  // G4ThreeVector theMomentum = theAbsMom*theMomUnit;
215  // theOne->SetMomentum(theMomentum);
216 
217  theOne->SetMomentum(aMomentum);
218  theResult.Get()->AddSecondary(theOne);
219  }
220 
221  // Now fill in the gammas.
222  for (i = 0; i < nPhotons; i++) {
223  // back to lab system
224  G4DynamicParticle* theOne = new G4DynamicParticle;
225  theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
226  theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
227  theResult.Get()->AddSecondary(theOne);
228  delete thePhotons->operator[](i);
229  }
230  delete thePhotons;
231 
232  // 101203TK
233  G4bool residual = false;
234  G4ParticleDefinition* aRecoil = G4IonTable::GetIonTable()->GetIon(
235  static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);
236  for (G4int j = 0; j != theResult.Get()->GetNumberOfSecondaries(); j++) {
237  if (theResult.Get()->GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil)
238  residual = true;
239  }
240 
241  if (residual == false) {
242  G4int nNonZero = 0;
243  G4LorentzVector p_photons(0, 0, 0, 0);
244  for (G4int j = 0; j != theResult.Get()->GetNumberOfSecondaries(); j++) {
245  p_photons += theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum();
246  // To many 0 momentum photons -> Check PhotonDist
247  if (theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0)
248  nNonZero++;
249  }
250 
251  // Can we include kinetic energy here?
252  G4double deltaE = (theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy()) -
253  (p_photons.e() + aRecoil->GetPDGMass());
254 
255  // Add photons
256  if (nPhotons - nNonZero > 0) {
257  // G4cout << "TKDB ArParticleHPCaptureFS::ApplyYourself we will create additional " <<
258  // nPhotons - nNonZero << " photons" << G4endl;
259  std::vector<G4double> vRand;
260  vRand.push_back(0.0);
261  for (G4int j = 0; j != nPhotons - nNonZero - 1; j++) {
262  vRand.push_back(G4UniformRand());
263  }
264  vRand.push_back(1.0);
265  std::sort(vRand.begin(), vRand.end());
266 
267  std::vector<G4double> vEPhoton;
268  for (G4int j = 0; j < (G4int)vRand.size() - 1; j++) {
269  vEPhoton.push_back(deltaE * (vRand[j + 1] - vRand[j]));
270  }
271  std::sort(vEPhoton.begin(), vEPhoton.end());
272 
273  for (G4int j = 0; j < nPhotons - nNonZero - 1; j++) {
274  // Isotopic in LAB OK?
275  // Bug # 1745 DHW G4double theta = pi*G4UniformRand();
276  G4double costheta = 2. * G4UniformRand() - 1.;
277  G4double theta = std::acos(costheta);
278  G4double phi = twopi * G4UniformRand();
279  G4double sinth = std::sin(theta);
280  G4double en = vEPhoton[j];
281  G4ThreeVector tempVector(
282  en * sinth * std::cos(phi), en * sinth * std::sin(phi), en * costheta);
283 
284  p_photons += G4LorentzVector(tempVector, tempVector.mag());
285  G4DynamicParticle* theOne = new G4DynamicParticle;
286  theOne->SetDefinition(G4Gamma::Gamma());
287  theOne->SetMomentum(tempVector);
288  theResult.Get()->AddSecondary(theOne);
289  }
290 
291  // Add last photon
292  G4DynamicParticle* theOne = new G4DynamicParticle;
293  theOne->SetDefinition(G4Gamma::Gamma());
294  // For better momentum conservation
295  G4ThreeVector lastPhoton = -p_photons.vect().unit() * vEPhoton.back();
296  p_photons += G4LorentzVector(lastPhoton, lastPhoton.mag());
297  theOne->SetMomentum(lastPhoton);
298  theResult.Get()->AddSecondary(theOne);
299  }
300 
301  // Add residual
302  G4DynamicParticle* theOne = new G4DynamicParticle;
303  G4ThreeVector aMomentum =
304  theTrack.Get4Momentum().vect() + theTarget.GetMomentum() - p_photons.vect();
305  theOne->SetDefinition(aRecoil);
306  theOne->SetMomentum(aMomentum);
307  theResult.Get()->AddSecondary(theOne);
308  }
309  // 101203TK END
310 
311  // clean up the primary neutron
312  theResult.Get()->SetStatusChange(stopAndKill);
313 
314  //<--for(int i=0;i<100;i++){
315  //<-- std::cout<<"I am ArParticleHPCaptureFS"<<std::endl;
316  //<--}
317 
318  return theResult.Get();
319 }
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
ArCaptureGammas::GetGammas
G4ReactionProductVector * GetGammas()
Definition: ArCaptureGammas.cc:30
theTarget
G4ErrorTarget * theTarget
Definition: errprop.cc:58
ArParticleHPCaptureFS::theMF6FinalState
G4ParticleHPEnAngCorrelation theMF6FinalState
Definition: ArParticleHPCaptureFS.hh:73
ArParticleHPCaptureFS::theFinalStatePhotons
G4ParticleHPPhotonDist theFinalStatePhotons
Definition: ArParticleHPCaptureFS.hh:70
void ArParticleHPCaptureFS::Init ( G4double  A,
G4double  Z,
G4int  M,
G4String &  dirName,
G4String &  aFSType,
G4ParticleDefinition *   
)
override

Definition at line 323 of file ArParticleHPCaptureFS.cc.

References util::abs(), hasExactMF6, in, targetMass, theFinalStatePhotons, theMF6FinalState, theNames, and Z.

Referenced by ArParticleHPCaptureFS().

329 {
330 
331  // TK110430 BEGIN
332  std::stringstream ss;
333  ss << static_cast<G4int>(Z);
334  G4String sZ;
335  ss >> sZ;
336  ss.clear();
337  ss << static_cast<G4int>(A);
338  G4String sA;
339  ss >> sA;
340 
341  ss.clear();
342  G4String sM;
343  if (M > 0) {
344  ss << "m";
345  ss << M;
346  ss >> sM;
347  ss.clear();
348  }
349 
350  G4String element_name = theNames.GetName(static_cast<G4int>(Z) - 1);
351  G4String filenameMF6 = dirName + "/FSMF6/" + sZ + "_" + sA + sM + "_" + element_name;
352  // std::ifstream dummyIFS(filenameMF6, std::ios::in);
353  // if ( dummyIFS.good() == true ) hasExactMF6=true;
354  std::istringstream theData(std::ios::in);
355  G4ParticleHPManager::GetInstance()->GetDataStream(filenameMF6, theData);
356 
357  // TK110430 Only use MF6MT102 which has exactly same A and Z
358  // Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0
359  if (theData.good() == true) {
360  hasExactMF6 = true;
361  theMF6FinalState.Init(theData);
362  // theData.close();
363  return;
364  }
365  // TK110430 END
366 
367  G4String tString = "/FS";
368  G4bool dbool;
369  G4ParticleHPDataUsed aFile =
370  theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
371 
372  G4String filename = aFile.GetName();
373  SetAZMs(A, Z, M, aFile);
374  // theBaseA = A;
375  // theBaseZ = G4int(Z+.5);
376  if (!dbool || (Z < 2.5 && (std::abs(theBaseZ - Z) > 0.0001 || std::abs(theBaseA - A) > 0.0001))) {
377  hasAnyData = false;
378  hasFSData = false;
379  hasXsec = false;
380  return;
381  }
382  // std::ifstream theData(filename, std::ios::in);
383  // std::istringstream theData(std::ios::in);
384  theData.clear();
385  G4ParticleHPManager::GetInstance()->GetDataStream(filename, theData);
386  hasFSData = theFinalStatePhotons.InitMean(theData);
387  if (hasFSData) {
388  targetMass = theFinalStatePhotons.GetTargetMass();
389  theFinalStatePhotons.InitAngular(theData);
390  theFinalStatePhotons.InitEnergies(theData);
391  }
392  // theData.close();
393 }
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
Z
Float_t Z
Definition: plot.C:37
ArParticleHPCaptureFS::theMF6FinalState
G4ParticleHPEnAngCorrelation theMF6FinalState
Definition: ArParticleHPCaptureFS.hh:73
in
ifstream in
Definition: comparison.C:7
ArParticleHPCaptureFS::theFinalStatePhotons
G4ParticleHPPhotonDist theFinalStatePhotons
Definition: ArParticleHPCaptureFS.hh:70
G4ParticleHPFinalState* ArParticleHPCaptureFS::New ( )
inlineoverride

Definition at line 62 of file ArParticleHPCaptureFS.hh.

References ArParticleHPCaptureFS().

63  {
64  return new ArParticleHPCaptureFS;
65  }

Member Data Documentation

G4bool ArParticleHPCaptureFS::hasExactMF6 {false}
private

Definition at line 74 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself(), and Init().

G4double ArParticleHPCaptureFS::targetMass {0}
private

Definition at line 68 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself(), and Init().

ArCaptureGammas ArParticleHPCaptureFS::theFinalgammas
private

Definition at line 71 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself().

G4ParticleHPPhotonDist ArParticleHPCaptureFS::theFinalStatePhotons
private

Definition at line 70 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself(), and Init().

G4ParticleHPEnAngCorrelation ArParticleHPCaptureFS::theMF6FinalState
private

Definition at line 73 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself(), and Init().

G4ParticleHPNames ArParticleHPCaptureFS::theNames
private

Definition at line 76 of file ArParticleHPCaptureFS.hh.

Referenced by Init().

bool ArParticleHPCaptureFS::useArCapGamma = true

Definition at line 47 of file ArParticleHPCaptureFS.hh.

Referenced by ApplyYourself().

The documentation for this class was generated from the following files: