source: trunk/source/processes/hadronic/stopping/src/G4PionMinusAbsorptionAtRest.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 14.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// G4PionMinusAbsorptionAtRest physics process
27// Larry Felawka (TRIUMF), April 1998
28//---------------------------------------------------------------------
29
30#include "G4PionMinusAbsorptionAtRest.hh"
31#include "G4DynamicParticle.hh"
32#include "G4ParticleTypes.hh"
33#include "Randomize.hh"
34#include <string.h>
35#include <cmath>
36#include <stdio.h>
37
38#define MAX_SECONDARIES 100
39
40// constructor
41
42G4PionMinusAbsorptionAtRest::G4PionMinusAbsorptionAtRest(const G4String& processName,
43 G4ProcessType aType ) :
44 G4VRestProcess (processName, aType), // initialization
45 massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
46 pdefGamma(G4Gamma::Gamma()),
47 pdefPionZero(G4PionZero::PionZero()),
48 pdefPionMinus(G4PionMinus::PionMinus()),
49 pdefProton(G4Proton::Proton()),
50 pdefNeutron(G4Neutron::Neutron()),
51 pdefDeuteron(G4Deuteron::Deuteron()),
52 pdefTriton(G4Triton::Triton()),
53 pdefAlpha(G4Alpha::Alpha())
54{
55 if (verboseLevel>0) {
56 G4cout << GetProcessName() << " is created "<< G4endl;
57 }
58 SetProcessSubType(fHadronAtRest);
59 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
60 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
61 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
62
63}
64
65// destructor
66
67G4PionMinusAbsorptionAtRest::~G4PionMinusAbsorptionAtRest()
68{
69 delete [] pv;
70 delete [] eve;
71 delete [] gkin;
72}
73
74
75// methods.............................................................................
76
77G4bool G4PionMinusAbsorptionAtRest::IsApplicable(
78 const G4ParticleDefinition& particle
79 )
80{
81 return ( &particle == pdefPionMinus );
82
83}
84
85// Warning - this method may be optimized away if made "inline"
86G4int G4PionMinusAbsorptionAtRest::GetNumberOfSecondaries()
87{
88 return ( ngkine );
89
90}
91
92// Warning - this method may be optimized away if made "inline"
93G4GHEKinematicsVector* G4PionMinusAbsorptionAtRest::GetSecondaryKinematics()
94{
95 return ( &gkin[0] );
96
97}
98
99G4double G4PionMinusAbsorptionAtRest::AtRestGetPhysicalInteractionLength(
100 const G4Track& track,
101 G4ForceCondition* condition
102 )
103{
104 // beggining of tracking
105 ResetNumberOfInteractionLengthLeft();
106
107 // condition is set to "Not Forced"
108 *condition = NotForced;
109
110 // get mean life time
111 currentInteractionLength = GetMeanLifeTime(track, condition);
112
113 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
114 G4cout << "G4PionMinusAbsorptionAtRestProcess::AtRestGetPhysicalInteractionLength ";
115 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
116 track.GetDynamicParticle()->DumpInfo();
117 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
118 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
119 }
120
121 return theNumberOfInteractionLengthLeft * currentInteractionLength;
122
123}
124
125G4VParticleChange* G4PionMinusAbsorptionAtRest::AtRestDoIt(
126 const G4Track& track,
127 const G4Step&
128 )
129//
130// Handles PionMinuss at rest; a PionMinus can either create secondaries or
131// do nothing (in which case it should be sent back to decay-handling
132// section
133//
134{
135
136// Initialize ParticleChange
137// all members of G4VParticleChange are set to equal to
138// corresponding member in G4Track
139
140 aParticleChange.Initialize(track);
141
142// Store some global quantities that depend on current material and particle
143
144 globalTime = track.GetGlobalTime()/s;
145 G4Material * aMaterial = track.GetMaterial();
146 const G4int numberOfElements = aMaterial->GetNumberOfElements();
147 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
148
149 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
150 G4double normalization = 0;
151 for ( G4int i1=0; i1 < numberOfElements; i1++ )
152 {
153 normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
154 // probabilities are included.
155 }
156 G4double runningSum= 0.;
157 G4double random = G4UniformRand()*normalization;
158 for ( G4int i2=0; i2 < numberOfElements; i2++ )
159 {
160 runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
161 // probabilities are included.
162 if (random<=runningSum)
163 {
164 targetCharge = G4double((*theElementVector)[i2]->GetZ());
165 targetAtomicMass = (*theElementVector)[i2]->GetN();
166 }
167 }
168 if (random>runningSum)
169 {
170 targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
171 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
172
173 }
174
175 if (verboseLevel>1) {
176 G4cout << "G4PionMinusAbsorptionAtRest::AtRestDoIt is invoked " <<G4endl;
177 }
178
179 G4ParticleMomentum momentum;
180 G4float localtime;
181
182 G4ThreeVector position = track.GetPosition();
183
184 GenerateSecondaries(); // Generate secondaries
185
186 aParticleChange.SetNumberOfSecondaries( ngkine );
187
188 for ( G4int isec = 0; isec < ngkine; isec++ ) {
189 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
190 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
191 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
192
193 localtime = globalTime + gkin[isec].GetTOF();
194
195 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
196 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
197 aParticleChange.AddSecondary( aNewTrack );
198
199 }
200
201 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
202
203 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident PionMinus
204
205// clear InteractionLengthLeft
206
207 ResetNumberOfInteractionLengthLeft();
208
209 return &aParticleChange;
210
211}
212
213
214void G4PionMinusAbsorptionAtRest::GenerateSecondaries()
215{
216 static G4int index;
217 static G4int l;
218 static G4int nopt;
219 static G4int i;
220 static G4ParticleDefinition* jnd;
221
222 for (i = 1; i <= MAX_SECONDARIES; ++i) {
223 pv[i].SetZero();
224 }
225
226 ngkine = 0; // number of generated secondary particles
227 ntot = 0;
228 result.SetZero();
229 result.SetMass( massPionMinus );
230 result.SetKineticEnergyAndUpdate( 0. );
231 result.SetTOF( 0. );
232 result.SetParticleDef( pdefPionMinus );
233
234 PionMinusAbsorption(&nopt);
235
236 // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
237 if (ntot != 0 || result.GetParticleDef() != pdefPionMinus) {
238 // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
239 // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
240
241 // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
242 // --- THE GEANT TEMPORARY STACK ---
243
244 // --- PUT PARTICLE ON THE STACK ---
245 gkin[0] = result;
246 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
247 ngkine = 1;
248
249 // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
250 // --- CONVENTION IS THE FOLLOWING ---
251
252 // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
253 for (l = 1; l <= ntot; ++l) {
254 index = l - 1;
255 jnd = eve[index].GetParticleDef();
256
257 // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
258 if (ngkine < MAX_SECONDARIES) {
259 gkin[ngkine] = eve[index];
260 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
261 ++ngkine;
262 }
263 }
264 }
265 else {
266 // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
267 // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
268 ngkine = 0;
269 ntot = 0;
270 globalTime += result.GetTOF() * G4float(5e-11);
271 }
272
273 // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
274 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
275
276} // GenerateSecondaries
277
278
279void G4PionMinusAbsorptionAtRest::PionMinusAbsorption(G4int *nopt)
280{
281 static G4int i;
282 static G4int nt, nbl;
283 static G4float ran, tex;
284 static G4int isw;
285 static G4float ran2, tof1, ekin;
286 static G4float ekin1, ekin2, black;
287 static G4float pnrat;
288 static G4ParticleDefinition* ipa1;
289 static G4ParticleDefinition* inve;
290
291 // *** CHARGED PION ABSORPTION BY A NUCLEUS ***
292 // *** NVE 04-MAR-1988 CERN GENEVA ***
293
294 // ORIGIN : H.FESEFELDT (09-JULY-1987)
295
296 // PANOFSKY RATIO (PI- P --> N PI0/PI- P --> N GAMMA) = 3/2
297 // FOR CAPTURE ON PROTON (HYDROGEN),
298 // STAR PRODUCTION FOR HEAVIER ELEMENTS
299
300 pv[1].SetZero();
301 pv[1].SetMass( massPionMinus );
302 pv[1].SetKineticEnergyAndUpdate( 0. );
303 pv[1].SetTOF( result.GetTOF() );
304 pv[1].SetParticleDef( result.GetParticleDef() );
305 if (targetAtomicMass <= G4float(1.5)) {
306 ran = G4UniformRand();
307 isw = 1;
308 if (ran < G4float(.33)) {
309 isw = 2;
310 }
311 *nopt = isw;
312 ran = G4UniformRand();
313 tof1 = std::log(ran) * G4float(-25.);
314 tof1 *= G4float(20.);
315 if (isw != 1) {
316 pv[2].SetZero();
317 pv[2].SetMass( 0. );
318 pv[2].SetKineticEnergyAndUpdate( .02 );
319 pv[2].SetTOF( result.GetTOF() + tof1 );
320 pv[2].SetParticleDef( pdefGamma );
321 }
322 else {
323 pv[2] = pv[1];
324 pv[2].SetTOF( result.GetTOF() + tof1 );
325 pv[2].SetParticleDef( pdefPionZero );
326 }
327 result = pv[2];
328 }
329 else {
330 // **
331 // ** STAR PRODUCTION FOR PION ABSORPTION IN HEAVY ELEMENTS
332 // **
333 evapEnergy1 = G4float(.0135);
334 evapEnergy3 = G4float(.0058);
335 nt = 1;
336 tex = evapEnergy1;
337 black = std::log(targetAtomicMass) * G4float(.5);
338 Poisso(black, &nbl);
339 if (nbl <= 0) {
340 nbl = 1;
341 }
342 if (nt + nbl > (MAX_SECONDARIES - 2)) {
343 nbl = (MAX_SECONDARIES - 2) - nt;
344 }
345 ekin = tex / nbl;
346 ekin2 = G4float(0.);
347 for (i = 1; i <= nbl; ++i) {
348 if (nt == (MAX_SECONDARIES - 2)) {
349 continue;
350 }
351 ran2 = G4UniformRand();
352 ekin1 = -G4double(ekin) * std::log(ran2);
353 ekin2 += ekin1;
354 ipa1 = pdefNeutron;
355 pnrat = G4float(1.) - targetCharge / targetAtomicMass;
356 if (G4UniformRand() > pnrat) {
357 ipa1 = pdefProton;
358 }
359 ++nt;
360 pv[nt].SetZero();
361 pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
362 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
363 pv[nt].SetTOF( 2. );
364 pv[nt].SetParticleDef( ipa1 );
365 if (ekin2 > tex) {
366 break;
367 }
368 }
369 tex = evapEnergy3;
370 black = std::log(targetAtomicMass) * G4float(.5);
371 Poisso(black, &nbl);
372 if (nt + nbl > (MAX_SECONDARIES - 2)) {
373 nbl = (MAX_SECONDARIES - 2) - nt;
374 }
375 if (nbl <= 0) {
376 nbl = 1;
377 }
378 ekin = tex / nbl;
379 ekin2 = G4float(0.);
380 for (i = 1; i <= nbl; ++i) {
381 if (nt == (MAX_SECONDARIES - 2)) {
382 continue;
383 }
384 ran2 = G4UniformRand();
385 ekin1 = -G4double(ekin) * std::log(ran2);
386 ekin2 += ekin1;
387 ++nt;
388 ran = G4UniformRand();
389 inve= pdefDeuteron;
390 if (ran > G4float(.6)) {
391 inve = pdefTriton;
392 }
393 if (ran > G4float(.9)) {
394 inve = pdefAlpha;
395 }
396 pv[nt].SetZero();
397 pv[nt].SetMass( inve->GetPDGMass()/GeV );
398 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
399 pv[nt].SetTOF( 2. );
400 pv[nt].SetParticleDef( inve );
401 if (ekin2 > tex) {
402 break;
403 }
404 }
405 // **
406 // ** STORE ON EVENT COMMON
407 // **
408 ran = G4UniformRand();
409 tof1 = std::log(ran) * G4float(-25.);
410 tof1 *= G4float(20.);
411 for (i = 2; i <= nt; ++i) {
412 pv[i].SetTOF( result.GetTOF() + tof1 );
413 }
414 result = pv[2];
415 for (i = 3; i <= nt; ++i) {
416 if (ntot >= MAX_SECONDARIES) {
417 break;
418 }
419 eve[ntot++] = pv[i];
420 }
421 }
422
423} // PionMinusAbsorption
424
425
426void G4PionMinusAbsorptionAtRest::Poisso(G4float xav, G4int *iran)
427{
428 static G4int i;
429 static G4float r, p1, p2, p3;
430 static G4int mm;
431 static G4float rr, ran, rrr, ran1;
432
433 // *** GENERATION OF POISSON DISTRIBUTION ***
434 // *** NVE 16-MAR-1988 CERN GENEVA ***
435 // ORIGIN : H.FESEFELDT (27-OCT-1983)
436
437 // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
438 if (xav > G4float(9.9)) {
439 // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
440 Normal(&ran1);
441 ran1 = xav + ran1 * std::sqrt(xav);
442 *iran = G4int(ran1);
443 if (*iran < 0) {
444 *iran = 0;
445 }
446 }
447 else {
448 mm = G4int(xav * G4float(5.));
449 *iran = 0;
450 if (mm > 0) {
451 r = std::exp(-G4double(xav));
452 ran1 = G4UniformRand();
453 if (ran1 > r) {
454 rr = r;
455 for (i = 1; i <= mm; ++i) {
456 ++(*iran);
457 if (i <= 5) {
458 rrr = std::pow(xav, G4float(i)) / NFac(i);
459 }
460 // ** STIRLING' S FORMULA FOR LARGE NUMBERS
461 if (i > 5) {
462 rrr = std::exp(i * std::log(xav) -
463 (i + G4float(.5)) * std::log(i * G4float(1.)) +
464 i - G4float(.9189385));
465 }
466 rr += r * rrr;
467 if (ran1 <= rr) {
468 break;
469 }
470 }
471 }
472 }
473 else {
474 // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
475 p1 = xav * std::exp(-G4double(xav));
476 p2 = xav * p1 / G4float(2.);
477 p3 = xav * p2 / G4float(3.);
478 ran = G4UniformRand();
479 if (ran >= p3) {
480 if (ran >= p2) {
481 if (ran >= p1) {
482 *iran = 0;
483 }
484 else {
485 *iran = 1;
486 }
487 }
488 else {
489 *iran = 2;
490 }
491 }
492 else {
493 *iran = 3;
494 }
495 }
496 }
497
498} // Poisso
499
500
501G4int G4PionMinusAbsorptionAtRest::NFac(G4int n)
502{
503 G4int ret_val;
504
505 static G4int i, m;
506
507 // *** NVE 16-MAR-1988 CERN GENEVA ***
508 // ORIGIN : H.FESEFELDT (27-OCT-1983)
509
510 ret_val = 1;
511 m = n;
512 if (m > 1) {
513 if (m > 10) {
514 m = 10;
515 }
516 for (i = 2; i <= m; ++i) {
517 ret_val *= i;
518 }
519 }
520 return ret_val;
521
522} // NFac
523
524
525void G4PionMinusAbsorptionAtRest::Normal(G4float *ran)
526{
527 static G4int i;
528
529 // *** NVE 14-APR-1988 CERN GENEVA ***
530 // ORIGIN : H.FESEFELDT (27-OCT-1983)
531
532 *ran = G4float(-6.);
533 for (i = 1; i <= 12; ++i) {
534 *ran += G4UniformRand();
535 }
536
537} // Normal
Note: See TracBrowser for help on using the repository browser.