Changeset 1055 for trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QCollision.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QCollision.cc (modified) (60 diffs)

trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QCollision.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4QCollision.cc,v 1.28 2008/10/02 21:10:07 dennis Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4QCollision.cc,v 1.32 2009/04/17 15:22:31 mkossov Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //      ---------------- G4QCollision class -----------------
@@ -34 +34 @@
 // ********** This CLASS is temporary moved from the photolepton_hadron directory *********
 // ****************************************************************************************
-
+// Short description: This is a universal class for the incoherent (inelastic)
+// nuclear interactions in the CHIPS model.
+// ---------------------------------------------------------------------------
 //#define debug
 //#define pdebug
+//#define ldebug
 //#define ppdebug
 //#define qedebug
@@ -250 +253 @@
     G4cout<<"G4QCollision::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
 #endif
-    if(isoSize)                             // The Element has non-trivial abumdance set
+    if(isoSize)                             // The Element has non-trivial abundance set
     {
       indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
@@ -262 +265 @@
           G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
           if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QCollision::GetMeanFreePath"
-                                                                                                                                                                                                                                                                        <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
+                                 <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
           G4double abund=abuVector[j];
-                                                                  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
+          std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
 #ifdef debug
           G4cout<<"G4QCollision::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
 #endif
           newAbund->push_back(pr);
-                                                  }
+        }
 #ifdef debug
         G4cout<<"G4QCollision::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
@@ -384 +387 @@
   static const G4double mPi0s= mPi0*mPi0;
   static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
+  //static const G4double mDeut2=mDeut*mDeut;                      // SquaredMassOfDeuteron
   static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
   static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
@@ -432 +436 @@
   static G4bool CWinit = true;              // CHIPS Warld needs to be initted
   if(CWinit)
