Geant4  10.02.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 91837 2015-08-07 07:27:08Z 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 
110  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
111  while(0==particlesFromProjectile)
112  {
113  do
114  {
115  x = 2*G4UniformRand() - 1;
116  y = 2*G4UniformRand() - 1;
117  }
118  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
119  while(x*x + y*y > 1);
120  impactParameter = std::sqrt(x*x+y*y)*(targetOuterRadius+projectileOuterRadius);
121  ++totalTries;
122  area = pi*(targetOuterRadius+projectileOuterRadius)*
123  (targetOuterRadius+projectileOuterRadius);
124  G4double projectileHorizon = impactParameter-targetOuterRadius;
125 
126  // Empirical boundary transparency.
127  G4double empirical = G4UniformRand();
128  if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
129 
130  // Calculate the number of nucleons involved in collision
131  // From projectile
132  aPrim.StartLoop();
133 
134  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
135  while((pNucleon = aPrim.GetNextNucleon()))
136  {
137  if(pNucleon->GetPosition().y()>projectileHorizon)
138  {
139  // We have one
140  ++particlesFromProjectile;
141  if(pNucleon->GetParticleType() == proton)
142  {
143  ++chargedFromProjectile;
144  }
145  }
146  }
147  }
148  ++hits;
149 
150  // From target:
151  G4double targetHorizon = impactParameter-projectileOuterRadius;
152  G4int chargedFromTarget = 0;
153  G4int particlesFromTarget = 0;
154  aTarg.StartLoop();
155  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
156  while((pNucleon = aTarg.GetNextNucleon()))
157  {
158  if(pNucleon->GetPosition().y()>targetHorizon)
159  {
160  // We have one
161  ++particlesFromTarget;
162  if(pNucleon->GetParticleType() == proton)
163  {
164  ++chargedFromTarget;
165  }
166  }
167  }
168 
169  // Energy sharing between projectile and target.
170  // Note that this is a quite simplistic kinetically.
171  G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
172  G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
173 
174  G4double projTotEnergy = thePrimary.GetTotalEnergy();
175  G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
176  G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
177 
178  // take the nucleons and fill the Fragments
179  G4Fragment anInitialState(aTargetA+particlesFromProjectile,
180  aTargetZ+chargedFromProjectile,
181  fragment4Momentum);
182  // M.A. Cortes fix
183  //anInitialState.SetNumberOfParticles(particlesFromProjectile);
184  anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile
185  + particlesFromTarget,
186  chargedFromProjectile
187  + chargedFromTarget);
188  anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
189  chargedFromProjectile + chargedFromTarget);
190  G4double time = thePrimary.GetGlobalTime();
191  anInitialState.SetCreationTime(time);
192 
193  // Fragment the Fragment using Pre-compound
194  G4ReactionProductVector* thePreCompoundResult =
195  theModel->DeExcite(anInitialState);
196 
197  // De-excite the projectile using ExcitationHandler
198  G4ReactionProductVector * theExcitationResult = 0;
199  if(particlesFromProjectile < aProjectileA)
200  {
201  G4LorentzVector residual4Momentum(momentum*(1.0-w), projTotEnergy*(1.0-w));
202 
203  G4Fragment initialState2(aProjectileA-particlesFromProjectile,
204  aProjectileZ-chargedFromProjectile,
205  residual4Momentum );
206 
207  // half of particles are excited (?!)
208  G4int pinit = (aProjectileA-particlesFromProjectile)/2;
209  G4int cinit = (aProjectileZ-chargedFromProjectile)/2;
210 
211  initialState2.SetNumberOfExcitedParticle(pinit,cinit);
212  initialState2.SetNumberOfHoles(pinit,cinit);
213  initialState2.SetCreationTime(time);
214 
215  theExcitationResult = theHandler->BreakItUp(initialState2);
216  }
217 
218  // Fill the particle change and clear intermediate vectors
219  G4int nexc = 0;
220  G4int npre = 0;
221  if(theExcitationResult) { nexc = theExcitationResult->size(); }
222  if(thePreCompoundResult) { npre = thePreCompoundResult->size();}
223 
224  if(nexc > 0) {
225  for(G4int k=0; k<nexc; ++k) {
226  G4ReactionProduct* p = (*theExcitationResult)[k];
228  p->GetMomentum()));
229  delete p;
230  }
231  }
232 
233  if(npre > 0) {
234  for(G4int k=0; k<npre; ++k) {
235  G4ReactionProduct* p = (*thePreCompoundResult)[k];
237  p->GetMomentum()));
238  delete p;
239  }
240  }
241 
242  delete thePreCompoundResult;
243  delete theExcitationResult;
244 
245  // return the particle change
246  return &theResult;
247 }
G4PreCompoundModel * theModel
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
const G4ParticleDefinition * GetParticleType() const
Definition: G4Nucleon.hh:84
static G4double GetNuclearMass(const G4double A, const G4double Z)
const G4ParticleDefinition * proton
CLHEP::Hep3Vector G4ThreeVector
const G4double w[NPOINTSGL]
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:357
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
int G4int
Definition: G4Types.hh:78
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:338
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:97
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:425
static const double pi
Definition: G4SIunits.hh:74
int G4lrint(double ad)
Definition: templates.hh:163
const G4double x[NPOINTSGL]
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ThreeVector GetMomentum() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4double GetOuterRadius()
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
G4Nucleon * GetNextNucleon()
G4ExcitationHandler * theHandler
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector