Geant4-11
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ChipsKaonMinusElasticXS.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//
28//
29// G4 Physics class: G4ChipsKaonMinusElasticXS for pA elastic cross sections
30// Created: M.V. Kossov, CERN/ITEP(Moscow), 5-Feb-2010
31// The last update: M.V. Kossov, CERN/ITEP (Moscow) 5-Feb-2010
32//
33// -------------------------------------------------------------------------------
34// Short description: Interaction cross-sections for the elastic process.
35// Class extracted from CHIPS and integrated in Geant4 by W.Pokorski
36// -------------------------------------------------------------------------------
37//
38
40#include "G4SystemOfUnits.hh"
41#include "G4DynamicParticle.hh"
43#include "G4KaonMinus.hh"
44#include "G4Nucleus.hh"
45#include "G4ParticleTable.hh"
46#include "G4NucleiProperties.hh"
47#include "G4IonTable.hh"
48#include "G4AutoLock.hh"
49
50// factory
52//
54
55
56namespace {
57 G4double mK;//= G4KaonMinus::KaonMinus()->GetPDGMass()*.001; // MeV to GeV//Cannot initialize here, needs particles
58 G4double mK2;//= mK*mK;
61 const G4double third=1./3.;
62 const G4double fifth=1./5.;
63 const G4double sevth=1./7.;
64 const G4double pwd=2727;
65 const G4int n_kmpel=36; // #of parameters for pp-elastic ( < nPoints=128)
66 // -0- -1- -2- -3- -4- -5- -6- -7- -8- -9--10- -11--12-
67 const G4double kmp_el[n_kmpel]={5.2,.0557,3.5,2.23,.7,.075,.004,.39,.000156,.15,1.,.0156,5.,
68 74.,3.,3.4,.2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,
69 1.2e6,3.5e6,5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
70 // -13-14--15-16--17--18--19- -20- -21- -22- -23- -24- -25-
71 // -26- -27- -28- -29- -30- -31- -32- -33- -34- -35-
73}
74
76{
77 G4AutoLock l(&initM);
79 mK2 = mK*mK;
80 l.unlock();
81 lPMin=-8.; //Min tabulatedLogarithmMomentum/D
82 lPMax= 8.; //Max tabulatedLogarithmMomentum/D
83 dlnP=(lPMax-lPMin)/nLast;// LogStep inTable /D
84 onlyCS=true;//Flag toCalculOnlyCS(not Si/Bi)/L
85 lastSIG=0.; //Last calculated cross section /L
86 lastLP=-10.;//LastLog(mom_of IncidentHadron)/L
87 lastTM=0.; //Last t_maximum /L
88 theSS=0.; //TheLastSqSlope of 1st difr.Max/L
89 theS1=0.; //TheLastMantissa of 1st difrMax/L
90 theB1=0.; //TheLastSlope of 1st difructMax/L
91 theS2=0.; //TheLastMantissa of 2nd difrMax/L
92 theB2=0.; //TheLastSlope of 2nd difructMax/L
93 theS3=0.; //TheLastMantissa of 3d difr.Max/L
94 theB3=0.; //TheLastSlope of 3d difruct.Max/L
95 theS4=0.; //TheLastMantissa of 4th difrMax/L
96 theB4=0.; //TheLastSlope of 4th difructMax/L
97 lastTZ=0; // Last atomic number of theTarget
98 lastTN=0; // Last # of neutrons in theTarget
99 lastPIN=0.;// Last initialized max momentum
100 lastCST=0; // Elastic cross-section table
101 lastPAR=0; // ParametersForFunctionCalculation
102 lastSST=0; // E-dep ofSqardSlope of 1st difMax
103 lastS1T=0; // E-dep of mantissa of 1st dif.Max
104 lastB1T=0; // E-dep of the slope of 1st difMax
105 lastS2T=0; // E-dep of mantissa of 2nd difrMax
106 lastB2T=0; // E-dep of the slope of 2nd difMax
107 lastS3T=0; // E-dep of mantissa of 3d difr.Max
108 lastB3T=0; // E-dep of the slope of 3d difrMax
109 lastS4T=0; // E-dep of mantissa of 4th difrMax
110 lastB4T=0; // E-dep of the slope of 4th difMax
111 lastN=0; // The last N of calculated nucleus
112 lastZ=0; // The last Z of calculated nucleus
113 lastP=0.; // LastUsed inCrossSection Momentum
114 lastTH=0.; // Last threshold momentum
115 lastCS=0.; // Last value of the Cross Section
116 lastI=0; // The last position in the DAMDB
117}
118
120{
121 std::vector<G4double*>::iterator pos;
122 for (pos=CST.begin(); pos<CST.end(); pos++)
123 { delete [] *pos; }
124 CST.clear();
125 for (pos=PAR.begin(); pos<PAR.end(); pos++)
126 { delete [] *pos; }
127 PAR.clear();
128 for (pos=SST.begin(); pos<SST.end(); pos++)
129 { delete [] *pos; }
130 SST.clear();
131 for (pos=S1T.begin(); pos<S1T.end(); pos++)
132 { delete [] *pos; }
133 S1T.clear();
134 for (pos=B1T.begin(); pos<B1T.end(); pos++)
135 { delete [] *pos; }
136 B1T.clear();
137 for (pos=S2T.begin(); pos<S2T.end(); pos++)
138 { delete [] *pos; }
139 S2T.clear();
140 for (pos=B2T.begin(); pos<B2T.end(); pos++)
141 { delete [] *pos; }
142 B2T.clear();
143 for (pos=S3T.begin(); pos<S3T.end(); pos++)
144 { delete [] *pos; }
145 S3T.clear();
146 for (pos=B3T.begin(); pos<B3T.end(); pos++)
147 { delete [] *pos; }
148 B3T.clear();
149 for (pos=S4T.begin(); pos<S4T.end(); pos++)
150 { delete [] *pos; }
151 S4T.clear();
152 for (pos=B4T.begin(); pos<B4T.end(); pos++)
153 { delete [] *pos; }
154 B4T.clear();
155}
156
157void
159{
160 outFile << "G4ChipsKaonMinusElasticXS provides the elastic cross\n"
161 << "section for K- nucleus scattering as a function of incident\n"
162 << "momentum. The cross section is calculated using M. Kossov's\n"
163 << "CHIPS parameterization of cross section data.\n";
164}
165
167 const G4Element*,
168 const G4Material*)
169{
170 return true;
171}
172
173// The main member function giving the collision cross section (P is in IU, CS is in mb)
174// Make pMom in independent units ! (Now it is MeV)
176 const G4Isotope*,
177 const G4Element*,
178 const G4Material*)
179{
180 G4double pMom=Pt->GetTotalMomentum();
181 G4int tgN = A - tgZ;
182
183 return GetChipsCrossSection(pMom, tgZ, tgN, -321);
184}
185
187{
188
189 G4bool fCS = false;
190
191 G4double pEn=pMom;
192 onlyCS=fCS;
193
194 G4bool in=false; // By default the isotope must be found in the AMDB
195 lastP = 0.; // New momentum history (nothing to compare with)
196 lastN = tgN; // The last N of the calculated nucleus
197 lastZ = tgZ; // The last Z of the calculated nucleus
198 lastI = colN.size(); // Size of the Associative Memory DB in the heap
199 if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
200 { // The nucleus with projPDG is found in AMDB
201 if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
202 {
203 lastI=i;
204 lastTH =colTH[i]; // Last THreshold (A-dependent)
205 if(pEn<=lastTH)
206 {
207 return 0.; // Energy is below the Threshold value
208 }
209 lastP =colP [i]; // Last Momentum (A-dependent)
210 lastCS =colCS[i]; // Last CrossSect (A-dependent)
211 // if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance)
212 if(lastP == pMom) // Do not recalculate
213 {
214 CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom); // Update param's only
215 return lastCS*millibarn; // Use theLastCS
216 }
217 in = true; // This is the case when the isotop is found in DB
218 // Momentum pMom is in IU ! @@ Units
219 lastCS=CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom); // read & update
220 if(lastCS<=0. && pEn>lastTH) // Correct the threshold
221 {
222 lastTH=pEn;
223 }
224 break; // Go out of the LOOP with found lastI
225 }
226 } // End of attampt to find the nucleus in DB
227 if(!in) // This nucleus has not been calculated previously
228 {
230 lastCS=CalculateCrossSection(fCS,0,lastI,-321,lastZ,lastN,pMom);//calculate&create
231 if(lastCS<=0.)
232 {
233 lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
234 if(pEn>lastTH)
235 {
236 lastTH=pEn;
237 }
238 }
239 colN.push_back(tgN);
240 colZ.push_back(tgZ);
241 colP.push_back(pMom);
242 colTH.push_back(lastTH);
243 colCS.push_back(lastCS);
244 return lastCS*millibarn;
245 } // End of creation of the new set of parameters
246 else
247 {
248 colP[lastI]=pMom;
250 }
251 return lastCS*millibarn;
252}
253
254// Calculation of total elastic cross section (p in IU, CS in mb) @@ Units (?)
255// F=0 - create AMDB, F=-1 - read&update AMDB, F=1 - update AMDB (sinchro with higher AMDB)
257 G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
258{
259 G4double pMom=pIU/GeV; // All calculations are in GeV
260 onlyCS=CS; // Flag to calculate only CS (not Si/Bi)
261 lastLP=std::log(pMom); // Make a logarithm of the momentum for calculation
262 if(F) // This isotope was found in AMDB =>RETRIEVE/UPDATE
263 {
264 if(F<0) // the AMDB must be loded
265 {
266 lastPIN = PIN[I]; // Max log(P) initialised for this table set
267 lastPAR = PAR[I]; // Pointer to the parameter set
268 lastCST = CST[I]; // Pointer to the total sross-section table
269 lastSST = SST[I]; // Pointer to the first squared slope
270 lastS1T = S1T[I]; // Pointer to the first mantissa
271 lastB1T = B1T[I]; // Pointer to the first slope
272 lastS2T = S2T[I]; // Pointer to the second mantissa
273 lastB2T = B2T[I]; // Pointer to the second slope
274 lastS3T = S3T[I]; // Pointer to the third mantissa
275 lastB3T = B3T[I]; // Pointer to the rhird slope
276 lastS4T = S4T[I]; // Pointer to the 4-th mantissa
277 lastB4T = B4T[I]; // Pointer to the 4-th slope
278 }
280 {
281 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
282 PIN[I]=lastPIN; // Remember the new P-Limit of the tables
283 }
284 }
285 else // This isotope wasn't initialized => CREATE
286 {
287 lastPAR = new G4double[nPoints]; // Allocate memory for parameters of CS function
288 lastPAR[nLast]=0; // Initialization for VALGRIND
289 lastCST = new G4double[nPoints]; // Allocate memory for Tabulated CS function
290 lastSST = new G4double[nPoints]; // Allocate memory for Tabulated first sqaredSlope
291 lastS1T = new G4double[nPoints]; // Allocate memory for Tabulated first mantissa
292 lastB1T = new G4double[nPoints]; // Allocate memory for Tabulated first slope
293 lastS2T = new G4double[nPoints]; // Allocate memory for Tabulated second mantissa
294 lastB2T = new G4double[nPoints]; // Allocate memory for Tabulated second slope
295 lastS3T = new G4double[nPoints]; // Allocate memory for Tabulated third mantissa
296 lastB3T = new G4double[nPoints]; // Allocate memory for Tabulated third slope
297 lastS4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th mantissa
298 lastB4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th slope
299 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
300 PIN.push_back(lastPIN); // Fill parameters of CS function to AMDB
301 PAR.push_back(lastPAR); // Fill parameters of CS function to AMDB
302 CST.push_back(lastCST); // Fill Tabulated CS function to AMDB
303 SST.push_back(lastSST); // Fill Tabulated first sq.slope to AMDB
304 S1T.push_back(lastS1T); // Fill Tabulated first mantissa to AMDB
305 B1T.push_back(lastB1T); // Fill Tabulated first slope to AMDB
306 S2T.push_back(lastS2T); // Fill Tabulated second mantissa to AMDB
307 B2T.push_back(lastB2T); // Fill Tabulated second slope to AMDB
308 S3T.push_back(lastS3T); // Fill Tabulated third mantissa to AMDB
309 B3T.push_back(lastB3T); // Fill Tabulated third slope to AMDB
310 S4T.push_back(lastS4T); // Fill Tabulated 4-th mantissa to AMDB
311 B4T.push_back(lastB4T); // Fill Tabulated 4-th slope to AMDB
312 } // End of creation/update of the new set of parameters and tables
313 // =----------= NOW Update (if necessary) and Calculate the Cross Section =-----------=
315 {
316 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
317 }
318 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
319 if(lastLP>lPMin && lastLP<=lastPIN) // Linear fit is made using precalculated tables
320 {
321 if(lastLP==lastPIN)
322 {
323 G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
324 G4int blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
325 if(blast<0 || blast>=nLast) G4cout<<"G4QKMElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl;
326 lastSIG = lastCST[blast];
327 if(!onlyCS) // Skip the differential cross-section parameters
328 {
329 theSS = lastSST[blast];
330 theS1 = lastS1T[blast];
331 theB1 = lastB1T[blast];
332 theS2 = lastS2T[blast];
333 theB2 = lastB2T[blast];
334 theS3 = lastS3T[blast];
335 theB3 = lastB3T[blast];
336 theS4 = lastS4T[blast];
337 theB4 = lastB4T[blast];
338 }
339 }
340 else
341 {
342 G4double shift=(lastLP-lPMin)/dlnP; // a shift from the beginning of the table
343 G4int blast=static_cast<int>(shift); // the lower bin number
344 if(blast<0) blast=0;
345 if(blast>=nLast) blast=nLast-1; // low edge of the last bin
346 shift-=blast; // step inside the unit bin
347 G4int lastL=blast+1; // the upper bin number
348 G4double SIGL=lastCST[blast]; // the basic value of the cross-section
349 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
350 if(!onlyCS) // Skip the differential cross-section parameters
351 {
352 G4double SSTL=lastSST[blast]; // the low bin of the first squared slope
353 theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
354 G4double S1TL=lastS1T[blast]; // the low bin of the first mantissa
355 theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
356 G4double B1TL=lastB1T[blast]; // the low bin of the first slope
357 theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
358 G4double S2TL=lastS2T[blast]; // the low bin of the second mantissa
359 theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
360 G4double B2TL=lastB2T[blast]; // the low bin of the second slope
361 theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
362 G4double S3TL=lastS3T[blast]; // the low bin of the third mantissa
363 theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
364 G4double B3TL=lastB3T[blast]; // the low bin of the third slope
365 theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
366 G4double S4TL=lastS4T[blast]; // the low bin of the 4-th mantissa
367 theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
368 G4double B4TL=lastB4T[blast]; // the low bin of the 4-th slope
369 theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
370 }
371 }
372 }
373 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
374 if(lastSIG<0.) lastSIG = 0.; // @@ a Warning print can be added
375 return lastSIG;
376}
377
378// It has parameter sets for all tZ/tN/PDG, using them the tables can be created/updated
380 G4int tgZ, G4int tgN)
381{
382 // @@ At present all nA==pA ---------> Each neucleus can have not more than 51 parameters
383 if(PDG == -321)
384 {
385 // -- Total pp elastic cross section cs & s1/b1 (main), s2/b2 (tail1), s3/b3 (tail2) --
386 //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|
387 //CS=2.865/p2s/(1+.0022/p2s)+(18.9+.6461*dl1*dl1+9./p)/(1.+.425*lp)/(1.+.4276/p4);
388 // par(0) par(7) par(1) par(2) par(4) par(5) par(6)
389 //dl2=lp-5., s1=(74.+3.*dl2*dl2)/(1+3.4/p4/p)+(.2/p2+17.*p)/(p4+.001*sp),
390 // par(8) par(9) par(10) par(11) par(12)par(13) par(14)
391 // b1=8.*p**.055/(1.+3.64/p3); s2=5.e-5+4000./(p4+1500.*p); b2=.46+1.2e6/(p4+3.5e6/sp);
392 // par(15) par(16) par(17) par(18) par(19) par(20) par(21) par(22) par(23)
393 // s3=5.e-5+1.e10/(p4*p4+8.5e8*p2+1.e10); b3=1.1+3.4e6/(p4+6.8e6); ss=0.
394 // par(24) par(25) par(26) par(27) par(28) par(29) par(30) par(31)
395 //
396 if(lastPAR[nLast]!=pwd) // A unique flag to avoid the repeatable definition
397 {
398 if ( tgZ == 1 && tgN == 0 )
399 {
400 for (G4int ip=0; ip<n_kmpel; ip++) lastPAR[ip]=kmp_el[ip]; // PiMinus+P
401 }
402 else
403 {
404 G4double a=tgZ+tgN;
405 G4double sa=std::sqrt(a);
406 G4double ssa=std::sqrt(sa);
407 G4double asa=a*sa;
408 G4double a2=a*a;
409 G4double a3=a2*a;
410 G4double a4=a3*a;
411 G4double a5=a4*a;
412 G4double a6=a4*a2;
413 G4double a7=a6*a;
414 G4double a8=a7*a;
415 G4double a9=a8*a;
416 G4double a10=a5*a5;
417 G4double a12=a6*a6;
418 G4double a14=a7*a7;
419 G4double a16=a8*a8;
420 G4double a17=a16*a;
421 //G4double a20=a16*a4;
422 G4double a32=a16*a16;
423 // Reaction cross-section parameters (kmael_fit.f)
424 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa)); // p1
425 lastPAR[1]=.75*asa/(1.+.009*a); // p2
426 lastPAR[2]=.1*a2*ssa/(1.+.0015*a2/ssa); // p3
427 lastPAR[3]=1./(1.+500./a2); // p4
428 lastPAR[4]=4.2; // p5
429 lastPAR[5]=0.; // p6 not used
430 lastPAR[6]=0.; // p7 not used
431 lastPAR[7]=0.; // p8 not used
432 lastPAR[8]=0.; // p9 not used
433 // @@ the differential cross-section is parameterized separately for A>6 & A<7
434 if(a<6.5)
435 {
436 G4double a28=a16*a12;
437 // The main pre-exponent (pel_sg)
438 lastPAR[ 9]=4000*a; // p1
439 lastPAR[10]=1.2e7*a8+380*a17; // p2
440 lastPAR[11]=.7/(1.+4.e-12*a16); // p3
441 lastPAR[12]=2.5/a8/(a4+1.e-16*a32); // p4
442 lastPAR[13]=.28*a; // p5
443 lastPAR[14]=1.2*a2+2.3; // p6
444 lastPAR[15]=3.8/a; // p7
445 // The main slope (pel_sl)
446 lastPAR[16]=.01/(1.+.0024*a5); // p1
447 lastPAR[17]=.2*a; // p2
448 lastPAR[18]=9.e-7/(1.+.035*a5); // p3
449 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a); // p4
450 // The main quadratic (pel_sh)
451 lastPAR[20]=2.25*a3; // p1
452 lastPAR[21]=18.; // p2
453 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7); // p3
454 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a); // p4
455 // The 1st max pre-exponent (pel_qq)
456 lastPAR[24]=1.e5/(a8+2.5e12/a16); // p1
457 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28); // p2
458 lastPAR[26]=.0006*a3; // p3
459 // The 1st max slope (pel_qs)
460 lastPAR[27]=10.+4.e-8*a12*a; // p1
461 lastPAR[28]=.114; // p2
462 lastPAR[29]=.003; // p3
463 lastPAR[30]=2.e-23; // p4
464 // The effective pre-exponent (pel_ss)
465 lastPAR[31]=1./(1.+.0001*a8); // p1
466 lastPAR[32]=1.5e-4/(1.+5.e-6*a12); // p2
467 lastPAR[33]=.03; // p3
468 // The effective slope (pel_sb)
469 lastPAR[34]=a/2; // p1
470 lastPAR[35]=2.e-7*a4; // p2
471 lastPAR[36]=4.; // p3
472 lastPAR[37]=64./a3; // p4
473 // The gloria pre-exponent (pel_us)
474 lastPAR[38]=1.e8*std::exp(.32*asa); // p1
475 lastPAR[39]=20.*std::exp(.45*asa); // p2
476 lastPAR[40]=7.e3+2.4e6/a5; // p3
477 lastPAR[41]=2.5e5*std::exp(.085*a3); // p4
478 lastPAR[42]=2.5*a; // p5
479 // The gloria slope (pel_ub)
480 lastPAR[43]=920.+.03*a8*a3; // p1
481 lastPAR[44]=93.+.0023*a12; // p2
482 }
483 else
484 {
485 G4double p1a10=2.2e-28*a10;
486 G4double r4a16=6.e14/a16;
487 G4double s4a16=r4a16*r4a16;
488 // a24
489 // a36
490 // The main pre-exponent (peh_sg)
491 lastPAR[ 9]=4.5*std::pow(a,1.15); // p1
492 lastPAR[10]=.06*std::pow(a,.6); // p2
493 lastPAR[11]=.6*a/(1.+2.e15/a16); // p3
494 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32); // p4
495 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5); // p5
496 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12); // p6
497 // The main slope (peh_sl)
498 lastPAR[15]=400./a12+2.e-22*a9; // p1
499 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14); // p2
500 lastPAR[17]=1000./a2+9.5*sa*ssa; // p3
501 lastPAR[18]=4.e-6*a*asa+1.e11/a16; // p4
502 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16); // p5
503 lastPAR[20]=9.+100./a; // p6
504 // The main quadratic (peh_sh)
505 lastPAR[21]=.002*a3+3.e7/a6; // p1
506 lastPAR[22]=7.e-15*a4*asa; // p2
507 lastPAR[23]=9000./a4; // p3
508 // The 1st max pre-exponent (peh_qq)
509 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4); // p1
510 lastPAR[25]=1.e-5*a2+2.e14/a16; // p2
511 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12); // p3
512 lastPAR[27]=.016*asa/(1.+5.e16/a16); // p4
513 // The 1st max slope (peh_qs)
514 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14)); // p1
515 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11); // p2
516 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8); // p3
517 lastPAR[31]=100./asa; // p4
518 // The 2nd max pre-exponent (peh_ss)
519 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4); // p1
520 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8); // p2
521 lastPAR[34]=1.3+3.e5/a4; // p3
522 lastPAR[35]=500./(a2+50.)+3; // p4
523 lastPAR[36]=1.e-9/a+s4a16*s4a16; // p5
524 // The 2nd max slope (peh_sb)
525 lastPAR[37]=.4*asa+3.e-9*a6; // p1
526 lastPAR[38]=.0005*a5; // p2
527 lastPAR[39]=.002*a5; // p3
528 lastPAR[40]=10.; // p4
529 // The effective pre-exponent (peh_us)
530 lastPAR[41]=.05+.005*a; // p1
531 lastPAR[42]=7.e-8/sa; // p2
532 lastPAR[43]=.8*sa; // p3
533 lastPAR[44]=.02*sa; // p4
534 lastPAR[45]=1.e8/a3; // p5
535 lastPAR[46]=3.e32/(a32+1.e32); // p6
536 // The effective slope (peh_ub)
537 lastPAR[47]=24.; // p1
538 lastPAR[48]=20./sa; // p2
539 lastPAR[49]=7.e3*a/(sa+1.); // p3
540 lastPAR[50]=900.*sa/(1.+500./a3); // p4
541 }
542 // Parameter for lowEnergyNeutrons
543 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
544 }
546 // and initialize the zero element of the table
547 G4double lp=lPMin; // ln(momentum)
548 G4bool memCS=onlyCS; // ??
549 onlyCS=false;
550 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables
551 onlyCS=memCS;
552 lastSST[0]=theSS;
553 lastS1T[0]=theS1;
554 lastB1T[0]=theB1;
555 lastS2T[0]=theS2;
556 lastB2T[0]=theB2;
557 lastS3T[0]=theS3;
558 lastB3T[0]=theB3;
559 lastS4T[0]=theS4;
560 lastB4T[0]=theB4;
561 }
562 if(LP>ILP)
563 {
564 G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1; // already inited till this
565 if(ini<0) ini=0;
566 if(ini<nPoints)
567 {
568 G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1; // final bin of initialization
569 if(fin>=nPoints) fin=nLast; // Limit of the tabular initialization
570 if(fin>=ini)
571 {
572 G4double lp=0.;
573 for(G4int ip=ini; ip<=fin; ip++) // Calculate tabular CS,S1,B1,S2,B2,S3,B3
574 {
575 lp=lPMin+ip*dlnP; // ln(momentum)
576 G4bool memCS=onlyCS;
577 onlyCS=false;
578 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables (ret CS)
579 onlyCS=memCS;
580 lastSST[ip]=theSS;
581 lastS1T[ip]=theS1;
582 lastB1T[ip]=theB1;
583 lastS2T[ip]=theS2;
584 lastB2T[ip]=theB2;
585 lastS3T[ip]=theS3;
586 lastB3T[ip]=theB3;
587 lastS4T[ip]=theS4;
588 lastB4T[ip]=theB4;
589 }
590 return lp;
591 }
592 else G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetPTables: PDG="<<PDG
593 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
594 <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
595 }
596 else G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetPTables: PDG="<<PDG
597 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
598 <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
599 }
600 }
601 else
602 {
603 // G4cout<<"*Error*G4ChipsKaonMinusElasticXS::GetPTables: PDG="<<PDG<<", Z="<<tgZ
604 // <<", N="<<tgN<<", while it is defined only for PDG=-321"<<G4endl;
605 // throw G4QException("G4ChipsKaonMinusElasticXS::GetPTables:onlyK-'s implemented");
607 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
608 << ", while it is defined only for PDG=-321 (K-) " << G4endl;
609 G4Exception("G4ChipsKaonMinusElasticXS::GetPTables()", "HAD_CHPS_0000",
610 FatalException, ed);
611 }
612 return ILP;
613}
614
615// Returns Q2=-t in independent units (MeV^2) (all internal calculations are in GeV)
617{
618 if(PDG==310 || PDG==130) PDG=-321;
619 if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetET:PDG="<<PDG<<G4endl;
620 if(onlyCS) G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetExT: onlyCS=1"<<G4endl;
621 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
622 G4double q2=0.;
623 if(tgZ==1 && tgN==0) // ===> p+p=p+p
624 {
626 G4double R1=(1.-std::exp(-E1));
628 G4double R2=(1.-std::exp(-E2*E2*E2));
630 G4double R3=(1.-std::exp(-E3));
631 G4double I1=R1*theS1/theB1;
632 G4double I2=R2*theS2;
633 G4double I3=R3*theS3;
634 G4double I12=I1+I2;
635 G4double rand=(I12+I3)*G4UniformRand();
636 if (rand<I1 )
637 {
638 G4double ran=R1*G4UniformRand();
639 if(ran>1.) ran=1.;
640 q2=-std::log(1.-ran)/theB1;
641 }
642 else if(rand<I12)
643 {
644 G4double ran=R2*G4UniformRand();
645 if(ran>1.) ran=1.;
646 q2=-std::log(1.-ran);
647 if(q2<0.) q2=0.;
648 q2=std::pow(q2,third)/theB2;
649 }
650 else
651 {
652 G4double ran=R3*G4UniformRand();
653 if(ran>1.) ran=1.;
654 q2=-std::log(1.-ran)/theB3;
655 }
656 }
657 else
658 {
659 G4double a=tgZ+tgN;
661 G4double R1=(1.-std::exp(-E1));
662 G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
664 G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st)
665 if(a>6.5)E2*=tm2; // for heavy nuclei
666 G4double R2=(1.-std::exp(-E2));
668 if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA
669 G4double R3=(1.-std::exp(-E3));
671 G4double R4=(1.-std::exp(-E4));
672 G4double I1=R1*theS1;
673 G4double I2=R2*theS2;
674 G4double I3=R3*theS3;
675 G4double I4=R4*theS4;
676 G4double I12=I1+I2;
677 G4double I13=I12+I3;
678 G4double rand=(I13+I4)*G4UniformRand();
679 if(rand<I1)
680 {
681 G4double ran=R1*G4UniformRand();
682 if(ran>1.) ran=1.;
683 q2=-std::log(1.-ran)/theB1;
684 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
685 }
686 else if(rand<I12)
687 {
688 G4double ran=R2*G4UniformRand();
689 if(ran>1.) ran=1.;
690 q2=-std::log(1.-ran)/theB2;
691 if(q2<0.) q2=0.;
692 if(a<6.5) q2=std::pow(q2,third);
693 else q2=std::pow(q2,fifth);
694 }
695 else if(rand<I13)
696 {
697 G4double ran=R3*G4UniformRand();
698 if(ran>1.) ran=1.;
699 q2=-std::log(1.-ran)/theB3;
700 if(q2<0.) q2=0.;
701 if(a>6.5) q2=std::pow(q2,sevth);
702 }
703 else
704 {
705 G4double ran=R4*G4UniformRand();
706 if(ran>1.) ran=1.;
707 q2=-std::log(1.-ran)/theB4;
708 if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0)
709 }
710 }
711 if(q2<0.) q2=0.;
712 if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<G4endl;
713 if(q2>lastTM)
714 {
715 q2=lastTM;
716 }
717 return q2*GeVSQ;
718}
719
720// Returns B in independent units (MeV^-2) (all internal calculations are in GeV) see ExT
722{
723 if(onlyCS)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetSl:onlCS=true"<<G4endl;
724 if(lastLP<-4.3) return 0.; // S-wave for p<14 MeV/c (kinE<.1MeV)
725 if(PDG != -321)
726 {
727 // G4cout<<"*Error*G4ChipsKaonMinusElasticXS::GetSlope: PDG="<<PDG<<", Z="<<tgZ
728 // <<", N="<<tgN<<", while it is defined only for PDG=-321"<<G4endl;
729 // throw G4QException("G4ChipsKaonMinusElasticXS::GetSlope:Only K- is implemented");
731 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
732 << ", while it is defined only for PDG=-321 (K-)" << G4endl;
733 }
734 if(theB1<0.) theB1=0.;
735 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<G4endl;
736 return theB1/GeVSQ;
737}
738
739// Returns half max(Q2=-t) in independent units (MeV^2)
741{
742 return lastTM*HGeVSQ;
743}
744
745// lastLP is used, so calculating tables, one need to remember and then recover lastLP
747 G4int tgN)
748{
749 if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetTV:PDG="<<PDG<<G4endl;
750
751 //AR-24Apr2018 Switch to allow transuranic elements
752 const G4bool isHeavyElementAllowed = true;
753 if(tgZ<0 || ( !isHeavyElementAllowed && tgZ>92))
754 {
755 G4cout<<"*Warning*G4QKaonMinusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
756 return 0.;
757 }
758 G4int iZ=tgZ-1; // Z index
759 if(iZ<0)
760 {
761 iZ=0; // conversion of the neutron target to the proton target
762 tgZ=1;
763 tgN=0;
764 }
765 G4double p=std::exp(lp); // momentum
766 G4double sp=std::sqrt(p); // sqrt(p)
767 G4double psp=p*sp; // p*sqrt(p)
768 G4double p2=p*p;
769 G4double p3=p2*p;
770 G4double p4=p3*p;
771 if ( tgZ == 1 && tgN == 0 ) // KaonMinus+P
772 {
773 G4double dl2=lp-lastPAR[12];
774 theSS=lastPAR[35];
775 theS1=(lastPAR[13]+lastPAR[14]*dl2*dl2)/(1.+lastPAR[15]/p4/p)+
776 (lastPAR[16]/p2+lastPAR[17]*p)/(p4+lastPAR[18]*sp);
777 theB1=lastPAR[19]*std::pow(p,lastPAR[20])/(1.+lastPAR[21]/p3);
778 theS2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]*p);
779 theB2=lastPAR[25]+lastPAR[26]/(p4+lastPAR[27]/sp);
780 theS3=lastPAR[28]+lastPAR[29]/(p4*p4+lastPAR[30]*p2+lastPAR[31]);
781 theB3=lastPAR[32]+lastPAR[33]/(p4+lastPAR[34]);
782 theS4=0.;
783 theB4=0.;
784 // Returns the total elastic pim-p cross-section (to avoid spoiling lastSIG)
785 G4double dp=lp-lastPAR[2];
786 return lastPAR[0]/psp+(lastPAR[1]*dp*dp+lastPAR[3])/(1.-lastPAR[4]/sp+lastPAR[5]/p4)+
787 lastPAR[6]/(sqr(p-lastPAR[7])+lastPAR[8])+lastPAR[9]/(sqr(p-lastPAR[10])+lastPAR[11]);
788 }
789 else
790 {
791 G4double p5=p4*p;
792 G4double p6=p5*p;
793 G4double p8=p6*p2;
794 G4double p10=p8*p2;
795 G4double p12=p10*p2;
796 G4double p16=p8*p8;
797 //G4double p24=p16*p8;
798 G4double dl=lp-5.;
799 G4double a=tgZ+tgN;
800 G4double pah=std::pow(p,a/2);
801 G4double pa=pah*pah;
802 G4double pa2=pa*pa;
803 if(a<6.5)
804 {
805 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
806 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
807 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
808 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
809 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
810 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
811 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
812 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
813 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
814 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
815 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
816 }
817 else
818 {
819 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
820 lastPAR[13]/(p5+lastPAR[14]/p16);
821 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
822 lastPAR[17]/(1.+lastPAR[18]/p4);
823 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
824 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
825 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
826 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
827 lastPAR[33]/(1.+lastPAR[34]/p6);
828 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
829 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
830 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
831 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
832 }
833 // Returns the total elastic (n/p)A cross-section (to avoid spoiling lastSIG)
834 G4double dlp=lp-lastPAR[4]; // ax
835 // p1 p2 p3 p4
836 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p3)/(1.+lastPAR[3]/p2/sp);
837 }
838 return 0.;
839} // End of GetTableValues
840
841// Returns max -t=Q2 (GeV^2) for the momentum pP(GeV) and the target nucleus (tgN,tgZ)
843 G4double pP)
844{
845 G4double pP2=pP*pP; // squared momentum of the projectile
846 if(tgZ || tgN>-1) // ---> pipA
847 {
848 G4double mt=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(tgZ,tgZ+tgN,0)->GetPDGMass()*.001; // Target mass in GeV
849
850 G4double dmt=mt+mt;
851 G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt; // Mondelstam mds
852 return dmt*dmt*pP2/mds;
853 }
854 else
855 {
857 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
858 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
859 G4Exception("G4ChipsKaonMinusElasticXS::GetQ2max()", "HAD_CHPS_0000",
860 FatalException, ed);
861 return 0;
862 }
863}
G4_DECLARE_XS_FACTORY(G4ChipsKaonMinusElasticXS)
static const G4double pos
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
static constexpr double gigaelectronvolt
Definition: G4SIunits.hh:194
static constexpr double millibarn
Definition: G4SIunits.hh:86
static constexpr double GeV
Definition: G4SIunits.hh:203
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:85
std::mutex G4Mutex
Definition: G4Threading.hh:81
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4double A[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
virtual void CrossSectionDescription(std::ostream &) const
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int tgZ, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int pPDG, G4int Z, G4int N, G4double pP)
G4double GetQ2max(G4int pPDG, G4int tgZ, G4int tgN, G4double pP)
G4double GetTabValues(G4double lp, G4int pPDG, G4int tgZ, G4int tgN)
G4double GetSlope(G4int tZ, G4int tN, G4int pPDG)
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
G4double GetPTables(G4double lpP, G4double lPm, G4int PDG, G4int tZ, G4int tN)
virtual G4bool IsIsoApplicable(const G4DynamicParticle *Pt, G4int Z, G4int A, const G4Element *elm, const G4Material *mat)
G4double GetTotalMomentum() const
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:112
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
T sqr(const T &x)
Definition: templates.hh:128