#include <G4QDiffractionRatio.hh>
Public Member Functions | |
| ~G4QDiffractionRatio () | |
| G4double | GetRatio (G4double pIU, G4int prPDG, G4int tgZ, G4int tgN) |
| G4QHadronVector * | ProjFragment (G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN) |
| G4QHadronVector * | TargFragment (G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN) |
| G4double | GetTargSingDiffXS (G4double pIU, G4int prPDG, G4int tgZ, G4int tgN) |
Static Public Member Functions | |
| static G4QDiffractionRatio * | GetPointer () |
Protected Member Functions | |
| G4QDiffractionRatio () | |
Definition at line 53 of file G4QDiffractionRatio.hh.
| G4QDiffractionRatio::G4QDiffractionRatio | ( | ) | [inline, protected] |
| G4QDiffractionRatio::~G4QDiffractionRatio | ( | ) | [inline] |
| G4QDiffractionRatio * G4QDiffractionRatio::GetPointer | ( | ) | [static] |
Definition at line 48 of file G4QDiffractionRatio.cc.
Referenced by G4ProjectileDiffractiveChannel::G4ProjectileDiffractiveChannel(), G4QFragmentation::G4QFragmentation(), and G4QDiffraction::PostStepDoIt().
00049 { 00050 static G4QDiffractionRatio theRatios; // *** Static body of the Diffraction Ratio *** 00051 return &theRatios; 00052 }
Definition at line 55 of file G4QDiffractionRatio.cc.
References G4cout, G4endl, and CLHEP::detail::n.
Referenced by G4QFragmentation::G4QFragmentation(), and G4ProjectileDiffractiveChannel::GetFraction().
00056 { 00057 static const G4double mNeut= G4QPDGCode(2112).GetMass()/GeV; // in GeV 00058 static const G4double mProt= G4QPDGCode(2212).GetMass()/GeV; // in GeV 00059 static const G4double mN=.5*(mNeut+mProt); // mean nucleon mass in GeV 00060 static const G4double dmN=mN+mN; // doubled nuc. mass in GeV 00061 static const G4double mN2=mN*mN; // squared nuc. mass in GeV^2 00062 // Table parameters 00063 static const G4int nps=100; // Number of steps in the R(s) LinTable 00064 static const G4int mps=nps+1; // Number of elements in the R(s) LinTable 00065 static const G4double sma=6.; // The first LinTabEl(sqs=0)=1., sqs>sma -> logTab 00066 static const G4double ds=sma/nps; // Step of the linear Table 00067 static const G4int nls=150; // Number of steps in the R(lns) logTable 00068 static const G4int mls=nls+1; // Number of elements in the R(lns) logTable 00069 static const G4double lsi=1.79; // The min ln(sqs) logTabEl(sqs=5.99 < sma=6.) 00070 static const G4double lsa=8.; // The max ln(sqs) logTabEl(sqs=5.99 - 2981 GeV) 00071 static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 5.99 GeV) 00072 static const G4double max_s=std::exp(lsa);// The max s of logTabEl(~ 2981 GeV) 00073 static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table 00074 static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table 00075 static const G4double toler=.0001; // Tolarence (GeV) defining the same sqs 00076 static G4double lastS=0.; // Last sqs value for which R was calculated 00077 static G4double lastR=0.; // Last ratio R which was calculated 00078 // Local Associative Data Base: 00079 static std::vector<G4int> vA; // Vector of calculated A 00080 //static std::vector<G4double> vH; // Vector of max sqs initialized in the LinTable 00081 //static std::vector<G4int> vN; // Vector of topBin number initialized in LinTab 00082 //static std::vector<G4double> vM; // Vector of relMax ln(sqs) initialized in LogTab 00083 //static std::vector<G4int> vK; // Vector of topBin number initialized in LogTab 00084 static std::vector<G4double*> vT; // Vector of pointers to LinTable in C++ heap 00085 static std::vector<G4double*> vL; // Vector of pointers to LogTable in C++ heap 00086 // Last values of the Associative Data Base: 00087 //static G4int lastPDG=0; // Last PDG for which R was calculated (now fake) 00088 static G4int lastA=0; // theLast of calculated A 00089 //static G4double lastH=0.; // theLast max sqs initialized in the LinTable 00090 //static G4int lastN=0; // theLast of topBin number initialized in LinTab 00091 //static G4double lastM=0.; // theLast relMax ln(sqs) initialized in LogTab 00092 //static G4int lastK=0; // theLast of topBin number initialized in LogTab 00093 static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap 00094 static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap 00095 // LogTable is created only if necessary. R(sqs>2981GeV) calcul by formula for any nuclei 00096 G4int A=tgN+tgZ; 00097 if(pIU<toler || A<1) return 1.; // Fake use of toler as non zero number 00098 if(A>238) 00099 { 00100 G4cout<<"-*-Warning-*-G4QuasiFreeRatio::GetRatio: A="<<A<<">238, return zero"<<G4endl; 00101 return 0.; 00102 } 00103 //lastPDG=pPDG; // @@ at present ratio is PDG independent @@ 00104 // Calculate sqs 00105 G4double pM=G4QPDGCode(pPDG).GetMass()/GeV; // Projectile mass in GeV 00106 G4double pM2=pM*pM; 00107 G4double mom=pIU/GeV; // Projectile momentum in GeV 00108 G4double s_value=std::sqrt(mN2+pM2+dmN*std::sqrt(pM2+mom*mom)); // in GeV 00109 G4int nDB=vA.size(); // A number of nuclei already initialized in AMDB 00110 if(nDB && lastA==A && std::fabs(s_value-lastS)<toler) return lastR; 00111 if(s_value>max_s) 00112 { 00113 lastR=CalcDiff2Prod_Ratio(s_value,A); // @@ Probably user ought to be notified about bigS 00114 return lastR; 00115 } 00116 G4bool found=false; 00117 G4int i=-1; 00118 if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB 00119 { 00120 found=true; // The A value is found 00121 break; 00122 } 00123 if(!nDB || !found) // Create new line in the AMDB 00124 { 00125 lastA = A; 00126 lastT = new G4double[mps]; // Create the linear Table 00127 //lastN = static_cast<int>(s_value/ds)+1; // MaxBin to be initialized 00128 //if(lastN>nps) // ===> Now initialize all lin table 00129 //{ 00130 // lastN=nps; 00131 // lastH=sma; 00132 //} 00133 //else lastH = lastN*ds; // Calculate max initialized s for LinTab 00134 G4double sv=0; 00135 lastT[0]=1.; 00136 //for(G4int j=1; j<=lastN; j++) // Calculate LinTab values 00137 for(G4int j=1; j<=nps; j++) // Calculate LinTab values 00138 { 00139 sv+=ds; 00140 lastT[j]=CalcDiff2Prod_Ratio(sv,A); 00141 } 00142 lastL = new G4double[mls]; // Create the logarithmic Table 00143 //G4double ls=std::log(s_value); 00144 //lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB 00145 //if(lastK>nls) // ===> Now initialize all lin table 00146 //{ 00147 // lastK=nls; 00148 // lastM=lsa-lsi; 00149 //} 00150 //else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab 00151 sv=mi; 00152 //for(G4int j=0; j<=lastK; j++) // Calculate LogTab values 00153 for(G4int j=0; j<=nls; j++) // Calculate LogTab values 00154 { 00155 lastL[j]=CalcDiff2Prod_Ratio(sv,A); 00156 //if(j!=lastK) sv*=edl; 00157 sv*=edl; 00158 } 00159 i++; // Make a new record to AMDB and position on it 00160 vA.push_back(lastA); 00161 //vH.push_back(lastH); 00162 //vN.push_back(lastN); 00163 //vM.push_back(lastM); 00164 //vK.push_back(lastK); 00165 vT.push_back(lastT); 00166 vL.push_back(lastL); 00167 } 00168 else // The A value was found in AMDB 00169 { 00170 lastA=vA[i]; 00171 //lastH=vH[i]; 00172 //lastN=vN[i]; 00173 //lastM=vM[i]; 00174 //lastK=vK[i]; 00175 lastT=vT[i]; 00176 lastL=vL[i]; 00177 // ==> Now all bins of the tables are initialized immediately for the A 00178 //if(s_value>lastH) // At least LinTab must be updated 00179 //{ 00180 // G4int nextN=lastN+1; // The next bin to be initialized 00181 // if(lastN<nps) 00182 // { 00183 // lastN = static_cast<int>(s_value/ds)+1;// MaxBin to be initialized 00184 // G4double sv=lastH; 00185 // if(lastN>nps) 00186 // { 00187 // lastN=nps; 00188 // lastH=sma; 00189 // } 00190 // else lastH = lastN*ds; // Calculate max initialized s for LinTab 00191 // for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values 00192 // { 00193 // sv+=ds; 00194 // lastT[j]=CalcDiff2Prod_Ratio(sv,A); 00195 // } 00196 // } // End of LinTab update 00197 // if(lastN>=nextN) 00198 // { 00199 // vH[i]=lastH; 00200 // vN[i]=lastN; 00201 // } 00202 // G4int nextK=lastK+1; 00203 // if(s_value>sma && lastK<nls) // LogTab must be updated 00204 // { 00205 // G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed) 00206 // G4double ls=std::log(s_value); 00207 // lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB 00208 // if(lastK>nls) 00209 // { 00210 // lastK=nls; 00211 // lastM=lsa-lsi; 00212 // } 00213 // else lastM = lastK*dl; // Calcul. max initialized ln(s)-lsi for LogTab 00214 // for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values 00215 // { 00216 // sv*=edl; 00217 // lastL[j]=CalcDiff2Prod_Ratio(sv,A); 00218 // } 00219 // } // End of LogTab update 00220 // if(lastK>=nextK) 00221 // { 00222 // vM[i]=lastM; 00223 // vK[i]=lastK; 00224 // } 00225 //} 00226 } 00227 // Now one can use tabeles to calculate the value 00228 if(s_value<sma) // Use linear table 00229 { 00230 G4int n=static_cast<int>(s_value/ds); // Low edge number of the bin 00231 G4double d=s_value-n*ds; // Linear shift 00232 G4double v=lastT[n]; // Base 00233 lastR=v+d*(lastT[n+1]-v)/ds; // Result 00234 } 00235 else // Use log table 00236 { 00237 G4double ls=std::log(s_value)-lsi; // ln(s)-l_min 00238 G4int n=static_cast<int>(ls/dl); // Low edge number of the bin 00239 G4double d=ls-n*dl; // Log shift 00240 G4double v=lastL[n]; // Base 00241 lastR=v+d*(lastL[n+1]-v)/dl; // Result 00242 } 00243 if(lastR<0.) lastR=0.; 00244 if(lastR>1.) lastR=1.; 00245 return lastR; 00246 } // End of GetRatio
Definition at line 1290 of file G4QDiffractionRatio.cc.
References G4cerr, and G4endl.
01291 { 01292 G4double mom=pIU/gigaelectronvolt; // Projectile momentum in GeV 01293 if ( mom < 1. || (pPDG != 2212 && pPDG != 2112) ) 01294 G4cerr<<"G4QDiffractionRatio::GetTargSingDiffXS isn't applicable p="<<mom<<" GeV, PDG=" 01295 <<pPDG<<G4endl; 01296 G4double A=Z+N; // A of the target 01297 //return 4.5*std::pow(A,.364)*millibarn; // Result 01298 return 3.7*std::pow(A,.364)*millibarn; // Result after mpi0 correction 01299 01300 } // End of ProjFragment
| G4QHadronVector * G4QDiffractionRatio::ProjFragment | ( | G4int | pPDG, | |
| G4LorentzVector | p4M, | |||
| G4int | tgZ, | |||
| G4int | tgN | |||
| ) |
Definition at line 484 of file G4QDiffractionRatio.cc.
References FatalException, G4Quasmon::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4RandomDirection(), G4UniformRand, and G4QNucleus::GetMZNS().
Referenced by G4QFragmentation::G4QFragmentation(), and G4ProjectileDiffractiveChannel::Scatter().
00486 { 00487 static const G4double pFm= 250.; // Fermi momentum in MeV (delta function) 00488 static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function) 00489 static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction) 00490 static const G4double mPi= G4QPDGCode(211).GetMass(); // pi+- mass (MeV) 00491 static const G4double mNeut= G4QPDGCode(2112).GetMass(); 00492 static const G4double mNeut2=mNeut*mNeut; 00493 static const G4double dmNeut=mNeut+mNeut; 00494 static const G4double mProt= G4QPDGCode(2212).GetMass(); 00495 static const G4double mProt2=mProt*mProt; 00496 static const G4double dmProt=mProt+mProt; 00497 static const G4double maxDM=mProt*12.; 00498 static const G4double mLamb= G4QPDGCode(3122).GetMass(); 00499 static const G4double mSigZ= G4QPDGCode(3212).GetMass(); 00500 static const G4double mSigM= G4QPDGCode(3112).GetMass(); 00501 static const G4double mSigP= G4QPDGCode(3222).GetMass(); 00502 static const G4double eps=.003; 00503 // 00504 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) 00505 // prepare the DONOTHING answer 00506 G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !! 00507 G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile 00508 ResHV->push_back(hadron); // It must be cleaned up for real scattering sec 00509 // @@ diffraction is simulated as noncoherent (coherent is small) 00510 G4int tgA=tgZ+tgN; // A of the target 00511 G4int tPDG=90000000+tgZ*1000+tgN; // PDG code of the targetNucleus/recoilNucleus 00512 G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus 00513 G4int rPDG=2112; // prototype of PDG code of the recoiled nucleon 00514 if(tgA*G4UniformRand()>tgN) // Substitute by a proton 00515 { 00516 rPDG=2212; // PDG code of the recoiled QF nucleon 00517 tPDG-=1000; // PDG code of the recoiled nucleus 00518 } 00519 else tPDG-=1; // PDG code of the recoiled nucleus 00520 G4double tM=G4QPDGCode(tPDG).GetMass(); // Mass of the recoiled nucleus 00521 G4double tE=std::sqrt(tM*tM+pFm2); 00522 G4ThreeVector tP=pFm*G4RandomDirection(); 00523 G4LorentzVector t4M(tP,tE); // 4M of the recoil nucleus 00524 G4LorentzVector tg4M(0.,0.,0.,tgM); 00525 G4LorentzVector N4M=tg4M-t4M; // Quasi-free target nucleon 00526 G4LorentzVector tot4M=N4M+p4M; // total momentum of quasi-free diffraction 00527 G4double mT=mNeut; 00528 G4double mT2=mNeut2; // Squared mass of the free nucleon spectator 00529 G4double dmT=dmNeut; 00530 //G4int Z=0; 00531 //G4int N=1; 00532 if(rPDG==2212) 00533 { 00534 mT=mProt; 00535 mT2=mProt2; 00536 dmT=dmProt; 00537 //Z=1; 00538 //N=0; 00539 } 00540 G4double mP2=pr4M.m2(); // Squared mass of the projectile 00541 if(mP2<0.) mP2=0.; // A possible problem for photon (m_min = 2*m_pi0) 00542 G4double s_value=tot4M.m2(); // @@ Check <0 ... 00543 G4double E=(s_value-mT2-mP2)/dmT; // Effective interactin energy (virt. nucl. target) 00544 G4double E2=E*E; 00545 if(E<0. || E2<mP2) 00546 { 00547 #ifdef pdebug 00548 G4cerr<<"-Warning-G4DifR::PFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; 00549 #endif 00550 return ResHV; // Do Nothing Action 00551 } 00552 G4double mP=std::sqrt(mP2); 00553 if(mP<.1)mP=mPi0; // For photons min diffraction is gamma+P->Pi0+Pi0 00554 G4double mMin=mP+mPi0; // Minimum diffractive mass 00555 G4double ss=std::sqrt(s_value); // CM compound mass (sqrt(s)) 00556 G4double mMax=ss-mT; // Maximum diffraction mass 00557 if(mMax>maxDM) mMax=maxDM; // Restriction to avoid too big masses 00558 if(mMin>=mMax) 00559 { 00560 #ifdef pdebug 00561 G4cerr<<"-Warning-G4DifR::PFra:ZeroDiffractionMRange, mi="<<mMin<<",ma="<<mMax<<G4endl; 00562 #endif 00563 return ResHV; // Do Nothing Action 00564 } 00565 G4double R = G4UniformRand(); 00566 G4double mDif=std::exp(R*std::log(mMax)+(1.-R)*std::log(mMin)); // LowMassApproximation 00567 G4double mDif2=mDif*mDif; 00568 G4double ds=s_value-mT2-mDif2; 00569 //G4double e=ds/dmT; 00570 //G4double P=std::sqrt(e*e-mDif2); // Momentum in pseudo laboratory system 00571 #ifdef debug 00572 G4cout<<"G4QDiffR::PFra: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; 00573 #endif 00574 // @@ Temporary NN t-dependence for all hadrons 00575 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl; 00576 G4double tsl=140000.; // slope in MeV^2 00577 G4double t=-std::log(G4UniformRand())*tsl; 00578 G4double maxt=(ds*ds-4*mT2*mDif2)/s_value; // maximum possible -t 00579 #ifdef pdebug 00580 G4cout<<"G4QDifR::PFra:ph="<<pPDG<<",P="<<P<<",t="<<mint<<"<"<<maxt<<G4endl; 00581 #endif 00582 #ifdef nandebug 00583 if(mint>-.0000001); // To make the Warning for NAN 00584 else G4cout<<"******G4QDiffractionRatio::ProjFragment: -t="<<mint<<G4endl; 00585 #endif 00586 G4double rt=t/maxt; 00587 G4double cost=1.-rt-rt; // cos(theta) in CMS 00588 #ifdef ppdebug 00589 G4cout<<"G4QDiffRatio::ProjFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl; 00590 #endif 00591 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) 00592 { 00593 if (cost>1.) cost=1.; 00594 else if(cost<-1.) cost=-1.; 00595 else 00596 { 00597 G4cerr<<"G4QDiffRat::ProjFragm: *NAN* cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl; 00598 return ResHV; // Do Nothing Action 00599 } 00600 } 00601 G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon 00602 G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif); // 4mom of the diffract. Quasmon 00603 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); 00604 if(!G4QHadron(tot4M).RelDecayIn2(d4M, r4M, dir4M, cost, cost)) 00605 { 00606 G4cerr<<"G4QDifR::PFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl; 00607 //G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp"); 00608 return ResHV; // Do Nothing Action 00609 } 00610 #ifdef debug 00611 G4cout<<"G4QDiffR::ProjFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl; 00612 #endif 00613 // Now everything is ready for fragmentation and DoNothing projHadron must be wiped out 00614 delete hadron; // Delete the fake (doNothing) projectile hadron 00615 ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile 00616 hadron = new G4QHadron(tPDG,t4M); // Hadron for the recoil neucleus 00617 ResHV->push_back(hadron); // Fill the recoil nucleus 00618 #ifdef debug 00619 G4cout<<"G4QDiffractionRatio::ProjFragment: *Filled* RecNucleus="<<t4M<<tPDG<<G4endl; 00620 #endif 00621 hadron = new G4QHadron(rPDG,r4M); // Hadron for the recoil nucleon 00622 ResHV->push_back(hadron); // Fill the recoil nucleon 00623 #ifdef debug 00624 G4cout<<"G4QDiffractionRatio::ProjFragment: *Filled* RecNucleon="<<r4M<<rPDG<<G4endl; 00625 #endif 00626 G4LorentzVector sum4M(0.,0.,0.,0.); 00627 // Now the (pPdg,d4M) Quasmon must be fragmented 00628 G4QHadronVector* leadhs = 0; // Prototype of QuasmOutput G4QHadronVector 00629 G4QContent dQC=G4QPDGCode(pPDG).GetQuarkContent(); // Quark Content of the projectile 00630 G4Quasmon* pan= new G4Quasmon(dQC,d4M); // --->---->---->----->-----> DELETED -->---* 00631 try // | 00632 { // | 00633 G4QNucleus vac(90000000); // | 00634 leadhs=pan->Fragment(vac,1); // DELETED after it is copied to ResHV vector -->---+-* 00635 } // | | 00636 catch (G4QException& error) // | | 00637 { // | | 00638 G4cerr<<"***G4QDiffractionRatio::ProjFragment: G4Quasmon Exception"<<G4endl; //| | 00639 // G4Exception("G4QDiffractionRatio::ProjFragment","72",FatalException,"*Quasmon");//| | 00640 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0072", 00641 FatalException, "*Quasmon"); 00642 } // | | 00643 delete pan; // Delete the Nuclear Environment <----<---* | 00644 G4int qNH=leadhs->size(); // A#of collected hadrons from diff.frag. | 00645 if(qNH) for(G4int iq=0; iq<qNH; iq++) // Loop over hadrons to fill the result | 00646 { // | 00647 G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron | 00648 G4int nL=loh->GetStrangeness(); // A number of Lambdas in the Hypernucleus | 00649 G4int nB=loh->GetBaryonNumber(); // Total Baryon Number of the Hypernucleus | 00650 G4int nC = loh->GetCharge(); // Charge of the Hypernucleus | 00651 G4int oPDG = loh->GetPDGCode(); // Original CHIPS PDG Code of the hadron | 00652 //if((nC>nB || nC<0) && nB>0 && nL>=0 && nL<=nB && oPDG>80000000) // Iso-nucleus | 00653 if(2>3) // Closed because "G4QDR::F:90002999,M=-7.768507e-04,B=2,S=0,C=3" is found | 00654 { 00655 G4LorentzVector q4M = loh->Get4Momentum(); // Get 4-momentum of the Isonucleus | 00656 G4double qM=q4M.m(); // Real mass of the Isonucleus 00657 #ifdef fdebug 00658 G4cout<<"G4QDR::PF:"<<oPDG<<",M="<<qM<<",B="<<nB<<",S="<<nL<<",C="<<nC<<G4endl;// | 00659 #endif 00660 G4int qPN=nC-nB; // Number of pions in the Isonucleus | 00661 G4int fPDG = 2212; // Prototype for nP+(Pi+) case | 00662 G4int sPDG = 211; 00663 tPDG = 3122; // @@ Sigma0 (?) | 00664 G4double fMass= mProt; 00665 G4double sMass= mPi; 00666 G4double tMass= mLamb; // @@ Sigma0 (?) | 00667 G4bool cont=true; // Continue flag | 00668 // =--------= Negative state =---------= 00669 if(nC<0) // =----= Only Pi- can help | 00670 { 00671 if(nL&&nB==nL) // --- n*Lamb + k*(Pi-) State --- | 00672 { 00673 sPDG = -211; 00674 if(-nC==nL && nL==1) // Only one Sigma- like (nB=1) | 00675 { 00676 if(std::fabs(qM-mSigM)<eps) 00677 { 00678 loh->SetQPDG(G4QPDGCode(3112)); // This is Sigma- | 00679 cont=false; // Skip decay | 00680 } 00681 else if(qM>mLamb+mPi) //(2) Sigma- => Lambda + Pi- decay | 00682 { 00683 fPDG = 3122; 00684 fMass= mLamb; 00685 } 00686 else if(qM>mSigM) //(2) Sigma+=>Sigma++gamma decay | 00687 { 00688 fPDG = 3112; 00689 fMass= mSigM; 00690 sPDG = 22; 00691 sMass= 0.; 00692 } 00693 else //(2) Sigma-=>Neutron+Pi- decay | 00694 { 00695 fPDG = 2112; 00696 fMass= mNeut; 00697 } 00698 qPN = 1; // #of (Pi+ or gamma)'s = 1 | 00699 } 00700 else if(-nC==nL) //(2) a few Sigma- like | 00701 { 00702 qPN = 1; // One separated Sigma- | 00703 fPDG = 3112; 00704 sPDG = 3112; 00705 sMass= mSigM; 00706 nB--; 00707 fMass= mSigM; 00708 } 00709 else if(-nC>nL) //(2) n*(Sigma-)+m*(Pi-) | 00710 { 00711 qPN = -nC-nL; // #of Pi-'s | 00712 fPDG = 3112; 00713 fMass= mSigM; 00714 } 00715 else //(2) n*(Sigma-)+m*Lambda(-nC<nL) | 00716 { 00717 nB += nC; // #of Lambda's | 00718 fPDG = 3122; 00719 fMass= mLamb; 00720 qPN = -nC; // #of Sigma+'s | 00721 sPDG = 3112; 00722 sMass= mSigM; 00723 } 00724 nL = 0; // Only decays in two are above | 00725 } 00726 else if(nL) // ->n*Lamb+m*Neut+k*(Pi-) State (nL<nB) | 00727 { 00728 nB -= nL; // #of neutrons | 00729 fPDG = 2112; 00730 fMass= mNeut; 00731 G4int nPin = -nC; // #of Pi-'s 00732 if(nL==nPin) //(2) m*Neut+n*Sigma- | 00733 { 00734 qPN = nL; // #of Sigma- | 00735 sPDG = 3112; 00736 sMass= mSigM; 00737 nL = 0; 00738 } 00739 else if(nL>nPin) //(3) m*P+n*(Sigma+)+k*Lambda | 00740 { 00741 nL-=nPin; // #of Lambdas | 00742 qPN = nPin; // #of Sigma+ | 00743 sPDG = 3112; 00744 sMass= mSigM; 00745 } 00746 else //(3) m*N+n*(Sigma-)+k*(Pi-) (nL<nPin) | 00747 { 00748 qPN = nPin-nL; // #of Pi- | 00749 sPDG = -211; 00750 tPDG = 3112; 00751 tMass= mSigM; 00752 } 00753 } 00754 else //(2) n*N+m*(Pi-) (nL=0) | 00755 { 00756 sPDG = -211; 00757 qPN = -nC; 00758 fPDG = 2112; 00759 fMass= mNeut; 00760 } 00761 } 00762 else if(!nC) // *** Should not be here *** | 00763 { 00764 if(nL && nL<nB) //(2) n*Lamb+m*N ***Should not be here*** | 00765 { 00766 qPN = nL; 00767 fPDG = 2112; // mN+nL case | 00768 sPDG = 3122; 00769 sMass= mLamb; 00770 nB -= nL; 00771 fMass= mNeut; 00772 nL = 0; 00773 } 00774 else if(nL>1 && nB==nL) //(2) m*Lamb(m>1) ***Should not be here*** | 00775 { 00776 qPN = 1; 00777 fPDG = 3122; 00778 sPDG = 3122; 00779 sMass= mLamb; 00780 nB--; 00781 fMass= mLamb; 00782 } 00783 else if(!nL && nB>1) //(2) n*Neut(n>1) ***Should not be here*** | 00784 { 00785 qPN = 1; 00786 fPDG = 2112; 00787 sPDG = 2112; 00788 sMass= mNeut; 00789 nB--; 00790 fMass= mNeut; 00791 } 00792 else G4cout<<"*?*G4QDiffractionRatio::ProjFragment: (1) oPDG="<<oPDG<<G4endl;// | 00793 } 00794 else if(nC>0) // n*Lamb+(m*P)+(k*Pi+) | 00795 { 00796 if(nL && nL+nC==nB) //(2) n*Lamb+m*P ***Should not be here*** | 00797 { 00798 qPN = nL; 00799 nL = 0; 00800 fPDG = 2212; 00801 sPDG = 3122; 00802 sMass= mLamb; 00803 nB = nC; 00804 fMass= mProt; 00805 } 00806 else if(nL && nC<nB-nL) //(3)n*L+m*P+k*N ***Should not be here*** | 00807 { 00808 qPN = nC; // #of protons | 00809 fPDG = 2112; // mP+nL case | 00810 sPDG = 2212; 00811 sMass= mProt; 00812 nB -= nL+nC; // #of neutrons | 00813 fMass= mNeut; 00814 } 00815 else if(nL && nB==nL) // ---> n*L+m*Pi+ State | 00816 { 00817 if(nC==nL && nL==1) // Only one Sigma+ like State | 00818 { 00819 if(std::fabs(qM-mSigP)<eps) 00820 { 00821 loh->SetQPDG(G4QPDGCode(3222)); // This is GS Sigma+ | 00822 cont=false; // Skip decay | 00823 } 00824 else if(qM>mLamb+mPi) //(2) Sigma+=>Lambda+Pi+ decay | 00825 { 00826 fPDG = 3122; 00827 fMass= mLamb; 00828 } 00829 else if(qM>mNeut+mPi) //(2) Sigma+=>Neutron+Pi+ decay | 00830 { 00831 fPDG = 2112; 00832 fMass= mNeut; 00833 } 00834 else if(qM>mSigP) //(2) Sigma+=>Sigma++gamma decay | 00835 { 00836 fPDG = 3222; 00837 fMass= mSigP; 00838 sPDG = 22; 00839 sMass= 0.; 00840 } 00841 else //(2) Sigma+=>Proton+gamma decay | 00842 { 00843 fPDG = 2212; 00844 fMass= mProt; 00845 sPDG = 22; 00846 sMass= 0.; 00847 } 00848 qPN = 1; // #of (Pi+ or gamma)'s = 1 | 00849 } 00850 else if(nC==nL) //(2) a few Sigma+ like hyperons | 00851 { 00852 qPN = 1; 00853 fPDG = 3222; 00854 sPDG = 3222; 00855 sMass= mSigP; 00856 nB--; 00857 fMass= mSigP; 00858 } 00859 else if(nC>nL) //(2) n*(Sigma+)+m*(Pi+) | 00860 { 00861 qPN = nC-nL; // #of Pi+'s | 00862 fPDG = 3222; 00863 nB = nL; // #of Sigma+'s | 00864 fMass= mSigP; 00865 } 00866 else //(2) n*(Sigma+)+m*Lambda | 00867 { 00868 nB -= nC; // #of Lambda's | 00869 fPDG = 3122; 00870 fMass= mLamb; 00871 qPN = nC; // #of Sigma+'s | 00872 sPDG = 3222; 00873 sMass= mSigP; 00874 } 00875 nL = 0; // All above are decays in 2 | 00876 } 00877 else if(nL && nC>nB-nL) // n*Lamb+m*P+k*Pi+ | 00878 { 00879 nB -= nL; // #of protons | 00880 G4int nPip = nC-nB; // #of Pi+'s | 00881 if(nL==nPip) //(2) m*P+n*Sigma+ | 00882 { 00883 qPN = nL; // #of Sigma+ | 00884 sPDG = 3222; 00885 sMass= mSigP; 00886 nL = 0; 00887 } 00888 else if(nL>nPip) //(3) m*P+n*(Sigma+)+k*Lambda | 00889 { 00890 nL -= nPip; // #of Lambdas | 00891 qPN = nPip; // #of Sigma+ | 00892 sPDG = 3222; 00893 sMass= mSigP; 00894 } 00895 else //(3) m*P+n*(Sigma+)+k*(Pi+) | 00896 { 00897 qPN = nPip-nL; // #of Pi+ | 00898 tPDG = 3222; 00899 tMass= mSigP; 00900 } 00901 } 00902 if(nC<nB) //(2) n*P+m*N ***Should not be here*** | 00903 { 00904 fPDG = 2112; 00905 fMass= mNeut; 00906 qPN = nC; 00907 sPDG = 2212; 00908 sMass= mProt; 00909 } 00910 else if(nB==nC && nC>1) //(2) m*Prot(m>1) ***Should not be here*** | 00911 { 00912 qPN = 1; 00913 fPDG = 2212; 00914 sPDG = 2212; 00915 sMass= mProt; 00916 nB--; 00917 fMass= mProt; 00918 } 00919 else if(nC<=nB||!nB) G4cout<<"*?*G4QDR::ProjFragm: (2) oPDG="<<oPDG<<G4endl; // | 00920 // !nL && nC>nB //(2) Default condition n*P+m*(Pi+) | 00921 } 00922 if(cont) // Make a decay | 00923 { 00924 G4double tfM=nB*fMass; 00925 G4double tsM=qPN*sMass; 00926 G4double ttM=0.; 00927 if(nL) ttM=nL*tMass; 00928 G4LorentzVector f4Mom(0.,0.,0.,tfM); 00929 G4LorentzVector s4Mom(0.,0.,0.,tsM); 00930 G4LorentzVector t4Mom(0.,0.,0.,ttM); 00931 G4double sum=tfM+tsM+ttM; 00932 if(std::fabs(qM-sum)<eps) 00933 { 00934 f4Mom=q4M*(tfM/sum); 00935 s4Mom=q4M*(tsM/sum); 00936 if(nL) t4Mom=q4M*(ttM/sum); 00937 } 00938 else if(!nL && (qM<sum || !G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))) // Error | 00939 { 00940 //#ifdef fdebug 00941 G4cout<<"***G4QDR::PrFragm:fPDG="<<fPDG<<"*"<<nB<<"(fM="<<fMass<<")+sPDG="<<sPDG 00942 <<"*"<<qPN<<"(sM="<<sMass<<")"<<"="<<sum<<" > TM="<<qM<<q4M<<oPDG<<G4endl; 00943 //#endif 00944 // throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 2"); // | 00945 G4ExceptionDescription ed; 00946 ed << "***G4QDR::PrFragm:fPDG=" << fPDG << "*" << nB << "(fM=" 00947 << fMass << ")+sPDG=" << sPDG << "*" << qPN << "(sM=" << sMass 00948 << ")" << "=" << sum << " > TM=" << qM << q4M << oPDG << G4endl; 00949 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0002", 00950 FatalException, ed); 00951 } 00952 else if(nL && (qM<sum || !G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom)))// Error| 00953 { 00954 //#ifdef fdebug 00955 G4cout<<"***G4DF::PrFrag: "<<fPDG<<"*"<<nB<<"("<<fMass<<")+"<<sPDG<<"*"<<qPN<<"(" 00956 <<sMass<<")+Lamb*"<<nL<<"="<<sum<<" > TotM="<<qM<<q4M<<oPDG<<G4endl; 00957 //#endif 00958 // throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 3"); // | 00959 G4ExceptionDescription ed; 00960 ed << "***G4DF::PrFrag: " << fPDG << "*" << nB << "(" << fMass << ")+" 00961 << sPDG << "*" << qPN << "(" << sMass << ")+Lamb*" << nL << "=" 00962 << sum << " > TotM=" << qM << q4M << oPDG << G4endl; 00963 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0003", 00964 FatalException, ed); 00965 } 00966 #ifdef fdebug 00967 G4cout<<"G4QDF::ProjFragm: *DONE* n="<<nB<<f4Mom<<fPDG<<", m="<<qPN<<s4Mom<<sPDG 00968 <<", l="<<nL<<t4Mom<<G4endl; 00969 #endif 00970 G4bool notused=true; 00971 if(nB) // There are baryons | 00972 { 00973 f4Mom/=nB; 00974 loh->Set4Momentum(f4Mom); // ! Update the Hadron ! | 00975 loh->SetQPDG(G4QPDGCode(fPDG)); // Baryons | 00976 notused=false; // Loh was used | 00977 if(nB>1) for(G4int ih=1; ih<nB; ih++) // Loop over the rest of baryons | 00978 { 00979 G4QHadron* Hi = new G4QHadron(fPDG,f4Mom); // Create a Hadron for Baryon | 00980 ResHV->push_back(Hi); // Fill in the additional nucleon | 00981 #ifdef fdebug 00982 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 00983 G4cout<<"G4QDR::ProjFrag: *additional Nucleon*="<<f4Mom<<fPDG<<G4endl; // | 00984 #endif 00985 } 00986 } 00987 if(qPN) // There are pions | 00988 { 00989 s4Mom/=qPN; 00990 G4int min=0; 00991 if(notused) 00992 { 00993 loh->Set4Momentum(s4Mom); // ! Update the Hadron 4M ! | 00994 loh->SetQPDG(G4QPDGCode(sPDG)); // Update PDG | 00995 notused=false; // loh was used | 00996 min=1; // start value | 00997 } 00998 if(qPN>min) for(G4int ip=min; ip<qPN; ip++) // Loop over pions | 00999 { 01000 G4QHadron* Hj = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the meson | 01001 ResHV->push_back(Hj); // Fill in the additional pion | 01002 #ifdef fdebug 01003 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01004 G4cout<<"G4QDR::ProjFragm: *additional Pion*="<<f4Mom<<fPDG<<G4endl; // | 01005 #endif 01006 } 01007 } 01008 if(nL) // There are Hyperons | 01009 { 01010 t4Mom/=nL; 01011 G4int min=0; 01012 if(notused) 01013 { 01014 loh->Set4Momentum(t4Mom); // ! Update the Hadron 4M ! | 01015 loh->SetQPDG(G4QPDGCode(tPDG));// Update PDG | 01016 notused=false; // loh was used | 01017 min=1; // 01018 } 01019 if(nL>min) for(G4int il=min; il<nL; il++) // Loop over Hyperons | 01020 { 01021 G4QHadron* Hk = new G4QHadron(tPDG,t4Mom); // Create a Hadron for Lambda | 01022 ResHV->push_back(Hk); // Fill in the additional pion | 01023 #ifdef fdebug 01024 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01025 G4cout<<"G4QDR::ProjFragm: *additional Hyperon*="<<f4Mom<<fPDG<<G4endl; // | 01026 #endif 01027 } 01028 } 01029 } // --> End of decay | 01030 } // -> End of Iso-nuclear treatment | 01031 else if( (nL > 0 && nB > 1) || (nL < 0 && nB < -1) ) 01032 { // Hypernucleus is found | 01033 G4bool anti=false; // Default=Nucleus (true=antinucleus | 01034 if(nB<0) // Anti-nucleus | 01035 { 01036 anti=true; // Flag of anti-hypernucleus | 01037 nB=-nB; // Reverse the baryon number | 01038 nC=-nC; // Reverse the charge | 01039 nL=-nL; // Reverse the strangeness | 01040 } 01041 G4int hPDG = 90000000+nL*999999+nC*999+nB; // CHIPS PDG Code for Hypernucleus | 01042 G4int nSM=0; // A#0f unavoidable Sigma- | 01043 G4int nSP=0; // A#0f unavoidable Sigma+ | 01044 if(nC<0) // Negative hypernucleus | 01045 { 01046 if(-nC<=nL) // Partial compensation by Sigma- | 01047 { 01048 nSM=-nC; // Can be compensated by Sigma- | 01049 nL+=nC; // Reduce the residual strangeness | 01050 } 01051 else // All Charge is compensated by Sigma- | 01052 { 01053 nSM=nL; // The maximum number of Sigma- | 01054 nL=0; // Kill the residual strangeness | 01055 } 01056 } 01057 else if(nC>nB-nL) // Extra positive hypernucleus | 01058 { 01059 if(nC<=nB) // Partial compensation by Sigma+ | 01060 { 01061 G4int dH=nB-nC; // Isotopic shift | 01062 nSP=nL-dH; // Can be compensated by Sigma+ | 01063 nL=dH; // Reduce the residual strangeness | 01064 } 01065 else // All Charge is compensated by Sigma+ | 01066 { 01067 nSP=nL; // The maximum number of Sigma+ | 01068 nL=0; // Kill the residual strangeness | 01069 } 01070 } 01071 r4M=loh->Get4Momentum(); // Real 4-momentum of the hypernucleus ! 01072 G4double reM=r4M.m(); // Real mass of the hypernucleus | 01073 #ifdef fdebug 01074 G4cout<<"G4QDiffRatio::PrFrag:oPDG=="<<oPDG<<",hPDG="<<hPDG<<",M="<<reM<<G4endl;//| 01075 #endif 01076 G4int rlPDG=hPDG-nL*1000000-nSP*1000999-nSM*999001;// Subtract Lamb/Sig from Nucl.| 01077 G4int sPDG=3122; // Prototype for the Hyperon PDG (Lambda)| 01078 G4double MLa=mLamb; // Prototype for one Hyperon decay | 01079 #ifdef fdebug 01080 G4cout<<"G4QDiffRatio::PrFrag:*G4*nS+="<<nSP<<",nS-="<<nSM<<",nL="<<nL<<G4endl;// | 01081 #endif 01082 if(nSP||nSM) // Sigma+/- improvement | 01083 { 01084 if(nL) // By mistake Lambda improvement is found | 01085 { 01086 G4cout<<"***G4QDR::PFr:HypN="<<hPDG<<": bothSigm&Lamb -> ImproveIt"<<G4endl;//| 01087 //throw G4QException("*G4QDiffractionRatio::Fragment:BothLambda&SigmaInHN");//| 01088 // @@ Correction, which does not conserv the charge !! (-> add decay in 3) | 01089 if(nSP) nL+=nSP; // Convert Sigma+ to Lambda | 01090 else nL+=nSM; // Convert Sigma- to Lambda | 01091 } 01092 if(nSP) // Sibma+ should be decayed | 01093 { 01094 nL=nSP; // #of decaying hyperons | 01095 sPDG=3222; // PDG code of decaying hyperons | 01096 MLa=mSigP; // Mass of decaying hyperons | 01097 } 01098 else // Sibma+ should be decayed | 01099 { 01100 nL=nSM; // #of decaying hyperons | 01101 sPDG=3112; // PDG code of decaying hyperons | 01102 MLa=mSigM; // Mass of decaying hyperons | 01103 } 01104 } 01105 #ifdef fdebug 01106 G4cout<<"G4QDiffRat::ProjFrag:*G4*mS="<<MLa<<",sPDG="<<sPDG<<",nL="<<nL<<G4endl;//| 01107 #endif 01108 if(nL>1) MLa*=nL; // Total mass of the decaying hyperons | 01109 G4double rlM=G4QNucleus(rlPDG).GetMZNS();// Mass of the NonstrangeNucleus | 01110 if(!nSP&&!nSM&&nL==1&&reM>rlM+mSigZ&&G4UniformRand()>.5) // Conv Lambda->Sigma0 | 01111 { 01112 sPDG=3212; // PDG code of a decaying hyperon | 01113 MLa=mSigZ; // Mass of the decaying hyperon | 01114 } 01115 G4int rnPDG = hPDG-nL*999999; // Convert Lambdas to neutrons (for convInN) | 01116 G4QNucleus rnN(rnPDG); // New nonstrange nucleus | 01117 G4double rnM=rnN.GetMZNS(); // Mass of the new nonstrange nucleus | 01118 // @@ In future take into account Iso-Hypernucleus (Add PI+,R & Pi-,R decays) | 01119 if(rlPDG==90000000) // Multy Hyperon (HyperNuc of only hyperons) | 01120 { 01121 if(nL>1) r4M=r4M/nL; // split the 4-mom for the MultyLambda | 01122 for(G4int il=0; il<nL; il++) // loop over Lambdas | 01123 { 01124 if(anti) sPDG=-sPDG; // For anti-nucleus case | 01125 G4QHadron* theLam = new G4QHadron(sPDG,r4M); // Make NewHadr for the Hyperon | 01126 ResHV->push_back(theLam); // Fill in the Lambda | 01127 #ifdef fdebug 01128 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01129 G4cout<<"G4QDR::ProjFrag: *additional Lambda*="<<r4M<<sPDG<<G4endl; // | 01130 #endif 01131 } 01132 } 01133 else if(reM>rlM+MLa-eps) // Lambda (or Sigma) can be split | 01134 { 01135 G4LorentzVector n4M(0.,0.,0.,rlM); // 4-mom of the residual nucleus | 01136 G4LorentzVector h4M(0.,0.,0.,MLa); // 4-mom of the Hyperon | 01137 G4double sum=rlM+MLa; // Safety sum | 01138 if(std::fabs(reM-sum)<eps) // At rest in CMS | 01139 { 01140 n4M=r4M*(rlM/sum); // Split tot 4-mom for resNuc | 01141 h4M=r4M*(MLa/sum); // Split tot 4-mom for Hyperon | 01142 } 01143 else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay | 01144 { 01145 G4cerr<<"***G4QDF::PF:HypN,M="<<reM<<"<A+n*L="<<sum<<",d="<<sum-reM<<G4endl;//| 01146 // throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclusDecay");//| 01147 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0100", 01148 FatalException, "Error in hypernuclus decay"); 01149 } 01150 #ifdef fdebug 01151 G4cout<<"*G4QDR::PF:HypN="<<r4M<<"->A="<<rlPDG<<n4M<<",n*L="<<nL<<h4M<<G4endl;//| 01152 #endif 01153 loh->Set4Momentum(n4M); // ! Update the Hadron ! | 01154 if(anti && rlPDG==90000001) rlPDG=-2112; // Convert to anti-neutron | 01155 if(anti && rlPDG==90001000) rlPDG=-2212; // Convert to anti-proton | 01156 loh->SetQPDG(G4QPDGCode(rlPDG)); // ConvertedHypernucleus to nonstrange(@anti)| 01157 if(rlPDG==90000002) // Additional action with loH changed to 2n | 01158 { 01159 G4LorentzVector newLV=n4M/2.; // Split 4-momentum | 01160 loh->Set4Momentum(newLV); // Reupdate the hadron | 01161 if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG | 01162 else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG | 01163 G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron | 01164 ResHV->push_back(secHadr); // Fill in the additional neutron | 01165 #ifdef fdebug 01166 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01167 G4cout<<"G4QDR::ProgFrag: *additional Neutron*="<<r4M<<sPDG<<G4endl; // | 01168 #endif 01169 } 01170 else if(rlPDG==90002000) // Additional action with loH change to 2p | 01171 { 01172 G4LorentzVector newLV=n4M/2.; // Split 4-momentum | 01173 loh->Set4Momentum(newLV); // Reupdate the hadron | 01174 if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG | 01175 else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG | 01176 G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton | 01177 ResHV->push_back(secHadr); // Fill in the additional neutron | 01178 #ifdef fdebug 01179 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01180 G4cout<<"G4QDR::ProjFrag: *additional Proton*="<<r4M<<sPDG<<G4endl; // | 01181 #endif 01182 } 01183 // @@(?) Add multybaryon decays if necessary (Now it anyhow is made later) | 01184 #ifdef fdebug 01185 G4cout<<"*G4QDiffractionRatio::PrFrag:resNucPDG="<<loh->GetPDGCode()<<G4endl;// | 01186 #endif 01187 if(nL>1) h4M=h4M/nL; // split the lambda's 4-mom if necessary | 01188 for(G4int il=0; il<nL; il++) // A loop over excessive hyperons | 01189 { 01190 if(anti) sPDG=-sPDG; // For anti-nucleus case | 01191 G4QHadron* theLamb = new G4QHadron(sPDG,h4M); // Make NewHadr for the Hyperon | 01192 ResHV->push_back(theLamb); // Fill in the additional neutron | 01193 #ifdef fdebug 01194 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01195 G4cout<<"G4QDR::ProjFrag: *additional Hyperon*="<<r4M<<sPDG<<G4endl; // | 01196 #endif 01197 } 01198 } 01199 else if(reM>rnM+mPi0-eps&&!nSP&&!nSM)// Lambda->N only if Sigmas are absent | 01200 { 01201 G4int nPi=static_cast<G4int>((reM-rnM)/mPi0); // Calc. pion multiplicity | 01202 if (nPi>nL) nPi=nL; // Cut the pion multiplicity | 01203 G4double npiM=nPi*mPi0; // Total pion mass | 01204 G4LorentzVector n4M(0.,0.,0.,rnM); // Residual nucleus 4-momentum | 01205 G4LorentzVector h4M(0.,0.,0.,npiM);// 4-momentum of pions | 01206 G4double sum=rnM+npiM; // Safety sum | 01207 if(std::fabs(reM-sum)<eps) // At rest | 01208 { 01209 n4M=r4M*(rnM/sum); // The residual nucleus part | 01210 h4M=r4M*(npiM/sum); // The pion part | 01211 } 01212 else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay | 01213 { 01214 G4cerr<<"*G4QDR::PF:HypN,M="<<reM<<"<A+n*Pi0="<<sum<<",d="<<sum-reM<<G4endl;//| 01215 // throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclDecay"); // | 01216 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0101", 01217 FatalException, "Error in HypernuclDecay"); 01218 } 01219 loh->Set4Momentum(n4M); // ! Update the Hadron ! | 01220 if(anti && rnPDG==90000001) rnPDG=-2112; // Convert to anti-neutron | 01221 if(anti && rnPDG==90001000) rnPDG=-2212; // Convert to anti-proton | 01222 loh->SetQPDG(G4QPDGCode(rnPDG)); // convert hyperNuc to nonstrangeNuc(@@anti) | 01223 #ifdef fdebug 01224 G4cout<<"*G4QDR::PF:R="<<r4M<<"->A="<<rnPDG<<n4M<<",n*Pi0="<<nPi<<h4M<<G4endl;//| 01225 #endif 01226 if(nPi>1) h4M=h4M/nPi; // Split the 4-mom if necessary | 01227 for(G4int ihn=0; ihn<nPi; ihn++) // A loop over additional pions | 01228 { 01229 G4QHadron* thePion = new G4QHadron(111,h4M); // Make a New Hadr for the pi0 | 01230 ResHV->push_back(thePion); // Fill in the Pion | 01231 #ifdef fdebug 01232 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01233 G4cout<<"G4QDR::ProjFrag: *additional Pion*="<<r4M<<sPDG<<G4endl; // | 01234 #endif 01235 } 01236 if(rnPDG==90000002) // Additional action with loH change to 2n | 01237 { 01238 G4LorentzVector newLV=n4M/2.; // Split 4-momentum | 01239 loh->Set4Momentum(newLV); // Reupdate the hadron | 01240 if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG | 01241 else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG | 01242 G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron | 01243 ResHV->push_back(secHadr); // Fill in the additional neutron | 01244 #ifdef fdebug 01245 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01246 G4cout<<"G4QDR::ProjFrag: *additional Neutron*="<<r4M<<sPDG<<G4endl; // | 01247 #endif 01248 } 01249 else if(rnPDG==90002000) // Additional action with loH change to 2p | 01250 { 01251 G4LorentzVector newLV=n4M/2.; // Split 4-momentum | 01252 loh->Set4Momentum(newLV); // Reupdate the hadron | 01253 if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG | 01254 else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG | 01255 G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton | 01256 ResHV->push_back(secHadr); // Fill in the additional neutron | 01257 #ifdef fdebug 01258 sum4M+=r4M; // Sum 4-momenta for the EnMom check | 01259 G4cout<<"G4QDR::ProjFrag: *additional Proton*="<<r4M<<sPDG<<G4endl; // | 01260 #endif 01261 } 01262 // @@ Add multybaryon decays if necessary | 01263 } 01264 else // If this Excepton shows up (lowProbable appearance) => include gamma decay | 01265 { 01266 G4double d=rlM+MLa-reM; // Hyperon Excessive energy | 01267 G4cerr<<"G4QDR::PF:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d="<<d<<G4endl; 01268 d=rnM+mPi0-reM; // Pion Excessive energy | 01269 G4cerr<<"G4QDR::PF:"<<oPDG<<","<<hPDG<<",M="<<reM<<"<"<<rnM+mPi0<<",d="<<d<<G4endl; 01270 // throw G4QException("G4QDiffractionRatio::ProjFragment: Hypernuclear conver");// | 01271 G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0102", 01272 FatalException, "Excessive hypernuclear energy"); 01273 } 01274 } // => End of G4 Hypernuclear decay | 01275 ResHV->push_back(loh); // Fill in the result | 01276 #ifdef debug 01277 sum4M+=loh->Get4Momentum(); // Sum 4-momenta for the EnMom check | 01278 G4cout<<"G4QDR::PrFra:#"<<iq<<","<<loh->Get4Momentum()<<loh->GetPDGCode()<<G4endl;//| 01279 #endif 01280 } // | 01281 leadhs->clear();// | 01282 delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<--* 01283 #ifdef debug 01284 G4cout<<"G4QDiffractionRatio::ProjFragment: *End* Sum="<<sum4M<<" =?= d4M="<<d4M<<G4endl; 01285 #endif 01286 return ResHV; // Result 01287 } // End of ProjFragment
| G4QHadronVector * G4QDiffractionRatio::TargFragment | ( | G4int | pPDG, | |
| G4LorentzVector | p4M, | |||
| G4int | tgZ, | |||
| G4int | tgN | |||
| ) |
Definition at line 303 of file G4QDiffractionRatio.cc.
References G4QEnvironment::AddQuasmon(), FatalException, G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4RandomDirection(), and G4UniformRand.
Referenced by G4QDiffraction::PostStepDoIt().
00305 { 00306 static const G4double pFm= 0.; // Fermi momentum in MeV (delta function) 00307 //static const G4double pFm= 250.; // Fermi momentum in MeV (delta function) 00308 static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function) 00309 static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction) 00310 //static const G4double mPi= G4QPDGCode(211).GetMass(); // pi+- mass (MeV) 00311 static const G4double mNeut= G4QPDGCode(2112).GetMass(); 00312 static const G4double mNeut2=mNeut*mNeut; 00313 static const G4double dmNeut=mNeut+mNeut; 00314 static const G4double mProt= G4QPDGCode(2212).GetMass(); 00315 static const G4double mProt2=mProt*mProt; 00316 static const G4double dmProt=mProt+mProt; 00317 static const G4double maxDM=mProt*12.; 00318 //static const G4double mLamb= G4QPDGCode(3122).GetMass(); 00319 //static const G4double mSigZ= G4QPDGCode(3212).GetMass(); 00320 //static const G4double mSigM= G4QPDGCode(3112).GetMass(); 00321 //static const G4double mSigP= G4QPDGCode(3222).GetMass(); 00322 //static const G4double eps=.003; 00323 static const G4double third=1./3.; 00324 // 00325 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) 00326 // prepare the DONOTHING answer 00327 G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !! 00328 G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile 00329 ResHV->push_back(hadron); // It must be cleaned up for real scattering sec 00330 // @@ diffraction is simulated as noncoherent (coherent is small) 00331 G4int tgA=tgZ+tgN; // A of the target 00332 G4int tPDG=90000000+tgZ*1000+tgN; // PDG code of the targetNucleus/recoilNucleus 00333 G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus 00334 G4int rPDG=2112; // prototype of PDG code of the recoiled nucleon 00335 if(tgA*G4UniformRand()>tgN) // Substitute by a proton 00336 { 00337 rPDG=2212; // PDG code of the recoiled QF nucleon 00338 tPDG-=1000; // PDG code of the recoiled nucleus 00339 } 00340 else tPDG-=1; // PDG code of the recoiled nucleus 00341 G4double tM=G4QPDGCode(tPDG).GetMass(); // Mass of the recoiled nucleus 00342 G4double tE=std::sqrt(tM*tM+pFm2); // Free energy of the recoil nucleus 00343 G4ThreeVector tP=pFm*G4RandomDirection();// 3-mom of the recoiled nucleus 00344 G4LorentzVector t4M(tP,tE); // 4M of the recoil nucleus 00345 G4LorentzVector tg4M(0.,0.,0.,tgM); // Full target 4-momentum 00346 G4LorentzVector N4M=tg4M-t4M; // 4-mom of Quasi-free target nucleon 00347 G4LorentzVector tot4M=N4M+p4M; // total momentum of quasi-free diffraction 00348 G4double mT=mNeut; // Prototype of mass of QF nucleon 00349 G4double mT2=mNeut2; // Squared mass of a free nucleon to be excited 00350 G4double dmT=dmNeut; // Doubled mass 00351 //G4int Z=0; // Prototype of the isotope Z 00352 //G4int N=1; // Prototype of the Isotope N 00353 if(rPDG==2212) // Correct it, if this is a proton 00354 { 00355 mT=mProt; // Prototype of mass of QF nucleon to be excited 00356 mT2=mProt2; // Squared mass of the free nucleon 00357 dmT=dmProt; // Doubled mass 00358 //Z=1; // Z of the isotope 00359 //N=0; // N of the Isotope 00360 } 00361 G4double mP2=pr4M.m2(); // Squared mass of the projectile 00362 if(mP2<0.) mP2=0.; // Can be a problem for photon (m_min = 2*m_pi0) 00363 G4double s_value=tot4M.m2(); // @@ Check <0 ... 00364 G4double E=(s_value-mT2-mP2)/dmT; // Effective interactionEnergy (virtNucl target) 00365 G4double E2=E*E; 00366 if(E<0. || E2<mP2) // Impossible to fragment: return projectile 00367 { 00368 #ifdef pdebug 00369 G4cerr<<"-Warning-G4DifR::TFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; 00370 #endif 00371 return ResHV; // *** Do Nothing Action *** 00372 } 00373 G4double mP=std::sqrt(mP2); // Calculate mass of the projectile (to be exc.) 00374 if(mP<.1) mP=mPi0; // For photons minDiffraction is gam+P->P+Pi0 00375 //G4double dmP=mP+mP; // Doubled mass of the projectile 00376 G4double mMin=mP+mPi0; // Minimum diffractive mass 00377 G4double tA=tgA; // Real A of the target 00378 G4double sA=5./std::pow(tA,third); // Mass-screaning 00379 //mMin+=mPi0+G4UniformRand()*(mP*sA+mPi0); // *Experimental* 00380 mMin+=G4UniformRand()*(mP*sA+mPi0); // *Experimental* 00381 G4double ss=std::sqrt(s_value); // CM compound mass (sqrt(s)) 00382 G4double mMax=ss-mP; // Maximum diffraction mass of the projectile 00383 if(mMax>maxDM) mMax=maxDM; // Restriction to avoid too big masses 00384 if(mMin>=mMax) 00385 { 00386 #ifdef pdebug 00387 G4cerr<<"-Warning-G4DifR::TFra:ZeroDiffractionMRange, mi="<<mMin<<",ma="<<mMax<<G4endl; 00388 #endif 00389 return ResHV; // Do Nothing Action 00390 } 00391 G4double R = G4UniformRand(); 00392 G4double mDif=std::exp(R*std::log(mMax)+(1.-R)*std::log(mMin)); // Low Mass Approximation 00393 G4double mDif2=mDif*mDif; 00394 G4double ds=s_value-mP2-mDif2; 00395 //G4double e=ds/dmP; 00396 //G4double P=std::sqrt(e*e-mDif2); // Momentum in pseudo laboratory system 00397 #ifdef debug 00398 G4cout<<"G4QDiffR::TargFrag:Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl; 00399 #endif 00400 // @@ Temporary NN t-dependence for all hadrons 00401 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl; 00402 G4double maxt=(ds*ds-4*mP2*mDif2)/s_value; // maximum possible -t 00403 G4double tsl=140000.; // slope in MeV^2 00404 G4double t=-std::log(G4UniformRand())*tsl; 00405 #ifdef pdebug 00406 G4cout<<"G4QDifR::TFra:ph="<<pPDG<<",P="<<P<<",t="<<t<<"<"<<maxt<<G4endl; 00407 #endif 00408 #ifdef nandebug 00409 if(mint>-.0000001); // To make the Warning for NAN 00410 else G4cout<<"******G4QDiffractionRatio::TargFragment: -t="<<mint<<G4endl; 00411 #endif 00412 G4double rt=t/maxt; 00413 G4double cost=1.-rt-rt; // cos(theta) in CMS 00414 #ifdef ppdebug 00415 G4cout<<"G4QDiffraRatio::TargFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl; 00416 #endif 00417 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) 00418 { 00419 if (cost>1.) cost=1.; 00420 else if(cost<-1.) cost=-1.; 00421 else 00422 { 00423 G4cerr<<"G4QDiffRat::TargFragm: *NAN* cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl; 00424 return ResHV; // Do Nothing Action 00425 } 00426 } 00427 G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mP); // 4mom of the leading nucleon 00428 G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif); // 4mom of the diffract. Quasmon 00429 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); 00430 if(!G4QHadron(tot4M).RelDecayIn2(r4M, d4M, dir4M, cost, cost)) 00431 { 00432 G4cerr<<"G4QDifR::TFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl; 00433 //G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp"); 00434 return ResHV; // Do Nothing Action 00435 } 00436 #ifdef debug 00437 G4cout<<"G4QDifRat::TargFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl; 00438 #endif 00439 // Now everything is ready for fragmentation and DoNothing projHadron must be wiped out 00440 delete hadron; // Delete the fake (doNothing) projectile hadron 00441 ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile 00442 hadron = new G4QHadron(pPDG,r4M); // Hadron for the recoil nucleon 00443 ResHV->push_back(hadron); // Fill the recoil nucleon 00444 #ifdef debug 00445 G4cout<<"G4QDiffractionRatio::TargFragm: *Filled* LeadingNuc="<<r4M<<pPDG<<G4endl; 00446 #endif 00447 G4QHadronVector* leadhs = 0; // Prototype of Quasmon Output G4QHadronVector ---->---* 00448 G4QContent dQC=G4QPDGCode(rPDG).GetQuarkContent(); // QuarkContent of quasiFreeNucleon | 00449 G4Quasmon* quasm = new G4Quasmon(dQC,d4M); // Quasmon=DiffractionExcitationQuasmon-* | 00450 #ifdef debug 00451 G4cout<<"G4QDiffRatio::TgFrag:tPDG="<<tPDG<<",rPDG="<<rPDG<<",d4M="<<d4M<<G4endl;//| | 00452 #endif 00453 G4QEnvironment* pan= new G4QEnvironment(G4QNucleus(tPDG));// --> DELETED --->---* | | 00454 pan->AddQuasmon(quasm); // Add diffractiveQuasmon to Environ.| | | 00455 #ifdef debug 00456 G4cout<<"G4QDiffractionRatio::TargFragment: EnvPDG="<<tPDG<<G4endl; // | | | 00457 #endif 00458 try // | | | 00459 { // | | | 00460 leadhs = pan->Fragment();// DESTROYED in the end of the LOOP work space | | <-| 00461 } // | | | 00462 catch (G4QException& error)// | | | 00463 { // | | | 00464 //#ifdef pdebug 00465 G4cerr<<"***G4QDiffractionRatio::TargFrag: G4QException is catched"<<G4endl;//| | | 00466 //#endif 00467 // G4Exception("G4QDiffractionRatio::TargFragm:","27",FatalException,"*Nucl");// | | | 00468 G4Exception("G4QDiffractionRatio::TargFragment()","HAD_CHPS_0027", 00469 FatalException, "Nucl"); 00470 } // | | | 00471 delete pan; // Delete the Nuclear Environment <-<--*--* | 00472 G4int qNH=leadhs->size(); // A#of collected hadrons from diff.frag. | 00473 if(qNH) for(G4int iq=0; iq<qNH; iq++) // Loop over hadrons to fill the result | 00474 { // | 00475 G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron | 00476 ResHV->push_back(loh); // Fill in the result | 00477 } // | 00478 leadhs->clear();// | 00479 delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<---* 00480 return ResHV; // Result 00481 } // End of TargFragment
1.4.7