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

Last change on this file since 975 was 962, checked in by garnier, 17 years ago

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