G4QDiffractionRatio Class Reference

#include <G4QDiffractionRatio.hh>


Public Member Functions

 ~G4QDiffractionRatio ()
G4double GetRatio (G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
G4QHadronVectorProjFragment (G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
G4QHadronVectorTargFragment (G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
G4double GetTargSingDiffXS (G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)

Static Public Member Functions

static G4QDiffractionRatioGetPointer ()

Protected Member Functions

 G4QDiffractionRatio ()


Detailed Description

Definition at line 53 of file G4QDiffractionRatio.hh.


Constructor & Destructor Documentation

G4QDiffractionRatio::G4QDiffractionRatio (  )  [inline, protected]

Definition at line 57 of file G4QDiffractionRatio.hh.

00057 {}                 // Constructor

G4QDiffractionRatio::~G4QDiffractionRatio (  )  [inline]

Definition at line 61 of file G4QDiffractionRatio.hh.

00061 {}                 // Destructor


Member Function Documentation

G4QDiffractionRatio * G4QDiffractionRatio::GetPointer (  )  [static]

Definition at line 48 of file G4QDiffractionRatio.cc.

Referenced by G4ProjectileDiffractiveChannel::G4ProjectileDiffractiveChannel(), G4QFragmentation::G4QFragmentation(), and G4QDiffraction::PostStepDoIt().

00049 {
00050   static G4QDiffractionRatio theRatios;   // *** Static body of the Diffraction Ratio ***
00051   return &theRatios;
00052 }

G4double G4QDiffractionRatio::GetRatio ( G4double  pIU,
G4int  prPDG,
G4int  tgZ,
G4int  tgN 
)

Definition at line 55 of file G4QDiffractionRatio.cc.

References G4cout, G4endl, and CLHEP::detail::n.

Referenced by G4QFragmentation::G4QFragmentation(), and G4ProjectileDiffractiveChannel::GetFraction().

00056 {
00057   static const G4double mNeut= G4QPDGCode(2112).GetMass()/GeV; // in GeV
00058   static const G4double mProt= G4QPDGCode(2212).GetMass()/GeV; // in GeV
00059   static const G4double mN=.5*(mNeut+mProt);  // mean nucleon mass in GeV
00060   static const G4double dmN=mN+mN;            // doubled nuc. mass in GeV
00061   static const G4double mN2=mN*mN;            // squared nuc. mass in GeV^2
00062   // Table parameters
00063   static const G4int    nps=100;        // Number of steps in the R(s) LinTable
00064   static const G4int    mps=nps+1;      // Number of elements in the R(s) LinTable
00065   static const G4double sma=6.;         // The first LinTabEl(sqs=0)=1., sqs>sma -> logTab
00066   static const G4double ds=sma/nps;     // Step of the linear Table
00067   static const G4int    nls=150;        // Number of steps in the R(lns) logTable
00068   static const G4int    mls=nls+1;      // Number of elements in the R(lns) logTable
00069   static const G4double lsi=1.79;       // The min ln(sqs) logTabEl(sqs=5.99 < sma=6.)
00070   static const G4double lsa=8.;         // The max ln(sqs) logTabEl(sqs=5.99 - 2981 GeV)
00071   static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 5.99 GeV)
00072   static const G4double max_s=std::exp(lsa);// The max s of logTabEl(~ 2981 GeV)
00073   static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table
00074   static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
00075   static const G4double toler=.0001;    // Tolarence (GeV) defining the same sqs
00076   static G4double lastS=0.;             // Last sqs value for which R was calculated
00077   static G4double lastR=0.;             // Last ratio R which was calculated
00078   // Local Associative Data Base:
00079   static std::vector<G4int>     vA;     // Vector of calculated A
00080   //static std::vector<G4double>  vH;     // Vector of max sqs initialized in the LinTable
00081   //static std::vector<G4int>     vN;     // Vector of topBin number initialized in LinTab
00082   //static std::vector<G4double>  vM;     // Vector of relMax ln(sqs) initialized in LogTab
00083   //static std::vector<G4int>     vK;     // Vector of topBin number initialized in LogTab
00084   static std::vector<G4double*> vT;     // Vector of pointers to LinTable in C++ heap
00085   static std::vector<G4double*> vL;     // Vector of pointers to LogTable in C++ heap
00086   // Last values of the Associative Data Base:
00087   //static G4int     lastPDG=0;           // Last PDG for which R was calculated (now fake)
00088   static G4int     lastA=0;             // theLast of calculated A
00089   //static G4double  lastH=0.;            // theLast max sqs initialized in the LinTable
00090   //static G4int     lastN=0;             // theLast of topBin number initialized in LinTab
00091   //static G4double  lastM=0.;            // theLast relMax ln(sqs) initialized in LogTab
00092   //static G4int     lastK=0;             // theLast of topBin number initialized in LogTab
00093   static G4double* lastT=0;             // theLast of pointer to LinTable in the C++ heap
00094   static G4double* lastL=0;             // theLast of pointer to LogTable in the C++ heap
00095   // LogTable is created only if necessary. R(sqs>2981GeV) calcul by formula for any nuclei
00096   G4int A=tgN+tgZ;
00097   if(pIU<toler || A<1) return 1.;       // Fake use of toler as non zero number
00098   if(A>238)
00099   {
00100     G4cout<<"-*-Warning-*-G4QuasiFreeRatio::GetRatio: A="<<A<<">238, return zero"<<G4endl;
00101     return 0.;
00102   }
00103   //lastPDG=pPDG;                         // @@ at present ratio is PDG independent @@
00104   // Calculate sqs
00105   G4double pM=G4QPDGCode(pPDG).GetMass()/GeV; // Projectile mass in GeV
00106   G4double pM2=pM*pM;
00107   G4double mom=pIU/GeV;                 // Projectile momentum in GeV
00108   G4double s_value=std::sqrt(mN2+pM2+dmN*std::sqrt(pM2+mom*mom)); // in GeV
00109   G4int nDB=vA.size();                  // A number of nuclei already initialized in AMDB
00110   if(nDB && lastA==A && std::fabs(s_value-lastS)<toler) return lastR;
00111   if(s_value>max_s)
00112   {
00113     lastR=CalcDiff2Prod_Ratio(s_value,A);     // @@ Probably user ought to be notified about bigS
00114     return lastR;
00115   }
00116   G4bool found=false;
00117   G4int i=-1;
00118   if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB
00119   {
00120     found=true;                         // The A value is found
00121     break;
00122   }
00123   if(!nDB || !found)                    // Create new line in the AMDB
00124   {
00125     lastA = A;
00126     lastT = new G4double[mps];          // Create the linear Table
00127     //lastN = static_cast<int>(s_value/ds)+1;   // MaxBin to be initialized
00128     //if(lastN>nps)                     // ===> Now initialize all lin table
00129     //{
00130     //  lastN=nps;
00131     //  lastH=sma;
00132     //}
00133     //else lastH = lastN*ds;              // Calculate max initialized s for LinTab
00134     G4double sv=0;
00135     lastT[0]=1.;
00136     //for(G4int j=1; j<=lastN; j++)       // Calculate LinTab values
00137     for(G4int j=1; j<=nps; j++)       // Calculate LinTab values
00138     {
00139       sv+=ds;
00140       lastT[j]=CalcDiff2Prod_Ratio(sv,A);
00141     }
00142     lastL = new G4double[mls];          // Create the logarithmic Table
00143     //G4double ls=std::log(s_value);
00144     //lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
00145     //if(lastK>nls)                     // ===> Now initialize all lin table
00146     //{
00147     //  lastK=nls;
00148     //  lastM=lsa-lsi;
00149     //}
00150     //else lastM = lastK*dl;              // Calculate max initialized ln(s)-lsi for LogTab
00151     sv=mi;
00152     //for(G4int j=0; j<=lastK; j++)     // Calculate LogTab values
00153     for(G4int j=0; j<=nls; j++)     // Calculate LogTab values
00154     {
00155       lastL[j]=CalcDiff2Prod_Ratio(sv,A);
00156       //if(j!=lastK) sv*=edl;
00157       sv*=edl;
00158     }
00159     i++;                                // Make a new record to AMDB and position on it
00160     vA.push_back(lastA);
00161     //vH.push_back(lastH);
00162     //vN.push_back(lastN);
00163     //vM.push_back(lastM);
00164     //vK.push_back(lastK);
00165     vT.push_back(lastT);
00166     vL.push_back(lastL);
00167   }
00168   else                                  // The A value was found in AMDB
00169   {
00170     lastA=vA[i];
00171     //lastH=vH[i];
00172     //lastN=vN[i];
00173     //lastM=vM[i];
00174     //lastK=vK[i];
00175     lastT=vT[i];
00176     lastL=vL[i];
00177     // ==> Now all bins of the tables are initialized immediately for the A
00178     //if(s_value>lastH)                    // At least LinTab must be updated
00179     //{
00180     //  G4int nextN=lastN+1;               // The next bin to be initialized
00181     //  if(lastN<nps)
00182     //  {
00183     //    lastN = static_cast<int>(s_value/ds)+1;// MaxBin to be initialized
00184     //    G4double sv=lastH;
00185     //    if(lastN>nps)
00186     //    {
00187     //      lastN=nps;
00188     //      lastH=sma;
00189     //    }
00190     //    else lastH = lastN*ds;           // Calculate max initialized s for LinTab
00191     //    for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values
00192     //    {
00193     //      sv+=ds;
00194     //      lastT[j]=CalcDiff2Prod_Ratio(sv,A);
00195     //    }
00196     //  } // End of LinTab update
00197     //  if(lastN>=nextN)
00198     //  {
00199     //    vH[i]=lastH;
00200     //    vN[i]=lastN;
00201     //  }
00202     //  G4int nextK=lastK+1;
00203     //  if(s_value>sma && lastK<nls)             // LogTab must be updated
00204     //  {
00205     //    G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed)
00206     //    G4double ls=std::log(s_value);
00207     //    lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
00208     //    if(lastK>nls)
00209     //    {
00210     //      lastK=nls;
00211     //      lastM=lsa-lsi;
00212     //    }
00213     //    else lastM = lastK*dl;           // Calcul. max initialized ln(s)-lsi for LogTab
00214     //    for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
00215     //    {
00216     //      sv*=edl;
00217     //      lastL[j]=CalcDiff2Prod_Ratio(sv,A);
00218     //    }
00219     //  } // End of LogTab update
00220     //  if(lastK>=nextK)
00221     //  {
00222     //    vM[i]=lastM;
00223     //    vK[i]=lastK;
00224     //  }
00225     //}
00226   }
00227   // Now one can use tabeles to calculate the value
00228   if(s_value<sma)                             // Use linear table
00229   {
00230     G4int n=static_cast<int>(s_value/ds);     // Low edge number of the bin
00231     G4double d=s_value-n*ds;                  // Linear shift
00232     G4double v=lastT[n];                      // Base
00233     lastR=v+d*(lastT[n+1]-v)/ds;              // Result
00234   }
00235   else                                  // Use log table
00236   {
00237     G4double ls=std::log(s_value)-lsi;  // ln(s)-l_min
00238     G4int n=static_cast<int>(ls/dl);    // Low edge number of the bin
00239     G4double d=ls-n*dl;                 // Log shift
00240     G4double v=lastL[n];                // Base
00241     lastR=v+d*(lastL[n+1]-v)/dl;        // Result
00242   }
00243   if(lastR<0.) lastR=0.;
00244   if(lastR>1.) lastR=1.;
00245   return lastR;
00246 } // End of GetRatio

G4double G4QDiffractionRatio::GetTargSingDiffXS ( G4double  pIU,
G4int  prPDG,
G4int  tgZ,
G4int  tgN 
)

Definition at line 1290 of file G4QDiffractionRatio.cc.

References G4cerr, and G4endl.

01291 {
01292   G4double mom=pIU/gigaelectronvolt;    // Projectile momentum in GeV
01293   if ( mom < 1. || (pPDG != 2212 && pPDG != 2112) )
01294     G4cerr<<"G4QDiffractionRatio::GetTargSingDiffXS isn't applicable p="<<mom<<" GeV, PDG="
01295          <<pPDG<<G4endl;
01296   G4double A=Z+N;                        // A of the target
01297   //return 4.5*std::pow(A,.364)*millibarn; // Result
01298   return 3.7*std::pow(A,.364)*millibarn; // Result after mpi0 correction
01299 
01300 } // End of ProjFragment

G4QHadronVector * G4QDiffractionRatio::ProjFragment ( G4int  pPDG,
G4LorentzVector  p4M,
G4int  tgZ,
G4int  tgN 
)

Definition at line 484 of file G4QDiffractionRatio.cc.

References FatalException, G4Quasmon::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4RandomDirection(), G4UniformRand, and G4QNucleus::GetMZNS().

Referenced by G4QFragmentation::G4QFragmentation(), and G4ProjectileDiffractiveChannel::Scatter().

00486 {
00487   static const G4double pFm= 250.; // Fermi momentum in MeV (delta function)
00488   static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function)
00489   static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction)
00490   static const G4double mPi= G4QPDGCode(211).GetMass();  // pi+- mass (MeV)
00491   static const G4double mNeut= G4QPDGCode(2112).GetMass();
00492   static const G4double mNeut2=mNeut*mNeut;
00493   static const G4double dmNeut=mNeut+mNeut;
00494   static const G4double mProt= G4QPDGCode(2212).GetMass();
00495   static const G4double mProt2=mProt*mProt;
00496   static const G4double dmProt=mProt+mProt;
00497   static const G4double maxDM=mProt*12.;
00498   static const G4double mLamb= G4QPDGCode(3122).GetMass();
00499   static const G4double mSigZ= G4QPDGCode(3212).GetMass();
00500   static const G4double mSigM= G4QPDGCode(3112).GetMass();
00501   static const G4double mSigP= G4QPDGCode(3222).GetMass();
00502   static const G4double eps=.003;
00503   //
00504   G4LorentzVector pr4M=p4M/megaelectronvolt;   // Convert 4-momenta in MeV (keep p4M)
00505   // prepare the DONOTHING answer
00506   G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !!
00507   G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile
00508   ResHV->push_back(hadron);                // It must be cleaned up for real scattering sec
00509   // @@ diffraction is simulated as noncoherent (coherent is small)
00510   G4int tgA=tgZ+tgN;                       // A of the target
00511   G4int tPDG=90000000+tgZ*1000+tgN;        // PDG code of the targetNucleus/recoilNucleus
00512   G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus
00513   G4int rPDG=2112;                         // prototype of PDG code of the recoiled nucleon
00514   if(tgA*G4UniformRand()>tgN)              // Substitute by a proton
00515   {
00516     rPDG=2212;                             // PDG code of the recoiled QF nucleon
00517     tPDG-=1000;                            // PDG code of the recoiled nucleus
00518   }
00519   else tPDG-=1;                            // PDG code of the recoiled nucleus
00520   G4double tM=G4QPDGCode(tPDG).GetMass();  // Mass of the recoiled nucleus
00521   G4double tE=std::sqrt(tM*tM+pFm2);
00522   G4ThreeVector tP=pFm*G4RandomDirection();
00523   G4LorentzVector t4M(tP,tE);              // 4M of the recoil nucleus
00524   G4LorentzVector tg4M(0.,0.,0.,tgM);
00525   G4LorentzVector N4M=tg4M-t4M;            // Quasi-free target nucleon
00526   G4LorentzVector tot4M=N4M+p4M;           // total momentum of quasi-free diffraction
00527   G4double mT=mNeut;
00528   G4double mT2=mNeut2;                 // Squared mass of the free nucleon spectator
00529   G4double dmT=dmNeut;
00530   //G4int Z=0;
00531   //G4int N=1;
00532   if(rPDG==2212)
00533   {
00534     mT=mProt;
00535     mT2=mProt2;
00536     dmT=dmProt;
00537     //Z=1;
00538     //N=0;
00539   }
00540   G4double mP2=pr4M.m2();               // Squared mass of the projectile
00541   if(mP2<0.) mP2=0.;                    // A possible problem for photon (m_min = 2*m_pi0)
00542   G4double s_value=tot4M.m2();          // @@ Check <0 ...
00543   G4double E=(s_value-mT2-mP2)/dmT;     // Effective interactin energy (virt. nucl. target)
00544   G4double E2=E*E;
00545   if(E<0. || E2<mP2)
00546   {
00547 #ifdef pdebug
00548     G4cerr<<"-Warning-G4DifR::PFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
00549 #endif
00550     return ResHV; // Do Nothing Action
00551   }
00552   G4double mP=std::sqrt(mP2);
00553   if(mP<.1)mP=mPi0;                      // For photons min diffraction is gamma+P->Pi0+Pi0
00554   G4double mMin=mP+mPi0;                 // Minimum diffractive mass
00555   G4double ss=std::sqrt(s_value);        // CM compound mass (sqrt(s))
00556   G4double mMax=ss-mT;                   // Maximum diffraction mass
00557   if(mMax>maxDM) mMax=maxDM;             // Restriction to avoid too big masses
00558   if(mMin>=mMax)
00559   {
00560 #ifdef pdebug
00561     G4cerr<<"-Warning-G4DifR::PFra:ZeroDiffractionMRange, mi="<<mMin<<",ma="<<mMax<<G4endl;
00562 #endif
00563     return ResHV; // Do Nothing Action
00564   }
00565   G4double R = G4UniformRand();
00566   G4double mDif=std::exp(R*std::log(mMax)+(1.-R)*std::log(mMin)); // LowMassApproximation
00567   G4double mDif2=mDif*mDif;
00568   G4double ds=s_value-mT2-mDif2;
00569   //G4double e=ds/dmT;
00570   //G4double P=std::sqrt(e*e-mDif2);          // Momentum in pseudo laboratory system
00571 #ifdef debug
00572   G4cout<<"G4QDiffR::PFra: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
00573 #endif
00574   // @@ Temporary NN t-dependence for all hadrons
00575   if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl;
00576   G4double tsl=140000.;                        // slope in MeV^2 
00577   G4double t=-std::log(G4UniformRand())*tsl;
00578   G4double maxt=(ds*ds-4*mT2*mDif2)/s_value;   // maximum possible -t
00579 #ifdef pdebug
00580   G4cout<<"G4QDifR::PFra:ph="<<pPDG<<",P="<<P<<",t="<<mint<<"<"<<maxt<<G4endl;
00581 #endif
00582 #ifdef nandebug
00583   if(mint>-.0000001);                          // To make the Warning for NAN
00584   else  G4cout<<"******G4QDiffractionRatio::ProjFragment: -t="<<mint<<G4endl;
00585 #endif
00586   G4double rt=t/maxt;
00587   G4double cost=1.-rt-rt;                      // cos(theta) in CMS
00588 #ifdef ppdebug
00589   G4cout<<"G4QDiffRatio::ProjFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
00590 #endif
00591   if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00592   {
00593     if     (cost>1.)  cost=1.;
00594     else if(cost<-1.) cost=-1.;
00595     else
00596     {
00597       G4cerr<<"G4QDiffRat::ProjFragm: *NAN* cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl;
00598       return ResHV; // Do Nothing Action
00599     }
00600   }
00601   G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
00602   G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif);    // 4mom of the diffract. Quasmon
00603   G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00604   if(!G4QHadron(tot4M).RelDecayIn2(d4M, r4M, dir4M, cost, cost))
00605   {
00606     G4cerr<<"G4QDifR::PFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl;
00607     //G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp");
00608     return ResHV; // Do Nothing Action
00609   }
00610 #ifdef debug
00611   G4cout<<"G4QDiffR::ProjFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl;
00612 #endif
00613   // Now everything is ready for fragmentation and DoNothing projHadron must be wiped out
00614   delete hadron;     // Delete the fake (doNothing) projectile hadron
00615   ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile
00616   hadron = new G4QHadron(tPDG,t4M);  // Hadron for the recoil neucleus
00617   ResHV->push_back(hadron);          // Fill the recoil nucleus
00618 #ifdef debug
00619   G4cout<<"G4QDiffractionRatio::ProjFragment: *Filled* RecNucleus="<<t4M<<tPDG<<G4endl;
00620 #endif
00621   hadron = new G4QHadron(rPDG,r4M);  // Hadron for the recoil nucleon
00622   ResHV->push_back(hadron);          // Fill the recoil nucleon
00623 #ifdef debug
00624   G4cout<<"G4QDiffractionRatio::ProjFragment: *Filled* RecNucleon="<<r4M<<rPDG<<G4endl;
00625 #endif
00626   G4LorentzVector sum4M(0.,0.,0.,0.);
00627   // Now the (pPdg,d4M) Quasmon must be fragmented
00628   G4QHadronVector* leadhs = 0;       // Prototype of QuasmOutput G4QHadronVector
00629   G4QContent dQC=G4QPDGCode(pPDG).GetQuarkContent(); // Quark Content of the projectile
00630   G4Quasmon* pan= new G4Quasmon(dQC,d4M); // --->---->---->----->-----> DELETED -->---*
00631   try                                                           //                    |
00632   {                                                             //                    |
00633     G4QNucleus vac(90000000);                                   //                    |
00634     leadhs=pan->Fragment(vac,1);  // DELETED after it is copied to ResHV vector -->---+-*
00635   }                                                             //                    | |
00636   catch (G4QException& error)                                   //                    | |
00637   {                                                             //                    | |
00638     G4cerr<<"***G4QDiffractionRatio::ProjFragment: G4Quasmon Exception"<<G4endl;    //| |
00639     // G4Exception("G4QDiffractionRatio::ProjFragment","72",FatalException,"*Quasmon");//| |
00640     G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0072",
00641                 FatalException, "*Quasmon");
00642   }                                                             //                    | |
00643   delete pan;                              // Delete the Nuclear Environment <----<---* |
00644   G4int qNH=leadhs->size();                // A#of collected hadrons from diff.frag.    |
00645   if(qNH) for(G4int iq=0; iq<qNH; iq++)    // Loop over hadrons to fill the result      |
00646   {                                        //                                           |
00647     G4QHadron* loh=(*leadhs)[iq];          // Pointer to the output hadron              |
00648     G4int nL=loh->GetStrangeness();        // A number of Lambdas in the Hypernucleus   |
00649     G4int nB=loh->GetBaryonNumber();       // Total Baryon Number of the Hypernucleus   |
00650     G4int nC = loh->GetCharge();           // Charge of the Hypernucleus                |
00651     G4int oPDG = loh->GetPDGCode();        // Original CHIPS PDG Code of the hadron     |
00652     //if((nC>nB || nC<0) && nB>0 && nL>=0 && nL<=nB && oPDG>80000000) // Iso-nucleus    |
00653     if(2>3) // Closed because "G4QDR::F:90002999,M=-7.768507e-04,B=2,S=0,C=3" is found  |
00654     {
00655       G4LorentzVector q4M = loh->Get4Momentum(); // Get 4-momentum of the Isonucleus    |
00656       G4double qM=q4M.m();                 // Real mass of the Isonucleus
00657 #ifdef fdebug
00658       G4cout<<"G4QDR::PF:"<<oPDG<<",M="<<qM<<",B="<<nB<<",S="<<nL<<",C="<<nC<<G4endl;// |
00659 #endif
00660       G4int    qPN=nC-nB;                  // Number of pions in the Isonucleus         |
00661       G4int    fPDG = 2212;                // Prototype for nP+(Pi+) case               |
00662       G4int    sPDG = 211;
00663       tPDG = 3122;                         // @@ Sigma0 (?)                             |
00664       G4double fMass= mProt;
00665       G4double sMass= mPi;
00666       G4double tMass= mLamb;               // @@ Sigma0 (?)                             |
00667       G4bool   cont=true;                  // Continue flag                             |
00668       // =--------= Negative state =---------=
00669       if(nC<0)                             // =----= Only Pi- can help                  |
00670       {
00671         if(nL&&nB==nL)                     // --- n*Lamb + k*(Pi-) State ---            |
00672         {
00673           sPDG = -211;
00674           if(-nC==nL && nL==1)             // Only one Sigma- like (nB=1)               |
00675           {
00676             if(std::fabs(qM-mSigM)<eps)
00677             {
00678               loh->SetQPDG(G4QPDGCode(3112));  // This is Sigma-                        |
00679               cont=false;                  // Skip decay                                |
00680             }
00681             else if(qM>mLamb+mPi)          //(2) Sigma- => Lambda + Pi- decay           |
00682             {
00683               fPDG = 3122;
00684               fMass= mLamb;
00685             }
00686             else if(qM>mSigM)              //(2) Sigma+=>Sigma++gamma decay             |
00687             {
00688               fPDG = 3112;
00689               fMass= mSigM;
00690               sPDG = 22;
00691               sMass= 0.;
00692             }
00693             else                           //(2) Sigma-=>Neutron+Pi- decay              |
00694             {
00695               fPDG = 2112;
00696               fMass= mNeut;
00697             }
00698             qPN  = 1;                      // #of (Pi+ or gamma)'s = 1                  |
00699           }
00700           else if(-nC==nL)                 //(2) a few Sigma- like                      |
00701           {
00702             qPN  = 1;                      // One separated Sigma-                      |
00703             fPDG = 3112;
00704             sPDG = 3112;
00705             sMass= mSigM;
00706             nB--;
00707             fMass= mSigM;
00708           }
00709           else if(-nC>nL)                  //(2) n*(Sigma-)+m*(Pi-)                     |
00710           {
00711             qPN  = -nC-nL;                 // #of Pi-'s                                 |
00712             fPDG = 3112;
00713             fMass= mSigM;
00714           }
00715           else                             //(2) n*(Sigma-)+m*Lambda(-nC<nL)            |
00716           {
00717             nB += nC;                      // #of Lambda's                              |
00718             fPDG = 3122;
00719             fMass= mLamb;
00720             qPN  = -nC;                    // #of Sigma+'s                              |
00721             sPDG = 3112;
00722             sMass= mSigM;
00723           }
00724           nL   = 0;                        // Only decays in two are above              |
00725         }
00726         else if(nL)                        // ->n*Lamb+m*Neut+k*(Pi-) State (nL<nB)     |
00727         {
00728           nB -= nL;                        // #of neutrons                              |
00729           fPDG = 2112;
00730           fMass= mNeut;
00731           G4int nPin = -nC;                           // #of Pi-'s                    
00732           if(nL==nPin)                                //(2) m*Neut+n*Sigma-             |
00733           {
00734             qPN  = nL;                                // #of Sigma-                     |
00735             sPDG = 3112;
00736             sMass= mSigM;
00737             nL   = 0;
00738           }
00739           else if(nL>nPin)                            //(3) m*P+n*(Sigma+)+k*Lambda     |
00740           {
00741             nL-=nPin;                                 // #of Lambdas                    |
00742             qPN  = nPin;                              // #of Sigma+                     |
00743             sPDG = 3112;
00744             sMass= mSigM;
00745           }
00746           else                                 //(3) m*N+n*(Sigma-)+k*(Pi-) (nL<nPin)   |
00747           {
00748             qPN  = nPin-nL;                           // #of Pi-                        |
00749             sPDG = -211;
00750             tPDG = 3112;
00751             tMass= mSigM;
00752           }
00753         }
00754         else                                          //(2) n*N+m*(Pi-)   (nL=0)        |
00755         {
00756           sPDG = -211;
00757           qPN  = -nC;
00758           fPDG = 2112;
00759           fMass= mNeut;
00760         }
00761       }
00762       else if(!nC)                                   // *** Should not be here ***      |
00763       {
00764         if(nL && nL<nB)          //(2) n*Lamb+m*N ***Should not be here***              |
00765         {
00766           qPN  = nL;
00767           fPDG = 2112;                               // mN+nL case                      |
00768           sPDG = 3122;
00769           sMass= mLamb;
00770           nB -= nL;
00771           fMass= mNeut;
00772           nL   = 0;
00773         }
00774         else if(nL>1 && nB==nL)  //(2) m*Lamb(m>1) ***Should not be here***             |
00775         {
00776           qPN  = 1;
00777           fPDG = 3122;
00778           sPDG = 3122;
00779           sMass= mLamb;
00780           nB--;
00781           fMass= mLamb;
00782         }
00783         else if(!nL && nB>1)     //(2) n*Neut(n>1) ***Should not be here***             |
00784         {
00785           qPN  = 1;
00786           fPDG = 2112;
00787           sPDG = 2112;
00788           sMass= mNeut;
00789           nB--;
00790           fMass= mNeut;
00791         }
00792         else G4cout<<"*?*G4QDiffractionRatio::ProjFragment: (1) oPDG="<<oPDG<<G4endl;// |
00793       }
00794       else if(nC>0)              // n*Lamb+(m*P)+(k*Pi+)                                |
00795       {
00796         if(nL && nL+nC==nB)      //(2) n*Lamb+m*P ***Should not be here***              |
00797         {
00798           qPN  = nL;
00799           nL   = 0;
00800           fPDG = 2212;
00801           sPDG = 3122;
00802           sMass= mLamb;
00803           nB  = nC;
00804           fMass= mProt;
00805         }
00806         else if(nL  && nC<nB-nL) //(3)n*L+m*P+k*N ***Should not be here***              |
00807         {
00808           qPN  = nC;                                  // #of protons                    |
00809           fPDG = 2112;                                // mP+nL case                     |
00810           sPDG = 2212;
00811           sMass= mProt;
00812           nB -= nL+nC;                                // #of neutrons                   |
00813           fMass= mNeut;
00814         }
00815         else if(nL  && nB==nL)                        // ---> n*L+m*Pi+ State           |
00816         {
00817           if(nC==nL && nL==1)                         // Only one Sigma+ like State     |
00818           {
00819             if(std::fabs(qM-mSigP)<eps)
00820             {
00821               loh->SetQPDG(G4QPDGCode(3222));         // This is GS Sigma+              |
00822               cont=false;                  // Skip decay                                |
00823             }
00824             else if(qM>mLamb+mPi)                     //(2) Sigma+=>Lambda+Pi+ decay    |
00825             {
00826               fPDG = 3122;
00827               fMass= mLamb;
00828             }
00829             else if(qM>mNeut+mPi)                     //(2) Sigma+=>Neutron+Pi+ decay   |
00830             {
00831               fPDG = 2112;
00832               fMass= mNeut;
00833             }
00834             else if(qM>mSigP)                         //(2) Sigma+=>Sigma++gamma decay  |
00835             {
00836               fPDG = 3222;
00837               fMass= mSigP;
00838               sPDG = 22;
00839               sMass= 0.;
00840             }
00841             else                                      //(2) Sigma+=>Proton+gamma decay  |
00842             {
00843               fPDG = 2212;
00844               fMass= mProt;
00845               sPDG = 22;
00846               sMass= 0.;
00847             }
00848             qPN  = 1;                                 // #of (Pi+ or gamma)'s = 1       |
00849           }
00850           else if(nC==nL)                             //(2) a few Sigma+ like hyperons  |
00851           {
00852             qPN  = 1;
00853             fPDG = 3222;
00854             sPDG = 3222;
00855             sMass= mSigP;
00856             nB--;
00857             fMass= mSigP;
00858           }
00859           else if(nC>nL)                              //(2) n*(Sigma+)+m*(Pi+)          |
00860           {
00861             qPN  = nC-nL;                             // #of Pi+'s                      |
00862             fPDG = 3222;
00863             nB  = nL;                                 // #of Sigma+'s                   |
00864             fMass= mSigP;
00865           }
00866           else                                        //(2) n*(Sigma+)+m*Lambda         |
00867           {
00868             nB -= nC;                                 // #of Lambda's                   |
00869             fPDG = 3122;
00870             fMass= mLamb;
00871             qPN  = nC;                                // #of Sigma+'s                   |
00872             sPDG = 3222;
00873             sMass= mSigP;
00874           }
00875           nL   = 0;                                   // All above are decays in 2      |
00876         }
00877         else if(nL && nC>nB-nL)                       // n*Lamb+m*P+k*Pi+               |
00878         {
00879           nB -= nL;                                   // #of protons                    |
00880           G4int nPip = nC-nB;                         // #of Pi+'s                      |
00881           if(nL==nPip)                                //(2) m*P+n*Sigma+                |
00882           {
00883             qPN  = nL;                                // #of Sigma+                     |
00884             sPDG = 3222;
00885             sMass= mSigP;
00886             nL   = 0;
00887           }
00888           else if(nL>nPip)                            //(3) m*P+n*(Sigma+)+k*Lambda     |
00889           {
00890             nL  -= nPip;                              // #of Lambdas                    |
00891             qPN  = nPip;                              // #of Sigma+                     |
00892             sPDG = 3222;
00893             sMass= mSigP;
00894           }
00895           else                                        //(3) m*P+n*(Sigma+)+k*(Pi+)      |
00896           {
00897             qPN  = nPip-nL;                           // #of Pi+                        |
00898             tPDG = 3222;
00899             tMass= mSigP;
00900           }
00901         }
00902         if(nC<nB)                 //(2) n*P+m*N ***Should not be here***                |
00903         {
00904           fPDG = 2112;
00905           fMass= mNeut;
00906           qPN  = nC;
00907           sPDG = 2212;
00908           sMass= mProt;
00909         }
00910         else if(nB==nC && nC>1)   //(2) m*Prot(m>1) ***Should not be here***            |
00911         {
00912           qPN  = 1;
00913           fPDG = 2212;
00914           sPDG = 2212;
00915           sMass= mProt;
00916           nB--;
00917           fMass= mProt;
00918         }
00919         else if(nC<=nB||!nB) G4cout<<"*?*G4QDR::ProjFragm: (2) oPDG="<<oPDG<<G4endl; // |
00920         // !nL && nC>nB                             //(2) Default condition n*P+m*(Pi+) |
00921       }
00922       if(cont)                                      // Make a decay                     |
00923       {
00924         G4double tfM=nB*fMass;
00925         G4double tsM=qPN*sMass;
00926         G4double ttM=0.;
00927         if(nL) ttM=nL*tMass;
00928         G4LorentzVector f4Mom(0.,0.,0.,tfM);
00929         G4LorentzVector s4Mom(0.,0.,0.,tsM);
00930         G4LorentzVector t4Mom(0.,0.,0.,ttM);
00931         G4double sum=tfM+tsM+ttM;
00932         if(std::fabs(qM-sum)<eps)
00933         {
00934           f4Mom=q4M*(tfM/sum);
00935           s4Mom=q4M*(tsM/sum);
00936           if(nL) t4Mom=q4M*(ttM/sum);
00937         }
00938         else if(!nL && (qM<sum || !G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))) // Error     |
00939         {
00940           //#ifdef fdebug
00941           G4cout<<"***G4QDR::PrFragm:fPDG="<<fPDG<<"*"<<nB<<"(fM="<<fMass<<")+sPDG="<<sPDG
00942                 <<"*"<<qPN<<"(sM="<<sMass<<")"<<"="<<sum<<" > TM="<<qM<<q4M<<oPDG<<G4endl;
00943           //#endif
00944           // throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 2"); //  |
00945           G4ExceptionDescription ed;
00946           ed << "***G4QDR::PrFragm:fPDG=" << fPDG << "*" << nB << "(fM="
00947              << fMass << ")+sPDG=" << sPDG << "*" << qPN << "(sM=" << sMass
00948              << ")" << "=" << sum << " > TM=" << qM << q4M << oPDG << G4endl;
00949           G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0002",
00950                       FatalException, ed);
00951         }
00952         else if(nL && (qM<sum || !G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom)))// Error|
00953         {
00954           //#ifdef fdebug
00955           G4cout<<"***G4DF::PrFrag: "<<fPDG<<"*"<<nB<<"("<<fMass<<")+"<<sPDG<<"*"<<qPN<<"("
00956                 <<sMass<<")+Lamb*"<<nL<<"="<<sum<<" > TotM="<<qM<<q4M<<oPDG<<G4endl;
00957           //#endif
00958           // throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 3"); //  |
00959           G4ExceptionDescription ed;
00960           ed << "***G4DF::PrFrag: " << fPDG << "*" << nB << "(" << fMass << ")+"
00961              << sPDG << "*" << qPN << "(" << sMass << ")+Lamb*" << nL << "="
00962              << sum << " > TotM=" << qM << q4M << oPDG << G4endl;
00963           G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0003",
00964                       FatalException, ed);
00965         }
00966 #ifdef fdebug
00967         G4cout<<"G4QDF::ProjFragm: *DONE* n="<<nB<<f4Mom<<fPDG<<", m="<<qPN<<s4Mom<<sPDG
00968               <<", l="<<nL<<t4Mom<<G4endl;
00969 #endif
00970         G4bool notused=true;
00971         if(nB)                               // There are baryons                       |
00972         {
00973           f4Mom/=nB;
00974           loh->Set4Momentum(f4Mom);          // ! Update the Hadron !                   |
00975           loh->SetQPDG(G4QPDGCode(fPDG));    // Baryons                                 |
00976           notused=false;                     // Loh was used                            |
00977           if(nB>1) for(G4int ih=1; ih<nB; ih++) // Loop over the rest of baryons        |
00978           {
00979             G4QHadron* Hi = new G4QHadron(fPDG,f4Mom); // Create a Hadron for Baryon    |
00980             ResHV->push_back(Hi);            // Fill in the additional nucleon          |
00981 #ifdef fdebug
00982             sum4M+=r4M;                      // Sum 4-momenta for the EnMom check       |
00983             G4cout<<"G4QDR::ProjFrag: *additional Nucleon*="<<f4Mom<<fPDG<<G4endl; //   |
00984 #endif
00985           }
00986         }
00987         if(qPN)                              // There are pions                         |
00988         {
00989           s4Mom/=qPN;
00990           G4int min=0;
00991           if(notused)
00992           {
00993             loh->Set4Momentum(s4Mom);        // ! Update the Hadron 4M !                |
00994             loh->SetQPDG(G4QPDGCode(sPDG));  // Update PDG                              |
00995             notused=false;                   // loh was used                            |
00996             min=1;                           // start value                             |
00997           }
00998           if(qPN>min) for(G4int ip=min; ip<qPN; ip++) // Loop over pions                |
00999           {
01000             G4QHadron* Hj = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the meson |
01001             ResHV->push_back(Hj);            // Fill in the additional pion             |
01002 #ifdef fdebug
01003             sum4M+=r4M;                      // Sum 4-momenta for the EnMom check       |
01004             G4cout<<"G4QDR::ProjFragm: *additional Pion*="<<f4Mom<<fPDG<<G4endl; //     |
01005 #endif
01006           }
01007         }
01008         if(nL)                               // There are Hyperons                      |
01009         {
01010           t4Mom/=nL;
01011           G4int min=0;
01012           if(notused)
01013           {
01014             loh->Set4Momentum(t4Mom);      // ! Update the Hadron 4M !                  |
01015             loh->SetQPDG(G4QPDGCode(tPDG));// Update PDG                                |
01016             notused=false;                 // loh was used                              |
01017             min=1;                         //   
01018           }
01019           if(nL>min) for(G4int il=min; il<nL; il++) // Loop over Hyperons               |
01020           {
01021             G4QHadron* Hk = new G4QHadron(tPDG,t4Mom); // Create a Hadron for Lambda    |
01022             ResHV->push_back(Hk);          // Fill in the additional pion               |
01023 #ifdef fdebug
01024             sum4M+=r4M;                    // Sum 4-momenta for the EnMom check         |
01025             G4cout<<"G4QDR::ProjFragm: *additional Hyperon*="<<f4Mom<<fPDG<<G4endl; //  |
01026 #endif
01027           }
01028         }
01029       }                                    // --> End of decay                          |
01030     }                                      // -> End of Iso-nuclear treatment           |
01031     else if( (nL > 0 && nB > 1) || (nL < 0 && nB < -1) ) 
01032     {     // Hypernucleus is found                                                      |
01033       G4bool anti=false;                   // Default=Nucleus (true=antinucleus         |
01034       if(nB<0)                             // Anti-nucleus                              |
01035       {
01036         anti=true;                         // Flag of anti-hypernucleus                 |
01037         nB=-nB;                            // Reverse the baryon number                 |
01038         nC=-nC;                            // Reverse the charge                        |
01039         nL=-nL;                            // Reverse the strangeness                   |
01040       }
01041       G4int hPDG = 90000000+nL*999999+nC*999+nB; // CHIPS PDG Code for Hypernucleus     |
01042       G4int nSM=0;                         // A#0f unavoidable Sigma-                   |
01043       G4int nSP=0;                         // A#0f unavoidable Sigma+                   |
01044       if(nC<0)                             // Negative hypernucleus                     |
01045       {
01046         if(-nC<=nL)                        // Partial compensation by Sigma-            |
01047         {
01048           nSM=-nC;                         // Can be compensated by Sigma-              |
01049           nL+=nC;                          // Reduce the residual strangeness           |
01050         }
01051         else                               // All Charge is compensated by Sigma-       |
01052         {
01053           nSM=nL;                          // The maximum number of Sigma-              |
01054           nL=0;                            // Kill the residual strangeness             |
01055         }
01056       }
01057       else if(nC>nB-nL)                    // Extra positive hypernucleus               |
01058       {
01059         if(nC<=nB)                         // Partial compensation by Sigma+            |
01060         {
01061           G4int dH=nB-nC;                  // Isotopic shift                            |
01062           nSP=nL-dH;                       // Can be compensated by Sigma+              |
01063           nL=dH;                           // Reduce the residual strangeness           |
01064         }
01065         else                               // All Charge is compensated by Sigma+       |
01066         {
01067           nSP=nL;                          // The maximum number of Sigma+              |
01068           nL=0;                            // Kill the residual strangeness             |
01069         }
01070       }
01071       r4M=loh->Get4Momentum();                 // Real 4-momentum of the hypernucleus   !
01072       G4double reM=r4M.m();                    // Real mass of the hypernucleus         |
01073 #ifdef fdebug
01074       G4cout<<"G4QDiffRatio::PrFrag:oPDG=="<<oPDG<<",hPDG="<<hPDG<<",M="<<reM<<G4endl;//|
01075 #endif
01076       G4int rlPDG=hPDG-nL*1000000-nSP*1000999-nSM*999001;// Subtract Lamb/Sig from Nucl.|
01077       G4int    sPDG=3122;                      // Prototype for the Hyperon PDG (Lambda)|
01078       G4double MLa=mLamb;                      // Prototype for one Hyperon decay       |
01079 #ifdef fdebug
01080       G4cout<<"G4QDiffRatio::PrFrag:*G4*nS+="<<nSP<<",nS-="<<nSM<<",nL="<<nL<<G4endl;// |
01081 #endif
01082       if(nSP||nSM)                         // Sigma+/- improvement                      |
01083       {
01084         if(nL)                             // By mistake Lambda improvement is found    |
01085         {
01086           G4cout<<"***G4QDR::PFr:HypN="<<hPDG<<": bothSigm&Lamb -> ImproveIt"<<G4endl;//|
01087           //throw G4QException("*G4QDiffractionRatio::Fragment:BothLambda&SigmaInHN");//|
01088           // @@ Correction, which does not conserv the charge !! (-> add decay in 3)    |
01089           if(nSP) nL+=nSP;                 // Convert Sigma+ to Lambda                  |
01090           else    nL+=nSM;                 // Convert Sigma- to Lambda                  |
01091         }
01092         if(nSP)                            // Sibma+ should be decayed                  |
01093         {
01094           nL=nSP;                          // #of decaying hyperons                     |
01095           sPDG=3222;                       // PDG code of decaying hyperons             |
01096           MLa=mSigP;                       // Mass of decaying hyperons                 |
01097         }
01098         else                               // Sibma+ should be decayed                  |
01099         {
01100           nL=nSM;                          // #of decaying hyperons                     |
01101           sPDG=3112;                       // PDG code of decaying hyperons             |
01102           MLa=mSigM;                       // Mass of decaying hyperons                 |
01103         }
01104       }
01105 #ifdef fdebug
01106       G4cout<<"G4QDiffRat::ProjFrag:*G4*mS="<<MLa<<",sPDG="<<sPDG<<",nL="<<nL<<G4endl;//|
01107 #endif
01108       if(nL>1) MLa*=nL;                    // Total mass of the decaying hyperons       |
01109       G4double rlM=G4QNucleus(rlPDG).GetMZNS();// Mass of the NonstrangeNucleus         |
01110       if(!nSP&&!nSM&&nL==1&&reM>rlM+mSigZ&&G4UniformRand()>.5) // Conv Lambda->Sigma0   |
01111       {
01112         sPDG=3212;                         // PDG code of a decaying hyperon            |
01113         MLa=mSigZ;                         // Mass of the decaying hyperon              |
01114       }
01115       G4int rnPDG = hPDG-nL*999999;        // Convert Lambdas to neutrons (for convInN) |
01116       G4QNucleus rnN(rnPDG);               // New nonstrange nucleus                    |
01117       G4double rnM=rnN.GetMZNS();          // Mass of the new nonstrange nucleus        |
01118       // @@ In future take into account Iso-Hypernucleus (Add PI+,R & Pi-,R decays)     |
01119       if(rlPDG==90000000)                  // Multy Hyperon (HyperNuc of only hyperons) |
01120       {
01121         if(nL>1) r4M=r4M/nL;               // split the 4-mom for the MultyLambda       |
01122         for(G4int il=0; il<nL; il++)       // loop over Lambdas                         |
01123         {
01124           if(anti) sPDG=-sPDG;             // For anti-nucleus case                     |
01125           G4QHadron* theLam = new G4QHadron(sPDG,r4M); // Make NewHadr for the Hyperon  |
01126           ResHV->push_back(theLam);        // Fill in the Lambda                        |
01127 #ifdef fdebug
01128           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01129           G4cout<<"G4QDR::ProjFrag: *additional Lambda*="<<r4M<<sPDG<<G4endl; //        |
01130 #endif
01131         }
01132       }
01133       else if(reM>rlM+MLa-eps)              // Lambda (or Sigma) can be split           |
01134       {
01135         G4LorentzVector n4M(0.,0.,0.,rlM);  // 4-mom of the residual nucleus            |
01136         G4LorentzVector h4M(0.,0.,0.,MLa);  // 4-mom of the Hyperon                     |
01137         G4double sum=rlM+MLa;               // Safety sum                               |
01138         if(std::fabs(reM-sum)<eps)          // At rest in CMS                           |
01139         {
01140           n4M=r4M*(rlM/sum);                // Split tot 4-mom for resNuc               |
01141           h4M=r4M*(MLa/sum);                // Split tot 4-mom for Hyperon              |
01142         }
01143         else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay         |
01144         {
01145           G4cerr<<"***G4QDF::PF:HypN,M="<<reM<<"<A+n*L="<<sum<<",d="<<sum-reM<<G4endl;//|
01146           // throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclusDecay");//|
01147           G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0100",
01148                       FatalException, "Error in hypernuclus decay");
01149         }
01150 #ifdef fdebug
01151         G4cout<<"*G4QDR::PF:HypN="<<r4M<<"->A="<<rlPDG<<n4M<<",n*L="<<nL<<h4M<<G4endl;//|
01152 #endif
01153         loh->Set4Momentum(n4M);            // ! Update the Hadron !                     |
01154         if(anti && rlPDG==90000001) rlPDG=-2112; // Convert to anti-neutron             |
01155         if(anti && rlPDG==90001000) rlPDG=-2212; // Convert to anti-proton              |
01156         loh->SetQPDG(G4QPDGCode(rlPDG));   // ConvertedHypernucleus to nonstrange(@anti)|
01157         if(rlPDG==90000002)                // Additional action with loH changed to 2n  |
01158         {
01159           G4LorentzVector newLV=n4M/2.;    // Split 4-momentum                          |
01160           loh->Set4Momentum(newLV);        // Reupdate the hadron                       |
01161           if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG            |
01162           else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG                      |
01163           G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron             |
01164           ResHV->push_back(secHadr);       // Fill in the additional neutron            |
01165 #ifdef fdebug
01166           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01167           G4cout<<"G4QDR::ProgFrag: *additional Neutron*="<<r4M<<sPDG<<G4endl; //       |
01168 #endif
01169         }
01170         else if(rlPDG==90002000)           // Additional action with loH change to 2p   |
01171         {
01172           G4LorentzVector newLV=n4M/2.;    // Split 4-momentum                          |
01173           loh->Set4Momentum(newLV);        // Reupdate the hadron                       |
01174           if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG            |
01175           else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG                      |
01176           G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton              |
01177           ResHV->push_back(secHadr);       // Fill in the additional neutron            |
01178 #ifdef fdebug
01179           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01180           G4cout<<"G4QDR::ProjFrag: *additional Proton*="<<r4M<<sPDG<<G4endl; //        |
01181 #endif
01182         }
01183         // @@(?) Add multybaryon decays if necessary (Now it anyhow is made later)      |
01184 #ifdef fdebug
01185         G4cout<<"*G4QDiffractionRatio::PrFrag:resNucPDG="<<loh->GetPDGCode()<<G4endl;// |
01186 #endif
01187         if(nL>1) h4M=h4M/nL;               // split the lambda's 4-mom if necessary     |
01188         for(G4int il=0; il<nL; il++)       // A loop over excessive hyperons            |
01189         {
01190           if(anti) sPDG=-sPDG;             // For anti-nucleus case                     |
01191           G4QHadron* theLamb = new G4QHadron(sPDG,h4M); // Make NewHadr for the Hyperon |
01192           ResHV->push_back(theLamb);       // Fill in the additional neutron            |
01193 #ifdef fdebug
01194           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01195           G4cout<<"G4QDR::ProjFrag: *additional Hyperon*="<<r4M<<sPDG<<G4endl; //       |
01196 #endif
01197         }
01198       }
01199       else if(reM>rnM+mPi0-eps&&!nSP&&!nSM)// Lambda->N only if Sigmas are absent       |
01200       {
01201         G4int nPi=static_cast<G4int>((reM-rnM)/mPi0); // Calc. pion multiplicity        |
01202         if (nPi>nL) nPi=nL;                // Cut the pion multiplicity                 |
01203         G4double npiM=nPi*mPi0;            // Total pion mass                           |
01204         G4LorentzVector n4M(0.,0.,0.,rnM); // Residual nucleus 4-momentum               |
01205         G4LorentzVector h4M(0.,0.,0.,npiM);// 4-momentum of pions                       |
01206         G4double sum=rnM+npiM;             // Safety sum                                |
01207         if(std::fabs(reM-sum)<eps)         // At rest                                   |
01208         {
01209           n4M=r4M*(rnM/sum);               // The residual nucleus part                 |
01210           h4M=r4M*(npiM/sum);              // The pion part                             |
01211         }
01212         else if(reM<sum || !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay         |
01213         {
01214           G4cerr<<"*G4QDR::PF:HypN,M="<<reM<<"<A+n*Pi0="<<sum<<",d="<<sum-reM<<G4endl;//|
01215           // throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclDecay"); // |
01216           G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0101",
01217                        FatalException, "Error in HypernuclDecay"); 
01218         }
01219         loh->Set4Momentum(n4M);            // ! Update the Hadron !                     |
01220         if(anti && rnPDG==90000001) rnPDG=-2112; // Convert to anti-neutron             |
01221         if(anti && rnPDG==90001000) rnPDG=-2212; // Convert to anti-proton              |
01222         loh->SetQPDG(G4QPDGCode(rnPDG));   // convert hyperNuc to nonstrangeNuc(@@anti) |
01223 #ifdef fdebug
01224         G4cout<<"*G4QDR::PF:R="<<r4M<<"->A="<<rnPDG<<n4M<<",n*Pi0="<<nPi<<h4M<<G4endl;//|
01225 #endif
01226         if(nPi>1) h4M=h4M/nPi;             // Split the 4-mom if necessary              |
01227         for(G4int ihn=0; ihn<nPi; ihn++)   // A loop over additional pions              |
01228         {
01229           G4QHadron* thePion = new G4QHadron(111,h4M); // Make a New Hadr for the pi0   |
01230           ResHV->push_back(thePion);       // Fill in the Pion                          |
01231 #ifdef fdebug
01232           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01233           G4cout<<"G4QDR::ProjFrag: *additional Pion*="<<r4M<<sPDG<<G4endl; //          |
01234 #endif
01235         }
01236         if(rnPDG==90000002)                // Additional action with loH change to 2n   |
01237         {
01238           G4LorentzVector newLV=n4M/2.;    // Split 4-momentum                          |
01239           loh->Set4Momentum(newLV);        // Reupdate the hadron                       |
01240           if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG            |
01241           else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG                      |
01242           G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron             |
01243           ResHV->push_back(secHadr);       // Fill in the additional neutron            |
01244 #ifdef fdebug
01245           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01246           G4cout<<"G4QDR::ProjFrag: *additional Neutron*="<<r4M<<sPDG<<G4endl; //       |
01247 #endif
01248         }
01249         else if(rnPDG==90002000)           // Additional action with loH change to 2p   |
01250         {
01251           G4LorentzVector newLV=n4M/2.;    // Split 4-momentum                          |
01252           loh->Set4Momentum(newLV);        // Reupdate the hadron                       |
01253           if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG            |
01254           else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG                      |
01255           G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton              |
01256           ResHV->push_back(secHadr);       // Fill in the additional neutron            |
01257 #ifdef fdebug
01258           sum4M+=r4M;                      // Sum 4-momenta for the EnMom check         |
01259           G4cout<<"G4QDR::ProjFrag: *additional Proton*="<<r4M<<sPDG<<G4endl; //        |
01260 #endif
01261         }
01262         // @@ Add multybaryon decays if necessary                                       |
01263       }
01264       else // If this Excepton shows up (lowProbable appearance) => include gamma decay |
01265       {
01266         G4double d=rlM+MLa-reM;            // Hyperon Excessive energy                  |
01267         G4cerr<<"G4QDR::PF:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d="<<d<<G4endl;
01268         d=rnM+mPi0-reM;                    // Pion Excessive energy                     |
01269         G4cerr<<"G4QDR::PF:"<<oPDG<<","<<hPDG<<",M="<<reM<<"<"<<rnM+mPi0<<",d="<<d<<G4endl;
01270         // throw G4QException("G4QDiffractionRatio::ProjFragment: Hypernuclear conver");// |
01271         G4Exception("G4QDiffractionRatio::ProjFragment()", "HAD_CHPS_0102",
01272                     FatalException, "Excessive hypernuclear energy");
01273       }
01274     }                                      // => End of G4 Hypernuclear decay           |
01275     ResHV->push_back(loh);                 // Fill in the result                        |
01276 #ifdef debug
01277     sum4M+=loh->Get4Momentum();            // Sum 4-momenta for the EnMom check         |
01278     G4cout<<"G4QDR::PrFra:#"<<iq<<","<<loh->Get4Momentum()<<loh->GetPDGCode()<<G4endl;//|
01279 #endif
01280   }                                        //                                           |
01281   leadhs->clear();//                                                                    |
01282   delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<--*
01283 #ifdef debug
01284   G4cout<<"G4QDiffractionRatio::ProjFragment: *End* Sum="<<sum4M<<" =?= d4M="<<d4M<<G4endl;
01285 #endif
01286   return ResHV; // Result
01287 } // End of ProjFragment

