Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QHyperonElasticCrossSection.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id$
28//
29//
30// G4 Physics class: G4QHyperonElasticCrossSection for pA elastic cross sections
31// Created: M.V. Kossov, CERN/ITEP(Moscow), 5-Feb-2010
32// The last update: M.V. Kossov, CERN/ITEP (Moscow) 5-Feb-2010
33//
34// -------------------------------------------------------------------------------
35// Short description: Interaction cross-sections for the G4QElastic process
36// -------------------------------------------------------------------------------
37
38//#define debug
39//#define isodebug
40//#define pdebug
41//#define ppdebug
42//#define tdebug
43//#define sdebug
44
46#include "G4SystemOfUnits.hh"
47
48// Initialization of the static parameters
49const G4int G4QHyperonElasticCrossSection::nPoints=128;//#ofPt in AMDB table(>anyPar)(D)
50const G4int G4QHyperonElasticCrossSection::nLast=nPoints-1;// theLastElement inTable (D)
51G4double G4QHyperonElasticCrossSection::lPMin=-8.; //Min tabulatedLogarithmMomentum(D)
52G4double G4QHyperonElasticCrossSection::lPMax= 8.; //Max tabulatedLogarithmMomentum(D)
53G4double G4QHyperonElasticCrossSection::dlnP=(lPMax-lPMin)/nLast;// LogStep inTable (D)
54G4bool G4QHyperonElasticCrossSection::onlyCS=true;//Flag toCalculOnlyCS(not Si/Bi)(L)
55G4double G4QHyperonElasticCrossSection::lastSIG=0.; //Last calculated cross section (L)
56G4double G4QHyperonElasticCrossSection::lastLP=-10.;//LastLog(mom_of IncidentHadron)(L)
57G4double G4QHyperonElasticCrossSection::lastTM=0.; //Last t_maximum (L)
58G4double G4QHyperonElasticCrossSection::theSS=0.; //TheLastSqSlope of 1st difr.Max(L)
59G4double G4QHyperonElasticCrossSection::theS1=0.; //TheLastMantissa of 1st difrMax(L)
60G4double G4QHyperonElasticCrossSection::theB1=0.; //TheLastSlope of 1st difructMax(L)
61G4double G4QHyperonElasticCrossSection::theS2=0.; //TheLastMantissa of 2nd difrMax(L)
62G4double G4QHyperonElasticCrossSection::theB2=0.; //TheLastSlope of 2nd difructMax(L)
63G4double G4QHyperonElasticCrossSection::theS3=0.; //TheLastMantissa of 3d difr.Max(L)
64G4double G4QHyperonElasticCrossSection::theB3=0.; //TheLastSlope of 3d difruct.Max(L)
65G4double G4QHyperonElasticCrossSection::theS4=0.; //TheLastMantissa of 4th difrMax(L)
66G4double G4QHyperonElasticCrossSection::theB4=0.; //TheLastSlope of 4th difructMax(L)
67G4int G4QHyperonElasticCrossSection::lastTZ=0; // Last atomic number of the target
68G4int G4QHyperonElasticCrossSection::lastTN=0; // Last # of neutrons in the target
69G4double G4QHyperonElasticCrossSection::lastPIN=0.; // Last initialized max momentum
70G4double* G4QHyperonElasticCrossSection::lastCST=0; // Elastic cross-section table
71G4double* G4QHyperonElasticCrossSection::lastPAR=0; // ParametersForFunctionCalculation
72G4double* G4QHyperonElasticCrossSection::lastSST=0; // E-dep ofSqardSlope of 1st difMax
73G4double* G4QHyperonElasticCrossSection::lastS1T=0; // E-dep of mantissa of 1st dif.Max
74G4double* G4QHyperonElasticCrossSection::lastB1T=0; // E-dep of the slope of 1st difMax
75G4double* G4QHyperonElasticCrossSection::lastS2T=0; // E-dep of mantissa of 2nd difrMax
76G4double* G4QHyperonElasticCrossSection::lastB2T=0; // E-dep of the slope of 2nd difMax
77G4double* G4QHyperonElasticCrossSection::lastS3T=0; // E-dep of mantissa of 3d difr.Max
78G4double* G4QHyperonElasticCrossSection::lastB3T=0; // E-dep of the slope of 3d difrMax
79G4double* G4QHyperonElasticCrossSection::lastS4T=0; // E-dep of mantissa of 4th difrMax
80G4double* G4QHyperonElasticCrossSection::lastB4T=0; // E-dep of the slope of 4th difMax
81G4int G4QHyperonElasticCrossSection::lastN=0; // The last N of calculated nucleus
82G4int G4QHyperonElasticCrossSection::lastZ=0; // The last Z of calculated nucleus
83G4double G4QHyperonElasticCrossSection::lastP=0.; // LastUsed inCrossSection Momentum
84G4double G4QHyperonElasticCrossSection::lastTH=0.; // Last threshold momentum
85G4double G4QHyperonElasticCrossSection::lastCS=0.; // Last value of the Cross Section
86G4int G4QHyperonElasticCrossSection::lastI=0; // The last position in the DAMDB
87
88std::vector<G4double*> G4QHyperonElasticCrossSection::PAR;// Vector ofParsForFunctCalcul
89std::vector<G4double*> G4QHyperonElasticCrossSection::CST;// Vector ofCrossSection table
90std::vector<G4double*> G4QHyperonElasticCrossSection::SST;// Vector ofThe1st SquardSlope
91std::vector<G4double*> G4QHyperonElasticCrossSection::S1T;// Vector of the 1st mantissa
92std::vector<G4double*> G4QHyperonElasticCrossSection::B1T;// Vector of the1st slope
93std::vector<G4double*> G4QHyperonElasticCrossSection::S2T;// Vector of the2nd mantissa
94std::vector<G4double*> G4QHyperonElasticCrossSection::B2T;// Vector of the2nd slope
95std::vector<G4double*> G4QHyperonElasticCrossSection::S3T;// Vector of the3d mantissa
96std::vector<G4double*> G4QHyperonElasticCrossSection::B3T;// Vector of the3d slope
97std::vector<G4double*> G4QHyperonElasticCrossSection::S4T;// Vector ofThe4th mantissa(g)
98std::vector<G4double*> G4QHyperonElasticCrossSection::B4T;// Vector ofThe4th slope(glor)
99
101{
102}
103
105{
106 std::vector<G4double*>::iterator pos;
107 for (pos=CST.begin(); pos<CST.end(); pos++)
108 { delete [] *pos; }
109 CST.clear();
110 for (pos=PAR.begin(); pos<PAR.end(); pos++)
111 { delete [] *pos; }
112 PAR.clear();
113 for (pos=SST.begin(); pos<SST.end(); pos++)
114 { delete [] *pos; }
115 SST.clear();
116 for (pos=S1T.begin(); pos<S1T.end(); pos++)
117 { delete [] *pos; }
118 S1T.clear();
119 for (pos=B1T.begin(); pos<B1T.end(); pos++)
120 { delete [] *pos; }
121 B1T.clear();
122 for (pos=S2T.begin(); pos<S2T.end(); pos++)
123 { delete [] *pos; }
124 S2T.clear();
125 for (pos=B2T.begin(); pos<B2T.end(); pos++)
126 { delete [] *pos; }
127 B2T.clear();
128 for (pos=S3T.begin(); pos<S3T.end(); pos++)
129 { delete [] *pos; }
130 S3T.clear();
131 for (pos=B3T.begin(); pos<B3T.end(); pos++)
132 { delete [] *pos; }
133 B3T.clear();
134 for (pos=S4T.begin(); pos<S4T.end(); pos++)
135 { delete [] *pos; }
136 S4T.clear();
137 for (pos=B4T.begin(); pos<B4T.end(); pos++)
138 { delete [] *pos; }
139 B4T.clear();
140}
141
142// Returns Pointer to the G4VQCrossSection class
144{
145 static G4QHyperonElasticCrossSection theCrossSection;//StaticBody ofTheQElCrossSection
146 return &theCrossSection;
147}
148
149// The main member function giving the collision cross section (P is in IU, CS is in mb)
150// Make pMom in independent units ! (Now it is MeV)
152 G4int tgZ, G4int tgN, G4int pPDG)
153{
154 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops)
155 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops)
156 static std::vector <G4double> colP; // Vector of last momenta for the reaction
157 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
158 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
159 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
160 G4double pEn=pMom;
161 onlyCS=fCS;
162#ifdef pdebug
163 G4cout<<"G4QHyperElCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="
164 <<tgN<<"("<<lastN<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="<<colN.size()<<G4endl;
165 //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
166#endif
167 if(pPDG==3222 || pPDG<3000 || pPDG>3334)
168 {
169 G4cout<<"*Warning*G4QHyperonElaCS::GetCS:**> Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
170 //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
171 return 0.; // projectile PDG=0 is a mistake (?!) @@
172 }
173 G4bool in=false; // By default the isotope must be found in the AMDB
174 lastP = 0.; // New momentum history (nothing to compare with)
175 lastN = tgN; // The last N of the calculated nucleus
176 lastZ = tgZ; // The last Z of the calculated nucleus
177 lastI = colN.size(); // Size of the Associative Memory DB in the heap
178 if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
179 { // The nucleus with projPDG is found in AMDB
180 if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
181 {
182 lastI=i;
183 lastTH =colTH[i]; // Last THreshold (A-dependent)
184#ifdef pdebug
185 G4cout<<"G4QHyElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
186 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
187#endif
188 if(pEn<=lastTH)
189 {
190#ifdef pdebug
191 G4cout<<"G4QHyElCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
192 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
193#endif
194 return 0.; // Energy is below the Threshold value
195 }
196 lastP =colP [i]; // Last Momentum (A-dependent)
197 lastCS =colCS[i]; // Last CrossSect (A-dependent)
198 // if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance)
199 if(lastP == pMom) // Do not recalculate
200 {
201#ifdef pdebug
202 G4cout<<"G4QHyperonElasticCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
203#endif
204 CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // Update param's only
205 return lastCS*millibarn; // Use theLastCS
206 }
207 in = true; // This is the case when the isotop is found in DB
208 // Momentum pMom is in IU ! @@ Units
209#ifdef pdebug
210 G4cout<<"G4QHyElCS::G:UpdateDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",i="<<i<<G4endl;
211#endif
212 lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // read & update
213#ifdef pdebug
214 G4cout<<"G4QHyperElCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
215 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
216#endif
217 if(lastCS<=0. && pEn>lastTH) // Correct the threshold
218 {
219#ifdef pdebug
220 G4cout<<"G4QHyperonElCS::GetCS: NewT="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
221#endif
222 lastTH=pEn;
223 }
224 break; // Go out of the LOOP with found lastI
225 }
226#ifdef pdebug
227 G4cout<<"---G4QHyperonElCrossSection::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N="<<colN[i]
228 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
229 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
230#endif
231 } // End of attampt to find the nucleus in DB
232 if(!in) // This nucleus has not been calculated previously
233 {
234#ifdef pdebug
235 G4cout<<"G4QHyElCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
236#endif
237 //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
238 lastCS=CalculateCrossSection(fCS,0,lastI,pPDG,lastZ,lastN,pMom);//calculate&create
239 if(lastCS<=0.)
240 {
241 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
242#ifdef pdebug
243 G4cout<<"G4QHyperonElasticCrossSection::GetCS:NewThres="<<lastTH<<",T="<<pEn<<G4endl;
244#endif
245 if(pEn>lastTH)
246 {
247#ifdef pdebug
248 G4cout<<"G4QHyperonElCS::GetCS:1st T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
249#endif
250 lastTH=pEn;
251 }
252 }
253#ifdef pdebug
254 G4cout<<"G4QHyElCS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
255 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
256#endif
257 colN.push_back(tgN);
258 colZ.push_back(tgZ);
259 colP.push_back(pMom);
260 colTH.push_back(lastTH);
261 colCS.push_back(lastCS);
262#ifdef pdebug
263 G4cout<<"G4QHyElCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
264 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
265#endif
266 return lastCS*millibarn;
267 } // End of creation of the new set of parameters
268 else
269 {
270#ifdef pdebug
271 G4cout<<"G4QHyperonElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl;
272#endif
273 colP[lastI]=pMom;
274 colCS[lastI]=lastCS;
275 }
276#ifdef pdebug
277 G4cout<<"G4QHyElCS::GetCSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
278 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
279 G4cout<<"G4QHyperonElasticCrossSection::GetCrosSec:***End***,onlyCS="<<onlyCS<<G4endl;
280#endif
281 return lastCS*millibarn;
282}
283
284// Calculation of total elastic cross section (p in IU, CS in mb) @@ Units (?)
285// F=0 - create AMDB, F=-1 - read&update AMDB, F=1 - update AMDB (sinchro with higher AMDB)
287 G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
288{
289 // *** Begin of Associative Memory DB for acceleration of the cross section calculations
290 static std::vector <G4double> PIN; // Vector of max initialized log(P) in the table
291 // *** End of Static Definitions (Associative Memory Data Base) ***
292 G4double pMom=pIU/GeV; // All calculations are in GeV
293 onlyCS=CS; // Flag to calculate only CS (not Si/Bi)
294#ifdef pdebug
295 G4cout<<"G4QHyperonElasticCS::CalcCS:->onlyCS="<<onlyCS<<",F="<<F<<",p="<<pIU<<G4endl;
296#endif
297 lastLP=std::log(pMom); // Make a logarithm of the momentum for calculation
298 if(F) // This isotope was found in AMDB =>RETRIEVE/UPDATE
299 {
300 if(F<0) // the AMDB must be loded
301 {
302 lastPIN = PIN[I]; // Max log(P) initialised for this table set
303 lastPAR = PAR[I]; // Pointer to the parameter set
304 lastCST = CST[I]; // Pointer to the total sross-section table
305 lastSST = SST[I]; // Pointer to the first squared slope
306 lastS1T = S1T[I]; // Pointer to the first mantissa
307 lastB1T = B1T[I]; // Pointer to the first slope
308 lastS2T = S2T[I]; // Pointer to the second mantissa
309 lastB2T = B2T[I]; // Pointer to the second slope
310 lastS3T = S3T[I]; // Pointer to the third mantissa
311 lastB3T = B3T[I]; // Pointer to the rhird slope
312 lastS4T = S4T[I]; // Pointer to the 4-th mantissa
313 lastB4T = B4T[I]; // Pointer to the 4-th slope
314#ifdef pdebug
315 G4cout<<"G4QHyperonElCS::CalcCS: DB's updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
316#endif
317 }
318#ifdef pdebug
319 G4cout<<"G4QHyperonElasticCS::CalcCS:*read*, LP="<<lastLP<<",PIN="<<lastPIN<<G4endl;
320#endif
321 if(lastLP>lastPIN && lastLP<lPMax)
322 {
323 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
324#ifdef pdebug
325 G4cout<<"G4QHyElCS::CalcCS:updated(I),LP="<<lastLP<<"<IN["<<I<<"]="<<lastPIN<<G4endl;
326#endif
327 PIN[I]=lastPIN; // Remember the new P-Limit of the tables
328 }
329 }
330 else // This isotope wasn't initialized => CREATE
331 {
332 lastPAR = new G4double[nPoints]; // Allocate memory for parameters of CS function
333 lastPAR[nLast]=0; // Initialization for VALGRIND
334 lastCST = new G4double[nPoints]; // Allocate memory for Tabulated CS function
335 lastSST = new G4double[nPoints]; // Allocate memory for Tabulated first sqaredSlope
336 lastS1T = new G4double[nPoints]; // Allocate memory for Tabulated first mantissa
337 lastB1T = new G4double[nPoints]; // Allocate memory for Tabulated first slope
338 lastS2T = new G4double[nPoints]; // Allocate memory for Tabulated second mantissa
339 lastB2T = new G4double[nPoints]; // Allocate memory for Tabulated second slope
340 lastS3T = new G4double[nPoints]; // Allocate memory for Tabulated third mantissa
341 lastB3T = new G4double[nPoints]; // Allocate memory for Tabulated third slope
342 lastS4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th mantissa
343 lastB4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th slope
344#ifdef pdebug
345 G4cout<<"G4QHyperonElasticCrosS::CalcCS:*ini*,lastLP="<<lastLP<<",min="<<lPMin<<G4endl;
346#endif
347 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
348#ifdef pdebug
349 G4cout<<"G4QHypElCS::CCS:i, Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",LP"<<lastPIN<<G4endl;
350#endif
351 PIN.push_back(lastPIN); // Fill parameters of CS function to AMDB
352 PAR.push_back(lastPAR); // Fill parameters of CS function to AMDB
353 CST.push_back(lastCST); // Fill Tabulated CS function to AMDB
354 SST.push_back(lastSST); // Fill Tabulated first sq.slope to AMDB
355 S1T.push_back(lastS1T); // Fill Tabulated first mantissa to AMDB
356 B1T.push_back(lastB1T); // Fill Tabulated first slope to AMDB
357 S2T.push_back(lastS2T); // Fill Tabulated second mantissa to AMDB
358 B2T.push_back(lastB2T); // Fill Tabulated second slope to AMDB
359 S3T.push_back(lastS3T); // Fill Tabulated third mantissa to AMDB
360 B3T.push_back(lastB3T); // Fill Tabulated third slope to AMDB
361 S4T.push_back(lastS4T); // Fill Tabulated 4-th mantissa to AMDB
362 B4T.push_back(lastB4T); // Fill Tabulated 4-th slope to AMDB
363 } // End of creation/update of the new set of parameters and tables
364 // =-----------= NOW Update (if necessary) and Calculate the Cross Section =-----------=
365#ifdef pdebug
366 G4cout<<"G4QHypElCS::CalCS:?update?,LP="<<lastLP<<",IN="<<lastPIN<<",ML="<<lPMax<<G4endl;
367#endif
368 if(lastLP>lastPIN && lastLP<lPMax)
369 {
370 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
371#ifdef pdebug
372 G4cout<<"G4QHyperionElCS::CalcCS:*updated(O)*, LP="<<lastLP<<" < IN="<<lastPIN<<G4endl;
373#endif
374 }
375#ifdef pdebug
376 G4cout<<"G4QHypElCS::CalcCS: lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
377#endif
378 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
379#ifdef pdebug
380 G4cout<<"G4QHyperElCrosSec::CalcCS:oCS="<<onlyCS<<",-t="<<lastTM<<", p="<<lastLP<<G4endl;
381#endif
382 if(lastLP>lPMin && lastLP<=lastPIN) // Linear fit is made using precalculated tables
383 {
384 if(lastLP==lastPIN)
385 {
386 G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
387 G4int blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
388 if(blast<0 || blast>=nLast)G4cout<<"G4QHyperElCS::CCS:b="<<blast<<","<<nLast<<G4endl;
389 lastSIG = lastCST[blast];
390 if(!onlyCS) // Skip the differential cross-section parameters
391 {
392 theSS = lastSST[blast];
393 theS1 = lastS1T[blast];
394 theB1 = lastB1T[blast];
395 theS2 = lastS2T[blast];
396 theB2 = lastB2T[blast];
397 theS3 = lastS3T[blast];
398 theB3 = lastB3T[blast];
399 theS4 = lastS4T[blast];
400 theB4 = lastB4T[blast];
401 }
402#ifdef pdebug
403 G4cout<<"G4QHyperonElasticCS::CalculateCS:(E) S1="<<theS1<<", B1="<<theB1<<G4endl;
404#endif
405 }
406 else
407 {
408 G4double shift=(lastLP-lPMin)/dlnP; // a shift from the beginning of the table
409 G4int blast=static_cast<int>(shift); // the lower bin number
410 if(blast<0) blast=0;
411 if(blast>=nLast) blast=nLast-1; // low edge of the last bin
412 shift-=blast; // step inside the unit bin
413 G4int lastL=blast+1; // the upper bin number
414 G4double SIGL=lastCST[blast]; // the basic value of the cross-section
415 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
416#ifdef pdebug
417 G4cout<<"G4QHyperonElasticCrossSection::CalcCrossSection:Sig="<<lastSIG<<",P="
418 <<pMom<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",onlyCS="<<onlyCS<<G4endl;
419#endif
420 if(!onlyCS) // Skip the differential cross-section parameters
421 {
422 G4double SSTL=lastSST[blast]; // the low bin of the first squared slope
423 theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
424 G4double S1TL=lastS1T[blast]; // the low bin of the first mantissa
425 theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
426 G4double B1TL=lastB1T[blast]; // the low bin of the first slope
427#ifdef pdebug
428 G4cout<<"G4QHyperonElasticCS::CalcCrossSection:b="<<blast<<",ls="<<lastL<<",SL="
429 <<S1TL<<",SU="<<lastS1T[lastL]<<",BL="<<B1TL<<",BU="<<lastB1T[lastL]<<G4endl;
430#endif
431 theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
432 G4double S2TL=lastS2T[blast]; // the low bin of the second mantissa
433 theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
434 G4double B2TL=lastB2T[blast]; // the low bin of the second slope
435 theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
436 G4double S3TL=lastS3T[blast]; // the low bin of the third mantissa
437 theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
438#ifdef pdebug
439 G4cout<<"G4QHyperonElasticCrossSection::CCS:s3l="<<S3TL<<",sh3="<<shift<<",s3h="
440 <<lastS3T[lastL]<<",b="<<blast<<",l="<<lastL<<G4endl;
441#endif
442 G4double B3TL=lastB3T[blast]; // the low bin of the third slope
443 theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
444 G4double S4TL=lastS4T[blast]; // the low bin of the 4-th mantissa
445 theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
446#ifdef pdebug
447 G4cout<<"G4QHyperonElasticCrossSection::CCS:s4l="<<S4TL<<",sh4="<<shift<<",s4h="
448 <<lastS4T[lastL]<<",b="<<blast<<",l="<<lastL<<G4endl;
449#endif
450 G4double B4TL=lastB4T[blast]; // the low bin of the 4-th slope
451 theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
452 }
453#ifdef pdebug
454 G4cout<<"G4QHyperonElasticCS::CalculateCS:(I) S1="<<theS1<<", B1="<<theB1<<G4endl;
455#endif
456 }
457 }
458 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
459 if(lastSIG<0.) lastSIG = 0.; // @@ a Warning print can be added
460#ifdef pdebug
461 G4cout<<"G4QHyperonElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
462#endif
463 return lastSIG;
464}
465
466// It has parameter sets for all tZ/tN/PDG, using them the tables can be created/updated
467G4double G4QHyperonElasticCrossSection::GetPTables(G4double LP, G4double ILP, G4int PDG,
468 G4int tgZ, G4int tgN)
469{
470 // @@ At present all nA==pA ---------> Each neucleus can have not more than 51 parameters
471 static const G4double pwd=2727;
472 const G4int n_hypel=33; // #of parameters for pp-elastic (<nPoints=128)
473 // -0- -1- -2- -3- -4- -5- -6--7--8--9--10--11--12-13--14-
474 G4double hyp_el[n_hypel]={1.,.002,.12,.0557,3.5,6.72,99.,2.,3.,5.,74.,3.,3.4,.2,.17,
475 .001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,5.e-5,
476 1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
477 // -15--16- -17- -18- -19- -20- -21- -22- -23- -24- -25-
478 // -26- -27- -28- -29- -30- -31- -32-
479 if(PDG!=3222 && PDG>3000 && PDG<3335)
480 {
481 // -- Total pp elastic cross section cs & s1/b1 (main), s2/b2 (tail1), s3/b3 (tail2) --
482 //p2=p*p;p3=p2*p;sp=sqrt(p);p2s=p2*sp;lp=log(p);dl1=lp-(3.=par(3));p4=p2*p2; p=|3-mom|
483 //CS=2.865/p2s/(1+.0022/p2s)+(18.9+.6461*dl1*dl1+9./p)/(1.+.425*lp)/(1.+.4276/p4);
484 // par(0) par(7) par(1) par(2) par(4) par(5) par(6)
485 //dl2=lp-5., s1=(74.+3.*dl2*dl2)/(1+3.4/p4/p)+(.2/p2+17.*p)/(p4+.001*sp),
486 // par(8) par(9) par(10) par(11) par(12)par(13) par(14)
487 // b1=8.*p**.055/(1.+3.64/p3); s2=5.e-5+4000./(p4+1500.*p); b2=.46+1.2e6/(p4+3.5e6/sp);
488 // par(15) par(16) par(17) par(18) par(19) par(20) par(21) par(22) par(23)
489 // s3=5.e-5+1.e10/(p4*p4+8.5e8*p2+1.e10); b3=1.1+3.4e6/(p4+6.8e6); ss=0.
490 // par(24) par(25) par(26) par(27) par(28) par(29) par(30) par(31)
491 //
492 if(lastPAR[nLast]!=pwd) // A unique flag to avoid the repeatable definition
493 {
494 if ( tgZ == 1 && tgN == 0 )
495 {
496 for (G4int ip=0; ip<n_hypel; ip++) lastPAR[ip]=hyp_el[ip]; // Hyperon+P
497 }
498 else
499 {
500 G4double a=tgZ+tgN;
501 G4double sa=std::sqrt(a);
502 G4double ssa=std::sqrt(sa);
503 G4double asa=a*sa;
504 G4double a2=a*a;
505 G4double a3=a2*a;
506 G4double a4=a3*a;
507 G4double a5=a4*a;
508 G4double a6=a4*a2;
509 G4double a7=a6*a;
510 G4double a8=a7*a;
511 G4double a9=a8*a;
512 G4double a10=a5*a5;
513 G4double a12=a6*a6;
514 G4double a14=a7*a7;
515 G4double a16=a8*a8;
516 G4double a17=a16*a;
517 //G4double a20=a16*a4;
518 G4double a32=a16*a16;
519 // Reaction cross-section parameters (pel=peh_fit.f)
520 lastPAR[0]=4./(1.+22/asa); // p1
521 lastPAR[1]=2.36*asa/(1.+a*.055/ssa); // p2
522 lastPAR[2]=(1.+.00007*a3/ssa)/(1.+.0026*a2); // p3
523 lastPAR[3]=1.76*a/ssa+.00003*a3; // p4
524 lastPAR[4]=(.03+200./a3)/(1.+1.E5/a3/sa); // p5
525 lastPAR[5]=5.; // p6
526 lastPAR[6]=0.; // p7 not used
527 lastPAR[7]=0.; // p8 not used
528 lastPAR[8]=0.; // p9 not used
529 // @@ the differential cross-section is parameterized separately for A>6 & A<7
530 if(a<6.5)
531 {
532 G4double a28=a16*a12;
533 // The main pre-exponent (pel_sg)
534 lastPAR[ 9]=4000*a; // p1
535 lastPAR[10]=1.2e7*a8+380*a17; // p2
536 lastPAR[11]=.7/(1.+4.e-12*a16); // p3
537 lastPAR[12]=2.5/a8/(a4+1.e-16*a32); // p4
538 lastPAR[13]=.28*a; // p5
539 lastPAR[14]=1.2*a2+2.3; // p6
540 lastPAR[15]=3.8/a; // p7
541 // The main slope (pel_sl)
542 lastPAR[16]=.01/(1.+.0024*a5); // p1
543 lastPAR[17]=.2*a; // p2
544 lastPAR[18]=9.e-7/(1.+.035*a5); // p3
545 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a); // p4
546 // The main quadratic (pel_sh)
547 lastPAR[20]=2.25*a3; // p1
548 lastPAR[21]=18.; // p2
549 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7); // p3
550 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a); // p4
551 // The 1st max pre-exponent (pel_qq)
552 lastPAR[24]=1.e5/(a8+2.5e12/a16); // p1
553 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28); // p2
554 lastPAR[26]=.0006*a3; // p3
555 // The 1st max slope (pel_qs)
556 lastPAR[27]=10.+4.e-8*a12*a; // p1
557 lastPAR[28]=.114; // p2
558 lastPAR[29]=.003; // p3
559 lastPAR[30]=2.e-23; // p4
560 // The effective pre-exponent (pel_ss)
561 lastPAR[31]=1./(1.+.0001*a8); // p1
562 lastPAR[32]=1.5e-4/(1.+5.e-6*a12); // p2
563 lastPAR[33]=.03; // p3
564 // The effective slope (pel_sb)
565 lastPAR[34]=a/2; // p1
566 lastPAR[35]=2.e-7*a4; // p2
567 lastPAR[36]=4.; // p3
568 lastPAR[37]=64./a3; // p4
569 // The gloria pre-exponent (pel_us)
570 lastPAR[38]=1.e8*std::exp(.32*asa); // p1
571 lastPAR[39]=20.*std::exp(.45*asa); // p2
572 lastPAR[40]=7.e3+2.4e6/a5; // p3
573 lastPAR[41]=2.5e5*std::exp(.085*a3); // p4
574 lastPAR[42]=2.5*a; // p5
575 // The gloria slope (pel_ub)
576 lastPAR[43]=920.+.03*a8*a3; // p1
577 lastPAR[44]=93.+.0023*a12; // p2
578#ifdef pdebug
579 G4cout<<"G4QHyElCS::CalcCS:la "<<lastPAR[38]<<", "<<lastPAR[39]<<", "<<lastPAR[40]
580 <<", "<<lastPAR[42]<<", "<<lastPAR[43]<<", "<<lastPAR[44]<<G4endl;
581#endif
582 }
583 else
584 {
585 G4double p1a10=2.2e-28*a10;
586 G4double r4a16=6.e14/a16;
587 G4double s4a16=r4a16*r4a16;
588 // a24
589 // a36
590 // The main pre-exponent (peh_sg)
591 lastPAR[ 9]=4.5*std::pow(a,1.15); // p1
592 lastPAR[10]=.06*std::pow(a,.6); // p2
593 lastPAR[11]=.6*a/(1.+2.e15/a16); // p3
594 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32); // p4
595 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5); // p5
596 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12); // p6
597 // The main slope (peh_sl)
598 lastPAR[15]=400./a12+2.e-22*a9; // p1
599 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14); // p2
600 lastPAR[17]=1000./a2+9.5*sa*ssa; // p3
601 lastPAR[18]=4.e-6*a*asa+1.e11/a16; // p4
602 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16); // p5
603 lastPAR[20]=9.+100./a; // p6
604 // The main quadratic (peh_sh)
605 lastPAR[21]=.002*a3+3.e7/a6; // p1
606 lastPAR[22]=7.e-15*a4*asa; // p2
607 lastPAR[23]=9000./a4; // p3
608 // The 1st max pre-exponent (peh_qq)
609 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4); // p1
610 lastPAR[25]=1.e-5*a2+2.e14/a16; // p2
611 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12); // p3
612 lastPAR[27]=.016*asa/(1.+5.e16/a16); // p4
613 // The 1st max slope (peh_qs)
614 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14)); // p1
615 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11); // p2
616 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8); // p3
617 lastPAR[31]=100./asa; // p4
618 // The 2nd max pre-exponent (peh_ss)
619 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4); // p1
620 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8); // p2
621 lastPAR[34]=1.3+3.e5/a4; // p3
622 lastPAR[35]=500./(a2+50.)+3; // p4
623 lastPAR[36]=1.e-9/a+s4a16*s4a16; // p5
624 // The 2nd max slope (peh_sb)
625 lastPAR[37]=.4*asa+3.e-9*a6; // p1
626 lastPAR[38]=.0005*a5; // p2
627 lastPAR[39]=.002*a5; // p3
628 lastPAR[40]=10.; // p4
629 // The effective pre-exponent (peh_us)
630 lastPAR[41]=.05+.005*a; // p1
631 lastPAR[42]=7.e-8/sa; // p2
632 lastPAR[43]=.8*sa; // p3
633 lastPAR[44]=.02*sa; // p4
634 lastPAR[45]=1.e8/a3; // p5
635 lastPAR[46]=3.e32/(a32+1.e32); // p6
636 // The effective slope (peh_ub)
637 lastPAR[47]=24.; // p1
638 lastPAR[48]=20./sa; // p2
639 lastPAR[49]=7.e3*a/(sa+1.); // p3
640 lastPAR[50]=900.*sa/(1.+500./a3); // p4
641#ifdef pdebug
642 G4cout<<"G4QHyElCS::CalcCS:ha "<<lastPAR[41]<<", "<<lastPAR[42]<<", "<<lastPAR[43]
643 <<", "<<lastPAR[44]<<", "<<lastPAR[45]<<", "<<lastPAR[46]<<G4endl;
644#endif
645 }
646 // Parameter for lowEnergyNeutrons
647 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
648 }
649 lastPAR[nLast]=pwd;
650 // and initialize the zero element of the table
651 G4double lp=lPMin; // ln(momentum)
652 G4bool memCS=onlyCS; // ??
653 onlyCS=false;
654 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables
655 onlyCS=memCS;
656 lastSST[0]=theSS;
657 lastS1T[0]=theS1;
658 lastB1T[0]=theB1;
659 lastS2T[0]=theS2;
660 lastB2T[0]=theB2;
661 lastS3T[0]=theS3;
662 lastB3T[0]=theB3;
663 lastS4T[0]=theS4;
664 lastB4T[0]=theB4;
665#ifdef pdebug
666 G4cout<<"G4QHyperonElasticCrossSection::GetPTables:ip=0(init), lp="<<lp<<",S1="
667 <<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB3<<",S3="<<theS3
668 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
669#endif
670 }
671 if(LP>ILP)
672 {
673 G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1; // already inited till this
674 if(ini<0) ini=0;
675 if(ini<nPoints)
676 {
677 G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1; // final bin of initialization
678 if(fin>=nPoints) fin=nLast; // Limit of the tabular initialization
679 if(fin>=ini)
680 {
681 G4double lp=0.;
682 for(G4int ip=ini; ip<=fin; ip++) // Calculate tabular CS,S1,B1,S2,B2,S3,B3
683 {
684 lp=lPMin+ip*dlnP; // ln(momentum)
685 G4bool memCS=onlyCS;
686 onlyCS=false;
687 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables (ret CS)
688 onlyCS=memCS;
689 lastSST[ip]=theSS;
690 lastS1T[ip]=theS1;
691 lastB1T[ip]=theB1;
692 lastS2T[ip]=theS2;
693 lastB2T[ip]=theB2;
694 lastS3T[ip]=theS3;
695 lastB3T[ip]=theB3;
696 lastS4T[ip]=theS4;
697 lastB4T[ip]=theB4;
698#ifdef pdebug
699 G4cout<<"G4QHyperonElasticCrossSection::GetPTables:ip="<<ip<<",lp="<<lp
700 <<",S1="<<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="
701 <<theS3<<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
702#endif
703 }
704 return lp;
705 }
706 else G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetPTables: PDG="<<PDG
707 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
708 <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
709 }
710 else G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetPTables: PDG="<<PDG
711 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
712 <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
713 }
714#ifdef pdebug
715 else G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetPTab:PDG="<<PDG<<",Z="<<tgZ
716 <<", N="<<tgN<<", LP="<<LP<<" <= ILP="<<ILP<<" nothing is done!"<<G4endl;
717#endif
718 } else {
719 // G4cout<<"*Error*G4QHyperonElasticCrossSection::GetPTables: PDG="<<PDG<<", Z="<<tgZ
720 // <<", N="<<tgN<<", while it is defined only for Hyperons"<<G4endl;
721 // throw G4QException("G4QHyperonElasticCrossSection::GetPTables:onlyaBA implemented");
723 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
724 << ", while it is defined only for Hyperons" << G4endl;
725 G4Exception("G4QHyperonElasticCrossSection::GetPTables()", "HAD_CHPS_0000",
726 FatalException, ed);
727 }
728 return ILP;
729}
730
731// Returns Q2=-t in independent units (MeV^2) (all internal calculations are in GeV)
733{
734 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
735 static const G4double third=1./3.;
736 static const G4double fifth=1./5.;
737 static const G4double sevth=1./7.;
738#ifdef tdebug
739 G4cout<<"G4QHyperElCS::GetExcT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
740#endif
741 if(PDG==3222 || PDG<3000 || PDG>3334)G4cout<<"*Warning*G4QHyElCS::GET:PDG="<<PDG<<G4endl;
742 if(onlyCS)G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExchanT: onlyCS=1"<<G4endl;
743 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
744 G4double q2=0.;
745 if(tgZ==1 && tgN==0) // ===> p+p=p+p
746 {
747#ifdef tdebug
748 G4cout<<"G4QHyperElCS::GetExchangeT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",S2="
749 <<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",GeV2="<<GeVSQ<<G4endl;
750#endif
751 G4double E1=lastTM*theB1;
752 G4double R1=(1.-std::exp(-E1));
753#ifdef tdebug
754 G4double ts1=-std::log(1.-R1)/theB1;
755 G4double ds1=std::fabs(ts1-lastTM)/lastTM;
756 if(ds1>.0001)
757 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:1p "<<ts1<<"#"<<lastTM
758 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
759#endif
760 G4double E2=lastTM*theB2;
761 G4double R2=(1.-std::exp(-E2*E2*E2));
762#ifdef tdebug
763 G4double ts2=std::pow(-std::log(1.-R2),.333333333)/theB2;
764 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
765 if(ds2>.0001)
766 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:2p "<<ts2<<"#"<<lastTM
767 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
768#endif
769 G4double E3=lastTM*theB3;
770 G4double R3=(1.-std::exp(-E3));
771#ifdef tdebug
772 G4double ts3=-std::log(1.-R3)/theB3;
773 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
774 if(ds3>.0001)
775 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:3p "<<ts3<<"#"<<lastTM
776 <<",d="<<ds3<<",R3="<<R1<<",E3="<<E3<<G4endl;
777#endif
778 G4double I1=R1*theS1/theB1;
779 G4double I2=R2*theS2;
780 G4double I3=R3*theS3;
781 G4double I12=I1+I2;
782 G4double rand=(I12+I3)*G4UniformRand();
783 if (rand<I1 )
784 {
785 G4double ran=R1*G4UniformRand();
786 if(ran>1.) ran=1.;
787 q2=-std::log(1.-ran)/theB1;
788 }
789 else if(rand<I12)
790 {
791 G4double ran=R2*G4UniformRand();
792 if(ran>1.) ran=1.;
793 q2=-std::log(1.-ran);
794 if(q2<0.) q2=0.;
795 q2=std::pow(q2,third)/theB2;
796 }
797 else
798 {
799 G4double ran=R3*G4UniformRand();
800 if(ran>1.) ran=1.;
801 q2=-std::log(1.-ran)/theB3;
802 }
803 }
804 else
805 {
806 G4double a=tgZ+tgN;
807#ifdef tdebug
808 G4cout<<"G4QHyperonElasticCrossSection::GetExT:a="<<a<<",t="<<lastTM<<",S1="<<theS1
809 <<",B1="<<theB1<<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3
810 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
811#endif
812 G4double E1=lastTM*(theB1+lastTM*theSS);
813 G4double R1=(1.-std::exp(-E1));
814 G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
815#ifdef tdebug
816 G4double ts1=-std::log(1.-R1)/theB1;
817 if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss;
818 G4double ds1=(ts1-lastTM)/lastTM;
819 if(ds1>.0001)
820 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:1a "<<ts1<<"#"<<lastTM
821 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
822#endif
823 G4double tm2=lastTM*lastTM;
824 G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st)
825 if(a>6.5)E2*=tm2; // for heavy nuclei
826 G4double R2=(1.-std::exp(-E2));
827#ifdef tdebug
828 G4double ts2=-std::log(1.-R2)/theB2;
829 if(a<6.5)ts2=std::pow(ts2,third);
830 else ts2=std::pow(ts2,fifth);
831 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
832 if(ds2>.0001)
833 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:2a "<<ts2<<"#"<<lastTM
834 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
835#endif
836 G4double E3=lastTM*theB3;
837 if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA
838 G4double R3=(1.-std::exp(-E3));
839#ifdef tdebug
840 G4double ts3=-std::log(1.-R3)/theB3;
841 if(a>6.5)ts3=std::pow(ts3,sevth);
842 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
843 if(ds3>.0001)
844 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:3a "<<ts3<<"#"<<lastTM
845 <<",d="<<ds3<<",R3="<<R3<<",E3="<<E3<<G4endl;
846#endif
847 G4double E4=lastTM*theB4;
848 G4double R4=(1.-std::exp(-E4));
849#ifdef tdebug
850 G4double ts4=-std::log(1.-R4)/theB4;
851 G4double ds4=std::fabs(ts4-lastTM)/lastTM;
852 if(ds4>.0001)
853 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:4a "<<ts4<<"#"<<lastTM
854 <<",d="<<ds4<<",R4="<<R4<<",E4="<<E4<<G4endl;
855#endif
856 G4double I1=R1*theS1;
857 G4double I2=R2*theS2;
858 G4double I3=R3*theS3;
859 G4double I4=R4*theS4;
860 G4double I12=I1+I2;
861 G4double I13=I12+I3;
862 G4double rand=(I13+I4)*G4UniformRand();
863#ifdef tdebug
864 G4cout<<"G4QHyElCS::GExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl;
865#endif
866 if(rand<I1)
867 {
868 G4double ran=R1*G4UniformRand();
869 if(ran>1.) ran=1.;
870 q2=-std::log(1.-ran)/theB1;
871 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
872#ifdef tdebug
873 G4cout<<"G4QHyElCS::GExT:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl;
874#endif
875 }
876 else if(rand<I12)
877 {
878 G4double ran=R2*G4UniformRand();
879 if(ran>1.) ran=1.;
880 q2=-std::log(1.-ran)/theB2;
881 if(q2<0.) q2=0.;
882 if(a<6.5) q2=std::pow(q2,third);
883 else q2=std::pow(q2,fifth);
884#ifdef tdebug
885 G4cout<<"G4QHyElCS::GetExT: Q2="<<q2<<",r2="<<R2<<",b2="<<theB2<<",t2="<<ts2<<G4endl;
886#endif
887 }
888 else if(rand<I13)
889 {
890 G4double ran=R3*G4UniformRand();
891 if(ran>1.) ran=1.;
892 q2=-std::log(1.-ran)/theB3;
893 if(q2<0.) q2=0.;
894 if(a>6.5) q2=std::pow(q2,sevth);
895#ifdef tdebug
896 G4cout<<"G4QHyElCS::GetExT: Q2="<<q2<<",r3="<<R2<<",b3="<<theB2<<",t3="<<ts2<<G4endl;
897#endif
898 }
899 else
900 {
901 G4double ran=R4*G4UniformRand();
902 if(ran>1.) ran=1.;
903 q2=-std::log(1.-ran)/theB4;
904 if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0)
905#ifdef tdebug
906 G4cout<<"G4QHyElCS::GExT:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl;
907#endif
908 }
909 }
910 if(q2<0.) q2=0.;
911 if(!(q2>=-1.||q2<=1.))G4cout<<"*NAN*G4QHyElasticCrossSect::GetExchangeT:-t="<<q2<<G4endl;
912 if(q2>lastTM)
913 {
914#ifdef tdebug
915 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GExT: -t="<<q2<<" > "<<lastTM<<G4endl;
916#endif
917 q2=lastTM;
918 }
919 return q2*GeVSQ;
920}
921
922// Returns B in independent units (MeV^-2) (all internal calculations are in GeV) see ExT
924{
925 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
926#ifdef tdebug
927 G4cout<<"G4QHyperElCS::GetSlope:"<<onlyCS<<", Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
928#endif
929 if(onlyCS)G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetSlope: onlCS=true"<<G4endl;
930 if(lastLP<-4.3) return 0.; // S-wave for p<14 MeV/c (kinE<.1MeV)
931 if(PDG==3222 || PDG<3000 || PDG>3334)
932 {
933 // G4cout<<"*Error*G4QHyperonElasticCrossSection::GetSlope: PDG="<<PDG<<", Z="<<tgZ
934 // <<", N="<<tgN<<", while it is defined only for Hyperons"<<G4endl;
935 // throw G4QException("G4QHyperonElasticCrossSection::GetSlope: HypA are implemented");
937 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
938 << ", while it is defined only for Hyperons" << G4endl;
939 G4Exception("G4QHyperonElasticCrossSection::GetSlope()", "HAD_CHPS_0000",
940 FatalException, ed);
941 }
942 if(theB1<0.) theB1=0.;
943 if(!(theB1>=-1.||theB1<=1.)) G4cout<<"*NAN*G4QHyElasticCrossS::Getslope:"<<theB1<<G4endl;
944 return theB1/GeVSQ;
945}
946
947// Returns half max(Q2=-t) in independent units (MeV^2)
949{
950 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
951 return lastTM*HGeVSQ;
952}
953
954// lastLP is used, so calculating tables, one need to remember and then recover lastLP
955G4double G4QHyperonElasticCrossSection::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
956 G4int tgN)
957{
958 if(PDG==3222 || PDG<3000 || PDG>3334) G4cout<<"*Warning*G4QHypElCS::GTV:P="<<PDG<<G4endl;
959 if(tgZ<0 || tgZ>92)
960 {
961 G4cout<<"*Warning*G4QHyperonElastCS::GetTabValue:(1-92) NoIsotopesFor Z="<<tgZ<<G4endl;
962 return 0.;
963 }
964 G4int iZ=tgZ-1; // Z index
965 if(iZ<0)
966 {
967 iZ=0; // conversion of the neutron target to the proton target
968 tgZ=1;
969 tgN=0;
970 }
971 //if(nN[iZ][0] < 0)
972 //{
973#ifdef isodebug
974 // G4cout<<"*Warning*G4QHyperonElasticCS::GetTabValue: No isotopes for Z="<<tgZ<<G4endl;
975#endif
976 // return 0.;
977 //}
978#ifdef pdebug
979 G4cout<<"G4QHyElasticCS::GetTabVal:l="<<lp<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
980#endif
981 G4double p=std::exp(lp); // momentum
982 G4double sp=std::sqrt(p); // sqrt(p)
983 G4double p2=p*p;
984 G4double p3=p2*p;
985 G4double p4=p3*p;
986 if ( tgZ == 1 && tgN == 0 ) // Hyperon+P
987 {
988 G4double dl2=lp-lastPAR[9];
989 theSS=lastPAR[32];
990 theS1=(lastPAR[10]+lastPAR[11]*dl2*dl2)/(1.+lastPAR[12]/p4/p)+
991 (lastPAR[13]/p2+lastPAR[14]*p)/(p4+lastPAR[15]*sp);
992 theB1=lastPAR[16]*std::pow(p,lastPAR[17])/(1.+lastPAR[18]/p3);
993 theS2=lastPAR[19]+lastPAR[20]/(p4+lastPAR[21]*p);
994 theB2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]/sp);
995 theS3=lastPAR[25]+lastPAR[26]/(p4*p4+lastPAR[27]*p2+lastPAR[28]);
996 theB3=lastPAR[29]+lastPAR[30]/(p4+lastPAR[31]);
997 theS4=0.;
998 theB4=0.;
999#ifdef tdebug
1000 G4cout<<"G4QHyperonElasticCrossSection::GetTableVal:TM="<<lastTM<<",S1="<<theS1<<",B1="
1001 <<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS1<<",B3="<<theB1<<G4endl;
1002#endif
1003 // Returns the total elastic pim-p cross-section (to avoid spoiling lastSIG)
1004 G4double dp=lp-lastPAR[4];
1005 return lastPAR[0]/(lastPAR[1]+p2*(lastPAR[2]+p2))+(lastPAR[3]*dp*dp+lastPAR[5]+
1006 lastPAR[6]/p2)/(1.+lastPAR[7]/sp+lastPAR[8]/p4);
1007 }
1008 else
1009 {
1010 G4double p5=p4*p;
1011 G4double p6=p5*p;
1012 G4double p8=p6*p2;
1013 G4double p10=p8*p2;
1014 G4double p12=p10*p2;
1015 G4double p16=p8*p8;
1016 //G4double p24=p16*p8;
1017 G4double dl=lp-5.;
1018 G4double a=tgZ+tgN;
1019 G4double pah=std::pow(p,a/2);
1020 G4double pa=pah*pah;
1021 G4double pa2=pa*pa;
1022 if(a<6.5)
1023 {
1024 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
1025 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
1026 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
1027 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
1028 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
1029 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
1030 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
1031 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
1032 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
1033 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
1034 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
1035#ifdef tdebug
1036 G4cout<<"G4QHyperonElasticCS::GetTabV: lA, p="<<p<<",S1="<<theS1<<",B1="<<theB1
1037 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
1038 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
1039#endif
1040 }
1041 else
1042 {
1043 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
1044 lastPAR[13]/(p5+lastPAR[14]/p16);
1045 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
1046 lastPAR[17]/(1.+lastPAR[18]/p4);
1047 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
1048 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
1049 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
1050 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
1051 lastPAR[33]/(1.+lastPAR[34]/p6);
1052 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
1053 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
1054 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
1055 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
1056#ifdef tdebug
1057 G4cout<<"G4QHyperonElasticCS::GetTabV:hA,p="<<p<<",S1="<<theS1<<",B1="<<theB1<<",SS="
1058 <<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",S4="
1059 <<theS4<<",B4="<<theB4<<G4endl;
1060#endif
1061 }
1062 // Returns the total elastic (n/p)A cross-section (to avoid spoiling lastSIG)
1063#ifdef tdebug
1064 G4cout<<"G4QHyperonElCS::GetTabV: PDG="<<PDG<<",P="<<p<<",N="<<tgN<<",Z="<<tgZ<<G4endl;
1065#endif
1066 G4double dlp=lp-lastPAR[5]; // ax
1067 // p1 p2 p3 p4 p5
1068 return (lastPAR[0]*dlp*dlp+lastPAR[1])/(1.+lastPAR[2]/p)+lastPAR[3]/(p3+lastPAR[4]);
1069 }
1070 return 0.;
1071} // End of GetTableValues
1072
1073// Returns max -t=Q2 (GeV^2) for the momentum pP(GeV) and the target nucleus (tgN,tgZ)
1074G4double G4QHyperonElasticCrossSection::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
1075 G4double pP)
1076{
1077 //static const G4double mNeut= G4QPDGCode(2112).GetMass()*.001; // MeV to GeV
1078 //static const G4double mPi= G4QPDGCode(211).GetMass()*.001; // pion mass MeV to GeV
1079 //static const G4double mProt= G4QPDGCode(2212).GetMass()*.001; // MeV to GeV
1080 static const G4double mLamb= G4QPDGCode(3122).GetMass()*.001; // MeV to GeV
1081 //static const G4double mHe3 = G4QPDGCode(2112).GetNuclMass(2,1,0)*.001; // MeV to GeV
1082 //static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0)*.001; // MeV to GeV
1083 //static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0)*.001; // MeV to GeV
1084 //static const G4double mPi2= mPi*mPi;
1085 static const G4double mLa2= mLamb*mLamb;
1086 //static const G4double mProt2= mProt*mProt;
1087 //static const G4double mNeut2= mNeut*mNeut;
1088 //static const G4double mDeut2= mDeut*mDeut;
1089 G4double pP2=pP*pP; // squared momentum of the projectile
1090 if(tgZ || tgN>-1) // --> Hyperon-A
1091 {
1092 G4double mt=G4QPDGCode(90000000+tgZ*1000+tgN).GetMass()*.001; // Target mass in GeV
1093 G4double dmt=mt+mt;
1094 G4double s_value=dmt*std::sqrt(pP2+mLa2)+mLa2+mt*mt; // Mondelstam s (@@ other hyperons?)
1095 return dmt*dmt*pP2/s_value;
1096 }
1097 else
1098 {
1099 // G4cout<<"*Error*G4QHyperonElasticCrossSection::GetQ2ma:PDG="<<PDG<<",Z="<<tgZ<<",N="
1100 // <<tgN<<", while it is defined only for p projectiles & Z_target>0"<<G4endl;
1101 // throw G4QException("G4QHyperonElasticCrossSection::GetQ2max: only HyperA implemented");
1103 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
1104 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
1105 G4Exception("G4QHyperonElasticCrossSection::GetQ2max()", "HAD_CHPS_0000",
1106 FatalException, ed);
1107 return 0;
1108 }
1109}
@ FatalException
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int pPDG, G4int Z, G4int N, G4double pP)
virtual G4double GetCrossSection(G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=3122)
G4double GetSlope(G4int tZ, G4int tN, G4int pPDG)
G4double GetMass()
Definition: G4QPDGCode.cc:693
virtual G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=0)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76