#include <G4ChipsProtonElasticXS.hh>
Inheritance diagram for G4ChipsProtonElasticXS:
Public Member Functions | |
G4ChipsProtonElasticXS () | |
~G4ChipsProtonElasticXS () | |
virtual G4bool | IsIsoApplicable (const G4DynamicParticle *Pt, G4int Z, G4int A, const G4Element *elm, const G4Material *mat) |
virtual G4double | GetIsoCrossSection (const G4DynamicParticle *, G4int tgZ, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0) |
virtual G4double | GetChipsCrossSection (G4double momentum, G4int Z, G4int N, G4int pdg) |
G4double | GetExchangeT (G4int tZ, G4int tN, G4int pPDG) |
G4double | GetHMaxT () |
Static Public Member Functions | |
static const char * | Default_Name () |
Definition at line 46 of file G4ChipsProtonElasticXS.hh.
G4ChipsProtonElasticXS::G4ChipsProtonElasticXS | ( | ) |
Definition at line 55 of file G4ChipsProtonElasticXS.cc.
00055 :G4VCrossSectionDataSet(Default_Name()), nPoints(128), nLast(nPoints-1) 00056 { 00057 // Initialization of the parameters 00058 lPMin=-8.; // Min tabulated logarithmicMomentum(D) 00059 lPMax= 8.; // Max tabulated logarithmicMomentum(D) 00060 dlnP=(lPMax-lPMin)/nLast;// LogStep in the table(D) 00061 onlyCS=false;// Flag toCalculateOnlyCS(not Si/Bi)(L) 00062 lastSIG=0.; // Last calculated cross section (L) 00063 lastLP=-10.;// Last log(mom_ofTheIncidentHadron)(L) 00064 lastTM=0.; // Last t_maximum (L) 00065 theSS=0.; // The Last sq.slope of 1st difr.Max(L) 00066 theS1=0.; // The Last mantissa of 1st difr.Max(L) 00067 theB1=0.; // The Last slope of 1st difruct.Max(L) 00068 theS2=0.; // The Last mantissa of 2nd difr.Max(L) 00069 theB2=0.; // The Last slope of 2nd difruct.Max(L) 00070 theS3=0.; // The Last mantissa of 3d difr. Max(L) 00071 theB3=0.; // The Last slope of 3d difruct. Max(L) 00072 theS4=0.; // The Last mantissa of 4th difr.Max(L) 00073 theB4=0.; // The Last slope of 4th difruct.Max(L) 00074 lastTZ=0; // Last atomic number of the target 00075 lastTN=0; // Last # of neutrons in the target 00076 lastPIN=0.; // Last initialized max momentum 00077 lastCST=0; // Elastic cross-section table 00078 lastPAR=0; // Parameters for FunctionalCalculation 00079 lastSST=0; // E-dep of sq.slope of the 1st dif.Max 00080 lastS1T=0; // E-dep of mantissa of the 1st dif.Max 00081 lastB1T=0; // E-dep of the slope of the 1st difMax 00082 lastS2T=0; // E-dep of mantissa of the 2nd difrMax 00083 lastB2T=0; // E-dep of the slope of the 2nd difMax 00084 lastS3T=0; // E-dep of mantissa of the 3d difr.Max 00085 lastB3T=0; // E-dep of the slope of the 3d difrMax 00086 lastS4T=0; // E-dep of mantissa of the 4th difrMax 00087 lastB4T=0; // E-dep of the slope of the 4th difMax 00088 lastN=0; // The last N of calculated nucleus 00089 lastZ=0; // The last Z of calculated nucleus 00090 lastP=0.; // Last used in cross section Momentum 00091 lastTH=0.; // Last threshold momentum 00092 lastCS=0.; // Last value of the Cross Section 00093 lastI=0; // The last position in the DAMDB 00094 }
G4ChipsProtonElasticXS::~G4ChipsProtonElasticXS | ( | ) |
Definition at line 97 of file G4ChipsProtonElasticXS.cc.
00098 { 00099 std::vector<G4double*>::iterator pos; 00100 for (pos=CST.begin(); pos<CST.end(); pos++) 00101 { delete [] *pos; } 00102 CST.clear(); 00103 for (pos=PAR.begin(); pos<PAR.end(); pos++) 00104 { delete [] *pos; } 00105 PAR.clear(); 00106 for (pos=SST.begin(); pos<SST.end(); pos++) 00107 { delete [] *pos; } 00108 SST.clear(); 00109 for (pos=S1T.begin(); pos<S1T.end(); pos++) 00110 { delete [] *pos; } 00111 S1T.clear(); 00112 for (pos=B1T.begin(); pos<B1T.end(); pos++) 00113 { delete [] *pos; } 00114 B1T.clear(); 00115 for (pos=S2T.begin(); pos<S2T.end(); pos++) 00116 { delete [] *pos; } 00117 S2T.clear(); 00118 for (pos=B2T.begin(); pos<B2T.end(); pos++) 00119 { delete [] *pos; } 00120 B2T.clear(); 00121 for (pos=S3T.begin(); pos<S3T.end(); pos++) 00122 { delete [] *pos; } 00123 S3T.clear(); 00124 for (pos=B3T.begin(); pos<B3T.end(); pos++) 00125 { delete [] *pos; } 00126 B3T.clear(); 00127 for (pos=S4T.begin(); pos<S4T.end(); pos++) 00128 { delete [] *pos; } 00129 S4T.clear(); 00130 for (pos=B4T.begin(); pos<B4T.end(); pos++) 00131 { delete [] *pos; } 00132 B4T.clear(); 00133 00134 }
static const char* G4ChipsProtonElasticXS::Default_Name | ( | ) | [inline, static] |
Definition at line 54 of file G4ChipsProtonElasticXS.hh.
Referenced by G4HadronElasticPhysics::ConstructProcess(), G4ChipsComponentXS::G4ChipsComponentXS(), G4ChipsElasticModel::G4ChipsElasticModel(), and G4QuasiElRatios::G4QuasiElRatios().
G4double G4ChipsProtonElasticXS::GetChipsCrossSection | ( | G4double | momentum, | |
G4int | Z, | |||
G4int | N, | |||
G4int | pdg | |||
) | [virtual] |
Definition at line 160 of file G4ChipsProtonElasticXS.cc.
Referenced by G4QuasiElRatios::ChExer(), GetIsoCrossSection(), G4ChipsElasticModel::SampleInvariantT(), and G4QuasiElRatios::Scatter().
00161 { 00162 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00163 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00164 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00165 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00166 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00167 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00168 00169 G4double pEn=pMom; 00170 onlyCS=false; 00171 00172 G4bool in=false; // By default the isotope must be found in the AMDB 00173 lastP = 0.; // New momentum history (nothing to compare with) 00174 lastN = tgN; // The last N of the calculated nucleus 00175 lastZ = tgZ; // The last Z of the calculated nucleus 00176 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00177 if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB 00178 { // The nucleus with projPDG is found in AMDB 00179 if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB 00180 { 00181 lastI=i; 00182 lastTH =colTH[i]; // Last THreshold (A-dependent) 00183 if(pEn<=lastTH) 00184 { 00185 return 0.; // Energy is below the Threshold value 00186 } 00187 lastP =colP [i]; // Last Momentum (A-dependent) 00188 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00189 if(lastP == pMom) // Do not recalculate 00190 { 00191 CalculateCrossSection(onlyCS,-1,i,2212,lastZ,lastN,pMom); // Update param's only 00192 return lastCS*millibarn; // Use theLastCS 00193 } 00194 in = true; // This is the case when the isotop is found in DB 00195 // Momentum pMom is in IU ! @@ Units 00196 lastCS=CalculateCrossSection(onlyCS,-1,i,2212,lastZ,lastN,pMom); // read & update 00197 if(lastCS<=0. && pEn>lastTH) // Correct the threshold 00198 { 00199 lastTH=pEn; 00200 } 00201 break; // Go out of the LOOP with found lastI 00202 } 00203 } // End of attampt to find the nucleus in DB 00204 if(!in) // This nucleus has not been calculated previously 00205 { 00207 lastCS=CalculateCrossSection(onlyCS,0,lastI,2212,lastZ,lastN,pMom);//calculate&create 00208 if(lastCS<=0.) 00209 { 00210 lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00211 if(pEn>lastTH) 00212 { 00213 lastTH=pEn; 00214 } 00215 } 00216 colN.push_back(tgN); 00217 colZ.push_back(tgZ); 00218 colP.push_back(pMom); 00219 colTH.push_back(lastTH); 00220 colCS.push_back(lastCS); 00221 return lastCS*millibarn; 00222 } // End of creation of the new set of parameters 00223 else 00224 { 00225 colP[lastI]=pMom; 00226 colCS[lastI]=lastCS; 00227 } 00228 return lastCS*millibarn; 00229 }
Definition at line 614 of file G4ChipsProtonElasticXS.cc.
References G4cout, and G4UniformRand.
Referenced by G4QuasiElRatios::ChExer(), G4ChipsElasticModel::SampleInvariantT(), and G4QuasiElRatios::Scatter().
00615 { 00616 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt; 00617 static const G4double third=1./3.; 00618 static const G4double fifth=1./5.; 00619 static const G4double sevth=1./7.; 00620 if(PDG!=2212) G4cout<<"**Warning*G4ChipsProtonElasticXS::GetExT:PDG="<<PDG<<G4endl; 00621 if(onlyCS) G4cout<<"**Warning*G4ChipsProtonElasticXS::GetExchanT:onlyCS=1"<<G4endl; 00622 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV) 00623 G4double q2=0.; 00624 if(tgZ==1 && tgN==0) // ===> p+p=p+p 00625 { 00626 G4double E1=lastTM*theB1; 00627 G4double R1=(1.-std::exp(-E1)); 00628 G4double E2=lastTM*theB2; 00629 G4double R2=(1.-std::exp(-E2*E2*E2)); 00630 G4double E3=lastTM*theB3; 00631 G4double R3=(1.-std::exp(-E3)); 00632 G4double I1=R1*theS1/theB1; 00633 G4double I2=R2*theS2; 00634 G4double I3=R3*theS3; 00635 G4double I12=I1+I2; 00636 G4double rand=(I12+I3)*G4UniformRand(); 00637 if (rand<I1 ) 00638 { 00639 G4double ran=R1*G4UniformRand(); 00640 if(ran>1.) ran=1.; 00641 q2=-std::log(1.-ran)/theB1; 00642 } 00643 else if(rand<I12) 00644 { 00645 G4double ran=R2*G4UniformRand(); 00646 if(ran>1.) ran=1.; 00647 q2=-std::log(1.-ran); 00648 if(q2<0.) q2=0.; 00649 q2=std::pow(q2,third)/theB2; 00650 } 00651 else 00652 { 00653 G4double ran=R3*G4UniformRand(); 00654 if(ran>1.) ran=1.; 00655 q2=-std::log(1.-ran)/theB3; 00656 } 00657 } 00658 else 00659 { 00660 G4double a=tgZ+tgN; 00661 G4double E1=lastTM*(theB1+lastTM*theSS); 00662 G4double R1=(1.-std::exp(-E1)); 00663 G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check) 00664 G4double tm2=lastTM*lastTM; 00665 G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st) 00666 if(a>6.5)E2*=tm2; // for heavy nuclei 00667 G4double R2=(1.-std::exp(-E2)); 00668 G4double E3=lastTM*theB3; 00669 if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA 00670 G4double R3=(1.-std::exp(-E3)); 00671 G4double E4=lastTM*theB4; 00672 G4double R4=(1.-std::exp(-E4)); 00673 G4double I1=R1*theS1; 00674 G4double I2=R2*theS2; 00675 G4double I3=R3*theS3; 00676 G4double I4=R4*theS4; 00677 G4double I12=I1+I2; 00678 G4double I13=I12+I3; 00679 G4double rand=(I13+I4)*G4UniformRand(); 00680 if(rand<I1) 00681 { 00682 G4double ran=R1*G4UniformRand(); 00683 if(ran>1.) ran=1.; 00684 q2=-std::log(1.-ran)/theB1; 00685 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss; 00686 } 00687 else if(rand<I12) 00688 { 00689 G4double ran=R2*G4UniformRand(); 00690 if(ran>1.) ran=1.; 00691 q2=-std::log(1.-ran)/theB2; 00692 if(q2<0.) q2=0.; 00693 if(a<6.5) q2=std::pow(q2,third); 00694 else q2=std::pow(q2,fifth); 00695 } 00696 else if(rand<I13) 00697 { 00698 G4double ran=R3*G4UniformRand(); 00699 if(ran>1.) ran=1.; 00700 q2=-std::log(1.-ran)/theB3; 00701 if(q2<0.) q2=0.; 00702 if(a>6.5) q2=std::pow(q2,sevth); 00703 } 00704 else 00705 { 00706 G4double ran=R4*G4UniformRand(); 00707 if(ran>1.) ran=1.; 00708 q2=-std::log(1.-ran)/theB4; 00709 if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0) 00710 } 00711 } 00712 if(q2<0.) q2=0.; 00713 if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QElasticCrossSect::GetExchangeT: -t="<<q2<<G4endl; 00714 if(q2>lastTM) 00715 { 00716 q2=lastTM; 00717 } 00718 return q2*GeVSQ; 00719 }
G4double G4ChipsProtonElasticXS::GetHMaxT | ( | ) |
Definition at line 744 of file G4ChipsProtonElasticXS.cc.
Referenced by G4QuasiElRatios::ChExer(), and G4QuasiElRatios::Scatter().
00745 { 00746 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.; 00747 return lastTM*HGeVSQ; 00748 }
G4double G4ChipsProtonElasticXS::GetIsoCrossSection | ( | const G4DynamicParticle * | , | |
G4int | tgZ, | |||
G4int | A, | |||
const G4Isotope * | iso = 0 , |
|||
const G4Element * | elm = 0 , |
|||
const G4Material * | mat = 0 | |||
) | [virtual] |
Reimplemented from G4VCrossSectionDataSet.
Definition at line 146 of file G4ChipsProtonElasticXS.cc.
References GetChipsCrossSection(), and G4DynamicParticle::GetTotalMomentum().
00150 { 00151 G4double pMom=Pt->GetTotalMomentum(); 00152 G4int tgN = A - tgZ; 00153 00154 return GetChipsCrossSection(pMom, tgZ, tgN, 2212); 00155 }
G4bool G4ChipsProtonElasticXS::IsIsoApplicable | ( | const G4DynamicParticle * | Pt, | |
G4int | Z, | |||
G4int | A, | |||
const G4Element * | elm, | |||
const G4Material * | mat | |||
) | [virtual] |
Reimplemented from G4VCrossSectionDataSet.
Definition at line 136 of file G4ChipsProtonElasticXS.cc.
References G4DynamicParticle::GetDefinition(), and G4Proton::Proton().
00139 { 00140 G4ParticleDefinition* particle = Pt->GetDefinition(); 00141 if (particle == G4Proton::Proton() ) return true; 00142 return false; 00143 }