#include <G4QFreeScattering.hh>
Public Member Functions | |
~G4QFreeScattering () | |
std::pair< G4LorentzVector, G4LorentzVector > | Scatter (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
G4QHadronVector * | InElF (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
std::pair< G4double, G4double > | GetElTotMean (G4double pIU, G4int hPDG, G4int Z, G4int N) |
std::pair< G4double, G4double > | GetElTotXS (G4double Mom, G4int PDG, G4bool F) |
std::pair< G4double, G4double > | FetchElTot (G4double pGeV, G4int PDG, G4bool F) |
std::pair< G4double, G4double > | CalcElTot (G4double pGeV, G4int Index) |
Static Public Member Functions | |
static G4QFreeScattering * | GetPointer () |
Protected Member Functions | |
G4QFreeScattering () |
Definition at line 49 of file G4QFreeScattering.hh.
G4QFreeScattering::G4QFreeScattering | ( | ) | [protected] |
G4QFreeScattering::~G4QFreeScattering | ( | ) |
Definition at line 56 of file G4QFreeScattering.cc.
00057 { 00058 std::vector<std::pair<G4double,G4double>*>::iterator pos2; 00059 for(pos2=vX.begin(); pos2<vX.end(); pos2++) delete [] * pos2; 00060 vX.clear(); 00061 }
Definition at line 71 of file G4QFreeScattering.cc.
References FatalException, G4cout, G4endl, G4Exception(), LE, and G4InuclParticleNames::sp.
Referenced by FetchElTot(), G4QuasiFreeRatios::GetElTotXS(), and GetElTotXS().
00072 { 00073 // ---------> Each parameter set can have not more than nPoints=128 parameters 00074 static const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola 00075 static const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient 00076 static const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient 00077 static const G4double pmi=.1; // Below that fast LE calculation is made 00078 static const G4double pma=1000.; // Above that fast HE calculation is made 00079 G4double El=0.; // prototype of the elastic hN cross-section 00080 G4double To=0.; // prototype of the total hN cross-section 00081 if(p<=0.) 00082 { 00083 //G4cout<<"-Warning-G4QFreeScattering::CalcElTot: p="<<p<<" is 0 or negative"<<G4endl; 00084 return std::make_pair(El,To); 00085 } 00086 if (!I) // pp/nn 00087 { 00088 #ifdef debug 00089 G4cout<<"G4QFreeScatter::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl; 00090 #endif 00091 if(p<pmi) 00092 { 00093 G4double p2=p*p; 00094 El=1./(.00012+p2*.2); 00095 To=El; 00096 #ifdef debug 00097 G4cout<<"G4QFreeScatter::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl; 00098 #endif 00099 } 00100 else if(p>pma) 00101 { 00102 G4double lp=std::log(p)-lmi; 00103 G4double lp2=lp*lp; 00104 El=pbe*lp2+6.72; 00105 To=pbt*lp2+38.2; 00106 #ifdef debug 00107 G4cout<<"G4QFreeScat::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl; 00108 #endif 00109 } 00110 else 00111 { 00112 G4double p2=p*p; 00113 G4double LE=1./(.00012+p2*.2); 00114 G4double lp=std::log(p)-lmi; 00115 G4double lp2=lp*lp; 00116 G4double rp2=1./p2; 00117 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p); 00118 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2); 00119 #ifdef debug 00120 G4cout<<"G4QFreeScat::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl; 00121 #endif 00122 } 00123 } 00124 else if(I==1) // np/pn 00125 { 00126 if(p<pmi) 00127 { 00128 G4double p2=p*p; 00129 El=1./(.00012+p2*(.051+.1*p2)); 00130 To=El; 00131 } 00132 else if(p>pma) 00133 { 00134 G4double lp=std::log(p)-lmi; 00135 G4double lp2=lp*lp; 00136 El=pbe*lp2+6.72; 00137 To=pbt*lp2+38.2; 00138 } 00139 else 00140 { 00141 G4double p2=p*p; 00142 G4double LE=1./(.00012+p2*(.051+.1*p2)); 00143 G4double lp=std::log(p)-lmi; 00144 G4double lp2=lp*lp; 00145 G4double rp2=1./p2; 00146 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p); 00147 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2); 00148 } 00149 } 00150 else if(I==2) // pimp/pipn 00151 { 00152 G4double lp=std::log(p); 00153 if(p<pmi) 00154 { 00155 G4double lr=lp+1.27; 00156 El=1.53/(lr*lr+.0676); 00157 To=El*3; 00158 } 00159 else if(p>pma) 00160 { 00161 G4double ld=lp-lmi; 00162 G4double ld2=ld*ld; 00163 G4double sp=std::sqrt(p); 00164 El=pbe*ld2+2.4+7./sp; 00165 To=pbt*ld2+22.3+12./sp; 00166 } 00167 else 00168 { 00169 G4double lr=lp+1.27; // p1 00170 G4double LE=1.53/(lr*lr+.0676); // p2, p3 00171 G4double ld=lp-lmi; // p4 (lmi=3.5) 00172 G4double ld2=ld*ld; 00173 G4double p2=p*p; 00174 G4double p4=p2*p2; 00175 G4double sp=std::sqrt(p); 00176 G4double lm=lp+.36; // p5 00177 G4double md=lm*lm+.04; // p6 00178 G4double lh=lp-.017; // p7 00179 G4double hd=lh*lh+.0025; // p8 00180 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14 00181 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd; 00182 } 00183 } 00184 else if(I==3) // pipp/pimn 00185 { 00186 G4double lp=std::log(p); 00187 if(p<pmi) 00188 { 00189 G4double lr=lp+1.27; 00190 G4double lr2=lr*lr; 00191 El=13./(lr2+lr2*lr2+.0676); 00192 To=El; 00193 } 00194 else if(p>pma) 00195 { 00196 G4double ld=lp-lmi; 00197 G4double ld2=ld*ld; 00198 G4double sp=std::sqrt(p); 00199 El=pbe*ld2+2.4+6./sp; 00200 To=pbt*ld2+22.3+5./sp; 00201 } 00202 else 00203 { 00204 G4double lr=lp+1.27; // p1 00205 G4double lr2=lr*lr; 00206 G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3 00207 G4double ld=lp-lmi; // p4 (lmi=3.5) 00208 G4double ld2=ld*ld; 00209 G4double p2=p*p; 00210 G4double p4=p2*p2; 00211 G4double sp=std::sqrt(p); 00212 G4double lm=lp-.32; // p5 00213 G4double md=lm*lm+.0576; // p6 00214 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11 00215 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md; 00216 } 00217 } 00218 else if(I==4) // Kmp/Kmn/K0p/K0n 00219 { 00220 00221 if(p<pmi) 00222 { 00223 G4double psp=p*std::sqrt(p); 00224 El=5.2/psp; 00225 To=14./psp; 00226 } 00227 else if(p>pma) 00228 { 00229 G4double ld=std::log(p)-lmi; 00230 G4double ld2=ld*ld; 00231 El=pbe*ld2+2.23; 00232 To=pbt*ld2+19.5; 00233 } 00234 else 00235 { 00236 G4double ld=std::log(p)-lmi; 00237 G4double ld2=ld*ld; 00238 G4double sp=std::sqrt(p); 00239 G4double psp=p*sp; 00240 G4double p2=p*p; 00241 G4double p4=p2*p2; 00242 G4double lm=p-.39; 00243 G4double md=lm*lm+.000156; 00244 G4double lh=p-1.; 00245 G4double hd=lh*lh+.0156; 00246 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd; 00247 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd; 00248 } 00249 } 00250 else if(I==5) // Kpp/Kpn/aKp/aKn 00251 { 00252 if(p<pmi) 00253 { 00254 G4double lr=p-.38; 00255 G4double lm=p-1.; 00256 G4double md=lm*lm+.372; 00257 El=.7/(lr*lr+.0676)+2./md; 00258 To=El+.6/md; 00259 } 00260 else if(p>pma) 00261 { 00262 G4double ld=std::log(p)-lmi; 00263 G4double ld2=ld*ld; 00264 El=pbe*ld2+2.23; 00265 To=pbt*ld2+19.5; 00266 } 00267 else 00268 { 00269 G4double ld=std::log(p)-lmi; 00270 G4double ld2=ld*ld; 00271 G4double lr=p-.38; 00272 G4double LE=.7/(lr*lr+.0676); 00273 G4double sp=std::sqrt(p); 00274 G4double p2=p*p; 00275 G4double p4=p2*p2; 00276 G4double lm=p-1.; 00277 G4double md=lm*lm+.372; 00278 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md; 00279 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md; 00280 } 00281 } 00282 else if(I==6) // hyperon-N 00283 { 00284 if(p<pmi) 00285 { 00286 G4double p2=p*p; 00287 El=1./(.002+p2*(.12+p2)); 00288 To=El; 00289 } 00290 else if(p>pma) 00291 { 00292 G4double lp=std::log(p)-lmi; 00293 G4double lp2=lp*lp; 00294 G4double sp=std::sqrt(p); 00295 El=(pbe*lp2+6.72)/(1.+2./sp); 00296 To=(pbt*lp2+38.2+900./sp)/(1.+27./sp); 00297 } 00298 else 00299 { 00300 G4double p2=p*p; 00301 G4double LE=1./(.002+p2*(.12+p2)); 00302 G4double lp=std::log(p)-lmi; 00303 G4double lp2=lp*lp; 00304 G4double p4=p2*p2; 00305 G4double sp=std::sqrt(p); 00306 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4); 00307 To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4); 00308 } 00309 } 00310 else if(I==7) // antibaryon-N 00311 { 00312 if(p>pma) 00313 { 00314 G4double lp=std::log(p)-lmi; 00315 G4double lp2=lp*lp; 00316 El=pbe*lp2+6.72; 00317 To=pbt*lp2+38.2; 00318 } 00319 else 00320 { 00321 G4double ye=std::pow(p,1.25); 00322 G4double yt=std::pow(p,.35); 00323 G4double lp=std::log(p)-lmi; 00324 G4double lp2=lp*lp; 00325 El=80./(ye+1.)+pbe*lp2+6.72; 00326 To=(80./yt+.3)/yt+pbt*lp2+38.2; 00327 } 00328 } 00329 else 00330 { 00331 G4cout<<"*Error*G4QFreeScattering::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl; 00332 G4Exception("G4QFreeScattering::CalcElTot:","23",FatalException,"CHIPScrash"); 00333 } 00334 if(El>To) El=To; 00335 return std::make_pair(El,To); 00336 } // End of CalcElTot
std::pair< G4double, G4double > G4QFreeScattering::FetchElTot | ( | G4double | pGeV, | |
G4int | PDG, | |||
G4bool | F | |||
) |
Definition at line 373 of file G4QFreeScattering.cc.
References CalcElTot(), G4cout, G4endl, G4Exception(), G4UniformRand, JustWarning, and CLHEP::detail::n.
Referenced by G4QuasiFreeRatios::GetChExFactor(), G4QuasiFreeRatios::GetElTot(), and GetElTotMean().
00374 { 00375 static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step) 00376 static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable 00377 static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) 00378 static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) 00379 static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c) 00380 static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV) 00381 static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table 00382 static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table 00383 //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum" 00384 static G4double lastP=0.; // The last momentum for which XS was calculated 00385 static G4int lastH=0; // The last projPDG for which XS was calculated 00386 static G4bool lastF=true; // The last nucleon for which XS was calculated 00387 static std::pair<G4double,G4double> lastR(0.,0.); // The last result 00388 // Local Associative Data Base: 00389 static std::vector<G4int> vI; // Vector of index for which XS was calculated 00390 static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable 00391 static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable 00392 // Last values of the Associative Data Base: 00393 static G4int lastI=0; // The Last index for which XS was calculated 00394 static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable 00395 static G4int lastK=0; // The Last topBin number initialized in LogTable 00396 static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap 00397 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei 00398 G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB 00399 if(PDG==90001000) PDG=2212; 00400 if(PDG==90000001) PDG=2112; 00401 if(PDG==91000000) PDG=3122; 00402 #ifdef pdebug 00403 G4cout<<"G4QFreeScatter::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl; 00404 #endif 00405 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler. 00406 // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR; 00407 lastH=PDG; 00408 lastF=F; 00409 G4int ind=-1; // Prototipe of the index of the PDG/F combination 00410 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p), 00411 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) 00412 G4bool kfl=true; // Flag of K0/aK0 oscillation 00413 G4bool kf=false; 00414 if(PDG==130 || PDG==310) 00415 { 00416 kf=true; 00417 if(G4UniformRand()>.5) kfl=false; 00418 } 00419 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn 00420 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn 00421 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn 00422 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn 00423 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N 00424 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N 00425 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda 00426 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) 00427 else 00428 { 00429 // VI changed Fatal to Warning, because this method cannot do better 00430 // source of problem in classes which make this call 00431 ind = 0; 00432 G4cout<<"*Error*G4QFreeScattering::FetchElTot: PDG="<<PDG 00433 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; 00434 G4Exception("G4QFreeScattering::FetchElTot:","22",JustWarning,"CHIPS problem"); 00435 } 00436 if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler 00437 // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR; 00438 if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?) 00439 G4bool found=false; 00440 G4int i=-1; 00441 if(nDB) for (i=0; i<nDB; ++i) if(ind==vI[i]) // Sirch for this index in AMDB 00442 { 00443 found=true; // The index is found 00444 break; 00445 } 00446 G4double lp=std::log(p); 00447 #ifdef pdebug 00448 G4cout<<"G4QFreeScat::FetchElTot: I="<<ind<<", i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl; 00449 #endif 00450 if(!nDB || !found) // Create new line in the AMDB 00451 { 00452 #ifdef pdebug 00453 G4cout<<"G4QFreeScattering::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl; 00454 #endif 00455 lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot 00456 lastI = ind; // Remember the initialized inex 00457 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB 00458 if(lastK>nlp) 00459 { 00460 lastK=nlp; 00461 lastM=lpa-lpi; 00462 } 00463 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab 00464 G4double pv=mi; 00465 for(G4int j=0; j<=lastK; ++j) // Calculate LogTab values 00466 { 00467 lastX[j]=CalcElTot(pv,ind); 00468 #ifdef pdebug 00469 G4cout<<"G4QFreeScat::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T=" 00470 <<lastX[j].second<<G4endl; 00471 #endif 00472 if(j!=lastK) pv*=edl; 00473 } 00474 i++; // Make a new record to AMDB and position on it 00475 vI.push_back(lastI); 00476 vM.push_back(lastM); 00477 vK.push_back(lastK); 00478 vX.push_back(lastX); 00479 } 00480 else // The A value was found in AMDB 00481 { 00482 lastI=vI[i]; 00483 lastM=vM[i]; 00484 lastK=vK[i]; 00485 lastX=vX[i]; 00486 G4int nextK=lastK+1; 00487 G4double lpM=lastM+lpi; 00488 #ifdef pdebug 00489 G4cout<<"G4QFreeScatt::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl; 00490 #endif 00491 if(lp>lpM && lastK<nlp) // LogTab must be updated 00492 { 00493 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab 00494 #ifdef pdebug 00495 G4cout<<"G4QFreeScat::FetET: K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl; 00496 #endif 00497 if(lastK>nlp) 00498 { 00499 lastK=nlp; 00500 lastM=lpa-lpi; 00501 } 00502 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab 00503 G4double pv=std::exp(lpM); // momentum of the last calculated beam 00504 for(G4int j=nextK; j<=lastK; ++j) // Calculate LogTab values 00505 { 00506 pv*=edl; 00507 lastX[j]=CalcElTot(pv,ind); 00508 #ifdef pdebug 00509 G4cout<<"G4QFreeScat::FetchElTot: U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T=" 00510 <<lastX[j].second<<G4endl; 00511 #endif 00512 } 00513 } // End of LogTab update 00514 if(lastK >= nextK) // The AMDB was apdated 00515 { 00516 vM[i]=lastM; 00517 vK[i]=lastK; 00518 } 00519 } 00520 // Now one can use tabeles to calculate the value 00521 G4double dlp=lp-lpi; // Shifted log(p) value 00522 G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin 00523 G4double d=dlp-n*dl; // Log shift 00524 G4double e=lastX[n].first; // E-Base 00525 lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result 00526 if(lastR.first<0.) lastR.first = 0.; 00527 G4double t=lastX[n].second; // T-Base 00528 lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result 00529 if(lastR.second<0.) lastR.second= 0.; 00530 #ifdef pdebug 00531 G4cout<<"=O=>G4QFreeScat::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl; 00532 #endif 00533 if(lastR.first>lastR.second) lastR.first = lastR.second; 00534 return lastR; 00535 } // End of FetchElTot
std::pair< G4double, G4double > G4QFreeScattering::GetElTotMean | ( | G4double | pIU, | |
G4int | hPDG, | |||
G4int | Z, | |||
G4int | N | |||
) |
Definition at line 538 of file G4QFreeScattering.cc.
References FetchElTot(), G4cout, and G4endl.
00540 { 00541 G4double pGeV=pIU/gigaelectronvolt; 00542 if(hPDG==90001000) hPDG=2212; 00543 if(hPDG==90000001) hPDG=2112; 00544 if(hPDG==91000000) hPDG=3122; 00545 #ifdef pdebug 00546 G4cout<<"->G4QFreeSc::GetElTotMean: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl; 00547 #endif 00548 if(Z<1 && N<1) 00549 { 00550 G4cout<<"-Warning-G4QFreeScat::GetElTotMean: Z="<<Z<<",N="<<N<<", return zero"<<G4endl; 00551 return std::make_pair(0.,0.); 00552 } 00553 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); 00554 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); 00555 #ifdef pdebug 00556 G4cout<<"-OUT->G4QFreeScat::GetElTotMean: hp("<<hp.first<<","<<hp.second<<"), hn(" 00557 <<hn.first<<","<<hn.second<<")"<<G4endl; 00558 #endif 00559 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU) 00560 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A); 00561 } // End of GetElTotMean
Definition at line 340 of file G4QFreeScattering.cc.
References CalcElTot(), G4cout, G4endl, G4Exception(), G4UniformRand, and JustWarning.
00341 { 00342 G4int ind=0; // Prototype of the reaction index 00343 G4bool kfl=true; // Flag of K0/aK0 oscillation 00344 G4bool kf=false; 00345 if(PDG==90001000) PDG=2212; 00346 if(PDG==90000001) PDG=2112; 00347 if(PDG==91000000) PDG=3122; 00348 if(PDG==130 || PDG==310) 00349 { 00350 kf=true; 00351 if(G4UniformRand()>.5) kfl=false; 00352 } 00353 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn 00354 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn 00355 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn 00356 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn 00357 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N 00358 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N 00359 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda 00360 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) 00361 else { 00362 // VI changed Fatal to Warning, because this method cannot do better 00363 // source of problem in classes which make this call 00364 ind = 0; 00365 G4cout<<"*Error*G4QFreeScattering::CalcElTotXS: PDG="<<PDG 00366 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; 00367 G4Exception("G4QFreeScattering::CalcElTotXS:","22",JustWarning,"CHIPS_crash"); 00368 } 00369 return CalcElTot(p,ind); // This is a slow direct method, better use FetchElTot 00370 }
G4QFreeScattering * G4QFreeScattering::GetPointer | ( | ) | [static] |
Definition at line 64 of file G4QFreeScattering.cc.
Referenced by G4QEnvironment::G4QEnvironment(), G4QuasiFreeRatios::G4QuasiFreeRatios(), and G4QEnvironment::operator=().
00065 { 00066 static G4QFreeScattering theQFS; // *** Static body of the Quasi-Free Scattering *** 00067 return &theQFS; 00068 }
G4QHadronVector * G4QFreeScattering::InElF | ( | G4int | NPDG, | |
G4LorentzVector | N4M, | |||
G4int | pPDG, | |||
G4LorentzVector | p4M | |||
) |
Definition at line 667 of file G4QFreeScattering.cc.
References FatalException, G4cout, G4endl, G4Exception(), and G4UniformRand.
00669 { 00670 static const G4double mPi0 = G4QPDGCode(111).GetMass(); 00671 static const G4double mPi = G4QPDGCode(211).GetMass(); // Pi+ (Pi-: -211) 00672 static const G4double mK = G4QPDGCode(321).GetMass(); // K+ (K- : -321) 00673 static const G4double mK0 = G4QPDGCode(311).GetMass(); // aK0 (K0 : -311) 00674 static const G4double mNeut= G4QPDGCode(2112).GetMass(); 00675 static const G4double mProt= G4QPDGCode(2212).GetMass(); 00676 static const G4double mLamb= G4QPDGCode(3122).GetMass(); 00677 //static const G4double mSigZ= G4QPDGCode(3212).GetMass(); 00678 static const G4double mSigM= G4QPDGCode(3112).GetMass(); 00679 static const G4double mSigP= G4QPDGCode(3222).GetMass(); 00680 static const G4double mXiM = G4QPDGCode(3312).GetMass(); 00681 static const G4double mXiZ = G4QPDGCode(3322).GetMass(); 00682 //static const G4double mOmM = G4QPDGCode(3334).GetMass(); 00683 static const G4double third =1./3.; 00684 static const G4double twothd =2./3.; 00685 00686 p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units) 00687 N4M/=megaelectronvolt; 00688 G4LorentzVector c4M=N4M+p4M; 00689 if(pPDG==90001000) pPDG=2212; 00690 if(pPDG==90000001) pPDG=2112; 00691 if(pPDG==91000000) pPDG=3122; 00692 #ifdef ppdebug 00693 G4cout<<"->G4QFR::InElF: p4M="<<p4M<<",N4M="<<N4M<<", c4M="<<c4M<<",NPDG="<<NPDG<<G4endl; 00694 #endif 00695 G4double mT=mNeut; 00696 G4int Z=0; 00697 //G4int N=1; 00698 if(NPDG==2212 || NPDG==90001000) 00699 { 00700 NPDG=2212; 00701 mT=mProt; 00702 Z=1; 00703 //N=0; 00704 } 00705 else if(NPDG!=2112 && NPDG!=90000001) 00706 { 00707 G4cout<<"Error:G4QFreeScattering::InElF:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; 00708 G4Exception("G4QFreeScattering::InElF:","21",FatalException,"CHIPS complain"); 00709 } 00710 else NPDG=2112; 00711 G4double mC2=c4M.m2(); 00712 G4double mC=0.; 00713 if( mC2 < -0.0000001 ) 00714 { 00715 #ifdef ppdebug 00716 G4cout<<"-Warning-G4QFS::InElF: Negative compoundMass ="<<mC2<<", c4M="<<c4M<<G4endl; 00717 #endif 00718 return 0; // Do Nothing Action 00719 } 00720 else if ( mC2 > 0.0000001) mC=std::sqrt(mC2); 00721 #ifdef ppdebug 00722 G4cout<<"G4QFS::InElF: mC ="<<mC<<G4endl; 00723 #endif 00724 G4double mP=0.; 00725 G4double mP2=p4M.m2(); // a projectile squared mass 00726 if( mP2 < -0.0000001 ) 00727 { 00728 #ifdef ppdebug 00729 G4cout<<"-Warning-G4QFS::InElF: Negative projectileMass ="<<mC2<<", c4M="<<c4M<<G4endl; 00730 #endif 00731 return 0; // Do Nothing Action 00732 } 00733 else if ( mP2 > 0.0000001) mP=std::sqrt(mP2); 00734 #ifdef pdebug 00735 G4cout<<"G4QFS::InElF:mT("<<mT<<")+mP("<<mP<<")+mPi0="<<mP+mT+mPi0<<"<? mC="<<mC<<G4endl; 00736 #endif 00737 if(pPDG > 3334 || pPDG < -321) // Annihilation/Inelastic of anti-barions not implemented 00738 { 00739 G4cout<<"-Warning-G4QFreeScat::InElF: pPDG="<<pPDG<<G4endl; 00740 return 0; // Do Nothing Action 00741 } 00742 if(pPDG==130 || pPDG==310) 00743 { 00744 if(G4UniformRand()>.5) pPDG = 311; // K0 00745 else pPDG =-311; // aK0 00746 } 00747 G4int mPDG=111; // Additional newtral pion by default 00748 G4double mM=mPi0; // Default Pi0 mass 00749 G4double r=G4UniformRand(); // A random number to split pi0/pi 00750 if (pPDG == 2212) // p-N 00751 { 00752 if(Z) // p-p 00753 { 00754 if(r<twothd) // -> n + Pi+ + p 00755 { 00756 pPDG=2112; // n 00757 mP=mNeut; 00758 mPDG= 211; // Pi+ 00759 mM=mPi; 00760 } 00761 } 00762 else // p-n 00763 { 00764 if(r<third) // -> n + Pi+ + n 00765 { 00766 pPDG=2112; // n 00767 mP=mNeut; 00768 mPDG= 211; // Pi+ 00769 mM=mPi; 00770 } 00771 else if(r<twothd) // -> p + Pi- + p 00772 { 00773 mPDG=-211; // Pi- 00774 mM=mPi; 00775 NPDG=2212; // p 00776 mT=mProt; 00777 } 00778 } 00779 } 00780 else if (pPDG == 2112) // n-N 00781 { 00782 if(Z) // n-p 00783 { 00784 if(r<third) // -> n + Pi+ + n 00785 { 00786 mPDG= 211; // Pi+ 00787 mM=mPi; 00788 NPDG=2112; // n 00789 mT=mNeut; 00790 } 00791 else if(r<twothd) // -> p + Pi- + p 00792 { 00793 pPDG=2212; // p 00794 mP=mProt; 00795 mPDG=-211; // Pi- 00796 mM=mPi; 00797 } 00798 } 00799 else // n-n 00800 { 00801 if(r<twothd) // -> p + Pi- + n 00802 { 00803 pPDG=2212; // p 00804 mP=mProt; 00805 mPDG= -211; // Pi- 00806 mM=mPi; 00807 } 00808 } 00809 } 00810 else if (pPDG == 211) // pip-N 00811 { 00812 if(Z) // pip-p 00813 { 00814 if(r<0.5) // -> Pi+ + Pi+ + n 00815 { 00816 mPDG= 211; // Pi+ 00817 mM=mPi; 00818 NPDG=2112; // n 00819 mT=mNeut; 00820 } 00821 } 00822 else // pip-n 00823 { 00824 if(r<third) // -> Pi0 + Pi0 + p 00825 { 00826 pPDG= 111; // Pi0 00827 mP=mPi0; 00828 NPDG=2212; // p 00829 mT=mProt; 00830 } 00831 else if(r<twothd) // -> Pi+ + Pi- + p 00832 { 00833 mPDG=-211; // Pi- 00834 mM=mPi; 00835 NPDG=2212; // p 00836 mT=mProt; 00837 } 00838 } 00839 } 00840 else if (pPDG == -211) // pim-N 00841 { 00842 if(Z) // pim-p 00843 { 00844 if(r<third) // -> Pi0 + Pi0 + n 00845 { 00846 pPDG= 111; // Pi0 00847 mP=mPi0; 00848 NPDG=2112; // n 00849 mT=mNeut; 00850 } 00851 else if(r<twothd) // -> Pi- + Pi+ + n 00852 { 00853 mPDG= 211; // Pi+ 00854 mM=mPi; 00855 NPDG=2112; // n 00856 mT=mNeut; 00857 } 00858 } 00859 else // pim-n 00860 { 00861 if(r<0.5) // -> Pi- + Pi- + p 00862 { 00863 mPDG=-211; // Pi- 00864 mM=mPi; 00865 NPDG=2212; // p 00866 mT=mProt; 00867 } 00868 } 00869 } 00870 else if (pPDG == -321) // Km-N 00871 { 00872 if(Z) // Km-p 00873 { 00874 if(r<0.25) // K0 + Pi0 + n 00875 { 00876 pPDG=-311; // K0 00877 mP=mK0; 00878 NPDG=2112; // n 00879 mT=mNeut; 00880 } 00881 else if(r<0.5) // -> K- + Pi+ + n 00882 { 00883 mPDG= 211; // Pi+ 00884 mM=mPi; 00885 NPDG=2112; // n 00886 mT=mNeut; 00887 } 00888 else if(r<0.75) // -> K0 + Pi- + p 00889 { 00890 pPDG=-311; // K0 00891 mP=mK0; 00892 mPDG=-211; // Pi- 00893 mM=mPi; 00894 } 00895 } 00896 else // Km-n 00897 { 00898 if(r<third) // -> K- + Pi- + p 00899 { 00900 mPDG=-211; // Pi- 00901 mM=mPi; 00902 NPDG=2212; // p 00903 mT=mProt; 00904 } 00905 else if(r<twothd) // -> K0 + Pi- + n 00906 { 00907 pPDG=-311; // K0 00908 mP=mK0; 00909 mPDG=-211; // Pi- 00910 mM=mPi; 00911 } 00912 } 00913 } 00914 else if (pPDG == -311) // K0-N 00915 { 00916 if(Z) // K0-p 00917 { 00918 if(r<third) // K- + Pi+ + p 00919 { 00920 pPDG=-321; // K- 00921 mP=mK; 00922 NPDG=2212; // p 00923 mT=mProt; 00924 } 00925 else if(r<twothd) // -> K0 + Pi+ + n 00926 { 00927 mPDG= 211; // Pi+ 00928 mM=mPi; 00929 NPDG=2112; // n 00930 mT=mNeut; 00931 } 00932 } 00933 else // K0-n 00934 { 00935 if(r<0.25) // -> K- + Pi+ + n 00936 { 00937 pPDG=-321; // K- 00938 mP=mK; 00939 mPDG= 211; // Pi+ 00940 mM=mPi; 00941 } 00942 else if(r<0.5) // -> K- + Pi0 + p 00943 { 00944 pPDG=-321; // K- 00945 mP=mK; 00946 NPDG=2212; // p 00947 mT=mProt; 00948 } 00949 else if(r<0.75) // -> K0 + Pi- + p 00950 { 00951 mPDG=-211; // Pi- 00952 mM=mPi; 00953 NPDG=2212; // p 00954 mT=mProt; 00955 } 00956 } 00957 } 00958 else if (pPDG == 321) // Kp-N 00959 { 00960 if(Z) // Kp-p 00961 { 00962 if(r<third) // -> K+ + Pi+ + n 00963 { 00964 mPDG= 211; // Pi+ 00965 mM=mPi; 00966 NPDG=2112; // n 00967 mT=mNeut; 00968 } 00969 else if(r<twothd) // -> aK0 + Pi+ + p 00970 { 00971 pPDG= 311; // aK0 00972 mP=mK0; 00973 mPDG= 211; // Pi+ 00974 mM=mPi; 00975 } 00976 } 00977 else // Kp-n 00978 { 00979 if(r<0.25) // -> aK0 + Pi+ + n 00980 { 00981 pPDG= 311; // aK0 00982 mP=mK0; 00983 mPDG= 211; // Pi+ 00984 mM=mPi; 00985 } 00986 else if(r<0.5) // -> Kp + Pi- + p 00987 { 00988 mPDG=-211; // Pi- 00989 mM=mPi; 00990 NPDG=2212; // p 00991 mT=mProt; 00992 } 00993 else if(r<0.75) // -> aK0 + Pi0 + p 00994 { 00995 pPDG= 311; // aK0 00996 mP=mK0; 00997 NPDG=2212; // p 00998 mT=mProt; 00999 } 01000 } 01001 } 01002 else if (pPDG == 311) // aK0-N 01003 { 01004 if(Z) // aK0-p 01005 { 01006 if(r<0.25) // -> K+ + Pi- + p 01007 { 01008 pPDG= 321; // K+ 01009 mP=mK; 01010 mPDG=-211; // Pi- 01011 mM=mPi; 01012 } 01013 else if(r<0.5) // -> aK0 + Pi+ + n 01014 { 01015 mPDG= 211; // Pi+ 01016 mM=mPi; 01017 NPDG=2112; // n 01018 mT=mNeut; 01019 } 01020 else if(r<0.75) // -> K+ + Pi0 + n 01021 { 01022 pPDG= 321; // K+ 01023 mP=mK; 01024 NPDG=2112; // n 01025 mT=mNeut; 01026 } 01027 } 01028 else // aK0-n 01029 { 01030 if(r<third) // -> aK0 + Pi- + p 01031 { 01032 mPDG=-211; // Pi- 01033 mM=mPi; 01034 NPDG=2212; // p 01035 mT=mProt; 01036 } 01037 else if(r<twothd) // -> K+ + Pi- + n 01038 { 01039 pPDG= 321; // K+ 01040 mP=mK; 01041 mPDG=-211; // Pi- 01042 mM=mPi; 01043 } 01044 } 01045 } 01046 else if (pPDG == 3122 || pPDG== 3212) // Lambda/Sigma0-N 01047 { 01048 if(pPDG == 3212) 01049 { 01050 pPDG=3122; 01051 mP=mLamb; 01052 } 01053 if(Z) // L/S0-p 01054 { 01055 if(r<0.2) // -> SigP + Pi0 + n 01056 { 01057 pPDG=3222; // SigP 01058 mP=mSigP; 01059 NPDG=2112; // n 01060 mT=mNeut; 01061 } 01062 else if(r<0.4) // -> SigP + Pi- + p 01063 { 01064 pPDG=3222; // SigP 01065 mP=mSigP; 01066 mPDG=-211; // Pi- 01067 mM=mPi; 01068 } 01069 else if(r<0.6) // -> SigM + Pi+ + p 01070 { 01071 pPDG=3112; // SigM 01072 mP=mSigM; 01073 mPDG= 211; // Pi+ 01074 mM=mPi; 01075 } 01076 else if(r<0.8) // -> Lamb + Pi+ + n 01077 { 01078 mPDG= 211; // Pi+ 01079 mM=mPi; 01080 NPDG=2112; // n 01081 mT=mNeut; 01082 } 01083 } 01084 else // L/S0-n 01085 { 01086 if(r<0.2) // -> SigM + Pi0 + p 01087 { 01088 pPDG=3112; // SigM 01089 mP=mSigM; 01090 NPDG=2212; // p 01091 mT=mProt; 01092 } 01093 else if(r<0.4) // -> SigP + Pi- + n 01094 { 01095 pPDG=3222; // SigP 01096 mP=mSigP; 01097 mPDG=-211; // Pi- 01098 mM=mPi; 01099 } 01100 else if(r<0.6) // -> SigM + Pi+ + n 01101 { 01102 pPDG=3112; // SigM 01103 mP=mSigM; 01104 mPDG= 211; // Pi+ 01105 mM=mPi; 01106 } 01107 else if(r<0.8) // -> Lamb + Pi- + p 01108 { 01109 mPDG=-211; // Pi- 01110 mM=mPi; 01111 NPDG=2212; // p 01112 mT=mProt; 01113 } 01114 } 01115 } 01116 else if (pPDG == 3112) // Sigma- -N 01117 { 01118 if(Z) // SigM-p 01119 { 01120 if(r<0.25) // -> Lamb + Pi0 + n 01121 { 01122 pPDG=3122; // Lamb 01123 mP=mLamb; 01124 NPDG=2112; // n 01125 mT=mNeut; 01126 } 01127 else if(r<0.5) // -> Lamb + Pi- + p 01128 { 01129 pPDG=3122; // Lamb 01130 mP=mLamb; 01131 mPDG=-211; // Pi- 01132 mM=mPi; 01133 } 01134 else if(r<0.75) // -> SigM + Pi+ + n 01135 { 01136 mPDG= 211; // Pi+ 01137 mM=mPi; 01138 NPDG=2112; // n 01139 mT=mNeut; 01140 } 01141 } 01142 else // SigM-n 01143 { 01144 if(r<third) // -> Lamb + Pi- + n 01145 { 01146 pPDG=3122; // Lamb 01147 mP=mLamb; 01148 mPDG=-211; // Pi- 01149 mM=mPi; 01150 } 01151 else if(r<twothd) // -> SigM + Pi- + p 01152 { 01153 mPDG=-211; // Pi- 01154 mM=mPi; 01155 NPDG=2212; // p 01156 mT=mProt; 01157 } 01158 } 01159 } 01160 else if (pPDG == 3222) // Sigma+ -N 01161 { 01162 if(Z) // SigP-p 01163 { 01164 if(r<third) // -> Lamb + Pi+ + p 01165 { 01166 pPDG=3122; // Lamb 01167 mP=mLamb; 01168 mPDG= 211; // Pi+ 01169 mM=mPi; 01170 } 01171 else if(r<twothd) // -> SigP + Pi+ + n 01172 { 01173 mPDG= 211; // Pi+ 01174 mM=mPi; 01175 NPDG=2112; // n 01176 mT=mNeut; 01177 } 01178 } 01179 else // SigP-n 01180 { 01181 if(r<0.25) // -> Lamb + Pi0 + p 01182 { 01183 pPDG=3122; // Lamb 01184 mP=mLamb; 01185 NPDG=2212; // p 01186 mT=mProt; 01187 } 01188 else if(r<0.5) // -> Lamb + Pi+ + n 01189 { 01190 pPDG=3122; // Lamb 01191 mP=mLamb; 01192 mPDG= 211; // Pi+ 01193 mM=mPi; 01194 } 01195 else if(r<0.75) // -> SigP + Pi- + p 01196 { 01197 mPDG=-211; // Pi- 01198 mM=mPi; 01199 NPDG=2212; // p 01200 mT=mProt; 01201 } 01202 } 01203 } 01204 else if (pPDG == 3312) // Xi- -N 01205 { 01206 if(Z) // XiM-p 01207 { 01208 if(r<0.25) // -> Xi0 + Pi0 + n 01209 { 01210 pPDG=3322; // Xi0 01211 mP=mXiZ; 01212 NPDG=2112; // n 01213 mT=mNeut; 01214 } 01215 else if(r<0.5) // -> Xi0 + Pi- + p 01216 { 01217 pPDG=3322; // Xi0 01218 mP=mXiZ; 01219 mPDG=-211; // Pi- 01220 mM=mPi; 01221 } 01222 else if(r<0.75) // -> XiM + Pi+ + n 01223 { 01224 mPDG= 211; // Pi+ 01225 mM=mPi; 01226 NPDG=2112; // n 01227 mT=mNeut; 01228 } 01229 } 01230 else // XiM-n 01231 { 01232 if(r<third) // -> Xi0 + Pi- + n 01233 { 01234 pPDG=3322; // Xi0 01235 mP=mXiZ; 01236 mPDG=-211; // Pi- 01237 mM=mPi; 01238 } 01239 else if(r<twothd) // -> XiM + Pi- + p 01240 { 01241 mPDG=-211; // Pi- 01242 mM=mPi; 01243 NPDG=2212; // p 01244 mT=mProt; 01245 } 01246 } 01247 } 01248 else if (pPDG == 3322) // Xi0 -N 01249 { 01250 if(Z) // Xi0-p 01251 { 01252 if(r<third) // -> Xi- + Pi+ + p 01253 { 01254 pPDG=3312; // Xi- 01255 mP=mXiM; 01256 mPDG= 211; // Pi+ 01257 mM=mPi; 01258 } 01259 else if(r<twothd) // -> Xi0 + Pi+ + n 01260 { 01261 mPDG= 211; // Pi+ 01262 mM=mPi; 01263 NPDG=2112; // n 01264 mT=mNeut; 01265 } 01266 } 01267 else // Xi0-n 01268 { 01269 if(r<0.25) // -> Xi- + Pi0 + p 01270 { 01271 pPDG=3312; // Xi- 01272 mP=mXiM; 01273 NPDG=2212; // p 01274 mT=mProt; 01275 } 01276 else if(r<0.5) // -> Xi- + Pi+ + n 01277 { 01278 pPDG=3312; // Xi- 01279 mP=mXiM; 01280 mPDG= 211; // Pi+ 01281 mM=mPi; 01282 } 01283 else if(r<0.75) // -> Xi0 + Pi- + p 01284 { 01285 mPDG=-211; // Pi- 01286 mM=mPi; 01287 NPDG=2212; // p 01288 mT=mProt; 01289 } 01290 } 01291 } 01292 else if (pPDG == 3334) // Om- -N 01293 { 01294 if(Z) // OmM-p 01295 { 01296 if(r<0.5) // -> Om- + Pi+ + n 01297 { 01298 mPDG= 211; // Pi+ 01299 mM=mPi; 01300 NPDG=2112; // n 01301 mT=mNeut; 01302 } 01303 } 01304 else // OmM-n 01305 { 01306 if(r<0.5) // -> Om- + Pi- + p 01307 { 01308 mPDG=-211; // Pi- 01309 mM=mPi; 01310 NPDG=2212; // p 01311 mT=mProt; 01312 } 01313 } 01314 } 01315 else 01316 { 01317 G4cout<<"*Error*G4QFreeScattering::InElF: PDG="<<pPDG 01318 <<", while it is defined only for p,n,hyperons(not Omega),pi,K/antiK"<<G4endl; 01319 G4Exception("G4QFreeScattering::InElF:","22",FatalException,"CHIPS_crash"); 01320 } 01321 if (mC-mP-mM-mT <-0.000001) return 0; 01322 G4QHadronVector* TripQH = new G4QHadronVector; // Proto of the Result 01323 G4LorentzVector m4M(0.,0.,0.,mM); 01324 if(mC-mP-mM-mT < 0.000001) // Equal share 01325 { 01326 p4M=(mP/mC)*c4M; 01327 m4M=(mM/mC)*c4M; 01328 N4M=(mT/mC)*c4M; 01329 } 01330 else 01331 { 01332 p4M=G4LorentzVector(0.,0.,0.,mP); 01333 N4M=G4LorentzVector(0.,0.,0.,mT); 01334 if(!G4QHadron(c4M).DecayIn3(p4M,m4M,N4M)) 01335 { 01336 G4ExceptionDescription ed; 01337 ed << "DecayIn3, TotM=" << mC /* << " <? " << mT + mP + mM */ << G4endl; 01338 G4Exception("G4QFreeScattering::InElF()","HAD_CHPS_0027", FatalException, ed); 01339 } 01340 G4QHadron* h1 = new G4QHadron(pPDG, p4M); 01341 TripQH->push_back(h1); // (delete equivalent, responsibility of users) 01342 #ifdef debug 01343 G4cout << "G4QFreeScat::InElF: H1=" << pPDG << p4M << G4endl; 01344 #endif 01345 G4QHadron* h2 = new G4QHadron(mPDG, m4M); 01346 TripQH->push_back(h2); // (delete equivalent, responsibility of users) 01347 #ifdef debug 01348 G4cout << "G4QFreeScat::InElF: H2=" << mPDG << m4M << G4endl; 01349 #endif 01350 G4QHadron* h3 = new G4QHadron(NPDG, N4M); 01351 TripQH->push_back(h3); // (delete equivalent, responsibility of users) 01352 #ifdef debug 01353 G4cout << "G4QFreeScat::InElF: H3=" << NPDG << N4M << G4endl; 01354 #endif 01355 } 01356 return TripQH; // The Result 01357 } // End of Scatter
std::pair< G4LorentzVector, G4LorentzVector > G4QFreeScattering::Scatter | ( | G4int | NPDG, | |
G4LorentzVector | N4M, | |||
G4int | pPDG, | |||
G4LorentzVector | p4M | |||
) |
Definition at line 565 of file G4QFreeScattering.cc.
References FatalException, G4cerr, G4cout, G4endl, G4Exception(), and G4UniformRand.
00567 { 00568 static const G4double mNeut= G4QPDGCode(2112).GetMass(); 00569 static const G4double mProt= G4QPDGCode(2212).GetMass(); 00570 p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units) 00571 N4M/=megaelectronvolt; 00572 G4LorentzVector tot4M=N4M+p4M; 00573 if(pPDG==90001000) pPDG=2212; 00574 if(pPDG==90000001) pPDG=2112; 00575 if(pPDG==91000000) pPDG=3122; 00576 #ifdef ppdebug 00577 G4cout<<"->G4QFR::Scat:p4M="<<p4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl; 00578 #endif 00579 G4double mT=mNeut; 00580 #ifdef ppdebug 00581 G4int Z=0; 00582 G4int N=1; 00583 #endif 00584 if(NPDG==2212 || NPDG==90001000) 00585 { 00586 mT=mProt; 00587 #ifdef ppdebug 00588 Z=1; 00589 N=0; 00590 #endif 00591 } 00592 else if(NPDG!=2112 && NPDG!=90000001) 00593 { 00594 G4cout<<"Error:G4QFreeScattering::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; 00595 G4Exception("G4QFreeScattering::Scatter:","21",FatalException,"CHIPScomplain"); 00596 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Use this if not exception 00597 } 00598 G4double mT2=mT*mT; // a squared mass of the free scattered nuclead cluster (FSNC) 00599 G4double mP2=p4M.m2(); // a projectile squared mass 00600 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); // a projectile energy in the CMS of FSNC 00601 #ifdef pdebug 00602 G4cout<<"G4QFS::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl; 00603 #endif 00604 G4double E2=E*E; 00605 if( E < 0. || E2 < mP2) 00606 { 00607 #ifdef ppdebug 00608 G4cout<<"-Warning-G4QFS::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; 00609 #endif 00610 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 00611 } 00612 G4double pP2=E2-mP2; // Squared Momentum in pseudo laboratory system (final particles) 00613 G4double P=std::sqrt(pP2); // Momentum in pseudo laboratory system (final particles) ? 00614 #ifdef ppdebug 00615 G4cout<<"G4QFreeS::Scatter: Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl; 00616 #endif 00617 // @@ Temporary NN t-dependence for all hadrons 00618 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QFreeScat::Scatter: pPDG="<<pPDG<<G4endl; 00619 // @@ check a possibility to separate p, n, or alpha (!) 00620 G4double dmT=mT+mT; 00621 G4double s_value=dmT*E2+mP2+mT2; // Mondelstam s (?) 00622 G4double maxt=dmT*dmT*pP2/s_value; // max possible |t| 00623 G4double mint=0.; // Prototype of functional rand -t (MeV^2) 00624 if(P < 14.) mint=maxt*G4UniformRand(); // S-wave 00625 else // Calculate slopes (no cash !) 00626 { 00627 G4double p4=pP2*pP2/1000000.; // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ 00628 G4double theB1=(7.2+4.32/p4/(p4+12./P))/(1.+2.5/p4); // p4 in GeV, P in MeV 00629 mint=-std::log(G4UniformRand())/theB1; // t-chan only 00630 } 00631 #ifdef ppdebug 00632 G4cout<<"G4QFS::Scat:PDG="<<pPDG<<",P="<<P<<",-t="<<mint<<"<"<<maxt<<", Z="<<Z<<",N="<<N 00633 <<G4endl; 00634 #endif 00635 #ifdef nandebug 00636 if(mint>-.0000001); 00637 else G4cout<<"*Warning*G4QFreeScattering::Scatter: -t="<<mint<<G4endl; 00638 #endif 00639 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS 00640 #ifdef ppdebug 00641 G4cout<<"G4QFS::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl; 00642 #endif 00643 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) 00644 { 00645 if (cost > 1.) cost = 1.; 00646 else if(cost <-1.) cost =-1.; 00647 else 00648 { 00649 G4cerr<<"G4QFreeScatter::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<maxt<<G4endl; 00650 return std::make_pair(G4LorentzVector(0.,0.,0.,0.), p4M); // Do Nothing Action 00651 } 00652 } 00653 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon 00654 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); 00655 if(!G4QHadron(tot4M).RelDecayIn2(p4M, reco4M, dir4M, cost, cost)) 00656 { 00657 G4cerr<<"G4QFS::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl; 00658 //G4Exception("G4QFS::Scat:","009",FatalException,"Decay of ElasticComp"); 00659 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 00660 } 00661 #ifdef ppdebug 00662 G4cout<<"G4QFS::Scat:p4M="<<p4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl; 00663 #endif 00664 return std::make_pair( reco4M*megaelectronvolt, p4M*megaelectronvolt ); // The Result 00665 } // End of Scatter