source: trunk/source/processes/hadronic/models/rpg/src/G4RPGProtonInelastic.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 10.0 KB
Line 
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// $Id: G4RPGProtonInelastic.cc,v 1.4 2008/05/05 21:21:55 dennis Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
28//
29 
30#include "G4RPGProtonInelastic.hh"
31#include "Randomize.hh"
32 
33G4HadFinalState*
34G4RPGProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
35                                    G4Nucleus& targetNucleus )
36{
37  theParticleChange.Clear();
38  const G4HadProjectile *originalIncident = &aTrack;
39  if (originalIncident->GetKineticEnergy()<= 0.1) 
40  {
41    theParticleChange.SetStatusChange(isAlive);
42    theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
43    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
44    return &theParticleChange;     
45  }
46
47  //
48  // create the target particle
49  //
50  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
51
52  if (originalIncident->GetKineticEnergy()/GeV < 0.01+2.*G4UniformRand()/9. )
53  {
54    SlowProton( originalIncident, targetNucleus );
55    delete originalTarget;
56    return &theParticleChange;
57  }
58
59  // Fermi motion and evaporation
60  // As of Geant3, the Fermi energy calculation had not been Done
61
62  G4double ek = originalIncident->GetKineticEnergy();
63  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
64  G4ReactionProduct modifiedOriginal;
65  modifiedOriginal = *originalIncident;
66   
67  G4double tkin = targetNucleus.Cinema( ek );
68  ek += tkin;
69  modifiedOriginal.SetKineticEnergy(ek);
70  G4double et = ek + amas;
71  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
72  G4double pp = modifiedOriginal.GetMomentum().mag();
73  if (pp > 0.0) {
74    G4ThreeVector momentum = modifiedOriginal.GetMomentum();
75    modifiedOriginal.SetMomentum( momentum * (p/pp) );
76  }
77  //
78  // calculate black track energies
79  //
80  tkin = targetNucleus.EvaporationEffects(ek);
81  ek -= tkin;
82  modifiedOriginal.SetKineticEnergy(ek);
83  et = ek + amas;
84  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
85  pp = modifiedOriginal.GetMomentum().mag();
86  if (pp > 0.0) {
87    G4ThreeVector momentum = modifiedOriginal.GetMomentum();
88    modifiedOriginal.SetMomentum( momentum * (p/pp) );
89  }
90  const G4double cutOff = 0.1;
91  if (modifiedOriginal.GetKineticEnergy() < cutOff) {
92    SlowProton( originalIncident, targetNucleus );
93    delete originalTarget;
94    return &theParticleChange;
95  }
96
97  G4ReactionProduct currentParticle = modifiedOriginal;
98  G4ReactionProduct targetParticle;
99  targetParticle = *originalTarget;
100  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
101  targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
102  G4bool incidentHasChanged = false;
103  G4bool targetHasChanged = false;
104  G4bool quasiElastic = false;
105  G4FastVector<G4ReactionProduct,256> vec;  // vec will contain the sec. particles
106  G4int vecLen = 0;
107  vec.Initialize( 0 );
108
109  InitialCollision(vec, vecLen, currentParticle, targetParticle,
110                   incidentHasChanged, targetHasChanged);
111
112  CalculateMomenta(vec, vecLen,
113                   originalIncident, originalTarget, modifiedOriginal,
114                   targetNucleus, currentParticle, targetParticle,
115                   incidentHasChanged, targetHasChanged, quasiElastic);
116
117  SetUpChange( vec, vecLen,
118               currentParticle, targetParticle,
119               incidentHasChanged );
120
121  delete originalTarget;
122  return &theParticleChange;
123}
124
125
126void
127G4RPGProtonInelastic::SlowProton(const G4HadProjectile *originalIncident,
128                                 G4Nucleus &targetNucleus )
129{
130  const G4double A = targetNucleus.GetN();    // atomic weight
131  const G4double Z = targetNucleus.GetZ();    // atomic number
132//  G4double currentKinetic = originalIncident->GetKineticEnergy();
133  //
134  // calculate Q-value of reactions
135  //
136  G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
137  G4double massVec[9];
138  massVec[0] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
139  massVec[1] = 0.;
140  if (A > Z+1.0)
141     massVec[1] = targetNucleus.AtomicMass( A    , Z+1.0 );
142  massVec[2] = theAtomicMass;
143  massVec[3] = 0.;
144  if (A > 1.0 && A-1.0 > Z) 
145     massVec[3] = targetNucleus.AtomicMass( A-1.0, Z );
146  massVec[4] = 0.;
147  if (A > 2.0 && A-2.0 > Z) 
148     massVec[4] = targetNucleus.AtomicMass( A-2.0, Z     );
149  massVec[5] = 0.;
150  if (A > 3.0 && Z > 1.0 && A-3.0 > Z-1.0) 
151     massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-1.0 );
152  massVec[6] = 0.;
153  if (A > 1.0 && A-1.0 > Z+1.0) 
154     massVec[6] = targetNucleus.AtomicMass( A-1.0, Z+1.0 );
155  massVec[7] = massVec[3];
156  massVec[8] = 0.;
157  if (A > 1.0 && Z > 1.0) 
158     massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
159   
160  G4FastVector<G4ReactionProduct,4> vec;  // vec will contain the secondary particles
161  G4int vecLen = 0;
162  vec.Initialize( 0 );
163   
164  twoBody.NuclearReaction( vec, vecLen, originalIncident,
165                           targetNucleus, theAtomicMass, massVec );
166   
167  theParticleChange.SetStatusChange( stopAndKill );
168  theParticleChange.SetEnergyChange( 0.0 );
169   
170  G4DynamicParticle *pd;
171  for( G4int i=0; i<vecLen; ++i )
172  {
173    pd = new G4DynamicParticle();
174    pd->SetDefinition( vec[i]->GetDefinition() );
175    pd->SetMomentum( vec[i]->GetMomentum() );
176    theParticleChange.AddSecondary( pd );
177    delete vec[i];
178  }
179}
180
181
182// Initial Collision
183//   selects the particle types arising from the initial collision of
184//   the proton and target nucleon.  Secondaries are assigned to forward
185//   and backward reaction hemispheres, but final state energies and
186//   momenta are not calculated here.
187
188void 
189G4RPGProtonInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
190                                  G4int& vecLen,
191                                  G4ReactionProduct& currentParticle,
192                                  G4ReactionProduct& targetParticle,
193                                  G4bool& incidentHasChanged,
194                                  G4bool& targetHasChanged)
195{
196  G4double KE = currentParticle.GetKineticEnergy()/GeV;
197
198  G4int mult;
199  G4int partType;
200  std::vector<G4int> fsTypes;
201  G4int part1; 
202  G4int part2;
203
204  G4double testCharge;
205  G4double testBaryon;
206  G4double testStrange;
207 
208  // Get particle types according to incident and target types
209 
210  if (targetParticle.GetDefinition() == particleDef[pro]) {
211    mult = GetMultiplicityT1(KE);
212    fsTypes = GetFSPartTypesForPP(mult, KE);
213
214    part1 = fsTypes[0];
215    part2 = fsTypes[1];
216    currentParticle.SetDefinition(particleDef[part1]);
217    targetParticle.SetDefinition(particleDef[part2]);
218    if (part1 == pro) {
219      if (part2 == neu) {
220        if (G4UniformRand() > 0.5) {
221          incidentHasChanged = true;
222          targetParticle.SetDefinition(particleDef[part1]);
223          currentParticle.SetDefinition(particleDef[part2]);
224        } else {
225          targetHasChanged = true;
226        }
227      } else if (part2 > neu && part2 < xi0) {
228        targetHasChanged = true;
229      }
230
231    } else {  // neutron
232      targetHasChanged = true;
233      incidentHasChanged = true;
234    }
235
236    testCharge = 2.0;
237    testBaryon = 2.0;
238    testStrange = 0.0;
239 
240  } else {   // target was a neutron
241    mult = GetMultiplicityT0(KE);
242    fsTypes = GetFSPartTypesForPN(mult, KE);
243
244    part1 = fsTypes[0];
245    part2 = fsTypes[1];
246    currentParticle.SetDefinition(particleDef[part1]);
247    targetParticle.SetDefinition(particleDef[part2]);
248    if (part1 == pro) {
249      if (part2 == pro) {
250        targetHasChanged = true;
251      } else if (part2 == neu) {
252        if (G4UniformRand() > 0.5) {
253          incidentHasChanged = true;
254          targetHasChanged = true;
255          targetParticle.SetDefinition(particleDef[part1]);
256          currentParticle.SetDefinition(particleDef[part2]);
257        }
258      } else {  // hyperon
259        targetHasChanged = true;
260      }
261
262    } else {  // neutron
263      incidentHasChanged = true;
264      if (part2 > neu && part2 < xi0) targetHasChanged = true;
265    }
266
267    testCharge = 1.0;
268    testBaryon = 2.0;
269    testStrange = 0.0;
270  }
271
272  // Remove incident and target from fsTypes
273 
274  fsTypes.erase(fsTypes.begin());
275  fsTypes.erase(fsTypes.begin());
276
277  // Remaining particles are secondaries.  Put them into vec.
278 
279  G4ReactionProduct* rp(0);
280  for(G4int i=0; i < mult-2; ++i ) {
281    partType = fsTypes[i];
282    rp = new G4ReactionProduct();
283    rp->SetDefinition(particleDef[partType]);
284    (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
285    vec.SetElement(vecLen++, rp);
286  }
287
288  // Check conservation of charge, strangeness, baryon number
289 
290  CheckQnums(vec, vecLen, currentParticle, targetParticle,
291             testCharge, testBaryon, testStrange);
292 
293  return;
294}
Note: See TracBrowser for help on using the repository browser.