// // ******************************************************************** // * 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. * // ******************************************************************** // // 1 2 3 4 5 6 7 8 9 //34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901 // // // $Id: G4Quasmon.cc,v 1.125 2010/06/10 08:37:27 mkossov Exp $ // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ // // ---------------- G4Quasmon ---------------- // by Mikhail Kossov, July 1999. // class for an excited hadronic state used by the CHIPS Model // --------------------------------------------------------------------------- // Short description: The Quasmon is the main object of the CHIPS model, where // hadronic states are defined by the quark content and the 4-momentum (mass). // In this sense a Quasmon is a generalised hadron - hadrons are the low // mass discrete states of quasmons. Here the hadron can be not only an // elementary particle, consisting of 2 or 3 quarks, but a nucleonic // cluster, consisting of 6, 9, ... , 3n, ... quarks. In the CHIPS model // the nucleus is considered as a group of the nucleons and the nucleonic // clusters. In hadronic reactions a Quasmon is constructed in vacuum as // a result of the collision (G4QCollision process). In this case only // the G4Quasmon class can be used for the reaction. In nuclear reactions // one or a few Quasmons can be created as a result of the colision of // the projectile hadrons with the nucleons and the nucleonic clusters // of the nuclear target. In this case the Quasmons are created and // fragmented in the nuclear environment (G4QNucleus) and the G4QEnvironment // class must be used for the reaction. For nuclear-nuclear reactions // two G4QEnvironments are used with the common (which can fragment in both // nuclear environments - mostly at low energies), individual (which can // fragment in only one of two G4QEnvironments) and vacuum (which can // fragment in vacuum being too far by rapidity from both nuclei) Quasmons. // -------------------------------------------------------------------------- // //#define debug //#define pdebug //#define pardeb //#define psdebug //#define rdebug //#define ppdebug //#define chdebug //#define tdebug //#define sdebug #include "G4Quasmon.hh" #include #include using namespace std; G4Quasmon::G4Quasmon(G4QContent qQCont, G4LorentzVector q4M, G4LorentzVector ph4M): q4Mom(q4M), valQ(qQCont), phot4M(ph4M), f2all(0), rEP(0.), rMo(0.) { #ifdef debug G4cout<<"G4Quasmon:Constructor:QC="<G4Q:Con:(1),T="<0.) q4Mom+=phot4M; // InCaseOf CaptureByQuark it will be subtracted back valQ.DecQAQ(-1); status=4; } G4Quasmon::G4Quasmon(const G4Quasmon& right) { q4Mom = right.q4Mom; valQ = right.valQ; //theEnvironment = right.theEnvironment; status = right.status; //theQHadrons (Vector) G4int nQH = right.theQHadrons.size(); if(nQH) for(G4int ih=0; ih>>DONE<<<"<>>START QUASMON HADRONIZATION<<<=*, aP="< s) cQ.DecQAQ((tQ-s)/2); // Reduce QC to minimum QC #ifdef debug G4int rsPDG = cQ.GetSPDGCode(); // PDG for the lowest residual Quasmon state G4cout<<"G4Q::HQ:eN="< diPiM) qCountMax=(G4int)(excE/mPi0); // Try more for big excess // //G4int kCountMax=27; //G4int kCountMax=9; //G4int kCountMax=3; // Try different k if they are below minK G4int kCountMax=1; // "No reson to increase it" //G4int kCountMax=qCountMax; // "No reson to increase it" //G4int kCountMax=0; //@@ *** Close search for the minimum k *** // //G4int pCountMax=27; //Try differentHadrons(Parents) forBetterRecoil //G4int pCountMax=9; //Try differentHadrons(Parents) forBetterRecoil //G4int pCountMax=3; //Try differentHadrons(Parents) forBetterRecoil G4int pCountMax=1; //Try differentHadrons(Parents) forBetterRecoil //if(envA>0) pCountMax=3; if(envA>0&&envA<19) pCountMax=36/envA; //if(envA>0&&envA<31) pCountMax=60/envA; //if(envA>0&&envA<61) pCountMax=120/envA; G4bool gintFlag=false; // Flag of gamma interaction with one quark while(kCountmPi0) miM2=mPi02; else miM2=mP2; } else // "Env. exists" case - find k_min & k_max { minK=100000.; // @@ ??? May be for the light nuclei ??? //if(piF&&quasM>iniQM) maxK=mNeut/2.+(quasM-iniQM)*iniQM/(iniQM+mNeut); //else if(quasM>iniQM) //{ // G4double limK=envM/2.+(quasM-iniQM)*iniQM/(iniQM+envM); // if(limK0&&iniBN>1) // "Neutron-hadron" estimate //{ // G4double iK=1000000.; // G4double dqm=quasM+quasM; // if(envA>0) // { // G4QContent rtQC=valQ-neutQC+envQC; // Total Residual Quark Content // G4QNucleus rtN(rtQC); // Create pseudo-nucleus for the TotalResidual // G4double rtM =rtN.GetMZNS(); // Min Mass of total residual Nucleus // G4double bnRQ=rtM-envM; // Bound mass of residual Quasmon // G4double sm2=qM2+m2N-bnRQ*bnRQ; // G4double fqm=dqm+dqm; // G4double aK=sm2/fqm; #ifdef debug // G4double kts=.135; // Test value of k // G4double dkts=kts+kts; // G4double fu=dkts*(dkts*quasM-sm2)+quasM*m2N; // G4cout<<"G4Q::HQ:M="<=mNeut+bnRQ) // { // G4double srm=sqrt(sm2*sm2-4.*qM2*m2N)/fqm; // iK=aK-srm; // aK+=srm; // } // else iK=aK; // if(aK0&&totBN>1) // "Neutron-cluster" estimate { G4QNucleus tmpNN(totZ,totN-1,0); G4double delN=tmpNN.GetMZNS()+mNeut-totM; if(envN>0&&envA>1) { G4QNucleus envNN(envZ,envN-1,0); G4double delEN=envNN.GetMZNS()+mNeut-envM; if(delEN>delN) delN=delEN; } delN*=qurF; if(delN0&&totBN>1) // "Proton-cluster" estimate { G4double proCB=theEnvironment.CoulombBarrier(1,1); G4QNucleus tmpPN(totZ-1,totN,0); G4double delP=tmpPN.GetMZNS()+mProt-totM+proCB; if(envZ>0&&envA>1) { G4QNucleus envPN(envZ-1,envN,0); G4double delEP=envPN.GetMZNS()+mProt-envM+proCB; if(delEP>delP) delP=delEP; } delP*=qurF; if(delP0&&totZ>0&&totBN>2) // "Deuteron-cluster" estimate { G4double proCB=theEnvironment.CoulombBarrier(1,2); G4QNucleus tmpDN(totZ-1,totN-1,0); G4double delD=tmpDN.GetMZNS()+mDeut-totM+proCB; if(envN>0&&envZ>0&&envA>2) { G4QNucleus envDN(envZ-1,envN-1,0); G4double delED=envDN.GetMZNS()+mDeut-envM+proCB; if(delED>delD) delD=delED; } delD*=qurF; if(delD1&&totZ>0&&totBN>3) // "Triton-cluster" estimate { G4double proCB=theEnvironment.CoulombBarrier(1,3); G4QNucleus tmpTN(totZ-1,totN-2,0); G4double delT=tmpTN.GetMZNS()+mTrit-totM+proCB; if(envN>1&&envZ>0&&envA>3) { G4QNucleus envTN(envZ-1,envN-2,0); G4double delET=envTN.GetMZNS()+mTrit-envM+proCB; if(delET>delT) delT=delET; } delT*=qurF; if(delT0&&totZ>1&&totBN>3) // "He3-cluster" estimate { G4double proCB=theEnvironment.CoulombBarrier(2,3); G4QNucleus tmpRN(totZ-2,totN-1,0); G4double delR=tmpRN.GetMZNS()+mHel3-totM+proCB; if(envN>0&&envZ>1&&envA>3) { G4QNucleus envRN(envZ-2,envN-1,0); G4double delER=envRN.GetMZNS()+mHel3-envM+proCB; if(delER>delR) delR=delER; } delR*=qurF; if(delR1&&totZ>1&&totBN>4) // "Alpha-cluster" estimate { G4double proCB=theEnvironment.CoulombBarrier(2,4); G4QNucleus tmpAN(totZ-2,totN-2,0); G4double delA=tmpAN.GetMZNS()+mAlph-totM+proCB; if(envN>1&&envZ>1&&envA>4) { G4QNucleus envAN(envZ-2,envN-2,0); G4double delEA=envAN.GetMZNS()+mAlph-envM+proCB; if(delEA>delA) delA=delEA; } delA*=qurF; if(delA*qurF2&&totZ>2&&totBN>6) // "Li6-cluster" estimate ?????? //{ // G4double proCB=theEnvironment.CoulombBarrier(3,6); // G4QNucleus tmpLN(totZ-3,totN-3,0); // G4double delL=tmpLN.GetMZNS()+mLit6-totM+proCB; // if(envN>2&&envZ>2&&envA>6) // { // G4QNucleus envLN(envZ-3,envN-3,0); // G4double delEL=envLN.GetMZNS()+mLit6-envM+proCB; // if(delEL>delL) delL=delEL; // } // delL*=qurF; // if(delLquasM) minK=0.; } //if(addPhoton>0.&&quasM<1500.&&G4UniformRand()0.&&iniBN<2)// PhotonAbsorbDiagramContrib(@@HiddenPar) //if(addPhoton>0.&&(G4UniformRand()<.7||iniBN<2))//PhotoAbsorbDiagramContrib //if(addPhoton>0.&&(G4UniformRand()<.5||iniBN<2))//PhotonAbsorbDiagramContrib //if(addPhoton>0.&&(G4UniformRand()<.27||iniBN<2))// PhotonAbsorbDiagramContrib //if(addPhoton>0.&&iniBN<2)// PhotonAbsorbDiagramContrib(@@HiddenPar) //if(addPhoton>0.&&iniBN>1)// PhotonAbsorbDiagramContrib(@@HiddenPar) //if(addPhoton>0.&&G4UniformRand()<0.5)// PhotonAbsorbDiagramContrib(@@HiddenPar) //if(addPhoton>0.&&G4UniformRand()<0.75)// PhotonAbsorbDiagramContrib(@@HiddenPar) ///if(addPhoton>0.&&G4UniformRand()<0.8)// PhotonAbsorbDiagramContrib(@@HiddenPar) //if(addPhoton>0.&&iniBN>1)//PhotonAbsorbDiagramContrib(Photon's captured by N quark) if(addPhoton>0.)// PhotonAbsorbDiagramContrib(Photon is always captured by quarks) //if(2>3) // Photon capture by quark is closed { //nOfQ=valQ.GetQ()-valQ.GetAQ(); // Recover nOfQ for the not excited cluster // @@ 1/(p_g,p_q) interaction probability (? (p_g,p_q)=0 always) gintFlag=true; q4Mom-= phot4M; // recover Quasmon-Cluster without the Photon qM2 = q4Mom.m2(); // Update the Current squared mass of Quasmon quasM = sqrt(qM2); // Update the Current Mass of Quasmon G4double kpow=static_cast(nOfQ-2); // n-3+1 because of integration if(kpow<1.) kpow=1.; G4double xmi=(momPhoton-addPhoton)*quasM; // Minimum residual mass 2kM if(xmi<0.) xmi=0.; // While must be commented from here and down __________________________ //G4bool frat=true;//[k/(k+p) factor of QuarkExch]*[p/k factor of fixed ct]=p/(p+k) //G4int cMax=27;//For m_gam=0:*[k/(k+p) factor 1/s{m<kMom/(addPhoton+kMom); // //frat=G4UniformRand()>kMom*kMom/(addPhoton+kMom)/(addPhoton+kMom); // //frat=G4UniformRand()>addPhoton/(addPhoton+kMom); // //frat=G4UniformRand()>addPhoton*addPhoton/(addPhoton+kMom)/(addPhoton+kMom); // pCount++; //} if(cost>1.) cost=1.; if(cost<-1.) cost=-1.; #ifdef debug G4cout<<"G4Q::HQ:**PHOTON out of Q**k="<1.||r2<.0001) // | //{ // x = G4UniformRand(); // | // y = G4UniformRand(); // | // z = G4UniformRand(); // | // r2=x*x+y*y+z*z; // | //} //G4double r=140./sqrt(r2); // | //G4double xs=x*r; // | //G4double ys=y*r+kMom*sqrt(1.-cost*cost); // | //G4double zs=z*r+kMom*cost; // | //kMom=sqrt(xs*xs+ys*ys+zs*zs); // | //cost=zs/kMom; // | //if(kMom>=quasM/2.) kMom=quasM/2.-.001; // | // --- End of the pseudo Fermi-motion corection-----------------------------------^ // --- Virtual gluoCompton correction starts here --------------------------------* //G4double hms=32400.; // T_c^2 | ////G4double hms=16200.; // T_c^2/2 (OK) | ////G4double hms=8100.; // T_c^2/4 | //G4double x=kMom*kMom/hms; // | //G4double dc=(pow(x,G4UniformRand())-1.)/x; // smearing delta-function absorbtion| //if(dc<0.)dc=0.; // only positive smearing | //cost-=dc; // for virtPhotons smear in both dir.| // --- Quark mass correction ends here (?) ---------------------------------------* } else { gaF=false; // GammaFirstAct flag is only for the gamma-q gintFlag=false; // ==== Probabiliti proportional to k (without 1/k factor of hadronization) if(!miM2) miM2=(minK+minK)*quasM; // Make minimum mass for randomization if(qM2<.0001) kMom=0.; else kMom = GetQPartonMomentum(maxK,miM2); // Calculate value of primary qParton // ==== Direct calculation of the quark spectrum //G4double kpow=static_cast(nOfQ-2); //G4double kst=0.; //if(maxK+maxKminK) //{ // G4double rn=(pow((1.-(minK+minK)/quasM),kpow)-kst)*G4UniformRand(); // kMom=(1.-pow(kst+rn,1./kpow))*quasM/2.; //} //else kMom=(minK+maxK); //^^^ Direct calculation G4double rnc=G4UniformRand(); cost = 1.-rnc-rnc; } G4double cQM2=qM2-(kMom+kMom)*quasM; if(cQM2<0.) { //if(cQM2<-.0001)G4cerr<<"--Warning--G4Q::HQ:(PhBack) cQM2="<ColResQ+k_part { #ifdef debug G4cerr<<"*G4Q::HQ:PB,M="<0.) quasM = sqrt(qM2); // Update the Current Mass of Quasmon else { if(qM2<-.0001)G4cerr<<"--Warning-- G4Q::HQ:Phot.M2="< The decay succeeded { if(addPhoton&&gintFlag) // Make it as if the phaton was a part of the Q { q4Mom+= phot4M; // Recover the Full Quasmon with the Photon k4Mom+= phot4M; // add photon energy to a quark of the Cluster qM2 = q4Mom.m2(); // Update the Current squared mass of Quasmon quasM = sqrt(qM2); // Update the Current Mass of Quasmon kMom=k4Mom.e(); // Update the k_parton momentum gintFlag=false; } #ifdef debug G4cout<<"G4Q::HQ:(PhBack) k="<>>K-ITERATION #"<kMin=" <kMin) { kCond=false; break; } } } } } kCount++; } // End of search for "k": *** Here we forget about the initial photon forever ***???*** #ifdef debug if(addPhoton) G4cout<<"G4Q::HQ:***PHOTON OK***k="<GetIntegProbability(); #ifdef debug G4cout<<"G4Q::HQ:***RANDOMIZE CANDIDATEs***a#OfCand="<mPi0+iniQM) { gamPDG=111; gamM=mPi0; } G4LorentzVector r4Mom(0.,0.,0.,gamM); // mass of the photon/Pi0 G4LorentzVector s4Mom(0.,0.,0.,iniQM); // mass of the GSQuasmon G4double sum=gamM+iniQM; if(sum>0. && fabs(quasM-sum)g/pi0(M="<HadrVec, Q="<g/pi0("<minK="<New p-Attempt#"<GetIntegProbability() < totP) i++; if (i>=nCandid) { G4cerr<<"***G4Q::HQ: Cand#"<= Tot#"<1&&totMass>totM&&totS>=0) { #ifdef pdebug G4cout<<"G4Quas::HQ:NEED-EVAP-1:Q="< "Chipolino" case { fchipo=true; G4QChipolino chipQ(valQ); G4QPDGCode QPDG1=chipQ.GetQPDG1(); sPDG = QPDG1.GetPDGCode(); sMass= QPDG1.GetMass(); G4QPDGCode QPDG2=chipQ.GetQPDG2(); rPDG = QPDG2.GetPDGCode(); rMass= QPDG2.GetMass(); if(sMass+rMass>quasM) //@@ Evaporate but try 3Pt decay with Environ { if(totBN>1&&totMass>totM&&totS>=0) { #ifdef pdebug G4cout<<"G4Q::HQ:NEED-EVAP-2:Q="<minM="<1) { #ifdef pdebug G4cout<<"G4Quas::HQ:NEED-EVAP-0:Q="<@@ Final decay in MinHadr+(PI0 or GAM)?? { sMass=G4QPDGCode(sPDG).GetMass(); rPDG=envPDG; if (rPDG>MINPDG&&rPDG!=NUCPDG) { rMass=theEnvironment.GetMZNS(); q4Mom+=G4LorentzVector(0.,0.,0.,rMass); valQ +=theEnvironment.GetQC(); quasM=q4Mom.m(); KillEnvironment(); fred=true; // Flag of environment reduction } else if(rPDG!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0) // ===> "Evaporation" case { #ifdef pdebug G4QContent nTotQC=totQC-neutQC; G4QNucleus nTotN(nTotQC); // PseudoNucleus for TotalResidual to neutron G4double nTotM =nTotN.GetMZNS(); // MinMass of the Total Residual to neutron G4double dMnT=totMass-nTotM-mNeut;// Energy excess for the neutron G4QContent pTotQC=totQC-protQC; G4QNucleus pTotN(pTotQC); // PseudoNucleus for TotalResidual to proton G4double pTotM =pTotN.GetMZNS(); // MinMass of the Total Residual to proton G4double dMpT=totMass-pTotM-mProt;// Energy excess for neutron if(dMpT>dMnT)dMnT=dMpT; G4cerr<<"G4Q::HQ:NEED-EVAP3:s="< "Something was selected" case { G4QCandidate* curCand = theQCandidates[i];// Pointer toSelectedCandidate:hadr/fragm sPDG = curCand->GetPDGCode(); // PDG of the selected candidate //////////////////G4double prpr=curCand->GetPreProbability(); #ifdef debug G4cout<<"G4Q::HQ:hsfl="<GetPClustEntries(); G4double sppm = curCand->TakeParClust(nParCandid-1)->GetProbability(); if (sppm<=0) // Impossible to find a parent cluster { G4cerr<<"***G4Quasmon::HadronizeQ:P="<GetIntegProbability() <<",nC="<TakeParClust(ip)->GetProbability(); G4cerr< "Parent cluster was found" case { G4double totPP = sppm * G4UniformRand(); while(curCand->TakeParClust(ip)->GetProbability() < totPP) ip++; #ifdef debug G4cout<<"G4Q::HQ:p#ip="<TakeParClust(ip); pPDG = parCluster->GetPDGCode(); G4QPDGCode pQPDG(pPDG); pQC = pQPDG.GetQuarkContent(); pBaryn= pQC.GetBaryonNumber(); pMass = parCluster->GetEBMass(); // Environment Bounded Mass pNMass = parCluster->GetNBMass(); // NuclBoundedMass @@Env.Reduce dosn't work transQC = parCluster->GetTransQC(); delta = parCluster->GetEBind(); // Environmental Binding deltaN = parCluster->GetNBind(); // Nuclear Binding loli = parCluster->GetLow(); hili = parCluster->GetHigh(); //G4double dhil=.0001*(hili-loli); // Safety factor //loli += dhil; //hili -= dhil; npqp2 = parCluster->GetNQPart2(); // @@ One can get otherUsefulParameters of the parent cluster for hadronization G4QPDGCode sQPDG(curCand->GetPDGCode()); sQC = sQPDG.GetQuarkContent(); //if(sPDG==90001001 && G4UniformRand()>0.75) sMass=np; //@@ n-p pair if(sPDG==90001001 && G4UniformRand()>1.0) sMass=np; //@@ no n-p pair (close) else sMass = sQPDG.GetMass(); sM2 = sMass*sMass; // Squared mass of the fragment curQ += transQC; // Subtract ExchangeMesonQC from QuasmonQC #ifdef debug G4cout<<"G4Q::HQ:valQ="< "Hadron" case { pBaryn=0; // @@ ? sQC=theQCandidates[i]->GetQC(); sMass = theQCandidates[i]->GetNBMass();// Mass is randomized on probability level sM2=sMass*sMass; // SqMass is randomized on probability level curQ-= sQC; // Subtract outHadron QC from QC of Quasmon #ifdef debug G4cout<<"G4Q::HQ: hsfl="< Calculate the SquaredGroundStateMass of theResidualQuasmon+ResidualEnvironment //if(envPDG>pPDG) if(envA>=pBaryn) { // *** LIM *** G4QContent RNQC=curQ+envQC; if(sPDG>MINPDG&&sPDG!=NUCPDG) RNQC-=pQC; // ==> "Nuclear Fragment Radiation" case if(envA-pBaryn>bEn) RNQC=curQ+bEnQC; // Leave the minimum environment G4int RNPDG = RNQC.GetSPDGCode(); // Total kinematically involved nuclear mass if(RNPDG==10) minSqN=G4QChipolino(RNQC).GetMass2();// MinSqM of DiHadron of Chipo else if(!RNPDG) // It never should happen as curQ is real and env/bEn QC is > 0 { //#ifdef debug G4cout<<"**G4Q::HQ:*KinematicTotal*, PDG="< MINPDG && sPDG != NUCPDG) {// ==> NuclearFragmentCandidate hadronization #ifdef debug G4cout<<"G4Q::HQ: BoundM="< colouredCluster (cc) G4LorentzVector cl4Mom(0.,0.,0.,pMass);// 4-momentum prototype for parent cluster G4LorentzVector tot4Mom=q4Mom+cl4Mom; // @@ Just for checking #ifdef debug G4cout<<"G4Q::HQ:Q("<k("< frM2="< sM="< rTM="<bEn) tmpEQ=bEnQC; // Leave the minimum environment G4QNucleus tmpN(tmpEQ); // Pseudo nucleus for Residual Environment G4double tmpNM=tmpN.GetMZNS(); // Mass of Residual Environment #ifdef debug G4cout<<"G4Q::HQ:eQC="<1.0001)G4cerr<<"Warn-G4Q::HQ: cost_max="<1.CorHadrProb"<cta+.0001)G4cerr<<"**G4Q::HQ:ci="<ca="< 1.)ctkk= 1.; if(ctkk<-1.)ctkk=-1.; } G4double cen=kLS+pMass; // LS Energy of k+parentCluster CompSystem G4double ctc=(cen*ctkk-kLS)/(cen-kLS*ctkk);//cos(theta_k,kap) in k+pClastSyst if(abs(ctc)>1.) { //G4cout<<"***G4Quasm:HadrQ: e="<ColDec>>>CF("<F("<=minSqB && fr4Mom.e()>=sCBE) // minSqB = SqGSMass of BoundResidQ // ***VFQ*** //G4double rQM2=rQ4Mom.m2(); // TMP (Before the "kt" is defined) //if(rQM2>=minSqT && fr4Mom.e()>=sCBE) // minSqB = SqGSMass of BoundResidQ // ***VQU*** fQ4Mom=rQ4Mom+G4LorentzVector(0.,0.,0.,tmpBE); // Free Quasmon 4-mom G4double fQM2=fQ4Mom.m2(); // TMP (Before the "kt" is defined) if(fQM2>=minSqT && reTNM2>=tmpTM2 && fr4Mom.e()>=sCBE) // minSqT = SqGSMass { qCond=false; // Ok, the appropriate q is found //ffdc=false; #ifdef debug // ***VTN*** //G4cout<<"G4Q::HQ:Attemp#"<"< CB+M="<"<"<MEMrQM2 && fr4Mom.e()>=sCBE)//tM>minTM,maxRQM,CB // ***VBQ***VFQ*** //if(reTNM2MEMrQM2 && fr4Mom.e()>=sCBE)//tM>minTM,maxRQM,CB // ***VTN*** //if(reTNM2MEMrQM2 && fr4Mom.e()>=sCBE) { // ***VQU*** MEMrQM2=fQM2; // ***VBQ***VFQ*** //MEMrQM2=rQM2; // ***VTN*** //MEMrQM2=reTNM2; //------------ MEMkp4M=kp4Mom; MEMfr4M=fr4Mom; MEMrQ4M=rQ4Mom; MEMreM2=reTNM2; } // ***VQU*** else if(fr4Mom.e()MEMsCBE&&reTNM2>=tmpTM2&&fQM2>MEMrQM2) // ***VBQ***VFQ*** //else if(fr4Mom.e()MEMsCBE&&reTNM2>tmpTM2&&rQM2>MEMrQM2) // ***VTN*** //else if(fr4Mom.e()MEMsCBE && reTNM2>=tmpTM2) { MEMsCBE=fr4Mom.e(); // Remember the best choice MEMkp4M=kp4Mom; MEMfr4M=fr4Mom; MEMrQ4M=rQ4Mom; MEMreM2=reTNM2; } else if(!qCount) //@@ Should not be here { // ***VQU*** MEMrQM2=fQM2; // ***VBQ***VFQ*** //MEMrQM2=rQM2; // ***VTN*** //MEMrQM2=reTNM2; //------------ MEMsCBE=fr4Mom.e(); MEMkp4M=kp4Mom; MEMfr4M=fr4Mom; MEMrQ4M=rQ4Mom; MEMreM2=reTNM2; } else { // ***VQU*** fQM2=MEMrQM2; // ***VBQ***VFQ*** //rQM2=MEMrQM2; // ***VTN*** //reTNM2=MEMrQM2; //----------- kp4Mom=MEMkp4M; // Make the best choice actual instead of the last fr4Mom=MEMfr4M; rQ4Mom=MEMrQ4M; reTNM2=MEMreM2; } } qCount++; } // End of the WHILE of the q-choice for the fixed parent // If q-choice is exhosted, then get the best, not the last quexf=true; // Quark Exchange is successfully done #ifdef debug G4cout<<"G4Q::HadQ:RQ("<bEn) totC4Mom+=G4LorentzVector(0.,0.,0.,mbEn); else totC4Mom+=G4LorentzVector(0.,0.,0.,envM-pMass); kn=totC4Mom.m2(); #ifdef debug G4cout<<"G4Q::HQ:A="<"<"<"<"< "HadronicCandidate hadronization" case { kt = (quasM-dk)*(quasM-sM2/dk); // squared mass of the RecoilQuasmon G4double rQM=0.; if(kt>0.) rQM=sqrt(kt); // Mass of the residual quasmon fr4Mom=G4LorentzVector(0.,0.,0.,sMass); // 4-mom update for the fragment rQ4Mom=G4LorentzVector(0.,0.,0.,rQM); // 4-mom update for the RecoilQuasmon if(!G4QHadron(q4Mom).DecayIn2(fr4Mom, rQ4Mom)) { G4cerr<<"*G4Quasm::HadrQ: q4M="<MINPDG&&envPDG!=NUCPDG) { // *** LIM *** G4LorentzVector TCRN=rQ4Mom; if(envA>bEn) TCRN+=bEn4M; else TCRN+=env4M; kn=TCRN.m2(); // tot4M - fr4Mom } else kn=kt; } // *** LIM *** G4LorentzVector tL=rQ4Mom; // @@ Is rQ4Mom calculated for hadrons?? tL+=G4LorentzVector(0.,0.,0.,reM); G4double tM=tL.m(); // Real Residual Total Nucleus Mass (hadr/frag) // Residual S+S limit (absoluteLowLimit for corresponding P-res.) for R->S+S Decay #ifdef debug G4cout<<"G4Q::HQ:k="<"<"<rtM="< T="<"<>>rM="<.5) rPDG=221; //if(rPDG==221&&sPDG!=221&&sPDG!=331&&rrn<.5) rPDG=111; G4int aPDG = abs(rPDG); G4int rb = abs(curQ.GetBaryonNumber()); // BaryNum of residual hadronic state G4double rcMass=-BIG; //@@ just BIG number // Prototype of minimalMass for residual if (!rPDG) { G4cerr<<"***G4Quasmon::HadronizeQuasmon: rQ ="<rMass)rFl=true;// @@ Kills Resonance+Resonance //if(sPDGH("<rM="< HadronVec, Q="< s4M="<minSqN && kt>minSqT) // ***VQU*** if (kn > minSqN && ku > minSqT) // ***VBQ*** //if(kn>minSqN && kt>minSqB) // ***VTN*** //if(kt>minSqB&&sPDGMINPDG&&kn>minSqN) { pCond=false; // Ok, the appropriate parent cluster is found #ifdef debug // ***VTN***VBQ*** //G4cout<<"G4Q::HQ:P-Attempt#"<" // <"<" <"<" // <"<=minSqB || sPDG>MINPDG&&kn>minSqN) // ***VBQ*** //if(kn=minSqB) // ***VFQ*** //if(kn=minSqT) // ***VQU*** if (kn < minSqN && ku < minSqT) { // ***VTN*** (former default) //if(ktPMEMktM2 || // knMINPDG && kn>PMEMknM2) // ***VBQ*** //if(ktPMEMktM2) // ***VFQ*** //if(ktPMEMktM2) // ***VQU*** if(ku < minSqT && ku > PMEMktM2) { // ***VQU*** PMEMktM2=ku; // ***VFQ***VBQ***VTN*** //PMEMktM2=kt; // --------- PMEMknM2=kn; PMEMfr4M=fr4Mom; PMEMrQ4M=rQ4Mom; PMEMreM2=reTNM2; PMEMrMas=rMass; PMEMpMas=pMass; PMEMsMas=sMass; PMEMdMas=dMass; PMEMmiSN=minSqN; PMEMmiST=minSqT; PMEMmiSB=minSqB; PMEMrPDG=rPDG; PMEMsPDG=sPDG; PMEMpPDG=pPDG; PMEMpQC =pQC; PMEMsQC =sQC; PMEMtQC =transQC; PMEMcQC =curQ; PMEMhsfl=hsflag; PMEMnucf=nucflag; #ifdef debug G4cout<<"G4Q::HQ:RemTheBest rPDG="<rPDG="<.5) rPDG=221; if(rPDG==221&&sPDG!=221&&sPDG!=331&&rrr<.5) rPDG=111; } //G4double reMass=sqrt(minSqT); // Min ResidQuasmon Mass after decay G4double reMass=sqrt(minSqB); // Min ResidQuasmon Mass after decay if (!rPDG) { G4cerr<<"***G4Q::HQ:Q="<MINPDG&&(sPDG MINPDG && envPDG != NUCPDG) || (sPDG > MINPDG && sPDG!=NUCPDG && envPDG > pPDG) ) && iniBN > 0 ) || iniBN > 1 || rPDG == 10 ) aMass=0.; // No Pi0 cond.(eg in NucE) #ifdef debug G4cout <<"G4Q::HQ:Is hsfl="<pPDG="< MINPDG && envPDG > pPDG && reTNM2 < tmpTM2) || fdul ) { // >>>> Decay Q->S+H or Q/C->H1+H2 or suck in or evaporate or slow down or decay etc. // ========> Decide what to do, if fragmentation in this Candidate is impossible === #ifdef debug G4cout<<"G4Q::HQ: Yes(No), hsf="<retNM) tmpT=G4QNucleus(tmpTQ,retN4Mom); G4QPDGCode sQPDG(sPDG); // tmpNM - residualenvironmentm (M),retNM - mass of TotalNucl [MN=sqrt((E+M)^2-p^2] // tmpRM - ResidQuasmongsm (m_GS), rQ4Mom - 4-momentum of ResidQuasmon (E,p,m) ////G4double m2=rQ4Mom.m2(); // Real Squared Mass of the ResidualQuasmon // potE=[sqrt(E^2*M^2-m^2*M^2+m_GS^2*MN*2)-m^2-E*M]/MN=-U (bindEn should be cutOff) //G4double pEc=2*(tmpRM+tmpNM-retNM); // DoubledBindingEnergy (virial theorem) /////////G4double pEc=tmpRM+tmpNM-retNM; // BindingEnergy (relativistic effect) //G4int rB=rQPDG.GetBaryNum(); // Baryon number of residQ //G4double rCB=theEnvironment.CoulombBarrier(rQPDG.GetCharge(),rB);// CB for residQ rMass=rQPDG.GetMass(); //G4int sB=sQPDG.GetBaryNum(); // Baryon number of residQ //G4double sCB=theEnvironment.CoulombBarrier(sQPDG.GetCharge(),sB);// CB of Fragm. #ifdef debug //G4int bSplit=tmpT.SplitBaryon(); // Possibility to split baryon from TotResN //G4double EQ=rQ4Mom.e(); // EnergyOfResidualQuasmon (E) //G4double em=tmpNM*EQ; // ResEnvM * EnergyOfResidualQuasmon (M*E) //G4double mM=retNM*tmpRM; // TotResNuclM*GSMassResidQuasmon (m_GS*MN) //G4double pEn=(sqrt(em*em-m2*tmpNM*tmpNM+mM*mM)-m2-em)/retNM; //Real BindingEnergy //G4double pEt=tmpNM*(EQ+tmpNM-retNM)/retNM; // Energy Transfer to nucleus //G4double PQ=rQ4Mom.rho(); // mod3MomentumOfResidualQuasmon //G4double pPt=tmpNM*PQ/retNM; // mod3Momentum Transfer to nucleus //G4cout<<"G4Q::HQ:tM="<RE="<3) //*** Attempt GoOutDropExcNucl is closed *** { G4LorentzVector fr4M = G4LorentzVector(0.,0.,0.,sMass);//GSM of Fragment G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,tmpNM);//GSM of ResidualEnviron G4LorentzVector rq4M = G4LorentzVector(0.,0.,0.,rMass);//GSM of ResidualQuasmon #ifdef debug G4double cfM=fr4Mom.m(); // @@ ? G4double ctM=tot4M.m(); // @@ ? G4cout<<"G4Q::HQ: *YES*,tM="< GSM="<tM="<totalMass"); #endif #ifdef debug G4cout<<"G4Q::HQ:Q="<mLamb+mPi0) { nucM = mLamb; nucPDG= 3122; piM = mPi0; piPDG = 111; } else if(abs(totZ)==1&&totMass>mLamb+mPi) { nucM = mLamb; nucPDG= 3122; piM = mPi; if(totZ>0) piPDG = 211; else piPDG =-211; } else { G4cerr<<"***G4Q::HQ:Z="<mNeut+mK0) { nucM = mNeut; nucPDG= 2112; piM = mK0; piPDG = 311; } else if(totZ==2&&totMass>mProt+mK) { piM = mK; piPDG = 321; } else if(totZ==1&&totMass>mProt+mK0&&G4UniformRand()>0.5) { piM = mK0; piPDG = 311; } else if(totZ==1&&totMass>=mNeut+mK) { nucM = mNeut; nucPDG= 2112; piM = mK; piPDG = 321; } else { G4cerr<<"***G4Q::HQ:Z="<PiNM&&!totS) // Decay in nucleon & pion { if(!totZ&&totMass>mProt+mPi&&G4UniformRand()<0.5) { piM = mPi; piPDG = -211; } else if(!totZ&&totMass>mNeut+mPi0) { nucM = mNeut; nucPDG= 2112; piM = mPi0; piPDG = 111; } else if(totZ==1&&totMass>mNeut+mPi&&G4UniformRand()<0.5) { nucM = mNeut; nucPDG= 2112; piM = mPi; piPDG = 211; } else if(totZ==1&&totMass>mProt+mPi0) { piM = mPi0; piPDG = 111; } else if(totZ==-1) { nucM = mNeut; nucPDG= 2112; piM = mPi; piPDG = -211; } else if(totZ==2) { piM = mPi; piPDG = 211; } else { G4cerr<<"*G4Q::HQ:Z="<1) { G4cerr<<"*G4Q::HQ:Z="<gam/pi/K("<GPK="< M="<dm&&!rWi) // Try to use the h-resonance width or reduce its spin { G4double sWi=G4QPDGCode(sPDG).GetWidth(); G4double sMM=G4QPDGCode(sPDG).GetMass(); if(sWi) // Hadron is a resonance { G4double mmm=theWorld->GetQParticle(G4QPDGCode(sPDG))->MinMassOfFragm(); G4double ddm=quasM-rMass; // Minimum mass of the sHadron if(fabs(sMM-ddm)<1.5*sWi-.001 && ddm>mmm) { #ifdef debug G4double msm=sMass; #endif sMass=GetRandomMass(sPDG,ddm); // Randomize mass of the Reson-Hadron if(fabs(sMass)<.001) { #ifdef debug G4cerr<<"***G4Q::HQ:ChangeToM=0, "<mPi0) { rPDG=111; reMass=mPi0; } else { rPDG=22; reMass=0.; } } G4double freeRQM=rQPDG.GetMass(); G4int RQB = rQPDG.GetBaryNum(); G4double fRQW= 3*rQPDG.GetWidth(); if(fRQW<.001) fRQW=.001; G4QPDGCode sQPDG(sPDG); G4int sChg=sQPDG.GetCharge(); G4int sBaryn=sQPDG.GetBaryNum(); G4double sCB=theEnvironment.CoulombBarrier(sChg,sBaryn); #ifdef debug G4cout<<"G4Q::HQ:h="<QM="<tM="<quasM) { if(mNeut+mK<=quasM+.001) { reMass=mNeut; rPDG =2112; rQPDG=G4QPDGCode(rPDG); rChg=rQPDG.GetCharge(); rBaryn=rQPDG.GetBaryNum(); rCB=theEnvironment.CoulombBarrier(rChg,rBaryn); #ifdef debug G4cout<<"G4Q::HQ:NCB="<RQ+QEX s="<rPDG="<rPDG="<rPDG="<mK0+nreM) { sMass=mK0; sPDG=311; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass curQ+=KpQC-K0QC; reMass=nreM; rPDG=nNuc.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<rPDG="<mPi0+nreM) { sMass=mPi0; sPDG=111; curQ+=PiQC; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreM; rPDG=nNuc.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<nreM) { sMass=0.; sPDG=22; curQ+=PiQC; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreM; rPDG=nNuc.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<rPDG="<0) // Switch From pi- To pi0 (or gamma) { G4QContent crQC=tmpQPDG.GetQuarkContent()-Pi0QC+PiMQC; // new hadron's QC G4QNucleus nNuc(crQC); // New neucleus for the hadron/fragment G4double nreM=nNuc.GetGSMass(); // Mass of the new isotope if(tmM>mPi0+nreM) { sMass=mPi0; sPDG=111; curQ+=PiMQC; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreM; rPDG=nNuc.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<nreM) { sMass=0.; sPDG=22; curQ+=PiMQC; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreM; rPDG=nNuc.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<mPi0+reMass) { sMass=mPi0; sPDG=111; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to pi0 mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<reMass) { sMass=0.; sPDG=22; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to gamma sum=reMass+sMass; if(fabs(tmM-reMass)rPDG="<0 && iniS>0) // Force Lamb->p+PiM (2/3) or Lamb->n+Pi0 decays @@ tot { G4QContent tmpSQC=G4QPDGCode(sPDG).GetQuarkContent();//QuarkContent of the hadron G4QContent lanQC=tmpQPDG.GetQuarkContent()+tmpSQC+K0QC;// switch from Lambda to n G4QNucleus nucM(lanQC-PiMQC); // New neucleus for the residual for Pi- G4double nreM=nucM.GetGSMass(); // Mass of the residual for Pi- G4QNucleus nucZ(lanQC-Pi0QC); // New neucleus for the residual for Pi- G4double nreZ=nucZ.GetGSMass(); // Mass of the residual for Pi- #ifdef debug G4cout<<"G4Q::HQ:LsPDG="<tmM) && mPi0+nreZ n+Pi0 case { sMass=mPi0; sPDG=111; curQ+=tmpSQC+K0QC; // LToN correction for curQC s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreZ; rPDG=nucZ.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<n)+Pi0 DecayIn2"); } } else if(mPi+nreM p+Pi- case { sMass=mPi; sPDG=-211; curQ+=tmpSQC+K0QC-PiMQC; // LToN correction for curQC (-QC_PIM) s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreM; rPDG=nucM.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<n)+Pi- DecayIn2"); } } else if(nreM N+gamma case { sMass=0.; sPDG=22; curQ+=tmpSQC+K0QC; // LToN correction for curQC s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to new hadron mass reMass=nreZ; rPDG=nucZ.GetPDG(); r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-sum)rPDG="<n)+Gamma DecayIn2"); } } else { G4cerr<<"***G4Q::HQ: LamToN M="<rM="<iniQM) { G4QContent tmpSQC=G4QPDGCode(sPDG).GetQuarkContent();//QuarkContent of the hadron sMass=0.; sPDG=22; s4Mom=G4LorentzVector(0.,0.,0.,sMass); // Switch to gamma curQ+=tmpSQC; // totQC correction for curQC reMass=iniQM; rPDG=iniPDG; r4Mom=G4LorentzVector(0.,0.,0.,reMass);// Switch to other isotope mass sum=reMass+sMass; if(fabs(tmM-reMass)rPDG="<resTNM+sMass) // Just decay in sPDG and total residual nucleus { G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,resTNM); //GSM of ResidTotEnvir G4LorentzVector rs4M = G4LorentzVector(0.,0.,0.,sMass); //GSM of a Hadron #ifdef debug G4cout<<"G4Q::HQ:EMERGENCY,rEM="<R="<totM) // Just decay in minimal total nucleus and gamma { G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,totM); // GSM of ResidTotEnvir G4LorentzVector rs4M = G4LorentzVector(0.,0.,0.,0.); // GSM of a Photon #ifdef debug G4cout<<"G4Q::HQ:EMERGENSY,minM="<g+M="<rPDG="<QHVect s4M="< "KinEn is below CB, try once more" case { #ifdef pdebug G4cout<<"****G4Q::HQ:E-9: sKE="<MINPDG) q4Mom-=G4LorentzVector(0.,0.,0.,pMass); qEnv=theEnvironment; return theQHadrons; } else if(abs(reMass-freeRQM)"ResidQ is a GSHadron/Frag" case { G4QHadron* curHadr1 = new G4QHadron(rPDG,r4Mom);// Create RealHadron for the ResidQ FillHadronVector(curHadr1); // Fill "new curHadr1" (del. eq.) G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);// Creation Hadron for theCandidate FillHadronVector(curHadr2); // Fill "new curHadr2" (del. eq.) #ifdef rdebug G4cout<<"G4Q::HQ:DecayQuasmon "<1&&totMass>totM&&totS>=0) //@@ ?? #ifdef pdebug G4cout<<"G4Q::HQ:*NEEDS-EVAP-10* M="<MINPDG) q4Mom-=G4LorentzVector(0.,0.,0.,pMass); qEnv=theEnvironment; return theQHadrons; } else // Only theTotResidNucl can evaporate { G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);// Create Hadron for theOutHadron FillHadronVector(curHadr2); // Fill "new curHadr2" (del.equiv.) q4Mom = r4Mom; if(sPDG>MINPDG) { theEnvironment.Reduce(pPDG); // Update NuclEnviron after Q->RQ+QEXF valQ += transQC; // Update the Quark Content of Quasmon } else valQ = curQ; // Update the Quark Content of Quasmon #ifdef rdebug G4cout<<"OK***>G4Q::HQ:S="<G4Q::HQ:After,S="<>"<>"Hadron candidate in Vacuum" case //if(sPDG> "Hadron candidate with the only Quasmon" case if(sPDG rK0M + mK0) { rPDG = rK0PDG; // PDG of the Residual System to K0 rMass = rK0M; sPDG = 311; sMass = mK0; } else { rPDG = rKPPDG; // PDG of the Residual System to K+ rMass = rKPM; sPDG = 321; sMass = mK; } G4double ctM=tot4M.m(); G4LorentzVector r4Mom(0.,0.,0.,rMass); G4LorentzVector s4Mom(0.,0.,0.,sMass); // Mass is random since probab. time G4double sum=rMass+sMass; if(fabs(ctM-sum) rPDG="<HadrVec s="<=iniBN) // Try to decay in K+ { sPDG = 321; sMass= mK; G4QNucleus totQN(valQ+KpQC); // Nucleus Residual after Kp sub. rPDG = totQN.GetPDG(); rMass= totQN.GetMZNS(); } else if(iniS<0) // Try to decay in K0 { sPDG = 311; sMass= mK0; G4QNucleus totQN(valQ+K0QC); // Nucleus Residual after K0 sub. rPDG = totQN.GetPDG(); rMass= totQN.GetMZNS(); } else if(iniQChg>iniBN-iniS) // Try to decay in Pi+ { sPDG = 211; sMass= mPi; G4QNucleus totQN(valQ-PiQC); // Nucleus Residual after Pi+ sub. rPDG = totQN.GetPDG(); rMass= totQN.GetMZNS(); } else if(iniQChg<0) // Try to decay in Pi- { sPDG = -211; sMass= mPi; G4QNucleus totQN(valQ+PiQC); // Nucleus Residual after Pi- sub. rPDG = totQN.GetPDG(); rMass= totQN.GetMZNS(); } else if(quasM>iniQM+mPi0) // Try to decay in Pi0 { sPDG = 111; sMass= mPi0; rPDG = iniPDG; rMass= iniQM; } else // Decay in gamma as a final decision { sPDG = 22; sMass= 0.; rPDG = iniPDG; rMass= iniQM; } ffin = true; } #ifdef debug G4cout<<"G4Q::HQ:MQ="<sPDG="< rPDG="<s="<NUCPDG) { G4int outBar = outQC.GetBaryonNumber(); G4double outCB = theEnvironment.CoulombBarrier(outChg,outBar);//ChrgIsNeglected // Now the CoulBar reflection should be taken into account G4double outT = s4Mom.e()-s4Mom.m(); outProb = theEnvironment.CoulBarPenProb(outCB,outT,outChg,outBar); } G4double rnd=G4UniformRand(); #ifdef debug G4cout<<"G4Q::HQ: for "<"HadronInNuclMedia or NuclCand" { // Now the CoulBar reflection should be taken into account G4QContent outQC = G4QPDGCode(sPDG).GetQuarkContent(); G4int outBar = outQC.GetBaryonNumber(); G4int outChg = outQC.GetCharge(); G4double outCB = theEnvironment.CoulombBarrier(outChg,outBar); // Now theCoulombPotential should be taken into account //@@ How to do this ?? if(nucflag) rQ4Mom+=dMass; // Make a correction of ResidQuasmon4M G4QHadron tmpRQH(valQ+transQC,rQ4Mom); // Tmp Hadron for the Residual Quasmon // Now theCoulBar reflection should be taken into account G4double outT = fr4Mom.e()-fr4Mom.m(); G4double outProb = theEnvironment.CoulBarPenProb(outCB,outT,outChg,outBar); if(G4UniformRand()RQ+F G4LorentzVector sumL=theEnvironment.Get4Momentum()+q4Mom; //@@ Check Print Only check += fr4Mom; //@@ Just for checking ccheck+=G4QPDGCode(sPDG).GetCharge(); //@@ Just for checking // --- @@ --- Potential recovery of the secondary --- Bad experience //G4double bindE=sMass-pMass; //G4double fE=fr4Mom.e(); //G4double nfE=fE+bindE; //G4double frM2=sMass*sMass; // MinSuaredMass of OutgoingFragment //G4double rpf=sqrt((nfE*nfE-frM2)/(fE*fE-frM2))-1.; //G4LorentzVector bind4M(rpf*fr4Mom.vect(),bindE); //G4LorentzVector ren4M=tot4M-fr4Mom-bind4M; //G4double rnM2= ren4M.m2(); // ResidNucleusMass after separation //G4QContent renQC=theEnvironment.GetQCZNS()+valQ; //G4QNucleus renTot(renQC); // PseudoNucleus for TotResidualNucl //G4double renTotM=renTot.GetMZNS(); // GSMass of Total Residual Nucleus //if(rnM2>renTotM*renTotM) //{ // fr4Mom+=bind4M; // rQ4Mom-=bind4M; //} // --- @@ --- End of Potential recovery of the secondary q4Mom = rQ4Mom; // Update the 4Mom of the Quasmon if(sPDG>MINPDG) valQ += transQC; // Update the Quark Content of Quasmon G4QHadron* candHadr = new G4QHadron(sPDG,fr4Mom);// Createa Hadron for Candidate FillHadronVector(candHadr); // Fill the RadiatedHadron(del.equiv.) #ifdef debug G4cout<<"G4Q::HQ:QuarkExchHadronizThroughCB Q="<>>>>>>>> NuclearMatter SUBCHECK >>>>>>"<envA*envA/400) // ?? M.K. ?? //{ status=3; if(gaF) { phot4M=zeroLV; gaF=false; } //} } } } // End of skip // === Check of boundary to escape not existing state of residual nucleus === if(CheckGroundState()) //if(CheckGroundState(true)) { ClearQuasmon(); // This Quasmon is done qEnv=theEnvironment; return theQHadrons; } G4LorentzVector sumLor=theEnvironment.Get4Momentum()+q4Mom+check; #ifdef debug G4int eZ = theEnvironment.GetZ(); G4int sumC = eZ+valQ.GetCharge()+ccheck; G4int curPDG=valQ.GetSPDGCode(); G4cout<<"G4Q::HQ:Z="<FinalCHECK***>>4M="< ResidualQ 4M="<GetPDGCode(); // Get PDG code of the Hadron to switch G4LorentzVector t = qH->Get4Momentum(); // 4-Mom of Chipolino #ifdef psdebug if(thePDG==113 && fabs(t.m()-770.)<.001) { G4cerr<<"G4Q::FillHadronVector: PDG="<Get4Momentum(); G4LorentzVector s4Mom = sHadr->Get4Momentum(); if(!qH->DecayIn2(f4Mom,s4Mom)) { delete fHadr; // Delete "new fHadr" delete sHadr; // Delete "new sHadr" G4cerr<<"***G4Q::FillHadrV:ChipQC"<Set4Momentum(f4Mom); // Put the randomized 4Mom to 1-st Hadron FillHadronVector(fHadr); // Fill 1st Hadron (delete equivalent) sHadr->Set4Momentum(s4Mom); // Put the randomized 4Mom to 2-nd Hadron FillHadronVector(sHadr); // Fill 2nd Hadron (delete equivalent) } } else if(thePDG>80000000&&thePDG!=90000000) //==Decay-Evaporation of theBarionicFragment== { G4double fragMas=t.m(); // Real Mass of the nuclear fragment //G4double fragMas=qH->GetMass(); // GrStMass of the nuclear fragment (wrong?) G4QNucleus qNuc(t,thePDG); // Make a Nucleus out of the Hadron // @@ Probably, when nucleus is initialized, the mass is not initialized ? Was OK! Why? //G4double GSMass =qNuc.GetGSMass(); // GrState Mass of the nuclear fragment (?) G4double GSMass = G4QPDGCode(thePDG).GetMass(); // More robust definition G4QContent totQC=qNuc.GetQCZNS(); // Total Quark Content of Residual Nucleus G4int nN =qNuc.GetN(); // A#of neutrons in the Nucleus G4int nZ =qNuc.GetZ(); // A#of protons in the Nucleus G4int nS =qNuc.GetS(); // A#of protons in the Nucleus G4int bA =qNuc.GetA(); // A#of baryons in the Nucleus #ifdef pdebug G4cout<<"G4Quasm::FillHadrVect:Nucl="<0) // => "Anti-strangeness or ISOBAR" case { G4double m1=mPi; // Prototypes for the nZ<0 case G4int PDG1=-211; G4QNucleus newNpm(totQC+PiQC); G4int newS=newNpm.GetStrangeness(); if(newS>0) newNpm=G4QNucleus(totQC+PiQC+newS*K0QC); G4int PDG2=newNpm.GetPDG(); G4double m2=newNpm.GetMZNS(); if(nS<0) { m1 =mK; PDG1 =321; G4QNucleus newNp(totQC-KpQC); PDG2 =newNp.GetPDG(); m2 =newNp.GetMZNS(); G4QNucleus newN0(totQC-K0QC); G4double m3=newN0.GetMZNS(); if (m3+mK0 "aK0+ResA is better" case { m1 =mK0; PDG1=311; m2 =m3; PDG2=newN0.GetPDG(); } } else if(nS>0&&nZ+nN>0) { if(nN<0) { m1 =mSigP; PDG1 =3222; G4QNucleus newNp(totQC-sigpQC); PDG2 =newNp.GetPDG(); m2 =newNp.GetMZNS(); } else { m1 =mSigM; PDG1 =3112; G4QNucleus newNp(totQC-sigmQC); PDG2 =newNp.GetPDG(); m2 =newNp.GetMZNS(); } } else if(nN<0) { PDG1 =211; G4QNucleus newNpp(totQC-PiQC); PDG2 =newNpp.GetPDG(); m2 =newNpp.GetMZNS(); } if(fragMas>m1+m2) // => "can decay" case { G4LorentzVector fq4M(0.,0.,0.,m1); G4LorentzVector qe4M(0.,0.,0.,m2); if(!qH->DecayIn2(fq4M,qe4M)) { G4cerr<<"***G4Quasm::FillHadrV: QM="< Mes="<SetNFragments(2); // Put a#of Fragments=2 //theQHadrons.push_back(qH); // Instead delete qH; // G4QHadron* H1 = new G4QHadron(PDG1,r1*t); theQHadrons.push_back(H1); // (delete equivalent) G4QHadron* H2 = new G4QHadron(PDG2,r2*t); FillHadronVector(H2); // (delete equivalent) } else { #ifdef debug G4cerr<<"-Warning-G4Q::FillHVec:PDG="<0&&bA>1) // It's nucleus and there is the neutron { 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 residualMass for the proton G4int pResPDG=0; // Prototype of PDGCode of the proton if(nZ>0&&bA>1) // It's nucleus and there is theroton { 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 residualMass for the Lambda G4int lResPDG=0; // Prototype of PDGCode of the Lambda if(nS>0&&bA>1) // It's nucleus and there is the Lambda { 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 } #ifdef debug G4cout<<"G4Quasm::FillHadrVec:rP="<DecayIn2(a4Mom,b4Mom)) { theQHadrons.push_back(qH); // No decay (delete equivalent) G4cerr<<"---Warning---G4Q::FillHadronVector: Be8 decay did not succeed"<SetNFragments(2); // Fill a#of fragments to decaying Hadron //theQHadrons.push_back(qH); // Fill hadron with nf=2 (del. eq.) // Instead delete qH; // G4QHadron* HadrB = new G4QHadron(barPDG,a4Mom); FillHadronVector(HadrB); // Fill 1st Hadron (delete equivalent) G4QHadron* HadrR = new G4QHadron(resPDG,b4Mom); FillHadronVector(HadrR); // Fill 2nd Hadron (delete equivalent) } } else { #ifdef debug G4cout<<"G4Quasm::FillHadrVect: Leave as it is"<GSM="<>>G4Quasm::FillHadrVect: Leave as it is Instead of Exception"<GSMass) { G4int gamPDG=22; G4double gamM=0.; if(fragMas>mPi0+GSMass) { gamPDG=111; gamM=mPi0; } G4LorentzVector a4Mom(0.,0.,0.,gamM); G4LorentzVector b4Mom(0.,0.,0.,GSMass); if(!qH->DecayIn2(a4Mom,b4Mom)) { theQHadrons.push_back(qH); // No decay (delete equivalent) G4cerr<<"---Warning---G4Q::FillHadrVect:N*->gamma/pi0+N decay error"<Set4Momentum(b4Mom); // Put the new 4-mom in the residual GS fragment theQHadrons.push_back(qH); // Fill corrected baryon in the HadronVector } } else // ===> Evaporation of excited system { #ifdef ppdebug G4cout<<"G4Quasm::FillHadrVect:Evaporate "< GS="<Products are added to QHV, nQHV="<clear(); delete tmpQHadVec; // That who calls DecayQHadron is responsible for clear & delete #ifdef pdebug G4cout<<"G4Q::FillHV: TemporaryQHV of DecayProducts is deleted"<kMin QParton Momentum for the current Quasmon #ifdef debug //G4cout<<"G4Quasmon::GetQPartonMom:**called**mR="<=kLim) vRndm = vRndm*pow((1.-xMin),n)*(1.+n*xMin); // "xMin - 1." Range else { G4double xMax=kMax/kLim; G4double vRmin = pow((1.-xMin),n)*(1.+n*xMin); G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax); vRndm = vRmax + vRndm*(vRmin-vRmax); // Randomization in the "xMin - xMax" Range } } else if (kMax1.) { //G4cout<<"-Warning-G4Quasmon::GetQPM: R="<.999999999) x=.999999999; G4double r = 1.-x; G4double p = r; if (n>2) p = pow(r,n-1); // (1-x)**(n-1) G4double nx = n*x; G4double c = p*r*(1.+nx); // vRndm=(1-x)**n*(1+n*x) is the equation too be solved G4double od = c - vRndm; // the old Residual Difference #ifdef debug G4cout<<"G4Q::GetQPMom:>>>First x="<.999999999) x=.999999999; r = 1.-x; if (n>2) p = pow(r,n-1); else p = r; nx = n*x; c = p*r*(1.+nx); d = c - vRndm; } else { if (f>1.0001&&f<27.) f=1.+(f-1.)*2; // Make a regular correction of OptStepFactor if (f<0.99999999999) f=1.+(1.-f)/2.; if (f>=27.) f=27.; } #ifdef debug G4cout<<"G4Q::GetQPMom: Iter#"< nitMax="<xMax) x=xMax; if(x=kMax)kCand=kMax-.001; if(kCand<=kMin)kCand=kMin+.001; return kCand; } // ********** Possible Performance Improvement (corrupt physics!) ********** // This is a simplified algorithm, which takes an empirical reduction 1/k of HadrProbab // As a result the equation to be solved is just R=(1-x)**n=n*INT_x^1[(1-z)**(n-1)*dz] //if (kMin>0.) // ==> There is a minimum cut for the QuarkPartonMomentum //{ // G4double xMin=kMin/kLim; // Minimal value for "x" // if (kMax>=kLim) vRndm = vRndm*pow((1.-xMin),n); // Shrink to the "xMin - 1." Range // else // { // G4double xMax=kMax/kLim; // Maximum value for "x" // G4double vRmin = pow((1.-xMin),n); // G4double vRmax = pow((1.-xMax),n); // vRndm = vRmax + vRndm*(vRmin-vRmax); // Randomization in the "xMin - xMax" Range // } //} //else if (kMax at this point kMin<=0 -> kMin=0 //{ // G4double xMax=kMax/kLim; // Maximum value for "x" // G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax); // vRndm = vRmax + vRndm*(1.-vRmax); // Shrink to the "0 - xMax" Range //} //// For the kMin=0, kMax=1 the normalization is not needed //if (n==1) return kLim*(1.-vRndm); // Direct solution for nOfQ==3 //else if (n==2) return kLim*(1.-sqrt(vRndm)); // Direct solution for nOfQ==4 //else return kLim*(1.-pow(vRndm,1./fn)); // Direct solution for nOfQ>4 } // End of "GetQPartonMomentum" // For the given quasmon mass calculate a number of quark-partons in the system G4int G4Quasmon::CalculateNumberOfQPartons(G4double qMass) // ==================================================== { static const G4double mK0 = G4QPDGCode(311).GetMass(); // @@ Temporary here. To have 3 quarks in Nucleon Temperature should be < M_N/4 (234 MeV) // M^2=4*n*(n-1)*T^2 => n = M/2T + 1/2 + T/4M + o(T^3/16M^3) // @@ Genius (better than 10**(-3) even for n=2!) but useless //G4double qMOver2T = qMass/(Temperature+Temperature); //G4double est = qMOver2T+1.+0.125/qMOver2T; // @@ Longer but exact G4double qMOverT = qMass/Temperature; G4int valc = valQ.GetTot(); // ................................. // --- Exponent, Double Split, Poisson 1 ============ ///G4int b = valQ.GetBaryonNumber(); ///G4int mq= 3*b; ///if (!b) mq=2; ///G4double mean = ((1.+sqrt(1.+qMOverT*qMOverT))/2. - mq)/2.; ///if(mean<0.) nOfQ=mq; // --- Uncomment up to here ================^^^^^^^^^ // Exponent ------ //else nOfQ=mq-2*mean*log(G4UniformRand()); // Poisson 1 ------ ///else nOfQ=mq+2*RandomPoisson(mean); // Double Split ------ //else //{ // G4int imean = static_cast(mean); // G4double dm = mean - imean; // if(G4UniformRand()>dm) nOfQ=mq+imean+imean; // else nOfQ=mq+imean+imean+2; //} // ......... // Poisson 2 ============= //if(valc%2==0)nOfQ = 2*RandomPoisson((1.+sqrt(1.+qMOverT*qMOverT))/4.);// abs(b) is even //else nOfQ = 1+2*RandomPoisson((1.+sqrt(1.+qMOverT*qMOverT))/4.-0.5);// abs(b) is odd // Poisson 3 ============= nOfQ = RandomPoisson((1.+sqrt(1.+qMOverT*qMOverT))/2.); G4int ev = valc%2; if (!ev && nOfQ<2) nOfQ=2; // #of valence quarks is even else if ( ev && nOfQ<3) nOfQ=3; // #of valence quarks is odd // #ifdef debug G4cout<<"G4Q::Calc#ofQP:QM="<stran) stran=astra; G4int nMaxStrangeSea=static_cast((qMass-stran*mK0)/(mK0+mK0));//KK is min for s-sea if (absb) nMaxStrangeSea=static_cast((qMass-absb)/672.); //LambdaK is min for s-sea #ifdef debug G4cout<<"G4Q::Calc#ofQP:"<0)valQ.IncQAQ(nSeaPairs,SSin2Gluons); else morDec=valQ.DecQAQ(-nSeaPairs); #ifdef debug if(morDec) G4cout<<"G4Q::Calc#ofQP: "<nMaxStrangeSea) // @@@@@@@ Too many strange sea ?????? { #ifdef debug G4cout<<"G4Q::Calc#ofQP:**Reduce** S="<GetPDGCode(); // PDGC of the Candidate G4bool poss = curCand->GetPossibility(); // Possibility for the Candidate G4QContent tmpTQ=totQC-curCand->GetQC(); G4QNucleus tmpT(tmpTQ); // Nucleus for TotResidNucleus for Fragment G4double tmpTM=tmpT.GetMZNS(); // GSMass of TotalResidNucleus for Fragment G4QPDGCode cQPDG(cPDG); // QPDG for the candidate G4double frM=cQPDG.GetMass(); // Vacuum mass of the candidate if(cPDG>80000000&&cPDG!=90000000) // Modify Fragments toTakeIntoAccount CurNuc { if(totMasstM=" <SetPossibility(false); } G4QNucleus cNuc(cPDG); // Fake nucleus for the Candidate G4int cP = cNuc.GetZ(); // A#of protons in the Current Environment G4int cN = cNuc.GetN(); // A#of neutrons in the Current Environment G4int cL = cNuc.GetS(); // A#of lambdas in the Current Environment G4QPDGCode cQPDG(cPDG); // QPDG of the Current Cluster #ifdef debug if(cPDG==90001000) G4cout<<"G4Q::MIM:>>>cPDG=90001000<<<,possibility="<=cP&&eN>=cN&&eL>=cL&&poss) // Cluster exists & possible { G4double clME = 0.; // Prototype of the BoundClMass in Environ G4double clMN = 0.; // Prototype of the BoundClMass in TotNucl G4double renvM = 0; // Prototype of the residual Environ mass if(cP==eP&&cN==eN&&cL==eL)clME=cQPDG.GetMass();// The only notBoundCluster of Envir else // Bound Cluster in the Environment { renvM = cQPDG.GetNuclMass(eP-cP,eN-cN,eL-cL); // Mass of residual for Environment clME = envM-renvM; } if(cP==tP&&cN==tN&&cL==tL)clMN=cQPDG.GetMass(); // The only NotBoundCluster of TotN else // Bound Cluster in Total Nucleus { renvM = cQPDG.GetNuclMass(tP-cP,tN-cN,tL-cL); // TotalResidualNucleus Mass clMN = totM-renvM; } curCand->SetParPossibility(true); curCand->SetEBMass(clME); curCand->SetNBMass(clMN); #ifdef sdebug G4int envPDGC = theEnvironment.GetPDGCode(); // PDG Code of Current Environment G4cout<<"G4Q:ModInMatCand:C="<SetParPossibility(false); } // @@ Modification of hadron masses in nuclear matter are not implemented yet } } // End of "ModifyInMatterCandidates" // Randomize the Resonance masses and calculate probabilities of hadronization for them void G4Quasmon::CalculateHadronizationProbabilities (G4double E, G4double kVal, G4LorentzVector k4M,G4bool piF, G4bool gaF, G4bool first) // =====================================================================E is not used==== { // ^ static const G4double mPi0 = G4QPDGCode(111).GetMass(); // | /////////static const G4double mEta = G4QPDGCode(221).GetMass(); // | G4double kLS=E; // | kLS=k4M.e(); // Temporary trick to avoid worning G4int vap = nOfQ-3; // Vacuum power //G4double kLSi= kLS; // Initial (without photon) kLS //if(addPhoton) kLSi=kLS-addPhoton; // @@ probabilities for k+gam can be wrong G4double mQ2 = q4Mom.m2(); // Squared Mass of the Quasmon G4double eQ = q4Mom.e(); // LS Energy of the Quasmon G4double mQ = sqrt(mQ2); // Mass of the decaying Quasmon G4double dk = kVal + kVal; // Double momentu of quark-parton in QCM G4double rQ2 = mQ2-dk*mQ; // Min Residual Colored Quasmon Squared Mass //////////////G4double rQ = sqrt(rQ2); G4double mQk = mQ-dk; // For acceleration G4double var = theEnvironment.GetProbability();// Vacuum to medium ratio G4double vaf = 0; //@@ !! Vacuum factor if(vap>0)vaf = var*mQk/kVal/vap; //@@ //if(vap>0)vaf = mQk/kVal/vap; //@@ VacuumFactor(instead of in G4QNucleus) G4double accumulatedProbability = 0.; G4double secondAccumProbability = 0.; G4int qBar =valQ.GetBaryonNumber(); // BaryNum of Quasmon G4int nofU = valQ.GetU()- valQ.GetAU(); // A#of u-quarks G4int dofU = nofU+nofU; G4int nofD = valQ.GetD()- valQ.GetAD(); // A#of d-quarks G4int dofD = nofD+nofD; G4int qChg = valQ.GetCharge(); G4int qIso = qBar-qChg-qChg; // Charge of Quasmon ////////////////G4int aQuI = abs(qIso); // Abs value for estimate ///////////////G4int oeQB = qBar%2; // odd(1)/even(0) QBaryn flag G4int maxC = theQCandidates.size(); // A#of candidates G4double totZ = theEnvironment.GetZ() + valQ.GetCharge(); // Z of the Nucleus ///////////////G4double totA = theEnvironment.GetA() + valQ.GetBaryonNumber(); //A of Nuc G4double envM = theEnvironment.GetMass(); // Mass of the Current Environment G4int envPDGC = theEnvironment.GetPDGCode();// PDG Code of Current Environment G4int envN = theEnvironment.GetN(); // N of current Nuclear Environment G4int envZ = theEnvironment.GetZ(); // Z of current Nuclear Environment ////////////////G4int envS = theEnvironment.GetS(); // S of CurrentNuclearEnvironment G4int envA = theEnvironment.GetA(); // A of current Nuclear Environment G4QContent envQC=theEnvironment.GetQCZNS(); // QuarkContent of the CurrentNuclearEnviron. //////G4double addPhoton=phot4M.e(); // PhotonEn for capture by quark-parton #ifdef debug G4int absb = abs(qBar); // Abs BaryNum of Quasmon G4int maxB = theEnvironment.GetMaxClust(); // Maximum BaryNum for clusters G4cout<<"G4Q::CalcHadronizationProbab:Q="<>>c="< tmpTM+frM="<80000000&&cPDG!=90000000&&dUD<1)))// 1 ** never try if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000&&dUD<2)))//2 ***The best*** //if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000&&dUD<3))) // 3 *** good *** //if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000&&dUD<4)))//4 almost the same //if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000))) // no restrictions { G4QContent curQ = valQ; // Make current copy of theQuasmonQuarkCont G4int baryn= candQC.GetBaryonNumber(); // Baryon number of the Candidate G4int cC = candQC.GetCharge(); // Charge of the Candidate G4double CB=0.; if(envA) CB=theEnvironment.CoulombBarrier(cC,baryn); /////////////G4int cI = baryn-cC-cC; //G4int cNQ= candQC.GetTot()-1-baryn; // A#of quarks/diquarksInTheCandidate - 2 G4int cNQ= candQC.GetTot()-2; // #of quark-partonsInTheCandidate - 2 (OK) //G4int cNQ= candQC.GetTot()+baryn-2; // A#of q-partons+b_Sj-2 (string junction) //G4int cNQ= candQC.GetTot()+3*baryn-4; // A#of q-partons+b+2*(b-1)-2 (choc q-link) G4double resM=0.; // Prototype for minMass of residual hadron #ifdef debug if(pPrint)G4cout<<"G4Q::CHP:B="<80000000&&cPDG!=90000000&&baryn<=envA)//==>Nuclear Fragment (QUarkEXchange) { G4int pc=0; // Parent counter counter G4double pcomb=0.; // Summed probability of parent clusters G4double frM2=frM*frM; // Squared mass of the nuclear fragment G4double qMax=frM+CB-kLS; // ParClustM-qmax value (k-q>frM-prM+CB) //G4double qMax=frM-kLS; // ParClustM-qmax value(k-q>frM-prM+) //////////G4double qM2=qMax+qMax; G4int iQmin=0; // IncomingToCluster quarks are d=0,u=1,s=2 G4int iQmax=3; // 3 is bigger than s-quark (2, iq<3) G4int oQmin=0; // Returning from cluster quarks: d=0,u=1 G4int oQmax=2; // 2 is bigger than u-quark (1, oq<2) @@ 3? if (dofU<=nofD) iQmax=1; // Too many Dquarks (in-Uquark is forbiden) else if(dofD<=nofU) iQmin=1; // Too many Uquarks (in-Dquark is forbiden) // @@ This is how netrons are increased for the pion capture at rest case @@ if(piF) // force Pi- transfer its charge to a quark { iQmin=0; iQmax=1; } #ifdef debug if(pPrint)G4cout<<"G4Q::CHP:***!!!***>>F="<"<0.) // Kinematical analysis of hadronization { #ifdef debug if(pPrint) G4cout<<"G4Q::CHP:fM="<="<"< 0. && nz < 1. && nz < ne) { ne=nz; newh=pow(nz,cNQ); if(newh < kf) kf=newh; } else if(nz <= 0.) kf=0.; } // *** VBQ *** CHECK Residual Quazmon (Never use: reduces PS) //if(minBM2>rQ2&&!piF&&!gaF&&baryn>3) // ==>Check ResidVirtualQuasm //if(minBM2>rQ2&&!piF&&!gaF&&baryn>2) // ==>Check ResidVirtualQuasm //if(minBM2>rQ2&&!piF&&!gaF) // ==> Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&baryn>2)//==>Check of ResidualVirtualQuasmon**OK** //if(minBM2>rQ2&&baryn>1)//==>Check ResidualVirtualQuasm **Better** //if(minBM2>rQ2&&piF&&(cPDG==90000001||cPDG==90002002))//CheckRVirQ //if(minBM2>rQ2&&piF&&(cPDG==90000001||baryn>3))//CheckResidVirtQua //if(minBM2>rQ2&&(piF&&cPDG==90000001||baryn>2))//CheckResidVirtQua //if(minBM2>rQ2&&baryn>2) // ==> Check of Residual Virtual Quasmon //if(minBM2>rQ2&&!piF&&baryn>2&&cPDG!=90001002)//ResidVirtQuasmon //if(minBM2>rQ2&&!piF&&baryn>2)//==>Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&!piF&&baryn>3)//==>Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&!piF&&baryn>1)//==>Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&!piF&&!gaF)// ==> Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&!piF)// ==> Check of ResidualVirtualQuasmon ALWAYS //if(minBM2>rQ2&&piF&&baryn>3)//==>Check of ResidualVirtualQuasmon //if(minBM2>rQ2&&piF&&(baryn==1||baryn>2))//==>Check ResidVirtQ //if(minBM2>rQ2&&(!piF||piF&&(cPDG!=90000001||G4UniformRand()<.5))) //if(minBM2>rQ2&&(!piF||piF&&(cPDG!=90000001))) //if(minBM2>rQ2&&(!piF&&baryn>4 || piF && cPDG!=90000001 && if (minBM2 > rQ2 && ( !piF || ( piF && cPDG != 90000001 && cPDG != 90001001 && cPDG != 90001002 ) ) ) //if(minBM2>rQ2) // ==> Check of Residual (Virtual?) Quasmon //if(2>3) { G4double nz=0.; if(atrest) nz=1.-(minBM2-rQ2+pmk*dked)/(boundM*(rEP+pmk)); else nz=1.-(minBM2-rQ2)/(boundM*rEP); //if(atrest) nz=1.-(minBM2-rQ2+pmk*dkd)/(nucBM*(rEP+pmk)); //else nz=1.-(minBM2-rQ2)/(nucBM*rEP); #ifdef debug if(pPrint) G4cout<<"G4Q::CHP:q="<"<0.&&nz<1.&&nzrQ2&&baryn>3) // ==> Check of Residual Quasmon //if(minM2>rQ2&&!piF&&!gaF&&baryn>3)// ==> Check of ResidualQuasmon //if(minM2>rQ2&&!piF&&baryn>1) // ==> Check of ResidualQuasmon //if(minM2>rQ2&&!piF&&baryn>2) // ==> Check of ResidualQuasmon //if(minM2>rQ2&&!piF&&baryn>2&&cPDG!=90001002)//=>CheckResidQuasmon //if(minM2>rQ2&&!piF&&baryn>3) // ==> Check Residual Quasmon **OK** //if(minM2>rQ2&&!piF&&!gaF) // ==> Check of Residual Quasmon //if(minM2>rQ2&&!piF) // ==> Check of Residual Quasmon //if(minM2>rQ2&&piF) // ==> Check of Residual Quasmon //if(minM2>rQ2&&baryn>1) // ==> Check Residual Quasmon **Better** //if(minM2>rQ2&&(baryn>1||!piF))//==>CheckResidualQuasmon**Better** //if(minM2>rQ2&&baryn>1&&cPDG!=90002002) //==> CheckResidualQuasmon //if(minM2>rQ2&&!piF) // ==> Check of Residual Quasmon //if(minM2>rQ2&&baryn>3) //=>CheckResidQuasmon *** The Best *** //if(minM2>rQ2 && (!piF || piF && //if(minM2>rQ2 && (!piF&&baryn>3 || piF && if ( minM2 > rQ2 && ( (!piF && baryn > 4) || (piF && (cPDG != 90000001 || G4UniformRand() > .3333333) && cPDG != 90001001) ) ) //cPDG!=90001001) ) //if(minM2>rQ2) // ==> Check of Residual Quasmon //if(2>3) { G4double nz=0.; if(atrest) nz=1.-(minM2-rQ2+pmk*dked)/(boundM*(rEP+pmk)); else nz=1.-(minM2-rQ2)/(boundM*rEP); //if(atrest) nz=1.-(minM2-rQ2+pmk*dkd)/(nucBM*(rEP+pmk)); //else nz=1.-(minM2-rQ2)/(nucBM*rEP); #ifdef debug if(pPrint) G4cout<<"G4Q::CHP:q="<"<0.&&nz<1.&&nz1.)kf=1.; G4double high = kf; // after this kf can be changed #ifdef debug if(pPrint) G4cout<<"G4Q::CHP:"<0.&&lz<1.)low=pow(lz,cNQ); else if(lz>=1.)low=1.; #ifdef debug G4double pl=low; // Just for printing if(pPrint) G4cout<<"G4Q::CHP:="<="<0.&&nz>lz) { newl=pow(nz,cNQ); if(newl>low) low=newl; } else if(nz>1.) low=10.; } // ***** End of restrictions ***** kf-=low; #ifdef debug if(pPrint) G4cout<<"G4Q::CHP:>>"< End of Quark Exchange "iq" Test LOOP } // ---> End of Nuclear Case of fragmentation else if(cPDG<80000000) // ===> Hadron case (QUark FUsion mechanism) { // Calculation of the existing hadrons G4int curnh=theQHadrons.size(); G4int npip=0; G4int npin=0; G4int npiz=0; for (G4int ind=0; indGetPDGCode(); // PDG Code of the hadron if (curhPDG== 111) npiz++; if (curhPDG== 211) npip++; if (curhPDG==-211) npin++; } // End of the hadron counting comb = valQ.NOfCombinations(candQC); if(!comb) { if ( (aPDG==111)|(aPDG==211) ) comb=1.; // Permit pions @@ ? else if ( (aPDG==311)|(aPDG==321) ) comb=SSin2Gluons; // Permit kaons @@ ? } if(cPDG== 211&&npip>0) comb*=(npip+1); // Bose multyplication for pi+ if(cPDG==-211&&npip>0) comb*=(npin+1); // Bose multyplication for pi- if(cPDG==111||cPDG==221||cPDG==331||cPDG==113||cPDG==223||cPDG==333||cPDG==115|| cPDG==225||cPDG==335||cPDG==117||cPDG==227||cPDG==337||cPDG==110||cPDG==220|| cPDG==330) // @@ Can it be shorter if? { G4QContent tQCd(1,0,0,1,0,0); G4QContent tQCu(0,1,0,0,1,0); G4QContent tQCs(0,0,1,0,0,1); G4double cmd=valQ.NOfCombinations(tQCd); G4double cmu=valQ.NOfCombinations(tQCu); G4double cms=valQ.NOfCombinations(tQCs); if(cPDG!=333&&cPDG!=335&&cPDG!=337) comb=(cmd+cmu)/2.; //if(cPDG==331||cPDG==221) comb =(comb + cms)/2.; //eta,eta' if(cPDG==331||cPDG==221) comb =(comb + cms)/4.; //eta,eta'(factor2 suppression) if(cPDG==113) comb*=4.; //@@ if(cPDG==223) comb*=2.; //@@ if(cPDG==111&&npiz>0) comb*=(npiz+1); // Bose multyplication #ifdef debug if(abs(cPDG)<3) G4cout<<"G4Q::CHP:comb="<>cPDG="<"< cM="<0. && vap>0 && zMax>zMin) { probability = vaf*(pow(zMax, vap)-pow(zMin, vap)); #ifdef debug if(priCon) G4cout<<"G4Q::CHP:#"< LowBaryonNumber case (tuned on p-ap) //{ if(cPDG==110||cPDG==220||cPDG==330) probability*=comb; // f0 has spin 0 else probability*=comb*(abs(cPDG)%10); // Spin of resonance G4int BarRQC=curQ.GetBaryonNumber(); // Res Quasmon BaryonNumber G4int StrRQC=curQ.GetStrangeness(); // Res Quasmon Strangeness if(BarRQC==2 && !StrRQC) // --> DiBaryon Correction { G4int ChgRQC=curQ.GetCharge(); // Res Quasmon Charge if(ChgRQC==1) probability/=2; // Only one S else probability*=2; // One S + three P } //} } } else { #ifdef debug if(priCon) G4cout<<"G4Q::CHP:cM=0[cPDG"<0"< 0 && corFlag) { // *** CORRECTION *** G4QHadron* theLast = theQHadrons[nOfOUT-1]; if(!(theLast->GetNFragments()) && theLast->GetPDGCode()!=22)//NotDecayedHadron & NotGam { G4LorentzVector hadr4M=theLast->Get4Momentum(); G4double hadrMa=hadr4M.m(); // Mass of the Last hadron (==GSMass) G4LorentzVector tmpTLV=reTLV+hadr4M;// Tot (ResidNucl+LastHadron) 4-Mom #ifdef debug G4cout<<"G4Q::CheckGS:YES,T="<rM+hM="<resSMa+hadrMa) // resMa contains 2 Hadrons: resQ and Environ { if(resEMa) // => "Non vacuum Environment exists" case { G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa); // GS Mass of Quasmon if(!G4QHadron(tmpTLV).DecayIn3(hadr4M,quas4M,enva4M)) { G4cerr<<"---Warning---G4Q::CheckGS:DecIn Fragm+ResQ+ResEnv Error"<Set4Momentum(quas4M); // @@ FillHadronVector(quasH); // Fill ResidQuasm Hadron (del.equiv.) theLast->Set4Momentum(hadr4M); } } else //=>"The Env is vacuum" case (DecayIn2) { G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa); // GS Mass of Quasmon G4QHadron* quasH = new G4QHadron(valQ, quas4M); if(!G4QHadron(tmpTLV).DecayIn2(hadr4M,quas4M)) { delete quasH; // Delete "new Quasmon Hadron" G4cerr<<"---Warning---G4Q::CheckGS: Decay in Fragm+ResQ Error"<Set4Momentum(hadr4M); quasH->Set4Momentum(quas4M); FillHadronVector(quasH); // Fill ResidQuasmHadron (del.equivalent) } } } else // "CORRECTION" !!!! { if(nOfOUT>1 && corFlag) { G4QHadron* thePrev = theQHadrons[nOfOUT-2];// Get pointer to prev-before-lastHad if(thePrev->GetNFragments()||thePrev->GetPDGCode()==22)return false;//DecH or Gam G4LorentzVector prev4M=thePrev->Get4Momentum(); // 4M of thePreviousButLastHadron G4double prevMa=prev4M.m(); // PreviousHadronMass (==HadrGSM) tmpTLV+=prev4M; // IncrementTotal4M of TotResNucl G4QContent totQC=valQ+envaQC; // QCont of theResidNucl=ResQ+Env G4int totPDG=totQC.GetSPDGCode(); // PDG Code of TotResidualNucleus G4double totQMa=G4QPDGCode(totPDG).GetMass(); // GS Mass of the ResidualNucleus G4double totNMa=tmpTLV.m(); // RealMass of TotResidualNucleus #ifdef debug G4cout<<"G4Q::CheckGS:NO, M="<"<totQMa+hadrMa+prevMa) { G4LorentzVector nuc4M = G4LorentzVector(0.,0.,0.,totQMa); // ResNuclAtRest 4Mom if(!G4QHadron(tmpTLV).DecayIn3(hadr4M,prev4M,nuc4M)) { G4cerr<<"---Warning---G4Q::CheckGS:DecIn3 ResN+Last+Prev Error"<Set4Momentum(hadr4M); thePrev->Set4Momentum(prev4M); #ifdef debug G4cout<<"G4Q::CheckGS: Yes, Check D4M="<GetQC(); G4int resLPDG =tmpLNQ.GetSPDGCode(); G4double resLastM=G4QPDGCode(resLPDG).GetMass();//GSM of ResidNucl for theLastH G4QContent tmpPNQ=totQC+theLast->GetQC(); G4int resPPDG =tmpPNQ.GetSPDGCode(); G4double resPrevM=G4QPDGCode(resPPDG).GetMass();//GSM of ResidNucl for thePrevH //////G4bool which = true; // LastH is absorbed, PrevH is radiated #ifdef debug G4cout<<"G4Quasm::CheckGS: NO, tM="< rp+l="< rl+p="<resLastM+hadrMa) // "Just exclude the Prev" case { theQHadrons.pop_back(); // theLast* is excluded from OUTPUT HV theQHadrons.pop_back(); // thePrev* is excluded from OUTPUT HV theQHadrons.push_back(theLast); // theLast substitutes thePrev in OUTPUT delete thePrev; // thePrev QHadron is destructed thePrev=theLast; resPPDG=resLPDG; resPrevM=resLastM; prev4M = hadr4M; } else if (totNMa>resPrevM+prevMa) // "Just exclude the Last" case { theQHadrons.pop_back(); delete theLast; } else return false; G4LorentzVector nuc4M = G4LorentzVector(0.,0.,0.,resPrevM);//ResNucl4m to PrevH if(!G4QHadron(tmpTLV).DecayIn2(prev4M,nuc4M)) { G4cerr<<"---Warning---G4Q::CheckGS:DecIn2 (ResN+Last)+Prev Error"<Set4Momentum(prev4M); #ifdef debug G4cout<<"G4Q::CheckGS:Yes, Check D4M="<GetQPDG(); G4int thePDG = theQPDG.GetPDGCode(); // Get the PDG code of decaying hadron G4int pap = 0; // --- particle if(thePDG<0) pap = 1; // --- anti-particle G4LorentzVector t = qH->Get4Momentum(); // Get 4-momentum of decaying hadron G4double m = t.m(); // Get the mass value of decaying Hadron // --- Randomize a channel of decay G4QDecayChanVector decV = theWorld->GetQParticle(theQPDG)->GetDecayVector(); G4int nChan = decV.size(); #ifdef debug G4cout<<"G4Quasm::DecQHadron: PDG="<1) { G4double rnd = G4UniformRand(); // Random value to select a Decay Channel for(i=0; iGetPDGCode() <<","<GetPDGCode(); if(nPart>2) G4cout<<","<GetPDGCode(); G4cout<3) { G4cerr<<"---Warning---G4Q::DecayQHadr:n="<push_back(qH); // Fill as it is (del.equiv.) return theFragments; } #ifdef debug G4cout<<"G4Q::DecQH:Decay("<GetMass()<<")->"<push_back(qH); // Fill as it is (del.equiv.) return theFragments; } else { //qH->SetNFragments(2); //theFragments.push_back(qH); // Fill with NFr=2 (del.equiv.) // Instead delete qH; // Delete it (without History) // fHadr->Set4Momentum(f4Mom); // Put the randomized 4Mom to 1-st Hadron G4QHadronVector* theTmpQHV=DecayQHadron(fHadr); // Try to decay G4int nProd=theTmpQHV->size(); #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH1="<push_back((*theTmpQHV)[0]);// Final = no Further Decay else for(G4int ip1=0; ip1size(); if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);// Final = no Further Decay else { theFragments->resize(tmpS+theFragments->size());// Resize theFragments length copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS); } #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec11 nProd="<clear(); delete intTmpQHV; } theTmpQHV->clear(); delete theTmpQHV; sHadr->Set4Momentum(s4Mom); // Put the randomized 4Mom to 2-nd Hadron theTmpQHV=DecayQHadron(sHadr); // Try to decay nProd=theTmpQHV->size(); #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH2="<push_back((*theTmpQHV)[0]);// Final = no Further Decay else for(G4int ip1=0; ip1size(); if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);// Final = no Further Decay else { theFragments->resize(tmpS+theFragments->size());// Resize theFragments length copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS); } #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec12 nProd="<clear(); delete intTmpQHV; } theTmpQHV->clear(); delete theTmpQHV; } #ifdef debug G4cout<<"G4Q::DecQHadr: DecayIn2 is made with nH="<size()<push_back(qH); // Fill hadron as it is (del.equivalent) } } else if(nPart==3) { G4QHadron* fHadr; G4QHadron* sHadr; G4QHadron* tHadr; G4int fPDG=cV[0]->GetPDGCode(); G4int sPDG=cV[1]->GetPDGCode(); G4int tPDG=cV[2]->GetPDGCode(); //The radiative decays of the GS hadrons In3 are closed if ElMaDecays=false if ( (fPDG != 22 && sPDG != 22 && tPDG != 22) || ElMaDecays) { #ifdef debug G4cout<<"G4Q::DQH:Y,f="<GetWidth()==0.) // Don't randomize theFirstHardon { fHadr = new G4QHadron(fPDG); // theFirst Hadron is created *1* if(cV[1]->GetWidth()==0.) { sHadr = new G4QHadron(sPDG); // theSecond Hadron is created *2* if(cV[2]->GetWidth()==0.)tHadr = new G4QHadron(tPDG);//theThirdHadron isCreated else { G4QParticle* tPart=theWorld->GetQParticle(cV[2]);// Pt for the3-rdH G4double tdm = m-fHadr->GetMass()-sHadr->GetMass();// MaxMass for the 2d Hadr tHadr = new G4QHadron(tPart,tdm); //the3rdHadron is created if(tPDG<0) tHadr->MakeAntiHadron(); } } else // Randomize 2nd & 3rd Hadrons { m-=fHadr->GetMass(); // Reduce the residual MaxMass G4QParticle* tPart=theWorld->GetQParticle(cV[2]);// Pt for the 3-rd Hadron G4double mim = tPart->MinMassOfFragm(); // MinMassLimit for the 3rd Hd G4double sdm = m - mim; // MaxMassLimit for the 2nd Hd G4QParticle* sPart=theWorld->GetQParticle(cV[1]);// Pt for the 2-nd Hadron sHadr = new G4QHadron(sPart,sdm); // theSecondHadron is created if(sPDG<0) sHadr->MakeAntiHadron(); G4double tdm = m - sHadr->GetMass(); // MaxMassLimit for the 3-rd H tHadr = new G4QHadron(tPart,tdm); // the Third Hadron is created if(tPDG<0) tHadr->MakeAntiHadron(); } } else // Randomize masses of all three Hadrons { G4QParticle* sPart=theWorld->GetQParticle(cV[1]); // Pt for theSecondHadr G4double smim = sPart->MinMassOfFragm(); // MinMassLim for SecondHadron G4QParticle* tPart=theWorld->GetQParticle(cV[2]); // Pt for the Third Hadron G4double tmim = tPart->MinMassOfFragm(); // MinMassLimit for theThirdHd G4double fdm = m - smim - tmim; // MaxMassLimit for theFirstHd G4QParticle* fPart=theWorld->GetQParticle(cV[0]); // Pt for the First Hadron fHadr = new G4QHadron(fPart,fdm); // the First Hadron is created if(fPDG<0) fHadr->MakeAntiHadron(); m-=fHadr->GetMass(); // Reduce the residual MaxMass G4double sdm = m - tmim; // MaxMassLimit for theSecondH sHadr = new G4QHadron(sPart,sdm); // theSecondHadron is created if(sPDG<0) sHadr->MakeAntiHadron(); G4double tdm = m - sHadr->GetMass(); // MaxMassLimit for theThird H tHadr = new G4QHadron(tPart,tdm); // the Third Hadron is created if(tPDG<0) tHadr->MakeAntiHadron(); } #ifdef debug G4cout<<"G4Quasmon::DecayQHadron:3Dec. m1="<GetMass() <<",m2="<GetMass()<<",m3="<GetMass()<MakeAntiHadron(); sHadr->MakeAntiHadron(); tHadr->MakeAntiHadron(); } G4LorentzVector f4Mom = fHadr->Get4Momentum(); // Get 4M of the First Hadron (mass) G4LorentzVector s4Mom = sHadr->Get4Momentum(); // Get 4M of the SecondHadron (mass) G4LorentzVector t4Mom = tHadr->Get4Momentum(); // Get 4M of the Third Hadron (mass) if(!qH->DecayIn3(f4Mom,s4Mom,t4Mom)) { delete fHadr; // Delete "new fHadr" delete sHadr; // Delete "new sHadr" delete tHadr; // Delete "new tHadr" G4cerr<<"---Warning---G4Q::DecayQHadron:in3,PDGC="<push_back(qH); // Fill as it is (delete equivalent) return theFragments; } else { //qH->SetNFragments(3); //theFragments.push_back(q); // Fill with NFr=3 (del.equiv.) // Instead delete qH; // fHadr->Set4Momentum(f4Mom); // Put the randomized 4Mom to 1-st Hadron G4QHadronVector* theTmpQHV=DecayQHadron(fHadr); // Try to decay G4int nProd=theTmpQHV->size(); #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH1="<push_back((*theTmpQHV)[0]);// Final = no Further Decay else for(G4int ip1=0; ip1size(); if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);// Final = no Further Decay else { theFragments->resize(tmpS+theFragments->size());// Resize theFragments length copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS); } #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec11 nProd="<clear(); delete intTmpQHV; } theTmpQHV->clear(); delete theTmpQHV; sHadr->Set4Momentum(s4Mom); // Put the randomized 4Mom to 2-nd Hadron theTmpQHV=DecayQHadron(sHadr); // Try to decay nProd=theTmpQHV->size(); #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH2="<push_back((*theTmpQHV)[0]);// Final = no Further Decay else for(G4int ip1=0; ip1size(); if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);// Final = no Further Decay else { theFragments->resize(tmpS+theFragments->size());// Resize theFragments length copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS); } #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec12 nProd="<clear(); delete intTmpQHV; } theTmpQHV->clear(); delete theTmpQHV; tHadr->Set4Momentum(t4Mom); // Put the randomized 4Mom to 3-rd Hadron theTmpQHV=DecayQHadron(tHadr); // Try to decay nProd=theTmpQHV->size(); #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH3="<push_back((*theTmpQHV)[0]);// Final = no Further Decay else for(G4int ip1=0; ip1size(); if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);// Final = no Further Decay else { theFragments->resize(tmpS+theFragments->size());// Resize theFragments length copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS); } #ifdef debug G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec13 nProd="<clear(); delete intTmpQHV; } theTmpQHV->clear(); delete theTmpQHV; } #ifdef debug G4cout<<"G4Q::DecQHadr: DecayIn3 is made with nH="<size()<push_back(qH); // Fill hadron as it is (del.equivalent) } } else { #ifdef debug G4cout<<"G4Quas::DecQHadr:Fill PDG= "<>>>"<push_back(qH); // Fill as it is (delete equivalent) } #ifdef debug G4cout<<"G4Q::DecQHadr:====HADRON IS DECAYED==== with nH="<size()<