Geant4  10.00.p02
G4LowEIonFragmentation.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 // $Id: G4LowEIonFragmentation.cc 68717 2013-04-05 09:15:15Z gcosmo $
27 //
28 //---------------------------------------------------------------------------
29 //
30 // ClassName: G4LowEIonFragmentation
31 //
32 // Author: H.P. Wellisch
33 //
34 // Modified:
35 // 02 Jun 2010 M. A. Cortes Giraldo fix: particlesFromTarget must be
36 // accounted for as particles of initial compound nucleus
37 // 28 Oct 2010 V.Ivanchenko complete migration to integer Z and A;
38 // use updated G4Fragment methods
39 
40 #include <algorithm>
41 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Fancy3DNucleus.hh"
46 #include "G4Proton.hh"
47 #include "G4NucleiProperties.hh"
48 
50 {
51  theHandler = value;
54  hits = 0;
55  totalTries = 1;
56  area = 0.0;
57 }
58 
60 {
64  hits = 0;
65  totalTries = 1;
66  area = 0.0;
67 }
68 
70 {
71  delete theModel;
72 }
73 
75 ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus & theNucleus)
76 {
77  area = 0.0;
78  // initialize the particle change
79  theResult.Clear();
82 
83  // Get Target A, Z
84  G4int aTargetA = theNucleus.GetA_asInt();
85  G4int aTargetZ = theNucleus.GetZ_asInt();
86 
87  // Get Projectile A, Z
88  G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
89  G4int aProjectileZ = G4lrint(thePrimary.GetDefinition()->GetPDGCharge()/eplus);
90 
91  // Get Maximum radius of both
92 
93  G4Fancy3DNucleus aPrim;
94  aPrim.Init(aProjectileA, aProjectileZ);
95  G4double projectileOuterRadius = aPrim.GetOuterRadius();
96 
97  G4Fancy3DNucleus aTarg;
98  aTarg.Init(aTargetA, aTargetZ);
99  G4double targetOuterRadius = aTarg.GetOuterRadius();
100 
101  // Get the Impact parameter
102  G4int particlesFromProjectile = 0;
103  G4int chargedFromProjectile = 0;
104  G4double impactParameter = 0;
105  G4double x,y;
106  G4Nucleon * pNucleon;
107  // need at lease one particle from the projectile model beyond the
108  // projectileHorizon.
109  while(0==particlesFromProjectile)
110  {
111  do
112  {
113  x = 2*G4UniformRand() - 1;
114  y = 2*G4UniformRand() - 1;
115  }
116  while(x*x + y*y > 1);
117  impactParameter = std::sqrt(x*x+y*y)*(targetOuterRadius+projectileOuterRadius);
118  ++totalTries;
119  area = pi*(targetOuterRadius+projectileOuterRadius)*
120  (targetOuterRadius+projectileOuterRadius);
121  G4double projectileHorizon = impactParameter-targetOuterRadius;
122 
123  // Empirical boundary transparency.
124  G4double empirical = G4UniformRand();
125  if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
126 
127  // Calculate the number of nucleons involved in collision
128  // From projectile
129  aPrim.StartLoop();
130  while((pNucleon = aPrim.GetNextNucleon()))
131  {
132  if(pNucleon->GetPosition().y()>projectileHorizon)
133  {
134  // We have one
135  ++particlesFromProjectile;
136  if(pNucleon->GetParticleType() == proton)
137  {
138  ++chargedFromProjectile;
139  }
140  }
141  }
142  }
143  ++hits;
144 
145  // From target:
146  G4double targetHorizon = impactParameter-projectileOuterRadius;
147  G4int chargedFromTarget = 0;
148  G4int particlesFromTarget = 0;
149  aTarg.StartLoop();
150  while((pNucleon = aTarg.GetNextNucleon()))
151  {
152  if(pNucleon->GetPosition().y()>targetHorizon)
153  {
154  // We have one
155  ++particlesFromTarget;
156  if(pNucleon->GetParticleType() == proton)
157  {
158  ++chargedFromTarget;
159  }
160  }
161  }
162 
163  // Energy sharing between projectile and target.
164  // Note that this is a quite simplistic kinetically.
165  G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
166  G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
167 
168  G4double projTotEnergy = thePrimary.GetTotalEnergy();
169  G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
170  G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
171 
172  // take the nucleons and fill the Fragments
173  G4Fragment anInitialState(aTargetA+particlesFromProjectile,
174  aTargetZ+chargedFromProjectile,
175  fragment4Momentum);
176  // M.A. Cortes fix
177  //anInitialState.SetNumberOfParticles(particlesFromProjectile);
178  anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile+particlesFromTarget,
179  chargedFromProjectile + chargedFromTarget);
180  anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
181  chargedFromProjectile + chargedFromTarget);
182  G4double time = thePrimary.GetGlobalTime();
183  anInitialState.SetCreationTime(time);
184 
185  // Fragment the Fragment using Pre-compound
186  G4ReactionProductVector* thePreCompoundResult = theModel->DeExcite(anInitialState);
187 
188  // De-excite the projectile using ExcitationHandler
189  G4ReactionProductVector * theExcitationResult = 0;
190  if(particlesFromProjectile < aProjectileA)
191  {
192  G4LorentzVector residual4Momentum(momentum*(1.0-w), projTotEnergy*(1.0-w));
193 
194  G4Fragment initialState2(aProjectileA-particlesFromProjectile,
195  aProjectileZ-chargedFromProjectile,
196  residual4Momentum );
197 
198  // half of particles are excited (?!)
199  G4int pinit = (aProjectileA-particlesFromProjectile)/2;
200  G4int cinit = (aProjectileZ-chargedFromProjectile)/2;
201 
202  initialState2.SetNumberOfExcitedParticle(pinit,cinit);
203  initialState2.SetNumberOfHoles(pinit,cinit);
204  initialState2.SetCreationTime(time);
205 
206  theExcitationResult = theHandler->BreakItUp(initialState2);
207  }
208 
209  // Fill the particle change and clear intermediate vectors
210  G4int nexc = 0;
211  G4int npre = 0;
212  if(theExcitationResult) { nexc = theExcitationResult->size(); }
213  if(thePreCompoundResult) { npre = thePreCompoundResult->size();}
214 
215  if(nexc > 0) {
216  for(G4int k=0; k<nexc; ++k) {
217  G4ReactionProduct* p = (*theExcitationResult)[k];
219  delete p;
220  }
221  }
222 
223  if(npre > 0) {
224  for(G4int k=0; k<npre; ++k) {
225  G4ReactionProduct* p = (*thePreCompoundResult)[k];
227  delete p;
228  }
229  }
230 
231  delete thePreCompoundResult;
232  delete theExcitationResult;
233 
234  // return the particle change
235  return &theResult;
236 
237 }
G4PreCompoundModel * theModel
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static G4double GetNuclearMass(const G4double A, const G4double Z)
const G4ParticleDefinition * proton
CLHEP::Hep3Vector G4ThreeVector
const G4double pi
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:355
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
int G4int
Definition: G4Types.hh:78
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const
G4ParticleDefinition * GetDefinition() const
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:336
#define G4UniformRand()
Definition: Randomize.hh:87
const G4ParticleDefinition * GetDefinition() const
void Init(G4int theA, G4int theZ)
G4double GetGlobalTime() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
const G4LorentzVector & Get4Momentum() const
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)
void SetEnergyChange(G4double anEnergy)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:403
int G4lrint(double ad)
Definition: templates.hh:163
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ThreeVector GetMomentum() const
G4ParticleDefinition * GetParticleType() const
Definition: G4Nucleon.hh:84
G4double GetOuterRadius()
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
G4Nucleon * GetNextNucleon()
void AddSecondary(G4DynamicParticle *aP)
G4ExcitationHandler * theHandler
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector