source: trunk/source/processes/hadronic/models/high_energy/src/G4HEKaonZeroShortInelastic.cc @ 1340

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

update ti head
File size: 29.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// $Id: G4HEKaonZeroShortInelastic.cc,v 1.11 2010/02/09 21:58:57 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29//
30
31#include "globals.hh"
32#include "G4ios.hh"
33
34//
35// G4 Process: Gheisha High Energy Collision model.
36// This includes the high energy cascading model, the two-body-resonance model
37// and the low energy two-body model. Not included are the low energy stuff like
38// nuclear reactions, nuclear fission without any cascading and all processes for
39// particles at rest.
40// 
41// New version by D.H. Wright (SLAC) to fix seg fault in old version
42// 21 January 2010
43
44 
45#include "G4HEKaonZeroShortInelastic.hh"
46
47G4HadFinalState* G4HEKaonZeroShortInelastic::
48ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
49{
50  G4HEVector * pv = new G4HEVector[MAXPART];
51  const G4HadProjectile *aParticle = &aTrack;
52  const G4double atomicWeight = targetNucleus.GetN();
53  const G4double atomicNumber = targetNucleus.GetZ();
54  G4HEVector incidentParticle(aParticle);
55
56  G4int    incidentCode          = incidentParticle.getCode();
57  G4double incidentMass          = incidentParticle.getMass();
58  G4double incidentTotalEnergy   = incidentParticle.getEnergy();
59  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
60  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
61
62  if(incidentKineticEnergy < 1)
63    G4cout << "GHEKaonZeroShortInelastic: incident energy < 1 GeV" << G4endl;
64
65  if(verboseLevel > 1) {
66    G4cout << "G4HEKaonZeroShortInelastic::ApplyYourself" << G4endl;
67    G4cout << "incident particle " << incidentParticle.getName()
68           << "mass "              << incidentMass
69           << "kinetic energy "    << incidentKineticEnergy
70           << G4endl;
71    G4cout << "target material with (A,Z) = (" 
72           << atomicWeight << "," << atomicNumber << ")" << G4endl;
73  }
74   
75  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
76                                              atomicWeight, atomicNumber);
77  if(verboseLevel > 1)
78    G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
79   
80  incidentKineticEnergy -= inelasticity;
81   
82  G4double excitationEnergyGNP = 0.;
83  G4double excitationEnergyDTA = 0.; 
84
85  G4double excitation = NuclearExcitation(incidentKineticEnergy,
86                                          atomicWeight, atomicNumber,
87                                          excitationEnergyGNP,
88                                          excitationEnergyDTA);
89  if(verboseLevel > 1)
90    G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
91           << excitationEnergyDTA << G4endl;             
92
93  incidentKineticEnergy -= excitation;
94  incidentTotalEnergy    = incidentKineticEnergy + incidentMass;
95  incidentTotalMomentum  = std::sqrt( (incidentTotalEnergy-incidentMass)                   
96                                  *(incidentTotalEnergy+incidentMass));
97
98
99  G4HEVector targetParticle;
100  if(G4UniformRand() < atomicNumber/atomicWeight) { 
101    targetParticle.setDefinition("Proton");
102  } else { 
103    targetParticle.setDefinition("Neutron");
104  }
105
106  G4double targetMass = targetParticle.getMass();
107  G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
108                                       + 2.0*targetMass*incidentTotalEnergy);
109  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
110                                    // this was the meaning of inElastic in the
111                                    // original Gheisha stand-alone version.
112//  G4bool inElastic = InElasticCrossSectionInFirstInt
113//                      (availableEnergy, incidentCode, incidentTotalMomentum); 
114// for unknown reasons, it has been replaced by the following code in Geant???
115
116  G4bool inElastic = true;
117//    if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;   
118
119  vecLength = 0;           
120       
121  if(verboseLevel > 1)
122    G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
123           << incidentCode << G4endl;
124
125  G4bool successful = false; 
126   
127  if(inElastic || (!inElastic && atomicWeight < 1.5)) {
128 
129    // Split K0L into K0 and K0bar
130    if (G4UniformRand() < 0.5)
131      FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
132                                incidentParticle, targetParticle );
133    else
134      FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength,
135                            incidentParticle, targetParticle, atomicWeight );
136
137    // Do nuclear interaction with either K0 or K0bar
138    if ((vecLength > 0) && (availableEnergy > 1.)) 
139      StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
140                                    pv, vecLength,
141                                    incidentParticle, targetParticle);
142    HighEnergyCascading(successful, pv, vecLength,
143                        excitationEnergyGNP, excitationEnergyDTA,
144                        incidentParticle, targetParticle,
145                        atomicWeight, atomicNumber);
146    if (!successful)
147      HighEnergyClusterProduction(successful, pv, vecLength,
148                                  excitationEnergyGNP, excitationEnergyDTA,
149                                  incidentParticle, targetParticle,
150                                  atomicWeight, atomicNumber);
151    if (!successful) 
152      MediumEnergyCascading(successful, pv, vecLength, 
153                            excitationEnergyGNP, excitationEnergyDTA, 
154                            incidentParticle, targetParticle,
155                            atomicWeight, atomicNumber);
156
157    if (!successful)
158      MediumEnergyClusterProduction(successful, pv, vecLength,
159                                    excitationEnergyGNP, excitationEnergyDTA,       
160                                    incidentParticle, targetParticle,
161                                    atomicWeight, atomicNumber);
162    if (!successful)
163      QuasiElasticScattering(successful, pv, vecLength,
164                             excitationEnergyGNP, excitationEnergyDTA,
165                             incidentParticle, targetParticle, 
166                             atomicWeight, atomicNumber);
167  }
168
169  if (!successful) 
170    ElasticScattering(successful, pv, vecLength,
171                      incidentParticle,   
172                      atomicWeight, atomicNumber);
173
174  if (!successful) 
175    G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
176           << G4endl;
177
178  // Check for K0, K0bar and change particle types to K0L, K0S if necessary
179  G4int kcode;
180  for (G4int i = 0; i < vecLength; i++) {
181    kcode = pv[i].getCode();
182    if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) {
183      if (G4UniformRand() < 0.5) 
184        pv[i] = KaonZeroShort; 
185      else
186        pv[i] = KaonZeroLong;
187    }
188  } 
189
190  //      ................
191 
192  FillParticleChange(pv,  vecLength);
193  delete [] pv;
194  theParticleChange.SetStatusChange(stopAndKill);
195  return & theParticleChange;
196}
197
198
199void
200G4HEKaonZeroShortInelastic::FirstIntInCasKaonZero(G4bool &inElastic,
201                                                 const G4double availableEnergy,
202                                                 G4HEVector pv[],
203                                                 G4int &vecLen,
204                                                 G4HEVector incidentParticle,
205                                                 G4HEVector targetParticle,
206                                                 const G4double atomicWeight)
207
208// Kaon0 undergoes interaction with nucleon within a nucleus.  Check if it is
209// energetically possible to produce pions/kaons.  In not, assume nuclear excitation
210// occurs and input particle is degraded in energy. No other particles are produced.
211// If reaction is possible, find the correct number of pions/protons/neutrons
212// produced using an interpolation to multiplicity data.  Replace some pions or
213// protons/neutrons by kaons or strange baryons according to the average
214// multiplicity per inelastic reaction.
215
216{
217  static const G4double expxu =  std::log(MAXFLOAT); // upper bound for arg. of exp
218  static const G4double expxl = -expxu;         // lower bound for arg. of exp
219
220  static const G4double protb = 0.7;
221  static const G4double neutb = 0.7;
222  static const G4double     c = 1.25;
223
224  static const G4int   numMul = 1200;
225  static const G4int   numSec = 60;
226
227  G4int neutronCode = Neutron.getCode();
228  G4int protonCode  = Proton.getCode();
229
230  G4int targetCode = targetParticle.getCode();
231  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
232
233  static G4bool first = true;
234  static G4double protmul[numMul], protnorm[numSec];  // proton constants
235  static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
236
237  // misc. local variables
238  // np = number of pi+,  nm = number of pi-,  nz = number of pi0
239
240  G4int i, counter, nt, np, nm, nz;
241
242  if (first) {
243    // compute normalization constants, this will only be done once
244    first = false;
245    for( i=0; i<numMul; i++ )protmul[i]  = 0.0;
246    for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
247    counter = -1;
248    for (np=0; np<(numSec/3); np++) {
249      for (nm=std::max(0,np-1); nm<=(np+1); nm++) {
250        for (nz=0; nz<numSec/3; nz++) {
251          if (++counter < numMul) {
252            nt = np+nm+nz;
253            if( (nt>0) && (nt<=numSec) ) {
254              protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ;
255              protnorm[nt-1] += protmul[counter];
256            }
257          }
258        }
259      }
260    }
261
262    for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
263    for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
264    counter = -1;
265    for (np=0; np<numSec/3; np++) {
266      for (nm=np; nm<=(np+2); nm++) {
267        for (nz=0; nz<numSec/3; nz++) {
268          if (++counter < numMul) {
269            nt = np+nm+nz;
270            if( (nt>0) && (nt<=numSec) ) {
271              neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c);
272              neutnorm[nt-1] += neutmul[counter];
273            }
274          }
275        }
276      }
277    }
278
279    for (i=0; i<numSec; i++) {
280      if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
281      if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
282    }
283  }    // end of initialization
284
285
286  // Initialize the first two particles
287  // the same as beam and target
288  pv[0] = incidentParticle;
289  pv[1] = targetParticle;
290  vecLen = 2;
291
292  if( !inElastic ) {
293    // quasi-elastic scattering, no pions produced
294    if( targetCode == protonCode) {
295      G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
296      G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) );
297      if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) {
298        // charge exchange  K+ n -> K0 p
299        pv[0] = KaonPlus;
300        pv[1] = Neutron;
301      }
302    }
303    return;
304  } else if (availableEnergy <= PionPlus.getMass()) {
305    return;
306  }
307
308  // Inelastic scattering
309
310  np = 0, nm = 0, nz = 0;
311  G4double eab = availableEnergy;
312  G4int ieab = G4int( eab*5.0 );
313
314  G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
315  if( (ieab <= 9) && (G4UniformRand() >=  supp[ieab])) {
316    // Suppress high multiplicity events at low momentum
317    // only one additional pion will be produced
318    G4double w0, wp, wm, wt, ran;
319    if (targetCode == neutronCode) {
320      // target is a neutron
321      w0 = - sqr(1.+protb)/(2.*c*c);
322      w0 = std::exp(w0);
323      wm = - sqr(-1.+protb)/(2.*c*c);
324      wm = std::exp(wm);
325      w0 = w0/2.;
326      wm = wm*1.5;
327      if (G4UniformRand() < w0/(w0+wm) ) {
328        np = 0;
329        nm = 0;
330        nz = 1;
331      } else {
332        np = 0;
333        nm = 1;
334        nz = 0;
335      }
336
337    } else {
338      // target is a proton
339      w0 = -sqr(1.+neutb)/(2.*c*c);
340      wp = w0 = std::exp(w0);
341      wm = -sqr(-1.+neutb)/(2.*c*c);
342      wm = std::exp(wm);
343      wt = w0+wp+wm;
344      wp = w0+wp;
345      ran = G4UniformRand();
346      if ( ran < w0/wt) {
347        np = 0;
348        nm = 0;
349        nz = 1;
350      } else if (ran < wp/wt) {
351        np = 1;
352        nm = 0;
353        nz = 0;
354      } else {
355        np = 0;
356        nm = 1;
357        nz = 0;
358      }
359    }
360  } else {
361    // number of total particles vs. centre of mass Energy - 2*proton mass
362
363    G4double aleab = std::log(availableEnergy);
364    G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
365                 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
366
367    // Normalization constant for kno-distribution.
368    // Calculate first the sum of all constants, check for numerical problems.
369    G4double test, dum, anpn = 0.0;
370
371       for (nt=1; nt<=numSec; nt++) {
372         test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
373         dum = pi*nt/(2.0*n*n);
374         if (std::fabs(dum) < 1.0) {
375           if( test >= 1.0e-10 )anpn += dum*test;
376         } else {
377           anpn += dum*test;
378         }
379       }
380
381       G4double ran = G4UniformRand();
382       G4double excs = 0.0;
383       if( targetCode == protonCode )
384         {
385           counter = -1;
386           for( np=0; np<numSec/3; np++ )
387              {
388                for( nm=std::max(0,np-1); nm<=(np+1); nm++ )
389                   {
390                     for (nz=0; nz<numSec/3; nz++) {
391                       if (++counter < numMul) {
392                         nt = np+nm+nz;
393                         if ( (nt>0) && (nt<=numSec) ) {
394                           test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
395                           dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
396                           if (std::fabs(dum) < 1.0) {
397                             if( test >= 1.0e-10 )excs += dum*test;
398                           } else {
399                             excs += dum*test;
400                           }
401                           if (ran < excs) goto outOfLoop;      //----------------------->
402                         }
403                       }
404                     }
405                   }
406              }
407
408                                     // 3 previous loops continued to the end
409           inElastic = false;                 // quasi-elastic scattering
410           return;
411         }
412       else
413         {                                         // target must be a neutron
414           counter = -1;
415           for( np=0; np<numSec/3; np++ )
416              {
417                for( nm=np; nm<=(np+2); nm++ )
418                   {
419                     for (nz=0; nz<numSec/3; nz++) {
420                       if (++counter < numMul) {
421                         nt = np+nm+nz;
422                         if ( (nt>=1) && (nt<=numSec) ) {
423                           test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
424                           dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
425                           if (std::fabs(dum) < 1.0) {
426                             if( test >= 1.0e-10 )excs += dum*test;
427                           } else {
428                             excs += dum*test;
429                           }
430                           if (ran < excs) goto outOfLoop;       // -------------------------->
431                         }
432                       }
433                     }
434                   }
435              }
436                                                  // 3 previous loops continued to the end
437           inElastic = false;                     // quasi-elastic scattering.
438           return;
439         }
440     }
441   outOfLoop:           //  <-----------------------------------------------
442
443   if( targetCode == neutronCode)
444     {
445       if( np == nm)
446         {
447         }
448       else if (np == (nm-1))
449         {
450           if( G4UniformRand() < 0.5)
451             {
452               pv[0] = KaonPlus;
453             }
454           else
455             {
456               pv[1] = Proton;
457             }
458         }
459       else
460         {
461           pv[0] = KaonPlus;
462           pv[1] = Proton;
463         }
464     }
465   else
466     {
467       if( np == nm )
468         {
469           if( G4UniformRand() < 0.25)
470             {
471               pv[0] = KaonPlus;
472               pv[1] = Neutron;
473             }
474           else
475             {
476             }
477         }
478       else if ( np == (nm+1))
479         {
480           pv[1] = Neutron;
481         }
482       else
483         {
484           pv[0] = KaonPlus;
485         }
486     }
487
488  nt = np + nm + nz;
489  while (nt > 0) {
490    G4double ran = G4UniformRand();
491    if (ran < (G4double)np/nt) {
492      if (np > 0) {
493        pv[vecLen++] = PionPlus;
494        np--;
495      }
496    } else if ( ran < (G4double)(np+nm)/nt) {
497      if (nm > 0) {
498        pv[vecLen++] = PionMinus;
499        nm--;
500      }
501    } else {
502      if (nz > 0) {
503        pv[vecLen++] = PionZero;
504        nz--;
505      }
506    }
507    nt = np + nm + nz;
508  }
509
510  if (verboseLevel > 1) {
511    G4cout << "Particles produced: " ;
512    G4cout << pv[0].getName() << " " ;
513    G4cout << pv[1].getName() << " " ;
514    for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
515    G4cout << G4endl;
516  }
517
518  return;
519}
520
521
522void
523G4HEKaonZeroShortInelastic::FirstIntInCasAntiKaonZero(G4bool &inElastic,
524                                                     const G4double availableEnergy,
525                                                     G4HEVector pv[],
526                                                     G4int &vecLen,
527                                                     G4HEVector incidentParticle,
528                                                     G4HEVector targetParticle )
529
530// AntiKaon0 undergoes interaction with nucleon within a nucleus.  Check if it is
531// energetically possible to produce pions/kaons.  In not, assume nuclear excitation
532// occurs and input particle is degraded in energy. No other particles are produced.
533// If reaction is possible, find the correct number of pions/protons/neutrons
534// produced using an interpolation to multiplicity data.  Replace some pions or
535// protons/neutrons by kaons or strange baryons according to the average
536// multiplicity per inelastic reaction.
537
538{
539  static const G4double expxu =  std::log(MAXFLOAT); // upper bound for arg. of exp
540  static const G4double expxl = -expxu;         // lower bound for arg. of exp
541
542  static const G4double protb = 0.7;
543  static const G4double neutb = 0.7;
544  static const G4double     c = 1.25;
545
546  static const G4int   numMul = 1200;
547  static const G4int   numSec = 60;
548
549  G4int neutronCode = Neutron.getCode();
550  G4int protonCode = Proton.getCode();
551  G4int kaonMinusCode = KaonMinus.getCode();
552  G4int kaonZeroCode = KaonZero.getCode();
553  G4int antiKaonZeroCode = AntiKaonZero.getCode(); 
554
555  G4int targetCode = targetParticle.getCode();
556  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
557
558  static G4bool first = true;
559  static G4double protmul[numMul], protnorm[numSec];  // proton constants
560  static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
561
562//  misc. local variables
563//  np = number of pi+,  nm = number of pi-,  nz = number of pi0
564
565  G4int i, counter, nt, np, nm, nz;
566
567  if(first) {
568    // compute normalization constants, this will only be done once
569    first = false;
570    for( i=0; i<numMul; i++ )protmul[i]  = 0.0;
571    for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
572    counter = -1;
573    for(np=0; np<(numSec/3); np++) {
574      for(nm=std::max(0,np-2); nm<=np; nm++) {
575        for(nz=0; nz<numSec/3; nz++) {
576          if(++counter < numMul) {
577            nt = np+nm+nz;
578            if( (nt>0) && (nt<=numSec) ) {
579              protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ;
580              protnorm[nt-1] += protmul[counter];
581            }
582          }
583        }
584      }
585    }
586
587    for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
588    for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
589    counter = -1;
590    for(np=0; np<numSec/3; np++) {
591      for(nm=std::max(0,np-1); nm<=(np+1); nm++) {
592        for(nz=0; nz<numSec/3; nz++) {
593          if(++counter < numMul) {
594            nt = np+nm+nz;
595            if( (nt>0) && (nt<=numSec) ) {
596              neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c);
597              neutnorm[nt-1] += neutmul[counter];
598            }
599          }
600        }
601      }
602    }
603
604    for(i=0; i<numSec; i++) {
605      if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
606      if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
607    }
608  }                                // end of initialization
609
610  // initialize the first two particles
611  // the same as beam and target                                   
612  pv[0] = incidentParticle;
613  pv[1] = targetParticle;
614  vecLen = 2;
615
616  if (!inElastic || (availableEnergy <= PionPlus.getMass())) 
617    return;
618                                       
619  // Inelastic scattering
620
621  np = 0, nm = 0, nz = 0;
622  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
623  G4int iplab = G4int( incidentTotalMomentum*5.);
624  if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
625    G4int     iplab = std::min(19, G4int( incidentTotalMomentum*5.));
626    G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
627                       0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
628    if(G4UniformRand() < cnk0[iplab]) {
629      if(targetCode == protonCode) {
630        return;
631      } else {
632        pv[0] = KaonMinus;
633        pv[1] = Proton;
634        return;
635      }
636    }
637
638    G4double ran = G4UniformRand();
639    if(targetCode == protonCode) {
640
641      // target is a proton
642      if( ran < 0.25 ) {
643        ; 
644      } else if (ran < 0.50) {
645        pv[0] = PionPlus;
646        pv[1] = SigmaZero;
647      } else if (ran < 0.75) {
648        ;
649      } else {
650        pv[0] = PionPlus;
651        pv[1] = Lambda;
652      }
653    } else {
654
655      // target is a neutron
656      if( ran < 0.25 ) { 
657        pv[0] = PionMinus;
658        pv[1] = SigmaPlus;
659      } else if (ran < 0.50) {
660        pv[0] = PionZero;
661        pv[1] = SigmaZero;
662      } else if (ran < 0.75) { 
663        pv[0] = PionPlus;
664        pv[1] = SigmaMinus;
665      } else {
666        pv[0] = PionZero;
667        pv[1] = Lambda;
668      }
669    }
670    return;
671
672  } else {
673    // number of total particles vs. centre of mass Energy - 2*proton mass
674   
675    G4double aleab = std::log(availableEnergy);
676    G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
677                 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
678   
679    // Normalization constant for kno-distribution.
680    // Calculate first the sum of all constants, check for numerical problems.   
681    G4double test, dum, anpn = 0.0;
682
683    for (nt=1; nt<=numSec; nt++) {
684      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
685      dum = pi*nt/(2.0*n*n);
686      if (std::fabs(dum) < 1.0) { 
687        if( test >= 1.0e-10 )anpn += dum*test;
688      } else { 
689        anpn += dum*test;
690      }
691    }
692   
693    G4double ran = G4UniformRand();
694    G4double excs = 0.0;
695    if (targetCode == protonCode) {
696      counter = -1;
697      for (np=0; np<numSec/3; np++) {
698        for (nm=std::max(0,np-2); nm<=np; nm++) {
699          for (nz=0; nz<numSec/3; nz++) {
700            if (++counter < numMul) {
701              nt = np+nm+nz;
702              if( (nt>0) && (nt<=numSec) ) {
703                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
704                dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
705
706                if (std::fabs(dum) < 1.0) { 
707                  if( test >= 1.0e-10 )excs += dum*test;
708                } else { 
709                  excs += dum*test;
710                }
711
712                if (ran < excs) goto outOfLoop;      //----------------------->
713              }
714            }
715          }
716        }
717      }
718                            // 3 previous loops continued to the end
719      inElastic = false;    // quasi-elastic scattering   
720      return;
721
722    } else {         // target must be a neutron
723      counter = -1;
724      for (np=0; np<numSec/3; np++) {
725        for (nm=std::max(0,np-1); nm<=(np+1); nm++) {
726          for (nz=0; nz<numSec/3; nz++) {
727            if (++counter < numMul) {
728              nt = np+nm+nz;
729              if( (nt>=1) && (nt<=numSec) ) {
730                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
731                dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
732
733                if (std::fabs(dum) < 1.0) { 
734                  if( test >= 1.0e-10 )excs += dum*test;
735                } else { 
736                  excs += dum*test;
737                }
738
739                if (ran < excs) goto outOfLoop;   // -------------------------->
740              }
741            }
742          }
743        }
744      }
745                              // 3 previous loops continued to the end
746      inElastic = false;      // quasi-elastic scattering.
747      return;
748    }
749  } 
750  outOfLoop:   //  <------------------------------------------------------------------------   
751   
752  if( targetCode == protonCode)
753     {
754       if( np == nm)
755         {
756         }
757       else if (np == (1+nm))
758         {
759           if( G4UniformRand() < 0.5)
760             {
761               pv[0] = KaonMinus;
762             }
763           else
764             {
765               pv[1] = Neutron;
766             }
767         }
768       else     
769         {
770           pv[0] = KaonMinus;
771           pv[1] = Neutron;
772         } 
773     } 
774  else
775     {
776       if( np == nm)
777         {
778           if( G4UniformRand() < 0.75)
779             {
780             }
781           else
782             {
783               pv[0] = KaonMinus;
784               pv[1] = Proton;
785             }
786         } 
787       else if ( np == (1+nm))
788         {
789           pv[0] = KaonMinus;
790         }
791       else
792         {
793           pv[1] = Proton;
794         }
795     }     
796
797
798  if( G4UniformRand() < 0.5 )   
799     {
800       if(    (    (pv[0].getCode() == kaonMinusCode)
801                && (pv[1].getCode() == neutronCode)  )
802           || (    (pv[0].getCode() == kaonZeroCode)
803                && (pv[1].getCode() == protonCode)   )
804           || (    (pv[0].getCode() == antiKaonZeroCode)
805                && (pv[1].getCode() == protonCode)   )   )
806         {
807           G4double ran = G4UniformRand();
808           if( pv[1].getCode() == protonCode)
809             { 
810               if(ran < 0.68)
811                 {
812                   pv[0] = PionPlus;
813                   pv[1] = Lambda;
814                 }
815               else if (ran < 0.84)
816                 {
817                   pv[0] = PionZero;
818                   pv[1] = SigmaPlus;
819                 }
820               else
821                 {
822                   pv[0] = PionPlus;
823                   pv[1] = SigmaZero;
824                 }
825             }
826           else
827             {
828               if(ran < 0.68)
829                 {
830                   pv[0] = PionMinus;
831                   pv[1] = Lambda;
832                 }
833               else if (ran < 0.84)
834                 {
835                   pv[0] = PionMinus;
836                   pv[1] = SigmaZero;
837                 }
838               else
839                 {
840                   pv[0] = PionZero;
841                   pv[1] = SigmaMinus;
842                 }
843             }
844         } 
845       else
846         {
847           G4double ran = G4UniformRand();
848           if (ran < 0.67)
849              {
850                pv[0] = PionZero;
851                pv[1] = Lambda;
852              }
853           else if (ran < 0.78)
854              {
855                pv[0] = PionMinus;
856                pv[1] = SigmaPlus;
857              }
858           else if (ran < 0.89)
859              {
860                pv[0] = PionZero;
861                pv[1] = SigmaZero;
862              }
863           else
864              {
865                pv[0] = PionPlus;
866                pv[1] = SigmaMinus;
867              }
868         }
869    }
870
871  nt = np + nm + nz;
872  while ( nt > 0) {
873    G4double ran = G4UniformRand();
874    if ( ran < (G4double)np/nt) { 
875      if( np > 0 ) {
876        pv[vecLen++] = PionPlus;
877        np--;
878      }
879    } else if (ran < (G4double)(np+nm)/nt) {   
880      if( nm > 0 ) { 
881        pv[vecLen++] = PionMinus;
882        nm--;
883      }
884    } else {
885      if( nz > 0 ) { 
886        pv[vecLen++] = PionZero;
887        nz--;
888      }
889    }
890    nt = np + nm + nz;
891  }
892 
893  if (verboseLevel > 1) {
894    G4cout << "Particles produced: " ;
895    G4cout << pv[0].getName() << " " ;
896    G4cout << pv[1].getName() << " " ;
897    for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
898    G4cout << G4endl;
899  }
900
901  return;
902}
Note: See TracBrowser for help on using the repository browser.