Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/rpg/src/G4RPGNeutronInelastic.cc

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4RPGNeutronInelastic.cc,v 1.1 2007/07/18 21:04:20 dennis Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4RPGNeutronInelastic.cc,v 1.4 2008/05/05 21:21:55 dennis Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929 
    3030#include "G4RPGNeutronInelastic.hh"
    3131#include "Randomize.hh"
    32 #include "G4Electron.hh"
    33 // #include "DumpFrame.hh"
    34 
    35  G4HadFinalState *
    36   G4RPGNeutronInelastic::ApplyYourself( const G4HadProjectile &aTrack,
    37                                        G4Nucleus &targetNucleus )
    38   {
    39     theParticleChange.Clear();
    40     const G4HadProjectile *originalIncident = &aTrack;
    41     //
    42     // create the target particle
    43     //
    44     G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
    45    
    46     if( verboseLevel > 1 )
    47     {
    48       const G4Material *targetMaterial = aTrack.GetMaterial();
    49       G4cout << "G4RPGNeutronInelastic::ApplyYourself called" << G4endl;
    50       G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
    51       G4cout << "target material = " << targetMaterial->GetName() << ", ";
    52       G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
    53            << G4endl;
    54     }
    55 /* not true, for example for Fe56, etc..
    56     if( originalIncident->GetKineticEnergy()/MeV < 0.000001 )
    57       throw G4HadronicException(__FILE__, __LINE__, "G4RPGNeutronInelastic: should be capture process!");
    58     if( originalIncident->Get4Momentum().vect().mag()/MeV < 0.000001 )
    59       throw G4HadronicException(__FILE__, __LINE__, "G4RPGNeutronInelastic: should be capture process!");
    60 */
    61    
    62     G4ReactionProduct modifiedOriginal;
    63     modifiedOriginal = *originalIncident;
    64     G4ReactionProduct targetParticle;
    65     targetParticle = *originalTarget;
    66     if( originalIncident->GetKineticEnergy()/GeV < 0.01 + 2.*G4UniformRand()/9. )
    67     {
    68       SlowNeutron( originalIncident, modifiedOriginal, targetParticle, targetNucleus );
    69       delete originalTarget;
    70       return &theParticleChange;
    71     }
    72     //
    73     // Fermi motion and evaporation
    74     // As of Geant3, the Fermi energy calculation had not been Done
    75     //
    76     G4double ek = originalIncident->GetKineticEnergy()/MeV;
    77     G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
    78    
    79     G4double tkin = targetNucleus.Cinema( ek );
    80     ek += tkin;
    81     modifiedOriginal.SetKineticEnergy( ek*MeV );
    82     G4double et = ek + amas;
    83     G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
    84     G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
    85     if( pp > 0.0 )
    86     {
    87       G4ThreeVector momentum = modifiedOriginal.GetMomentum();
    88       modifiedOriginal.SetMomentum( momentum * (p/pp) );
    89     }
    90     //
    91     // calculate black track energies
    92     //
    93     tkin = targetNucleus.EvaporationEffects( ek );
    94     ek -= tkin;
    95     modifiedOriginal.SetKineticEnergy( ek*MeV );
    96     et = ek + amas;
    97     p = std::sqrt( std::abs((et-amas)*(et+amas)) );
    98     pp = modifiedOriginal.GetMomentum().mag()/MeV;
    99     if( pp > 0.0 )
    100     {
    101       G4ThreeVector momentum = modifiedOriginal.GetMomentum();
    102       modifiedOriginal.SetMomentum( momentum * (p/pp) );
    103     }
    104     const G4double cutOff = 0.1;
    105     if( modifiedOriginal.GetKineticEnergy()/MeV <= cutOff )
    106     {
    107       SlowNeutron( originalIncident, modifiedOriginal, targetParticle, targetNucleus );
    108       delete originalTarget;
    109       return &theParticleChange;
    110     }
    111     G4ReactionProduct currentParticle = modifiedOriginal;
    112     currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
    113     targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
    114     G4bool incidentHasChanged = false;
    115     G4bool targetHasChanged = false;
    116     G4bool quasiElastic = false;
    117     G4FastVector<G4ReactionProduct,256> vec;  // vec will contain the secondary particles
    118     G4int vecLen = 0;
    119     vec.Initialize( 0 );
    120    
    121     Cascade( vec, vecLen,
    122              originalIncident, currentParticle, targetParticle,
    123              incidentHasChanged, targetHasChanged, quasiElastic );
    124    
    125     CalculateMomenta( vec, vecLen,
    126                       originalIncident, originalTarget, modifiedOriginal,
    127                       targetNucleus, currentParticle, targetParticle,
    128                       incidentHasChanged, targetHasChanged, quasiElastic );
    129    
    130     SetUpChange( vec, vecLen,
    131                  currentParticle, targetParticle,
    132                  incidentHasChanged );
    133    
     32
     33
     34G4HadFinalState*
     35G4RPGNeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
     36                                      G4Nucleus& targetNucleus)
     37{
     38  theParticleChange.Clear();
     39  const G4HadProjectile* originalIncident = &aTrack;
     40
     41  //
     42  // create the target particle
     43  //
     44  G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
     45 
     46  G4ReactionProduct modifiedOriginal;
     47  modifiedOriginal = *originalIncident;
     48  G4ReactionProduct targetParticle;
     49  targetParticle = *originalTarget;
     50  if( originalIncident->GetKineticEnergy()/GeV < 0.01 + 2.*G4UniformRand()/9. )
     51  {
     52    SlowNeutron(originalIncident,modifiedOriginal,targetParticle,targetNucleus );
    13453    delete originalTarget;
    13554    return &theParticleChange;
    13655  }
    137  
    138  void
    139   G4RPGNeutronInelastic::SlowNeutron(
    140    const G4HadProjectile *originalIncident,
    141    G4ReactionProduct &modifiedOriginal,
    142    G4ReactionProduct &targetParticle,
    143    G4Nucleus &targetNucleus )
    144   {       
    145     const G4double A = targetNucleus.GetN();    // atomic weight
    146     const G4double Z = targetNucleus.GetZ();    // atomic number
    147    
    148     G4double currentKinetic = modifiedOriginal.GetKineticEnergy()/MeV;
    149     G4double currentMass = modifiedOriginal.GetMass()/MeV;
    150     if( A < 1.5 )   // Hydrogen
     56
     57  //
     58  // Fermi motion and evaporation
     59  // As of Geant3, the Fermi energy calculation had not been Done
     60  //
     61  G4double ek = originalIncident->GetKineticEnergy()/MeV;
     62  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
     63   
     64  G4double tkin = targetNucleus.Cinema( ek );
     65  ek += tkin;
     66  modifiedOriginal.SetKineticEnergy( ek*MeV );
     67  G4double et = ek + amas;
     68  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
     69  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
     70  if( pp > 0.0 )
     71  {
     72    G4ThreeVector momentum = modifiedOriginal.GetMomentum();
     73    modifiedOriginal.SetMomentum( momentum * (p/pp) );
     74  }
     75  //
     76  // calculate black track energies
     77  //
     78  tkin = targetNucleus.EvaporationEffects( ek );
     79  ek -= tkin;
     80  modifiedOriginal.SetKineticEnergy(ek);
     81  et = ek + amas;
     82  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
     83  pp = modifiedOriginal.GetMomentum().mag();
     84  if( pp > 0.0 )
     85  {
     86    G4ThreeVector momentum = modifiedOriginal.GetMomentum();
     87    modifiedOriginal.SetMomentum( momentum * (p/pp) );
     88  }
     89  const G4double cutOff = 0.1;
     90  if( modifiedOriginal.GetKineticEnergy()/MeV <= cutOff )
     91  {
     92    SlowNeutron( originalIncident, modifiedOriginal, targetParticle, targetNucleus );
     93    delete originalTarget;
     94    return &theParticleChange;
     95  }
     96
     97  G4ReactionProduct currentParticle = modifiedOriginal;
     98  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
     99  targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
     100  G4bool incidentHasChanged = false;
     101  G4bool targetHasChanged = false;
     102  G4bool quasiElastic = false;
     103  G4FastVector<G4ReactionProduct,256> vec;  // vec will contain sec. particles
     104  G4int vecLen = 0;
     105  vec.Initialize( 0 );
     106   
     107  InitialCollision(vec, vecLen, currentParticle, targetParticle,
     108                   incidentHasChanged, targetHasChanged);
     109   
     110  CalculateMomenta(vec, vecLen,
     111                   originalIncident, originalTarget, modifiedOriginal,
     112                   targetNucleus, currentParticle, targetParticle,
     113                   incidentHasChanged, targetHasChanged, quasiElastic);
     114   
     115  SetUpChange(vec, vecLen,
     116              currentParticle, targetParticle,
     117              incidentHasChanged);
     118   
     119  delete originalTarget;
     120  return &theParticleChange;
     121}
     122 
     123void
     124G4RPGNeutronInelastic::SlowNeutron(const G4HadProjectile* originalIncident,
     125                              G4ReactionProduct& modifiedOriginal,
     126                              G4ReactionProduct& targetParticle,
     127                              G4Nucleus& targetNucleus)
     128{       
     129  const G4double A = targetNucleus.GetN();    // atomic weight
     130  const G4double Z = targetNucleus.GetZ();    // atomic number
     131   
     132  G4double currentKinetic = modifiedOriginal.GetKineticEnergy()/MeV;
     133  G4double currentMass = modifiedOriginal.GetMass()/MeV;
     134  if( A < 1.5 )   // Hydrogen
     135  {
     136    //
     137    // very simple simulation of scattering angle and energy
     138    // nonrelativistic approximation with isotropic angular
     139    // distribution in the cms system
     140    //
     141    G4double cost1, eka = 0.0;
     142    while (eka <= 0.0)
    151143    {
    152       //
    153       // very simple simulation of scattering angle and energy
    154       // nonrelativistic approximation with isotropic angular
    155       // distribution in the cms system
    156       //
    157       G4double cost1, eka = 0.0;
    158       while (eka <= 0.0)
    159       {
    160         cost1 = -1.0 + 2.0*G4UniformRand();
    161         eka = 1.0 + 2.0*cost1*A + A*A;
     144      cost1 = -1.0 + 2.0*G4UniformRand();
     145      eka = 1.0 + 2.0*cost1*A + A*A;
     146    }
     147    G4double cost = std::min( 1.0, std::max( -1.0, (A*cost1+1.0)/std::sqrt(eka) ) );
     148    eka /= (1.0+A)*(1.0+A);
     149    G4double ek = currentKinetic*MeV/GeV;
     150    G4double amas = currentMass*MeV/GeV;
     151    ek *= eka;
     152    G4double en = ek + amas;
     153    G4double p = std::sqrt(std::abs(en*en-amas*amas));
     154    G4double sint = std::sqrt(std::abs(1.0-cost*cost));
     155    G4double phi = G4UniformRand()*twopi;
     156    G4double px = sint*std::sin(phi);
     157    G4double py = sint*std::cos(phi);
     158    G4double pz = cost;
     159    targetParticle.SetMomentum( px*GeV, py*GeV, pz*GeV );
     160    G4double pxO = originalIncident->Get4Momentum().x()/GeV;
     161    G4double pyO = originalIncident->Get4Momentum().y()/GeV;
     162    G4double pzO = originalIncident->Get4Momentum().z()/GeV;
     163    G4double ptO = pxO*pxO + pyO+pyO;
     164    if( ptO > 0.0 )
     165    {
     166      G4double pO = std::sqrt(pxO*pxO+pyO*pyO+pzO*pzO);
     167      cost = pzO/pO;
     168      sint = 0.5*(std::sqrt(std::abs((1.0-cost)*(1.0+cost)))+std::sqrt(ptO)/pO);
     169      G4double ph = pi/2.0;
     170      if( pyO < 0.0 )ph = ph*1.5;
     171      if( std::abs(pxO) > 0.000001 )ph = std::atan2(pyO,pxO);
     172      G4double cosp = std::cos(ph);
     173      G4double sinp = std::sin(ph);
     174      px = cost*cosp*px - sinp*py+sint*cosp*pz;
     175      py = cost*sinp*px + cosp*py+sint*sinp*pz;
     176      pz = -sint*px     + cost*pz;
     177    }
     178    else
     179    {
     180      if( pz < 0.0 )pz *= -1.0;
     181    }
     182    G4double pu = std::sqrt(px*px+py*py+pz*pz);
     183    modifiedOriginal.SetMomentum( targetParticle.GetMomentum() * (p/pu) );
     184    modifiedOriginal.SetKineticEnergy( ek*GeV );
     185     
     186    targetParticle.SetMomentum(
     187     originalIncident->Get4Momentum().vect() - modifiedOriginal.GetMomentum() );
     188    G4double pp = targetParticle.GetMomentum().mag();
     189    G4double tarmas = targetParticle.GetMass();
     190    targetParticle.SetTotalEnergy( std::sqrt( pp*pp + tarmas*tarmas ) );
     191     
     192    theParticleChange.SetEnergyChange( modifiedOriginal.GetKineticEnergy() );
     193    G4DynamicParticle *pd = new G4DynamicParticle;
     194    pd->SetDefinition( targetParticle.GetDefinition() );
     195    pd->SetMomentum( targetParticle.GetMomentum() );
     196    theParticleChange.AddSecondary( pd );
     197    return;
     198  }
     199
     200  G4FastVector<G4ReactionProduct,4> vec;  // vec will contain the secondary particles
     201  G4int vecLen = 0;
     202  vec.Initialize( 0 );
     203   
     204  G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
     205  G4double massVec[9];
     206  massVec[0] = targetNucleus.AtomicMass( A+1.0, Z     );
     207  massVec[1] = theAtomicMass;
     208  massVec[2] = 0.;
     209  if (Z > 1.0) massVec[2] = targetNucleus.AtomicMass(A, Z-1.0);
     210  massVec[3] = 0.;
     211  if (Z > 1.0 && A > 1.0) massVec[3] = targetNucleus.AtomicMass(A-1.0, Z-1.0 );
     212  massVec[4] = 0.;
     213  if (Z > 1.0 && A > 2.0 && A-2.0 > Z-1.0)
     214    massVec[4] = targetNucleus.AtomicMass( A-2.0, Z-1.0 );
     215  massVec[5] = 0.;
     216  if (Z > 2.0 && A > 3.0 && A-3.0 > Z-2.0)
     217    massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-2.0 );
     218  massVec[6] = 0.;
     219  if (A > 1.0 && A-1.0 > Z) massVec[6] = targetNucleus.AtomicMass(A-1.0, Z);
     220  massVec[7] = massVec[3];
     221  massVec[8] = 0.;
     222  if (Z > 2.0 && A > 1.0) massVec[8] = targetNucleus.AtomicMass( A-1.0,Z-2.0 );
     223   
     224  twoBody.NuclearReaction(vec, vecLen, originalIncident,
     225                          targetNucleus, theAtomicMass, massVec );
     226   
     227  theParticleChange.SetStatusChange( stopAndKill );
     228  theParticleChange.SetEnergyChange( 0.0 );
     229   
     230  G4DynamicParticle* pd;
     231  for( G4int i=0; i<vecLen; ++i ) {
     232    pd = new G4DynamicParticle();
     233    pd->SetDefinition( vec[i]->GetDefinition() );
     234    pd->SetMomentum( vec[i]->GetMomentum() );
     235    theParticleChange.AddSecondary( pd );
     236    delete vec[i];
     237  }
     238}
     239
     240
     241// Initial Collision
     242//   selects the particle types arising from the initial collision of
     243//   the neutron and target nucleon.  Secondaries are assigned to
     244//   forward and backward reaction hemispheres, but final state energies
     245//   and momenta are not calculated here.
     246
     247void
     248G4RPGNeutronInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
     249                                   G4int& vecLen,
     250                                   G4ReactionProduct& currentParticle,
     251                                   G4ReactionProduct& targetParticle,
     252                                   G4bool& incidentHasChanged,
     253                                   G4bool& targetHasChanged)
     254{
     255  G4double KE = currentParticle.GetKineticEnergy()/GeV;
     256 
     257  G4int mult;
     258  G4int partType;
     259  std::vector<G4int> fsTypes;
     260  G4int part1;
     261  G4int part2;
     262 
     263  G4double testCharge;
     264  G4double testBaryon;
     265  G4double testStrange;
     266 
     267  // Get particle types according to incident and target types
     268
     269  if (targetParticle.GetDefinition() == particleDef[neu]) {
     270    mult = GetMultiplicityT1(KE);
     271    fsTypes = GetFSPartTypesForNN(mult, KE);
     272
     273    part1 = fsTypes[0];
     274    part2 = fsTypes[1];
     275    currentParticle.SetDefinition(particleDef[part1]);
     276    targetParticle.SetDefinition(particleDef[part2]);
     277    if (part1 == pro) {     
     278      if (part2 == neu) {
     279        if (G4UniformRand() > 0.5) {
     280          incidentHasChanged = true;
     281        } else {
     282          targetHasChanged = true;
     283          currentParticle.SetDefinition(particleDef[part2]);
     284          targetParticle.SetDefinition(particleDef[part1]);
     285        }
     286      } else {
     287        targetHasChanged = true;
     288        incidentHasChanged = true;
    162289      }
    163       G4double cost = std::min( 1.0, std::max( -1.0, (A*cost1+1.0)/std::sqrt(eka) ) );
    164       eka /= (1.0+A)*(1.0+A);
    165       G4double ek = currentKinetic*MeV/GeV;
    166       G4double amas = currentMass*MeV/GeV;
    167       ek *= eka;
    168       G4double en = ek + amas;
    169       G4double p = std::sqrt(std::abs(en*en-amas*amas));
    170       G4double sint = std::sqrt(std::abs(1.0-cost*cost));
    171       G4double phi = G4UniformRand()*twopi;
    172       G4double px = sint*std::sin(phi);
    173       G4double py = sint*std::cos(phi);
    174       G4double pz = cost;
    175       targetParticle.SetMomentum( px*GeV, py*GeV, pz*GeV );
    176       G4double pxO = originalIncident->Get4Momentum().x()/GeV;
    177       G4double pyO = originalIncident->Get4Momentum().y()/GeV;
    178       G4double pzO = originalIncident->Get4Momentum().z()/GeV;
    179       G4double ptO = pxO*pxO + pyO+pyO;
    180       if( ptO > 0.0 )
    181       {
    182         G4double pO = std::sqrt(pxO*pxO+pyO*pyO+pzO*pzO);
    183         cost = pzO/pO;
    184         sint = 0.5*(std::sqrt(std::abs((1.0-cost)*(1.0+cost)))+std::sqrt(ptO)/pO);
    185         G4double ph = pi/2.0;
    186         if( pyO < 0.0 )ph = ph*1.5;
    187         if( std::abs(pxO) > 0.000001 )ph = std::atan2(pyO,pxO);
    188         G4double cosp = std::cos(ph);
    189         G4double sinp = std::sin(ph);
    190         px = cost*cosp*px - sinp*py+sint*cosp*pz;
    191         py = cost*sinp*px + cosp*py+sint*sinp*pz;
    192         pz = -sint*px     + cost*pz;
     290
     291    } else { // neutron
     292      if (part2 > neu && part2 < xi0) targetHasChanged = true;
     293    }
     294
     295    testCharge = 0.0;
     296    testBaryon = 2.0;
     297    testStrange = 0.0;
     298
     299  } else { // target was a proton
     300    mult = GetMultiplicityT0(KE);
     301    fsTypes = GetFSPartTypesForNP(mult, KE);
     302 
     303    part1 = fsTypes[0];
     304    part2 = fsTypes[1];
     305    currentParticle.SetDefinition(particleDef[part1]);
     306    targetParticle.SetDefinition(particleDef[part2]);
     307    if (part1 == pro) {
     308      if (part2 == pro) {
     309        incidentHasChanged = true;
     310      } else if (part2 == neu) {
     311        if (G4UniformRand() > 0.5) {
     312          incidentHasChanged = true;
     313          targetHasChanged = true;
     314        } else {
     315          currentParticle.SetDefinition(particleDef[part2]);
     316          targetParticle.SetDefinition(particleDef[part1]);
     317        }
     318       
     319      } else if (part2 > neu && part2 < xi0) {
     320        incidentHasChanged = true;
     321        targetHasChanged = true;
    193322      }
    194       else
    195       {
    196         if( pz < 0.0 )pz *= -1.0;
    197       }
    198       G4double pu = std::sqrt(px*px+py*py+pz*pz);
    199       modifiedOriginal.SetMomentum( targetParticle.GetMomentum() * (p/pu) );
    200       modifiedOriginal.SetKineticEnergy( ek*GeV );
    201      
    202       targetParticle.SetMomentum(
    203        originalIncident->Get4Momentum().vect() - modifiedOriginal.GetMomentum() );
    204       G4double pp = targetParticle.GetMomentum().mag();
    205       G4double tarmas = targetParticle.GetMass();
    206       targetParticle.SetTotalEnergy( std::sqrt( pp*pp + tarmas*tarmas ) );
    207      
    208       theParticleChange.SetEnergyChange( modifiedOriginal.GetKineticEnergy() );
    209       G4DynamicParticle *pd = new G4DynamicParticle;
    210       pd->SetDefinition( targetParticle.GetDefinition() );
    211       pd->SetMomentum( targetParticle.GetMomentum() );
    212       theParticleChange.AddSecondary( pd );
    213       return;
    214     }
    215     G4FastVector<G4ReactionProduct,4> vec;  // vec will contain the secondary particles
    216     G4int vecLen = 0;
    217     vec.Initialize( 0 );
    218    
    219     G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
    220     G4double massVec[9];
    221     massVec[0] = targetNucleus.AtomicMass( A+1.0, Z     );
    222     massVec[1] = theAtomicMass;
    223     massVec[2] = 0.;
    224     if (Z > 1.0)
    225         massVec[2] = targetNucleus.AtomicMass( A    , Z-1.0 );
    226     massVec[3] = 0.;
    227     if (Z > 1.0 && A > 1.0)
    228         massVec[3] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
    229     massVec[4] = 0.;
    230     if (Z > 1.0 && A > 2.0 && A-2.0 > Z-1.0)
    231         massVec[4] = targetNucleus.AtomicMass( A-2.0, Z-1.0 );
    232     massVec[5] = 0.;
    233     if (Z > 2.0 && A > 3.0 && A-3.0 > Z-2.0)
    234         massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-2.0 );
    235     massVec[6] = 0.;
    236     if (A > 1.0 && A-1.0 > Z)
    237         massVec[6] = targetNucleus.AtomicMass( A-1.0, Z     );
    238     massVec[7] = massVec[3];
    239     massVec[8] = 0.;
    240     if (Z > 2.0 && A > 1.0)
    241         massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-2.0 );
    242    
    243     twoBody.NuclearReaction(vec, vecLen, originalIncident,
    244                             targetNucleus, theAtomicMass, massVec );
    245    
    246     theParticleChange.SetStatusChange( stopAndKill );
    247     theParticleChange.SetEnergyChange( 0.0 );
    248    
    249     G4DynamicParticle * pd;
    250     for( G4int i=0; i<vecLen; ++i )
    251     {
    252       pd = new G4DynamicParticle();
    253       pd->SetDefinition( vec[i]->GetDefinition() );
    254       pd->SetMomentum( vec[i]->GetMomentum() );
    255       theParticleChange.AddSecondary( pd );
    256       delete vec[i];
    257     }
    258   }
    259  
    260  void
    261   G4RPGNeutronInelastic::Cascade(
    262    G4FastVector<G4ReactionProduct,256> &vec,
    263    G4int& vecLen,
    264    const G4HadProjectile *originalIncident,
    265    G4ReactionProduct &currentParticle,
    266    G4ReactionProduct &targetParticle,
    267    G4bool &incidentHasChanged,
    268    G4bool &targetHasChanged,
    269    G4bool &quasiElastic )
    270   {
    271     // derived from original FORTRAN code CASN by H. Fesefeldt (13-Sep-1987)
    272     //
    273     // Neutron undergoes interaction with nucleon within a nucleus.  Check if it is
    274     // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
    275     // occurs and input particle is degraded in energy. No other particles are produced.
    276     // If reaction is possible, find the correct number of pions/protons/neutrons
    277     // produced using an interpolation to multiplicity data.  Replace some pions or
    278     // protons/neutrons by kaons or strange baryons according to the average
    279     // multiplicity per Inelastic reaction.
    280     //
    281     const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
    282     const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
    283     const G4double targetMass = targetParticle.GetMass()/MeV;
    284     G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
    285                                         targetMass*targetMass +
    286                                         2.0*targetMass*etOriginal );
    287     G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
    288     if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
    289     {
    290       quasiElastic = true;
    291       return;
    292     }
    293     static G4bool first = true;
    294     const G4int numMul = 1200;
    295     const G4int numSec = 60;
    296     static G4double protmul[numMul], protnorm[numSec]; // proton constants
    297     static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
    298     // np = number of pi+, nm = number of pi-, nz = number of pi0
    299     G4int counter, nt=0, np=0, nm=0, nz=0;
    300     const G4double c = 1.25;   
    301     const G4double b[] = { 0.35, 0.0 };
    302     if( first )      // compute normalization constants, this will only be Done once
    303     {
    304       first = false;
    305       G4int i;
    306       for( i=0; i<numMul; ++i )protmul[i] = 0.0;
    307       for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
    308       counter = -1;
    309       for( np=0; np<numSec/3; ++np )
    310       {
    311         for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
    312         {
    313           for( nz=0; nz<numSec/3; ++nz )
    314           {
    315             if( ++counter < numMul )
    316             {
    317               nt = np+nm+nz;
    318               if( nt > 0 )
    319               {
    320                 protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c) /
    321                   (Factorial(1-np+nm)*Factorial(1+np-nm) );
    322                 protnorm[nt-1] += protmul[counter];
    323               }
    324             }
    325           }
    326         }
    327       }
    328       for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
    329       for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
    330       counter = -1;
    331       for( np=0; np<(numSec/3); ++np )
    332       {
    333         for( nm=np; nm<=(np+2); ++nm )
    334         {
    335           for( nz=0; nz<numSec/3; ++nz )
    336           {
    337             if( ++counter < numMul )
    338             {
    339               nt = np+nm+nz;
    340               if( (nt>0) && (nt<=numSec) )
    341               {
    342                 neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c) /
    343                   (Factorial(nm-np)*Factorial(2-nm+np) );
    344                 neutnorm[nt-1] += neutmul[counter];
    345               }
    346             }
    347           }
    348         }
    349       }
    350       for( i=0; i<numSec; ++i )
    351       {
    352         if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
    353         if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
    354       }
    355     }   // end of initialization
    356    
    357     const G4double expxu = 82.;      // upper bound for arg. of exp
    358     const G4double expxl = -expxu;        // lower bound for arg. of exp
    359     G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    360     G4ParticleDefinition *aProton = G4Proton::Proton();
    361     G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
    362     const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
    363     G4double test, w0, wp, wt, wm;
    364     if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
    365     {
    366       // suppress high multiplicity events at low momentum
    367       // only one pion will be produced
    368 
    369       nm = np = nz = 0;
    370       if( targetParticle.GetDefinition() == aNeutron )
    371       {
    372         test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
    373         w0 = test/2.0;
    374         wm = test;
    375         if( G4UniformRand() < w0/(w0+wm) )
    376           nz = 1;
    377         else
    378           nm = 1;
    379       } else { // target is a proton
    380         test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
    381         w0 = test;
    382         wp = test/2.0;       
    383         test = std::exp( std::min( expxu, std::max( expxl, -(-1.0+b[0])*(-1.0+b[0])/(2.0*c*c) ) ) );
    384         wm = test/2.0;
    385         wt = w0+wp+wm;
    386         wp += w0;
    387         G4double ran = G4UniformRand();
    388         if( ran < w0/wt )
    389           nz = 1;
    390         else if( ran < wp/wt )
    391           np = 1;
    392         else
    393           nm = 1;
    394       }
    395     } else {  // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
    396       G4double n, anpn;
    397       GetNormalizationConstant( availableEnergy, n, anpn );
    398       G4double ran = G4UniformRand();
    399       G4double dum, excs = 0.0;
    400       if( targetParticle.GetDefinition() == aProton )
    401       {
    402         counter = -1;
    403         for( np=0; np<numSec/3 && ran>=excs; ++np )
    404         {
    405           for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
    406           {
    407             for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
    408             {
    409               if( ++counter < numMul )
    410               {
    411                 nt = np+nm+nz;
    412                 if( nt > 0 )
    413                 {
    414                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
    415                   dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
    416                   if( std::fabs(dum) < 1.0 ) {
    417                     if( test >= 1.0e-10 )excs += dum*test;
    418                   } else {
    419                     excs += dum*test;
    420                   }
    421                 }
    422               }
    423             }
    424           }
    425         }
    426         if( ran >= excs )  // 3 previous loops continued to the end
    427         {
    428           quasiElastic = true;
    429           return;
    430         }
    431         np--; nm--; nz--;
    432       } else { // target must be a neutron
    433         counter = -1;
    434         for( np=0; np<numSec/3 && ran>=excs; ++np )
    435         {
    436           for( nm=np; nm<=(np+2) && ran>=excs; ++nm )
    437           {
    438             for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
    439             {
    440               if( ++counter < numMul )
    441               {
    442                 nt = np+nm+nz;
    443                 if( (nt>=1) && (nt<=numSec) )
    444                 {
    445                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
    446                   dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
    447                   if( std::fabs(dum) < 1.0 ) {
    448                     if( test >= 1.0e-10 )excs += dum*test;
    449                   } else {
    450                     excs += dum*test;
    451                   }
    452                 }
    453               }
    454             }
    455           }
    456         }
    457         if( ran >= excs )  // 3 previous loops continued to the end
    458         {
    459           quasiElastic = true;
    460           return;
    461         }
    462         np--; nm--; nz--;
    463       }
    464     }
    465     if( targetParticle.GetDefinition() == aProton )
    466     {
    467       switch( np-nm )
    468       {
    469        case 0:
    470          if( G4UniformRand() < 0.33 )
    471          {
    472            currentParticle.SetDefinitionAndUpdateE( aProton );
    473            targetParticle.SetDefinitionAndUpdateE( aNeutron );
    474            incidentHasChanged = true;
    475            targetHasChanged = true;
    476          }
    477          break;
    478        case 1:
    479          targetParticle.SetDefinitionAndUpdateE( aNeutron );
    480          targetHasChanged = true;
    481          break;
    482        default:
    483          currentParticle.SetDefinitionAndUpdateE( aProton );
    484          incidentHasChanged = true;
    485          break;
    486       }
    487     } else { // target must be a neutron
    488       switch( np-nm )
    489       {
    490        case -1:                       // changed from +1 by JLC, 7Jul97
    491          if( G4UniformRand() < 0.5 )
    492          {
    493            currentParticle.SetDefinitionAndUpdateE( aProton );
    494            incidentHasChanged = true;
    495          } else {
    496            targetParticle.SetDefinitionAndUpdateE( aProton );
    497            targetHasChanged = true;
    498          }
    499          break;
    500        case 0:
    501          break;
    502        default:
    503          currentParticle.SetDefinitionAndUpdateE( aProton );
    504          targetParticle.SetDefinitionAndUpdateE( aProton );
    505          incidentHasChanged = true;
    506          targetHasChanged = true;
    507          break;
    508       }
    509     }
    510     SetUpPions( np, nm, nz, vec, vecLen );
    511 // DEBUG -->    DumpFrames::DumpFrame(vec, vecLen);
    512     return;
    513   }
    514 
    515  /* end of file */
    516  
     323
     324    } else { // neutron
     325      targetHasChanged = true;
     326    }
     327
     328    testCharge = 1.0;
     329    testBaryon = 2.0;
     330    testStrange = 0.0;
     331  }
     332
     333  //  if (mult == 2 && !incidentHasChanged && !targetHasChanged)
     334  //                                              quasiElastic = true;
     335 
     336  // Remove incident and target from fsTypes
     337 
     338  fsTypes.erase(fsTypes.begin());
     339  fsTypes.erase(fsTypes.begin());
     340
     341  // Remaining particles are secondaries.  Put them into vec.
     342 
     343  G4ReactionProduct* rp(0);
     344  for(G4int i=0; i < mult-2; ++i ) {
     345    partType = fsTypes[i];
     346    rp = new G4ReactionProduct();
     347    rp->SetDefinition(particleDef[partType]);
     348    (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
     349    vec.SetElement(vecLen++, rp);
     350  }
     351
     352  // Check conservation of charge, strangeness, baryon number
     353 
     354  CheckQnums(vec, vecLen, currentParticle, targetParticle,
     355             testCharge, testBaryon, testStrange);
     356 
     357  return;
     358}
Note: See TracChangeset for help on using the changeset viewer.