Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4QFreeScattering Class Reference

#include <G4QFreeScattering.hh>

Public Member Functions

 ~G4QFreeScattering ()
 
std::pair< G4LorentzVector, G4LorentzVector > Scatter (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
 
G4QHadronVector * InElF (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
 
std::pair< G4double, G4double > GetElTotMean (G4double pIU, G4int hPDG, G4int Z, G4int N)
 
std::pair< G4double, G4double > GetElTotXS (G4double Mom, G4int PDG, G4bool F)
 
std::pair< G4double, G4double > FetchElTot (G4double pGeV, G4int PDG, G4bool F)
 
std::pair< G4double, G4double > CalcElTot (G4double pGeV, G4int Index)
 

Static Public Member Functions

static G4QFreeScattering * GetPointer ()
 

Protected Member Functions

 G4QFreeScattering ()
 

Detailed Description

Definition at line 49 of file G4QFreeScattering.hh.

Constructor & Destructor Documentation

◆ G4QFreeScattering()

G4QFreeScattering::G4QFreeScattering ( )
protected

Definition at line 49 of file G4QFreeScattering.cc.

50{
51#ifdef pdebug
52 G4cout<<"***^^^*** G4QFreeScattering singletone is created ***^^^***"<<G4endl;
53#endif
54}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout

◆ ~G4QFreeScattering()

G4QFreeScattering::~G4QFreeScattering ( )

Definition at line 56 of file G4QFreeScattering.cc.

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}

Member Function Documentation

◆ CalcElTot()

std::pair< G4double, G4double > G4QFreeScattering::CalcElTot ( G4double  pGeV,
G4int  Index 
)

Definition at line 71 of file G4QFreeScattering.cc.

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
LE
@ LE
Definition: Evaluator.cc:66
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:64
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

Referenced by FetchElTot(), GetElTotXS(), and G4QuasiFreeRatios::GetElTotXS().

◆ FetchElTot()

std::pair< G4double, G4double > G4QFreeScattering::FetchElTot ( G4double  pGeV,
G4int  PDG,
G4bool  F 
)

Definition at line 373 of file G4QFreeScattering.cc.

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
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4QFreeScattering::CalcElTot
std::pair< G4double, G4double > CalcElTot(G4double pGeV, G4int Index)
Definition: G4QFreeScattering.cc:71

Referenced by G4QuasiFreeRatios::GetChExFactor(), G4QuasiFreeRatios::GetElTot(), and GetElTotMean().

◆ GetElTotMean()

std::pair< G4double, G4double > G4QFreeScattering::GetElTotMean ( G4double  pIU,
G4int  hPDG,
G4int  Z,
G4int  N 
)

Definition at line 538 of file G4QFreeScattering.cc.

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
G4QFreeScattering::FetchElTot
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
Definition: G4QFreeScattering.cc:373

◆ GetElTotXS()

std::pair< G4double, G4double > G4QFreeScattering::GetElTotXS ( G4double  Mom,
G4int  PDG,
G4bool  F 
)

Definition at line 340 of file G4QFreeScattering.cc.

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}

◆ GetPointer()

G4QFreeScattering * G4QFreeScattering::GetPointer ( )
static

Definition at line 64 of file G4QFreeScattering.cc.

65{
66 static G4QFreeScattering theQFS; // *** Static body of the Quasi-Free Scattering ***
67 return &theQFS;
68}

Referenced by G4QEnvironment::G4QEnvironment(), G4QuasiFreeRatios::G4QuasiFreeRatios(), and G4QEnvironment::operator=().

◆ InElF()

G4QHadronVector * G4QFreeScattering::InElF ( G4int  NPDG,
G4LorentzVector  N4M,
G4int  pPDG,
G4LorentzVector  p4M 
)

Definition at line 667 of file G4QFreeScattering.cc.

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
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76

◆ Scatter()

std::pair< G4LorentzVector, G4LorentzVector > G4QFreeScattering::Scatter ( G4int  NPDG,
G4LorentzVector  N4M,
G4int  pPDG,
G4LorentzVector  p4M 
)

Definition at line 565 of file G4QFreeScattering.cc.

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
G4cerr
G4DLLIMPORT std::ostream G4cerr
The documentation for this class was generated from the following files: