| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // The lust update: M.V. Kossov, CERN/ITEP(Mos
28 //
29 //
30 // G4 Physics class: G4ChipsKaonMinusInelastic
31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20
32 // The last update: M.V. Kossov, CERN/ITEP (Mo
33 // -------------------------------------------
34 // Short description: Cross-sections extracted
35 // kaon(minus)-nuclear interactions. Author: M
36 // -------------------------------------------
37 //
38
39 #include "G4ChipsKaonMinusInelasticXS.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DynamicParticle.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4KaonMinus.hh"
44
45 // factory
46 #include "G4CrossSectionFactory.hh"
47 //
48 G4_DECLARE_XS_FACTORY(G4ChipsKaonMinusInelasti
49
50 namespace {
51 const G4double THmin=27.; // default m
52 const G4double THmiG=THmin*.001; // minimu
53 const G4double dP=10.; // step for
54 const G4double dPG=dP*.001; // step for
55 const G4int nL=105; // A#of LEN
56 const G4double Pmin=THmin+(nL-1)*dP; // mi
57 const G4double Pmax=227000.; // maxP for
58 const G4int nH=224; // A#of HEN
59 const G4double milP=std::log(Pmin);// Low
60 const G4double malP=std::log(Pmax);// High
61 const G4double dlP=(malP-milP)/(nH-1); //
62 const G4double milPG=std::log(.001*Pmin);/
63 }
64 // Initialization of the
65
66 G4ChipsKaonMinusInelasticXS::G4ChipsKaonMinusI
67 {
68 lastLEN=0; // Pointer to lastArray of LowEn
69 lastHEN=0; // Pointer to lastArray of HighEn
70 lastN=0; // The last N of calculated nucle
71 lastZ=0; // The last Z of calculated nucle
72 lastP=0.; // Last used in CrossSection Mome
73 lastTH=0.; // Last threshold momentum
74 lastCS=0.; // Last value of the Cross Sectio
75 lastI=0; // The last position in the DAMDB
76 LEN = new std::vector<G4double*>;
77 HEN = new std::vector<G4double*>;
78 }
79
80 G4ChipsKaonMinusInelasticXS::~G4ChipsKaonMinus
81 {
82 std::size_t lens=LEN->size();
83 for(std::size_t i=0; i<lens; ++i) delete[] (
84 delete LEN;
85
86 std::size_t hens=HEN->size();
87 for(std::size_t i=0; i<hens; ++i) delete[] (
88 delete HEN;
89 }
90
91 void
92 G4ChipsKaonMinusInelasticXS::CrossSectionDescr
93 {
94 outFile << "G4ChipsKaonMinusInelasticXS pr
95 << "section for K- nucleus scatter
96 << "momentum. The cross section is
97 << "CHIPS parameterization of cros
98 }
99
100 G4bool G4ChipsKaonMinusInelasticXS::IsIsoAppli
101 const G4Element*,
102 const G4Material*)
103 {
104 return true;
105 }
106
107
108 // The main member function giving the collisi
109 // Make pMom in independent units ! (Now it is
110 G4double G4ChipsKaonMinusInelasticXS::GetIsoCr
111 const G4Isotope*,
112 const G4Element*,
113 const G4Material*)
114 {
115 G4double pMom=Pt->GetTotalMomentum();
116 G4int tgN = A - tgZ;
117
118 return GetChipsCrossSection(pMom, tgZ, tgN,
119 }
120
121 G4double G4ChipsKaonMinusInelasticXS::GetChips
122 {
123 G4bool in=false; // By d
124 if(tgN!=lastN || tgZ!=lastZ) // The
125 {
126 in = false; // By d
127 lastP = 0.; // New
128 lastN = tgN; // The
129 lastZ = tgZ; // The
130 lastI = (G4int)colN.size(); // Size
131 j = 0; // A#0f
132 if(lastI) for(G4int i=0; i<lastI; ++i) //
133 {
134 if(colN[i]==tgN && colZ[i]==tgZ) // Try
135 {
136 lastI=i; // Reme
137 lastTH =colTH[i]; // The
138 if(pMom<=lastTH)
139 {
140 return 0.; // Ener
141 }
142 lastP =colP [i]; // Last
143 lastCS =colCS[i]; // Last
144 in = true; // This
145 // Momentum pMom is in IU ! @@ Units
146 lastCS=CalculateCrossSection(-1,j,-321
147 if(lastCS<=0. && pMom>lastTH) // Corr
148 {
149 lastCS=0.;
150 lastTH=pMom;
151 }
152 break; // Go o
153 }
154 j++; // Incr
155 }
156 if(!in) // This
157 {
158 //!!The slave functions must provide cro
159 lastCS=CalculateCrossSection(0,j,-321,la
160 //if(lastCS>0.) // It
161 //{
162
163 // lastTH = ThresholdEnergy(tgZ,
164
165 lastTH = 0; // WP - to be checked!!!
166 colN.push_back(tgN);
167 colZ.push_back(tgZ);
168 colP.push_back(pMom);
169 colTH.push_back(lastTH);
170 colCS.push_back(lastCS);
171 //} // M.K. Presence of H1 with high thr
172 return lastCS*millibarn;
173 } // End of creation of the new set of par
174 else
175 {
176 colP[lastI]=pMom;
177 colCS[lastI]=lastCS;
178 }
179 } // End of parameters udate
180 else if(pMom<=lastTH)
181 {
182 return 0.; // Mome
183 }
184 else // It i
185 {
186 lastCS=CalculateCrossSection(1,j,-321,last
187 lastP=pMom;
188 }
189 return lastCS*millibarn;
190 }
191
192 // The main member function giving the gamma-A
193 G4double G4ChipsKaonMinusInelasticXS::Calculat
194 G4int,
195 {
196 if(F<=0) // This
197 {
198 if(F<0) // This
199 {
200 G4int sync=(G4int)LEN->size();
201 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS
202 lastLEN=(*LEN)[I]; // Poi
203 lastHEN=(*HEN)[I]; // Poi
204 }
205 else // This
206 {
207 lastLEN = new G4double[nL]; // Allo
208 lastHEN = new G4double[nH]; // Allo
209 // --- Instead of making a separate func
210 G4double P=THmiG; // Tabl
211 for(G4int k=0; k<nL; k++)
212 {
213 lastLEN[k] = CrossSectionLin(targZ, ta
214 P+=dPG;
215 }
216 G4double lP=milPG;
217 for(G4int n=0; n<nH; n++)
218 {
219 lastHEN[n] = CrossSectionLog(targZ, ta
220 lP+=dlP;
221 }
222 // --- End of possible separate function
223 // *** The synchronization check ***
224 G4int sync=(G4int)LEN->size();
225 if(sync!=I)
226 {
227 G4cerr<<"***G4ChipsKaonMinusCS::CalcCr
228 <<", N="<<targN<<", F="<<F<<G4en
229 //G4Exception("G4PiMinusNuclearCS::Cal
230 }
231 LEN->push_back(lastLEN); // reme
232 HEN->push_back(lastHEN); // reme
233 } // End of creation of the new set of par
234 } // End of parameters update
235 // =------------------= NOW the Magic Formul
236 G4double sigma;
237 if (Momentum<lastTH) return 0.; // It m
238 else if (Momentum<Pmin) // High
239 {
240 sigma=EquLinearFit(Momentum,nL,THmin,dP,la
241 }
242 else if (Momentum<Pmax) // High
243 {
244 G4double lP=std::log(Momentum);
245 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN)
246 }
247 else // UHE
248 {
249 G4double P=0.001*Momentum; // Appr
250 sigma=CrossSectionFormula(targZ, targN, P,
251 }
252 if(sigma<0.) return 0.;
253 return sigma;
254 }
255
256 // Calculation formula for piMinus-nuclear ine
257 G4double G4ChipsKaonMinusInelasticXS::CrossSec
258 {
259 G4double lP=std::log(P);
260 return CrossSectionFormula(tZ, tN, P, lP);
261 }
262
263 // Calculation formula for piMinus-nuclear ine
264 G4double G4ChipsKaonMinusInelasticXS::CrossSec
265 {
266 G4double P=std::exp(lP);
267 return CrossSectionFormula(tZ, tN, P, lP);
268 }
269 // Calculation formula for piMinus-nuclear ine
270 G4double G4ChipsKaonMinusInelasticXS::CrossSec
271
272 {
273 G4double sigma=0.;
274 if(tZ==1 && !tN) // P
275 {
276 G4double ld=lP-3.5;
277 G4double ld2=ld*ld;
278 G4double p2=P*P;
279 G4double p4=p2*p2;
280 G4double sp=std::sqrt(P);
281 G4double psp=P*sp;
282 G4double lm=P-.39;
283 G4double md=lm*lm+.000156;
284 G4double lh=P-1.;
285 G4double hd=lh*lh+.0156;
286 G4double El=(.0557*ld2+2.23)/(1.-.7/sp+.07
287 G4double To=(.3*ld2+19.5)/(1.-.21/sp+.52/p
288 sigma=8.8/psp+(To-El)+.002/md+.15/hd;
289 }
290 else if(tZ==1 && tN==1) //
291 {
292 G4double p2=P*P;
293 G4double dX=lP-3.7;
294 G4double dR=P-.94;
295 G4double sp=std::sqrt(P);
296 sigma=(.6*dX*dX+36.)/(1.-.11/sp+.52/p2/p2)
297 }
298 else if(tZ<97 && tN<152) // G
299 {
300 G4double d=lP-4.2;
301 G4double sp=std::sqrt(P);
302 G4double p2=P*P;
303 G4double a=tN+tZ; //
304 G4double sa=std::sqrt(a);
305 G4double al=std::log(a);
306 G4double a2=a*a;
307 G4double c=52.*std::exp(al*0.6)*(1.+97./a2
308 G4double gg=-.2-.003*a;
309 G4double h=.5+.07*a;
310 G4double v=P-1.;
311 G4double f=.6*a*sa/(1.+.00002*a2);
312 G4double u=.125+.127*al;
313 sigma=(c+d*d)/(1.+gg/sp+h/p2/p2)+f/(v*v+u*
314 }
315 else
316 {
317 G4cerr<<"-Warning-G4ChipsKMinusNuclearCroS
318 sigma=0.;
319 }
320 if(sigma<0.) return 0.;
321 return sigma;
322 }
323
324 G4double G4ChipsKaonMinusInelasticXS::EquLinea
325 {
326 if(DX<=0. || N<2)
327 {
328 G4cerr<<"***G4ChipsKaonMinusInelasticXS:
329 return Y[0];
330 }
331
332 G4int N2=N-2;
333 G4double d=(X-X0)/DX;
334 G4int jj=static_cast<int>(d);
335 if (jj<0) jj=0;
336 else if(jj>N2) jj=N2;
337 d-=jj; // excess
338 G4double yi=Y[jj];
339 G4double sigma=yi+(Y[jj+1]-yi)*d;
340
341 return sigma;
342 }
343