-                {
+  {
     CWinit=false;
     G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
@@ -457 +461 @@
 #ifdef debug
   G4double mp=proj4M.m();
-  G4cout<<"G4QCollis::PostStepDoIt:called, P="<<Momentum<<"="<<momentum<<",m="<<mp<<G4endl;
+  G4cout<<"-->G4QCollis::PostStDoIt:called,P="<<Momentum<<"="<<momentum<<",m="<<mp<<G4endl;
 #endif
   if (!IsApplicable(*particle))  // Check applicability
@@ -507 +511 @@
 #ifdef debug
   G4int prPDG=particle->GetPDGEncoding();
-                G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
+  G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
 #endif
   if(!projPDG)
@@ -516 +520 @@
   G4int EPIM=ElProbInMat.size();
 #ifdef debug
-                G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
+  G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
 #endif
   G4int i=0;
@@ -523 +527 @@
     G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
     for(i=0; i<nE; ++i)
-                  {
-#ifdef debug
-                                  G4cout<<"G4QCollision::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
+    {
+#ifdef debug
+      G4cout<<"G4QCollision::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
 #endif
       if (rnd<ElProbInMat[i]) break;
@@ -534 +538 @@
   Z=static_cast<G4int>(pElement->GetZ());
 #ifdef debug
-                                G4cout<<"G4QCollision::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
+    G4cout<<"G4QCollision::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
 #endif
   if(Z<=0)
@@ -545 +549 @@
   G4int nofIsot=SPI->size();               // #of isotopes in the element i
 #ifdef debug
-                G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
+  G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
 #endif
   G4int j=0;
@@ -554 +558 @@
     {
 #ifdef debug
-                                  G4cout<<"G4QCollision::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
+      G4cout<<"G4QCollision::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
 #endif
       if(rndI < (*SPI)[j]) break;
@@ -562 +566 @@
   G4int N =(*IsN)[j]; ;                    // Randomized number of neutrons
 #ifdef debug
-                G4cout<<"G4QCollision::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<G4endl;
+  G4cout<<"G4QCollision::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl;
 #endif
   if(N<0)
@@ -576 +580 @@
   if(dsr<dd)dsr=dd;
   if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // ManualP
-                else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
+  else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
   else if(projPDG==-211)  G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);//Pi- ClustPars
 #ifdef debug
@@ -607 +611 @@
   // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
   //
-  if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15) {
-
+  if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15)
+  {
 #ifdef debug
     G4cout<<"G4QCollision::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
@@ -646 +650 @@
     G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
 #ifdef debug
-                                G4cout<<"G4QCol::PStDoIt: kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
+    G4cout<<"G4QCol::PStDoIt: kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
 #endif
     if( kinEnergy < photonEnergy || photonEnergy < 0.)
@@ -708 +712 @@
     G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
 #ifdef pdebug
-                  G4cout<<"G4QC::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
+    G4cout<<"G4QC::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
 #endif
     if(cost>1.) cost=1.;                 // To avoid the accuracy of calculation problem
@@ -717 +721 @@
     G4double absEn = std::pow(am,third)*GeV;  // @@(b) Mean Energy Absorbed by a Nucleus
     //if(am>1 && absEn < photonEnergy)     // --> the absorption of energy can happen
-                                if(absEn < photonEnergy)     // --> the absorption of energy can happen
+    if(absEn < photonEnergy)     // --> the absorption of energy can happen
     {
       G4double abtEn = absEn+hdM;        // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
@@ -727 +731 @@
       absMom         = std::sqrt(abMo2); // Absorbed Momentum
       if(absMom < phMo)                  // --> the absorption of momentum can happen
-                                  {
+      {
         G4double dEn = photonEnergy - absEn; // Leading energy
         G4double dMo = phMo - absMom;    // Leading momentum
         G4double sF  = dEn*dEn - dMo*dMo;// s of leading particle
 #ifdef ppdebug
-                                    G4cout<<"-PhotoAbsorption-G4QCol::PStDoIt:sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
+        G4cout<<"-PhotoAbsorption-G4QCol::PStDoIt:sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
 #endif
         if(sF > stmPi)                   // --> Leading fragmentation is possible
-                                                                {
+        {
           photonEnergy = absEn;          // New value of the photon energy
           photonQ2=abMo2-absEn*absEn;    // New value of the photon Q2
@@ -757 +761 @@
     aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
 #ifdef pdebug
-                  G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
+    G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
           <<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
 #endif
@@ -765 +769 @@
       G4double ptm=photon3M.mag();                    // 3M of the virtual photon
 #ifdef ppdebug
-                    G4cout<<"-Absorption-G4QCollision::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
+      G4cout<<"-Absorption-G4QCollision::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
             <<iniP*dir<<", eFin3M="<<finP*findir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
 #endif
@@ -772 +776 @@
       proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
 #ifdef ppdebug
-                    G4cout<<"-->G4QC::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
+      G4cout<<"-->G4QC::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
 #endif
       lead4M=proj4M;                        // Remember 4-mom for the total 4-momentum
       G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
       try                                                           //                    |
-             {                                                             //                    |
-               delete leadhs;                                              //                    |
+      {                                                             //                    |
+        delete leadhs;                                              //                    |
         G4QNucleus vac(90000000);                                   //                    |
         leadhs=pan->Fragment(vac,1);  // DELETED after it is copied to output vector      |
       }                                                             //                    |
       catch (G4QException& error)                                   //                    |
-             {                                                             //                    |
+      {                                                             //                    |
         G4cerr<<"***G4QCollision::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//|
         G4Exception("G4QCollision::PostStepDoIt:","72",FatalException,"QuasmonCrash");  //|
@@ -789 +793 @@
       delete pan;                              // Delete the Nuclear Environment <----<---+
 #ifdef ppdebug
-                                                G4cout<<"G4QCol::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
+      G4cout<<"G4QCol::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
             <<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
 #endif
@@ -796 +800 @@
     proj4M=G4LorentzVector(photon3M,photonEnergy);
 #ifdef debug
-                  G4cout<<"G4QCollision::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
+    G4cout<<"G4QCollision::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
           <<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
 #endif
@@ -803 +807 @@
   // neutrinoNuclear interactions (nu_e, nu_mu only)
   //
-  } else if (aProjPDG == 12 || aProjPDG == 14) {
-
+  }
+  else if (aProjPDG == 12 || aProjPDG == 14)
+  {
     G4double kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
     G4double dKinE=kinEnergy+kinEnergy;  // doubled energy for s calculation
 #ifdef debug
-                  G4cout<<"G4QCollision::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
+    G4cout<<"G4QCollision::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
 #endif
     G4ParticleMomentum dir = projHadron->GetMomentumDirection(); // unit vector
     G4double ml  = mu;
-    G4double ml2 =      mu2;
-    //G4double mlN =    muN;
-    G4double mlN2=      muN2;
-    G4double fmlN=      fmuN;
-    G4double mlsN=      musN;
-    //G4double mlP =    muP;
-    G4double mlP2=      muP2;
-    G4double fmlP=      fmuP;
-    G4double mlsP=      musP;
-    G4double mldM=      mudM;
-    //G4double mlD =    muD;
-    G4double mlD2=      muD2;
+    G4double ml2 = mu2;
+    //G4double mlN = muN;
+    G4double mlN2= muN2;
+    G4double fmlN= fmuN;
+    G4double mlsN= musN;
+    //G4double mlP = muP;
+    G4double mlP2= muP2;
+    G4double fmlP= fmuP;
+    G4double mlsP= musP;
+    G4double mldM= mudM;
+    //G4double mlD = muD;
+    G4double mlD2= muD2;
     if(aProjPDG==12)
     {
       ml  = me;
-      ml2 =     me2;
-      //mlN =   meN;
-      mlN2=     meN2;
-      fmlN=     fmeN;
-      mlsN=     mesN;
-      //mlP =   meP;
-      mlP2=     meP2;
-      fmlP=     fmeP;
-      mlsP=     mesP;
-      mldM=     medM;
-      //mlD =   meD;
-      mlD2=     meD2;
+      ml2 = me2;
+      //mlN = meN;
+      mlN2= meN2;
+      fmlN= fmeN;
+      mlsN= mesN;
+      //mlP = meP;
+      mlP2= meP2;
+      fmlP= fmeP;
+      mlsP= mesP;
+      mldM= medM;
+      //mlD = meD;
+      mlD2= meD2;
     }
     G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
@@ -881 +886 @@
     G4bool secnu=false;
     if(xSec*G4UniformRand()>xSec1)               // recover neutrino/antineutrino
-                                {
+    {
       if(scatPDG>0) scatPDG++;
       else          scatPDG--;
@@ -900 +905 @@
       totCS1 = qelCS1;
       totCS2 = qelCS2;
-                                }
+    }
     // make different definitions for neutrino and antineutrino
     G4double mIN=mProt;                          // Just a prototype (for anu, Z=1, N=0)
@@ -956 +961 @@
     {
       if(Z>1||N>0)                               // Calculate the splitted mass
-                                                {
+      {
         targPDG-=1000;                           // Anti-Neutrino -> subtract proton
         G4QPDGCode targQPDG(targPDG);
@@ -973 +978 @@
     G4double s=mIN*(mIN+dKinE);                  // s=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
 #ifdef debug
-                  G4cout<<"G4QCollision::PostStDoIt: s="<<s<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
+    G4cout<<"G4QCollision::PostStDoIt: s="<<s<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
 #endif
     if(s<=OT)                                    // *** Do nothing solution ***
@@ -986 +991 @@
     }
 #ifdef debug
-                G4cout<<"G4QCollision::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
+    G4cout<<"G4QCollision::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
 #endif
     aParticleChange.ProposeEnergy(0.);
@@ -998 +1003 @@
       else      Q2=CSmanager->GetQEL_ExchangeQ2();
 #ifdef debug
-                  G4cout<<"G4QCollision::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s<<",Q2="<<Q2<<G4endl;
+      G4cout<<"G4QCollision::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s<<",Q2="<<Q2<<G4endl;
 #endif
       //G4double ds=s+s;                          // doubled s
@@ -1032 +1037 @@
       else      Q2=CSmanager2->GetNQE_ExchangeQ2();
 #ifdef debug
-                  G4cout<<"G4QColl::PStDoIt: MultiPeriferal s="<<s<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
+      G4cout<<"G4QColl::PStDoIt: MultiPeriferal s="<<s<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
 #endif
       if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
@@ -1055 +1060 @@
       {
 #ifdef debug
-                    G4cout<<"*G4QCollision::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
+        G4cout<<"*G4QCollision::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
               <<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
 #endif
@@ -1077 +1082 @@
       proj4M-=scat4M;                           // 4mom of the W/Z (effective pion/gamma)
 #ifdef debug
-                  G4cout<<"G4QCollision::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
+      G4cout<<"G4QCollision::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
 #endif
       // Check that the en/mom transfer is possible, if not -> elastic
@@ -1162 +1167 @@
       }
     }
-
   //
-  // quasi-elastic for p+A(Z,N)
+  // quasi-elastic (& pickup process) for p+A(Z,N)
   //
-  } else if (aProjPDG == 2212 && Z > 0 && N > 0) {
-    //else if(2>3)
-
+  }
+  else if (aProjPDG == 2212 && Z > 0 && N > 0)
+  //else if(2>3) 
+  {
     G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
     std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, aProjPDG, Z, N);
     G4double qepart=fief.first*fief.second;
-#ifdef qedebug
-    G4cout<<"G4QCol::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<",)="<<qepart<<G4endl;
-#endif
-    if(G4UniformRand()<qepart) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim
-    {
-      // First decay a nucleus in a nucleon and a residual (A-1) nucleus
-      G4double dmom=91.; // Fermi momentum (proto default for a deuteron)
-      if(Z>1||N>1) dmom=286.2*std::pow(-std::log(G4UniformRand()),third);// p_max=250 MeV/c
-
+    // Keep QE part for the quasi-free nucleons
+    const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
+    G4double clProb[lCl]={.65,.85,.95};// N/P,D,t/He3,He4, integratedProb for .65,.2,.1,.05
+    if(qepart>0.45) qepart=.45; // Quasielastic part is too large - shrink
+    else qepart/=clProb[0];     // Add corresponding number of 2N, 3N, & 4N clusters
+    G4double pickup=1.-qepart;  // Estimate the rest of the cross-section
+    G4double thresh=100.;
+    if(momentum > thresh) pickup*=50./momentum/std::pow(Z+N,third); // 50. is a par(!)
+    // pickup = 0.;               // To exclude the pickup process
+    if (N) pickup+=qepart;
+    else   pickup =qepart;
+    G4double rnd=G4UniformRand();
+#ifdef ldebug
+    G4cout<<"-->G4QCol::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<")]="<<qepart
+          <<", pickup="<<pickup<<G4endl;
+#endif
+    if(rnd<pickup) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim,CoulBar,PauliBl
+    {
+     G4double dmom=91.;  // Fermi momentum (proto default for a deuteron)
+     G4double fmm=287.;  // @@ Can be A-dependent
+     if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
+     // Values to be redefined for quasi-elastic
+     G4LorentzVector r4M(0.,0.,0.,0.);      // Prototype of 4mom of the residual nucleus
+     G4LorentzVector n4M(0.,0.,0.,0.);      // Prototype of 4mom of the quasi-cluster
+     G4int nPDG=90000001;                   // Prototype for quasi-cluster mass calculation
+     G4int restPDG=targPDG;                 // Prototype should be reduced by quasi-cluster
+     G4int rA=Z+N-1;                        // Prototype for the residualNucl definition
+     G4int rZ=Z;                            // residZ: OK for the quasi-free neutron
+     G4int nA=1;                            // Prototype for the quasi-cluster definition
+     G4int nZ=0;                            // nA=1,nZ=0: OK for the quasi-free neutron
+     G4double qM=mNeut;                     // Free mass of the quasi-free cluster
+     G4int A = Z + N;                       // Baryon number of the nucleus
+     if(rnd<qepart)
+     {
+      // ===> First decay a nucleus in a nucleon and a residual (A-1) nucleus
       // Calculate cluster probabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
-      const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
-      G4double clProb[lCl]={0.6,0.7,0.8}; // N/P,D,t/He3,Alpha, integrated prob for .6,.1,.1,.2
       G4double base=1.;  // Base for randomization (can be reduced by totZ & totN)
-      G4int max=lCl;   // Number of boundaries (can be reduced by totZ & totN)
-
+      G4int max=lCl;     // Number of boundaries (can be reduced by totZ & totN)
       // Take into account that at least one nucleon must be left !
-      // Change max-- to --max - DHW 05/08
-      G4int A = Z + N;       // Baryon number of the nucleus
       if (Z<2 || N<2 || A<5) base = clProb[--max]; // Alpha cluster is impossible
       if ( (Z > 1 && N < 2) || (Z < 2 && N > 1) ) 
         base=(clProb[max]+clProb[max-1])/2; // t or He3 is impossible
-
-      if ( (Z < 2 && N < 2) || A < 4) base=clProb[--max]; // Both He3 and t clusters are impossible
-
+      if ( (Z < 2 && N < 2) || A < 4) base=clProb[--max];// He3 & t clusters are impossible
       if(A<3)           base=clProb[--max]; // Deuteron cluster is impossible
       G4int cln=0;                          // Cluster#0 (Default for the selected nucleon)
@@ -1207 +1231 @@
       G4ParticleDefinition* theDefinition;  // Prototype for qfNucleon
       G4bool cp1 = cln+2==A;                // A=ClusterBN+1 condition
-      // Values to be defined in the following IF/ELSE
-      G4LorentzVector r4M(0.,0.,0.,0.);     // Prototype of 4mom of the residual nucleus
-      G4LorentzVector n4M(0.,0.,0.,0.);     // Prototype of 4mom of the quasi-cluster
-      G4int nPDG=90000001;                  // Prototype for quasi-cluster mass calculation
-      G4int restPDG=targPDG;                // Prototype should be reduced by quasi-cluster
-      G4int rA=Z+N-1;                       // Prototype for the residualNucl definition
-      G4int rZ=Z;                           // residZ: OK for the quasi-free neutron
-      G4int nA=1;                           // Prototype for the quasi-cluster definition
-      G4int nZ=0;                           // nA=1,nZ=0: OK for the quasi-free neutron
-      G4double qM=mNeut;                    // Free mass of the quasi-free cluster
       if(!cln || cp1)                       // Split in nucleon + (A-1) with Fermi momentum
       {
@@ -1239 +1253 @@
         n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
 #ifdef qedebug
-                      G4cout<<"G4QCollis::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
+        G4cout<<"G4QCollis::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
 #endif
         if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
@@ -1247 +1261 @@
         }
 #ifdef qedebug
-                      G4cout<<"G4QCol::PStDoIt:QE-N,RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
+        G4cout<<"G4QCol::PStDoIt:QE-N,RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
 #endif
         if(cp1 && cln)                           // Quasi-cluster case: swap the output
@@ -1266 +1280 @@
         }
       }
-      else // Split a cluster (w or w/o "Fermi motion" and "Fermi decay")
+      else // Split the cluster (with or w/o "Fermi motion" and "Fermi decay")
       {
         if(cln==1)
@@ -1275 +1289 @@
           nZ=1;
           restPDG-=1001;
+          // Recalculate dmom
+          G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
+                      fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+          dmom=sumv.mag();
         }
         else if(cln==2)
         {
           nA=3;
-                                                                                if(G4UniformRand()*(A-2)>(N-1)) // He3
+          if(G4UniformRand()*(A-2)>(N-1)) // He3
           {
             nPDG=90002001;
@@ -1293 +1311 @@
             restPDG-=1002;
           }
+          // Recalculate dmom
+          G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+          dmom=sumv.mag();
         }
         else
@@ -1301 +1324 @@
           nZ=2;
           restPDG-=2002;
+          G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+          dmom=sumv.mag();
         }
         rA=A-nA;
@@ -1317 +1345 @@
         n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
 #ifdef qedebug
-                      G4cout<<"G4QCollis::PStDoIt:QEC,p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
+        G4cout<<"G4QCollis::PStDoIt:QEC,p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
 #endif
         if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
@@ -1326 +1354 @@
         // --- End of the moving cluster implementation ---
 #ifdef qedebug
-                      G4cout<<"G4QCol::PStDoIt:QEC,RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
+        G4cout<<"G4QCol::PStDoIt:QEC,RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
 #endif
       }
       G4LorentzVector s4M=n4M+proj4M;             // Tot 4-momentum for scattering
-      G4double prjM2 = proj4M.m2();
+      G4double prjM2 = proj4M.m2();               // Squared mass of the projectile
       G4double prjM = std::sqrt(prjM2);           // @@ Get from pPDG (?)
       G4double minM = prjM+qM;                    // Min mass sum for the final products
@@ -1337 +1365 @@
       {
 #ifdef qedebug
-                      G4cout<<"G4QCol::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl;
+        G4cout<<"G4QCol::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl;
 #endif
         // Estimate and randomize charge-exchange with quasi-free cluster
         G4bool chex=false;                        // Flag of the charge exchange scattering
         G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
+        // This is reserved for the charge-exchange scattering (to help neutron spectras)
         //if(cln&&!cp1 &&(projPDG==2212&&rA>rZ || projPDG==2112&&rZ>1))// @@ Use proj chex
-               if(2>3)
+        if(2>3)
         {
 #ifdef qedebug
-                        G4cout<<"G4QCol::PStDoIt:-Enter,P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
+          G4cout<<"G4QCol::PStDoIt:-Enter,P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
 #endif
           G4double tprM=mProt;
@@ -1365 +1394 @@
           G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
 #ifdef qedebug
-                        G4cout<<"G4QCol::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
+          G4cout<<"G4QCol::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
 #endif
           if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl))     // minM is redefined
@@ -1381 +1410 @@
                                                                          projPDG, proj4M);
 #ifdef qedebug
-                      G4cout<<"G4QCollis::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
+        G4cout<<"G4QCollis::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
               <<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
 #endif
@@ -1431 +1460 @@
       else G4cout<<"G4QCol::PSD: OUT, M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
 #endif
-    }
+     } // end of proton quasi-elastic (including QE on NF)
+     else // @@ make cost-condition @@ Pickup process (pickup 1 or 2 n and make d or t)
+     {
+       restPDG--;                // Res. nucleus PDG is one neutron less than targ. PDG
+       G4double mPUF=mDeut;      // Prototype of the mass of the created pickup fragment
+       G4ParticleDefinition* theDefinition = G4Deuteron::Deuteron();
+       //theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
+       G4bool din=false;         // Flag of picking up 2 neutrons to create t (tritium)
+       //G4bool hin=false;         // Flag of picking up 1 n + 1 p to create He3 (NOT IMPL)
+       G4bool tin=false;         // Flag of picking up 2 n + 1 p to create alpha (He4)
+       G4double frn=G4UniformRand();
+       if(N>2 && frn > 0.95)      // .95 is a parameter to pickup 2n (cteate t or +1p=alph)
+       {
+         theDefinition = G4Triton::Triton();
+         restPDG--;              // Res. nucleus PDG is two neutrons less than target PDG
+         rA--;                   // Reduce the baryon number of the residual nucleus
+         mPUF=mTrit;             // The mass of the created pickup fragment (triton)
+         // Recalculate dmom
+         G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
+                       fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+         dmom=sumv.mag();
+         din=true;
+         if(Z>1 && frn > 0.99)  // 0.99 is a parameter to pickup 2n + 1p and create alpha
+         {
+           theDefinition = G4Alpha::Alpha();
+           restPDG-=1000;        // Res. nucleus PDG is (2 n + 1 p) less than target PDG
+           rA--;                 // Reduce the baryon number of the residual nucleus
+           rZ--;                 // Reduce the charge of the residual nucleus
+           mPUF=mAlph;           // The mass of the created pickup fragment (triton)
+           // Recalculate dmom
+           sumv += fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+           dmom=sumv.mag();
+           tin=true;
+         }
+       }
+       G4double mPUF2=mPUF*mPUF;
+       G4LorentzVector t4M(0.,0.,0.,tM);         // 4m of the target nucleus to be decayed
+       G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
+       G4double rM2=rM*rM;                       // Squared mass of the residual nucleus
+       G4double nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom);// M2 of boundQF-nucleon(2N)
+       if(nM2 < 0) nM2=0.;
+       G4double nM=std::sqrt(nM2);               // M of bQF-nucleon
+       G4double den2=(dmom*dmom+nM2);            // squared energy of bQF-nucleon
+       G4double den=std::sqrt(den2);             // energy of bQF-nucleon
+#ifdef qedebug
+       G4cout<<"G4QCollis::PStDoIt:PiU,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
+#endif
+       G4double qp=momentum*dmom;
+       G4double srm=proj4M.e()*den;
+       G4double cost=(nM2+mProt2+srm+srm-mPUF2)/(qp+qp);
+       G4int ict=0;
+       while(std::fabs(cost)>1. && ict<3)
+       {
+         dmom=91.;  // Fermi momentum (proto default for a deuteron)
+         if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
+         if(din)          // Recalculate dmom for the final triton
+         {
+          G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
+                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+          if(tin) sumv+=fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
+          dmom=sumv.mag();
+         }
+         nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom); // M2 of bQF-nucleon
+         if(nM2 < 0) nM2=0.;
+         nM=std::sqrt(nM2);                      // M of bQF-nucleon
+         den2=(dmom*dmom+nM2);                   // squared energy of bQF-nucleon
+         den=std::sqrt(den2);                    // energy of bQF-nucleon
+         qp=momentum*dmom;
+         srm=proj4M.e()*den;
+         cost=(nM2+mProt2+srm+srm-mPUF2)/(qp+qp);
+         ict++;
+#ifdef ldebug
+         if(ict>2)G4cout<<"G4QCollis::PStDoIt:i="<<ict<<",d="<<dmom<<",ct="<<cost<<G4endl;
+#endif
+       }
+       if(std::fabs(cost)<=1.)
+       {// ---- @@ From here can be a MF for QF nucleon extraction (if used by others)
+        G4ThreeVector vdir = proj4M.vect();  // 3-Vector in the projectile direction
+        G4ThreeVector vx(0.,0.,1.);          // ProtoOrt in theDirection of the projectile
+        G4ThreeVector vy(0.,1.,0.);          // First ProtoOrt orthogonal to the direction
+        G4ThreeVector vz(1.,0.,0.);          // SecondProtoOrt orthoganal to the direction
+        if(vdir.mag2() > 0.)                 // the projectile isn't at rest
+        {
+          vx = vdir.unit();                  // Ort in the direction of the projectile
+          G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
+          vy = vv.unit();                    // First ort orthogonal to the proj direction
+          vz = vx.cross(vy);                 // Second ort orthoganal to the projDirection
+        }
+        // ---- @@ End of possible MF (Similar is in G4QEnvironment)
+        G4double tdom=dmom*std::sqrt(1-cost*cost);// Transversal component of the Fermi-Mom
+        G4double phi=twopi*G4UniformRand();  // Phi of the Fermi-Mom
+        G4ThreeVector vedm=vx*dmom*cost+vy*tdom*std::sin(phi)+vz*tdom*std::cos(phi);// bQFN
+        G4LorentzVector bqf(vedm,den);      // 4-mom of the bQF nucleon
+        r4M=t4M-bqf;                         // 4mom of the residual nucleus
+        if(std::fabs(r4M.m()-rM)>.001) G4cout<<"G4QCol::PSD:rM="<<rM<<"#"<<r4M.m()<<G4endl;
+        G4LorentzVector f4M=proj4M+bqf;      // Prototype of 4-mom of Deuteron
+        if(std::fabs(f4M.m()-mPUF)>.001)
+                              G4cout<<"G4QCol::PSD:mDeut="<<mPUF<<" # "<<f4M.m()<<G4endl;
+        aParticleChange.ProposeEnergy(0.);   // @@ ??
+        aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed
+        aParticleChange.SetNumberOfSecondaries(2); 
+        // Fill pickup hadron
+        G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,f4M);
+        G4Track* scatQFN = new G4Track(theQFN, localtime, position ); //    pickup
+        scatQFN->SetWeight(weight);                                   //    weighted
+        scatQFN->SetTouchableHandle(trTouchable);                     //    nuclear
+        aParticleChange.AddSecondary(scatQFN);                        //    cluster
+        // ----------------------------------------------------
+        // Fill residual nucleus
+        if     (restPDG==90000001) theDefinition = G4Neutron::Neutron();
+        else if(restPDG==90001000) theDefinition = G4Proton::Proton();
+        else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
+        G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
+        G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
+        scatReN->SetWeight(weight);                                   //    weighted
+        scatReN->SetTouchableHandle(trTouchable);                     //    residual
+        aParticleChange.AddSecondary(scatReN);                        //    nucleus
+        return G4VDiscreteProcess::PostStepDoIt(track, step);
+         }
+     }
+    } // end of decoupled processes for the proton projectile
   }  // end lepto-nuclear, neutrino-nuclear, proton quasi-elastic
