Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4ChipsHyperonInelasticXS.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/cross_sections/src/G4ChipsHyperonInelasticXS.cc (Version 11.3.0) and /processes/hadronic/cross_sections/src/G4ChipsHyperonInelasticXS.cc (Version 10.0.p1)


  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 //                                                 26 //
 27 // The lust update: M.V. Kossov, CERN/ITEP(Mos     27 // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
 28 //                                                 28 //
 29 // *******************************************     29 // ****************************************************************************************
 30 // Short description: Cross-sections extracted     30 // Short description: Cross-sections extracted (by W.Pokorski) from the CHIPS package for 
 31 // Hyperon-nuclear  interactions. Original aut     31 // Hyperon-nuclear  interactions. Original author: M. Kossov
 32 // -------------------------------------------     32 // -------------------------------------------------------------------------------------
 33 //                                                 33 //
 34                                                    34 
 35 #include "G4ChipsHyperonInelasticXS.hh"            35 #include "G4ChipsHyperonInelasticXS.hh"
 36 #include "G4SystemOfUnits.hh"                      36 #include "G4SystemOfUnits.hh"
 37 #include "G4DynamicParticle.hh"                    37 #include "G4DynamicParticle.hh"
 38 #include "G4ParticleDefinition.hh"                 38 #include "G4ParticleDefinition.hh"
 39 #include "G4Lambda.hh"                             39 #include "G4Lambda.hh"
 40 #include "G4SigmaPlus.hh"                          40 #include "G4SigmaPlus.hh"
 41 #include "G4SigmaMinus.hh"                         41 #include "G4SigmaMinus.hh"
 42 #include "G4SigmaZero.hh"                          42 #include "G4SigmaZero.hh"
 43 #include "G4XiMinus.hh"                            43 #include "G4XiMinus.hh"
 44 #include "G4XiZero.hh"                             44 #include "G4XiZero.hh"
 45 #include "G4OmegaMinus.hh"                         45 #include "G4OmegaMinus.hh"
 46 #include "G4Log.hh"                            << 
 47 #include "G4Exp.hh"                            << 
 48                                                    46 
 49 // factory                                         47 // factory
 50 #include "G4CrossSectionFactory.hh"                48 #include "G4CrossSectionFactory.hh"
 51 //                                                 49 //
 52 G4_DECLARE_XS_FACTORY(G4ChipsHyperonInelasticX     50 G4_DECLARE_XS_FACTORY(G4ChipsHyperonInelasticXS);
 53                                                    51 
 54 G4ChipsHyperonInelasticXS::G4ChipsHyperonInela     52 G4ChipsHyperonInelasticXS::G4ChipsHyperonInelasticXS():G4VCrossSectionDataSet(Default_Name())
 55 {                                                  53 {
 56   // Initialization of the                         54   // Initialization of the
 57   lastLEN=0; // Pointer to the lastArray of Lo     55   lastLEN=0; // Pointer to the lastArray of LowEn CS
 58   lastHEN=0; // Pointer to the lastArray of Hi     56   lastHEN=0; // Pointer to the lastArray of HighEn CS
 59   lastN=0;   // The last N of calculated nucle     57   lastN=0;   // The last N of calculated nucleus
 60   lastZ=0;   // The last Z of calculated nucle     58   lastZ=0;   // The last Z of calculated nucleus
 61   lastP=0.;  // Last used in cross section Mom     59   lastP=0.;  // Last used in cross section Momentum
 62   lastTH=0.; // Last threshold momentum            60   lastTH=0.; // Last threshold momentum
 63   lastCS=0.; // Last value of the Cross Sectio     61   lastCS=0.; // Last value of the Cross Section
 64   lastI=0;   // The last position in the DAMDB     62   lastI=0;   // The last position in the DAMDB
 65   LEN = new std::vector<G4double*>;                63   LEN = new std::vector<G4double*>;
 66   HEN = new std::vector<G4double*>;                64   HEN = new std::vector<G4double*>;
 67 }                                                  65 }
 68                                                    66 
 69 G4ChipsHyperonInelasticXS::~G4ChipsHyperonInel     67 G4ChipsHyperonInelasticXS::~G4ChipsHyperonInelasticXS()
 70 {                                                  68 {
 71     std::size_t lens=LEN->size();              <<  69     G4int lens=LEN->size();
 72     for(std::size_t i=0; i<lens; ++i) delete[] <<  70     for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
 73     delete LEN;                                    71     delete LEN;
 74                                                    72 
 75     std::size_t hens=HEN->size();              <<  73     G4int hens=HEN->size();
 76     for(std::size_t i=0; i<hens; ++i) delete[] <<  74     for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
 77     delete HEN;                                    75     delete HEN;
 78 }                                                  76 }
 79                                                    77 
 80 void G4ChipsHyperonInelasticXS::CrossSectionDe <<  78 G4bool G4ChipsHyperonInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,    
 81 {                                              << 
 82   outFile << "G4ChipsHyperonInelasticXS provid << 
 83           << "section for hyperon nucleus scat << 
 84           << "momentum. The cross section is c << 
 85           << "CHIPS parameterization of cross  << 
 86 }                                              << 
 87                                                << 
 88 G4bool G4ChipsHyperonInelasticXS::IsIsoApplica << 
 89          const G4Element*,                         79          const G4Element*,
 90          const G4Material*)                        80          const G4Material*)
 91 {                                                  81 {
 92   return true;                                 <<  82   G4ParticleDefinition* particle = Pt->GetDefinition();
                                                   >>  83   if (particle == G4Lambda::Lambda()) 
                                                   >>  84     {
                                                   >>  85       return true;
                                                   >>  86     }
                                                   >>  87   else if(particle == G4SigmaPlus::SigmaPlus())
                                                   >>  88     {
                                                   >>  89     return true;
                                                   >>  90     }
                                                   >>  91   else if(particle == G4SigmaMinus::SigmaMinus())
                                                   >>  92     {
                                                   >>  93     return true;
                                                   >>  94     }
                                                   >>  95   else if(particle == G4SigmaZero::SigmaZero())
                                                   >>  96     {
                                                   >>  97       return true;
                                                   >>  98     }
                                                   >>  99   else if(particle == G4XiMinus::XiMinus())
                                                   >> 100     {
                                                   >> 101       return true;
                                                   >> 102     }
                                                   >> 103   else if(particle == G4XiZero::XiZero())
                                                   >> 104     {
                                                   >> 105       return true;
                                                   >> 106     }
                                                   >> 107   else if(particle == G4OmegaMinus::OmegaMinus())
                                                   >> 108     {
                                                   >> 109       return true;
                                                   >> 110     }
                                                   >> 111   return false;
 93 }                                                 112 }
 94                                                   113 
 95 // The main member function giving the collisi    114 // The main member function giving the collision cross section (P is in IU, CS is in mb)
 96 // Make pMom in independent units ! (Now it is    115 // Make pMom in independent units ! (Now it is MeV)
 97 G4double G4ChipsHyperonInelasticXS::GetIsoCros    116 G4double G4ChipsHyperonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,  
 98                const G4Isotope*,                  117                const G4Isotope*,
 99                const G4Element*,                  118                const G4Element*,
