source: trunk/source/processes/hadronic/models/rpg/src/G4RPGStrangeProduction.cc @ 1055

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

update to geant4.9.2
File size: 15.5 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: G4RPGStrangeProduction.cc,v 1.1 2007/07/18 21:04:21 dennis Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29 
30#include "G4RPGStrangeProduction.hh"
31// #include "G4AntiProton.hh"
32// #include "G4AntiNeutron.hh"
33#include "Randomize.hh"
34#include <iostream>
35#include "G4HadReentrentException.hh"
36#include <signal.h>
37
38
39G4RPGStrangeProduction::G4RPGStrangeProduction()
40  : G4RPGReaction() {}
41
42
43G4bool G4RPGStrangeProduction::
44ReactionStage(const G4HadProjectile* /*originalIncident*/,
45              G4ReactionProduct& modifiedOriginal,
46              G4bool& incidentHasChanged,
47              const G4DynamicParticle* originalTarget,
48              G4ReactionProduct& targetParticle,
49              G4bool& targetHasChanged,
50              const G4Nucleus& /*targetNucleus*/,
51              G4ReactionProduct& currentParticle,
52              G4FastVector<G4ReactionProduct,256>& vec,
53              G4int& vecLen,
54              G4bool /*leadFlag*/,
55              G4ReactionProduct& /*leadingStrangeParticle*/)
56{
57  // Derived from H. Fesefeldt's original FORTRAN code STPAIR
58  //
59  // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
60  //                            K+ Y0, K0 Y+,  K0 Y-
61  // For antibaryon induced reactions half of the cross sections KB YB
62  // pairs are produced.  Charge is not conserved, no experimental data available
63  // for exclusive reactions, therefore some average behaviour assumed.
64  // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
65  //
66
67  if( vecLen == 0 )return true;
68  //
69  // the following protects against annihilation processes
70  //
71  if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return true;
72   
73  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
74  const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
75  G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
76  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
77                                      targetMass*targetMass +
78                                      2.0*targetMass*etOriginal );  // GeV
79  G4double currentMass = currentParticle.GetMass()/GeV;
80  G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
81  if( availableEnergy <= 1.0 )return true;
82   
83  G4ParticleDefinition *aProton = G4Proton::Proton();
84  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
85  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
86  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
87  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
88  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
89  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
90  G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
91  G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
92  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
93  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
94  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
95  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
96  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
97  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
98  G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
99   
100  const G4double protonMass = aProton->GetPDGMass()/GeV;
101  const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
102  //
103  // determine the center of mass energy bin
104  //
105  const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
106
107  G4int ibin, i3, i4;
108  G4double avk, avy, avn, ran;
109  G4int i = 1;
110  while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
111  if( i == 12 )
112    ibin = 11;
113  else
114    ibin = i;
115  //
116  // the fortran code chooses a random replacement of produced kaons
117  //  but does not take into account charge conservation
118  //
119  if( vecLen == 1 )  // we know that vecLen > 0
120  {
121    i3 = 0;
122    i4 = 1;   // note that we will be adding a new secondary particle in this case only
123  }
124  else               // otherwise  0 <= i3,i4 < vecLen
125  {
126    G4double ran = G4UniformRand();
127    while( ran == 1.0 )ran = G4UniformRand();
128    i4 = i3 = G4int( vecLen*ran );
129    while( i3 == i4 )
130    {
131      ran = G4UniformRand();
132      while( ran == 1.0 )ran = G4UniformRand();
133      i4 = G4int( vecLen*ran );
134    }
135  }
136
137  //
138  // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
139  //
140  const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
141                             0.0975, 0.123, 0.28,  0.398,  0.495,  0.573 };
142  const G4double avky[] = { 0.005, 0.03,  0.064, 0.095, 0.115, 0.13,
143                            0.145, 0.155, 0.20,  0.205, 0.210, 0.212 };
144  const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
145                             0.04,    0.05,   0.12,   0.15,   0.18, 0.20 };
146   
147  avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
148    /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
149  avk = std::exp(avk);
150   
151  avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
152    /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
153  avy = std::exp(avy);
154   
155  avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
156    /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
157  avn = std::exp(avn);
158   
159  if( avk+avy+avn <= 0.0 )return true;
160   
161  if( currentMass < protonMass )avy /= 2.0;
162  if( targetMass < protonMass )avy = 0.0;
163  avy += avk+avn;
164  avk += avn;
165  ran = G4UniformRand();
166  if(  ran < avn )
167  {
168    if( availableEnergy < 2.0 )return true;
169    if( vecLen == 1 )                              // add a new secondary
170    {
171      G4ReactionProduct *p1 = new G4ReactionProduct;
172      if( G4UniformRand() < 0.5 )
173      {
174        vec[0]->SetDefinition( aNeutron );
175        p1->SetDefinition( anAntiNeutron );
176        (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
177        vec[0]->SetMayBeKilled(false);
178        p1->SetMayBeKilled(false);
179      }
180      else
181      {
182        vec[0]->SetDefinition( aProton );
183        p1->SetDefinition( anAntiProton );
184        (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
185        vec[0]->SetMayBeKilled(false);
186        p1->SetMayBeKilled(false);
187      }
188      vec.SetElement( vecLen++, p1 );
189    // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
190    }
191    else
192    {                                             // replace two secondaries
193      if( G4UniformRand() < 0.5 )
194      {
195        vec[i3]->SetDefinition( aNeutron );
196        vec[i4]->SetDefinition( anAntiNeutron );
197        vec[i3]->SetMayBeKilled(false);
198        vec[i4]->SetMayBeKilled(false);
199      }
200      else
201      {
202        vec[i3]->SetDefinition( aProton );
203        vec[i4]->SetDefinition( anAntiProton );
204        vec[i3]->SetMayBeKilled(false);
205        vec[i4]->SetMayBeKilled(false);
206      }
207    }
208  }
209  else if( ran < avk )
210  {
211    if( availableEnergy < 1.0 )return true;
212     
213    const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
214                               0.6875, 0.7500, 0.8750, 1.000 };
215    const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
216    const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
217    ran = G4UniformRand();
218    i = 0;
219    while( (i<9) && (ran>=kkb[i]) )++i;
220    if( i == 9 )return true;
221    //
222    // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
223    // charge       +  -   +  0   +  0   0  0   0  0   0  0   0  0   0  -   0  -
224    //
225    switch( ipakkb1[i] )
226    {
227     case 10:
228       vec[i3]->SetDefinition( aKaonPlus );
229       vec[i3]->SetMayBeKilled(false);
230       break;
231     case 11:
232       vec[i3]->SetDefinition( aKaonZS );
233       vec[i3]->SetMayBeKilled(false);
234       break;
235     case 12:
236       vec[i3]->SetDefinition( aKaonZL );
237       vec[i3]->SetMayBeKilled(false);
238       break;
239    }
240
241    if( vecLen == 1 )                          // add a secondary
242    {
243      G4ReactionProduct *p1 = new G4ReactionProduct;
244      switch( ipakkb2[i] )
245      {
246       case 11:
247         p1->SetDefinition( aKaonZS );
248         p1->SetMayBeKilled(false);
249         break;
250       case 12:
251         p1->SetDefinition( aKaonZL );
252         p1->SetMayBeKilled(false);
253         break;
254       case 13:
255         p1->SetDefinition( aKaonMinus );
256         p1->SetMayBeKilled(false);
257         break;
258      }
259      (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
260      vec.SetElement( vecLen++, p1 );
261
262    }
263    else                                        // replace
264    {
265      switch( ipakkb2[i] )
266      {
267       case 11:
268         vec[i4]->SetDefinition( aKaonZS );
269         vec[i4]->SetMayBeKilled(false);
270         break;
271       case 12:
272         vec[i4]->SetDefinition( aKaonZL );
273         vec[i4]->SetMayBeKilled(false);
274         break;
275       case 13:
276         vec[i4]->SetDefinition( aKaonMinus );
277         vec[i4]->SetMayBeKilled(false);
278         break;
279      }
280    }
281  }
282  else if( ran < avy )
283  {
284    if( availableEnergy < 1.6 )return true;
285     
286    const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
287                            0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
288    const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
289    const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
290    const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
291    const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
292    ran = G4UniformRand();
293    i = 0;
294    while( (i<12) && (ran>ky[i]) )++i;
295    if( i == 12 )return true;
296      if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
297      {
298        // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
299        //             0  +   0  0   0  0   +  +   +  0   +  0
300        //
301        //             21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
302        //             0  +   0  0   0  0   -  +   -  0   -  0
303        switch( ipaky1[i] )
304        {
305         case 18:
306           targetParticle.SetDefinition( aLambda );
307           break;
308         case 20:
309           targetParticle.SetDefinition( aSigmaPlus );
310           break;
311         case 21:
312           targetParticle.SetDefinition( aSigmaZero );
313           break;
314         case 22:
315           targetParticle.SetDefinition( aSigmaMinus );
316           break;
317        }
318        targetHasChanged = true;
319        switch( ipaky2[i] )
320        {
321         case 10:
322           vec[i3]->SetDefinition( aKaonPlus ); 
323           vec[i3]->SetMayBeKilled(false);
324           break;
325         case 11:
326           vec[i3]->SetDefinition( aKaonZS );
327           vec[i3]->SetMayBeKilled(false);
328           break;
329         case 12:
330           vec[i3]->SetDefinition( aKaonZL );
331           vec[i3]->SetMayBeKilled(false);
332           break;
333        }
334      }
335      else  // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
336      {
337        // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
338        //              24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
339        if( (currentParticle.GetDefinition() == anAntiProton) ||
340            (currentParticle.GetDefinition() == anAntiNeutron) ||
341            (currentParticle.GetDefinition() == anAntiLambda) ||
342            (currentMass > sigmaMinusMass) )
343        {
344          switch( ipakyb1[i] )
345          {
346           case 19:
347             currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
348             break;
349           case 23:
350             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
351             break;
352           case 24:
353             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
354             break;
355           case 25:
356             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
357             break;
358          }
359          incidentHasChanged = true;
360          switch( ipakyb2[i] )
361          {
362           case 11:
363             vec[i3]->SetDefinition( aKaonZS ); 
364             vec[i3]->SetMayBeKilled(false);
365             break;
366           case 12:
367             vec[i3]->SetDefinition( aKaonZL );
368             vec[i3]->SetMayBeKilled(false);
369             break;
370           case 13:
371             vec[i3]->SetDefinition( aKaonMinus );
372             vec[i3]->SetMayBeKilled(false);
373             break;
374          }
375        }
376        else
377        {
378          switch( ipaky1[i] )
379          {
380           case 18:
381             currentParticle.SetDefinitionAndUpdateE( aLambda );
382             break;
383           case 20:
384             currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
385             break;
386           case 21:
387             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
388             break;
389           case 22:
390             currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
391             break;
392          }
393          incidentHasChanged = true;
394          switch( ipaky2[i] )
395          {
396           case 10:
397             vec[i3]->SetDefinition( aKaonPlus ); 
398             vec[i3]->SetMayBeKilled(false);
399             break;
400           case 11:
401             vec[i3]->SetDefinition( aKaonZS );
402             vec[i3]->SetMayBeKilled(false);
403             break;
404           case 12:
405             vec[i3]->SetDefinition( aKaonZL );
406             vec[i3]->SetMayBeKilled(false);
407             break;
408          }
409        }
410      }
411  }
412  else return true;
413
414  //
415  //  check the available energy
416  //   if there is not enough energy for kkb/ky pair production
417  //   then reduce the number of secondary particles
418  //  NOTE:
419  //        the number of secondaries may have been changed
420  //        the incident and/or target particles may have changed
421  //        charge conservation is ignored (as well as strangness conservation)
422  //
423  currentMass = currentParticle.GetMass()/GeV;
424  targetMass = targetParticle.GetMass()/GeV;
425   
426  G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
427  for( i=0; i<vecLen; ++i )
428  {
429    energyCheck -= vec[i]->GetMass()/GeV;
430    if( energyCheck < 0.0 )      // chop off the secondary List
431    {
432      vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
433      G4int j;
434      for(j=i; j<vecLen; j++) delete vec[j];
435      break;
436    }
437  }
438
439  return true;
440}
441
442 
443 /* end of file */
Note: See TracBrowser for help on using the repository browser.