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