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

Last change on this file since 1330 was 1055, checked in by garnier, 17 years ago

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