Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4ChipsPionMinusInelasticXS.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/G4ChipsPionMinusInelasticXS.cc (Version 11.3.0) and /processes/hadronic/cross_sections/src/G4ChipsPionMinusInelasticXS.cc (Version 11.1.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 //                                                 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 // G4 Physics class: G4ChipsPionMinusInelastic     30 // G4 Physics class: G4ChipsPionMinusInelasticXS for gamma+A cross sections
 31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20     31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03
 32 // The last update: M.V. Kossov, CERN/ITEP (Mo     32 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04
 33 //                                                 33 //
 34 // -------------------------------------------     34 // -------------------------------------------------------------------------------------
 35 // Short description: Cross-sections extracted     35 // Short description: Cross-sections extracted (by W.Pokorski) from the CHIPS package for 
 36 // pion interactions. Original author: M. Koss     36 // pion interactions. Original author: M. Kossov
 37 // -------------------------------------------     37 // -------------------------------------------------------------------------------------
 38 //                                                 38 //
 39                                                    39 
 40 #include "G4ChipsPionMinusInelasticXS.hh"          40 #include "G4ChipsPionMinusInelasticXS.hh"
 41 #include "G4ChipsPionPlusInelasticXS.hh"           41 #include "G4ChipsPionPlusInelasticXS.hh"
 42 #include "G4SystemOfUnits.hh"                      42 #include "G4SystemOfUnits.hh"
 43 #include "G4DynamicParticle.hh"                    43 #include "G4DynamicParticle.hh"
 44 #include "G4ParticleDefinition.hh"                 44 #include "G4ParticleDefinition.hh"
 45 #include "G4PionMinus.hh"                          45 #include "G4PionMinus.hh"
 46                                                    46 
 47 #include "G4Log.hh"                                47 #include "G4Log.hh"
 48 #include "G4Exp.hh"                                48 #include "G4Exp.hh"
 49                                                    49 
 50 // factory                                         50 // factory
 51 #include "G4CrossSectionFactory.hh"                51 #include "G4CrossSectionFactory.hh"
 52 //                                                 52 //
 53 G4_DECLARE_XS_FACTORY(G4ChipsPionMinusInelasti     53 G4_DECLARE_XS_FACTORY(G4ChipsPionMinusInelasticXS);
 54                                                    54 
 55 G4ChipsPionMinusInelasticXS::G4ChipsPionMinusI     55 G4ChipsPionMinusInelasticXS::G4ChipsPionMinusInelasticXS():G4VCrossSectionDataSet("ChipsPionMinusInelasticXS")
 56 {                                                  56 {
 57   // Initialization of the                         57   // Initialization of the
 58   lastLEN=0; // Pointer to lastArray of LowEn      58   lastLEN=0; // Pointer to lastArray of LowEn CS
 59   lastHEN=0; // Pointer to lastArray of HighEn     59   lastHEN=0; // Pointer to lastArray of HighEn CS
 60   lastN=0;   // The last N of calculated nucle     60   lastN=0;   // The last N of calculated nucleus
 61   lastZ=0;   // The last Z of calculated nucle     61   lastZ=0;   // The last Z of calculated nucleus
 62   lastP=0.;  // Last used cross section Moment     62   lastP=0.;  // Last used cross section Momentum
 63   lastTH=0.; // Last threshold momentum            63   lastTH=0.; // Last threshold momentum
 64   lastCS=0.; // Last value of the Cross Sectio     64   lastCS=0.; // Last value of the Cross Section
 65   lastI=0;   // The last position in the DAMDB     65   lastI=0;   // The last position in the DAMDB
 66   LEN = new std::vector<G4double*>;                66   LEN = new std::vector<G4double*>;
 67   HEN = new std::vector<G4double*>;                67   HEN = new std::vector<G4double*>;
 68 }                                                  68 }
 69                                                    69 
 70 G4ChipsPionMinusInelasticXS::~G4ChipsPionMinus     70 G4ChipsPionMinusInelasticXS::~G4ChipsPionMinusInelasticXS()
 71 {                                                  71 {
 72   std::size_t lens=LEN->size();                    72   std::size_t lens=LEN->size();
 73   for(std::size_t i=0; i<lens; ++i) delete[] (     73   for(std::size_t i=0; i<lens; ++i) delete[] (*LEN)[i];
 74   delete LEN;                                      74   delete LEN;
 75   std::size_t hens=HEN->size();                    75   std::size_t hens=HEN->size();
 76   for(std::size_t i=0; i<hens; ++i) delete[] (     76   for(std::size_t i=0; i<hens; ++i) delete[] (*HEN)[i];
 77   delete HEN;                                      77   delete HEN;
 78 }                                                  78 }
 79                                                    79 
 80 void                                               80 void
 81 G4ChipsPionMinusInelasticXS::CrossSectionDescr     81 G4ChipsPionMinusInelasticXS::CrossSectionDescription(std::ostream& outFile) const
 82 {                                                  82 {
 83     outFile << "G4ChipsPionMinusInelasticXS pr     83     outFile << "G4ChipsPionMinusInelasticXS provides the inelastic cross\n"
 84             << "section for pion- nucleus scat     84             << "section for pion- nucleus scattering as a function of incident\n"
 85             << "momentum. The cross section is     85             << "momentum. The cross section is calculated using M. Kossov's\n"
 86             << "CHIPS parameterization of cros     86             << "CHIPS parameterization of cross section data.\n";
 87 }                                                  87 }
 88                                                    88 
 89 G4bool G4ChipsPionMinusInelasticXS::IsIsoAppli     89 G4bool G4ChipsPionMinusInelasticXS::IsIsoApplicable(const G4DynamicParticle*, G4int, G4int,    
 90          const G4Element*,                         90          const G4Element*,
 91          const G4Material*)                        91          const G4Material*)
 92 {                                                  92 {
 93   return true;                                     93   return true;
 94 }                                                  94 }
 95                                                    95 
 96 // The main member function giving the collisi     96 // The main member function giving the collision cross section (P is in IU, CS is in mb)
 97 // Make pMom in independent units ! (Now it is     97 // Make pMom in independent units ! (Now it is MeV)
 98 G4double G4ChipsPionMinusInelasticXS::GetIsoCr     98 G4double G4ChipsPionMinusInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,  
 99               const G4Isotope*,                    99               const G4Isotope*,
100               const G4Element*,                   100               const G4Element*,
101               const G4Material*)                  101               const G4Material*)
102 {                                                 102 {
103   G4double pMom=Pt->GetTotalMomentum();           103   G4double pMom=Pt->GetTotalMomentum();
104   G4int tgN = A - tgZ;                            104   G4int tgN = A - tgZ;
105                                                   105   
106   return GetChipsCrossSection(pMom, tgZ, tgN,     106   return GetChipsCrossSection(pMom, tgZ, tgN, -211);
107 }                                                 107 }
108                                                   108 
109 G4double G4ChipsPionMinusInelasticXS::GetChips    109 G4double G4ChipsPionMinusInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
110 {                                                 110 {
111                                                   111 
112   G4bool in=false;                     // By d    112   G4bool in=false;                     // By default the isotope must be found in the AMDB
113   if(tgN!=lastN || tgZ!=lastZ)         // The     113   if(tgN!=lastN || tgZ!=lastZ)         // The nucleus was not the last used isotope
114   {                                               114   {
115     in = false;                        // By d    115     in = false;                        // By default the isotope haven't be found in AMDB  
116     lastP   = 0.;                      // New     116     lastP   = 0.;                      // New momentum history (nothing to compare with)
117     lastN   = tgN;                     // The     117     lastN   = tgN;                     // The last N of the calculated nucleus
118     lastZ   = tgZ;                     // The     118     lastZ   = tgZ;                     // The last Z of the calculated nucleus
119     lastI   = (G4int)colN.size();      // Size    119     lastI   = (G4int)colN.size();      // Size of the Associative Memory DB in the heap
120     j  = 0;                            // A#0f    120     j  = 0;                            // A#0f records found in DB for this projectile
121     if(lastI) for(G4int i=0; i<lastI; ++i) //     121     if(lastI) for(G4int i=0; i<lastI; ++i) // AMDB exists, try to find the (Z,N) isotope
122     {                                             122     {
123       if(colN[i]==tgN && colZ[i]==tgZ) // Try     123       if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
124       {                                           124       {
125         lastI=i;                       // Reme    125         lastI=i;                       // Remember the index for future fast/last use
126         lastTH =colTH[i];              // The     126         lastTH =colTH[i];              // The last THreshold (A-dependent)
127         if(pMom<=lastTH)                          127         if(pMom<=lastTH)
128         {                                         128         {
129           return 0.;                   // Ener    129           return 0.;                   // Energy is below the Threshold value
130         }                                         130         }
131         lastP  =colP [i];              // Last    131         lastP  =colP [i];              // Last Momentum  (A-dependent)
132         lastCS =colCS[i];              // Last    132         lastCS =colCS[i];              // Last CrossSect (A-dependent)
133         in = true;                     // This    133         in = true;                     // This is the case when the isotop is found in DB
134         // Momentum pMom is in IU ! @@ Units      134         // Momentum pMom is in IU ! @@ Units
135         lastCS=CalculateCrossSection(-1,j,-211    135         lastCS=CalculateCrossSection(-1,j,-211,lastZ,lastN,pMom); // read & update
136         if(lastCS<=0. && pMom>lastTH)  // Corr    136         if(lastCS<=0. && pMom>lastTH)  // Correct the threshold (@@ No intermediate Zeros)
137         {                                         137         {
138           lastCS=0.;                              138           lastCS=0.;
139           lastTH=pMom;                            139           lastTH=pMom;
140         }                                         140         }
141         break;                         // Go o    141         break;                         // Go out of the LOOP
142       }                                           142       }
143       j++;                             // Incr    143       j++;                             // Increment a#0f records found in DB
144     }                                             144     }
145     if(!in)                            // This    145     if(!in)                            // This isotope has not been calculated previously
146     {                                             146     {
147       //!!The slave functions must provide cro    147       //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
148       lastCS=CalculateCrossSection(0,j,-211,la    148       lastCS=CalculateCrossSection(0,j,-211,lastZ,lastN,pMom); //calculate & create
149       //if(lastCS>0.)                   // It     149       //if(lastCS>0.)                   // It means that the AMBD was initialized
150       //{                                         150       //{
151                                                   151 
152       lastTH = 0; //ThresholdEnergy(tgZ, tgN);    152       lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
153         colN.push_back(tgN);                      153         colN.push_back(tgN);
154         colZ.push_back(tgZ);                      154         colZ.push_back(tgZ);
155         colP.push_back(pMom);                     155         colP.push_back(pMom);
156         colTH.push_back(lastTH);                  156         colTH.push_back(lastTH);
157         colCS.push_back(lastCS);                  157         colCS.push_back(lastCS);
158       //} // M.K. Presence of H1 with high thr    158       //} // M.K. Presence of H1 with high threshold breaks the syncronization
159       return lastCS*millibarn;                    159       return lastCS*millibarn;
160     } // End of creation of the new set of par    160     } // End of creation of the new set of parameters
161     else                                          161     else
162     {                                             162     {
163       colP[lastI]=pMom;                           163       colP[lastI]=pMom;
164       colCS[lastI]=lastCS;                        164       colCS[lastI]=lastCS;
165     }                                             165     }
166   } // End of parameters udate                    166   } // End of parameters udate
167   else if(pMom<=lastTH)                           167   else if(pMom<=lastTH)
168   {                                               168   {
169     return 0.;                         // Mome    169     return 0.;                         // Momentum is below the Threshold Value -> CS=0
170   }                                               170   }
171   else                                 // It i    171   else                                 // It is the last used -> use the current tables
172   {                                               172   {
173     lastCS=CalculateCrossSection(1,j,-211,last    173     lastCS=CalculateCrossSection(1,j,-211,lastZ,lastN,pMom); // Only read and UpdateDB
174     lastP=pMom;                                   174     lastP=pMom;
175   }                                               175   }
176   return lastCS*millibarn;                        176   return lastCS*millibarn;
177 }                                                 177 }
178                                                   178 
179 // The main member function giving the gamma-A    179 // The main member function giving the gamma-A cross section (E in GeV, CS in mb)
180 G4double G4ChipsPionMinusInelasticXS::Calculat    180 G4double G4ChipsPionMinusInelasticXS::CalculateCrossSection(G4int F, G4int I,
181                                         G4int,    181                                         G4int, G4int targZ, G4int targN, G4double Momentum)
182 {                                                 182 {
183   static const G4double THmin=27.;     // defa    183   static const G4double THmin=27.;     // default minimum Momentum (MeV/c) Threshold
184   static const G4double THmiG=THmin*.001; // m    184   static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
185   static const G4double dP=10.;        // step    185   static const G4double dP=10.;        // step for the LEN (Low ENergy) table MeV/c
186   static const G4double dPG=dP*.001;   // step    186   static const G4double dPG=dP*.001;   // step for the LEN (Low ENergy) table GeV/c
187   static const G4int    nL=105;        // A#of    187   static const G4int    nL=105;        // A#of LEN points in E (step 10 MeV/c)
188   static const G4double Pmin=THmin+(nL-1)*dP;     188   static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
189   static const G4double Pmax=227000.;  // maxP    189   static const G4double Pmax=227000.;  // maxP for the HEN (High ENergy) part 227 GeV
190   static const G4int    nH=224;        // A#of    190   static const G4int    nH=224;        // A#of HEN points in lnE
191   static const G4double milP=G4Log(Pmin);// Lo    191   static const G4double milP=G4Log(Pmin);// Low logarithm energy for the HEN part
192   static const G4double malP=G4Log(Pmax);// Hi    192   static const G4double malP=G4Log(Pmax);// High logarithm energy (each 2.75 percent)
193   static const G4double dlP=(malP-milP)/(nH-1)    193   static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
194   static const G4double milPG=G4Log(.001*Pmin)    194   static const G4double milPG=G4Log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
                                                   >> 195   G4double sigma=0.;
                                                   >> 196   if(F&&I) sigma=0.;                   // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
                                                   >> 197   //G4double A=targN+targZ;              // A of the target
195   if(F<=0)                             // This    198   if(F<=0)                             // This isotope was not the last used isotop
196   {                                               199   {
197     if(F<0)                            // This    200     if(F<0)                            // This isotope was found in DAMDB =-----=> RETRIEVE
198     {                                             201     {
199       G4int sync=(G4int)LEN->size();              202       G4int sync=(G4int)LEN->size();
200       if(sync<=I) G4cerr<<"*!*G4ChipsPiMinusNu    203       if(sync<=I) G4cerr<<"*!*G4ChipsPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
201       lastLEN=(*LEN)[I];               // Poin    204       lastLEN=(*LEN)[I];               // Pointer to prepared LowEnergy cross sections
202       lastHEN=(*HEN)[I];               // Poin    205       lastHEN=(*HEN)[I];               // Pointer to prepared High Energy cross sections
203     }                                             206     }
204     else                               // This    207     else                               // This isotope wasn't calculated before => CREATE
205     {                                             208     {
206       lastLEN = new G4double[nL];      // Allo    209       lastLEN = new G4double[nL];      // Allocate memory for the new LEN cross sections
207       lastHEN = new G4double[nH];      // Allo    210       lastHEN = new G4double[nH];      // Allocate memory for the new HEN cross sections
208       // --- Instead of making a separate func    211       // --- Instead of making a separate function ---
209       G4double P=THmiG;                // Tabl    212       G4double P=THmiG;                // Table threshold in GeV/c
210      for(G4int k=0; k<nL; k++)                    213      for(G4int k=0; k<nL; k++)
211       {                                           214       {
212         lastLEN[k] = CrossSectionLin(targZ, ta    215         lastLEN[k] = CrossSectionLin(targZ, targN, P);
213         P+=dPG;                                   216         P+=dPG;
214       }                                           217       }
215       G4double lP=milPG;                          218       G4double lP=milPG;
216       for(G4int n=0; n<nH; n++)                   219       for(G4int n=0; n<nH; n++)
217       {                                           220       {
218         lastHEN[n] = CrossSectionLog(targZ, ta    221         lastHEN[n] = CrossSectionLog(targZ, targN, lP);
219         lP+=dlP;                                  222         lP+=dlP;
220       }                                           223       }
221       // --- End of possible separate function    224       // --- End of possible separate function
222       // *** The synchronization check ***        225       // *** The synchronization check ***
223       G4int sync=(G4int)LEN->size();              226       G4int sync=(G4int)LEN->size();
224       if(sync!=I)                                 227       if(sync!=I)
225       {                                           228       {
226         G4cerr<<"***G4ChipsPiMinusNuclCS::Calc    229         G4cerr<<"***G4ChipsPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
227               <<", N="<<targN<<", F="<<F<<G4en    230               <<", N="<<targN<<", F="<<F<<G4endl;
228         //G4Exception("G4PiMinusNuclearCS::Cal    231         //G4Exception("G4PiMinusNuclearCS::CalculateCS:","39",FatalException,"DBoverflow");
229       }                                           232       }
230       LEN->push_back(lastLEN);         // reme    233       LEN->push_back(lastLEN);         // remember the Low Energy Table
231       HEN->push_back(lastHEN);         // reme    234       HEN->push_back(lastHEN);         // remember the High Energy Table
232     } // End of creation of the new set of par    235     } // End of creation of the new set of parameters
233   } // End of parameters udate                    236   } // End of parameters udate
234   // =---------------------= NOW the Magic For    237   // =---------------------= NOW the Magic Formula =---------------------------=
235   G4double sigma;                              << 
236   if (Momentum<lastTH) return 0.;      // It m    238   if (Momentum<lastTH) return 0.;      // It must be already checked in the interface class
237   else if (Momentum<Pmin)              // High    239   else if (Momentum<Pmin)              // High Energy region
238   {                                               240   {
239     sigma=EquLinearFit(Momentum,nL,THmin,dP,la    241     sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
240   }                                               242   }
241   else if (Momentum<Pmax)              // High    243   else if (Momentum<Pmax)              // High Energy region
242   {                                               244   {
243     G4double lP=G4Log(Momentum);                  245     G4double lP=G4Log(Momentum);
244     sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN)    246     sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
245   }                                               247   }
246   else                                 // UHE     248   else                                 // UHE region (calculation, not frequent)
247   {                                               249   {
248     G4double P=0.001*Momentum;         // Appr    250     G4double P=0.001*Momentum;         // Approximation formula is for P in GeV/c
249     sigma=CrossSectionFormula(targZ, targN, P,    251     sigma=CrossSectionFormula(targZ, targN, P, G4Log(P));
250   }                                               252   }
251   if(sigma<0.) return 0.;                         253   if(sigma<0.) return 0.;
252   return sigma;                                   254   return sigma;
253 }                                                 255 }
254                                                   256 
255 // Calculation formula for piMinus-nuclear ine    257 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) (P in GeV/c)
256 G4double G4ChipsPionMinusInelasticXS::CrossSec    258 G4double G4ChipsPionMinusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
257 {                                                 259 {
258   G4double lP=G4Log(P);                           260   G4double lP=G4Log(P);
259   return CrossSectionFormula(tZ, tN, P, lP);      261   return CrossSectionFormula(tZ, tN, P, lP);
260 }                                                 262 }
261                                                   263 
262 // Calculation formula for piMinus-nuclear ine    264 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
263 G4double G4ChipsPionMinusInelasticXS::CrossSec    265 G4double G4ChipsPionMinusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
264 {                                                 266 {
265   G4double P=G4Exp(lP);                           267   G4double P=G4Exp(lP);
266   return CrossSectionFormula(tZ, tN, P, lP);      268   return CrossSectionFormula(tZ, tN, P, lP);
267 }                                                 269 }
268 // Calculation formula for piMinus-nuclear ine    270 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
269 G4double G4ChipsPionMinusInelasticXS::CrossSec    271 G4double G4ChipsPionMinusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
270                                                   272                                                               G4double P, G4double lP)
271 {                                                 273 {
272   G4double sigma=0.;                              274   G4double sigma=0.;
273   if(tZ==1 && !tN)                        // P    275   if(tZ==1 && !tN)                        // PiMin-Proton interaction from G4QuasiElRatios
274   {                                               276   {
275     G4double lr=lP+1.27;                    //    277     G4double lr=lP+1.27;                    // From G4QuasiFreeRatios.cc Uzhi
276     G4double LE=1.53/(lr*lr+.0676);         //    278     G4double LE=1.53/(lr*lr+.0676);         // From G4QuasiFreeRatios.cc Uzhi
277     G4double ld=lP-3.5;                           279     G4double ld=lP-3.5;
278     G4double ld2=ld*ld;                           280     G4double ld2=ld*ld;
279     G4double p2=P*P;                              281     G4double p2=P*P;
280     G4double p4=p2*p2;                            282     G4double p4=p2*p2;
281     G4double sp=std::sqrt(P);                     283     G4double sp=std::sqrt(P);
282     G4double lm=lP+.36;                           284     G4double lm=lP+.36;
283     G4double md=lm*lm+.04;                        285     G4double md=lm*lm+.04;
284     G4double lh=lP-.017;                          286     G4double lh=lP-.017;
285     G4double hd=lh*lh+.0025;                      287     G4double hd=lh*lh+.0025;
286     G4double El=(.0557*ld2+2.4+7./sp)/(1.+.7/p    288     G4double El=(.0557*ld2+2.4+7./sp)/(1.+.7/p4);
287     G4double To=(.3*ld2+22.3+12./sp)/(1.+.4/p4    289     G4double To=(.3*ld2+22.3+12./sp)/(1.+.4/p4);
288     sigma=(To-El)+.4/md+.01/hd;                   290     sigma=(To-El)+.4/md+.01/hd;
289     sigma+=LE*2;                            //    291     sigma+=LE*2;                            //  Uzhi       
290   }                                               292   }
291   else if(tZ==1 && tN==1)                   //    293   else if(tZ==1 && tN==1)                   // pimp_tot
292   {                                               294   {
293     G4double p2=P*P;                              295     G4double p2=P*P;
294     G4double d=lP-2.7;                            296     G4double d=lP-2.7;
295     G4double f=lP+1.25;                           297     G4double f=lP+1.25;
296     G4double gg=lP-.017;                          298     G4double gg=lP-.017;
297     sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.    299     sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.3/p2/p2)+18./(f*f+.1089)+.02/(gg*gg+.0025);
298   }                                               300   }
299   else if(tZ<97 && tN<152)                // G    301   else if(tZ<97 && tN<152)                // General solution
300   {                                               302   {
301     G4double d=lP-4.2;                            303     G4double d=lP-4.2;
302     G4double p2=P*P;                              304     G4double p2=P*P;
303     G4double p4=p2*p2;                            305     G4double p4=p2*p2;
304     G4double a=tN+tZ;                     // A    306     G4double a=tN+tZ;                     // A of the target
305     G4double al=G4Log(a);                         307     G4double al=G4Log(a);
306     G4double sa=std::sqrt(a);                     308     G4double sa=std::sqrt(a);
307     G4double ssa=std::sqrt(sa);                   309     G4double ssa=std::sqrt(sa);
308     G4double a2=a*a;                              310     G4double a2=a*a;
309     G4double c=41.*G4Exp(al*.68)*(1.+44./a2)/(    311     G4double c=41.*G4Exp(al*.68)*(1.+44./a2)/(1.+8./a)/(1.+200./a2/a2);
310     G4double f=120*sa/(1.+24./a/ssa);             312     G4double f=120*sa/(1.+24./a/ssa);
311     G4double gg=-1.32-al*.043;                    313     G4double gg=-1.32-al*.043;
312     G4double u=lP-gg;                             314     G4double u=lP-gg;
313     G4double h=al*(.388-.046*al);                 315     G4double h=al*(.388-.046*al);
314     sigma=(c+d*d)/(1.+.17/p4)+f/(u*u+h*h);        316     sigma=(c+d*d)/(1.+.17/p4)+f/(u*u+h*h);
315   }                                               317   }
316   else                                            318   else
317   {                                               319   {
318     G4cerr<<"-Warning-G4ChipsPiMinusNuclearCro    320     G4cerr<<"-Warning-G4ChipsPiMinusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
319     sigma=0.;                                     321     sigma=0.;
320   }                                               322   }
321   if(sigma<0.) return 0.;                         323   if(sigma<0.) return 0.;
322   return sigma;                                   324   return sigma;  
323 }                                                 325 }
324                                                   326 
325 G4double G4ChipsPionMinusInelasticXS::EquLinea    327 G4double G4ChipsPionMinusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
326 {                                                 328 {
327   if(DX<=0. || N<2)                               329   if(DX<=0. || N<2)
328     {                                             330     {
329       G4cerr<<"***G4ChipsPionMinusInelasticXS:    331       G4cerr<<"***G4ChipsPionMinusInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
330       return Y[0];                                332       return Y[0];
331     }                                             333     }
332                                                   334   
333   G4int    N2=N-2;                                335   G4int    N2=N-2;
334   G4double d=(X-X0)/DX;                           336   G4double d=(X-X0)/DX;
335   G4int         jj=static_cast<int>(d);           337   G4int         jj=static_cast<int>(d);
336   if     (jj<0)  jj=0;                            338   if     (jj<0)  jj=0;
337   else if(jj>N2) jj=N2;                           339   else if(jj>N2) jj=N2;
338   d-=jj; // excess                                340   d-=jj; // excess
339   G4double yi=Y[jj];                              341   G4double yi=Y[jj];
340   G4double sigma=yi+(Y[jj+1]-yi)*d;               342   G4double sigma=yi+(Y[jj+1]-yi)*d;
341                                                   343   
342   return sigma;                                   344   return sigma;
343 }                                                 345 }
344                                                   346