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