// // ******************************************************************** // * 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: G4QEnvironment.cc,v 1.171 2010/06/25 14:03:44 mkossov Exp $ // GEANT4 tag $Name: hadr-chips-V09-03-08 $ // // ---------------- G4QEnvironment ---------------- // by Mikhail Kossov, August 2000. // class for Multy Quasmon Environment used by the CHIPS Model // ------------------------------------------------------------ // Short description: The G4QEnvironment class corresponds to the nuclear // environment, containing excited nucleons or nuclear clusters // (G4Quasmons). In the constructer the nucleus (G4QNucleus) is clusterized // and then the projectile hadron (or hadrons) can create one (or a few) // Quasmons, which are fragmented by the Fragment member function. As a // result a vector of G4QHadrons is created, which can include the residual // nucleus in the ground state. //--------------------------------------------------------------------- //#define debug //#define pdebug //#define qdebug //#define chdebug //#define sdebug //#define ppdebug //#define cdebug //#define cldebug //#define edebug //#define rdebug //#define fdebug //#define ffdebug //#define pcdebug //#define mudebug #include "G4QEnvironment.hh" #include #include using namespace std; G4QEnvironment::G4QEnvironment(const G4QNucleus theEnv) : theEnvironment(theEnv) { G4int envPDG=theEnv.GetPDG(); G4QPDGCode envQPDG(envPDG); G4int envA=envQPDG.GetBaryNum(); G4double envM=envQPDG.GetMass(); totCharge=envQPDG.GetCharge(); totBaryoN=envA; tot4Mom=G4LorentzVector(0.,0.,0.,envM); #ifdef debug G4cout<<"G4QEnviron::Const: t4M="<>>>>>G4QE::Const: Called targPDG="<GetQC()<Get4Momentum()<GetPDGCode() == 10) // Chipolino is found in the input -> Decay { G4QContent chQC=curQH->GetQC(); // Quark content of the Hadron-Chipolino G4QChipolino QCh(chQC); // Define a Chipolino instance for the Hadron G4LorentzVector ch4M=curQH->Get4Momentum(); // 4Mom of the Hadron-Chipolino G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron G4double h1M =h1QPDG.GetMass();// Mass of the first hadron G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron G4double h2M =h2QPDG.GetMass();// Mass of the second hadron G4double chM2 =ch4M.m2(); // Squared Mass of the Chipolino if( sqr(h1M+h2M) < chM2 ) // Decay is possible { G4LorentzVector h14M(0.,0.,0.,h1M); G4LorentzVector h24M(0.,0.,0.,h2M); if(!G4QHadron(ch4M).DecayIn2(h14M,h24M)) { G4cerr<<"***G4QE::Constr: CM="< h1="< H1="< H2="<GetQC()<Get4Momentum() <<", chipoM="< tC="<GetQPEntries(); // A#of init'ed particles in CHIPS World G4int nCl=nP-90; // A#of init'ed clusters in CHIPS World #ifdef debug G4cout<<"G4QEnv:Const:Before NCI:n="<GetNFragments()<<",tC=" < mPi0 ?"<999.&&!exEnviron.SplitBaryon())//Nucleus is below SplitFragmentThreshold { #ifdef debug G4cout<<"G4QEnv::Const:Photon's added to Output, Env="<Output, Env="<Get4Momentum(); //G4double qpen=-180.*log(G4UniformRand()); // Energy of target-quark-parton(T=180) G4double qpen=465.*sqrt(sqrt(G4UniformRand())); // UniformDistr for 3-q nucleon G4double qpct=2*G4UniformRand()-1.; // Cos(thet) of target-quark-parton G4double qpst=sqrt(1.-qpct*qpct); // Sin(theta) of target-quark-parton G4double qppt=qpen*qpst; // PT of target-quark-parton G4double qphi=twopi*G4UniformRand(); // Phi of target-quark-parton G4LorentzVector qi4m(qppt*sin(qphi),qppt*cos(qphi),qpen*qpct,qpen); // quark-parton G4LorentzVector qt4m=mu4m+qi4m; // Total 4mom (iniQP+lepton) G4LorentzVector nu4m(0.,0.,0.,0.); // Prototype of secondary neutrino 4mom G4LorentzVector qf4m(0.,0.,0.,0.); // Prototype for secondary quark-parton G4QContent targQC=targQPDG.GetQuarkContent(); // QC of the target nucleus (local!) targQC+=G4QContent(1,0,0,0,1,0); // Make iso-shift with fake pi- is added G4LorentzVector fn4m=G4LorentzVector(0.,0.,0.,0.); // Prototype of the residual 4M G4QNucleus fnN(targQC,fn4m); // Define the final state nucleus G4double fnm=fnN.GetMZNS(); // GS Mass of the final state nucleus //G4QContent resiQC=targQC-neutQC; // QC of resid nucleus (-neutron) //G4QNucleus rsN(resiQC,fn4m); // Define the final state nucleus //G4double rsm=rsN.GetMZNS()+mNeut; // GS Mass of residual nucleus w/o neutron G4double tm=0.; // Prototype of RealMass of the final nucleus G4LorentzVector tg4m=G4LorentzVector(0.,0.,0.,targM); // 4mom of all target nucleus G4LorentzVector fd4m=tg4m-qi4m; // 4mom of the residual coloured nuclear sys. #ifdef debug //G4cout<<">>>G4QEnv::Const:rM="<n="<nu,q) reaction succeded and Neutrino can be pushed to Output G4QHadron* neutrino = 0; // NeutrinoPrototype to be filled to Output #ifdef mudebug G4cout<<"G4QEnv::Const:fM="<nu"<SetQPDG(pimQPDG); //Convert (mu-)u->d to (virt pi-)u->d capture fake=false; // normal pi- for q-muon scattering //fake=true; // fake pi- for q-muon scattering ***** //if(G4UniformRand()>.5) fake=false; // normal pi- for q-muon scattering ***** opHad->Set4Momentum(fn4m); } } for(G4int ih=0; ih The main LOOP over projQHadrons { G4QHadron* curHadr=projHadrons[ih]; // Pointer to current projectile Hadron G4int hNFrag = curHadr->GetNFragments();// #0 means intermediate (skip) G4LorentzVector ch4M=curHadr->Get4Momentum(); // 4-momenyum of the current projectile #ifdef debug G4cout<<"G4QE:C:"<GetNFragments()<0.) // => "Final hadron" case { G4int envPDG=theEnvironment.GetPDG(); if(envPDG==90000000||(theEnvironment.Get4Momentum().m2())<1.) // ==>"Vacuum" { G4int hPDG = curHadr->GetPDGCode();// A PDG Code of the projQHadron //if(!hPDG||hPDG==10) // Check for the validity of the QHadron (@@ 10 OK?) if(!hPDG) // Check for the validity of the QHadron { //G4cerr<<"--Warning--G4QEnvironment::Constructor: wrong PDG("<GetQCode(); // One more check for valid of the QHadron if(hQ<0) { //G4cerr<<"--Warning--G4QEnv::Constructor:Q<0, PDG=("<GetPDGCode(); // A PDG Code of the projQHadron if(!hPDG||hPDG==10) // Check for the validity of the QHadron { G4cout<<"---Warning---G4QEnvironment::Constructor:Vacuum,1st Hadron wrong PDG="<GetQCode(); // One more check for valid of the QHadron if(hQ<0) { G4cout<<"---Warning---G4QEnviron::Constructor:Vacuum,Q<0, 1st HPDG="<Get4Momentum(); G4QContent hQC = curHadr->GetQC(); if(!targPDG||targPDG==10) G4cout<<"G4QEnv::CreateQ; (1) PDG="<GetQC()<Get4Momentum()<GetCharge(); finBaryoN+=theQHadrons[ih]->GetBaryonNumber(); } G4int nQuas=theQuasmons.size(); if(nQuas) for(G4int iq=0; iqGetCharge(); finBaryoN+=theQuasmons[iq]->GetBaryonNumber(); } if(finCharge!=totCharge || finBaryoN!=totBaryoN) { G4cout<<"***G4QEnv::C:tC="<GetQC()<Get4Momentum()<theQHadrons.size(); if(nQH) for(G4int ih=0; ihtheQHadrons[ih]); #ifdef debug G4cout<<"G4QE::CopyByPtr:cH#"<GetQC()<Get4Momentum()<theProjectiles.size(); if(nP) for(G4int ip=0; iptheProjectiles[ip]); theProjectiles.push_back(curP); // (delete equivalent) } theTargetPDG = right->theTargetPDG; theWorld = right->theWorld; nBarClust = right->nBarClust; f2all = right->f2all; tot4Mom = right->tot4Mom; totCharge = right->totCharge; totBaryoN = right->totBaryoN; // theQuasmons (Vector) G4int nQ = right->theQuasmons.size(); if(nQ) for(G4int iq=0; iqtheQuasmons[iq]); theQuasmons.push_back(curQ); // (delete equivalent) } // theQCandidates (Vector) G4int nQC = right->theQCandidates.size(); if(nQC) for(G4int ic=0; ictheQCandidates[ic]); theQCandidates.push_back(curQC); // (delete equivalent) } theEnvironment = right->theEnvironment; } G4QEnvironment::~G4QEnvironment() { #ifdef debug G4cout<<"~G4QEnvironment: before theQCandidates nC="<GetQC()<Get4Momentum()<>AnnihilationIsDone,C="<size(); // For the selection LOOP ^ ^ G4ThreeVector dir = G4RandomDirection(); // For the selection in LOOP (@@ at rest)^ ^ G4double ra=std::pow(G4double(totBaryoN),third); // ^ ^ #ifdef debug G4cout<<"G4QE::CQ:N="<operator[](ind); // Pointer to theCurrentHadron ^ ^ G4QHadron* curHadr = (*output)[ind]; // Pointer to theCurrentHadron ^ ^ G4int shDFL= curHadr->GetNFragments();// A#of decFragments for proj. ^ ^ G4LorentzVector sh4m = curHadr->Get4Momentum(); // 4Mom for the projectile ^ ^ G4ThreeVector shDIR= sh4m.vect().unit(); // unitVector in projMomDirect ^ ^ G4int shPDG= curHadr->GetPDGCode(); // PDG Code of the projectile ^ ^ G4int shCHG= curHadr->GetCharge(); // Charge of the projectile ^ ^ G4double shMOM= sh4m.rho(); // Momentum of the projectile ^ ^ #ifdef debug G4cout<<"G4QE::CrQ:"< no absorption ^ ^ else if(shMOM<.1) // Meson at rest ^ ^ { // ^ ^ if(shCHG<=0.) solAnCut=-3.; // Always totally absorbed ^ ^ else solAnCut=3.; // Positive are repelled from A ^ ^ } // ^ ^ else solAnCut+=1000*shCHG/shMOM/ra; // ChargeDepSolAngle(Normal) ^ ^ //G4double solAnCut=SolidAngle+20*shCHG*sqrt(1.*envZ)/shMOM;//ChargeDepSolAngle ^ ^ #ifdef debug G4cout<<"G4QE::CrQ: PDG="<H="<"<solAnCut||shMOM<120.) && abs(shPDG)>99) // Absorb mesons ^ ^ if(dir.dot(shDIR)>solAnCut && abs(shPDG)>99) // Absorb mesons ^ ^ { #ifdef debug G4cout<<"G4QE::CQ:>H="<"< Case of "Energy Flux approach" ^ ^ { G4QContent shQC = curHadr->GetQC();// QuarkContent of the Current Hadron ^ ^ ef4Mom+=sh4m; EnFlQC+=shQC; efCounter++; #ifdef debug G4int hPDG=curHadr->GetPDGCode(); // Only for gebug printing ^ ^ G4LorentzVector h4M = curHadr->Get4Momentum(); // Only for gebug printing^ ^ G4cout<<"G4QE::CrQ:#"<GetPDGCode(); // Only for debug printing ^ ^ G4LorentzVector h4M = curHadr->Get4Momentum(); // Only for gebug printing ^ ^ G4cout<<"G4QE::CrQ:Absorb#"<begin(), output->end(), DeleteQHadron()); //DESTRUCT output>-^-^ output->clear(); // ^ ^ delete output; // ==================================^=^ if(!efFlag) // =>NotEnergyFlux=MultyQuasmon Case ^ { G4int noh = theQHadrons.size(); // a#oh hadrons in Output UpToNow ^ if(noh) for(G4int kh=0; khGet4Momentum()< m="< M="<.001)G4cout<<"-W-G4QE::CQA:M="<.1) G4cout<<"-W-G4QE::CrQ:m2="<3) {allB=AP; allB=pPDG;} // A trick to use not used AP(3M) & pPDG for (unsigned index=0; indexGetPDGCode(); // PDG Code of the Candidate G4int cST = curCand->GetStrangeness(); // Strangeness of the candidate G4int cBN = curCand->GetBaryonNumber(); // Baryon Number of the candidate G4int cCH = curCand->GetCharge(); // Charge of the candidate #ifdef sdebug G4cout<<"G4QE::PIP:=====> #"<GetPreProbability(); // A number of clusters of the type G4double dOfCl=curCand->GetDenseProbability();// A number of clusters in dense region #ifdef sdebug G4cout<<"G4QE::PIP:Z="<GetQC(); // Quark Content of the candidate ////////////G4int pC = projQC.GetCharge(); // Charge of the projectile G4QContent qQC=pQC+projQC; // Total Quark content of the Compound G4QPDGCode qQPDG(qQC); G4int qC = qQPDG.GetQCode(); G4double d = abs(zc-nc); G4double fact=1./pow(2.,d); if (qC<-1) probab=0.; //else if(pPDG==-211&&AP<152.&&cBN<2) probab=0.; // PionCaptureByCluster //else if(pPDG==22&&AP<152. || pPDG>90000000) else if(pPDG==22&&AP<152.) { if(cBN<2)probab=nOfCl*cBN*fact; //Gamma Under Pi Threshold (QuarkCapture) else probab=0.; } else if(pPDG==2212) { //if(cBN<2)probab=nOfCl*cBN*fact; // Moving nucleons hits only nucleons //else probab=0.; probab=nOfCl*cBN*fact; // Moving nucleons hits composed clusters //probab=nOfCl*fact; // Moving nucleons hits compact clusters } ////////////////////////else if((pPDG==-211&&AP<10.)&&cBN<2) probab=0;//PiCapAtRst(D) //else if((pPDG==-211||pPDG==-13)&&AP<27.)probab=dOfCl*cBN*fact;//Pi/Mu-CaptureAtRest //else if(pPDG==-211&&AP<10.) probab=nOfCl*fact;// special PiCaptureAtRest //else if(pPDG==-211&&AP<10.) probab=nOfCl*cBN*(cBN-1)*fact; //else probab=nOfCl*fact; else probab=nOfCl*cBN*fact; //else probab=dOfCl*cBN*fact; //if(cBN>1) probab=0.; // Suppress clusters //if(cBN>2) probab=0.; // Suppress heavy clusters #ifdef sdebug G4int pPDG = projQC.GetSPDGCode(); // PDG code of the projectile particle G4int rPDG = qQC.GetSPDGCode(); G4double baryn = qQC.GetBaryonNumber(); G4double charge= qQC.GetCharge(); G4double dq= abs(baryn-charge-charge); G4cout<<"G4QE::PIP:P="< Fill QEnviron { G4int nPDG = theEnvironment.GetPDG(); // PDG code of the residual Nucl.Environ. #ifdef debug G4cout<<"G4QE::HQE:***NO QUASMONS***Env="<80000000) { G4QHadron* rNucleus = new G4QHadron(theEnvironment); // Create HadronEnvironment theQHadrons.push_back(rNucleus); // Fill GS - no further decay (del. equiv.) #ifdef fdebug G4cout<<"G4QEnv::HadrQE: >>>> Fill Environment="< "Environment is Vacuum" case { #ifdef rdebug G4cout<<"G4QEnv::HadrQE: ***Vacuum*** #ofQ="<Get4Momentum(); totIn4M += Q4M; totInC += pQ->GetQC().GetCharge(); } // End of TotInitial4Momentum summation LOOP over Quasmons G4int nsHadr = theQHadrons.size(); // Update the value of OUTPUT entries if(nsHadr) for(G4int jso=0; jsoGetNFragments(); // A#of secondary fragments if(!hsNF) // Add only final hadrons { G4LorentzVector hs4Mom = theQHadrons[jso]->Get4Momentum(); totIn4M += hs4Mom; totInC += theQHadrons[jso]->GetCharge(); } } #endif G4QNucleus vE(90000000); G4int nlq = 0; // Prototype of a#of Living Quasmons if(nQuasmons) for(G4int lq=0; lqGetStatus())nlq++; if(nQuasmons) for(G4int iq=0; iqGetCharge(); f1BaryoN+=theQHadrons[ih]->GetBaryonNumber(); } G4int nQuas=theQuasmons.size(); if(nQuas)for(G4int iqs=0; iqsGetCharge(); f1BaryoN+=theQuasmons[iqs]->GetBaryonNumber(); } if(f1Charge!=totCharge || f1BaryoN!=totBaryoN) { G4cout<<"**G4QE::HQ:q#"<GetQC();// ^ G4int tQBN=totQC.GetBaryonNumber();// Baryon Number of not decayed Quasmon ^ G4QNucleus tqN(totQC); // Define the quasmon as a nucleus ^ G4double gsM=tqN.GetMZNS(); // GS Mass ^ G4LorentzVector tot4M=theQuasmons[iq]->Get4Momentum(); G4double totQM=tot4M.m(); // Real Mass of Quasmon ^ if(tQBN>0&&totQM>gsM) // => "Try Quasmon evaporation" case ^ { // ^ G4QHadron* nuclQ = new G4QHadron(totQC,tot4M); // ^ #ifdef fdebug G4cout<<"G4QEnv::HadrQE:Vac,tQC"<KillQuasmon(); // Kill evaporated Quasmon ^ nlq--; // ^ } else if(iq+11) // => "Try to merge with next" case ^ { G4int s=theQuasmons[iq+1]->GetStatus();//Status of the next Quasmon ^ theQuasmons[iq+1]->IncreaseBy(theQuasmons[iq]);// Merge with the next Quasmon ^ theQuasmons[iq]->KillQuasmon(); // Kill the week Quasmon ^ if(s) nlq--; // Reduce a number of "living Quasmons" ^ } else if(iq+1==nQuasmons&&iq&&nlq>1) // => "Quasmon stack is exhosted" case ^ { G4int s=theQuasmons[0]->GetStatus(); // Status of the first Quasmon ^ theQuasmons[0]->IncreaseBy(theQuasmons[iq]);// Merge with the first Quasmon ^ theQuasmons[iq]->KillQuasmon(); // Kill the week Quasmon ^ if(s) nlq--; // Reduce a number of "living Quasmons" ^ } else // "Have a chance to recover" case ^ { // ^ #ifdef debug G4cout<<"***G4QE::HQE:"<InitQuasmon(totQC,tot4M);// Update the week Quasmon ^ G4QHadronVector* curout=theQuasmons[iq]->Fragment(vE,1);//!DESTROY! <---+ ^ G4int ast=theQuasmons[iq]->GetStatus(); // Status of the Quasmon ^ ^ if(!ast) nlq--; // Reduce nlq if Quasmon decayed ^ ^ G4int nHadrons=curout->size(); // A#of outputQHadrons in theDecayedQ ^ ^ #ifdef debug G4cout<<"G4QEnv::HadrQE:VacuumRecoverQ#"<0) // => "QHadrons from Quasmon to Output"^ ^ { // ^ ^ for (G4int ih=0; ihoperator[](ih)); // ^ ^ G4QHadron* curH = new G4QHadron((*curout)[ih]); // ^ ^ #ifdef debug G4cout<<"G4QEnv::HadrQE:Recovered, H#"<begin(), output->end(), DeleteQHadron()); // >------------------^ output->clear(); // ^ delete output; // >=====================================^ } // End of check for the already decayed Quasmon } // End of the LOOP over Quasmons } else // ==> "Nuclear environment" case { #ifdef rdebug G4cout<<"G4QEnv::HadrQE:FRAGMENTATION IN NUCLEAR ENVIRONMENT nQ="<Get4Momentum(); totIn4M += Q4M; totInC += pQ->GetQC().GetCharge(); } // End of TotInitial4Momentum summation LOOP over Quasmons G4int nsHadr = theQHadrons.size(); // Update the value of #of OUTPUT entries if(nsHadr) for(G4int jso=0; jsoGetNFragments(); // A#of secondary fragments if(!hsNF) // Add only final hadrons { G4LorentzVector hs4Mom = theQHadrons[jso]->Get4Momentum(); totIn4M += hs4Mom; totInC += theQHadrons[jso]->GetCharge(); } } #endif // @@ Experimental calculations G4QContent totInQC=theEnvironment.GetQCZNS(); G4LorentzVector totIn4M=theEnvironment.Get4Momentum(); for (G4int is=0; isGet4Momentum(); totInQC += pQ->GetQC(); } // End of TotInitial4Momentum/QC summation LOOP over Quasmons G4double totMass=totIn4M.m(); G4QNucleus totN(totInQC); G4double totM=totN.GetMZNS(); G4double excE = totMass-totM; // @@ End of experimental calculations G4int envA=theEnvironment.GetA(); //G4int c3Max = 27; // Big number (and any #0) slowes dow a lot //G4int c3Max = 9; // Max#of "no hadrons" steps (reduced below?) G4int c3Max = 3; if(excE > fPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess //G4int c3Max = 1; //if(excE > dPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess //G4int c3Max = 0; //if(excE > mPi0) c3Max=(G4int)(excE/mPi0); // Try more for big excess //G4int c3Max = 0; // It closes the force decay of Quasmon at all! // //G4int premC = 27; //G4int premC = 3; G4int premC = 1; //if(envA>1&&envA<13) premC = 24/envA; if(envA>1&&envA<19) premC = 36/envA; //if(envA>1&&envA<31) premC = 60/envA; //if(envA>1&&envA<61) premC = 120/envA; G4int sumstat= 2; // Sum of statuses of all Quasmons G4bool force = false; // Prototype of the Force Major Flag G4int cbR =0; // Counter of the "Stoped by Coulomb Barrier" // //G4int cbRM =0; //** MaxCounter of "StopedByCoulombBarrier" ** //G4int cbRM =1; // MaxCounter of the "StopedByCoulombBarrier" //G4int cbRM =3; // MaxCounter of the "StopedByCoulombBarrier" //G4int cbRM =9; // MaxCounter of the "Stoped byCoulombBarrier" G4int cbRM = c3Max; // MaxCounter of the "StopedByCoulombBarrier" G4int totC = 0; // Counter to break the "infinit" loop //G4int totCM = 227; // Limit for the "infinit" loop counter G4int totCM = envA; // Limit for the "infinit" loop counter //G4int totCM = 27; // Limit for this counter //G4int nCnMax = 1; // MaxCounterOfHadrFolts for shortCutSolutions //G4int nCnMax = 3; // MaxCounterOfHadrFolts for shortCutSolutions //G4int nCnMax = 9; // MaxCounterOfHadrFolts for shortCutSolutions G4int nCnMax = c3Max; // MaxCounterOfHadrFolts for shortCutSolutions G4bool first=true; // Flag of the first interaction (only NucMedia) G4int cAN=0; // Counter of the nucleon absorptions //G4int mcAN=27; // Max for the counter of nucleon absorptions G4int mcAN=9; // Max for the counter of nucleon absorptions //G4int mcAN=3; // Max for the counter of nucleon absorptions //G4int mcAN=1; // Max for the counter of nucleon absorptions while (sumstat||totC0.) { f2Charge=theEnvironment.GetCharge(); f2BaryoN=theEnvironment.GetA(); } G4int nHad=theQHadrons.size(); if(nHad) for(G4int ih=0; ihGetCharge(); f2BaryoN+=theQHadrons[ih]->GetBaryonNumber(); } G4int nQuas=theQuasmons.size(); if(nQuas)for(G4int iqs=0; iqsGetCharge(); f2BaryoN+=theQuasmons[iqs]->GetBaryonNumber(); } if(f2Charge!=totCharge || f2BaryoN!=totBaryoN) { G4cout<<"**G4QE::HQ:(NucEnv)i#"<"Force Quasmons decay" { if(!fCount) premC--; // Reduce premium efforts counter if(nQuasmons) for (G4int jq=0; jqGetStatus();// Old status of the Quasmon ^ #ifdef debug G4cout<<"G4QE::HQE:Status of Q#"<operator[](kpo);// Pointer to theOutputQHadron ^ ^ G4QHadron* insH = (*output)[kpo];// Pointer to the Output QHadron ^ ^ G4int qhsNF = insH->GetNFragments(); // A#of secondary fragments ^ ^ if(!qhsNF) // Subtract only final hadrons ^ ^ { // ^ ^ G4LorentzVector qhs4Mom = insH->Get4Momentum();// 4M of theOutputQHadron^ ^ G4int hPDG = insH->GetPDGCode(); // PDG Code of the Q-output Hadron ^ ^ G4cout<<"G4QE::HQE:SUM-4-Mom qh("<GetCharge(); // ^ ^ } // ^ ^ } // End of the LOOP over output QHadrons ^ ^ G4cout<<"G4QEnv::HadrQE:|||||4-MomCHECK||||d4M="<0) // Transfer QHadrons from Quasm to Output ^ ^ { // ^ ^ for (G4int ih=0; ihoperator[](ih); // ^ ^ G4QHadron* inpH = (*output)[ih]; // ^ ^ G4int hC=inpH->GetCharge(); // Charge of the Hadron ^ ^ G4int hF=inpH->GetNFragments();// Number of fragments ^ ^ G4double hCB=0.; // Coulomb Barrier ^ ^ G4double hKE=0.; // Kinetic Energy of the Hadron ^ ^ G4LorentzVector hLV=inpH->Get4Momentum(); // ^ ^ #ifdef debug G4cout<<"G4QEnv::HadrQE:H#"< "Suck into existing Quasmon" case ^ ^ { // ^ ^ G4QContent tQC=inpH->GetQC()+pQ->GetQC(); // ^ ^ G4LorentzVector tLV=hLV+pQ->Get4Momentum();// ^ ^ pQ->InitQuasmon(tQC,tLV); // Reinitialize the current Quasmon ^ ^ #ifdef debug G4cout<<"G4QE::HQE:Medium, H#"< "Suck in the existing Quasmon" case ^ ^ { // ^ ^ G4QContent tQC=inpH->GetQC()+theEnvironment.GetQCZNS();// ^ ^ G4LorentzVector tLV=hLV+theEnvironment.Get4Momentum(); // ^ ^ theEnvironment=G4QNucleus(tQC,tLV); // Reinit currentEnvironment ^ ^ #ifdef debug G4cout<<"G4QE::HQE:Med,H#"<GetQPDG()<<",4M="//^ ^ <Get4Momentum()<<" is suckedInEnvironment"< "Hadron can go out" case ^ ^ { // ^ ^ G4QHadron* curH = new G4QHadron(inpH); // ^ ^ #ifdef debug G4LorentzVector ph4M=curH->Get4Momentum(); // 4-mom of the hadron ^ ^ G4double phX=ph4M.x(); // p_x of the hadron ^ ^ G4double phY=ph4M.y(); // p_y of the hadron ^ ^ G4double phZ=ph4M.z(); // p_x of the hadron ^ ^ G4double phCost=phZ/sqrt(phX*phX+phY*phY+phZ*phZ);// Hadr cos(theta)^ ^ G4cout<<"G4QEnv::HadrQE:Medium, H#"<GetQPDG()//^ ^ <<", 4M="<500.)// Add N from E to Q^ { // ^ #ifdef debug G4cout<<"G4QE::HQE:E="< envZ ) // Change to a proton ^ { // ^ resPDG=envPDG-1000; // new PDG for the Environment ^ nucQC=protQC; // proton QContent ^ } // ^ #ifdef debug G4cout<<"G4QE::HQE:P,eZ="< 0.001) // the Environment is not at rest ^ { // ^ res4M=(resM/eM)*env4M; // Proportional 4M for residEnv ^ nuc4M=(nucM/eM)*env4M; // Proportional 4M for effNucleon ^ } // ^ theEnvironment=G4QNucleus(res4M,resPDG);// Update the Environment ^ theQuasmons[0]->IncreaseBy(nucQC,nuc4M);// Update the Only Quasmon ^ #ifdef debug G4cout<<"G4QE::HQE:P,Q="<nCnMax)// Treat PANIC for stat=2 (NothingWasDone) ^ { // ^ #ifdef debug G4cout<<"G4QE::HQE:PANIC,nC="<"<GetQC(); // QuarkContent of the Quasmon ^ G4int pqC=qQC.GetCharge(); // Charge (nP) of the Current Quasmon ^ G4int pqS=qQC.GetStrangeness(); // Strangeness (nL) of theCurrQuasmon^ G4int pqB=qQC.GetBaryonNumber(); // BaryNumber of the CurrentQuasmon ^ G4LorentzVector cq4M=pQ->Get4Momentum(); // 4Mom of the Current Quasmon ^ G4double cqMass=cq4M.m(); // Real Mass of the current Quasmon ^ G4double fqMass=G4QPDGCode(22).GetNuclMass(pqC,pqB-pqC-pqS,pqS);//CQ FreeM^ #ifdef edebug G4cout<<"G4QEnv::HQE:M="<fM="<KillQuasmon(); // If BackFusion succeeded, kill theQuasmon ^ #ifdef edebug G4cout<<"G4QEnv::HQE:Status after kill (#"<GetStatus()// ^ <<", nH="< "NuclearEnvironment" case ^ { // ^ if(eCount>1) // ^ { // ^ #ifdef fdebug G4cout<<"G4QE::HQE:TOTEVAP tPDG="<begin(), output->end(), DeleteQHadron());// >--------^ output->clear(); // ^ delete output; // >==========>===========>=================^ return theQHadrons; // ^ } // ^ G4LorentzVector t4M=cq4M+theEnvironment.Get4Momentum(); // Q+E tot4Mom ^ G4double tM=t4M.m(); // Real total (Quasmon+Environment) mass ^ G4QContent envQC=theEnvironment.GetQCZNS(); // QuarkCont of NucEnviron ^ G4QContent curQC=envQC+qQC; // Total Quark Content ^ G4QNucleus curE(curQC); // Pseudo nucleus for the Total System ^ G4double curM=curE.GetGSMass();// min mass of the Total System ^ #ifdef edebug G4cout<<"G4QEnv::HQE:Q#"<gsM="< "Quasmon-Chipolino or Environment-Dibaryon" case ^ if(qPDG==10 || qPDG==92000000 || qPDG==90002000 || qPDG==90000002)//^ { // ^ G4QPDGCode h1QPDG=nQPDG; // QPDG of the first hadron ^ G4double h1M =mNeut; // Mass of the first hadron ^ G4QPDGCode h2QPDG=h1QPDG;// QPDG of the second hadron ^ G4double h2M =mNeut; // Mass of the second hadron ^ if(qPDG==10) // CHIPOLINO decay case ^ { // ^ G4QChipolino QChip(qQC);// define the Quasmon as a Chipolino ^ h1QPDG=QChip.GetQPDG1();// QPDG of the first hadron ^ h1M =h1QPDG.GetMass();// Mass of the first hadron ^ h2QPDG=QChip.GetQPDG2();// QPDG of the second hadron ^ h2M =h2QPDG.GetMass();// Mass of the second hadron ^ } // ^ else if(qPDG==90002000) // DiProton decay case ^ { // ^ h1QPDG=pQPDG; // QPDG of the first hadron ^ h1M =mProt; // Mass of the first hadron ^ h2QPDG=h1QPDG; // QPDG of the second hadron ^ h2M =mProt; // Mass of the second hadron ^ } // ^ else if(qPDG==92000000) // Two lambdas case ^ { // ^ h1QPDG=lQPDG; // QPDG of the first hadron ^ h1M =mLamb; // Mass of the first hadron ^ h2QPDG=h1QPDG; // QPDG of the second hadron ^ h2M =mLamb; // Mass of the second hadron ^ G4double ddMass=totMass-envM; // Free CM energy ^ if(ddMass>mSigZ+mSigZ) // Sigma0+Sigma0 is possible ^ { // @@ Only two particles PS is used ^ G4double dd2=ddMass*ddMass; // Squared free energy ^ G4double sma=mLamb+mLamb; // Lambda+Lambda sum ^ G4double pr1=0.; // Prototype to avoid sqrt(-) ^ if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Lamb+Lamb PS ^ sma=mLamb+mSigZ; // Lambda+Sigma0 sum ^ G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference ^ G4double pr2=pr1; // Prototype of +L+S0 PS ^ if(ddMass>sma&&ddMass>smi) // ^ pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); // ^ sma=mSigZ+mSigZ; // Sigma0+Sigma0 sum ^ G4double pr3=pr2; // Prototype of +Sigma0+Sigma0 PS ^ if(ddMass>sma) pr3+=sqrt((dd2-sma*sma)*dd2); // ^ G4double hhRND=pr3*G4UniformRand(); // Randomize PS ^ if(hhRND>pr2) // --> "ENnv+Sigma0+Sigma0" case ^ { // ^ h1QPDG=s0QPDG; // QPDG of the first hadron ^ h1M =mSigZ; // Mass of the first hadron ^ h2QPDG=h1QPDG; // QPDG of the second hadron ^ h2M =mSigZ; // Mass of the second hadron ^ } // ^ else if(hhRND>pr1)// --> "ENnv+Sigma0+Lambda" case ^ { // ^ h1QPDG=s0QPDG; // QPDG of the first hadron ^ h1M =mSigZ; // Mass of the first hadron ^ } // ^ } // ^ else if(ddMass>mSigZ+mLamb) // Lambda+Sigma0 is possible ^ { // @@ Only two particles PS is used ^ G4double dd2=ddMass*ddMass; // Squared free energy ^ G4double sma=mLamb+mLamb; // Lambda+Lambda sum ^ G4double pr1=0.; // Prototype to avoid sqrt(-) ^ if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Lamb+Lamb PS ^ sma=mLamb+mSigZ; // Lambda+Sigma0 sum ^ G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference ^ G4double pr2=pr1; //+L+S0 PS ^ if(ddMass>sma && ddMass>smi) // ^ pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); // ^ if(pr2*G4UniformRand()>pr1) // --> "ENnv+Sigma0+Lambda" case ^ { // ^ h1QPDG=s0QPDG; // QPDG of the first hadron ^ h1M =mSigZ; // Mass of the first hadron ^ } // ^ } // ^ } // ^ if(h1M+h2M+envM "Three parts decay" case ^ { // ^ G4LorentzVector h14M(0.,0.,0.,h1M); // ^ G4LorentzVector h24M(0.,0.,0.,h2M); // ^ G4LorentzVector e4M(0.,0.,0.,envM); // ^ if(!G4QHadron(tot4M).DecayIn3(h14M,h24M,e4M)) // ^ { // ^ G4cerr<<"***G4QE::HQE:(0)tM="<h1="< totMass");//^ #endif CleanUp(); // ^ G4QHadron* evH = new G4QHadron(totQC,tot4M);// ResidualNuclHadr ^ EvaporateResidual(evH); // Evaporate residual (del. equiv.) ^ return theQHadrons; // ^ } // ^ } // ^ else // => "Two particles decay" case ^ { // ^ G4LorentzVector fq4M(0.,0.,0.,qMass); // ^ G4LorentzVector qe4M(0.,0.,0.,envM); // ^ if(!G4QHadron(tot4M).RelDecayIn2(fq4M,qe4M,cq4M,1.,1.))//Q ch.dir.^ { // ^ G4cerr<<"***G4QEnv::HadQE:(0)tM="< qPDG="<begin(),output->end(),DeleteQHadron()); //->->-^ output->clear(); // -->-->-->-->-->-->-->-->-->-->-->-->-->--^ delete output; // >========================================^ throw G4QException("***G4QEnv::HadrQEnv: Q+Env DecIn2 error");//^ } // ^ G4QHadron* qH = new G4QHadron(qPDG,fq4M);// the out going Quasmon ^ theQHadrons.push_back(qH); // (delete equivalent) ^ #ifdef debug G4cout<<"G4QE::HQE:QuasmH="<h1="< H1="<H2="<h1="<GetQPDG(); // QPDG Code for the Quasmon ^ G4int PDGQ=QPDGQ.GetPDGCode(); // PDG Code of the QUASMON ^ #ifdef edebug G4cout<<"G4QEnv::HadrQEnv: vacuum PDGQ="<mPi+mLamb) // "Lambda + Pi- is possible" case ^ { // ^ hyM = mLamb; // Lambda mass ^ hyPDG = 3122; // Lambda PDG Code ^ } // ^ else if(cqMass>mSigM) // "Sigma- gamma decay" case ^ { // ^ hyM=mSigM; // Sigma- mass ^ hyPDG=3112; // Sigma- PDG Code ^ pigM=0.; // Gamma mass ^ pigPDG=22; // Gamma PDG Code ^ } // ^ } // ^ else if(PDGQ==3222||PDGQ==91000999) // --> "Sigma+" case ^ { // ^ pigPDG= 211; // Pi+ PDG Code ^ if(cqMass>mPi+mLamb) // --- "Lambda + Pi+ is possible" case ^ { // ^ hyM = mLamb; // Lambda mass ^ hyPDG = 3122; // Lambda PDG Code ^ pigM = mPi; // Pi+ mass ^ pigPDG= 211; // Pi+ PDG Code ^ } // ^ else if(cqMass>mSigP) // "Sigma- gamma decay" case ^ { // ^ hyM=mSigP; // Sigma+ mass ^ hyPDG=3222; // Sigma+ PDG Code ^ pigM=0.; // Gamma mass ^ pigPDG=22; // Gamma PDG Code ^ } // ^ else if(cqMass>mPi0+mProt&&G4UniformRand()>.5) // "P + Pi0" case ^ { // ^ hyM = mProt; // Proton mass ^ hyPDG = 2212; // Proton PDG Code ^ pigM = mPi0; // Pi0 mass ^ pigPDG= 111; // Pi0 PDG Code ^ } // ^ else if(cqMass"< Hyperon="// ^ < "DiNeutron+Pi-" case ^ { // ^ dinFlag = true; // For the final decay ^ hyM=mNeut+mNeut; // Di-neutron ^ hyPDG=2112; // Neutron PDG Code ^ pigM=mPi; // Pi- mass ^ pigPDG=-211; // Pi- PDG Code ^ } // ^ } // ^ else if(PDGQ==91001999) // --> "p+Sigma+" case ^ { // ^ hyM=mSigP; // Sigma+ ^ hyPDG=3222; // PDG Code of Sigma+ ^ pigM=mProt; // Proton mass ^ pigPDG=2212; // PDG Code of proton ^ if(cqMass "Proton+Proton" case ^ { // ^ hyM=mProt; // Proton mass ^ hyPDG=2212; // Proton PDG Code ^ pigM=mProt; // Proton mass ^ pigPDG=2212; // Proton PDG Code ^ } // ^ } // ^ G4LorentzVector b4Mom(0.,0.,0.,hyM); // Hyperon (di-nucleon) mass ^ G4LorentzVector m4Mom(0.,0.,0.,pigM);// Nucleon (pion) mass ^ if(!G4QHadron(cq4M).DecayIn2(b4Mom, m4Mom)) // "DecayIn2 failed" case ^ { // ^ G4cout<<"--Warning--G4QE::HQE:S/D="<"< Sigma/dN="//^ < "DiNeutron+Pi-" case ^ { // ^ dinFlag = true; // For the final decay ^ pigM=mNeut+mNeut+mNeut;// Tri-neutron ^ pigPDG=2112; // Neutron PDG Code ^ hyM=mPi; // Pi- mass ^ hyPDG=-211; // Pi- PDG Code ^ } // ^ } // ^ else if(PDGQ==91002999) // --> "p+Sigma+" case ^ { // ^ hyM=mSigP; // Sigma+ ^ hyPDG=3222; // PDG Code of Sigma+ ^ pigM=mProt+mProt; // Di-Proton mass ^ pigPDG=2212; // PDG Code of proton ^ if(cqMass "DiProton+Pi+" case ^ { // ^ dinFlag = true; // For the final decay ^ pigM=mProt+mProt+mProt;// Tri-proton ^ pigPDG=2212; // Neutron PDG Code ^ hyM=mPi; // Pi+ mass ^ hyPDG=211; // Pi+ PDG Code ^ } // ^ } // ^ G4LorentzVector b4Mom(0.,0.,0.,hyM); // Hyperon (di-nucleon) mass ^ G4LorentzVector m4Mom(0.,0.,0.,pigM);// Nucleon (pion) mass ^ if(!G4QHadron(cq4M).DecayIn2(b4Mom, m4Mom)) // "DecayIn2 failed" case ^ { // ^ G4cout<<"--Warning--G4QE::HQE:S/Pi="<"< Sigma/Pi="//^ <gsM="<>** CGS Correction **>>**"<begin(), output->end(), DeleteQHadron()); //-->-->-^ output->clear(); // ^ delete output; // >========================================^ pQ->KillQuasmon(); // If BackFusion -> kill theQuasmon ^ eCount--; // Reduce a#of the living Quasmons ^ return theQHadrons; // ^ } // ^ else if(qMGetQC().GetSPDGCode()==1114 // ^ || pQ->GetQC().GetSPDGCode()==2224) // ^ &&qM>theWorld->GetQParticle(QPDGQ)->MinMassOfFragm())//Del&Kill ^ { // ^ #ifdef edebug G4cout<<"G4QEnv::HadrQEnv:**||** Copy&Decay **||**"<DecayQuasmon();//Dec.Quasm & fill decHV=*^ CopyAndDeleteHadronVector(decHV);// Copy output to QHadrV of G4Env ^ tot4M-=pQ->Get4Momentum(); // tot4M recalculation ^ totQC-=pQ->GetQC(); // totQC recalculation ^ pQ->KillQuasmon(); // Make done the current Quasmon ^ eCount--; // Reduce a#of the living Quasmons ^ } // ^ } // ^ } // End of the Medium/Vacuum IF ^ } // End of the status ELSE IF ^ else if(status==3) count3++; // ^ if(status<0) // Panic: Quasmon is below theMassShell ^ { // ^ //if(eCount==1 && DecayInEnvQ(pQ)) // ^ //{ // ^ // for_each(output->begin(), output->end(), DeleteQHadron());//-->-->-->-^ // output->clear(); // ^ // delete output; // >======================================^ // eCount--; // Reduce a#of the living Quasmons ^ // pQ->KillQuasmon(); // ^ // return theQHadrons; // ^ //} // ^ G4int ppm=jq; // Initialized by PANIC Quasmon pointer ^ G4int nRQ=0; // Prot. of a#of additionalRealQuasmons ^ #ifdef edebug G4cout<<"G4QEnv::HadrQEnv: ***PANIC*** for jq="<Get4Momentum().vect(); // PANICQuasmon momentum ^ G4double dpm=1.e+30; // Big number (dot product of momenta) ^ if(nQuasmons>1) for(G4int ir=0; irGetStatus();// Status of a Quasmon ^ #ifdef edebug G4cout<<"G4QEnv::HadrQEnv: ir="<0) // Skip the dead Quasmon ^ { nRQ++; // Increment real-Quasmon-counter ^ G4double dp=vp.dot(rQ->Get4Momentum().vect()); // ^ if(dpGetQC(); // ^ G4LorentzVector r4M= rQ->Get4Momentum(); // ^ rQC += pQ->GetQC(); // ^ r4M += pQ->Get4Momentum(); // ^ rQ->InitQuasmon(rQC, r4M); // Make new Quasmon ^ #ifdef edebug G4cout<<"G4QE::HQE:"<GetQC()<<"+"<GetQC()<<"="<KillQuasmon(); // Delete old Quasmon ^ eCount--; // Decrement counter of living Quasmons ^ } // ^ else // No candidate to resolve PANIC was found ^ { // ^ #ifdef edebug G4cout<<"G4QEnv::HQE: No Q-cand. nRQ="<begin(), output->end(), DeleteQHadron()); //-->-->---^ output->clear(); // ^ delete output; // >======================================^ pQ->KillQuasmon(); // ^ eCount--; // Reduce a#of the living Quasmons ^ return theQHadrons; // ^ } // ^ #ifdef fdebug G4cout<<"G4QEnv::HadrQEnv:NO PANICsolution,t="<GetStatus();// Status of a Quasmon ^ if(jr==jq) // ^ { // ^ totQC+=rQ->GetQC(); // QuarkContent of the Quasmon ^ tot4M+=rQ->Get4Momentum();// QuarkContent of the Quasmon ^ } // ^ else if(Qst) // Skip dead Quasmons ^ { // ^ totQC+=rQ->GetQC(); // QuarkContent of the Quasmon ^ tot4M+=rQ->Get4Momentum();// QuarkContent of the Quasmon ^ } // ^ } // End of the "No candidate to resolve PANIC" ELSE ^ pQ->KillQuasmon(); // Kill the only Quasmon ^ eCount--; // Reduce a#of the living Quasmons ^ CleanUp(); // Clean up THIS Quasmon and Environment ^ G4QHadron* evH = new G4QHadron(totQC,tot4M); // Create ResidNuclHadron ^ EvaporateResidual(evH); // Try to evaporate residual (del.equiv.) ^ for_each(output->begin(), output->end(), DeleteQHadron()); //-->-->-->--^ output->clear(); // ^ delete output; // >======================================^ force=true; // Make the force decision ^ break; // Out of the fragmentation loop >->+ ^ } // | ^ } // | ^ } // | ^ for_each(output->begin(), output->end(), DeleteQHadron()); // ->-->-->--|-----^ output->clear(); // | ^ delete output; // >================================|=====^ } // End of skip of the dead Quasmons | #ifdef debug G4cout<<"G4QE::HQE:QStat("<Get4Momentum()<totM+.001) // ==> "Try Evaporate or decay" case { #ifdef edebug G4cout<<"G4QEnv::HadrQE: NQ="<SetNFragments(2); // Put a#of Fragments=2 //theQHadrons.push_back(delta); // Fill the residual DELTA (del.Eq.) // Instead //delete delta; // G4LorentzVector b4Mom(0.,0.,0.,mBar); G4LorentzVector m4Mom(0.,0.,0.,mMes); if(!G4QHadron(tot4M).DecayIn2(b4Mom, m4Mom)) { G4cout<<"---Warning---G4QEnv::HadronizeQE:B="< mDel="<Bar+Mes error"); } delete quasH; return theQHadrons; } #ifdef debug G4cout<<"G4QEnv::HadronizeQEnv: DELTA="< Bar=" < mChipo="<1+2 decay failed"); } delete quasH; return theQHadrons; } #ifdef debug G4cout<<"G4QEnv::HadronizeQEnv: Chipo="< h1=" < mTot="<h="< mTot="<h="<Fragment(vacuum,1);//**!!DESTROY!!** <-<-+ ^ G4int rest = resid->GetStatus(); // New status after fragm attempt ^ ^ if(!rest) eCount--; // Dec ExistingQuasmonsCounter ^ ^ delete resid; //_________________________________^___^ G4int nHadrons = curout->size(); // a#of Hadrons in the outHV ^ if(nHadrons>0) // Transfer QHadrons to Output ^ { for (G4int ih=0; ihclear(); // The instances are filled above ^ delete curout; // >===============================^ return theQHadrons; } } else { G4QContent tQC =totQC; // Not subtracted copy for error prints G4int NSi =0; // a#of additional Sigma G4int SiPDG =0; // PDG of additional Sigma G4double MSi =0.; // TotalMass of additional Sigma G4int NaK =0; // a#of additional Kaons/anti-Kaons G4int aKPDG =0; // PDG of additional Kaons/anti-Kaons G4double MaK =0.; // TotalMass of additionalKaons/anti-Kaons G4int NPi =0; // a#of additional pions G4int PiPDG =0; // PDG of additional pions G4double MPi =0.; // Total Mass of additional pions if (totBN>0&&totS<0&&totChg+totChg>=totBN)// => "additional K+" case { aKPDG=321; NaK=-totS; MaK=mK*NaK; tQC+=totS*KpQC; totChg+=totS; // Charge reduction (totS<0!) totS=0; // Anti-strangness goes to anti-Kaons } else if (totBN>0&&totS<0) // => "additional aK0" case { aKPDG=311; NaK=-totS; MaK=mK0*NaK; tQC+=totS*K0QC; totS=0; // Anti-strangness goes to anti-Kaons } else if (totBN>1&&totS>0&&(totChg<0||totChg>totBN-totS))//=>"additional Sigma" { NSi=totS; // Prototype of a#of Sigmas if(totChg<0) // Negative Sigmas { SiPDG=3112; if(-totChg0&&totS>totBN&&totBN "additional K0" case {// @@ Here Ksi0 check should be added totS=2>totBN=1&&totBN=10&&totS>totBN&&totChg<0)// => "additional K-" case {// @@ Here Ksi- check should be added totS=2>totBN=1&&totChg=-1<0 aKPDG=-321; NaK=totS-totBN; MaK=mK0*NaK; tQC+=NaK*KpQC; totChg+=NaK; // Increase residual charge totS-=NaK; // Reduce residual strangeness } // === Now residual DELTAS should be subtracted === if (totBN>0&&totChg>totBN-totS) // => "additional PI+" case {// @@ Here Sigma+ check should be added totChg=1>totBn=1-totS=1 PiPDG=211; NPi=totChg-totBN+totS; MPi=mPi*NPi; tQC-=NPi*PiQC; totChg-=NPi; } else if (totBN>0&&totChg<0) // => "additional PI-" case {// @@ Here Sigma- check should be added totChg<0 PiPDG=-211; NPi=-totChg; MPi=mPi*NPi; tQC+=NPi*PiQC; // Now anti-Pions must be subtracted totChg+=NPi; } else if (!totBN&&totChg>1-totS) // => "additional PI+" case {// @@ Here Sigma+ check should be added totChg=1>totBn=1-totS=1 PiPDG=211; NPi=totChg+totS-1; MPi=mPi*NPi; tQC-=NPi*PiQC; totChg-=NPi; } else if (!totBN&&totChg<-1-totS) // => "additional PI-" case {// @@ Here Sigma- check should be added totChg<0 PiPDG=-211; NPi-=totChg+totS+1; MPi=mPi*NPi; tQC+=NPi*PiQC; // Now anti-Pions must be subtracted totChg+=NPi; } G4double totRM=0.; // min (GS) Mass of the Residual System if(totBN<2) // Calculate totPDG & totRM { totPDG=tQC.GetSPDGCode(); // MinPDGCode for the Residual compound if(totPDG==10&&tQC.GetBaryonNumber()>0) totPDG=tQC.GetZNSPDGCode(); if(totPDG) totRM=G4QPDGCode(totPDG).GetMass(); // minMass of theResidualSystem else { G4cerr<<"***G4QEnvironment::HadronizeQEnv: totPDG=0"< "Decay in K0 or K+ + NPi" case {//@@ Can (must) be moved to EvaporateResidual ?? if(!NPi) // ==> "Only anti-strange K" case { G4LorentzVector m4Mom(0.,0.,0.,MaK); G4LorentzVector n4Mom(0.,0.,0.,totRM); G4double sum=MaK+totRM; if(fabs(totMass-sum) mSN="< M=" <tM="< nK="< "Only Sigma's" case { G4LorentzVector m4Mom(0.,0.,0.,MSi); G4LorentzVector n4Mom(0.,0.,0.,totRM); G4double sum=MSi+totRM; if(fabs(totMass-sum) mSN="< Sig=" <tM="< nS="< "One isobar" case { G4LorentzVector m4Mom(0.,0.,0.,MPi); G4LorentzVector n4Mom(0.,0.,0.,totRM); if(!G4QHadron(tot4M).DecayIn2(m4Mom, n4Mom)) { G4cout<<"---Warning---G4QEnv::HadronizeQEnv:M="< mSN="< M="<(?)SN="< N*PI="<0,tBN="< "Decay Quasmon or Q+Environ" case { G4int envPDG = theEnvironment.GetPDG();// PDGCode of the NuclQEnvironment #ifdef debug G4cout<<"G4QEnv::HadrQEnv: nQ=1, QPDG=="< "Quasmon-Chipolino or Environment-Dibaryon" case else if(QPDG==10||QPDG==92000000||QPDG==90002000||QPDG==90000002) { G4QContent QQC = pQ->GetQC(); // Quark Content of the Quasmon G4QPDGCode h1QPDG=nQPDG; // QPDG of the first hadron G4double h1M =mNeut; // Mass of the first hadron G4QPDGCode h2QPDG=h1QPDG; // QPDG of the second hadron G4double h2M =mNeut; // Mass of the second hadron if(QPDG==10) { G4QChipolino QChip(QQC); // define the Quasmon as a Chipolino h1QPDG=QChip.GetQPDG1(); // QPDG of the first hadron h1M =h1QPDG.GetMass(); // Mass of the first hadron h2QPDG=QChip.GetQPDG2(); // QPDG of the second hadron h2M =h2QPDG.GetMass(); // Mass of the second hadron } else if(QPDG==90002000) { h1QPDG=pQPDG; // QPDG of the first hadron h1M =mProt; // Mass of the first hadron h2QPDG=h1QPDG; // QPDG of the second hadron h2M =mProt; // Mass of the second hadron } else if(QPDG==92000000) { h1QPDG=lQPDG; // QPDG of the first hadron h1M =mLamb; // Mass of the first hadron h2QPDG=h1QPDG; // QPDG of the second hadron h2M =mLamb; // Mass of the second hadron G4double ddMass=totMass-envM; // Free CM energy if(ddMass>mSigZ+mSigZ) // Sigma0+Sigma0 is possible { // @@ Only two particles PS is used G4double dd2=ddMass*ddMass; // Squared free energy G4double sma=mLamb+mLamb; // Lambda+Lambda sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Lamb+Lamb PS sma=mLamb+mSigZ; // Lambda+Sigma0 sum G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference G4double pr2=pr1; // Prototype of +L+S0 PS if(ddMass>sma && ddMass>smi) pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); sma=mSigZ+mSigZ; // Sigma0+Sigma0 sum G4double pr3=pr2; // Prototype of +Sigma0+Sigma0 PS if(ddMass>sma) pr3+=sqrt((dd2-sma*sma)*dd2); G4double hhRND=pr3*G4UniformRand(); // Randomize PS if(hhRND>pr2) // --> "ENnv+Sigma0+Sigma0" case { // h1QPDG=s0QPDG; // QPDG of the first hadron h1M =mSigZ; // Mass of the first hadron h2QPDG=h1QPDG; // QPDG of the second hadron h2M =mSigZ; // Mass of the second hadron } // else if(hhRND>pr1) // --> "ENnv+Sigma0+Lambda" case { // h1QPDG=s0QPDG; // QPDG of the first hadron h1M =mSigZ; // Mass of the first hadron } // } // else if(ddMass>mSigZ+mLamb) // Lambda+Sigma0 is possible { // @@ Only two particles PS is used G4double dd2=ddMass*ddMass; // Squared free energy G4double sma=mLamb+mLamb; // Lambda+Lambda sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Lamb+Lamb PS sma=mLamb+mSigZ; // Lambda+Sigma0 sum G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference G4double pr2=pr1; // Prototype of +L+S0 PS if(ddMass>sma && ddMass>smi) pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); if(pr2*G4UniformRand()>pr1) // --> "ENnv+Sigma0+Lambda" case { // h1QPDG=s0QPDG; // QPDG of the first hadron h1M =mSigZ; // Mass of the first hadron } // } // } // if(h1M+h2M+envM "Three parts decay" case { G4LorentzVector h14M(0.,0.,0.,h1M); G4LorentzVector h24M(0.,0.,0.,h2M); G4LorentzVector e4M(0.,0.,0.,envM); if(!G4QHadron(tot4M).DecayIn3(h14M,h24M,e4M)) { G4cout<<"Warning->G4QE::HQE:M="<"< than decaying mass"); #endif CleanUp(); G4QHadron* evH = new G4QHadron(totQC,tot4M);// Create a Hadron for ResidualNucl EvaporateResidual(evH); // Try to evaporate residual (del. equiv.) return theQHadrons; } else // No environment { G4int nHadrs=theQHadrons.size(); // #of available hadrons for(G4int ih=0; ihGet4Momentum(); G4double chM=ch4M.m(); G4LorentzVector tch4M=ch4M+tot4M; if(tch4M.m() > chM + totM) // Can be corrected { G4LorentzVector h14M(0.,0.,0.,chM); G4LorentzVector h24M(0.,0.,0.,totM); if(!G4QHadron(tch4M).DecayIn2(h14M,h24M)) { G4cout<<"-Warning->G4QE::HQE:M="<"<Set4Momentum(h14M); // Change 4M of the current hadron break; // Quit the loop } } } G4QHadron* rH = new G4QHadron(totQC,tot4M);// Create a Hadron for ResidualNucl theQHadrons.push_back(rH); } } else if(spbRN)// => "Join all quasmons to the residual compound and evaporate" case { #ifdef fdebug G4cout<<"***G4QEnv::HadQEnv: Evaporate the total residual tRN="<3) // "Try to correct" case (change condition) { #ifdef debug G4cout<<"***G4QEnv::HadrQE: M="<Get4Momentum(); G4QContent lastQC = theLast->GetQC(); G4int lastS = lastQC.GetStrangeness(); G4int totS = totQC.GetStrangeness(); G4int nFr = theLast->GetNFragments(); G4int gam = theLast->GetPDGCode(); if(gam!=22&&!nFr&&lastS<0&&lastS+totS<0&&nOfOUT>1) // => "Skip K,gam & decayed" { G4QHadron* thePrev = theQHadrons[nOfOUT-2]; theQHadrons.pop_back(); // the last QHadron is excluded from OUTPUT theQHadrons.pop_back(); // the prev QHadron is excluded from OUTPUT theQHadrons.push_back(thePrev); // thePast becomes theLast as an instance delete theLast; // theLast QHadron is deleated as an instance theLast = thePrev; // Update parameters(thePrev becomes theLast) last4M = theLast->Get4Momentum(); lastQC = theLast->GetQC(); } else { theQHadrons.pop_back(); // the last QHadron is excluded from OUTPUT delete theLast; // theLastQHadron is deleated as an instance } totQC+=lastQC; // Update (increase) the total QC tot4M+=last4M; // Update (increase) the total 4-momentum totMass=tot4M.m(); // Calculate new real total mass G4int bn=totQC.GetBaryonNumber(); // The BaryNum after addition totPDG=totQC.GetSPDGCode(); if(totPDG==10&&totQC.GetBaryonNumber()>0) totPDG=totQC.GetZNSPDGCode(); if(bn>1) { totS =totQC.GetStrangeness(); // Total Strangeness of this System if(totS>=0) // => "This is a normal nucleus" case { G4QNucleus newN(totQC,tot4M); totPDG=newN.GetPDG(); totM =newN.GetMZNS(); // Calculate new minimum (GS) mass } else if(totS==-1) // => "Try to decay in K+/aK0 and finish" { G4double m1=mK; G4int PDG1=321; G4QNucleus newNp(totQC-KpQC); G4int PDG2=newNp.GetPDG(); G4double 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(); } if(totMass>m1+m2) // => "can decay" case { G4LorentzVector fq4M(0.,0.,0.,m1); G4LorentzVector qe4M(0.,0.,0.,m2); if(!G4QHadron(tot4M).DecayIn2(fq4M,qe4M)) { G4cout<<"---Warning---G4QE::HadQE:tM="<aK="<"Try to decay in 2(K+/aK0) and finish" { G4double m1=mK; G4int PDG1=321; G4double m2=mK0; G4int PDG2=311; G4QNucleus newNp0(totQC-KpQC-K0QC); G4int PDG3=newNp0.GetPDG(); G4double m3=newNp0.GetMZNS(); // M-K^0-K^+ G4QNucleus newN00(totQC-K0QC-K0QC); G4double m4=newN00.GetMZNS(); // M-2*K^0 G4QNucleus newNpp(totQC-KpQC-KpQC); G4double m5=newNpp.GetMZNS(); // M-2*K^+ if (m4+mK0+mK0"2K0+ResA is theBest" { m1 =mK0; PDG1=311; m3 =m4; PDG3=newN00.GetPDG(); } else if(m5+mK+mK"2Kp+ResA isTheBest" { m2 =mK; PDG1=321; m3 =m5; PDG3=newNpp.GetPDG(); } if(totMass>m1+m2+m3) // => "can decay" case { G4LorentzVector k14M(0.,0.,0.,m1); G4LorentzVector k24M(0.,0.,0.,m2); G4LorentzVector ra4M(0.,0.,0.,m3); if(!G4QHadron(tot4M).DecayIn3(k14M,k24M,ra4M)) { G4cout<<"--Warning--G4QE::HQE:tM="<aK="< "Continue reversion" case } else { if (totPDG==1114||totPDG==2224||totPDG==10) // Decay right now and finish { G4double m1=mNeut; G4int PDG1=2112; G4double m2=mPi; G4int PDG2=-211; if(totPDG==2224) { m1=mProt; PDG1=2212; m2=mPi; PDG2=211; } else if(totPDG==10) // "Chipolino" case { G4QChipolino resChip(totQC); // define the residual as a Chipolino G4QPDGCode h1=resChip.GetQPDG1(); PDG1=h1.GetPDGCode(); // PDG code of the first hadron m1 =h1.GetMass(); // Mass of the first hadron G4QPDGCode h2=resChip.GetQPDG2(); PDG2=h2.GetPDGCode(); // PDG code of the second hadron m2 =h2.GetMass(); // Mass of the second hadron } if(totMass>m1+m2) { G4LorentzVector fq4M(0.,0.,0.,m1); G4LorentzVector qe4M(0.,0.,0.,m2); if(!G4QHadron(tot4M).DecayIn2(fq4M,qe4M)) { G4cout<<"---Warning---G4QE::HaQE:tM="< h1="<KillQuasmon(); } // End of CleanUp //Evaporate Residual Nucleus void G4QEnvironment::EvaporateResidual(G4QHadron* qH, G4bool fCGS) {// ============================================================== static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0); static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0); static const G4double mNeut = G4QPDGCode(2112).GetMass(); static const G4double mProt = G4QPDGCode(2212).GetMass(); static const G4double mAlPr = mAlph+mProt; static const G4double mAlNt = mAlph+mNeut; static const G4double dProt = mProt+mProt; static const G4double dNeut = mNeut+mNeut; static const G4double dAlph = mAlph+mAlph; static const G4double eps=.003; G4int thePDG = qH->GetPDGCode(); // Get PDG code of the Residual Nucleus G4int theBN = qH->GetBaryonNumber(); // A (Baryon number of the nucleus) G4QContent theQC = qH->GetQC(); // Quark Content of the hadron G4int theS=theQC.GetStrangeness(); // S (Strangeness of the nucleus) #ifdef debug G4cout<<"G4QE::EvaporateRes:Called for PDG="<Get4Momentum()<GetQC(); // Quark content of the Hadron-Chipolino G4QChipolino QCh(chQC); // Define a Chipolino instance for the Hadron G4LorentzVector ch4M=qH->Get4Momentum(); // 4Mom of the Hadron-Chipolino G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron G4double h1M =h1QPDG.GetMass(); // Mass of the first hadron G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron G4double h2M =h2QPDG.GetMass(); // Mass of the second hadron G4double chM2 =ch4M.m2(); // Squared Mass of the Chipolino if( sqr(h1M+h2M) < chM2 ) // Decay is possible { G4LorentzVector h14M(0.,0.,0.,h1M); G4LorentzVector h24M(0.,0.,0.,h2M); if(!G4QHadron(ch4M).DecayIn2(h14M,h24M)) { G4cerr<<"***G4QE::EvaporateRes: CM="< h1="< H1="< H2="<Get4Momentum(); // Get 4-momentum of theTotalResidNucleus G4double totMass = q4M.m(); // Get theRealMass of theTotalResidNucleus if(fabs(totMass-totGSM)totGSM) { theEnvironment.EvaporateNucleus(qH,&theQHadrons); #ifdef qdebug qH=0; #endif } else // Correction must be done { #ifdef debug G4cout<<"G4QE::EvaRes: *Correct* "<G4QE::HQ:h#"<GetQC()<<",PDG="<GetPDGCode()<G4QE::HQ:q#"<GetCharge()<<",QCont="<GetQC()<push_back(curQ); // deleted after the "while LOOP" } } G4QHadronVector* reQHadrons = new G4QHadronVector; // deleted after the "while LOOP" G4int nH = theQHadrons.size(); if(nH) { for(G4int ih=0; ihpush_back(curH); // deleted after the "while LOOP" } } G4QNucleus reEnvironment=theEnvironment; G4LorentzVector rem4M=tot4Mom; while (RepFlag && ExCountsize(); //intQHadrons.resize(tmpS+intQHadrons.size()); //copy(theFragments->begin(), theFragments->end(), intQHadrons.end()-tmpS); //tmpS=intQHadrons.size(); //theFragments->resize(tmpS); // Resize theFragments //copy(intQHadrons.begin(), intQHadrons.end(), theFragments->begin()); // Can be like this, but by performance it is closer to FNAL (but better than it) //copy(theFragments->begin(), theFragments->end(), back_inserter(intQHadrons)); //theFragments->resize(intQHadrons.size()); // Resize theFragments //copy(intQHadrons.begin(), intQHadrons.end(), theFragments->begin()); // The following (istead of all above) has worse performance ! theFragments->insert(theFragments->begin(), intQHadrons.begin(), intQHadrons.end() ); intQHadrons.clear(); } return theFragments; } // End of the Fragmentation member function //The Final State Interaction Filter for the resulting output of ::HadronizeQEnvironment() G4QHadronVector* G4QEnvironment::FSInteraction() {// =============================== static const G4QPDGCode gQPDG(22); static const G4QPDGCode pizQPDG(111); static const G4QPDGCode pipQPDG(211); static const G4QPDGCode pimQPDG(-211); static const G4QPDGCode nQPDG(2112); static const G4QPDGCode pQPDG(2212); static const G4QPDGCode lQPDG(3122); static const G4QPDGCode s0QPDG(3212); //static const G4QPDGCode dQPDG(90001001); static const G4QPDGCode tQPDG(90001002); static const G4QPDGCode he3QPDG(90002001); static const G4QPDGCode aQPDG(90002002); static const G4QPDGCode a6QPDG(90002004); static const G4QPDGCode be6QPDG(90004002); //static const G4QPDGCode b7QPDG(90005002); //static const G4QPDGCode he7QPDG(90002005); static const G4QPDGCode c8QPDG(90006002); static const G4QPDGCode a8QPDG(90002006); static const G4QPDGCode c10QPDG(90006004); static const G4QPDGCode o14QPDG(90008006); static const G4QPDGCode o15QPDG(90008007); static const G4QContent K0QC(1,0,0,0,0,1); static const G4QContent KpQC(0,1,0,0,0,1); static const G4double mPi = G4QPDGCode(211).GetMass(); static const G4double mPi0 = G4QPDGCode(111).GetMass(); static const G4double mK = G4QPDGCode(321).GetMass(); static const G4double mK0 = G4QPDGCode(311).GetMass(); static const G4double mNeut= G4QPDGCode(2112).GetMass(); static const G4double mProt= G4QPDGCode(2212).GetMass(); static const G4double mLamb= G4QPDGCode(3122).GetMass(); static const G4double mSigZ= G4QPDGCode(3212).GetMass(); static const G4double mSigM= G4QPDGCode(3112).GetMass(); static const G4double mSigP= G4QPDGCode(3222).GetMass(); static const G4double mXiZ = G4QPDGCode(3322).GetMass(); static const G4double mXiM = G4QPDGCode(3312).GetMass(); static const G4double mOmM = G4QPDGCode(3334).GetMass(); static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0); static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0); static const G4double mHe3 = G4QPDGCode(2112).GetNuclMass(2,1,0); static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); static const G4double mHe6 = G4QPDGCode(2112).GetNuclMass(2,4,0); static const G4double mBe6 = G4QPDGCode(2112).GetNuclMass(4,2,0); //static const G4double mHe7 = G4QPDGCode(2112).GetNuclMass(2,5,0); //static const G4double mB7 = G4QPDGCode(2112).GetNuclMass(5,2,0); static const G4double mHe8 = G4QPDGCode(2112).GetNuclMass(2,6,0); static const G4double mC8 = G4QPDGCode(2112).GetNuclMass(6,2,0); static const G4double mC10 = G4QPDGCode(2112).GetNuclMass(6,4,0); static const G4double mO14 = G4QPDGCode(2112).GetNuclMass(8,6,0); static const G4double mO15 = G4QPDGCode(2112).GetNuclMass(8,7,0); static const G4double mKmP = mK+mProt; static const G4double mKmN = mK+mNeut; static const G4double mK0mP = mK0+mProt; static const G4double mK0mN = mK0+mNeut; static const G4QNucleus vacuum(90000000); static const G4double eps=0.003; ///////////////static const G4double third=1./3.; ///////////////static const G4double nPDG=90000001; G4int envA=theEnvironment.GetBaryonNumber(); ///////////////G4int envC=theEnvironment.GetCharge(); #ifdef rdebug G4int totInC=theEnvironment.GetZ(); G4LorentzVector totIn4M=theEnvironment.Get4Momentum(); G4cout<<"G4QEnvironment(G4QE)::FSInter(FSI): ***called*** envA="<Get4Momentum(); G4cout<<"G4QE::FSI: Quasmon ("<GetBaryonNumber(); #ifdef debug G4cout<<"G4QE::FSI:after HQE,nH="<1) // TheLastHadron is nucleus:try to decay/evap/cor it { G4QHadron* theLast = theQHadrons[nHadr-1]; G4QHadron* curHadr = new G4QHadron(theLast); G4LorentzVector lh4M=curHadr->Get4Momentum(); // Actual mass of the last fragment G4double lhM=lh4M.m(); // Actual mass of the last fragment G4int lhPDG=curHadr->GetPDGCode(); // PDG code of the last fragment G4double lhGSM=G4QPDGCode(lhPDG).GetMass(); // GroundStateMass of the last fragment #ifdef debug G4cout<<"G4QE::FSI:lastHadr 4M/M="<lhGSM+eps) // ==> Try to evaporate the residual nucleus { theQHadrons.pop_back(); // the last QHadron-Nucleus is excluded from OUTPUT delete theLast;//*!!When kill,DON'T forget to delete theLastQHadron as an instance!* EvaporateResidual(curHadr); // Try to evaporate Hadr-Nucl (@@DecDib)(delete eq.) nHadr=theQHadrons.size(); #ifdef debug G4cout<<"G4QE::FSI:After nH="<GetPDGCode()< Try to make the HadronicSteck FSI correction { theQHadrons.pop_back(); //exclude LastHadronPointer from OUTPUT delete theLast; //*!! When killing, DON'T forget to delete the last QHadron !!* G4Quasmon* quasH = new G4Quasmon(curHadr->GetQC(),lh4M); // Fake Quasmon ctreation if(!CheckGroundState(quasH,true))// Try to correct with other hadrons { #ifdef debug // M.K. Fake complain in the low energy nHe/pHe reactions, while everything is OK G4cout<<"---Warning---G4QEnv::FSI:Correction error LeaveAsItIs h4m="< Leave the nucleus as it is (close to the GSM) } #ifdef debug G4LorentzVector ccs4M(0.,0.,0.,0.); // CurrentControlSum of outgoing Hadrons #endif // *** Initial Charge Control Sum Calculation *** G4int chContSum=0; // ChargeControlSum to keepTrack FSI transformations G4int bnContSum=0; // BaryoNControlSum to keepTrack FSI transformations if(nHadr)for(unsigned ich=0; ichGetNFragments())) { chContSum+=theQHadrons[ich]->GetCharge(); bnContSum+=theQHadrons[ich]->GetBaryonNumber(); } #ifdef chdebug if(chContSum!=totCharge || bnContSum!=totBaryoN) { G4cout<<"***G4QE::Fr:(1)tC="<"<SetQPDG(fQPDG); theQHadrons[ipo]->Set4Momentum(f4M); G4QHadron* sH = new G4QHadron(sPDG,s4M); theQHadrons.push_back(sH); // (delete equivalent) } } else { G4bool fOK=true; G4LorentzVector f4M(0.,0.,0.,fM); G4LorentzVector s4M(0.,0.,0.,sM); G4LorentzVector t4M(0.,0.,0.,tM); G4double sum=fM+sM+tM; if(fabs(hM-sum)<=eps) { f4M=h4Mom*(fM/sum); s4M=h4Mom*(sM/sum); t4M=h4Mom*(tM/sum); } else if(hM"<SetQPDG(fQPDG); theQHadrons[ipo]->Set4Momentum(f4M); G4QHadron* sH = new G4QHadron(sPDG,s4M); theQHadrons.push_back(sH); // (delete equivalent) G4QHadron* tH = new G4QHadron(tPDG,t4M); theQHadrons.push_back(tH); // (delete equivalent) } } nHadr=theQHadrons.size(); } else if(hPDG==89999003||hPDG==90002999||hPDG==90000003||hPDG==90003000|| hPDG==90999002||hPDG==91001999) // "3-particles decays of dibaryons and 3N" { #ifdef debug G4cout<<"G4QE::FSI:***nD-/pD++/nnn/ppp***i="<mSigZ+mNeut+mPi) { G4double ddMass=hM-mPi; // Free CM energy G4double dd2=ddMass*ddMass; // Squared free energy G4double sma=mLamb+mNeut; // Neutron+Lambda sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Neut+Lamb PS sma=mNeut+mSigZ; // Neutron+Sigma0 sum G4double smi=mSigZ-mNeut; // Sigma0-Neutron difference G4double pr2=pr1; // Prototype of +N+S0 PS if(ddMass>sma && ddMass>smi) pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); if(pr2*G4UniformRand()>pr1) // --> "ENnv+Sigma0+Lambda" case { barPDG = 3212; // Substitute Sigma0 for the second n barM = mSigZ; } else { barPDG = 3122; // Substitute Lambda for the second n barM = mLamb; } } else if(hM>mLamb+mNeut+mPi) { barPDG = 3122; // Substitute Lambda for the second n barM = mLamb; } else if(hM>mSigM+mNeut)// @@ Decay in 2 { barPDG = 3112; // Substitute Sigma- for the second n barM = mSigM; tPDG = 22; tM = 0; } } else if(hPDG==91001999) // p,Lambda,pi+/p,Sigma0,pi+/p,Sigma+,gamma(@@) { nuQPDG = pQPDG; // Substitute p for the first n nucM = mProt; tPDG = 211; // Substitute pi+ for the first pi- if(hM>mSigZ+mProt+mPi) { G4double ddMass=hM-mPi; // Free CM energy G4double dd2=ddMass*ddMass; // Squared free energy G4double sma=mLamb+mProt; // Lambda+Proton sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(ddMass>sma) pr1=sqrt((dd2-sma*sma)*dd2); // Lamb+Prot PS sma=mProt+mSigZ; // Proton+Sigma0 sum G4double smi=mSigZ-mProt; // Sigma0-Proton difference G4double pr2=pr1; // Prototype of +P+S0 PS if(ddMass>sma && ddMass>smi) pr2+=sqrt((dd2-sma*sma)*(dd2-smi*smi)); if(pr2*G4UniformRand()>pr1) // --> "ENnv+Sigma0+Lambda" case { barPDG = 3212; // Substitute Sigma0 for the second n barM = mSigZ; } else { barPDG = 3122; // Substitute Lambda for the second n barM = mLamb; } } if(hM>mLamb+mProt+mPi) { barPDG = 3122; // Substitute Lambda for the second n barM = mLamb; } else if(hM>mSigP+mProt) // @@ Decay in 2 { barPDG = 3222; // Substitute Sigma- for the second n barM = mSigP; tPDG = 22; tM = 0; } } else if(hPDG==90000003) // 3n { tPDG = 2112; // Substitute n for pi- tM = mNeut; } G4bool fOK=true; G4LorentzVector nu4M(0.,0.,0.,nucM); G4LorentzVector ba4M(0.,0.,0.,barM); G4LorentzVector pi4M(0.,0.,0.,tM); G4double sum=nucM+barM+tM; if(fabs(hM-sum)<=eps) { nu4M=h4Mom*(nucM/sum); ba4M=h4Mom*(barM/sum); pi4M=h4Mom*(tM/sum); } if(hM N="<SetQPDG(nuQPDG); theQHadrons[ipo]->Set4Momentum(nu4M); G4QHadron* baH = new G4QHadron(barPDG,ba4M); theQHadrons.push_back(baH); // (delete equivalent) G4QHadron* piH = new G4QHadron(tPDG,pi4M); theQHadrons.push_back(piH); // (delete equivalent) nHadr=theQHadrons.size(); } } else if(hBN>1 && !sBN && (cBN<0 || cBN>hBN)) // "nN+mD- or nP+mD++ decay" { #ifdef debug G4cout<<"G4QE::FSI:nNmD-/nPmD++ #"<hBN) // reinitialization for the "(n+m)*P,m*pi+" case { nuQPDG = pQPDG; // Substitute p for the first n nucM = mProt; barPDG = 2212; // Substitute p for the second n barM = mProt; tPDG = 211; // Substitute pi+ for the first pi- nPi = cBN-hBN; // a#0f Pi+'s } if(nPi>1) tM*=nPi; if(nN >1) barM*=nN; G4bool fOK=true; G4LorentzVector nu4M(0.,0.,0.,nucM); G4LorentzVector ba4M(0.,0.,0.,barM); G4LorentzVector pi4M(0.,0.,0.,tM); G4double sum=nucM+barM+tM; if(fabs(hM-sum)<=eps) { nu4M=h4Mom*(nucM/sum); ba4M=h4Mom*(barM/sum); pi4M=h4Mom*(tM/sum); } if(hMN="<SetQPDG(nuQPDG); theQHadrons[ipo]->Set4Momentum(nu4M); if(nN>1) ba4M=ba4M/nN; for(G4int ib=0; ib1) pi4M=pi4M/nPi; for(G4int im=0; imGetNFragments(); //Recover after swapping G4QContent hQC = theQHadrons[ipo]->GetQC(); // ... hPDG = theQHadrons[ipo]->GetPDGCode(); // ... h4Mom = theQHadrons[ipo]->Get4Momentum(); // ... ccs4M+=h4Mom; // Calculate CurSum of Hadrs G4cout<<"G4QE::FSI:#"<>>CurrentControlSumOf4MomOfHadrons="<GetNFragments())) { ccContSum+=theQHadrons[ic1]->GetCharge(); cbContSum+=theQHadrons[ic1]->GetBaryonNumber(); } if(ccContSum-chContSum || cbContSum-bnContSum) { G4cout<<"*G4QE::FSI:(1) dC="<0) p2cut/=envA*envA; //G4double p2cut2=0.; //cut for the alpha creation // G4int bfCountM=3; if(envA>10) bfCountM*=(envA-1)/3; G4bool bfAct = true; G4int bfCount= 0; G4LorentzVector tmp4Mom=tot4Mom; G4LorentzVector postp4M(0.,0.,0.,0.); G4int nPost=intQHadrons.size(); if(nPost) for(G4int psp=0; pspGetNFragments())) postp4M+=intQHadrons[psp]->Get4Momentum(); while(bfAct&&bfCountGetPDGCode(); G4QPDGCode hQPDG(hPDG); G4double hGSM = hQPDG.GetMass(); // Ground State Mass of the first fragment #ifdef debug G4cout<<"G4QE::FSI:LOOP START,h#"<Get4Momentum()<Get4Momentum(); G4double hM=h4Mom.m(); G4QPDGCode fQPDG=nQPDG; G4double fM =mNeut; G4int sPDG =2112; G4double sM =mNeut; G4int tPDG =-211; G4double tM =mPi; if(hPDG==90002999) { fQPDG=pQPDG; fM =mProt; sPDG =2212; sM =mProt; tPDG =211; } G4bool fOK=true; G4LorentzVector f4M(0.,0.,0.,fM); G4LorentzVector s4M(0.,0.,0.,sM); G4LorentzVector t4M(0.,0.,0.,tM); G4double sum=fM+sM+tM; if(fabs(hM-sum)<=eps) { f4M=h4Mom*(fM/sum); s4M=h4Mom*(sM/sum); t4M=h4Mom*(tM/sum); } else if(hM"<SetQPDG(fQPDG); curHadr->Set4Momentum(f4M); G4QHadron* sH = new G4QHadron(sPDG,s4M); theQHadrons.push_back(sH); // (delete equivalent) G4QHadron* tH = new G4QHadron(tPDG,t4M); theQHadrons.push_back(tH); // (delete equivalent) } hPDG = curHadr->GetPDGCode(); // Change PDG Code of theFirstFragment hQPDG= G4QPDGCode(hPDG); hGSM = hQPDG.GetMass(); // Change GroundStateMass of theFirstFragm } nHadr=theQHadrons.size(); if(hPDG==89001001||hPDG==89002000||hPDG==89000002) { G4cout<<"---WARNING---G4QE::FSI:***(K+N)*** (2),PDG="<Get4Momentum(); G4double hM=h4Mom.m(); G4QPDGCode fQPDG=nQPDG; G4double fM =mNeut; G4int sPDG =311; G4double sM =mK0; if(hPDG==89000002) { fQPDG=pQPDG; fM =mProt; sPDG =321; sM =mK; } if(hPDG==89001001) { if(hM.5) { sPDG =321; sM =mK; } else { fQPDG=pQPDG; fM =mProt; } } G4bool fOK=true; G4LorentzVector f4M(0.,0.,0.,fM); G4LorentzVector s4M(0.,0.,0.,sM); G4double sum=fM+sM; if(fabs(hM-sum)<=eps) { f4M=h4Mom*(fM/sum); s4M=h4Mom*(sM/sum); } else if(hM"<SetQPDG(fQPDG); curHadr->Set4Momentum(f4M); G4QHadron* sH = new G4QHadron(sPDG,s4M); theQHadrons.push_back(sH); // (delete equivalent) } hPDG = curHadr->GetPDGCode(); // Change PDG Code of theFirstFragment hQPDG= G4QPDGCode(hPDG); hGSM = hQPDG.GetMass(); // Change GroundStateMass of theFirstFragm } nHadr=theQHadrons.size(); G4int hS = curHadr->GetStrangeness(); G4int hF = curHadr->GetNFragments(); G4LorentzVector h4m= curHadr->Get4Momentum(); G4double hM = h4m.m(); // Real Mass of the first fragment G4int hB = curHadr->GetBaryonNumber(); //////////////////////G4int hC = curHadr->GetCharge(); #ifdef debug if(!hF&&(hPDG>80000000&&hPDG<90000000||hPDG==90000000|| hPDG>90000000&&(hPDG%1000000>200000||hPDG%1000>300))) G4cout<<"**G4QEnv::FSInteraction: PDG("<0&&!hS&&(nHadr>3||hB<2)) // ThermoBackFus (VIMP for gamA TotCS) //if(hadron&&!hF&&hB>0&&!hS&&(nHadr>2||hB<4)) // ThermoBackFus (VIMP for gamA TotCS) if(hadron&&!hF&&hB>0&&!hS) // ThermoBackFusion cond. (VIMP for gamA TotCS) //if(hadron&&!hF&&hB>0&&hB<4&&!hS) // ThermoBackFusion cond. (VIMP for gamA TotCS) //if(hadron&&!hF&&hB>0&&!hS&&nHadr>2)//ThermoBackFusion MAX condition (VIMP for gamA) //if(2>3) // Close the ThermoBackFusion (VIMP for gamA TotCS) { #ifdef debug //if(nHadr==3) G4cout<<"G4QE::FSI: h="<tM="<tM+.001&&pCM20||hB==1&&bB>0)&&bM+hM>tM+.00001 // && (pCM23||pCM2tM+.001&&(pCM23 || // pCM290000000)) //if(!bF&&bB<4&&bM+hM>tM+.001&&pCM20&&bM+hM>tM+.001&&pCM2tM+.001&&(pCM23||bPDG>90000000)&&bM+hM>tM+.001&&pCM2tM+.001&&pCM2tM+.001&&sM-bM-hMtM+.001&&sM-bM-hMGetPDGCode(); if(sPDG==89999003||sPDG==90002999||sPDG==89999002||sPDG==90001999) G4cout<<"G4QE::FSI:**nD-/pD++**,h="<["<mSigZ+mSigZ) // Sigma0+Sigma0 is possible { // @@ Only two particles PS is used G4double sma=mLamb+mLamb; // Lambda+Lambda sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(sM>sma) pr1=sqrt((sM2-sma*sma)*sM2); // Lamb+Lamb PS sma=mLamb+mSigZ; // Lambda+Sigma0 sum G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference G4double pr2=pr1; // Prototype of +L +S0 PS if(sM>sma && sM>smi) pr2+=sqrt((sM2-sma*sma)*(sM2-smi*smi)); sma=mSigZ+mSigZ; // Sigma0+Sigma0 sum G4double pr3=pr2; // Prototype of +Sigma0+Sigma0 PS if(sM>sma) pr3+=sqrt((sM2-sma*sma)*sM2); G4double hhRND=pr3*G4UniformRand(); // Randomize PS if(hhRND>pr2) // --> "ENnv+Sigma0+Sigma0" case { // fQPDG=s0QPDG; f4Mom=G4LorentzVector(0.,0.,0.,mSigZ); rQPDG=s0QPDG; g4Mom=f4Mom; } // else if(hhRND>pr1) // --> "ENnv+Sigma0+Lambda" case { // fQPDG=s0QPDG; f4Mom=G4LorentzVector(0.,0.,0.,mSigZ); } // } else if(sM>mSigZ+mLamb) // Lambda+Sigma0 is possible { // @@ Only two particles PS is used G4double sma=mLamb+mLamb; // Lambda+Lambda sum G4double pr1=0.; // Prototype to avoid sqrt(-) if(sM>sma) pr1=sqrt((sM2-sma*sma)*sM2); // Lamb+Lamb PS sma=mLamb+mSigZ; // Lambda+Sigma0 sum G4double smi=mSigZ-mLamb; // Sigma0-Lambda difference G4double pr2=pr1; // Prototype of +L +S0 PS if(sM>sma && sM>smi) pr2+=sqrt((sM2-sma*sma)*(sM2-smi*smi)); if(pr2*G4UniformRand()>pr1) // --> "ENnv+Sigma0+Lambda" case { // fQPDG=s0QPDG; f4Mom=G4LorentzVector(0.,0.,0.,mSigZ); } // } // } else if(sPDG==89999003) // A=2 { hQPDG=nQPDG; rQPDG=nQPDG; fQPDG=pimQPDG; t4Mom=G4LorentzVector(0.,0.,0.,mNeut); g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mPi); three=true; } else if(sPDG==90002999) { hQPDG=pQPDG; rQPDG=pQPDG; fQPDG=pipQPDG; t4Mom=G4LorentzVector(0.,0.,0.,mProt); g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mPi); three=true; } else if(sPDG==90000003) // A=3 { hQPDG=nQPDG; rQPDG=nQPDG; fQPDG=nQPDG; t4Mom=G4LorentzVector(0.,0.,0.,mNeut); g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mNeut); three=true; } else if(sPDG==90003000) { hQPDG=pQPDG; rQPDG=pQPDG; fQPDG=pQPDG; t4Mom=G4LorentzVector(0.,0.,0.,mProt); g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mProt); three=true; } else if(sPDG==90001003) // A=4 { rQPDG=nQPDG; fQPDG=tQPDG; g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mTrit); } else if(sPDG==90003001) { rQPDG=pQPDG; fQPDG=he3QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mHe3); } else if(sPDG==90002003) // A=5 { rQPDG=nQPDG; fQPDG=aQPDG; g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mAlph); } else if(sPDG==90003002) { rQPDG=pQPDG; fQPDG=aQPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mAlph); } else if(sPDG==90004002) // A=6 { hQPDG=pQPDG; rQPDG=pQPDG; fQPDG=aQPDG; t4Mom=G4LorentzVector(0.,0.,0.,mProt); g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mAlph); three=true; } else if(sPDG==90002005) // A=7 { rQPDG=nQPDG; fQPDG=a6QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mHe6); } else if(sPDG==90005002) { rQPDG=pQPDG; fQPDG=be6QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mBe6); } else if(sPDG==90004004) // A=8 { fQPDG=aQPDG; rQPDG=aQPDG; f4Mom=G4LorentzVector(0.,0.,0.,mAlph); g4Mom=f4Mom; } //else if(sPDG==90006002) //{ // hQPDG=pQPDG; // rQPDG=pQPDG; // fQPDG=be6QPDG; // t4Mom=G4LorentzVector(0.,0.,0.,mProt); // g4Mom=G4LorentzVector(0.,0.,0.,mProt); // f4Mom=G4LorentzVector(0.,0.,0.,mBe6); // three=true; //} //else if(sPDG==90002006) //{ // hQPDG=nQPDG; // rQPDG=nQPDG; // fQPDG=a6QPDG; // t4Mom=G4LorentzVector(0.,0.,0.,mNeut); // g4Mom=G4LorentzVector(0.,0.,0.,mNeut); // f4Mom=G4LorentzVector(0.,0.,0.,mHe6); // three=true; //} else if(sPDG==90002007) // A=9 { rQPDG=nQPDG; fQPDG=a8QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mNeut); f4Mom=G4LorentzVector(0.,0.,0.,mHe8); } else if(sPDG==90005004) // A=9 { rQPDG=pQPDG; fQPDG=aQPDG; hQPDG=aQPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mAlph); t4Mom=G4LorentzVector(0.,0.,0.,mAlph); three=true; } else if(sPDG==90007002) // A=9 { rQPDG=pQPDG; fQPDG=c8QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mC8); } else if(sPDG==90008004) // A=12 { hQPDG=pQPDG; rQPDG=pQPDG; fQPDG=c10QPDG; t4Mom=G4LorentzVector(0.,0.,0.,mProt); g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mC10); three=true; } else if(sPDG==90009006) // A=15 { rQPDG=pQPDG; fQPDG=o14QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mO14); } else if(sPDG==90009007) // A=16 { rQPDG=pQPDG; fQPDG=o15QPDG; g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mO15); } else if(sPDG==90010006) // A=16 { hQPDG=pQPDG; rQPDG=pQPDG; fQPDG=o14QPDG; t4Mom=G4LorentzVector(0.,0.,0.,mProt); g4Mom=G4LorentzVector(0.,0.,0.,mProt); f4Mom=G4LorentzVector(0.,0.,0.,mO14); three=true; } #ifdef debug G4cout<<"G4QE::FSI: "<sM="<sM="<Get4Momentum(); // 4-mom of the current hadron tot4Mom-=cH4Mom; // Reduce theTotal 4mom by theCurrent 4mom totCharge-=curHadr->GetCharge(); // @!@ totBaryoN-=curHadr->GetBaryonNumber(); // @!@ #ifdef pdebug G4cout<<"G4QE::FSI:Cur4M="< 0. && !totCharge && !totBaryoN && sdm/re2 < .0001)) // @!@ { #ifdef pdebug G4cout<<"...G4QE::FSI:E/M conservation is corrected by a photon"<900.) // MomentumIsConserved - recoverMissing { if(fabs(re-mNeut)<.01) { #ifdef pdebug G4cout<<"...G4QE::FSI:E/M conservation is corrected by neutron"<Get4Momentum(); // 4mom of thePrevHadron G4double cHM = curHadr->GetMass(); // Mass of the current hadron G4double pHM = prevHadr->GetMass(); // Mass of the previous hadron tot4Mom+=cH4Mom+pH4Mom; G4double totRM=tot4Mom.m(); if(cHM+pHM<=totRM) // *** Make the final correction *** { if(!G4QHadron(tot4Mom).DecayIn2(pH4Mom,cH4Mom)) { G4cout<<"***G4QE::FSI:**Correction**tot4M="<sM=" <"<GetMass(); // Mass of the current hadron G4int ch=0; for(ch=nHadr-2; ch>-1; --ch) { G4QHadron* prevHadr = theQHadrons[ch]; // GetPointer to Hadr prev to theLast G4LorentzVector pH4Mom = prevHadr->Get4Momentum();// 4-mom of thePreviousHadron G4double pHM = prevHadr->GetMass(); // Mass of the current hadron tot4Mom+=cH4Mom+pH4Mom; G4double totRM=tot4Mom.m(); if(cHM+pHM<=totRM) // *** Make the final correction *** { if(!G4QHadron(tot4Mom).DecayIn2(pH4Mom,cH4Mom)) { G4cout<<"***G4QEnv::FSI:**Correction**,tot4M="< sM=" <| 2 | 2 | 2 | 2 | 2 | 2 - old | 1. If gamma: add to sum4Mom // |>0->sum| 3 | 3 | 3 | 3 | 3 - old | 2. Compare BN with the Last // | 5 |>5 | 4 | 4 | 4 | 4 - old | 3. Swap if larger, del theLast // | 0 | 0 |>0->sum| 5<-sum| 5->evap| 2 - new | 4. If the Last: add the sum // | 4 | 4 | 5 | ex | |(0 - gamma?)| 5. Decay/Eporate the Last // | 3 | ex | | 3 - new #ifdef chdebug // *** (3) Charge Control Sum Calculation for the Charge Conservation Check *** ccContSum=0; // Intermediate ChargeControlSum cbContSum=0; // Intermediate BaryonNumberControlSum if(nHadr)for(unsigned ic3=0; ic3GetNFragments())) { ccContSum+=theQHadrons[ic3]->GetCharge(); cbContSum+=theQHadrons[ic3]->GetBaryonNumber(); } if(ccContSum-chContSum || cbContSum-bnContSum) { G4cout<<"*G4QE::FSI:(3) dC="<2)for(unsigned f=0; fGetBaryonNumber(); // Baryon number of the fragment #ifdef debug G4int fPDG=theQHadrons[f]->GetPDGCode(); // PDG code of the fragment G4LorentzVector fLV=theQHadrons[f]->Get4Momentum(); // 4Mom of the fragment G4cout<<"G4QE::FSI:"<1) frag=true; // The fragment with A>1 is found } #ifdef debug G4cout<<"G4QE::FSI:===Before Gamma Compression===, nH="<2 && frag) for(G4int h=nHadr-1; h>=0; h--)//Collect gammas & kill DecayedHadrs { G4QHadron* curHadr = theQHadrons[h]; // Get a pointer to the current Hadron G4int hF = curHadr->GetNFragments(); G4int hPDG = curHadr->GetPDGCode(); if(hPDG==89999003||hPDG==90002999) G4cout<<"---Warning---G4QEnv::FSI:nD-/pD++(1)="<1)for(unsigned hdr=0; hdrGetBaryonNumber()>0) // "Absorb photons & evaporate/decay" case { G4QHadron* theNew = new G4QHadron(theLast); // Make New Hadron of the Last Hadron #ifdef ffdebug G4cout<<"G4QE::FSI:BeforeLastSub,n="<GetPDGCode()<GetPDGCode(); G4LorentzVector new4M=theNew->Get4Momentum(); // 4-mom of the fragment #ifdef debug G4cout<<"G4QE::FSI:gSum4M="<3) //CloseTheFirstPriorityResN+gamSumEvaporation { theNew->Set4Momentum(exRes4M); // Icrease 4Mom of theLast by "sum" to Evapor #ifdef ffdebug G4cout<<"G4QE::FSI:BeforeE(1),n="<GetPDGCode()<GetPDGCode()==90002002&&exRes4M.m()>mHe3+mNeut&&G4UniformRand()>.5) { theNew->Set4Momentum(exRes4M); // Icrease 4Mom of theLast by "sum" to Evapor G4LorentzVector n4M(0.,0.,0.,mNeut); G4LorentzVector h4M(0.,0.,0.,mHe3); if(!theNew->DecayIn2(n4M,h4M)) { G4cerr<<"***G4QE::FSI:GamSup, tM="<GetNFragments())) { ccContSum+=theQHadrons[ic6]->GetCharge(); cbContSum+=theQHadrons[ic6]->GetBaryonNumber(); } if(ccContSum-chContSum || cbContSum-bnContSum) { G4cout<<"*G4QE::FSI:(6) dC="<GetBaryonNumber(); G4int tHadr=lHadr; // Total Baryon number if(nHadr>1)for(unsigned ipo=0; ipoGetPDGCode(); if(hPDG==22) gamcnt++; else { G4int hBN = theQHadrons[ipo]->GetBaryonNumber(); tHadr+=hBN; #ifdef debug G4cout<<"G4QE::FSI:h#"<lHadr) { lHadr=hBN; maxB=ipo; } // the current biggest nuclear fragment } } #ifdef debug G4cout<<"G4QE::FSI:max#"<GetNFragments())) { ccContSum+=theQHadrons[ic7]->GetCharge(); cbContSum+=theQHadrons[ic7]->GetBaryonNumber(); } if(ccContSum-chContSum || cbContSum-bnContSum) { G4cout<<"*G4QE::FSI:(7) dC="<1) // Only if there are gammas one should act { if(maxB+1GetQPDG(); // The QPDG of the last if(lQP.GetPDGCode()!=10) theCurr->SetQPDG(lQP); //CurHadr instead of PrevHadr else theCurr->SetQC(theLast->GetQC()); // CurHadrPDG instead of LastHadrPDG theCurr->Set4Momentum(theLast->Get4Momentum()); // ... continue substitution theQHadrons.pop_back(); // Rnt to theLastHadron is excluded from HV delete theLast;//*!!When kill,DON'T forget to delete the last QHadron as an inst. !!* theQHadrons.push_back(curHadr); // The current Hadron, which is the Biggest } nHadr=theQHadrons.size(); // Must be the same // Now it is necessary to absorb the photon (photons) and try to radiate alpha or evap. G4LorentzVector gamSum(0.,0.,0.,0.); if(nHadr>1)for(unsigned gp=0; gpGetPDGCode(); #ifdef debug G4cout<<"G4QE::FSI:gp#"<Get4Momentum(); // Accumulate the 4Momenta of the photon #ifdef debug G4cout<<"G4QE::FSI:Photon gp#"<gp) { gamSum=gamSum+theLast->Get4Momentum();// Accumulate 4-momentum of theLastPhoton #ifdef debug G4cout<<"G4QE::FSI:TheLastPhotonIsFound #"<GetPDGCode(); if(hPDG==22 && fabs(curHadr->Get4Momentum().e())<.00001) // E=0, gamma in the OUTPUT { unsigned lin=theQHadrons.size()-1; G4QHadron* theLast = theQHadrons[lin];// Pointer to theLastHadron in theQHadrVector if(lin>hd) { G4QPDGCode lQP=theLast->GetQPDG(); // The QPDG of the last if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP); //CurHadr instead of PrevHadr else curHadr->SetQC(theLast->GetQC()); // CurHadrPDG instead of LastHadrPDG curHadr->Set4Momentum(theLast->Get4Momentum()); // ... continue substitution (4Mom) } //else break; // It's not necessary to delete: not copy to theFragments is enough delete theLast;//*!!When kill, DON'T forget to delete theLastQHadron asAnInstance!! theQHadrons.pop_back(); //pointer to theLastHadron is excluded from HV hPDG=curHadr->GetPDGCode(); // Redefine thePDGCode of theCurrentHadron if(lin==hd) break; } #ifdef fdebug G4cout<<"G4QE::FSI: LOOP starts nH="<GetQPDG(); // The QPDG of the last if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP); //CurHadr instead of PrevHadr else curHadr->SetQC(theLast->GetQC()); // CurHadrPDG instead of LastHadrPDG curHadr->Set4Momentum(theLast->Get4Momentum()); // ... 4Momentum substitution } theQHadrons.pop_back(); // exclude theLastHadron pointer from the OUTPUT delete theLast; //*!! When killing, DON'T forget to delete the last QHadron !!* G4Quasmon* quasH = new G4Quasmon(qH->GetQC(),qH->Get4Momentum());//Fake Quasmon if(!CheckGroundState(quasH,true)) // Try to correct by other hadrons { qH->SetNFragments(-1); //@@ To avoid looping G4cout<<"---Warning---G4QE::FSI:PiNSig Failed, LeaveAsItIs 4m="<GetPDGCode(); G4LorentzVector hpLV=cpHadr->Get4Momentum(); G4cout<<"G4QE::FSI:h#"<"<N="<"<N="<"<N="<1) n4M/=hBN; // Calculate the 4mom per one nucleon curHadr->Set4Momentum(n4M); curHadr->SetQPDG(G4QPDGCode(PDGnu)); // At least one nucleun always exists if(hBN>1) for(G4int ih=1; ihpush_back(Hi); // Fill nucleons (delete equivalent) } if(nPi>1) k4M/=nPi; for(G4int ip=0; ip2 neutrons conversion (***Check it***)"<Get4Momentum())/hBN; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(hPDG)); for(G4int bd=1; bdpush_back(secHadr);// (del.equiv. - user is responsible for that) } } else if(hST<0 && hBN>0) // AntistrangeNucleus (@@ see above, already done) { G4LorentzVector r4M=curHadr->Get4Momentum(); // Real 4-mom of theAntiStrangeNucleus G4double reM=r4M.m(); // Real mass of the antiStrange nucleus G4QContent nQC=curHadr->GetQC(); // QCont of the Antistrange nucleus G4QNucleus newN0(nQC-K0QC); G4int RPDG=newN0.GetPDG(); G4double mR=newN0.GetMZNS(); // Residual mass for K0 G4double mKaon=mK0; // Prototype is mass of K0 G4int kPDG =311; // Prototype is PDG of K0 G4QNucleus newNp(nQC-KpQC); G4double mp=newNp.GetMZNS(); // Residual mass for K+ if(mp+mKreM) // for GEANT4 (Can not decay in kaon and residual) { if(kPDG==321) // *** Very seldom *** "K+->pi+ conversion" { kPDG=211; mKaon=mPi; } else // *** Very seldom *** "K0(S=-1)->pi- conversion" { kPDG=111; mKaon=mPi0; } sum=mR+mKaon; } G4LorentzVector n4M(0.,0.,0.,mR); G4LorentzVector k4M(0.,0.,0.,mKaon); if(fabs(reM-sum)"<A="<Set4Momentum(k4M); // Kaon curHadr->SetQPDG(G4QPDGCode(kPDG)); G4QHadron* theRes = new G4QHadron(RPDG,n4M);// Make a New Hadr for the residual EvaporateResidual(theRes); // Try to evaporate theResNuc (del.eq.) } } // =Lambda(Sigma0)= = Sigma+ = = Sigma- = else if(hPDG>90500000 && hPDG!=91000000 && hPDG!=91000999 && hPDG!=90999001 && hPDG!=91999000 && hPDG!=91999999 && hPDG!=92998999 ) // == Xi- == == Xi0 == === Omega- === { #ifdef pdebug G4cout<<"***G4QEnvironment::FSI:*G4* Hypernucleus PDG="<GetStrangeness(); // A number of Lambdas in the Hypernucleus G4int nB=curHadr->GetBaryonNumber(); // Total Baryon Number of the Hypernucleus G4int nC=curHadr->GetCharge(); // Charge of the Hypernucleus G4int nSM=0; // A#0f unavoidable Sigma- G4int nSP=0; // A#0f unavoidable Sigma+ G4int nXZ=0; // A#0f unavoidable Xi0 G4int nXM=0; // A#0f unavoidable Xi- G4int nOM=0; // A#0f unavoidable Omega- // @@ Now it does not include more complicated clusters as Omega-,Xi- (Improve) if(nC < 0) // Negative hypernucleus { if(-nC <= nL) // Negative Charge is smaller than Strange { if(nL <= nB) // ==> Big Hyper-nucleus { nSM=-nC; // Can be compensated by Sigma- nL+=nC; // Reduce the residual strangeness } else // ==> Lack of BaryonNumber (better use Xi-) { G4int nX=nL-nB; if(-nC <= nX && nL >= nX+nX) // Part or all Xi are Xi- (e.g. Xi-,Lambda) { nXM=-nC; if(nXM != nX) nXZ=nX-nXM; // The rest are Xi0 nL-=nX+nX; // Xi reduce srangeness twice faster } else if(nL >= nX-nC) // Sigma- should be used in addition to Xi- { // e.g. Xi-, Sigma- nXM=nX; nSM=-nC-nX; nL-=nX-nC; } // @@ Don't close this warning as it costs days to find out this growing point! else G4cout<<"-Warning-G4QEn::FSI:*Improve*,-nC<=nL,nL>nB, PDG="<nB-nL) // Extra positive hypernucleus (nC>=0 here) { // Can't compensate theCharge by onlyProtons if(nC<=nB) { if(nL <= nB) { G4int dH=nB-nC; nSP=nL-dH; // Can be compensated by Sigma+ nL=dH; // Reduce the residual strangeness } else if(nL nB && ((nL-nB)/2)*2 == nL-nB && nC+nL <= nB+nB)// Omega- can be used { nOM=(nL-nB)/2; // This is how many Omega- can be made nL=nB+nB-nL-nC; // a#of Lambdas should be 0 ! nSP=nOM-nC; } // @@ Do not close this warning as it costs days to find out this growing point ! else G4cout<<"-Warning-G4QEn::FSI:*Improve*,nC>nB-nL,nC<=nB, PDG="<Get4Momentum(); // Real 4-momentum of the hypernucleus G4double reM=r4M.m(); // Real mass of the hypernucleus // Subtract Lamb/Sig/Xi/Omega from the Nucleus and decay G4int SS=nL+nSP+nSM+nXZ+nXM; if(SS add decay in 3) if(nSP) nL+=nSP; else nL+=nSM; } if(nSP) { nL=nSP; sPDG=3222; MLa=mSigP; } else { nL=nSM; sPDG=3112; MLa=mSigM; } } #ifdef pdebug G4cout<<"G4QEnvironment::FSI:*G4* mS="<A="<Set4Momentum(n4M); curHadr->SetQPDG(G4QPDGCode(rlPDG)); // converted hypernucleus if(rlPDG==90000002) // Additional action with curHadr change { G4LorentzVector newLV=n4M/2.; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90000001)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) } else if(rlPDG==90002000) // Additional action with curHadr change { G4LorentzVector newLV=n4M/2.; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90001000)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) } // @@(?) Add multybaryon decays if necessary (Now it anyhow is made later) #ifdef pdebug G4cout<<"*G4QE::FSI: resNucPDG="<GetPDGCode()<1) h4M=h4M/nL; // split the lambda's 4-mom if necessary for(G4int il=0; ilrnM+mPi0-eps&&!nSP&&!nSM) // Lambda->N only if Sigmas are absent { G4int nPi=static_cast((reM-rnM)/mPi0); if (nPi>nL) nPi=nL; G4double npiM=nPi*mPi0; G4LorentzVector n4M(0.,0.,0.,rnM); G4LorentzVector h4M(0.,0.,0.,npiM); G4double sum=rnM+npiM; if(fabs(reM-sum)Set4Momentum(n4M); curHadr->SetQPDG(G4QPDGCode(rnPDG)); // converted hypernucleus #ifdef pdebug G4cout<<"*G4QE::FSI:HypN "<A="<1) h4M=h4M/nPi; // split the 4-mom if necessary for(G4int ihn=0; ihnpush_back(thePion); // (delete equivalent for the pion) } if(rnPDG==90000002) // Additional action with curHadr change { G4LorentzVector newLV=n4M/2.; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90000001)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) //theFragments->push_back(secHadr); // (del.eq.- user is responsible for del) } else if(rnPDG==90002000) // Additional action with curHadr change { G4LorentzVector newLV=n4M/2.; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90001000)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) //theFragments->push_back(secHadr); // (del.eq.- user is responsible for del) } // @@ Add multybaryon decays if necessary } else if(reM>rnM-eps) // Decay in nonstrange and gamma { G4LorentzVector n4M(0.,0.,0.,rnM); G4LorentzVector h4M(0.,0.,0.,0.); G4double sum=rnM; if(fabs(reM-sum)Set4Momentum(n4M); curHadr->SetQPDG(G4QPDGCode(rnPDG)); // converted hypernucleus #ifdef pdebug G4cout<<"*G4QE::FSI:HypNg "<A="<Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90000001)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) //theFragments->push_back(secHadr); // (del.eq.- user is responsible for del) } else if(rnPDG==90002000) // Additional action with curHadr change { G4LorentzVector newLV=n4M/2.; curHadr->Set4Momentum(newLV); curHadr->SetQPDG(G4QPDGCode(90001000)); G4QHadron* secHadr = new G4QHadron(curHadr); theQHadrons.push_back(secHadr); // (del.eq.- user is responsible for del) //theFragments->push_back(secHadr); // (del.eq.- user is responsible for del) } // @@ Add multybaryon decays if necessary } #ifdef pdebug else // If this Error shows up (lowProbable appearance) => now it is left as is { G4cout<<"-Warning-G4QEnv::F:R="<First="<1) { f4M/=fS; // split the Hyperons 4-mom if necessary for(G4int il=1; ilSet4Momentum(f4M); curHadr->SetQPDG(G4QPDGCode(fPDG)); // converted hypernucleus if(sS>1) s4M/=sS; // split the Hyperon 4-mom if necessary for(G4int il=0; ilFirst="<1) { f4M/=fS; // split the Hyperons 4-mom if necessary for(G4int il=1; ilSet4Momentum(f4M); curHadr->SetQPDG(G4QPDGCode(fPDG)); // converted hypernucleus if(sS>1) s4M/=sS; // split the Hyperon 4-mom if necessary for(G4int il=0; ilB="<OUT,nH="<size()<size(); nHadr=theQHadrons.size(); G4int cfContSum=0; // Final ChargeControlSum G4int bfContSum=0; // Final BaryonNumberControlSum if(nHadr)for(unsigned f=0; fGet4Momentum()<GetPDGCode()<GetNFragments())) { G4int chg=curHadr->GetCharge(); cfContSum+=chg; G4int brn=curHadr->GetBaryonNumber(); bfContSum+=brn; G4int str=curHadr->GetStrangeness(); if ( brn > 1 && ( (!str && (chg == brn || !chg)) || (!chg && str == brn) ) ) // Check for multibaryon(split) { G4int bPDG=90000001; // Prototype: multineutron if (chg==brn) bPDG=90001000; // Multyproton else if(str==brn) bPDG=91000000; // Multilambda G4LorentzVector sp4m=curHadr->Get4Momentum()/brn; // Split 4-mom of the Multibaryon curHadr->Set4Momentum(sp4m); curHadr->SetQPDG(G4QPDGCode(bPDG)); // Substitute Multibaryon by a baryon for (G4int ib=1; ibpush_back(bH); //(del eq. - user is responsible for deleting) } } theFragments->push_back(curHadr); //(del eq. - user is responsible for deleting) } } #ifdef chdebug if(cfContSum-chContSum || bfContSum-bnContSum) { G4cerr<<"*G4QE::FSI:(9) Ch="< fill as it is"<GetCharge(); // Chsrge of the Baryon G4int theLS= qH->GetStrangeness(); // Strangness of the Baryon //G4int qPDG = qH->GetPDGCode(); // PDG Code of the decaying baryon G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the Baryon G4double qM = q4M.m(); // Mass of the Baryon #ifdef debug G4cout<<"G4QEnv::DecayBaryon: *Called* S="<mPPi) // Both n+pi0 (p=2/3) & p+Pi- (p=1/3) are possible { if(G4UniformRand()>.333333333) // Clebsh value for the Isobar decay { fQPDG=nQPDG; // Baryon is neutron (meson is Pi0) fMass=mNeut; } else { sQPDG=pimQPDG; // Meson is Pi- (baryon is proton) sMass=mPi; } } else // Only n+Pi0 decay is possible { fQPDG=nQPDG; // Baryon is neutron fMass=mNeut; } } else if(theLC==1) // Proton like: p+gam, p+Pi0, n+Pi+ are possible { if(qMmNPi) // Both p+pi0 (p=2/3) & n+Pi+ (p=1/3) are possible { if(G4UniformRand()<.333333333) // Clebsh value for the Isobar decay { fQPDG=nQPDG; // Baryon is neutron fMass=mNeut; sQPDG=pipQPDG; // Meson is Pi+ sMass=mPi; } // p+Pi0 is a default case } // p+Pi0 is a default case } else if(theLC==2) // Delta++ like: only p+PiP is possible { if(qM>mPPi) // Only p+gamma decay is possible { sQPDG=pipQPDG; // Meson is Pi+ (baryon is proton) sMass=mPi; } else // @@ Can be aReason to search for anError in Fragmentation { #ifdef debug G4cout<<"-Warning-G4QE::DecBary:*AsIs* DEL++ M="<mNPi) // Only p+gamma decay is possible { fQPDG=nQPDG; // Baryon is neutron fMass=mNeut; sQPDG=pimQPDG; // Meson is Pi- sMass=mPi; } else // @@ Can be aReason to search for anError in Fragmentation { #ifdef debug G4cout<<"-Warning-G4QE::DecBary:*AsIs* DEL++ M="<.6) // @@ Relative probability (take into account Phase Space) { fQPDG=szQPDG; // Baryon is Sigma0 fMass=mSigZ; } else { fQPDG=lQPDG; // Baryon is Lambda fMass=mLamb; } } else if(qM=mSpPi) // SigmaPlus+PiPlus decay is possible { fQPDG=spQPDG; // Baryon is SigmaP fMass=mSigP; sQPDG=pipQPDG; // Pi+ Meson sMass=mPi; } if(theLC==-2 && qM>=mSmPi)// SigmaPlus+PiPlus decay is possible { fQPDG=smQPDG; // Baryon is SigmaP fMass=mSigM; sQPDG=pimQPDG; // Pi- Meson sMass=mPi; } else { #ifdef debug G4cout<<"-Warning-G4QE::DecBary:*AsIs* Baryon(S=1,|C|=2),QC="<GetQC()<GetQC()< TotM="<GetQC(),qH->Get4Momentum()); if(!CheckGroundState(quasH,true)) theQHadrons.push_back(qH); // Cor or fill asItIs else delete qH; delete quasH; return; } else { delete qH; G4cerr<<"***G4QEnv::DecayBaryon: Can't Correct, *EmptyEnv*="<GetStrangeness(); // Strangness of the Nucleus if(theLS>=0) { G4cerr<<"***G4QEnvironment::DecayAntistrange: S="<Improve"<GetBaryonNumber(); // Baryon number of the Nucleus G4int theLC= qH->GetCharge(); // Chsrge of the Nucleus G4int qPDG = qH->GetPDGCode(); // PDG Code of the decaying Nucleus G4int K0PDG= qPDG+astr*999999; // Residual nonStrange nucleus for S*antiK0 G4QPDGCode K0QPDG(K0PDG); // QPDG of the nuclear residual for S*antiK0 G4double rK0M=K0QPDG.GetMass(); // Mass of the nuclear residual for S*antiK0 G4int KpPDG= qPDG+astr*999000; // Residual nonStrange nucleus for S*K+ G4QPDGCode KpQPDG(KpPDG); // QPDG of the nuclear residual for S*K+ G4double rKpM=KpQPDG.GetMass(); // Mass of the nuclear residual for S*K+ G4LorentzVector q4M = qH->Get4Momentum(); // Get 4-momentum of the Nucleus G4double qM = q4M.m(); // Mass of the Nucleus #ifdef debug G4cout<<"G4QE::DecayAntistrang:S="<qM) // Can not be K0 { if(astr*mK+rKpM>qM) // Can not be K+ too { #ifdef debug // @@ Survices, but... G4cout<<"*Warning*G4QEnvironment::DecayAntistrange: Too low mass, keep AsIs"<theLB-theLC) // Switch to K+ if Z>N { fQPDG = kpQPDG; // Positive Kaon fMass= mK; sQPDG = KpQPDG; // Residual nucleus to K+ sMass= rKpM; } G4double afMass=fMass; if(astr>1) afMass*=astr; G4LorentzVector f4Mom(0.,0.,0.,afMass); G4LorentzVector s4Mom(0.,0.,0.,sMass); G4double sum=afMass+sMass; if(fabs(qM-sum) TotM="<GetQC(); // Quark content of the Quasmon G4int resQPDG=valQ.GetSPDGCode(); // Reachable in a member function G4int resB=valQ.GetBaryonNumber(); // Baryon number of the Quasmon G4int resC=valQ.GetCharge(); // Charge of the Quasmon G4int resS=valQ.GetStrangeness(); // Strangeness of the Quasmon if(resQPDG==10 && resB>0) resQPDG=valQ.GetZNSPDGCode(); G4double resQMa=G4QPDGCode(resQPDG).GetMass();// GS Mass of the Residual Quasmon G4double resEMa=0.; // GS Mass of the EmptyResidualEnvironment G4bool bsCond=false; // FragSeparatCondition for QuasmonInVacuum G4LorentzVector enva4M=G4LorentzVector(0.,0.,0.,0.); // Prototype of the environment Mass G4QContent reTQC=valQ; // Prototype QuarkContent of the ResidNucl G4LorentzVector reTLV=quasm->Get4Momentum(); // Prototyoe 4-Mom of the Residual Nucleus G4double reTM=reTLV.m(); // Real mass of the Quasmon G4int envPDG=theEnvironment.GetPDG(); if ( resB > 1 && ( ( !resS && ( (resC == resB && reTM > resC*mProt) || (!resC && reTM > resB*mNeut) ) ) || (resS == resB && reTM > resS*mLamb) ) ) // Immediate Split(@@Decay) MultiBaryon { #ifdef chdebug G4cout<<"G4QE::CGS:*MultyBar*E="< "Can split neutron?" { G4QContent resQC=valQ-neutQC; // QC of Residual for the Neutron G4int resPDG=resQC.GetSPDGCode(); // PDG of Residual for the Neutron if(resPDG==10&&resQC.GetBaryonNumber()>0) resPDG=resQC.GetZNSPDGCode(); G4double resMas=G4QPDGCode(resPDG).GetMass(); // GS Mass of the Residual if(resMas+mNeut "Can split proton?" { G4QContent resQC=valQ-protQC; // QC of Residual for the Proton G4int resPDG=resQC.GetSPDGCode(); // PDG of Residual for the Proton if(resPDG==10&&resQC.GetBaryonNumber()>0) resPDG=resQC.GetZNSPDGCode(); G4double resMas=G4QPDGCode(resPDG).GetMass(); // GS Mass of the Residual if(resMas+mProt "Can split Lambda?" { G4QContent resQC=valQ-lambQC; // QC of Residual for the Lambda G4int resPDG=resQC.GetSPDGCode(); // PDG of Residual for the Lambda if(resPDG==10&&resQC.GetBaryonNumber()>0) resPDG=resQC.GetZNSPDGCode(); G4double resMas=G4QPDGCode(resPDG).GetMass(); // GS Mass of the Residual if(resMas+mLambresSMa && (resEMa || bsCond)) return true;// Why not ?? @@ (See G4Q the same) G4int nOfOUT = theQHadrons.size(); // Total #of QHadrons at this point #ifdef cdebug G4int reTPDG=reTQC.GetSPDGCode(); if(reTPDG==10&&reTQC.GetBaryonNumber()>0) reTPDG=reTQC.GetZNSPDGCode(); G4cout<<"G4QEnv::CheckGS:(tM="<0 & F="<rSM+hM="<resSMa+hadrMa &&!cNf && crPDG!=22) // Q(E)L contain QM(+EM)+lM ***Last CORR*** { if(resEMa) // => "NonVacuumEnvironment exists" case { G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa); // GS Mass of Quasmon G4LorentzVector enva4M = G4LorentzVector(0.,0.,0.,resEMa); // GS Mass of ResidEnvir if(tmpTM>=resQMa+resEMa+hadrMa && G4QHadron(tmpTLV).DecayIn3(hadr4M,quas4M,enva4M)) { //@@CHECK CoulBar (only for ResQuasmon in respect to ResEnv) and may be evaporate #ifdef debug G4cout<<"G4QE::CGS: Modify the Last 4-momentum to "<Set4Momentum(hadr4M); G4QHadron* quasH = new G4QHadron(valQ, quas4M); G4QContent theEQC=theEnvironment.GetQCZNS(); G4QHadron* envaH = new G4QHadron(theEQC,enva4M); #ifdef debug G4cout<<"G4QE::CGS: Fill Quasm "<Set4Momentum(enva4M); EvaporateResidual(envaH); // Try to evaporate residual environment(del.eq.) // Kill environment as it is already included in the decay theEnvironment=G4QNucleus(G4QContent(0,0,0,0,0,0), G4LorentzVector(0.,0.,0.,0.)); } else { #ifdef cdebug G4cout<<"***G4QEnv::CheckGroundState: Decay in Frag+ResQ+ResE error"< "Environ is vacuum" case (DecayIn2) { G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa); // GS Mass of Quasmon G4QHadron* quasH = new G4QHadron(valQ, quas4M); if(tmpTM>=resQMa+hadrMa && G4QHadron(tmpTLV).DecayIn2(hadr4M,quas4M))//DeIn2(noEnv) { //@@CHECK CoulBar (only for ResQuasmon in respect to ResEnv) & may be evaporate theLast->Set4Momentum(hadr4M); #ifdef debug G4cout<<"G4QE::CGS: modify theLast 4M="<Set4Momentum(quas4M); #ifdef debug G4cout<<"G4QE::CGS: fill newQH "<t+h+p="<tQMa+hadrMa+prevMa) { G4LorentzVector nuc4M = G4LorentzVector(0.,0.,0.,tQMa);// 4-Mom of ResidNucAtRest G4QHadron* nucH = new G4QHadron(reTQC, nuc4M); if(!G4QHadron(tmpTLV).DecayIn3(hadr4M,prev4M,nuc4M)) { delete nucH; // Delete "Residual Nucleus Hadron" G4cout<<"---Warning---G4QE::CGS:DecayIn ResNuc+LastH+PrevH Error"<Set4Momentum(hadr4M); thePrev->Set4Momentum(prev4M); nucH->Set4Momentum(nuc4M); #ifdef cdebug G4cout<<"G4QE::CGS:*SUCCESS**>CHECK, D4M="< Try to use any hadron from the output { #ifdef cdebug G4cout<<"G4QE::CGS:Prev&Last didn't help,nH="<2, MQ="<=0; id--) // Search for photons and pi+, and pi- { G4QHadron* curHadr = theQHadrons[id]; G4int hPDG=curHadr->GetPDGCode(); if(hPDG == 22) nphot=id; // Get the earliest gamma else if(hPDG==211 && npip<1) npip=id; // Get the latest Pi+ else if(hPDG==111) { G4double piE=curHadr->Get4Momentum().e(); #ifdef chdebug G4cout<<"G4QE::CheckGroundState:"<"<emaz) { npiz=id; emaz=piE; } } else if(hPDG==-211 && npim<1) npim=id; // Get the latest Pi- } if(nphot>=0) // Photon is found, try to use it to resolve PANIC { G4QHadron* curHadr = theQHadrons[nphot]; // Pointer to the photon G4LorentzVector ch4M=curHadr->Get4Momentum(); // 4-Mom of the Photon G4LorentzVector tt4M=ch4M+reTLV;// (resQMa=GSMass of the ResidQuasmon(+Env.)) G4double ttM=tt4M.m(); // Mass of the Phot+ResidQm compaund system if(resQMa1) fn4M/=rB; for(G4int ib=0; ibSet4Momentum(pi4M);// Change 4M of the Pion (reduced by decay) curHadr->SetQPDG(piQPDG); // Change Charge of thePion #ifdef cdebug if(resQPDG==89998005) G4cout<<"G4QE::CGS:1="<1&&npiz>=0&&(resC<-1||resC-resB>1)) { #ifdef chdebug G4cout<<"***G4QE::CGS:Pi0, rPDG="<(k+m)*k+(pi-) G4int piPD=-211; G4int nuPD=2112; G4double nuM=mNeut; if(resC!=-2) // k*(Delta++)+m*p+pi0->(k+m)*p+k*(pi+) { npi=resC-resB; piPD=211; nuPD=2212; nuM=mProt; } G4QPDGCode piQPDG(piPD); G4double suB=resB*nuM; // Total mass of nucleons G4double suM=npi*mPi; // Total mass of pions G4double sum=suB+suM; // Total mass of secondaries G4QHadron* curHadr = theQHadrons[npiz]; // Pointer to pi0 G4LorentzVector ch4M=curHadr->Get4Momentum(); // 4-Mom of the Pion G4LorentzVector tt4M=ch4M+reTLV;// (resQMa=GSMass of the ResidQuasmon(+Env.)) G4double ttM=tt4M.m(); // Mass of the Pion+ResidQm compaund system #ifdef chdebug G4cout<<"G4QE::CGS:sum="<1) pi4M/=npi; curHadr->Set4Momentum(pi4M);// Change 4M of the Pion (reduced by decay) curHadr->SetQPDG(piQPDG); // Change Charge of thePion if(npi>1) for(G4int ip=1; ip1) fn4M/=resB; for(G4int ib=0; ib Photons did not help, try to find an appropriate partner to join and decay G4int reTBN=reTQC.GetBaryonNumber(); // Baryon number of theHadronicState G4int reTCH=reTQC.GetCharge(); // Charge of theHadronicState G4bool isoN = reTCH-reTBN>0 || reTCH<0; // UnavoidableIsonucleus (Delta cond.) G4bool norN = reTCH<=reTBN || reTCH>=0; // "Regular nucleus" condition G4double nnM=resSMa; // Fake prototype of the NormalNucleusMass G4QContent ipiQC=pipQC; // Prototype of QCont for the Residual Pion+ G4QContent nnQC=reTQC+pimQC; // Prototype of theIsoReduceNucleus(Delta++) G4int nnPDG=nnQC.GetSPDGCode(); // Prot. PDGCode of the ResidNormalNucleus if((!nnPDG||nnPDG==10)&&nnQC.GetBaryonNumber()>0) nnPDG=nnQC.GetZNSPDGCode(); #ifdef cdebug G4cout<<"G4QE::CGS: nnPDR="<0) nnPDG=nnQC.GetZNSPDGCode(); } G4QPDGCode nnQPDG(nnPDG); // Now can even have Q-code ! if(nnPDG<80000000) nnM=nnQPDG.GetMass(); // Mass for the Fundamental Hadron else nnM=nnQPDG.GetNuclMass(nnPDG); // Mass for the Nucleus } G4bool chx2g=true; G4bool force=false; // Force-flag to initiate gamma decays (ChEx=>"Pi0"->2gamma) while(chx2g) { if(force) chx2g=false; for(G4int hd=nOfOUT-1; hd>=0; hd--)// Try to use any hadron to resolve PANIC { G4QHadron* curHadr = theQHadrons[hd]; G4int chNF=curHadr->GetNFragments(); G4int chCH=curHadr->GetCharge(); G4int chBN=curHadr->GetBaryonNumber(); //G4int chS=curHadr->GetStrangeness(); G4LorentzVector ch4M=curHadr->Get4Momentum(); // 4Mom of the Current Hadron #ifdef cdebug G4cout<<"G4QE::CGS:#"<2 gamma { if(!G4QHadron(tt4M).DecayIn3(ch4M,quas4M,gam4M)) { delete rqH; // Delete tmp "Residual Quasmon Hadron" #ifdef cdebug G4cout<<"***G4QE::CGS:DecayIn3 CurH+2Gammas,d="<"<Set4Momentum(quas4M); // Fill 4M of the GS Residual Quasmon rqH->SetQPDG(gamQPDG); // Change QPDG of the ResidualQuasmon theQHadrons.push_back(rqH); // Fill Gamma 1 as QHadron (del. eq.) #ifdef debug G4cout<<"G4QE::CGS:Fill (SubRQ) Gamma 1,(22)4M="<2 gammas) { if(!G4QHadron(tt4M).DecayIn2(ch4M,quas4M)) { delete rqH; // Delete tmp "Residual Quasmon Hadron" #ifdef cdebug G4cout<<"**G4QE::CGS:DecayIn2 CurH+ResQ d="<"<1) { G4QContent tcQC=chQC+nnQC; //QuarkCont for theTotalCompound nucleus G4QPDGCode tcQPDG(tcQC); // QPDG for the Total Compound G4double tcM=tcQPDG.GetMass(); // GS Mass of the TotalCompound G4QHadron* tcH = new G4QHadron(tcQPDG,tt4M);// Hadron=TotalCompound G4int tcS=tcQC.GetStrangeness(); // Total Strangeness G4int tcC=tcQC.GetCharge(); // Total Charge G4int tcBN=tcQC.GetBaryonNumber(); // Total Baryon Number //Try to decay in DiBaryon or MultyBaryon if ( tcBN == 2 || (!tcS && !tcC) || tcS == tcBN || tcC == tcBN ) { if(tcBN==2) theEnvironment.DecayDibaryon(tcH,&theQHadrons); // DB else theEnvironment.DecayMultyBaryon(tcH,&theQHadrons); // MB G4QHadron* theLast = theQHadrons[theQHadrons.size()-1]; curHadr->Set4Momentum(theLast->Get4Momentum());//4-Mom of CurHadr G4QPDGCode lQP=theLast->GetQPDG(); if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP); else curHadr->SetQC(theLast->GetQC()); theQHadrons.pop_back(); // theLastQHadron is excluded from OUTPUT delete theLast;//*!!When kill, delete theLastQHadr asAnInstance!* } else if(tcMSet4Momentum(tt4M);// Change 4M of the Current Hadron curHadr->SetQPDG(tcQPDG); // Change QPDG of theCurrentHadron #ifdef cdebug G4cout<<"**G4QE::CGS:DecayIn2 TotComp+gam d="<Set4Momentum(gc4M);//Change 4Mom of the Current Hadron curHadr->SetQPDG(gamQPDG);//Change PDG of theCurHadron to gamma tcH->Set4Momentum(tc4M); // Fill 4-Mom of the GS Total Compound #ifdef cdebug G4cout<<"G4QE::CGS:#"<"<Set4Momentum(tt4M); // Change 4Mom of the Current Hadron curHadr->SetQPDG(tcQPDG); // Change QPDG of the Current Hadron } return true; } // Check that the residual is a nucleus, not just a nucleon } // Check if pion can still be radiated or must be absorbed } // Check that the iso-exchange could help } // End of the IF: KINEMATIC CHECK FOR THE CURRENT HADRON(Isonuclear Case) else if(norN) { if(resSMa+chM1) { G4int hind=-1; if(npiz>-1) hind=npiz; else { if(resC+resC>resB && npim>-1) hind=npim; else hind=npip; } if(hind>-1) { G4QHadron* curHadr = theQHadrons[hind]; G4int chNF=curHadr->GetNFragments(); if(!chNF) { G4LorentzVector ch4M=curHadr->Get4Momentum(); // 4Mom of the Current Hadron G4LorentzVector tt4M=ch4M+reTLV; // resSMa=GSMass of the ResidQuasmon(+Env) G4double ttM=tt4M.m(); // TotalMass of CurHadr+Residual compaund G4QContent ttQC=valQ+curHadr->GetQC(); // total Quark Content G4int ttPDG=ttQC.GetZNSPDGCode(); G4QPDGCode ttQPDG(ttPDG); G4double ttGSM=ttQPDG.GetMass(); if(ttM>ttGSM && ttQC.GetStrangeness()>0)//Hypernucleus can be L->N degraded { ttPDG-=ttQC.GetStrangeness()*999999; // S Neutrons instead of S Lambdas ttQPDG=G4QPDGCode(ttPDG); // Update QPDGcode defining fragment ttGSM=ttQPDG.GetMass(); // Update the degraded mass value #ifdef debug G4cout<<"G4QEnv::CheckGS:Hypernucleus degraded to QPDG="<ttGSM) // Decay of Pi0 in 2 gammas with the residual nucleus { #ifdef cdebug G4cout<<"G4QEnv::CheckGS: Decay in ResQ="<Set4Momentum(fgam4M); // Change 4M of the Pion to Gamma curHadr->SetQPDG(gamQPDG); // Change QPDG of the Pion to Gamma G4QHadron* sgH = new G4QHadron(gamQPDG,sgam4M); // Gamma 2 theQHadrons.push_back(sgH); // Fill Gamma 2 as QHadron (del. eq.) G4QHadron* rqH = new G4QHadron(ttQPDG,quas4M); // GSResidualQuasmon theQHadrons.push_back(rqH); // Fill GSResidQuasmon as QHadron (del.eq.) return true; } } #ifdef debug else G4cout<<"-W-G4QEn::CheckGS:M="<size(); if(nHadrons) { for(G4int ih=0; ihoperator[](ih); // Pointer to the i-th QHadron G4QHadron* inH = (*HV)[ih]; // Pointer to the i-th QHadron G4int hNF = inH->GetNFragments(); // A#of secondary fragments if(!hNF) // Fill only final hadrons { #ifdef debug G4cout<<"G4QEnv::Copy&DeleteHV:#"<GetQC(); // Quark content of the Quasmon G4int resQPDG=valQ.GetSPDGCode(); // Reachable in a member function if(resQPDG==10&&valQ.GetBaryonNumber()>0) resQPDG=valQ.GetZNSPDGCode(); G4double resQMa=G4QPDGCode(resQPDG).GetMass(); // GS Mass of the Residual Quasmon G4LorentzVector enva4M=G4LorentzVector(0.,0.,0.,0.); G4LorentzVector reTLV=quasm->Get4Momentum(); // Prototyoe of the 4-Mom of the ResidNuc G4double resSMa=resQMa; // Prototype of MinSplitMass of ResidNucl G4int envPDG=theEnvironment.GetPDG(); // PDG Code of the Environment if(envPDG!=90000000) // => "Environment is not vacuum" case { G4double resEMa=theEnvironment.GetMZNS(); // GSMass of the Residual Environment enva4M=G4LorentzVector(0.,0.,0.,resEMa); // 4-Mom of the Residual Environment reTLV+=enva4M; // 4-Mom of Residual Nucleus resSMa+=resEMa; // Minimal Split Mass of Residual Nucleus G4double resTMa=reTLV.m(); // CM Mass of theResidNucleus (Quasm+Env) //#ifdef debug G4cout<<"G4QEnv::DecayInEnvQ: totM="< rQM+rEM="<resSMa) { G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa); // GS Mass of Quasmon G4QHadron* quasH = new G4QHadron(valQ, quas4M); G4QHadron* envaH = new G4QHadron(envPDG, enva4M); if(!G4QHadron(reTLV).DecayIn2(enva4M,quas4M)) { delete quasH; // Delete "Quasmon Hadron" delete envaH; // Delete "Environment Hadron" G4cout<<"---Warning---G4Q::DecInEnvQ:Decay in Environment+ResidualQuasmon"<Set4Momentum(quas4M); EvaporateResidual(quasH); // Try to evaporate quasmon (del. equiv.) envaH->Set4Momentum(enva4M); EvaporateResidual(envaH); // Try to evaporate residual (del. equiv.) } } else return false; } else return false; // => "Environment is vacuum" case return true; } // End of "DecayInEnvQ" // Add a Quasmon to the Environment void G4QEnvironment::AddQuasmon(G4Quasmon* Q) { // ======================================== theQuasmons.push_back(Q); totCharge+=Q->GetCharge(); totBaryoN+=Q->GetBaryonNumber(); tot4Mom +=Q->Get4Momentum(); #ifdef debug G4cout<<"G4QEnv::AddQuasmon:t4M="< now G4RandomDirection() is used //G4ThreeVector RndmDir() //{// -------- ========= // G4double x = G4UniformRand(), y = G4UniformRand(), z = G4UniformRand(); // G4double r2= x*x+y*y+z*z; // while(r2>1.||r2<.000001) // { // x = G4UniformRand(); y = G4UniformRand(); z = G4UniformRand(); // r2=x*x+y*y+z*z; // } // G4double r=sqrt(r2), quad=G4UniformRand(); // if(quad>0.5) // { // if(quad>0.75) // { // if(quad>0.875) return G4ThreeVector(-x/r,-y/r,-z/r); // else return G4ThreeVector(-x/r,-y/r, z/r); // } // else // { // if(quad>0.625) return G4ThreeVector(-x/r, y/r,-z/r); // else return G4ThreeVector(-x/r, y/r, z/r); // } // } // else // { // if(quad>0.25) // { // if(quad>0.375) return G4ThreeVector( x/r,-y/r,-z/r); // else return G4ThreeVector( x/r,-y/r, z/r); // } // else if(quad>0.125) return G4ThreeVector( x/r, y/r,-z/r); // } // return G4ThreeVector( x/r, y/r, z/r); //} // End of RndmDir