65 const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
67 const G4int n_kppel=35;
69 G4double kpp_el[n_kppel]={.7,.38,.0676,.0557,3.5,2.23,.7,.1,2.,1.,.372,5.,74.,3.,3.4,
70 .2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,
71 5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
86 dlnP=(lPMax-lPMin)/nLast;
124 std::vector<G4double*>::iterator pos;
125 for (pos=CST.begin(); pos<CST.end(); pos++)
128 for (pos=PAR.begin(); pos<PAR.end(); pos++)
131 for (pos=SST.begin(); pos<SST.end(); pos++)
134 for (pos=S1T.begin(); pos<S1T.end(); pos++)
137 for (pos=B1T.begin(); pos<B1T.end(); pos++)
140 for (pos=S2T.begin(); pos<S2T.end(); pos++)
143 for (pos=B2T.begin(); pos<B2T.end(); pos++)
146 for (pos=S3T.begin(); pos<S3T.end(); pos++)
149 for (pos=B3T.begin(); pos<B3T.end(); pos++)
152 for (pos=S4T.begin(); pos<S4T.end(); pos++)
155 for (pos=B4T.begin(); pos<B4T.end(); pos++)
184 static G4ThreadLocal std::vector <G4int> *colN_G4MT_TLS_ = 0 ;
if (!colN_G4MT_TLS_) colN_G4MT_TLS_ =
new std::vector <G4int> ; std::vector <G4int> &colN = *colN_G4MT_TLS_;
185 static G4ThreadLocal std::vector <G4int> *colZ_G4MT_TLS_ = 0 ;
if (!colZ_G4MT_TLS_) colZ_G4MT_TLS_ =
new std::vector <G4int> ; std::vector <G4int> &colZ = *colZ_G4MT_TLS_;
186 static G4ThreadLocal std::vector <G4double> *colP_G4MT_TLS_ = 0 ;
if (!colP_G4MT_TLS_) colP_G4MT_TLS_ =
new std::vector <G4double> ; std::vector <G4double> &colP = *colP_G4MT_TLS_;
187 static G4ThreadLocal std::vector <G4double> *colTH_G4MT_TLS_ = 0 ;
if (!colTH_G4MT_TLS_) colTH_G4MT_TLS_ =
new std::vector <G4double> ; std::vector <G4double> &colTH = *colTH_G4MT_TLS_;
188 static G4ThreadLocal std::vector <G4double> *colCS_G4MT_TLS_ = 0 ;
if (!colCS_G4MT_TLS_) colCS_G4MT_TLS_ =
new std::vector <G4double> ; std::vector <G4double> &colCS = *colCS_G4MT_TLS_;
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);
261 static G4ThreadLocal std::vector <G4double> *PIN_G4MT_TLS_ = 0 ;
if (!PIN_G4MT_TLS_) PIN_G4MT_TLS_ =
new std::vector <G4double> ; std::vector <G4double> &PIN = *PIN_G4MT_TLS_;
265 lastLP=std::log(pMom);
283 if(lastLP>lastPIN && lastLP<lPMax)
285 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
303 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
304 PIN.push_back(lastPIN);
305 PAR.push_back(lastPAR);
306 CST.push_back(lastCST);
307 SST.push_back(lastSST);
308 S1T.push_back(lastS1T);
309 B1T.push_back(lastB1T);
310 S2T.push_back(lastS2T);
311 B2T.push_back(lastB2T);
312 S3T.push_back(lastS3T);
313 B3T.push_back(lastB3T);
314 S4T.push_back(lastS4T);
315 B4T.push_back(lastB4T);
318 if(lastLP>lastPIN && lastLP<lPMax)
320 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
322 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
323 if(lastLP>lPMin && lastLP<=lastPIN)
327 G4double shift=(lastLP-lPMin)/dlnP+.000001;
328 G4int blast=
static_cast<int>(shift);
329 if(blast<0 || blast>=nLast)
G4cout<<
"G4QKPElCS::CCS:b="<<blast<<
",n="<<nLast<<
G4endl;
330 lastSIG = lastCST[blast];
333 theSS = lastSST[blast];
334 theS1 = lastS1T[blast];
335 theB1 = lastB1T[blast];
336 theS2 = lastS2T[blast];
337 theB2 = lastB2T[blast];
338 theS3 = lastS3T[blast];
339 theB3 = lastB3T[blast];
340 theS4 = lastS4T[blast];
341 theB4 = lastB4T[blast];
347 G4int blast=
static_cast<int>(shift);
349 if(blast>=nLast) blast=nLast-1;
353 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
357 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
359 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
361 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
363 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
365 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
367 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
369 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
371 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
373 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
377 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
378 if(lastSIG<0.) lastSIG = 0.;
399 if(lastPAR[nLast]!=pwd)
401 if ( tgZ == 1 && tgN == 0 )
403 for (
G4int ip=0; ip<n_kppel; ip++) lastPAR[ip]=kpp_el[ip];
427 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa));
428 lastPAR[1]=.75*asa/(1.+.009*
a);
429 lastPAR[2]=.9*asa*ssa/(1.+.03*
a);
442 lastPAR[10]=1.2e7*a8+380*a17;
443 lastPAR[11]=.7/(1.+4.e-12*a16);
444 lastPAR[12]=2.5/a8/(a4+1.e-16*a32);
446 lastPAR[14]=1.2*a2+2.3;
449 lastPAR[16]=.01/(1.+.0024*a5);
451 lastPAR[18]=9.e-7/(1.+.035*a5);
452 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a);
456 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7);
457 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/
a);
459 lastPAR[24]=1.e5/(a8+2.5e12/a16);
460 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28);
461 lastPAR[26]=.0006*a3;
463 lastPAR[27]=10.+4.e-8*a12*
a;
468 lastPAR[31]=1./(1.+.0001*a8);
469 lastPAR[32]=1.5e-4/(1.+5.e-6*a12);
473 lastPAR[35]=2.e-7*a4;
477 lastPAR[38]=1.e8*std::exp(.32*asa);
478 lastPAR[39]=20.*std::exp(.45*asa);
479 lastPAR[40]=7.e3+2.4e6/a5;
480 lastPAR[41]=2.5e5*std::exp(.085*a3);
483 lastPAR[43]=920.+.03*a8*a3;
484 lastPAR[44]=93.+.0023*a12;
494 lastPAR[ 9]=4.5*std::pow(a,1.15);
495 lastPAR[10]=.06*std::pow(a,.6);
496 lastPAR[11]=.6*a/(1.+2.e15/a16);
497 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32);
498 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5);
499 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12);
501 lastPAR[15]=400./a12+2.e-22*a9;
502 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14);
503 lastPAR[17]=1000./a2+9.5*sa*ssa;
504 lastPAR[18]=4.e-6*a*asa+1.e11/a16;
505 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16);
506 lastPAR[20]=9.+100./
a;
508 lastPAR[21]=.002*a3+3.e7/a6;
509 lastPAR[22]=7.e-15*a4*asa;
510 lastPAR[23]=9000./a4;
512 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4);
513 lastPAR[25]=1.e-5*a2+2.e14/a16;
514 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12);
515 lastPAR[27]=.016*asa/(1.+5.e16/a16);
517 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14));
518 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11);
519 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8);
520 lastPAR[31]=100./asa;
522 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4);
523 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);
524 lastPAR[34]=1.3+3.e5/a4;
525 lastPAR[35]=500./(a2+50.)+3;
526 lastPAR[36]=1.e-9/a+s4a16*s4a16;
528 lastPAR[37]=.4*asa+3.e-9*a6;
529 lastPAR[38]=.0005*a5;
533 lastPAR[41]=.05+.005*
a;
534 lastPAR[42]=7.e-8/sa;
538 lastPAR[46]=3.e32/(a32+1.e32);
542 lastPAR[49]=7.e3*a/(sa+1.);
543 lastPAR[50]=900.*sa/(1.+500./a3);
546 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
553 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN);
567 G4int ini =
static_cast<int>((ILP-lPMin+.000001)/dlnP)+1;
571 G4int fin =
static_cast<int>((LP-lPMin)/dlnP)+1;
572 if(fin>=nPoints) fin=nLast;
581 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN);
595 else G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
596 <<
", Z="<<tgZ<<
", N="<<tgN<<
", i="<<ini<<
" > fin="<<fin<<
", LP="<<LP
597 <<
" > ILP="<<ILP<<
" nothing is done!"<<
G4endl;
599 else G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
600 <<
", Z="<<tgZ<<
", N="<<tgN<<
", i="<<ini<<
">= max="<<nPoints<<
", LP="<<LP
601 <<
" > ILP="<<ILP<<
", lPMax="<<lPMax<<
" nothing is done!"<<
G4endl;
610 ed <<
"PDG = " << PDG <<
", Z = " << tgZ <<
", N = " << tgN
611 <<
", while it is defined only for PDG=321 (K+) " <<
G4endl;
612 G4Exception(
"G4ChipsKaonPlusElasticXS::GetPTables()",
"HAD_CHPS_0000",
621 if(PDG!=321)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetExT:PDG="<<PDG<<
G4endl;
622 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetExT: onlyCS=1"<<
G4endl;
630 G4double R2=(1.-std::exp(-E2*E2*E2));
642 q2=-std::log(1.-ran)/theB1;
648 q2=-std::log(1.-ran);
650 q2=std::pow(q2,third)/theB2;
656 q2=-std::log(1.-ran)/theB3;
662 G4double E1=lastTM*(theB1+lastTM*theSS);
670 if(a>6.5)E3*=tm2*tm2*tm2;
685 q2=-std::log(1.-ran)/theB1;
686 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
692 q2=-std::log(1.-ran)/theB2;
694 if(a<6.5) q2=std::pow(q2,third);
695 else q2=std::pow(q2,fifth);
701 q2=-std::log(1.-ran)/theB3;
703 if(a>6.5) q2=std::pow(q2,sevth);
709 q2=-std::log(1.-ran)/theB4;
710 if(a<6.5) q2=lastTM-q2;
714 if(!(q2>=-1.||q2<=1.))
G4cout<<
"*NAN*G4QKaonPlusElasticCS::GetExchT: -t="<<q2<<
G4endl;
725 if(onlyCS)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetSl:onlCS=true"<<
G4endl;
726 if(lastLP<-4.3)
return 0.;
730 ed <<
"PDG = " << PDG <<
", Z = " << tgZ <<
", N = " << tgN
731 <<
", while it is defined only for PDG=321 (K+)" <<
G4endl;
732 G4Exception(
"G4ChipsKaonPlusElasticXS::GetSlope()",
"HAD_CHPS_0000",
735 if(theB1<0.) theB1=0.;
736 if(!(theB1>=-1.||theB1<=1.))
G4cout<<
"*NAN*G4QKaonPlusElCS::GetSlope:B1="<<theB1<<
G4endl;
741 G4double G4ChipsKaonPlusElasticXS::GetHMaxT()
743 return lastTM*HGeVSQ;
750 if(PDG!=321)
G4cout<<
"*Warning*G4ChipsKaonPlusElasticXS::GetTaV:PDG="<<PDG<<
G4endl;
753 G4cout<<
"*Warning*G4QKaonPlusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<
G4endl;
768 if ( tgZ == 1 && tgN == 0 )
772 theS1=(lastPAR[12]+lastPAR[13]*dl2*dl2)/(1.+lastPAR[14]/p4/p)+
773 (lastPAR[15]/p2+lastPAR[16]*
p)/(p4+lastPAR[17]*sp);
774 theB1=lastPAR[18]*std::pow(p,lastPAR[19])/(1.+lastPAR[20]/p3);
775 theS2=lastPAR[21]+lastPAR[22]/(p4+lastPAR[23]*
p);
776 theB2=lastPAR[24]+lastPAR[25]/(p4+lastPAR[26]/
sp);
777 theS3=lastPAR[27]+lastPAR[28]/(p4*p4+lastPAR[29]*p2+lastPAR[30]);
778 theB3=lastPAR[31]+lastPAR[32]/(p4+lastPAR[33]);
784 return lastPAR[0]/(lastPAR[2]+
sqr(p-lastPAR[1]))+(lastPAR[3]*dp*dp+lastPAR[5])/
785 (1.-lastPAR[6]/sp+lastPAR[7]/p4)
786 +lastPAR[8]/(
sqr(p-lastPAR[9])+lastPAR[10]);
805 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
806 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
807 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
808 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
809 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
810 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
811 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
812 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
813 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
814 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
815 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
819 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
820 lastPAR[13]/(p5+lastPAR[14]/p16);
821 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
822 lastPAR[17]/(1.+lastPAR[18]/p4);
823 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
824 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
825 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
826 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
827 lastPAR[33]/(1.+lastPAR[34]/p6);
828 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
829 theS4=(lastPAR[41]/p4+lastPAR[46]/
p)/(1.+lastPAR[42]/p10)+
830 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
831 theB4=lastPAR[47]/(1.+lastPAR[48]/
p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
836 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p2)/(1.+lastPAR[3]/p2/
sp);
851 G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
852 return dmt*dmt*pP2/mds;
857 ed <<
"PDG = " << PDG <<
",Z = " << tgZ <<
", N = " << tgN
858 <<
", while it is defined only for p projectiles & Z_target>0" <<
G4endl;
859 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)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4ParticleDefinition * GetDefinition() const
#define G4MUTEX_INITIALIZER
G4double GetTotalMomentum() const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
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()
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
~G4ChipsKaonPlusElasticXS()
static G4KaonPlus * KaonPlus()