Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPInelasticBaseFS.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 // 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 // 101111 Add Special treatment for Be9(n,2n)Be8(2a) case by T. Koi
34 //
35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 //
38 #include "G4ParticleHPManager.hh"
39 #include "G4Nucleus.hh"
40 #include "G4NucleiProperties.hh"
41 #include "G4He3.hh"
42 #include "G4Alpha.hh"
43 #include "G4Electron.hh"
44 #include "G4ParticleHPDataUsed.hh"
45 
46 #include "G4IonTable.hh"
47 
49 {
50  // char the[100] = {""};
51  // std::ostrstream ost(the, 100, std::ios::out);
52  // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
53  // G4String * aName = new G4String(the);
54  // std::ifstream from(*aName, std::ios::in);
55 
56  std::ostringstream ost;
57  ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
58  G4String aName = ost.str();
59  std::ifstream from(aName, std::ios::in);
60 
61  if(!from) return; // no data found for this isotope
62  // std::ifstream theGammaData(*aName, std::ios::in);
63  std::ifstream theGammaData(aName, std::ios::in);
64 
65  G4double eps = 0.001;
67  G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(AR+eps),static_cast<G4int>(ZR+eps)) -
68  G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps));
69  theGammas.Init(theGammaData);
70  // delete aName;
71 
72 }
73 
75 {
76  gammaPath = "/Inelastic/Gammas/";
77  if(!getenv("G4NEUTRONHPDATA"))
78  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
79  G4String tBase = getenv("G4NEUTRONHPDATA");
80  gammaPath = tBase+gammaPath;
81  G4String tString = dirName;
82  G4bool dbool;
83  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M,tString, bit, dbool);
84  G4String filename = aFile.GetName();
85 #ifdef G4PHPDEBUG
86  if( getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPInelasticBaseFS::Init FILE " << filename << G4endl;
87 #endif
88  SetAZMs( A, Z, M, aFile);
89  //theBaseA = aFile.GetA();
90  //theBaseZ = aFile.GetZ();
91  // theNDLDataA = (int)aFile.GetA();
92  // theNDLDataZ = aFile.GetZ();
93  //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
94  if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
95  {
96 #ifdef G4PHPDEBUG
97  if(getenv("G4ParticleHPDebug_NamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
98 #endif
99  hasAnyData = false;
100  hasFSData = false;
101  hasXsec = false;
102  return;
103  }
104  // theBaseA = A;
105  // theBaseZ = G4int(Z+.5);
106  //std::ifstream theData(filename, std::ios::in);
107  //130205 For compressed data files
108  std::istringstream theData(std::ios::in);
109  G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
110  //130205 END
111  if(!theData) //"!" is a operator of ios
112  {
113  hasAnyData = false;
114  hasFSData = false;
115  hasXsec = false;
116  // theData.close();
117  return; // no data for exactly this isotope and FS
118  }
119  // here we go
120  G4int infoType, dataType, dummy=INT_MAX;
121  hasFSData = false;
122  while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
123  {
124  theData >> dataType;
125 
126  if(dummy==INT_MAX) theData >> dummy >> dummy;
127  if(dataType==3)
128  {
129  G4int total;
130  theData >> total;
131  theXsection->Init(theData, total, CLHEP::eV);
132  }
133  else if(dataType==4)
134  {
136  theAngularDistribution->Init(theData);
137  hasFSData = true;
138  }
139  else if(dataType==5)
140  {
142  theEnergyDistribution->Init(theData);
143  hasFSData = true;
144  }
145  else if(dataType==6)
146  {
148  //- G4cout << this << " BaseFS theEnergyAngData " << theEnergyAngData << G4endl; //GDEB
149  theEnergyAngData->Init(theData);
150  hasFSData = true;
151  }
152  else if(dataType==12)
153  {
155  theFinalStatePhotons->InitMean(theData);
156  hasFSData = true;
157  }
158  else if(dataType==13)
159  {
162  hasFSData = true;
163  }
164  else if(dataType==14)
165  {
167  hasFSData = true;
168  }
169  else if(dataType==15)
170  {
172  hasFSData = true;
173  }
174  else
175  {
176  throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4ParticleHPInelasticBaseFS");
177  }
178  }
179  // theData.close();
180 }
181 
183  G4ParticleDefinition ** theDefs,
184  G4int nDef)
185 {
186 
187 // prepare neutron
188  if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
189  theResult.Get()->Clear();
190  G4double eKinetic = theTrack.GetKineticEnergy();
191  const G4HadProjectile *hadProjectile = &theTrack;
192  G4ReactionProduct incidReactionProduct( const_cast<G4ParticleDefinition *>(hadProjectile->GetDefinition()) );
193  incidReactionProduct.SetMomentum( hadProjectile->Get4Momentum().vect() );
194  incidReactionProduct.SetKineticEnergy( eKinetic );
195 
196 // prepare target
197  G4double targetMass;
198  G4double eps = 0.0001;
199  targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
200  //theProjectile->GetPDGMass();
202 
203  //give priority to ENDF vales for target mass
204  if(theEnergyAngData!=0)
205  { targetMass = theEnergyAngData->GetTargetMass(); }
207  { targetMass = theAngularDistribution->GetTargetMass(); }
208 
209 //080731a
210 //110512 TKDB ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np) targetMass is not properly recorded.
211 //if ( 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;
212  if ( targetMass == 0 )
213  {
214  //G4cout << "TKDB targetMass = 0; ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np) targetMass is not properly recorded. This could be a similar situation." << G4endl;
215  //targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) / theProjectile->GetPDGMass();
216  targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) / G4Neutron::Neutron()->GetPDGMass();
217  }
218 
219  G4Nucleus aNucleus;
221  //G4ThreeVector neuVelo = (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum();
222  //theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
223  //G4Nucleus::GetBiasedThermalNucleus requests normalization of mass and velocity in neutron mass
224  G4ThreeVector neuVelo = (1./G4Neutron::Neutron()->GetPDGMass())*incidReactionProduct.GetMomentum();
225  theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
226 
227  theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) );
228 
229 // prepare energy in target rest frame
230  G4ReactionProduct boosted;
231  boosted.Lorentz(incidReactionProduct, theTarget);
232  eKinetic = boosted.GetKineticEnergy();
233  G4double orgMomentum = boosted.GetMomentum().mag();
234 
235 // Take N-body phase-space distribution, if no other data present.
236  if(!HasFSData()) // adding the residual is trivial here @@@
237  {
238  G4ParticleHPNBodyPhaseSpace thePhaseSpaceDistribution;
239  G4double aPhaseMass=0;
240  G4int ii;
241  for(ii=0; ii<nDef; ii++)
242  {
243  aPhaseMass+=theDefs[ii]->GetPDGMass();
244  }
245  thePhaseSpaceDistribution.Init(aPhaseMass, nDef);
246  thePhaseSpaceDistribution.SetProjectileRP(&incidReactionProduct);
247  thePhaseSpaceDistribution.SetTarget(&theTarget);
248  for(ii=0; ii<nDef; ii++)
249  {
250  G4double massCode = 1000.*std::abs(theDefs[ii]->GetPDGCharge());
251  massCode += theDefs[ii]->GetBaryonNumber();
252  G4double dummy = 0;
253  G4ReactionProduct * aSec = thePhaseSpaceDistribution.Sample(eKinetic, massCode, dummy);
254  aSec->Lorentz(*aSec, -1.*theTarget);
255  G4DynamicParticle * aPart = new G4DynamicParticle();
256  aPart->SetDefinition(aSec->GetDefinition());
257  aPart->SetMomentum(aSec->GetMomentum());
258  delete aSec;
259  theResult.Get()->AddSecondary(aPart);
260 #ifdef G4PHPDEBUG
261  if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply NoFSData add secondary " << aPart->GetParticleDefinition()->GetParticleName() << " E= " << aPart->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
262 #endif
263  }
265  //TK120607
266  //Final momentum check should be done before return
268  G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
269  G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
270  G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
271  adjust_final_state ( init_4p_lab );
272 
273  return;
274  }
275 
276 // set target and neutron in the relevant exit channel
277  if(theAngularDistribution!=0)
278  {
279  theAngularDistribution->SetTarget(theTarget);
280  theAngularDistribution->SetProjectileRP(incidReactionProduct);
281  }
282  else if(theEnergyAngData!=0)
283  {
284  theEnergyAngData->SetTarget(theTarget);
285  theEnergyAngData->SetProjectileRP(incidReactionProduct);
286  }
287 
288  G4ReactionProductVector * tmpHadrons = 0;
289 #ifdef G4PHPDEBUG
290  //To avoid compilation error around line 532.
291  G4int ii(0);
292 #endif
293  G4int dummy;
294  unsigned int i;
295 
296  if(theEnergyAngData != 0)
297  {
298  tmpHadrons = theEnergyAngData->Sample(eKinetic);
299 
300  if ( !getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) {
301  //141017 Fix BEGIN
302  //Adjust A and Z in the case of miss much between selected data and target nucleus
303  if ( tmpHadrons != NULL ) {
304  G4int sumA = 0;
305  G4int sumZ = 0;
306  G4int maxA = 0;
307  G4int jAtMaxA = 0;
308  for ( G4int j = 0 ; j != (G4int)tmpHadrons->size() ; j++ ) {
309  //G4cout << __FILE__ << " " << __LINE__ << "th line: tmpHadrons->at(j)->GetDefinition()->GetParticleName() = " << tmpHadrons->at(j)->GetDefinition()->GetParticleName() << G4endl;
310  if ( tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
311  maxA = tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
312  jAtMaxA = j;
313  }
314  sumA += tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
315  sumZ += G4int( tmpHadrons->at(j)->GetDefinition()->GetPDGCharge() + eps );
316  }
317  G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
318  G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
319  if ( dA < 0 || dZ < 0 ) {
320  G4int newA = tmpHadrons->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
321  G4int newZ = G4int( tmpHadrons->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
322  G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
323  tmpHadrons->at( jAtMaxA )->SetDefinition( pd );
324  }
325  }
326  //141017 Fix END
327  }
328  }
329  else if(theAngularDistribution!= 0)
330  {
331  G4bool * Done = new G4bool[nDef];
332  G4int i0;
333  for(i0=0; i0<nDef; i0++) Done[i0] = false;
334 
335 //Following lines are commented out to fix coverity defeat 58622
336 // if(tmpHadrons == 0)
337 // {
338  tmpHadrons = new G4ReactionProductVector;
339 // }
340 // else
341 // {
342 // for(i=0; i<tmpHadrons->size(); i++)
343 // {
344 // for(ii=0; ii<nDef; ii++)
345 // if(!Done[ii] && tmpHadrons->operator[](i)->GetDefinition() == theDefs[ii])
346 // Done[ii] = true;
347 // }
348 #ifdef G4PHPDEBUG
349 // if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply secondary previously added " << tmpHadrons->operator[](i)->GetDefinition()->GetParticleName() << " E= " << tmpHadrons->operator[](i)->GetKineticEnergy() << G4endl;
350 #endif
351 // }
352  G4ReactionProduct * aHadron;
353  G4double localMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps)));
354  G4ThreeVector bufferedDirection(0,0,0);
355  for(i0=0; i0<nDef; i0++)
356  {
357  if(!Done[i0])
358  {
359  aHadron = new G4ReactionProduct;
360  if(theEnergyDistribution!=0)
361  {
362  aHadron->SetDefinition(theDefs[i0]);
363  aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
364  }
365  else if(nDef == 1)
366  {
367  aHadron->SetDefinition(theDefs[i0]);
368  aHadron->SetKineticEnergy(eKinetic);
369  }
370  else if(nDef == 2)
371  {
372  aHadron->SetDefinition(theDefs[i0]);
373  aHadron->SetKineticEnergy(50*CLHEP::MeV);
374  }
375  else
376  {
377  throw G4HadronicException(__FILE__, __LINE__, "No energy distribution to sample from in InelasticBaseFS::BaseApply");
378  }
380  if(theEnergyDistribution==0 && nDef == 2)
381  {
382  if(i0==0)
383  {
384  G4double mass1 = theDefs[0]->GetPDGMass();
385  G4double mass2 = theDefs[1]->GetPDGMass();
386  G4double massn = theProjectile->GetPDGMass();
387  G4int z1 = static_cast<G4int>(theBaseZ+eps-theDefs[0]->GetPDGCharge()-theDefs[1]->GetPDGCharge());
388  G4int a1 = static_cast<G4int>(theBaseA+eps)-theDefs[0]->GetBaryonNumber()-theDefs[1]->GetBaryonNumber();
389  G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
390  G4double availableEnergy = eKinetic+massn+localMass-mass1-mass2-concreteMass;
391  // available kinetic energy in CMS (non relativistic)
392  G4double emin = availableEnergy+mass1+mass2 - std::sqrt((mass1+mass2)*(mass1+mass2)+orgMomentum*orgMomentum);
393  G4double p1=std::sqrt(2.*mass2*emin);
394  bufferedDirection = p1*aHadron->GetMomentum().unit();
395 #ifdef G4PHPDEBUG
396  if(getenv("G4ParticleHPDebug")) // @@@@@ verify the nucleon counting...
397  {
398  G4cout << "G4ParticleHPInelasticBaseFS "<<z1<<" "<<theBaseZ<<" "<<a1<<" "<<theBaseA<<" "<<availableEnergy<<" "
399  << emin<<G4endl;
400  }
401 #endif
402  }
403  else
404  {
405  bufferedDirection = -bufferedDirection;
406  }
407  // boost from cms to lab
408 #ifdef G4PHPDEBUG
409  if(getenv("G4ParticleHPDebug"))
410  {
411  G4cout << " G4ParticleHPInelasticBaseFS "<<bufferedDirection.mag2()<<G4endl;
412  }
413 #endif
414  aHadron->SetTotalEnergy( std::sqrt(aHadron->GetMass()*aHadron->GetMass()
415  +bufferedDirection.mag2()) );
416  aHadron->SetMomentum(bufferedDirection);
417  aHadron->Lorentz(*aHadron, -1.*(theTarget+incidReactionProduct));
418 #ifdef G4PHPDEBUG
419  if(getenv("G4ParticleHPDebug"))
420  {
421  G4cout << " G4ParticleHPInelasticBaseFS "<<aHadron->GetTotalEnergy()<<" "<<aHadron->GetMomentum()<<G4endl;
422  }
423 #endif
424  }
425  tmpHadrons->push_back(aHadron);
426 #ifdef G4PHPDEBUG
427  if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply FSData add secondary " << aHadron->GetDefinition()->GetParticleName() << " E= " << aHadron->GetKineticEnergy() << G4endl;
428 #endif
429  }
430  }
431  delete [] Done;
432  }
433  else
434  {
435  throw G4HadronicException(__FILE__, __LINE__, "No data to create the neutrons in NInelasticFS");
436  }
437 
438  G4ReactionProductVector * thePhotons = 0;
439  if(theFinalStatePhotons!=0)
440  {
441  // the photon distributions are in the Nucleus rest frame.
442  G4ReactionProduct boosted_tmp;
443  boosted_tmp.Lorentz(incidReactionProduct, theTarget);
444  G4double anEnergy = boosted_tmp.GetKineticEnergy();
445  thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
446  if(thePhotons!=0)
447  {
448  for(i=0; i<thePhotons->size(); i++)
449  {
450  // back to lab
451  thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1.*theTarget);
452  }
453  }
454  }
455  else if(theEnergyAngData!=0)
456  {
457 
458  // PA130927: do not create photons to adjust binding energy
459  G4bool bAdjustPhotons = true;
460 #ifdef PHP_AS_HP
461  bAdjustPhotons = true;
462 #else
463  if ( getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) bAdjustPhotons = false;
464 #endif
465 
466  if( bAdjustPhotons ) {
467  G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
468  G4double anEnergy = boosted.GetKineticEnergy();
469  theGammaEnergy = anEnergy-theGammaEnergy;
470  theGammaEnergy += theNuclearMassDifference;
471  G4double eBindProducts = 0;
472  G4double eBindN = 0;
473  G4double eBindP = 0;
478  G4int ia=0;
479  for(i=0; i<tmpHadrons->size(); i++)
480  {
481  if(tmpHadrons->operator[](i)->GetDefinition() == G4Neutron::Neutron())
482  {
483  eBindProducts+=eBindN;
484  }
485  else if(tmpHadrons->operator[](i)->GetDefinition() == G4Proton::Proton())
486  {
487  eBindProducts+=eBindP;
488  }
489  else if(tmpHadrons->operator[](i)->GetDefinition() == G4Deuteron::Deuteron())
490  {
491  eBindProducts+=eBindD;
492  }
493  else if(tmpHadrons->operator[](i)->GetDefinition() == G4Triton::Triton())
494  {
495  eBindProducts+=eBindT;
496  }
497  else if(tmpHadrons->operator[](i)->GetDefinition() == G4He3::He3())
498  {
499  eBindProducts+=eBindHe3;
500  }
501  else if(tmpHadrons->operator[](i)->GetDefinition() == G4Alpha::Alpha())
502  {
503  eBindProducts+=eBindA;
504  ia++;
505  }
506  }
507 
508  theGammaEnergy += eBindProducts;
509 
510 #ifdef G4PHPDEBUG
511  if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply gamma Energy " << theGammaEnergy << " eBindProducts " << eBindProducts << G4endl;
512 #endif
513 
514  //101111
515  //Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
516  if ( (G4int)(theBaseZ+eps) == 4 && (G4int)(theBaseA+eps) == 9 )
517  {
518  // This only valid for G4NDL3.13,,,
519  if ( std::abs( theNuclearMassDifference -
522  && ia == 2 )
523  {
524  theGammaEnergy -= (2*eBindA);
525  }
526  }
527 
528  G4ReactionProductVector * theOtherPhotons = 0;
529  G4int iLevel;
530  while(theGammaEnergy>=theGammas.GetLevelEnergy(0)) // Loop checking, 11.05.2015, T. Koi
531  {
532  for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
533  {
534  if(theGammas.GetLevelEnergy(iLevel)<theGammaEnergy) break;
535  }
536  if(iLevel==0||iLevel==theGammas.GetNumberOfLevels()-1)
537  {
538  theOtherPhotons = theGammas.GetDecayGammas(iLevel);
539 #ifdef G4PHPDEBUG
540  if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply adding gamma from level " << iLevel << theOtherPhotons->operator[](ii)->GetKineticEnergy() << G4endl;
541 #endif
542  }
543  else
544  {
545  G4double random = G4UniformRand();
546  G4double eLow = theGammas.GetLevelEnergy(iLevel);
547  G4double eHigh = theGammas.GetLevelEnergy(iLevel+1);
548  if(random > (eHigh-eLow)/(theGammaEnergy-eLow)) iLevel++;
549  theOtherPhotons = theGammas.GetDecayGammas(iLevel);
550  }
551  if(thePhotons==0) thePhotons = new G4ReactionProductVector;
552  if(theOtherPhotons != 0)
553  {
554  for(unsigned int iii=0; iii<theOtherPhotons->size(); iii++)
555  {
556  thePhotons->push_back(theOtherPhotons->operator[](iii));
557 #ifdef G4PHPDEBUG
558  if( getenv("G4ParticleHPDebug"))
559  G4cout << iii << " G4ParticleHPInelasticBaseFS::BaseApply adding gamma " << theOtherPhotons->operator[](iii)->GetKineticEnergy() << G4endl;
560 #endif
561  }
562  delete theOtherPhotons;
563  }
564  theGammaEnergy -= theGammas.GetLevelEnergy(iLevel);
565  if(iLevel == -1) break;
566  }
567  }
568  }
569 
570  // fill the result
571  unsigned int nSecondaries = tmpHadrons->size();
572  unsigned int nPhotons = 0;
573  if(thePhotons!=0) { nPhotons = thePhotons->size(); }
574  nSecondaries += nPhotons;
575  G4DynamicParticle * theSec;
576 #ifdef G4PHPDEBUG
577  if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N hadrons " << nSecondaries-nPhotons << G4endl;
578 #endif
579 
580  for(i=0; i<nSecondaries-nPhotons; i++)
581  {
582  theSec = new G4DynamicParticle;
583  theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
584  theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
585  theResult.Get()->AddSecondary(theSec);
586 #ifdef G4PHPDEBUG
587  if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
588 #endif
589  if( getenv("G4PHPTEST") ) G4cout << " InelasticBaseFS COS THETA " << std::cos(theSec->GetMomentum().theta()) << " " << (theSec->GetMomentum().theta()) << " " << theSec->GetMomentum() << " E "<< theSec->GetKineticEnergy() << " " << theSec->GetDefinition()->GetParticleName() << G4endl; //GDEB
590  delete tmpHadrons->operator[](i);
591  }
592 #ifdef G4PHPDEBUG
593  if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N photons " << nPhotons << G4endl;
594 #endif
595  if(thePhotons != 0)
596  {
597  for(i=0; i<nPhotons; i++)
598  {
599  theSec = new G4DynamicParticle;
600  theSec->SetDefinition(thePhotons->operator[](i)->GetDefinition());
601  theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
602  theResult.Get()->AddSecondary(theSec);
603 #ifdef G4PHPDEBUG
604  if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply add secondary3 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
605 #endif
606  delete thePhotons->operator[](i);
607  }
608  }
609 
610 // some garbage collection
611  delete thePhotons;
612  delete tmpHadrons;
613 
614 //080721
616  G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
617  G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
618  G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
619  adjust_final_state ( init_4p_lab );
620 
621 // clean up the primary neutron
623 }
static G4He3 * He3()
Definition: G4He3.cc:94
G4double total(Particle const *const p1, Particle const *const p2)
static G4ParticleHPManager * GetInstance()
G4double Sample(G4double anEnergy, G4int &it)
void Init(std::istream &aDataFile)
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * theProjectile
G4ReactionProductVector * Sample(G4double anEnergy)
const G4String & GetParticleName() const
void SetTarget(G4ReactionProduct &aTarget)
G4double GetPDGCharge() const
void SetTarget(G4ReactionProduct *aTarget)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
const G4ParticleDefinition * GetDefinition() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ParticleHPPhotonDist * theFinalStatePhotons
G4double GetPDGMass() const
double theta() const
void SetTotalEnergy(const G4double en)
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
void SetTarget(const G4ReactionProduct &aTarget)
G4ReactionProductVector * GetDecayGammas(G4int aLevel)
Float_t Z
void GetDataStream(G4String, std::istringstream &iss)
G4double GetTotalEnergy() const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void SetKineticEnergy(const G4double en)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
G4double GetKineticEnergy() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &bit, G4ParticleDefinition *)
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4ParticleDefinition * GetDefinition() const
const G4Material * GetMaterial() const
G4double GetLevelEnergy(G4int aLevel)
void SetProjectileRP(G4ReactionProduct *aIncidentParticleRP)
void InitGammas(G4double AR, G4double ZR)
void SetProjectileRP(const G4ReactionProduct &anIncidentParticleRP)
void Init(G4double aMass, G4int aCount)
G4Cache< G4HadFinalState * > theResult
static constexpr double MeV
G4bool InitMean(std::istream &aDataFile)
#define G4UniformRand()
Definition: Randomize.hh:53
G4ErrorTarget * theTarget
Definition: errprop.cc:59
double A(double temperature)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static const G4int maxA
static G4double GetNuclearMass(const G4double A, const G4double Z)
double mag2() const
value_type & Get() const
Definition: G4Cache.hh:314
static constexpr double eV
G4ThreeVector GetMomentum() const
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
#define INT_MAX
Definition: templates.hh:111
void Put(const value_type &val) const
Definition: G4Cache.hh:318
const G4LorentzVector & Get4Momentum() const
G4int GetNumberOfSecondaries() const
G4ParticleHPEnAngCorrelation * theEnergyAngData
G4ReactionProductVector * GetPhotons(G4double anEnergy)
Hep3Vector unit() const
G4ParticleHPAngular * theAngularDistribution
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
void Init(std::istream &aDataFile)
void adjust_final_state(G4LorentzVector)
G4double GetMass() const
const G4ParticleDefinition * GetDefinition() const
double mag() const
void SetProjectileRP(G4ReactionProduct &aIncidentPart)
int G4int
Definition: G4Types.hh:78
void SampleAndUpdate(G4ReactionProduct &anIncidentParticle)
static G4Triton * Triton()
Definition: G4Triton.cc:95
ifstream in
Definition: comparison.C:7
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void InitPartials(std::istream &aDataFile)
void InitAngular(std::istream &aDataFile)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetKineticEnergy() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4GLOB_DLL std::ostream G4cout
void Init(std::istream &aDataFile)
Hep3Vector vect() const
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:504
const G4ParticleDefinition * GetParticleDefinition() const
static const G4double eps
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass)
static constexpr double keV
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:113
void BaseApply(const G4HadProjectile &theTrack, G4ParticleDefinition **theDefs, G4int nDef)
void SetStatusChange(G4HadFinalStateStatus aS)
G4double GetTemperature() const
Definition: G4Material.hh:183
void InitEnergies(std::istream &aDataFile)
G4ParticleHPEnergyDistribution * theEnergyDistribution