| 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 //
28 //
29 // G4 Physics class: G4QuasiElRatios for N+A e
30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10
31 // The last update: M.V. Kossov, CERN/ITEP (Mo
32 //
33 // -------------------------------------------
34 // This class has been extracted from the CHIP
35 // All the dependencies on CHIPS classes have
36 // Short description: Provides percentage of q
37 // reactions in the inelastic reactions.
38 // -------------------------------------------
39
40
41 #include "G4QuasiElRatios.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Proton.hh"
45 #include "G4Neutron.hh"
46 #include "G4Deuteron.hh"
47 #include "G4Triton.hh"
48 #include "G4He3.hh"
49 #include "G4Alpha.hh"
50 #include "G4ThreeVector.hh"
51 #include "G4CrossSectionDataSetRegistry.hh"
52 #include "G4Pow.hh"
53 #include "G4Log.hh"
54 #include "G4Exp.hh"
55
56 namespace {
57 const G4int nps=150; // Number o
58 const G4int mps=nps+1; // Number o
59 const G4double sma=150.; // The firs
60 const G4double ds=sma/nps; // Step of
61 const G4int nls=100; // Number o
62 const G4int mls=nls+1; // Number o
63 const G4double lsi=5.; // The min
64 const G4double lsa=9.; // The max
65 const G4double mi=G4Exp(lsi);// The min s
66 const G4double min_s=G4Exp(lsa);// The max
67 const G4double dls=(lsa-lsi)/nls;// Step o
68 const G4double edls=G4Exp(dls);// Multipli
69 const G4double toler=.01; // The tola
70 const G4double C=1.246;
71 const G4double lmi=3.5; // min of (l
72 const G4double pbe=.0557; // elastic (
73 const G4double pbt=.3; // total (ln
74 const G4double pmi=.1; // Below tha
75 const G4double pma=1000.; // Above tha
76 const G4int nlp=300; // Number
77 const G4int mlp=nlp+1; // Number
78 const G4double lpi=-5.; // The min
79 const G4double lpa=10.; // The max
80 const G4double mip=G4Exp(lpi);// The min p
81 const G4double map=G4Exp(lpa);// The max p
82 const G4double dlp=(lpa-lpi)/nlp;// Step o
83 const G4double edlp=G4Exp(dlp);// Multipli
84 }
85
86
87 G4QuasiElRatios::G4QuasiElRatios()
88 {
89 vT = new std::vector<G4double*>;
90 vL = new std::vector<G4double*>;
91 vX = new std::vector<std::pair<G4double,G4
92
93 lastSRatio=0.; // The last sig
94 lastRRatio=0.; // The last rat
95 lastARatio=0; // theLast of ca
96 lastHRatio=0.; // theLast of ma
97 lastNRatio=0; // theLast of to
98 lastMRatio=0.; // theLast of re
99 lastKRatio=0; // theLast of to
100 lastTRatio=0; // theLast of po
101 lastLRatio=0; // theLast of po
102 lastPtot=0.; // The last mome
103 lastHtot=0; // The last proj
104 lastFtot=true; // The last nucl
105 lastItot=0; // The Last index
106 lastMtot=0.; // The Last rel m
107 lastKtot=0; // The Last topBin
108 lastXtot = nullptr; // The Last ETPoin
109
110 PCSmanager=(G4ChipsProtonElasticXS*)G4Cros
111
112 NCSmanager=(G4ChipsNeutronElasticXS*)G4Cro
113 }
114
115 G4QuasiElRatios::~G4QuasiElRatios()
116 {
117 std::vector<G4double*>::iterator pos;
118 for(pos=vT->begin(); pos<vT->end(); pos++)
119 { delete [] *pos; }
120 vT->clear();
121 delete vT;
122
123 for(pos=vL->begin(); pos<vL->end(); pos++)
124 { delete [] *pos; }
125 vL->clear();
126 delete vL;
127
128 std::vector<std::pair<G4double,G4double>*>
129 for(pos2=vX->begin(); pos2<vX->end(); pos2
130 { delete [] *pos2; }
131 vX->clear();
132 delete vX;
133 }
134
135 // Calculation of pair(QuasiFree/Inelastic,Qua
136 std::pair<G4double,G4double> G4QuasiElRatios::
137
138 {
139 G4double R=0.;
140 G4double QF2In=1.;
141 G4int tgA=tgZ+tgN;
142 if(tgA<2) return std::make_pair(QF2In,R);
143 std::pair<G4double,G4double> ElTot=GetElTo
144 //if( ( (pPDG>999 && pIU<227.) || pIU<27.)
145 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.
146 else if(ElTot.second>0.)
147 {
148 R=ElTot.first/ElTot.second;
149 QF2In=GetQF2IN_Ratio(ElTot.second/mill
150 }
151 return std::make_pair(QF2In,R);
152 }
153
154 // Calculatio QasiFree/Inelastic Ratio as a fu
155 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4dou
156 {
157 // LogTable is created only if necessary.
158 if(m_s<toler || A<2) return 1.;
159 if(m_s>min_s) return 0.;
160
161 //if(A>238)
162 //{
163 // G4cout<<"-Warning-G4QuasiElRatio::Ge
164 // return 0.;
165 //}
166 G4int nDB=(G4int)vARatio.size(); // A num
167 if(nDB && lastARatio==A && m_s==lastSRatio
168 G4bool found=false;
169 G4int i=-1;
170 if(nDB) for (i=0; i<nDB; i++) if(A==vARati
171 {
172 found=true; //
173 break;
174 }
175 if(!nDB || !found) // C
176 {
177 lastARatio = A;
178 lastTRatio = new G4double[mps];
179 lastNRatio = static_cast<int>(m_s/ds)+
180 if(lastNRatio>nps)
181 {
182 lastNRatio=nps;
183 lastHRatio=sma;
184 }
185 else lastHRatio = lastNRatio*ds;
186 G4double sv=0;
187 lastTRatio[0]=1.;
188 for(G4int j=1; j<=lastNRatio; j++)
189 {
190 sv+=ds;
191 lastTRatio[j]=CalcQF2IN_Ratio(sv,A
192 }
193 lastLRatio=new G4double[mls];
194 // Initialize the logarithmic Table
195 for(G4int j=0; j<mls; ++j) lastLRatio[j]=0.0
196 if(m_s>sma)
197 {
198 G4double ls=G4Log(m_s);
199 lastKRatio = static_cast<int>((ls-
200 if(lastKRatio>nls)
201 {
202 lastKRatio=nls;
203 lastMRatio=lsa-lsi;
204 }
205 else lastMRatio = lastKRatio*dls;
206 sv=mi;
207 for(G4int j=0; j<=lastKRatio; j++)
208 {
209 lastLRatio[j]=CalcQF2IN_Ratio(
210 if(j!=lastKRatio) sv*=edls;
211 }
212 }
213 else //
214 {
215 lastKRatio = 0;
216 lastMRatio = 0.;
217 }
218 i++; //
219 vARatio.push_back(lastARatio);
220 vHRatio.push_back(lastHRatio);
221 vNRatio.push_back(lastNRatio);
222 vMRatio.push_back(lastMRatio);
223 vKRatio.push_back(lastKRatio);
224 vT->push_back(lastTRatio);
225 vL->push_back(lastLRatio);
226 }
227 else // T
228 {
229 lastARatio=vARatio[i];
230 lastHRatio=vHRatio[i];
231 lastNRatio=vNRatio[i];
232 lastMRatio=vMRatio[i];
233 lastKRatio=vKRatio[i];
234 lastTRatio=(*vT)[i];
235 lastLRatio=(*vL)[i];
236 if(m_s>lastHRatio)
237 {
238 G4int nextN=lastNRatio+1;
239 if(lastNRatio<nps)
240 {
241 G4double sv=lastHRatio; // bug fix by
242
243 lastNRatio = static_cast<int>(
244 if(lastNRatio>nps)
245 {
246 lastNRatio=nps;
247 lastHRatio=sma;
248 }
249 else lastHRatio = lastNRatio*d
250
251 for(G4int j=nextN; j<=lastNRat
252 {
253 sv+=ds;
254 lastTRatio[j]=CalcQF2IN_Ra
255 }
256 } // End of LinTab update
257 if(lastNRatio>=nextN)
258 {
259 vHRatio[i]=lastHRatio;
260 vNRatio[i]=lastNRatio;
261 }
262 G4int nextK=lastKRatio+1;
263 if(!lastKRatio) nextK=0;
264 if(m_s>sma && lastKRatio<nls)
265 {
266 G4double sv=G4Exp(lastMRatio+l
267 G4double ls=G4Log(m_s);
268 lastKRatio = static_cast<int>(
269 if(lastKRatio>nls)
270 {
271 lastKRatio=nls;
272 lastMRatio=lsa-lsi;
273 }
274 else lastMRatio = lastKRatio*d
275 for(G4int j=nextK; j<=lastKRat
276 {
277 sv*=edls;
278 lastLRatio[j]=CalcQF2IN_Ra
279 }
280 } // End of LogTab update
281 if(lastKRatio>=nextK)
282 {
283 vMRatio[i]=lastMRatio;
284 vKRatio[i]=lastKRatio;
285 }
286 }
287 }
288 // Now one can use tabeles to calculate th
289 if(m_s<sma) //
290 {
291 G4int n=static_cast<int>(m_s/ds);
292 G4double d=m_s-n*ds;
293 G4double v=lastTRatio[n];
294 lastRRatio=v+d*(lastTRatio[n+1]-v)/ds;
295 }
296 else // U
297 {
298 G4double ls=G4Log(m_s)-lsi; //
299 G4int n=static_cast<int>(ls/dls); /
300 G4double d=ls-n*dls; /
301 G4double v=lastLRatio[n];
302 lastRRatio=v+d*(lastLRatio[n+1]-v)/dls
303 }
304 if(lastRRatio<0.) lastRRatio=0.;
305 if(lastRRatio>1.) lastRRatio=1.;
306 return lastRRatio;
307 } // End of CalcQF2IN_Ratio
308
309 // Calculatio QasiFree/Inelastic Ratio as a fu
310 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4do
311 {
312 G4double s2=m_s*m_s;
313 G4double s4=s2*s2;
314 G4double ss=std::sqrt(std::sqrt(m_s));
315 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2
316 G4double E=.2644+.016/(1.+G4Exp((29.54-m_s
317 G4double F=ss*.1526*G4Exp(-s2*ss*.0000859)
318 return C*G4Exp(-E*G4Pow::GetInstance()->po
319 } // End of CalcQF2IN_Ratio
320
321 // Calculatio pair(hN_el,hN_tot) (mb): p in Ge
322 std::pair<G4double,G4double> G4QuasiElRatios::
323 {
324 // ---------> Each parameter set can have
325
326 G4double El=0.; // pr
327 G4double To=0.; // pr
328 if(p<=0.)
329 {
330 G4cout<<"-Warning-G4QuasiElRatios::Cal
331 return std::make_pair(El,To);
332 }
333 if (!I) // pp
334 {
335 if(p<pmi)
336 {
337 G4double p2=p*p;
338 El=1./(.00012+p2*.2);
339 To=El;
340 }
341 else if(p>pma)
342 {
343 G4double lp=G4Log(p)-lmi;
344 G4double lp2=lp*lp;
345 El=pbe*lp2+6.72;
346 To=pbt*lp2+38.2;
347 }
348 else
349 {
350 G4double p2=p*p;
351 G4double LE=1./(.00012+p2*.2);
352 G4double lp=G4Log(p)-lmi;
353 G4double lp2=lp*lp;
354 G4double rp2=1./p2;
355 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp
356 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+
357 }
358 }
359 else if(I==1) // np
360 {
361 if(p<pmi)
362 {
363 G4double p2=p*p;
364 El=1./(.00012+p2*(.051+.1*p2));
365 To=El;
366 }
367 else if(p>pma)
368 {
369 G4double lp=G4Log(p)-lmi;
370 G4double lp2=lp*lp;
371 El=pbe*lp2+6.72;
372 To=pbt*lp2+38.2;
373 }
374 else
375 {
376 G4double p2=p*p;
377 G4double LE=1./(.00012+p2*(.051+.1
378 G4double lp=G4Log(p)-lmi;
379 G4double lp2=lp*lp;
380 G4double rp2=1./p2;
381 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49
382 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*r
383 }
384 }
385 else if(I==2) // pi
386 {
387 G4double lp=G4Log(p);
388 if(p<pmi)
389 {
390 G4double lr=lp+1.27;
391 El=1.53/(lr*lr+.0676);
392 To=El*3;
393 }
394 else if(p>pma)
395 {
396 G4double ld=lp-lmi;
397 G4double ld2=ld*ld;
398 G4double sp=std::sqrt(p);
399 El=pbe*ld2+2.4+7./sp;
400 To=pbt*ld2+22.3+12./sp;
401 }
402 else
403 {
404 G4double lr=lp+1.27;
405 G4double LE=1.53/(lr*lr+.0676);
406 G4double ld=lp-lmi;
407 G4double ld2=ld*ld;
408 G4double p2=p*p;
409 G4double p4=p2*p2;
410 G4double sp=std::sqrt(p);
411 G4double lm=lp+.36;
412 G4double md=lm*lm+.04;
413 G4double lh=lp-.017;
414 G4double hd=lh*lh+.0025;
415 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p
416 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+
417 }
418 }
419 else if(I==3) // pi
420 {
421 G4double lp=G4Log(p);
422 if(p<pmi)
423 {
424 G4double lr=lp+1.27;
425 G4double lr2=lr*lr;
426 El=13./(lr2+lr2*lr2+.0676);
427 To=El;
428 }
429 else if(p>pma)
430 {
431 G4double ld=lp-lmi;
432 G4double ld2=ld*ld;
433 G4double sp=std::sqrt(p);
434 El=pbe*ld2+2.4+6./sp;
435 To=pbt*ld2+22.3+5./sp;
436 }
437 else
438 {
439 G4double lr=lp+1.27;
440 G4double lr2=lr*lr;
441 G4double LE=13./(lr2+lr2*lr2+.0676
442 G4double ld=lp-lmi;
443 G4double ld2=ld*ld;
444 G4double p2=p*p;
445 G4double p4=p2*p2;
446 G4double sp=std::sqrt(p);
447 G4double lm=lp-.32;
448 G4double md=lm*lm+.0576;
449 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p
450 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./
451 }
452 }
453 else if(I==4) // Km
454 {
455
456 if(p<pmi)
457 {
458 G4double psp=p*std::sqrt(p);
459 El=5.2/psp;
460 To=14./psp;
461 }
462 else if(p>pma)
463 {
464 G4double ld=G4Log(p)-lmi;
465 G4double ld2=ld*ld;
466 El=pbe*ld2+2.23;
467 To=pbt*ld2+19.5;
468 }
469 else
470 {
471 G4double ld=G4Log(p)-lmi;
472 G4double ld2=ld*ld;
473 G4double sp=std::sqrt(p);
474 G4double psp=p*sp;
475 G4double p2=p*p;
476 G4double p4=p2*p2;
477 G4double lm=p-.39;
478 G4double md=lm*lm+.000156;
479 G4double lh=p-1.;
480 G4double hd=lh*lh+.0156;
481 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/s
482 To=14./psp+(pbt*ld2+19.5)/(1.-.21/
483 }
484 }
485 else if(I==5) // Kp
486 {
487 if(p<pmi)
488 {
489 G4double lr=p-.38;
490 G4double lm=p-1.;
491 G4double md=lm*lm+.372;
492 El=.7/(lr*lr+.0676)+2./md;
493 To=El+.6/md;
494 }
495 else if(p>pma)
496 {
497 G4double ld=G4Log(p)-lmi;
498 G4double ld2=ld*ld;
499 El=pbe*ld2+2.23;
500 To=pbt*ld2+19.5;
501 }
502 else
503 {
504 G4double ld=G4Log(p)-lmi;
505 G4double ld2=ld*ld;
506 G4double lr=p-.38;
507 G4double LE=.7/(lr*lr+.0676);
508 G4double sp=std::sqrt(p);
509 G4double p2=p*p;
510 G4double p4=p2*p2;
511 G4double lm=p-1.;
512 G4double md=lm*lm+.372;
513 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/
514 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.
515 }
516 }
517 else if(I==6) // hy
518 {
519 if(p<pmi)
520 {
521 G4double p2=p*p;
522 El=1./(.002+p2*(.12+p2));
523 To=El;
524 }
525 else if(p>pma)
526 {
527 G4double lp=G4Log(p)-lmi;
528 G4double lp2=lp*lp;
529 G4double sp=std::sqrt(p);
530 El=(pbe*lp2+6.72)/(1.+2./sp);
531 To=(pbt*lp2+38.2+900./sp)/(1.+27./
532 }
533 else
534 {
535 G4double p2=p*p;
536 G4double LE=1./(.002+p2*(.12+p2));
537 G4double lp=G4Log(p)-lmi;
538 G4double lp2=lp*lp;
539 G4double p4=p2*p2;
540 G4double sp=std::sqrt(p);
541 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2.
542 To=LE+(pbt*lp2+38.2+900./sp)/(1.+2
543 }
544 }
545 else if(I==7) // an
546 {
547 if(p>pma)
548 {
549 G4double lp=G4Log(p)-lmi;
550 G4double lp2=lp*lp;
551 El=pbe*lp2+6.72;
552 To=pbt*lp2+38.2;
553 }
554 else
555 {
556 G4double ye=G4Pow::GetInstance()->
557 G4double yt=G4Pow::GetInstance()->
558 G4double lp=G4Log(p)-lmi;
559 G4double lp2=lp*lp;
560 El=80./(ye+1.)+pbe*lp2+6.72;
561 To=(80./yt+.3)/yt+pbt*lp2+38.2;
562 }
563 }
564 else
565 {
566 G4cout<<"*Error*G4QuasiElRatios::CalcE
567 G4Exception("G4QuasiElRatios::CalcElTo
568 }
569 if(El>To) El=To;
570 return std::make_pair(El,To);
571 } // End of CalcElTot
572
573 // For hadron PDG with momentum Mom (GeV/c) on
574 std::pair<G4double,G4double> G4QuasiElRatios::
575 {
576 G4int ind=0;
577 G4bool kfl=true;
578 G4bool kf=false;
579 if(PDG==130||PDG==310)
580 {
581 kf=true;
582 if(G4UniformRand()>.5) kfl=false;
583 }
584 if ( (PDG == 2212 && F) || (PDG == 21
585 else if ( (PDG == 2112 && F) || (PDG == 22
586 else if ( (PDG == -211 && F) || (PDG == 21
587 else if ( (PDG == 211 && F) || (PDG == -21
588 //AR-Jul2020: Extended to charmed and bott
589 // - treat mesons with constituent c quark
590 // - treat mesons with constituent cbar an
591 // - treat all heavy baryons (i.e. hyperon
592 // - treat all heavy anti-baryons (i.e. an
593 else if ( PDG == -321 || PDG == -311 || (k
594 PDG == 411 || PDG == 421 || PD
595 PDG == -521 || PDG == -511 || PDG == -
596 else if ( PDG == 321 || PDG == 311 || (k
597 PDG == -411 || PDG == -421 || PD
598 PDG == 521 || PDG == 511 || PDG ==
599 else if ( PDG > 3000 && PDG < 5333 ) ind
600 else if ( PDG > -5333 && PDG < -2000 ) ind
601 else {
602 G4cout<<"*Error*G4QuasiElRatios::CalcE
603 <<", while it is defined only for p,n,
604 G4Exception("G4QuasiElRatio::CalcElTot
605 }
606 return CalcElTot(p,ind);
607 }
608
609 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV
610 std::pair<G4double,G4double> G4QuasiElRatios::
611 {
612 // LogTable is created only if necessary.
613 G4int nDB=(G4int)vItot.size(); // A numbe
614 if(nDB && lastHtot==PDG && lastFtot==F &&
615 // if(nDB && lastH==PDG && lastF==F && p>
616 lastHtot=PDG;
617 lastFtot=F;
618 G4int ind=-1; //
619 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pip
620 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i
621 G4bool kfl=true;
622 G4bool kf=false;
623 if(PDG==130||PDG==310)
624 {
625 kf=true;
626 if(G4UniformRand()>.5) kfl=false;
627 }
628 if ( (PDG == 2212 && F) || (PDG == 21
629 else if ( (PDG == 2112 && F) || (PDG == 22
630 else if ( (PDG == -211 && F) || (PDG == 21
631 else if ( (PDG == 211 && F) || (PDG == -21
632 //AR-Jul2020: Extended to charmed and bott
633 // - treat mesons with constituent c quark
634 // - treat mesons with constituent cbar an
635 // - treat all heavy baryons (i.e. hyperon
636 // - treat all heavy anti-baryons (i.e. an
637 else if ( PDG == -321 || PDG == -311 || (k
638 PDG == 411 || PDG == 421 || PD
639 PDG == -521 || PDG == -511 || PDG == -
640 else if ( PDG == 321 || PDG == 311 || (k
641 PDG == -411 || PDG == -421 || PD
642 PDG == 521 || PDG == 511 || PDG ==
643 else if ( PDG > 3000 && PDG < 5333 ) ind
644 else if ( PDG > -5333 && PDG < -2000 ) ind
645 else {
646 G4cout<<"*Error*G4QuasiElRatios::Fetch
647 <<", while it is defined only for p,n,
648 G4Exception("G4QuasiELRatio::FetchElTo
649 }
650 if(nDB && lastItot==ind && p>0. && p==last
651 // if(nDB && lastI==ind && p>0. && std::f
652 if(p<=mip || p>=map) return CalcElTot(p,in
653 G4bool found=false;
654 G4int i=-1;
655 if(nDB) for (i=0; i<nDB; i++) if(ind==vIto
656 {
657 found=true;
658 break;
659 }
660 G4double lp=G4Log(p);
661 if(!nDB || !found)
662 {
663 lastXtot = new std::pair<G4double,G4do
664 lastItot = ind;
665 lastKtot = static_cast<int>((lp-lpi)/d
666 if(lastKtot>nlp)
667 {
668 lastKtot=nlp;
669 lastMtot=lpa-lpi;
670 }
671 else lastMtot = lastKtot*dlp;
672 G4double pv=mip;
673 for(G4int j=0; j<=lastKtot; j++)
674 {
675 lastXtot[j]=CalcElTot(pv,ind);
676 if(j!=lastKtot) pv*=edlp;
677 }
678 i++; /
679 vItot.push_back(lastItot);
680 vMtot.push_back(lastMtot);
681 vKtot.push_back(lastKtot);
682 vX->push_back(lastXtot);
683 }
684 else //
685 {
686 lastItot=vItot[i];
687 lastMtot=vMtot[i];
688 lastKtot=vKtot[i];
689 lastXtot=(*vX)[i];
690 G4int nextK=lastKtot+1;
691 G4double lpM=lastMtot+lpi;
692 if(lp>lpM && lastKtot<nlp)
693 {
694 lastKtot = static_cast<int>((lp-lp
695 if(lastKtot>nlp)
696 {
697 lastKtot=nlp;
698 lastMtot=lpa-lpi;
699 }
700 else lastMtot = lastKtot*dlp;
701 G4double pv=G4Exp(lpM); // m
702 for(G4int j=nextK; j<=lastKtot; j+
703 {
704 pv*=edlp;
705 lastXtot[j]=CalcElTot(pv,ind);
706 }
707 } // End of LogTab update
708 if(lastKtot>=nextK)
709 {
710 vMtot[i]=lastMtot;
711 vKtot[i]=lastKtot;
712 }
713 }
714 // Now one can use tabeles to calculate th
715 G4double dlpp=lp-lpi;
716 G4int n=static_cast<int>(dlpp/dlp);
717 G4double d=dlpp-n*dlp;
718 G4double e=lastXtot[n].first;
719 lastRtot.first=e+d*(lastXtot[n+1].first-e)
720 if(lastRtot.first<0.) lastRtot.first = 0.
721 G4double t=lastXtot[n].second;
722 lastRtot.second=t+d*(lastXtot[n+1].second-
723 if(lastRtot.second<0.) lastRtot.second= 0.
724 if(lastRtot.first>lastRtot.second) lastRto
725 return lastRtot;
726 } // End of FetchElTot
727
728 // (Mean Elastic and Mean Total) Cross-Section
729 std::pair<G4double,G4double> G4QuasiElRatios::
730
731 {
732 G4double pGeV=pIU/gigaelectronvolt;
733 if(Z<1 && N<1)
734 {
735 G4cout<<"-Warning-G4QuasiElRatio::GetE
736 return std::make_pair(0.,0.);
737 }
738 std::pair<G4double,G4double> hp=FetchElTot
739 std::pair<G4double,G4double> hn=FetchElTot
740 G4double A=(Z+N)/millibarn;
741 return std::make_pair((Z*hp.first+N*hn.fir
742 } // End of GetElTot
743
744 // (Mean Elastic and Mean Total) Cross-Section
745 std::pair<G4double,G4double> G4QuasiElRatios::
746
747 {
748 G4double pGeV=pIU/gigaelectronvolt;
749 G4double resP=0.;
750 G4double resN=0.;
751 if(Z<1 && N<1)
752 {
753 G4cout<<"-Warning-G4QuasiElRatio::GetC
754 return std::make_pair(resP,resN);
755 }
756 G4double A=Z+N;
757 G4double pf=0.;
758 G4double nf=0.;
759 if (hPDG==-211||hPDG==-321||hPDG==3112||
760 else if(hPDG==211||hPDG==321||hPDG==3222||
761 else if(hPDG==-311||hPDG==311||hPDG==130||
762 {
763 G4double dA=A+A;
764 pf=Z/(dA+N+N);
765 nf=N/(dA+Z+Z);
766 }
767 G4double mult=1.; // Factor of increasing
768 if(pGeV>.5)
769 {
770 mult=1./(1.+G4Log(pGeV+pGeV))/pGeV;
771 if(mult>1.) mult=1.;
772 }
773 if(pf)
774 {
775 std::pair<G4double,G4double> hp=FetchE
776 resP=pf*(hp.second/hp.first-1.)*mult;
777 }
778 if(nf)
779 {
780 std::pair<G4double,G4double> hn=FetchE
781 resN=nf*(hn.second/hn.first-1.)*mult;
782 }
783 return std::make_pair(resP,resN);
784 } // End of GetChExFactor
785
786 // scatter (pPDG,p4M) on a virtual nucleon (NP
787 // if(newN4M.e()==0.) - below threshold, XS=0,
788 std::pair<G4LorentzVector,G4LorentzVector> G4Q
789
790 {
791 static const G4double mNeut= G4Neutron::Ne
792 static const G4double mProt= G4Proton::Pro
793 static const G4double mDeut= G4Deuteron::D
794 static const G4double mTrit= G4Triton::Tri
795 static const G4double mHel3= G4He3::He3()-
796 static const G4double mAlph= G4Alpha::Alph
797
798 G4LorentzVector pr4M=p4M/megaelectronvolt;
799 N4M/=megaelectronvolt;
800 G4LorentzVector tot4M=N4M+p4M;
801 G4double mT=mNeut;
802 G4int Z=0;
803 G4int N=1;
804 if(NPDG==2212||NPDG==90001000)
805 {
806 mT=mProt;
807 Z=1;
808 N=0;
809 }
810 else if(NPDG==90001001)
811 {
812 mT=mDeut;
813 Z=1;
814 N=1;
815 }
816 else if(NPDG==90002001)
817 {
818 mT=mHel3;
819 Z=2;
820 N=1;
821 }
822 else if(NPDG==90001002)
823 {
824 mT=mTrit;
825 Z=1;
826 N=2;
827 }
828 else if(NPDG==90002002)
829 {
830 mT=mAlph;
831 Z=2;
832 N=2;
833 }
834 else if(NPDG!=2112&&NPDG!=90000001)
835 {
836 G4cout<<"Error:G4QuasiElRatios::Scatte
837 G4Exception("G4QuasiElRatios::Scatter:
838 //return std::make_pair(G4LorentzVecto
839 }
840 G4double mT2=mT*mT;
841 G4double mP2=pr4M.m2();
842 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
843 G4double E2=E*E;
844 if(E<0. || E2<mP2)
845 {
846 return std::make_pair(G4LorentzVector(
847 }
848 G4double P=std::sqrt(E2-mP2);
849 // @@ Temporary NN t-dependence for all ha
850 //if(pPDG>3400 || pPDG<-3400) G4cout<<"-Wa
851 G4int PDG=2212;
852 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG
853 if(!Z && N==1) // Change f
854 {
855 Z=1;
856 N=0;
857 if (PDG==2212) PDG=2112;
858 else if(PDG==2112) PDG=2212;
859 }
860 G4double xSec=0.; /
861 if(PDG==2212) xSec=PCSmanager->GetChipsCro
862 else xSec=NCSmanager->GetChipsCro
863 // @@ check a possibility to separate p, n
864 if(xSec <= 0.)
865 {
866 return std::make_pair(G4LorentzVector(
867 }
868 G4double mint=0.; /
869 if(PDG==2212) mint=PCSmanager->GetExchange
870 else mint=NCSmanager->GetExchange
871 G4double maxt=0.;
872 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
873 else maxt=NCSmanager->GetHMaxT();
874 G4double cost=1.-(mint+mint)/maxt; // cos(
875 if(cost>1. || cost<-1. || !(cost>-1. || co
876 {
877 if (cost>1.) cost=1.;
878 else if(cost<-1.) cost=-1.;
879 else
880 {
881 G4double tm=0.;
882 if(PDG==2212) tm=PCSmanager->GetHM
883 else tm=NCSmanager->GetHM
884 G4cerr<<"G4QuasiFreeRatio::Scat:*N
885 return std::make_pair(G4LorentzVec
886 }
887 }
888 G4LorentzVector reco4M=G4LorentzVector(0.,
889 G4LorentzVector dir4M=tot4M-G4LorentzVecto
890 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M
891 {
892 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M
893 //G4Exception("G4QFR::Scat:","009",Fat
894 return std::make_pair(G4LorentzVector(
895 }
896 return std::make_pair(reco4M*megaelectronv
897 } // End of Scatter
898
899 // scatter (pPDG,p4M) on a virtual nucleon (NP
900 // if(newN4M.e()==0.) - below threshold, XS=0,
901 // User should himself change the charge (PDG)
902 //AR-Jul2020: No need to change this method in
903 // charmed and bottom mesons, becau
904 std::pair<G4LorentzVector,G4LorentzVector> G4Q
905
906 {
907 static const G4double mNeut= G4Neutron::Ne
908 static const G4double mProt= G4Proton::Pro
909 G4LorentzVector pr4M=p4M/megaelectronvolt;
910 N4M/=megaelectronvolt;
911 G4LorentzVector tot4M=N4M+p4M;
912 G4int Z=0;
913 G4int N=1;
914 G4int sPDG=0;
915 G4double mS=0.;
916 G4double mT=mProt;
917 if(NPDG==2212)
918 {
919 mT=mNeut;
920 Z=1;
921 N=0;
922 if(pPDG==-211) sPDG=111;
923 else if(pPDG==-321)
924 {
925 sPDG=310;
926 if(G4UniformRand()>.5) sPDG=130;
927 }
928 else if(pPDG==-311||pPDG==311||pPDG==1
929 else if(pPDG==3112) sPDG=3212;
930 else if(pPDG==3212) sPDG=3222;
931 else if(pPDG==3312) sPDG=3322;
932 }
933 else if(NPDG==2112) // Default
934 {
935 if(pPDG==211) sPDG=111;
936 else if(pPDG==321)
937 {
938 sPDG=310;
939 if(G4UniformRand()>.5) sPDG=130;
940 }
941 else if(pPDG==-311||pPDG==311||pPDG==1
942 else if(pPDG==3222) sPDG=3212;
943 else if(pPDG==3212) sPDG=3112;
944 else if(pPDG==3322) sPDG=3312;
945 }
946 else
947 {
948 G4cout<<"Error:G4QuasiElRatios::ChExer
949 G4Exception("G4QuasiElRatios::ChExer:"
950 //return std::make_pair(G4LorentzVecto
951 }
952 if(sPDG) mS=mNeut;
953 else
954 {
955 G4cout<<"Error:G4QuasiElRatios::ChExer
956 G4Exception("G4QuasiElRatios::ChExer:"
957 //return std::make_pair(G4LorentzVecto
958 }
959 G4double mT2=mT*mT;
960 G4double mS2=mS*mS;
961 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
962 G4double E2=E*E;
963 if(E<0. || E2<mS2)
964 {
965 return std::make_pair(G4LorentzVector(
966 }
967 G4double P=std::sqrt(E2-mS2);
968 // @@ Temporary NN t-dependence for all ha
969 G4int PDG=2212;
970 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG
971 if(!Z && N==1) // Change f
972 {
973 Z=1;
974 N=0;
975 if (PDG==2212) PDG=2112;
976 else if(PDG==2112) PDG=2212;
977 }
978 G4double xSec=0.; /
979 if(PDG==2212) xSec=PCSmanager->GetChipsCro
980 else xSec=NCSmanager->GetChipsCro
981 // @@ check a possibility to separate p, n
982 if(xSec <= 0.) // The cross-section iz 0 -
983 {
984 return std::make_pair(G4LorentzVector(
985 }
986 G4double mint=0.; /
987 if(PDG==2212) mint=PCSmanager->GetExchange
988 else mint=NCSmanager->GetExchange
989 G4double maxt=0.;
990 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
991 else maxt=NCSmanager->GetHMaxT();
992 G4double cost=1.-mint/maxt;
993 if(cost>1. || cost<-1. || !(cost>-1. || co
994 {
995 if (cost>1.) cost=1.;
996 else if(cost<-1.) cost=-1.;
997 else
998 {
999 G4cerr<<"G4QuasiFreeRatio::ChExer:
1000 return std::make_pair(G4LorentzVe
1001 }
1002 }
1003 G4LorentzVector reco4M=G4LorentzVector(0.
1004 pr4M=G4LorentzVector(0.,0.,0.,mS);
1005 G4LorentzVector dir4M=tot4M-G4LorentzVect
1006 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4
1007 {
1008 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4
1009 //G4Exception("G4QFR::ChExer:","009",
1010 return std::make_pair(G4LorentzVector
1011 }
1012 return std::make_pair(reco4M*megaelectron
1013 } // End of ChExer
1014
1015 // Calculate ChEx/El ratio (p is in independe
1016 G4double G4QuasiElRatios::ChExElCoef(G4double
1017 {
1018
1019 p/=MeV; //
1020 G4double A=Z+N;
1021 if(A<1.5) return 0.;
1022 G4double Cex=0.;
1023 if (pPDG==2212) Cex=N/(A+Z);
1024 else if(pPDG==2112) Cex=Z/(A+N);
1025 else G4cout<<"*Warning*G4CohChrgExchange:
1026 Cex*=Cex; // Cohe
1027 // @@ This is true only for nucleons: oth
1028 G4double sp=std::sqrt(p);
1029 G4double p2=p*p;
1030 G4double p4=p2*p2;
1031 G4double dl1=G4Log(p)-5.;
1032 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.
1033 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)
1034 G4double R=U/T;
1035 return Cex*R*R;
1036 }
1037
1038 // Decay of Hadron In2Particles f&s, f is in
1039 G4bool G4QuasiElRatios::RelDecayIn2(G4Lorentz
1040 G4Lorent
1041 {
1042 G4double fM2 = f4Mom.m2();
1043 G4double fM = std::sqrt(fM2);
1044 G4double sM2 = s4Mom.m2();
1045 G4double sM = std::sqrt(sM2);
1046 G4double iM2 = theMomentum.m2();
1047 G4double iM = std::sqrt(iM2);
1048 G4double vP = theMomentum.rho(); //
1049 G4double dE = theMomentum.e(); //
1050 if(dE<vP)
1051 {
1052 G4cerr<<"***G4QHad::RelDecIn2: Tachio
1053 G4double accuracy=.000001*vP;
1054 G4double emodif=std::fabs(dE-vP);
1055 //if(emodif<accuracy)
1056 //{
1057 G4cerr<<"G4QHadron::RelDecIn2: *Boost
1058 theMomentum.setE(vP+emodif+.01*accura
1059 //}
1060 }
1061 G4ThreeVector ltb = theMomentum.boostVect
1062 G4ThreeVector ltf = -ltb; //
1063 G4LorentzVector cdir = dir; //
1064 cdir.boost(ltf); //
1065 G4ThreeVector vdir = cdir.vect(); //
1066 G4ThreeVector vx(0.,0.,1.); //
1067 G4ThreeVector vy(0.,1.,0.); //
1068 G4ThreeVector vz(1.,0.,0.); //
1069 if(vdir.mag2() > 0.) //
1070 {
1071 vx = vdir.unit();
1072 G4ThreeVector vv= vx.orthogonal();
1073 vy = vv.unit();
1074 vz = vx.cross(vy);
1075 }
1076 if(maxCost> 1.) maxCost= 1.;
1077 if(minCost<-1.) minCost=-1.;
1078 if(maxCost<-1.) maxCost=-1.;
1079 if(minCost> 1.) minCost= 1.;
1080 if(minCost> maxCost) minCost=maxCost;
1081 if(std::fabs(iM-fM-sM)<.00000001)
1082 {
1083 G4double fR=fM/iM;
1084 G4double sR=sM/iM;
1085 f4Mom=fR*theMomentum;
1086 s4Mom=sR*theMomentum;
1087 return true;
1088 }
1089 else if (iM+.001<fM+sM || iM==0.)
1090 {//@@ Later on make a quark content check
1091 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM
1092 return false;
1093 }
1094 G4double d2 = iM2-fM2-sM2;
1095 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
1096 if(p2<0.)
1097 {
1098 p2=0.;
1099 }
1100 G4double p = std::sqrt(p2);
1101 G4double ct = maxCost;
1102 if(maxCost>minCost)
1103 {
1104 G4double dcost=maxCost-minCost;
1105 ct = minCost+dcost*G4UniformRand();
1106 }
1107 G4double phi= twopi*G4UniformRand(); //
1108 G4double pss=0.;
1109 if(std::fabs(ct)<1.) pss = p * std::sqrt(
1110 else
1111 {
1112 if(ct>1.) ct=1.;
1113 if(ct<-1.) ct=-1.;
1114 }
1115 G4ThreeVector pVect=(pss*std::sin(phi))*v
1116
1117 f4Mom.setVect(pVect);
1118 f4Mom.setE(std::sqrt(fM2+p2));
1119 s4Mom.setVect((-1)*pVect);
1120 s4Mom.setE(std::sqrt(sM2+p2));
1121
1122 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G
1123 <<f4Mom.e()-f4Mom.rho()<<G4endl;
1124 f4Mom.boost(ltb);
1125 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G
1126 <<s4Mom.e()-s4Mom.rho()<<G4endl;
1127 s4Mom.boost(ltb);
1128 return true;
1129 } // End of "RelDecayIn2"
1130
1131
1132
1133
1134
1135
1136