-
   EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM);    // Total 4-mom of the reaction
   if(absMom) EnMomConservation+=lead4M;         // Add E/M of leading System
@@ -1444 +1592 @@
   //G4double eWei=1.;
   // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End  
-#ifdef debug
-  G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
-#endif
-  G4QHadron* pH = new G4QHadron(projPDG,proj4M);                // ---> DELETED -->--  -+
+#ifdef ldebug
+  G4cout<<"^^G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
+#endif
+  G4QHadron* pH = new G4QHadron(projPDG,proj4M);                // ---> DELETED -->----+*
   //if(momentum<1000.) // Condition for using G4QEnvironment (not G4QuasmonString)      |
   { //                                                                                  |
@@ -1465 +1613 @@
     catch (G4QException& error)//                                                   |   .
     {                                                             //                |   .
-             //#ifdef pdebug
+      //#ifdef pdebug
       G4cerr<<"***G4QCollision::PostStepDoIt: G4QE Exception is catched"<<G4endl;// |   .
-             //#endif
+      //#endif
       G4Exception("G4QCollision::PostStepDoIt:","27",FatalException,"CHIPSCrash");//|   .
     }                                                             //                |   .
@@ -1483 +1631 @@
   //  try                                                           //                 |
   //  {                                                             //                 |
-  //              delete output;                                  //                   |
+  //    delete output;                                  //                   |
   //    output = pan->Fragment();// DESTROYED in the end of the LOOP work space        |
   //    // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin                 |
   //    //tNH=pan->GetNOfHadrons();     // For the test purposes of the String         |
   //    //if(tNH==2)                    // At least 2 hadrons are in the Constr.Output |
-  //    //{//                                                                          |
+  // //{//                                                                          |
   //    //  elF=true;                   // Just put a flag for the ellastic Scattering |
   //    //  delete output;              // Delete a prototype of dummy G4QHadronVector |
@@ -1503 +1651 @@
   //    //#ifdef pdebug
   //    G4cerr<<"***G4QCollision::PostStepDoIt: GEN Exception is catched"<<G4endl; //  |
-         //    //#endif
+  //    //#endif
   //    G4Exception("G4QCollision::PostStDoIt:","27",FatalException,"QString Excep");//|
   //  }                                                             //                 |
@@ -1524 +1672 @@
     delete leadhs;
   }
