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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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