Geant4_10
G4NeutronHPFinalState.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 //
27 //080721 Create adjust_final_state method by T. Koi
28 //080801 Residual reconstruction with theNDLDataA,Z (A, Z, and momentum are adjusted) by T. Koi
29 //101110 Set lower limit for gamma energy(1keV) by T. Koi
30 
31 #include "G4NeutronHPFinalState.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4IonTable.hh"
35 #include "G4Gamma.hh"
36 #include "G4Neutron.hh"
37 
39 {
40 
41  G4double minimum_energy = 1*keV;
42 
43  if ( adjustResult != true ) return;
44 
45  G4int nSecondaries = theResult.GetNumberOfSecondaries();
46 
47  G4int sum_Z = 0;
48  G4int sum_A = 0;
49  G4int max_SecZ = 0;
50  G4int max_SecA = 0;
51  G4int imaxA = -1;
52  for ( int i = 0 ; i < nSecondaries ; i++ )
53  {
55  max_SecZ = std::max ( max_SecZ , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
57  max_SecA = std::max ( max_SecA , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
58  if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
59  }
60 
61  G4ParticleDefinition* resi_pd = NULL;
62  G4bool needOneMoreSec = false;
63  G4ParticleDefinition* oneMoreSec_pd = NULL;
64  if ( (int)(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) == 0 )
65  {
66  //All secondaries are already created;
67  resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition();
68  }
69  else
70  {
71  if ( max_SecA > int(theBaseA + 1 - sum_A) )
72  {
73  //Most heavy secondary is interpreted as residual
74  resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition();
75  needOneMoreSec = true;
76  }
77  else
78  {
79  //creation of residual is requierd
80  resi_pd = G4IonTable::GetIonTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
81  }
82 
83  if ( needOneMoreSec )
84  {
85  if ( int(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) > 0 )
86  {
87  if ( int(theBaseA + 1 - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
88  oneMoreSec_pd = G4Neutron::Neutron();
89  }
90  else
91  {
92  oneMoreSec_pd = G4IonTable::GetIonTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
93  }
94  }
95 
96  if ( resi_pd == NULL )
97  {
98  // theNDLDataZ,A has the Z and A of used NDL file
99  if ( (int)(theNDLDataZ - sum_Z) == 0 && (int)(theNDLDataA + 1 - sum_A) == 0 )
100  {
101  G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
102  G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
103  resi_pd = G4IonTable::GetIonTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
104  for ( int i = 0 ; i < nSecondaries ; i++ )
105  {
106  if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ
107  && theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
108  {
110  p = p * resi_pd->GetPDGMass()/ G4IonTable::GetIonTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass();
111  theResult.GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
113  }
114  }
115  }
116  }
117  }
118 
119 
120  G4LorentzVector secs_4p_lab( 0.0 );
121 
123  G4double fast = 0;
124  G4double slow = 1;
125  G4int ifast = 0;
126  G4int islow = 0;
127  G4int ires = -1;
128 
129  for ( G4int i = 0 ; i < n_sec ; i++ )
130  {
131 
132  //G4cout << "HP_DB " << i
133  // << " " << theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName()
134  // << " 4p " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum()
135  // << " ke " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e() - theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetPDGMass()
136  // << G4endl;
137 
138  secs_4p_lab += theResult.GetSecondary( i )->GetParticle()->Get4Momentum();
139 
140  G4double beta = 0;
142  {
143  beta = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().beta();
144  }
145  else
146  {
147  beta = 1;
148  }
149 
150  if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i;
151 
152  if ( slow > beta && beta != 0 )
153  {
154  slow = beta;
155  islow = i;
156  }
157 
158  if ( fast <= beta )
159  {
160  if ( fast != 1 )
161  {
162  fast = beta;
163  ifast = i;
164  }
165  else
166  {
167 // fast is already photon then check E
169  if ( e > theResult.GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )
170  {
171 // among photons, the highest E becomes the fastest
172  ifast = i;
173  }
174  }
175  }
176  }
177 
178 
179  G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
180 
181  //G4cout << "HP_DB dif_4p " << init_4p_lab - secs_4p_lab << G4endl;
182  //G4cout << "HP_DB dif_3p mag " << ( dif_4p.v() ).mag() << G4endl;
183  //G4cout << "HP_DB dif_e " << dif_4p.e() - ( dif_4p.v() ).mag()<< G4endl;
184 
185  G4LorentzVector p4(0);
186  if ( ires == -1 )
187  {
188 // Create and Add Residual Nucleus
189  ires = nSecondaries;
190  nSecondaries += 1;
191 
192  G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );
193  theResult.AddSecondary ( res );
194 
195  p4 = res->Get4Momentum();
196  if ( slow > p4.beta() )
197  {
198  slow = p4.beta();
199  islow = ires;
200  }
201  dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
202  }
203 
204  if ( needOneMoreSec )
205  {
206  nSecondaries += 1;
207  G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );
208  theResult.AddSecondary ( one );
209  p4 = one->Get4Momentum();
210  if ( slow > p4.beta() )
211  {
212  slow = p4.beta();
213  islow = nSecondaries-1; //Because the first is 0th, so the last becomes "nSecondaries-1"
214  }
215  dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
216  }
217 
218  //Which is bigger dif_p or dif_e
219 
220  if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
221  {
222 
223  // Adjust p
224  //if ( dif_4p.v().mag() < 1*MeV )
225  if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
226  {
227 
228  nSecondaries += 1;
229  theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );
230 
231  }
232  else
233  {
234  //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
235  }
236 
237  }
238  else
239  {
240 
241  // dif_p > dif_e
242  // at first momentum
243  // Move residual momentum
244 
245  p4 = theResult.GetSecondary( ires )->GetParticle()->Get4Momentum();
246  theResult.GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() );
247  dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.GetSecondary( ires )->GetParticle()->Get4Momentum() );
248 
249  //G4cout << "HP_DB new residual kinetic energy " << theResult.GetSecondary( ires )->GetParticle()->GetKineticEnergy() << G4endl;
250 
251  }
252 
253  G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
254  //G4cout << "HP_DB dif_e " << dif_e << G4endl;
255 
256  if ( dif_e > 0 )
257  {
258 
259 // create 2 gamma
260 
261  nSecondaries += 2;
262  G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
263 
264  if ( minimum_energy < e1 )
265  {
266  G4double costh = 2.*G4UniformRand()-1.;
267  G4double phi = twopi*G4UniformRand();
268  G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi),
269  std::sin(std::acos(costh))*std::sin(phi),
270  costh);
271  theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );
272  theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );
273  }
274  else
275  {
276  //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
277  }
278 
279  }
280  else //dif_e < 0
281  {
282 
283 // At first reduce KE of the fastest secondary;
286  G4ThreeVector dir = ( theResult.GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
287 
288  //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
289 
290  if ( ke0 + dif_e > 0 )
291  {
292  theResult.GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e );
293  G4ThreeVector dp = p0 - theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
294 
296  //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir );
297  theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p + dp );
298  }
299  else
300  {
301  //G4cout << "HP_DB Difference in dif_e too large ( <0MeV ) to adjust, so that give up tuning" << G4endl;
302  }
303 
304  }
305 
306 }
void SetMomentum(const G4ThreeVector &momentum)
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
Hep3Vector v() const
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
G4int GetAtomicNumber() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4int GetAtomicMass() const
void SetKineticEnergy(G4double aEnergy)
G4LorentzVector Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:80
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4DynamicParticle * GetParticle()
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void adjust_final_state(G4LorentzVector)
double mag() const
G4int GetNumberOfSecondaries() const
void AddSecondary(G4DynamicParticle *aP)
TDirectory * dir
Definition: macro.C:5
G4ThreeVector GetMomentum() const