Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPProduct.cc
Go to the documentation of this file.
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 // 080718 As for secondary photons, if its mean value has a value of integer,
31 // then a sampling of multiplicity that based on Poisson Distribution
32 // is not carried out and the mean is used as a multiplicity.
33 // modified by T. Koi.
34 // 080721 Using ClearHistories() methodl for limiting the sum of secondary energies
35 // modified by T. Koi.
36 // 080901 bug fix of too many secnodaries production in nd reactinos by T. Koi
37 //
38 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
39 //
40 #include "G4ParticleHPProduct.hh"
41 #include "G4Poisson.hh"
42 #include "G4Proton.hh"
43 
45 {
46  //if(theDist == 0) { return 0; }
47  //151120 TK Modified for solving reproducibility problem
48  if ( theDist == 0 ) {
49  fCache.Get().theCurrentMultiplicity = 0;
50  return 0;
51  }
52 
53  G4double mean = theYield.GetY(anEnergy);
54  //g G4cout << "G4ParticleHPProduct MEAN NUMBER OF PARTICLES " << mean << " for " << theMass << G4endl;
55  //if( mean <= 0. ) return 0;
56  //151120 TK Modified for solving reproducibility problem
57  //This is also a real fix
58  if ( mean <= 0. ) {
59  fCache.Get().theCurrentMultiplicity = 0;
60  return 0;
61  }
62 
63  G4int multi;
64  multi = G4int(mean+0.0001);
65  //if(theMassCode==0) multi = G4Poisson(mean); // @@@@gammas. please X-check this
66  //080718
67 #ifdef PHP_AS_HP
68  if ( theMassCode == 0 ) // DELETE THIS: IT MUST BE DONE FOR ALL PARTICLES
69 #endif
70  {
71  if ( G4int ( mean ) == mean )
72  {
73  multi = (G4int) mean;
74  }
75  else
76  {
77 #ifdef PHP_AS_HP
78  multi = G4Poisson ( mean );
79 #else
80  if( theMultiplicityMethod == G4HPMultiPoisson ) {
81  multi = G4Poisson ( mean );
82  if( getenv("G4PHPTEST") ) G4cout << " MULTIPLICITY MULTIPLIED " << multi << " " << theMassCode << G4endl;
83  } else { // if( theMultiplicityMethod == G4HPMultiBetweenInts ) {
85  G4int imulti = G4int(mean);
86  multi = imulti + G4int(radnf < mean-imulti);
87  // G4cout << theMass << " multi " << multi << " mean " << mean
88  // << " radnf " << radnf << " mean-imulti " << mean-imulti << G4endl;
89  }
90 #endif
91  // multi = int(mean);
92  // if( CLHEP::RandFlat::shoot() > mean-multi ) multi++;
93  }
94 #ifdef G4PHPDEBUG
95  if( getenv("G4ParticleHPDebug") ) G4cout << "G4ParticleHPProduct::GetMultiplicity " << theMassCode << " " << theMass << " multi " << multi << " mean " << mean << G4endl;
96 #endif
97  }
98 
99  fCache.Get().theCurrentMultiplicity = static_cast<G4int>(mean);
100 
101  return multi;
102 }
103 
104 
106 {
107  if(theDist == 0) { return 0; }
109 
110  theDist->SetTarget(fCache.Get().theTarget);
111  theDist->SetProjectileRP(fCache.Get().theProjectileRP);
112  G4int i;
113 // G4double eMax = GetTarget()->GetMass()+GetNeutron()->GetMass()
114 // - theActualStateQValue;
115  G4ReactionProduct * tmp;
116  theDist->ClearHistories();
117 
118  for(i=0;i<multi;i++)
119  {
120 #ifdef G4PHPDEBUG
121  if( getenv("G4PHPTEST") )
122  if( getenv("G4ParticleHPDebug") && tmp != 0 ) G4cout << multi << " " << i << " @@@ G4ParticleHPProduct::Sample " << anEnergy << " Mass " << theMassCode << " " << theMass << G4endl;
123 #endif
124  tmp = theDist->Sample(anEnergy, theMassCode, theMass);
125  if(tmp != 0) { result->push_back(tmp); }
126 #ifndef G4PHPDEBUG //GDEB
127  if( getenv("G4ParticleHPDebug") && tmp != 0 ) G4cout << multi << " " << i << " @@@ G4ParticleHPProduct::Sample " << tmp->GetDefinition()->GetParticleName() << " E= " << tmp->GetKineticEnergy() << G4endl;
128 #endif
129  }
130  if(multi == 0)
131  {
132  tmp = theDist->Sample(anEnergy, theMassCode, theMass);
133  delete tmp;
134  }
135  /*
136  //080901 TK Comment out, too many secondaries are produced in deuteron reactions
137  if(theTarget->GetMass()<2*GeV) // @@@ take care of residuals in all cases
138  {
139  tmp = theDist->Sample(anEnergy, theMassCode, theMass);
140  tmp->SetDefinition(G4Proton::Proton());
141  if(tmp != 0) { result->push_back(tmp); }
142  }
143  */
144 
145  return result;
146 }
G4double G4ParticleHPJENDLHEData::G4double result
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
value_type & Get() const
Definition: G4Cache.hh:282
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
std::vector< G4ReactionProduct * > G4ReactionProductVector
static double shoot()
Definition: RandFlat.cc:59
const G4ParticleDefinition * GetDefinition() const
G4GLOB_DLL std::ostream G4cout
G4double GetY(G4double x)
void SetProjectileRP(G4ReactionProduct *aIncidentParticleRP)
G4double GetKineticEnergy() const
G4int GetMultiplicity(G4double anEnergy)
void SetTarget(G4ReactionProduct *aTarget)
virtual G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass)=0
G4ReactionProductVector * Sample(G4double anEnergy, G4int nParticles)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76