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