G4QHadronVector * G4QDiffractionRatio::TargFragment ( G4int  pPDG,
G4LorentzVector  p4M,
G4int  tgZ,
G4int  tgN 
)

Definition at line 303 of file G4QDiffractionRatio.cc.

References G4QEnvironment::AddQuasmon(), FatalException, G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4RandomDirection(), and G4UniformRand.

Referenced by G4QDiffraction::PostStepDoIt().

00305 {
00306   static const G4double pFm= 0.; // Fermi momentum in MeV (delta function)
00307   //static const G4double pFm= 250.; // Fermi momentum in MeV (delta function)
00308   static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function)
00309   static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction)
00310   //static const G4double mPi= G4QPDGCode(211).GetMass();  // pi+- mass (MeV)
00311   static const G4double mNeut= G4QPDGCode(2112).GetMass();
00312   static const G4double mNeut2=mNeut*mNeut;
00313   static const G4double dmNeut=mNeut+mNeut;
00314   static const G4double mProt= G4QPDGCode(2212).GetMass();
00315   static const G4double mProt2=mProt*mProt;
00316   static const G4double dmProt=mProt+mProt;
00317   static const G4double maxDM=mProt*12.;
00318   //static const G4double mLamb= G4QPDGCode(3122).GetMass();
00319   //static const G4double mSigZ= G4QPDGCode(3212).GetMass();
00320   //static const G4double mSigM= G4QPDGCode(3112).GetMass();
00321   //static const G4double mSigP= G4QPDGCode(3222).GetMass();
00322   //static const G4double eps=.003;
00323   static const G4double third=1./3.;
00324   //
00325   G4LorentzVector pr4M=p4M/megaelectronvolt;   // Convert 4-momenta in MeV (keep p4M)
00326   // prepare the DONOTHING answer
00327   G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !!
00328   G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile
00329   ResHV->push_back(hadron);                // It must be cleaned up for real scattering sec
00330   // @@ diffraction is simulated as noncoherent (coherent is small)
00331   G4int tgA=tgZ+tgN;                       // A of the target
00332   G4int tPDG=90000000+tgZ*1000+tgN;        // PDG code of the targetNucleus/recoilNucleus
00333   G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus
00334   G4int rPDG=2112;                         // prototype of PDG code of the recoiled nucleon
00335   if(tgA*G4UniformRand()>tgN)              // Substitute by a proton
00336   {
00337     rPDG=2212;                             // PDG code of the recoiled QF nucleon
00338     tPDG-=1000;                            // PDG code of the recoiled nucleus
00339   }
00340   else tPDG-=1;                            // PDG code of the recoiled nucleus
00341   G4double tM=G4QPDGCode(tPDG).GetMass();  // Mass of the recoiled nucleus
00342   G4double tE=std::sqrt(tM*tM+pFm2);       // Free energy of the recoil nucleus
00343   G4ThreeVector tP=pFm*G4RandomDirection();// 3-mom of the recoiled nucleus
00344   G4LorentzVector t4M(tP,tE);              // 4M of the recoil nucleus
00345   G4LorentzVector tg4M(0.,0.,0.,tgM);      // Full target 4-momentum
00346   G4LorentzVector N4M=tg4M-t4M;            // 4-mom of Quasi-free target nucleon
00347   G4LorentzVector tot4M=N4M+p4M;           // total momentum of quasi-free diffraction
00348   G4double mT=mNeut;                       // Prototype of mass of QF nucleon
00349   G4double mT2=mNeut2;                     // Squared mass of a free nucleon to be excited
00350   G4double dmT=dmNeut;                     // Doubled mass              
00351   //G4int Z=0;                               // Prototype of the isotope Z
00352   //G4int N=1;                               // Prototype of the Isotope N
00353   if(rPDG==2212)                           // Correct it, if this is a proton
00354   {
00355     mT=mProt;                              // Prototype of mass of QF nucleon to be excited
00356     mT2=mProt2;                            // Squared mass of the free nucleon
00357     dmT=dmProt;                            // Doubled mass              
00358     //Z=1;                                   // Z of the isotope
00359     //N=0;                                   // N of the Isotope
00360   }
00361   G4double mP2=pr4M.m2();                  // Squared mass of the projectile
00362   if(mP2<0.) mP2=0.;                       // Can be a problem for photon (m_min = 2*m_pi0)
00363   G4double s_value=tot4M.m2();             // @@ Check <0 ...
00364   G4double E=(s_value-mT2-mP2)/dmT;        // Effective interactionEnergy (virtNucl target)
00365   G4double E2=E*E;
00366   if(E<0. || E2<mP2)                       // Impossible to fragment: return projectile
00367   {
00368 #ifdef pdebug
00369     G4cerr<<"-Warning-G4DifR::TFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
00370 #endif
00371     return ResHV;                          // *** Do Nothing Action ***
00372   }
00373   G4double mP=std::sqrt(mP2);              // Calculate mass of the projectile (to be exc.)
00374   if(mP<.1) mP=mPi0;                       // For photons minDiffraction is gam+P->P+Pi0
00375   //G4double dmP=mP+mP;                      // Doubled mass of the projectile
00376   G4double mMin=mP+mPi0;                   // Minimum diffractive mass
00377   G4double tA=tgA;                         // Real A of the target
00378   G4double sA=5./std::pow(tA,third);       // Mass-screaning
00379   //mMin+=mPi0+G4UniformRand()*(mP*sA+mPi0); // *Experimental*
00380   mMin+=G4UniformRand()*(mP*sA+mPi0);      // *Experimental*
00381   G4double ss=std::sqrt(s_value);          // CM compound mass (sqrt(s))
00382   G4double mMax=ss-mP;                     // Maximum diffraction mass of the projectile
00383   if(mMax>maxDM) mMax=maxDM;               // Restriction to avoid too big masses
00384   if(mMin>=mMax)
00385   {
00386 #ifdef pdebug
00387     G4cerr<<"-Warning-G4DifR::TFra:ZeroDiffractionMRange, mi="<<mMin<<",ma="<<mMax<<G4endl;
00388 #endif
00389     return ResHV;                          // Do Nothing Action
00390   }
00391   G4double R = G4UniformRand();
00392   G4double mDif=std::exp(R*std::log(mMax)+(1.-R)*std::log(mMin)); // Low Mass Approximation
00393   G4double mDif2=mDif*mDif;
00394   G4double ds=s_value-mP2-mDif2;               
00395   //G4double e=ds/dmP;
00396   //G4double P=std::sqrt(e*e-mDif2);      // Momentum in pseudo laboratory system
00397 #ifdef debug
00398   G4cout<<"G4QDiffR::TargFrag:Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
00399 #endif
00400   // @@ Temporary NN t-dependence for all hadrons
00401   if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl;
00402   G4double maxt=(ds*ds-4*mP2*mDif2)/s_value;  // maximum possible -t
00403   G4double tsl=140000.;                 // slope in MeV^2 
00404   G4double t=-std::log(G4UniformRand())*tsl;
00405 #ifdef pdebug
00406   G4cout<<"G4QDifR::TFra:ph="<<pPDG<<",P="<<P<<",t="<<t<<"<"<<maxt<<G4endl;
00407 #endif
00408 #ifdef nandebug
00409   if(mint>-.0000001);                   // To make the Warning for NAN
00410   else  G4cout<<"******G4QDiffractionRatio::TargFragment: -t="<<mint<<G4endl;
00411 #endif
00412   G4double rt=t/maxt;
00413   G4double cost=1.-rt-rt;               // cos(theta) in CMS
00414 #ifdef ppdebug
00415   G4cout<<"G4QDiffraRatio::TargFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
00416 #endif
00417   if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00418   {
00419     if     (cost>1.)  cost=1.;
00420     else if(cost<-1.) cost=-1.;
00421     else
00422     {
00423       G4cerr<<"G4QDiffRat::TargFragm: *NAN* cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl;
00424       return ResHV;                     // Do Nothing Action
00425     }
00426   }
00427   G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mP);      // 4mom of the leading nucleon
00428   G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif);    // 4mom of the diffract. Quasmon
00429   G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00430   if(!G4QHadron(tot4M).RelDecayIn2(r4M, d4M, dir4M, cost, cost))
00431   {
00432     G4cerr<<"G4QDifR::TFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl;
00433     //G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp");
00434     return ResHV; // Do Nothing Action
00435   }
00436 #ifdef debug
00437   G4cout<<"G4QDifRat::TargFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl;
00438 #endif
00439   // Now everything is ready for fragmentation and DoNothing projHadron must be wiped out
00440   delete hadron;     // Delete the fake (doNothing) projectile hadron
00441   ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile
00442   hadron = new G4QHadron(pPDG,r4M);     // Hadron for the recoil nucleon
00443   ResHV->push_back(hadron);             // Fill the recoil nucleon
00444 #ifdef debug
00445   G4cout<<"G4QDiffractionRatio::TargFragm: *Filled* LeadingNuc="<<r4M<<pPDG<<G4endl;
00446 #endif
00447   G4QHadronVector* leadhs = 0;   // Prototype of Quasmon Output G4QHadronVector  ---->---*
00448   G4QContent dQC=G4QPDGCode(rPDG).GetQuarkContent(); // QuarkContent of quasiFreeNucleon | 
00449   G4Quasmon* quasm = new G4Quasmon(dQC,d4M); // Quasmon=DiffractionExcitationQuasmon-*   |
00450 #ifdef debug
00451   G4cout<<"G4QDiffRatio::TgFrag:tPDG="<<tPDG<<",rPDG="<<rPDG<<",d4M="<<d4M<<G4endl;//|   |
00452 #endif
00453   G4QEnvironment* pan= new G4QEnvironment(G4QNucleus(tPDG));// --> DELETED --->---*  |   |
00454   pan->AddQuasmon(quasm);                    // Add diffractiveQuasmon to Environ.|  |   |
00455 #ifdef debug
00456   G4cout<<"G4QDiffractionRatio::TargFragment: EnvPDG="<<tPDG<<G4endl; //          |  |   |
00457 #endif
00458   try                                                           //                |  |   |
00459   {                                                             //                |  |   |
00460     leadhs = pan->Fragment();// DESTROYED in the end of the LOOP work space       |  | <-|
00461   }                                                             //                |  |   |
00462   catch (G4QException& error)//                                                   |  |   |
00463   {                                                             //                |  |   |
00464     //#ifdef pdebug
00465     G4cerr<<"***G4QDiffractionRatio::TargFrag: G4QException is catched"<<G4endl;//|  |   |
00466     //#endif
00467     //  G4Exception("G4QDiffractionRatio::TargFragm:","27",FatalException,"*Nucl");// |  |   |
00468     G4Exception("G4QDiffractionRatio::TargFragment()","HAD_CHPS_0027",
00469                 FatalException, "Nucl");
00470   }                                                             //                |  |   |
00471   delete pan;                              // Delete the Nuclear Environment <-<--*--*   |
00472   G4int qNH=leadhs->size();                // A#of collected hadrons from diff.frag.     |
00473   if(qNH) for(G4int iq=0; iq<qNH; iq++)    // Loop over hadrons to fill the result       |
00474   {                                        //                                            |
00475     G4QHadron* loh=(*leadhs)[iq];          // Pointer to the output hadron               |
00476     ResHV->push_back(loh);                 // Fill in the result                         |
00477   }                                        //                                            |
00478   leadhs->clear();//                                                                     |
00479   delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<---*
00480   return ResHV; // Result
00481 } // End of TargFragment


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:05 2013 for Geant4 by  doxygen 1.4.7