Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc (Version 9.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 // $Id: G4ChargeExchange.cc,v 1.11 2007/05/25 17:46:52 dennis Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01 $
                                                   >>  29 //
                                                   >>  30 //
 26 // G4 Model: Charge and strangness exchange ba     31 // G4 Model: Charge and strangness exchange based on G4LightMedia model
 27 //           28 May 2006 V.Ivanchenko              32 //           28 May 2006 V.Ivanchenko
 28 //                                                 33 //
 29 // Modified:                                       34 // Modified:
 30 // 07-Jun-06 V.Ivanchenko fix problem of rotat     35 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
 31 // 25-Jul-06 V.Ivanchenko add 19 MeV low energ     36 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
 32 // 12-Jun-12 A.Ribon fix warnings of shadowed  << 
 33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log << 
 34 //                                                 37 //
 35                                                    38 
 36 #include "G4ChargeExchange.hh"                     39 #include "G4ChargeExchange.hh"
 37 #include "G4ChargeExchangeXS.hh"               << 
 38 #include "G4PhysicalConstants.hh"              << 
 39 #include "G4SystemOfUnits.hh"                  << 
 40 #include "G4ParticleTable.hh"                      40 #include "G4ParticleTable.hh"
 41 #include "G4ParticleDefinition.hh"                 41 #include "G4ParticleDefinition.hh"
 42 #include "G4IonTable.hh"                           42 #include "G4IonTable.hh"
                                                   >>  43 #include "G4QElasticCrossSection.hh"
                                                   >>  44 #include "G4VQCrossSection.hh"
                                                   >>  45 #include "G4ElasticHadrNucleusHE.hh"
 43 #include "Randomize.hh"                            46 #include "Randomize.hh"
 44 #include "G4NucleiProperties.hh"               <<  47 #include "G4HadronElastic.hh"
 45 #include "G4DecayTable.hh"                     << 
 46 #include "G4VDecayChannel.hh"                  << 
 47 #include "G4DecayProducts.hh"                  << 
 48 #include "G4NistManager.hh"                    << 
 49 #include "G4Fragment.hh"                       << 
 50 #include "G4ExcitationHandler.hh"              << 
 51 #include "G4ReactionProductVector.hh"          << 
 52                                                << 
 53 #include "G4Exp.hh"                            << 
 54 #include "G4Log.hh"                            << 
 55 #include "G4Pow.hh"                            << 
 56                                                    48 
 57 #include "G4HadronicParameters.hh"             << 
 58 #include "G4PhysicsModelCatalog.hh"            << 
 59                                                << 
 60 namespace                                      << 
 61 {                                              << 
 62   constexpr G4int maxN = 1000;                 << 
 63   constexpr G4double emin = 2*136.9*CLHEP::MeV << 
 64 }                                              << 
 65                                                    49 
 66 G4ChargeExchange::G4ChargeExchange(G4ChargeExc <<  50 G4ChargeExchange::G4ChargeExchange(G4HadronElastic* hel, G4double elim,
 67   : G4HadronicInteraction("ChargeExchange"),   <<  51                                    G4double ehigh)
 68     fXSection(ptr), fXSWeightFactor(1.0)       <<  52 : G4HadronicInteraction("G4ChargeExchange"),
                                                   >>  53   fElastic(hel),
                                                   >>  54   native(false),
                                                   >>  55   ekinlim(elim),
                                                   >>  56   ekinhigh(ehigh)
 69 {                                                  57 {
 70   lowEnergyLimit = 1.*CLHEP::MeV;              <<  58   SetMinEnergy( 0.0*GeV );
 71   secID = G4PhysicsModelCatalog::GetModelID( " <<  59   SetMaxEnergy( 100.*TeV );
 72   nist = G4NistManager::Instance();            <<  60   ekinlow = 19.0*MeV;
 73   fHandler = new G4ExcitationHandler();        <<  61 
 74   if (nullptr != fXSection) {                  <<  62   verboseLevel= 0;
 75     fXSWeightFactor = 1.0/fXSection->GetCrossS <<  63   if(!fElastic) {
 76   }                                            <<  64     native = true;
                                                   >>  65     fElastic = new G4HadronElastic();
                                                   >>  66   }
                                                   >>  67   qCManager   = fElastic->GetCS();
                                                   >>  68   hElastic    = fElastic->GetHElastic();
                                                   >>  69 
                                                   >>  70   theProton   = G4Proton::Proton();
                                                   >>  71   theNeutron  = G4Neutron::Neutron();
                                                   >>  72   theAProton  = G4AntiProton::AntiProton();
                                                   >>  73   theANeutron = G4AntiNeutron::AntiNeutron();
                                                   >>  74   thePiPlus   = G4PionPlus::PionPlus();
                                                   >>  75   thePiMinus  = G4PionMinus::PionMinus();
                                                   >>  76   thePiZero   = G4PionZero::PionZero();
                                                   >>  77   theKPlus    = G4KaonPlus::KaonPlus();
                                                   >>  78   theKMinus   = G4KaonMinus::KaonMinus();
                                                   >>  79   theK0S      = G4KaonZeroShort::KaonZeroShort();
                                                   >>  80   theK0L      = G4KaonZeroLong::KaonZeroLong();
                                                   >>  81   theL        = G4Lambda::Lambda();
                                                   >>  82   theAntiL    = G4AntiLambda::AntiLambda();
                                                   >>  83   theSPlus    = G4SigmaPlus::SigmaPlus();
                                                   >>  84   theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
                                                   >>  85   theSMinus   = G4SigmaMinus::SigmaMinus();
                                                   >>  86   theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
                                                   >>  87   theS0       = G4SigmaZero::SigmaZero();
                                                   >>  88   theAS0      = G4AntiSigmaZero::AntiSigmaZero();
                                                   >>  89   theXiMinus  = G4XiMinus::XiMinus();
                                                   >>  90   theXi0      = G4XiZero::XiZero();
                                                   >>  91   theAXiMinus = G4AntiXiMinus::AntiXiMinus();
                                                   >>  92   theAXi0     = G4AntiXiZero::AntiXiZero();
                                                   >>  93   theOmega    = G4OmegaMinus::OmegaMinus();
                                                   >>  94   theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
                                                   >>  95   theD        = G4Deuteron::Deuteron();
                                                   >>  96   theT        = G4Triton::Triton();
                                                   >>  97   theA        = G4Alpha::Alpha();
                                                   >>  98   theA        = G4He3::He3();
 77 }                                                  99 }
 78                                                   100 
 79 G4ChargeExchange::~G4ChargeExchange()             101 G4ChargeExchange::~G4ChargeExchange()
 80 {                                                 102 {
 81   delete fHandler;                             << 103   if(native) delete fElastic;
 82 }                                                 104 }
 83                                                   105 
 84 G4HadFinalState* G4ChargeExchange::ApplyYourse    106 G4HadFinalState* G4ChargeExchange::ApplyYourself(
 85      const G4HadProjectile& aTrack, G4Nucleus&    107      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 86 {                                                 108 {
 87   theParticleChange.Clear();                      109   theParticleChange.Clear();
 88   auto part = aTrack.GetDefinition();          << 110   const G4HadProjectile* aParticle = &aTrack;
 89   G4double ekin = aTrack.GetKineticEnergy();   << 111   G4double ekin = aParticle->GetKineticEnergy();
                                                   >> 112 
                                                   >> 113   G4double aTarget = targetNucleus.GetN();
                                                   >> 114   G4double zTarget = targetNucleus.GetZ();
 90                                                   115 
 91   G4int A = targetNucleus.GetA_asInt();        << 116   G4int Z = static_cast<G4int>(zTarget+0.5);
 92   G4int Z = targetNucleus.GetZ_asInt();        << 117   G4int A = static_cast<G4int>(aTarget+0.5);
 93                                                   118 
 94   if (ekin <= lowEnergyLimit) {                << 119   if(ekin == 0.0 || A < 3) {
                                                   >> 120     theParticleChange.SetEnergyChange(ekin);
                                                   >> 121     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
 95     return &theParticleChange;                    122     return &theParticleChange;
 96   }                                               123   }
 97   theParticleChange.SetWeightChange(fXSWeightF << 
 98                                                   124 
 99   G4int projPDG = part->GetPDGEncoding();      << 125   G4double plab = aParticle->GetTotalMomentum();
                                                   >> 126 
                                                   >> 127   if (verboseLevel > 1)
                                                   >> 128     G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
                                                   >> 129      << plab/GeV << " GeV/c "
                                                   >> 130      << " ekin(MeV) = " << ekin/MeV << "  "
                                                   >> 131      << aParticle->GetDefinition()->GetParticleName() << G4endl;
                                                   >> 132 
                                                   >> 133   // Scattered particle referred to axis of incident particle
                                                   >> 134   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
                                                   >> 135   G4double m1 = theParticle->GetPDGMass();
100                                                   136 
101   // for hydrogen targets and positive project << 137   G4int N = A - Z;
102   // is not possible on proton, only on deuter << 138   G4int projPDG = theParticle->GetPDGEncoding();
103   if (1 == Z && (211 == projPDG || 321 == proj << 
104                                                << 
105   if (verboseLevel > 1)                           139   if (verboseLevel > 1)
106     G4cout << "G4ChargeExchange for " << part- << 140     G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
107      << " PDGcode= " << projPDG << " on nucleu    141      << " PDGcode= " << projPDG << " on nucleus Z= " << Z
108      << " A= " << A << " N= " << A - Z         << 142      << " A= " << A << " N= " << N
109      << G4endl;                                   143      << G4endl;
110                                                   144 
111   G4double mass1 = G4NucleiProperties::GetNucl << 145   G4ParticleDefinition * theDef = 0;
112   G4LorentzVector lv0 = aTrack.Get4Momentum(); << 146 
113   G4double etot = mass1 + lv0.e();             << 147   if (Z == 1 && A == 3) theDef = theT;
114                                                << 148   else if (Z == 2 && A == 3) theDef = theHe3;
115   // select final state                        << 149   else if (Z == 2 && A == 4) theDef = theA;
116   const G4ParticleDefinition* theSecondary =   << 150   else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
117     fXSection->SampleSecondaryType(part, Z, A) << 151 
118   G4int pdg = theSecondary->GetPDGEncoding();  << 152   G4double m2 = theDef->GetPDGMass();
119                                                << 153   G4LorentzVector lv1 = aParticle->Get4Momentum();
120   // omega(782) and f2(1270)                   << 154   G4LorentzVector lv0(0.0,0.0,0.0,m2);
121   G4bool isShortLived = (pdg == 223 || pdg ==  << 155 
122                                                << 156   G4LorentzVector lv  = lv0 + lv1;
123   // atomic number of the recoil nucleus       << 157   G4ThreeVector bst = lv.boostVector();
124   if (projPDG == -211) { --Z; }                << 158   lv1.boost(-bst);
125   else if (projPDG == 211) { ++Z; }            << 159   lv0.boost(-bst);
126   else if (projPDG == -321) { --Z; }           << 
127   else if (projPDG == 321) { ++Z; }            << 
128   else if (projPDG == 130) {                   << 
129     if (theSecondary->GetPDGCharge() > 0.0) {  << 
130     else { ++Z; }                              << 
131   } else {                                     << 
132     // not ready for other projectile          << 
133     return &theParticleChange;                 << 
134   }                                            << 
135                                                   160 
136   // recoil nucleus                            << 161   // Sample final particles
137   const G4ParticleDefinition* theRecoil = null << 162   G4bool theHyperon = false;
138   if (Z == 0 && A == 1) { theRecoil = G4Neutro << 163   G4ParticleDefinition* theRecoil = 0;
139   else if (Z == 1 && A == 1) { theRecoil = G4P << 164   G4ParticleDefinition* theSecondary = 0;
140   else if (Z == 1 && A == 2) { theRecoil = G4D << 165 
141   else if (Z == 1 && A == 3) { theRecoil = G4T << 166   if(theParticle == theProton) {
142   else if (Z == 2 && A == 3) { theRecoil = G4H << 167     theSecondary = theNeutron;
143   else if (Z == 2 && A == 4) { theRecoil = G4A << 168     Z++;
144   else if (nist->GetIsotopeAbundance(Z, A) > 0 << 169   } else if(theParticle == theNeutron) {
145     theRecoil = G4ParticleTable::GetParticleTa << 170     theSecondary = theProton;
146       ->GetIonTable()->GetIon(Z, A, 0.0);      << 171     Z--;
                                                   >> 172   } else if(theParticle == thePiPlus) {
                                                   >> 173     theSecondary = thePiZero;
                                                   >> 174     Z++;
                                                   >> 175   } else if(theParticle == thePiMinus) {
                                                   >> 176     theSecondary = thePiZero;
                                                   >> 177     Z--;
                                                   >> 178   } else if(theParticle == theKPlus) {
                                                   >> 179     if(G4UniformRand()<0.5) theSecondary = theK0S;
                                                   >> 180     else  theSecondary = theK0L;
                                                   >> 181     Z++;
                                                   >> 182   } else if(theParticle == theKMinus) {
                                                   >> 183     if(G4UniformRand()<0.5) theSecondary = theK0S;
                                                   >> 184     else  theSecondary = theK0L;
                                                   >> 185     Z--;
                                                   >> 186   } else if(theParticle == theK0S || theParticle == theK0L) {
                                                   >> 187     if(G4UniformRand()*aTarget < zTarget) {
                                                   >> 188       theSecondary = theKPlus;
                                                   >> 189       Z--;
                                                   >> 190     } else {
                                                   >> 191       theSecondary = theKMinus;
                                                   >> 192       Z++;
                                                   >> 193     }
                                                   >> 194   } else if(theParticle == theANeutron) {
                                                   >> 195     theSecondary = theAProton;
                                                   >> 196     Z++;
                                                   >> 197   } else if(theParticle == theAProton) {
                                                   >> 198     theSecondary = theANeutron;
                                                   >> 199     Z--;
                                                   >> 200   } else if(theParticle == theL) {
                                                   >> 201     G4double x = G4UniformRand();
                                                   >> 202     if(G4UniformRand()*aTarget < zTarget) {
                                                   >> 203       if(x < 0.2) {
                                                   >> 204         theSecondary = theS0;
                                                   >> 205       } else if (x < 0.4) {
                                                   >> 206         theSecondary = theSPlus;
                                                   >> 207         Z--;
                                                   >> 208       } else if (x < 0.6) {
                                                   >> 209         theSecondary = theProton;
                                                   >> 210   theRecoil = theL;
                                                   >> 211         theHyperon = true;
                                                   >> 212   A--;
                                                   >> 213       } else if (x < 0.8) {
                                                   >> 214         theSecondary = theProton;
                                                   >> 215   theRecoil = theS0;
                                                   >> 216         theHyperon = true;
                                                   >> 217         A--;
                                                   >> 218       } else {
                                                   >> 219         theSecondary = theNeutron;
                                                   >> 220   theRecoil = theSPlus;
                                                   >> 221         theHyperon = true;
                                                   >> 222         A--;
                                                   >> 223       }
                                                   >> 224     } else {
                                                   >> 225       if(x < 0.2) {
                                                   >> 226         theSecondary = theS0;
                                                   >> 227       } else if (x < 0.4) {
                                                   >> 228         theSecondary = theSMinus;
                                                   >> 229         Z++;
                                                   >> 230       } else if (x < 0.6) {
                                                   >> 231         theSecondary = theNeutron;
                                                   >> 232   theRecoil = theL;
                                                   >> 233         A--;
                                                   >> 234         theHyperon = true;
                                                   >> 235       } else if (x < 0.8) {
                                                   >> 236         theSecondary = theNeutron;
                                                   >> 237   theRecoil = theS0;
                                                   >> 238         theHyperon = true;
                                                   >> 239         A--;
                                                   >> 240       } else {
                                                   >> 241         theSecondary = theProton;
                                                   >> 242   theRecoil = theSMinus;
                                                   >> 243         theHyperon = true;
                                                   >> 244         A--;
                                                   >> 245       }
                                                   >> 246     }
147   }                                               247   }
148                                                   248 
149   // check if there is enough energy for the f << 249   if (Z == 1 && A == 2) theDef = theD;
150   // and sample mass of produced state         << 250   else if (Z == 1 && A == 3) theDef = theT;
151   const G4double mass0 = theSecondary->GetPDGM << 251   else if (Z == 2 && A == 3) theDef = theHe3;
152   G4double mass3 = (nullptr == theRecoil) ?    << 252   else if (Z == 2 && A == 4) theDef = theA;
153     G4NucleiProperties::GetNuclearMass(A, Z) : << 253   else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
154   G4double mass2 = mass0;                      << 254 
155   if (isShortLived &&                          << 255   G4double m11 = theSecondary->GetPDGMass();
156       !SampleMass(mass2, theSecondary->GetPDGW << 256   G4double m21 = theDef->GetPDGMass();
157     return &theParticleChange;                 << 257   if(theRecoil)  m21 += theRecoil->GetPDGMass();
                                                   >> 258   else           theRecoil = theDef;
                                                   >> 259 
                                                   >> 260   G4double etot = lv0.e() + lv1.e();
                                                   >> 261   if(etot < m11 + m21) return &theParticleChange;
                                                   >> 262 
                                                   >> 263   G4ThreeVector p1 = lv1.vect();
                                                   >> 264   G4double e1 = 0.5*etot*(1.0 + (m21*m21 - m11*m11)/(etot*etot));
                                                   >> 265   //  G4double e2 = etot - e1;
                                                   >> 266   G4double ptot = std::sqrt(e1*e1 - m11*m11);
                                                   >> 267 
                                                   >> 268   G4double tmax = 4.0*ptot*ptot;
                                                   >> 269   G4double t = 0.0;
                                                   >> 270 
                                                   >> 271   // Choose generator
                                                   >> 272   G4ElasticGenerator gtype = fLElastic;
                                                   >> 273   if ((theParticle == theProton || theParticle == theNeutron) && 
                                                   >> 274       Z <= 2 && ekin >= ekinlow) {
                                                   >> 275     gtype = fQElastic;
                                                   >> 276   } else {
                                                   >> 277     if(ekin >= ekinlow)       gtype = fSWave;
                                                   >> 278     else if(ekin >= ekinhigh) gtype = fHElastic;
158   }                                               279   }
159                                                   280 
160   // not possible kinematically                << 281   // Sample t
161   if (etot <= mass2 + mass3) {                 << 282   if(gtype == fQElastic) {
162     return &theParticleChange;                 << 283     if (verboseLevel > 1)
                                                   >> 284       G4cout << "G4ChargeExchange: Z= " << Z << " N= "
                                                   >> 285        << N << " pdg= " <<  projPDG
                                                   >> 286        << " mom(GeV)= " << plab/GeV << "  " << qCManager << G4endl;
                                                   >> 287     if(Z == 1 && N == 2) N = 1;
                                                   >> 288     else if (Z == 2 && N == 1) N = 2;
                                                   >> 289     G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
                                                   >> 290     if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
                                                   >> 291     else gtype = fSWave;
                                                   >> 292   }
                                                   >> 293 
                                                   >> 294   if(gtype == fHElastic) {
                                                   >> 295     t = hElastic->SampleT(theParticle,plab,Z,A);
                                                   >> 296     if(t > tmax) gtype = fSWave;
                                                   >> 297   }
                                                   >> 298 
                                                   >> 299   if(gtype == fLElastic) {
                                                   >> 300     t = GeV*GeV*fElastic->SampleT(ptot,m1,m2,aTarget);
                                                   >> 301     if(t > tmax) gtype = fSWave;
                                                   >> 302   }
                                                   >> 303 
                                                   >> 304   // NaN finder
                                                   >> 305   if(!(t < 0.0 || t >= 0.0)) {
                                                   >> 306     if (verboseLevel > -1) {
                                                   >> 307       G4cout << "G4ChargeExchange:WARNING: Z= " << Z << " N= " 
                                                   >> 308        << N << " pdg= " <<  projPDG
                                                   >> 309        << " mom(GeV)= " << plab/GeV 
                                                   >> 310        << " the model type " << gtype;
                                                   >> 311       if(gtype ==  fQElastic) G4cout << " CHIPS ";
                                                   >> 312       else if(gtype ==  fLElastic) G4cout << " LElastic ";
                                                   >> 313       else if(gtype ==  fHElastic) G4cout << " HElastic ";
                                                   >> 314       G4cout << " t= " << t
                                                   >> 315        << " S-wave will be sampled" 
                                                   >> 316        << G4endl; 
                                                   >> 317     }
                                                   >> 318     gtype = fSWave;
163   }                                               319   }
164                                                   320 
165   // sample kinematics                         << 321   if(gtype == fSWave) t = G4UniformRand()*tmax;
166   G4LorentzVector lv1(0.0, 0.0, 0.0, mass1);   << 
167   G4LorentzVector lv = lv0 + lv1;              << 
168   G4ThreeVector bst = lv.boostVector();        << 
169   G4double ss = lv.mag2();                     << 
170                                                   322 
171   // tmax = 4*momCMS^2                         << 323   if(verboseLevel>1)
172   G4double e2 = ss + mass2*mass2 - mass3*mass3 << 324     G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax
173   G4double tmax = e2*e2/ss - 4*mass2*mass2;    << 325      << " ptot= " << ptot << G4endl;
174                                                << 326 
175   G4double t = SampleT(theSecondary, A, tmax); << 327   // Sampling in CM system
176                                                << 328   G4double phi  = G4UniformRand()*twopi;
177   G4double phi  = G4UniformRand()*CLHEP::twopi << 
178   G4double cost = 1. - 2.0*t/tmax;                329   G4double cost = 1. - 2.0*t/tmax;
                                                   >> 330   if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
                                                   >> 331   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
179                                                   332 
180   if (cost > 1.0) { cost = 1.0; }              << 333   if (verboseLevel > 1)
181   else if(cost < -1.0) { cost = -1.0; }        << 334     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
182                                                   335 
183   G4double sint = std::sqrt((1.0-cost)*(1.0+co << 336   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
                                                   >> 337   v1 *= ptot;
                                                   >> 338   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
                                                   >> 339   G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
184                                                   340 
185   if (verboseLevel>1) {                        << 341   nlv0.boost(bst);
186     G4cout << " t= " << t << " tmax(GeV^2)= "  << 342   nlv1.boost(bst);
187      << " cos(t)=" << cost << " sin(t)=" << si << 
188   }                                            << 
189   G4double momentumCMS = 0.5*std::sqrt(tmax);  << 
190   G4LorentzVector lv2(momentumCMS*sint*std::co << 
191           momentumCMS*sint*std::sin(phi),      << 
192           momentumCMS*cost,                    << 
193           std::sqrt(momentumCMS*momentumCMS +  << 
194                                                << 
195   // kinematics in the final state, may be a w << 
196   lv2.boost(bst);                              << 
197   if (lv2.e() < mass2) {                       << 
198     lv2.setE(mass2);                           << 
199   }                                            << 
200   lv -= lv2;                                   << 
201   if (lv.e() < mass3) {                        << 
202     lv.setE(mass3);                            << 
203   }                                            << 
204                                                   343 
205   // prepare secondary particles               << 
206   theParticleChange.SetStatusChange(stopAndKil    344   theParticleChange.SetStatusChange(stopAndKill);
207   theParticleChange.SetEnergyChange(0.0);      << 345   G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
208                                                << 346   theParticleChange.AddSecondary(aSec);
209   if (!isShortLived) {                         << 
210     auto aSec = new G4DynamicParticle(theSecon << 
211     theParticleChange.AddSecondary(aSec, secID << 
212   } else {                                     << 
213     auto channel = theSecondary->GetDecayTable << 
214     auto products = channel->DecayIt(mass2);   << 
215     G4ThreeVector bst1 = lv2.boostVector();    << 
216     G4int N = products->entries();             << 
217     for (G4int i=0; i<N; ++i) {                << 
218       auto p = (*products)[i];                 << 
219       auto lvp = p->Get4Momentum();            << 
220       lvp.boost(bst1);                         << 
221       p->Set4Momentum(lvp);                    << 
222       theParticleChange.AddSecondary(p, secID) << 
223     }                                          << 
224     delete products;                           << 
225   }                                            << 
226                                                   347 
227   // recoil is a stable isotope                << 348   G4double erec =  nlv0.e() - m21;
228   if (nullptr != theRecoil) {                  << 349   if(theHyperon) {
229     auto aRec = new G4DynamicParticle(theRecoi << 350     theParticleChange.SetLocalEnergyDeposit(erec);
230     theParticleChange.AddSecondary(aRec, secID << 351     aSec = new G4DynamicParticle();
                                                   >> 352     aSec->SetDefinition(theRecoil);
                                                   >> 353     aSec->SetKineticEnergy(0.0);
                                                   >> 354   } else if(erec > ekinlim) {
                                                   >> 355     aSec = new G4DynamicParticle(theRecoil, nlv0);
                                                   >> 356     theParticleChange.AddSecondary(aSec);
231   } else {                                        357   } else {
232     // recoil is an unstable fragment          << 358     theParticleChange.SetLocalEnergyDeposit(erec);
233     G4Fragment frag(A, Z, lv);                 << 
234     auto products = fHandler->BreakItUp(frag); << 
235     for (auto & prod : *products) {            << 
236       auto dp = new G4DynamicParticle(prod->Ge << 
237       theParticleChange.AddSecondary(dp, secID << 
238       delete prod;                             << 
239     }                                          << 
240     delete products;                           << 
241   }                                               359   }
242   return &theParticleChange;                      360   return &theParticleChange;
243 }                                              << 
244                                                << 
245 G4double G4ChargeExchange::SampleT(const G4Par << 
246                                    const G4int << 
247 {                                              << 
248   G4double aa, bb, cc, dd;                     << 
249   G4Pow* g4pow = G4Pow::GetInstance();         << 
250   if (A <= 62.) {                              << 
251     aa = g4pow->powZ(A, 1.63);                 << 
252     bb = 14.5*g4pow->powZ(A, 0.66);            << 
253     cc = 1.4*g4pow->powZ(A, 0.33);             << 
254     dd = 10.;                                  << 
255   } else {                                     << 
256     aa = g4pow->powZ(A, 1.33);                 << 
257     bb = 60.*g4pow->powZ(A, 0.33);             << 
258     cc = 0.4*g4pow->powZ(A, 0.40);             << 
259     dd = 10.;                                  << 
260   }                                            << 
261   G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb; << 
262   G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd; << 
263                                                << 
264   G4double t;                                  << 
265   G4double y = bb;                             << 
266   if(G4UniformRand()*(x1 + x2) < x2) y = dd;   << 
267                                                << 
268   for (G4int i=0; i<maxN; ++i) {               << 
269     t = -G4Log(G4UniformRand())/y;             << 
270     if (t <= tmax) { return t; }               << 
271   }                                            << 
272   return 0.0;                                  << 
273 }                                              << 
274                                                << 
275 G4bool G4ChargeExchange::SampleMass(G4double&  << 
276 {                                              << 
277   // +- 4 width but above 2 pion mass          << 
278   const G4double e1 = std::max(M - 4*G, emin); << 
279   const G4double e2 = std::min(M + 4*G, elim)  << 
280   if (e2 <= 0.0) { return false; }             << 
281   const G4double M2 = M*M;                     << 
282   const G4double MG2 = M2*G*G;                 << 
283                                                << 
284   // sampling Breit-Wigner function            << 
285   for (G4int i=0; i<maxN; ++i) {               << 
286     G4double e = e1 + e2*G4UniformRand();      << 
287     G4double x = e*e - M2;                     << 
288     G4double y = MG2/(x*x + MG2);              << 
289     if (y >= G4UniformRand()) {                << 
290       M = e;                                   << 
291       return true;                             << 
292     }                                          << 
293   }                                            << 
294   return false;                                << 
295 }                                                 361 }
296                                                   362