Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPFinalState.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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 //
32 
34 
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4IonTable.hh"
38 #include "G4Gamma.hh"
39 #include "G4Neutron.hh"
40 #include "G4Proton.hh"
41 #include "G4Deuteron.hh"
42 #include "G4Triton.hh"
43 #include "G4He3.hh"
44 #include "G4Alpha.hh"
45 
47 {
48 
49  G4double minimum_energy = 1*keV;
50 
51  if ( G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState() ) return;
52 
53  G4int nSecondaries = theResult.Get()->GetNumberOfSecondaries();
54 
55  G4int sum_Z = 0;
56  G4int sum_A = 0;
57  G4int max_SecZ = 0;
58  G4int max_SecA = 0;
59  G4int imaxA = -1;
60  for ( int i = 0 ; i < nSecondaries ; i++ )
61  {
62  //G4cout << "G4ParticleHPFinalState::adjust_final_state theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() = " << theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() << G4endl;
64  max_SecZ = std::max ( max_SecZ , theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
66  max_SecA = std::max ( max_SecA , theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
67  if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
68 #ifdef G4PHPDEBUG
69  if( getenv("G4ParticleHPDebug")) G4cout << "G4ParticleHPFinalState::adjust_final_stat SECO " << i << " " <<theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() << G4endl;
70 #endif
71 
72  }
73 
74  G4ParticleDefinition* resi_pd = NULL;
75 
76  G4double baseZNew = theBaseZ;
77  G4double baseANew = theBaseA;
79  baseANew ++;
80  } else if( theProjectile == G4Proton::Proton() ) {
81  baseZNew ++;
82  baseANew ++;
83  } else if( theProjectile == G4Deuteron::Deuteron() ) {
84  baseZNew ++;
85  baseANew += 2;
86  } else if( theProjectile == G4Triton::Triton() ) {
87  baseZNew ++;
88  baseANew += 3;
89  } else if( theProjectile == G4He3::He3() ) {
90  baseZNew += 2;
91  baseANew += 3;
92  } else if( theProjectile == G4Alpha::Alpha() ) {
93  baseZNew += 2;
94  baseANew += 4;
95  }
96 
97 #ifdef G4PHPDEBUG
98  if( getenv("G4ParticleHPDebug")) G4cout << "G4ParticleHPFinalState::adjust_final_stat BaseZ " << baseZNew << " BaseA " << baseANew << " sum_Z " << sum_Z << " sum_A " << sum_A << G4endl;
99 #endif
100 
101  G4bool needOneMoreSec = false;
102  G4ParticleDefinition* oneMoreSec_pd = NULL;
103  if ( (int)(baseZNew - sum_Z) == 0 && (int)(baseANew - sum_A) == 0 )
104  {
105  //All secondaries are already created;
106  resi_pd = theResult.Get()->GetSecondary( imaxA )->GetParticle()->GetDefinition();
107  }
108  else
109  {
110  if ( max_SecA >= int(baseANew - sum_A) )
111  {
112  //Most heavy secondary is interpreted as residual
113  resi_pd = theResult.Get()->GetSecondary( imaxA )->GetParticle()->GetDefinition();
114  needOneMoreSec = true;
115  }
116  else
117  {
118  //creation of residual is requierd
119  resi_pd = G4IonTable::GetIonTable()->GetIon ( int(baseZNew - sum_Z) , (int)(baseANew - sum_A) , 0.0 );
120  }
121 
122  if ( needOneMoreSec )
123  {
124  if ( int(baseZNew - sum_Z) == 0 && (int)(baseANew - sum_A) > 0 )
125  {
126  //In this case, one neutron is added to secondaries
127  if ( int(baseANew - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
128  oneMoreSec_pd = G4Neutron::Neutron();
129  }
130  else
131  {
132 #ifdef G4PHPDEBUG
133  if( getenv("G4ParticleHPDebug")) G4cout << this << "G4ParticleHPFinalState oneMoreSec_pd Z " << baseZNew << " - " << sum_Z << " A " << baseANew << " - " << sum_A << " projectile " << theProjectile->GetParticleName() << G4endl;
134 #endif
135  oneMoreSec_pd = G4IonTable::GetIonTable()->GetIon ( int(baseZNew - sum_Z) , (int)(baseANew - sum_A) , 0.0 );
136  if( !oneMoreSec_pd ) {
137  G4cerr << this << "G4ParticleHPFinalState oneMoreSec_pd Z " << baseZNew << " - " << sum_Z << " A " << baseANew << " - " << sum_A << " projectile " << theProjectile->GetParticleName() << G4endl;
138  G4Exception("G4ParticleHPFinalState:adjust_final_state",
139  "Warning",
140  JustWarning,
141  "No adjustment will be done!");
142  return;
143  }
144  }
145  }
146 
147  if ( resi_pd == NULL )
148  {
149  // theNDLDataZ,A has the Z and A of used NDL file
150  G4double ndlZNew = theNDLDataZ;
151  G4double ndlANew = theNDLDataA;
152  if( theProjectile == G4Neutron::Neutron() ) {
153  ndlANew ++;
154  } else if( theProjectile == G4Proton::Proton() ) {
155  ndlZNew ++;
156  ndlANew ++;
157  } else if( theProjectile == G4Deuteron::Deuteron() ) {
158  ndlZNew ++;
159  ndlANew += 2;
160  } else if( theProjectile == G4Triton::Triton() ) {
161  ndlZNew ++;
162  ndlANew += 3;
163  } else if( theProjectile == G4He3::He3() ) {
164  ndlZNew += 2;
165  ndlANew += 3;
166  } else if( theProjectile == G4Alpha::Alpha() ) {
167  ndlZNew += 2;
168  ndlANew += 4;
169  }
170  // theNDLDataZ,A has the Z and A of used NDL file
171  if ( (int)(ndlZNew - sum_Z) == 0 && (int)(ndlANew - sum_A) == 0 )
172  {
173  G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
174  G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
175  resi_pd = G4IonTable::GetIonTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
176  if( !resi_pd ) {
177  G4cerr << "G4ParticleHPFinalState resi_pd Z " << max_SecZ << " - " << dif_Z << " A " << max_SecA << " - " << dif_A << " projectile " << theProjectile->GetParticleName() << G4endl;
178  G4Exception("G4ParticleHPFinalState:adjust_final_state",
179  "Warning",
180  JustWarning,
181  "No adjustment will be done!");
182  return;
183  }
184 
185  for ( int i = 0 ; i < nSecondaries ; i++ )
186  {
187  if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ
188  && theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
189  {
191  p = p * resi_pd->GetPDGMass()/ G4IonTable::GetIonTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass();
192  theResult.Get()->GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
194  }
195  }
196  }
197  }
198  }
199 
200 
201  G4LorentzVector secs_4p_lab( 0.0 );
202 
204  G4double fast = 0;
205  G4double slow = 1;
206  G4int ifast = 0;
207  G4int islow = 0;
208  G4int ires = -1;
209 
210  for ( G4int i = 0 ; i < n_sec ; i++ )
211  {
212 
213  //G4cout << "HP_DB " << i
214  // << " " << theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName()
215  // << " 4p " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum()
216  // << " ke " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e() - theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetPDGMass()
217  // << G4endl;
218 
219  secs_4p_lab += theResult.Get()->GetSecondary( i )->GetParticle()->Get4Momentum();
220 
221  G4double beta = 0;
223  {
224  beta = theResult.Get()->GetSecondary( i )->GetParticle()->Get4Momentum().beta();
225  }
226  else
227  {
228  beta = 1;
229  }
230 
231  if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i;
232 
233  if ( slow > beta && beta != 0 )
234  {
235  slow = beta;
236  islow = i;
237  }
238 
239  if ( fast <= beta )
240  {
241  if ( fast != 1 )
242  {
243  fast = beta;
244  ifast = i;
245  }
246  else
247  {
248 // fast is already photon then check E
250  if ( e > theResult.Get()->GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )
251  {
252 // among photons, the highest E becomes the fastest
253  ifast = i;
254  }
255  }
256  }
257  }
258 
259 
260  G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
261 
262  //G4cout << "HP_DB dif_4p " << init_4p_lab - secs_4p_lab << G4endl;
263  //G4cout << "HP_DB dif_3p mag " << ( dif_4p.v() ).mag() << G4endl;
264  //G4cout << "HP_DB dif_e " << dif_4p.e() - ( dif_4p.v() ).mag()<< G4endl;
265 
266  G4LorentzVector p4(0);
267  if ( ires == -1 )
268  {
269 // Create and Add Residual Nucleus
270  ires = nSecondaries;
271  nSecondaries += 1;
272 
273  G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );
274  theResult.Get()->AddSecondary ( res );
275 
276  p4 = res->Get4Momentum();
277  if ( slow > p4.beta() )
278  {
279  slow = p4.beta();
280  islow = ires;
281  }
282  dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
283  }
284 
285  if ( needOneMoreSec && oneMoreSec_pd)
286  //
287  // fca: this is not a fix, this is a crash avoidance...
288  // fca: the baryon number is still wrong, most probably because it
289  // fca: should have been decreased, but since we could not create a particle
290  // fca: we just do not add it
291  //
292  {
293  nSecondaries += 1;
294  G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );
295  theResult.Get()->AddSecondary ( one );
296  p4 = one->Get4Momentum();
297  if ( slow > p4.beta() )
298  {
299  slow = p4.beta();
300  islow = nSecondaries-1; //Because the first is 0th, so the last becomes "nSecondaries-1"
301  }
302  dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
303  }
304 
305  //Which is bigger dif_p or dif_e
306 
307  if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
308  {
309 
310  // Adjust p
311  //if ( dif_4p.v().mag() < 1*MeV )
312  if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
313  {
314 
315  nSecondaries += 1;
316  theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );
317 
318  }
319  else
320  {
321  //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
322  }
323 
324  }
325  else
326  {
327 
328  // dif_p > dif_e
329  // at first momentum
330  // Move residual momentum
331 
332  p4 = theResult.Get()->GetSecondary( ires )->GetParticle()->Get4Momentum();
333  theResult.Get()->GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() );
334  dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.Get()->GetSecondary( ires )->GetParticle()->Get4Momentum() );
335 
336  //G4cout << "HP_DB new residual kinetic energy " << theResult.GetSecondary( ires )->GetParticle()->GetKineticEnergy() << G4endl;
337 
338  }
339 
340  G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
341  //G4cout << "HP_DB dif_e " << dif_e << G4endl;
342 
343  if ( dif_e > 0 )
344  {
345 
346 // create 2 gamma
347 
348  nSecondaries += 2;
349  G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
350 
351  if ( minimum_energy < e1 )
352  {
353  G4double costh = 2.*G4UniformRand()-1.;
354  G4double phi = twopi*G4UniformRand();
355  G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi),
356  std::sin(std::acos(costh))*std::sin(phi),
357  costh);
358  theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );
359  theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );
360  }
361  else
362  {
363  //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
364  }
365 
366  }
367  else //dif_e < 0
368  {
369 
370 // At first reduce KE of the fastest secondary;
373  G4ThreeVector dir = ( theResult.Get()->GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
374 
375  //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
376 
377  if ( ke0 + dif_e > 0 )
378  {
379  theResult.Get()->GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e );
380  G4ThreeVector dp = p0 - theResult.Get()->GetSecondary( ifast )->GetParticle()->GetMomentum();
381 
383  //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir );
384  theResult.Get()->GetSecondary( islow )->GetParticle()->SetMomentum( p + dp );
385  }
386  else
387  {
388  //G4cout << "HP_DB Difference in dif_e too large ( <0MeV ) to adjust, so that give up tuning" << G4endl;
389  }
390 
391  }
392 
393 }
394 
static G4ParticleHPManager * GetInstance()
G4Cache< G4HadFinalState * > theResult
void SetMomentum(const G4ThreeVector &momentum)
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
value_type & Get() const
Definition: G4Cache.hh:282
Hep3Vector v() const
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
G4ParticleDefinition * theProjectile
static constexpr double twopi
Definition: G4SIunits.hh:76
void adjust_final_state(G4LorentzVector)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4int GetAtomicMass() const
void SetKineticEnergy(G4double aEnergy)
G4LorentzVector Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
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
static constexpr double MeV
Definition: G4SIunits.hh:214
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
G4int GetNumberOfSecondaries() const
G4ThreeVector GetMomentum() const
G4GLOB_DLL std::ostream G4cerr