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 10.7)


  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 //
                                                   >>  28 //
 26 // G4 Model: Charge and strangness exchange ba     29 // G4 Model: Charge and strangness exchange based on G4LightMedia model
 27 //           28 May 2006 V.Ivanchenko              30 //           28 May 2006 V.Ivanchenko
 28 //                                                 31 //
 29 // Modified:                                       32 // Modified:
 30 // 07-Jun-06 V.Ivanchenko fix problem of rotat     33 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
 31 // 25-Jul-06 V.Ivanchenko add 19 MeV low energ     34 // 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      35 // 12-Jun-12 A.Ribon fix warnings of shadowed variables
 33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log     36 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow
 34 //                                                 37 //
 35                                                    38 
 36 #include "G4ChargeExchange.hh"                     39 #include "G4ChargeExchange.hh"
 37 #include "G4ChargeExchangeXS.hh"               << 
 38 #include "G4PhysicalConstants.hh"                  40 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      41 #include "G4SystemOfUnits.hh"
 40 #include "G4ParticleTable.hh"                      42 #include "G4ParticleTable.hh"
 41 #include "G4ParticleDefinition.hh"                 43 #include "G4ParticleDefinition.hh"
 42 #include "G4IonTable.hh"                           44 #include "G4IonTable.hh"
 43 #include "Randomize.hh"                            45 #include "Randomize.hh"
 44 #include "G4NucleiProperties.hh"                   46 #include "G4NucleiProperties.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                                                    47 
 53 #include "G4Exp.hh"                                48 #include "G4Exp.hh"
 54 #include "G4Log.hh"                                49 #include "G4Log.hh"
 55 #include "G4Pow.hh"                                50 #include "G4Pow.hh"
 56                                                    51 
 57 #include "G4HadronicParameters.hh"                 52 #include "G4HadronicParameters.hh"
 58 #include "G4PhysicsModelCatalog.hh"            << 
 59                                                    53 
 60 namespace                                      << 
 61 {                                              << 
 62   constexpr G4int maxN = 1000;                 << 
 63   constexpr G4double emin = 2*136.9*CLHEP::MeV << 
 64 }                                              << 
 65                                                    54 
 66 G4ChargeExchange::G4ChargeExchange(G4ChargeExc <<  55 G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange")
 67   : G4HadronicInteraction("ChargeExchange"),   << 
 68     fXSection(ptr), fXSWeightFactor(1.0)       << 
 69 {                                                  56 {
 70   lowEnergyLimit = 1.*CLHEP::MeV;              <<  57   SetMinEnergy( 0.0*GeV );
 71   secID = G4PhysicsModelCatalog::GetModelID( " <<  58   SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
 72   nist = G4NistManager::Instance();            <<  59 
 73   fHandler = new G4ExcitationHandler();        <<  60   lowestEnergyLimit    = 1.*MeV;  
 74   if (nullptr != fXSection) {                  <<  61 
 75     fXSWeightFactor = 1.0/fXSection->GetCrossS <<  62   theProton   = G4Proton::Proton();
 76   }                                            <<  63   theNeutron  = G4Neutron::Neutron();
                                                   >>  64   theAProton  = G4AntiProton::AntiProton();
                                                   >>  65   theANeutron = G4AntiNeutron::AntiNeutron();
                                                   >>  66   thePiPlus   = G4PionPlus::PionPlus();
                                                   >>  67   thePiMinus  = G4PionMinus::PionMinus();
                                                   >>  68   thePiZero   = G4PionZero::PionZero();
                                                   >>  69   theKPlus    = G4KaonPlus::KaonPlus();
                                                   >>  70   theKMinus   = G4KaonMinus::KaonMinus();
                                                   >>  71   theK0S      = G4KaonZeroShort::KaonZeroShort();
                                                   >>  72   theK0L      = G4KaonZeroLong::KaonZeroLong();
                                                   >>  73   theL        = G4Lambda::Lambda();
                                                   >>  74   theAntiL    = G4AntiLambda::AntiLambda();
                                                   >>  75   theSPlus    = G4SigmaPlus::SigmaPlus();
                                                   >>  76   theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
                                                   >>  77   theSMinus   = G4SigmaMinus::SigmaMinus();
                                                   >>  78   theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
                                                   >>  79   theS0       = G4SigmaZero::SigmaZero();
                                                   >>  80   theAS0      = G4AntiSigmaZero::AntiSigmaZero();
                                                   >>  81   theXiMinus  = G4XiMinus::XiMinus();
                                                   >>  82   theXi0      = G4XiZero::XiZero();
                                                   >>  83   theAXiMinus = G4AntiXiMinus::AntiXiMinus();
                                                   >>  84   theAXi0     = G4AntiXiZero::AntiXiZero();
                                                   >>  85   theOmega    = G4OmegaMinus::OmegaMinus();
                                                   >>  86   theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
                                                   >>  87   theD        = G4Deuteron::Deuteron();
                                                   >>  88   theT        = G4Triton::Triton();
                                                   >>  89   theA        = G4Alpha::Alpha();
                                                   >>  90   theHe3      = G4He3::He3();
 77 }                                                  91 }
 78                                                    92 
 79 G4ChargeExchange::~G4ChargeExchange()              93 G4ChargeExchange::~G4ChargeExchange()
 80 {                                              <<  94 {}
 81   delete fHandler;                             << 
 82 }                                              << 
 83                                                    95 
 84 G4HadFinalState* G4ChargeExchange::ApplyYourse     96 G4HadFinalState* G4ChargeExchange::ApplyYourself(
 85      const G4HadProjectile& aTrack, G4Nucleus&     97      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 86 {                                                  98 {
 87   theParticleChange.Clear();                       99   theParticleChange.Clear();
 88   auto part = aTrack.GetDefinition();          << 100   const G4HadProjectile* aParticle = &aTrack;
 89   G4double ekin = aTrack.GetKineticEnergy();   << 101   G4double ekin = aParticle->GetKineticEnergy();
 90                                                   102 
 91   G4int A = targetNucleus.GetA_asInt();           103   G4int A = targetNucleus.GetA_asInt();
 92   G4int Z = targetNucleus.GetZ_asInt();           104   G4int Z = targetNucleus.GetZ_asInt();
 93                                                   105 
 94   if (ekin <= lowEnergyLimit) {                << 106   if(ekin <= lowestEnergyLimit || A < 3) {
                                                   >> 107     theParticleChange.SetEnergyChange(ekin);
                                                   >> 108     theParticleChange.SetMomentumChange(0.0,0.0,1.0);
 95     return &theParticleChange;                    109     return &theParticleChange;
 96   }                                               110   }
 97   theParticleChange.SetWeightChange(fXSWeightF << 
 98                                                   111 
 99   G4int projPDG = part->GetPDGEncoding();      << 112   G4double plab = aParticle->GetTotalMomentum();
100                                                   113 
101   // for hydrogen targets and positive project << 
102   // is not possible on proton, only on deuter << 
103   if (1 == Z && (211 == projPDG || 321 == proj << 
104                                                << 
105   if (verboseLevel > 1)                           114   if (verboseLevel > 1)
106     G4cout << "G4ChargeExchange for " << part- << 115     G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
                                                   >> 116      << plab/GeV << " GeV/c "
                                                   >> 117      << " ekin(MeV) = " << ekin/MeV << "  "
                                                   >> 118      << aParticle->GetDefinition()->GetParticleName() << G4endl;
                                                   >> 119 
                                                   >> 120   // Scattered particle referred to axis of incident particle
                                                   >> 121   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
                                                   >> 122 
                                                   >> 123   G4int N = A - Z;
                                                   >> 124   G4int projPDG = theParticle->GetPDGEncoding();
                                                   >> 125   if (verboseLevel > 1)
                                                   >> 126     G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
107      << " PDGcode= " << projPDG << " on nucleu    127      << " PDGcode= " << projPDG << " on nucleus Z= " << Z
108      << " A= " << A << " N= " << A - Z         << 128      << " A= " << A << " N= " << N
109      << G4endl;                                   129      << G4endl;
110                                                   130 
111   G4double mass1 = G4NucleiProperties::GetNucl << 131   const G4ParticleDefinition* theDef = nullptr;
112   G4LorentzVector lv0 = aTrack.Get4Momentum(); << 
113   G4double etot = mass1 + lv0.e();             << 
114                                                << 
115   // select final state                        << 
116   const G4ParticleDefinition* theSecondary =   << 
117     fXSection->SampleSecondaryType(part, Z, A) << 
118   G4int pdg = theSecondary->GetPDGEncoding();  << 
119                                                << 
120   // omega(782) and f2(1270)                   << 
121   G4bool isShortLived = (pdg == 223 || pdg ==  << 
122                                                << 
123   // atomic number of the recoil nucleus       << 
124   if (projPDG == -211) { --Z; }                << 
125   else if (projPDG == 211) { ++Z; }            << 
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                                                   132 
136   // recoil nucleus                            << 133   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
                                                   >> 134   G4LorentzVector lv1 = aParticle->Get4Momentum();
                                                   >> 135   G4LorentzVector lv0(0.0,0.0,0.0,mass2);
                                                   >> 136 
                                                   >> 137   G4LorentzVector lv  = lv0 + lv1;
                                                   >> 138   G4ThreeVector bst = lv.boostVector();
                                                   >> 139   lv1.boost(-bst);
                                                   >> 140   lv0.boost(-bst);
                                                   >> 141 
                                                   >> 142   // Sample final particles
                                                   >> 143   G4bool theHyperon = false;
137   const G4ParticleDefinition* theRecoil = null    144   const G4ParticleDefinition* theRecoil = nullptr;
138   if (Z == 0 && A == 1) { theRecoil = G4Neutro << 145   const G4ParticleDefinition* theSecondary = nullptr;
139   else if (Z == 1 && A == 1) { theRecoil = G4P << 
140   else if (Z == 1 && A == 2) { theRecoil = G4D << 
141   else if (Z == 1 && A == 3) { theRecoil = G4T << 
142   else if (Z == 2 && A == 3) { theRecoil = G4H << 
143   else if (Z == 2 && A == 4) { theRecoil = G4A << 
144   else if (nist->GetIsotopeAbundance(Z, A) > 0 << 
145     theRecoil = G4ParticleTable::GetParticleTa << 
146       ->GetIonTable()->GetIon(Z, A, 0.0);      << 
147   }                                            << 
148                                                   146 
149   // check if there is enough energy for the f << 147   if(theParticle == theProton) {
150   // and sample mass of produced state         << 148     theSecondary = theNeutron;
151   const G4double mass0 = theSecondary->GetPDGM << 149     Z++;
152   G4double mass3 = (nullptr == theRecoil) ?    << 150   } else if(theParticle == theNeutron) {
153     G4NucleiProperties::GetNuclearMass(A, Z) : << 151     theSecondary = theProton;
154   G4double mass2 = mass0;                      << 152     Z--;
155   if (isShortLived &&                          << 153   } else if(theParticle == thePiPlus) {
156       !SampleMass(mass2, theSecondary->GetPDGW << 154     theSecondary = thePiZero;
157     return &theParticleChange;                 << 155     Z++;
                                                   >> 156   } else if(theParticle == thePiMinus) {
                                                   >> 157     theSecondary = thePiZero;
                                                   >> 158     Z--;
                                                   >> 159   } else if(theParticle == theKPlus) {
                                                   >> 160     if(G4UniformRand()<0.5) theSecondary = theK0S;
                                                   >> 161     else  theSecondary = theK0L;
                                                   >> 162     Z++;
                                                   >> 163   } else if(theParticle == theKMinus) {
                                                   >> 164     if(G4UniformRand()<0.5) theSecondary = theK0S;
                                                   >> 165     else  theSecondary = theK0L;
                                                   >> 166     Z--;
                                                   >> 167   } else if(theParticle == theK0S || theParticle == theK0L) {
                                                   >> 168     if(G4UniformRand()*A < G4double(Z)) {
                                                   >> 169       theSecondary = theKPlus;
                                                   >> 170       Z--;
                                                   >> 171     } else {
                                                   >> 172       theSecondary = theKMinus;
                                                   >> 173       Z++;
                                                   >> 174     }
                                                   >> 175   } else if(theParticle == theANeutron) {
                                                   >> 176     theSecondary = theAProton;
                                                   >> 177     Z++;
                                                   >> 178   } else if(theParticle == theAProton) {
                                                   >> 179     theSecondary = theANeutron;
                                                   >> 180     Z--;
                                                   >> 181   } else if(theParticle == theL) {
                                                   >> 182     G4double x = G4UniformRand();
                                                   >> 183     if(G4UniformRand()*A < G4double(Z)) {
                                                   >> 184       if(x < 0.2) {
                                                   >> 185         theSecondary = theS0;
                                                   >> 186       } else if (x < 0.4) {
                                                   >> 187         theSecondary = theSPlus;
                                                   >> 188         Z--;
                                                   >> 189       } else if (x < 0.6) {
                                                   >> 190         theSecondary = theProton;
                                                   >> 191   theRecoil = theL;
                                                   >> 192         theHyperon = true;
                                                   >> 193   A--;
                                                   >> 194       } else if (x < 0.8) {
                                                   >> 195         theSecondary = theProton;
                                                   >> 196   theRecoil = theS0;
                                                   >> 197         theHyperon = true;
                                                   >> 198         A--;
                                                   >> 199       } else {
                                                   >> 200         theSecondary = theNeutron;
                                                   >> 201   theRecoil = theSPlus;
                                                   >> 202         theHyperon = true;
                                                   >> 203         A--;
                                                   >> 204       }
                                                   >> 205     } else {
                                                   >> 206       if(x < 0.2) {
                                                   >> 207         theSecondary = theS0;
                                                   >> 208       } else if (x < 0.4) {
                                                   >> 209         theSecondary = theSMinus;
                                                   >> 210         Z++;
                                                   >> 211       } else if (x < 0.6) {
                                                   >> 212         theSecondary = theNeutron;
                                                   >> 213   theRecoil = theL;
                                                   >> 214         A--;
                                                   >> 215         theHyperon = true;
                                                   >> 216       } else if (x < 0.8) {
                                                   >> 217         theSecondary = theNeutron;
                                                   >> 218   theRecoil = theS0;
                                                   >> 219         theHyperon = true;
                                                   >> 220         A--;
                                                   >> 221       } else {
                                                   >> 222         theSecondary = theProton;
                                                   >> 223   theRecoil = theSMinus;
                                                   >> 224         theHyperon = true;
                                                   >> 225         A--;
                                                   >> 226       }
                                                   >> 227     }
158   }                                               228   }
159                                                   229 
160   // not possible kinematically                << 230   if (Z == 1 && A == 2) theDef = theD;
161   if (etot <= mass2 + mass3) {                 << 231   else if (Z == 1 && A == 3) theDef = theT;
                                                   >> 232   else if (Z == 2 && A == 3) theDef = theHe3;
                                                   >> 233   else if (Z == 2 && A == 4) theDef = theA;
                                                   >> 234   else {
                                                   >> 235     theDef = 
                                                   >> 236       G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
                                                   >> 237   }
                                                   >> 238   if(!theSecondary) { return &theParticleChange; }
                                                   >> 239 
                                                   >> 240   G4double m11 = theSecondary->GetPDGMass();
                                                   >> 241   G4double m21 = theDef->GetPDGMass();
                                                   >> 242   if(theRecoil)  { m21 += theRecoil->GetPDGMass(); }
                                                   >> 243   else           { theRecoil = theDef; }
                                                   >> 244 
                                                   >> 245   G4double etot = lv0.e() + lv1.e();
                                                   >> 246 
                                                   >> 247   // kinematiacally impossible
                                                   >> 248   if(etot < m11 + m21) {
                                                   >> 249     theParticleChange.SetEnergyChange(ekin);
                                                   >> 250     theParticleChange.SetMomentumChange(0.0,0.0,1.0);
162     return &theParticleChange;                    251     return &theParticleChange;
163   }                                               252   }
164                                                   253 
165   // sample kinematics                         << 254   G4ThreeVector p1 = lv1.vect();
166   G4LorentzVector lv1(0.0, 0.0, 0.0, mass1);   << 255   G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
167   G4LorentzVector lv = lv0 + lv1;              << 256   //  G4double e2 = etot - e1;
168   G4ThreeVector bst = lv.boostVector();        << 257   G4double ptot = std::sqrt(e1*e1 - m11*m11);
169   G4double ss = lv.mag2();                     << 258 
170                                                << 259   G4double tmax = 4.0*ptot*ptot;
171   // tmax = 4*momCMS^2                         << 260   G4double g2 = GeV*GeV; 
172   G4double e2 = ss + mass2*mass2 - mass3*mass3 << 261 
173   G4double tmax = e2*e2/ss - 4*mass2*mass2;    << 262   G4double t = g2*SampleT(tmax/g2, A);
174                                                << 263 
175   G4double t = SampleT(theSecondary, A, tmax); << 264   if(verboseLevel>1) {
176                                                << 265     G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
177   G4double phi  = G4UniformRand()*CLHEP::twopi << 266      << " ptot= " << ptot << G4endl;
                                                   >> 267   }
                                                   >> 268   // Sampling in CM system
                                                   >> 269   G4double phi  = G4UniformRand()*twopi;
178   G4double cost = 1. - 2.0*t/tmax;                270   G4double cost = 1. - 2.0*t/tmax;
                                                   >> 271   if(std::abs(cost) > 1.0) cost = 1.0;
                                                   >> 272   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
179                                                   273 
180   if (cost > 1.0) { cost = 1.0; }              << 274   //if (verboseLevel > 1)
181   else if(cost < -1.0) { cost = -1.0; }        << 275   //  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
182                                                   276 
183   G4double sint = std::sqrt((1.0-cost)*(1.0+co << 277   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
                                                   >> 278   v1 *= ptot;
                                                   >> 279   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
                                                   >> 280   G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
184                                                   281 
185   if (verboseLevel>1) {                        << 282   nlv0.boost(bst);
186     G4cout << " t= " << t << " tmax(GeV^2)= "  << 283   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                                                   284 
205   // prepare secondary particles               << 
206   theParticleChange.SetStatusChange(stopAndKil    285   theParticleChange.SetStatusChange(stopAndKill);
207   theParticleChange.SetEnergyChange(0.0);         286   theParticleChange.SetEnergyChange(0.0);
                                                   >> 287   G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
                                                   >> 288   theParticleChange.AddSecondary(aSec);
208                                                   289 
209   if (!isShortLived) {                         << 290   G4double erec = std::max(nlv0.e() - m21, 0.0);
210     auto aSec = new G4DynamicParticle(theSecon << 291 
211     theParticleChange.AddSecondary(aSec, secID << 292   //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;  
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                                                   293 
227   // recoil is a stable isotope                << 294   if(theHyperon) {
228   if (nullptr != theRecoil) {                  << 295     theParticleChange.SetLocalEnergyDeposit(erec);
229     auto aRec = new G4DynamicParticle(theRecoi << 296     aSec = new G4DynamicParticle();
230     theParticleChange.AddSecondary(aRec, secID << 297     aSec->SetDefinition(theRecoil);
                                                   >> 298     aSec->SetKineticEnergy(0.0);
                                                   >> 299   } else if(erec > GetRecoilEnergyThreshold()) {
                                                   >> 300     aSec = new G4DynamicParticle(theRecoil, nlv0);
                                                   >> 301     theParticleChange.AddSecondary(aSec);
231   } else {                                        302   } else {
232     // recoil is an unstable fragment          << 303     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   }                                               304   }
242   return &theParticleChange;                      305   return &theParticleChange;
243 }                                                 306 }
244                                                   307 
245 G4double G4ChargeExchange::SampleT(const G4Par << 308 G4double G4ChargeExchange::SampleT(G4double tmax, G4int A)
246                                    const G4int << 
247 {                                                 309 {
248   G4double aa, bb, cc, dd;                        310   G4double aa, bb, cc, dd;
249   G4Pow* g4pow = G4Pow::GetInstance();            311   G4Pow* g4pow = G4Pow::GetInstance();
250   if (A <= 62.) {                                 312   if (A <= 62.) {
251     aa = g4pow->powZ(A, 1.63);                    313     aa = g4pow->powZ(A, 1.63);
252     bb = 14.5*g4pow->powZ(A, 0.66);               314     bb = 14.5*g4pow->powZ(A, 0.66);
253     cc = 1.4*g4pow->powZ(A, 0.33);                315     cc = 1.4*g4pow->powZ(A, 0.33);
254     dd = 10.;                                     316     dd = 10.;
255   } else {                                        317   } else {
256     aa = g4pow->powZ(A, 1.33);                    318     aa = g4pow->powZ(A, 1.33);
257     bb = 60.*g4pow->powZ(A, 0.33);                319     bb = 60.*g4pow->powZ(A, 0.33);
258     cc = 0.4*g4pow->powZ(A, 0.40);                320     cc = 0.4*g4pow->powZ(A, 0.40);
259     dd = 10.;                                     321     dd = 10.;
260   }                                               322   }
261   G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb;    323   G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb;
262   G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd;    324   G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd;
263                                                   325   
264   G4double t;                                     326   G4double t;
265   G4double y = bb;                                327   G4double y = bb;
266   if(G4UniformRand()*(x1 + x2) < x2) y = dd;      328   if(G4UniformRand()*(x1 + x2) < x2) y = dd;
267                                                   329 
268   for (G4int i=0; i<maxN; ++i) {               << 330   const G4int maxNumberOfLoops = 10000;
                                                   >> 331   G4int loopCounter = 0;
                                                   >> 332   do {
269     t = -G4Log(G4UniformRand())/y;                333     t = -G4Log(G4UniformRand())/y;
270     if (t <= tmax) { return t; }               << 334   } while ( (t > tmax) &&
                                                   >> 335             ++loopCounter < maxNumberOfLoops );  /* Loop checking, 10.08.2015, A.Ribon */
                                                   >> 336   if ( loopCounter >= maxNumberOfLoops ) {
                                                   >> 337     t = 0.0;
271   }                                               338   }
272   return 0.0;                                  << 339   return t;
273 }                                                 340 }
274                                                   341 
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 }                                              << 
296                                                   342