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

Last change on this file since 1307 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 15.2 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 "G4HadronicProcessStore.hh"
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  }
59  SetProcessSubType(fHadronAtRest);
60  pv   = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
61  eve  = new G4GHEKinematicsVector [MAX_SECONDARIES];
62  gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
63
64  G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
65}
66 
67// destructor
68 
69G4PionMinusAbsorptionAtRest::~G4PionMinusAbsorptionAtRest()
70{
71  G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
72  delete [] pv;
73  delete [] eve;
74  delete [] gkin;
75}
76 
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}
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.