source: trunk/source/processes/hadronic/stopping/src/G4PiMinusAbsorptionAtRest.cc@ 1331

Last change on this file since 1331 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 9.6 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// File name: G4PiMinusAbsorptionAtRest.hh
27//
28// Author: Maria Grazia Pia (pia@genova.infn.it)
29//
30// Creation date: 8 May 1998
31//
32// Modifications:
33// MGP 4 Jul 1998 Changed excitation energy calculation
34// MGP 14 Sep 1998 Fixed excitation energy calculation
35//
36// -------------------------------------------------------------------
37
38#include "G4ios.hh"
39
40#include "G4PiMinusAbsorptionAtRest.hh"
41
42#include "G4PiMinusStopLi.hh"
43#include "G4PiMinusStopC.hh"
44#include "G4PiMinusStopN.hh"
45#include "G4PiMinusStopO.hh"
46#include "G4PiMinusStopAl.hh"
47#include "G4PiMinusStopCu.hh"
48#include "G4PiMinusStopCo.hh"
49#include "G4PiMinusStopTa.hh"
50#include "G4PiMinusStopPb.hh"
51#include "G4StopTheoDeexcitation.hh"
52#include "G4StopDummyDeexcitation.hh"
53#include "G4DynamicParticle.hh"
54#include "G4DynamicParticleVector.hh"
55#include "Randomize.hh"
56#include "G4ThreeVector.hh"
57#include "G4LorentzVector.hh"
58#include "G4HadronicProcessStore.hh"
59
60// Constructor
61
62G4PiMinusAbsorptionAtRest::G4PiMinusAbsorptionAtRest(const G4String& processName,
63 G4ProcessType aType) :
64 G4VRestProcess (processName, aType)
65{
66 SetProcessSubType(fHadronAtRest);
67
68 _indexDeexcitation = 0;
69
70 if (verboseLevel>0)
71 { G4cout << GetProcessName() << " is created "<< G4endl; }
72 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
73}
74
75
76// Destructor
77
78G4PiMinusAbsorptionAtRest::~G4PiMinusAbsorptionAtRest()
79{
80 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
81}
82
83void G4PiMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
84{
85 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
86}
87
88void G4PiMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
89{
90 G4HadronicProcessStore::Instance()->PrintInfo(&p);
91}
92
93G4VParticleChange* G4PiMinusAbsorptionAtRest::AtRestDoIt(const G4Track& track, const G4Step& )
94{
95 const G4DynamicParticle* stoppedHadron = track.GetDynamicParticle();
96
97 // Check applicability
98 if (! IsApplicable(*(stoppedHadron->GetDefinition())))
99 {
100 G4cerr
101 << "G4PiMinusAbsorptionAtRest: ERROR, particle must be a pion minus!"
102 << G4endl;
103 return NULL;
104 }
105
106 // Get the current material
107 const G4Material* material = track.GetMaterial();
108
109 G4double A=-1;
110 G4double Z=-1;
111 G4double random = G4UniformRand();
112 const G4ElementVector* theElementVector = material->GetElementVector();
113 unsigned int i;
114 G4double sum = 0;
115 G4double totalsum=0;
116 for(i=0; i<material->GetNumberOfElements(); ++i)
117 {
118 if((*theElementVector)[i]->GetZ()!=1) totalsum+=material->GetFractionVector()[i];
119 }
120 for (i = 0; i<material->GetNumberOfElements(); ++i)
121 {
122 if((*theElementVector)[i]->GetZ()!=1) sum += material->GetFractionVector()[i];
123 if ( sum/totalsum > random )
124 {
125 A = (*theElementVector)[i]->GetA()*mole/g;
126 Z = (*theElementVector)[i]->GetZ();
127 break;
128 }
129 }
130
131 // Do the interaction with the nucleon cluster
132
133 G4PiMinusStopMaterial* algorithm = LoadAlgorithm(static_cast<G4int>(Z));
134 G4PiMinusStopAbsorption stopAbsorption(algorithm,Z,A);
135 stopAbsorption.SetVerboseLevel(verboseLevel);
136
137 G4DynamicParticleVector* absorptionProducts = stopAbsorption.DoAbsorption();
138
139 // Deal with the leftover nucleus
140
141 G4double pionEnergy = stoppedHadron->GetTotalEnergy();
142 G4double excitation = pionEnergy - stopAbsorption.Energy();
143 if (excitation < 0.)
144 {
145 G4Exception("G4PiMinusAbsorptionAtRest", "007", FatalException,
146 "AtRestDoIt -- excitation energy < 0");
147 }
148 if (verboseLevel>0) { G4cout << " excitation " << excitation << G4endl; }
149
150 G4StopDeexcitationAlgorithm* nucleusAlgorithm = LoadNucleusAlgorithm();
151 G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
152
153 G4double newZ = Z - stopAbsorption.NProtons();
154 G4double newN = A - Z - stopAbsorption.NNeutrons();
155 G4double newA = newZ + newN;
156 G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,excitation,stopAbsorption.RecoilMomentum());
157
158 unsigned int nAbsorptionProducts = 0;
159 if (absorptionProducts != 0)
160 { nAbsorptionProducts = absorptionProducts->size(); }
161
162 unsigned int nFragmentationProducts = 0;
163 if (fragmentationProducts != 0)
164 { nFragmentationProducts = fragmentationProducts->size(); }
165
166 if (verboseLevel>0)
167 {
168 G4cout << "nAbsorptionProducts = " << nAbsorptionProducts
169 << " nFragmentationProducts = " << nFragmentationProducts
170 << G4endl;
171 }
172
173 // Deal with ParticleChange final state
174
175 aParticleChange.Initialize(track);
176 aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts + nFragmentationProducts));
177
178 for (i = 0; i<nAbsorptionProducts; i++)
179 { aParticleChange.AddSecondary((*absorptionProducts)[i]); }
180
181// for (i = 0; i<nFragmentationProducts; i++)
182// { aParticleChange.AddSecondary(fragmentationProducts->at(i)); }
183 for(i=0; i<nFragmentationProducts; i++)
184 {
185 G4DynamicParticle * aNew =
186 new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
187 (*fragmentationProducts)[i]->GetMomentum());
188 G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
189 aParticleChange.AddSecondary(aNew, newTime);
190 delete (*fragmentationProducts)[i];
191 }
192
193 if (fragmentationProducts != 0) delete fragmentationProducts;
194
195 if (_indexDeexcitation == 1) aParticleChange.ProposeLocalEnergyDeposit(excitation);
196
197 // Kill the absorbed pion
198 aParticleChange.ProposeTrackStatus(fStopAndKill);
199
200 return &aParticleChange;
201
202}
203
204G4PiMinusStopMaterial* G4PiMinusAbsorptionAtRest::LoadAlgorithm(int Z)
205{
206 if (verboseLevel>0)
207 {
208 G4cout << "Load material algorithm " << Z << G4endl;
209 }
210
211 G4int index = 3;
212 if (Z <= 3) { index = 3;}
213 if (Z > 3 && Z<= 6) {index = 6;}
214 if (Z == 7) {index = 7;}
215 if (Z >= 8 && Z<= 11) {index = 8;}
216 if (Z >= 12 && Z<= 18) {index = 13;}
217 if (Z >=19 && Z<= 27) {index = 27;}
218 if (Z >= 28 && Z<= 51) {index = 29;}
219 if (Z >=52 ) {index = 73;}
220
221 switch (index)
222 {
223 case 3:
224 if (verboseLevel>0)
225 { G4cout << " =================== Load Li algorithm " << G4endl; }
226 return new G4PiMinusStopLi();
227 case 6:
228 if (verboseLevel>0)
229 { G4cout << " =================== Load C algorithm " << G4endl; }
230 return new G4PiMinusStopC();
231 case 7:
232 if (verboseLevel>0)
233 { G4cout << " =================== Load N algorithm " << G4endl; }
234 return new G4PiMinusStopN();
235 case 8:
236 if (verboseLevel>0)
237 { G4cout << " =================== Load O algorithm " << G4endl; }
238 return new G4PiMinusStopO();
239 case 13:
240 if (verboseLevel>0)
241 { G4cout << " =================== Load Al algorithm " << G4endl; }
242 return new G4PiMinusStopAl();
243 case 27:
244 if (verboseLevel>0)
245 { G4cout << " =================== Load Cu algorithm " << G4endl; }
246 return new G4PiMinusStopCu();
247 case 29:
248 if (verboseLevel>0)
249 { G4cout << " =================== Load Co algorithm " << G4endl; }
250 return new G4PiMinusStopCo();
251 case 73:
252 if (verboseLevel>0)
253 { G4cout << " =================== Load Ta algorithm " << G4endl; }
254 return new G4PiMinusStopTa();
255 default:
256 if (verboseLevel>0)
257 { G4cout << " =================== Load default material algorithm " << G4endl; }
258 return new G4PiMinusStopC();
259 }
260}
261
262G4StopDeexcitationAlgorithm* G4PiMinusAbsorptionAtRest::LoadNucleusAlgorithm()
263{
264
265 switch (_indexDeexcitation)
266 {
267 case 0:
268 if (verboseLevel>0)
269 { G4cout << " =================== Load Theo deexcitation " << G4endl; }
270 return new G4StopTheoDeexcitation();
271 case 1:
272 if (verboseLevel>0)
273 { G4cout << " =================== Load Dummy deexcitation " << G4endl; }
274 return new G4StopDummyDeexcitation();
275 default:
276 if (verboseLevel>0)
277 { G4cout << " =================== Load default deexcitation " << G4endl; }
278 return new G4StopTheoDeexcitation();
279 }
280}
281
282void G4PiMinusAbsorptionAtRest::SetDeexcitationAlgorithm(G4int index)
283{
284 _indexDeexcitation = index;
285}
Note: See TracBrowser for help on using the repository browser.