Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4QFreeScattering.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: G4QFreeScattering for quasi-free scattering
31// Created: M.V. Kossov, CERN/ITEP(Moscow), 26-OCT-11
32// The last update: M.V. Kossov, CERN/ITEP (Moscow) 29-Oct-11
33//
34// ----------------------------------------------------------------------
35// Short description: Provides quasi-free scattering on nuclear nucleons.
36// ----------------------------------------------------------------------
37
38//#define debug
39//#define pdebug
40//#define ppdebug
41//#define nandebug
42
43#include "G4QFreeScattering.hh"
44#include "G4SystemOfUnits.hh"
45
46// initialisation of statics
47std::vector<std::pair<G4double,G4double>*> G4QFreeScattering::vX; // ETPointers to LogTable
48
49G4QFreeScattering::G4QFreeScattering()
50{
51#ifdef pdebug
52 G4cout<<"***^^^*** G4QFreeScattering singletone is created ***^^^***"<<G4endl;
53#endif
54}
55
56G4QFreeScattering::~G4QFreeScattering()
57{
58 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
59 for(pos2=vX.begin(); pos2<vX.end(); pos2++) delete [] * pos2;
60 vX.clear();
61}
62
63// Returns Pointer to the G4VQCrossSection class
64G4QFreeScattering* G4QFreeScattering::GetPointer()
65{
66 static G4QFreeScattering theQFS; // *** Static body of the Quasi-Free Scattering ***
67 return &theQFS;
68}
69
70// Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot)
71std::pair<G4double,G4double> G4QFreeScattering::CalcElTot(G4double p, G4int I)
72{
73 // ---------> Each parameter set can have not more than nPoints=128 parameters
74 static const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola
75 static const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient
76 static const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient
77 static const G4double pmi=.1; // Below that fast LE calculation is made
78 static const G4double pma=1000.; // Above that fast HE calculation is made
79 G4double El=0.; // prototype of the elastic hN cross-section
80 G4double To=0.; // prototype of the total hN cross-section
81 if(p<=0.)
82 {
83 //G4cout<<"-Warning-G4QFreeScattering::CalcElTot: p="<<p<<" is 0 or negative"<<G4endl;
84 return std::make_pair(El,To);
85 }
86 if (!I) // pp/nn
87 {
88#ifdef debug
89 G4cout<<"G4QFreeScatter::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl;
90#endif
91 if(p<pmi)
92 {
93 G4double p2=p*p;
94 El=1./(.00012+p2*.2);
95 To=El;
96#ifdef debug
97 G4cout<<"G4QFreeScatter::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl;
98#endif
99 }
100 else if(p>pma)
101 {
102 G4double lp=std::log(p)-lmi;
103 G4double lp2=lp*lp;
104 El=pbe*lp2+6.72;
105 To=pbt*lp2+38.2;
106#ifdef debug
107 G4cout<<"G4QFreeScat::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl;
108#endif
109 }
110 else
111 {
112 G4double p2=p*p;
113 G4double LE=1./(.00012+p2*.2);
114 G4double lp=std::log(p)-lmi;
115 G4double lp2=lp*lp;
116 G4double rp2=1./p2;
117 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
118 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
119#ifdef debug
120 G4cout<<"G4QFreeScat::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl;
121#endif
122 }
123 }
124 else if(I==1) // np/pn
125 {
126 if(p<pmi)
127 {
128 G4double p2=p*p;
129 El=1./(.00012+p2*(.051+.1*p2));
130 To=El;
131 }
132 else if(p>pma)
133 {
134 G4double lp=std::log(p)-lmi;
135 G4double lp2=lp*lp;
136 El=pbe*lp2+6.72;
137 To=pbt*lp2+38.2;
138 }
139 else
140 {
141 G4double p2=p*p;
142 G4double LE=1./(.00012+p2*(.051+.1*p2));
143 G4double lp=std::log(p)-lmi;
144 G4double lp2=lp*lp;
145 G4double rp2=1./p2;
146 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
147 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
148 }
149 }
150 else if(I==2) // pimp/pipn
151 {
152 G4double lp=std::log(p);
153 if(p<pmi)
154 {
155 G4double lr=lp+1.27;
156 El=1.53/(lr*lr+.0676);
157 To=El*3;
158 }
159 else if(p>pma)
160 {
161 G4double ld=lp-lmi;
162 G4double ld2=ld*ld;
163 G4double sp=std::sqrt(p);
164 El=pbe*ld2+2.4+7./sp;
165 To=pbt*ld2+22.3+12./sp;
166 }
167 else
168 {
169 G4double lr=lp+1.27; // p1
170 G4double LE=1.53/(lr*lr+.0676); // p2, p3
171 G4double ld=lp-lmi; // p4 (lmi=3.5)
172 G4double ld2=ld*ld;
173 G4double p2=p*p;
174 G4double p4=p2*p2;
175 G4double sp=std::sqrt(p);
176 G4double lm=lp+.36; // p5
177 G4double md=lm*lm+.04; // p6
178 G4double lh=lp-.017; // p7
179 G4double hd=lh*lh+.0025; // p8
180 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
181 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
182 }
183 }
184 else if(I==3) // pipp/pimn
185 {
186 G4double lp=std::log(p);
187 if(p<pmi)
188 {
189 G4double lr=lp+1.27;
190 G4double lr2=lr*lr;
191 El=13./(lr2+lr2*lr2+.0676);
192 To=El;
193 }
194 else if(p>pma)
195 {
196 G4double ld=lp-lmi;
197 G4double ld2=ld*ld;
198 G4double sp=std::sqrt(p);
199 El=pbe*ld2+2.4+6./sp;
200 To=pbt*ld2+22.3+5./sp;
201 }
202 else
203 {
204 G4double lr=lp+1.27; // p1
205 G4double lr2=lr*lr;
206 G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3
207 G4double ld=lp-lmi; // p4 (lmi=3.5)
208 G4double ld2=ld*ld;
209 G4double p2=p*p;
210 G4double p4=p2*p2;
211 G4double sp=std::sqrt(p);
212 G4double lm=lp-.32; // p5
213 G4double md=lm*lm+.0576; // p6
214 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
215 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
216 }
217 }
218 else if(I==4) // Kmp/Kmn/K0p/K0n
219 {
220
221 if(p<pmi)
222 {
223 G4double psp=p*std::sqrt(p);
224 El=5.2/psp;
225 To=14./psp;
226 }
227 else if(p>pma)
228 {
229 G4double ld=std::log(p)-lmi;
230 G4double ld2=ld*ld;
231 El=pbe*ld2+2.23;
232 To=pbt*ld2+19.5;
233 }
234 else
235 {
236 G4double ld=std::log(p)-lmi;
237 G4double ld2=ld*ld;
238 G4double sp=std::sqrt(p);
239 G4double psp=p*sp;
240 G4double p2=p*p;
241 G4double p4=p2*p2;
242 G4double lm=p-.39;
243 G4double md=lm*lm+.000156;
244 G4double lh=p-1.;
245 G4double hd=lh*lh+.0156;
246 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
247 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
248 }
249 }
250 else if(I==5) // Kpp/Kpn/aKp/aKn
251 {
252 if(p<pmi)
253 {
254 G4double lr=p-.38;
255 G4double lm=p-1.;
256 G4double md=lm*lm+.372;
257 El=.7/(lr*lr+.0676)+2./md;
258 To=El+.6/md;
259 }
260 else if(p>pma)
261 {
262 G4double ld=std::log(p)-lmi;
263 G4double ld2=ld*ld;
264 El=pbe*ld2+2.23;
265 To=pbt*ld2+19.5;
266 }
267 else
268 {
269 G4double ld=std::log(p)-lmi;
270 G4double ld2=ld*ld;
271 G4double lr=p-.38;
272 G4double LE=.7/(lr*lr+.0676);
273 G4double sp=std::sqrt(p);
274 G4double p2=p*p;
275 G4double p4=p2*p2;
276 G4double lm=p-1.;
277 G4double md=lm*lm+.372;
278 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
279 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
280 }
281 }
282 else if(I==6) // hyperon-N
283 {
284 if(p<pmi)
285 {
286 G4double p2=p*p;
287 El=1./(.002+p2*(.12+p2));
288 To=El;
289 }
290 else if(p>pma)
291 {
292 G4double lp=std::log(p)-lmi;
293 G4double lp2=lp*lp;
294 G4double sp=std::sqrt(p);
295 El=(pbe*lp2+6.72)/(1.+2./sp);
296 To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
297 }
298 else
299 {
300 G4double p2=p*p;
301 G4double LE=1./(.002+p2*(.12+p2));
302 G4double lp=std::log(p)-lmi;
303 G4double lp2=lp*lp;
304 G4double p4=p2*p2;
305 G4double sp=std::sqrt(p);
306 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
307 To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
308 }
309 }
310 else if(I==7) // antibaryon-N
311 {
312 if(p>pma)
313 {
314 G4double lp=std::log(p)-lmi;
315 G4double lp2=lp*lp;
316 El=pbe*lp2+6.72;
317 To=pbt*lp2+38.2;
318 }
319 else
320 {
321 G4double ye=std::pow(p,1.25);
322 G4double yt=std::pow(p,.35);
323 G4double lp=std::log(p)-lmi;
324 G4double lp2=lp*lp;
325 El=80./(ye+1.)+pbe*lp2+6.72;
326 To=(80./yt+.3)/yt+pbt*lp2+38.2;
327 }
328 }
329 else
330 {
331 G4cout<<"*Error*G4QFreeScattering::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
332 G4Exception("G4QFreeScattering::CalcElTot:","23",FatalException,"CHIPScrash");
333 }
334 if(El>To) El=To;
335 return std::make_pair(El,To);
336} // End of CalcElTot
337
338// For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
339// *** Only for external use. For simulation it is better to use GetElTot(...) [AMDB] ***
340std::pair<G4double,G4double> G4QFreeScattering::GetElTotXS(G4double p, G4int PDG, G4bool F)
341{
342 G4int ind=0; // Prototype of the reaction index
343 G4bool kfl=true; // Flag of K0/aK0 oscillation
344 G4bool kf=false;
345 if(PDG==90001000) PDG=2212;
346 if(PDG==90000001) PDG=2112;
347 if(PDG==91000000) PDG=3122;
348 if(PDG==130 || PDG==310)
349 {
350 kf=true;
351 if(G4UniformRand()>.5) kfl=false;
352 }
353 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
354 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
355 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
356 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
357 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
358 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
359 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
360 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
361 else {
362 // VI changed Fatal to Warning, because this method cannot do better
363 // source of problem in classes which make this call
364 ind = 0;
365 G4cout<<"*Error*G4QFreeScattering::CalcElTotXS: PDG="<<PDG
366 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
367 G4Exception("G4QFreeScattering::CalcElTotXS:","22",JustWarning,"CHIPS_crash");
368 }
369 return CalcElTot(p,ind); // This is a slow direct method, better use FetchElTot
370}
371
372// Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron
373std::pair<G4double,G4double> G4QFreeScattering::FetchElTot(G4double p, G4int PDG, G4bool F)
374{
375 static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step)
376 static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable
377 static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
378 static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
379 static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
380 static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
381 static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
382 static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
383 //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum"
384 static G4double lastP=0.; // The last momentum for which XS was calculated
385 static G4int lastH=0; // The last projPDG for which XS was calculated
386 static G4bool lastF=true; // The last nucleon for which XS was calculated
387 static std::pair<G4double,G4double> lastR(0.,0.); // The last result
388 // Local Associative Data Base:
389 static std::vector<G4int> vI; // Vector of index for which XS was calculated
390 static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable
391 static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable
392 // Last values of the Associative Data Base:
393 static G4int lastI=0; // The Last index for which XS was calculated
394 static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable
395 static G4int lastK=0; // The Last topBin number initialized in LogTable
396 static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
397 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
398 G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB
399 if(PDG==90001000) PDG=2212;
400 if(PDG==90000001) PDG=2112;
401 if(PDG==91000000) PDG=3122;
402#ifdef pdebug
403 G4cout<<"G4QFreeScatter::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl;
404#endif
405 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
406 // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
407 lastH=PDG;
408 lastF=F;
409 G4int ind=-1; // Prototipe of the index of the PDG/F combination
410 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
411 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp)
412 G4bool kfl=true; // Flag of K0/aK0 oscillation
413 G4bool kf=false;
414 if(PDG==130 || PDG==310)
415 {
416 kf=true;
417 if(G4UniformRand()>.5) kfl=false;
418 }
419 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
420 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
421 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
422 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
423 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
424 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
425 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
426 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
427 else
428 {
429 // VI changed Fatal to Warning, because this method cannot do better
430 // source of problem in classes which make this call
431 ind = 0;
432 G4cout<<"*Error*G4QFreeScattering::FetchElTot: PDG="<<PDG
433 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
434 G4Exception("G4QFreeScattering::FetchElTot:","22",JustWarning,"CHIPS problem");
435 }
436 if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler
437 // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
438 if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?)
439 G4bool found=false;
440 G4int i=-1;
441 if(nDB) for (i=0; i<nDB; ++i) if(ind==vI[i]) // Sirch for this index in AMDB
442 {
443 found=true; // The index is found
444 break;
445 }
446 G4double lp=std::log(p);
447#ifdef pdebug
448 G4cout<<"G4QFreeScat::FetchElTot: I="<<ind<<", i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl;
449#endif
450 if(!nDB || !found) // Create new line in the AMDB
451 {
452#ifdef pdebug
453 G4cout<<"G4QFreeScattering::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl;
454#endif
455 lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
456 lastI = ind; // Remember the initialized inex
457 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB
458 if(lastK>nlp)
459 {
460 lastK=nlp;
461 lastM=lpa-lpi;
462 }
463 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
464 G4double pv=mi;
465 for(G4int j=0; j<=lastK; ++j) // Calculate LogTab values
466 {
467 lastX[j]=CalcElTot(pv,ind);
468#ifdef pdebug
469 G4cout<<"G4QFreeScat::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T="
470 <<lastX[j].second<<G4endl;
471#endif
472 if(j!=lastK) pv*=edl;
473 }
474 i++; // Make a new record to AMDB and position on it
475 vI.push_back(lastI);
476 vM.push_back(lastM);
477 vK.push_back(lastK);
478 vX.push_back(lastX);
479 }
480 else // The A value was found in AMDB
481 {
482 lastI=vI[i];
483 lastM=vM[i];
484 lastK=vK[i];
485 lastX=vX[i];
486 G4int nextK=lastK+1;
487 G4double lpM=lastM+lpi;
488#ifdef pdebug
489 G4cout<<"G4QFreeScatt::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl;
490#endif
491 if(lp>lpM && lastK<nlp) // LogTab must be updated
492 {
493 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
494#ifdef pdebug
495 G4cout<<"G4QFreeScat::FetET: K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl;
496#endif
497 if(lastK>nlp)
498 {
499 lastK=nlp;
500 lastM=lpa-lpi;
501 }
502 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
503 G4double pv=std::exp(lpM); // momentum of the last calculated beam
504 for(G4int j=nextK; j<=lastK; ++j) // Calculate LogTab values
505 {
506 pv*=edl;
507 lastX[j]=CalcElTot(pv,ind);
508#ifdef pdebug
509 G4cout<<"G4QFreeScat::FetchElTot: U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T="
510 <<lastX[j].second<<G4endl;
511#endif
512 }
513 } // End of LogTab update
514 if(lastK >= nextK) // The AMDB was apdated
515 {
516 vM[i]=lastM;
517 vK[i]=lastK;
518 }
519 }
520 // Now one can use tabeles to calculate the value
521 G4double dlp=lp-lpi; // Shifted log(p) value
522 G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin
523 G4double d=dlp-n*dl; // Log shift
524 G4double e=lastX[n].first; // E-Base
525 lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result
526 if(lastR.first<0.) lastR.first = 0.;
527 G4double t=lastX[n].second; // T-Base
528 lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
529 if(lastR.second<0.) lastR.second= 0.;
530#ifdef pdebug
531 G4cout<<"=O=>G4QFreeScat::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl;
532#endif
533 if(lastR.first>lastR.second) lastR.first = lastR.second;
534 return lastR;
535} // End of FetchElTot
536
537// (Mean Elastic and Mean Total over (Z,N)) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
538std::pair<G4double,G4double> G4QFreeScattering::GetElTotMean(G4double pIU, G4int hPDG,
539 G4int Z, G4int N)
540{
541 G4double pGeV=pIU/gigaelectronvolt;
542 if(hPDG==90001000) hPDG=2212;
543 if(hPDG==90000001) hPDG=2112;
544 if(hPDG==91000000) hPDG=3122;
545#ifdef pdebug
546 G4cout<<"->G4QFreeSc::GetElTotMean: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
547#endif
548 if(Z<1 && N<1)
549 {
550 G4cout<<"-Warning-G4QFreeScat::GetElTotMean: Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
551 return std::make_pair(0.,0.);
552 }
553 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
554 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
555#ifdef pdebug
556 G4cout<<"-OUT->G4QFreeScat::GetElTotMean: hp("<<hp.first<<","<<hp.second<<"), hn("
557 <<hn.first<<","<<hn.second<<")"<<G4endl;
558#endif
559 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
560 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
561} // End of GetElTotMean
562
563// scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
564// if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
565std::pair<G4LorentzVector,G4LorentzVector> G4QFreeScattering::Scatter(G4int NPDG,
566 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
567{
568 static const G4double mNeut= G4QPDGCode(2112).GetMass();
569 static const G4double mProt= G4QPDGCode(2212).GetMass();
570 p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units)
571 N4M/=megaelectronvolt;
572 G4LorentzVector tot4M=N4M+p4M;
573 if(pPDG==90001000) pPDG=2212;
574 if(pPDG==90000001) pPDG=2112;
575 if(pPDG==91000000) pPDG=3122;
576#ifdef ppdebug
577 G4cout<<"->G4QFR::Scat:p4M="<<p4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
578#endif
579 G4double mT=mNeut;
580#ifdef ppdebug
581 G4int Z=0;
582 G4int N=1;
583#endif
584 if(NPDG==2212 || NPDG==90001000)
585 {
586 mT=mProt;
587#ifdef ppdebug
588 Z=1;
589 N=0;
590#endif
591 }
592 else if(NPDG!=2112 && NPDG!=90000001)
593 {
594 G4cout<<"Error:G4QFreeScattering::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
595 G4Exception("G4QFreeScattering::Scatter:","21",FatalException,"CHIPScomplain");
596 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Use this if not exception
597 }
598 G4double mT2=mT*mT; // a squared mass of the free scattered nuclead cluster (FSNC)
599 G4double mP2=p4M.m2(); // a projectile squared mass
600 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); // a projectile energy in the CMS of FSNC
601#ifdef pdebug
602 G4cout<<"G4QFS::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
603#endif
604 G4double E2=E*E;
605 if( E < 0. || E2 < mP2)
606 {
607#ifdef ppdebug
608 G4cout<<"-Warning-G4QFS::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
609#endif
610 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
611 }
612 G4double pP2=E2-mP2; // Squared Momentum in pseudo laboratory system (final particles)
613 G4double P=std::sqrt(pP2); // Momentum in pseudo laboratory system (final particles) ?
614#ifdef ppdebug
615 G4cout<<"G4QFreeS::Scatter: Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
616#endif
617 // @@ Temporary NN t-dependence for all hadrons
618 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QFreeScat::Scatter: pPDG="<<pPDG<<G4endl;
619 // @@ check a possibility to separate p, n, or alpha (!)
620 G4double dmT=mT+mT;
621 G4double s_value=dmT*E2+mP2+mT2; // Mondelstam s (?)
622 G4double maxt=dmT*dmT*pP2/s_value; // max possible |t|
623 G4double mint=0.; // Prototype of functional rand -t (MeV^2)
624 if(P < 14.) mint=maxt*G4UniformRand(); // S-wave
625 else // Calculate slopes (no cash !)
626 {
627 G4double p4=pP2*pP2/1000000.; // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
628 G4double theB1=(7.2+4.32/p4/(p4+12./P))/(1.+2.5/p4); // p4 in GeV, P in MeV
629 mint=-std::log(G4UniformRand())/theB1; // t-chan only
630 }
631#ifdef ppdebug
632 G4cout<<"G4QFS::Scat:PDG="<<pPDG<<",P="<<P<<",-t="<<mint<<"<"<<maxt<<", Z="<<Z<<",N="<<N
633 <<G4endl;
634#endif
635#ifdef nandebug
636 if(mint>-.0000001);
637 else G4cout<<"*Warning*G4QFreeScattering::Scatter: -t="<<mint<<G4endl;
638#endif
639 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
640#ifdef ppdebug
641 G4cout<<"G4QFS::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
642#endif
643 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
644 {
645 if (cost > 1.) cost = 1.;
646 else if(cost <-1.) cost =-1.;
647 else
648 {
649 G4cerr<<"G4QFreeScatter::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<maxt<<G4endl;
650 return std::make_pair(G4LorentzVector(0.,0.,0.,0.), p4M); // Do Nothing Action
651 }
652 }
653 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
654 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
655 if(!G4QHadron(tot4M).RelDecayIn2(p4M, reco4M, dir4M, cost, cost))
656 {
657 G4cerr<<"G4QFS::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
658 //G4Exception("G4QFS::Scat:","009",FatalException,"Decay of ElasticComp");
659 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
660 }
661#ifdef ppdebug
662 G4cout<<"G4QFS::Scat:p4M="<<p4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
663#endif
664 return std::make_pair( reco4M*megaelectronvolt, p4M*megaelectronvolt ); // The Result
665} // End of Scatter
666
667G4QHadronVector* G4QFreeScattering::InElF(G4int NPDG, G4LorentzVector N4M,
668 G4int pPDG, G4LorentzVector p4M)
669{
670 static const G4double mPi0 = G4QPDGCode(111).GetMass();
671 static const G4double mPi = G4QPDGCode(211).GetMass(); // Pi+ (Pi-: -211)
672 static const G4double mK = G4QPDGCode(321).GetMass(); // K+ (K- : -321)
673 static const G4double mK0 = G4QPDGCode(311).GetMass(); // aK0 (K0 : -311)
674 static const G4double mNeut= G4QPDGCode(2112).GetMass();
675 static const G4double mProt= G4QPDGCode(2212).GetMass();
676 static const G4double mLamb= G4QPDGCode(3122).GetMass();
677 //static const G4double mSigZ= G4QPDGCode(3212).GetMass();
678 static const G4double mSigM= G4QPDGCode(3112).GetMass();
679 static const G4double mSigP= G4QPDGCode(3222).GetMass();
680 static const G4double mXiM = G4QPDGCode(3312).GetMass();
681 static const G4double mXiZ = G4QPDGCode(3322).GetMass();
682 //static const G4double mOmM = G4QPDGCode(3334).GetMass();
683 static const G4double third =1./3.;
684 static const G4double twothd =2./3.;
685
686 p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units)
687 N4M/=megaelectronvolt;
688 G4LorentzVector c4M=N4M+p4M;
689 if(pPDG==90001000) pPDG=2212;
690 if(pPDG==90000001) pPDG=2112;
691 if(pPDG==91000000) pPDG=3122;
692#ifdef ppdebug
693 G4cout<<"->G4QFR::InElF: p4M="<<p4M<<",N4M="<<N4M<<", c4M="<<c4M<<",NPDG="<<NPDG<<G4endl;
694#endif
695 G4double mT=mNeut;
696 G4int Z=0;
697 //G4int N=1;
698 if(NPDG==2212 || NPDG==90001000)
699 {
700 NPDG=2212;
701 mT=mProt;
702 Z=1;
703 //N=0;
704 }
705 else if(NPDG!=2112 && NPDG!=90000001)
706 {
707 G4cout<<"Error:G4QFreeScattering::InElF:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
708 G4Exception("G4QFreeScattering::InElF:","21",FatalException,"CHIPS complain");
709 }
710 else NPDG=2112;
711 G4double mC2=c4M.m2();
712 G4double mC=0.;
713 if( mC2 < -0.0000001 )
714 {
715#ifdef ppdebug
716 G4cout<<"-Warning-G4QFS::InElF: Negative compoundMass ="<<mC2<<", c4M="<<c4M<<G4endl;
717#endif
718 return 0; // Do Nothing Action
719 }
720 else if ( mC2 > 0.0000001) mC=std::sqrt(mC2);
721#ifdef ppdebug
722 G4cout<<"G4QFS::InElF: mC ="<<mC<<G4endl;
723#endif
724 G4double mP=0.;
725 G4double mP2=p4M.m2(); // a projectile squared mass
726 if( mP2 < -0.0000001 )
727 {
728#ifdef ppdebug
729 G4cout<<"-Warning-G4QFS::InElF: Negative projectileMass ="<<mC2<<", c4M="<<c4M<<G4endl;
730#endif
731 return 0; // Do Nothing Action
732 }
733 else if ( mP2 > 0.0000001) mP=std::sqrt(mP2);
734#ifdef pdebug
735 G4cout<<"G4QFS::InElF:mT("<<mT<<")+mP("<<mP<<")+mPi0="<<mP+mT+mPi0<<"<? mC="<<mC<<G4endl;
736#endif
737 if(pPDG > 3334 || pPDG < -321) // Annihilation/Inelastic of anti-barions not implemented
738 {
739 G4cout<<"-Warning-G4QFreeScat::InElF: pPDG="<<pPDG<<G4endl;
740 return 0; // Do Nothing Action
741 }
742 if(pPDG==130 || pPDG==310)
743 {
744 if(G4UniformRand()>.5) pPDG = 311; // K0
745 else pPDG =-311; // aK0
746 }
747 G4int mPDG=111; // Additional newtral pion by default
748 G4double mM=mPi0; // Default Pi0 mass
749 G4double r=G4UniformRand(); // A random number to split pi0/pi
750 if (pPDG == 2212) // p-N
751 {
752 if(Z) // p-p
753 {
754 if(r<twothd) // -> n + Pi+ + p
755 {
756 pPDG=2112; // n
757 mP=mNeut;
758 mPDG= 211; // Pi+
759 mM=mPi;
760 }
761 }
762 else // p-n
763 {
764 if(r<third) // -> n + Pi+ + n
765 {
766 pPDG=2112; // n
767 mP=mNeut;
768 mPDG= 211; // Pi+
769 mM=mPi;
770 }
771 else if(r<twothd) // -> p + Pi- + p
772 {
773 mPDG=-211; // Pi-
774 mM=mPi;
775 NPDG=2212; // p
776 mT=mProt;
777 }
778 }
779 }
780 else if (pPDG == 2112) // n-N
781 {
782 if(Z) // n-p
783 {
784 if(r<third) // -> n + Pi+ + n
785 {
786 mPDG= 211; // Pi+
787 mM=mPi;
788 NPDG=2112; // n
789 mT=mNeut;
790 }
791 else if(r<twothd) // -> p + Pi- + p
792 {
793 pPDG=2212; // p
794 mP=mProt;
795 mPDG=-211; // Pi-
796 mM=mPi;
797 }
798 }
799 else // n-n
800 {
801 if(r<twothd) // -> p + Pi- + n
802 {
803 pPDG=2212; // p
804 mP=mProt;
805 mPDG= -211; // Pi-
806 mM=mPi;
807 }
808 }
809 }
810 else if (pPDG == 211) // pip-N
811 {
812 if(Z) // pip-p
813 {
814 if(r<0.5) // -> Pi+ + Pi+ + n
815 {
816 mPDG= 211; // Pi+
817 mM=mPi;
818 NPDG=2112; // n
819 mT=mNeut;
820 }
821 }
822 else // pip-n
823 {
824 if(r<third) // -> Pi0 + Pi0 + p
825 {
826 pPDG= 111; // Pi0
827 mP=mPi0;
828 NPDG=2212; // p
829 mT=mProt;
830 }
831 else if(r<twothd) // -> Pi+ + Pi- + p
832 {
833 mPDG=-211; // Pi-
834 mM=mPi;
835 NPDG=2212; // p
836 mT=mProt;
837 }
838 }
839 }
840 else if (pPDG == -211) // pim-N
841 {
842 if(Z) // pim-p
843 {
844 if(r<third) // -> Pi0 + Pi0 + n
845 {
846 pPDG= 111; // Pi0
847 mP=mPi0;
848 NPDG=2112; // n
849 mT=mNeut;
850 }
851 else if(r<twothd) // -> Pi- + Pi+ + n
852 {
853 mPDG= 211; // Pi+
854 mM=mPi;
855 NPDG=2112; // n
856 mT=mNeut;
857 }
858 }
859 else // pim-n
860 {
861 if(r<0.5) // -> Pi- + Pi- + p
862 {
863 mPDG=-211; // Pi-
864 mM=mPi;
865 NPDG=2212; // p
866 mT=mProt;
867 }
868 }
869 }
870 else if (pPDG == -321) // Km-N
871 {
872 if(Z) // Km-p
873 {
874 if(r<0.25) // K0 + Pi0 + n
875 {
876 pPDG=-311; // K0
877 mP=mK0;
878 NPDG=2112; // n
879 mT=mNeut;
880 }
881 else if(r<0.5) // -> K- + Pi+ + n
882 {
883 mPDG= 211; // Pi+
884 mM=mPi;
885 NPDG=2112; // n
886 mT=mNeut;
887 }
888 else if(r<0.75) // -> K0 + Pi- + p
889 {
890 pPDG=-311; // K0
891 mP=mK0;
892 mPDG=-211; // Pi-
893 mM=mPi;
894 }
895 }
896 else // Km-n
897 {
898 if(r<third) // -> K- + Pi- + p
899 {
900 mPDG=-211; // Pi-
901 mM=mPi;
902 NPDG=2212; // p
903 mT=mProt;
904 }
905 else if(r<twothd) // -> K0 + Pi- + n
906 {
907 pPDG=-311; // K0
908 mP=mK0;
909 mPDG=-211; // Pi-
910 mM=mPi;
911 }
912 }
913 }
914 else if (pPDG == -311) // K0-N
915 {
916 if(Z) // K0-p
917 {
918 if(r<third) // K- + Pi+ + p
919 {
920 pPDG=-321; // K-
921 mP=mK;
922 NPDG=2212; // p
923 mT=mProt;
924 }
925 else if(r<twothd) // -> K0 + Pi+ + n
926 {
927 mPDG= 211; // Pi+
928 mM=mPi;
929 NPDG=2112; // n
930 mT=mNeut;
931 }
932 }
933 else // K0-n
934 {
935 if(r<0.25) // -> K- + Pi+ + n
936 {
937 pPDG=-321; // K-
938 mP=mK;
939 mPDG= 211; // Pi+
940 mM=mPi;
941 }
942 else if(r<0.5) // -> K- + Pi0 + p
943 {
944 pPDG=-321; // K-
945 mP=mK;
946 NPDG=2212; // p
947 mT=mProt;
948 }
949 else if(r<0.75) // -> K0 + Pi- + p
950 {
951 mPDG=-211; // Pi-
952 mM=mPi;
953 NPDG=2212; // p
954 mT=mProt;
955 }
956 }
957 }
958 else if (pPDG == 321) // Kp-N
959 {
960 if(Z) // Kp-p
961 {
962 if(r<third) // -> K+ + Pi+ + n
963 {
964 mPDG= 211; // Pi+
965 mM=mPi;
966 NPDG=2112; // n
967 mT=mNeut;
968 }
969 else if(r<twothd) // -> aK0 + Pi+ + p
970 {
971 pPDG= 311; // aK0
972 mP=mK0;
973 mPDG= 211; // Pi+
974 mM=mPi;
975 }
976 }
977 else // Kp-n
978 {
979 if(r<0.25) // -> aK0 + Pi+ + n
980 {
981 pPDG= 311; // aK0
982 mP=mK0;
983 mPDG= 211; // Pi+
984 mM=mPi;
985 }
986 else if(r<0.5) // -> Kp + Pi- + p
987 {
988 mPDG=-211; // Pi-
989 mM=mPi;
990 NPDG=2212; // p
991 mT=mProt;
992 }
993 else if(r<0.75) // -> aK0 + Pi0 + p
994 {
995 pPDG= 311; // aK0
996 mP=mK0;
997 NPDG=2212; // p
998 mT=mProt;
999 }
1000 }
1001 }
1002 else if (pPDG == 311) // aK0-N
1003 {
1004 if(Z) // aK0-p
1005 {
1006 if(r<0.25) // -> K+ + Pi- + p
1007 {
1008 pPDG= 321; // K+
1009 mP=mK;
1010 mPDG=-211; // Pi-
1011 mM=mPi;
1012 }
1013 else if(r<0.5) // -> aK0 + Pi+ + n
1014 {
1015 mPDG= 211; // Pi+
1016 mM=mPi;
1017 NPDG=2112; // n
1018 mT=mNeut;
1019 }
1020 else if(r<0.75) // -> K+ + Pi0 + n
1021 {
1022 pPDG= 321; // K+
1023 mP=mK;
1024 NPDG=2112; // n
1025 mT=mNeut;
1026 }
1027 }
1028 else // aK0-n
1029 {
1030 if(r<third) // -> aK0 + Pi- + p
1031 {
1032 mPDG=-211; // Pi-
1033 mM=mPi;
1034 NPDG=2212; // p
1035 mT=mProt;
1036 }
1037 else if(r<twothd) // -> K+ + Pi- + n
1038 {
1039 pPDG= 321; // K+
1040 mP=mK;
1041 mPDG=-211; // Pi-
1042 mM=mPi;
1043 }
1044 }
1045 }
1046 else if (pPDG == 3122 || pPDG== 3212) // Lambda/Sigma0-N
1047 {
1048 if(pPDG == 3212)
1049 {
1050 pPDG=3122;
1051 mP=mLamb;
1052 }
1053 if(Z) // L/S0-p
1054 {
1055 if(r<0.2) // -> SigP + Pi0 + n
1056 {
1057 pPDG=3222; // SigP
1058 mP=mSigP;
1059 NPDG=2112; // n
1060 mT=mNeut;
1061 }
1062 else if(r<0.4) // -> SigP + Pi- + p
1063 {
1064 pPDG=3222; // SigP
1065 mP=mSigP;
1066 mPDG=-211; // Pi-
1067 mM=mPi;
1068 }
1069 else if(r<0.6) // -> SigM + Pi+ + p
1070 {
1071 pPDG=3112; // SigM
1072 mP=mSigM;
1073 mPDG= 211; // Pi+
1074 mM=mPi;
1075 }
1076 else if(r<0.8) // -> Lamb + Pi+ + n
1077 {
1078 mPDG= 211; // Pi+
1079 mM=mPi;
1080 NPDG=2112; // n
1081 mT=mNeut;
1082 }
1083 }
1084 else // L/S0-n
1085 {
1086 if(r<0.2) // -> SigM + Pi0 + p
1087 {
1088 pPDG=3112; // SigM
1089 mP=mSigM;
1090 NPDG=2212; // p
1091 mT=mProt;
1092 }
1093 else if(r<0.4) // -> SigP + Pi- + n
1094 {
1095 pPDG=3222; // SigP
1096 mP=mSigP;
1097 mPDG=-211; // Pi-
1098 mM=mPi;
1099 }
1100 else if(r<0.6) // -> SigM + Pi+ + n
1101 {
1102 pPDG=3112; // SigM
1103 mP=mSigM;
1104 mPDG= 211; // Pi+
1105 mM=mPi;
1106 }
1107 else if(r<0.8) // -> Lamb + Pi- + p
1108 {
1109 mPDG=-211; // Pi-
1110 mM=mPi;
1111 NPDG=2212; // p
1112 mT=mProt;
1113 }
1114 }
1115 }
1116 else if (pPDG == 3112) // Sigma- -N
1117 {
1118 if(Z) // SigM-p
1119 {
1120 if(r<0.25) // -> Lamb + Pi0 + n
1121 {
1122 pPDG=3122; // Lamb
1123 mP=mLamb;
1124 NPDG=2112; // n
1125 mT=mNeut;
1126 }
1127 else if(r<0.5) // -> Lamb + Pi- + p
1128 {
1129 pPDG=3122; // Lamb
1130 mP=mLamb;
1131 mPDG=-211; // Pi-
1132 mM=mPi;
1133 }
1134 else if(r<0.75) // -> SigM + Pi+ + n
1135 {
1136 mPDG= 211; // Pi+
1137 mM=mPi;
1138 NPDG=2112; // n
1139 mT=mNeut;
1140 }
1141 }
1142 else // SigM-n
1143 {
1144 if(r<third) // -> Lamb + Pi- + n
1145 {
1146 pPDG=3122; // Lamb
1147 mP=mLamb;
1148 mPDG=-211; // Pi-
1149 mM=mPi;
1150 }
1151 else if(r<twothd) // -> SigM + Pi- + p
1152 {
1153 mPDG=-211; // Pi-
1154 mM=mPi;
1155 NPDG=2212; // p
1156 mT=mProt;
1157 }
1158 }
1159 }
1160 else if (pPDG == 3222) // Sigma+ -N
1161 {
1162 if(Z) // SigP-p
1163 {
1164 if(r<third) // -> Lamb + Pi+ + p
1165 {
1166 pPDG=3122; // Lamb
1167 mP=mLamb;
1168 mPDG= 211; // Pi+
1169 mM=mPi;
1170 }
1171 else if(r<twothd) // -> SigP + Pi+ + n
1172 {
1173 mPDG= 211; // Pi+
1174 mM=mPi;
1175 NPDG=2112; // n
1176 mT=mNeut;
1177 }
1178 }
1179 else // SigP-n
1180 {
1181 if(r<0.25) // -> Lamb + Pi0 + p
1182 {
1183 pPDG=3122; // Lamb
1184 mP=mLamb;
1185 NPDG=2212; // p
1186 mT=mProt;
1187 }
1188 else if(r<0.5) // -> Lamb + Pi+ + n
1189 {
1190 pPDG=3122; // Lamb
1191 mP=mLamb;
1192 mPDG= 211; // Pi+
1193 mM=mPi;
1194 }
1195 else if(r<0.75) // -> SigP + Pi- + p
1196 {
1197 mPDG=-211; // Pi-
1198 mM=mPi;
1199 NPDG=2212; // p
1200 mT=mProt;
1201 }
1202 }
1203 }
1204 else if (pPDG == 3312) // Xi- -N
1205 {
1206 if(Z) // XiM-p
1207 {
1208 if(r<0.25) // -> Xi0 + Pi0 + n
1209 {
1210 pPDG=3322; // Xi0
1211 mP=mXiZ;
1212 NPDG=2112; // n
1213 mT=mNeut;
1214 }
1215 else if(r<0.5) // -> Xi0 + Pi- + p
1216 {
1217 pPDG=3322; // Xi0
1218 mP=mXiZ;
1219 mPDG=-211; // Pi-
1220 mM=mPi;
1221 }
1222 else if(r<0.75) // -> XiM + Pi+ + n
1223 {
1224 mPDG= 211; // Pi+
1225 mM=mPi;
1226 NPDG=2112; // n
1227 mT=mNeut;
1228 }
1229 }
1230 else // XiM-n
1231 {
1232 if(r<third) // -> Xi0 + Pi- + n
1233 {
1234 pPDG=3322; // Xi0
1235 mP=mXiZ;
1236 mPDG=-211; // Pi-
1237 mM=mPi;
1238 }
1239 else if(r<twothd) // -> XiM + Pi- + p
1240 {
1241 mPDG=-211; // Pi-
1242 mM=mPi;
1243 NPDG=2212; // p
1244 mT=mProt;
1245 }
1246 }
1247 }
1248 else if (pPDG == 3322) // Xi0 -N
1249 {
1250 if(Z) // Xi0-p
1251 {
1252 if(r<third) // -> Xi- + Pi+ + p
1253 {
1254 pPDG=3312; // Xi-
1255 mP=mXiM;
1256 mPDG= 211; // Pi+
1257 mM=mPi;
1258 }
1259 else if(r<twothd) // -> Xi0 + Pi+ + n
1260 {
1261 mPDG= 211; // Pi+
1262 mM=mPi;
1263 NPDG=2112; // n
1264 mT=mNeut;
1265 }
1266 }
1267 else // Xi0-n
1268 {
1269 if(r<0.25) // -> Xi- + Pi0 + p
1270 {
1271 pPDG=3312; // Xi-
1272 mP=mXiM;
1273 NPDG=2212; // p
1274 mT=mProt;
1275 }
1276 else if(r<0.5) // -> Xi- + Pi+ + n
1277 {
1278 pPDG=3312; // Xi-
1279 mP=mXiM;
1280 mPDG= 211; // Pi+
1281 mM=mPi;
1282 }
1283 else if(r<0.75) // -> Xi0 + Pi- + p
1284 {
1285 mPDG=-211; // Pi-
1286 mM=mPi;
1287 NPDG=2212; // p
1288 mT=mProt;
1289 }
1290 }
1291 }
1292 else if (pPDG == 3334) // Om- -N
1293 {
1294 if(Z) // OmM-p
1295 {
1296 if(r<0.5) // -> Om- + Pi+ + n
1297 {
1298 mPDG= 211; // Pi+
1299 mM=mPi;
1300 NPDG=2112; // n
1301 mT=mNeut;
1302 }
1303 }
1304 else // OmM-n
1305 {
1306 if(r<0.5) // -> Om- + Pi- + p
1307 {
1308 mPDG=-211; // Pi-
1309 mM=mPi;
1310 NPDG=2212; // p
1311 mT=mProt;
1312 }
1313 }
1314 }
1315 else
1316 {
1317 G4cout<<"*Error*G4QFreeScattering::InElF: PDG="<<pPDG
1318 <<", while it is defined only for p,n,hyperons(not Omega),pi,K/antiK"<<G4endl;
1319 G4Exception("G4QFreeScattering::InElF:","22",FatalException,"CHIPS_crash");
1320 }
1321 if (mC-mP-mM-mT <-0.000001) return 0;
1322 G4QHadronVector* TripQH = new G4QHadronVector; // Proto of the Result
1323 G4LorentzVector m4M(0.,0.,0.,mM);
1324 if(mC-mP-mM-mT < 0.000001) // Equal share
1325 {
1326 p4M=(mP/mC)*c4M;
1327 m4M=(mM/mC)*c4M;
1328 N4M=(mT/mC)*c4M;
1329 }
1330 else
1331 {
1332 p4M=G4LorentzVector(0.,0.,0.,mP);
1333 N4M=G4LorentzVector(0.,0.,0.,mT);
1334 if(!G4QHadron(c4M).DecayIn3(p4M,m4M,N4M))
1335 {
1336 G4ExceptionDescription ed;
1337 ed << "DecayIn3, TotM=" << mC /* << " <? " << mT + mP + mM */ << G4endl;
1338 G4Exception("G4QFreeScattering::InElF()","HAD_CHPS_0027", FatalException, ed);
1339 }
1340 G4QHadron* h1 = new G4QHadron(pPDG, p4M);
1341 TripQH->push_back(h1); // (delete equivalent, responsibility of users)
1342#ifdef debug
1343 G4cout << "G4QFreeScat::InElF: H1=" << pPDG << p4M << G4endl;
1344#endif
1345 G4QHadron* h2 = new G4QHadron(mPDG, m4M);
1346 TripQH->push_back(h2); // (delete equivalent, responsibility of users)
1347#ifdef debug
1348 G4cout << "G4QFreeScat::InElF: H2=" << mPDG << m4M << G4endl;
1349#endif
1350 G4QHadron* h3 = new G4QHadron(NPDG, N4M);
1351 TripQH->push_back(h3); // (delete equivalent, responsibility of users)
1352#ifdef debug
1353 G4cout << "G4QFreeScat::InElF: H3=" << NPDG << N4M << G4endl;
1354#endif
1355 }
1356 return TripQH; // The Result
1357} // End of Scatter
LE
@ LE
Definition: Evaluator.cc:66
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4QFreeScattering::FetchElTot
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
Definition: G4QFreeScattering.cc:373
G4QFreeScattering::CalcElTot
std::pair< G4double, G4double > CalcElTot(G4double pGeV, G4int Index)
Definition: G4QFreeScattering.cc:71
G4QFreeScattering::GetElTotXS
std::pair< G4double, G4double > GetElTotXS(G4double Mom, G4int PDG, G4bool F)
Definition: G4QFreeScattering.cc:340
G4QFreeScattering::GetPointer
static G4QFreeScattering * GetPointer()
Definition: G4QFreeScattering.cc:64
G4QFreeScattering::GetElTotMean
std::pair< G4double, G4double > GetElTotMean(G4double pIU, G4int hPDG, G4int Z, G4int N)
Definition: G4QFreeScattering.cc:538
G4QFreeScattering::InElF
G4QHadronVector * InElF(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
Definition: G4QFreeScattering.cc:667
G4QFreeScattering::Scatter
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
Definition: G4QFreeScattering.cc:565
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76