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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 14.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//   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
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.