100                const G4Material*)                 119                const G4Material*)
101 {                                                 120 {
102   G4double pMom=Pt->GetTotalMomentum();           121   G4double pMom=Pt->GetTotalMomentum();
103   G4int tgN = A - tgZ;                            122   G4int tgN = A - tgZ;
104   G4int pdg = Pt->GetDefinition()->GetPDGEncod    123   G4int pdg = Pt->GetDefinition()->GetPDGEncoding();
105                                                   124   
106   return GetChipsCrossSection(pMom, tgZ, tgN,     125   return GetChipsCrossSection(pMom, tgZ, tgN, pdg);
107 }                                                 126 }
108                                                   127 
109 G4double G4ChipsHyperonInelasticXS::GetChipsCr    128 G4double G4ChipsHyperonInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int PDG)
110 {                                                 129 {
                                                   >> 130   static G4ThreadLocal G4int j;                      // A#0f Z/N-records already tested in AMDB
                                                   >> 131   static G4ThreadLocal std::vector <G4int>    *colN_G4MT_TLS_ = 0 ; if (!colN_G4MT_TLS_) colN_G4MT_TLS_ = new  std::vector <G4int>     ;  std::vector <G4int>    &colN = *colN_G4MT_TLS_;  // Vector of N for calculated nuclei (isotops)
                                                   >> 132   static G4ThreadLocal std::vector <G4int>    *colZ_G4MT_TLS_ = 0 ; if (!colZ_G4MT_TLS_) colZ_G4MT_TLS_ = new  std::vector <G4int>     ;  std::vector <G4int>    &colZ = *colZ_G4MT_TLS_;  // Vector of Z for calculated nuclei (isotops)
                                                   >> 133   static G4ThreadLocal std::vector <G4double> *colP_G4MT_TLS_ = 0 ; if (!colP_G4MT_TLS_) colP_G4MT_TLS_ = new  std::vector <G4double>  ;  std::vector <G4double> &colP = *colP_G4MT_TLS_;  // Vector of last momenta for the reaction
                                                   >> 134   static G4ThreadLocal std::vector <G4double> *colTH_G4MT_TLS_ = 0 ; if (!colTH_G4MT_TLS_) colTH_G4MT_TLS_ = new  std::vector <G4double>  ;  std::vector <G4double> &colTH = *colTH_G4MT_TLS_; // Vector of energy thresholds for the reaction
                                                   >> 135   static G4ThreadLocal std::vector <G4double> *colCS_G4MT_TLS_ = 0 ; if (!colCS_G4MT_TLS_) colCS_G4MT_TLS_ = new  std::vector <G4double>  ;  std::vector <G4double> &colCS = *colCS_G4MT_TLS_; // Vector of last cross sections for the reaction
                                                   >> 136   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
111                                                   137  
112   G4bool in=false;                     // By d    138   G4bool in=false;                     // By default the isotope must be found in the AMDB
113   if(tgN!=lastN || tgZ!=lastZ)         // The     139   if(tgN!=lastN || tgZ!=lastZ)         // The nucleus was not the last used isotope
114   {                                               140   {
115     in = false;                        // By d    141     in = false;                        // By default the isotope haven't be found in AMDB  
116     lastP   = 0.;                      // New     142     lastP   = 0.;                      // New momentum history (nothing to compare with)
117     lastN   = tgN;                     // The     143     lastN   = tgN;                     // The last N of the calculated nucleus
118     lastZ   = tgZ;                     // The     144     lastZ   = tgZ;                     // The last Z of the calculated nucleus
119     lastI   = (G4int)colN.size();      // Size << 145     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
120     j  = 0;                            // A#0f    146     j  = 0;                            // A#0f records found in DB for this projectile
121                                                   147 
122     if(lastI) for(G4int i=0; i<lastI; ++i) //  << 148     if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope
123     {                                             149     {
124       if(colN[i]==tgN && colZ[i]==tgZ) // Try     150       if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
125       {                                           151       {
126         lastI=i;                       // Reme    152         lastI=i;                       // Remember the index for future fast/last use
127         lastTH =colTH[i];              // The     153         lastTH =colTH[i];              // The last THreshold (A-dependent)
128                                                   154 
129         if(pMom<=lastTH)                          155         if(pMom<=lastTH)
130         {                                         156         {
131           return 0.;                   // Ener    157           return 0.;                   // Energy is below the Threshold value
132         }                                         158         }
133         lastP  =colP [i];              // Last    159         lastP  =colP [i];              // Last Momentum  (A-dependent)
134         lastCS =colCS[i];              // Last    160         lastCS =colCS[i];              // Last CrossSect (A-dependent)
135         in = true;                     // This    161         in = true;                     // This is the case when the isotop is found in DB
136         // Momentum pMom is in IU ! @@ Units      162         // Momentum pMom is in IU ! @@ Units
137         lastCS=CalculateCrossSection(-1,j,PDG,    163         lastCS=CalculateCrossSection(-1,j,PDG,lastZ,lastN,pMom); // read & update
138                                                   164 
139         if(lastCS<=0. && pMom>lastTH)  // Corr    165         if(lastCS<=0. && pMom>lastTH)  // Correct the threshold (@@ No intermediate Zeros)
140         {                                         166         {
141           lastCS=0.;                              167           lastCS=0.;
142           lastTH=pMom;                            168           lastTH=pMom;
143         }                                         169         }
144         break;                         // Go o    170         break;                         // Go out of the LOOP
145       }                                           171       }
146       j++;                             // Incr    172       j++;                             // Increment a#0f records found in DB
147     }                                             173     }
148     if(!in)                            // This    174     if(!in)                            // This isotope has not been calculated previously
149     {                                             175     {
150       //!!The slave functions must provide cro    176       //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
151       lastCS=CalculateCrossSection(0,j,PDG,las    177       lastCS=CalculateCrossSection(0,j,PDG,lastZ,lastN,pMom); //calculate & create
                                                   >> 178       //if(lastCS>0.)                   // It means that the AMBD was initialized
                                                   >> 179       //{
152                                                   180 
153       lastTH = 0; //ThresholdEnergy(tgZ, tgN);    181       lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
154         colN.push_back(tgN);                      182         colN.push_back(tgN);
155         colZ.push_back(tgZ);                      183         colZ.push_back(tgZ);
156         colP.push_back(pMom);                     184         colP.push_back(pMom);
157         colTH.push_back(lastTH);                  185         colTH.push_back(lastTH);
158         colCS.push_back(lastCS);                  186         colCS.push_back(lastCS);
159   //} // M.K. Presence of H1 with high thresho    187   //} // M.K. Presence of H1 with high threshold breaks the syncronization
160       return lastCS*millibarn;                    188       return lastCS*millibarn;
161     } // End of creation of the new set of par    189     } // End of creation of the new set of parameters
162     else                                          190     else
163     {                                             191     {
164       colP[lastI]=pMom;                           192       colP[lastI]=pMom;
165       colCS[lastI]=lastCS;                        193       colCS[lastI]=lastCS;
166     }                                             194     }
167   } // End of parameters udate                    195   } // End of parameters udate
168   else if(pMom<=lastTH)                           196   else if(pMom<=lastTH)
169   {                                               197   {
170     return 0.;                         // Mome    198     return 0.;                         // Momentum is below the Threshold Value -> CS=0
171   }                                               199   }
172   else                                 // It i    200   else                                 // It is the last used -> use the current tables
173   {                                               201   {
174     lastCS=CalculateCrossSection(1,j,PDG,lastZ    202     lastCS=CalculateCrossSection(1,j,PDG,lastZ,lastN,pMom); // Only read and UpdateDB
175     lastP=pMom;                                   203     lastP=pMom;
176   }                                               204   }
177   return lastCS*millibarn;                        205   return lastCS*millibarn;
178 }                                                 206 }
179                                                   207 
180 // The main member function giving the gamma-A    208 // The main member function giving the gamma-A cross section (E in GeV, CS in mb)
181 G4double G4ChipsHyperonInelasticXS::CalculateC    209 G4double G4ChipsHyperonInelasticXS::CalculateCrossSection(G4int F, G4int I,
182                 G4int, G4int targZ, G4int targ    210                 G4int, G4int targZ, G4int targN, G4double Momentum)
183 {                                                 211 {
184   static const G4double THmin=27.;     // defa    212   static const G4double THmin=27.;     // default minimum Momentum (MeV/c) Threshold
185   static const G4double THmiG=THmin*.001; // m    213   static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
186   static const G4double dP=10.;        // step    214   static const G4double dP=10.;        // step for the LEN (Low ENergy) table MeV/c
187   static const G4double dPG=dP*.001;   // step    215   static const G4double dPG=dP*.001;   // step for the LEN (Low ENergy) table GeV/c
188   static const G4int    nL=105;        // A#of    216   static const G4int    nL=105;        // A#of LEN points in E (step 10 MeV/c)
189   static const G4double Pmin=THmin+(nL-1)*dP;     217   static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
190   static const G4double Pmax=227000.;  // maxP    218   static const G4double Pmax=227000.;  // maxP for the HEN (High ENergy) part 227 GeV
191   static const G4int    nH=224;        // A#of    219   static const G4int    nH=224;        // A#of HEN points in lnE
192   static const G4double milP=G4Log(Pmin);// Lo << 220   static const G4double milP=std::log(Pmin);// Low logarithm energy for the HEN part
193   static const G4double malP=G4Log(Pmax);// Hi << 221   static const G4double malP=std::log(Pmax);// High logarithm energy (each 2.75 percent)
194   static const G4double dlP=(malP-milP)/(nH-1)    222   static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
195   static const G4double milPG=G4Log(.001*Pmin) << 223   static const G4double milPG=std::log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
196                                                << 224   G4double sigma=0.;
                                                   >> 225   if(F&&I) sigma=0.;                   // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
                                                   >> 226   //G4double A=targN+targZ;              // A of the target
197   if(F<=0)                             // This    227   if(F<=0)                             // This isotope was not the last used isotop
198   {                                               228   {
199     if(F<0)                            // This    229     if(F<0)                            // This isotope was found in DAMDB =-----=> RETRIEVE
200     {                                             230     {
201       G4int sync=(G4int)LEN->size();           << 231       G4int sync=LEN->size();
202       if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS    232       if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
203       lastLEN=(*LEN)[I];               // Poin    233       lastLEN=(*LEN)[I];               // Pointer to prepared LowEnergy cross sections
204       lastHEN=(*HEN)[I];               // Poin    234       lastHEN=(*HEN)[I];               // Pointer to prepared High Energy cross sections
205     }                                             235     }
206     else                               // This    236     else                               // This isotope wasn't calculated before => CREATE
207     {                                             237     {
208       lastLEN = new G4double[nL];      // Allo    238       lastLEN = new G4double[nL];      // Allocate memory for the new LEN cross sections
209       lastHEN = new G4double[nH];      // Allo    239       lastHEN = new G4double[nH];      // Allocate memory for the new HEN cross sections
210       // --- Instead of making a separate func    240       // --- Instead of making a separate function ---
211       G4double P=THmiG;                // Tabl    241       G4double P=THmiG;                // Table threshold in GeV/c
212       for(G4int k=0; k<nL; k++)                   242       for(G4int k=0; k<nL; k++)
213       {                                           243       {
214         lastLEN[k] = CrossSectionLin(targZ, ta    244         lastLEN[k] = CrossSectionLin(targZ, targN, P);
215         P+=dPG;                                   245         P+=dPG;
216       }                                           246       }
217       G4double lP=milPG;                          247       G4double lP=milPG;
218       for(G4int n=0; n<nH; n++)                   248       for(G4int n=0; n<nH; n++)
219       {                                           249       {
220         lastHEN[n] = CrossSectionLog(targZ, ta    250         lastHEN[n] = CrossSectionLog(targZ, targN, lP);
221         lP+=dlP;                                  251         lP+=dlP;
222       }                                           252       }
223       // --- End of possible separate function    253       // --- End of possible separate function
224       // *** The synchronization check ***        254       // *** The synchronization check ***
225       G4int sync=(G4int)LEN->size();           << 255       G4int sync=LEN->size();
226       if(sync!=I)                                 256       if(sync!=I)
227       {                                           257       {
228         G4cerr<<"***G4QHyperNuclCS::CalcCrossS    258         G4cerr<<"***G4QHyperNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
229               <<", N="<<targN<<", F="<<F<<G4en    259               <<", N="<<targN<<", F="<<F<<G4endl;
230         //G4Exception("G4PiMinusNuclearCS::Cal    260         //G4Exception("G4PiMinusNuclearCS::CalculateCS:","39",FatalException,"DBoverflow");
231       }                                           261       }
232       LEN->push_back(lastLEN);         // reme    262       LEN->push_back(lastLEN);         // remember the Low Energy Table
233       HEN->push_back(lastHEN);         // reme    263       HEN->push_back(lastHEN);         // remember the High Energy Table
234     } // End of creation of the new set of par    264     } // End of creation of the new set of parameters
235   } // End of parameters udate                    265   } // End of parameters udate
236   // =--------------------------= NOW the Magi    266   // =--------------------------= NOW the Magic Formula =------------------------------=
237   G4double sigma;                              << 
238   if (Momentum<lastTH) return 0.;      // It m    267   if (Momentum<lastTH) return 0.;      // It must be already checked in the interface class
239   else if (Momentum<Pmin)              // High    268   else if (Momentum<Pmin)              // High Energy region
240   {                                               269   {
241     sigma=EquLinearFit(Momentum,nL,THmin,dP,la    270     sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
242   }                                               271   }
243   else if (Momentum<Pmax)              // High    272   else if (Momentum<Pmax)              // High Energy region
244   {                                               273   {
245     G4double lP=G4Log(Momentum);               << 274     G4double lP=std::log(Momentum);
246     sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN)    275     sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
247   }                                               276   }
248   else                                 // UHE     277   else                                 // UHE region (calculation, not frequent)
249   {                                               278   {
250     G4double P=0.001*Momentum;         // Appr    279     G4double P=0.001*Momentum;         // Approximation formula is for P in GeV/c
251     sigma=CrossSectionFormula(targZ, targN, P, << 280     sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
252   }                                               281   }
253   if (sigma<0.) return 0.;                     << 282   if(sigma<0.) return 0.;
254   return sigma;                                   283   return sigma;
255 }                                                 284 }
256                                                   285 
257 // Calculation formula for piMinus-nuclear ine    286 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) (P in GeV/c)
258 G4double G4ChipsHyperonInelasticXS::CrossSecti    287 G4double G4ChipsHyperonInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
259 {                                                 288 {
260   G4double lP=G4Log(P);                        << 289   G4double lP=std::log(P);
261   return CrossSectionFormula(tZ, tN, P, lP);      290   return CrossSectionFormula(tZ, tN, P, lP);
262 }                                                 291 }
263                                                   292 
264 // Calculation formula for piMinus-nuclear ine    293 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
265 G4double G4ChipsHyperonInelasticXS::CrossSecti    294 G4double G4ChipsHyperonInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
266 {                                                 295 {
267   G4double P=G4Exp(lP);                        << 296   G4double P=std::exp(lP);
268   return CrossSectionFormula(tZ, tN, P, lP);      297   return CrossSectionFormula(tZ, tN, P, lP);
269 }                                                 298 }
270 // Calculation formula for piMinus-nuclear ine    299 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
271 G4double G4ChipsHyperonInelasticXS::CrossSecti    300 G4double G4ChipsHyperonInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
272                                                   301                                                               G4double P, G4double lP)
273 {                                                 302 {
274   G4double sigma=0.;                              303   G4double sigma=0.;
275                                                << 
276   //AR-24Apr2018 Switch to allow transuranic e << 
277   const G4bool isHeavyElementAllowed = true;   << 
278                                                << 
279   if(tZ==1 && !tN)                        // H    304   if(tZ==1 && !tN)                        // Hyperon-P interaction from G4QuasiElastRatios
280   {                                               305   {
281     G4double ld=lP-3.5;                           306     G4double ld=lP-3.5;
282     G4double ld2=ld*ld;                           307     G4double ld2=ld*ld;
283     G4double p2=P*P;                              308     G4double p2=P*P;
284     G4double p4=p2*p2;                            309     G4double p4=p2*p2;
285     G4double sp=std::sqrt(P);                     310     G4double sp=std::sqrt(P);
286     G4double El=(.0557*ld2+6.72+99./p2)/(1.+2.    311     G4double El=(.0557*ld2+6.72+99./p2)/(1.+2./sp+2./p4);
287     G4double To=(.3*ld2+38.2+900./sp)/(1.+27./    312     G4double To=(.3*ld2+38.2+900./sp)/(1.+27./sp+3./p4);
288     sigma=To-El;                                  313     sigma=To-El;
289   }                                               314   }
290   else if((tZ<97 && tN<152) || isHeavyElementA << 315   else if(tZ<97 && tN<152)                // General solution
291   {                                               316   {
292     G4double d=lP-4.2;                            317     G4double d=lP-4.2;
293     G4double p2=P*P;                              318     G4double p2=P*P;
294     G4double p4=p2*p2;                            319     G4double p4=p2*p2;
295     G4double sp=std::sqrt(P);                     320     G4double sp=std::sqrt(P);
296     G4double ssp=std::sqrt(sp);                   321     G4double ssp=std::sqrt(sp);
297     G4double a=tN+tZ;                      //     322     G4double a=tN+tZ;                      // A of the target
298     G4double al=G4Log(a);                      << 323     G4double al=std::log(a);
299     G4double sa=std::sqrt(a);                     324     G4double sa=std::sqrt(a);
300     G4double a2=a*a;                              325     G4double a2=a*a;
301     G4double a2s=a2*sa;                           326     G4double a2s=a2*sa;
302     G4double a4=a2*a2;                            327     G4double a4=a2*a2;
303     G4double a8=a4*a4;                            328     G4double a8=a4*a4;
304     G4double c=(170.+3600./a2s)/(1.+65./a2s);     329     G4double c=(170.+3600./a2s)/(1.+65./a2s);
305     G4double gg=42.*(G4Exp(al*0.8)+4.E-8*a4)/( << 330     G4double gg=42.*(std::exp(al*0.8)+4.E-8*a4)/(1.+28./a)/(1.+5.E-5*a2);
306     G4double e=390.;                       //     331     G4double e=390.;                       // Defolt values for deutrons
307     G4double r=0.27;                              332     G4double r=0.27;
308     G4double h=2.E-7;                             333     G4double h=2.E-7;
309     G4double t=0.3;                               334     G4double t=0.3;
310     if(tZ>1 || tN>1)                              335     if(tZ>1 || tN>1)
311     {                                             336     {
312       e=380.+18.*a2/(1.+a2/60.)/(1.+2.E-19*a8)    337       e=380.+18.*a2/(1.+a2/60.)/(1.+2.E-19*a8);
313       r=0.15;                                     338       r=0.15;
314       h=1.E-8*a2/(1.+a2/17.)/(1.+3.E-20*a8);      339       h=1.E-8*a2/(1.+a2/17.)/(1.+3.E-20*a8);
315       t=(.2+.00056*a2)/(1.+a2*.0006);             340       t=(.2+.00056*a2)/(1.+a2*.0006);
316     }                                             341     }
317     sigma=(c+d*d)/(1.+t/ssp+r/p4)+(gg+e*G4Exp( << 342     sigma=(c+d*d)/(1.+t/ssp+r/p4)+(gg+e*std::exp(-6.*P))/(1.+h/p4/p4);
318 #ifdef pdebug                                     343 #ifdef pdebug
319     G4cout<<"G4QHyperonNucCS::CSForm: A="<<a<<    344     G4cout<<"G4QHyperonNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<gg
320           <<",d="<<d<<",r="<<r<<",e="<<e<<",h=    345           <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
321 #endif                                            346 #endif
322   }                                               347   }
323   else                                            348   else
324   {                                               349   {
325     G4cerr<<"-Warning-G4QHyperonNuclearCroSect    350     G4cerr<<"-Warning-G4QHyperonNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
326     sigma=0.;                                     351     sigma=0.;
327   }                                               352   }
328   if(sigma<0.) return 0.;                         353   if(sigma<0.) return 0.;
329   return sigma;                                   354   return sigma;  
330 }                                                 355 }
331                                                   356 
332 G4double G4ChipsHyperonInelasticXS::EquLinearF    357 G4double G4ChipsHyperonInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
333 {                                                 358 {
334   if(DX<=0. || N<2)                               359   if(DX<=0. || N<2)
335     {                                             360     {
336       G4cerr<<"***G4ChipsHyperonInelasticXS::E    361       G4cerr<<"***G4ChipsHyperonInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
337       return Y[0];                                362       return Y[0];
338     }                                             363     }
339                                                   364   
340   G4int    N2=N-2;                                365   G4int    N2=N-2;
341   G4double d=(X-X0)/DX;                           366   G4double d=(X-X0)/DX;
342   G4int         jj=static_cast<int>(d);        << 367   G4int         j=static_cast<int>(d);
343   if     (jj<0)  jj=0;                         << 368   if     (j<0)  j=0;
344   else if(jj>N2) jj=N2;                        << 369   else if(j>N2) j=N2;
345   d-=jj; // excess                             << 370   d-=j; // excess
346   G4double yi=Y[jj];                           << 371   G4double yi=Y[j];
347   G4double sigma=yi+(Y[jj+1]-yi)*d;            << 372   G4double sigma=yi+(Y[j+1]-yi)*d;
348                                                   373   
349   return sigma;                                   374   return sigma;
350 }                                                 375 }
351                                                   376