Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/quasi_elastic/src/G4QuasiElRatios.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/quasi_elastic/src/G4QuasiElRatios.cc (Version 11.3.0) and /processes/hadronic/models/quasi_elastic/src/G4QuasiElRatios.cc (Version 9.6.p2)


  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 // $Id$
 27 //                                                 28 //
 28 //                                                 29 //
 29 // G4 Physics class: G4QuasiElRatios for N+A e     30 // G4 Physics class: G4QuasiElRatios for N+A elastic cross sections
 30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10     31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
 31 // The last update: M.V. Kossov, CERN/ITEP (Mo     32 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06
 32 //                                                 33 // 
 33 // -------------------------------------------     34 // ----------------------------------------------------------------------
 34 // This class has been extracted from the CHIP     35 // This class has been extracted from the CHIPS model. 
 35 // All the dependencies on CHIPS classes have      36 // All the dependencies on CHIPS classes have been removed.
 36 // Short description: Provides percentage of q     37 // Short description: Provides percentage of quasi-free and quasi-elastic
 37 // reactions in the inelastic reactions.           38 // reactions in the inelastic reactions.
 38 // -------------------------------------------     39 // ----------------------------------------------------------------------
 39                                                    40 
 40                                                    41 
 41 #include "G4QuasiElRatios.hh"                      42 #include "G4QuasiElRatios.hh"
 42 #include "G4PhysicalConstants.hh"                  43 #include "G4PhysicalConstants.hh"
 43 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 44 #include "G4Proton.hh"                             45 #include "G4Proton.hh"
 45 #include "G4Neutron.hh"                            46 #include "G4Neutron.hh"
 46 #include "G4Deuteron.hh"                           47 #include "G4Deuteron.hh"
 47 #include "G4Triton.hh"                             48 #include "G4Triton.hh"
 48 #include "G4He3.hh"                                49 #include "G4He3.hh"
 49 #include "G4Alpha.hh"                              50 #include "G4Alpha.hh"
 50 #include "G4ThreeVector.hh"                        51 #include "G4ThreeVector.hh"
 51 #include "G4CrossSectionDataSetRegistry.hh"        52 #include "G4CrossSectionDataSetRegistry.hh"
 52 #include "G4Pow.hh"                            << 
 53 #include "G4Log.hh"                            << 
 54 #include "G4Exp.hh"                            << 
 55                                                << 
 56 namespace  {                                   << 
 57     const G4int    nps=150;        // Number o << 
 58     const G4int    mps=nps+1;      // Number o << 
 59     const G4double sma=150.;       // The firs << 
 60     const G4double ds=sma/nps;     // Step of  << 
 61     const G4int    nls=100;        // Number o << 
 62     const G4int    mls=nls+1;      // Number o << 
 63     const G4double lsi=5.;         // The min  << 
 64     const G4double lsa=9.;         // The max  << 
 65     const G4double mi=G4Exp(lsi);// The min s  << 
 66     const G4double min_s=G4Exp(lsa);// The max << 
 67     const G4double dls=(lsa-lsi)/nls;// Step o << 
 68     const G4double edls=G4Exp(dls);// Multipli << 
 69     const G4double toler=.01;      // The tola << 
 70     const G4double C=1.246;                    << 
 71     const G4double lmi=3.5;       // min of (l << 
 72     const G4double pbe=.0557;     // elastic ( << 
 73     const G4double pbt=.3;        // total (ln << 
 74     const G4double pmi=.1;        // Below tha << 
 75     const G4double pma=1000.;     // Above tha << 
 76     const G4int    nlp=300;         // Number  << 
 77     const G4int    mlp=nlp+1;       // Number  << 
 78     const G4double lpi=-5.;         // The min << 
 79     const G4double lpa=10.;         // The max << 
 80     const G4double mip=G4Exp(lpi);// The min p << 
 81     const G4double map=G4Exp(lpa);// The max p << 
 82     const G4double dlp=(lpa-lpi)/nlp;// Step o << 
 83     const G4double edlp=G4Exp(dlp);// Multipli << 
 84 }                                              << 
 85                                                    53 
 86                                                    54 
                                                   >>  55 // initialisation of statics
                                                   >>  56 std::vector<G4double*> G4QuasiElRatios::vT; // Vector of pointers to LinTable in C++ heap
                                                   >>  57 std::vector<G4double*> G4QuasiElRatios::vL; // Vector of pointers to LogTable in C++ heap
                                                   >>  58 std::vector<std::pair<G4double,G4double>*> G4QuasiElRatios::vX; // ETPointers to LogTable
                                                   >>  59 
 87 G4QuasiElRatios::G4QuasiElRatios()                 60 G4QuasiElRatios::G4QuasiElRatios()
 88 {                                                  61 {
 89     vT = new std::vector<G4double*>;           << 
 90     vL = new std::vector<G4double*>;           << 
 91     vX = new std::vector<std::pair<G4double,G4 << 
 92                                                << 
 93     lastSRatio=0.;             // The last sig << 
 94     lastRRatio=0.;             // The last rat << 
 95     lastARatio=0;             // theLast of ca << 
 96     lastHRatio=0.;            // theLast of ma << 
 97     lastNRatio=0;             // theLast of to << 
 98     lastMRatio=0.;            // theLast of re << 
 99     lastKRatio=0;             // theLast of to << 
100     lastTRatio=0;             // theLast of po << 
101     lastLRatio=0;             // theLast of po << 
102     lastPtot=0.;              // The last mome << 
103     lastHtot=0;               // The last proj << 
104     lastFtot=true;            // The last nucl << 
105     lastItot=0;              // The Last index << 
106     lastMtot=0.;             // The Last rel m << 
107     lastKtot=0;             // The Last topBin << 
108     lastXtot = nullptr;     // The Last ETPoin << 
109                                                    62     
110     PCSmanager=(G4ChipsProtonElasticXS*)G4Cros     63     PCSmanager=(G4ChipsProtonElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name());
111                                                    64     
112     NCSmanager=(G4ChipsNeutronElasticXS*)G4Cro     65     NCSmanager=(G4ChipsNeutronElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name());
113 }                                                  66 }
114                                                    67 
115 G4QuasiElRatios::~G4QuasiElRatios()                68 G4QuasiElRatios::~G4QuasiElRatios()
116 {                                                  69 {
117     std::vector<G4double*>::iterator pos;          70     std::vector<G4double*>::iterator pos;
118     for(pos=vT->begin(); pos<vT->end(); pos++) <<  71     for(pos=vT.begin(); pos<vT.end(); pos++)
119     { delete [] *pos; }                            72     { delete [] *pos; }
120     vT->clear();                               <<  73     vT.clear();
121     delete vT;                                 <<  74     for(pos=vL.begin(); pos<vL.end(); pos++)
122                                                << 
123     for(pos=vL->begin(); pos<vL->end(); pos++) << 
124     { delete [] *pos; }                            75     { delete [] *pos; }
125     vL->clear();                               <<  76     vL.clear();
126     delete vL;                                 << 
127                                                    77     
128     std::vector<std::pair<G4double,G4double>*>     78     std::vector<std::pair<G4double,G4double>*>::iterator pos2;
129     for(pos2=vX->begin(); pos2<vX->end(); pos2 <<  79     for(pos2=vX.begin(); pos2<vX.end(); pos2++)
130     { delete [] *pos2; }                           80     { delete [] *pos2; }
131     vX->clear();                               <<  81     vX.clear();
132     delete vX;                                 <<  82 }
                                                   >>  83 
                                                   >>  84 // Returns Pointer to the G4VQCrossSection class
                                                   >>  85 G4QuasiElRatios* G4QuasiElRatios::GetPointer()
                                                   >>  86 {
                                                   >>  87     static G4QuasiElRatios theRatios;   // *** Static body of the QEl Cross Section ***
                                                   >>  88     return &theRatios;
133 }                                                  89 }
134                                                    90 
135 // Calculation of pair(QuasiFree/Inelastic,Qua     91 // Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree)
136 std::pair<G4double,G4double> G4QuasiElRatios::     92 std::pair<G4double,G4double> G4QuasiElRatios::GetRatios(G4double pIU, G4int pPDG,
137                                                    93                                                          G4int tgZ,    G4int tgN)
138 {                                                  94 {
139     G4double R=0.;                                 95     G4double R=0.;
140     G4double QF2In=1.;                             96     G4double QF2In=1.;                        // Prototype of QuasiFree/Inel ratio for hN_tot
141     G4int tgA=tgZ+tgN;                             97     G4int tgA=tgZ+tgN;
142     if(tgA<2) return std::make_pair(QF2In,R);      98     if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
143     std::pair<G4double,G4double> ElTot=GetElTo     99     std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU)
144     //if( ( (pPDG>999 && pIU<227.) || pIU<27.)    100     //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
145     if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.    101     if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;                // To accelerate @lowE
146     else if(ElTot.second>0.)                      102     else if(ElTot.second>0.)
147     {                                             103     {
148         R=ElTot.first/ElTot.second;               104         R=ElTot.first/ElTot.second;             // El/Total ratio (does not depend on units
149         QF2In=GetQF2IN_Ratio(ElTot.second/mill    105         QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN);   // QuasiFree/Inelastic ratio
150     }                                             106     }
151     return std::make_pair(QF2In,R);               107     return std::make_pair(QF2In,R);
152 }                                                 108 }
153                                                   109 
154 // Calculatio QasiFree/Inelastic Ratio as a fu    110 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A
155 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4dou    111 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4double m_s, G4int A)
156 {                                                 112 {
                                                   >> 113     static const G4int    nps=150;        // Number of steps in the R(s) LinTable
                                                   >> 114     static const G4int    mps=nps+1;      // Number of elements in the R(s) LinTable
                                                   >> 115     static const G4double sma=150.;       // The first LinTabEl(s=0)=1., s>sma -> logTab
                                                   >> 116     static const G4double ds=sma/nps;     // Step of the linear Table
                                                   >> 117     static const G4int    nls=100;        // Number of steps in the R(lns) logTable
                                                   >> 118     static const G4int    mls=nls+1;      // Number of elements in the R(lns) logTable
                                                   >> 119     static const G4double lsi=5.;         // The min ln(s) logTabEl(s=148.4 < sma=150.)
                                                   >> 120     static const G4double lsa=9.;         // The max ln(s) logTabEl(s=148.4 - 8103. mb)
                                                   >> 121     static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 148.4 mb)
                                                   >> 122     static const G4double min_s=std::exp(lsa);// The max s of logTabEl(~ 8103. mb)
                                                   >> 123     static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table
                                                   >> 124     static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
                                                   >> 125     static const G4double toler=.01;      // The tolarence mb defining the same cross-section
                                                   >> 126     static G4double lastS=0.;             // The last sigma value for which R was calculated
                                                   >> 127     static G4double lastR=0.;             // The last ratio R which was calculated
                                                   >> 128     // Local Associative Data Base:
                                                   >> 129     static std::vector<G4int>     vA;     // Vector of calculated A
                                                   >> 130     static std::vector<G4double>  vH;     // Vector of max s initialized in the LinTable
                                                   >> 131     static std::vector<G4int>     vN;     // Vector of topBin number initialized in LinTable
                                                   >> 132     static std::vector<G4double>  vM;     // Vector of rel max ln(s) initialized in LogTable
                                                   >> 133     static std::vector<G4int>     vK;     // Vector of topBin number initialized in LogTable
                                                   >> 134     // Last values of the Associative Data Base:
                                                   >> 135     static G4int     lastA=0;             // theLast of calculated A
                                                   >> 136     static G4double  lastH=0.;            // theLast of max s initialized in the LinTable
                                                   >> 137     static G4int     lastN=0;             // theLast of topBin number initialized in LinTable
                                                   >> 138     static G4double  lastM=0.;            // theLast of rel max ln(s) initialized in LogTable
                                                   >> 139     static G4int     lastK=0;             // theLast of topBin number initialized in LogTable
                                                   >> 140     static G4double* lastT=0;             // theLast of pointer to LinTable in the C++ heap
                                                   >> 141     static G4double* lastL=0;             // theLast of pointer to LogTable in the C++ heap
157     // LogTable is created only if necessary.     142     // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
158     if(m_s<toler || A<2) return 1.;               143     if(m_s<toler || A<2) return 1.;
159     if(m_s>min_s) return 0.;                      144     if(m_s>min_s) return 0.;
160                                                << 145     if(A>238)
161     //if(A>238)                                << 146     {
162     //{                                        << 147         G4cout<<"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
163     //    G4cout<<"-Warning-G4QuasiElRatio::Ge << 148         return 0.;
164     //    return 0.;                           << 149     }
165     //}                                        << 150     G4int nDB=vA.size();                  // A number of nuclei already initialized in AMDB
166     G4int nDB=(G4int)vARatio.size();  // A num << 151     if(nDB && lastA==A && m_s==lastS) return lastR;  // VI do not use tolerance
167     if(nDB && lastARatio==A && m_s==lastSRatio << 
168     G4bool found=false;                           152     G4bool found=false;
169     G4int i=-1;                                   153     G4int i=-1;
170     if(nDB) for (i=0; i<nDB; i++) if(A==vARati << 154     if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Search for this A in AMDB
171     {                                             155     {
172         found=true;                         //    156         found=true;                         // The A value is found
173         break;                                    157         break;
174     }                                             158     }
175     if(!nDB || !found)                    // C    159     if(!nDB || !found)                    // Create new line in the AMDB
176     {                                             160     {
177         lastARatio = A;                        << 161         lastA = A;
178         lastTRatio = new G4double[mps];        << 162         lastT = new G4double[mps];          // Create the linear Table
179         lastNRatio = static_cast<int>(m_s/ds)+ << 163         lastN = static_cast<int>(m_s/ds)+1;   // MaxBin to be initialized
180         if(lastNRatio>nps)                     << 164         if(lastN>nps)
181         {                                         165         {
182             lastNRatio=nps;                    << 166             lastN=nps;
183             lastHRatio=sma;                    << 167             lastH=sma;
184         }                                         168         }
185         else lastHRatio = lastNRatio*ds;       << 169         else lastH = lastN*ds;              // Calculate max initialized s for LinTab
186         G4double sv=0;                            170         G4double sv=0;
187         lastTRatio[0]=1.;                      << 171         lastT[0]=1.;
188         for(G4int j=1; j<=lastNRatio; j++)     << 172         for(G4int j=1; j<=lastN; j++)       // Calculate LogTab values
189         {                                         173         {
190             sv+=ds;                               174             sv+=ds;
191             lastTRatio[j]=CalcQF2IN_Ratio(sv,A << 175             lastT[j]=CalcQF2IN_Ratio(sv,A);
192         }                                         176         }
193         lastLRatio=new G4double[mls];          << 177         lastL=new G4double[mls];            // Create the logarithmic Table
194         // Initialize the logarithmic Table    << 178         if(m_s>sma)                           // Initialize the logarithmic Table
195   for(G4int j=0; j<mls; ++j) lastLRatio[j]=0.0 << 179         {
196         if(m_s>sma)                            << 180     G4double ls=std::log(m_s);
197         {                                      << 181             lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
198     G4double ls=G4Log(m_s);                    << 182             if(lastK>nls)
199             lastKRatio = static_cast<int>((ls- << 
200             if(lastKRatio>nls)                 << 
201             {                                     183             {
202                 lastKRatio=nls;                << 184                 lastK=nls;
203                 lastMRatio=lsa-lsi;            << 185                 lastM=lsa-lsi;
204             }                                     186             }
205             else lastMRatio = lastKRatio*dls;  << 187             else lastM = lastK*dl;            // Calculate max initialized ln(s)-lsi for LogTab
206             sv=mi;                                188             sv=mi;
207             for(G4int j=0; j<=lastKRatio; j++) << 189             for(G4int j=0; j<=lastK; j++)     // Calculate LogTab values
208             {                                     190             {
209                 lastLRatio[j]=CalcQF2IN_Ratio( << 191                 lastL[j]=CalcQF2IN_Ratio(sv,A);
210                 if(j!=lastKRatio) sv*=edls;    << 192                 if(j!=lastK) sv*=edl;
211             }                                     193             }
212         }                                         194         }
213         else                                //    195         else                                // LogTab is not initialized
214         {                                         196         {
215             lastKRatio = 0;                    << 197             lastK = 0;
216             lastMRatio = 0.;                   << 198             lastM = 0.;
217         }                                         199         }
218         i++;                                //    200         i++;                                // Make a new record to AMDB and position on it
219         vARatio.push_back(lastARatio);         << 201         vA.push_back(lastA);
220         vHRatio.push_back(lastHRatio);         << 202         vH.push_back(lastH);
221         vNRatio.push_back(lastNRatio);         << 203         vN.push_back(lastN);
222         vMRatio.push_back(lastMRatio);         << 204         vM.push_back(lastM);
223         vKRatio.push_back(lastKRatio);         << 205         vK.push_back(lastK);
224         vT->push_back(lastTRatio);             << 206         vT.push_back(lastT);
225         vL->push_back(lastLRatio);             << 207         vL.push_back(lastL);
226     }                                             208     }
227     else                                  // T    209     else                                  // The A value was found in AMDB
228     {                                             210     {
229         lastARatio=vARatio[i];                 << 211         lastA=vA[i];
230         lastHRatio=vHRatio[i];                 << 212         lastH=vH[i];
231         lastNRatio=vNRatio[i];                 << 213         lastN=vN[i];
232         lastMRatio=vMRatio[i];                 << 214         lastM=vM[i];
233         lastKRatio=vKRatio[i];                 << 215         lastK=vK[i];
234         lastTRatio=(*vT)[i];                   << 216         lastT=vT[i];
235         lastLRatio=(*vL)[i];                   << 217         lastL=vL[i];
236         if(m_s>lastHRatio)                     << 218         if(m_s>lastH)                          // At least LinTab must be updated
237         {                                         219         {
238             G4int nextN=lastNRatio+1;          << 220             G4int nextN=lastN+1;               // The next bin to be initialized
239             if(lastNRatio<nps)                 << 221             if(lastN<nps)
240             {                                     222             {
241         G4double sv=lastHRatio; // bug fix by  << 223         G4double sv=lastH; // bug fix by WP
242                                                   224 
243                 lastNRatio = static_cast<int>( << 225                 lastN = static_cast<int>(m_s/ds)+1;// MaxBin to be initialized
244                 if(lastNRatio>nps)             << 226                 if(lastN>nps)
245                 {                                 227                 {
246                     lastNRatio=nps;            << 228                     lastN=nps;
247                     lastHRatio=sma;            << 229                     lastH=sma;
248                 }                                 230                 }
249                 else lastHRatio = lastNRatio*d << 231                 else lastH = lastN*ds;           // Calculate max initialized s for LinTab
250                                                   232 
251                 for(G4int j=nextN; j<=lastNRat << 233                 for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values
252                 {                                 234                 {
253                     sv+=ds;                       235                     sv+=ds;
254                     lastTRatio[j]=CalcQF2IN_Ra << 236                     lastT[j]=CalcQF2IN_Ratio(sv,A);
255                 }                                 237                 }
256             } // End of LinTab update             238             } // End of LinTab update
257             if(lastNRatio>=nextN)              << 239             if(lastN>=nextN)
258             {                                     240             {
259                 vHRatio[i]=lastHRatio;         << 241                 vH[i]=lastH;
260                 vNRatio[i]=lastNRatio;         << 242                 vN[i]=lastN;
261             }                                     243             }
262             G4int nextK=lastKRatio+1;          << 244             G4int nextK=lastK+1;
263             if(!lastKRatio) nextK=0;           << 245             if(!lastK) nextK=0;
264             if(m_s>sma && lastKRatio<nls)      << 246             if(m_s>sma && lastK<nls)             // LogTab must be updated
265             {                                     247             {
266                 G4double sv=G4Exp(lastMRatio+l << 248                 G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed)
267                 G4double ls=G4Log(m_s);        << 249                 G4double ls=std::log(m_s);
268                 lastKRatio = static_cast<int>( << 250                 lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
269                 if(lastKRatio>nls)             << 251                 if(lastK>nls)
270                 {                                 252                 {
271                     lastKRatio=nls;            << 253                     lastK=nls;
272                     lastMRatio=lsa-lsi;        << 254                     lastM=lsa-lsi;
273                 }                                 255                 }
274                 else lastMRatio = lastKRatio*d << 256                 else lastM = lastK*dl;           // Calculate max initialized ln(s)-lsi for LogTab
275                 for(G4int j=nextK; j<=lastKRat << 257                 for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
276                 {                                 258                 {
277                     sv*=edls;                  << 259                     sv*=edl;
278                     lastLRatio[j]=CalcQF2IN_Ra << 260                     lastL[j]=CalcQF2IN_Ratio(sv,A);
279                 }                                 261                 }
280             } // End of LogTab update             262             } // End of LogTab update
281             if(lastKRatio>=nextK)              << 263             if(lastK>=nextK)
282             {                                     264             {
283                 vMRatio[i]=lastMRatio;         << 265                 vM[i]=lastM;
284                 vKRatio[i]=lastKRatio;         << 266                 vK[i]=lastK;
285             }                                     267             }
286         }                                         268         }
287     }                                             269     }
288     // Now one can use tabeles to calculate th    270     // Now one can use tabeles to calculate the value
289     if(m_s<sma)                             //    271     if(m_s<sma)                             // Use linear table
290     {                                             272     {
291         G4int n=static_cast<int>(m_s/ds);         273         G4int n=static_cast<int>(m_s/ds);     // Low edge number of the bin
292         G4double d=m_s-n*ds;                      274         G4double d=m_s-n*ds;                  // Linear shift
293         G4double v=lastTRatio[n];              << 275         G4double v=lastT[n];                // Base
294         lastRRatio=v+d*(lastTRatio[n+1]-v)/ds; << 276         lastR=v+d*(lastT[n+1]-v)/ds;        // Result
295     }                                             277     }
296     else                                  // U    278     else                                  // Use log table
297     {                                             279     {
298         G4double ls=G4Log(m_s)-lsi;        //  << 280         G4double ls=std::log(m_s)-lsi;        // ln(s)-l_min
299         G4int n=static_cast<int>(ls/dls);    / << 281         G4int n=static_cast<int>(ls/dl);    // Low edge number of the bin
300         G4double d=ls-n*dls;                 / << 282         G4double d=ls-n*dl;                 // Log shift
301         G4double v=lastLRatio[n];              << 283         G4double v=lastL[n];                // Base
302         lastRRatio=v+d*(lastLRatio[n+1]-v)/dls << 284         lastR=v+d*(lastL[n+1]-v)/dl;        // Result
303     }                                          << 285     }
304     if(lastRRatio<0.) lastRRatio=0.;           << 286     if(lastR<0.) lastR=0.;
305     if(lastRRatio>1.) lastRRatio=1.;           << 287     if(lastR>1.) lastR=1.;
306     return lastRRatio;                         << 288     return lastR;
307 } // End of CalcQF2IN_Ratio                       289 } // End of CalcQF2IN_Ratio
308                                                   290 
309 // Calculatio QasiFree/Inelastic Ratio as a fu    291 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A
310 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4do    292 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4double m_s, G4int A)
311 {                                                 293 {
                                                   >> 294     static const G4double C=1.246;
312     G4double s2=m_s*m_s;                          295     G4double s2=m_s*m_s;
313     G4double s4=s2*s2;                            296     G4double s4=s2*s2;
314     G4double ss=std::sqrt(std::sqrt(m_s));        297     G4double ss=std::sqrt(std::sqrt(m_s));
315     G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2    298     G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
316     G4double E=.2644+.016/(1.+G4Exp((29.54-m_s << 299     G4double E=.2644+.016/(1.+std::exp((29.54-m_s)/2.49));
317     G4double F=ss*.1526*G4Exp(-s2*ss*.0000859) << 300     G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
318     return C*G4Exp(-E*G4Pow::GetInstance()->po << 301     return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P);
319 } // End of CalcQF2IN_Ratio                       302 } // End of CalcQF2IN_Ratio
320                                                   303 
321 // Calculatio pair(hN_el,hN_tot) (mb): p in Ge    304 // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot)
322 std::pair<G4double,G4double> G4QuasiElRatios::    305 std::pair<G4double,G4double> G4QuasiElRatios::CalcElTot(G4double p, G4int I)
323 {                                                 306 {
324     // ---------> Each parameter set can have     307     // ---------> Each parameter set can have not more than nPoints=128 parameters
325                                                << 308     static const G4double lmi=3.5;       // min of (lnP-lmi)^2 parabola
                                                   >> 309     static const G4double pbe=.0557;     // elastic (lnP-lmi)^2 parabola coefficient
                                                   >> 310     static const G4double pbt=.3;        // total (lnP-lmi)^2 parabola coefficient
                                                   >> 311     static const G4double pmi=.1;        // Below that fast LE calculation is made
                                                   >> 312     static const G4double pma=1000.;     // Above that fast HE calculation is made
326     G4double El=0.;                      // pr    313     G4double El=0.;                      // prototype of the elastic hN cross-section
327     G4double To=0.;                      // pr    314     G4double To=0.;                      // prototype of the total hN cross-section
328     if(p<=0.)                                     315     if(p<=0.)
329     {                                             316     {
330         G4cout<<"-Warning-G4QuasiElRatios::Cal    317         G4cout<<"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl;
331         return std::make_pair(El,To);             318         return std::make_pair(El,To);
332     }                                             319     }
333     if     (!I)                          // pp    320     if     (!I)                          // pp/nn
334     {                                             321     {
335         if(p<pmi)                                 322         if(p<pmi)
336         {                                         323         {
337             G4double p2=p*p;                      324             G4double p2=p*p;
338             El=1./(.00012+p2*.2);                 325             El=1./(.00012+p2*.2);
339             To=El;                                326             To=El;
340         }                                         327         }
341         else if(p>pma)                            328         else if(p>pma)
342         {                                         329         {
343             G4double lp=G4Log(p)-lmi;          << 330             G4double lp=std::log(p)-lmi;
344             G4double lp2=lp*lp;                   331             G4double lp2=lp*lp;
345             El=pbe*lp2+6.72;                      332             El=pbe*lp2+6.72;
346             To=pbt*lp2+38.2;                      333             To=pbt*lp2+38.2;
347         }                                         334         }
348         else                                      335         else
349         {                                         336         {
350             G4double p2=p*p;                      337             G4double p2=p*p;
351             G4double LE=1./(.00012+p2*.2);        338             G4double LE=1./(.00012+p2*.2);
352             G4double lp=G4Log(p)-lmi;          << 339             G4double lp=std::log(p)-lmi;
353             G4double lp2=lp*lp;                   340             G4double lp2=lp*lp;
354             G4double rp2=1./p2;                   341             G4double rp2=1./p2;
355             El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp    342             El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
356             To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+    343             To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
357         }                                         344         }
358     }                                             345     }
359     else if(I==1)                        // np    346     else if(I==1)                        // np/pn
360     {                                             347     {
361         if(p<pmi)                                 348         if(p<pmi)
362         {                                         349         {
363             G4double p2=p*p;                      350             G4double p2=p*p;
364             El=1./(.00012+p2*(.051+.1*p2));       351             El=1./(.00012+p2*(.051+.1*p2));
365             To=El;                                352             To=El;
366         }                                         353         }
367         else if(p>pma)                            354         else if(p>pma)
368         {                                         355         {
369             G4double lp=G4Log(p)-lmi;          << 356             G4double lp=std::log(p)-lmi;
370             G4double lp2=lp*lp;                   357             G4double lp2=lp*lp;
371             El=pbe*lp2+6.72;                      358             El=pbe*lp2+6.72;
372             To=pbt*lp2+38.2;                      359             To=pbt*lp2+38.2;
373         }                                         360         }
374         else                                      361         else
375         {                                         362         {
376             G4double p2=p*p;                      363             G4double p2=p*p;
377             G4double LE=1./(.00012+p2*(.051+.1    364             G4double LE=1./(.00012+p2*(.051+.1*p2));
378             G4double lp=G4Log(p)-lmi;          << 365             G4double lp=std::log(p)-lmi;
379             G4double lp2=lp*lp;                   366             G4double lp2=lp*lp;
380             G4double rp2=1./p2;                   367             G4double rp2=1./p2;
381             El=LE+(pbe*lp2+6.72+30./p)/(1.+.49    368             El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
382             To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*r    369             To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
383         }                                         370         }
384     }                                             371     }
385     else if(I==2)                        // pi    372     else if(I==2)                        // pimp/pipn
386     {                                             373     {
387         G4double lp=G4Log(p);                  << 374         G4double lp=std::log(p);
388         if(p<pmi)                                 375         if(p<pmi)
389         {                                         376         {
390             G4double lr=lp+1.27;                  377             G4double lr=lp+1.27;
391             El=1.53/(lr*lr+.0676);                378             El=1.53/(lr*lr+.0676);
392             To=El*3;                              379             To=El*3;
393         }                                         380         }
394         else if(p>pma)                            381         else if(p>pma)
395         {                                         382         {
396             G4double ld=lp-lmi;                   383             G4double ld=lp-lmi;
397             G4double ld2=ld*ld;                   384             G4double ld2=ld*ld;
398             G4double sp=std::sqrt(p);             385             G4double sp=std::sqrt(p);
399             El=pbe*ld2+2.4+7./sp;                 386             El=pbe*ld2+2.4+7./sp;
400             To=pbt*ld2+22.3+12./sp;               387             To=pbt*ld2+22.3+12./sp;
401         }                                         388         }
402         else                                      389         else
403         {                                         390         {
404             G4double lr=lp+1.27;                  391             G4double lr=lp+1.27;                    // p1
405             G4double LE=1.53/(lr*lr+.0676);       392             G4double LE=1.53/(lr*lr+.0676);         // p2, p3        
406             G4double ld=lp-lmi;                   393             G4double ld=lp-lmi;                     // p4 (lmi=3.5)
407             G4double ld2=ld*ld;                   394             G4double ld2=ld*ld;
408             G4double p2=p*p;                      395             G4double p2=p*p;
409             G4double p4=p2*p2;                    396             G4double p4=p2*p2;
410             G4double sp=std::sqrt(p);             397             G4double sp=std::sqrt(p);
411             G4double lm=lp+.36;                   398             G4double lm=lp+.36;                     // p5
412             G4double md=lm*lm+.04;                399             G4double md=lm*lm+.04;                  // p6
413             G4double lh=lp-.017;                  400             G4double lh=lp-.017;                    // p7
414             G4double hd=lh*lh+.0025;              401             G4double hd=lh*lh+.0025;                // p8
415             El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p    402             El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
416             To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+    403             To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
417         }                                         404         }
418     }                                             405     }
419     else if(I==3)                        // pi    406     else if(I==3)                        // pipp/pimn
420     {                                             407     {
421         G4double lp=G4Log(p);                  << 408         G4double lp=std::log(p);
422         if(p<pmi)                                 409         if(p<pmi)
423         {                                         410         {
424             G4double lr=lp+1.27;                  411             G4double lr=lp+1.27;
425             G4double lr2=lr*lr;                   412             G4double lr2=lr*lr;
426             El=13./(lr2+lr2*lr2+.0676);           413             El=13./(lr2+lr2*lr2+.0676);
427             To=El;                                414             To=El;
428         }                                         415         }
429         else if(p>pma)                            416         else if(p>pma)
430         {                                         417         {
431             G4double ld=lp-lmi;                   418             G4double ld=lp-lmi;
432             G4double ld2=ld*ld;                   419             G4double ld2=ld*ld;
433             G4double sp=std::sqrt(p);             420             G4double sp=std::sqrt(p);
434             El=pbe*ld2+2.4+6./sp;                 421             El=pbe*ld2+2.4+6./sp;
435             To=pbt*ld2+22.3+5./sp;                422             To=pbt*ld2+22.3+5./sp;
436         }                                         423         }
437         else                                      424         else
438         {                                         425         {
439             G4double lr=lp+1.27;                  426             G4double lr=lp+1.27;                   // p1
440             G4double lr2=lr*lr;                   427             G4double lr2=lr*lr;
441             G4double LE=13./(lr2+lr2*lr2+.0676    428             G4double LE=13./(lr2+lr2*lr2+.0676);   // p2, p3
442             G4double ld=lp-lmi;                   429             G4double ld=lp-lmi;                    // p4 (lmi=3.5)
443             G4double ld2=ld*ld;                   430             G4double ld2=ld*ld;
444             G4double p2=p*p;                      431             G4double p2=p*p;
445             G4double p4=p2*p2;                    432             G4double p4=p2*p2;
446             G4double sp=std::sqrt(p);             433             G4double sp=std::sqrt(p);
447             G4double lm=lp-.32;                   434             G4double lm=lp-.32;                    // p5
448             G4double md=lm*lm+.0576;              435             G4double md=lm*lm+.0576;               // p6
449             El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p    436             El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
450             To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./    437             To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
451         }                                         438         }
452     }                                             439     }
453     else if(I==4)                        // Km    440     else if(I==4)                        // Kmp/Kmn/K0p/K0n
454     {                                             441     {
455                                                   442         
456         if(p<pmi)                                 443         if(p<pmi)
457         {                                         444         {
458             G4double psp=p*std::sqrt(p);          445             G4double psp=p*std::sqrt(p);
459             El=5.2/psp;                           446             El=5.2/psp;
460             To=14./psp;                           447             To=14./psp;
461         }                                         448         }
462         else if(p>pma)                            449         else if(p>pma)
463         {                                         450         {
464             G4double ld=G4Log(p)-lmi;          << 451             G4double ld=std::log(p)-lmi;
465             G4double ld2=ld*ld;                   452             G4double ld2=ld*ld;
466             El=pbe*ld2+2.23;                      453             El=pbe*ld2+2.23;
467             To=pbt*ld2+19.5;                      454             To=pbt*ld2+19.5;
468         }                                         455         }
469         else                                      456         else
470         {                                         457         {
471             G4double ld=G4Log(p)-lmi;          << 458             G4double ld=std::log(p)-lmi;
472             G4double ld2=ld*ld;                   459             G4double ld2=ld*ld;
473             G4double sp=std::sqrt(p);             460             G4double sp=std::sqrt(p);
474             G4double psp=p*sp;                    461             G4double psp=p*sp;
475             G4double p2=p*p;                      462             G4double p2=p*p;
476             G4double p4=p2*p2;                    463             G4double p4=p2*p2;
477             G4double lm=p-.39;                    464             G4double lm=p-.39;
478             G4double md=lm*lm+.000156;            465             G4double md=lm*lm+.000156;
479             G4double lh=p-1.;                     466             G4double lh=p-1.;
480             G4double hd=lh*lh+.0156;              467             G4double hd=lh*lh+.0156;
481             El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/s    468             El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
482             To=14./psp+(pbt*ld2+19.5)/(1.-.21/    469             To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
483         }                                         470         }
484     }                                             471     }
485     else if(I==5)                        // Kp    472     else if(I==5)                        // Kpp/Kpn/aKp/aKn
486     {                                             473     {
487         if(p<pmi)                                 474         if(p<pmi)
488         {                                         475         {
489             G4double lr=p-.38;                    476             G4double lr=p-.38;
490             G4double lm=p-1.;                     477             G4double lm=p-1.;
491             G4double md=lm*lm+.372;               478             G4double md=lm*lm+.372;   
492             El=.7/(lr*lr+.0676)+2./md;            479             El=.7/(lr*lr+.0676)+2./md;
493             To=El+.6/md;                          480             To=El+.6/md;
494         }                                         481         }
495         else if(p>pma)                            482         else if(p>pma)
496         {                                         483         {
497             G4double ld=G4Log(p)-lmi;          << 484             G4double ld=std::log(p)-lmi;
498             G4double ld2=ld*ld;                   485             G4double ld2=ld*ld;
499             El=pbe*ld2+2.23;                      486             El=pbe*ld2+2.23;
500             To=pbt*ld2+19.5;                      487             To=pbt*ld2+19.5;
501         }                                         488         }
502         else                                      489         else
503         {                                         490         {
504             G4double ld=G4Log(p)-lmi;          << 491             G4double ld=std::log(p)-lmi;
505             G4double ld2=ld*ld;                   492             G4double ld2=ld*ld;
506             G4double lr=p-.38;                    493             G4double lr=p-.38;
507             G4double LE=.7/(lr*lr+.0676);         494             G4double LE=.7/(lr*lr+.0676);
508             G4double sp=std::sqrt(p);             495             G4double sp=std::sqrt(p);
509             G4double p2=p*p;                      496             G4double p2=p*p;
510             G4double p4=p2*p2;                    497             G4double p4=p2*p2;
511             G4double lm=p-1.;                     498             G4double lm=p-1.;
512             G4double md=lm*lm+.372;               499             G4double md=lm*lm+.372;
513             El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/    500             El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
514             To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.    501             To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
515         }                                         502         }
516     }                                             503     }
517     else if(I==6)                        // hy    504     else if(I==6)                        // hyperon-N
518     {                                             505     {
519         if(p<pmi)                                 506         if(p<pmi)
520         {                                         507         {
521             G4double p2=p*p;                      508             G4double p2=p*p;
522             El=1./(.002+p2*(.12+p2));             509             El=1./(.002+p2*(.12+p2));
523             To=El;                                510             To=El;
524         }                                         511         }
525         else if(p>pma)                            512         else if(p>pma)
526         {                                         513         {
527             G4double lp=G4Log(p)-lmi;          << 514             G4double lp=std::log(p)-lmi;
528             G4double lp2=lp*lp;                   515             G4double lp2=lp*lp;
529             G4double sp=std::sqrt(p);             516             G4double sp=std::sqrt(p);
530             El=(pbe*lp2+6.72)/(1.+2./sp);         517             El=(pbe*lp2+6.72)/(1.+2./sp);
531             To=(pbt*lp2+38.2+900./sp)/(1.+27./    518             To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
532         }                                         519         }
533         else                                      520         else
534         {                                         521         {
535             G4double p2=p*p;                      522             G4double p2=p*p;
536             G4double LE=1./(.002+p2*(.12+p2));    523             G4double LE=1./(.002+p2*(.12+p2));
537             G4double lp=G4Log(p)-lmi;          << 524             G4double lp=std::log(p)-lmi;
538             G4double lp2=lp*lp;                   525             G4double lp2=lp*lp;
539             G4double p4=p2*p2;                    526             G4double p4=p2*p2;
540             G4double sp=std::sqrt(p);             527             G4double sp=std::sqrt(p);
541             El=LE+(pbe*lp2+6.72+99./p2)/(1.+2.    528             El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
542             To=LE+(pbt*lp2+38.2+900./sp)/(1.+2    529             To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
543         }                                         530         }
544     }                                             531     }
545     else if(I==7)                        // an    532     else if(I==7)                        // antibaryon-N
546     {                                             533     {
547         if(p>pma)                                 534         if(p>pma)
548         {                                         535         {
549             G4double lp=G4Log(p)-lmi;          << 536             G4double lp=std::log(p)-lmi;
550             G4double lp2=lp*lp;                   537             G4double lp2=lp*lp;
551             El=pbe*lp2+6.72;                      538             El=pbe*lp2+6.72;
552             To=pbt*lp2+38.2;                      539             To=pbt*lp2+38.2;
553         }                                         540         }
554         else                                      541         else
555         {                                         542         {
556             G4double ye=G4Pow::GetInstance()-> << 543             G4double ye=std::pow(p,1.25);
557             G4double yt=G4Pow::GetInstance()-> << 544             G4double yt=std::pow(p,.35);
558             G4double lp=G4Log(p)-lmi;          << 545             G4double lp=std::log(p)-lmi;
559             G4double lp2=lp*lp;                   546             G4double lp2=lp*lp;
560             El=80./(ye+1.)+pbe*lp2+6.72;          547             El=80./(ye+1.)+pbe*lp2+6.72;
561             To=(80./yt+.3)/yt+pbt*lp2+38.2;       548             To=(80./yt+.3)/yt+pbt*lp2+38.2;
562         }                                         549         }
563     }                                             550     }
564     else                                          551     else
565     {                                             552     {
566         G4cout<<"*Error*G4QuasiElRatios::CalcE    553         G4cout<<"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
567         G4Exception("G4QuasiElRatios::CalcElTo    554         G4Exception("G4QuasiElRatios::CalcElTot:","23",FatalException,"QEcrash");
568     }                                             555     }
569     if(El>To) El=To;                              556     if(El>To) El=To;
570     return std::make_pair(El,To);                 557     return std::make_pair(El,To);
571 } // End of CalcElTot                             558 } // End of CalcElTot
572                                                   559 
573 // For hadron PDG with momentum Mom (GeV/c) on    560 // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
574 std::pair<G4double,G4double> G4QuasiElRatios::    561 std::pair<G4double,G4double> G4QuasiElRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
575 {                                                 562 {
576     G4int ind=0;                                  563     G4int ind=0;                                 // Prototype of the reaction index
577     G4bool kfl=true;                              564     G4bool kfl=true;                             // Flag of K0/aK0 oscillation
578     G4bool kf=false;                              565     G4bool kf=false;
579     if(PDG==130||PDG==310)                        566     if(PDG==130||PDG==310)
580     {                                             567     {
581         kf=true;                                  568         kf=true;
582         if(G4UniformRand()>.5) kfl=false;         569         if(G4UniformRand()>.5) kfl=false;
583     }                                             570     }
584     if      ( (PDG == 2212 && F) || (PDG == 21    571     if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
585     else if ( (PDG == 2112 && F) || (PDG == 22    572     else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
586     else if ( (PDG == -211 && F) || (PDG == 21    573     else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
587     else if ( (PDG == 211 && F) || (PDG == -21    574     else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
588     //AR-Jul2020: Extended to charmed and bott << 575     else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
589     // - treat mesons with constituent c quark << 576     else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
590     // - treat mesons with constituent cbar an << 577     else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
591     // - treat all heavy baryons (i.e. hyperon << 578     else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
592     // - treat all heavy anti-baryons (i.e. an << 
593     else if ( PDG == -321 || PDG == -311 || (k << 
594               PDG ==  411 || PDG ==  421 || PD << 
595         PDG == -521 || PDG == -511 || PDG == - << 
596     else if ( PDG ==  321 || PDG ==  311 || (k << 
597               PDG == -411 || PDG == -421 || PD << 
598         PDG ==  521 || PDG ==  511 || PDG ==   << 
599     else if ( PDG >  3000 && PDG <  5333 ) ind << 
600     else if ( PDG > -5333 && PDG < -2000 ) ind << 
601     else {                                        579     else {
602         G4cout<<"*Error*G4QuasiElRatios::CalcE    580         G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
603         <<", while it is defined only for p,n,    581         <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
604         G4Exception("G4QuasiElRatio::CalcElTot    582         G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash");
605     }                                             583     }
606     return CalcElTot(p,ind);                      584     return CalcElTot(p,ind);
607 }                                                 585 }
608                                                   586 
609 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV    587 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron
610 std::pair<G4double,G4double> G4QuasiElRatios::    588 std::pair<G4double,G4double> G4QuasiElRatios::FetchElTot(G4double p, G4int PDG, G4bool F)
611 {                                                 589 {
                                                   >> 590     static const G4int    nlp=300;         // Number of steps in the S(lnp) logTable(5% step)
                                                   >> 591     static const G4int    mlp=nlp+1;       // Number of elements in the S(lnp) logTable
                                                   >> 592     static const G4double lpi=-5.;         // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
                                                   >> 593     static const G4double lpa=10.;         // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
                                                   >> 594     static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
                                                   >> 595     static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
                                                   >> 596     static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
                                                   >> 597     static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
                                                   >> 598     //static const G4double toler=.001;      // Relative Tolarence defining "theSameMomentum"
                                                   >> 599     static G4double lastP=0.;              // The last momentum for which XS was calculated
                                                   >> 600     static G4int    lastH=0;               // The last projPDG for which XS was calculated
                                                   >> 601     static G4bool   lastF=true;            // The last nucleon for which XS was calculated
                                                   >> 602     static std::pair<G4double,G4double> lastR(0.,0.); // The last result
                                                   >> 603     // Local Associative Data Base:
                                                   >> 604     static std::vector<G4int>     vI;      // Vector of index for which XS was calculated
                                                   >> 605     static std::vector<G4double>  vM;      // Vector of rel max ln(p) initialized in LogTable
                                                   >> 606     static std::vector<G4int>     vK;      // Vector of topBin number initialized in LogTable
                                                   >> 607     // Last values of the Associative Data Base:
                                                   >> 608     static G4int     lastI=0;              // The Last index for which XS was calculated
                                                   >> 609     static G4double  lastM=0.;             // The Last rel max ln(p) initialized in LogTable
                                                   >> 610     static G4int     lastK=0;             // The Last topBin number initialized in LogTable
                                                   >> 611     static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
612     // LogTable is created only if necessary.     612     // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
613     G4int nDB=(G4int)vItot.size();  // A numbe << 613     G4int nDB=vI.size();                   // A number of hadrons already initialized in AMDB
614     if(nDB && lastHtot==PDG && lastFtot==F &&  << 614     if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
615     //  if(nDB && lastH==PDG && lastF==F && p>    615     //  if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
616     lastHtot=PDG;                              << 616     lastH=PDG;
617     lastFtot=F;                                << 617     lastF=F;
618     G4int ind=-1;                          //     618     G4int ind=-1;                          // Prototipe of the index of the PDG/F combination
619     // i=0: pp(nn), i=1: np(pn), i=2: pimp(pip    619     // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
620     // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i    620     // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) 
621     G4bool kfl=true;                              621     G4bool kfl=true;                             // Flag of K0/aK0 oscillation
622     G4bool kf=false;                              622     G4bool kf=false;
623     if(PDG==130||PDG==310)                        623     if(PDG==130||PDG==310)
624     {                                             624     {
625         kf=true;                                  625         kf=true;
626         if(G4UniformRand()>.5) kfl=false;         626         if(G4UniformRand()>.5) kfl=false;
627     }                                             627     }
628     if      ( (PDG == 2212 && F) || (PDG == 21    628     if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
629     else if ( (PDG == 2112 && F) || (PDG == 22    629     else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
630     else if ( (PDG == -211 && F) || (PDG == 21    630     else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
631     else if ( (PDG == 211 && F) || (PDG == -21    631     else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
632     //AR-Jul2020: Extended to charmed and bott << 632     else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
633     // - treat mesons with constituent c quark << 633     else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
634     // - treat mesons with constituent cbar an << 634     else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
635     // - treat all heavy baryons (i.e. hyperon << 635     else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
636     // - treat all heavy anti-baryons (i.e. an << 
637     else if ( PDG == -321 || PDG == -311 || (k << 
638               PDG ==  411 || PDG ==  421 || PD << 
639         PDG == -521 || PDG == -511 || PDG == - << 
640     else if ( PDG ==  321 || PDG ==  311 || (k << 
641               PDG == -411 || PDG == -421 || PD << 
642         PDG ==  521 || PDG ==  511 || PDG ==   << 
643     else if ( PDG >  3000 && PDG <  5333 ) ind << 
644     else if ( PDG > -5333 && PDG < -2000 ) ind << 
645     else {                                        636     else {
646         G4cout<<"*Error*G4QuasiElRatios::Fetch    637         G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
647         <<", while it is defined only for p,n,    638         <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
648         G4Exception("G4QuasiELRatio::FetchElTo    639         G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash");
649     }                                             640     }
650     if(nDB && lastItot==ind && p>0. && p==last << 641     if(nDB && lastI==ind && p>0. && p==lastP) return lastR;  // VI do not use toler
651     //  if(nDB && lastI==ind && p>0. && std::f    642     //  if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
652     if(p<=mip || p>=map) return CalcElTot(p,in << 643     if(p<=mi || p>=ma) return CalcElTot(p,ind);   // @@ Slow calculation ! (Warning?)
653     G4bool found=false;                           644     G4bool found=false;
654     G4int i=-1;                                   645     G4int i=-1;
655     if(nDB) for (i=0; i<nDB; i++) if(ind==vIto << 646     if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i])  // Sirch for this index in AMDB
656     {                                             647     {
657         found=true;                               648         found=true;                                 // The index is found
658         break;                                    649         break;
659     }                                             650     }
660     G4double lp=G4Log(p);                      << 651     G4double lp=std::log(p);
661     if(!nDB || !found)                            652     if(!nDB || !found)                            // Create new line in the AMDB
662     {                                             653     {
663         lastXtot = new std::pair<G4double,G4do << 654         lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
664         lastItot = ind;                        << 655         lastI = ind;                                // Remember the initialized inex
665         lastKtot = static_cast<int>((lp-lpi)/d << 656         lastK = static_cast<int>((lp-lpi)/dl)+1;    // MaxBin to be initialized in LogTaB
666         if(lastKtot>nlp)                       << 657         if(lastK>nlp)
667         {                                      << 658         {
668             lastKtot=nlp;                      << 659             lastK=nlp;
669             lastMtot=lpa-lpi;                  << 660             lastM=lpa-lpi;
670         }                                      << 661         }
671         else lastMtot = lastKtot*dlp;          << 662         else lastM = lastK*dl;               // Calculate max initialized ln(p)-lpi for LogTab
672         G4double pv=mip;                       << 663         G4double pv=mi;
673         for(G4int j=0; j<=lastKtot; j++)       << 664         for(G4int j=0; j<=lastK; j++)        // Calculate LogTab values
674         {                                         665         {
675             lastXtot[j]=CalcElTot(pv,ind);     << 666             lastX[j]=CalcElTot(pv,ind);
676             if(j!=lastKtot) pv*=edlp;          << 667             if(j!=lastK) pv*=edl;
677         }                                         668         }
678         i++;                                 /    669         i++;                                 // Make a new record to AMDB and position on it
679         vItot.push_back(lastItot);             << 670         vI.push_back(lastI);
680         vMtot.push_back(lastMtot);             << 671         vM.push_back(lastM);
681         vKtot.push_back(lastKtot);             << 672         vK.push_back(lastK);
682         vX->push_back(lastXtot);               << 673         vX.push_back(lastX);
683     }                                             674     }
684     else                                   //     675     else                                   // The A value was found in AMDB
685     {                                             676     {
686         lastItot=vItot[i];                     << 677         lastI=vI[i];
687         lastMtot=vMtot[i];                     << 678         lastM=vM[i];
688         lastKtot=vKtot[i];                     << 679         lastK=vK[i];
689         lastXtot=(*vX)[i];                     << 680         lastX=vX[i];
690         G4int nextK=lastKtot+1;                << 681         G4int nextK=lastK+1;
691         G4double lpM=lastMtot+lpi;             << 682         G4double lpM=lastM+lpi;
692         if(lp>lpM && lastKtot<nlp)             << 683         if(lp>lpM && lastK<nlp)              // LogTab must be updated
693         {                                         684         {
694             lastKtot = static_cast<int>((lp-lp << 685             lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
695             if(lastKtot>nlp)                   << 686             if(lastK>nlp)
696             {                                     687             {
697                 lastKtot=nlp;                  << 688                 lastK=nlp;
698                 lastMtot=lpa-lpi;              << 689                 lastM=lpa-lpi;
699             }                                     690             }
700             else lastMtot = lastKtot*dlp;      << 691             else lastM = lastK*dl;           // Calculate max initialized ln(p)-lpi for LogTab
701             G4double pv=G4Exp(lpM);       // m << 692             G4double pv=std::exp(lpM);       // momentum of the last calculated beam
702             for(G4int j=nextK; j<=lastKtot; j+ << 693             for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
703             {                                     694             {
704                 pv*=edlp;                      << 695                 pv*=edl;
705                 lastXtot[j]=CalcElTot(pv,ind); << 696                 lastX[j]=CalcElTot(pv,ind);
706             }                                     697             }
707         } // End of LogTab update                 698         } // End of LogTab update
708         if(lastKtot>=nextK)                    << 699         if(lastK>=nextK)                   // The AMDB was apdated
709         {                                         700         {
710             vMtot[i]=lastMtot;                 << 701             vM[i]=lastM;
711             vKtot[i]=lastKtot;                 << 702             vK[i]=lastK;
712         }                                         703         }
713     }                                             704     }
714     // Now one can use tabeles to calculate th    705     // Now one can use tabeles to calculate the value
715     G4double dlpp=lp-lpi;                      << 706     G4double dlp=lp-lpi;                       // Shifted log(p) value
716     G4int n=static_cast<int>(dlpp/dlp);        << 707     G4int n=static_cast<int>(dlp/dl);          // Low edge number of the bin
717     G4double d=dlpp-n*dlp;                     << 708     G4double d=dlp-n*dl;                       // Log shift
718     G4double e=lastXtot[n].first;              << 709     G4double e=lastX[n].first;                 // E-Base
719     lastRtot.first=e+d*(lastXtot[n+1].first-e) << 710     lastR.first=e+d*(lastX[n+1].first-e)/dl;   // E-Result
720     if(lastRtot.first<0.)  lastRtot.first = 0. << 711     if(lastR.first<0.)  lastR.first = 0.;
721     G4double t=lastXtot[n].second;             << 712     G4double t=lastX[n].second;                // T-Base
722     lastRtot.second=t+d*(lastXtot[n+1].second- << 713     lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
723     if(lastRtot.second<0.) lastRtot.second= 0. << 714     if(lastR.second<0.) lastR.second= 0.;
724     if(lastRtot.first>lastRtot.second) lastRto << 715     if(lastR.first>lastR.second) lastR.first = lastR.second;
725     return lastRtot;                           << 716     return lastR;
726 } // End of FetchElTot                            717 } // End of FetchElTot
727                                                   718 
728 // (Mean Elastic and Mean Total) Cross-Section    719 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
729 std::pair<G4double,G4double> G4QuasiElRatios::    720 std::pair<G4double,G4double> G4QuasiElRatios::GetElTot(G4double pIU, G4int hPDG,
730                                                   721                                                         G4int Z,       G4int N)
731 {                                                 722 {
732     G4double pGeV=pIU/gigaelectronvolt;           723     G4double pGeV=pIU/gigaelectronvolt;
733     if(Z<1 && N<1)                                724     if(Z<1 && N<1)
734     {                                             725     {
735         G4cout<<"-Warning-G4QuasiElRatio::GetE    726         G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
736         return std::make_pair(0.,0.);             727         return std::make_pair(0.,0.);
737     }                                             728     }
738     std::pair<G4double,G4double> hp=FetchElTot    729     std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
739     std::pair<G4double,G4double> hn=FetchElTot    730     std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
740     G4double A=(Z+N)/millibarn;                   731     G4double A=(Z+N)/millibarn;                // To make the result in independent units(IU)
741     return std::make_pair((Z*hp.first+N*hn.fir    732     return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
742 } // End of GetElTot                              733 } // End of GetElTot
743                                                   734 
744 // (Mean Elastic and Mean Total) Cross-Section    735 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
745 std::pair<G4double,G4double> G4QuasiElRatios::    736 std::pair<G4double,G4double> G4QuasiElRatios::GetChExFactor(G4double pIU, G4int hPDG,
746                                                   737                                                              G4int Z, G4int N)
747 {                                                 738 {
748     G4double pGeV=pIU/gigaelectronvolt;           739     G4double pGeV=pIU/gigaelectronvolt;
749     G4double resP=0.;                             740     G4double resP=0.;
750     G4double resN=0.;                             741     G4double resN=0.;
751     if(Z<1 && N<1)                                742     if(Z<1 && N<1)
752     {                                             743     {
753         G4cout<<"-Warning-G4QuasiElRatio::GetC    744         G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
754         return std::make_pair(resP,resN);         745         return std::make_pair(resP,resN);
755     }                                             746     }
756     G4double A=Z+N;                               747     G4double A=Z+N;
757     G4double pf=0.;                               748     G4double pf=0.;                              // Possibility to interact with a proton
758     G4double nf=0.;                               749     G4double nf=0.;                              // Possibility to interact with a neutron
759     if   (hPDG==-211||hPDG==-321||hPDG==3112||    750     if   (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
760     else if(hPDG==211||hPDG==321||hPDG==3222||    751     else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
761     else if(hPDG==-311||hPDG==311||hPDG==130||    752     else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
762     {                                             753     {
763         G4double dA=A+A;                          754         G4double dA=A+A;
764         pf=Z/(dA+N+N);                            755         pf=Z/(dA+N+N);
765         nf=N/(dA+Z+Z);                            756         nf=N/(dA+Z+Z);
766     }                                             757     }
767     G4double mult=1.;  // Factor of increasing    758     G4double mult=1.;  // Factor of increasing multiplicity ( ? @@)
768     if(pGeV>.5)                                   759     if(pGeV>.5)
769     {                                             760     {
770         mult=1./(1.+G4Log(pGeV+pGeV))/pGeV;    << 761         mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
771         if(mult>1.) mult=1.;                      762         if(mult>1.) mult=1.;
772     }                                             763     }
773     if(pf)                                        764     if(pf)
774     {                                             765     {
775         std::pair<G4double,G4double> hp=FetchE    766         std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
776         resP=pf*(hp.second/hp.first-1.)*mult;     767         resP=pf*(hp.second/hp.first-1.)*mult;
777     }                                             768     }
778     if(nf)                                        769     if(nf)
779     {                                             770     {
780         std::pair<G4double,G4double> hn=FetchE    771         std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
781         resN=nf*(hn.second/hn.first-1.)*mult;     772         resN=nf*(hn.second/hn.first-1.)*mult;
782     }                                             773     }
783     return std::make_pair(resP,resN);             774     return std::make_pair(resP,resN);
784 } // End of GetChExFactor                         775 } // End of GetChExFactor
785                                                   776 
786 // scatter (pPDG,p4M) on a virtual nucleon (NP    777 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
787 // if(newN4M.e()==0.) - below threshold, XS=0,    778 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
788 std::pair<G4LorentzVector,G4LorentzVector> G4Q    779 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::Scatter(G4int NPDG,
789                                                   780                                                                      G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
790 {                                              << 781 {    
791     static const G4double mNeut= G4Neutron::Ne    782     static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
792     static const G4double mProt= G4Proton::Pro    783     static const G4double mProt= G4Proton::Proton()->GetPDGMass();
793     static const G4double mDeut= G4Deuteron::D    784     static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass();
794     static const G4double mTrit= G4Triton::Tri    785     static const G4double mTrit= G4Triton::Triton()->GetPDGMass();
795     static const G4double mHel3= G4He3::He3()-    786     static const G4double mHel3= G4He3::He3()->GetPDGMass();
796     static const G4double mAlph= G4Alpha::Alph    787     static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass();
797                                                   788     
798     G4LorentzVector pr4M=p4M/megaelectronvolt;    789     G4LorentzVector pr4M=p4M/megaelectronvolt;   // Convert 4-momenta in MeV (keep p4M)
799     N4M/=megaelectronvolt;                        790     N4M/=megaelectronvolt;
800     G4LorentzVector tot4M=N4M+p4M;                791     G4LorentzVector tot4M=N4M+p4M;
801     G4double mT=mNeut;                            792     G4double mT=mNeut;
802     G4int Z=0;                                    793     G4int Z=0;
803     G4int N=1;                                    794     G4int N=1;
804     if(NPDG==2212||NPDG==90001000)                795     if(NPDG==2212||NPDG==90001000)
805     {                                             796     {
806         mT=mProt;                                 797         mT=mProt;
807         Z=1;                                      798         Z=1;
808         N=0;                                      799         N=0;
809     }                                             800     }
810     else if(NPDG==90001001)                       801     else if(NPDG==90001001)
811     {                                             802     {
812         mT=mDeut;                                 803         mT=mDeut;
813         Z=1;                                      804         Z=1;
814         N=1;                                      805         N=1;
815     }                                             806     }
816     else if(NPDG==90002001)                       807     else if(NPDG==90002001)
817     {                                             808     {
818         mT=mHel3;                                 809         mT=mHel3;
819         Z=2;                                      810         Z=2;
820         N=1;                                      811         N=1;
821     }                                             812     }
822     else if(NPDG==90001002)                       813     else if(NPDG==90001002)
823     {                                             814     {
824         mT=mTrit;                                 815         mT=mTrit;
825         Z=1;                                      816         Z=1;
826         N=2;                                      817         N=2;
827     }                                             818     }
828     else if(NPDG==90002002)                       819     else if(NPDG==90002002)
829     {                                             820     {
830         mT=mAlph;                                 821         mT=mAlph;
831         Z=2;                                      822         Z=2;
832         N=2;                                      823         N=2;
833     }                                             824     }
834     else if(NPDG!=2112&&NPDG!=90000001)           825     else if(NPDG!=2112&&NPDG!=90000001)
835     {                                             826     {
836         G4cout<<"Error:G4QuasiElRatios::Scatte    827         G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
837         G4Exception("G4QuasiElRatios::Scatter:    828         G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain");
838         //return std::make_pair(G4LorentzVecto    829         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
839     }                                             830     }
840     G4double mT2=mT*mT;                           831     G4double mT2=mT*mT;
841     G4double mP2=pr4M.m2();                       832     G4double mP2=pr4M.m2();
842     G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);      833     G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
843     G4double E2=E*E;                              834     G4double E2=E*E;
844     if(E<0. || E2<mP2)                            835     if(E<0. || E2<mP2)
845     {                                             836     {
846         return std::make_pair(G4LorentzVector(    837         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
847     }                                             838     }
848     G4double P=std::sqrt(E2-mP2);                 839     G4double P=std::sqrt(E2-mP2);                   // Momentum in pseudo laboratory system
849     // @@ Temporary NN t-dependence for all ha    840     // @@ Temporary NN t-dependence for all hadrons
850     //if(pPDG>3400 || pPDG<-3400) G4cout<<"-Wa << 841     if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl;
851     G4int PDG=2212;                               842     G4int PDG=2212;                                                // *TMP* instead of pPDG
852     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG    843     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
853     if(!Z && N==1)                 // Change f    844     if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
854     {                                             845     {
855         Z=1;                                      846         Z=1;
856         N=0;                                      847         N=0;
857         if     (PDG==2212) PDG=2112;              848         if     (PDG==2212) PDG=2112;
858         else if(PDG==2112) PDG=2212;              849         else if(PDG==2112) PDG=2212;
859     }                                             850     }
860     G4double xSec=0.;                        /    851     G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
861     if(PDG==2212) xSec=PCSmanager->GetChipsCro    852     if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
862     else          xSec=NCSmanager->GetChipsCro    853     else          xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
863     // @@ check a possibility to separate p, n    854     // @@ check a possibility to separate p, n, or alpha (!)
864     if(xSec <= 0.)                                855     if(xSec <= 0.)                                    // The cross-section iz 0 -> Do Nothing
865     {                                             856     {
866         return std::make_pair(G4LorentzVector(    857         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
867     }                                             858     }
868     G4double mint=0.;                        /    859     G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
869     if(PDG==2212) mint=PCSmanager->GetExchange    860     if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
870     else          mint=NCSmanager->GetExchange    861     else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
871     G4double maxt=0.;                             862     G4double maxt=0.;                                    // Prototype of max possible -t
872     if(PDG==2212) maxt=PCSmanager->GetHMaxT();    863     if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
873     else          maxt=NCSmanager->GetHMaxT();    864     else          maxt=NCSmanager->GetHMaxT();           // max possible -t
874     G4double cost=1.-(mint+mint)/maxt; // cos(    865     G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
875     if(cost>1. || cost<-1. || !(cost>-1. || co    866     if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
876     {                                             867     {
877         if     (cost>1.)  cost=1.;                868         if     (cost>1.)  cost=1.;
878         else if(cost<-1.) cost=-1.;               869         else if(cost<-1.) cost=-1.;
879         else                                      870         else
880         {                                         871         {
881             G4double tm=0.;                       872             G4double tm=0.;
882             if(PDG==2212) tm=PCSmanager->GetHM    873             if(PDG==2212) tm=PCSmanager->GetHMaxT();
883             else          tm=NCSmanager->GetHM    874             else          tm=NCSmanager->GetHMaxT();
884             G4cerr<<"G4QuasiFreeRatio::Scat:*N    875             G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
885             return std::make_pair(G4LorentzVec    876             return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
886         }                                         877         }
887     }                                             878     }
888     G4LorentzVector reco4M=G4LorentzVector(0.,    879     G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
889     G4LorentzVector dir4M=tot4M-G4LorentzVecto    880     G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
890     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M    881     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
891     {                                             882     {
892         G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M    883         G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
893         //G4Exception("G4QFR::Scat:","009",Fat    884         //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
894         return std::make_pair(G4LorentzVector(    885         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
895     }                                             886     }
896     return std::make_pair(reco4M*megaelectronv    887     return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
897 } // End of Scatter                               888 } // End of Scatter
898                                                   889 
899 // scatter (pPDG,p4M) on a virtual nucleon (NP    890 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
900 // if(newN4M.e()==0.) - below threshold, XS=0,    891 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
901 // User should himself change the charge (PDG)    892 // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.)
902 //AR-Jul2020: No need to change this method in << 
903 //            charmed and bottom mesons, becau << 
904 std::pair<G4LorentzVector,G4LorentzVector> G4Q    893 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::ChExer(G4int NPDG,
905                                                   894                                                                     G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
906 {                                                 895 {
907     static const G4double mNeut= G4Neutron::Ne    896     static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
908     static const G4double mProt= G4Proton::Pro    897     static const G4double mProt= G4Proton::Proton()->GetPDGMass();
909     G4LorentzVector pr4M=p4M/megaelectronvolt;    898     G4LorentzVector pr4M=p4M/megaelectronvolt;          // Convert 4-momenta in MeV(keep p4M)
910     N4M/=megaelectronvolt;                        899     N4M/=megaelectronvolt;
911     G4LorentzVector tot4M=N4M+p4M;                900     G4LorentzVector tot4M=N4M+p4M;
912     G4int Z=0;                                    901     G4int Z=0;
913     G4int N=1;                                    902     G4int N=1;
914     G4int sPDG=0;                                 903     G4int sPDG=0;                                        // PDG code of the scattered hadron
915     G4double mS=0.;                               904     G4double mS=0.;                                      // proto of mass of scattered hadron
916     G4double mT=mProt;                            905     G4double mT=mProt;                                   // mass of the recoil nucleon
917     if(NPDG==2212)                                906     if(NPDG==2212)
918     {                                             907     {
919         mT=mNeut;                                 908         mT=mNeut;
920         Z=1;                                      909         Z=1;
921         N=0;                                      910         N=0;
922         if(pPDG==-211) sPDG=111;                  911         if(pPDG==-211) sPDG=111;                           // pi+    -> pi0
923         else if(pPDG==-321)                       912         else if(pPDG==-321)
924         {                                         913         {
925             sPDG=310;                             914             sPDG=310;                                        // K+     -> K0S
926             if(G4UniformRand()>.5) sPDG=130;      915             if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
927         }                                         916         }
928         else if(pPDG==-311||pPDG==311||pPDG==1    917         else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;  // K0     -> K+ (?)
929         else if(pPDG==3112) sPDG=3212;            918         else if(pPDG==3112) sPDG=3212;                     // Sigma- -> Sigma0
930         else if(pPDG==3212) sPDG=3222;            919         else if(pPDG==3212) sPDG=3222;                     // Sigma0 -> Sigma+
931         else if(pPDG==3312) sPDG=3322;            920         else if(pPDG==3312) sPDG=3322;                     // Xi-    -> Xi0
932     }                                             921     }
933     else if(NPDG==2112) // Default                922     else if(NPDG==2112) // Default
934     {                                             923     {
935         if(pPDG==211)  sPDG=111;                  924         if(pPDG==211)  sPDG=111;                           // pi+    -> pi0
936         else if(pPDG==321)                        925         else if(pPDG==321)
937         {                                         926         {
938             sPDG=310;                             927             sPDG=310;                                        // K+     -> K0S
939             if(G4UniformRand()>.5) sPDG=130;      928             if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
940         }                                         929         }
941         else if(pPDG==-311||pPDG==311||pPDG==1    930         else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0     -> K- (?)
942         else if(pPDG==3222) sPDG=3212;            931         else if(pPDG==3222) sPDG=3212;                     // Sigma+ -> Sigma0
943         else if(pPDG==3212) sPDG=3112;            932         else if(pPDG==3212) sPDG=3112;                     // Sigma0 -> Sigma-
944         else if(pPDG==3322) sPDG=3312;            933         else if(pPDG==3322) sPDG=3312;                     // Xi0    -> Xi-
945     }                                             934     }
946     else                                          935     else
947     {                                             936     {
948         G4cout<<"Error:G4QuasiElRatios::ChExer    937         G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
949         G4Exception("G4QuasiElRatios::ChExer:"    938         G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
950         //return std::make_pair(G4LorentzVecto    939         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
951     }                                             940     }
952     if(sPDG) mS=mNeut;                            941     if(sPDG) mS=mNeut;
953     else                                          942     else
954     {                                             943     {
955         G4cout<<"Error:G4QuasiElRatios::ChExer    944         G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
956         G4Exception("G4QuasiElRatios::ChExer:"    945         G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
957         //return std::make_pair(G4LorentzVecto    946         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
958     }                                             947     }
959     G4double mT2=mT*mT;                           948     G4double mT2=mT*mT;
960     G4double mS2=mS*mS;                           949     G4double mS2=mS*mS;
961     G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);      950     G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
962     G4double E2=E*E;                              951     G4double E2=E*E;
963     if(E<0. || E2<mS2)                            952     if(E<0. || E2<mS2)
964     {                                             953     {
965         return std::make_pair(G4LorentzVector(    954         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
966     }                                             955     }
967     G4double P=std::sqrt(E2-mS2);                 956     G4double P=std::sqrt(E2-mS2);                   // Momentum in pseudo laboratory system
968     // @@ Temporary NN t-dependence for all ha    957     // @@ Temporary NN t-dependence for all hadrons
969     G4int PDG=2212;                               958     G4int PDG=2212;                                                // *TMP* instead of pPDG
970     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG    959     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
971     if(!Z && N==1)                 // Change f    960     if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
972     {                                             961     {
973         Z=1;                                      962         Z=1;
974         N=0;                                      963         N=0;
975         if     (PDG==2212) PDG=2112;              964         if     (PDG==2212) PDG=2112;
976         else if(PDG==2112) PDG=2212;              965         else if(PDG==2112) PDG=2212;
977     }                                             966     }
978     G4double xSec=0.;                        /    967     G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
979     if(PDG==2212) xSec=PCSmanager->GetChipsCro    968     if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
980     else          xSec=NCSmanager->GetChipsCro    969     else          xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
981     // @@ check a possibility to separate p, n    970     // @@ check a possibility to separate p, n, or alpha (!)
982     if(xSec <= 0.) // The cross-section iz 0 -    971     if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
983     {                                             972     {
984         return std::make_pair(G4LorentzVector(    973         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
985     }                                             974     }
986     G4double mint=0.;                        /    975     G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
987     if(PDG==2212) mint=PCSmanager->GetExchange    976     if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
988     else          mint=NCSmanager->GetExchange    977     else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
989     G4double maxt=0.;                             978     G4double maxt=0.;                                    // Prototype of max possible -t
990     if(PDG==2212) maxt=PCSmanager->GetHMaxT();    979     if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
991     else          maxt=NCSmanager->GetHMaxT();    980     else          maxt=NCSmanager->GetHMaxT();           // max possible -t
992     G4double cost=1.-mint/maxt;                   981     G4double cost=1.-mint/maxt;                          // cos(theta) in CMS
993     if(cost>1. || cost<-1. || !(cost>-1. || co    982     if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
994     {                                             983     {
995         if     (cost>1.)  cost=1.;                984         if     (cost>1.)  cost=1.;
996         else if(cost<-1.) cost=-1.;               985         else if(cost<-1.) cost=-1.;
997         else                                      986         else
998         {                                         987         {
999             G4cerr<<"G4QuasiFreeRatio::ChExer:    988             G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
1000             return std::make_pair(G4LorentzVe    989             return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1001         }                                        990         }
1002     }                                            991     }
1003     G4LorentzVector reco4M=G4LorentzVector(0.    992     G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
1004     pr4M=G4LorentzVector(0.,0.,0.,mS);           993     pr4M=G4LorentzVector(0.,0.,0.,mS);                        // 4mom of the scattered hadron
1005     G4LorentzVector dir4M=tot4M-G4LorentzVect    994     G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
1006     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4    995     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
1007     {                                            996     {
1008         G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4    997         G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
1009         //G4Exception("G4QFR::ChExer:","009",    998         //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
1010         return std::make_pair(G4LorentzVector    999         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1011     }                                            1000     }
1012     return std::make_pair(reco4M*megaelectron    1001     return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
1013 } // End of ChExer                               1002 } // End of ChExer
1014                                                  1003 
1015 // Calculate ChEx/El ratio (p is in independe    1004 // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile)
1016 G4double G4QuasiElRatios::ChExElCoef(G4double    1005 G4double G4QuasiElRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) 
1017 {                                                1006 {
1018                                               << 
1019     p/=MeV;                                //    1007     p/=MeV;                                // Converted from independent units
1020     G4double A=Z+N;                              1008     G4double A=Z+N;
1021     if(A<1.5) return 0.;                         1009     if(A<1.5) return 0.;
1022     G4double Cex=0.;                          << 1010     G4double C=0.;
1023     if     (pPDG==2212) Cex=N/(A+Z);          << 1011     if     (pPDG==2212) C=N/(A+Z);
1024     else if(pPDG==2112) Cex=Z/(A+N);          << 1012     else if(pPDG==2112) C=Z/(A+N);
1025     else G4cout<<"*Warning*G4CohChrgExchange:    1013     else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
1026     Cex*=Cex;                         // Cohe << 1014     C*=C;                         // Coherent processes squares the amplitude
1027     // @@ This is true only for nucleons: oth    1015     // @@ This is true only for nucleons: other projectiles must be treated differently
1028     G4double sp=std::sqrt(p);                    1016     G4double sp=std::sqrt(p);
1029     G4double p2=p*p;                             1017     G4double p2=p*p;            
1030     G4double p4=p2*p2;                           1018     G4double p4=p2*p2;
1031     G4double dl1=G4Log(p)-5.;                 << 1019     G4double dl1=std::log(p)-5.;
1032     G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.    1020     G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1033     G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)    1021     G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; 
1034     G4double R=U/T;                              1022     G4double R=U/T;
1035     return Cex*R*R;                           << 1023     return C*R*R;
1036 }                                                1024 }
1037                                                  1025 
1038 // Decay of Hadron In2Particles f&s, f is in     1026 // Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir
1039 G4bool G4QuasiElRatios::RelDecayIn2(G4Lorentz    1027 G4bool G4QuasiElRatios::RelDecayIn2(G4LorentzVector& theMomentum, G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
1040                                      G4Lorent    1028                                      G4LorentzVector& dir, G4double maxCost, G4double minCost)
1041 {                                                1029 {
1042     G4double fM2 = f4Mom.m2();                   1030     G4double fM2 = f4Mom.m2();
1043     G4double fM  = std::sqrt(fM2);               1031     G4double fM  = std::sqrt(fM2);              // Mass of the 1st Hadron
1044     G4double sM2 = s4Mom.m2();                   1032     G4double sM2 = s4Mom.m2();
1045     G4double sM  = std::sqrt(sM2);               1033     G4double sM  = std::sqrt(sM2);              // Mass of the 2nd Hadron
1046     G4double iM2 = theMomentum.m2();             1034     G4double iM2 = theMomentum.m2();
1047     G4double iM  = std::sqrt(iM2);               1035     G4double iM  = std::sqrt(iM2);              // Mass of the decaying hadron
1048     G4double vP  = theMomentum.rho();      //    1036     G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
1049     G4double dE  = theMomentum.e();        //    1037     G4double dE  = theMomentum.e();        // Energy of the decaying hadron
1050     if(dE<vP)                                    1038     if(dE<vP)
1051     {                                            1039     {
1052         G4cerr<<"***G4QHad::RelDecIn2: Tachio    1040         G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
1053         G4double accuracy=.000001*vP;            1041         G4double accuracy=.000001*vP;
1054         G4double emodif=std::fabs(dE-vP);        1042         G4double emodif=std::fabs(dE-vP);
1055         //if(emodif<accuracy)                    1043         //if(emodif<accuracy)
1056         //{                                      1044         //{
1057         G4cerr<<"G4QHadron::RelDecIn2: *Boost    1045         G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
1058         theMomentum.setE(vP+emodif+.01*accura    1046         theMomentum.setE(vP+emodif+.01*accuracy);
1059         //}                                      1047         //}
1060     }                                            1048     }
1061     G4ThreeVector ltb = theMomentum.boostVect    1049     G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
1062     G4ThreeVector ltf = -ltb;              //    1050     G4ThreeVector ltf = -ltb;              // Boost vector for forward Lorentz Trans.
1063     G4LorentzVector cdir = dir;            //    1051     G4LorentzVector cdir = dir;            // A copy to make a transformation to CMS
1064     cdir.boost(ltf);                       //    1052     cdir.boost(ltf);                       // Direction transpormed to CMS of the Momentum
1065     G4ThreeVector vdir = cdir.vect();      //    1053     G4ThreeVector vdir = cdir.vect();      // 3-Vector of the direction-particle
1066     G4ThreeVector vx(0.,0.,1.);            //    1054     G4ThreeVector vx(0.,0.,1.);            // Ort in the direction of the reference particle
1067     G4ThreeVector vy(0.,1.,0.);            //    1055     G4ThreeVector vy(0.,1.,0.);            // First ort orthogonal to the direction
1068     G4ThreeVector vz(1.,0.,0.);            //    1056     G4ThreeVector vz(1.,0.,0.);            // Second ort orthoganal to the direction
1069     if(vdir.mag2() > 0.)                   //    1057     if(vdir.mag2() > 0.)                   // the refference particle isn't at rest in CMS
1070     {                                            1058     {
1071         vx = vdir.unit();                        1059         vx = vdir.unit();                    // Ort in the direction of the reference particle
1072         G4ThreeVector vv= vx.orthogonal();       1060         G4ThreeVector vv= vx.orthogonal();   // Not normed orthogonal vector (!)
1073         vy = vv.unit();                          1061         vy = vv.unit();                      // First ort orthogonal to the direction
1074         vz = vx.cross(vy);                       1062         vz = vx.cross(vy);                   // Second ort orthoganal to the direction
1075     }                                            1063     }
1076     if(maxCost> 1.) maxCost= 1.;                 1064     if(maxCost> 1.) maxCost= 1.;
1077     if(minCost<-1.) minCost=-1.;                 1065     if(minCost<-1.) minCost=-1.;
1078     if(maxCost<-1.) maxCost=-1.;                 1066     if(maxCost<-1.) maxCost=-1.;
1079     if(minCost> 1.) minCost= 1.;                 1067     if(minCost> 1.) minCost= 1.;
1080     if(minCost> maxCost) minCost=maxCost;        1068     if(minCost> maxCost) minCost=maxCost;
1081     if(std::fabs(iM-fM-sM)<.00000001)            1069     if(std::fabs(iM-fM-sM)<.00000001)
1082     {                                            1070     {
1083         G4double fR=fM/iM;                       1071         G4double fR=fM/iM;
1084         G4double sR=sM/iM;                       1072         G4double sR=sM/iM;
1085         f4Mom=fR*theMomentum;                    1073         f4Mom=fR*theMomentum;
1086         s4Mom=sR*theMomentum;                    1074         s4Mom=sR*theMomentum;
1087         return true;                             1075         return true;
1088     }                                            1076     }
1089     else if (iM+.001<fM+sM || iM==0.)            1077     else if (iM+.001<fM+sM || iM==0.)
1090     {//@@ Later on make a quark content check    1078     {//@@ Later on make a quark content check for the decay
1091         G4cerr<<"***G4QH::RelDecIn2: fM="<<fM    1079         G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
1092         return false;                            1080         return false;
1093     }                                            1081     }
1094     G4double d2 = iM2-fM2-sM2;                   1082     G4double d2 = iM2-fM2-sM2;
1095     G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;        1083     G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
1096     if(p2<0.)                                    1084     if(p2<0.)
1097     {                                            1085     {
1098         p2=0.;                                   1086         p2=0.;
1099     }                                            1087     }
1100     G4double p  = std::sqrt(p2);                 1088     G4double p  = std::sqrt(p2);
1101     G4double ct = maxCost;                       1089     G4double ct = maxCost;
1102     if(maxCost>minCost)                          1090     if(maxCost>minCost)
1103     {                                            1091     {
1104         G4double dcost=maxCost-minCost;          1092         G4double dcost=maxCost-minCost;
1105         ct = minCost+dcost*G4UniformRand();      1093         ct = minCost+dcost*G4UniformRand();
1106     }                                            1094     }
1107     G4double phi= twopi*G4UniformRand();  //     1095     G4double phi= twopi*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
1108     G4double pss=0.;                          << 1096     G4double ps=0.;
1109     if(std::fabs(ct)<1.) pss = p * std::sqrt( << 1097     if(std::fabs(ct)<1.) ps = p * std::sqrt(1.-ct*ct);
1110     else                                         1098     else
1111     {                                            1099     {
1112         if(ct>1.) ct=1.;                         1100         if(ct>1.) ct=1.;
1113         if(ct<-1.) ct=-1.;                       1101         if(ct<-1.) ct=-1.;
1114     }                                            1102     }
1115     G4ThreeVector pVect=(pss*std::sin(phi))*v << 1103     G4ThreeVector pVect=(ps*std::sin(phi))*vz+(ps*std::cos(phi))*vy+p*ct*vx;
1116                                                  1104     
1117     f4Mom.setVect(pVect);                        1105     f4Mom.setVect(pVect);
1118     f4Mom.setE(std::sqrt(fM2+p2));               1106     f4Mom.setE(std::sqrt(fM2+p2));
1119     s4Mom.setVect((-1)*pVect);                   1107     s4Mom.setVect((-1)*pVect);
1120     s4Mom.setE(std::sqrt(sM2+p2));               1108     s4Mom.setE(std::sqrt(sM2+p2));
1121                                                  1109     
1122     if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G    1110     if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
1123         <<f4Mom.e()-f4Mom.rho()<<G4endl;         1111         <<f4Mom.e()-f4Mom.rho()<<G4endl;
1124     f4Mom.boost(ltb);                            1112     f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
1125     if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G    1113     if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
1126         <<s4Mom.e()-s4Mom.rho()<<G4endl;         1114         <<s4Mom.e()-s4Mom.rho()<<G4endl;
1127     s4Mom.boost(ltb);                            1115     s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
1128     return true;                                 1116     return true;
1129 } // End of "RelDecayIn2"                        1117 } // End of "RelDecayIn2"
1130                                                  1118 
1131                                                  1119 
1132                                                  1120 
1133                                                  1121 
1134                                                  1122 
1135                                                  1123 
1136                                                  1124