| Geant4 Cross Reference |
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
41 #include "G4LowEIonFragmentation.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Fancy3DNucleus.hh"
45 #include "G4Proton.hh"
46 #include "G4NucleiProperties.hh"
47 #include "G4PhysicsModelCatalog.hh"
48
49 G4LowEIonFragmentation::G4LowEIonFragmentation(G4ExcitationHandler * const value)
50 : G4HadronicInteraction("LowEIonPreco")
51 {
52 theHandler = value;
53 theModel = new G4PreCompoundModel(theHandler);
54 proton = G4Proton::Proton();
55 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
56 }
57
58 G4LowEIonFragmentation::~G4LowEIonFragmentation()
59 {
60 theResult.Clear();
61 }
62
63 G4HadFinalState * G4LowEIonFragmentation::
64 ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus & theNucleus)
65 {
66 area = 0.0;
67 // initialize the particle change
68 theResult.Clear();
69 theResult.SetStatusChange( stopAndKill );
70 theResult.SetEnergyChange( 0.0 );
71
72 // Get Target A, Z
73 G4int aTargetA = theNucleus.GetA_asInt();
74 G4int aTargetZ = theNucleus.GetZ_asInt();
75
76 // Get Projectile A, Z
77 G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
78 G4int aProjectileZ =
79 G4lrint(thePrimary.GetDefinition()->GetPDGCharge()/eplus);
80
81 // Get Maximum radius of both
82
83 G4Fancy3DNucleus aPrim;
84 aPrim.Init(aProjectileA, aProjectileZ);
85 G4double projectileOuterRadius = aPrim.GetOuterRadius();
86
87 G4Fancy3DNucleus aTarg;
88 aTarg.Init(aTargetA, aTargetZ);
89 G4double targetOuterRadius = aTarg.GetOuterRadius();
90
91 // Get the Impact parameter
92 G4int particlesFromProjectile = 0;
93 G4int chargedFromProjectile = 0;
94 G4double impactParameter = 0;
95 G4double x,y;
96 G4Nucleon * pNucleon;
97 // need at lease one particle from the projectile model beyond the
98 // projectileHorizon.
99
100 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
101 while(0==particlesFromProjectile)
102 {
103 do
104 {
105 x = 2*G4UniformRand() - 1;
106 y = 2*G4UniformRand() - 1;
107 }
108 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
109 while(x*x + y*y > 1);
110 impactParameter = std::sqrt(x*x+y*y)*
111 (targetOuterRadius+projectileOuterRadius);
112 ++totalTries;
113 area = pi*(targetOuterRadius+projectileOuterRadius)*
114 (targetOuterRadius+projectileOuterRadius);
115 G4double projectileHorizon = impactParameter-targetOuterRadius;
116
117 // Empirical boundary transparency.
118 G4double empirical = G4UniformRand();
119 if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
120
121 // Calculate the number of nucleons involved in collision
122 // From projectile
123 aPrim.StartLoop();
124
125 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
126 while((pNucleon = aPrim.GetNextNucleon()))
127 {
128 if(pNucleon->GetPosition().y()>projectileHorizon)
129 {
130 // We have one
131 ++particlesFromProjectile;
132 if(pNucleon->GetParticleType() == proton)
133 {
134 ++chargedFromProjectile;
135 }
136 }
137 }
138 }
139 ++hits;
140
141 // From target:
142 G4double targetHorizon = impactParameter-projectileOuterRadius;
143 G4int chargedFromTarget = 0;
144 G4int particlesFromTarget = 0;
145 aTarg.StartLoop();
146 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
147 while((pNucleon = aTarg.GetNextNucleon()))
148 {
149 if(pNucleon->GetPosition().y()>targetHorizon)
150 {
151 // We have one
152 ++particlesFromTarget;
153 if(pNucleon->GetParticleType() == proton)
154 {
155 ++chargedFromTarget;
156 }
157 }
158 }
159
160 // Energy sharing between projectile and target.
161 // Note that this is a quite simplistic kinetically.
162 G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
163 G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
164
165 G4double projTotEnergy = thePrimary.GetTotalEnergy();
166 G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
167 G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
168
169 // take the nucleons and fill the Fragments
170 G4Fragment anInitialState(aTargetA+particlesFromProjectile,
171 aTargetZ+chargedFromProjectile,
172 fragment4Momentum);
173 // M.A. Cortes fix
174 anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile
175 + particlesFromTarget,
176 chargedFromProjectile
177 + chargedFromTarget);
178 anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
179 chargedFromProjectile + chargedFromTarget);
180 G4double time = thePrimary.GetGlobalTime();
181 anInitialState.SetCreationTime(time);
182 anInitialState.SetCreatorModelID(secID);
183
184 // Fragment the Fragment using Pre-compound
185 G4ReactionProductVector* thePreCompoundResult =
186 theModel->DeExcite(anInitialState);
187
188 // De-excite the projectile using ExcitationHandler
189 G4ReactionProductVector * theExcitationResult = nullptr;
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 initialState2.SetCreatorModelID(secID);
206
207 theExcitationResult = theHandler->BreakItUp(initialState2);
208 }
209
210 // Fill the particle change and clear intermediate vectors
211 std::size_t nexc = (nullptr != theExcitationResult) ?
212 theExcitationResult->size() : 0;
213 std::size_t npre = (nullptr != thePreCompoundResult) ?
214 thePreCompoundResult->size() : 0;
215
216 for(std::size_t k=0; k<nexc; ++k) {
217 G4ReactionProduct* p = (*theExcitationResult)[k];
218 G4HadSecondary secondary(new G4DynamicParticle(p->GetDefinition(), p->GetMomentum()));
219 secondary.SetTime(p->GetTOF());
220 secondary.SetCreatorModelID(secID);
221 theResult.AddSecondary(secondary);
222 delete p;
223 }
224 for(std::size_t k=0; k<npre; ++k) {
225 G4ReactionProduct* p = (*thePreCompoundResult)[k];
226 G4HadSecondary secondary(new G4DynamicParticle(p->GetDefinition(), p->GetMomentum()));
227 secondary.SetTime(p->GetTOF());
228 secondary.SetCreatorModelID(secID);
229 theResult.AddSecondary(secondary);
230 delete p;
231 }
232
233 delete thePreCompoundResult;
234 delete theExcitationResult;
235
236 // return the particle change
237 return &theResult;
238 }
239