ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEIonFragmentation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LowEIonFragmentation.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 //---------------------------------------------------------------------------
28 //
29 // ClassName: G4LowEIonFragmentation
30 //
31 // Author: H.P. Wellisch
32 //
33 // Modified:
34 // 02 Jun 2010 M. A. Cortes Giraldo fix: particlesFromTarget must be
35 // accounted for as particles of initial compound nucleus
36 // 28 Oct 2010 V.Ivanchenko complete migration to integer Z and A;
37 // use updated G4Fragment methods
38 
39 #include <algorithm>
40 
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Fancy3DNucleus.hh"
45 #include "G4Proton.hh"
46 #include "G4NucleiProperties.hh"
47 
49 {
50  theHandler = value;
53  hits = 0;
54  totalTries = 1;
55  area = 0.0;
56 }
57 
59 {
63  hits = 0;
64  totalTries = 1;
65  area = 0.0;
66 }
67 
69 {
70  delete theModel;
71 }
72 
74 ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus & theNucleus)
75 {
76  area = 0.0;
77  // initialize the particle change
78  theResult.Clear();
81 
82  // Get Target A, Z
83  G4int aTargetA = theNucleus.GetA_asInt();
84  G4int aTargetZ = theNucleus.GetZ_asInt();
85 
86  // Get Projectile A, Z
87  G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
88  G4int aProjectileZ = G4lrint(thePrimary.GetDefinition()->GetPDGCharge()/eplus);
89 
90  // Get Maximum radius of both
91 
92  G4Fancy3DNucleus aPrim;
93  aPrim.Init(aProjectileA, aProjectileZ);
94  G4double projectileOuterRadius = aPrim.GetOuterRadius();
95 
96  G4Fancy3DNucleus aTarg;
97  aTarg.Init(aTargetA, aTargetZ);
98  G4double targetOuterRadius = aTarg.GetOuterRadius();
99 
100  // Get the Impact parameter
101  G4int particlesFromProjectile = 0;
102  G4int chargedFromProjectile = 0;
103  G4double impactParameter = 0;
104  G4double x,y;
105  G4Nucleon * pNucleon;
106  // need at lease one particle from the projectile model beyond the
107  // projectileHorizon.
108 
109  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
110  while(0==particlesFromProjectile)
111  {
112  do
113  {
114  x = 2*G4UniformRand() - 1;
115  y = 2*G4UniformRand() - 1;
116  }
117  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
118  while(x*x + y*y > 1);
119  impactParameter = std::sqrt(x*x+y*y)*(targetOuterRadius+projectileOuterRadius);
120  ++totalTries;
121  area = pi*(targetOuterRadius+projectileOuterRadius)*
122  (targetOuterRadius+projectileOuterRadius);
123  G4double projectileHorizon = impactParameter-targetOuterRadius;
124 
125  // Empirical boundary transparency.
126  G4double empirical = G4UniformRand();
127  if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
128 
129  // Calculate the number of nucleons involved in collision
130  // From projectile
131  aPrim.StartLoop();
132 
133  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
134  while((pNucleon = aPrim.GetNextNucleon()))
135  {
136  if(pNucleon->GetPosition().y()>projectileHorizon)
137  {
138  // We have one
139  ++particlesFromProjectile;
140  if(pNucleon->GetParticleType() == proton)
141  {
142  ++chargedFromProjectile;
143  }
144  }
145  }
146  }
147  ++hits;
148 
149  // From target:
150  G4double targetHorizon = impactParameter-projectileOuterRadius;
151  G4int chargedFromTarget = 0;
152  G4int particlesFromTarget = 0;
153  aTarg.StartLoop();
154  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
155  while((pNucleon = aTarg.GetNextNucleon()))
156  {
157  if(pNucleon->GetPosition().y()>targetHorizon)
158  {
159  // We have one
160  ++particlesFromTarget;
161  if(pNucleon->GetParticleType() == proton)
162  {
163  ++chargedFromTarget;
164  }
165  }
166  }
167 
168  // Energy sharing between projectile and target.
169  // Note that this is a quite simplistic kinetically.
170  G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
171  G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
172 
173  G4double projTotEnergy = thePrimary.GetTotalEnergy();
174  G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
175  G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
176 
177  // take the nucleons and fill the Fragments
178  G4Fragment anInitialState(aTargetA+particlesFromProjectile,
179  aTargetZ+chargedFromProjectile,
180  fragment4Momentum);
181  // M.A. Cortes fix
182  //anInitialState.SetNumberOfParticles(particlesFromProjectile);
183  anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile
184  + particlesFromTarget,
185  chargedFromProjectile
186  + chargedFromTarget);
187  anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
188  chargedFromProjectile + chargedFromTarget);
189  G4double time = thePrimary.GetGlobalTime();
190  anInitialState.SetCreationTime(time);
191 
192  // Fragment the Fragment using Pre-compound
193  G4ReactionProductVector* thePreCompoundResult =
194  theModel->DeExcite(anInitialState);
195 
196  // De-excite the projectile using ExcitationHandler
197  G4ReactionProductVector * theExcitationResult = 0;
198  if(particlesFromProjectile < aProjectileA)
199  {
200  G4LorentzVector residual4Momentum(momentum*(1.0-w), projTotEnergy*(1.0-w));
201 
202  G4Fragment initialState2(aProjectileA-particlesFromProjectile,
203  aProjectileZ-chargedFromProjectile,
204  residual4Momentum );
205 
206  // half of particles are excited (?!)
207  G4int pinit = (aProjectileA-particlesFromProjectile)/2;
208  G4int cinit = (aProjectileZ-chargedFromProjectile)/2;
209 
210  initialState2.SetNumberOfExcitedParticle(pinit,cinit);
211  initialState2.SetNumberOfHoles(pinit,cinit);
212  initialState2.SetCreationTime(time);
213 
214  theExcitationResult = theHandler->BreakItUp(initialState2);
215  }
216 
217  // Fill the particle change and clear intermediate vectors
218  G4int nexc = 0;
219  G4int npre = 0;
220  if(theExcitationResult) { nexc = theExcitationResult->size(); }
221  if(thePreCompoundResult) { npre = thePreCompoundResult->size();}
222 
223  if(nexc > 0) {
224  for(G4int k=0; k<nexc; ++k) {
225  G4ReactionProduct* p = (*theExcitationResult)[k];
227  p->GetMomentum()));
228  delete p;
229  }
230  }
231 
232  if(npre > 0) {
233  for(G4int k=0; k<npre; ++k) {
234  G4ReactionProduct* p = (*thePreCompoundResult)[k];
236  p->GetMomentum()));
237  delete p;
238  }
239  }
240 
241  delete thePreCompoundResult;
242  delete theExcitationResult;
243 
244  // return the particle change
245  return &theResult;
246 }