Loading...
Searching...
No Matches
G4QNucleus Class Reference
#include <G4QNucleus.hh>
Inheritance diagram for G4QNucleus:
Static Public Member Functions | |
| static void | SetParameters (G4double fN=.1, G4double fD=.05, G4double cP=4., G4double mR=1., G4double nD=.8 *CLHEP::fermi) |
Additional Inherited Members | |
| Protected Attributes inherited from G4QHadron | |
| G4LorentzVector | theMomentum |
Detailed Description
Definition at line 51 of file G4QNucleus.hh.
Constructor & Destructor Documentation
◆ G4QNucleus() [1/9]
| G4QNucleus::G4QNucleus | ( | ) |
Definition at line 65 of file G4QNucleus.cc.
65 : G4QHadron(), Z(0), N(0), S(0), dZ(0), dN(0), dS(0), maxClust(0),
66 theNucleons(),currentNucleon(-1),
67 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
68{
69 probVect[0]=mediRatio;
71#ifdef pardeb
74#endif
75}
G4DLLIMPORT std::ostream G4cout
Referenced by BindingEnergy(), CoulBarPenProb(), DecayAntiStrange(), EvaporateNucleus(), PrepareCandidates(), and SubtractNucleon().
◆ G4QNucleus() [2/9]
| G4QNucleus::G4QNucleus | ( | G4int | nucPDG | ) |
Definition at line 109 of file G4QNucleus.cc.
109 :
110 G4QHadron(nucPDG), maxClust(0), theNucleons(),
111 currentNucleon(-1), rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
112{
113 if(nucPDG==22) nucPDG=90000000;
114 InitByPDG(nucPDG);
116 Set4Momentum(p);
117#ifdef pardeb
120#endif
121}
Definition: LorentzVector.h:72
◆ G4QNucleus() [3/9]
| G4QNucleus::G4QNucleus | ( | G4LorentzVector | p, |
| G4int | nucPDG | ||
| ) |
Definition at line 123 of file G4QNucleus.cc.
123 :
124 G4QHadron(nucPDG, p), maxClust(0), theNucleons(),
125 currentNucleon(-1), rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
126{
127 InitByPDG(nucPDG);
128 Set4Momentum(p);
129#ifdef pardeb
132#endif
133}
◆ G4QNucleus() [4/9]
| G4QNucleus::G4QNucleus | ( | G4QContent | nucQC | ) |
Definition at line 155 of file G4QNucleus.cc.
155 :
156 G4QHadron(nucQC), dZ(0), dN(0), dS(0), maxClust(0), theNucleons(), currentNucleon(-1),
157 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
158{
160#ifdef debug
162#endif
163 probVect[0]=mediRatio;
170 Z = (b-S-du)/2; // protons
171 N = Z+du; // neutrons
172 SetQC(nucQC);
173#ifdef debug
175#endif
176 G4int nucPDG=90000000+S*1000000+Z*1000+N;
177 G4QPDGCode nQPDG(nucPDG);
178#ifdef debug
180#endif
181 G4double mass=nQPDG.GetNuclMass(Z,N,S);
182 if(nucPDG==90000000)
183 {
185 else mass=0.;
186 }
187#ifdef debug
189#endif
190 SetQPDG(nQPDG);
191#ifdef debug
193#endif
194 G4LorentzVector p(0.,0.,0.,mass);
195 Set4Momentum(p);
196 SetNFragments(0);
197#ifdef pardeb
200#endif
201}
Definition: G4QPDGCode.hh:142
◆ G4QNucleus() [5/9]
| G4QNucleus::G4QNucleus | ( | G4QContent | nucQC, |
| G4LorentzVector | p | ||
| ) |
Definition at line 203 of file G4QNucleus.cc.
203 :
204 G4QHadron(nucQC,p), dZ(0), dN(0), dS(0), maxClust(0), theNucleons(), currentNucleon(-1),
205 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
206{
207#ifdef debug
209#endif
210 probVect[0]=mediRatio;
212 Set4Momentum(p);
218 Z = (b-S-du)/2; // protons
219 N = Z+du; // neutrons
220 SetQC(nucQC);
221#ifdef debug
223#endif
224 G4QPDGCode nPDG(90000000+S*1000000+Z*1000+N);
225 SetQPDG(nPDG);
226 SetNFragments(0);
227#ifdef pardeb
230#endif
231}
◆ G4QNucleus() [6/9]
Definition at line 77 of file G4QNucleus.cc.
77 :
78 G4QHadron(90000000+s_value*1000000+z*1000+n), Z(z),N(n),S(s_value), dZ(0),dN(0),dS(0), maxClust(0),
79 theNucleons(), currentNucleon(-1), rho0(1.), radius(1.),
80 Tb(), TbActive(false), RhoActive(false)
81{
82 probVect[0]=mediRatio;
84#ifdef debug
86#endif
87 SetZNSQC(z,n,s_value);
89#ifdef debug
91#endif
92 G4double mass=nQPDG.GetNuclMass(Z,N,S);
93#ifdef debug
95#endif
97#ifdef debug
99#endif
100 G4LorentzVector p(0.,0.,0.,mass);
101 Set4Momentum(p);
102 SetNFragments(0);
103#ifdef debug
106#endif
107}
◆ G4QNucleus() [7/9]
| G4QNucleus::G4QNucleus | ( | G4int | z, |
| G4int | n, | ||
| G4int | s, | ||
| G4LorentzVector | p | ||
| ) |
Definition at line 135 of file G4QNucleus.cc.
135 :
136 G4QHadron(90000000+s_value*1000000+z*1000+n,p), Z(z),N(n),S(s_value), dZ(0),dN(0),dS(0), maxClust(0),
137 theNucleons(), currentNucleon(-1), rho0(1.),radius(1.),
138 Tb(), TbActive(false), RhoActive(false)
139{
140 probVect[0]=mediRatio;
142 Set4Momentum(p);
143 SetNFragments(0);
144 G4int ZNS=Z+N+S;
145 G4QPDGCode nPDG(90000000+S*1000000+Z*1000+N);
146 SetQPDG(nPDG);
147 G4QContent nQC(N+ZNS,Z+ZNS,S,0,0,0);
148 SetZNSQC(z,n,s_value);
149#ifdef pardeb
152#endif
153}
Definition: G4QContent.hh:53
◆ G4QNucleus() [8/9]
| G4QNucleus::G4QNucleus | ( | G4QNucleus * | right, |
| G4bool | cop3D = false |
||
| ) |
Definition at line 233 of file G4QNucleus.cc.
233 : currentNucleon(-1)
234{
235 Z = right->Z;
236 N = right->N;
237 S = right->S;
238 dZ = right->dZ;
239 dN = right->dN;
240 dS = right->dS;
241 maxClust = right->maxClust;
243 probVect[254] = right->probVect[254];
244 probVect[255] = right->probVect[255];
245 Tb = right->Tb;
246 TbActive = right->TbActive;
247 RhoActive = right->RhoActive;
252 rho0 = right->rho0;
253 radius = right->radius;
254 if(cop3D)
255 {
258 {
260 theNucleons.push_back(nucleon);
261 }
262 }
263#ifdef pardeb
266#endif
267}
Definition: G4QHadron.hh:53
◆ G4QNucleus() [9/9]
| G4QNucleus::G4QNucleus | ( | const G4QNucleus & | right, |
| G4bool | cop3D = false |
||
| ) |
Definition at line 269 of file G4QNucleus.cc.
269 :
270 G4QHadron(), currentNucleon(-1)
271{
272 Z = right.Z;
273 N = right.N;
274 S = right.S;
275 dZ = right.dZ;
276 dN = right.dN;
277 dS = right.dS;
278 maxClust = right.maxClust;
280 probVect[254] = right.probVect[254];
281 probVect[255] = right.probVect[255];
282 Tb = right.Tb;
283 TbActive = right.TbActive;
284 RhoActive = right.RhoActive;
289 rho0 = right.rho0;
290 radius = right.radius;
291 if(cop3D)
292 {
295 {
297 theNucleons.push_back(nucleon);
298 }
299 }
300#ifdef pardeb
303#endif
304}
◆ ~G4QNucleus()
| G4QNucleus::~G4QNucleus | ( | ) |
Definition at line 341 of file G4QNucleus.cc.
342{
343 for_each(theNucleons.begin(),theNucleons.end(),DeleteQHadron());
344}
Definition: G4QHadronVector.hh:57
Member Function Documentation
◆ ActivateBThickness()
| void G4QNucleus::ActivateBThickness | ( | ) |
Definition at line 3991 of file G4QNucleus.cc.
3992{
3997 // @@ make better approximation for light nuclei
4004 mT/=1000.; // Min b-thickness (@@ make 1000 a parameter)
4005 G4double b = 0.;
4006 while(T>mT)
4007 {
4008 //Tb->push_back(T); // Fill the thickness vector starting with b=0
4009 Tb.push_back(T); // Fill the thickness vector starting with b=0
4010 b+=db; // increment impact parameter
4012 T=C*E/(1.+E); // T(b) in fm^-2
4013 }
4014 TbActive=true; // Flag of activation
4015} // End of "ActivateBThickness"
Referenced by GetBThickness(), GetTbIntegral(), and GetThickness().
◆ BindingEnergy()
| G4double G4QNucleus::BindingEnergy | ( | const G4double & | cZ = 0, |
| const G4double & | cA = 0, |
||
| G4double | dZ = 0., |
||
| G4double | dA = 0. |
||
| ) |
Definition at line 3418 of file G4QNucleus.cc.
3420{
3427 G4int sZ=iZ;
3428 G4int sN=cN;
3429 if(delZ||dA)
3430 {
3433 GSM=G4QNucleus(Z-dz,N-dn,S).GetGSMass();
3434 sZ-=dz;
3435 sN-=dn;
3436 }
3438} // End of "BindingEnergy"
Referenced by Init3D().
◆ CalculateMass()
|
inline |
Definition at line 142 of file G4QNucleus.hh.
◆ ChooseFermiMomenta()
| void G4QNucleus::ChooseFermiMomenta | ( | ) |
Definition at line 3773 of file G4QNucleus.cc.
3774{
3777 //static const G4double mNeut= G4QPDGCode(2112).GetMass(); // Mass of neutron
3779 G4int i=0;
3785#ifdef debug
3787#endif
3788 for(i=0; i<theA; i++) // momenta for all, including the last, in case we swap nucleons
3789 {
3795 if( i == am1) dir=-sumMom.unit(); // try to compensate the mom noncons.
3798 {
3800 if(eMax>mProt) mom=sqrt(eMax*eMax - mProt2)*rn3*dir; // 3D proton momentum
3801#ifdef debug
3803#endif
3804 }
3805 else mom=ferm*rn3*dir; // 3-vector for the neutron momentum
3806 momentum[i]= mom;
3807 sumMom+= mom;
3808#ifdef debug
3810#endif
3811 }
3813 //G4double bindEn=BindingEnergy()/theA;
3817 G4double rMp2=rMp*rMp;
3818 G4double rMn2=rMn*rMn;
3819 //G4double rM=rMn;
3820 G4double rM2=rMn2;
3822#ifdef debug
3823 G4LorentzVector sum(0.,0.,0.,0.);
3824#endif
3825 for(i=0; i< theA ; i++ )
3826 {
3828 {
3829 //rM=rMp;
3830 rM2=rMp2;
3831 }
3832 else
3833 {
3834 //rM=rMn;
3835 rM2=rMn2;
3836 }
3837 G4ThreeVector curMom = momentum[i];
3839 G4LorentzVector tempV(curMom,energy);
3840#ifdef debug
3842 sum+=tempV;
3843#endif
3844 theNucleons[i]->Set4Momentum(tempV);
3845 }
3846#ifdef debug
3848#endif
3849 delete [] momentum;
3850} // End of ChooseFermiMomenta
G4DLLIMPORT std::ostream G4cerr
Definition: ThreeVector.h:41
double mag2() const
void SimpleSumReduction(G4ThreeVector *vectors, G4ThreeVector sum)
Definition: G4QNucleus.cc:3853
G4double GetFermiMomentum(G4double density)
Definition: G4QNucleus.cc:3765
G4double CoulombBarrier(const G4double &cZ=1, const G4double &cA=1, G4double dZ=0., G4double dA=0.)
Definition: G4QNucleus.cc:3386
Referenced by Init3D().
◆ ChooseImpactXandY()
Definition at line 3551 of file G4QNucleus.cc.
3552{
3553 G4double x=1.;
3554 G4double y=1.;
3555 do
3556 {
3557 x = G4UniformRand();
3558 x = x+x-1.;
3559 y = G4UniformRand();
3560 y = y+y-1.;
3561 }
3562 while(x*x+y*y > 1.);
3563 std::pair<G4double, G4double> theImpactParameter;
3564 theImpactParameter.first = x*maxImpact;
3565 theImpactParameter.second = y*maxImpact;
3566 return theImpactParameter;
3567} // End of "ChooseImpactXandY"
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ ChooseNucleons()
| void G4QNucleus::ChooseNucleons | ( | ) |
Definition at line 3570 of file G4QNucleus.cc.
3571{
3572#ifdef debug
3574#endif
3575 G4int protons =0;
3576 G4int nucleons=0;
3578 while (nucleons < theA)
3579 {
3581 {
3582 protons++;
3583 nucleons++;
3585 theNucleons.push_back(proton);
3586 }
3587 else if ( (nucleons-protons) < N )
3588 {
3589 nucleons++;
3591 theNucleons.push_back(neutron);
3592 }
3594 }
3595 return;
3596} // End of ChooseNucleons
Referenced by Init3D().
◆ ChoosePositions()
| void G4QNucleus::ChoosePositions | ( | ) |
Definition at line 3599 of file G4QNucleus.cc.
3600{
3604#ifdef debug
3606#endif
3622 while( i < theA && maxR > 0.) // LOOP over all nucleons
3623 {
3628#ifdef debug
3630#endif
3632 {
3633 // @@ Gaussian oscilator distribution is good only up to He4 (s-wave). Above: p-wave
3634 // (1+k*(r^2/R^2)]*exp[-r^2/R^2]. A=s+p=4+3*4=16 (M.K.) So Li,Be,C,N,O are wrong
3635 if(i==lastN) aPos=-rPos*sumPos.unit(); // TheLast tries toCompensate CenterOfGravity
3637 freeplace = true; // Imply that there is a free space
3639 {
3640 delta = places[j] - aPos; // Distance to nucleon j
3642 }
3643 // protons must at least have binding energy of CoulombBarrier (@@ ? M.K.), so
3644 // assuming Fermi Energy corresponds to Potential, we must place protons such
3645 // that the Fermi Energy > CoulombBarrier (?)
3646 // @@ M.K.: 1. CoulBar depends on aPos; 2. Makes Isotopic assymetry (!); 3. Perform.
3647 G4int nucPDG= theNucleons[i]->GetPDGCode();
3648#ifdef debug
3650#endif
3651 if(freeplace && nucPDG == 2212) // Free Space Protons
3652 {
3656 }
3657 if(rPos<mirPos)
3658 {
3659 mirPos=rPos;
3660 minPos=aPos;
3661 }
3662 if( freeplace || failCNT > maxCNT )
3663 {
3664 if( failCNT > maxCNT ) aPos=minPos;
3665#ifdef debug
3667#endif
3668 places[i]=aPos;
3669 sumPos+=aPos;
3670 ++i;
3671 failCNT=0;
3672 mirPos=maxR;
3673 }
3674 else ++failCNT;
3675 }
3676 }
3677#ifdef debug
3679#endif
3681#ifdef debug
3683#endif
3685 delete [] places;
3686#ifdef debug
3688#endif
3689} // End of ChoosePositions
G4double GetRadius(const G4double maxRelativeDenisty=0.5)
Definition: G4QNucleus.cc:3746
G4bool ReduceSum(G4ThreeVector *vectors, G4ThreeVector sum)
Definition: G4QNucleus.cc:3861
Referenced by Init3D().
◆ CoulBarPenProb()
| G4double G4QNucleus::CoulBarPenProb | ( | const G4double & | CB, |
| const G4double & | E, | ||
| const G4int & | C, | ||
| const G4int & | B | ||
| ) |
Definition at line 3441 of file G4QNucleus.cc.
3443{// = A.Lepretre et al, Nucl.Phys., A390 (1982) 221-239 =
3449 //static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0); // Mass of tritium
3450 //static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0); // Mass of Helium 3
3451 //static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0); // Mass of alpha
3453 // @@ --- Temporary 1 ---> close the OverBarrierReflection for all
3454 //return 1.;
3455 // ^^^^^^^---> End of Themporary 1
3456 // @@ --- Temporary 2 ---> close the OverBarrierReflection for fragments and mesons
3457 //if(E<CB) return 0.;
3458 //if(B!=1) return 1.;
3459 if(B<1 || B>2) return 1.;
3460 if(C>B+1)
3461 {
3462#ifdef debug
3464#endif
3465 return 1.;
3466 }
3467 // ^^^^^^^---> End of Themporary 2
3468 //G4double nA=GetA();
3469 //G4double nA=GetA()-B;
3470 //if(nA==40) G4cout<<"G4QN::CBPP:Z="<<GetZ()<<",C="<<C<<",B="<<B<<G4endl;
3471 //else return 1.; // @@@@@ Over barrier reflection is closed @@@ !!! @@@
3472 // Li6 C12 Al27
3473 //else if(nA<7||nA>8&&nA<12||nA>16&&nA<40) return 1.;// "OverBarrierReflection is closed"
3474 //else if(nA>8&&nA<12||nA>16&&nA<40) return 1.; // "OverBarrierReflection is closed" Cond
3475 //else if(nA<12||nA>16&&nA<40) return 1.; // "OverBarrierReflection is closed" Condition
3476 //else if(nA<12||nA>16) return 1.; // "OverBarrierReflection is closed" Condition
3477 //else if(nA<12) return 1.; // @@@@@ Over barrier reflection is closed @@@ !!! @@@
3478 //if(B+B>Z+N+S) return 1.;
3479 //G4double wD=wellDebth*B;
3480 G4double wD=wellDebth;
3481 //G4double wD=0.;
3482 // @@ --- Temporary 3 ---> close the OverBarrierReflection for mesons
3483 //if(!B) wD=0.;
3484 // ^^^^^^^---> End of Themporary 3
3486#ifdef debug
3488#endif
3489 if(2>3);
3490 // @@ Temporary "Mass Barrier for mesons" @@ __________________
3491 //else if(!B) wD=40.;
3492 // @@ End of Temporary^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3493 ////else if(nA<7&&B>0) wD=0.; // Only Coulomb Barrier can reflect !!!
3494 ////else if((nA<12||nA>16)&&B>0) wD=0.;// Only CoulombB can reflect !!! O16 E-dep of gamA
3495 ////else if((nA<12||nA>27)&&B>0) wD=0.;// Only CoulombB can reflect !!! O16 E-dep of gamA
3496 ////else if(nA<9&&B>0) return 1.;// Only CoulombBarrier can reflect !!! O16 E-dep of gamA
3497 ////else if(B>0) wD=0.; // Only Coulomb Barrier can reflect !!!
3498 ////else if(B==1) wD=0.;
3501 ////else if(B==1&&C==0) wD=0.;
3502 ////else if(B>1) return 1.;
3503 ////else if(B>1) wD=0.;
3504 ////else if(B==2) wD=0.;
3505 else if(B==2)
3506 {
3510 // @@ Temporary "Local B=2 Anti Virial factor" @@ __________________
3511 wD=80.; // 40 MeV per each nucleon
3512 //wD=wD/2;
3513 // @@ End of Temporary^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3514 }
3515 ////else if(B>2) wD=0.;
3516 ////else if(B>2) return 1.;
3517 ////else if(B==3) wD=0.;
3518 //else if(B==3&&C==1) wD=G4QNucleus(Z-1,N-2,S).GetGSMass()+mTrit-GSM;
3519 ////else if(B==3&&C==1) wD=0.;
3520 //else if(B==3&&C==2) wD=G4QNucleus(Z-2,N-1,S).GetGSMass()+mHel3-GSM;
3521 ////else if(B>3) wD=0.;
3522 ////else if(B==4) wD=0.;
3523 //else if(B==4&&C==2) wD=G4QNucleus(Z-2,N-2,S).GetGSMass()+mAlph-GSM;
3524 ////else if(B>4) wD=0.;
3525 ////else if(B>4) return 1.;
3526 //else if(B>4)wD=G4QNucleus(Z-C,N-B+C,S).GetGSMass()+G4QNucleus(C,B-C,S).GetGSMass()-GSM;
3527 if(wD<0.) wD=0.;
3528#ifdef debug
3530#endif
3531 // @@ Temporary "Virial factor" @@ __________________
3532 wD=wD+wD;
3533 // @@ End of Temporary^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3534 G4double sR=0.;
3535 G4double CBD=CB+wD;
3536 G4double ED=E+wD;
3537 if(CBD<0.) return 1.;
3538 if(ED<=0.) return 0.;
3539 //if(nA<27) sR=sqrt(wD/ED);
3540 //else sR=sqrt(CBD/ED);
3541 sR=sqrt(CBD/ED);
3542 //sR=sqrt(wD/ED);
3543#ifdef debug
3545#endif
3546 if(sR>=1.) return 0.;
3547 return 1.-sR*sR*sR;
3548} // End of "CoulBarPenProb"
Referenced by EvaporateBaryon().
◆ CoulombBarGen()
| G4double G4QNucleus::CoulombBarGen | ( | const G4double & | rZ, |
| const G4double & | rA, | ||
| const G4double & | cZ, | ||
| const G4double & | cA | ||
| ) |
Definition at line 3358 of file G4QNucleus.cc.
3360{
3363 if(cA < 0.) ca=-cA;
3364#ifdef debug
3366#endif
3368 if(rA < 0.) ra=-rA;
3369 G4double zz=rZ*cZ;
3370 // Naitive CHIPS radius: CB={1.46=200(MeV)/137}*z*Z/{R=1.13}*((a*z)**1/3+A**1/3) (?)
3371 //G4double cb=1.29*zz/(pow(rA,third)+pow(cA,third));
3372 //double cb=zz/(pow(rA,third)+pow(ca,third)+.1);
3374 // Geant4 solution for protons is practically the same:
3375 // G4double cb=1.263*Z/(1.0 + pow(rA,third));
3376 // @@ --- Temporary "Lambda/Delta barrier for mesons"
3377 //if(!cA) cb+=40.;
3378 // --- End of Temporary ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3379#ifdef debug
3381#endif
3382 return cb;
3383} // End of "CoulombBarier"
Referenced by CoulombBarrier(), and G4QLowEnergy::PostStepDoIt().
◆ CoulombBarrier()
| G4double G4QNucleus::CoulombBarrier | ( | const G4double & | cZ = 1, |
| const G4double & | cA = 1, |
||
| G4double | dZ = 0., |
||
| G4double | dA = 0. |
||
| ) |
Definition at line 3386 of file G4QNucleus.cc.
3388{
3390 if (dA != 0.) rA-=dA;
3391#ifdef debug
3393#endif
3394 G4double rZ=Z-cZ;
3395 if(delZ != 0.) rZ-=delZ;
3396#ifdef debug
3398#endif
3400} // End of "CoulombBarier"
G4double CoulombBarGen(const G4double &rZ, const G4double &rA, const G4double &cZ, const G4double &cA)
Definition: G4QNucleus.cc:3358
Referenced by ChooseFermiMomenta(), ChoosePositions(), EvaporateBaryon(), EvaporateNucleus(), G4QFragmentation::Fragment(), G4QFragmentation::G4QFragmentation(), and SplitBaryon().
◆ DecayAlphaAlpha()
| void G4QNucleus::DecayAlphaAlpha | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 7665 of file G4QNucleus.cc.
7666{
7671 if(qPDG!=90004004)
7672 {
7673 delete qH;
7674 // G4cerr<<"***G4QNucleus::DecayAlphaAlpha: qPDG="<<qPDG<<G4endl;
7675 // throw G4QException("***G4QNucleus::DecayAlphaAlpha: Not Be8 state decais in 2 alphas");
7676 G4ExceptionDescription ed;
7679 FatalException, ed);
7680 }
7683#ifdef debug
7685#endif
7686 //if(qM>aaGSM+.01) // @@ Be8*->gamma+Be8 (as in evaporation) @@ gamma cooling
7687 if(2>3)
7688 {
7689 G4int fPDG = 22;
7690 G4int sPDG = 90004004;
7691 G4double fMass= 0.;
7692 G4double sMass= aaGSM;
7693 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7695 G4double sum=fMass+sMass;
7696 if(fabs(qM-sum)<eps)
7697 {
7698 f4Mom=q4M*(fMass/sum);
7699 s4Mom=q4M*(sMass/sum);
7700 }
7702 {
7705 //G4cerr<<"***G4QNuc::DecayAlphAlph: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7706 //throw G4QException("G4QNucleus::DecayAlphaAlpha:g+diAlph DecayIn2 didn't succeed");
7707 evaHV->push_back(qH);
7708 return;
7709 }
7710#ifdef debug
7712#endif
7714 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7715 qH->Set4Momentum(s4Mom);
7716 q4M=s4Mom;
7717 }
7718 G4int fPDG = 90002002;
7719 G4int sPDG = 90002002;
7720 G4double fMass= mAlph;
7721 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7722 G4LorentzVector s4Mom(0.,0.,0.,fMass);
7723 G4double sum=fMass+fMass;
7724 if(fabs(qM-sum)<eps)
7725 {
7726 f4Mom=q4M*(fMass/sum);
7727 s4Mom=f4Mom;
7728 }
7730 {
7733 //G4cerr<<"***G4QNuc::DecayAlphAlph: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7734 //throw G4QException("G4QNucleus::DecayAlphaAlpha: diAlpha DecayIn2 didn't succeed");
7735 evaHV->push_back(qH);
7736 return;
7737 }
7738#ifdef debug
7740#endif
7741 delete qH;
7743 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7745 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7746} // End of DecayAlphaAlpha
double m() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
Referenced by EvaporateNucleus().
◆ DecayAlphaBar()
| void G4QNucleus::DecayAlphaBar | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 7414 of file G4QNucleus.cc.
7415{
7427#ifdef debug
7429#endif
7433
7434 if ( ( (!totS && !totC) || totC == totBN || totS == totBN)
7435 && totBN > 1) DecayMultyBaryon(qH,evaHV);
7436 else if(qPDG==92001002||qPDG==92002001||qPDG==91003001||qPDG==91001003||qPDG==93001001)
7437 evaHV->push_back(qH);
7438 else if(qPDG==92000003||qPDG==92003000||qPDG==93000002||qPDG==93002000)
7439 {
7441 G4double fMass= mLamb;
7442 G4int sPDG = 2112;
7443 G4double sMass= mNeut;
7444 if (qPDG==92003000) // "2L+3p" case
7445 {
7446 sPDG = 2212;
7447 sMass= mProt;
7448 }
7449 else if(qPDG==93000002) // "2n+3L" case
7450 {
7451 fPDG = 2112;
7452 fMass= mNeut;
7453 sPDG = 3122;
7454 sMass= mLamb;
7455 }
7456 else if(qPDG==93002000) // "2p+3L" case
7457 {
7458 fPDG = 2212;
7459 fMass= mProt;
7460 sPDG = 3122;
7461 sMass= mLamb;
7462 }
7463 G4double tfM=fMass+fMass;
7464 G4double tsM=sMass+sMass+sMass;
7465 G4LorentzVector f4Mom(0.,0.,0.,tfM);
7466 G4LorentzVector s4Mom(0.,0.,0.,tsM);
7467 G4double sum=tfM+tsM;
7468 if(fabs(qM-sum)<eps)
7469 {
7470 f4Mom=q4M*(tfM/sum);
7471 s4Mom=q4M*(tsM/sum);
7472 }
7474 {
7477 //G4cerr<<"***G4QN::DecayAlphaBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7478 //throw G4QException("G4QNucleus::DecayAlphaBar: DecayIn2 didn't succeed for 3/2");
7479 evaHV->push_back(qH);
7480 return;
7481 }
7482#ifdef debug
7484#endif
7485 delete qH;
7486 G4LorentzVector rf4Mom=f4Mom/2;
7488 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7489 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7490 G4LorentzVector rs4Mom=s4Mom/3;
7492 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7493 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7494 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7495 }
7496 else if(qPDG==90004001||qPDG==90001004)
7497 {
7499 G4double fMass= mHe3;
7500 G4int sPDG = 2212;
7501 G4double sMass= mProt;
7502 if (qPDG==90001004) // "t+2n" case
7503 {
7504 fPDG = 90001002;
7505 fMass= mTrit;
7506 sPDG = 2112;
7507 sMass= mNeut;
7508 }
7509 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7510 G4LorentzVector s4Mom(0.,0.,0.,sMass);
7511 G4LorentzVector t4Mom(0.,0.,0.,sMass);
7512 G4double sum=fMass+sMass+sMass;
7513 if(fabs(qM-sum)<eps)
7514 {
7515 f4Mom=q4M*(fMass/sum);
7516 s4Mom=q4M*(sMass/sum);
7517 t4Mom=s4Mom;
7518 }
7520 {
7523 //G4cerr<<"*G4QNuc::DecayAlphaBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7524 //throw G4QException("G4QNucleus::DecayAlphaBar: t/nn,He3/pp DecayIn3 didn't");
7525 evaHV->push_back(qH);
7526 return;
7527 }
7528#ifdef debug
7530#endif
7531 delete qH;
7533 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7535 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7537 evaHV->push_back(H3); // Fill "H3" (delete equivalent)
7538 }
7539 else if(qPDG==94000001||qPDG==94001000||qPDG==91000004||qPDG==91004000)
7540 {
7542 G4double fMass= mLamb+mLamb;
7543 G4int sPDG = 2112;
7544 G4double sMass= mNeut;
7545 if (qPDG==94001000) // "4L+p" case
7546 {
7547 sPDG = 2212;
7548 sMass= mProt;
7549 }
7550 else if(qPDG==91000004) // "4n+L" case
7551 {
7552 fPDG = 2112;
7553 fMass= mNeut+mNeut;
7554 sPDG = 3122;
7555 sMass= mLamb;
7556 }
7557 else if(qPDG==91004000) // "4p+L" case
7558 {
7559 fPDG = 2212;
7560 fMass= mProt+mProt;
7561 sPDG = 3122;
7562 sMass= mLamb;
7563 }
7564 G4LorentzVector f4Mom(0.,0.,0.,fMass+fMass);
7565 G4LorentzVector s4Mom(0.,0.,0.,sMass);
7566 G4double sum=fMass+fMass+sMass;
7567 if(fabs(qM-sum)<eps)
7568 {
7569 f4Mom=q4M*((fMass+fMass)/sum);
7570 s4Mom=q4M*(sMass/sum);
7571 }
7573 {
7576 //G4cerr<<"*G4QNuc::DecayAlphaBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7577 //throw G4QException("G4QNucl::DecayAlphaBar:QuintBaryon DecayIn2 didn't succeed");
7578 evaHV->push_back(qH);
7579 return;
7580 }
7581#ifdef debug
7583#endif
7584 delete qH;
7585 G4LorentzVector rf4Mom=f4Mom/4;
7587 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7588 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7589 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7590 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7592 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7593 }
7594 else if(qPDG==90003002||qPDG==90002003||qPDG==91002002)
7595 {
7597 G4int sPDG = 2112;
7598 G4double fMass= mAlph;
7599 G4double sMass= mNeut;
7600 if(qPDG==90003002) // "alpha+p" case
7601 {
7602 sPDG = 2212;
7603 sMass= mProt;
7604 }
7605 else if(qPDG==9100202) // "alpha+l" case
7606 {
7607 sPDG = 3122;
7608 sMass= mLamb;
7609 }
7610 else if(qPDG!=90002003)
7611 {
7612 evaHV->push_back(qH); // Fill hadron as it is (delete equivalent)
7613 //EvaporateNucleus(qH, evaHV); // Evaporate Nucleus (delete equivivalent)
7614 return;
7615 }
7616 G4double dM=fMass+sMass-qM;
7617 if(dM>0.&&dM<1.)
7618 {
7619#ifdef debug
7622#endif
7623 G4double hdM=dM/2;
7624 fMass-=hdM;
7625 sMass-=hdM;
7626 }
7627 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7629 G4double sum=fMass+sMass;
7630 if(fabs(qM-sum)<eps)
7631 {
7632 f4Mom=q4M*(fMass/sum);
7633 s4Mom=q4M*(sMass/sum);
7634 }
7636 {
7639 //G4cout<<"*G4QNuc::DecayAlphaBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7640 //throw G4QException("***G4QNucl::DecayAlphaBar:Alpha+Baryon DecIn2 didn't succeed");
7641 evaHV->push_back(qH);
7642 return;
7643 }
7644#ifdef debug
7646#endif
7647 delete qH;
7649 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7651 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7652 }
7654#ifdef qdebug
7655 if (qH)
7656 {
7658 << qH->GetPDGCode() << G4endl;
7659 delete qH;
7660 }
7661#endif
7662} // End of DecayAlphaBar
void DecayMultyBaryon(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:7190
Referenced by EvaporateNucleus().
◆ DecayAlphaDiN()
| void G4QNucleus::DecayAlphaDiN | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 7331 of file G4QNucleus.cc.
7332{
7341#ifdef debug
7343#endif
7345 G4double fMass= mProt;
7346 G4int sPDG = 90002002;
7347 G4double sMass= mAlph;
7348 if (qPDG==90002004) // "alpha+2neutrons" case
7349 {
7350 if(fabs(qM-mHel6)<eps)
7351 {
7352 evaHV->push_back(qH); // Fill as it is (delete equivalent)
7353 return;
7354 }
7355 else if(mNeut+mNeut+mAlph<qM)
7356 {
7357 fPDG = 2112;
7358 fMass= mNeut;
7359 }
7360 else
7361 {
7362 delete qH;
7363 // G4cerr<<"***G4QNu::DecAlDiN:M(He6="<<mHel6<<")="<<qM<<"<"<<mNeut+mNeut+mAlph<<G4endl;
7364 // throw G4QException("G4QNuc::DecayAlphaDiN: Cannot decay excited He6 with this mass");
7365 G4ExceptionDescription ed;
7366 ed << "Cannot decay excited He6 with this mass: M(He6=" << mHel6 << ")="
7369 FatalException, ed);
7370 }
7371 }
7372 else if(qPDG!=90004002) // "Bad call" case
7373 {
7374 delete qH;
7375 // G4cerr<<"***G4QNuc::DecayAlphaDiN: PDG="<<qPDG<<G4endl;
7376 // throw G4QException("G4QNuc::DecayAlphaDiN: Can not decay this PDG Code");
7377 G4ExceptionDescription ed;
7380 FatalException, ed);
7381 }
7382 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7383 G4LorentzVector s4Mom(0.,0.,0.,fMass);
7384 G4LorentzVector t4Mom(0.,0.,0.,sMass);
7385 G4double sum=fMass+fMass+sMass;
7386 if(fabs(qM-sum)<eps)
7387 {
7388 f4Mom=q4M*(fMass/sum);
7389 s4Mom=f4Mom;
7390 t4Mom=q4M*(sMass/sum);
7391 }
7393 {
7396 //G4cerr<<"***G4QNuc::DecayAlphaDiN: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7397 //throw G4QException("G4QNuc::DecayAlphaDiN: Alpha+N+N DecayIn3 error");
7398 evaHV->push_back(qH);
7399 return;
7400 }
7401#ifdef debug
7403#endif
7404 delete qH;
7406 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7408 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7410 evaHV->push_back(H3); // Fill "H3" (delete equivalent)
7411} // End of DecayAlphaDiN
Referenced by EvaporateNucleus().
◆ DecayAntiDibaryon()
| void G4QNucleus::DecayAntiDibaryon | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 6591 of file G4QNucleus.cc.
6592{
6626#ifdef debug
6628#endif
6629 // Select a chanel of the dibaryon decay (including Delta+Delta-> 4 particle decay
6631 G4int sPDG = -2212;
6633 G4double fMass= mProt;
6634 G4double sMass= mProt;
6635 G4double tMass= mPi;
6636 if (qPDG==89996002 && rM>=dmPiP) // "anti-diDelta++" case
6637 {
6638 sPDG = -211;
6639 sMass= mPi;
6640 four = true;
6641 }
6642 else if(qPDG==90001996 && rM>=dmPiN) // "diDelta--" case
6643 {
6644 sPDG = 211;
6645 fPDG = -2112;
6646 sMass= mPi;
6647 fMass= mNeut;
6648 four = true;
6649 }
6650 else if(qPDG==89999998 && rM>=dNeut) // "dineutron" case
6651 {
6652 fPDG = -2112;
6653 sPDG = -2112;
6654 fMass= mNeut;
6655 sMass= mNeut;
6656 }
6657 else if(qPDG==89998999 && rM>=mDeut) // "exited deutron" case
6658 {
6659 if(fabs(qM-mDeut)<eps)
6660 {
6661 evaHV->push_back(qH); // Fill as it is (delete equivalent)
6662 return;
6663 }
6664 else if(mProt+mNeut<rM)
6665 {
6666 fPDG = -2112;
6667 fMass= mNeut;
6668 }
6669 else
6670 {
6671 fPDG = 22;
6672 sPDG = 89998999; // Anti-deuteron
6673 fMass= 0.;
6674 sMass= mDeut;
6676 }
6677 }
6678 else if(qPDG==88999999 && rM>=dLaNe) // "Lambda-neutron" case
6679 {
6680 fPDG = -2112;
6681 sPDG = -3122;
6682 fMass= mNeut;
6683 sMass= mLamb;
6684 }
6685 else if(qPDG==88999999 && rM>=dLaPr) // "Lambda-proton" case
6686 {
6687 sPDG = -3122;
6688 sMass= mLamb;
6689 }
6690 else if(qPDG==90000997 && rM>=nnPi) // "neutron/Delta-" case
6691 {
6692 fPDG = -2112;
6693 sPDG = -2112;
6694 tPDG = 211;
6695 fMass= mNeut;
6696 sMass= mNeut;
6697 }
6698 else if(qPDG==89997001 && rM>=ppPi) // "proton/Delta++" case
6699 {
6700 tPDG = -211;
6701 }
6702 else if(qPDG==89000998 && rM>=lnPi) // "lambda/Delta-" case
6703 {
6704 fPDG = -2112;
6705 sPDG = -3122;
6706 tPDG = 211;
6707 fMass= mNeut;
6708 sMass= mLamb;
6709 }
6710 else if(qPDG==889998001 && rM>=lpPi) // "lambda/Delta+" case
6711 {
6712 sPDG = -3122;
6713 tPDG = -211;
6714 sMass= mLamb;
6715 }
6716 else if(qPDG==89000998 && rM>=dSiNe) // "Sigma-/neutron" case
6717 {
6718 fPDG = -2112;
6719 sPDG = -3112;
6720 fMass= mNeut;
6721 sMass= mSigM;
6722 }
6723 else if(qPDG==88998001 && rM>=dSiPr) // "Sigma+/proton" case
6724 {
6725 sPDG = -3222;
6726 sMass= mSigP;
6727 }
6728 else if(qPDG==88000000 && rM>=dLamb) // "diLambda" case
6729 {
6730 fPDG = -3122;
6731 sPDG = -3122;
6732 fMass= mLamb;
6733 sMass= mLamb;
6734 }
6735 else if(qPDG==88000999 && rM>=dKsNe) // "Ksi-/neutron" case
6736 {
6737 fPDG = -2112;
6738 sPDG = -3312;
6739 fMass= mNeut;
6740 sMass= mKsiM;
6741 }
6742 else if(qPDG==87999001 && rM>=dKsPr) // "Ksi0/proton" case
6743 {
6744 sPDG = -3322;
6745 sMass= mKsiZ;
6746 }
6747 else if(qPDG!=89998000|| rM<dProt) // Other possibilities (if not a default)
6748 {
6749 G4int qS = qH->GetStrangeness();
6750 G4int qB = qH->GetBaryonNumber();
6751 if(qB>0&&qS<0) // Antistrange diBarion
6752 {
6753 DecayAntiStrange(qH,evaHV);
6754 return;
6755 }
6756 else
6757 {
6758 delete qH;
6761 // @@ Nothing to do. Just 2 GeV disappears... Very rare! Just to avoid the exception.
6762 //throw G4QException("G4QNucleus::DecayDibar: Unknown PDG code or small Mass of DB");
6763 }
6764 }
6765 G4LorentzVector f4Mom(0.,0.,0.,fMass);
6766 G4LorentzVector s4Mom(0.,0.,0.,sMass);
6767 G4LorentzVector t4Mom(0.,0.,0.,tMass);
6768 if(!tPDG&&!four)
6769 {
6770 G4double sum=fMass+sMass;
6771 if(fabs(qM-sum)<eps)
6772 {
6773 f4Mom=q4M*(fMass/sum);
6774 s4Mom=q4M*(sMass/sum);
6775 }
6777 {
6780 //G4cerr<<"***G4QNucl::DecayDiBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6781 //throw G4QException("***G4QNucleus::DecayDibaryon: DiBaryon DecayIn2 error");
6782 evaHV->push_back(qH);
6783 return;
6784 }
6785#ifdef debug
6788#endif
6789 delete qH;
6791 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6793 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6794 }
6795 else if(four)
6796 {
6797 q4M=q4M/2.; // Divided in 2 !!!
6798 qM/=2.; // Divide the mass in 2 !
6799 G4double sum=fMass+sMass;
6800 if(fabs(qM-sum)<eps)
6801 {
6802 f4Mom=q4M*(fMass/sum);
6803 s4Mom=q4M*(sMass/sum);
6804 }
6806 {
6809 //G4cerr<<"***G4QN::DecayDibaryon: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6810 //throw G4QException("***G4QNucleus::DecDibaryon: Dibaryon DecayIn2 error");
6811 evaHV->push_back(qH);
6812 return;
6813 }
6814#ifdef debug
6817#endif
6819 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6821 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6822 // Now the second pair mus be decayed
6823 if(fabs(qM-sum)<eps)
6824 {
6825 f4Mom=q4M*(fMass/sum);
6826 s4Mom=q4M*(sMass/sum);
6827 }
6829 {
6830 // Should not be here as sum was already compared with qM above for the first delta
6831 delete qH;
6832 // G4cerr<<"**G4QNucl::DecAntiDibar:fPDG="<<fPDG<<"(fM="<<fMass<<")+sPDG="<<sPDG<<"(sM="
6833 // <<sMass<<")="<<sum<<" >? (DD2,Can't be here) TotM="<<q4M.m()<<q4M<<G4endl;
6834 // throw G4QException("G4QNucleus::DecayAntiDibaryon: General DecayIn2 error");
6835 G4ExceptionDescription ed;
6836 ed << " General DecayIn2 error: fPDG=" << fPDG << "(fM=" << fMass
6837 << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
6840 FatalException, ed);
6841 }
6842#ifdef debug
6845#endif
6847 evaHV->push_back(H3); // Fill "H1" (delete equivalent)
6849 evaHV->push_back(H4); // Fill "H2" (delete equivalent)
6850 delete qH;
6851 }
6852 else
6853 {
6854 G4double sum=fMass+sMass+tMass;
6855 if(fabs(qM-sum)<eps)
6856 {
6857 f4Mom=q4M*(fMass/sum);
6858 s4Mom=q4M*(sMass/sum);
6859 t4Mom=q4M*(tMass/sum);
6860 }
6862 {
6865 //G4cerr<<"***G4QNuc::DecayDibaryon:qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6866 //throw G4QException("G4QNucleus::DecayDibaryon: diBar DecayIn3 error");
6867 evaHV->push_back(qH);
6868 return;
6869 }
6870#ifdef debug
6873#endif
6874 //qH->SetNFragments(2); // Fill a#of fragments to decaying Dibaryon
6875 //evaHV->push_back(qH); // Fill hadron with nf=2 (delete equivalent)
6876 // Instead
6877 delete qH;
6878 //
6880 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6882 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6884 evaHV->push_back(H3); // Fill "H3" (delete equivalent)
6885 }
6886#ifdef qdebug
6887 if (qH)
6888 {
6890 delete qH;
6891 }
6892#endif
6893} // End of DecayAntiDibaryon
Definition: G4QNucleus.hh:52
void DecayAntiStrange(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:6896
Referenced by EvaporateNucleus().
◆ DecayAntiStrange()
| void G4QNucleus::DecayAntiStrange | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 6896 of file G4QNucleus.cc.
6897{
6907#ifdef debug
6909#endif
6911 if(qS>=0 || qB<1)
6912 {
6913 delete qH;
6914 // G4cerr<<"G4QNuc::DecayAntiStrange:QC="<<qQC<<",S="<<qS<<",B="<<qB<<",4M="<<q4M<<G4endl;
6915 // throw G4QException("G4QNucleus::DecayAntiStrange: not an Anti Strange Nucleus");
6916 G4ExceptionDescription ed;
6917 ed << "not an Anti Strange Nucleus: QC=" << qQC << ",S=" << qS << ",B="
6920 FatalException, ed);
6921 }
6923 G4double k1M=mK0;
6926 G4double k2M=mK;
6931 if(qP>0 && qP>qN) // a#of protons > a#of neutrons
6932 {
6933 if(qP>=bH) // => "Enough protons in nucleus" case
6934 {
6935 if(qN>=sH)
6936 {
6937 n1=sH;
6938 n2=bH;
6939 }
6940 else
6941 {
6942 G4int dPN=qP-qN;
6943 if(dPN>=aS)
6944 {
6945 n1=0;
6946 n2=aS;
6947 }
6948 else
6949 {
6950 G4int sS=(aS-dPN)/2;
6951 G4int bS=aS-dPN-sS;
6952 sS+=dPN;
6953 if(qP>=sS && qN>=bS)
6954 {
6955 n1=bS;
6956 n2=sS;
6957 }
6958 else if(qP<sS)
6959 {
6960 n1=aS-qP;
6961 n2=qP;
6962 }
6963 else
6964 {
6965 n1=qN;
6966 n2=aS-qN;
6967 }
6968 }
6969 }
6970 }
6971 }
6972 else if(qN>=bH)
6973 {
6974 if(qP>=sH)
6975 {
6976 n2=sH;
6977 n1=bH;
6978 }
6979 else
6980 {
6981 G4int dNP=qN-qP;
6982 if(dNP>=aS)
6983 {
6984 n1=aS;
6985 n2=0;
6986 }
6987 else
6988 {
6989 G4int sS=(aS-dNP)/2;
6990 G4int bS=aS-sS;
6991 if(qN>=bS && qP>=sS)
6992 {
6993 n1=bS;
6994 n2=sS;
6995 }
6996 else if(qN<bS)
6997 {
6998 n1=qN;
6999 n2=aS-qN;
7000 }
7001 else
7002 {
7003 n1=aS-qP;
7004 n2=qP;
7005 }
7006 }
7007 }
7008 }
7011#ifdef debug
7013#endif
7016 if(qP>=-qS) m2_value=-qS; // Enough charge for K+'s
7017 else if(qP>0) m1=-qS-qP; // Anti-Lambdas are partially compensated by neutrons
7020 if(mucM+m1*mK+m2_value*mK0<nucM+n1*mK+n2*mK0) // New is smaller
7021 {
7022 qPDG=sPDG;
7023 nucM=mucM;
7024 n1=m1;
7025 n2=m2_value;
7026 }
7027#ifdef debug
7029#endif
7030 if(!n1||!n2) // AntiKaons of only one sort are found
7031 {
7032 if(!n1) // No K0's only K+'s
7033 {
7034 if(n2==1 && mK+nucM>qM+.0001) // Mass limit: switch to K0
7035 {
7036 k1M=mK0;
7037 n1=1;
7038 qPDG=90000000+qP*1000+qN-1; // PDG of the Residual Nucleus
7040 }
7041 else
7042 {
7043 k1M=mK;
7044 k1PDG=321; // Only K+'s (default K0's)
7045 n1=n2; // only n1 is used
7046 }
7047 }
7048 else // No K+'s only K0's
7049 {
7050 if(n1==1 && mK0+nucM>qM+.0001) // Mass limit: switch to K+
7051 {
7052 k1M=mK;
7053 k1PDG=321; // K+ instead of K0
7054 qPDG=90000000+(qP-1)*1000+qN; // PDG of the Residual Nucleus
7056 }
7057 else k1M=mK0; // Only anti-K0's (default k1PDG)
7058 }
7059#ifdef debug
7063 if(k1PDG==321) // Calculate alternative to K+
7064 {
7065 naPDG=90000000+qP*1000+qN-1; // PDG of the Alternative Residual Nucleus
7067 kaM=mK0; // Prot Mass of the Alternative kaon (K0)
7068 }
7071#endif
7072 G4double n1M=n1*k1M;
7073 G4LorentzVector f4Mom(0.,0.,0.,n1M);
7074 G4LorentzVector s4Mom(0.,0.,0.,nucM);
7075 G4double sum=nucM+n1M;
7076 if(sum>qM+eps && n1==1) // Try to use another K if this is the only kaon
7077 {
7082 if(k1PDG==321) // Calculate alternative to the K+ meson
7083 {
7084 naPDG=90000000+qP*1000+qN-1; // PDG of the Alternative Residual Nucleus
7086 akPDG=311; // PDG Code of the Alternative kaon (K0)
7087 kaM=mK0; // Mass of the Alternative kaon (K0)
7088 }
7089 G4double asum=naM+kaM;
7090 if(asum<sum) // Make a KSwap correction
7091 {
7092 nucM=naM;
7093 n1M=kaM;
7094 k1M=kaM;
7095 k1PDG=akPDG;
7096 qPDG=naPDG;
7097 f4Mom=G4LorentzVector(0.,0.,0.,n1M);
7098 s4Mom=G4LorentzVector(0.,0.,0.,nucM);
7099 }
7100 }
7101 if(fabs(qM-sum)<eps)
7102 {
7103 f4Mom=q4M*(n1M/sum);
7104 s4Mom=q4M*(nucM/sum);
7105 }
7107 {
7108#ifdef debug
7111#endif
7112 evaHV->push_back(qH); // @@ Can cause problems with particle conversion in G4
7113 return;
7114 }
7115#ifdef debug
7117#endif
7118 delete qH;
7119 //
7120 f4Mom/=n1;
7122 {
7124 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7125 }
7127 //evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7129#ifdef debug
7131#endif
7132 }
7133 else
7134 {
7135 G4double n1M=n1*k1M;
7136 G4double n2M=n2*k2M;
7137 G4LorentzVector f4Mom(0.,0.,0.,n1M);
7138 G4LorentzVector s4Mom(0.,0.,0.,n2M);
7139 G4LorentzVector t4Mom(0.,0.,0.,nucM);
7140 G4double sum=nucM+n1M+n2M;
7141 if(fabs(qM-sum)<eps)
7142 {
7143 f4Mom=q4M*(n1M/sum);
7144 s4Mom=q4M*(n2M/sum);
7145 t4Mom=q4M*(nucM/sum);
7146 }
7148 {
7151 evaHV->push_back(qH); // @@ Can cause problems with particle conversion in G4
7152 return;
7153 }
7154#ifdef debug
7156 <<G4endl;
7157#endif
7158 delete qH;
7159 //
7160 f4Mom/=n1;
7162 {
7164 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7165 }
7166 s4Mom/=n2;
7168 {
7170 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7171 }
7173 //evaHV->push_back(H3); // Fill "H3" (delete equivalent)
7175 }
7176#ifdef qdebug
7177 if (qH)
7178 {
7180 << qH->GetPDGCode() << G4endl;
7181 delete qH;
7182 }
7183#endif
7184#ifdef debug
7186#endif
7187} // End of DecayAntiStrange
void EvaporateNucleus(G4QHadron *hA, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:4171
Referenced by DecayAntiDibaryon(), DecayDibaryon(), and EvaporateNucleus().
◆ DecayDibaryon()
| void G4QNucleus::DecayDibaryon | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 6286 of file G4QNucleus.cc.
6287{
6321#ifdef debug
6323#endif
6324 // Select a chanel of the dibaryon decay (including Delta+Delta-> 4 particle decay
6326 G4int sPDG = 2212;
6328 G4double fMass= mProt;
6329 G4double sMass= mProt;
6330 G4double tMass= mPi;
6331 if (qPDG==90003998 && rM>=dmPiP) // "diDelta++" case
6332 {
6333 sPDG = 211;
6334 sMass= mPi;
6335 four = true;
6336 }
6337 else if(qPDG==89998004 && rM>=dmPiN) // "diDelta--" case
6338 {
6339 sPDG = -211;
6340 fPDG = 2112;
6341 sMass= mPi;
6342 fMass= mNeut;
6343 four = true;
6344 }
6345 else if(qPDG==90000002 && rM>=dNeut) // "dineutron" case
6346 {
6347 fPDG = 2112;
6348 sPDG = 2112;
6349 fMass= mNeut;
6350 sMass= mNeut;
6351 }
6352 else if(qPDG==90001001 && rM>=mDeut) // "exited deutron" case
6353 {
6354 if(fabs(qM-mDeut)<eps)
6355 {
6356 evaHV->push_back(qH); // Fill as it is (delete equivalent)
6357 return;
6358 }
6359 else if(mProt+mNeut<rM)
6360 {
6361 fPDG = 2112;
6362 fMass= mNeut;
6363 }
6364 else
6365 {
6366 fPDG = 22;
6367 sPDG = 90001001;
6368 fMass= 0.;
6369 sMass= mDeut;
6370 }
6371 }
6372 else if(qPDG==91000001 && rM>=dLaNe) // "Lambda-neutron" case
6373 {
6374 fPDG = 2112;
6375 sPDG = 3122;
6376 fMass= mNeut;
6377 sMass= mLamb;
6378 }
6379 else if(qPDG==91001000 && rM>=dLaPr) // "Lambda-proton" case
6380 {
6381 sPDG = 3122;
6382 sMass= mLamb;
6383 }
6384 else if(qPDG==89999003 && rM>=nnPi) // "neutron/Delta-" case
6385 {
6386 fPDG = 2112;
6387 sPDG = 2112;
6388 tPDG = -211;
6389 fMass= mNeut;
6390 sMass= mNeut;
6391 }
6392 else if(qPDG==90002999 && rM>=ppPi) // "proton/Delta++" case
6393 {
6394 tPDG = 211;
6395 }
6396 else if(qPDG==90999002 && rM>=lnPi) // "lambda/Delta-" case
6397 {
6398 fPDG = 2112;
6399 sPDG = 3122;
6400 tPDG = -211;
6401 fMass= mNeut;
6402 sMass= mLamb;
6403 }
6404 else if(qPDG==91001999 && rM>=lpPi) // "lambda/Delta+" case
6405 {
6406 sPDG = 3122;
6407 tPDG = 211;
6408 sMass= mLamb;
6409 }
6410 else if(qPDG==90999002 && rM>=dSiNe) // "Sigma-/neutron" case
6411 {
6412 fPDG = 2112;
6413 sPDG = 3112;
6414 fMass= mNeut;
6415 sMass= mSigM;
6416 }
6417 else if(qPDG==91001999 && rM>=dSiPr) // "Sigma+/proton" case
6418 {
6419 sPDG = 3222;
6420 sMass= mSigP;
6421 }
6422 else if(qPDG==92000000 && rM>=dLamb) // "diLambda" case
6423 {
6424 fPDG = 3122;
6425 sPDG = 3122;
6426 fMass= mLamb;
6427 sMass= mLamb;
6428 }
6429 else if(qPDG==91999001 && rM>=dKsNe) // "Ksi-/neutron" case
6430 {
6431 fPDG = 2112;
6432 sPDG = 3312;
6433 fMass= mNeut;
6434 sMass= mKsiM;
6435 }
6436 else if(qPDG==92000999 && rM>=dKsPr) // "Ksi0/proton" case
6437 {
6438 sPDG = 3322;
6439 sMass= mKsiZ;
6440 }
6441 else if(qPDG!=90002000|| rM<dProt) // Other possibilities (if not a default)
6442 {
6443 G4int qS = qH->GetStrangeness();
6444 G4int qB = qH->GetBaryonNumber();
6445 if(qB>0&&qS<0) // Antistrange diBarion
6446 {
6447 DecayAntiStrange(qH,evaHV);
6448 return;
6449 }
6450 else
6451 {
6452 delete qH;
6455 // @@ Nothing to do. Just 2 GeV disappears... Very rare! Just to avoid the exception.
6456 //throw G4QException("G4QNucleus::DecayDibar: Unknown PDG code or small Mass of DB");
6457 }
6458 }
6459 G4LorentzVector f4Mom(0.,0.,0.,fMass);
6460 G4LorentzVector s4Mom(0.,0.,0.,sMass);
6461 G4LorentzVector t4Mom(0.,0.,0.,tMass);
6462 if(!tPDG&&!four)
6463 {
6464 G4double sum=fMass+sMass;
6465 if(fabs(qM-sum)<eps)
6466 {
6467 f4Mom=q4M*(fMass/sum);
6468 s4Mom=q4M*(sMass/sum);
6469 }
6471 {
6474 //G4cerr<<"***G4QNucl::DecayDiBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6475 //throw G4QException("***G4QNucleus::DecayDibaryon: DiBaryon DecayIn2 error");
6476 evaHV->push_back(qH);
6477 return;
6478 }
6479#ifdef debug
6482#endif
6483 delete qH;
6485 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6487 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6488 }
6489 else if(four)
6490 {
6491 q4M=q4M/2.; // Divided in 2 !!!
6492 qM/=2.; // Divide the mass in 2 !
6493 G4double sum=fMass+sMass;
6494 if(fabs(qM-sum)<eps)
6495 {
6496 f4Mom=q4M*(fMass/sum);
6497 s4Mom=q4M*(sMass/sum);
6498 }
6500 {
6503 //G4cerr<<"***G4QN::DecayDibaryon: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6504 //throw G4QException("***G4QNucleus::DecDibaryon: Dibaryon DecayIn2 error");
6505 evaHV->push_back(qH);
6506 return;
6507 }
6508#ifdef debug
6511#endif
6513 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6515 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6516 // Now the second pair mus be decayed
6517 if(fabs(qM-sum)<eps)
6518 {
6519 f4Mom=q4M*(fMass/sum);
6520 s4Mom=q4M*(sMass/sum);
6521 }
6523 {
6524 // Should not be here as sum was already compared with qM above for the first delta
6525 delete qH;
6526 // G4cerr<<"***G4QNucl::DecDibar:fPDG="<<fPDG<<"(fM="<<fMass<<") + sPDG="<<sPDG<<"(sM="
6527 // <<sMass<<")="<<sum<<" >? (DD2,Can't be here) TotM="<<q4M.m()<<q4M<<G4endl;
6528 // throw G4QException("G4QNucleus::DecayDibaryon: General DecayIn2 error");
6529 G4ExceptionDescription ed;
6530 ed << "General DecayIn2 error: fPDG=" << fPDG << "(fM=" << fMass
6531 << ") + sPDG=" << sPDG <<"(sM=" << sMass << ")=" << sum
6534 FatalException, ed);
6535 }
6536#ifdef debug
6539#endif
6541 evaHV->push_back(H3); // Fill "H1" (delete equivalent)
6543 evaHV->push_back(H4); // Fill "H2" (delete equivalent)
6544 delete qH;
6545 }
6546 else
6547 {
6548 G4double sum=fMass+sMass+tMass;
6549 if(fabs(qM-sum)<eps)
6550 {
6551 f4Mom=q4M*(fMass/sum);
6552 s4Mom=q4M*(sMass/sum);
6553 t4Mom=q4M*(tMass/sum);
6554 }
6556 {
6559 //G4cerr<<"***G4QNuc::DecayDibaryon:qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
6560 //throw G4QException("G4QNucleus::DecayDibaryon: diBar DecayIn3 error");
6561 evaHV->push_back(qH);
6562 return;
6563 }
6564#ifdef debug
6567#endif
6568 //qH->SetNFragments(2); // Fill a#of fragments to decaying Dibaryon
6569 //evaHV->push_back(qH); // Fill hadron with nf=2 (delete equivalent)
6570 // Instead
6571 delete qH;
6572 //
6574 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
6576 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
6578 evaHV->push_back(H3); // Fill "H3" (delete equivalent)
6579 }
6580#ifdef qdebug
6581 if (qH)
6582 {
6584 << qH->GetPDGCode() << G4endl;
6585 delete qH;
6586 }
6587#endif
6588} // End of DecayDibaryon
Referenced by EvaporateNucleus().
◆ DecayIsonucleus()
| void G4QNucleus::DecayIsonucleus | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 5925 of file G4QNucleus.cc.
5926{
5931 //static const G4double mSigZ= G4QPDGCode(3212).GetMass();
5941#ifdef debug
5943#endif
5944 if(qS<0||qS>qBN) // *** Should not be here ***
5945 {
5947 evaHV->push_back(qH); // fill as it is (delete equivalent)
5948 return;
5949 }
5952 G4int sPDG = 211;
5954 G4double fMass= mProt;
5955 G4double sMass= mPi;
5957 // =---------= Negative state =-----------=
5958 if(qC<0) // =-------= Only Pi- can help
5959 {
5960 if(qS&&qBN==qS) // --- n*Lamb + k*(Pi-) State ---
5961 {
5962 sPDG = -211;
5963 if(-qC==qS && qS==1) // Only one Sigma- like (qBN=1)
5964 {
5965 if(fabs(qM-mSigM)<eps)
5966 {
5967 evaHV->push_back(qH); // Fill Sigma- as it is
5968 return;
5969 }
5970 else if(qM>mLamb+mPi) //(2) Sigma- => Lambda + Pi- decay
5971 {
5972 fPDG = 3122;
5973 fMass= mLamb;
5974 }
5975 else if(qM>mSigM) //(2) Sigma+ => Sigma+ + gamma decay
5976 {
5977 fPDG = 3112;
5978 fMass= mSigM;
5979 sPDG = 22;
5980 sMass= 0.;
5981 }
5982 else //(2) Sigma- => Neutron + Pi- decay
5983 {
5984 fPDG = 2112;
5985 fMass= mNeut;
5986 }
5987 qPN = 1; // #of (Pi+ or gamma)'s = 1
5988 }
5989 else if(-qC==qS) //(2) a few Sigma- like
5990 {
5991 qPN = 1; // One separated Sigma-
5992 fPDG = 3112;
5993 sPDG = 3112;
5994 sMass= mSigM;
5995 qBN--;
5996 fMass= mSigM;
5997 }
5998 else if(-qC>qS) //(2) n*(Sigma-)+m*(Pi-)
5999 {
6000 qPN = -qC-qS; // #of Pi-'s
6001 fPDG = 3112;
6002 fMass= mSigM;
6003 }
6004 else //(2) n*(Sigma-)+m*Lambda (-qC<qS)
6005 {
6006 qBN += qC; // #of Lambda's
6007 fPDG = 3122;
6008 fMass= mLamb;
6009 qPN = -qC; // #of Sigma+'s
6010 sPDG = 3112;
6011 sMass= mSigM;
6012 }
6013 qS = 0; // Only decays in two are above
6014 }
6015 else if(qS) // ->n*Lamb+m*Neut+k*(Pi-) State (qS<qBN)
6016 {
6017 qBN -= qS; // #of neutrons
6018 fPDG = 2112;
6019 fMass= mNeut;
6021 if(qS==nPin) //(2) m*Neut+n*Sigma-
6022 {
6023 qPN = qS; // #of Sigma-
6024 sPDG = 3112;
6025 sMass= mSigM;
6026 qS = 0;
6027 }
6028 else if(qS>nPin) //(3) m*P+n*(Sigma+)+k*Lambda
6029 {
6030 qS-=nPin; // #of Lambdas
6031 qPN = nPin; // #of Sigma+
6032 sPDG = 3112;
6033 sMass= mSigM;
6034 }
6035 else //(3) m*N+n*(Sigma-)+k*(Pi-) (qS<nPin)
6036 {
6037 qPN = nPin-qS; // #of Pi-
6038 sPDG = -211;
6039 tPDG = 3112;
6040 tMass= mSigM;
6041 }
6042 }
6043 else //(2) n*N+m*(Pi-) (qS=0)
6044 {
6045 sPDG = -211;
6046 qPN = -qC;
6047 fPDG = 2112;
6048 fMass= mNeut;
6049 }
6050 }
6051 else if(!qC) // *** Should not be here ***
6052 {
6053 if(qS && qS<qBN) //(2) n*Lamb+m*N ***Should not be here***
6054 {
6055 qPN = qS;
6056 fPDG = 2112; // mN+nL case
6057 sPDG = 3122;
6058 sMass= mLamb;
6059 qBN -= qS;
6060 fMass= mNeut;
6061 qS = 0;
6062 }
6063 else if(qS>1 && qBN==qS) //(2) m*Lamb(m>1) ***Should not be here***
6064 {
6065 qPN = 1;
6066 fPDG = 3122;
6067 sPDG = 3122;
6068 sMass= mLamb;
6069 qBN--;
6070 fMass= mLamb;
6071 }
6072 else if(!qS && qBN>1) //(2) n*Neut(n>1) ***Should not be here***
6073 {
6074 qPN = 1;
6075 fPDG = 2112;
6076 sPDG = 2112;
6077 sMass= mNeut;
6078 qBN--;
6079 fMass= mNeut;
6080 }
6082 }
6083 else if(qC>0) // n*Lamb+(m*P)+(k*Pi+)
6084 {
6085 if(qS && qS+qC==qBN) //(2) n*Lamb+m*P ***Should not be here***
6086 {
6087 qPN = qS;
6088 qS = 0;
6089 fPDG = 2212;
6090 sPDG = 3122;
6091 sMass= mLamb;
6092 qBN = qC;
6093 fMass= mProt;
6094 }
6095 else if(qS && qC<qBN-qS) //(3)n*L+m*P+k*N ***Should not be here***
6096 {
6097 qPN = qC; // #of protons
6098 fPDG = 2112; // mP+nL case
6099 sPDG = 2212;
6100 sMass= mProt;
6101 qBN -= qS+qC; // #of neutrons
6102 fMass= mNeut;
6103 }
6104 else if(qS && qBN==qS) // ---> n*L+m*Pi+ State
6105 {
6106 if(qC==qS && qS==1) // Only one Sigma+ like State
6107 {
6108 if(fabs(qM-mSigP)<eps) // Fill Sigma+ as it is
6109 {
6110 evaHV->push_back(qH);
6111 return;
6112 }
6113 else if(qM>mLamb+mPi) //(2) Sigma+ => Lambda + Pi+ decay
6114 {
6115 fPDG = 3122;
6116 fMass= mLamb;
6117 }
6118 else if(qM>mNeut+mPi) //(2) Sigma+ => Neutron + Pi+ decay
6119 {
6120 fPDG = 2112;
6121 fMass= mNeut;
6122 }
6123 else if(qM>mSigP) //(2) Sigma+ => Sigma+ + gamma decay
6124 {
6125 fPDG = 3222;
6126 fMass= mSigP;
6127 sPDG = 22;
6128 sMass= 0.;
6129 }
6130 else //(2) Sigma+ => Proton + gamma decay
6131 {
6132 fPDG = 2212;
6133 fMass= mProt;
6134 sPDG = 22;
6135 sMass= 0.;
6136 }
6137 qPN = 1; // #of (Pi+ or gamma)'s = 1
6138 }
6139 else if(qC==qS) //(2) a few Sigma+ like hyperons
6140 {
6141 qPN = 1;
6142 fPDG = 3222;
6143 sPDG = 3222;
6144 sMass= mSigP;
6145 qBN--;
6146 fMass= mSigP;
6147 }
6148 else if(qC>qS) //(2) n*(Sigma+)+m*(Pi+)
6149 {
6150 qPN = qC-qS; // #of Pi+'s
6151 fPDG = 3222;
6152 qBN = qS; // #of Sigma+'s
6153 fMass= mSigP;
6154 }
6155 else //(2) n*(Sigma+)+m*Lambda
6156 {
6157 qBN -= qC; // #of Lambda's
6158 fPDG = 3122;
6159 fMass= mLamb;
6160 qPN = qC; // #of Sigma+'s
6161 sPDG = 3222;
6162 sMass= mSigP;
6163 }
6164 qS = 0; // All above are decays in 2
6165 }
6166 else if(qS && qC>qBN-qS) // n*Lamb+m*P+k*Pi+
6167 {
6168 qBN -= qS; // #of protons
6170 if(qS==nPip) //(2) m*P+n*Sigma+
6171 {
6172 qPN = qS; // #of Sigma+
6173 sPDG = 3222;
6174 sMass= mSigP;
6175 qS = 0;
6176 }
6177 else if(qS>nPip) //(3) m*P+n*(Sigma+)+k*Lambda
6178 {
6179 qS -= nPip; // #of Lambdas
6180 qPN = nPip; // #of Sigma+
6181 sPDG = 3222;
6182 sMass= mSigP;
6183 }
6184 else //(3) m*P+n*(Sigma+)+k*(Pi+)
6185 {
6186 qPN = nPip-qS; // #of Pi+
6187 tPDG = 3222;
6188 tMass= mSigP;
6189 }
6190 }
6191 if(qC<qBN) //(2) n*P+m*N ***Should not be here***
6192 {
6193 fPDG = 2112;
6194 fMass= mNeut;
6195 qPN = qC;
6196 sPDG = 2212;
6197 sMass= mProt;
6198 }
6199 else if(qBN==qC && qC>1) //(2) m*Prot(m>1) ***Should not be here***
6200 {
6201 qPN = 1;
6202 fPDG = 2212;
6203 sPDG = 2212;
6204 sMass= mProt;
6205 qBN--;
6206 fMass= mProt;
6207 }
6209 // !qS && qC>qBN //(2) Default condition n*P+m*(Pi+)
6210 }
6211 G4double tfM=qBN*fMass;
6212 G4double tsM=qPN*sMass;
6213 G4double ttM=0.;
6214 if(qS) ttM=qS*tMass;
6215 G4LorentzVector f4Mom(0.,0.,0.,tfM);
6216 G4LorentzVector s4Mom(0.,0.,0.,tsM);
6217 G4LorentzVector t4Mom(0.,0.,0.,ttM);
6218 G4double sum=tfM+tsM+ttM;
6219 if(fabs(qM-sum)<eps)
6220 {
6221 f4Mom=q4M*(tfM/sum);
6222 s4Mom=q4M*(tsM/sum);
6223 if(qS) t4Mom=q4M*(ttM/sum);
6224 }
6226 {
6227#ifdef debug
6230#endif
6231 evaHV->push_back(qH); // fill as it is (delete equivalent)
6232 return;
6233 }
6235 {
6236#ifdef debug
6239#endif
6240 evaHV->push_back(qH); // fill as it is (delete equivalent)
6241 return;
6242 }
6243#ifdef debug
6246#endif
6247 delete qH;
6248 if(qBN)
6249 {
6250 f4Mom/=qBN;
6252 {
6254 evaHV->push_back(Hi); // Fill "Hi" (delete equivalent)
6255 }
6256 }
6257 if(qPN)
6258 {
6259 s4Mom/=qPN;
6261 {
6263 evaHV->push_back(Hj); // Fill "Hj" (delete equivalent)
6264 }
6265 }
6266 if(qS)
6267 {
6268 t4Mom/=qS;
6270 {
6272 evaHV->push_back(Hk); // Fill "Hk" (delete equivalent)
6273 }
6274 }
6275#ifdef qdebug
6276 if (qH)
6277 {
6279 << qH->GetPDGCode() << G4endl;
6280 delete qH;
6281 }
6282#endif
6283} // End of DecayIsonucleus
Referenced by EvaporateNucleus().
◆ DecayMultyBaryon()
| void G4QNucleus::DecayMultyBaryon | ( | G4QHadron * | dB, |
| G4QHadronVector * | oHV | ||
| ) |
Definition at line 7190 of file G4QNucleus.cc.
7191{
7200#ifdef debug
7202#endif
7208 G4double fMass= mLamb;
7209 if (totN==totBN) // "A-neutron" case
7210 {
7211 fPDG = 2112;
7212 fMass= mNeut;
7213 }
7214 else if(totC==totBN) // "A-protons" case
7215 {
7216 fPDG = 2212;
7217 fMass= mProt;
7218 }
7219 else if(totS!=totBN) // "Bad call" case
7220 {
7221 delete qH;
7222 // G4cerr<<"***G4QNuc::DecayMultyBaryon: PDG="<<qPDG<<G4endl;
7223 // throw G4QException("***G4QNuc::DecayMultyBaryon: Can not decay this PDG Code");
7224 G4ExceptionDescription ed;
7227 FatalException, ed);
7228 }
7229#ifdef debug
7230 else
7231 {
7232 delete qH;
7233 // G4cerr<<"**G4QNucleus::DecayMultyBaryon: PDG="<<qPDG<<G4endl;
7234 // throw G4QException("***G4QNuc::DecayMultyBaryon: Unknown PDG code of the MultiBaryon");
7235 G4ExceptionDescription ed;
7238 FatalException, ed);
7239 }
7240#endif
7241 if(totBN==1) evaHV->push_back(qH);
7242 else if(totBN==2)
7243 {
7244 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7245 G4LorentzVector s4Mom(0.,0.,0.,fMass);
7246 G4double sum=fMass+fMass;
7247 if(fabs(qM-sum)<eps)
7248 {
7249 f4Mom=q4M/2.;
7250 s4Mom=f4Mom;
7251 }
7253 {
7256 //G4cerr<<"***G4QNuc::DecayMultyBaryon:qM="<<qM<<"<sum="<<sum<<",d="<<sum-qM<<G4endl;
7257 //throw G4QException("G4QNuc::DecayMultyBaryon:diBaryon DecayIn2 didn't succeed");
7258 evaHV->push_back(qH);
7259 return;
7260 }
7261#ifdef debug
7263#endif
7264 delete qH;
7266 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7268 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7269 }
7270 else if(totBN==3)
7271 {
7272 G4LorentzVector f4Mom(0.,0.,0.,fMass);
7273 G4LorentzVector s4Mom(0.,0.,0.,fMass);
7274 G4LorentzVector t4Mom(0.,0.,0.,fMass);
7275 G4double sum=fMass+fMass+fMass;
7276 if(fabs(qM-sum)<eps)
7277 {
7278 f4Mom=q4M/3.;
7279 s4Mom=f4Mom;
7280 t4Mom=f4Mom;
7281 }
7283 {
7286 //G4cerr<<"***G4QN::DecayMultyBar: qM="<<qM<<" < sum="<<sum<<",d="<<sum-qM<<G4endl;
7287 //throw G4QException("G4QNucleus::DecayMultyBar:ThreeBaryonDecayIn3 didn't succeed");
7288 evaHV->push_back(qH);
7289 return;
7290 }
7291#ifdef debug
7293 <<t4Mom<<G4endl;
7294#endif
7295 delete qH;
7297 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7299 evaHV->push_back(H2); // Fill "H2" (delete equivalent)
7301 evaHV->push_back(H3); // Fill "H3" (delete equivalent)
7302 }
7303 else
7304 {
7305 // @@It must be checked, that they are not under the mass shell
7306 // !! OK !! Checked by the warning print that they are mostly in the Ground State !!
7308#ifdef debug
7309 // Warning for the future development
7312#endif
7313 delete qH;
7315 {
7317 evaHV->push_back(H1); // Fill "H1" (delete equivalent)
7318 }
7319 }
7320#ifdef qdebug
7321 if (qH)
7322 {
7324 << qH->GetPDGCode() << G4endl;
7325 delete qH;
7326 }
7327#endif
7328} // End of DecayMultyBaryon
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
Definition: G4QHadron.cc:544
G4bool DecayIn3(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &t4Mom)
Definition: G4QHadron.cc:783
Referenced by DecayAlphaBar(), and EvaporateNucleus().
◆ DeleteNucleons()
| void G4QNucleus::DeleteNucleons | ( | ) |
Definition at line 697 of file G4QNucleus.cc.
698{
699 G4QHadronVector::iterator u; // iterator for the nucleons
700 for(u=theNucleons.begin(); u!=theNucleons.end(); u++) delete *u;
702}
Referenced by G4QIonIonCollision::G4QIonIonCollision().
◆ DoLorentzBoost() [1/2]
|
inline |
Definition at line 155 of file G4QNucleus.hh.
Referenced by G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), and Init3D().
◆ DoLorentzBoost() [2/2]
|
inline |
Definition at line 165 of file G4QNucleus.hh.
166 {
169 }
◆ DoLorentzContraction() [1/2]
|
inline |
Definition at line 170 of file G4QNucleus.hh.
Hep3Vector vect() const
double e() const
void DoLorentzContraction(const G4LorentzVector &B)
Definition: G4QNucleus.hh:170
Referenced by DoLorentzContraction().
◆ DoLorentzContraction() [2/2]
| void G4QNucleus::DoLorentzContraction | ( | const G4ThreeVector & | theBeta | ) |
Definition at line 3960 of file G4QNucleus.cc.
3961{
3964 G4int theA=theNucleons.size();
3966 {
3967 G4ThreeVector pos=theNucleons[i]->GetPosition();
3968 pos -= factor*(theBeta*pos)*theBeta;
3969 theNucleons[i]->SetPosition(pos);
3970 }
3971} // End of DoLorentzContraction(G4ThreeVector)
◆ DoLorentzRotation()
|
inline |
Definition at line 160 of file G4QNucleus.hh.
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ DoTranslation()
| void G4QNucleus::DoTranslation | ( | const G4ThreeVector & | theShift | ) |
Definition at line 3974 of file G4QNucleus.cc.
◆ EvaporateBaryon()
Definition at line 951 of file G4QNucleus.cc.
952{
953 //static const G4double uWell=2.7; // EffectiveDepth of potential well B
954 //static const G4double uWell=7.; // EffectiveDepth of potential well B
956 //static const G4double uWell=0.0; // EffectiveDepth of potential well B
957 //////////static const G4double gunA=80.; // Switch A-parameter for BaryonGun
958 //static const G4double gunB=exp(1)/gunA;
959 ///////////////////static const G4double gunB=exp(2)/4/gunA/gunA;
960 //////////////static const G4double gunP2=200000.; // Switch P2-parameter for BaryonGun
961 //////////////static const G4double maSht=1.2; // shift for maximal x approximation
962 ///////////static const G4double coSht=.19; // multiple for maximal x approximation
963 //////////////static const G4double third=1./3.;// power for maximal x approximation
1003 //static const G4double mNN = mNeut+mNeut; // 2 neutrons mass
1004 //static const G4double mPP = mProt+mProt; // 2 protons mass
1005 //static const G4double mNL = mNeut+mLamb; // neutron and Lambda mass
1006 //static const G4double mPL = mProt+mLamb; // proton and Lambda mass
1007 //static const G4double mLL = mLamb+mLamb; // 2 Lambdas mass
1009 G4double uW=uWell;
1012 //if(a>4.5) evalph=2.7/sqrt(a-4.); // Probability for alpha to evaporate
1013 //G4double evalph=clustProb*clustProb*clustProb;
1014#ifdef debug
1016#endif
1017 G4double a1= a-1;
1018 //////////G4double z = Z;
1019 //////////G4double zn= Z+N;
1026 //G4double PPPBarr= CoulombBarrier(1,1,2,2); // CoulombBarrier for proton (after 2 pr)
1027 //G4double AAABarr= CoulombBarrier(2,4,4,8); // CoulombBarrier for alpha (after 2alph)
1028 //////G4double APABarr= CoulombBarrier(2,4,3,5);// CoulombBarrier for alpha (after p+al)
1029 //G4double PPABarr= CoulombBarrier(1,1,3,5); // CoulombBarrier for proton (after p+al)
1032 //G4double SPPPBarr=SPPBarr+PPPBarr; // SummedCoulombBarrier for 3 protons
1033 //G4double SAAABarr=SAABarr+AAABarr; // SummedCoulombBarrier for 3 alphas
1036 if(DAPBarr>SAPBarr)SAPBarr=DAPBarr; // Get max to make possible BothSequences
1037 ///////G4double SAPABarr=APABarr+SAPBarr; // Summed Coulomb Barrier for alph+p+alph
1038 //G4double SPPABarr=PPABarr+SAPBarr; // Summed Coulomb Barrier for p+p+alpha
1039 G4LorentzVector h1mom;
1040 G4LorentzVector h2mom;
1041 G4LorentzVector h3mom;
1043#ifdef debug
1045#endif
1046 if(a==-2)
1047 {
1048 if(Z==1 || N==1)
1049 {
1050 G4int nucPDG = -2112;
1051 G4int piPDG = 211;
1052 G4double nucM = mNeut;
1053 G4QPDGCode del = aDmQPDG;
1055 if(N>0)
1056 {
1057 nucPDG = -2212;
1058 piPDG = -211;
1059 nucM = mProt;
1060 del = aDppQPDG;
1061 nuc = apQPDG;
1062 }
1063 if(totMass > mPi+nucM+nucM)
1064 {
1065 G4LorentzVector n14M(0.,0.,0.,nucM);
1066 G4LorentzVector n24M(0.,0.,0.,nucM);
1067 G4LorentzVector pi4M(0.,0.,0.,mPi);
1069 {
1072 //throw G4QException("G4QNucl::EvapBary:ISO-dibaryon DecayIn3 did not succeed");
1073 return false;
1074 }
1075 n14M+=pi4M;
1076 h1->SetQPDG(del);
1077 h2->SetQPDG(nuc);
1078 h1->Set4Momentum(n14M);
1079 h2->Set4Momentum(n24M);
1080 return true;
1081 }
1082 else
1083 {
1086 //throw G4QException("***G4QNucl::EvaporateBaryon: ISO-dibaryon under Mass Shell");
1087 return false;
1088 }
1089 }
1090 else if(Z==2 || N==2)
1091 {
1092 G4int nucPDG = -2112;
1093 G4int piPDG = 211;
1094 G4double nucM = mNeut;
1095 G4QPDGCode del = aDmQPDG;
1096 if(N==2)
1097 {
1098 nucPDG = -2212;
1099 piPDG = -211;
1100 nucM = mProt;
1101 del = aDppQPDG;
1102 }
1103 if(totMass > mPi+mPi+nucM+nucM)
1104 {
1105 G4LorentzVector n14M(0.,0.,0.,nucM);
1106 G4LorentzVector n24M(0.,0.,0.,nucM);
1107 G4LorentzVector pi4M(0.,0.,0.,mPi+mPi);
1109 {
1112 //throw G4QException("G4QNucl::EvapBary:ISO-dibaryon DecayIn3 did not succeed");
1113 return false;
1114 }
1115 G4LorentzVector hpi4M=pi4M/2.;
1116 n14M+=hpi4M;
1117 n24M+=hpi4M;
1118 h1->SetQPDG(del);
1119 h2->SetQPDG(del);
1120 h1->Set4Momentum(n14M);
1121 h2->Set4Momentum(n24M);
1122 return true;
1123 }
1124 else
1125 {
1128 //throw G4QException("***G4QNucl::EvaporateBaryon: ISO-dibaryon under Mass Shell");
1129 return false;
1130 }
1131 }
1132 else if(Z==-2)
1133 {
1134 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1135 h2mom=h1mom;
1136 h1->SetQPDG(apQPDG);
1137 h2->SetQPDG(apQPDG);
1139 }
1140 else if(N==-2)
1141 {
1142 h1mom=G4LorentzVector(0.,0.,0.,mNeut);
1143 h2mom=h1mom;
1144 h1->SetQPDG(anQPDG);
1145 h2->SetQPDG(anQPDG);
1147 }
1148 else if(N==-1 && Z==-1)
1149 {
1150 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1151 h2mom=G4LorentzVector(0.,0.,0.,mNeut);
1152 h1->SetQPDG(apQPDG);
1153 h2->SetQPDG(anQPDG);
1155 }
1156 else if(Z==-1 && S==-1)
1157 {
1158 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1159 h2mom=G4LorentzVector(0.,0.,0.,mLamb);
1160 h1->SetQPDG(apQPDG);
1161 h2->SetQPDG(alQPDG);
1163 }
1164 else
1165 {
1166 h1mom=G4LorentzVector(0.,0.,0.,mNeut);
1167 h2mom=G4LorentzVector(0.,0.,0.,mLamb);
1168 h1->SetQPDG(anQPDG);
1169 h2->SetQPDG(alQPDG);
1171 }
1172 h1->Set4Momentum(h1mom);
1173 h2->Set4Momentum(h2mom);
1174 return true;
1175 }
1176 else if(a==2)
1177 {
1178 if(Z<0||N<0)
1179 {
1180 G4int nucPDG = 2112;
1181 G4double nucM = mNeut;
1182 G4int piPDG = -211;
1183 G4QPDGCode db = NNQPDG;
1184 G4QPDGCode pi_value = PIMQPDG;
1185 if(N<0)
1186 {
1187 nucPDG = 2212;
1188 nucM = mProt;
1189 piPDG = 211;
1190 db = PPQPDG;
1191 pi_value = PIPQPDG;
1192 }
1193 if(totMass>mPi+nucM+nucM)
1194 {
1195 G4LorentzVector n14M(0.,0.,0.,nucM);
1196 G4LorentzVector n24M(0.,0.,0.,nucM);
1197 G4LorentzVector pi4M(0.,0.,0.,mPi);
1199 {
1202 //throw G4QException("G4QNucl::EvapBary:ISO-dibaryon DecayIn3 did not succeed");
1203 return false;
1204 }
1205 n14M+=n24M;
1206 h1->SetQPDG(db);
1207 h2->SetQPDG(pi_value);
1208 h1->Set4Momentum(n14M);
1209 h2->Set4Momentum(pi4M);
1210 return true;
1211 }
1212 else
1213 {
1216 //throw G4QException("***G4QNucl::EvaporateBaryon: ISO-dibaryon under Mass Shell");
1217 return false;
1218 }
1219 }
1220 else if(Z==2)
1221 {
1222 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1223 h2mom=h1mom;
1224 h1->SetQPDG(pQPDG);
1225 h2->SetQPDG(pQPDG);
1227 }
1228 else if(N==2)
1229 {
1230 h1mom=G4LorentzVector(0.,0.,0.,mNeut);
1231 h2mom=h1mom;
1232 h1->SetQPDG(nQPDG);
1233 h2->SetQPDG(nQPDG);
1235 }
1236 else if(N==1&&Z==1)
1237 {
1238 if(totMass<=mNP)
1239 {
1240#ifdef debug
1242#endif
1243 h1mom=G4LorentzVector(0.,0.,0.,0.);
1244 h2mom=G4LorentzVector(0.,0.,0.,mDeut);
1245 h1->SetQPDG(gQPDG);
1246 h2->SetQPDG(NPQPDG);
1247 }
1248 else
1249 {
1250 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1251 h2mom=G4LorentzVector(0.,0.,0.,mNeut);
1252 h1->SetQPDG(pQPDG);
1253 h2->SetQPDG(nQPDG);
1254 }
1256 }
1257 else if(Z==1&&S==1)
1258 {
1259 h1mom=G4LorentzVector(0.,0.,0.,mProt);
1260 h2mom=G4LorentzVector(0.,0.,0.,mLamb);
1261 h1->SetQPDG(pQPDG);
1262 h2->SetQPDG(lQPDG);
1264 }
1265 else
1266 {
1267 h1mom=G4LorentzVector(0.,0.,0.,mNeut);
1268 h2mom=G4LorentzVector(0.,0.,0.,mLamb);
1269 h1->SetQPDG(nQPDG);
1270 h2->SetQPDG(lQPDG);
1272 }
1273 h1->Set4Momentum(h1mom);
1274 h2->Set4Momentum(h2mom);
1275 return true;
1276 }
1277 else if(a>2)
1278 {
1293 //G4bool nnnF = false; // Evaporation brunch is closed
1294 //G4bool nnpF = false;
1295 //G4bool nppF = false;
1296 //G4bool pppF = false;
1297 //G4bool nnlF = false;
1298 //G4bool nplF = false;
1299 //G4bool pplF = false;
1300 //G4bool nllF = false;
1301 //G4bool pllF = false;
1302 //G4bool lllF = false;
1303 //G4bool nnaF = false;
1304 //G4bool npaF = false;
1305 //G4bool ppaF = false;
1306 //G4bool nlaF = false;
1307 //G4bool plaF = false;
1308 //G4bool llaF = false;
1309 //G4bool paaF = false;
1310 //G4bool naaF = false;
1311 //G4bool laaF = false;
1312 //G4bool aaaF = false;
1325 /*
1326 // DHW 16 June 2011 : these variables set but not used. Comment out to fix
1327 // compiler warnings
1328 G4double GSReNNN= GSMass; // Prototype of Residual Nuclear Mass for n+n+n
1329 G4double GSReNNP= GSMass; // Prototype of Residual Nuclear Mass for n+n+p
1330 G4double GSReNPP= GSMass; // Prototype of Residual Nuclear Mass for n+p+p
1331 G4double GSRePPP= GSMass; // Prototype of Residual Nuclear Mass for p+p+p
1332 G4double GSReNNL= GSMass; // Prototype of Residual Nuclear Mass for n+n+l
1333 G4double GSReNPL= GSMass; // Prototype of Residual Nuclear Mass for n+p+l
1334 G4double GSRePPL= GSMass; // Prototype of Residual Nuclear Mass for p+p+l
1335 G4double GSReNLL= GSMass; // Prototype of Residual Nuclear Mass for n+l+l
1336 G4double GSRePLL= GSMass; // Prototype of Residual Nuclear Mass for p+l+l
1337 G4double GSReLLL= GSMass; // Prototype of Residual Nuclear Mass for l+l+l
1338 G4double GSReNNA= GSMass; // Prototype of Residual Nuclear Mass for n+n+a
1339 G4double GSReNPA= GSMass; // Prototype of Residual Nuclear Mass for n+p+a
1340 G4double GSRePPA= GSMass; // Prototype of Residual Nuclear Mass for p+p+a
1341 G4double GSReNLA= GSMass; // Prototype of Residual Nuclear Mass for n+l+a
1342 G4double GSRePLA= GSMass; // Prototype of Residual Nuclear Mass for p+l+a
1343 G4double GSReLLA= GSMass; // Prototype of Residual Nuclear Mass for l+l+a
1344 G4double GSRePAA= GSMass; // Prototype of Residual Nuclear Mass for p+a+a
1345 G4double GSReNAA= GSMass; // Prototype of Residual Nuclear Mass for n+a+a
1346 G4double GSReLAA= GSMass; // Prototype of Residual Nuclear Mass for l+a+a
1347 G4double GSReAAA= GSMass; // Prototype of Residual Nuclear Mass for a+a+a
1348 */
1368#ifdef debug
1370 <<G4endl;
1371#endif
1372 G4double tM2 = totMass*totMass;
1373 G4double qtM2 = 4*tM2;
1382 if(Z>0)
1383 {
1384 PQPDG=G4QPDGCode(90000000+1000*(1000*S+Z-1)+N);
1385 GSResNp=PQPDG.GetMass();
1386 G4double mpls=GSResNp+mProt;
1387 G4double mmin=GSResNp-mProt;
1388 pp2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1389 if(pp2m>=0.000001)
1390 {
1391 pFlag=true;
1392 pBnd=mProt-GSMass+GSResNp; // Binding energy for proton
1393 G4double eMax=sqrt(mP2+pp2m);
1394#ifdef debug
1397#endif
1398 pExcess=eMax-mProt+pBnd; // Max Kin Energy from bottom
1399 }
1400 else pExcess=pBnd;
1401 if(Z>1)
1402 {
1403 ppQPDG=G4QPDGCode(90000000+1000*(1000*S+Z-2)+N);
1404 GSResPP=ppQPDG.GetMass();
1405#ifdef debug
1409#endif
1410 if(GSResPP+mProt+mProt+SPPBarr<totMass) ppFlag=true;
1411 if(Z>2&&a>3)
1412 {
1413 /*
1414 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1415 GSRePPP=G4QPDGCode().GetNuclMass(Z-3,N,S);
1416 */
1417 //if(GSRePPP+mProt+mProt+mProt+SPPPBarr<totMass) pppF=true;
1418 if(N>1&&a>5)
1419 {
1420 paQPDG =G4QPDGCode(90000000+1000*(1000*S+Z-3)+N-2);
1421 GSResPA=paQPDG.GetMass();
1422#ifdef debug
1423 G4double s_value=GSResPA+mAlph+mProt+SAPBarr;
1425#endif
1426 if(GSResPA+mProt+SAPBarr+mAlph<totMass) paFlag=true;
1427 }
1428 }
1429 if(N>0&&a>3)
1430 {
1431 /*
1432 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1433 GSReNPP=G4QPDGCode().GetNuclMass(Z-2,N-1,S);
1434 */
1435 //if(GSReNPP+mProt+mProt+SPPBarr+mNeut<totMass) nppF=true;
1436 }
1437 if(S>0&&a>3)
1438 {
1439 /*
1440 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1441 GSRePPL=G4QPDGCode().GetNuclMass(Z-2,N,S-1);
1442 */
1443 //if(GSRePPL+mProt+mProt+SPPBarr+mLamb<totMass) pplF=true;
1444 }
1445 if(N>1&&a>4)
1446 {
1447 if(a>6)
1448 {
1449 if(S>1)
1450 {
1451 /*
1452 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1453 GSReLLA=G4QPDGCode().GetNuclMass(Z-2,N-2,S-2);
1454 */
1455 //if(GSReLLA+mAlph+ABarr+mLamb+mLamb<totMass) llaF=true;
1456 }
1457 if(N>2&&S>0)
1458 {
1459 /*
1460 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1461 GSReNLA=G4QPDGCode().GetNuclMass(Z-2,N-3,S-1);
1462 */
1463 //if(GSReNLA+mAlph+ABarr+mNeut+mLamb<totMass) nlaF=true;
1464 }
1465 if(Z>2&&S>0)
1466 {
1467 /*
1468 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1469 GSRePLA=G4QPDGCode().GetNuclMass(Z-3,N-2,S-1);
1470 */
1471 //if(GSRePLA+mAlph+SAPBarr+mProt+mLamb<totMass) plaF=true;
1472 }
1473 if(N>3)
1474 {
1475 /*
1476 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1477 GSReNNA=G4QPDGCode().GetNuclMass(Z-2,N-4,S);
1478 */
1479 //if(GSReNNA+mAlph+ABarr+mNeut+mNeut<totMass) nnaF=true;
1480 }
1481 if(Z>2&&N>2)
1482 {
1483 /*
1484 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1485 GSReNPA=G4QPDGCode().GetNuclMass(Z-3,N-3,S);
1486 */
1487 //if(GSReNPA+mAlph+SAPBarr+mProt+mNeut<totMass) npaF=true;
1488 }
1489 if(N>3)
1490 {
1491 /*
1492 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1493 GSRePPA=G4QPDGCode().GetNuclMass(Z-4,N-2,S);
1494 */
1495 //if(GSRePPA+mAlph+SPPABarr+mProt+mProt<totMass) ppaF=true;
1496 }
1497 if(a>9)
1498 {
1499 if(Z>3&&N>3&&S>0)
1500 {
1501 /*
1502 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1503 GSReLAA=G4QPDGCode().GetNuclMass(Z-4,N-4,S-1);
1504 */
1505 //if(GSReLAA+mLamb+mAlph+mAlph+SAABarr<totMass) laaF=true;
1506 }
1507 if(Z>3&&N>4)
1508 {
1509 /*
1510 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1511 GSReNAA=G4QPDGCode().GetNuclMass(Z-4,N-5,S);
1512 */
1513 //if(GSReNAA+mNeut+mAlph+mAlph+SAABarr<totMass) naaF=true;
1514 }
1515 if(Z>4&&N>3)
1516 {
1517 /*
1518 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1519 GSRePAA=G4QPDGCode().GetNuclMass(Z-5,N-4,S);
1520 */
1521 //if(GSRePAA+mProt+mAlph+mAlph+SAABarr<totMass) paaF=true;
1522 }
1523 if(a>12&&N>5&&Z>5)
1524 {
1525 /*
1526 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1527 GSReAAA=G4QPDGCode().GetNuclMass(Z-6,N-6,S);
1528 */
1529 //if(GSReAAA+mAlph+mAlph+mAlph+SAAABarr<totMass) aaaF=true;
1530 }
1531 }
1532 }
1533 if(N>3&&Z>3&&a>8)
1534 {
1535 aaQPDG =G4QPDGCode(90000000+1000*(1000*S+Z-4)+N-4);
1536 GSResAA=aaQPDG.GetMass();
1537#ifdef debug
1538 G4double s_value=GSResAA+mAlph+mAlph+SAABarr;
1541#endif
1542 if(GSResAA+mAlph+mAlph+SAABarr<totMass) aaFlag=true;
1543 }
1544 if(N>2&&a>5)
1545 {
1546 naQPDG =G4QPDGCode(90000000+1000*(1000*S+Z-2)+N-3);
1547 GSResNA=naQPDG.GetMass();
1548#ifdef debug
1549 G4double s_value=GSResNA+mAlph+mNeut;
1552#endif
1553 if(GSResNA+mNeut+mAlph+ABarr<totMass) naFlag=true;
1554 }
1555 if(S>0&&a>5)
1556 {
1557 laQPDG =G4QPDGCode(90000000+1000*(1000*S+Z-1002)+N-2);
1558 GSResLA=laQPDG.GetMass();
1559 if(GSResLA+mLamb+mAlph+ABarr<totMass) laFlag=true;
1560 }
1561 AQPDG =G4QPDGCode(90000000+1000*(1000*S+Z-2)+N-2);
1562 GSResNa=AQPDG.GetMass();
1563 mpls=GSResNa+mAlph;
1564 mmin=GSResNa-mAlph;
1565 ap2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1566 if(ap2m>=0.000001)
1567 {
1568 aFlag=true;
1569 aBnd=mAlph-GSMass+GSResNa; // Binding energy for ALPHA
1570 G4double eMax=sqrt(mA2+ap2m);
1571#ifdef debug
1574#endif
1575 aExcess=eMax-mAlph+aBnd; // Max Kin Energy from bottom
1576 }
1577 else aExcess=pBnd;
1578 }
1579 }
1580 if(N>0)
1581 {
1582 if(Z>0)
1583 {
1584 npQPDG=G4QPDGCode(90000000+1000*(1000*S+Z-1)+N-1);
1585 GSResNP=npQPDG.GetMass();
1586#ifdef debug
1587 G4double s_value=GSResNP+mNeut+mProt;
1590#endif
1591 if(GSResNP+mNeut+mProt+PBarr<totMass) npFlag=true;
1592 }
1593 if(N>1)
1594 {
1595 /*
1596 // DHW 16 June 2011: variable set but not used. See note at line 1317;
1597 GSReNNP=G4QPDGCode().GetNuclMass(Z-1,N-2,S);
1598 */
1599 //if(GSReNNP+mProt+PBarr+mNeut+mNeut<totMass) nnpF=true;
1600 }
1601 if(S>0)
1602 {
1603 /*
1604 // DHW 16 June 2011: variable set but not used. See note at line 1317;
1605 GSReNPL=G4QPDGCode().GetNuclMass(Z-1,N-1,S-1);
1606 */
1607 //if(GSReNPL+mProt+PBarr+mNeut+mLamb<totMass) nplF=true;
1608 }
1609 }
1610 if(S>0)
1611 {
1612 if(Z>0)
1613 {
1614 plQPDG=G4QPDGCode(90000000+1000*(1000*(S-1)+Z-1)+N);
1615 GSResPL=plQPDG.GetMass();
1616 if(GSResPL+mProt+PBarr+mLamb<totMass) plFlag=true;
1617 }
1618 if(S>1)
1619 {
1620 /*
1621 // DHW 16 June 2011: variable set but not used. See note at line 1317;
1622 GSRePLL=G4QPDGCode().GetNuclMass(Z-1,N,S-2);
1623 */
1624 //if(GSRePLL+mProt+PBarr+mLamb+mLamb<totMass) pllF=true;
1625 }
1626 }
1627 }
1632 if(N>0)
1633 {
1634 NQPDG=G4QPDGCode(90000000+1000*(1000*S+Z)+N-1);
1635 GSResNn=NQPDG.GetMass();
1636#ifdef debug
1639#endif
1640 G4double mpls=GSResNn+mNeut;
1641 G4double mmin=GSResNn-mNeut;
1642 np2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1643 if(np2m>=0.000001)
1644 {
1645 nFlag=true;
1646 nBnd=mNeut-GSMass+GSResNn; // Binding energy for neutron
1647 G4double eMax=sqrt(mN2+np2m);
1648#ifdef debug
1651#endif
1652 nExcess=eMax-mNeut+nBnd;
1653 }
1654 else nExcess=nBnd;
1655 if(N>1)
1656 {
1657 nnQPDG=G4QPDGCode(90000000+1000*(1000*S+Z)+N-2);
1658 GSResNN=nnQPDG.GetMass();
1659 if(GSResNN+mNeut+mNeut<totMass) nnFlag=true;
1660 if(N>2)
1661 {
1662 /*
1663 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1664 GSReNNN=G4QPDGCode().GetNuclMass(Z,N-3,S);
1665 */
1666 //if(GSReNNN+mNeut*3<totMass) nnnF=true;
1667 }
1668 if(S>0)
1669 {
1670 /*
1671 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1672 GSReNNL=G4QPDGCode().GetNuclMass(Z,N-2,S-1);
1673 */
1674 //if(GSReNNL+mNeut+mNeut+mLamb<totMass) nnlF=true;
1675 }
1676 }
1677 if(S>0)
1678 {
1679 nlQPDG=G4QPDGCode(90000000+1000*(1000*(S-1)+Z)+N-1);
1680 GSResNL=nlQPDG.GetMass();
1681 if(GSResNL+mNeut+mLamb<totMass) nlFlag=true;
1682 if(S>1)
1683 {
1684 /*
1685 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1686 GSReNLL=G4QPDGCode().GetNuclMass(Z,N-1,S-2);
1687 */
1688 //if(GSReNLL+mNeut+mLamb+mLamb<totMass) nllF=true;
1689 }
1690 }
1691 }
1696 if(S>0)
1697 {
1698 LQPDG=G4QPDGCode(90000000+1000*(1000*(S-1)+Z)+N);
1699 GSResNl=LQPDG.GetMass();
1700 G4double mpls=GSResNl+mLamb;
1701 G4double mmin=GSResNl-mLamb;
1702 lp2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1703 if(lp2m>=0.000001)
1704 {
1705 lFlag=true;
1706 lBnd=mLamb-GSMass+GSResNl; // Binding energy for lambda
1707 G4double eMax=sqrt(mL2+lp2m);
1708#ifdef debug
1711#endif
1712 lExcess=eMax-mLamb+lBnd;
1713 }
1714 else lExcess=lBnd;
1715 if(S>1)
1716 {
1717 llQPDG=G4QPDGCode(90000000+1000*(1000*(S-2)+Z)+N);
1718 GSResLL=llQPDG.GetMass();
1719 if(GSResLL+mLamb+mLamb<totMass) llFlag=true;
1720 if(S>2)
1721 {
1722 /*
1723 // DHW 16 June 2011: variable set but not used. See note at line 1317.
1724 GSReLLL=G4QPDGCode().GetNuclMass(Z,N,S-3);
1725 */
1726 //if(GSReLLL+mLamb*3<totMass) lllF=true;
1727 }
1728 }
1729 }
1734 //G4bool nTrF=nnnF||nnpF||nppF||nnlF||nplF||nllF; //Pos of 3-d baryon radiation after n
1735 //G4bool pTrF=nnpF||nppF||pppF||nplF||pplF||pllF; //Pos of 3-d baryon radiation after p
1736 //G4bool lTrF=nnlF||nplF||pplF||nllF||pllF||lllF; //Pos of 3-d baryon radiation after l
1737 //G4bool aTrF=nnaF||npaF||ppaF||nlaF||plaF||llaF; //Pos of 3-d baryon radiation after a
1739 //G4bool thdB = nTrF||pTrF||lTrF||aTrF||naaF||paaF||laaF||aaaF;// Pos to radiate three
1740#ifdef debug
1742 <<nnFlag<<", np="<<npFlag<<",pp="<<ppFlag<<",pa="<<paFlag<<",na="<<naFlag<<",aa="
1743 <<aaFlag<<G4endl;
1744#endif
1745 G4QPDGCode bQPDG;
1746 G4QPDGCode rQPDG;
1747 if(secB) // Decay in two baryons is possible
1748 //if(thdB) //@@CHECK@@ Decay in three baryons is possible
1749 {
1750 if(!nSecF) nFlag=false;
1751 if(!pSecF) pFlag=false;
1752 if(!lSecF) lFlag=false;
1753 if(!aSecF) aFlag=false;
1754#ifdef debug
1756#endif
1758 if(nFlag&&nExcess>maxE) maxE=nExcess;
1759 if(pFlag&&pExcess>maxE) maxE=pExcess;
1760 if(lFlag&&lExcess>maxE) maxE=lExcess;
1761 if(lFlag&&aExcess>maxE) maxE=aExcess;
1763 if(pFlag)pMin+= PBarr; // Add Coulomb Barrier for protons
1767 if(aFlag)aMin+= ABarr; // Add Coulomb Barrier for alpha
1769 if(nFlag&&nMin<minE) minE=nMin;
1770 if(pFlag&&pMin<minE) minE=pMin;
1771 if(lFlag&&lMin<minE) minE=lMin;
1772 if(evalph&&aFlag&&aMin<minE) minE=aMin;
1773
1774#ifdef debug
1776 <<",sE="<<lExcess<<">"<<lMin<<",E="<<aExcess<<">"<<aMin<<",miE="<<minE<<"<maE="
1777 <<maxE<<G4endl;
1778#endif
1779 // @@ Here one can put a condition for the Baryon Gun
1780 G4int cntr= 0;
1781 //G4int cntm= 27;
1782 //G4int cntm= 72; // Important difference !!DOn't change
1783 //G4int cntm= 80; // Important difference !!DOn'tChange"IsoNuclei"
1784 //G4int cntm= 90; // Important difference !!DOn'tChange "Lept/Hyper"
1786 if( ( (pFlag && pExcess > pMin) ||
1787 (nFlag && nExcess > nMin) ||
1788 (lFlag && lExcess > lMin) ||
1789 (aFlag && aExcess > aMin) ) && minE<maxE )
1790 {
1794 if(mi<0.)
1795 {
1796 uW-=mi;
1797 mm_value-=mi;
1798 mi=0.;
1799 }
1800 //G4bool good=true;
1801 if(ma<mm_value)
1802 {
1803 ma=mm_value;
1804 //good=false;
1805 }
1806#ifdef debug
1808 <<ma<<G4endl;
1809#endif
1812 //G4double xCa=maSht-coSht*log(a); // Maximal value of x (approximation)
1813 //G4double xMa=xCa; // Maximal value of x
1814 //if(xMm<xMa) xMa=xMm;
1815 G4double xMa=xMm;
1816 if(xMa>1.)xMa=1.;
1817 if(xMi<0.)xMi=0.;
1818 if(xMi>xMa)
1819 {
1821 return false;
1822 }
1823 xMi=sqrt(xMi); // @@ ?
1824 xMa=sqrt(xMa); // @@ ?
1825#ifdef debug
1827#endif
1830#ifdef debug
1832#endif
1834 G4double maxR=pow(1.-xMi*xMi,powr);
1835#ifdef debug
1837#endif
1839 G4int PDG = 0;
1841 while(cond&&cntr<cntm)
1842 {
1844 //if(!good)R = maxR;
1845 G4double x2= 1.-pow(R,revP);
1846 G4double x = sqrt(x2);
1847 if(x<xMi||x>xMa)
1848 {
1849#ifdef debug
1851#endif
1852 if(x<xMi) x=xMi;
1853 else x=xMa;
1854 x2 = x*x;
1855 }
1857 //if(rn<x/xMa||!good)
1858 if(rn<x/xMa) // Randomization cut
1859 {
1860 tk= ma*x2-uW; // Kinetic energy of the fragment
1861 G4double psum =0.;
1863#ifdef debug
1865 <<aMin<<G4endl;
1866#endif
1867 if(pFlag&&tk>pMin)
1868 {
1869 G4double kin=tk-pBnd;
1870 //if(barf) kin-=PBarr; //@@ This is a mistake
1871#ifdef debug
1874#endif
1875 zCBPP=Z*CoulBarPenProb(PBarr,kin,1,1)*sqrt(kin);
1876 }
1877 psum+=zCBPP;
1879 if(nFlag&&tk>nMin)
1880 {
1881 G4double kin=tk-nBnd;
1882#ifdef debug
1884#endif
1885 nCBPP=N*CoulBarPenProb(0.,kin,0,1)*sqrt(kin);
1886 }
1887 psum+=nCBPP;
1888 nCBPP+=zCBPP;
1890 if(lFlag&&tk>lMin)
1891 {
1892 G4double kin=tk-lBnd;
1893#ifdef debug
1895#endif
1896 lCBPP=S*CoulBarPenProb(0.,kin,0,1)*sqrt(kin);
1897 }
1898 psum+=lCBPP;
1899 lCBPP+=nCBPP;
1900 if(evalph&&aFlag&&tk>aMin)
1901 {
1902 G4double kin=tk-aBnd;
1903 //if(barf) kin-=ABarr; //@@ This is a mistake
1904#ifdef debug
1907#endif
1908 psum+=CoulBarPenProb(ABarr,kin,2,4)*sqrt(kin)*evalph*Z*(Z-1)*N*(N-1)
1909 *6/a1/(a-2)/(a-3);
1910 }
1912#ifdef debug
1914 <<psum<<G4endl;
1915#endif
1916 cond = false;
1917 if (r&&r>lCBPP)
1918 {
1919#ifdef debug
1921#endif
1922 PDG=aPDG;
1923 }
1924 else if(r&&r>nCBPP&&r<=lCBPP)
1925 {
1926#ifdef debug
1928#endif
1929 PDG=lPDG;
1930 }
1931 else if(r&&r>zCBPP&&r<=nCBPP)
1932 {
1933#ifdef debug
1935#endif
1936 PDG=nPDG;
1937 }
1938 else if(r&&r<=zCBPP)
1939 {
1940#ifdef debug
1942#endif
1943 PDG=pPDG;
1944 }
1945 else cond=true;
1946 }
1947#ifdef debug
1950#endif
1951 cntr++;
1952 }
1953 if(cntr<cntm) // => Succeeded to find the evaporation channel
1954 {
1955 G4double p2=0.;
1956 if (PDG==aPDG)
1957 {
1958 tk-=aBnd-mAlph; // Pays for binding and convert to total energy
1959 p2=tk*tk-mA2;
1960 if(p2>ap2m)
1961 {
1962 p2=ap2m;
1963 tk=sqrt(p2+mA2);
1964 }
1965 eMass=mAlph;
1966 bQPDG=aQPDG;
1967 rQPDG=AQPDG;
1968 }
1969 else if(PDG==pPDG)
1970 {
1971 tk-=pBnd-mProt; // Pays for binding and convert to total energy
1972 p2=tk*tk-mP2;
1973 if(p2>pp2m)
1974 {
1975 p2=pp2m;
1976 tk=sqrt(p2+mP2);
1977 }
1978 eMass=mProt;
1979 bQPDG=pQPDG;
1980 rQPDG=PQPDG;
1981 }
1982 else if(PDG==nPDG)
1983 {
1984 tk-=nBnd-mNeut; // Pays for binding and convert to total energy
1985 p2=tk*tk-mN2;
1986#ifdef debug
1988#endif
1989 if(p2>np2m)
1990 {
1991 p2=np2m;
1992 tk=sqrt(p2+mN2);
1993 }
1994 eMass=mNeut;
1995 bQPDG=nQPDG;
1996 rQPDG=NQPDG;
1997 }
1998 else if(PDG==lPDG)
1999 {
2000 tk-=lBnd-mLamb; // Pays for binding and convert to total energy
2001 p2=tk*tk-mL2;
2002 if(p2>lp2m)
2003 {
2004 p2=lp2m;
2005 tk=sqrt(p2+mL2);
2006 }
2007 eMass=mLamb;
2008 bQPDG=lQPDG;
2009 rQPDG=LQPDG;
2010 }
2012 G4double rEn=totMass-tk;
2013 G4double rEn2=rEn*rEn;
2014 if (rEn2 > p2) rMass=sqrt(rEn2-p2); // Mass of the Residual Nucleus
2015 else rMass=0.0;
2016 // Find out if the ResidualNucleus is below of the SecondBaryonDecayLimit
2017 //@@ Calculate it depending on PDG !!!!!!!
2018 G4bool nnCond = !nnFlag || (nnFlag && GSResNN+mNeut > rMass);
2019 G4bool npCond = !npFlag || (npFlag && GSResNP+mProt+PBarr > rMass);
2020 G4bool nlCond = !nlFlag || (nlFlag && GSResNL+mLamb > rMass);
2021 G4bool naCond = !naFlag || (naFlag && GSResNA+mAlph+ABarr > rMass);
2022 G4bool pnCond = !npFlag || (npFlag && GSResNP+mNeut > rMass);
2023 if(barf) pnCond = !npFlag || (npFlag && GSResNP+mNeut+PBarr > rMass);
2024 G4bool ppCond = !ppFlag || (ppFlag && GSResPP+mProt+PPBarr > rMass);
2025 if(barf) ppCond = !ppFlag || (ppFlag && GSResPP+mProt+SPPBarr > rMass);
2026 G4bool plCond = !plFlag || (plFlag && GSResPL+mLamb > rMass);
2027 if(barf) plCond = !plFlag || (plFlag && GSResPL+mLamb+PBarr > rMass);
2028 G4bool paCond = !paFlag || (paFlag && GSResPA+mAlph+APBarr > rMass);
2029 if(barf) paCond = !paFlag || (paFlag && GSResPA+mAlph+SAPBarr > rMass);
2030 G4bool lnCond = !nlFlag || (nlFlag && GSResNL+mNeut > rMass);
2031 G4bool lpCond = !plFlag || (plFlag && GSResPL+mProt+PBarr > rMass);
2032 G4bool llCond = !llFlag || (llFlag && GSResLL+mLamb > rMass);
2033 G4bool laCond = !laFlag || (laFlag && GSResLA+mAlph+ABarr > rMass);
2034 G4bool anCond = !naFlag || (naFlag && GSResNA+mNeut > rMass);
2035 if(barf) anCond = !naFlag || (naFlag && GSResNA+mNeut+ABarr > rMass);
2036 G4bool apCond = !paFlag || (paFlag && GSResPA+mProt+PABarr > rMass);
2037 if(barf) apCond = !paFlag || (paFlag && GSResPA+mProt+SAPBarr > rMass);
2038 G4bool alCond = !laFlag || (laFlag && GSResLA+mLamb > rMass);
2039 if(barf) alCond = !laFlag || (laFlag && GSResLA+mLamb+ABarr > rMass);
2040 G4bool aaCond = !aaFlag || (aaFlag && GSResAA+mAlph+AABarr > rMass);
2041 if(barf) aaCond = !aaFlag || (aaFlag && GSResAA+mAlph+SAABarr > rMass);
2042#ifdef debug
2045 <<GSResPP+mProt+PPBarr<<"("<<ppCond<<"),PL="
2046 <<GSResPL+mLamb<<"("<<plCond<<"),PA="
2047 <<GSResPA+mAlph+APBarr<<"("<<paCond;
2049 <<GSResNP+mProt+PBarr<<"("<<npCond<<"),NL="
2050 <<GSResNL+mLamb<<"("<<nlCond<<"),NA="
2051 <<GSResNA+mAlph+ABarr<<"("<<naCond;
2053 <<GSResPL+mProt+PBarr<<"("<<lpCond<<"),LL="
2054 <<GSResLL+mLamb<<"("<<llCond<<"),LA="
2055 <<GSResLA+mAlph+ABarr<<"("<<laCond;
2057 <<GSResPA+mProt+PABarr<<"("<<apCond<<"),AL="
2058 <<GSResLA+mLamb<<"("<<alCond<<"),AA="
2059 <<GSResAA+mAlph+AABarr<<"("<<aaCond;
2061#endif
2062 three=false; // Flag of b+b+ResNuc decay
2063 //if(3>2)three=false; // @@@@@@@@@@@@@@@@@@
2064 //else if(PDG==pPDG&&(pnCond&&ppCond&&plCond&&paCond)) // @@@@@@@@@@@@@@@@@@@
2065 if(PDG==pPDG&&(pnCond&&ppCond&&plCond&&paCond))//p+RN decay, p+b+RN dec is closed
2066 {
2067#ifdef debug
2069 <<G4endl;
2070#endif
2071 fMass=mProt;
2072 fQPDG=pQPDG;
2073 G4double nLim=0.;
2074 if(N&&GSResNP!=GSMass&&fMass+PBarr+mNeut+GSResNP<totMass)
2075 {
2076 if(barf) nLim+=(N+N)*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2077 else nLim+=(N+N)*pow(totMass-mNeut-mProt-GSResNP,2);
2078 }
2079 G4double zLim=nLim;
2080 if(Z>1&&GSResPP!=GSMass&&fMass+mProt+SPPBarr+GSResPP<totMass)
2081 {
2082 if(barf) zLim+=(Z-1)*pow(totMass-mProt-mProt-SPPBarr-GSResPP,2);
2083 else zLim+=(Z-1)*pow(totMass-mProt-mProt-GSResPP,2);
2084 }
2085 G4double sLim=zLim;
2086 if(S&&GSResPL!=GSMass&&fMass+PBarr+mLamb+GSResPL<totMass)
2087 {
2088 if(barf) sLim+=(S+S)*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2089 else sLim+=(S+S)*pow(totMass-mProt-mLamb-GSResPL,2);
2090 }
2091 G4double aLim=sLim;
2092 if(evalph&&Z>2&&N>1&&a>4&&GSResPL!=GSMass&&fMass+SAPBarr+mAlph+GSResPA<totMass)
2093 {
2094 if(barf) aLim+=pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2)*evalph*
2095 (Z-1)*(Z-2)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2096 else aLim+=pow(totMass-mProt-mAlph-GSResPA,2)*evalph*(Z-1)*(Z-2)*N*(N-1)
2097 *12/(a-2)/(a-3)/(a-4);
2098 }
2100#ifdef debug
2102 <<aLim<<G4endl;
2103#endif
2104 three=true; // Flag of b+b+ResNuc decay
2105 if(!aLim) three=false;
2106 else if(r>sLim)
2107 {
2108 eMass = mAlph;
2109 dbQPDG= PAQPDG;
2110 rMass = GSResPA;
2111 rQPDG = paQPDG;
2112#ifdef debug
2114#endif
2115 }
2116 else if(zLim<sLim&&r>zLim&&r<=sLim)
2117 {
2118 eMass = mLamb;
2119 dbQPDG= PLQPDG;
2120 rMass = GSResPL;
2121 rQPDG = plQPDG;
2122#ifdef debug
2124#endif
2125 }
2126 else if(nLim<zLim&&r>nLim&&r<=zLim)
2127 {
2128 eMass = mProt;
2129 dbQPDG= PPQPDG;
2130 rMass = GSResPP;
2131 rQPDG = ppQPDG;
2132#ifdef debug
2134#endif
2135 }
2136 else if(r<=nLim)
2137 {
2138 eMass = mNeut;
2139 dbQPDG= NPQPDG;
2140 rMass = GSResNP;
2141 rQPDG = npQPDG;
2142#ifdef debug
2144#endif
2145 }
2146 else three=false;
2147 }
2148 else if(PDG==nPDG&&(nnCond&&npCond&&nlCond&&naCond)) // n+b+RN decay can't happen
2149 { //@@ Take into account Coulomb Barier Penetration Probability
2150#ifdef debug
2152 <<G4endl;
2153#endif
2154 fMass=mNeut;
2155 fQPDG=nQPDG;
2156 G4double nLim=0.;
2157 if(N>1&&GSResNN!=GSMass&&fMass+mNeut+GSResNN<totMass)
2158 nLim+=(N-1)*pow(totMass-mNeut-mNeut-GSResNN,2);
2159 G4double zLim=nLim;
2160 if(Z&&GSResNP!=GSMass&&fMass+mProt+PBarr+GSResNP<totMass)
2161 {
2162 if(barf) zLim+=(Z+Z)*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2163 else zLim+=(Z+Z)*pow(totMass-mNeut-mProt-GSResNP,2);
2164 }
2165 G4double sLim=zLim;
2166 if(S&&GSResNL!=GSMass&&fMass+mLamb+GSResNL<totMass)
2167 sLim+=(S+S)*pow(totMass-mNeut-mLamb-GSResNL,2);
2168 G4double aLim=sLim;
2169 if(evalph&&Z>1&&N>2&&a>4&&GSResNA!=GSMass&&fMass+mAlph+ABarr+GSResNA<totMass)
2170 {
2171 if(barf) aLim+=pow(totMass-mNeut-mAlph-ABarr-GSResNA,2)*
2172 evalph*Z*(Z-1)*(N-1)*(N-2)*12/(a-2)/(a-3)/(a-4);
2173 else aLim+=pow(totMass-mNeut-mAlph-GSResNA,2)*
2174 evalph*Z*(Z-1)*(N-1)*(N-2)*12/(a-2)/(a-3)/(a-4);
2175 }
2177#ifdef debug
2179 <<G4endl;
2180#endif
2181 three=true; // Flag of b+b+ResNuc decay
2182 if(!aLim) three=false;
2183 else if(r>sLim)
2184 {
2185 eMass = mAlph;
2186 dbQPDG= NAQPDG;
2187 rMass = GSResNA;
2188 rQPDG = naQPDG;
2189#ifdef debug
2191#endif
2192 }
2193 else if(zLim<sLim&&r>zLim&&r<=sLim)
2194 {
2195 eMass = mLamb;
2196 dbQPDG= NLQPDG;
2197 rMass = GSResNL;
2198 rQPDG = nlQPDG;
2199#ifdef debug
2201#endif
2202 }
2203 else if(nLim<zLim&&r>nLim&&r<=zLim)
2204 {
2205 eMass = mProt;
2206 dbQPDG= NPQPDG;
2207 rMass = GSResNP;
2208 rQPDG = npQPDG;
2209#ifdef debug
2211#endif
2212 }
2213 else if(r<=nLim)
2214 {
2215 eMass = mNeut;
2216 dbQPDG= NNQPDG;
2217 rMass = GSResNN;
2218 rQPDG = nnQPDG;
2219#ifdef debug
2221#endif
2222 }
2223 else three=false;
2224 }
2225 else if(PDG==lPDG&&(lnCond&&lpCond&&llCond&&laCond)) // l+b+RN decay can't happen
2226 { //@@ Take into account Coulomb Barier Penetration Probability
2227#ifdef debug
2229 <<G4endl;
2230#endif
2231 fMass=mLamb;
2232 fQPDG=lQPDG;
2233 G4double nLim=0.;
2234 if(N&&GSResNL!=GSMass&&fMass+mNeut+GSResNL<totMass)
2235 nLim+=(N+N)*pow(totMass-mNeut-mLamb-GSResNL,2);
2236 G4double zLim=nLim;
2237 if(Z&&GSResPL!=GSMass&&fMass+mProt+PBarr+GSResPL<totMass)
2238 {
2239 if(barf) zLim+=(Z+Z)*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2240 else zLim+=(Z+Z)*pow(totMass-mProt-mLamb-GSResPL,2);
2241 }
2242 G4double sLim=zLim;
2243 if(S>1&&GSResLL!=GSMass&&fMass+mLamb+GSResLL<totMass)
2244 sLim+=(S-1)*pow(totMass-mLamb-mLamb-GSResLL,2);
2245 G4double aLim=sLim;
2246 if(evalph&&Z>1&&N>1&&a>4&&GSResLA!=GSMass&&fMass+mAlph+ABarr+GSResLA<totMass)
2247 {
2248 if(barf) aLim+=pow(totMass-mLamb-mAlph-ABarr-GSResLA,2)*
2249 evalph*Z*(Z-1)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2250 else aLim+=pow(totMass-mLamb-mAlph-GSResLA,2)*
2251 evalph*Z*(Z-1)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2252 }
2254#ifdef debug
2256 <<G4endl;
2257#endif
2258 three=true; // Flag of b+b+ResNuc decay
2259 if(!aLim) three=false;
2260 else if(r>sLim)
2261 {
2262 eMass = mAlph;
2263 dbQPDG= LAQPDG;
2264 rMass = GSResLA;
2265 rQPDG = laQPDG;
2266#ifdef debug
2268#endif
2269 }
2270 else if(zLim<sLim&&r>zLim&&r<=sLim)
2271 {
2272 eMass = mLamb;
2273 dbQPDG= LLQPDG;
2274 rMass = GSResLL;
2275 rQPDG = llQPDG;
2276#ifdef debug
2278#endif
2279 }
2280 else if(nLim<zLim&&r>nLim&&r<=zLim)
2281 {
2282 eMass = mProt;
2283 dbQPDG= PLQPDG;
2284 rMass = GSResPL;
2285 rQPDG = plQPDG;
2286#ifdef debug
2288#endif
2289 }
2290 else if(r<=nLim)
2291 {
2292 eMass = mNeut;
2293 dbQPDG= NLQPDG;
2294 rMass = GSResNL;
2295 rQPDG = nlQPDG;
2296#ifdef debug
2298#endif
2299 }
2300 else three=false;
2301 }
2302 else if(PDG==aPDG&&(anCond&&apCond&&alCond&&aaCond)) // a+b+RN decay can't happen
2303 { //@@ Take into account Coulomb Barier Penetration Probability
2304#ifdef debug
2306 <<G4endl;
2307#endif
2308 fMass=mAlph;
2309 fQPDG=aQPDG;
2310 G4double nLim=0.;
2311 if(N>2&&GSResNA!=GSMass&&fMass+mNeut+ABarr+GSResNA<totMass)
2312 {
2313 if(barf) nLim+=(N-2)*pow(totMass-mNeut-mAlph-ABarr-GSResNA,2);
2314 else nLim+=(N-2)*pow(totMass-mNeut-mAlph-GSResNA,2);
2315 }
2316 G4double zLim=nLim;
2317 if(Z>2&&GSResPA!=GSMass&&fMass+mProt+SAPBarr+GSResPA<totMass)
2318 {
2319 if(barf) zLim+=(Z-2)*pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2);
2320 else zLim+=(Z-2)*pow(totMass-mProt-mAlph-GSResPA,2);
2321 }
2322 G4double sLim=zLim;
2323 if(S&&GSResLA!=GSMass&&fMass+mLamb+ABarr+GSResLA<totMass)
2324 {
2325 if(barf) sLim+=S*pow(totMass-mLamb-mAlph-ABarr-GSResLA,2);
2326 else sLim+=S*pow(totMass-mLamb-mAlph-GSResLA,2);
2327 }
2328 G4double aLim=sLim;
2329 if(evalph&&Z>3&&N>3&&a>7&&GSResAA!=GSMass&&fMass+mAlph+SAABarr+GSResAA<totMass)
2330 {
2331 if(barf) aLim+=pow(totMass-mAlph-mAlph-SAABarr-GSResAA,2)*
2332 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*12/(a-5)/(a-6)/(a-7);
2333 else aLim+=pow(totMass-mAlph-mAlph-GSResAA,2)*
2334 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*12/(a-5)/(a-6)/(a-7);
2335 }
2337#ifdef debug
2339 <<G4endl;
2340#endif
2341 three=true; // Flag of b+b+ResNuc decay
2342 if(!aLim) three=false;
2343 else if(r>sLim)
2344 {
2345 eMass = mAlph;
2346 dbQPDG= AAQPDG;
2347 rMass = GSResAA;
2348 rQPDG = aaQPDG;
2349#ifdef debug
2351#endif
2352 }
2353 else if(zLim<sLim&&r>zLim&&r<=sLim)
2354 {
2355 eMass = mLamb;
2356 dbQPDG= LAQPDG;
2357 rMass = GSResLA;
2358 rQPDG = laQPDG;
2359#ifdef debug
2361#endif
2362 }
2363 else if(nLim<zLim&&r>nLim&&r<=zLim)
2364 {
2365 eMass = mProt;
2366 dbQPDG= PAQPDG;
2367 rMass = GSResPA;
2368 rQPDG = paQPDG;
2369#ifdef debug
2371#endif
2372 }
2373 else if(r<=nLim)
2374 {
2375 eMass = mNeut;
2376 dbQPDG= NAQPDG;
2377 rMass = GSResNA;
2378 rQPDG = naQPDG;
2379#ifdef debug
2381#endif
2382 }
2383 else three=false;
2384 }
2385 else three=false;
2386 if(rMass<1600.)
2387 {
2388 if (rQPDG==pQPDG)rMass=mProt;
2389 else if(rQPDG==nQPDG)rMass=mNeut;
2390 else if(rQPDG==lQPDG)rMass=mLamb;
2391 }
2392#ifdef debug
2395#endif
2396 }
2397 }
2398 else // =-------------=> Just decay in a baryon and a residual (to avoid gamma-decay)
2399 { //@@ Take into account Coulomb Barier Penetration Probability (?? - Emergency)
2400 G4double nLim=0.;
2401 if(nFlag&&mNeut+GSResNn<totMass)
2402 {
2403 G4double ken=totMass-mNeut-GSResNn;
2404 nLim+=N*CoulBarPenProb(0.,ken,0,1)*sqrt(ken);
2405 }
2406 G4double zLim=nLim;
2407 if(pFlag&&mProt+PBarr+GSResNp<totMass)
2408 {
2409 G4double ken=totMass-mProt-GSResNp;
2410 if(barf) ken-=PBarr;
2411 zLim+=Z*CoulBarPenProb(PBarr,ken,1,1)*sqrt(ken);
2412 }
2413 G4double sLim=zLim;
2414 if(lFlag&&mLamb+GSResNl<totMass)
2415 {
2416 G4double ken=totMass-mLamb-GSResNl;
2417 sLim+=S*CoulBarPenProb(0.,ken,0,1)*sqrt(ken);
2418 }
2419 G4double aLim=sLim;
2420 if(evalph&&aFlag&&mAlph+GSResNa+ABarr<totMass)
2421 {
2422 G4double ken=totMass-mAlph-GSResNa;
2423 if(barf) ken-=ABarr;
2424 aLim+=CoulBarPenProb(ABarr,ken,2,4)*sqrt(ken)*evalph*Z*(Z-1)*N*(N-1)
2425 *6/a1/(a-2)/(a-3);
2426 }
2428#ifdef debug
2431#endif
2432 if (aFlag&&r>sLim)
2433 {
2434 bQPDG=aQPDG;
2435 eMass=mAlph;
2436 rQPDG=AQPDG;
2437 rMass=GSResNa;
2438 }
2439 else if(lFlag&&r>=zLim&&r<=sLim&&zLim<sLim)
2440 {
2441 bQPDG=lQPDG;
2442 eMass=mLamb;
2443 rQPDG=LQPDG;
2444 rMass=GSResNl;
2445 }
2446 else if(nFlag&&r>=nLim&&r<=zLim&&nLim<zLim)
2447 {
2448 bQPDG=pQPDG;
2449 eMass=mProt;
2450 rQPDG=PQPDG;
2451 rMass=GSResNp;
2452 }
2453 else if(pFlag&&r<nLim)
2454 {
2455 bQPDG=nQPDG;
2456 eMass=mNeut;
2457 rQPDG=NQPDG;
2458 rMass=GSResNn;
2459 }
2460 else
2461 {
2462#ifdef debug
2464#endif
2465 bQPDG=gQPDG;
2466 rQPDG=GetQPDG();
2467 eMass=0.;
2468 rMass=GSMass;
2469 }
2470#ifdef debug
2472#endif
2473 }
2474 if(three) // Decay in two baryons + Residual Nucleus
2475 {
2476#ifdef debug
2478#endif
2479 h1mom=G4LorentzVector(0.,0.,0.,eMass);
2480 h2mom=G4LorentzVector(0.,0.,0.,rMass);
2481 h3mom=G4LorentzVector(0.,0.,0.,fMass);
2483 {
2484#ifdef debug
2487#endif
2488 return false;
2489 }
2490 h1mom+=h3mom;
2491 bQPDG=dbQPDG;
2492 }
2493 else
2494 {
2495 if(eMass+rMass<totMass&&cntr<cntm)
2496 {
2497#ifdef debug
2500#endif
2501 if(rMass<1600.)
2502 {
2503 if (rQPDG==pQPDG)rMass=mProt;
2504 else if(rQPDG==nQPDG)rMass=mNeut;
2505 else if(rQPDG==lQPDG)rMass=mLamb;
2506 }
2507 h1mom=G4LorentzVector(0.,0.,0.,eMass);
2508 h2mom=G4LorentzVector(0.,0.,0.,rMass);
2509 }
2510 else if(totMass>mNeut+GSResNn) // Neutron if 2-Decay failed
2511 {
2512#ifdef debug
2514#endif
2515 bQPDG=nQPDG;
2516 rQPDG=NQPDG;
2517 h1mom=G4LorentzVector(0.,0.,0.,mNeut);
2518 h2mom=G4LorentzVector(0.,0.,0.,GSResNn);
2519 }
2520 else if(totMass>mProt+PBarr+GSResNp) // Proton if 2-Decay failed
2521 {
2522#ifdef debug
2524#endif
2525 bQPDG=pQPDG;
2526 rQPDG=PQPDG;
2527 h1mom=G4LorentzVector(0.,0.,0.,mProt);
2528 h2mom=G4LorentzVector(0.,0.,0.,GSResNp);
2529 }
2530 else if(totMass>mAlph+ABarr+GSResNa) // Alpha if 2-Decay failed
2531 {
2532#ifdef debug
2534#endif
2535 bQPDG=aQPDG;
2536 rQPDG=AQPDG;
2537 h1mom=G4LorentzVector(0.,0.,0.,mAlph);
2538 h2mom=G4LorentzVector(0.,0.,0.,GSResNa);
2539 }
2540 else if(totMass>GSMass) // Photon if 2-Decay failed
2541 {
2542#ifdef debug
2544#endif
2545 bQPDG=gQPDG;
2546 rQPDG=GetQPDG();
2547 h1mom=G4LorentzVector(0.,0.,0.,0.);
2548 h2mom=G4LorentzVector(0.,0.,0.,GSMass);
2549 }
2550 else
2551 {
2553 return false;
2554 }
2555 }
2556 }
2557 else // ==> Decay in 3 Baryons + Residual is impossible at this point
2558 {
2559 if(secB) // Decay in 2Baryons(2a,a+bary)+ResidN is possible
2560 //if(2>3)
2561 {
2562#ifdef debug
2564#endif
2566 //@@ Coulomb Barrier penetration can be added
2567 G4double alp=0.;
2568 if(aSecF)alp=evalph*Z*(Z-1)*N*(N-1)*10/(a-2)/(a-3)/(a-4);
2569 G4double nnLim=0.;
2570 if(nnFlag&&totMass>mNeut+mNeut+GSResNN)
2571 nnLim+=N*(N-1)*pow(totMass-mNeut-mNeut-GSResNN,2);
2572 G4double nzLim=nnLim;
2573 if(npFlag&&totMass>mNeut+mProt+PBarr+GSResNP)
2574 {
2575 if(barf) nzLim+=2*N*Z*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2576 else nzLim+=2*N*Z*pow(totMass-mNeut-mProt-GSResNP,2);
2577 }
2578 G4double zzLim=nzLim;
2579 if(ppFlag&&totMass>mProt+mProt+SPPBarr+GSResPP)
2580 {
2581 if(barf) zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-SPPBarr-GSResPP,2);
2582 else zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-GSResPP,2);
2583 }
2584 G4double nlLim=zzLim;
2585 if(nlFlag&&totMass>mNeut+mLamb+GSResNL)
2586 nlLim+=2*N*S*pow(totMass-mNeut-mLamb-GSResNL,2);
2587 G4double zlLim=nlLim;
2588 if(plFlag&&totMass>mProt+PBarr+mLamb+GSResPL)
2589 {
2590 if(barf) zlLim+=2*Z*S*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2591 else zlLim+=2*Z*S*pow(totMass-mProt-mLamb-GSResPL,2);
2592 }
2593 G4double llLim=zlLim;
2594 if(llFlag&&totMass>mLamb+mLamb+GSResLL)
2595 llLim+=S*(S-1)*pow(totMass-mLamb-mLamb-GSResLL,2);
2596 G4double naLim=llLim;
2597 if(naFlag&&totMass>mNeut+mAlph+ABarr+GSResNA)
2598 {
2599 if(barf) naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-ABarr-GSResNA,2);
2600 else naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-GSResNA,2);
2601 }
2602 G4double zaLim=naLim;
2603 if(paFlag&&totMass>mProt+SAPBarr+mAlph+GSResPA)
2604 {
2605 if(barf) zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2);
2606 else zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-GSResPA,2);
2607 }
2608 G4double laLim=zaLim;
2609 if(laFlag&&totMass>mLamb+mAlph+ABarr+GSResLA)
2610 {
2611 if(barf) laLim+=S*alp*pow(totMass-mLamb-mAlph-ABarr-GSResLA,2);
2612 else laLim+=S*alp*pow(totMass-mLamb-mAlph-GSResLA,2);
2613 }
2614 G4double aaLim=laLim;
2615 if(evalph&&aaFlag&&totMass>mAlph+mAlph+SAABarr+GSResAA)
2616 {
2617 if(barf) aaLim+=alp*pow(totMass-mAlph-mAlph-SAABarr-GSResAA,2)*
2618 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7);
2619 else aaLim+=alp*pow(totMass-mAlph-mAlph-GSResAA,2)*
2620 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7);
2621 }
2623 if (aaLim&&laLim<r)
2624 {
2625 dbQPDG= AAQPDG;
2626 eMass=mAlph;
2627 fMass=mAlph;
2628 rQPDG=aaQPDG;
2629 rMass=GSResAA;
2630#ifdef debug
2632#endif
2633 }
2634 else if(aaLim&&zaLim<r&&r<laLim)
2635 {
2636 dbQPDG= LAQPDG;
2637 eMass=mAlph;
2638 fMass=mLamb;
2639 rQPDG=laQPDG;
2640 rMass=GSResLA;
2641#ifdef debug
2643#endif
2644 }
2645 else if(aaLim&&naLim<r&&r<zaLim)
2646 {
2647 dbQPDG= PAQPDG;
2648 eMass=mAlph;
2649 fMass=mProt;
2650 rQPDG=paQPDG;
2651 rMass=GSResPA;
2652#ifdef debug
2654#endif
2655 }
2656 else if(aaLim&&llLim<r&&r<naLim)
2657 {
2658 dbQPDG= NAQPDG;
2659 eMass=mAlph;
2660 fMass=mNeut;
2661 rQPDG=naQPDG;
2662 rMass=GSResNA;
2663#ifdef debug
2665#endif
2666 }
2667 else if(aaLim&&zlLim<r&&r<llLim)
2668 {
2669 dbQPDG= LLQPDG;
2670 eMass=mLamb;
2671 fMass=mLamb;
2672 rQPDG=llQPDG;
2673 rMass=GSResLL;
2674#ifdef debug
2676#endif
2677 }
2678 else if(aaLim&&nlLim<r&&r<zlLim)
2679 {
2680 dbQPDG= PLQPDG;
2681 eMass=mLamb;
2682 fMass=mProt;
2683 rQPDG=plQPDG;
2684 rMass=GSResPL;
2685#ifdef debug
2687#endif
2688 }
2689 else if(aaLim&&zzLim<r&&r<nlLim)
2690 {
2691 dbQPDG= NLQPDG;
2692 eMass=mLamb;
2693 fMass=mNeut;
2694 rQPDG=nlQPDG;
2695 rMass=GSResNL;
2696#ifdef debug
2698#endif
2699
2700 }
2701 else if(aaLim&&nzLim<r&&r<zzLim)
2702 {
2703 dbQPDG= PPQPDG;
2704 eMass=mProt;
2705 fMass=mProt;
2706 rQPDG=ppQPDG;
2707 rMass=GSResPP;
2708#ifdef debug
2710#endif
2711 }
2712 else if(aaLim&&nnLim<r&&r<nzLim)
2713 {
2714 dbQPDG= NPQPDG;
2715 eMass=mNeut;
2716 fMass=mProt;
2717 rQPDG=npQPDG;
2718 rMass=GSResNP;
2719#ifdef debug
2721#endif
2722 }
2723 else if(aaLim&&r<nnLim)
2724 {
2725 dbQPDG= NNQPDG;
2726 eMass=mNeut;
2727 fMass=mNeut;
2728 rQPDG=nnQPDG;
2729 rMass=GSResNN;
2730#ifdef debug
2732#endif
2733 }
2734 //Two particle decay only possible (not frequent event!)
2735 else if(nFlag)
2736 {
2737#ifdef debug
2739#endif
2740 tpd=false;
2741 bQPDG=nQPDG;
2742 rQPDG=NQPDG;
2743 eMass=mNeut;
2744 rMass=GSResNn;
2745 }
2746 else if(pFlag)
2747 {
2748#ifdef debug
2750#endif
2751 tpd=false;
2752 bQPDG=pQPDG;
2753 rQPDG=PQPDG;
2754 eMass=mProt;
2755 rMass=GSResNp;
2756 }
2757 else if(aFlag)
2758 {
2759#ifdef debug
2761#endif
2762 tpd=false;
2763 bQPDG=aQPDG;
2764 rQPDG=AQPDG;
2765 eMass=mAlph;
2766 rMass=GSResNa;
2767 }
2768 else if(lFlag)
2769 {
2770#ifdef debug
2772#endif
2773 tpd=false;
2774 bQPDG=lQPDG;
2775 rQPDG=LQPDG;
2776 eMass=mLamb;
2777 rMass=GSResNl;
2778 }
2779 else
2780 {
2781#ifdef debug
2783#endif
2784 tpd=false;
2785 bQPDG=gQPDG;
2786 rQPDG=GetQPDG();
2787 eMass=0.;
2788 rMass=GSMass;
2789 }
2790 if(tpd)
2791 {
2792 h1mom=G4LorentzVector(0.,0.,0.,eMass);
2793 h2mom=G4LorentzVector(0.,0.,0.,rMass);
2794 h3mom=G4LorentzVector(0.,0.,0.,fMass);
2796 {
2797#ifdef debug
2800#endif
2801 return false;
2802 }
2803 h1mom+=h3mom;
2804 bQPDG=dbQPDG;
2805#ifdef debug
2807 G4double dma=sma-eMass-fMass;
2809 <<rMass<<G4endl;
2810#endif
2811 }
2812 else
2813 {
2814 if(rMass<1600.)
2815 {
2816 if (rQPDG==pQPDG)rMass=mProt;
2817 else if(rQPDG==nQPDG)rMass=mNeut;
2818 else if(rQPDG==lQPDG)rMass=mLamb;
2819 }
2820 h1mom=G4LorentzVector(0.,0.,0.,eMass);
2821 h2mom=G4LorentzVector(0.,0.,0.,rMass);
2823 {
2824#ifdef debug
2827#endif
2828 return false;
2829 }
2830 h1->SetQPDG(bQPDG);
2831 h2->SetQPDG(rQPDG);
2832 h1->Set4Momentum(h1mom);
2833 h2->Set4Momentum(h2mom);
2834#ifdef debug
2837#endif
2838 return true;
2839 }
2840 }
2841 else // Only decay in Baryon+Residual is possible
2842 {
2843#ifdef debug
2845#endif
2846 //@@ Take into account Coulomb Barier Penetration Probability
2847 G4double nLim=0.;
2848 if(nFlag&&mNeut+GSResNn<totMass)
2849 {
2850 G4double ken=totMass-mNeut-GSResNn;
2851 nLim+=N*CoulBarPenProb(0.,ken,0,1)*sqrt(ken);
2852 }
2853 G4double zLim=nLim;
2854 if(pFlag&&mProt+PBarr+GSResNp<totMass)
2855 {
2856 G4double ken=totMass-mProt-GSResNp;
2857 if(barf) ken-=PBarr;
2858 zLim+=Z*CoulBarPenProb(PBarr,ken,1,1)*sqrt(ken);
2859 }
2860 G4double sLim=zLim;
2861 if(lFlag&&mLamb+GSResNl<totMass)
2862 {
2863 G4double ken=totMass-mLamb-GSResNl;
2864 sLim+=S*CoulBarPenProb(0.,ken,0,1)*sqrt(ken);
2865 }
2866 G4double aLim=sLim;
2867 if(aFlag&&mAlph+GSResNa+ABarr<totMass)
2868 {
2869 G4double ken=totMass-mAlph-GSResNa;
2870 if(barf) ken-=ABarr;
2871 aLim+=CoulBarPenProb(ABarr,ken,2,4)*sqrt(ken)*evalph*Z*(Z-1)*N*(N-1)
2872 *6/a1/(a-2)/(a-3);
2873#ifdef debug
2876#endif
2877 }
2879#ifdef debug
2882#endif
2883 if (aFlag&&r>sLim)
2884 {
2885 bQPDG=aQPDG;
2886 eMass=mAlph;
2887 rQPDG=AQPDG;
2888 rMass=GSResNa;
2889#ifdef debug
2891#endif
2892 }
2893 else if(lFlag&&r>zLim&&r<sLim)
2894 {
2895 bQPDG=lQPDG;
2896 eMass=mLamb;
2897 rQPDG=LQPDG;
2898 rMass=GSResNl;
2899#ifdef debug
2901#endif
2902 }
2903 else if(pFlag&&r>nLim&&r<zLim)
2904 {
2905 bQPDG=pQPDG;
2906 eMass=mProt;
2907 rQPDG=PQPDG;
2908 rMass=GSResNp;
2909#ifdef debug
2911#endif
2912 }
2913 else if(nFlag&&r<nLim)
2914 {
2915 bQPDG=nQPDG;
2916 eMass=mNeut;
2917 rQPDG=NQPDG;
2918 rMass=GSResNn;
2919#ifdef debug
2921#endif
2922 }
2923 else if(mProt+GSResNp<totMass)
2924 {
2925#ifdef debug
2927#endif
2928 bQPDG=pQPDG;
2929 rQPDG=PQPDG;
2930 eMass=mProt;
2931 rMass=GSResNp;
2932 }
2933 else if(mAlph+GSResNa<totMass)
2934 {
2935#ifdef debug
2937#endif
2938 bQPDG=aQPDG;
2939 rQPDG=AQPDG;
2940 eMass=mAlph;
2941 rMass=GSResNa;
2942 }
2943 else
2944 {
2945#ifdef debug
2947#endif
2948 bQPDG=gQPDG;
2949 rQPDG=GetQPDG();
2950 eMass=0.;
2951 rMass=GSMass;
2952 }
2953 if(rMass<1600.)
2954 {
2955 if (rQPDG==pQPDG)rMass=mProt;
2956 else if(rQPDG==nQPDG)rMass=mNeut;
2957 else if(rQPDG==lQPDG)rMass=mLamb;
2958 }
2959 h1mom=G4LorentzVector(0.,0.,0.,eMass);
2960 h2mom=G4LorentzVector(0.,0.,0.,rMass);
2961 }
2962 }
2964 {
2965#ifdef debug
2968#endif
2969 return false;
2970 }
2971#ifdef debug
2973#endif
2974 h1->SetQPDG(bQPDG);
2975 h2->SetQPDG(rQPDG);
2976 h1->Set4Momentum(h1mom);
2977 h2->Set4Momentum(h2mom);
2978#ifdef debug
2981#endif
2982 return true;
2983 }
2984 else if(a==1)
2985 {
2986#ifdef debug
2988#endif
2989 h1->SetQPDG(gQPDG);
2991 if (Z>0)
2992 {
2993 h2->SetQPDG(pQPDG);
2995 }
2996 else if(N>0)
2997 {
2998 h2->SetQPDG(nQPDG);
3000 }
3001 else if(S>0)
3002 {
3003 h2->SetQPDG(lQPDG);
3005 }
3006 else return false;
3007 return true;
3008 }
3010 return false;
3011}
G4double CoulBarPenProb(const G4double &CB, const G4double &E, const G4int &C, const G4int &B)
Definition: G4QNucleus.cc:3441
Referenced by EvaporateNucleus().
◆ EvaporateNucleus()
| void G4QNucleus::EvaporateNucleus | ( | G4QHadron * | hA, |
| G4QHadronVector * | oHV | ||
| ) |
@@@@ *** TEMPORARY TO AVOID HYPERMUCLEI FOR GEANT4 *** @@@@
@@ *** ^^^ END OF TEMPORARY ^^^ *** @@
!When kill,DON'T forget to del. theLastQHadron as an instance!!
!When kill,DON'T forget to delete theLastQHadron asAnInstance!!
Definition at line 4171 of file G4QNucleus.cc.
4172{
4190#ifdef pdebug
4192#endif
4194#ifdef pdebug
4196#endif
4199 if(!theBN || thePDG<80000000 || thePDG==90000000) // Hadron, anti-nucleous, or vacuum
4200 {
4201#ifdef debug
4203#endif
4204 if(thePDG==90000000)
4205 {
4206#ifdef qdebug
4208#endif
4209 delete qH;
4211 }
4212 else // Put input to the output (delete equivalent)
4213 {
4215 evaHV->push_back(qH);
4216 }
4217 return;
4218 }
4219 /// @@@@@@@ *** TEMPORARY TO AVOID HYPERMUCLEI FOR GEANT4 *** @@@@@@@
4220 if(thePDG>91000000) //@@MadeForGeant4@@: If there is a Lambda, substitute it by A neutron
4221 {
4223 thePDG-=SSS*999999; // S Neutrons instead of S Lambdas
4224 qH->SetQPDG(G4QPDGCode(thePDG));
4225 theQC = qH->GetQC(); // Quark Content of the hadron
4226#ifdef debug
4228#endif
4229 }
4230 /// @@@ *** ^^^ END OF TEMPORARY ^^^ *** @@@
4231 if(thePDG<80000000)
4232 {
4233#ifdef debug
4235#endif
4236 evaHV->push_back(qH); // TheFundamentalParticles must be FilledAsTheyAre (del.eq)
4237 return;
4238 }
4244#ifdef debug
4247#endif
4248 if(theBN<-2)
4249 {
4251 <<thePDG<<G4endl;
4252 evaHV->push_back(qH);
4253 return;
4254 }
4255 else if(thePDG==91000000||thePDG==90001000||thePDG==90000001) // Excited Lambda* or N*
4256 //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
4257 {
4258 G4double gsM=mNeut;
4259 if(thePDG==90001000) gsM=mProt;
4260 else if(thePDG==91000000) gsM=mLamb;
4261 if(fabs(totMass-gsM)<.001)
4262 {
4263#ifdef debug
4265#endif
4266 evaHV->push_back(qH); // (delete equivalent)
4267 return;
4268 }
4269 else if(totMass<gsM)
4270 {
4271#ifdef qdebug
4273#endif
4274 delete qH;
4275 // G4cerr<<"***G4QN::EvaNuc:Baryon "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
4276 // throw G4QException("G4QNucleus::EvaporateNucleus: Baryon is below Mass Shell");
4277 G4ExceptionDescription ed;
4278 ed << "Baryon is below Mass Shell: Baryon " << thePDG
4281 FatalException, ed);
4282 }
4283 else // Decay in gamma or charged pion (@@ neutral)
4284 {
4285 G4double d=totMass-gsM;
4286#ifdef debug
4288#endif
4289 G4int decPDG=22;
4290 G4double decM=0.;
4291 if(d>142.) // @@ to avoid more precise calculations
4292 {
4293 if(thePDG==90001000) // p* -> n + pi+
4294 {
4295 gsM=mNeut;
4296 thePDG=90000001;
4297 decPDG=211;
4298 decM=mPi;
4299 }
4300 else if(thePDG==90000001) // n* -> p + pi-
4301 {
4302 gsM=mProt;
4303 thePDG=90001000;
4304 decPDG=-211;
4305 decM=mPi;
4306 }
4307 else // decay in Pi0
4308 {
4309 decPDG=111;
4310 decM=mPi0;
4311 }
4312 }
4314 G4LorentzVector g4Mom(0.,0.,0.,decM);
4316 {
4317#ifdef qdebug
4319#endif
4320 delete qH;
4321 // G4cerr<<"***G4QNuc::EvaNuc:h="<<thePDG<<"(GSM="<<gsM<<")+gam>tM="<<totMass<<G4endl;
4322 // throw G4QException("G4QNucleus::EvaporateNucleus:BaryonDecay In Baryon+Gam Error");
4323 G4ExceptionDescription ed;
4324 ed << "BaryonDecay In Baryon+Gam Error: h=" << thePDG << "(GSM="
4327 FatalException, ed);
4328 }
4329#ifdef debug
4331 <<evaHV->size()<<G4endl;
4332#endif
4334#ifdef debug
4336#endif
4337 evaHV->push_back(curH); // Fill Baryon (delete equiv.)
4339#ifdef debug
4341#endif
4342 evaHV->push_back(curG); // Fill gamma/pion (delete equivalent)
4343#ifdef qdebug
4345#endif
4346 delete qH;
4347 }
4348 }
4349 else if(thePDG==89000000||thePDG==89999000||thePDG==89999999) // anti-Lambda* or anti-N*
4350 //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
4351 {
4352 G4double gsM=mNeut;
4353 if(thePDG==89999000) gsM=mProt;
4354 else if(thePDG==89000000) gsM=mLamb;
4355 if(fabs(totMass-gsM)<.001)
4356 {
4357#ifdef debug
4359#endif
4360 evaHV->push_back(qH); // (delete equivalent)
4361 return;
4362 }
4363 else if(totMass<gsM)
4364 {
4365#ifdef qdebug
4367#endif
4368 delete qH;
4369 // G4cerr<<"*G4QN::EvaNuc:antiBaryon="<<thePDG<<" below MassShell, M="<<totMass<<G4endl;
4370 // throw G4QException("G4QNucleus::EvaporateNucleus: anti-Baryon is below Mass Shell");
4371 G4ExceptionDescription ed;
4372 ed << "anti-Baryon is below Mass Shell: antiBaryon=" << thePDG
4375 FatalException, ed);
4376 }
4377 else // Decay in gamma or charged pion (@@ neutral)
4378 {
4379 G4double d=totMass-gsM;
4380#ifdef debug
4382#endif
4383 G4int decPDG=22;
4384 G4double decM=0.;
4385 if(d>142.) // @@ to avoid more precise calculations
4386 {
4387 if(thePDG==89999000) // anti (p* -> n + pi+)
4388 {
4389 gsM=mNeut;
4390 thePDG=89999999;
4391 decPDG=-211;
4392 decM=mPi;
4393 }
4394 else if(thePDG==89999999) // anti (n* -> p + pi-)
4395 {
4396 gsM=mProt;
4397 thePDG=89999000;
4398 decPDG=211;
4399 decM=mPi;
4400 }
4401 else // decay in Pi0
4402 {
4403 decPDG=111;
4404 decM=mPi0;
4405 }
4406 }
4408 G4LorentzVector g4Mom(0.,0.,0.,decM);
4410 {
4411#ifdef qdebug
4413#endif
4414 delete qH;
4415 // G4cerr<<"**G4QNuc::EvaNuc:ah="<<thePDG<<"(GSM="<<gsM<<")+gam>tM="<<totMass<<G4endl;
4416 // throw G4QException("G4QNucleus::EvaporateNucleus:BaryonDecay In Baryon+Gam Error");
4417 G4ExceptionDescription ed;
4418 ed << "BaryonDecay In Baryon+Gam Error: ah=" << thePDG << "(GSM="
4421 FatalException, ed);
4422 }
4423#ifdef debug
4425 <<evaHV->size()<<G4endl;
4426#endif
4428#ifdef debug
4430#endif
4431 evaHV->push_back(curH); // Fill Baryon (delete equiv.)
4433#ifdef debug
4435#endif
4436 evaHV->push_back(curMes); // Fill gamma/pion (delete equivalent)
4437#ifdef qdebug
4439#endif
4440 delete qH;
4441 }
4442 }
4443 else if((thePDG==90001999||thePDG==89999002) && totMass>1080.) // @@ to avoid threshold
4444 //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
4445 {
4446 G4double gsM=mNeut;
4447 G4int barPDG=2112;
4448 G4int mesPDG=-211;
4449 if(thePDG==90001999)
4450 {
4451 gsM=mProt;
4452 barPDG=2212;
4453 mesPDG=211;
4454 }
4455 if(fabs(totMass-gsM-mPi)<.001)
4456 {
4457#ifdef debug
4459#endif
4462 evaHV->push_back(curB); // (delete equivalent)
4463 G4LorentzVector g4Mom=q4M*(mPi/totMass);
4465 evaHV->push_back(curM); // (delete equivalent)
4466#ifdef qdebug
4468#endif
4469 delete qH;
4470 return;
4471 }
4472 else if(totMass<gsM+mPi)
4473 {
4474#ifdef qdebug
4476#endif
4477 delete qH;
4478 // G4cerr<<"***G4QN::EvaNuc:Delta "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
4479 // throw G4QException("G4QNucleus::EvaporateNucleus: Delta is below Mass Shell");
4480 G4ExceptionDescription ed;
4481 ed << "Delta is below Mass Shell: Delta " << thePDG
4484 FatalException, ed);
4485 }
4486 else // Decay in gamma or charged pion (@@ neutral)
4487 {
4489 G4LorentzVector g4Mom(0.,0.,0.,mPi);
4491 {
4492#ifdef qdebug
4494#endif
4495 delete qH;
4496 // G4cerr<<"***G4QNuc::EvaNuc:Dh="<<thePDG<<"N+pi="<<gsM+mPi<<">tM="<<totMass<<G4endl;
4497 // throw G4QException("G4QNucleus::EvaporateNucleus: DeltaDecInBaryon+Pi Error");
4498 G4ExceptionDescription ed;
4499 ed << "DeltaDecInBaryon+Pi Error: Dh=" << thePDG << "N+pi=" << gsM+mPi
4502 FatalException, ed);
4503 }
4504#ifdef debug
4506 <<evaHV->size()<<G4endl;
4507#endif
4509#ifdef debug
4511#endif
4512 evaHV->push_back(curH); // Fill the nucleon (delete equiv.)
4514#ifdef debug
4516#endif
4517 evaHV->push_back(curG); // Fill the pion (delete equivalent)
4518#ifdef qdebug
4520#endif
4521 delete qH;
4522 }
4523 }
4524 else if((thePDG==89998001||thePDG==90000998) && totMass>1080.) // @@ to avoid threshold
4525 //else if(2>3)// One can easily close this decay as it will be done later (time of calc?)
4526 {
4527 G4double gsM=mNeut;
4528 G4int barPDG=-2112;
4529 G4int mesPDG=211;
4530 if(thePDG==89998001)
4531 {
4532 gsM=mProt;
4533 barPDG=-2212;
4534 mesPDG=-211;
4535 }
4536 if(fabs(totMass-gsM-mPi)<.001)
4537 {
4538#ifdef debug
4540#endif
4543 evaHV->push_back(curB); // (delete equivalent)
4544 G4LorentzVector g4Mom=q4M*(mPi/totMass);
4546 evaHV->push_back(curM); // (delete equivalent)
4547#ifdef qdebug
4549#endif
4550 delete qH;
4551 return;
4552 }
4553 else if(totMass<gsM+mPi)
4554 {
4555#ifdef qdebug
4557#endif
4558 delete qH;
4559 // G4cerr<<"***G4QN::EvaNuc:aDelta "<<thePDG<<" is belowMassShell, M="<<totMass<<G4endl;
4560 // throw G4QException("G4QNucleus::EvaporateNucleus: anti-Delta is below Mass Shell");
4561 G4ExceptionDescription ed;
4562 ed << "anti-Delta is below Mass Shell: aDelta " << thePDG
4565 FatalException, ed);
4566 }
4567 else // Decay in gamma or charged pion (@@ neutral)
4568 {
4570 G4LorentzVector g4Mom(0.,0.,0.,mPi);
4572 {
4573#ifdef qdebug
4575#endif
4576 delete qH;
4577 // G4cerr<<"***G4QNuc::EvaNuc:aD="<<thePDG<<"N+pi="<<gsM+mPi<<">tM="<<totMass<<G4endl;
4578 // throw G4QException("G4QNucleus::EvaporateNucleus:AntiDeltaDecayInBaryon+Pi Error");
4579 G4ExceptionDescription ed;
4580 ed << "AntiDeltaDecayInBaryon+Pi Error: aD=" << thePDG << "N+pi="
4583 FatalException, ed);
4584 }
4585#ifdef debug
4588#endif
4590#ifdef debug
4592#endif
4593 evaHV->push_back(curH); // Fill the nucleon (delete equiv.)
4595#ifdef debug
4597#endif
4598 evaHV->push_back(curMes); // Fill the pion (delete equivalent)
4599#ifdef qdebug
4601#endif
4602 delete qH;
4603 }
4604 }
4607 else if((thePDG==89999003 || thePDG==90002999) && totMass>2020.) //=> "ISO-dibarion"
4608 {
4609 // @@ Check that it never comes here !!
4610 G4int nucPDG = 2112;
4611 G4double nucM = mNeut;
4612 G4int piPDG = -211;
4613 if(thePDG==90002999)
4614 {
4615 nucPDG = 2212;
4616 nucM = mProt;
4617 piPDG = 211;
4618 }
4619 if(totMass>mPi+nucM+nucM)
4620 {
4621 G4LorentzVector n14M(0.,0.,0.,nucM);
4622 G4LorentzVector n24M(0.,0.,0.,nucM);
4623 G4LorentzVector pi4M(0.,0.,0.,mPi);
4625 {
4626#ifdef qdebug
4628#endif
4629 delete qH;
4630 // G4cerr<<"***G4QNucleus::EvaporateNucleus: tM="<<totMass<<"-> 2N="<<nucPDG<<"(M="
4631 // <<nucM<<") + pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
4632 // throw G4QException("G4QNucleus::EvaporateNucleus: ISO-Dibaryon DecayIn3 error");
4633 G4ExceptionDescription ed;
4634 ed << " ISO-Dibaryon DecayIn3 error: tM=" << totMass << "-> 2N="
4635 << nucPDG << "(M=" << nucM << ") + pi=" << piPDG << "(M="
4638 FatalException, ed);
4639 }
4640#ifdef qdebug
4642#endif
4643 delete qH;
4645#ifdef debug
4647#endif
4648 evaHV->push_back(h1H); // (delete equivalent)
4650#ifdef debug
4652#endif
4653 evaHV->push_back(h2H); // (delete equivalent)
4655#ifdef debug
4657#endif
4658 evaHV->push_back(piH); // (delete equivalent)
4659 }
4660 else
4661 {
4662#ifdef qdebug
4664#endif
4665 delete qH;
4666 // G4cerr<<"***G4QNucleus::EvapNucleus: IdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
4667 // <<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
4668 // throw G4QException("G4QNucleus::EvaporateNucleus:ISO-DiBaryon is under MassShell");
4669 G4ExceptionDescription ed;
4670 ed << "ISO-DiBaryon is under MassShell: IdPDG=" << thePDG << ", q4M="
4671 << q4M << ", M=" << totMass << " < M_2N+Pi, d=" << totMass-2*nucM-mPi
4672 << G4endl;
4674 FatalException, ed);
4675 }
4676 }
4677 else if((thePDG==90000002||thePDG==90001001||thePDG==90002000)&&totMass>2020.) //=> IsoBP
4678 {
4679 // @@ Pi0 can be taken into account !
4680 G4int n1PDG = 2212;
4681 G4int n2PDG = 2112;
4682 G4int piPDG = -211;
4683 G4double n1M = mProt;
4684 G4double n2M = mNeut;
4685 if (thePDG==90002000) piPDG = 211; // pp -> np + pi-
4686 else if(thePDG==90000002) piPDG = -211; // nn -> np + pi-
4687 else // np -> 50%(nnpi+) 50%(pppi-)
4688 {
4690 {
4691 n1PDG=2112;
4692 n1M=mNeut;
4693 piPDG = 211;
4694 }
4695 else
4696 {
4697 n2PDG=2212;
4698 n2M=mProt;
4699 piPDG = -211;
4700 }
4701 }
4702 if(totMass>mPi+n1M+n2M)
4703 {
4704 G4LorentzVector n14M(0.,0.,0.,n1M);
4705 G4LorentzVector n24M(0.,0.,0.,n2M);
4706 G4LorentzVector pi4M(0.,0.,0.,mPi);
4708 {
4709 // G4cerr<<"***G4QNucl::EvapNucleus:IsoDN, tM="<<totMass<<"-> N1="<<n1PDG<<"(M="<<n1M
4710 // <<") + N2="<<n2PDG<<"(M="<<n2M<<") + pi="<<piPDG<<" (Mpi="<<mPi<<")"<<G4endl;
4711 // throw G4QException("G4QNucl::EvaporateNucleus:ISO-dibaryon excit. DecayIn3 error");
4712 G4ExceptionDescription ed;
4713 ed << "ISO-dibaryon excit. DecayIn3 error: IsoDN, tM=" << totMass
4714 << "-> N1=" << n1PDG << "(M=" << n1M << ") + N2=" << n2PDG
4715 << "(M=" << n2M << ") + pi=" << piPDG << " (Mpi=" << mPi << ")"
4716 << G4endl;
4718 FatalException, ed);
4719 }
4720#ifdef qdebug
4722#endif
4723 delete qH;
4725#ifdef debug
4727#endif
4728 evaHV->push_back(h1H); // (delete equivalent)
4730#ifdef debug
4732#endif
4733 evaHV->push_back(h2H); // (delete equivalent)
4735#ifdef debug
4737#endif
4738 evaHV->push_back(piH); // (delete equivalent)
4739 }
4740 else
4741 {
4742#ifdef qdebug
4744#endif
4745 delete qH;
4746 // G4cerr<<"***G4QNuc::EvaporateNucleus: IdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
4747 // <<" < M1+M2+Pi, d="<<totMass-n1M-n2M-mPi<<G4endl;
4748 // throw G4QException("G4QNucleus::EvaporateNucleus:IsoDiBarState is under MassShell");
4749 G4ExceptionDescription ed;
4750 ed << "IsoDiBarState is under MassShell: IdPDG=" << thePDG << ", q4M="
4751 << q4M << ", M=" << totMass << " < M1+M2+Pi, d="
4752 << totMass-n1M-n2M-mPi << G4endl;
4754 FatalException, ed);
4755 }
4756 }
4758 else if((thePDG==90000997 || thePDG==89997001) && totMass>2020.) //=> "anti-ISO-dibarion"
4759 {
4760 // @@ Check that it never comes here !!
4761 G4int nucPDG = -2112;
4762 G4double nucM = mNeut;
4763 G4int piPDG = 211;
4764 if(thePDG==90002999)
4765 {
4766 nucPDG = -2212;
4767 nucM = mProt;
4768 piPDG = -211;
4769 }
4770 if(totMass>mPi+nucM+nucM)
4771 {
4772 G4LorentzVector n14M(0.,0.,0.,nucM);
4773 G4LorentzVector n24M(0.,0.,0.,nucM);
4774 G4LorentzVector pi4M(0.,0.,0.,mPi);
4776 {
4777#ifdef qdebug
4779#endif
4780 delete qH;
4781 // G4cerr<<"***G4QNucleus::EvaporateNucleus:antiM="<<totMass<<"-> 2aN="<<nucPDG<<"(M="
4782 // <<nucM<<") + pi="<<piPDG<<"(M="<<mPi<<")"<<G4endl;
4783 // throw G4QException("G4QNucleus::EvaporateNucleus:Anti-ISO-DibaryonDecayIn3 error");
4784 G4ExceptionDescription ed;
4785 ed << "Anti-ISO-DibaryonDecayIn3 error: antiM=" << totMass
4786 << "-> 2aN=" << nucPDG << "(M=" << nucM << ") + pi=" << piPDG
4789 FatalException, ed);
4790 }
4791#ifdef qdebug
4793#endif
4794 delete qH;
4796#ifdef debug
4798#endif
4799 evaHV->push_back(h1H); // (delete equivalent)
4801#ifdef debug
4803#endif
4804 evaHV->push_back(h2H); // (delete equivalent)
4806#ifdef debug
4808#endif
4809 evaHV->push_back(piH); // (delete equivalent)
4810 }
4811 else
4812 {
4813#ifdef qdebug
4815#endif
4816 delete qH;
4817 // G4cerr<<"***G4QNucleus::EvapNucleus: aIdPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
4818 // <<" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<G4endl;
4819 // throw G4QException("G4QNucleus::EvaporateNucleus:AntiISODiBaryon is underMassShell");
4820 G4ExceptionDescription ed;
4821 ed << "AntiISODiBaryon is underMassShell: aIdPDG=" << thePDG << ", q4M="
4822 << q4M << ", M=" << totMass << " < M_2N+Pi, d=" << totMass-2*nucM-mPi
4823 << G4endl;
4825 FatalException, ed);
4826 }
4827 }
4828 else if((thePDG==89999998||thePDG==89998999||thePDG==89998000)&&totMass>2020.)//=>AnIsoBP
4829 {
4830 // @@ Pi0 can be taken into account !
4831 G4int n1PDG = -2212;
4832 G4int n2PDG = -2112;
4834 G4double n1M = mProt;
4835 G4double n2M = mNeut;
4836 if (thePDG==89998000) piPDG = -211; // anti ( pp -> np + pi- )
4837 else if(thePDG==89999998) piPDG = 211; // anti ( nn -> np + pi- )
4838 else // anti ( np -> 50%(nnpi+) 50%(pppi-) )
4839 {
4841 {
4842 n1PDG=-2112;
4843 n1M=mNeut;
4844 piPDG = -211;
4845 }
4846 else
4847 {
4848 n2PDG=-2212;
4849 n2M=mProt;
4850 piPDG = 211;
4851 }
4852 }
4853 if(totMass>mPi+n1M+n2M)
4854 {
4855 G4LorentzVector n14M(0.,0.,0.,n1M);
4856 G4LorentzVector n24M(0.,0.,0.,n2M);
4857 G4LorentzVector pi4M(0.,0.,0.,mPi);
4859 {
4860 // G4cerr<<"**G4QNucl::EvapNucleus:IsoDN,antiM="<<totMass<<"-> N1="<<n1PDG<<"(M="<<n1M
4861 // <<") + N2="<<n2PDG<<"(M="<<n2M<<") + pi="<<piPDG<<" (Mpi="<<mPi<<")"<<G4endl;
4862 // throw G4QException("G4QNucl::EvaporateNucleus:AntiExcitedDibaryon DecayIn3 error");
4863 G4ExceptionDescription ed;
4864 ed << "AntiExcitedDibaryon DecayIn3 error: IsoDN,antiM=" << totMass
4865 << "-> N1=" << n1PDG << "(M=" << n1M << ") + N2=" << n2PDG << "(M="
4868 FatalException, ed);
4869 }
4870#ifdef qdebug
4872#endif
4873 delete qH;
4875#ifdef debug
4877#endif
4878 evaHV->push_back(h1H); // (delete equivalent)
4880#ifdef debug
4882#endif
4883 evaHV->push_back(h2H); // (delete equivalent)
4885#ifdef debug
4887#endif
4888 evaHV->push_back(piH); // (delete equivalent)
4889 }
4890 else
4891 {
4892#ifdef qdebug
4894#endif
4895 delete qH;
4896 // G4cerr<<"***G4QNuc::EvaporateNucleus:andPDG="<<thePDG<<", q4M="<<q4M<<", M="<<totMass
4897 // <<" < M1+M2+Pi, d="<<totMass-n1M-n2M-mPi<<G4endl;
4898 // throw G4QException("G4QNucleus::EvaporateNucleus:AntiDiBarState is under MassShell");
4899 G4ExceptionDescription ed;
4900 ed << "AntiDiBarState is under MassShell: andPDG=" << thePDG << ", q4M="
4901 << q4M << ", M=" << totMass << " < M1+M2+Pi, d="
4902 << totMass-n1M-n2M-mPi << G4endl;
4904 FatalException, ed);
4905 }
4906 }
4908 else if(fabs(totMass-totGSM)<.001) // Fill as it is or decay Be8, He5, Li5 (@@ add more)
4909 {
4910#ifdef debug
4912#endif
4913 if(thePDG==90004004)
4914 {
4915 DecayAlphaAlpha(qH,evaHV);
4916 } // "Alpha+Alpha Decay" case (del eq.)
4917 else if(thePDG==90004002)
4918 {
4919 DecayAlphaDiN(qH,evaHV);
4920 } // Decay alpha+2p(alpha+2n is stable)
4921 else if((theC==theBN||theN==theBN||theS==theBN)&&theBN>1)
4922 {
4923 DecayMultyBaryon(qH,evaHV);
4924 }
4925 else if(theBN==5)
4926 {
4927 DecayAlphaBar(qH,evaHV);
4928 } // Decay unstable A5 system (del eq.)
4929 else
4930 {
4931 evaHV->push_back(qH);
4932 } // Fill as it is (del eq.)
4933 }
4934 else if(theBN>1 && thePDG>88000000 && thePDG<89000000) //==> 2antiK in the nucleus
4935 {
4938 G4int bN=theBN-bZ;
4939 G4int k1PDG = 321;
4940 G4double mK1= mK;
4941 G4int k2PDG = 321;
4942 G4double mK2= mK;
4943 G4int nucPDG = thePDG;
4944 if(bZ>=bN) nucPDG+=999000;
4945 else
4946 {
4947 nucPDG+=999999;
4948 k1PDG = 311;
4949 mK1= mK0;
4950 }
4951 if(bZ>bN) nucPDG+=999000;
4952 else
4953 {
4954 nucPDG+=999999;
4955 k2PDG = 311;
4956 mK2= mK0;
4957 }
4960 G4LorentzVector n4M(0.,0.,0.,nucM);
4961 G4LorentzVector k14M(0.,0.,0.,mK1);
4962 G4LorentzVector k24M(0.,0.,0.,mK2);
4964 {
4965#ifdef qdebug
4967#endif
4968 delete qH;
4969 // G4cout<<"***G4QNucleus::EvaNuc:tM="<<totMass<<"-> N="<<nucPDG<<"(M="<<nucM<<") + k1="
4970 // <<k1PDG<<"(M="<<mK1<<") + k2="<<k2PDG<<"(M="<<mK2<<")"<<G4endl;
4971 // throw G4QException("G4QNucleus::EvaporateNucleus: Nucleus+2antiK DecayIn3 error");
4972 G4ExceptionDescription ed;
4973 ed << " Nucleus+2antiK DecayIn3 error: tM=" << totMass << "-> N="
4974 << nucPDG << "(M=" << nucM << ") + k1=" << k1PDG << "(M=" << mK1
4977 FatalException, ed);
4978 }
4979#ifdef qdebug
4981#endif
4982 delete qH;
4984#ifdef debug
4986#endif
4987 evaHV->push_back(k1H); // (delete equivalent)
4989#ifdef debug
4991#endif
4992 evaHV->push_back(k2H); // (delete equivalent)
4994#ifdef debug
4996#endif
4997 evaHV->push_back(nH); // (delete equivalent)
4998 }
4999 // ***>> From here the EVA code starts (baryons/hyperons can be excited) <<***
5000 else if ( (thePDG > 80000000 && thePDG != 90000000) ||
5001 thePDG == 2112 || thePDG == 2212 || thePDG == 3122)
5002 { // @@ Improve for Sigma+, Sigma-, Ksi0 & Ksi- content in the Total Np/Nn Nuclei
5003 if(thePDG<80000000) // Switch from QHadronCode to QNuclearCode
5004 {
5005 if (thePDG==2112) thePDG=90000001; // n
5006 else if(thePDG==2212) thePDG=90001000; // p
5007 else if(thePDG==3122) thePDG=91000000; // lambda
5008 }
5015#ifdef ppdebug
5017#endif
5019#ifdef debug
5024 G4double dM=totMass-GSMass;
5026 ////////if(dM>7.) throw G4QException("G4QNucleus::EvaporateNucleus: Big Excitation");
5027#endif
5030#ifdef debug
5032#endif
5033 if(fabs(totMass-GSMass)<.003&&!bsCond&&!dbsCond) // GS or can't split 1(2)B FillAsItIs
5034 {
5035#ifdef debug
5037#endif
5038 evaHV->push_back(qH);
5039 return;
5040 }
5041 else if ( ( bA == 1 || (!bsCond && !dbsCond) ) && totMass > GSMass+.003 )//=>Fuse&Decay
5042 //else if(2>3) // Close "Fuse&Decay Technology"***@@@***
5043 {
5044#ifdef debug
5046#endif
5049 while(nOfOUT) // Try BackFusionDecays till something is in
5050 {
5051 G4QHadron* theLast = (*evaHV)[nOfOUT-1];
5054 //////////////////G4int gam = theLast->GetPDGCode();
5055#ifdef debug
5057#endif
5058 while(nFrag) // => "Delete Last Decayed Hadron" case
5059 {
5060 G4QHadron* thePrev = (*evaHV)[nOfOUT-2];
5061 nFrag = thePrev->GetNFragments();
5063#ifdef debug
5066#endif
5067 evaHV->pop_back(); // the prev QHadron is excluded from OUTPUT
5068 delete theLast;//!!When kill,DON'T forget to del. theLastQHadron as an instance!!
5069 theLast = thePrev; // Update theLastPntr(Prev instead of Last)
5070 lastBN=prevBN;
5071 nOfOUT--; // Reduce the stack for the Last decayed hadron
5072 }
5073 if(nOfOUT)
5074 {
5075 if(lastBN<1&&nOfOUT>1) // Try Once To Skip Mesons/Gammas & Antibaryons
5076 {
5078 evaHV->pop_back(); // the last QHadron is excluded from OUTPUT
5079 evaHV->pop_back(); // the prev QHadron is excluded from OUTPUT
5080 evaHV->push_back(theLast); // the Last becomes the Prev (1st part of exch)
5081 evaHV->push_back(thePrev); // the Prev becomes the Last (2nd part of exch)
5082 theLast = thePrev; // Update theLastPointer (Prev instead of Last)
5083 }
5087 totQC+=lastQC; // Update (increase) the total QC (prototype)
5088 q4M+=last4M; // Update (increase) the total 4-momentum
5093 G4double dM=totMass-GSMass-lastM;
5094#ifdef debug
5096 <<dM<<G4endl;
5097#endif
5098 if(dM>-0.001) // Decay in the qH and the last is impossible
5099 {
5101 if(dM<=0.)
5102 {
5103 evaHV->pop_back(); // lastQHadron is excluded from QHadrV asIs in TRN
5104 delete theLast; //When kill, DON'T forget to delete lastQHadron asAnInstance!
5105#ifdef qdebug
5107#endif
5108 delete qH;
5109#ifdef debug
5111#endif
5113 else evaHV->push_back(evH);// Fill TRN to HVect asIs (delete equivalent)
5114 }
5115 else // Decay TotalResidualNucleus in GSM+Last
5116 {
5117 G4LorentzVector r4Mom(0.,0.,0.,GSMass);
5119 {
5120 evaHV->pop_back(); // lastQHadron is excluded from QHadrV as is in TRN
5121 delete theLast; //When kill,DON'T forget to delete lastQHadron asAnInstance
5122#ifdef qdebug
5124#endif
5125 delete qH;
5126#ifdef debug
5128#endif
5129 evaHV->push_back(evH);// Fill TRN to Vect as it is (delete equivalent)
5130#ifdef debug
5132#endif
5133 }
5134 else // Decay in GSM+theLast succeeded
5135 {
5136 delete evH;
5137#ifdef qdebug
5139#endif
5140 delete qH;
5143#ifdef debug
5145#endif
5146 // @@ What about others, not DB possibilities?
5147 if(thePDG==92000000||thePDG==90002000||thePDG==90000002)
5149 else evaHV->push_back(nucH);// Fill the Residual Nucleus (del.eq.)
5150 }
5151 }
5152 bnfound=false;
5153 break;
5154 }
5155 thePDG=totPDG; // Make ResidualNucleus outOf theTotResidualNucl
5157 evaHV->pop_back(); // the last QHadron is excluded from OUTPUT
5158 delete theLast;//!!When kill,DON'T forget to delete theLastQHadron asAnInstance!!
5159 nOfOUT--; // Update the value of OUTPUT entries
5160 }
5161 }
5162 if(bnfound)
5163 {
5165 G4LorentzVector g4Mom(0.,0.,0.,0.);
5167 {
5168#ifdef qdebug
5170#endif
5171 delete qH;
5172 // G4cerr<<"**G4QN::EvaNuc:h="<<thePDG<<"(GSM="<<GSMass<<")+g>tM="<<totMass<<G4endl;
5173 // throw G4QException("G4QNucleus::EvaporateNucleus: Decay in Gamma failed");
5174 G4ExceptionDescription ed;
5175 ed << " Decay in Gamma failed: h=" << thePDG << "(GSM=" << GSMass
5178 FatalException, ed);
5179 }
5180#ifdef debug
5182#endif
5184#ifdef debug
5186#endif
5188 else evaHV->push_back(curH); // Fill the TotalResidualNucleus (del.equiv.)
5190#ifdef debug
5192#endif
5193 evaHV->push_back(curG); // Fill the gamma (delete equivalent)
5194#ifdef qdebug
5196#endif
5197 delete qH;
5198 }
5199 }
5203 else if(totMass<GSMass+.003&&(bsCond||dbsCond))//==>" M<GSM but decay is possible" case
5204 {
5205#ifdef pdebug
5207#endif
5208 //G4double gResM =1000000.; // Prototype of mass of residual for a gamma
5210 if(bN==4&&bZ==2&&!bS) // It's He6 nucleus
5211 {
5212 gResPDG= thePDG; // PDG of the Residual Nucleus
5213 //gResM = mHel6; // min mass of the Residual Nucleus
5214 }
5217 if(bsCond==2112&&bN>0&&bA>1) // There's aNeutr in theNucl, which can be split
5218 {
5219 G4QContent resQC=totQC-neutQC;
5221 nResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5222 if (nResPDG==90000001) nResM=mNeut;
5223 else if(nResPDG==90001000) nResM=mProt;
5224 else if(nResPDG==91000000) nResM=mLamb;
5225 else nResM=resN.GetMZNS(); // min mass of the Residual Nucleus
5226 }
5229 if(bsCond==2212&&bZ>0&&bA>1) // There's aProton in Nucl, which can be split
5230 {
5231 G4QContent resQC=totQC-protQC;
5233 pResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5234 if (pResPDG==90000001) pResM=mNeut;
5235 else if(pResPDG==90001000) pResM=mProt;
5236 else if(pResPDG==91000000) pResM=mLamb;
5237 else pResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5238 }
5241 if(bsCond==3122&&bS>0&&bA>1) // There's aLambd in theNucl, which can be split
5242 {
5243 G4QContent resQC=totQC-lambQC;
5245 lResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5246 if (lResPDG==90000001) lResM=mNeut;
5247 else if(lResPDG==90001000) lResM=mProt;
5248 else if(lResPDG==91000000) lResM=mLamb;
5249 else lResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5250 }
5253 if(bsCond==90001001&&bN>0&&bZ>0&&bA>2)// There's aDeuter in Nucl, which canBeRadiated
5254 {
5255 G4QContent resQC=totQC-deutQC;
5257 dResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5258 if (dResPDG==90000001) dResM=mNeut;
5259 else if(dResPDG==90001000) dResM=mProt;
5260 else if(dResPDG==91000000) dResM=mLamb;
5261 else dResM =resN.GetMZNS(); // minMass of the Residual Nucleus
5262 }
5265 if(bsCond==90002002&&bN>1&&bZ>1&&bA>4)// There's Alpha in theNucl, which can be split
5266 {
5267 G4QContent resQC=totQC-alphQC;
5269 aResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5270 if (aResPDG==90000001) aResM=mNeut;
5271 else if(aResPDG==90001000) aResM=mProt;
5272 else if(aResPDG==91000000) aResM=mLamb;
5273 else aResM =resN.GetMZNS(); // minMass of the Residual Nucleus
5274 }
5277 if(dbsCond&&bN>1&&bA>2) // It's nucleus and there is a dineutron
5278 {
5279 G4QContent resQC=totQC-neutQC-neutQC;
5281 nnResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5282 if (nnResPDG==90000001) nnResM=mNeut;
5283 else if(nnResPDG==90001000) nnResM=mProt;
5284 else if(nnResPDG==91000000) nnResM=mLamb;
5285 else nnResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5286 }
5289 if(dbsCond&&bZ>1&&bA>2) // It's nucleus and there is a diproton
5290 {
5291 G4QContent resQC=totQC-protQC-protQC;
5293 ppResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5294 if (ppResPDG==90000001) ppResM=mNeut;
5295 else if(ppResPDG==90001000) ppResM=mProt;
5296 else if(ppResPDG==91000000) ppResM=mLamb;
5297 else ppResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5298 }
5301 if(dbsCond&&bN>0&&bZ>0&&bA>2) // It's nucleus and there is aProton & aNeutron
5302 {
5303 G4QContent resQC=totQC-neutQC-protQC;
5305 npResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5306 if (npResPDG==90000001) npResM=mNeut;
5307 else if(npResPDG==90001000) npResM=mProt;
5308 else if(npResPDG==91000000) npResM=mLamb;
5309 else npResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5310 }
5313 if(dbsCond&&bN>0&&bS>0&&bA>2) // It's nucleus and there is aLambda & aNeutron
5314 {
5315 G4QContent resQC=totQC-lambQC-protQC;
5317 lnResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5318 if (lnResPDG==90000001) lnResM=mNeut;
5319 else if(lnResPDG==90001000) lnResM=mProt;
5320 else if(lnResPDG==91000000) lnResM=mLamb;
5321 else lnResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5322 }
5325 if(dbsCond&&bS>0&&bZ>0&&bA>2) // It's nucleus and there is aProton and aLambda
5326 {
5327 G4QContent resQC=totQC-neutQC-protQC;
5329 lpResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5330 if (lpResPDG==90000001) lpResM=mNeut;
5331 else if(lpResPDG==90001000) lpResM=mProt;
5332 else if(lpResPDG==91000000) lpResM=mLamb;
5333 else lpResM =resN.GetMZNS(); // minMass of the Residual Nucleus
5334 }
5337 if(dbsCond&&bS>1&&bA>2) // It's nucleus and there is a di-lambda
5338 {
5339 G4QContent resQC=totQC-neutQC-protQC;
5341 llResPDG=resN.GetPDG(); // PDG of the Residual Nucleus
5342 if (llResPDG==90000001) llResM=mNeut;
5343 else if(llResPDG==90001000) llResM=mProt;
5344 else if(llResPDG==91000000) llResM=mLamb;
5345 else llResM =resN.GetMZNS(); // min mass of the Residual Nucleus
5346 }
5347#ifdef debug
5349 <<bN<<",Z="<<bZ<<",nL="<<bS<<",totM="<<totMass<<",n="<<totMass-nResM-mNeut
5351#endif
5352 if ( thePDG == 90004004 ||
5353 (thePDG == 90002004 && totMass > mHel6+.003) ||
5354 (bA > 4 && bsCond && bN > 1 && bZ > 1 && totMass > aResM+mAlph) ||
5355 (bA > 1 && bsCond && ( (bN > 0 && totMass > nResM+mNeut) ||
5356 (bZ > 0 && totMass > pResM+mProt) ||
5357 (bS > 0 && totMass > lResM+mLamb) ) ) ||
5358 (bA > 2 &&
5359 (( bN > 0 && bZ > 0 &&
5360 ( (bsCond && totMass > dResM+mDeut) || (dbsCond && totMass > dResM+mDeut) )
5361 ) || ( dbsCond && ( (bN > 1 && totMass > nnResM+mNeut+mNeut) ||
5362 (bZ > 1 && totMass > ppResM+mProt+mProt) ||
5363 (bS > 1 && totMass > llResM+mLamb+mLamb) ||
5364 (bN && bS && totMass > lnResM+mLamb+mNeut) ||
5365 (bZ && bS && totMass > lpResM+mLamb+mProt)
5366 )
5367 )
5368 )
5369 )
5370 )
5371 {
5373 G4int resPDG = 90002002;
5374 G4int thdPDG = 0;
5375 G4double barM= mAlph;
5376 G4double resM= mAlph;
5378 G4double tMC=totMass+.0002;
5379 if(gResPDG&&tMC>mHel6+.003) // Can make radiative decay of He6 (priority 0)
5380 {
5381 barPDG=90002004;
5382 resPDG=22;
5383 barM =mHel6;
5384 resM =0.;
5385 }
5386 else if(nResPDG&&tMC>nResM+mNeut) // Can radiate a neutron (priority 1)
5387 {
5388 barPDG=90000001;
5389 resPDG=nResPDG;
5390 barM =mNeut;
5391 resM =nResM;
5392 }
5393 else if(pResPDG&&totMass+.001>pResM+mProt) // Can radiate a proton (priority 2)
5394 {
5395 barPDG=90001000;
5396 resPDG=pResPDG;
5397 barM =mProt;
5398 resM =pResM;
5399 }
5400 else if(lResPDG&&tMC>lResM+mLamb) // Can radiate a Lambda (priority 3) @@ Sigma0
5401 {
5402 barPDG=91000000;
5403 resPDG=lResPDG;
5404 barM =mLamb;
5405 resM =lResM;
5406 }
5407 else if(thePDG!=90004004&&bN>1&&bZ>1&&bA>4&&tMC>aResM+mAlph)// Decay in alpha (p4)
5408 {
5409 barPDG=90002002;
5410 resPDG=aResPDG;
5411 barM =mAlph;
5412 resM =aResM;
5413 }
5414 else if(dResPDG&&tMC>dResM+mDeut) // Can radiate a Deuteron (priority 5)
5415 {
5416 barPDG=90001001;
5417 resPDG=dResPDG;
5418 barM =mDeut;
5419 resM =dResM;
5420 }
5421 else if(ppResPDG&&tMC>ppResM+mProt+mProt)// Can radiate a DiProton (priority 6)
5422 {
5423 barPDG=90001000;
5424 resPDG=ppResPDG;
5425 thdPDG=90001000;
5426 barM =mProt;
5427 resM =ppResM;
5428 thdM =mProt;
5429 }
5430 else if(nnResPDG&&tMC>nnResM+mNeut+mNeut)// Can radiate a DiNeutron (priority 7)
5431 {
5432 barPDG=90000001;
5433 resPDG=nnResPDG;
5434 thdPDG=90000001;
5435 barM =mNeut;
5436 resM =nnResM;
5437 }
5438 else if(npResPDG&&tMC>npResM+mProt+mNeut)// Can radiate a neutron+proton (prior 8)
5439 {
5440 barPDG=90001000;
5441 resPDG=npResPDG;
5442 thdPDG=90000001;
5443 barM =mProt;
5444 resM =npResM;
5445 }
5446 else if(lnResPDG&&tMC>lnResM+mLamb+mNeut)// Can radiate a Lambda+neutron (prior 9)
5447 {
5448 barPDG=91000000; // @@ Sigma0
5449 resPDG=lnResPDG;
5450 thdPDG=90000001;
5451 barM =mLamb;
5452 resM =lnResM;
5453 }
5454 else if(lpResPDG&&tMC>lpResM+mLamb+mProt)// Can radiate a Lambda+proton (prior 10)
5455 {
5456 barPDG=91000000; // @@ Sigma0
5457 resPDG=lpResPDG;
5458 thdPDG=90001000;
5459 barM =mLamb;
5460 resM =lpResM;
5461 thdM =mProt;
5462 }
5463 else if(llResPDG&&tMC>llResM+mLamb+mLamb)// Can radiate a DiLambda (priority 11)
5464 {
5465 barPDG=91000000; // @@ Sigma0
5466 resPDG=llResPDG;
5467 thdPDG=91000000; // @@ Sigma0
5468 barM =mLamb;
5469 resM =llResM;
5470 thdM =mLamb;
5471 }
5472 else if(thePDG!=90004004&&tMC>GSMass)// If it's not Be8 decay in gamma & GSM
5473 {
5474 barPDG=thePDG;
5475 resPDG=22;
5476 barM =GSMass;
5477 resM =0.;
5478 }
5479 else if(thePDG!=90004004)
5480 {
5481 // G4cerr<<"***G4QNuc::EvaN:PDG="<<thePDG<<",M="<<totMass<<"< GSM="<<GSMass<<G4endl;
5482 // throw G4QException("G4QNucleus::EvaporateNucleus: M<GSM & can't decayInPNL");
5483 G4ExceptionDescription ed;
5484 ed << "M<GSM & can't decayInPNL: PDG=" << thePDG << ",M=" << totMass
5487 FatalException, ed);
5488 }
5489 G4LorentzVector a4Mom(0.,0.,0.,barM);
5490 G4LorentzVector b4Mom(0.,0.,0.,resM);
5491 if(!thdPDG)
5492 {
5493 if(!qH->DecayIn2(a4Mom,b4Mom))
5494 {
5495 evaHV->push_back(qH); // Fill as it is (delete equivalent)
5497 <<lResPDG<<",N="<<bN<<",Z="<<bZ<<",L="<<bS<<",totM="<<totMass<<",n="
5498 <<totMass-nResM-mNeut<<",p="<<totMass-pResM-mProt<<",l="
5499 <<totMass-lResM-mLamb<<G4endl;
5501 return;
5502 }
5503 else
5504 {
5505#ifdef qdebug
5507#endif
5508 delete qH;
5510#ifdef debug
5512#endif
5513 evaHV->push_back(HadrB); // Fill the baryon (delete equivalent)
5515#ifdef debug
5517#endif
5518 // @@ Self-call !!
5520 else evaHV->push_back(HadrR); // Fill ResidNucl=Baryon to OutHadronVector
5521 }
5522 }
5523 else
5524 {
5525 G4LorentzVector c4Mom(0.,0.,0.,thdM);
5526 if(!qH->DecayIn3(a4Mom,b4Mom,c4Mom))
5527 {
5528 evaHV->push_back(qH); // Fill as it is (delete equivalent)
5530 <<ppResPDG<<",N="<<bN<<",Z="<<bZ<<",L="<<bS<<",tM="<<totMass<<",nn="
5531 <<totMass-nnResM-mNeut-mNeut<<",np="<<totMass-npResM-mProt-mNeut<<",pp="
5532 <<totMass-ppResM-mProt-mProt<<G4endl;
5534 return;
5535 }
5536 else
5537 {
5538#ifdef qdebug
5540#endif
5541 delete qH;
5543#ifdef debug
5545#endif
5546 evaHV->push_back(HadrB); // Fill the first baryon (del.equiv.)
5548#ifdef debug
5550#endif
5551 evaHV->push_back(HadrC); // Fill the second baryon (del.equiv.)
5553#ifdef debug
5555#endif
5556 // @@ Self-call !!
5558 else evaHV->push_back(HadrR); // Fill ResidNucl=Baryon to OutputHadrVector
5559 }
5560 }
5561 }
5562 else if (fabs(totMass-GSMass)<.003) // @@ Looks like a duplication of the prev. check
5563 {
5564#ifdef debug
5566#endif
5567 evaHV->push_back(qH); // FillAsItIs (del.eq.)
5568 return;
5569 }
5570 else // "System is below mass shell and can't decay" case
5571 {
5572#ifdef debug
5575#endif
5576 evaHV->push_back(qH); // Correct or fill as it is
5577 return;
5578 }
5579 }
5580 else // ===> Evaporation of the excited system
5581 {
5582#ifdef pdebug
5585#endif
5586 G4LorentzVector b4M;
5587 G4LorentzVector r4M;
5589 G4int bPDG=0;
5590 G4int rPDG=0;
5591 //G4double bM = 0.; // Prototype of Real Mass of the EvaporatedDibaryon
5597 G4int evcn=0;
5598 //G4int evcm=27;
5600 // @@ Does not look like it can make two attempts @@ Improve, Simplify @@
5601 while(evC&&evcn<evcm)
5602 {
5603 evC=true;
5604 evcn++;
5605 if(!qNuc.EvaporateBaryon(bHadron,rHadron)) // Evaporation did not succeed
5606 {
5607#ifdef debug
5609#endif
5610 delete bHadron;
5611 delete rHadron;
5612#ifdef debug
5614#endif
5615 evaHV->push_back(qH); // fill AsItIs
5616 return;
5617 }
5618 evC=false;
5619 b4M=bHadron->Get4Momentum();
5620 r4M=rHadron->Get4Momentum();
5621 //bM = b4M.m(); // Real mass of the evaporated dibaryon
5625 bPDG=bHadron->GetPDGCode();
5626 rPDG=rHadron->GetPDGCode();
5627#ifdef debug
5629 //G4int rC=rHadron->GetCharge(); // Baryon number of the residual nucleus
5630 G4double bCB=qNuc.CoulombBarrier(bC,bB);
5631 //G4double rCB=qNuc.CoulombBarrier(rC,rB);
5634#endif
5635 //if(b4M.e()-b4M.m()<bCB&&evcn<evcm) evC=true;
5636 } // End of while
5637#ifdef debug
5640#endif
5641#ifdef qdebug
5643#endif
5644 delete qH;
5645 if(bB<2) evaHV->push_back(bHadron); // Fill EvaporatedBaryon (del.equivalent)
5647 else if(bB==4) evaHV->push_back(bHadron); // "Alpha radiation" case (del.eq.)
5651 else
5652 {
5653 delete bHadron;
5654 // G4cerr<<"***G4QNuc::EvaNuc:bB="<<bB<<">2 - unexpected evaporated fragment"<<G4endl;
5655 // throw G4QException("G4QNucleus::EvaporateNucleus: Wrong evaporation act");
5656 G4ExceptionDescription ed;
5657 ed << "Wrong evaporation act: EvaNuc:bB=" << bB
5660 FatalException, ed);
5661 }
5663 else if(rB==2) // => "Dibaryon" case needs decay @@ DecayDibaryon
5664 {
5666#ifdef debug
5668#endif
5669 if(rM<=rGSM-0.01)
5670 {
5671 delete rHadron;
5672 // G4cerr<<"***G4QNucleus::EvaporNucl: <residual> M="<<rM<<" < GSM="<<rGSM<<G4endl;
5673 // throw G4QException("G4QNucleus::EvaporateNucleus: Evaporation below MassShell");
5674 G4ExceptionDescription ed;
5675 ed << "Evaporation below MassShell: <residual> M=" << rM << " < GSM="
5676 << rGSM << G4endl;
5678 FatalException, ed);
5679 }
5680 else if(fabs(rM-rGSM)<0.01&&rPDG==90001001) evaHV->push_back(rHadron); // (DE)
5682 }
5686 else evaHV->push_back(rHadron); // Fill ResidNucl=Baryon to OutputHadronVector
5687 } // End of Evaporation of excited system
5688#ifdef debug
5690#endif
5691 }
5692 else // => "Decay if impossible evaporate" case
5693 {
5694#ifdef debug
5696#endif
5697 if(thePDG)
5698 {
5699 if(thePDG==10) // "Chipolino decay" case
5700 {
5703 G4QPDGCode h1=resChip.GetQPDG1();
5705 G4QPDGCode h2=resChip.GetQPDG2();
5707 if(totMass+.0001>m1+m2_value)
5708 {
5709#ifdef qdebug
5711#endif
5712 delete qH; // Chipolino should not be in a sequence
5713 G4LorentzVector fq4M(0.,0.,0.,m1);
5714 G4LorentzVector qe4M(0.,0.,0.,m2_value);
5716 {
5717 // G4cerr<<"***G4QNuc::EvaNuc:tM="<<totMass<<"-> h1M="<<m1<<" + h2M="<<m2_value<<G4endl;
5718 // throw G4QException("G4QNucleus::EvaporateNucleus: Chipol->h1+h2 DecIn2 error");
5719 G4ExceptionDescription ed;
5720 ed << "Chipol->h1+h2 DecIn2 error: tM=" << totMass << "-> h1M="
5723 FatalException, ed);
5724 }
5726#ifdef debug
5728#endif
5729 evaHV->push_back(H2); // (delete equivalent)
5731#ifdef debug
5733#endif
5734 evaHV->push_back(H1); // (delete equivalent)
5735 }
5736 else
5737 {
5738#ifdef qdebug
5740#endif
5741 delete qH;
5742 // G4cerr<<"**G4QN::EN:M="<<totMass<<"<"<<m1<<"+"<<m2_value<<",d="<<m1+m2_value-totMass<<G4endl;
5743 // throw G4QException("G4QNucleus::EvaporateNucleus: Chipolino is under MassShell");
5744 G4ExceptionDescription ed;
5745 ed << "Chipolino is under MassShell: M=" << totMass << "<" << m1
5748 FatalException, ed);
5749 }
5750 }
5751 else // "Hadron" case
5752 {
5754 if(fabs(totMass-totM)<0.001||abs(thePDG)-10*(abs(thePDG)/10)>2)
5755 {
5756#ifdef debug
5758#endif
5759 evaHV->push_back(qH);
5760 }
5761 else if ((thePDG==221||thePDG==331)&&totMass>mPi+mPi) // "Decay in pipi" case
5762 {
5763#ifdef qdebug
5765#endif
5766 delete qH;
5767 G4LorentzVector fq4M(0.,0.,0.,mPi);
5768 G4LorentzVector qe4M(0.,0.,0.,mPi);
5770 {
5771 // G4cerr<<"***G4QNucleus::EvaporateNucleus:tM="<<totMass<<"-> pi+ + pi-"<<G4endl;
5772 // throw G4QException("G4QNucleus::EvaporateNucleus: H->Pi+Pi DecayIn2 error");
5773 G4ExceptionDescription ed;
5774 ed << "H->Pi+Pi DecayIn2 error: tM=" << totMass << "-> pi+ + pi-"
5775 << G4endl;
5777 FatalException, ed);
5778 }
5780#ifdef debug
5782#endif
5783 evaHV->push_back(H1); // (delete equivalent)
5785#ifdef debug
5787#endif
5788 evaHV->push_back(H2); // (delete equivalent)
5789 }
5790 else if ((thePDG==221||thePDG==331)&&totMass>mPi0+mPi0) // "Decay in 2pi0" case
5791 {
5792#ifdef qdebug
5794#endif
5795 delete qH;
5796 G4LorentzVector fq4M(0.,0.,0.,mPi0);
5797 G4LorentzVector qe4M(0.,0.,0.,mPi0);
5799 {
5800 // G4cerr<<"***G4QNucleus::EvaporateNucleus:tM="<<totMass<<"-> pi0 + pi0"<<G4endl;
5801 // throw G4QException("G4QNucleus::EvaporateNucleus: H->Pi+Pi DecayIn2 error");
5802 G4ExceptionDescription ed;
5803 ed << "H->Pi+Pi DecayIn2 error: tM=" << totMass << "-> pi0 + pi0"
5804 << G4endl;
5806 FatalException, ed);
5807 }
5809#ifdef debug
5811#endif
5812 evaHV->push_back(H1); // (delete equivalent)
5814#ifdef debug
5816#endif
5817 evaHV->push_back(H2); // (delete equivalent)
5818 }
5819 else if (totMass>totM) // "Radiative Hadron decay" case
5820 {
5821#ifdef qdebug
5823#endif
5824 delete qH;
5825 G4LorentzVector fq4M(0.,0.,0.,0.);
5826 G4LorentzVector qe4M(0.,0.,0.,totM);
5828 {
5829 // G4cerr<<"***G4QNuc::EvaporateNuc:tM="<<totMass<<"->h1M="<<totM<<"+gam"<<G4endl;
5830 // throw G4QException("G4QNucleus::EvaporateNucleus: H*->H+gamma DecIn2 error");
5831 G4ExceptionDescription ed;
5832 ed << "H*->H+gamma DecIn2 error: tM=" << totMass << "->h1M="
5835 FatalException, ed);
5836 }
5838#ifdef debug
5840#endif
5841 evaHV->push_back(H2); // (delete equivalent)
5843#ifdef debug
5845#endif
5846 evaHV->push_back(H1); // (delete equivalent)
5847 }
5848 else if (thePDG==111||thePDG==221||thePDG==331) // "Gamma+Gamma decay" case
5849 {
5850#ifdef qdebug
5852#endif
5853 delete qH;
5854 G4LorentzVector fq4M(0.,0.,0.,0.);
5855 G4LorentzVector qe4M(0.,0.,0.,0.);
5857 {
5858 // G4cerr<<"***G4QNucl::EvaporateNucleus:tM="<<totMass<<"->gamma + gamma"<<G4endl;
5859 // throw G4QException("G4QNucleus::EvaporateNucleus: pi/eta->g+g DecIn2 error");
5860 G4ExceptionDescription ed;
5861 ed << "pi/eta->g+g DecIn2 error: tM=" << totMass
5864 FatalException, ed);
5865 }
5867#ifdef debug
5869#endif
5870 evaHV->push_back(H2); // (delete equivalent)
5872#ifdef debug
5874#endif
5875 evaHV->push_back(H1); // (delete equivalent)
5876 }
5877 else
5878 {
5879#ifdef qdebug
5881#endif
5882 delete qH;
5883 // G4cerr<<"***G4QNucl::EvaNuc: Nuc="<<thePDG<<theQC<<", q4M="<<q4M<<", M="<<totMass
5884 // <<" < GSM="<<totM<<", 2Pi="<<mPi+mPi<<", 2Pi0="<<mPi0+mPi0<<G4endl;
5885 // throw G4QException("G4QNucleus::EvaporateNucleus: Hadron is under MassShell");
5886 G4ExceptionDescription ed;
5887 ed << "Hadron is under MassShell: Nuc=" << thePDG << theQC
5888 << ", q4M=" << q4M << ", M=" << totMass <<" < GSM=" << totM
5891 FatalException, ed);
5892 }
5893 }
5894 }
5895 else
5896 {
5897#ifdef qdebug
5899#endif
5900 delete qH;
5901 // G4cerr<<"**G4QNuc::EvaNuc:RN="<<thePDG<<theQC<<",q4M="<<q4M<<",qM="<<totMass<<G4endl;
5902 // throw G4QException("G4QNucleus::EvaporateNucleus: This is not aNucleus nor aHadron");
5903 G4ExceptionDescription ed;
5904 ed << "This is not aNucleus nor aHadron: RN=" << thePDG << theQC
5907 FatalException, ed);
5908 }
5909 }
5910#ifdef qdebug
5911 if (qH)
5912 {
5915 else delete qH;
5916 }
5917#endif
5918#ifdef debug
5920#endif
5921 return;
5922} // End of EvaporateNucleus
Definition: G4QChipolino.hh:47
void DecayAlphaBar(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:7414
void DecayAntiDibaryon(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:6591
void DecayIsonucleus(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:5925
void DecayAlphaAlpha(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:7665
void DecayDibaryon(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:6286
void DecayAlphaDiN(G4QHadron *dB, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:7331
Referenced by DecayAntiStrange(), EvaporateNucleus(), G4QFragmentation::EvaporateResidual(), and G4QLowEnergy::PostStepDoIt().
◆ FissionCoulombBarrier()
| G4double G4QNucleus::FissionCoulombBarrier | ( | const G4double & | cZ, |
| const G4double & | cA, | ||
| G4double | dZ = 0., |
||
| G4double | dA = 0. |
||
| ) |
Definition at line 3403 of file G4QNucleus.cc.
3405{
3407 if(cZ<=0.) return 0.;
3409 if(dA) rA-=dA; // Reduce rA f CB is calculated for wounded nucleus
3410 G4double rZ=Z-cZ;
3411 if(delZ) rZ-=delZ; // Reduce rZ f CB is calculated for wounded nucleus
3413 G4double r=(pow(rA,third)+pow(cA,third))*(1.51+.00921*zz)/(1.+.009443*zz);
3414 return 1.44*zz/r;
3415} // End of "FissionCoulombBarier"
◆ GetA()
|
inline |
Definition at line 73 of file G4QNucleus.hh.
73{return Z+N+S;} // Get A of the nucleus
Referenced by ActivateBThickness(), G4QFragmentation::Breeder(), G4QIonIonCollision::Breeder(), ChooseFermiMomenta(), ChooseNucleons(), ChoosePositions(), CoulombBarrier(), G4QEnvironment::DecayBaryon(), G4QEnvironment::DecayMeson(), EvaporateBaryon(), EvaporateNucleus(), FissionCoulombBarrier(), G4QEnvironment::Fragment(), G4QFragmentation::Fragment(), G4QIonIonCollision::Fragment(), G4QEnvironment::G4QEnvironment(), G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), GetDeriv(), GetFermiMomentum(), GetNextNucleon(), GetRadius(), GetRelativeDensity(), GetThickness(), Init3D(), InitDensity(), ReduceSum(), SimpleSumReduction(), Split2Baryons(), and SplitBaryon().
◆ GetBThickness() [1/2]
Definition at line 114 of file G4QNucleus.hh.
114{return &Tb;} // T(b) function, step .1 fm
Referenced by GetThickness().
◆ GetBThickness() [2/2]
◆ GetDA()
|
inline |
Definition at line 77 of file G4QNucleus.hh.
77{return dZ+dN+dS;} // Get A of the dense part of nucleus
◆ GetDensity()
|
inline |
Definition at line 90 of file G4QNucleus.hh.
G4double GetRelativeDensity(const G4ThreeVector &aPosition)
Definition: G4QNucleus.cc:3757
Referenced by ChooseFermiMomenta(), ChoosePositions(), and GetDeriv().
◆ GetDeriv()
| G4double G4QNucleus::GetDeriv | ( | const G4ThreeVector & | point | ) |
Definition at line 3735 of file G4QNucleus.cc.
3736{
3738 G4double rPos=aPosition.mag();
3740 // Wood-Saxon density distribution
3742 return -exp((rPos-radius)/WoodSaxonSurf)*dens*dens*rho0/WoodSaxonSurf;
3743} // End of GetDeriv
◆ GetDN()
|
inline |
Definition at line 75 of file G4QNucleus.hh.
75{return dN;} // Get a#of neutrons in dense region
◆ GetDS()
|
inline |
Definition at line 76 of file G4QNucleus.hh.
76{return dS;} // Get a#of lambdas in dense region
◆ GetDZ()
|
inline |
Definition at line 74 of file G4QNucleus.hh.
74{return dZ;} // Get a#of protons in dense region
◆ GetFermiMomentum()
Definition at line 3765 of file G4QNucleus.cc.
3766{
3770} // End of GetFermiMomentum
Referenced by ChooseFermiMomenta(), and ChoosePositions().
◆ GetGSMass()
|
inline |
Definition at line 82 of file G4QNucleus.hh.
Referenced by BindingEnergy(), G4QFragmentation::Breeder(), G4QIonIonCollision::Breeder(), CalculateMass(), ChooseFermiMomenta(), CoulBarPenProb(), EvaporateBaryon(), EvaporateNucleus(), G4QFragmentation::EvaporateResidual(), G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), G4QNucleus(), GetNucleons4Momentum(), operator<<(), and SubtractNucleon().
◆ GetMaxClust()
|
inline |
Definition at line 78 of file G4QNucleus.hh.
78{return maxClust;} // Get Max BarNum of Clusters
◆ GetMZNS()
|
inline |
Definition at line 80 of file G4QNucleus.hh.
Referenced by EvaporateNucleus(), G4QEnvironment::G4QEnvironment(), and G4QDiffractionRatio::ProjFragment().
◆ GetN()
|
inline |
Definition at line 71 of file G4QNucleus.hh.
71{return N;} // Get a#of neutrons
Referenced by EvaporateNucleus(), G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), and operator<<().
◆ GetNDefBaryonC()
|
inline |
Definition at line 89 of file G4QNucleus.hh.
89{return nDefBaryonC;};// max#of predefed baryonCandidates
◆ GetNDefMesonC()
|
inline |
Definition at line 88 of file G4QNucleus.hh.
88{return nDefMesonC;}; // max#of predefed mesonCandidates
◆ GetNextNucleon()
|
inline |
Definition at line 100 of file G4QNucleus.hh.
101 {
102 //G4cout<<"G4QNucleus::GetNextNucleon: cN="<<currentNucleon<<", A="<<GetA()<<G4endl;
104 }
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ GetNucleons4Momentum()
|
inline |
Definition at line 107 of file G4QNucleus.hh.
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ GetOuterRadius()
| G4double G4QNucleus::GetOuterRadius | ( | ) |
Definition at line 3947 of file G4QNucleus.cc.
3948{
3949 G4double maxradius2=0;
3950 G4int theA=theNucleons.size();
3952 {
3953 G4double nucr2=theNucleons[i]->GetPosition().mag2();
3954 if(nucr2 > maxradius2) maxradius2=nucr2;
3955 }
3956 return sqrt(maxradius2)+nucleonDistance;
3957} // End of GetOuterRadius
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ GetPDG()
|
inline |
Definition at line 69 of file G4QNucleus.hh.
69{return 90000000+1000*(1000*S+Z)+N;}// PDG Code of Nucleus
Referenced by G4QFragmentation::Breeder(), G4QIonIonCollision::Breeder(), ChooseFermiMomenta(), EvaporateNucleus(), G4QFragmentation::Fragment(), G4QIonIonCollision::Fragment(), G4QEnvironment::G4QEnvironment(), G4QFragmentation::G4QFragmentation(), Increase(), and Reduce().
◆ GetProbability()
Definition at line 79 of file G4QNucleus.hh.
79{return probVect[bn];} // clust(BarN)probabil
Referenced by G4Quasmon::Fragment().
◆ GetQCZNS()
|
inline |
Definition at line 83 of file G4QNucleus.hh.
84 {
87 }
Referenced by EvaporateNucleus(), G4QFragmentation::Fragment(), G4QIonIonCollision::Fragment(), Split2Baryons(), and SplitBaryon().
◆ GetRadius()
Definition at line 3746 of file G4QNucleus.cc.
3747{
3751 // Wood-Saxon density distribution
3752 return (maxRelDens>0 && maxRelDens <= 1. ) ? (radius + WoodSaxonSurf*
3753 log((1.-maxRelDens+exp(-radius/WoodSaxonSurf))/maxRelDens) ) : DBL_MAX;
3754} // End of GetRadius (check @@ radius=sqr0 (fm^2) for A<17, r0 (fm) for A>16 (units)
Referenced by ChoosePositions().
◆ GetRelativeDensity()
| G4double G4QNucleus::GetRelativeDensity | ( | const G4ThreeVector & | aPosition | ) |
Definition at line 3757 of file G4QNucleus.cc.
3758{
3762} // End of GetRelativeDensity
double mag() const
Referenced by GetDensity(), and GetDeriv().
◆ GetRelOMDensity()
Definition at line 95 of file G4QNucleus.hh.
95{return std::exp(-r2/radius);} // OscModelRelDens
Referenced by ChoosePositions(), and GetRelativeDensity().
◆ GetRelWSDensity()
Definition at line 93 of file G4QNucleus.hh.
94 {return 1./(1.+std::exp((r-radius)/WoodSaxonSurf));}
Referenced by ChoosePositions(), and GetRelativeDensity().
◆ GetRho0()
|
inline |
Definition at line 91 of file G4QNucleus.hh.
91{return rho0;} // One nucleon prob-density
◆ GetS()
|
inline |
Definition at line 72 of file G4QNucleus.hh.
72{return S;} // Get a#of lambdas
Referenced by EvaporateNucleus(), and operator<<().
◆ GetTbIntegral()
| G4double G4QNucleus::GetTbIntegral | ( | ) |
Definition at line 4018 of file G4QNucleus.cc.
◆ GetThickness()
Definition at line 4048 of file G4QNucleus.cc.
4049{
4051 if(tA<1)
4052 {
4054 return 0.;
4055 }
4056 else if(tA==1) return 0.;
4060}
std::vector< G4double > const * GetBThickness()
Definition: G4QNucleus.hh:114
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ GetZ()
|
inline |
Definition at line 70 of file G4QNucleus.hh.
70{return Z;} // Get a#of protons
Referenced by G4QFragmentation::Breeder(), G4QIonIonCollision::Breeder(), EvaporateNucleus(), G4QFragmentation::Fragment(), G4QIonIonCollision::Fragment(), G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), and operator<<().
◆ HadrToNucPDG()
Definition at line 4117 of file G4QNucleus.cc.
4118{
4119 G4int nPDG=hPDG;
4120 if (hPDG==2212) nPDG=90001000; // p
4121 else if(hPDG==2112) nPDG=90000001; // n
4122 else if(hPDG==3122||hPDG==3212) nPDG=91000000; // Lambda
4123 else if(hPDG== 211) nPDG=90000999; // pi+
4124 else if(hPDG==-211) nPDG=89999001; // pi-
4125 else if(hPDG== 311) nPDG=89000001; // K0 (anti-strange)
4126 else if(hPDG== 321) nPDG=89001000; // K+ (anti-strange)
4127 else if(hPDG==-311) nPDG=90999999; // anti-K0 (strange)
4128 else if(hPDG==-321) nPDG=90999000; // K- (strange)
4129 else if(hPDG==1114) nPDG=89999002; // Delta-
4130 else if(hPDG==2224) nPDG=90001999; // Delta++
4131 else if(hPDG==3112) nPDG=90999000; // Sigma-
4132 else if(hPDG==3222) nPDG=91000999; // Sigma+
4133 else if(hPDG==3312) nPDG=91999000; // Ksi-
4134 else if(hPDG==3322) nPDG=91999999; // Ksi0
4135 else if(hPDG==3334) nPDG=92998999; // Omega-
4136 else if(hPDG==-2212) nPDG=8999000; // anti-p
4137 else if(hPDG==-2112) nPDG=8999999; // anti-n
4138 else if(hPDG==-3122||hPDG==3212) nPDG=89000000; //anti- Lambda
4139 else if(hPDG==-3112) nPDG=89000999; // anti-Sigma-
4140 else if(hPDG==-3222) nPDG=88999001; // anti-Sigma+
4141 else if(hPDG==-3312) nPDG=88001000; // anti-Ksi-
4142 else if(hPDG==-3322) nPDG=88000001; // anti-Ksi0
4143 else if(hPDG==-3334) nPDG=87001001; // anti-Omega-
4144 return nPDG;
4145}
Referenced by InitByPDG().
◆ IncProbability()
| void G4QNucleus::IncProbability | ( | G4int | bn | ) |
◆ Increase() [1/2]
| void G4QNucleus::Increase | ( | G4int | PDG, |
| G4LorentzVector | LV = G4LorentzVector(0.,0.,0.,0.) |
||
| ) |
Definition at line 729 of file G4QNucleus.cc.
730{
732 if(cPDG>80000000&&cPDG!=NUCPDG)
733 {
737 {
739 t4M +=c4M;
740 Set4Momentum(t4M);
741 }
742 }
744}
◆ Increase() [2/2]
| void G4QNucleus::Increase | ( | G4QContent | QC, |
| G4LorentzVector | LV = G4LorentzVector(0.,0.,0.,0.) |
||
| ) |
Definition at line 747 of file G4QNucleus.cc.
748{
752 t4M +=q4M;
754}
◆ Init3D()
| void G4QNucleus::Init3D | ( | ) |
Definition at line 3910 of file G4QNucleus.cc.
3911{
3912#ifdef debug
3914#endif
3915 for_each(theNucleons.begin(),theNucleons.end(),DeleteQHadron());
3916 theNucleons.clear();
3918 ChooseNucleons();
3919#ifdef debug
3921#endif
3922 InitDensity();
3923#ifdef debug
3925#endif
3927#ifdef debug
3929#endif
3933 {
3936 }
3937#ifdef debug
3939#endif
3942 currentNucleon=0; // Automatically starts the LOOP
3943 return;
3944} // End of Init3D
double z() const
double x() const
double y() const
G4double BindingEnergy(const G4double &cZ=0, const G4double &cA=0, G4double dZ=0., G4double dA=0.)
Definition: G4QNucleus.cc:3418
void DoLorentzBoost(const G4LorentzVector &theBoost)
Definition: G4QNucleus.hh:155
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ InitByPDG()
| void G4QNucleus::InitByPDG | ( | G4int | newPDG | ) |
Definition at line 371 of file G4QNucleus.cc.
372{
374#ifdef debug
376#endif
377 dZ=0;
378 dN=0;
379 dS=0;
380 probVect[0]=mediRatio; // init Vacuum/Medium probability
382 //std::uninitialized_fill( probVect+1, probVect+256, 0.0 ); // Worse in performance!
384 G4int s_value=0;
385 G4int z=0;
387 if(nucPDG>80000000 && nucPDG<100000000) // Try to convert the NUCCoding to PDGCoding
388 {
390 Z =z;
391 N =n;
392 S =s_value;
393#ifdef debug
395#endif
396 SetZNSQC(Z,N,S); // @@ ??
397 G4QPDGCode nPDG(nucPDG);
398 G4double PDGMass=0.;
399 if(nucPDG!=NUCPDG) PDGMass=nPDG.GetMass();
400 SetQPDG(nPDG);
401 G4LorentzVector p(0.,0.,0.,PDGMass);
402 Set4Momentum(p);
403 SetNFragments(0);
404#ifdef debug
406#endif
407 }
408 else
409 {
411 //throw G4QException("G4QNucleus::InitByPDG:PDGCode can't be converted to NucPDGCode");
412 }
413}
void ConvertPDGToZNS(G4int PDG, G4int &z, G4int &n, G4int &s)
Definition: G4QPDGCode.cc:2377
Referenced by G4QEnvironment::G4QEnvironment(), G4QFragmentation::G4QFragmentation(), G4QIonIonCollision::G4QIonIonCollision(), G4QNucleus(), Increase(), InitByQC(), Reduce(), and SubtractNucleon().
◆ InitByQC()
|
inline |
◆ InitCandidateVector()
| void G4QNucleus::InitCandidateVector | ( | G4QCandidateVector & | theQCandidates, |
| G4int | nM = 45, |
||
| G4int | nB = 72, |
||
| G4int | nC = 117 |
||
| ) |
Definition at line 3037 of file G4QNucleus.cc.
3039{
3042 // Scalar resonances (0): Eta,Pi0,Pi+,APi-,Ka0,Ka+,AKa0,AKa-,Eta*
3044 // Vector resonances (1): omega,Rh0,Rh+,Rho-,K0*,K+*,AK0*,AK-*,Phi
3045 223,113,213,-213,313,323,-313,-323,333, // 9-18
3046 // Tensor D-resonances (2): f2 ,a20,a2+, a2-,K20,K2+,AK20,AK2-,f2'
3047 225,115,215,-215,315,325,-315,-325,335, // 19-27
3048 // Tensor F-resonances (3): om3,ro3,r3+,rh3-,K30,K3+,AK30,AK3-,Phi3
3049 227,117,217,-217,317,327,-317,-327,337, // 28-35
3050 // Tensor G-resonances (4): f4 ,a40,a4+, a4-,K40,K4+,AK40,AK4-,f4'
3051 229,119,219,-219,319,329,-319,-329,339}; // 36-44
3052 // Baryon octet (1/2): n , an , p , ap ,lamb,alamb, sig-,asig-
3054 // Hyperon octet (1/2): sig0,asig0,sig+,asig+,ksi-,aksi-,ksi0,aksi0
3055 3212,-3212,3222,-3222,3312,-3312,3322,-3322, // 53-60
3056 // Baryon decuplet (3/2): del-,adel-,del0,adel0,del+,adel+,dl++,adl++,sis-,asis-
3057 1114,-1114,2114,-2114,2214,-2214,2224,-2224,3114,-3114,//70
3058 // sis0,asis0,sis+,asis+,kss-,akss-,kss0,akss0,omeg,aomeg
3059 3214,-3214,3224,-3224,3314,-3314,3324,-3324,3334,-3334,//80
3060 // Baryon octet (5/2): n5/2,an5/2,p5/2,ap5/2,l5/2,al5/2,si5-,asi5-
3061 2116,-2116,2216,-2216,3126,-3126,3116,-3116, // 81-88
3062 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50
3063 3216,-3216,3226,-3226,3316,-3316,3326,-3326, // 89-96
3064 // Baryon decuplet (7/2): dl5-,adl5-,dl50,adl50,dl5+,adl5+,d5++,ad5++,si5-,asi5-
3065 1118,-1118,2118,-2118,2218,-2218,2228,-2228,3118,-3118, //106
3066 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50,ome5,aome5
3067 3218,-3218,3228,-3228,3318,-3318,3328,-3328,3338,-3338};//116
3068 G4int i=0;
3069#ifdef debug
3070 G4int ind=0;
3071#endif
3072 G4int iQC = theQCandidates.size();
3074 theQCandidates.clear();
3075 if(maxMes>nOfMesons) maxMes=nOfMesons;
3076 if(maxMes>=0) for (i=0; i<maxMes; i++)
3077 {
3079#ifdef debug
3082#endif
3083 }
3084 if(maxBar>nOfBaryons) maxBar=nOfBaryons;
3085 if(maxBar>=0) for (i=0; i<maxBar; i++)
3086 {
3087#ifdef debug
3089#endif
3091#ifdef debug
3093#endif
3094 theQCandidates.push_back(curBar); // delete equivalent
3095#ifdef debug
3098#endif
3099 }
3100 if(maxClst>=0) for (i=0; i<maxClst; i++)
3101 {
3103 G4QPDGCode clustQPDG;
3104 clustQPDG.InitByQCode(clustQCode);
3107#ifdef debug
3110#endif
3111 }
3112} // End of "InitCandidateVector"
Definition: G4QCandidate.hh:46
◆ InitDensity()
| void G4QNucleus::InitDensity | ( | ) |
Definition at line 3692 of file G4QNucleus.cc.
3693{
3697 G4double rA = iA;
3698 G4double At = pow(rA,third);
3699 G4double At2= At*At;
3700#ifdef debug
3702#endif
3703 if(iA<17) // Gaussian density distribution
3704 {
3705 radius = r0sq*At2; // R2 Mean Squared Radius (fm^2)
3706 if(radius<=0.)
3707 {
3709 radius=1.;
3710 }
3711 rho0 = pow(2*pi*radius, -1.5); // Central Density (M.K. 2 is added)
3712 // V=4pi*R2*sqrt(pi*R2/2)=(sqrt(2*pi*R2))^3
3713 }
3714 else // Wood-Saxon density distribution
3715 {
3717 radius = r0*At; // Half Density Radius (fm)
3718 if(radius<=0.)
3719 {
3721 radius=1.;
3722 }
3724 if(!(rd<=0.1) && !(rd>-0.1)) // NAN for rd
3725 {
3727 rd=1.;
3728 }
3730 }
3731 RhoActive=true;
3732} // End of InitDensity
Referenced by GetDeriv(), GetRadius(), GetRelativeDensity(), GetThickness(), and Init3D().
◆ NucToHadrPDG()
Definition at line 4148 of file G4QNucleus.cc.
4149{
4150 G4int hPDG=nPDG;
4151 if (nPDG==90001000) hPDG=2212; // p
4152 else if(nPDG==90000001) hPDG=2112; // n
4153 else if(nPDG==91000000) hPDG=3122; // Lambda
4154 else if(nPDG==90000999) hPDG= 211; // pi+
4155 else if(nPDG==89999001) hPDG=-211; // pi-
4156 else if(nPDG==89001000) hPDG= 213; // K0 (anti-strange)
4157 else if(nPDG==89000001) hPDG= 213; // K+ (anti-strange)
4158 else if(nPDG==90999000) hPDG=-213; // anti-K0 (strange)
4159 else if(nPDG==90999999) hPDG=-213; // K- (strange)
4160 else if(nPDG==90001999) hPDG=1114; // Delta-
4161 else if(nPDG==89999002) hPDG=2224; // Delta++
4162 else if(nPDG==91000999) hPDG=3112; // Sigma-
4163 else if(nPDG==90999001) hPDG=3222; // Sigma+
4164 else if(nPDG==91999999) hPDG=3312; // Ksi-
4165 else if(nPDG==91999000) hPDG=3322; // Ksi0
4166 else if(nPDG==92998999) hPDG=3334; // Omega-
4167 return hPDG;
4168}
◆ operator!=()
|
inline |
Definition at line 67 of file G4QNucleus.hh.
67{return this!=&right;}
◆ operator*=()
| G4QNucleus G4QNucleus::operator*= | ( | const G4int & | rhs | ) |
Definition at line 4099 of file G4QNucleus.cc.
4100{
4101 Z*=rhs;
4102 N*=rhs;
4103 S*=rhs;
4104 dZ*=rhs;
4105 dN*=rhs;
4106 dS*=rhs;
4107 // Atributes of aHadron
4109 SetQPDG (newPDG);
4111 SetQC (newQC);
4112 theMomentum *= rhs;
4113 return *this;
4114}
◆ operator+=()
| G4QNucleus G4QNucleus::operator+= | ( | const G4QNucleus & | rhs | ) |
Definition at line 4063 of file G4QNucleus.cc.
◆ operator-=()
| G4QNucleus G4QNucleus::operator-= | ( | const G4QNucleus & | rhs | ) |
Definition at line 4081 of file G4QNucleus.cc.
◆ operator=()
| const G4QNucleus & G4QNucleus::operator= | ( | const G4QNucleus & | right | ) |
Definition at line 307 of file G4QNucleus.cc.
308{
309 if(this != &right) // Beware of self assignment
310 {
311 currentNucleon= -1;
312 TbActive = right.TbActive;
313 Tb = right.Tb;
314 RhoActive = right.RhoActive;
315 rho0 = right.rho0;
316 radius = right.radius;
319 {
321 theNucleons.push_back(nucleon);
322 }
327 Z = right.Z;
328 N = right.N;
329 S = right.S;
330 dZ = right.dZ;
331 dN = right.dN;
332 dS = right.dS;
333 maxClust = right.maxClust;
335 probVect[254] = right.probVect[254];
336 probVect[255] = right.probVect[255];
337 }
338 return *this;
339}
◆ operator==()
|
inline |
Definition at line 66 of file G4QNucleus.hh.
66{return this==&right;}
◆ PrepareCandidates()
| void G4QNucleus::PrepareCandidates | ( | G4QCandidateVector & | theQCandidates, |
| G4bool | piF = false, |
||
| G4bool | gaF = false, |
||
| G4LorentzVector | LV = G4LorentzVector(0.,0.,0.,0.) |
||
| ) |
Definition at line 3115 of file G4QNucleus.cc.
3117{
3119 G4double ze = Z;
3120 G4double ne = N;
3121 G4double se = S;
3122 G4double ae = Z+N+S;
3123 G4double aea = ae*(ae-1);
3124 G4double ae2 = aea/2.;
3125 G4double ze1 = dZ + 1;
3126 G4double ne1 = dN + 1;
3127 G4double se1 = dS + 1;
3128 G4double ae0 = dZ + dN + dS;
3129 G4double ae1 = ae0 + 1;
3131#ifdef cldebug
3133#endif
3139 if(comb<=0.) comb=1.;
3140#ifdef cldebug
3144#endif
3146 {
3147 G4QCandidate* curCand=theQCandidates[index];
3151 // ***********************************************************************************
3152 // These are first 117 candidates which are defined in G4QNucleus::InitCandidateVector
3153 // ***!!!*** if they are changed there the corresponding change must be done here
3154 //static const G4int nOfMesons =45;//a#of S=0,1,2,3,4 Mesons, => candidates to hadrons
3155 //static const G4int nOfBaryons=72;//a#of 1/2,3/2,5/2,7/2 Baryons => cand's to hadrons
3156 // Scalar resonances (0): Eta,Pi0,Pi+,APi-,Ka0,Ka+,AKa0,AKa-,Eta*
3157 //static G4int mesonPDG[45] = {221,111,211,-211,311,321,-311,-321,331, // 0- 8
3158 // Vector resonances (1): omega,Rh0,Rh+,Rho-,K0*,K+*,AK0*,AK-*,Phi
3159 // 223,113,213,-213,313,323,-313,-323,333, // 9-18
3160 // Tensor D-resonances (2): f2 ,a20,a2+, a2-,K20,K2+,AK20,AK2-,f2'
3161 // 225,115,215,-215,315,325,-315,-325,335, // 19-27
3162 // Tensor F-resonances (3): om3,ro3,r3+,rh3-,K30,K3+,AK30,AK3-,Phi3
3163 // 227,117,217,-217,317,327,-317,-327,337, // 28-35
3164 // Tensor G-resonances (4): f4 ,a40,a4+, a4-,K40,K4+,AK40,AK4-,f4'
3165 // 229,119,219,-219,319,329,-319,-329,339}; // 36-44
3166 // Baryon octet (1/2): n , an , p , ap ,lamb,alamb, sig-,asig-
3167 //static G4int baryonPDG[72]={2112,-2112,2212,-2212,3122,-3122,3112,-3112, // 45-52
3168 // Hyperon octet (1/2): sig0,asig0,sig+,asig+,ksi-,aksi-,ksi0,aksi0
3169 // 3212,-3212,3222,-3222,3312,-3312,3322,-3322, // 53-60
3170 // Baryon decuplet (3/2): del-,adel-,del0,adel0,del+,adel+,dl++,adl++,sis-,asis-
3171 // 1114,-1114,2114,-2114,2214,-2214,2224,-2224,3114,-3114,//70
3172 // sis0,asis0,sis+,asis+,kss-,akss-,kss0,akss0,omeg,aomeg
3173 // 3214,-3214,3224,-3224,3314,-3314,3324,-3324,3334,-3334,//80
3174 // Baryon octet (5/2): n5/2,an5/2,p5/2,ap5/2,l5/2,al5/2,si5-,asi5-
3175 // 2116,-2116,2216,-2216,3126,-3126,3116,-3116, // 81-88
3176 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50
3177 // 3216,-3216,3226,-3226,3316,-3316,3326,-3326, // 89-96
3178 // Baryon decuplet (7/2): dl5-,adl5-,dl50,adl50,dl5+,adl5+,d5++,ad5++,si5-,asi5-
3179 // 1118,-1118,2118,-2118,2218,-2218,2228,-2228,3118,-3118, //106
3180 // si50,asi50,si5+,asi5+,ks5-,aks5-,ks50,aks50,ome5,aome5
3181 // 3218,-3218,3228,-3228,3318,-3318,3328,-3328,3338,-3338};//116
3182 // One should take into account, that #of mesons & baryons can be cut in G4Quas::HadrQE
3183 //G4int nP= theWorld->GetQPEntries(); // A#of initialized particles in CHIPS World
3184 ////@@ Make special parametyer to cut high resonances for nuclear fragmentation !!
3185 //G4int nMesons = 45;
3186 //if (nP<34) nMesons = 9;
3187 //else if(nP<51) nMesons = 18;
3188 //else if(nP<65) nMesons = 27;
3189 //else if(nP<82) nMesons = 36;
3190 //G4int nBaryons = 72;
3191 //if (nP<45) nBaryons = 16;
3192 //else if(nP<59) nBaryons = 36;
3193 //else if(nP<76) nBaryons = 52;
3194 // **********************************************************************************
3195 //G4int cS = curCand->GetStrangeness();
3196 //if(piF&&gaF&&cPDG!=90000001&&cPDG!=90001000) // Both flags, in case of pi-first-int
3197 //if(piF&&gaF&&cBN!=1&&cBN!=3) // Both flags, which is in case of pi-first-int
3198 if(piF&&gaF&&cBN!=1)// @@ Should be both, which is in case of pi-first-interaction @@ ?
3199 //if(piF&&gaF&&cBN!=1&&cBN!=4) // Should be both, in case of pi-first-interaction
3200 {
3201 curCand->SetPreProbability(0.);
3202 curCand->SetDenseProbability(0.);
3204#ifdef cldebug
3206#endif
3207 }
3208 // @@ in case of the Ksi or Omega- capture it can disturb the simulation
3209 else if(cPDG<80000000&&(abs(cPDG)%10>4||cST>2))// @@ PreClosed HighSpin/HighStrange
3210 {
3211 curCand->SetPreProbability(0.);
3212 curCand->SetDenseProbability(0.);
3214 }
3215 else
3216 {
3218 if(cPDG>80000000&&cPDG!=90000000) // ===> Cluster case
3219 {
3226 pos = probVect[ac]; // Cluster Probability NormalizationFactor
3228 {
3229
3230#ifdef cldebug
3232#endif
3233 G4int dac=ac+ac;
3234 if(ac && (piF || gaF)) // zc>=0
3235 {
3236 if (piF&&!gaF&&zc+ac) pos*=(zc+ac)/ac; // piF interaction (#of u-quarks)
3237 else if(gaF&&!piF&&zc+dac) pos*=(zc+dac)/ac; // gaF interaction (sum of Q_q^2)
3238 }
3239 G4double dense=1.;
3240 if (ac==1&&pos>0.)dense=probVect[254]/pos;
3241 else if(ac==2&&pos>0.)dense=probVect[255]/pos;
3242#ifdef cldebug
3245 G4double mp=pos;
3246#endif
3247 if (ac==1 && ae) // ae=0 protection (since here no /pos)
3248 {
3249 if (zc) pos*=ze/ae;
3250 else if(nc) pos*=ne/ae;
3251 else if(sc) pos*=se/ae;
3252 //if (zc) pos*=ze;
3253 //else if(nc) pos*=ne;
3254 //else if(sc) pos*=se;
3255 acm=1;
3256#ifdef cldebug
3257 if(pos)
3260 sZ+=pos*zc;
3261 sN+=pos*nc;
3262#endif
3263 }
3264 else if(ac==2)
3265 {
3266 if(ze<zc||ne<nc||se<sc) pos=0.;
3267 else if(aea) // Protection against aea=0.
3268 {
3269 if (zc==2) pos*=ze*(ze-1)/aea;
3270 else if(nc==2) pos*=ne*(ne-1)/aea;
3271 else if(sc==2) pos*=se*(se-1)/aea;
3272 else if(zc==1&&nc==1) pos*=ze*ne/ae2;
3273 else if(zc==1&&sc==1) pos*=ze*se/ae2;
3274 else if(sc==1&&nc==1) pos*=se*ne/ae2;
3275 //if (zc==2) pos*=ze*(ze-1)/2.;
3276 //else if(nc==2) pos*=ne*(ne-1)/2.;
3277 //else if(sc==2) pos*=se*(se-1)/2.;
3278 //else if(zc==1&&nc==1) pos*=ze*ne;
3279 //else if(zc==1&&sc==1) pos*=ze*se;
3280 //else if(sc==1&&nc==1) pos*=se*ne;
3282 // Normalization for only not strange matter
3283 }
3284 acm=2;
3285#ifdef cldebug
3286 if(pos)
3288 <<mac<<G4endl;
3289 sZ+=pos*zc;
3290 sN+=pos*nc;
3291#endif
3292 }
3293 else // ac>2
3294 {
3295 if(acm<ac) // first time that big cluster
3296 {
3297 if(ac<ae1 && ac>0) comb*=(ae1-ac)/ac;
3298 acm=ac;
3299 s_value=0.;
3300 cca=0;
3301#ifdef cldebug
3302 if(ac%2) mac=7; // @@Change it if cluster set is changed
3303 else mac=8; // @@ It is not yet automatic
3305 <<ae0<<G4endl;
3306#endif
3307 }
3308 G4double prod=1.;
3309 if(ze1<=zc||ne1<=nc||se1<=sc) prod=0.;
3310 else
3311 {
3312 if(zc>0) for(int iz=1; iz<=zc; iz++) prod*=(ze1-iz)/iz;
3313 if(nc>0) for(int in=1; in<=nc; in++) prod*=(ne1-in)/in;
3314 if(sc>0) for(int is=1; is<=sc; is++) prod*=(se1-is)/is;
3315 }
3316 s_value+=prod;
3317 pos*=prod;
3318 pos/=comb;
3319#ifdef cldebug
3322 sZ+=pos*zc;
3323 sN+=pos*nc;
3324#endif
3325 cca++;
3326 }
3327 curCand->SetPreProbability(pos);
3328 curCand->SetDenseProbability(pos*dense);
3329#ifdef cldebug
3331#endif
3332 }
3333 else // => "Cluster is too big" case
3334 {
3335 curCand->SetPreProbability(0.);
3336 curCand->SetDenseProbability(0.);
3338 }
3339 }
3340 else
3341 {
3342#ifdef cldebug
3344#endif
3347 }
3349 }
3350 } // End of the LOOP over Candidates
3351#ifdef cldebug
3353 //throw G4QException("G4QNucleus::PrepareCandidate: Temporary stop");
3354#endif
3355}// End of PrepareCandidates
void SetDenseProbability(G4double prep)
Definition: G4QCandidate.hh:132
◆ RandomizeBinom()
Definition at line 3015 of file G4QNucleus.cc.
3016{
3018 G4double d = 1.-p;
3019 if(d<=0.) return 0;
3020 G4double v = pow(d,aN);
3021 G4double s_value = v;
3022 if(r<s_value) return 0;
3023 G4int i=0;
3024 G4int j=aN+1;
3025 G4double f=p/d;
3026 while(r>s_value && i<aN)
3027 {
3028 j--;
3029 i++;
3030 v*=j*f/i;
3031 s_value+=v;
3032 }
3033 return i;
3034}
Referenced by UpdateClusters().
◆ Reduce()
| void G4QNucleus::Reduce | ( | G4int | PDG | ) |
Definition at line 705 of file G4QNucleus.cc.
706{
708 if(cPDG>80000000&&cPDG!=NUCPDG)
709 {
713 else
714 {
715 //if(abs(newPDG)<NUCPDG)
716 //{
717 // G4cerr<<"***G4QNucleus::Reduce:iPDG="<<curPDG<<"=newPDG="<<newPDG<<"+cPDG="<<cPDG
718 // <<G4endl;
719 // throw G4QException("*E*:::G4QNucleus::Reduce: Abnormal Nuclear Reduction");
720 //}
722 }
723 }
725 // in case of cPDG=90000000 - subtract nothing
726}
◆ ReduceSum()
| G4bool G4QNucleus::ReduceSum | ( | G4ThreeVector * | vectors, |
| G4ThreeVector | sum | ||
| ) |
Definition at line 3861 of file G4QNucleus.cc.
3862{
3864 if(theA<3) // Can not reduce for 1 or 2 nucleons
3865 {
3867 return false;
3868 }
3869 // The last vector must have the same direction as the SUM (do not take into account
3878 while(m_value)
3879 {
3880 m_value=0;
3882 {
3883 m_value++;
3884 G4double shift=fabs(dp[i]-hsum2);
3885 if(shift < minS)
3886 {
3887 minS=shift;
3888 minI=i;
3889 }
3890 }
3891 if(m_value) // There is a vector reducing the sum
3892 {
3894 vect[minI]-=x; // turn the minI-th vector
3895 sum-=x; // reduce the sum
3897 hsum2=sum2/2; // Current half squared sum
3898 }
3899 }
3900 if(sum2 > 0.)
3901 {
3902 sum/=theA;
3904 }
3905 delete[] dp;
3906 return true;
3907} // End of ReduceSum
double dot(const Hep3Vector &) const
Referenced by ChoosePositions().
◆ SetMaxClust()
|
inline |
Definition at line 143 of file G4QNucleus.hh.
143{maxClust=maxC;}// Set Max BarNum of Clusters
◆ SetParameters()
|
static |
Definition at line 347 of file G4QNucleus.cc.
348{
349 freeNuc=a;
350 freeDib=b;
351 clustProb=c;
352 mediRatio=d;
353 nucleonDistance=e;
354}
Referenced by G4ChiralInvariantPhaseSpace::ApplyYourself(), G4QCaptureAtRest::AtRestDoIt(), G4QAtomicElectronScattering::G4QAtomicElectronScattering(), G4QCaptureAtRest::G4QCaptureAtRest(), G4QInelastic::G4QInelastic(), G4QAtomicElectronScattering::PostStepDoIt(), G4QInelastic::PostStepDoIt(), G4QStringChipsParticleLevelInterface::Propagate(), G4StringChipsInterface::Propagate(), G4StringChipsParticleLevelInterface::Propagate(), G4QAtomicElectronScattering::SetParameters(), G4QCaptureAtRest::SetParameters(), and G4QInelastic::SetParameters().
◆ SimpleSumReduction()
| void G4QNucleus::SimpleSumReduction | ( | G4ThreeVector * | vectors, |
| G4ThreeVector | sum | ||
| ) |
Definition at line 3853 of file G4QNucleus.cc.
3854{
3856 sum/=theA;
3858}
Referenced by ChooseFermiMomenta().
◆ Split2Baryons()
| G4bool G4QNucleus::Split2Baryons | ( | ) |
Definition at line 863 of file G4QNucleus.cc.
864{
872 if(baryn<3) return false;
875#ifdef debug
877#endif
879 if(NQ>1) // ===> "Can try to split 2 neutrons" case
880 {
884 G4double sM=resMas+mNeut+mNeut;
885#ifdef debug
887#endif
888 if(sM<totM) return true;
889 }
891 if(PQ>1) // ===> "Can try to split 2 protons" case
892 {
896 G4double sM=resMas+mProt+mProt;
897#ifdef debug
899#endif
900 if(sM<totM) return true;
901 }
902 if(PQ&&NQ) // ===> "Can try to split proton+neutron"
903 {
907 G4double sM=resMas+mProt+mNeut;
908#ifdef debug
910#endif
911 if(sM<totM) return true;
912 }
914 if(LQ&&NQ) // ===> "Can try to split lambda+neutron"
915 {
919 G4double sM=resMas+mLamb+mNeut;
920#ifdef debug
922#endif
923 if(sM<totM) return true;
924 }
925 if(LQ&&PQ) // ===> "Can try to split lambda+proton"
926 {
930 G4double sM=resMas+mProt+mLamb;
931#ifdef debug
933#endif
934 if(sM<totM) return true;
935 }
936 if(LQ>1) // ===> "Can try to split 2 lambdas" case
937 {
941 G4double sM=resMas+mLamb+mLamb;
942#ifdef debug
944#endif
945 if(sM<totM) return true;
946 }
947 return false;
948}
Referenced by EvaporateNucleus().
◆ SplitBaryon()
| G4int G4QNucleus::SplitBaryon | ( | ) |
Definition at line 774 of file G4QNucleus.cc.
775{
787 if(baryn<2) return 0;
788 //G4double totM=GetGSMass(); // GS Mass value of the Nucleus
791#ifdef debug
793#endif
795 if(NQ) // ===> "Can try to split a neutron" case
796 {
800 G4double sM=resMas+mNeut;
801#ifdef debug
803#endif
804 if(sM<totM+.001) return 2112;
805 }
807 if(PQ) // ===> "Can try to split a proton" case
808 {
813 G4double sM=resMas+mProt+CB;
814 /////////G4double sM=resMas+mProt;
815#ifdef debug
817#endif
818 if(sM<totM+.001) return 2212;
819 }
821 if(LQ) // ===> "Can try to split a lambda" case
822 {
826 G4double sM=resMas+mLamb;
827#ifdef debug
829#endif
830 if(sM<totM+.001) return 3122;
831 }
832 G4int AQ=NQ+PQ+LQ;
833 if(NQ>0&&PQ>0&&AQ>2) // ===> "Can try to split deuteron" case
834 {
839 G4double sM=resMas+mDeut+CB;
840 //G4double sM=resMas+mDeut;
841#ifdef debug
843#endif
844 if(sM<totM+.001) return 90001001;
845 }
846 if(NQ>1&&PQ>1&&AQ>4) // ===> "Can try to split an alpha" case
847 {
852 G4double sM=resMas+mAlph;
853 if(NQ!=4||PQ!=4) sM+=CB;
854#ifdef debug
856#endif
857 if(sM<totM+.001) return 90002002;
858 }
859 return 0;
860}
Referenced by EvaporateNucleus(), and G4QEnvironment::G4QEnvironment().
◆ StartLoop()
| G4bool G4QNucleus::StartLoop | ( | ) |
Definition at line 3982 of file G4QNucleus.cc.
3983{
3984 G4int theA=theNucleons.size();
3985 if(theA) currentNucleon=0;
3987 return theA;
3988} // End of StartLoop
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ SubtractNucleon()
| void G4QNucleus::SubtractNucleon | ( | G4QHadron * | pNucleon | ) |
Definition at line 612 of file G4QNucleus.cc.
613{
615 G4QHadronVector::iterator u; // iterator of the used nucleon
616 for(u=theNucleons.begin(); u!=theNucleons.end(); u++)
617 {
618#ifdef debug
620#endif
621 if (uNuc==*u) // Find uNuceon-pointer
622 {
623 NotFound=false;
624 break;
625 }
626 }
627// if(NotFound) throw G4QException("G4QNucleus::SubtractNucleon: The nucleon isn't found");
630 else
631 {
634#ifdef debug
636 <<theNucleons.size()<<G4endl;
637#endif
640#ifdef debug
642#endif
643 delete *u; // Delete the nucleon as an object
644 theNucleons.erase(u); // exclude the nucleon pointer from the HV
645 --currentNucleon; // Continue selection from theSame position
646 t4M-=u4M; // Update the nucleus 4-momentum VALUE
647 if (uPDG==2212) tPDG-=1000; // Reduce the nucleus PDG Code by a proton
648 else if(uPDG==2112) tPDG--; // Reduce the nucleus PDG Code by a neutron
649 else
650 {
651 // G4cerr<<"***G4QNucleus::SubtractNucleon: Unexpected Nucleon PDGCode ="<<uPDG<<G4endl;
652 // throw G4QException("G4QNucleus::SubtractNucleon: Impossible nucleon PDG Code");
653 G4ExceptionDescription ed;
654 ed << "Impossible nucleon PDG Code: Unexpected Nucleon PDGCode ="
655 << uPDG << G4endl;
657 FatalException, ed);
658 }
659#ifdef debug
662#endif
665 //#ifdef debug
668#ifdef debug
671#endif
673 //#endif
677 for(u=theNucleons.begin(); u!=theNucleons.end(); u++) // Correct the nucleon's energies
678 {
682 if((*u)->GetPDGCode()==2212) m2_value=m2p;// change it to subResNucleusM2 for protons
684#ifdef debug
686#endif
688 (*u)->Set4Momentum(n4M); // Update the 4-momentum of the nucleon
689 }
690 }
691#ifdef debug
693#endif
694}
double m2() const
double mag2() const
Referenced by G4QFragmentation::G4QFragmentation(), and G4QIonIonCollision::G4QIonIonCollision().
◆ UpdateClusters()
Definition at line 417 of file G4QNucleus.cc.
418{
419 //static const G4double r0 = 1.1; // fm, for nuclear radius: r=r0*A^(1/3)
420 //static const G4double del= .55; // fm, for a difused surface of the nucleus
421 //static const G4double rCl= 2.0; // clusterization radius @@??
422 //static const G4double freeibuc = 0.10; // probab. of the quasi-free baryon on surface
423 //static const G4double freeDib = 0.05; // probab. of the quasi-free dibar. on surface
424 //static const G4double clustProb = 4.0; // clusterization probability in dense region
425 //static const G4double prQ = 1.0; // relative probability for a Quasmon
426 //static const G4double prQ = 0.; //@@for pi@@relative probability for Quasmon
427 //G4double probSInt[254]; // integratedStaticProbabilities @@ not used
428 probVect[0]=mediRatio;
430 //probSInt[0]=0; // integrated static probabilities
431 dZ=0;
432 dN=0;
433 dS=0;
435#ifdef debug
437#endif
438 G4double A=a;
439 if(A<=0.)
440 {
441#ifdef debug
443#endif
444 return 0;
445 }
450#ifdef debug
452#endif
454 if (din && dA<2 && a>2)
455 {
456 dA=2;
457 sA=a-2;
458 }
459#ifdef debug
461#endif
463 G4double pA=0.;
464 G4double uA=0.;
465 if(surf>0.)
466 {
467 pA=0.5*freeDib*sA/surf; //@@Randomize(?)// a#of quasi-free Nucleon Pairs on the surface
468 uA=sA-pA-pA; // a#of quasi-free nucleons on Nuclear Surface
469 }
470 uA=uA/A; // Normalization of probability
471 pA=pA/A;
472 G4double sum =0.;
473 if(dA<2) // There is no dense phase at all
474 {
475 //probVect[1]= dA/A; // a#of quasi-free nucleons (only dense)
476 //probVect[1]= (uA+dA)/A; // a#of quasi-free nucleons (different norm)
477 probVect[1]= uA+dA/A; // a#of quasi-free nucleons (correct)
478 sum = probVect[1];
479 //probSInt[1]=sum; // integrated static probabilities
480 maxi=2;
481 probVect[254]= 0; // a#of dense nucleons (correct)
482 if(A>1 && pA>0.)
483 {
484 //probVect[2]= (pA+pA)/A/(A-1); // a#of quasi-free "dibaryons" (correct)
485 probVect[2]= pA; // a#of quasi-free "dibaryons" (correct)
486 //probVect[2]= 0; // a#of quasi-free "dibaryons" (only dense)
487 sum+= probVect[2]+probVect[2];
488 //probSInt[2]=sum; // integrated static probabilities
489 maxi=3;
490 probVect[255]= 0; // a#of dense "dibaryons" (correct)
491 }
492#ifdef debug
494#endif
495 }
496 else
497 {
500 // dA=C*Sum_k=1-A[n*C^A_k*wrd^(k-1)]=C*dA*(1+wrd)^(dA-1) => C=1/sud, sud=(1+wrd)^(dA-1)
501 // =1
502 G4double rd= dA/sud/A;
503 //G4double comb=A;
504 //G4double prb=rd; // (only dense)
505 G4double prb=rd+uA;
506 sum =prb;
507#ifdef debug
509#endif
510 //probVect[1]= prb/comb; // a#of quasi-free nucleons (correct)
511 //probVect[254]= rd/comb; // a#of dense nucleons (correct)
512 probVect[1]= prb; // a#of quasi-free nucleons (correct)
513 probVect[254]= rd; // a#of dense nucleons (correct)
514 //probSInt[1]=sum; // integrated static probabilities
515 // =2
516 rd*=wrd*(dA-1.)/2;
517 //comb*=(A-1.)/2;
518 //prb=rd; // (only dense)
519 prb=rd+pA;
520 sum+=prb+prb;
521#ifdef debug
523#endif
524 //probVect[2]= prb/comb; // a#of quasi-free "dibaryons" (correct)
525 //probVect[255]= rd/comb; // a#of dense "dibaryons" (correct)
526 probVect[2]= prb; // a#of quasi-free "dibaryons" (correct)
527 probVect[255]= rd; // a#of dense "dibaryons" (correct)
528 //probSInt[2]=sum; // integrated static probabilities
529 // >2
530 maxi=3;
531#ifdef debug
534#endif
535 if(dA>2)
536 {
537 ///////////G4double itA=A+1.;
538 G4double idA=dA+1.;
539 G4int dLim=dA;
540 if(maxClust<dA) dLim=maxClust;
541 for (int i=3; i<=dLim; i++)
542 {
543 rd*=wrd*(idA-i)/i;
544 sum+=rd*i;
545#ifdef debug
547#endif
548 //comb*=(itA-i)/i;
549 //probVect[i]=rd/comb; // Divide by sum of combinations for N+Z+S
550 probVect[i]=rd; // Comb's for N,Z,S are canceled later(G4QNuc)
551 //probSInt[i]=sum; // integrated static probabilities
552 maxi=i+1;
553#ifdef debug
555#endif
556 }
557 }
558 dS = S; // @@ Lambdas are always in the dense region
559 dZ = static_cast<int>(static_cast<double>((dA-dS)*Z)/(Z+N) + 0.5);
560 dN = dA - dZ;
561 }
562#ifdef debug
564#endif
565 maxClust=maxi-1;
566 //for (G4int j=maxi; j<255; j++) probVect[j]=0.;//Make the rest to be 0 [preinited above]
567 // =----------------= From here probability randomization starts =---------------=
568 // G4int rA=a; // Residual number of nucleons
569 //#ifdef debug
570 //G4cout<<"G4QNuc::UpdateClust:A="<<A<<",M="<<k<<",P1="<<probVect[1]<<",P2="<<probVect[2]
571 // <<G4endl;
572 //#endif
573 //if (k>1) for (j=k; j>1; j--) // nucleons are not randomized
574 //{
575 // G4int jmax=rA/j; // Max number of this kind of clusters
576 // if (jmax)
577 // {
578 // G4double prob=probVect[j]/probSInt[j]; // Probab of the cluster in the dest nucleus
579 //#ifdef debug
580 // G4cout<<"G4QNucl::UpdateClusters: j="<<j<<",sP="<<probVect[j]<<",iP="<<probSInt[j]
581 // <<G4endl;
582 //#endif
583 // G4int m=RandomizeBinom(prob,jmax); // A#of clusters of this type
584 // if(m)
585 // {
586 // probVect[j]=m;
587 // rA-=m*j;
588 // }
589 // else
590 // {
591 // probVect[j]=0.;
592 // if(j==maxClust) maxClust--;
593 // }
594 //#ifdef debug
595 // G4cout<<"G4QNucl::UpdateClust:p="<<prob<<",r="<<rA<<",m="<<jmax<<",P="<<probVect[j]
596 // <<G4endl;
597 //#endif
598 // }
599 // else
600 // {
601 // probVect[j]=0.;
602 // if(j==maxClust) maxClust--;
603 // }
604 //}
605 //probVect[1]=rA;
606 // =------------------= From here probability randomization starts =-------------------=
607 return maxClust;
608}
The documentation for this class was generated from the following files:
