source: trunk/source/processes/hadronic/models/rpg/src/G4RPGInelastic.cc @ 1058

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

maj sur la beta de geant 4.9.3
File size: 20.4 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// $Id: G4RPGInelastic.cc,v 1.7 2009/05/25 19:07:15 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28//
29
30#include "G4RPGInelastic.hh"
31#include "Randomize.hh"
32#include "G4HadReentrentException.hh"
33#include "G4RPGStrangeProduction.hh"
34#include "G4RPGTwoBody.hh"
35
36
37G4RPGInelastic::G4RPGInelastic(const G4String& modelName)
38  : G4HadronicInteraction(modelName)
39{
40  cache = 0.0;
41  particleDef[0] = G4PionZero::PionZero();
42  particleDef[1] = G4PionPlus::PionPlus();
43  particleDef[2] = G4PionMinus::PionMinus();
44  particleDef[3] = G4KaonPlus::KaonPlus();
45  particleDef[4] = G4KaonMinus::KaonMinus();
46  particleDef[5] = G4KaonZero::KaonZero();
47  particleDef[6] = G4AntiKaonZero::AntiKaonZero();
48  particleDef[7] = G4Proton::Proton();
49  particleDef[8] = G4Neutron::Neutron();
50  particleDef[9] = G4Lambda::Lambda();
51  particleDef[10] = G4SigmaPlus::SigmaPlus();
52  particleDef[11] = G4SigmaZero::SigmaZero();
53  particleDef[12] = G4SigmaMinus::SigmaMinus();
54  particleDef[13] = G4XiZero::XiZero();
55  particleDef[14] = G4XiMinus::XiMinus();
56  particleDef[15] = G4OmegaMinus::OmegaMinus();
57  particleDef[16] = G4AntiProton::AntiProton();
58  particleDef[17] = G4AntiNeutron::AntiNeutron();
59}
60
61
62G4double G4RPGInelastic::Pmltpc(G4int np, G4int nm, G4int nz, 
63                                G4int n, G4double b, G4double c)
64{
65  const G4double expxu =  82.;           // upper bound for arg. of exp
66  const G4double expxl = -expxu;         // lower bound for arg. of exp
67  G4double npf = 0.0;
68  G4double nmf = 0.0;
69  G4double nzf = 0.0;
70  G4int i;
71  for( i=2; i<=np; i++ )npf += std::log((double)i);
72  for( i=2; i<=nm; i++ )nmf += std::log((double)i);
73  for( i=2; i<=nz; i++ )nzf += std::log((double)i);
74  G4double r;
75  r = std::min( expxu, std::max( expxl, -(np-nm+nz+b)*(np-nm+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
76  return std::exp(r);
77}
78
79
80G4int G4RPGInelastic::Factorial( G4int n )
81{
82  G4int m = std::min(n,10);
83  G4int result = 1;
84  if( m <= 1 )return result;
85  for( G4int i=2; i<=m; ++i )result *= i;
86  return result;
87}
88
89
90G4bool G4RPGInelastic::MarkLeadingStrangeParticle(
91     const G4ReactionProduct &currentParticle,
92     const G4ReactionProduct &targetParticle,
93     G4ReactionProduct &leadParticle )
94{
95  // The following was in GenerateXandPt and TwoCluster.
96  // Add a parameter to the GenerateXandPt function telling it about the
97  // strange particle.
98  //
99  // Assumes that the original particle was a strange particle
100  //
101  G4bool lead = false;
102  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
103      (currentParticle.GetDefinition() != G4Proton::Proton()) &&
104      (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
105  {
106    lead = true;
107    leadParticle = currentParticle;              //  set lead to the incident particle
108  }
109  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
110           (targetParticle.GetDefinition() != G4Proton::Proton()) &&
111           (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
112  {
113    lead = true;
114    leadParticle = targetParticle;              //   set lead to the target particle
115  }
116  return lead;
117}
118
119 
120 void G4RPGInelastic::SetUpPions(const G4int np, const G4int nm, 
121                                 const G4int nz,
122                             G4FastVector<G4ReactionProduct,256> &vec,
123                                 G4int &vecLen)
124 {
125   if( np+nm+nz == 0 )return;
126   G4int i;
127   G4ReactionProduct *p;
128   for( i=0; i<np; ++i )
129   {
130     p = new G4ReactionProduct;
131     p->SetDefinition( G4PionPlus::PionPlus() );
132     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
133     vec.SetElement( vecLen++, p );
134   }
135   for( i=np; i<np+nm; ++i )
136   {
137     p = new G4ReactionProduct;
138     p->SetDefinition( G4PionMinus::PionMinus() );
139     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
140     vec.SetElement( vecLen++, p );
141   }
142   for( i=np+nm; i<np+nm+nz; ++i )
143   {
144     p = new G4ReactionProduct;
145     p->SetDefinition( G4PionZero::PionZero() );
146     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
147     vec.SetElement( vecLen++, p );
148   }
149 }
150
151
152 void G4RPGInelastic::GetNormalizationConstant(
153   const G4double energy,  // MeV,  <0 means annihilation channels
154   G4double &n,
155   G4double &anpn )
156  {
157    const G4double expxu =  82.;          // upper bound for arg. of exp
158    const G4double expxl = -expxu;        // lower bound for arg. of exp
159    const G4int numSec = 60;
160    //
161    // the only difference between the calculation for annihilation channels
162    // and normal is the starting value, iBegin, for the loop below
163    //
164    G4int iBegin = 1;
165    G4double en = energy;
166    if( energy < 0.0 )
167    {
168      iBegin = 2;
169      en *= -1.0;
170    }
171    //
172    // number of total particles vs. centre of mass Energy - 2*proton mass
173    //
174    G4double aleab = std::log(en/GeV);
175    n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
176    n -= 2.0;
177    //
178    // normalization constant for kno-distribution
179    //
180    anpn = 0.0;
181    G4double test, temp;
182    for( G4int i=iBegin; i<=numSec; ++i )
183    {
184      temp = pi*i/(2.0*n*n);
185      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
186      if( temp < 1.0 )
187      {
188        if( test >= 1.0e-10 )anpn += temp*test;
189      }
190      else
191        anpn += temp*test;
192    }
193  }
194 
195void 
196G4RPGInelastic::CalculateMomenta(G4FastVector<G4ReactionProduct,256>& vec,
197                                 G4int& vecLen,
198                                 const G4HadProjectile* originalIncident,
199                                 const G4DynamicParticle* originalTarget,
200                                 G4ReactionProduct& modifiedOriginal,
201                                 G4Nucleus& targetNucleus,
202                                 G4ReactionProduct& currentParticle,
203                                 G4ReactionProduct& targetParticle,
204                                 G4bool& incidentHasChanged,
205                                 G4bool& targetHasChanged,
206                                 G4bool quasiElastic)
207{
208  cache = 0;
209  what = originalIncident->Get4Momentum().vect();
210
211  G4ReactionProduct leadingStrangeParticle;
212
213  //  strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
214  //                                  incidentHasChanged, originalTarget,
215  //                                  targetParticle, targetHasChanged,
216  //                                  targetNucleus, currentParticle,
217  //                                  vec, vecLen,
218  //                              false, leadingStrangeParticle);
219
220  if( quasiElastic )
221  {
222    twoBody.ReactionStage(originalIncident, modifiedOriginal,
223                          incidentHasChanged, originalTarget,
224                          targetParticle, targetHasChanged,
225                          targetNucleus, currentParticle,
226                          vec, vecLen,
227                          false, leadingStrangeParticle);
228    return;
229  }
230
231  G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
232                                               targetParticle,
233                                               leadingStrangeParticle );
234  //
235  // Note: the number of secondaries can be reduced in GenerateXandPt
236  // and TwoCluster
237  //
238  G4bool finishedGenXPt = false;
239  G4bool annihilation = false;
240  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
241      currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
242  {
243    // original was an anti-particle and annihilation has taken place
244    annihilation = true;
245    G4double ekcor = 1.0;
246    G4double ek = originalIncident->GetKineticEnergy();
247    G4double ekOrg = ek;
248     
249    const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
250    if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
251    const G4double atomicWeight = targetNucleus.GetN();
252    ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
253    G4double tkin = targetNucleus.Cinema(ek);
254    ek += tkin;
255    ekOrg += tkin;
256    //      modifiedOriginal.SetKineticEnergy( ekOrg );
257    //
258    // evaporation --  re-calculate black track energies
259    //                 this was Done already just before the cascade
260    //
261    tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
262    ekOrg -= tkin;
263    ekOrg = std::max( 0.0001*GeV, ekOrg );
264    modifiedOriginal.SetKineticEnergy( ekOrg );
265    G4double amas = originalIncident->GetDefinition()->GetPDGMass();
266    G4double et = ekOrg + amas;
267    G4double p = std::sqrt( std::abs(et*et-amas*amas) );
268    G4double pp = modifiedOriginal.GetMomentum().mag();
269    if( pp > 0.0 )
270    {
271      G4ThreeVector momentum = modifiedOriginal.GetMomentum();
272      modifiedOriginal.SetMomentum( momentum * (p/pp) );
273    }
274    if( ekOrg <= 0.0001 )
275    {
276      modifiedOriginal.SetKineticEnergy( 0.0 );
277      modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
278    }
279  }
280
281  // twsup gives percentage of time two-cluster model is called
282
283  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
284  G4double rand1 = G4UniformRand();
285  G4double rand2 = G4UniformRand();
286
287  // Cache current, target, and secondaries
288  G4ReactionProduct saveCurrent = currentParticle;
289  G4ReactionProduct saveTarget = targetParticle;
290  std::vector<G4ReactionProduct> savevec;
291  for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
292
293  // Call fragmentation code if
294  //   1) there is annihilation, or
295  //   2) there are more than 5 secondaries, or
296  //   3) incident KE is > 1 GeV AND
297  //        ( incident is a kaon AND rand < 0.5 OR twsup )
298  //
299
300  if( annihilation || vecLen > 5 ||
301      ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
302
303      (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
304         originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
305         originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
306         originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
307          rand1 < 0.5) 
308       || rand2 > twsup[vecLen]) ) )
309
310    finishedGenXPt =
311      fragmentation.ReactionStage(originalIncident, modifiedOriginal,
312                                  incidentHasChanged, originalTarget,
313                                  targetParticle, targetHasChanged,
314                                  targetNucleus, currentParticle,
315                                  vec, vecLen,
316                                  leadFlag, leadingStrangeParticle);
317
318  if (finishedGenXPt) return;
319
320  G4bool finishedTwoClu = false;
321
322  if (modifiedOriginal.GetTotalMomentum() < 1.0) {
323    for (G4int i = 0; i < vecLen; i++) delete vec[i];
324    vecLen = 0;
325
326  } else {
327    // Occaisionally, GenerateXandPt will fail in the annihilation channel.
328    // Restore current, target and secondaries to pre-GenerateXandPt state
329    // before trying annihilation in TwoCluster
330
331    if (!finishedGenXPt && annihilation) {
332      currentParticle = saveCurrent;
333      targetParticle = saveTarget;
334      for (G4int i = 0; i < vecLen; i++) delete vec[i];
335      vecLen = 0;
336      vec.Initialize( 0 );
337      for (G4int i = 0; i < G4int(savevec.size()); i++) {
338        G4ReactionProduct* p = new G4ReactionProduct;
339        *p = savevec[i];
340        vec.SetElement( vecLen++, p );
341      }
342    }
343
344    // Big violations of energy conservation in this method - don't use
345    //
346    //    pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
347    //                                  incidentHasChanged, originalTarget,
348    //                                  targetParticle, targetHasChanged,
349    //                                  targetNucleus, currentParticle,
350    //                                  vec, vecLen,
351    //                                  false, leadingStrangeParticle);
352
353    try
354    {
355      finishedTwoClu = 
356        twoCluster.ReactionStage(originalIncident, modifiedOriginal,
357                                 incidentHasChanged, originalTarget,
358                                 targetParticle, targetHasChanged,
359                                 targetNucleus, currentParticle,
360                                 vec, vecLen,
361                                 leadFlag, leadingStrangeParticle);
362    }
363     catch(G4HadReentrentException aC)
364    {
365       aC.Report(G4cout);
366       throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
367    }
368  }
369
370  if (finishedTwoClu) return;
371
372  twoBody.ReactionStage(originalIncident, modifiedOriginal,
373                        incidentHasChanged, originalTarget,
374                        targetParticle, targetHasChanged,
375                        targetNucleus, currentParticle,
376                        vec, vecLen,
377                        false, leadingStrangeParticle);
378}
379
380/*
381 void G4RPGInelastic::
382 Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
383 {
384   G4double rotation = 2.*pi*G4UniformRand();
385   cache = rotation;
386   G4int i;
387   for( i=0; i<vecLen; ++i )
388   {
389     G4ThreeVector momentum = vec[i]->GetMomentum();
390     momentum = momentum.rotate(rotation, what);
391     vec[i]->SetMomentum(momentum);
392   }
393 }     
394*/
395
396void 
397G4RPGInelastic::SetUpChange(G4FastVector<G4ReactionProduct,256>& vec,
398                            G4int& vecLen,
399                            G4ReactionProduct& currentParticle,
400                            G4ReactionProduct& targetParticle,
401                            G4bool& incidentHasChanged )
402{
403  theParticleChange.Clear();
404  G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
405  G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
406  G4int i;
407
408  if (currentParticle.GetDefinition() == particleDef[k0]) {
409    if (G4UniformRand() < 0.5) {
410      currentParticle.SetDefinitionAndUpdateE(aKaonZL);
411      incidentHasChanged = true;
412    } else {
413      currentParticle.SetDefinitionAndUpdateE(aKaonZS);
414    }
415  } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
416    if (G4UniformRand() < 0.5) {
417      currentParticle.SetDefinitionAndUpdateE(aKaonZL);
418    } else {
419      currentParticle.SetDefinitionAndUpdateE(aKaonZS);
420      incidentHasChanged = true;
421    }
422  }
423
424  if (targetParticle.GetDefinition() == particleDef[k0] || 
425      targetParticle.GetDefinition() == particleDef[k0b] ) {
426    if (G4UniformRand() < 0.5) {
427      targetParticle.SetDefinitionAndUpdateE(aKaonZL);
428    } else {
429      targetParticle.SetDefinitionAndUpdateE(aKaonZS);
430    }
431  }
432
433  for (i = 0; i < vecLen; ++i) {
434    if (vec[i]->GetDefinition() == particleDef[k0] ||
435        vec[i]->GetDefinition() == particleDef[k0b] ) {
436      if (G4UniformRand() < 0.5) {
437        vec[i]->SetDefinitionAndUpdateE(aKaonZL);
438      } else {
439        vec[i]->SetDefinitionAndUpdateE(aKaonZS);
440      }
441    }
442  }
443
444  if (incidentHasChanged) {
445    G4DynamicParticle* p0 = new G4DynamicParticle;
446    p0->SetDefinition(currentParticle.GetDefinition() );
447    p0->SetMomentum(currentParticle.GetMomentum() );
448    theParticleChange.AddSecondary( p0 );
449    theParticleChange.SetStatusChange( stopAndKill );
450    theParticleChange.SetEnergyChange( 0.0 );
451
452  } else {
453    G4double p = currentParticle.GetMomentum().mag()/MeV;
454    G4ThreeVector m = currentParticle.GetMomentum();
455    if (p > DBL_MIN)
456      theParticleChange.SetMomentumChange( m.x()/p, m.y()/p, m.z()/p );
457    else
458      theParticleChange.SetMomentumChange( 0.0, 0.0, 1.0 );
459
460    G4double aE = currentParticle.GetKineticEnergy();
461    if (std::fabs(aE)<.1*eV) aE=.1*eV;
462    theParticleChange.SetEnergyChange( aE );
463  }
464
465  if (targetParticle.GetMass() > 0.0)  // Tgt particle can be eliminated in TwoBody
466  {
467    G4ThreeVector momentum = targetParticle.GetMomentum();
468    momentum = momentum.rotate(cache, what);
469    G4double targKE = targetParticle.GetKineticEnergy();
470    G4ThreeVector dir(0.0, 0.0, 1.0);
471    if (targKE < DBL_MIN)
472      targKE = DBL_MIN;
473    else
474      dir = momentum/momentum.mag();
475
476    G4DynamicParticle* p1 = 
477      new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
478
479    theParticleChange.AddSecondary( p1 );
480  }
481
482  G4DynamicParticle* p;
483  for (i = 0; i < vecLen; ++i) {
484    G4double secKE = vec[i]->GetKineticEnergy();
485    G4ThreeVector momentum = vec[i]->GetMomentum();
486    G4ThreeVector dir(0.0, 0.0, 1.0);
487    if (secKE < DBL_MIN)
488      secKE = DBL_MIN;
489    else
490      dir = momentum/momentum.mag();
491
492    p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
493    theParticleChange.AddSecondary( p );
494    delete vec[i];
495  }
496}
497
498
499std::pair<G4int, G4double>
500G4RPGInelastic::interpolateEnergy(G4double e) const
501{
502  G4int index = 29;
503  G4double fraction = 0.0;
504
505  for (G4int i = 1; i < 30; i++) {
506    if (e < energyScale[i]) {
507      index = i-1;
508      fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
509      break;
510    }
511  }
512  return std::pair<G4int, G4double>(index, fraction);
513}
514
515
516G4int
517G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
518{
519  G4int i;
520  G4double sum(0.);
521  for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
522
523  G4double fsum = sum*G4UniformRand();
524  G4double partialSum = 0.0;
525  G4int channel = 0;
526
527  for (i = 0; i < G4int(sigma.size()); i++) {
528    partialSum += sigma[i];
529    if (fsum < partialSum) {
530      channel = i;
531      break;
532    }
533  }
534
535  return channel;
536}
537
538
539void G4RPGInelastic::CheckQnums(G4FastVector<G4ReactionProduct,256> &vec,
540                                G4int &vecLen,
541                                G4ReactionProduct &currentParticle,
542                                G4ReactionProduct &targetParticle,
543                                G4double Q, G4double B, G4double S)
544{
545  G4ParticleDefinition* projDef = currentParticle.GetDefinition();
546  G4ParticleDefinition* targDef = targetParticle.GetDefinition();
547  G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
548  G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
549  G4double strangenessSum = projDef->GetQuarkContent(3) - 
550                            projDef->GetAntiQuarkContent(3) + 
551                            targDef->GetQuarkContent(3) -
552                            targDef->GetAntiQuarkContent(3);
553
554  G4ParticleDefinition* secDef = 0;
555  for (G4int i = 0; i < vecLen; i++) {
556    secDef = vec[i]->GetDefinition();
557    chargeSum += secDef->GetPDGCharge();
558    baryonSum += secDef->GetBaryonNumber();
559    strangenessSum += secDef->GetQuarkContent(3) 
560                    - secDef->GetAntiQuarkContent(3);
561  }
562
563  G4bool OK = true;
564  if (chargeSum != Q) {
565    G4cout << " Charge not conserved " << G4endl;
566    OK = false;
567  }
568  if (baryonSum != B) {
569    G4cout << " Baryon number not conserved " << G4endl;
570    OK = false;
571  }
572  if (strangenessSum != S) {
573    G4cout << " Strangeness not conserved " << G4endl;
574    OK = false;
575  } 
576
577  if (!OK) {
578    G4cout << " projectile: " << projDef->GetParticleName() 
579           << "  target: " << targDef->GetParticleName() << G4endl;
580    for (G4int i = 0; i < vecLen; i++) {
581      secDef = vec[i]->GetDefinition();
582      G4cout << secDef->GetParticleName() << " " ;
583    }
584    G4cout << G4endl;
585  }
586
587}
588
589
590const G4double G4RPGInelastic::energyScale[30] = {
591  0.0,  0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
592  0.13, 0.18, 0.24,  0.32,  0.42,  0.56,  0.75,  1.0,   1.3,   1.8,
593  2.4,  3.2,  4.2,   5.6,   7.5,   10.0,  13.0,  18.0,  24.0, 32.0 };
594
595
596/* end of file */
Note: See TracBrowser for help on using the repository browser.