#include <G4QKaonMinusElasticCrossSection.hh>
Inheritance diagram for G4QKaonMinusElasticCrossSection:
Public Member Functions | |
~G4QKaonMinusElasticCrossSection () | |
virtual G4double | GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=-321) |
G4double | CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int pPDG, G4int Z, G4int N, G4double pP) |
G4double | GetSlope (G4int tZ, G4int tN, G4int pPDG) |
G4double | GetExchangeT (G4int tZ, G4int tN, G4int pPDG) |
G4double | GetHMaxT () |
Static Public Member Functions | |
static G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QKaonMinusElasticCrossSection () |
Definition at line 47 of file G4QKaonMinusElasticCrossSection.hh.
G4QKaonMinusElasticCrossSection::G4QKaonMinusElasticCrossSection | ( | ) | [protected] |
G4QKaonMinusElasticCrossSection::~G4QKaonMinusElasticCrossSection | ( | ) |
Definition at line 104 of file G4QKaonMinusElasticCrossSection.cc.
00105 { 00106 std::vector<G4double*>::iterator pos; 00107 for (pos=CST.begin(); pos<CST.end(); pos++) 00108 { delete [] *pos; } 00109 CST.clear(); 00110 for (pos=PAR.begin(); pos<PAR.end(); pos++) 00111 { delete [] *pos; } 00112 PAR.clear(); 00113 for (pos=SST.begin(); pos<SST.end(); pos++) 00114 { delete [] *pos; } 00115 SST.clear(); 00116 for (pos=S1T.begin(); pos<S1T.end(); pos++) 00117 { delete [] *pos; } 00118 S1T.clear(); 00119 for (pos=B1T.begin(); pos<B1T.end(); pos++) 00120 { delete [] *pos; } 00121 B1T.clear(); 00122 for (pos=S2T.begin(); pos<S2T.end(); pos++) 00123 { delete [] *pos; } 00124 S2T.clear(); 00125 for (pos=B2T.begin(); pos<B2T.end(); pos++) 00126 { delete [] *pos; } 00127 B2T.clear(); 00128 for (pos=S3T.begin(); pos<S3T.end(); pos++) 00129 { delete [] *pos; } 00130 S3T.clear(); 00131 for (pos=B3T.begin(); pos<B3T.end(); pos++) 00132 { delete [] *pos; } 00133 B3T.clear(); 00134 for (pos=S4T.begin(); pos<S4T.end(); pos++) 00135 { delete [] *pos; } 00136 S4T.clear(); 00137 for (pos=B4T.begin(); pos<B4T.end(); pos++) 00138 { delete [] *pos; } 00139 B4T.clear(); 00140 }
G4double G4QKaonMinusElasticCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | pPDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | pP | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 286 of file G4QKaonMinusElasticCrossSection.cc.
References G4cout, and G4endl.
Referenced by GetCrossSection().
00288 { 00289 // *** Begin of Associative Memory DB for acceleration of the cross section calculations 00290 static std::vector <G4double> PIN; // Vector of max initialized log(P) in the table 00291 // *** End of Static Definitions (Associative Memory Data Base) *** 00292 G4double pMom=pIU/GeV; // All calculations are in GeV 00293 onlyCS=CS; // Flag to calculate only CS (not Si/Bi) 00294 #ifdef pdebug 00295 G4cout<<"G4QKaonMinusElasticCS::CalcCS: onlyCS="<<onlyCS<<",F="<<F<<",p="<<pIU<<G4endl; 00296 #endif 00297 lastLP=std::log(pMom); // Make a logarithm of the momentum for calculation 00298 if(F) // This isotope was found in AMDB =>RETRIEVE/UPDATE 00299 { 00300 if(F<0) // the AMDB must be loded 00301 { 00302 lastPIN = PIN[I]; // Max log(P) initialised for this table set 00303 lastPAR = PAR[I]; // Pointer to the parameter set 00304 lastCST = CST[I]; // Pointer to the total sross-section table 00305 lastSST = SST[I]; // Pointer to the first squared slope 00306 lastS1T = S1T[I]; // Pointer to the first mantissa 00307 lastB1T = B1T[I]; // Pointer to the first slope 00308 lastS2T = S2T[I]; // Pointer to the second mantissa 00309 lastB2T = B2T[I]; // Pointer to the second slope 00310 lastS3T = S3T[I]; // Pointer to the third mantissa 00311 lastB3T = B3T[I]; // Pointer to the rhird slope 00312 lastS4T = S4T[I]; // Pointer to the 4-th mantissa 00313 lastB4T = B4T[I]; // Pointer to the 4-th slope 00314 #ifdef pdebug 00315 G4cout<<"G4QKaonMinusElCS::CalCS:DB updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl; 00316 #endif 00317 } 00318 #ifdef pdebug 00319 G4cout<<"G4QKaonMinusElasticCS::CalcCS:*read*,LP="<<lastLP<<",PIN="<<lastPIN<<G4endl; 00320 #endif 00321 if(lastLP>lastPIN && lastLP<lPMax) 00322 { 00323 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs 00324 #ifdef pdebug 00325 G4cout<<"G4QKMElCS::CalcCS:updated(I),LP="<<lastLP<<"<IN["<<I<<"]="<<lastPIN<<G4endl; 00326 #endif 00327 PIN[I]=lastPIN; // Remember the new P-Limit of the tables 00328 } 00329 } 00330 else // This isotope wasn't initialized => CREATE 00331 { 00332 lastPAR = new G4double[nPoints]; // Allocate memory for parameters of CS function 00333 lastPAR[nLast]=0; // Initialization for VALGRIND 00334 lastCST = new G4double[nPoints]; // Allocate memory for Tabulated CS function 00335 lastSST = new G4double[nPoints]; // Allocate memory for Tabulated first sqaredSlope 00336 lastS1T = new G4double[nPoints]; // Allocate memory for Tabulated first mantissa 00337 lastB1T = new G4double[nPoints]; // Allocate memory for Tabulated first slope 00338 lastS2T = new G4double[nPoints]; // Allocate memory for Tabulated second mantissa 00339 lastB2T = new G4double[nPoints]; // Allocate memory for Tabulated second slope 00340 lastS3T = new G4double[nPoints]; // Allocate memory for Tabulated third mantissa 00341 lastB3T = new G4double[nPoints]; // Allocate memory for Tabulated third slope 00342 lastS4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th mantissa 00343 lastB4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th slope 00344 #ifdef pdebug 00345 G4cout<<"G4QKaonMinusElasticCS::CalcCS:*ini*,lstLP="<<lastLP<<",min="<<lPMin<<G4endl; 00346 #endif 00347 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables 00348 #ifdef pdebug 00349 G4cout<<"G4QKaMiElCS::CCS:i,Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",LP"<<lastPIN<<G4endl; 00350 #endif 00351 PIN.push_back(lastPIN); // Fill parameters of CS function to AMDB 00352 PAR.push_back(lastPAR); // Fill parameters of CS function to AMDB 00353 CST.push_back(lastCST); // Fill Tabulated CS function to AMDB 00354 SST.push_back(lastSST); // Fill Tabulated first sq.slope to AMDB 00355 S1T.push_back(lastS1T); // Fill Tabulated first mantissa to AMDB 00356 B1T.push_back(lastB1T); // Fill Tabulated first slope to AMDB 00357 S2T.push_back(lastS2T); // Fill Tabulated second mantissa to AMDB 00358 B2T.push_back(lastB2T); // Fill Tabulated second slope to AMDB 00359 S3T.push_back(lastS3T); // Fill Tabulated third mantissa to AMDB 00360 B3T.push_back(lastB3T); // Fill Tabulated third slope to AMDB 00361 S4T.push_back(lastS4T); // Fill Tabulated 4-th mantissa to AMDB 00362 B4T.push_back(lastB4T); // Fill Tabulated 4-th slope to AMDB 00363 } // End of creation/update of the new set of parameters and tables 00364 // =----------= NOW Update (if necessary) and Calculate the Cross Section =-----------= 00365 #ifdef pdebug 00366 G4cout<<"G4QKMElCS::CalcCS:?update?,LP="<<lastLP<<",IN="<<lastPIN<<",ML="<<lPMax<<G4endl; 00367 #endif 00368 if(lastLP>lastPIN && lastLP<lPMax) 00369 { 00370 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN); 00371 #ifdef pdebug 00372 G4cout<<"G4QKaonMinusElCS::CalcCS:*updated(O)*,LP="<<lastLP<<"<IN="<<lastPIN<<G4endl; 00373 #endif 00374 } 00375 #ifdef pdebug 00376 G4cout<<"G4QKaMiElCS::CalcCS:lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl; 00377 #endif 00378 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2) 00379 #ifdef pdebug 00380 G4cout<<"G4QKaMiElCrosSec::CalcCS: oCS="<<onlyCS<<",-t="<<lastTM<<", p="<<lastLP<<G4endl; 00381 #endif 00382 if(lastLP>lPMin && lastLP<=lastPIN) // Linear fit is made using precalculated tables 00383 { 00384 if(lastLP==lastPIN) 00385 { 00386 G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin 00387 G4int blast=static_cast<int>(shift); // this is a bin number of the lower edge (0) 00388 if(blast<0 || blast>=nLast) G4cout<<"G4QKMElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl; 00389 lastSIG = lastCST[blast]; 00390 if(!onlyCS) // Skip the differential cross-section parameters 00391 { 00392 theSS = lastSST[blast]; 00393 theS1 = lastS1T[blast]; 00394 theB1 = lastB1T[blast]; 00395 theS2 = lastS2T[blast]; 00396 theB2 = lastB2T[blast]; 00397 theS3 = lastS3T[blast]; 00398 theB3 = lastB3T[blast]; 00399 theS4 = lastS4T[blast]; 00400 theB4 = lastB4T[blast]; 00401 } 00402 #ifdef pdebug 00403 G4cout<<"G4QKaonMinusElasticCS::CalculateCS:(E) S1="<<theS1<<",B1="<<theB1<<G4endl; 00404 #endif 00405 } 00406 else 00407 { 00408 G4double shift=(lastLP-lPMin)/dlnP; // a shift from the beginning of the table 00409 G4int blast=static_cast<int>(shift); // the lower bin number 00410 if(blast<0) blast=0; 00411 if(blast>=nLast) blast=nLast-1; // low edge of the last bin 00412 shift-=blast; // step inside the unit bin 00413 G4int lastL=blast+1; // the upper bin number 00414 G4double SIGL=lastCST[blast]; // the basic value of the cross-section 00415 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section 00416 #ifdef pdebug 00417 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCrossSection:Sig="<<lastSIG<<",P=" 00418 <<pMom<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",onlyCS="<<onlyCS<<G4endl; 00419 #endif 00420 if(!onlyCS) // Skip the differential cross-section parameters 00421 { 00422 G4double SSTL=lastSST[blast]; // the low bin of the first squared slope 00423 theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope 00424 G4double S1TL=lastS1T[blast]; // the low bin of the first mantissa 00425 theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa 00426 G4double B1TL=lastB1T[blast]; // the low bin of the first slope 00427 #ifdef pdebug 00428 G4cout<<"G4QKaonMinusElasticCS::CalcCrossSect: b="<<blast<<", ls="<<lastL<<",SL=" 00429 <<S1TL<<",SU="<<lastS1T[lastL]<<",BL="<<B1TL<<",BU="<<lastB1T[lastL]<<G4endl; 00430 #endif 00431 theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope 00432 G4double S2TL=lastS2T[blast]; // the low bin of the second mantissa 00433 theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa 00434 G4double B2TL=lastB2T[blast]; // the low bin of the second slope 00435 theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope 00436 G4double S3TL=lastS3T[blast]; // the low bin of the third mantissa 00437 theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa 00438 #ifdef pdebug 00439 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCS: s3l="<<S3TL<<", sh3="<<shift 00440 <<", s3h="<<lastS3T[lastL]<<", b="<<blast<<", l="<<lastL<<G4endl; 00441 #endif 00442 G4double B3TL=lastB3T[blast]; // the low bin of the third slope 00443 theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope 00444 G4double S4TL=lastS4T[blast]; // the low bin of the 4-th mantissa 00445 theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa 00446 #ifdef pdebug 00447 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCS: s4l="<<S4TL<<", sh4="<<shift 00448 <<",s4h="<<lastS4T[lastL]<<",b="<<blast<<",l="<<lastL<<G4endl; 00449 #endif 00450 G4double B4TL=lastB4T[blast]; // the low bin of the 4-th slope 00451 theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope 00452 } 00453 #ifdef pdebug 00454 G4cout<<"G4QKaonMinusElasticCS::CalculateCS:(I) S1="<<theS1<<",B1="<<theB1<<G4endl; 00455 #endif 00456 } 00457 } 00458 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table 00459 if(lastSIG<0.) lastSIG = 0.; // @@ a Warning print can be added 00460 #ifdef pdebug 00461 G4cout<<"G4QKaonMinusElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl; 00462 #endif 00463 return lastSIG; 00464 }
G4double G4QKaonMinusElasticCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = -321 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 151 of file G4QKaonMinusElasticCrossSection.cc.
References CalculateCrossSection(), G4cout, G4endl, and G4VQCrossSection::ThresholdEnergy().
00153 { 00154 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00155 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00156 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00157 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00158 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00159 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00160 G4double pEn=pMom; 00161 onlyCS=fCS; 00162 #ifdef pdebug 00163 G4cout<<"G4QKaMiusElCS::GetCS:>>>f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N=" 00164 <<tgN<<"("<<lastN<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="<<colN.size()<<G4endl; 00165 //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00166 #endif 00167 if(pPDG != -321) 00168 { 00169 G4cout<<"*Warning*G4QKaonMinusElaCS::GetCS:*>Found pPDG="<<pPDG<<" ==> CS=0"<<G4endl; 00170 //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00171 return 0.; // projectile PDG=0 is a mistake (?!) @@ 00172 } 00173 G4bool in=false; // By default the isotope must be found in the AMDB 00174 lastP = 0.; // New momentum history (nothing to compare with) 00175 lastN = tgN; // The last N of the calculated nucleus 00176 lastZ = tgZ; // The last Z of the calculated nucleus 00177 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00178 if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB 00179 { // The nucleus with projPDG is found in AMDB 00180 if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB 00181 { 00182 lastI=i; 00183 lastTH =colTH[i]; // Last THreshold (A-dependent) 00184 #ifdef pdebug 00185 G4cout<<"G4QKMElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl; 00186 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00187 #endif 00188 if(pEn<=lastTH) 00189 { 00190 #ifdef pdebug 00191 G4cout<<"G4QKMElCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl; 00192 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00193 #endif 00194 return 0.; // Energy is below the Threshold value 00195 } 00196 lastP =colP [i]; // Last Momentum (A-dependent) 00197 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00198 // if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance) 00199 if(lastP == pMom) // Do not recalculate 00200 { 00201 #ifdef pdebug 00202 G4cout<<"G4QKaonMinusElastCS::GetCS: P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00203 #endif 00204 CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // Update param's only 00205 return lastCS*millibarn; // Use theLastCS 00206 } 00207 in = true; // This is the case when the isotop is found in DB 00208 // Momentum pMom is in IU ! @@ Units 00209 #ifdef pdebug 00210 G4cout<<"G4QKMElCS::G:UpdateDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",i="<<i<<G4endl; 00211 #endif 00212 lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // read & update 00213 #ifdef pdebug 00214 G4cout<<"G4QKaonMinusElCS::GetCrosSec:****>New(inDB) Calculated CS="<<lastCS<<G4endl; 00215 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00216 #endif 00217 if(lastCS<=0. && pEn>lastTH) // Correct the threshold 00218 { 00219 #ifdef pdebug 00220 G4cout<<"G4QKaonMinusElCS::GetCS: T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00221 #endif 00222 lastTH=pEn; 00223 } 00224 break; // Go out of the LOOP with found lastI 00225 } 00226 #ifdef pdebug 00227 G4cout<<"---G4QKaonMinusElasticCrossSection::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N=" 00228 <<colN[i]<<",Z["<<i<<"]="<<colZ[i]<<G4endl; 00229 //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00230 #endif 00231 } // End of attampt to find the nucleus in DB 00232 if(!in) // This nucleus has not been calculated previously 00233 { 00234 #ifdef pdebug 00235 G4cout<<"G4QKaonMinusElCS::GetCS:CalNewP="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl; 00236 #endif 00238 lastCS=CalculateCrossSection(fCS,0,lastI,pPDG,lastZ,lastN,pMom);//calculate&create 00239 if(lastCS<=0.) 00240 { 00241 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00242 #ifdef pdebug 00243 G4cout<<"G4QKaonMinusElasticCrossSection::GetCS:NewThr="<<lastTH<<",T="<<pEn<<G4endl; 00244 #endif 00245 if(pEn>lastTH) 00246 { 00247 #ifdef pdebug 00248 G4cout<<"G4QKaonMinusElCS::GetCS:1st T="<<pEn<<"(CS=0)>Threshold="<<lastTH<<G4endl; 00249 #endif 00250 lastTH=pEn; 00251 } 00252 } 00253 #ifdef pdebug 00254 G4cout<<"G4QKaonMinusElCS::GetCS:NCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; 00255 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00256 #endif 00257 colN.push_back(tgN); 00258 colZ.push_back(tgZ); 00259 colP.push_back(pMom); 00260 colTH.push_back(lastTH); 00261 colCS.push_back(lastCS); 00262 #ifdef pdebug 00263 G4cout<<"G4QKMElCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00264 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00265 #endif 00266 return lastCS*millibarn; 00267 } // End of creation of the new set of parameters 00268 else 00269 { 00270 #ifdef pdebug 00271 G4cout<<"G4QKaonMinusElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl; 00272 #endif 00273 colP[lastI]=pMom; 00274 colCS[lastI]=lastCS; 00275 } 00276 #ifdef pdebug 00277 G4cout<<"G4QKMElCS::GetCSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00278 //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST 00279 G4cout<<"G4QKaonMinusElasticCrossSection::GetCrosSec:**End**, onlyCS="<<onlyCS<<G4endl; 00280 #endif 00281 return lastCS*millibarn; 00282 }
Reimplemented from G4VQCrossSection.
Definition at line 734 of file G4QKaonMinusElasticCrossSection.cc.
References G4cout, G4endl, and G4UniformRand.
00735 { 00736 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt; 00737 static const G4double third=1./3.; 00738 static const G4double fifth=1./5.; 00739 static const G4double sevth=1./7.; 00740 #ifdef tdebug 00741 G4cout<<"G4QKaMiElCS::GetExchT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl; 00742 #endif 00743 if(PDG==310 || PDG==130) PDG=-321; 00744 if(PDG!=-321)G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetET:PDG="<<PDG<<G4endl; 00745 if(onlyCS) G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT: onlyCS=1"<<G4endl; 00746 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV) 00747 G4double q2=0.; 00748 if(tgZ==1 && tgN==0) // ===> p+p=p+p 00749 { 00750 #ifdef tdebug 00751 G4cout<<"G4QKaonMinusElCS::GetExT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",S2=" 00752 <<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",GeV2="<<GeVSQ<<G4endl; 00753 #endif 00754 G4double E1=lastTM*theB1; 00755 G4double R1=(1.-std::exp(-E1)); 00756 #ifdef tdebug 00757 G4double ts1=-std::log(1.-R1)/theB1; 00758 G4double ds1=std::fabs(ts1-lastTM)/lastTM; 00759 if(ds1>.0001) 00760 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:1p "<<ts1<<"#"<<lastTM 00761 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl; 00762 #endif 00763 G4double E2=lastTM*theB2; 00764 G4double R2=(1.-std::exp(-E2*E2*E2)); 00765 #ifdef tdebug 00766 G4double ts2=std::pow(-std::log(1.-R2),.333333333)/theB2; 00767 G4double ds2=std::fabs(ts2-lastTM)/lastTM; 00768 if(ds2>.0001) 00769 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:2p "<<ts2<<"#"<<lastTM 00770 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl; 00771 #endif 00772 G4double E3=lastTM*theB3; 00773 G4double R3=(1.-std::exp(-E3)); 00774 #ifdef tdebug 00775 G4double ts3=-std::log(1.-R3)/theB3; 00776 G4double ds3=std::fabs(ts3-lastTM)/lastTM; 00777 if(ds3>.0001) 00778 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:3p "<<ts3<<"#"<<lastTM 00779 <<",d="<<ds3<<",R3="<<R1<<",E3="<<E3<<G4endl; 00780 #endif 00781 G4double I1=R1*theS1/theB1; 00782 G4double I2=R2*theS2; 00783 G4double I3=R3*theS3; 00784 G4double I12=I1+I2; 00785 G4double rand=(I12+I3)*G4UniformRand(); 00786 if (rand<I1 ) 00787 { 00788 G4double ran=R1*G4UniformRand(); 00789 if(ran>1.) ran=1.; 00790 q2=-std::log(1.-ran)/theB1; 00791 } 00792 else if(rand<I12) 00793 { 00794 G4double ran=R2*G4UniformRand(); 00795 if(ran>1.) ran=1.; 00796 q2=-std::log(1.-ran); 00797 if(q2<0.) q2=0.; 00798 q2=std::pow(q2,third)/theB2; 00799 } 00800 else 00801 { 00802 G4double ran=R3*G4UniformRand(); 00803 if(ran>1.) ran=1.; 00804 q2=-std::log(1.-ran)/theB3; 00805 } 00806 } 00807 else 00808 { 00809 G4double a=tgZ+tgN; 00810 #ifdef tdebug 00811 G4cout<<"G4QKaonMinusElasticCrossSection::GetExT:a="<<a<<",t="<<lastTM<<",S1="<<theS1 00812 <<",B1="<<theB1<<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3 00813 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl; 00814 #endif 00815 G4double E1=lastTM*(theB1+lastTM*theSS); 00816 G4double R1=(1.-std::exp(-E1)); 00817 G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check) 00818 #ifdef tdebug 00819 G4double ts1=-std::log(1.-R1)/theB1; 00820 if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss; 00821 G4double ds1=(ts1-lastTM)/lastTM; 00822 if(ds1>.0001) 00823 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:1a "<<ts1<<"#"<<lastTM 00824 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl; 00825 #endif 00826 G4double tm2=lastTM*lastTM; 00827 G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st) 00828 if(a>6.5)E2*=tm2; // for heavy nuclei 00829 G4double R2=(1.-std::exp(-E2)); 00830 #ifdef tdebug 00831 G4double ts2=-std::log(1.-R2)/theB2; 00832 if(a<6.5)ts2=std::pow(ts2,third); 00833 else ts2=std::pow(ts2,fifth); 00834 G4double ds2=std::fabs(ts2-lastTM)/lastTM; 00835 if(ds2>.0001) 00836 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:2a "<<ts2<<"#"<<lastTM 00837 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl; 00838 #endif 00839 G4double E3=lastTM*theB3; 00840 if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA 00841 G4double R3=(1.-std::exp(-E3)); 00842 #ifdef tdebug 00843 G4double ts3=-std::log(1.-R3)/theB3; 00844 if(a>6.5)ts3=std::pow(ts3,sevth); 00845 G4double ds3=std::fabs(ts3-lastTM)/lastTM; 00846 if(ds3>.0001) 00847 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:3a "<<ts3<<"#"<<lastTM 00848 <<",d="<<ds3<<",R3="<<R3<<",E3="<<E3<<G4endl; 00849 #endif 00850 G4double E4=lastTM*theB4; 00851 G4double R4=(1.-std::exp(-E4)); 00852 #ifdef tdebug 00853 G4double ts4=-std::log(1.-R4)/theB4; 00854 G4double ds4=std::fabs(ts4-lastTM)/lastTM; 00855 if(ds4>.0001) 00856 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:4a "<<ts4<<"#"<<lastTM 00857 <<",d="<<ds4<<",R4="<<R4<<",E4="<<E4<<G4endl; 00858 #endif 00859 G4double I1=R1*theS1; 00860 G4double I2=R2*theS2; 00861 G4double I3=R3*theS3; 00862 G4double I4=R4*theS4; 00863 G4double I12=I1+I2; 00864 G4double I13=I12+I3; 00865 G4double rand=(I13+I4)*G4UniformRand(); 00866 #ifdef tdebug 00867 G4cout<<"G4QKMElCS::GExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl; 00868 #endif 00869 if(rand<I1) 00870 { 00871 G4double ran=R1*G4UniformRand(); 00872 if(ran>1.) ran=1.; 00873 q2=-std::log(1.-ran)/theB1; 00874 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss; 00875 #ifdef tdebug 00876 G4cout<<"G4QKMElCS::GExT:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl; 00877 #endif 00878 } 00879 else if(rand<I12) 00880 { 00881 G4double ran=R2*G4UniformRand(); 00882 if(ran>1.) ran=1.; 00883 q2=-std::log(1.-ran)/theB2; 00884 if(q2<0.) q2=0.; 00885 if(a<6.5) q2=std::pow(q2,third); 00886 else q2=std::pow(q2,fifth); 00887 #ifdef tdebug 00888 G4cout<<"G4QKMElCS::GetExT: Q2="<<q2<<",r2="<<R2<<",b2="<<theB2<<",t2="<<ts2<<G4endl; 00889 #endif 00890 } 00891 else if(rand<I13) 00892 { 00893 G4double ran=R3*G4UniformRand(); 00894 if(ran>1.) ran=1.; 00895 q2=-std::log(1.-ran)/theB3; 00896 if(q2<0.) q2=0.; 00897 if(a>6.5) q2=std::pow(q2,sevth); 00898 #ifdef tdebug 00899 G4cout<<"G4QKMElCS::GetExT: Q2="<<q2<<",r3="<<R2<<",b3="<<theB2<<",t3="<<ts2<<G4endl; 00900 #endif 00901 } 00902 else 00903 { 00904 G4double ran=R4*G4UniformRand(); 00905 if(ran>1.) ran=1.; 00906 q2=-std::log(1.-ran)/theB4; 00907 if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0) 00908 #ifdef tdebug 00909 G4cout<<"G4QKMElCS::GET:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl; 00910 #endif 00911 } 00912 } 00913 if(q2<0.) q2=0.; 00914 if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<G4endl; 00915 if(q2>lastTM) 00916 { 00917 #ifdef tdebug 00918 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GET:-t="<<q2<<">"<<lastTM<<G4endl; 00919 #endif 00920 q2=lastTM; 00921 } 00922 return q2*GeVSQ; 00923 }
G4double G4QKaonMinusElasticCrossSection::GetHMaxT | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 949 of file G4QKaonMinusElasticCrossSection.cc.
00950 { 00951 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.; 00952 return lastTM*HGeVSQ; 00953 }
G4VQCrossSection * G4QKaonMinusElasticCrossSection::GetPointer | ( | ) | [static] |
Definition at line 143 of file G4QKaonMinusElasticCrossSection.cc.
Referenced by G4CHIPSElastic::G4CHIPSElastic(), G4CHIPSElasticXS::G4CHIPSElasticXS(), G4QHadronElasticDataSet::GetIsoCrossSection(), G4QElastic::GetMeanFreePath(), and G4QElastic::PostStepDoIt().
00144 { 00145 static G4QKaonMinusElasticCrossSection theCrossSection;//StaticBody of the QElCrossSect 00146 return &theCrossSection; 00147 }
Reimplemented from G4VQCrossSection.
Definition at line 926 of file G4QKaonMinusElasticCrossSection.cc.
References G4cout, and G4endl.
00927 { 00928 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt; 00929 #ifdef tdebug 00930 G4cout<<"G4QKaonMinusECS::GetS:"<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl; 00931 #endif 00932 if(onlyCS)G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetSl:onlCS=true"<<G4endl; 00933 if(lastLP<-4.3) return 0.; // S-wave for p<14 MeV/c (kinE<.1MeV) 00934 if(PDG != -321) 00935 { 00936 // G4cout<<"*Error*G4QKaonMinusElasticCrossSection::GetSlope: PDG="<<PDG<<", Z="<<tgZ 00937 // <<", N="<<tgN<<", while it is defined only for PDG=-321"<<G4endl; 00938 // throw G4QException("G4QKaonMinusElasticCrossSection::GetSlope:Only K- is implemented"); 00939 G4ExceptionDescription ed; 00940 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN 00941 << ", while it is defined only for PDG=-321 (K-)" << G4endl; 00942 } 00943 if(theB1<0.) theB1=0.; 00944 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<G4endl; 00945 return theB1/GeVSQ; 00946 }