64 const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
66 const G4int n_kppel=35;
68 G4double kpp_el[n_kppel]={.7,.38,.0676,.0557,3.5,2.23,.7,.1,2.,1.,.372,5.,74.,3.,3.4,
69 .2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,
70 5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
85 dlnP=(lPMax-lPMin)/nLast;
123 std::vector<G4double*>::iterator
pos;
124 for (pos=CST.begin(); pos<CST.end(); pos++)
127 for (pos=PAR.begin(); pos<PAR.end(); pos++)
130 for (pos=SST.begin(); pos<SST.end(); pos++)
133 for (pos=S1T.begin(); pos<S1T.end(); pos++)
136 for (pos=B1T.begin(); pos<B1T.end(); pos++)
139 for (pos=S2T.begin(); pos<S2T.end(); pos++)
142 for (pos=B2T.begin(); pos<B2T.end(); pos++)
145 for (pos=S3T.begin(); pos<S3T.end(); pos++)
148 for (pos=B3T.begin(); pos<B3T.end(); pos++)
151 for (pos=S4T.begin(); pos<S4T.end(); pos++)
154 for (pos=B4T.begin(); pos<B4T.end(); pos++)
162 outFile <<
"G4ChipsKaonPlusElasticXS provides the elastic cross\n"
163 <<
"section for K+ nucleus scattering as a function of incident\n"
164 <<
"momentum. The cross section is calculated using M. Kossov's\n"
165 <<
"CHIPS parameterization of cross section data.\n";
201 if(lastI)
for(
G4int i=0; i<lastI; i++)
203 if(colN[i]==tgN && colZ[i]==tgZ)
216 CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom);
221 lastCS=CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom);
222 if(lastCS<=0. && pEn>lastTH)
232 lastCS=CalculateCrossSection(fCS,0,lastI,321,lastZ,lastN,pMom);
243 colP.push_back(pMom);
244 colTH.push_back(lastTH);
245 colCS.push_back(lastCS);
263 lastLP=std::log(pMom);
281 if(lastLP>lastPIN && lastLP<lPMax)
283 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
301 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
302 PIN.push_back(lastPIN);
303 PAR.push_back(lastPAR);
304 CST.push_back(lastCST);
305 SST.push_back(lastSST);
306 S1T.push_back(lastS1T);
307 B1T.push_back(lastB1T);
308 S2T.push_back(lastS2T);
309 B2T.push_back(lastB2T);
310 S3T.push_back(lastS3T);
311 B3T.push_back(lastB3T);
312 S4T.push_back(lastS4T);
313 B4T.push_back(lastB4T);
316 if(lastLP>lastPIN && lastLP<lPMax)
318 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
320 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
321 if(lastLP>lPMin && lastLP<=lastPIN)
325 G4double shift=(lastLP-lPMin)/dlnP+.000001;
326 G4int blast=
static_cast<int>(shift);
327 if(blast<0 || blast>=nLast)
G4cout<<
"G4QKPElCS::CCS:b="<<blast<<
",n="<<nLast<<
G4endl;
328 lastSIG = lastCST[blast];
331 theSS = lastSST[blast];
332 theS1 = lastS1T[blast];
333 theB1 = lastB1T[blast];
334 theS2 = lastS2T[blast];
335 theB2 = lastB2T[blast];
336 theS3 = lastS3T[blast];
337 theB3 = lastB3T[blast];
338 theS4 = lastS4T[blast];
339 theB4 = lastB4T[blast];
345 G4int blast=
static_cast<int>(shift);
347 if(blast>=nLast) blast=nLast-1;
351 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
355 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
357 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
359 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
361 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
363 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
365 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
367 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
369 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
371 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
375 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
376 if(lastSIG<0.) lastSIG = 0.;
397 if(lastPAR[nLast]!=pwd)
399 if ( tgZ == 1 && tgN == 0 )
401 for (
G4int ip=0; ip<n_kppel; ip++) lastPAR[ip]=kpp_el[ip];
425 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa));
426 lastPAR[1]=.75*asa/(1.+.009*
a);
427 lastPAR[2]=.9*asa*ssa/(1.+.03*
a);
440 lastPAR[10]=1.2e7*a8+380*a17;
441 lastPAR[11]=.7/(1.+4.e-12*a16);
442 lastPAR[12]=2.5/a8/(a4+1.e-16*a32);
444 lastPAR[14]=1.2*a2+2.3;
447 lastPAR[16]=.01/(1.+.0024*a5);
449 lastPAR[18]=9.e-7/(1.+.035*a5);
450 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a);
454 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7);
455 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/
a);
457 lastPAR[24]=1.e5/(a8+2.5e12/a16);
458 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28);
459 lastPAR[26]=.0006*a3;
461 lastPAR[27]=10.+4.e-8*a12*
a;
466 lastPAR[31]=1./(1.+.0001*a8);
467 lastPAR[32]=1.5e-4/(1.+5.e-6*a12);
471 lastPAR[35]=2.e-7*a4;
475 lastPAR[38]=1.e8*std::exp(.32*asa);
476 lastPAR[39]=20.*std::exp(.45*asa);
477 lastPAR[40]=7.e3+2.4e6/a5;
478 lastPAR[41]=2.5e5*std::exp(.085*a3);
481 lastPAR[43]=920.+.03*a8*a3;
482 lastPAR[44]=93.+.0023*a12;
492 lastPAR[ 9]=4.5*std::pow(a,1.15);
493 lastPAR[10]=.06*std::pow(a,.6);
494 lastPAR[11]=.6*a/(1.+2.e15/a16);
495 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32);
496 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5);
497 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12);
499 lastPAR[15]=400./a12+2.e-22*a9;
500 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14);
501 lastPAR[17]=1000./a2+9.5*sa*ssa;
502 lastPAR[18]=4.e-6*a*asa+1.e11/a16;
503 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16);
504 lastPAR[20]=9.+100./
a;
506 lastPAR[21]=.002*a3+3.e7/a6;
507 lastPAR[22]=7.e-15*a4*asa;
508 lastPAR[23]=9000./a4;
510 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4);
511 lastPAR[25]=1.e-5*a2+2.e14/a16;
512 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12);
513 lastPAR[27]=.016*asa/(1.+5.e16/a16);
515 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14));
516 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11);
517 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8);
518 lastPAR[31]=100./asa;
520 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4);
521 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);
522 lastPAR[34]=1.3+3.e5/a4;
523 lastPAR[35]=500./(a2+50.)+3;
524 lastPAR[36]=1.e-9/a+s4a16*s4a16;
526 lastPAR[37]=.4*asa+3.e-9*a6;
527 lastPAR[38]=.0005*a5;
531 lastPAR[41]=.05+.005*
a;
532 lastPAR[42]=7.e-8/sa;
536 lastPAR[46]=3.e32/(a32+1.e32);
540 lastPAR[49]=7.e3*a/(sa+1.);
541 lastPAR[50]=900.*sa/(1.+500./a3);
544 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
551 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN);
565 G4int ini =
static_cast<int>((ILP-lPMin+.000001)/dlnP)+1;
569 G4int fin =
static_cast<int>((LP-lPMin)/dlnP)+1;
570 if(fin>=nPoints) fin=nLast;
574 for(
G4int ip=ini; ip<=fin; ip++)
579 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN);
593 else G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
594 <<
", Z="<<tgZ<<
", N="<<tgN<<
", i="<<ini<<
" > fin="<<fin<<
", LP="<<LP
595 <<
" > ILP="<<ILP<<
" nothing is done!"<<
G4endl;
597 else G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
598 <<
", Z="<<tgZ<<
", N="<<tgN<<
", i="<<ini<<
">= max="<<nPoints<<
", LP="<<LP
599 <<
" > ILP="<<ILP<<
", lPMax="<<lPMax<<
" nothing is done!"<<
G4endl;
608 ed <<
"PDG = " << PDG <<
", Z = " << tgZ <<
", N = " << tgN
609 <<
", while it is defined only for PDG=321 (K+) " <<
G4endl;
610 G4Exception(
"G4ChipsKaonPlusElasticXS::GetPTables()",
"HAD_CHPS_0000",
619 if(PDG!=321)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetExT:PDG="<<PDG<<
G4endl;
620 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetExT: onlyCS=1"<<
G4endl;
628 G4double R2=(1.-std::exp(-E2*E2*E2));
640 q2=-std::log(1.-ran)/theB1;
646 q2=-std::log(1.-ran);
648 q2=std::pow(q2,third)/theB2;
654 q2=-std::log(1.-ran)/theB3;
660 G4double E1=lastTM*(theB1+lastTM*theSS);
668 if(a>6.5)E3*=tm2*tm2*tm2;
683 q2=-std::log(1.-ran)/theB1;
684 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
690 q2=-std::log(1.-ran)/theB2;
692 if(a<6.5) q2=std::pow(q2,third);
693 else q2=std::pow(q2,fifth);
699 q2=-std::log(1.-ran)/theB3;
701 if(a>6.5) q2=std::pow(q2,sevth);
707 q2=-std::log(1.-ran)/theB4;
708 if(a<6.5) q2=lastTM-q2;
712 if(!(q2>=-1.||q2<=1.))
G4cout<<
"*NAN*G4QKaonPlusElasticCS::GetExchT: -t="<<q2<<
G4endl;
723 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetSl:onlCS=true"<<
G4endl;
724 if(lastLP<-4.3)
return 0.;
728 ed <<
"PDG = " << PDG <<
", Z = " << tgZ <<
", N = " << tgN
729 <<
", while it is defined only for PDG=321 (K+)" <<
G4endl;
730 G4Exception(
"G4ChipsKaonPlusElasticXS::GetSlope()",
"HAD_CHPS_0000",
733 if(theB1<0.) theB1=0.;
734 if(!(theB1>=-1.||theB1<=1.))
G4cout<<
"*NAN*G4QKaonPlusElCS::GetSlope:B1="<<theB1<<
G4endl;
739 G4double G4ChipsKaonPlusElasticXS::GetHMaxT()
741 return lastTM*HGeVSQ;
748 if(PDG!=321)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetTaV:PDG="<<PDG<<
G4endl;
751 G4cout<<
"*Warning*G4QKaonPlusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<
G4endl;
766 if ( tgZ == 1 && tgN == 0 )
770 theS1=(lastPAR[12]+lastPAR[13]*dl2*dl2)/(1.+lastPAR[14]/p4/p)+
771 (lastPAR[15]/p2+lastPAR[16]*
p)/(p4+lastPAR[17]*sp);
772 theB1=lastPAR[18]*std::pow(p,lastPAR[19])/(1.+lastPAR[20]/p3);
773 theS2=lastPAR[21]+lastPAR[22]/(p4+lastPAR[23]*
p);
774 theB2=lastPAR[24]+lastPAR[25]/(p4+lastPAR[26]/
sp);
775 theS3=lastPAR[27]+lastPAR[28]/(p4*p4+lastPAR[29]*p2+lastPAR[30]);
776 theB3=lastPAR[31]+lastPAR[32]/(p4+lastPAR[33]);
782 return lastPAR[0]/(lastPAR[2]+
sqr(p-lastPAR[1]))+(lastPAR[3]*dp*dp+lastPAR[5])/
783 (1.-lastPAR[6]/sp+lastPAR[7]/p4)
784 +lastPAR[8]/(
sqr(p-lastPAR[9])+lastPAR[10]);
803 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
804 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
805 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
806 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
807 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
808 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
809 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
810 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
811 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
812 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
813 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
817 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
818 lastPAR[13]/(p5+lastPAR[14]/p16);
819 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
820 lastPAR[17]/(1.+lastPAR[18]/p4);
821 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
822 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
823 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
824 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
825 lastPAR[33]/(1.+lastPAR[34]/p6);
826 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
827 theS4=(lastPAR[41]/p4+lastPAR[46]/
p)/(1.+lastPAR[42]/p10)+
828 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
829 theB4=lastPAR[47]/(1.+lastPAR[48]/
p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
834 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p2)/(1.+lastPAR[3]/p2/
sp);
849 G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
850 return dmt*dmt*pP2/mds;
855 ed <<
"PDG = " << PDG <<
",Z = " << tgZ <<
", N = " << tgN
856 <<
", while it is defined only for p projectiles & Z_target>0" <<
G4endl;
857 G4Exception(
"G4ChipsKaonPlusElasticXS::GetQ2max()",
"HAD_CHPS_0000",
G4ChipsKaonPlusElasticXS()
std::ostringstream G4ExceptionDescription
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int tgZ, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
std::vector< ExP01TrackerHit * > a
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
#define G4MUTEX_INITIALIZER
G4double GetTotalMomentum() const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
virtual G4bool IsIsoApplicable(const G4DynamicParticle *Pt, G4int Z, G4int A, const G4Element *elm, const G4Material *mat)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
#define G4_DECLARE_XS_FACTORY(cross_section)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static constexpr double gigaelectronvolt
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
static constexpr double GeV
~G4ChipsKaonPlusElasticXS()
static G4KaonPlus * KaonPlus()
static const G4double pos
virtual void CrossSectionDescription(std::ostream &) const
static constexpr double millibarn