-
-
   // ------------- From here the secondaries are filled -------------------------
   G4int tNH = output->size();       // A#of hadrons in the output
@@ -1565 +1711 @@
     {
 #ifdef debug
-             G4cout<<"G4QCollision::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
+      G4cout<<"G4QCollision::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
 #endif
       delete hadr;
@@ -1578 +1724 @@
     {
       if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong();   // K_L
-                                                else                   theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
+      else                   theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
     }
     else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
@@ -1585 +1731 @@
     else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
     else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
-           else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
+    else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
     {
       G4int aZ = hadr->GetCharge();
       G4int aA = hadr->GetBaryonNumber();
 #ifdef pdebug
-                                                G4cout<<"G4QCollision::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
+      G4cout<<"G4QCollision::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
 #endif
       theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
@@ -1598 +1744 @@
     {
 #ifdef pdebug
-                                                G4cout<<"G4QCollision::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
+      G4cout<<"G4QCollision::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
 #endif
       theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
 #ifdef pdebug
-                                                G4cout<<"G4QCollision::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
+      G4cout<<"G4QCollision::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
 #endif
     }
@@ -1643 +1789 @@
   delete output; // instances of the G4QHadrons from the output are already deleted above +
 #ifdef debug
-                G4cout<<"G4QCollision::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
+  G4cout<<"G4QCollision::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
 #endif
   if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
@@ -1649 +1795 @@
   //return &aParticleChange;                               // This is not enough (ClearILL)
 #ifdef pdebug
-                  G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()
+    G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()
           <<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
 #endif
-#ifdef debug
-                G4cout<<"G4QCollision::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St="
+#ifdef ldebug
+  G4cout<<"G4QCollision::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St="
         <<aParticleChange.GetTrackStatus()<<G4endl;
 #endif