// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // $Id: G4QNucleus.cc,v 1.120 2010/06/24 16:20:11 mkossov Exp $ // GEANT4 tag $Name: hadr-chips-V09-03-08 $ // // ---------------- G4QNucleus ---------------- // by Mikhail Kossov, Sept 1999. // class for Nuclei/Nuclear Environment used by CHIPS Model // --------------------------------------------------------------------- // Short description: a class describing properties of nuclei, which // are necessary for the CHIPS Model. // --------------------------------------------------------------------- //#define debug //#define pdebug //#define cldebug //#define qdebug //#define cldebug //#define pardeb //#define ppdebug #include "G4QNucleus.hh" #include "Randomize.hh" #include #include #include //#include using namespace std; // Static parameters definition G4double G4QNucleus::freeNuc=0.1; // probability to find quasiFreeBaryon on Surface G4double G4QNucleus::freeDib=.05; // probability to find quasiFreeDiBaryon on Surface G4double G4QNucleus::clustProb=4.; // clusterization probability in dense region G4double G4QNucleus::mediRatio=1.; // relative vacuum hadronization probability G4double G4QNucleus::nucleonDistance=.8*fermi; // Distance between nucleons (0.8 fm) (Body) G4double G4QNucleus::WoodSaxonSurf=.545*fermi; // WoodSaxon Surface Param (0.545 fm) (Body) G4QNucleus::G4QNucleus(): G4QHadron(), Z(0), N(0), S(0), dZ(0), dN(0), dS(0), maxClust(0), theNucleons(),currentNucleon(-1), rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false) { probVect[0]=mediRatio; for(G4int i=1; i<256; i++) {probVect[i] = 0.;} #ifdef pardeb G4cout<<"G4QNucleus::Constructor:(1) N="<Called< PDG="<80000000 && nucPDG<100000000) // Try to convert the NUCCoding to PDGCoding { G4QPDGCode(22).ConvertPDGToZNS(nucPDG, z, n, s); Z =z; N =n; S =s; #ifdef debug G4cout<<"G4QNucleus::InitByPDG:Z="<QPDG="<2 maxi=3; #ifdef debug G4cout<<"G4QNucleus::UpdateClusters:p1="<Get4Momentum(); // 4-mom of the current nucleon G4double srP2=(t4M-n4M).vect().mag2(); // p2 of the subResNucleus G4double m2=m2n; // default subResNucleusM2 (for neutrons) if((*u)->GetPDGCode()==2212) m2=m2p; // change it to subResNucleusM2 for protons G4double srE=std::sqrt(srP2+m2); // Energy of the subResNucleus #ifdef debug G4cout<<"G4QNucleus::SubtractNucleon:#"< 2N="<SetQPDG(del); h2->SetQPDG(nuc); h1->Set4Momentum(n14M); h2->Set4Momentum(n24M); return true; } else { G4cerr<<"***G4QNucleus::EvaporateBaryon: M="< mPi+mPi+nucM+nucM) { G4LorentzVector n14M(0.,0.,0.,nucM); G4LorentzVector n24M(0.,0.,0.,nucM); G4LorentzVector pi4M(0.,0.,0.,mPi+mPi); if(!DecayIn3(n14M, n24M, pi4M)) { G4cerr<<"***G4QNucl::EvapBary: (anti) tM="< 2N="<SetQPDG(del); h2->SetQPDG(del); h1->Set4Momentum(n14M); h2->Set4Momentum(n24M); return true; } else { G4cerr<<"***G4QNucleus::EvaporateBaryon: M="<SetQPDG(apQPDG); h2->SetQPDG(apQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else if(N==-2) { h1mom=G4LorentzVector(0.,0.,0.,mNeut); h2mom=h1mom; h1->SetQPDG(anQPDG); h2->SetQPDG(anQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else if(N==-1 && Z==-1) { h1mom=G4LorentzVector(0.,0.,0.,mProt); h2mom=G4LorentzVector(0.,0.,0.,mNeut); h1->SetQPDG(apQPDG); h2->SetQPDG(anQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else if(Z==-1 && S==-1) { h1mom=G4LorentzVector(0.,0.,0.,mProt); h2mom=G4LorentzVector(0.,0.,0.,mLamb); h1->SetQPDG(apQPDG); h2->SetQPDG(alQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else { h1mom=G4LorentzVector(0.,0.,0.,mNeut); h2mom=G4LorentzVector(0.,0.,0.,mLamb); h1->SetQPDG(anQPDG); h2->SetQPDG(alQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } h1->Set4Momentum(h1mom); h2->Set4Momentum(h2mom); return true; } else if(a==2) { if(Z<0||N<0) { G4int nucPDG = 2112; G4double nucM = mNeut; G4int piPDG = -211; G4QPDGCode db = NNQPDG; G4QPDGCode pi = PIMQPDG; if(N<0) { nucPDG = 2212; nucM = mProt; piPDG = 211; db = PPQPDG; pi = PIPQPDG; } if(totMass>mPi+nucM+nucM) { G4LorentzVector n14M(0.,0.,0.,nucM); G4LorentzVector n24M(0.,0.,0.,nucM); G4LorentzVector pi4M(0.,0.,0.,mPi); if(!DecayIn3(n14M,n24M,pi4M)) { G4cerr<<"***G4QNucl::EvapBary: tM="< 2N="<SetQPDG(db); h2->SetQPDG(pi); h1->Set4Momentum(n14M); h2->Set4Momentum(pi4M); return true; } else { G4cerr<<"***G4QNucleus::EvaporateBaryon: M="<SetQPDG(pQPDG); h2->SetQPDG(pQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else if(N==2) { h1mom=G4LorentzVector(0.,0.,0.,mNeut); h2mom=h1mom; h1->SetQPDG(nQPDG); h2->SetQPDG(nQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else if(N==1&&Z==1) { if(totMass<=mNP) { #ifdef debug G4cout<<"G4QNucl::EvaporateBaryon: Photon ### d+g ###, dM="<SetQPDG(gQPDG); h2->SetQPDG(NPQPDG); } else { h1mom=G4LorentzVector(0.,0.,0.,mProt); h2mom=G4LorentzVector(0.,0.,0.,mNeut); h1->SetQPDG(pQPDG); h2->SetQPDG(nQPDG); } if(!DecayIn2(h1mom,h2mom)) return false; } else if(Z==1&&S==1) { h1mom=G4LorentzVector(0.,0.,0.,mProt); h2mom=G4LorentzVector(0.,0.,0.,mLamb); h1->SetQPDG(pQPDG); h2->SetQPDG(lQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } else { h1mom=G4LorentzVector(0.,0.,0.,mNeut); h2mom=G4LorentzVector(0.,0.,0.,mLamb); h1->SetQPDG(nQPDG); h2->SetQPDG(lQPDG); if(!DecayIn2(h1mom,h2mom)) return false; } h1->Set4Momentum(h1mom); h2->Set4Momentum(h2mom); return true; } else if(a>2) { G4bool nFlag = false; // Flag of possibility to radiate neutron G4bool pFlag = false; // Flag of possibility to radiate proton G4bool lFlag = false; // Flag of possibility to radiate lambda G4bool aFlag = false; // Flag of possibility to radiate alpha G4bool nnFlag = false; // Flag of possibility to radiate 2 neutrons G4bool npFlag = false; // Flag of possibility to radiate neutron+proton G4bool nlFlag = false; // Flag of possibility to radiate neutron+lambda G4bool ppFlag = false; // Flag of possibility to radiate 2 protons G4bool plFlag = false; // Flag of possibility to radiate proton+lambda G4bool llFlag = false; // Flag of possibility to radiate 2 lambdas G4bool paFlag = false; // Flag of possibility to radiate proton+alpha G4bool naFlag = false; // Flag of possibility to radiate neutron+alpha G4bool laFlag = false; // Flag of possibility to radiate lambda+alpha G4bool aaFlag = false; // Flag of possibility to radiate alpha+alpha G4bool nnnF = false; G4bool nnpF = false; G4bool nppF = false; G4bool pppF = false; G4bool nnlF = false; G4bool nplF = false; G4bool pplF = false; G4bool nllF = false; G4bool pllF = false; G4bool lllF = false; G4bool nnaF = false; G4bool npaF = false; G4bool ppaF = false; G4bool nlaF = false; G4bool plaF = false; G4bool llaF = false; G4bool paaF = false; G4bool naaF = false; G4bool laaF = false; G4bool aaaF = false; G4double GSMass = GetGSMass(); // Ground State mass of the Nucleus G4double GSResNN= GSMass; // Prototype of Residual Nuclear Mass for n+n G4double GSResNP= GSMass; // Prototype of Residual Nuclear Mass for n+p G4double GSResNL= GSMass; // Prototype of Residual Nuclear Mass for n+l G4double GSResPP= GSMass; // Prototype of Residual Nuclear Mass for p+p G4double GSResPL= GSMass; // Prototype of Residual Nuclear Mass for p+l G4double GSResLL= GSMass; // Prototype of Residual Nuclear Mass for l+l G4double GSResNA= GSMass; // Prototype of Residual Nuclear Mass for n+alpha G4double GSResPA= GSMass; // Prototype of Residual Nuclear Mass for p+alpha G4double GSResLA= GSMass; // Prototype of Residual Nuclear Mass for l+alpha G4double GSResAA= GSMass; // Prototype of Residual Nuclear Mass for alpha+alpha G4double GSResNa= GSMass; // Prototype of Residual Nuclear Mass for alpha G4double GSReNNN= GSMass; // Prototype of Residual Nuclear Mass for n+n+n G4double GSReNNP= GSMass; // Prototype of Residual Nuclear Mass for n+n+p G4double GSReNPP= GSMass; // Prototype of Residual Nuclear Mass for n+p+p G4double GSRePPP= GSMass; // Prototype of Residual Nuclear Mass for p+p+p G4double GSReNNL= GSMass; // Prototype of Residual Nuclear Mass for n+n+l G4double GSReNPL= GSMass; // Prototype of Residual Nuclear Mass for n+p+l G4double GSRePPL= GSMass; // Prototype of Residual Nuclear Mass for p+p+l G4double GSReNLL= GSMass; // Prototype of Residual Nuclear Mass for n+l+l G4double GSRePLL= GSMass; // Prototype of Residual Nuclear Mass for p+l+l G4double GSReLLL= GSMass; // Prototype of Residual Nuclear Mass for l+l+l G4double GSReNNA= GSMass; // Prototype of Residual Nuclear Mass for n+n+a G4double GSReNPA= GSMass; // Prototype of Residual Nuclear Mass for n+p+a G4double GSRePPA= GSMass; // Prototype of Residual Nuclear Mass for p+p+a G4double GSReNLA= GSMass; // Prototype of Residual Nuclear Mass for n+l+a G4double GSRePLA= GSMass; // Prototype of Residual Nuclear Mass for p+l+a G4double GSReLLA= GSMass; // Prototype of Residual Nuclear Mass for l+l+a G4double GSRePAA= GSMass; // Prototype of Residual Nuclear Mass for p+a+a G4double GSReNAA= GSMass; // Prototype of Residual Nuclear Mass for n+a+a G4double GSReLAA= GSMass; // Prototype of Residual Nuclear Mass for l+a+a G4double GSReAAA= GSMass; // Prototype of Residual Nuclear Mass for a+a+a G4QPDGCode PQPDG(22); // Prototype of QPDG for ResidualNucleus to proton G4QPDGCode NQPDG(22); // Prototype of QPDG for ResidualNucleus to neutron G4QPDGCode LQPDG(22); // Prototype of QPDG for ResidualNucleus to lambda G4QPDGCode AQPDG(22); // Prototype of QPDG for ResidualNucleus to alpha G4QPDGCode nnQPDG(22); // Prototype of QPDG for ResidualNucleus to nn-dibar. G4QPDGCode npQPDG(22); // Prototype of QPDG for ResidualNucleus to np-dibar. G4QPDGCode nlQPDG(22); // Prototype of QPDG for ResidualNucleus to nl-dibar. G4QPDGCode ppQPDG(22); // Prototype of QPDG for ResidualNucleus to pp-dibar. G4QPDGCode plQPDG(22); // Prototype of QPDG for ResidualNucleus to pl-dibar. G4QPDGCode llQPDG(22); // Prototype of QPDG for ResidualNucleus to ll-dibar. G4QPDGCode naQPDG(22); // Prototype of QPDG for ResidualNucleus to n+alpha G4QPDGCode paQPDG(22); // Prototype of QPDG for ResidualNucleus to p+alpha G4QPDGCode laQPDG(22); // Prototype of QPDG for ResidualNucleus to l+alpha G4QPDGCode aaQPDG(22); // Prototype of QPDG for ResidualNucleus to alph+alph G4QPDGCode dbQPDG(22); // Prototype of chosen dibaryon QPDG G4QPDGCode fQPDG(22); // Prototype of QPDG of the Second Baryon G4double rMass = 0.; // Prototype of mass of Residual Nucleus G4double eMass = 0.; // Prototype of mass of Evaporated Baryon G4double fMass = 0.; // Prototype of mass of the Second Baryon #ifdef debug G4cout<<"G4QNuc::EvaB:a>2, totM="< GSMass="<0) { PQPDG=G4QPDGCode(90000000+1000*(1000*S+Z-1)+N); GSResNp=PQPDG.GetMass(); G4double mpls=GSResNp+mProt; G4double mmin=GSResNp-mProt; pp2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2; if(pp2m>=0.000001) { pFlag=true; pBnd=mProt-GSMass+GSResNp; // Binding energy for proton G4double eMax=sqrt(mP2+pp2m); #ifdef debug G4cout<<"G4QNuc::EvapBaryon:pm="<"<"<"<"< pMin) || (nFlag && nExcess > nMin) || (lFlag && lExcess > lMin) || (aFlag && aExcess > aMin) ) && minE"<nCBPP&&r<=lCBPP) { #ifdef debug G4cout<<"G4QNuc::EvaB:LAMBDA is selected for evap,r="<zCBPP&&r<=nCBPP) { #ifdef debug G4cout<<"G4QNuc::EvaBar: N is selected for evapor,r="<2)three=false; // @@@@@@@@@@@@@@@@@@ //else if(PDG==pPDG&&(pnCond&&ppCond&&plCond&&paCond)) // @@@@@@@@@@@@@@@@@@@ if(PDG==pPDG&&(pnCond&&ppCond&&plCond&&paCond))//p+RN decay, p+b+RN dec is closed { #ifdef debug G4cout<<"G4QN::EB:*p*: n="<1&&GSResPP!=GSMass&&fMass+mProt+SPPBarr+GSResPP2&&N>1&&a>4&&GSResPL!=GSMass&&fMass+SAPBarr+mAlph+GSResPAzLim&&r<=sLim) { eMass = mLamb; dbQPDG= PLQPDG; rMass = GSResPL; rQPDG = plQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: P+L"<nLim&&r<=zLim) { eMass = mProt; dbQPDG= PPQPDG; rMass = GSResPP; rQPDG = ppQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: P+P"<zLim&&r<=sLim) { eMass = mLamb; dbQPDG= NLQPDG; rMass = GSResNL; rQPDG = nlQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: N+L"<nLim&&r<=zLim) { eMass = mProt; dbQPDG= NPQPDG; rMass = GSResNP; rQPDG = npQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: N+P"<zLim&&r<=sLim) { eMass = mLamb; dbQPDG= LLQPDG; rMass = GSResLL; rQPDG = llQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: L+L"<nLim&&r<=zLim) { eMass = mProt; dbQPDG= PLQPDG; rMass = GSResPL; rQPDG = plQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: L+P"<zLim&&r<=sLim) { eMass = mLamb; dbQPDG= LAQPDG; rMass = GSResLA; rQPDG = laQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: A+L"<nLim&&r<=zLim) { eMass = mProt; dbQPDG= PAQPDG; rMass = GSResPA; rQPDG = paQPDG; #ifdef debug G4cout<<"G4QNucleus::EvaporateBary: A+P"< Decay in 3 Baryons + Residual is impossible at this point { if(secB) // Decay in 2Baryons(2a,a+bary)+ResidN is possible //if(2>3) { #ifdef debug G4cout<<"G4QNucleus::EvaporateBaryon: Decay in 2 baryons"<mNeut+mNeut+GSResNN) nnLim+=N*(N-1)*pow(totMass-mNeut-mNeut-GSResNN,2); G4double nzLim=nnLim; if(npFlag&&totMass>mNeut+mProt+PBarr+GSResNP) { if(barf) nzLim+=2*N*Z*pow(totMass-mNeut-mProt-PBarr-GSResNP,2); else nzLim+=2*N*Z*pow(totMass-mNeut-mProt-GSResNP,2); } G4double zzLim=nzLim; if(ppFlag&&totMass>mProt+mProt+SPPBarr+GSResPP) { if(barf) zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-SPPBarr-GSResPP,2); else zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-GSResPP,2); } G4double nlLim=zzLim; if(nlFlag&&totMass>mNeut+mLamb+GSResNL) nlLim+=2*N*S*pow(totMass-mNeut-mLamb-GSResNL,2); G4double zlLim=nlLim; if(plFlag&&totMass>mProt+PBarr+mLamb+GSResPL) { if(barf) zlLim+=2*Z*S*pow(totMass-mProt-mLamb-PBarr-GSResPL,2); else zlLim+=2*Z*S*pow(totMass-mProt-mLamb-GSResPL,2); } G4double llLim=zlLim; if(llFlag&&totMass>mLamb+mLamb+GSResLL) llLim+=S*(S-1)*pow(totMass-mLamb-mLamb-GSResLL,2); G4double naLim=llLim; if(naFlag&&totMass>mNeut+mAlph+ABarr+GSResNA) { if(barf) naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-ABarr-GSResNA,2); else naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-GSResNA,2); } G4double zaLim=naLim; if(paFlag&&totMass>mProt+SAPBarr+mAlph+GSResPA) { if(barf) zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2); else zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-GSResPA,2); } G4double laLim=zaLim; if(laFlag&&totMass>mLamb+mAlph+ABarr+GSResLA) { if(barf) laLim+=S*alp*pow(totMass-mLamb-mAlph-ABarr-GSResLA,2); else laLim+=S*alp*pow(totMass-mLamb-mAlph-GSResLA,2); } G4double aaLim=laLim; if(evalph&&aaFlag&&totMass>mAlph+mAlph+SAABarr+GSResAA) { if(barf) aaLim+=alp*pow(totMass-mAlph-mAlph-SAABarr-GSResAA,2)* evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7); else aaLim+=alp*pow(totMass-mAlph-mAlph-GSResAA,2)* evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7); } G4double r = aaLim*G4UniformRand(); if (aaLim&&laLimSetQPDG(gQPDG); h1->Set4Momentum(G4LorentzVector(0.,0.,0.,0.)); if (Z>0) { h2->SetQPDG(pQPDG); h2->Set4Momentum(G4LorentzVector(0.,0.,0.,mProt)); } else if(N>0) { h2->SetQPDG(nQPDG); h2->Set4Momentum(G4LorentzVector(0.,0.,0.,mNeut)); } else if(S>0) { h2->SetQPDG(lQPDG); h2->Set4Momentum(G4LorentzVector(0.,0.,0.,mLamb)); } else return false; return true; } if(a<-2) G4cerr<<"***G4QNucleus::EvaporateBaryon: A="<s && i candidates to hadrons static const G4int nOfBaryons=72; //a#of 1/2,3/2,5/2,7/2 Baryons => candidates to hadrons // Scalar resonances (0): Eta,Pi0,Pi+,APi-,Ka0,Ka+,AKa0,AKa-,Eta* static G4int mesonPDG[nOfMesons] = {221,111,211,-211,311,321,-311,-321,331, // 0- 8 // Vector resonances (1): omega,Rh0,Rh+,Rho-,K0*,K+*,AK0*,AK-*,Phi 223,113,213,-213,313,323,-313,-323,333, // 9-18 // Tensor D-resonances (2): f2 ,a20,a2+, a2-,K20,K2+,AK20,AK2-,f2' 225,115,215,-215,315,325,-315,-325,335, // 19-27 // Tensor F-resonances (3): om3,ro3,r3+,rh3-,K30,K3+,AK30,AK3-,Phi3 227,117,217,-217,317,327,-317,-327,337, // 28-35 // Tensor G-resonances (4): f4 ,a40,a4+, a4-,K40,K4+,AK40,AK4-,f4' 229,119,219,-219,319,329,-319,-329,339}; // 36-44 // Baryon octet (1/2): n , an , p , ap ,lamb,alamb, sig-,asig- static G4int baryonPDG[nOfBaryons]={2112,-2112,2212,-2212,3122,-3122,3112,-3112, // 45-52 // Hyperon octet (1/2): sig0,asig0,sig+,asig+,ksi-,aksi-,ksi0,aksi0 3212,-3212,3222,-3222,3312,-3312,3322,-3322, // 53-60 // Baryon decuplet (3/2): del-,adel-,del0,adel0,del+,adel+,dl++,adl++,sis-,asis- 1114,-1114,2114,-2114,2214,-2214,2224,-2224,3114,-3114,//70 // sis0,asis0,sis+,asis+,kss-,akss-,kss0,akss0,omeg,aomeg 3214,-3214,3224,-3224,3314,-3314,3324,-3324,3334,-3334,//80 // Baryon octet (5/2): n5/2,an5/2,p5/2,ap5/2,l5/2,al5/2,si5-,asi5- 2116,-2116,2216,-2216,3126,-3126,3116,-3116, // 81-88 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50 3216,-3216,3226,-3226,3316,-3316,3326,-3326, // 89-96 // Baryon decuplet (7/2): dl5-,adl5-,dl50,adl50,dl5+,adl5+,d5++,ad5++,si5-,asi5- 1118,-1118,2118,-2118,2218,-2218,2228,-2228,3118,-3118, //106 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50,ome5,aome5 3218,-3218,3228,-3228,3318,-3318,3328,-3328,3338,-3338};//116 G4int i=0; #ifdef debug G4int ind=0; #endif G4int iQC = theQCandidates.size(); if(iQC) for(G4int jq=0; jqnOfMesons) maxMes=nOfMesons; if(maxMes>=0) for (i=0; inOfBaryons) maxBar=nOfBaryons; if(maxBar>=0) for (i=0; i=0) for (i=0; i2) G4double comb=ae0*(ae0-1)/2; // Product up to ac=2 if(comb<=0.) comb=1.; #ifdef cldebug G4double sZ=0.; // Percent of protons G4double sN=0.; // Percent of neutrons G4cout<<"G4QN::PC:C#"<GetPDGCode(); G4int cBN = curCand->GetBaryonNumber(); G4int cST = curCand->GetStrangeness(); // *********************************************************************************** // These are first 117 candidates which are defined in G4QNucleus::InitCandidateVector // ***!!!*** if they are changed there the corresponding change must be done here //static const G4int nOfMesons =45;//a#of S=0,1,2,3,4 Mesons, => candidates to hadrons //static const G4int nOfBaryons=72;//a#of 1/2,3/2,5/2,7/2 Baryons => cand's to hadrons // Scalar resonances (0): Eta,Pi0,Pi+,APi-,Ka0,Ka+,AKa0,AKa-,Eta* //static G4int mesonPDG[45] = {221,111,211,-211,311,321,-311,-321,331, // 0- 8 // Vector resonances (1): omega,Rh0,Rh+,Rho-,K0*,K+*,AK0*,AK-*,Phi // 223,113,213,-213,313,323,-313,-323,333, // 9-18 // Tensor D-resonances (2): f2 ,a20,a2+, a2-,K20,K2+,AK20,AK2-,f2' // 225,115,215,-215,315,325,-315,-325,335, // 19-27 // Tensor F-resonances (3): om3,ro3,r3+,rh3-,K30,K3+,AK30,AK3-,Phi3 // 227,117,217,-217,317,327,-317,-327,337, // 28-35 // Tensor G-resonances (4): f4 ,a40,a4+, a4-,K40,K4+,AK40,AK4-,f4' // 229,119,219,-219,319,329,-319,-329,339}; // 36-44 // Baryon octet (1/2): n , an , p , ap ,lamb,alamb, sig-,asig- //static G4int baryonPDG[72]={2112,-2112,2212,-2212,3122,-3122,3112,-3112, // 45-52 // Hyperon octet (1/2): sig0,asig0,sig+,asig+,ksi-,aksi-,ksi0,aksi0 // 3212,-3212,3222,-3222,3312,-3312,3322,-3322, // 53-60 // Baryon decuplet (3/2): del-,adel-,del0,adel0,del+,adel+,dl++,adl++,sis-,asis- // 1114,-1114,2114,-2114,2214,-2214,2224,-2224,3114,-3114,//70 // sis0,asis0,sis+,asis+,kss-,akss-,kss0,akss0,omeg,aomeg // 3214,-3214,3224,-3224,3314,-3314,3324,-3324,3334,-3334,//80 // Baryon octet (5/2): n5/2,an5/2,p5/2,ap5/2,l5/2,al5/2,si5-,asi5- // 2116,-2116,2216,-2216,3126,-3126,3116,-3116, // 81-88 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50 // 3216,-3216,3226,-3226,3316,-3316,3326,-3326, // 89-96 // Baryon decuplet (7/2): dl5-,adl5-,dl50,adl50,dl5+,adl5+,d5++,ad5++,si5-,asi5- // 1118,-1118,2118,-2118,2218,-2218,2228,-2228,3118,-3118, //106 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50,ome5,aome5 // 3218,-3218,3228,-3228,3318,-3318,3328,-3328,3338,-3338};//116 // One should take into account, that #of mesons & baryons can be cut in G4Quas::HadrQE //G4int nP= theWorld->GetQPEntries(); // A#of initialized particles in CHIPS World ////@@ Make special parametyer to cut high resonances for nuclear fragmentation !! //G4int nMesons = 45; //if (nP<34) nMesons = 9; //else if(nP<51) nMesons = 18; //else if(nP<65) nMesons = 27; //else if(nP<82) nMesons = 36; //G4int nBaryons = 72; //if (nP<45) nBaryons = 16; //else if(nP<59) nBaryons = 36; //else if(nP<76) nBaryons = 52; // ********************************************************************************** //G4int cS = curCand->GetStrangeness(); //if(piF&&gaF&&cPDG!=90000001&&cPDG!=90001000) // Both flags, in case of pi-first-int //if(piF&&gaF&&cBN!=1&&cBN!=3) // Both flags, which is in case of pi-first-int if(piF&&gaF&&cBN!=1)// @@ Should be both, which is in case of pi-first-interaction @@ ? //if(piF&&gaF&&cBN!=1&&cBN!=4) // Should be both, in case of pi-first-interaction { curCand->SetPreProbability(0.); curCand->SetDenseProbability(0.); curCand->SetPossibility(false); #ifdef cldebug if(cPDG==90001001) G4cout<<"G4QNuc::PrepCand: piF/gaF fragments are blocked"<4||cST>2))// @@ PreClosed HighSpin/HighStrange { curCand->SetPreProbability(0.); curCand->SetDenseProbability(0.); curCand->SetPossibility(false); } else { G4double tnM=GetQPDG().GetMass(); // Total mass of this nucleus if(cPDG>80000000&&cPDG!=90000000) // ===> Cluster case { G4int sc = cST; // "S" of the cluster G4int zc = curCand->GetCharge(); // "Z" of the cluster G4int ac = cBN; // "A" of the cluster G4int nc = ac-zc-sc; // "N" of the cluster G4double cM=tnM-G4QNucleus(Z-zc,N-nc,S-sc).GetGSMass(); // BoundMass of the cluster G4LorentzVector intLV=pLV+G4LorentzVector(0.,0.,0.,cM); // 4-mom of the proj+clust pos = probVect[ac]; // Cluster Probability NormalizationFactor if(ac<=maxClust&&pos>0.&&(pLV==zeroLV||intLV.m()>.00001+cM)) { #ifdef cldebug G4cout<<"G4QNucleus::PrepareCand: ac="<2 { if(acm0) comb*=(ae1-ac)/ac; acm=ac; s=0.; cca=0; if(ac%2) mac=7; // @@Change it if cluster set is changed else mac=8; // @@ It is not yet automatic #ifdef cldebug G4cout<<"G4QNuc::PrepCl:c="<0) for(int iz=1; iz<=zc; iz++) prod*=(ze1-iz)/iz; if(nc>0) for(int in=1; in<=nc; in++) prod*=(ne1-in)/in; if(sc>0) for(int is=1; is<=sc; is++) prod*=(se1-is)/is; } s+=prod; pos*=prod; pos/=comb; #ifdef cldebug if(pos) G4cout<<"G4QN::PreC:c="<SetPreProbability(0.); curCand->SetDenseProbability(0.); curCand->SetPossibility(false); // This candidate is not possible } } else { #ifdef cldebug G4cout<<"G4QNucl::PrepCand:cPDG="<SetPreProbability(pos); // ===> Hadronic case in Vacuum curCand->SetDenseProbability(0.); // ===> Hadronic case in Vacuum } curCand->SetPossibility(true); // All candidates are possible at this point } } // End of the LOOP over Candidates #ifdef cldebug G4cout<<"G4QNucl::PrepCand:covP="<(cZ); G4int cN=static_cast(cA-cZ); G4int sZ=iZ; G4int sN=cN; if(delZ||dA) { G4int dz=static_cast(delZ); G4int dn=static_cast(dA-delZ); GSM=G4QNucleus(Z-dz,N-dn,S).GetGSMass(); sZ-=dz; sN-=dn; } return G4QNucleus(Z-sZ,N-sN,S).GetGSMass()+G4QNucleus(iZ,cN,0).GetGSMass()-GSM; } // End of "BindingEnergy" //Coulomb Barrier Reflection Probability (CB - Coulomb Barrier, E - Kinetic Energy) G4double G4QNucleus::CoulBarPenProb(const G4double& CB, const G4double& E, const G4int& C, const G4int& B) {// = A.Lepretre et al, Nucl.Phys., A390 (1982) 221-239 = static const G4double mNeut= G4QPDGCode(2112).GetMass(); // Mass of neutron static const G4double dNeut= mNeut+mNeut; // DiMass of neutron static const G4double mProt= G4QPDGCode(2212).GetMass(); // Mass of proton static const G4double dProt= mProt+mProt; // DiMass of proton static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0); // Mass of deuteron //static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0); // Mass of tritium //static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0); // Mass of Helium 3 //static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); // Mass of alpha static const G4double wellDebth=40.; //@@ Should be jus binding energy @@ done // @@ --- Temporary 1 ---> close the OverBarrierReflection for all //return 1.; // ^^^^^^^---> End of Themporary 1 // @@ --- Temporary 2 ---> close the OverBarrierReflection for fragments and mesons //if(E2) return 1.; if(C>B+1) { #ifdef debug G4cout<<"-Warning-G4QN::CBPP:SubtractedChrg="<SubtractedBaryonNmbr="< End of Themporary 2 //G4double nA=GetA(); //G4double nA=GetA()-B; //if(nA==40) G4cout<<"G4QN::CBPP:Z="<8&&nA<12||nA>16&&nA<40) return 1.;// "OverBarrierReflection is closed" //else if(nA>8&&nA<12||nA>16&&nA<40) return 1.; // "OverBarrierReflection is closed" Cond //else if(nA<12||nA>16&&nA<40) return 1.; // "OverBarrierReflection is closed" Condition //else if(nA<12||nA>16) return 1.; // "OverBarrierReflection is closed" Condition //else if(nA<12) return 1.; // @@@@@ Over barrier reflection is closed @@@ !!! @@@ //if(B+B>Z+N+S) return 1.; //G4double wD=wellDebth*B; G4double wD=wellDebth; //G4double wD=0.; // @@ --- Temporary 3 ---> close the OverBarrierReflection for mesons //if(!B) wD=0.; // ^^^^^^^---> End of Themporary 3 G4double GSM=GetGSMass(); #ifdef debug G4cout<<"G4QNucl::CBPenProb:GSM="<3); // @@ Temporary "Mass Barrier for mesons" @@ __________________ //else if(!B) wD=40.; // @@ End of Temporary^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ////else if(nA<7&&B>0) wD=0.; // Only Coulomb Barrier can reflect !!! ////else if((nA<12||nA>16)&&B>0) wD=0.;// Only CoulombB can reflect !!! O16 E-dep of gamA ////else if((nA<12||nA>27)&&B>0) wD=0.;// Only CoulombB can reflect !!! O16 E-dep of gamA ////else if(nA<9&&B>0) return 1.;// Only CoulombBarrier can reflect !!! O16 E-dep of gamA ////else if(B>0) wD=0.; // Only Coulomb Barrier can reflect !!! ////else if(B==1) wD=0.; else if(B==1&&C==1) wD=G4QNucleus(Z-1,N,S).GetGSMass()+mProt-GSM; else if(B==1&&C==0) wD=G4QNucleus(Z,N-1,S).GetGSMass()+mNeut-GSM; ////else if(B==1&&C==0) wD=0.; ////else if(B>1) return 1.; ////else if(B>1) wD=0.; ////else if(B==2) wD=0.; else if(B==2) { if (!C) wD=G4QNucleus(Z,N-2,S).GetGSMass()+dNeut-GSM; else if(C==1) wD=G4QNucleus(Z-1,N-1,S).GetGSMass()+mDeut-GSM; else if(C==2) wD=G4QNucleus(Z-2,N,S).GetGSMass()+dProt-GSM; // @@ Temporary "Local B=2 Anti Virial factor" @@ __________________ wD=80.; // 40 MeV per each nucleon //wD=wD/2; // @@ End of Temporary^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ } ////else if(B>2) wD=0.; ////else if(B>2) return 1.; ////else if(B==3) wD=0.; //else if(B==3&&C==1) wD=G4QNucleus(Z-1,N-2,S).GetGSMass()+mTrit-GSM; ////else if(B==3&&C==1) wD=0.; //else if(B==3&&C==2) wD=G4QNucleus(Z-2,N-1,S).GetGSMass()+mHel3-GSM; ////else if(B>3) wD=0.; ////else if(B==4) wD=0.; //else if(B==4&&C==2) wD=G4QNucleus(Z-2,N-2,S).GetGSMass()+mAlph-GSM; ////else if(B>4) wD=0.; ////else if(B>4) return 1.; //else if(B>4)wD=G4QNucleus(Z-C,N-B+C,S).GetGSMass()+G4QNucleus(C,B-C,S).GetGSMass()-GSM; if(wD<0.) wD=0.; #ifdef debug G4cout<<"G4QNucl::CBPenProb: wD="< G4QNucleus::ChooseImpactXandY(G4double maxImpact) { G4double x=1.; G4double y=1.; do { x = G4UniformRand(); x = x+x-1.; y = G4UniformRand(); y = y+y-1.; } while(x*x+y*y > 1.); std::pair theImpactParameter; theImpactParameter.first = x*maxImpact; theImpactParameter.second = y*maxImpact; return theImpactParameter; } // End of "ChooseImpactXandY" // Initializes 3D Nucleons in the nucleus using random sequence void G4QNucleus::ChooseNucleons() { #ifdef debug G4cout<<"G4QNucleus::ChooseNucleons: Nucleons search is started"< 0.) // LOOP over all nucleons { rPos = maxR*pow(G4UniformRand(),third); // Get random radius of the nucleon position G4double density=rPos*rPos; // Density at R (temporary squared radius) if(theA<17) density=GetRelOMDensity(density); // Oscilator model (M.K.?) else density=GetRelWSDensity(rPos); // Wood-Saxon model #ifdef debug G4cout<<"G4QNucl::ChoosePositions: i="<> fill N["< 2) ReduceSum(places,sumPos); // Reduce the CM shift (equal weights) #ifdef debug G4cout<<"G4QNucl::ChoosePositions: The reduced summ is made"<SetPosition(places[i]); delete [] places; #ifdef debug G4cout << "G4QNucleus::ChoosePositions: The positions are defined for A="<0 && dp[i] 0.) { sum/=theA; for(G4int i=0; iSetBindingEnergy(Ebind); // @@ ? M.K. currentNucleon=0; // Automatically starts the LOOP return; } // End of Init3D // Get radius of the most far nucleon (+ nucleon radius) G4double G4QNucleus::GetOuterRadius() { G4double maxradius2=0; G4int theA=theNucleons.size(); if(theA) for(G4int i=0; iGetPosition().mag2(); if(nucr2 > maxradius2) maxradius2=nucr2; } return sqrt(maxradius2)+nucleonDistance; } // End of GetOuterRadius // void G4QNucleus::DoLorentzContraction(const G4ThreeVector& theBeta) { G4double bet2=theBeta.mag2(); G4double factor=(1.-sqrt(1.-bet2))/bet2; // 1./(beta2*gamma2) G4int theA=theNucleons.size(); if(theA) for (G4int i=0; i< theA; i++) { G4ThreeVector pos=theNucleons[i]->GetPosition(); pos -= factor*(theBeta*pos)*theBeta; theNucleons[i]->SetPosition(pos); } } // End of DoLorentzContraction(G4ThreeVector) // Shift all nucleons of the 3D Nucleus (Used in GHAD-TFT) void G4QNucleus::DoTranslation(const G4ThreeVector& theShift) { G4int theA=theNucleons.size(); if(theA) for(G4int i=0; iSetPosition(theNucleons[i]->GetPosition() + theShift); } // End of DoTranslation // Initializes currentNucleon=0 returns size of theNucleons vector G4bool G4QNucleus::StartLoop() { G4int theA=theNucleons.size(); if(theA) currentNucleon=0; else G4cout<<"-Warning-G4QNucleus::StartLoop: LOOP starts for uninited nucleons"<mT) { //Tb->push_back(T); // Fill the thickness vector starting with b=0 Tb.push_back(T); // Fill the thickness vector starting with b=0 b+=db; // increment impact parameter G4double E=B*std::exp(-D*b*b); // b-dependent factor T=C*E/(1.+E); // T(b) in fm^-2 } TbActive=true; // Flag of activation } // End of "ActivateBThickness" // Calculate the integral of T(b) G4double G4QNucleus::GetTbIntegral() // Calculate the integral of T(b) //================================== { if(!TbActive) ActivateBThickness(); G4int nt = Tb.size(); G4double sum=0.; for(G4int i=0; i return 0"<GetPDGCode(); // Get PDG code of the Residual Nucleus G4int theBN = qH->GetBaryonNumber();// Baryon number of the nucleus G4QContent theQC = qH->GetQC(); // Quark Content of the hadron if((thePDG || thePDG==10) && theQC.GetBaryonNumber()>0) thePDG=theQC.GetZNSPDGCode(); G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of theTotalNucleus if(!theBN || thePDG<80000000 || thePDG==90000000) // Hadron, anti-nucleous, or vacuum { #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus: Nucleus="<Get4Momentum()<Hyper Change=>G4QNucleus::EvaporateNuceus: NewNucPDG="< PDG="<"<142.) // @@ to avoid more precise calculations { if(thePDG==90001000) // p* -> n + pi+ { gsM=mNeut; thePDG=90000001; decPDG=211; decM=mPi; } else if(thePDG==90000001) // n* -> p + pi- { gsM=mProt; thePDG=90001000; decPDG=-211; decM=mPi; } else // decay in Pi0 { decPDG=111; decM=mPi0; } } G4LorentzVector h4Mom(0.,0.,0.,gsM); // GSMass must be updated in previous while-LOOP G4LorentzVector g4Mom(0.,0.,0.,decM); if(!G4QHadron(q4M).DecayIn2(h4Mom, g4Mom)) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (3) qH=0"<tM="<"<size()<3)// One can easily close this decay as it will be done later (time of calc?) { G4double gsM=mNeut; if(thePDG==89999000) gsM=mProt; else if(thePDG==89000000) gsM=mLamb; if(fabs(totMass-gsM)<.001) { #ifdef debug G4cout<<"G4QNu::EvaNucl:(aB*),GSM="<"<142.) // @@ to avoid more precise calculations { if(thePDG==89999000) // anti (p* -> n + pi+) { gsM=mNeut; thePDG=89999999; decPDG=-211; decM=mPi; } else if(thePDG==89999999) // anti (n* -> p + pi-) { gsM=mProt; thePDG=89999000; decPDG=211; decM=mPi; } else // decay in Pi0 { decPDG=111; decM=mPi0; } } G4LorentzVector h4Mom(0.,0.,0.,gsM); // GSMass must be updated in previous while-LOOP G4LorentzVector g4Mom(0.,0.,0.,decM); if(!G4QHadron(q4M).DecayIn2(h4Mom, g4Mom)) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (3a) qH=0"<tM="<"<size()<1080.) // @@ to avoid threshold //else if(2>3)// One can easily close this decay as it will be done later (time of calc?) { G4double gsM=mNeut; G4int barPDG=2112; G4int mesPDG=-211; if(thePDG==90001999) { gsM=mProt; barPDG=2212; mesPDG=211; } if(fabs(totMass-gsM-mPi)<.001) { #ifdef debug G4cout<<"G4QN::EvaporateNuc:(D)GSM="<Get4Momentum()<push_back(curB); // (delete equivalent) G4LorentzVector g4Mom=q4M*(mPi/totMass); G4QHadron* curM = new G4QHadron(mesPDG,g4Mom); evaHV->push_back(curM); // (delete equivalent) #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (5) qH=0"<tM="<"<1080.) // @@ to avoid threshold //else if(2>3)// One can easily close this decay as it will be done later (time of calc?) { G4double gsM=mNeut; G4int barPDG=-2112; G4int mesPDG=211; if(thePDG==89998001) { gsM=mProt; barPDG=-2212; mesPDG=-211; } if(fabs(totMass-gsM-mPi)<.001) { #ifdef debug G4cout<<"G4QN::EvaporateNuc:(A)GSM="<Get4Momentum()<push_back(curB); // (delete equivalent) G4LorentzVector g4Mom=q4M*(mPi/totMass); G4QHadron* curM = new G4QHadron(mesPDG,g4Mom); evaHV->push_back(curM); // (delete equivalent) #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (5a) qH=0"<tM="<"<0&&theS<0) DecayAntiStrange(qH,evaHV); // "AntyStrangeNucleus" (del.eq.) else if(theBN>0&&(theC<0||theC>theBN-theS))DecayIsonucleus(qH,evaHV);//"Isonucleus"(d.e.) else if((thePDG==89999003 || thePDG==90002999) && totMass>2020.) //=> "ISO-dibarion" { // @@ Check that it never comes here !! G4int nucPDG = 2112; G4double nucM = mNeut; G4int piPDG = -211; if(thePDG==90002999) { nucPDG = 2212; nucM = mProt; piPDG = 211; } if(totMass>mPi+nucM+nucM) { G4LorentzVector n14M(0.,0.,0.,nucM); G4LorentzVector n24M(0.,0.,0.,nucM); G4LorentzVector pi4M(0.,0.,0.,mPi); if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M)) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (9) qH=0"< 2N="<2020.) //=> IsoBP { // @@ Pi0 can be taken into account ! G4int n1PDG = 2212; G4int n2PDG = 2112; G4int piPDG = -211; G4double n1M = mProt; G4double n2M = mNeut; if (thePDG==90002000) piPDG = 211; // pp -> np + pi- else if(thePDG==90000002) piPDG = -211; // nn -> np + pi- else // np -> 50%(nnpi+) 50%(pppi-) { if(G4UniformRand()>.5) { n1PDG=2112; n1M=mNeut; piPDG = 211; } else { n2PDG=2212; n2M=mProt; piPDG = -211; } } if(totMass>mPi+n1M+n2M) { G4LorentzVector n14M(0.,0.,0.,n1M); G4LorentzVector n24M(0.,0.,0.,n2M); G4LorentzVector pi4M(0.,0.,0.,mPi); if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M)) { G4cerr<<"***G4QNucl::EvapNucleus:IsoDN, tM="< N1="< "Dibaryon" case (del eq.) else if((thePDG==90000997 || thePDG==89997001) && totMass>2020.) //=> "anti-ISO-dibarion" { // @@ Check that it never comes here !! G4int nucPDG = -2112; G4double nucM = mNeut; G4int piPDG = 211; if(thePDG==90002999) { nucPDG = -2212; nucM = mProt; piPDG = -211; } if(totMass>mPi+nucM+nucM) { G4LorentzVector n14M(0.,0.,0.,nucM); G4LorentzVector n24M(0.,0.,0.,nucM); G4LorentzVector pi4M(0.,0.,0.,mPi); if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M)) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (9a) qH=0"< 2aN="<2020.)//=>AnIsoBP { // @@ Pi0 can be taken into account ! G4int n1PDG = -2212; G4int n2PDG = -2112; G4int piPDG = 211; // dummy initialization G4double n1M = mProt; G4double n2M = mNeut; if (thePDG==89998000) piPDG = -211; // anti ( pp -> np + pi- ) else if(thePDG==89999998) piPDG = 211; // anti ( nn -> np + pi- ) else // anti ( np -> 50%(nnpi+) 50%(pppi-) ) { if(G4UniformRand()>.5) { n1PDG=-2112; n1M=mNeut; piPDG = -211; } else { n2PDG=-2212; n2M=mProt; piPDG = 211; } } if(totMass>mPi+n1M+n2M) { G4LorentzVector n14M(0.,0.,0.,n1M); G4LorentzVector n24M(0.,0.,0.,n2M); G4LorentzVector pi4M(0.,0.,0.,mPi); if(!G4QHadron(q4M).DecayIn3(n14M,n24M,pi4M)) { G4cerr<<"**G4QNucl::EvapNucleus:IsoDN,antiM="< N1="< "Anti-Dibaryon" case (del eq.) else if(fabs(totMass-totGSM)<.001) // Fill as it is or decay Be8, He5, Li5 (@@ add more) { #ifdef debug G4cout<<"G4QNucleus::EvaNuc:GS "<GetQC()<Get4Momentum()<<" FillAsIs"<1) { DecayMultyBaryon(qH,evaHV); } else if(theBN==5) { DecayAlphaBar(qH,evaHV); } // Decay unstable A5 system (del eq.) else { evaHV->push_back(qH); } // Fill as it is (del eq.) } else if(theBN>1 && thePDG>88000000 && thePDG<89000000) //==> 2antiK in the nucleus { G4cout<<"---Warning---G4QNucl::EvaNuc:MustNotBeHere.PDG="<=bN) nucPDG+=999000; else { nucPDG+=999999; k1PDG = 311; mK1= mK0; } if(bZ>bN) nucPDG+=999000; else { nucPDG+=999999; k2PDG = 311; mK2= mK0; } G4double nucM = G4QNucleus(nucPDG).GetGSMass(); G4cout<<"-Warning-G4QN::EvN:M="< N="< " < " <GetQC()<Get4Momentum()<<" AsIs"<push_back(qH); return; } else if ( ( bA == 1 || (!bsCond && !dbsCond) ) && totMass > GSMass+.003 )//=>Fuse&Decay //else if(2>3) // Close "Fuse&Decay Technology"***@@@*** { #ifdef debug G4cout<<"G4QN::EvaN:SplitBar, s="< GSM="<size(); // Total#of QHadrons in Vector at this point G4bool bnfound=true; // Cure "back fusion fragment not found" while(nOfOUT) // Try BackFusionDecays till something is in { G4QHadron* theLast = (*evaHV)[nOfOUT-1]; G4int lastBN = theLast->GetBaryonNumber(); G4int nFrag = theLast->GetNFragments(); //////////////////G4int gam = theLast->GetPDGCode(); #ifdef debug G4cout<<"G4QN::EvaNuc:*BackFus*,BN="< "Delete Last Decayed Hadron" case { G4QHadron* thePrev = (*evaHV)[nOfOUT-2]; nFrag = thePrev->GetNFragments(); G4int prevBN = thePrev->GetBaryonNumber(); #ifdef debug G4int prevPDG = thePrev->GetPDGCode(); G4cout<<"G4QNucl::EvaNucl: DelTheLast, nFr="<-0.001) // Decay in the qH and the last is impossible { G4QHadron* evH = new G4QHadron(totPDG,q4M); // Create QHadron for CompResidNuc if(dM<=0.) { evaHV->pop_back(); // lastQHadron is excluded from QHadrV asIs in TRN delete theLast; //When kill, DON'T forget to delete lastQHadron asAnInstance! #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (16) qH=0"<push_back(evH);// Fill TRN to HVect asIs (delete equivalent) } else // Decay TotalResidualNucleus in GSM+Last { G4LorentzVector r4Mom(0.,0.,0.,GSMass); if(!G4QHadron(q4M).DecayIn2(last4M,r4Mom)) // Decay failed { evaHV->pop_back(); // lastQHadron is excluded from QHadrV as is in TRN delete theLast; //When kill,DON'T forget to delete lastQHadron asAnInstance #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (17) qH=0"<push_back(evH);// Fill TRN to Vect as it is (delete equivalent) #ifdef debug G4cout<<"***G4QN::EvaN:DecayIn L"<Set4Momentum(last4M);// Already exists:don't create&fill,->set4Mom G4QHadron* nucH = new G4QHadron(thePDG,r4Mom); // Create QHadron for qH-nuc #ifdef debug G4cout<<"G4QNucleus::EvaNuc:fill NH "<push_back(nucH);// Fill the Residual Nucleus (del.eq.) } } bnfound=false; break; } thePDG=totPDG; // Make ResidualNucleus outOf theTotResidualNucl GSMass=G4QPDGCode(thePDG).GetMass();// Update the Total Residual Nucleus mass evaHV->pop_back(); // the last QHadron is excluded from OUTPUT delete theLast;//!!When kill,DON'T forget to delete theLastQHadron asAnInstance!! nOfOUT--; // Update the value of OUTPUT entries } } if(bnfound) { G4LorentzVector h4Mom(0.,0.,0.,GSMass);//GSMass must be updated inPreviousWhileLOOP G4LorentzVector g4Mom(0.,0.,0.,0.); if(!G4QHadron(q4M).DecayIn2(h4Mom, g4Mom)) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (19) qH=0"<tM="<totResN="<0&&bS<0) DecayAntiStrange(qH,evaHV);// Decay nucleus with antistrangeness else if(bA==2) DecayDibaryon(qH,evaHV); // Decay the residual dibaryon (del.equivalent) else if(bA==-2) DecayAntiDibaryon(qH,evaHV); // Decay residual anti-dibaryon (del.eq) else if(totMass" M0&&bA>1) // There's aNeutr in theNucl, which can be split { G4QContent resQC=totQC-neutQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus nResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (nResPDG==90000001) nResM=mNeut; else if(nResPDG==90001000) nResM=mProt; else if(nResPDG==91000000) nResM=mLamb; else nResM=resN.GetMZNS(); // min mass of the Residual Nucleus } G4double pResM =1000000.; // Prototype of mass of residual for a proton G4int pResPDG=0; // Prototype of PDGCode of residual for a proton if(bsCond==2212&&bZ>0&&bA>1) // There's aProton in Nucl, which can be split { G4QContent resQC=totQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus pResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (pResPDG==90000001) pResM=mNeut; else if(pResPDG==90001000) pResM=mProt; else if(pResPDG==91000000) pResM=mLamb; else pResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double lResM =1000000.; // Prototype of mass of residual for a Lambda G4int lResPDG=0; // Prototype of PDGCode of residual for a Lambda if(bsCond==3122&&bS>0&&bA>1) // There's aLambd in theNucl, which can be split { G4QContent resQC=totQC-lambQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus lResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (lResPDG==90000001) lResM=mNeut; else if(lResPDG==90001000) lResM=mProt; else if(lResPDG==91000000) lResM=mLamb; else lResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double dResM =1000000.; // Prototype of mass of residual for a Alpha G4int dResPDG=0; // Prototype of PDGCode of residual for a Alpha if(bsCond==90001001&&bN>0&&bZ>0&&bA>2)// There's aDeuter in Nucl, which canBeRadiated { G4QContent resQC=totQC-deutQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus dResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (dResPDG==90000001) dResM=mNeut; else if(dResPDG==90001000) dResM=mProt; else if(dResPDG==91000000) dResM=mLamb; else dResM =resN.GetMZNS(); // minMass of the Residual Nucleus } G4double aResM =1000000.; // Prototype of mass of residual for a Alpha G4int aResPDG=0; // Prototype of PDGCode of residual for a Alpha if(bsCond==90002002&&bN>1&&bZ>1&&bA>4)// There's Alpha in theNucl, which can be split { G4QContent resQC=totQC-alphQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus aResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (aResPDG==90000001) aResM=mNeut; else if(aResPDG==90001000) aResM=mProt; else if(aResPDG==91000000) aResM=mLamb; else aResM =resN.GetMZNS(); // minMass of the Residual Nucleus } G4double nnResM =1000000.; // Prototype of mass of residual for a dineutron G4int nnResPDG=0; // Prototype of ResidualPDGCode for a dineutron if(dbsCond&&bN>1&&bA>2) // It's nucleus and there is a dineutron { G4QContent resQC=totQC-neutQC-neutQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus nnResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (nnResPDG==90000001) nnResM=mNeut; else if(nnResPDG==90001000) nnResM=mProt; else if(nnResPDG==91000000) nnResM=mLamb; else nnResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double ppResM =1000000.; // Prototype of mass of residual for a diproton G4int ppResPDG=0; // Prototype of ResidualPDGCode for a diproton if(dbsCond&&bZ>1&&bA>2) // It's nucleus and there is a diproton { G4QContent resQC=totQC-protQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus ppResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (ppResPDG==90000001) ppResM=mNeut; else if(ppResPDG==90001000) ppResM=mProt; else if(ppResPDG==91000000) ppResM=mLamb; else ppResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double npResM =1000000.; // Prototype of ResidualMass for proton+neutron G4int npResPDG=0; // Prototype of ResidualPDGCode for a prot+neut if(dbsCond&&bN>0&&bZ>0&&bA>2) // It's nucleus and there is aProton & aNeutron { G4QContent resQC=totQC-neutQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus npResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (npResPDG==90000001) npResM=mNeut; else if(npResPDG==90001000) npResM=mProt; else if(npResPDG==91000000) npResM=mLamb; else npResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double lnResM =1000000.; // Prototype of residualMass for lambda+neutron G4int lnResPDG=0; // Prototype of ResidualPDGCode for aLambda+Neut if(dbsCond&&bN>0&&bS>0&&bA>2) // It's nucleus and there is aLambda & aNeutron { G4QContent resQC=totQC-lambQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus lnResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (lnResPDG==90000001) lnResM=mNeut; else if(lnResPDG==90001000) lnResM=mProt; else if(lnResPDG==91000000) lnResM=mLamb; else lnResM =resN.GetMZNS(); // min mass of the Residual Nucleus } G4double lpResM =1000000.; // Prototype of residualMass for a proton+lambda G4int lpResPDG=0; // Prototype of ResidualPDGCode for theProt+lamb if(dbsCond&&bS>0&&bZ>0&&bA>2) // It's nucleus and there is aProton and aLambda { G4QContent resQC=totQC-neutQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus lpResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (lpResPDG==90000001) lpResM=mNeut; else if(lpResPDG==90001000) lpResM=mProt; else if(lpResPDG==91000000) lpResM=mLamb; else lpResM =resN.GetMZNS(); // minMass of the Residual Nucleus } G4double llResM =1000000.; // Prototype of mass of residual for a di-lambda G4int llResPDG=0; // Prototype of ResidPDGCode for the di-lambda if(dbsCond&&bS>1&&bA>2) // It's nucleus and there is a di-lambda { G4QContent resQC=totQC-neutQC-protQC; G4QNucleus resN(resQC); // Pseudo nucleus for the Residual Nucleus llResPDG=resN.GetPDG(); // PDG of the Residual Nucleus if (llResPDG==90000001) llResM=mNeut; else if(llResPDG==90001000) llResM=mProt; else if(llResPDG==91000000) llResM=mLamb; else llResM =resN.GetMZNS(); // min mass of the Residual Nucleus } #ifdef debug G4cout<<"G4QNucleus::EvaNucl: rP="<DecayIn2(a4Mom,b4Mom)) { evaHV->push_back(qH); // Fill as it is (delete equivalent) G4cout<<"---Warning---G4QNucleus::EvaNuc:rP="<push_back(qH); // FillAsItIs (del.eq.) return; } else // "System is below mass shell and can't decay" case { #ifdef debug G4cout<<"***G4QNucl::EvaNuc: tM="<push_back(qH); // Correct or fill as it is return; } } else // ===> Evaporation of the excited system { #ifdef pdebug G4cout<<"G4QN::EvaNuc:***EVA***tPDG="<GSM="< "< bCB="<push_back(bHadron); // Fill EvaporatedBaryon (del.equivalent) else if(bB==2) DecayDibaryon(bHadron,evaHV);// => "Dibaryon" case needs decay else if(bB==4) evaHV->push_back(bHadron); // "Alpha radiation" case (del.eq.) else if(bB==5) DecayAlphaBar(bHadron,evaHV);// "Alpha+Baryon Decay" case (del.equiv.) else if(bPDG==90004002) DecayAlphaDiN(bHadron,evaHV); // alph+2p(alph+2n is stable) else if(bPDG==90004004) DecayAlphaAlpha(bHadron,evaHV);// Alph+Alph Decay (del.eq.) else { delete bHadron; G4cerr<<"***G4QNuc::EvaNuc:bB="<2 - unexpected evaporated fragment"<2) EvaporateNucleus(rHadron,evaHV); // Continue evaporation (@@ Self-call) else if(rB==2) // => "Dibaryon" case needs decay @@ DecayDibaryon { G4double rGSM = rHadron->GetQPDG().GetMass(); // Ground State mass of the dibaryon #ifdef debug G4cout<<"G4QNuc::EvaNuc:ResidDibM="< M="<push_back(rHadron); // (DE) else DecayDibaryon(rHadron,evaHV); // => "Dibaryon Decay" case (del.equivalent) } else if(rB==5) DecayAlphaBar(rHadron,evaHV);// "Alpha+Baryon Decay" case (del.equiv.) else if(rPDG==90004002) DecayAlphaDiN(rHadron,evaHV);//alph+2p (alph+2n is stable) else if(rPDG==90004004) DecayAlphaAlpha(rHadron,evaHV);//Alpha+Alpha Decay (delEq) else evaHV->push_back(rHadron); // Fill ResidNucl=Baryon to OutputHadronVector } // End of Evaporation of excited system #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus: === End of the evaporation attempt"< "Decay if impossible evaporate" case { #ifdef debug G4cout<<"*G4QNucleus::EvaporateNucleus: InputHadron4M="<GetQC(); // Quark content of the hadron G4QChipolino resChip(totQC); // define the Residual as a Chipolino G4QPDGCode h1=resChip.GetQPDG1(); G4double m1 =h1.GetMass(); // Mass of the first hadron G4QPDGCode h2=resChip.GetQPDG2(); G4double m2 =h2.GetMass(); // Mass of the second hadron if(totMass+.0001>m1+m2) { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (24) qH=0"< h1M="<h1+h2 DecIn2 error"); } G4QHadron* H2 = new G4QHadron(h2.GetPDGCode(),qe4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(2) h2="<push_back(H2); // (delete equivalent) G4QHadron* H1 = new G4QHadron(h1.GetPDGCode(),fq4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(2) h1="<push_back(H1); // (delete equivalent) } else { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (25) qH=0"<2) { #ifdef debug G4cout<<"**G4QNuc::EvaNuc:EmerFill(2) "<GetQC()<Get4Momentum()<push_back(qH); } else if ((thePDG==221||thePDG==331)&&totMass>mPi+mPi) // "Decay in pipi" case { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (26) qH=0"< pi+ + pi-"<Pi+Pi DecayIn2 error"); } G4QHadron* H1 = new G4QHadron(211,fq4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(3) PiPlus="<push_back(H1); // (delete equivalent) G4QHadron* H2 = new G4QHadron(-211,qe4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(3) PiMinus="<push_back(H2); // (delete equivalent) } else if ((thePDG==221||thePDG==331)&&totMass>mPi0+mPi0) // "Decay in 2pi0" case { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (27) qH=0"< pi0 + pi0"<Pi+Pi DecayIn2 error"); } G4QHadron* H1 = new G4QHadron(111,fq4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(4) Pi01="<push_back(H1); // (delete equivalent) G4QHadron* H2 = new G4QHadron(111,qe4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(4) Pi02="<push_back(H2); // (delete equivalent) } else if (totMass>totM) // "Radiative Hadron decay" case { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (28) qH=0"<h1M="<H+gamma DecIn2 error"); } G4QHadron* H2 = new G4QHadron(thePDG,qe4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(5) tot="<push_back(H2); // (delete equivalent) G4QHadron* H1 = new G4QHadron(22,fq4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(5) GamFortot="<push_back(H1); // (delete equivalent) } else if (thePDG==111||thePDG==221||thePDG==331) // "Gamma+Gamma decay" case { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (29) qH=0"<gamma + gamma"<g+g DecIn2 error"); } G4QHadron* H2 = new G4QHadron(22,qe4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(6) gam1="<push_back(H2); // (delete equivalent) G4QHadron* H1 = new G4QHadron(22,fq4M); #ifdef debug G4cout<<"G4QNucleus::EvaporateNucleus:(6) gam2="<push_back(H1); // (delete equivalent) } else { #ifdef qdebug if(!qH) G4cout<<"-Warning-G4QNucleus::EvaporateNucleus: (30) qH=0"<>>> End. "<Get4Momentum(); // Get 4-momentum of the Isonucleus G4double qM=q4M.m(); // Real mass of the Isonucleus G4QContent qQC = qH->GetQC(); // Get QuarcContent of the Isonucleus G4int qBN=qQC.GetBaryonNumber(); // Baryon number of the IsoNucleus G4int qC=qQC.GetCharge(); // Charge of the IsoNucleus G4int qS=qQC.GetStrangeness(); // Strangness of the IsoNucleus #ifdef debug G4cout<<"G4QNuc:DecayIson:QC="<qBN) // *** Should not be here *** { G4cout<<"--Warning(Upgrade)--G4QNuc::DecIsonuc:FillAsIs,4M="<push_back(qH); // fill as it is (delete equivalent) return; } G4int qPN=qC-qBN; // Number of pions in the Isonucleus G4int fPDG = 2212; // Prototype for nP+(Pi+) case G4int sPDG = 211; G4int tPDG = 3122; // @@ Sigma0 (?) G4double fMass= mProt; G4double sMass= mPi; G4double tMass= mLamb; // @@ Sigma0 (?) // ========= Negative state ====== if(qC<0) // ====== Only Pi- can help { if(qS&&qBN==qS) // --- n*Lamb + k*(Pi-) State --- { sPDG = -211; if(-qC==qS && qS==1) // Only one Sigma- like (qBN=1) { if(fabs(qM-mSigM)push_back(qH); // Fill Sigma- as it is return; } else if(qM>mLamb+mPi) //(2) Sigma- => Lambda + Pi- decay { fPDG = 3122; fMass= mLamb; } else if(qM>mSigM) //(2) Sigma+ => Sigma+ + gamma decay { fPDG = 3112; fMass= mSigM; sPDG = 22; sMass= 0.; } else //(2) Sigma- => Neutron + Pi- decay { fPDG = 2112; fMass= mNeut; } qPN = 1; // #of (Pi+ or gamma)'s = 1 } else if(-qC==qS) //(2) a few Sigma- like { qPN = 1; // One separated Sigma- fPDG = 3112; sPDG = 3112; sMass= mSigM; qBN--; fMass= mSigM; } else if(-qC>qS) //(2) n*(Sigma-)+m*(Pi-) { qPN = -qC-qS; // #of Pi-'s fPDG = 3112; fMass= mSigM; } else //(2) n*(Sigma-)+m*Lambda (-qCn*Lamb+m*Neut+k*(Pi-) State (qSnPin) //(3) m*P+n*(Sigma+)+k*Lambda { qS-=nPin; // #of Lambdas qPN = nPin; // #of Sigma+ sPDG = 3112; sMass= mSigM; } else //(3) m*N+n*(Sigma-)+k*(Pi-) (qS1 && qBN==qS) //(2) m*Lamb(m>1) ***Should not be here*** { qPN = 1; fPDG = 3122; sPDG = 3122; sMass= mLamb; qBN--; fMass= mLamb; } else if(!qS && qBN>1) //(2) n*Neut(n>1) ***Should not be here*** { qPN = 1; fPDG = 2112; sPDG = 2112; sMass= mNeut; qBN--; fMass= mNeut; } else G4cout<<"*?*G4QNuc::DecayIsonucleus: (1) QC="<0) // n*Lamb+(m*P)+(k*Pi+) { if(qS && qS+qC==qBN) //(2) n*Lamb+m*P ***Should not be here*** { qPN = qS; qS = 0; fPDG = 2212; sPDG = 3122; sMass= mLamb; qBN = qC; fMass= mProt; } else if(qS && qC n*L+m*Pi+ State { if(qC==qS && qS==1) // Only one Sigma+ like State { if(fabs(qM-mSigP)push_back(qH); return; } else if(qM>mLamb+mPi) //(2) Sigma+ => Lambda + Pi+ decay { fPDG = 3122; fMass= mLamb; } else if(qM>mNeut+mPi) //(2) Sigma+ => Neutron + Pi+ decay { fPDG = 2112; fMass= mNeut; } else if(qM>mSigP) //(2) Sigma+ => Sigma+ + gamma decay { fPDG = 3222; fMass= mSigP; sPDG = 22; sMass= 0.; } else //(2) Sigma+ => Proton + gamma decay { fPDG = 2212; fMass= mProt; sPDG = 22; sMass= 0.; } qPN = 1; // #of (Pi+ or gamma)'s = 1 } else if(qC==qS) //(2) a few Sigma+ like hyperons { qPN = 1; fPDG = 3222; sPDG = 3222; sMass= mSigP; qBN--; fMass= mSigP; } else if(qC>qS) //(2) n*(Sigma+)+m*(Pi+) { qPN = qC-qS; // #of Pi+'s fPDG = 3222; qBN = qS; // #of Sigma+'s fMass= mSigP; } else //(2) n*(Sigma+)+m*Lambda { qBN -= qC; // #of Lambda's fPDG = 3122; fMass= mLamb; qPN = qC; // #of Sigma+'s sPDG = 3222; sMass= mSigP; } qS = 0; // All above are decays in 2 } else if(qS && qC>qBN-qS) // n*Lamb+m*P+k*Pi+ { qBN -= qS; // #of protons G4int nPip = qC-qBN; // #of Pi+'s if(qS==nPip) //(2) m*P+n*Sigma+ { qPN = qS; // #of Sigma+ sPDG = 3222; sMass= mSigP; qS = 0; } else if(qS>nPip) //(3) m*P+n*(Sigma+)+k*Lambda { qS -= nPip; // #of Lambdas qPN = nPip; // #of Sigma+ sPDG = 3222; sMass= mSigP; } else //(3) m*P+n*(Sigma+)+k*(Pi+) { qPN = nPip-qS; // #of Pi+ tPDG = 3222; tMass= mSigP; } } if(qC1) //(2) m*Prot(m>1) ***Should not be here*** { qPN = 1; fPDG = 2212; sPDG = 2212; sMass= mProt; qBN--; fMass= mProt; } else if(qC<=qBN||!qBN) G4cout<<"*?*G4QNuc::DecayIsonucleus: (2) QC="<qBN //(2) Default condition n*P+m*(Pi+) } G4double tfM=qBN*fMass; G4double tsM=qPN*sMass; G4double ttM=0.; if(qS) ttM=qS*tMass; G4LorentzVector f4Mom(0.,0.,0.,tfM); G4LorentzVector s4Mom(0.,0.,0.,tsM); G4LorentzVector t4Mom(0.,0.,0.,ttM); G4double sum=tfM+tsM+ttM; if(fabs(qM-sum) TotM="<push_back(qH); // fill as it is (delete equivalent) return; } else if(qS && (qM TotM="<push_back(qH); // fill as it is (delete equivalent) return; } #ifdef debug G4cout<<"G4QNuc::DecayIsonucleus: *DONE* n="<push_back(Hj); // Fill "Hj" (delete equivalent) } } if(qS) { t4Mom/=qS; for(G4int il=0; ilpush_back(Hk); // Fill "Hk" (delete equivalent) } } #ifdef qdebug if (qH) { G4cout << "G4QNucleus::DecayIsonucleus: deleted at end - PDG: " << qH->GetPDGCode() << G4endl; delete qH; } #endif } // End of DecayIsonucleus //Decay of the excited dibaryon in two baryons void G4QNucleus::DecayDibaryon(G4QHadron* qH, G4QHadronVector* evaHV) {// ================================================================ static const G4double mPi = G4QPDGCode(211).GetMass(); static const G4double mNeut= G4QPDGCode(2112).GetMass(); static const G4double mProt= G4QPDGCode(2212).GetMass(); static const G4double mSigM= G4QPDGCode(3112).GetMass(); static const G4double mLamb= G4QPDGCode(3122).GetMass(); static const G4double mSigP= G4QPDGCode(3222).GetMass(); static const G4double mKsiM= G4QPDGCode(3312).GetMass(); static const G4double mKsiZ= G4QPDGCode(3322).GetMass(); static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0); static const G4double mPiN = mPi+mNeut; static const G4double mPiP = mPi+mProt; static const G4double dmPiN= mPiN+mPiN; static const G4double dmPiP= mPiP+mPiP; static const G4double nnPi = mNeut+mPiN; static const G4double ppPi = mProt+mPiP; static const G4double lnPi = mLamb+mPiN; static const G4double lpPi = mLamb+mPiP; static const G4double dNeut= mNeut+mNeut; static const G4double dProt= mProt+mProt; static const G4double dLamb= mLamb+mLamb; static const G4double dLaNe= mLamb+mNeut; static const G4double dLaPr= mLamb+mProt; static const G4double dSiPr= mSigP+mProt; static const G4double dSiNe= mSigM+mNeut; static const G4double dKsPr= mKsiZ+mProt; static const G4double dKsNe= mKsiM+mNeut; static const G4double eps = 0.003; static const G4QNucleus vacuum(90000000); G4bool four=false; // defFALSE for 4-particle decay of diDelta G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the Dibaryon G4int qPDG = qH->GetPDGCode(); // PDG Code of the decaying dybaryon G4double qM = q4M.m(); G4double rM = qM+eps; // Just to avoid the computer accuracy #ifdef debug G4cout<<"G4QNucl::DecayDibaryon: *Called* PDG="<=dmPiN) // "diDelta--" case { sPDG = -211; fPDG = 2112; sMass= mPi; fMass= mNeut; four = true; } else if(qPDG==90000002 && rM>=dNeut) // "dineutron" case { fPDG = 2112; sPDG = 2112; fMass= mNeut; sMass= mNeut; } else if(qPDG==90001001 && rM>=mDeut) // "exited deutron" case { if(fabs(qM-mDeut)push_back(qH); // Fill as it is (delete equivalent) return; } else if(mProt+mNeut=dLaNe) // "Lambda-neutron" case { fPDG = 2112; sPDG = 3122; fMass= mNeut; sMass= mLamb; } else if(qPDG==91001000 && rM>=dLaPr) // "Lambda-proton" case { sPDG = 3122; sMass= mLamb; } else if(qPDG==89999003 && rM>=nnPi) // "neutron/Delta-" case { fPDG = 2112; sPDG = 2112; tPDG = -211; fMass= mNeut; sMass= mNeut; } else if(qPDG==90002999 && rM>=ppPi) // "proton/Delta++" case { tPDG = 211; } else if(qPDG==90999002 && rM>=lnPi) // "lambda/Delta-" case { fPDG = 2112; sPDG = 3122; tPDG = -211; fMass= mNeut; sMass= mLamb; } else if(qPDG==91001999 && rM>=lpPi) // "lambda/Delta+" case { sPDG = 3122; tPDG = 211; sMass= mLamb; } else if(qPDG==90999002 && rM>=dSiNe) // "Sigma-/neutron" case { fPDG = 2112; sPDG = 3112; fMass= mNeut; sMass= mSigM; } else if(qPDG==91001999 && rM>=dSiPr) // "Sigma+/proton" case { sPDG = 3222; sMass= mSigP; } else if(qPDG==92000000 && rM>=dLamb) // "diLambda" case { fPDG = 3122; sPDG = 3122; fMass= mLamb; sMass= mLamb; } else if(qPDG==91999001 && rM>=dKsNe) // "Ksi-/neutron" case { fPDG = 2112; sPDG = 3312; fMass= mNeut; sMass= mKsiM; } else if(qPDG==92000999 && rM>=dKsPr) // "Ksi0/proton" case { sPDG = 3322; sMass= mKsiZ; } else if(qPDG!=90002000|| rMGetStrangeness(); G4int qB = qH->GetBaryonNumber(); if(qB>0&&qS<0) // Antistrange diBarion { DecayAntiStrange(qH,evaHV); return; } else { delete qH; G4cerr<<"***G4QN::DecDiBar: badPDG="<? TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayDibaryon:(2) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) } else if(four) { q4M=q4M/2.; // Divided in 2 !!! qM/=2.; // Divide the mass in 2 ! G4double sum=fMass+sMass; if(fabs(qM-sum)tM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayDibaryon:(3) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) // Now the second pair mus be decayed if(fabs(qM-sum)? (DD2,Can't be here) TotM="<push_back(H4); // Fill "H2" (delete equivalent) delete qH; } else { G4double sum=fMass+sMass+tMass; if(fabs(qM-sum)TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNuc::DecayDibaryon:(5) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) G4QHadron* H3 = new G4QHadron(tPDG,t4Mom); // Create a Hadron for the meson evaHV->push_back(H3); // Fill "H3" (delete equivalent) } #ifdef qdebug if (qH) { G4cout << "G4QNucleus::DecayDiBaryon: deleted at end - PDG: " << qH->GetPDGCode() << G4endl; delete qH; } #endif } // End of DecayDibaryon //Decay of the excited anti-dibaryon in two anti-baryons void G4QNucleus::DecayAntiDibaryon(G4QHadron* qH, G4QHadronVector* evaHV) {// ==================================================================== static const G4double mPi = G4QPDGCode(211).GetMass(); static const G4double mNeut= G4QPDGCode(2112).GetMass(); static const G4double mProt= G4QPDGCode(2212).GetMass(); static const G4double mSigM= G4QPDGCode(3112).GetMass(); static const G4double mLamb= G4QPDGCode(3122).GetMass(); static const G4double mSigP= G4QPDGCode(3222).GetMass(); static const G4double mKsiM= G4QPDGCode(3312).GetMass(); static const G4double mKsiZ= G4QPDGCode(3322).GetMass(); static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0); static const G4double mPiN = mPi+mNeut; static const G4double mPiP = mPi+mProt; static const G4double dmPiN= mPiN+mPiN; static const G4double dmPiP= mPiP+mPiP; static const G4double nnPi = mNeut+mPiN; static const G4double ppPi = mProt+mPiP; static const G4double lnPi = mLamb+mPiN; static const G4double lpPi = mLamb+mPiP; static const G4double dNeut= mNeut+mNeut; static const G4double dProt= mProt+mProt; static const G4double dLamb= mLamb+mLamb; static const G4double dLaNe= mLamb+mNeut; static const G4double dLaPr= mLamb+mProt; static const G4double dSiPr= mSigP+mProt; static const G4double dSiNe= mSigM+mNeut; static const G4double dKsPr= mKsiZ+mProt; static const G4double dKsNe= mKsiM+mNeut; static const G4double eps = 0.003; static const G4QNucleus vacuum(90000000); G4bool four=false; // defFALSE for 4-particle decay of diDelta G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the Dibaryon G4int qPDG = qH->GetPDGCode(); // PDG Code of the decaying dybaryon G4double qM = q4M.m(); // Mass of the decaying anti-di-baryon G4double rM = qM+eps; // Just to avoid the computer accuracy #ifdef debug G4cout<<"G4QNucl::DecayAntiDibar:*Called* PDG="<=dmPiN) // "diDelta--" case { sPDG = 211; fPDG = -2112; sMass= mPi; fMass= mNeut; four = true; } else if(qPDG==89999998 && rM>=dNeut) // "dineutron" case { fPDG = -2112; sPDG = -2112; fMass= mNeut; sMass= mNeut; } else if(qPDG==89998999 && rM>=mDeut) // "exited deutron" case { if(fabs(qM-mDeut)push_back(qH); // Fill as it is (delete equivalent) return; } else if(mProt+mNeut=dLaPr) // "Lambda-proton" case { sPDG = -3122; sMass= mLamb; } else if(qPDG==90000997 && rM>=nnPi) // "neutron/Delta-" case { fPDG = -2112; sPDG = -2112; tPDG = 211; fMass= mNeut; sMass= mNeut; } else if(qPDG==89997001 && rM>=ppPi) // "proton/Delta++" case { tPDG = -211; } else if(qPDG==89000998 && rM>=lnPi) // "lambda/Delta-" case { fPDG = -2112; sPDG = -3122; tPDG = 211; fMass= mNeut; sMass= mLamb; } else if(qPDG==889998001 && rM>=lpPi) // "lambda/Delta+" case { sPDG = -3122; tPDG = -211; sMass= mLamb; } else if(qPDG==89000998 && rM>=dSiNe) // "Sigma-/neutron" case { fPDG = -2112; sPDG = -3112; fMass= mNeut; sMass= mSigM; } else if(qPDG==88998001 && rM>=dSiPr) // "Sigma+/proton" case { sPDG = -3222; sMass= mSigP; } else if(qPDG==88000000 && rM>=dLamb) // "diLambda" case { fPDG = -3122; sPDG = -3122; fMass= mLamb; sMass= mLamb; } else if(qPDG==88000999 && rM>=dKsNe) // "Ksi-/neutron" case { fPDG = -2112; sPDG = -3312; fMass= mNeut; sMass= mKsiM; } else if(qPDG==87999001 && rM>=dKsPr) // "Ksi0/proton" case { sPDG = -3322; sMass= mKsiZ; } else if(qPDG!=89998000|| rMGetStrangeness(); G4int qB = qH->GetBaryonNumber(); if(qB>0&&qS<0) // Antistrange diBarion { DecayAntiStrange(qH,evaHV); return; } else { delete qH; G4cerr<<"**G4QNuc::DecayAntiDiBar: badPDG="<? TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayAntiDibaryon:(2) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) } else if(four) { q4M=q4M/2.; // Divided in 2 !!! qM/=2.; // Divide the mass in 2 ! G4double sum=fMass+sMass; if(fabs(qM-sum)tM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayAntiDibaryon:(3) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) // Now the second pair mus be decayed if(fabs(qM-sum)? (DD2,Can't be here) TotM="<push_back(H4); // Fill "H2" (delete equivalent) delete qH; } else { G4double sum=fMass+sMass+tMass; if(fabs(qM-sum)TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNuc::DecayAbtiDibaryon:(5) *DONE* f4M="<push_back(H2); // Fill "H2" (delete equivalent) G4QHadron* H3 = new G4QHadron(tPDG,t4Mom); // Create a Hadron for the meson evaHV->push_back(H3); // Fill "H3" (delete equivalent) } #ifdef qdebug if (qH) { G4cout<<"G4QNucleus::DecayDiBaryon: deleted at end - PDG="<GetPDGCode()<Get4Momentum(); // Get 4-mom of the AntiStrangeNuclearState G4double qM = q4M.m(); // Real mass of the AntiStrangeNuclearState G4QContent qQC= qH->GetQC(); // PDGCode of theDecayingAntiStrangeNuclSt. G4int qS = qH->GetStrangeness(); // Strangness of the AntiStrangeNuclearState G4int qB = qH->GetBaryonNumber(); // BaryonNumber of the AntiStrangeNuclearSt. G4int qP = qH->GetCharge(); // Charge of the AntiStranNuclState (a#of p) #ifdef debug G4cout<<"G4QNuc::DecAntS:QC="<=0 || qB<1) { delete qH; G4cerr<<"G4QNuc::DecayAntiStrange:QC="<0 && qP>qN) // a#of protons > a#of neutrons { if(qP>=bH) // => "Enough protons in nucleus" case { if(qN>=sH) { n1=sH; n2=bH; } else { G4int dPN=qP-qN; if(dPN>=aS) { n1=0; n2=aS; } else { G4int sS=(aS-dPN)/2; G4int bS=aS-dPN-sS; sS+=dPN; if(qP>=sS && qN>=bS) { n1=bS; n2=sS; } else if(qP=bH) { if(qP>=sH) { n2=sH; n1=bH; } else { G4int dNP=qN-qP; if(dNP>=aS) { n1=aS; n2=0; } else { G4int sS=(aS-dNP)/2; G4int bS=aS-sS; if(qN>=bS && qP>=sS) { n1=bS; n2=sS; } else if(qNqM+eps && n1==1) // Try to use another K if this is the only kaon { G4int naPDG=90000000+(qP-1)*1000+qN; // Prot PDG of the Alternative Residual Nucleus G4double naM=G4QNucleus(naPDG).GetGSMass(); // Prot Mass of the Alt. Residual Nucleus G4int akPDG=321; // Prototype PDGCode of the AlternativeKaon (K+) G4double kaM=mK; // Prototype Mass of the AlternativeKaon (K+) if(k1PDG==321) // Calculate alternative to the K+ meson { naPDG=90000000+qP*1000+qN-1; // PDG of the Alternative Residual Nucleus naM=G4QNucleus(naPDG).GetGSMass(); // Mass of the Alt. Residual Nucleus akPDG=311; // PDG Code of the Alternative kaon (K0) kaM=mK0; // Mass of the Alternative kaon (K0) } G4double asum=naM+kaM; if(asumpush_back(qH); // @@ Can cause problems with particle conversion in G4 return; } #ifdef debug G4cout<<"G4QNuc::DecAntiS: nK+N "<push_back(H1); // Fill "H1" (delete equivalent) } G4QHadron* H2 = new G4QHadron(qPDG,s4Mom); // Create a Hadron for the Nucleus //evaHV->push_back(H2); // Fill "H2" (delete equivalent) EvaporateNucleus(H2,evaHV); // Fill "H2" (delete equivalent) #ifdef debug G4cout<<"G4QNucleus::DecAntiStr:*** After EvaporateNucleus nH="<size()<tM="<push_back(qH); // @@ Can cause problems with particle conversion in G4 return; } #ifdef debug G4cout<<"G4QNuc::DecAntiS:*DONE* nPDG="<push_back(H1); // Fill "H1" (delete equivalent) } s4Mom/=n2; for(G4int i2=0; i2push_back(H2); // Fill "H2" (delete equivalent) } G4QHadron* H3 = new G4QHadron(qPDG,t4Mom); // Create a Hadron for the nucleus //evaHV->push_back(H3); // Fill "H3" (delete equivalent) EvaporateNucleus(H3,evaHV); // Fill "H3" (delete equivalent) } #ifdef qdebug if (qH) { G4cout << "G4QNucleus::DecayAntiStrange: deleted at end - PDG: " << qH->GetPDGCode() << G4endl; delete qH; } #endif #ifdef debug G4cout<<"G4QNucleus::DecayAntiStrange: ===> End of DecayAntiStrangness"<Get4Momentum(); // Get 4-momentum of the MultyBaryon G4double qM = q4M.m(); // Mass of the Multybaryon G4int qPDG = qH->GetPDGCode(); // PDG Code of the decaying multybar G4QContent qQC = qH->GetQC(); // PDG Code of the decaying multibar #ifdef debug G4cout<<"G4QNuc::DecayMultyBaryon: *Called* PDG="<push_back(qH); else if(totBN==2) { G4LorentzVector f4Mom(0.,0.,0.,fMass); G4LorentzVector s4Mom(0.,0.,0.,fMass); G4double sum=fMass+fMass; if(fabs(qM-sum) TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecMulBar:*DONE* fPDG="<push_back(H1); // Fill "H1" (delete equivalent) G4QHadron* H2 = new G4QHadron(fPDG,s4Mom); // Create a Hadron for the 2-nd baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) } else if(totBN==3) { G4LorentzVector f4Mom(0.,0.,0.,fMass); G4LorentzVector s4Mom(0.,0.,0.,fMass); G4LorentzVector t4Mom(0.,0.,0.,fMass); G4double sum=fMass+fMass+fMass; if(fabs(qM-sum)? TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNuc::DecMBar:*DONE*, fPDG="<push_back(H1); // Fill "H1" (delete equivalent) G4QHadron* H2 = new G4QHadron(fPDG,s4Mom); // Create a Hadron for the 2-nd baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) G4QHadron* H3 = new G4QHadron(fPDG,t4Mom); // Create a Hadron for the 3-d baryon evaHV->push_back(H3); // Fill "H3" (delete equivalent) } else { // @@It must be checked, that they are not under the mass shell // !! OK !! Checked by the warning print that they are mostly in the Ground State !! G4LorentzVector f4Mom=q4M/totBN; // @@ Too simple solution (split in two parts!) #ifdef debug // Warning for the future development G4cout<<"*G4QNul::DecMulBar:SplitMultiBar inEqParts M="<GetPDGCode() << G4endl; delete qH; } #endif } // End of DecayMultyBaryon //Decay of the excited alpha+2p or alpha+2n systems void G4QNucleus::DecayAlphaDiN(G4QHadron* qH, G4QHadronVector* evaHV) {// ============================================ static const G4double mNeut= G4QPDGCode(2112).GetMass(); static const G4double mProt= G4QPDGCode(2212).GetMass(); static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); static const G4double mHel6= G4QPDGCode(2112).GetNuclMass(2,4,0); static const G4double eps=0.003; G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the AlphaDibaryon G4double qM = q4M.m(); // Real mass of the AlphaDibaryon G4int qPDG = qH->GetPDGCode(); // PDG Code of the decayin AlphaDybaryon #ifdef debug G4cout<<"G4QNuc::DecayAlphaDiN: *Called* PDG="<push_back(qH); // Fill as it is (delete equivalent) return; } else if(mNeut+mNeut+mAlph? TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNuc::DecAl2N: fPDG="<push_back(H1); // Fill "H1" (delete equivalent) G4QHadron* H2 = new G4QHadron(fPDG,s4Mom); // Create a Hadron for the 2-nd baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) G4QHadron* H3 = new G4QHadron(sPDG,t4Mom); // Create a Hadron for the alpha evaHV->push_back(H3); // Fill "H3" (delete equivalent) } // End of DecayAlphaDiN //Decay of the excited alpha+bayon state in alpha and baryons void G4QNucleus::DecayAlphaBar(G4QHadron* qH, G4QHadronVector* evaHV) {// ============================================ static const G4double mNeut= G4QPDGCode(2112).GetMass(); static const G4double mProt= G4QPDGCode(2212).GetMass(); static const G4double mLamb= G4QPDGCode(3122).GetMass(); static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0); static const G4double mHe3 = G4QPDGCode(2112).GetNuclMass(2,1,0); static const G4double eps=0.003; G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the Alpha-Baryon G4double qM = q4M.m(); // Mass of Alpha-Baryon G4int qPDG = qH->GetPDGCode(); // PDG Code of the decayin Alpha-Baryon G4QContent qQC = qH->GetQC(); // PDG Code of the decaying Alpha-Bar #ifdef debug G4cout<<"G4QNucleus::DecayAlphaBar: *Called* PDG="< 1) DecayMultyBaryon(qH,evaHV); else if(qPDG==92001002||qPDG==92002001||qPDG==91003001||qPDG==91001003||qPDG==93001001) evaHV->push_back(qH); else if(qPDG==92000003||qPDG==92003000||qPDG==93000002||qPDG==93002000) { G4int fPDG = 3122; // 1st Prototype for 2L+3n case G4double fMass= mLamb; G4int sPDG = 2112; G4double sMass= mNeut; if (qPDG==92003000) // "2L+3p" case { sPDG = 2212; sMass= mProt; } else if(qPDG==93000002) // "2n+3L" case { fPDG = 2112; fMass= mNeut; sPDG = 3122; sMass= mLamb; } else if(qPDG==93002000) // "2p+3L" case { fPDG = 2212; fMass= mProt; sPDG = 3122; sMass= mLamb; } G4double tfM=fMass+fMass; G4double tsM=sMass+sMass+sMass; G4LorentzVector f4Mom(0.,0.,0.,tfM); G4LorentzVector s4Mom(0.,0.,0.,tsM); G4double sum=tfM+tsM; if(fabs(qM-sum)M="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecAlB:*DONE*, fPDG="<push_back(H1); // Fill "H1" (delete equivalent) evaHV->push_back(H1); // Fill "H1" (delete equivalent) G4LorentzVector rs4Mom=s4Mom/3; G4QHadron* H2 = new G4QHadron(sPDG,rs4Mom); // Create a Hadron for the 2-nd baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) evaHV->push_back(H2); // Fill "H2" (delete equivalent) evaHV->push_back(H2); // Fill "H2" (delete equivalent) } else if(qPDG==90004001||qPDG==90001004) { G4int fPDG = 90002001; // Prototype for "He3+2p" case G4double fMass= mHe3; G4int sPDG = 2212; G4double sMass= mProt; if (qPDG==90001004) // "t+2n" case { fPDG = 90001002; fMass= mTrit; sPDG = 2112; sMass= mNeut; } G4LorentzVector f4Mom(0.,0.,0.,fMass); G4LorentzVector s4Mom(0.,0.,0.,sMass); G4LorentzVector t4Mom(0.,0.,0.,sMass); G4double sum=fMass+sMass+sMass; if(fabs(qM-sum) TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucl::DecAlB: *DONE*, f="<push_back(H2); // Fill "H2" (delete equivalent) G4QHadron* H3 = new G4QHadron(sPDG,t4Mom); // Create a Hadron for the 3-d baryon evaHV->push_back(H3); // Fill "H3" (delete equivalent) } else if(qPDG==94000001||qPDG==94001000||qPDG==91000004||qPDG==91004000) { G4int fPDG = 3122; // Prototype for "4L+n" case G4double fMass= mLamb+mLamb; G4int sPDG = 2112; G4double sMass= mNeut; if (qPDG==94001000) // "4L+p" case { sPDG = 2212; sMass= mProt; } else if(qPDG==91000004) // "4n+L" case { fPDG = 2112; fMass= mNeut+mNeut; sPDG = 3122; sMass= mLamb; } else if(qPDG==91004000) // "4p+L" case { fPDG = 2212; fMass= mProt+mProt; sPDG = 3122; sMass= mLamb; } G4LorentzVector f4Mom(0.,0.,0.,fMass+fMass); G4LorentzVector s4Mom(0.,0.,0.,sMass); G4double sum=fMass+fMass+sMass; if(fabs(qM-sum) TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNuc::DecAlphaB: *DONE*, fPDG="<push_back(H1); // Fill "H1" (delete equivalent) evaHV->push_back(H1); // Fill "H1" (delete equivalent) evaHV->push_back(H1); // Fill "H1" (delete equivalent) evaHV->push_back(H1); // Fill "H1" (delete equivalent) G4QHadron* H2 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the 2-nd baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) } else if(qPDG==90003002||qPDG==90002003||qPDG==91002002) { G4int fPDG = 90002002; // Prototype for "alpha+n" case G4int sPDG = 2112; G4double fMass= mAlph; G4double sMass= mNeut; if(qPDG==90003002) // "alpha+p" case { sPDG = 2212; sMass= mProt; } else if(qPDG==9100202) // "alpha+l" case { sPDG = 3122; sMass= mLamb; } else if(qPDG!=90002003) { evaHV->push_back(qH); // Fill hadron as it is (delete equivalent) //EvaporateNucleus(qH, evaHV); // Evaporate Nucleus (delete equivivalent) return; } G4double dM=fMass+sMass-qM; if(dM>0.&&dM<1.) { #ifdef debug G4cout<<"***Corrected***G4QNuc::DecayAlphaBar:fPDG="< TotM="< TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucl::DecAlBar:*DONE*a4M="<push_back(H1); // Fill "H1" (delete equivalent) G4QHadron* H2 = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the baryon evaHV->push_back(H2); // Fill "H2" (delete equivalent) } else G4cout<<"---Warning---G4QNucleus::DecayAlphaBar: Unknown PDG="<GetPDGCode() << G4endl; delete qH; } #endif } // End of DecayAlphaBar //Decay of the excited alpha+alpha state in 2 alphas void G4QNucleus::DecayAlphaAlpha(G4QHadron* qH, G4QHadronVector* evaHV) {// ============================================== static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); static const G4double aaGSM= G4QPDGCode(2112).GetNuclMass(4,4,0); static const G4double eps=0.003; G4int qPDG = qH->GetPDGCode(); // PDG Code of the decayin dialpha if(qPDG!=90004004) { delete qH; G4cerr<<"***G4QNucleus::DecayAlphaAlpha: qPDG="<Get4Momentum(); // Get 4-momentum of the Dibaryon G4double qM=q4M.m(); #ifdef debug G4cout<<"G4QNucleus::DecayAlAl: *Called* PDG="<"<aaGSM+.01) // @@ Be8*->gamma+Be8 (as in evaporation) @@ gamma cooling if(2>3) { G4int fPDG = 22; G4int sPDG = 90004004; G4double fMass= 0.; G4double sMass= aaGSM; G4LorentzVector f4Mom(0.,0.,0.,fMass); G4LorentzVector s4Mom(0.,0.,0.,sMass); // Mass is random since probab. time G4double sum=fMass+sMass; if(fabs(qM-sum) TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayAlphaAlpha: *DONE* gam4M="<Set4Momentum(s4Mom); q4M=s4Mom; } G4int fPDG = 90002002; G4int sPDG = 90002002; G4double fMass= mAlph; G4LorentzVector f4Mom(0.,0.,0.,fMass); G4LorentzVector s4Mom(0.,0.,0.,fMass); G4double sum=fMass+fMass; if(fabs(qM-sum) TotM="<push_back(qH); return; } #ifdef debug G4cout<<"G4QNucleus::DecayAlphaAlpha: *DONE* fal4M="<push_back(H2); // Fill "H2" (delete equivalent) } // End of DecayAlphaAlpha