66 const G4int n_kmpel=36;
68 const G4double kmp_el[n_kmpel]={5.2,.0557,3.5,2.23,.7,.075,.004,.39,.000156,.15,1.,.0156,5.,
69 74.,3.,3.4,.2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,
70 1.2e6,3.5e6,5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
73 const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
84 dlnP=(lPMax-lPMin)/nLast;
122 std::vector<G4double*>::iterator
pos;
123 for (pos=CST.begin(); pos<CST.end(); pos++)
126 for (pos=PAR.begin(); pos<PAR.end(); pos++)
129 for (pos=SST.begin(); pos<SST.end(); pos++)
132 for (pos=S1T.begin(); pos<S1T.end(); pos++)
135 for (pos=B1T.begin(); pos<B1T.end(); pos++)
138 for (pos=S2T.begin(); pos<S2T.end(); pos++)
141 for (pos=B2T.begin(); pos<B2T.end(); pos++)
144 for (pos=S3T.begin(); pos<S3T.end(); pos++)
147 for (pos=B3T.begin(); pos<B3T.end(); pos++)
150 for (pos=S4T.begin(); pos<S4T.end(); pos++)
153 for (pos=B4T.begin(); pos<B4T.end(); pos++)
161 outFile <<
"G4ChipsKaonMinusElasticXS provides the elastic cross\n"
162 <<
"section for K- nucleus scattering as a function of incident\n"
163 <<
"momentum. The cross section is calculated using M. Kossov's\n"
164 <<
"CHIPS parameterization of cross section data.\n";
200 if(lastI)
for(
G4int i=0; i<lastI; i++)
202 if(colN[i]==tgN && colZ[i]==tgZ)
215 CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom);
220 lastCS=CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom);
221 if(lastCS<=0. && pEn>lastTH)
231 lastCS=CalculateCrossSection(fCS,0,lastI,-321,lastZ,lastN,pMom);
242 colP.push_back(pMom);
243 colTH.push_back(lastTH);
244 colCS.push_back(lastCS);
262 lastLP=std::log(pMom);
280 if(lastLP>lastPIN && lastLP<lPMax)
282 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
300 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
301 PIN.push_back(lastPIN);
302 PAR.push_back(lastPAR);
303 CST.push_back(lastCST);
304 SST.push_back(lastSST);
305 S1T.push_back(lastS1T);
306 B1T.push_back(lastB1T);
307 S2T.push_back(lastS2T);
308 B2T.push_back(lastB2T);
309 S3T.push_back(lastS3T);
310 B3T.push_back(lastB3T);
311 S4T.push_back(lastS4T);
312 B4T.push_back(lastB4T);
315 if(lastLP>lastPIN && lastLP<lPMax)
317 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
319 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
320 if(lastLP>lPMin && lastLP<=lastPIN)
324 G4double shift=(lastLP-lPMin)/dlnP+.000001;
325 G4int blast=
static_cast<int>(shift);
326 if(blast<0 || blast>=nLast)
G4cout<<
"G4QKMElCS::CCS:b="<<blast<<
",n="<<nLast<<
G4endl;
327 lastSIG = lastCST[blast];
330 theSS = lastSST[blast];
331 theS1 = lastS1T[blast];
332 theB1 = lastB1T[blast];
333 theS2 = lastS2T[blast];
334 theB2 = lastB2T[blast];
335 theS3 = lastS3T[blast];
336 theB3 = lastB3T[blast];
337 theS4 = lastS4T[blast];
338 theB4 = lastB4T[blast];
344 G4int blast=
static_cast<int>(shift);
346 if(blast>=nLast) blast=nLast-1;
350 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
354 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
356 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
358 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
360 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
362 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
364 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
366 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
368 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
370 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
374 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
375 if(lastSIG<0.) lastSIG = 0.;
397 if(lastPAR[nLast]!=pwd)
399 if ( tgZ == 1 && tgN == 0 )
401 for (
G4int ip=0; ip<n_kmpel; ip++) lastPAR[ip]=kmp_el[ip];
425 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa));
426 lastPAR[1]=.75*asa/(1.+.009*
a);
427 lastPAR[2]=.1*a2*ssa/(1.+.0015*a2/ssa);
428 lastPAR[3]=1./(1.+500./a2);
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*G4ChipsKaonMinusElasticXS::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*G4ChipsKaonMinusElasticXS::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(
"G4ChipsKaonMinusElasticXS::GetPTables()",
"HAD_CHPS_0000",
619 if(PDG==310 || PDG==130) PDG=-321;
620 if(PDG!=-321)
G4cout<<
"*Warning*G4ChipsKaonMinusElasticXS::GetET:PDG="<<PDG<<
G4endl;
621 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonMinusElasticXS::GetExT: onlyCS=1"<<
G4endl;
629 G4double R2=(1.-std::exp(-E2*E2*E2));
641 q2=-std::log(1.-ran)/theB1;
647 q2=-std::log(1.-ran);
649 q2=std::pow(q2,third)/theB2;
655 q2=-std::log(1.-ran)/theB3;
661 G4double E1=lastTM*(theB1+lastTM*theSS);
669 if(a>6.5)E3*=tm2*tm2*tm2;
684 q2=-std::log(1.-ran)/theB1;
685 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
691 q2=-std::log(1.-ran)/theB2;
693 if(a<6.5) q2=std::pow(q2,third);
694 else q2=std::pow(q2,fifth);
700 q2=-std::log(1.-ran)/theB3;
702 if(a>6.5) q2=std::pow(q2,sevth);
708 q2=-std::log(1.-ran)/theB4;
709 if(a<6.5) q2=lastTM-q2;
713 if(!(q2>=-1.||q2<=1.))
G4cout<<
"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<
G4endl;
724 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonMinusElasticXS::GetSl:onlCS=true"<<
G4endl;
725 if(lastLP<-4.3)
return 0.;
732 ed <<
"PDG = " << PDG <<
", Z = " << tgZ <<
", N = " << tgN
733 <<
", while it is defined only for PDG=-321 (K-)" <<
G4endl;
735 if(theB1<0.) theB1=0.;
736 if(!(theB1>=-1.||theB1<=1.))
G4cout<<
"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<
G4endl;
741 G4double G4ChipsKaonMinusElasticXS::GetHMaxT()
743 return lastTM*HGeVSQ;
750 if(PDG!=-321)
G4cout<<
"*Warning*G4ChipsKaonMinusElasticXS::GetTV:PDG="<<PDG<<
G4endl;
753 G4cout<<
"*Warning*G4QKaonMinusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<
G4endl;
769 if ( tgZ == 1 && tgN == 0 )
773 theS1=(lastPAR[13]+lastPAR[14]*dl2*dl2)/(1.+lastPAR[15]/p4/p)+
774 (lastPAR[16]/p2+lastPAR[17]*
p)/(p4+lastPAR[18]*sp);
775 theB1=lastPAR[19]*std::pow(p,lastPAR[20])/(1.+lastPAR[21]/p3);
776 theS2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]*
p);
777 theB2=lastPAR[25]+lastPAR[26]/(p4+lastPAR[27]/
sp);
778 theS3=lastPAR[28]+lastPAR[29]/(p4*p4+lastPAR[30]*p2+lastPAR[31]);
779 theB3=lastPAR[32]+lastPAR[33]/(p4+lastPAR[34]);
784 return lastPAR[0]/psp+(lastPAR[1]*dp*dp+lastPAR[3])/(1.-lastPAR[4]/sp+lastPAR[5]/p4)+
785 lastPAR[6]/(
sqr(p-lastPAR[7])+lastPAR[8])+lastPAR[9]/(
sqr(p-lastPAR[10])+lastPAR[11]);
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]/p3)/(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(
"G4ChipsKaonMinusElasticXS::GetQ2max()",
"HAD_CHPS_0000",
~G4ChipsKaonMinusElasticXS()
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4ChipsKaonMinusElasticXS()
#define G4MUTEX_INITIALIZER
static G4KaonMinus * KaonMinus()
G4double GetTotalMomentum() const
virtual G4bool IsIsoApplicable(const G4DynamicParticle *Pt, G4int Z, G4int A, const G4Element *elm, const G4Material *mat)
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
virtual void CrossSectionDescription(std::ostream &) const
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
#define G4_DECLARE_XS_FACTORY(cross_section)
G4double GetPDGMass() const
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int tgZ, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
static G4ParticleTable * GetParticleTable()
static constexpr double gigaelectronvolt
static constexpr double GeV
static const G4double pos
static constexpr double millibarn
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)