ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPInelasticCompFS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPInelasticCompFS.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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 070606 bug fix and migrate to enable to Partial cases by T. Koi
32 // 080603 bug fix for Hadron Hyper News #932 by T. Koi
33 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34 // 080717 bug fix of calculation of residual momentum by T. Koi
35 // 080801 protect negative avalable energy by T. Koi
36 // introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
41 // add two_body_reaction
42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation
44 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
45 //
46 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
47 //
49 #include "G4ParticleHPManager.hh"
50 #include "G4Nucleus.hh"
51 #include "G4NucleiProperties.hh"
52 #include "G4He3.hh"
53 #include "G4Alpha.hh"
54 #include "G4Electron.hh"
55 #include "G4ParticleHPDataUsed.hh"
56 #include "G4IonTable.hh"
57 #include "G4Pow.hh"
58 #include "G4SystemOfUnits.hh"
59 
60 #include "G4NRESP71M03.hh" // nresp71_m03.hh and nresp71_m02.hh are alike. The only difference between m02 and m03 is in the total carbon cross section that is properly included in the latter. These data are not used in nresp71_m0*.hh.
61 
62 // June-2019 - E. Mendoza - re-build "two_body_reaction", to be used by incident charged particles (now isotropic emission in the CMS). Also restrict nresp use below 20 MeV (for future developments). Add photon emission when no data available.
63 
65 {
66  // char the[100] = {""};
67  // std::ostrstream ost(the, 100, std::ios::out);
68  // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
69  // G4String * aName = new G4String(the);
70  // std::ifstream from(*aName, std::ios::in);
71 
72  std::ostringstream ost;
73  ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
74  G4String aName = ost.str();
75  std::ifstream from(aName, std::ios::in);
76 
77  if(!from) return; // no data found for this isotope
78  // std::ifstream theGammaData(*aName, std::ios::in);
79  std::ifstream theGammaData(aName, std::ios::in);
80 
81  theGammas.Init(theGammaData);
82  // delete aName;
83 
84 }
85 
87 {
88  gammaPath = "/Inelastic/Gammas/"; //only in neutron data base
89  if(!std::getenv("G4NEUTRONHPDATA"))
90  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
91  G4String tBase = std::getenv("G4NEUTRONHPDATA");
92  gammaPath = tBase+gammaPath;
93  G4String tString = dirName;
94  G4bool dbool;
95  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
96  G4String filename = aFile.GetName();
97 #ifdef G4PHPDEBUG
98  if( std::getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
99 #endif
100 
101  SetAZMs( A, Z, M, aFile );
102  //theBaseA = aFile.GetA();
103  //theBaseZ = aFile.GetZ();
104  //theNDLDataA = (int)aFile.GetA();
105  //theNDLDataZ = aFile.GetZ();
106  //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
107  if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
108  {
109 #ifdef G4PHPDEBUG
110  if(std::getenv("G4ParticleHPDebug_NamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
111 #endif
112  hasAnyData = false;
113  hasFSData = false;
114  hasXsec = false;
115  return;
116  }
117  // theBaseA = A;
118  // theBaseZ = G4int(Z+.5);
119 //std::ifstream theData(filename, std::ios::in);
120  std::istringstream theData(std::ios::in);
121  G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
122  if(!theData) //"!" is a operator of ios
123  {
124  hasAnyData = false;
125  hasFSData = false;
126  hasXsec = false;
127  // theData.close();
128  return;
129  }
130  // here we go
131  G4int infoType, dataType, dummy;
132  G4int sfType, it;
133  hasFSData = false;
134  while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
135  {
136  hasFSData = true;
137  theData >> dataType;
138  theData >> sfType >> dummy;
139  it = 50;
140  if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
141  if(dataType==3)
142  {
143  //theData >> dummy >> dummy;
144  //TK110430
145  // QI and LR introudced since G4NDL3.15
146  G4double dqi;
147  G4int ilr;
148  theData >> dqi >> ilr;
149 
150  QI[ it ] = dqi*CLHEP::eV;
151  LR[ it ] = ilr;
153  G4int total;
154  theData >> total;
155  theXsection[it]->Init(theData, total, CLHEP::eV);
156  //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
157  }
158  else if(dataType==4)
159  {
161  theAngularDistribution[it]->Init(theData);
162  }
163  else if(dataType==5)
164  {
166  theEnergyDistribution[it]->Init(theData);
167  }
168  else if(dataType==6)
169  {
171  // G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl; //GDEB
172  theEnergyAngData[it]->Init(theData);
173  }
174  else if(dataType==12)
175  {
177  theFinalStatePhotons[it]->InitMean(theData);
178  }
179  else if(dataType==13)
180  {
183  }
184  else if(dataType==14)
185  {
187  }
188  else if(dataType==15)
189  {
191  }
192  else
193  {
194  throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4ParticleHPInelasticCompFS");
195  }
196  }
197  // theData.close();
198 }
199 
201 {
202 
203 // it = 0 has without Photon
204  G4double running[50];
205  running[0] = 0;
206  unsigned int i;
207  for(i=0; i<50; i++)
208  {
209  if(i!=0) running[i]=running[i-1];
210  if(theXsection[i] != 0)
211  {
212  running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
213  }
214  }
215  G4double random = G4UniformRand();
216  G4double sum = running[49];
217  G4int it = 50;
218  if(0!=sum)
219  {
220  G4int i0;
221  for(i0=0; i0<50; i0++)
222  {
223  it = i0;
224  // G4cout << " SelectExitChannel " << it << " " << random << " " << running[i0]/sum << " " << running[i0] << G4endl; //GDEB
225  if(random < running[i0]/sum) break;
226  }
227  }
228 //debug: it = 1;
229 // G4cout << " SelectExitChannel " << it << " " << sum << G4endl; //GDEB
230  return it;
231 }
232 
233 
234 // n,p,d,t,he3,a
236  G4ParticleDefinition* aDefinition)
237 {
238 
239 // prepare neutron
240  if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
241  theResult.Get()->Clear();
242  G4double eKinetic = theTrack.GetKineticEnergy();
243  const G4HadProjectile *hadProjectile = &theTrack;
244  G4ReactionProduct incidReactionProduct( const_cast<G4ParticleDefinition *>(hadProjectile->GetDefinition()) ); // incidReactionProduct
245  incidReactionProduct.SetMomentum( hadProjectile->Get4Momentum().vect() );
246  incidReactionProduct.SetKineticEnergy( eKinetic );
247 
248 // prepare target
249  G4int i;
250  for(i=0; i<50; i++)
251  { if(theXsection[i] != 0) { break; } }
252 
253  G4double targetMass=0;
254  G4double eps = 0.0001;
255  targetMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps));
256 #ifdef G4PHPDEBUG
257  if( std::getenv("G4ParticleHPDebug")) G4cout <<this <<" G4ParticleHPInelasticCompFS::CompositeApply A " <<theBaseA <<" Z " <<theBaseZ <<" incident " <<hadProjectile->GetDefinition()->GetParticleName() <<G4endl;
258 #endif
259 // if(theEnergyAngData[i]!=0)
260 // targetMass = theEnergyAngData[i]->GetTargetMass();
261 // else if(theAngularDistribution[i]!=0)
262 // targetMass = theAngularDistribution[i]->GetTargetMass();
263 // else if(theFinalStatePhotons[50]!=0)
264 // targetMass = theFinalStatePhotons[50]->GetTargetMass();
266  G4Nucleus aNucleus;
267  //G4ThreeVector neuVelo = (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum();
268  //theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() , neuVelo, theTrack.GetMaterial()->GetTemperature());
269  //G4Nucleus::GetBiasedThermalNucleus requests normalization of mass and velocity in neutron mass
270  G4ThreeVector neuVelo = ( 1./G4Neutron::Neutron()->GetPDGMass() )*incidReactionProduct.GetMomentum();
271  theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/G4Neutron::Neutron()->GetPDGMass()
272  , neuVelo, theTrack.GetMaterial()->GetTemperature() );
273 
274  theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //XX
275 
276 // prepare the residual mass
277  G4double residualMass=0;
278  G4double residualZ = theBaseZ + theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge();
279  G4double residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
280  residualMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps));
281 
282 // prepare energy in target rest frame
283  G4ReactionProduct boosted;
284  boosted.Lorentz(incidReactionProduct, theTarget);
285  eKinetic = boosted.GetKineticEnergy();
286 // G4double momentumInCMS = boosted.GetTotalMomentum();
287 
288 // select exit channel for composite FS class.
289  G4int it = SelectExitChannel( eKinetic );
290 
291  //E. Mendoza (2018) -- to use JENDL/AN-2005
292  if(theEnergyDistribution[it]==0 && theAngularDistribution[it]==0 && theEnergyAngData[it]==0){
293  if(theEnergyDistribution[50]!=0 || theAngularDistribution[50]!=0 || theEnergyAngData[50]!=0){
294  it=50;
295  }
296  }
297 
298  // set target and neutron in the relevant exit channel
299  InitDistributionInitialState(incidReactionProduct, theTarget, it);
300 
301  //---------------------------------------------------------------------//
302  //Hook for NRESP71MODEL
303  if ( G4ParticleHPManager::GetInstance()->GetUseNRESP71Model() && eKinetic<20*MeV) {
304  if ( (G4int)(theBaseZ+0.1) == 6 ) // If the reaction is with Carbon...
305  {
307  if ( use_nresp71_model( aDefinition , it , theTarget , boosted ) ) return;
308  }
309  }
310  }
311  //---------------------------------------------------------------------//
312 
313  G4ReactionProductVector * thePhotons = 0;
314  G4ReactionProductVector * theParticles = 0;
315  G4ReactionProduct aHadron;
316  aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
317  G4double availableEnergy = incidReactionProduct.GetKineticEnergy() + incidReactionProduct.GetMass() - aHadron.GetMass() +
318  (targetMass - residualMass);
319 //080730c
320  if ( availableEnergy < 0 )
321  {
322  //G4cout << "080730c Adjust availavleEnergy " << G4endl;
323  availableEnergy = 0;
324  }
325  G4int nothingWasKnownOnHadron = 0;
326  G4int dummy;
327  G4double eGamm = 0;
328  G4int iLevel=it-1;
329 
330 // TK without photon has it = 0
331  if( 50 == it )
332  {
333 
334 // TK Excitation level is not determined
335  iLevel=-1;
336  aHadron.SetKineticEnergy(availableEnergy*residualMass/
337  (aHadron.GetMass()+residualMass));
338 
339  //aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*
340  // std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
341  // aHadron.GetMass()*aHadron.GetMass()));
342 
343  //TK add safty 100909
344  G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
345  G4double p = 0.0;
346  if ( p2 > 0.0 ) p = std::sqrt( p );
347 
348  aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*p );
349 
350  }
351  else
352  {
353  while ( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; } // Loop checking, 11.05.2015, T. Koi
354  }
355 
356 
357  if ( theAngularDistribution[it] != 0 ) // MF4
358  {
359  if(theEnergyDistribution[it]!=0) // MF5
360  {
361  //************************************************************
362  /*
363  aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
364  G4double eSecN = aHadron.GetKineticEnergy();
365  */
366  //************************************************************
367  //EMendoza --> maximum allowable energy should be taken into account.
368  G4double dqi = 0.0;
369  if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
370  G4double MaxEne=eKinetic+dqi;
371  G4double eSecN;
372 
373  G4int icounter=0;
374  G4int icounter_max=1024;
375  do {
376  icounter++;
377  if ( icounter > icounter_max ) {
378  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
379  break;
380  }
381  eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
382  }while(eSecN>MaxEne); // Loop checking, 11.05.2015, T. Koi
383  aHadron.SetKineticEnergy(eSecN);
384  //************************************************************
385  eGamm = eKinetic-eSecN;
386  for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
387  {
388  if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
389  }
390  G4double random = 2*G4UniformRand();
391  iLevel+=G4int(random);
392  if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
393  }
394  else
395  {
396  G4double eExcitation = 0;
397  if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
398  while (eKinetic-eExcitation < 0 && iLevel>0) // Loop checking, 11.05.2015, T. Koi
399  {
400  iLevel--;
401  eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
402  }
403  //110610TK BEGIN
404  //Use QI value for calculating excitation energy of residual.
405  G4bool useQI=false;
406  G4double dqi = QI[it];
407  if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range
408 
409  if (useQI) {
410  // QI introudced since G4NDL3.15
411  // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
412  // eExcitation = QM-QI[it];
413  eExcitation = QI[0] - QI[it]; // Bug fix #1838
414  if(eExcitation < 20*CLHEP::keV) eExcitation = 0;
415 
416  // Re-evluate iLevel based on this eExcitation
417  iLevel = 0;
418  G4bool find = false;
419  G4int imaxEx = 0;
420  G4double level_tolerance = 1.0*CLHEP::keV;
421 
422  while( theGammas.GetLevel(iLevel+1) != 0 ) { // Loop checking, 11.05.2015, T. Koi
423  G4double maxEx = 0.0;
424  if (maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) {
425  maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();
426  imaxEx = iLevel;
427  }
428 
429  // Fix bug 1789 DHW - first if-branch added because gamma data come from ENSDF
430  // and do not necessarily match the excitations used in ENDF-B.VII
431  // Compromise solution: use 1 keV tolerance suggested by T. Koi
432  if (std::abs(eExcitation - theGammas.GetLevel(iLevel)->GetLevelEnergy() ) < level_tolerance) {
433  find = true;
434  break;
435 
436  } else if (eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) {
437  find = true;
438  iLevel--;
439  // very small eExcitation, iLevel becomes -1, this is protected below
440  if (theTrack.GetDefinition() == aDefinition) { // this line added as part of fix #1838
441  if (iLevel == -1) iLevel = 0;
442  }
443  break;
444  }
445  iLevel++;
446  }
447 
448  // If proper level cannot be found, use the maximum level
449  if (!find) iLevel = imaxEx;
450  }
451 
452  if(std::getenv("G4ParticleHPDebug") && eKinetic-eExcitation < 0)
453  {
454  throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
455  }
456  if(eKinetic-eExcitation < 0) eExcitation = 0;
457  if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
458 
459  }
461 
462  if( theFinalStatePhotons[it] == 0 )
463  {
464  //G4cout << "110610 USE Gamma Level" << G4endl;
465 // TK comment Most n,n* eneter to this
466  thePhotons = theGammas.GetDecayGammas(iLevel);
467  eGamm -= theGammas.GetLevelEnergy(iLevel);
468  if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
469  {
470  G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
471  theRestEnergy->SetDefinition(G4Gamma::Gamma());
472  theRestEnergy->SetKineticEnergy(eGamm);
473  G4double costh = 2.*G4UniformRand()-1.;
475  theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
476  eGamm*std::sin(std::acos(costh))*std::sin(phi),
477  eGamm*costh);
478  if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
479  thePhotons->push_back(theRestEnergy);
480  }
481  }
482  }
483  else if(theEnergyAngData[it] != 0) // MF6
484  {
485 
486  theParticles = theEnergyAngData[it]->Sample(eKinetic);
487 
488  //141017 Fix BEGIN
489  //Adjust A and Z in the case of miss much between selected data and target nucleus
490  if ( theParticles != NULL ) {
491  G4int sumA = 0;
492  G4int sumZ = 0;
493  G4int maxA = 0;
494  G4int jAtMaxA = 0;
495  for ( G4int j = 0 ; j != (G4int)theParticles->size() ; j++ ) {
496  if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
497  maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber();
498  jAtMaxA = j;
499  }
500  sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
501  sumZ += G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() + eps );
502  }
503  G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
504  G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
505  if ( dA < 0 || dZ < 0 ) {
506  G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
507  G4int newZ = G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
508  G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
509  theParticles->at( jAtMaxA )->SetDefinition( pd );
510  }
511  }
512  //141017 Fix END
513 
514  }
515  else
516  {
517  // @@@ what to do, if we have photon data, but no info on the hadron itself
518  nothingWasKnownOnHadron = 1;
519  }
520 
521  //G4cout << "theFinalStatePhotons it " << it << G4endl;
522  //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
523  //G4cout << "theFinalStatePhotons it " << it << G4endl;
524  //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
525  //G4cout << "thePhotons " << thePhotons << G4endl;
526 
527  if ( theFinalStatePhotons[it] != 0 )
528  {
529  // the photon distributions are in the Nucleus rest frame.
530  // TK residual rest frame
531  G4ReactionProduct boosted_tmp;
532  boosted_tmp.Lorentz(incidReactionProduct, theTarget);
533  G4double anEnergy = boosted_tmp.GetKineticEnergy();
534  thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
536  G4double testEnergy = 0;
537  if(thePhotons!=0 && thePhotons->size()!=0)
538  { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
539  if(theFinalStatePhotons[it]->NeedsCascade())
540  {
541  while(aBaseEnergy>0.01*CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
542  {
543  // cascade down the levels
544  G4bool foundMatchingLevel = false;
545  G4int closest = 2;
546  G4double deltaEold = -1;
547  for(G4int j=1; j<it; j++)
548  {
549  if(theFinalStatePhotons[j]!=0)
550  {
551  testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
552  }
553  else
554  {
555  testEnergy = 0;
556  }
557  G4double deltaE = std::abs(testEnergy-aBaseEnergy);
558  if(deltaE<0.1*CLHEP::keV)
559  {
560  G4ReactionProductVector * theNext =
561  theFinalStatePhotons[j]->GetPhotons(anEnergy);
562  if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
563  aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
564  delete theNext;
565  foundMatchingLevel = true;
566  break; // ===>
567  }
568  if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
569  {
570  closest = j;
571  deltaEold = deltaE;
572  }
573  } // <=== the break goes here.
574  if(!foundMatchingLevel)
575  {
576  G4ReactionProductVector * theNext =
577  theFinalStatePhotons[closest]->GetPhotons(anEnergy);
578  if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
579  aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
580  delete theNext;
581  }
582  }
583  }
584  }
585  unsigned int i0;
586  if(thePhotons!=0)
587  {
588  for(i0=0; i0<thePhotons->size(); i0++)
589  {
590  // back to lab
591  thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
592  }
593  }
594  //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
595  if(nothingWasKnownOnHadron)
596  {
597 // In this case, hadron should be isotropic in CM
598 // Next 12 lines are Emilio's replacement
599  // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
600  // G4double eExcitation = QM-QI[it];
601  G4double eExcitation = QI[0] - QI[it]; // Fix of bug #1838
602  if(eExcitation<20*CLHEP::keV){eExcitation=0;}
603  two_body_reaction(&incidReactionProduct,&theTarget,&aHadron,eExcitation);
604  if(thePhotons==0 && eExcitation>0){
605  for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
606  {
607  if(theGammas.GetLevelEnergy(iLevel)<eExcitation+5*keV) break; // 5 keV tolerance
608  }
609  thePhotons = theGammas.GetDecayGammas(iLevel);
610  }
611  }
612 // Emilio's replacement done
613 /*
614 // This code replaced by Emilio (previous 12 lines)
615 // mu and p should be correlated
616 //
617  //isotropic distribution in CM
618  G4double mu = 1.0 - 2.*G4UniformRand();
619 
620  // Need momenta in target rest frame
621  G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
622  G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
623  G4LorentzVector proj_in_LAB = hadProjectile->Get4Momentum();
624 
625  G4DynamicParticle* proj = new G4DynamicParticle(theProjectile, proj_in_LAB.boost(boostToTargetRest) );
626 // G4DynamicParticle* targ =
627 // new G4DynamicParticle(G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, totalPhotonEnergy), G4ThreeVector(0) );
628 // Fix bug 2166 (A. Zontikov): replace above two lines with next three lines
629  G4double excitationEnergy = theFinalStatePhotons[it] ? theFinalStatePhotons[it]->GetLevelEnergy() : 0.0;
630  G4DynamicParticle* targ =
631  new G4DynamicParticle(G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, excitationEnergy), G4ThreeVector(0) );
632  G4DynamicParticle* hadron =
633  new G4DynamicParticle(aHadron.GetDefinition(), G4ThreeVector(0) ); // Will fill in the momentum
634 
635  two_body_reaction ( proj , targ , hadron , mu );
636 
637  G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
638  G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
639  aHadron.SetMomentum( hadron_in_LAB.v() );
640  aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
641 
642  delete proj;
643  delete targ;
644  delete hadron;
645 
646  }
647 */
648 
649 // fill the result
650 // Beware - the recoil is not necessarily in the particles...
651 // Can be calculated from momentum conservation?
652 // The idea is that the particles ar emitted forst, and the gammas only once the
653 // recoil is on the residual; assumption is that gammas do not contribute to
654 // the recoil.
655 // This needs more design @@@
656 
657  G4int nSecondaries = 2; // the hadron and the recoil
658  G4bool needsSeparateRecoil = false;
659  G4int totalBaryonNumber = 0;
660  G4int totalCharge = 0;
661  G4ThreeVector totalMomentum(0);
662  if(theParticles != 0)
663  {
664  nSecondaries = theParticles->size();
665  const G4ParticleDefinition * aDef;
666  unsigned int ii0;
667  for(ii0=0; ii0<theParticles->size(); ii0++)
668  {
669  aDef = theParticles->operator[](ii0)->GetDefinition();
670  totalBaryonNumber+=aDef->GetBaryonNumber();
671  totalCharge+=G4int(aDef->GetPDGCharge()+eps);
672  totalMomentum += theParticles->operator[](ii0)->GetMomentum();
673  }
674  if(totalBaryonNumber!=G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()))
675  {
676  needsSeparateRecoil = true;
677  nSecondaries++;
678  residualA = G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()
679  -totalBaryonNumber);
680  residualZ = G4int(theBaseZ+eps+hadProjectile->GetDefinition()->GetPDGCharge()
681  -totalCharge);
682  }
683  }
684 
685  G4int nPhotons = 0;
686  if(thePhotons!=0) { nPhotons = thePhotons->size(); }
687  nSecondaries += nPhotons;
688 
689  G4DynamicParticle * theSec;
690 
691  if( theParticles==0 )
692  {
693  theSec = new G4DynamicParticle;
694  theSec->SetDefinition(aHadron.GetDefinition());
695  theSec->SetMomentum(aHadron.GetMomentum());
696  theResult.Get()->AddSecondary(theSec);
697 #ifdef G4PHPDEBUG
698  if( std::getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary1 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
699 #endif
700 
701  aHadron.Lorentz(aHadron, theTarget);
702  G4ReactionProduct theResidual;
704  ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
705  theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
706 
707  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
708  //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
709  G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
710  theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
711 
712  theResidual.Lorentz(theResidual, -1.*theTarget);
713  G4ThreeVector totalPhotonMomentum(0,0,0);
714  if(thePhotons!=0)
715  {
716  for(i=0; i<nPhotons; i++)
717  {
718  totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
719  }
720  }
721  theSec = new G4DynamicParticle;
722  theSec->SetDefinition(theResidual.GetDefinition());
723  theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
724  theResult.Get()->AddSecondary(theSec);
725 #ifdef G4PHPDEBUG
726  if( std::getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
727 #endif
728  }
729  else
730  {
731  for(i0=0; i0<theParticles->size(); i0++)
732  {
733  theSec = new G4DynamicParticle;
734  theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
735  theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
736  theResult.Get()->AddSecondary(theSec);
737 #ifdef G4PHPDEBUG
738  if( std::getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
739 #endif
740  delete theParticles->operator[](i0);
741  }
742  delete theParticles;
743  if(needsSeparateRecoil && residualZ!=0)
744  {
745  G4ReactionProduct theResidual;
747  ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
748  G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass();
749  resiualKineticEnergy += totalMomentum*totalMomentum;
750  resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
751 // cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
752  theResidual.SetKineticEnergy(resiualKineticEnergy);
753 
754  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
755  //theResidual.SetMomentum(-1.*totalMomentum);
756  //G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
757  //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
758 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
759  theResidual.SetMomentum( incidReactionProduct.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
760 
761  theSec = new G4DynamicParticle;
762  theSec->SetDefinition(theResidual.GetDefinition());
763  theSec->SetMomentum(theResidual.GetMomentum());
764  theResult.Get()->AddSecondary(theSec);
765 #ifdef G4PHPDEBUG
766  if( std::getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
767 #endif
768 
769  }
770  }
771  if(thePhotons!=0)
772  {
773  for(i=0; i<nPhotons; i++)
774  {
775  theSec = new G4DynamicParticle;
776  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
777  //theSec->SetDefinition(G4Gamma::Gamma());
778  theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
779  //But never cause real effect at least with G4NDL3.13 TK
780  theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
781  theResult.Get()->AddSecondary(theSec);
782 #ifdef G4PHPDEBUG
783  if( std::getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
784 #endif
785 
786  delete thePhotons->operator[](i);
787  }
788 // some garbage collection
789  delete thePhotons;
790  }
791 
792 //080721
794  G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
795  G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
796  G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
797  adjust_final_state ( init_4p_lab );
798 
799 // clean up the primary neutron
801 }
802 
803 
804 
805 //Re-implemented by E. Mendoza (2019). Isotropic emission in the CMS:
806 // proj: projectile in target-rest-frame (input)
807 // targ: target in target-rest-frame (input)
808 // product: secondary particle in target-rest-frame (output)
809 // resExcitationEnergy: excitation energy of the residual nucleus
810 
813  G4ReactionProduct* product,
814  G4double resExcitationEnergy)
815 {
816  //CMS system:
817  G4ReactionProduct theCMS= *proj+ *targ;
818 
819  //Residual definition:
820  G4int resZ=(G4int)(proj->GetDefinition()->GetPDGCharge()+targ->GetDefinition()->GetPDGCharge()-product->GetDefinition()->GetPDGCharge()+0.1);
821  G4int resA=proj->GetDefinition()->GetBaryonNumber()+targ->GetDefinition()->GetBaryonNumber()-product->GetDefinition()->GetBaryonNumber();
822  G4ReactionProduct theResidual;
823  theResidual.SetDefinition(G4IonTable::GetIonTable()->GetIon(resZ,resA,0.0));
824 
825  //CMS system:
826  G4ReactionProduct theCMSproj;
827  G4ReactionProduct theCMStarg;
828  theCMSproj.Lorentz(*proj,theCMS);
829  theCMStarg.Lorentz(*targ,theCMS);
830  //final Momentum in the CMS:
831  G4double totE=std::sqrt(theCMSproj.GetMass()*theCMSproj.GetMass()+theCMSproj.GetTotalMomentum()*theCMSproj.GetTotalMomentum())+std::sqrt(theCMStarg.GetMass()*theCMStarg.GetMass()+theCMStarg.GetTotalMomentum()*theCMStarg.GetTotalMomentum());
832  G4double prodmass=product->GetMass();
833  G4double resmass=theResidual.GetMass()+resExcitationEnergy;
834  G4double fmomsquared=1./4./totE/totE*(totE*totE-(prodmass-resmass)*(prodmass-resmass))*(totE*totE-(prodmass+resmass)*(prodmass+resmass));
835  G4double fmom=0;
836  if(fmomsquared>0){
837  fmom=std::sqrt(fmomsquared);
838  }
839 
840  //random (isotropic direction):
841  G4double cosTh = 2.*G4UniformRand()-1.;
843  G4double theta = std::acos(cosTh);
844  G4double sinth = std::sin(theta);
845  product->SetMomentum(fmom*sinth*std::cos(phi),fmom*sinth*std::sin(phi),fmom*cosTh); //CMS
846  product->SetTotalEnergy(std::sqrt(prodmass*prodmass+fmom*fmom)); //CMS
847  //Back to the LAB system:
848  product->Lorentz(*product,-1.*theCMS);
849 
850 }
851 
852 
854 {
855  if ( aDefinition == G4Neutron::Definition() ) // If the outgoing particle is a neutron...
856  {
857  // LR: flag LR in ENDF. It indicates whether there is breakup of the residual nucleus or not.
858  // it: exit channel (index of the carbon excited state)
859 
860  //if ( (G4int)(theBaseZ+0.1) == 6 ) G4cout << "LR[" << it << "] = " << LR[it] << G4endl;
861 
862  // Added by A. R. García (CIEMAT) to include the physics of C(N,N'3A) reactions from NRESP71.
863 
864  if ( LR[it] > 0 ) // If there is breakup of the residual nucleus LR(flag LR in ENDF)>0 (i.e. Z=6 MT=52-91 (it=MT-50)).
865  {
866  // Defining carbon as the target in the reference frame at rest.
867  G4ReactionProduct theCarbon(theTarget);
868 
869  theCarbon.SetMomentum(G4ThreeVector());
870  theCarbon.SetKineticEnergy(0.);
871 
872  // Creating four reaction products.
873  G4ReactionProduct theProds[4];
874 
875  // Applying C(N,N'3A) reaction mechanisms in the target rest frame.
876  if ( it == 41 )
877  {
878  // QI=QM=-7.275 MeV for C-0(N,N')C-C(3A) in ENDF/B-VII.1. This is not the value of the QI of the first step according to the model. So we don't take it. Instead, we set the one we have calculated: QI=(mn+m12C)-(ma+m9Be+Ex9Be)=-8.130 MeV.
879  nresp71_model.ApplyMechanismI_NBeA2A(boosted, theCarbon, theProds, -8.130/*QI[it]*/); // N+C --> A[0]+9BE* | 9BE* --> N[1]+8BE | 8BE --> 2*A[2,3].
880  //printf("- QI=%f\n", QI[it]);
881  }
882  else
883  {
884  nresp71_model.ApplyMechanismII_ACN2A(boosted, theCarbon, theProds, QI[it]); // N+C --> N'[0]+C* | C* --> A[1]+8BE | 8BE --> 2*A[2,3].
885  }
886 
887  //printf("it=%d qi=%f \n", it, QI[it]);
888 
889  // Returning to the reference frame where the target was in motion.
890  for ( G4int j=0; j<4; j++ )
891  {
892  theProds[j].Lorentz(theProds[j], -1.*theTarget);
893  theResult.Get()->AddSecondary(new G4DynamicParticle(theProds[j].GetDefinition(), theProds[j].GetMomentum()));
894  }
895 
896  /*G4double EN0 = theNeutron.GetKineticEnergy();
897  G4double EN1 = theProds[0].GetKineticEnergy();
898 
899  G4double EA1 = theProds[1].GetKineticEnergy();
900  G4double EA2 = theProds[2].GetKineticEnergy();
901  G4double EA3 = theProds[3].GetKineticEnergy();
902 
903  printf("Q=%f\n", EN1+EA1+EA2+EA3-EN0);*/
904 
905  // Killing the primary neutron.
907 
908  return true;
909  }
910  }
911  else if ( aDefinition == G4Alpha::Definition() ) // If the outgoing particle is an alpha, ...
912  {
913  // Added by A. R. García (CIEMAT) to include the physics of C(N,A)9BE reactions from NRESP71.
914 
915  if ( LR[it] == 0 ) // If Z=6, an alpha particle is emitted and there is no breakup of the residual nucleus LR(flag LR in ENDF)==0.
916  {
917  // Defining carbon as the target in the reference frame at rest.
918  G4ReactionProduct theCarbon(theTarget);
919 
920  theCarbon.SetMomentum(G4ThreeVector());
921  theCarbon.SetKineticEnergy(0.);
922 
923  // Creating four reaction products.
924  G4ReactionProduct theProds[2];
925 
926  // Applying C(N,A)9BE reaction mechanism.
927  nresp71_model.ApplyMechanismABE(boosted, theCarbon, theProds); // N+C --> A[0]+9BE[1].
928 
929  //G4DynamicParticle *theSec;
930  for ( G4int j=0; j<2; j++ )
931  {
932  // Returning to the system of reference where the target was in motion.
933  theProds[j].Lorentz(theProds[j], -1.*theTarget);
934  theResult.Get()->AddSecondary(new G4DynamicParticle(theProds[j].GetDefinition(), theProds[j].GetMomentum()));
935  }
936 
937  // Killing the primary neutron.
939 
940  return true;
941  }
942  else
943  {
944  G4Exception("G4ParticleHPInelasticCompFS::CompositeApply()", "G4ParticleInelasticCompFS.cc", FatalException, "Alpha production with LR!=0.");
945  }
946  }
947 
948  return false;
949 }