G4QLowEnergy.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 // GEANT4 tag $Name: geant4-09-04 $
00028 //
00029 //      ---------------- G4QLowEnergy class -----------------
00030 //                 by Mikhail Kossov, Aug 2007.
00031 // G4QLowEnergy class of the CHIPS Simulation Branch in GEANT4
00032 // ---------------------------------------------------------------
00033 // Short description: This is a fast low energy algorithm for the
00034 // inelastic interactions of nucleons and nuclei (ions) with nuclei.
00035 // This is a fase-space algorithm, but not quark level. Provides
00036 // nuclear fragments upto alpha only. Never was tumed (but can be).
00037 // ---------------------------------------------------------------
00038 //#define debug
00039 //#define pdebug
00040 //#define edebug
00041 //#define tdebug
00042 //#define nandebug
00043 //#define ppdebug
00044 
00045 #include "G4QLowEnergy.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4HadronicDeprecate.hh"
00048 
00049 
00050 // Initialization of static vectors
00051 //G4int    G4QLowEnergy::nPartCWorld=152;  // #of particles initialized in CHIPS
00052 //G4int    G4QLowEnergy::nPartCWorld=122;  // #of particles initialized in CHIPS
00053 G4int    G4QLowEnergy::nPartCWorld=85;     // #of particles initialized in CHIPS Reduced
00054 std::vector<G4int> G4QLowEnergy::ElementZ; // Z of element(i) in theLastCalc
00055 std::vector<G4double> G4QLowEnergy::ElProbInMat; // SumProbOfElem in Material
00056 std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;// N of isotope(j), E(i)
00057 std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;//SumProbIsotE(i)
00058 
00059 // Constructor
00060 G4QLowEnergy::G4QLowEnergy(const G4String& processName):
00061   G4VDiscreteProcess(processName, fHadronic), evaporate(true)
00062 {
00063   G4HadronicDeprecate("G4QLowEnergy");
00064 
00065 #ifdef debug
00066   G4cout<<"G4QLowEnergy::Constructor is called processName="<<processName<<G4endl;
00067 #endif
00068   if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created "<<G4endl;
00069   G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
00070 }
00071 
00072 // Destructor
00073 G4QLowEnergy::~G4QLowEnergy() {}
00074 
00075 
00076 G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;}
00077 
00078 G4int G4QLowEnergy::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
00079 
00080 // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
00081 // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
00082 // ********** All CHIPS cross sections are calculated in the surface units ************
00083 G4double G4QLowEnergy::GetMeanFreePath(const G4Track&Track, G4double, G4ForceCondition*F)
00084 {
00085   *F = NotForced;
00086   const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
00087   G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00088   if( !IsApplicable(*incidentParticleDefinition))
00089     G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl;
00090   // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
00091   G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
00092 #ifdef debug
00093   G4double KinEn = incidentParticle->GetKineticEnergy();
00094   G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
00095 #endif
00096   const G4Material* material = Track.GetMaterial();        // Get the current material
00097   const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00098   const G4ElementVector* theElementVector = material->GetElementVector();
00099   G4int nE=material->GetNumberOfElements();
00100 #ifdef debug
00101   G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl;
00102 #endif
00103   G4int pPDG=0;
00104   G4int Z=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge());
00105   G4int A=incidentParticleDefinition->GetBaryonNumber();
00106   if      ( incidentParticleDefinition == G4Proton::Proton()     ) pPDG = 2212;
00107   else if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020;
00108   else if ( incidentParticleDefinition == G4Alpha::Alpha()       ) pPDG = 1000020040;
00109   else if ( incidentParticleDefinition == G4Triton::Triton()     ) pPDG = 1000010030;
00110   else if ( incidentParticleDefinition == G4He3::He3()           ) pPDG = 1000020030;
00111   else if ( incidentParticleDefinition == G4GenericIon::GenericIon() || (Z > 0 && A > 0))
00112   {
00113     pPDG=incidentParticleDefinition->GetPDGEncoding();
00114 #ifdef debug
00115     G4int prPDG=1000000000+10000*A+10*Z;
00116     G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl;
00117 #endif
00118   }
00119   else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl;
00120   G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
00121   if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); // just to test
00122   Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A
00123   G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
00124   G4double sigma=0.;                        // Sums over elements for the material
00125   G4int IPIE=IsoProbInEl.size();            // How many old elements?
00126   if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
00127   {
00128     std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
00129     SPI->clear();
00130     delete SPI;
00131     std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
00132     IsN->clear();
00133     delete IsN;
00134   }
00135   ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
00136   ElementZ.clear();                         // Clear the body vector for Z of Elements
00137   IsoProbInEl.clear();                      // Clear the body vector for SPI
00138   ElIsoN.clear();                           // Clear the body vector for N of Isotopes
00139   for(G4int i=0; i<nE; ++i)
00140   {
00141     G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
00142     Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
00143     ElementZ.push_back(Z);                  // Remember Z of the Element
00144     G4int isoSize=0;                        // The default for the isoVectorLength is 0
00145     G4int indEl=0;                          // Index of non-natural element or 0(default)
00146     G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
00147     if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
00148 #ifdef debug
00149     G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
00150 #endif
00151     if(isoSize)                             // The Element has non-trivial abundance set
00152     {
00153       indEl=pElement->GetIndex()+1;         // Index of the non-trivial element is an order
00154 #ifdef debug
00155       G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00156 #endif
00157       if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
00158       {
00159         std::vector<std::pair<G4int,G4double>*>* newAbund =
00160                                                new std::vector<std::pair<G4int,G4double>*>;
00161         G4double* abuVector=pElement->GetRelativeAbundanceVector();
00162         for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
00163         {
00164           G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
00165           if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
00166                                          <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00167           G4double abund=abuVector[j];
00168           std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00169 #ifdef debug
00170           G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl;
00171 #endif
00172           newAbund->push_back(pr);
00173         }
00174 #ifdef debug
00175         G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
00176 #endif
00177         indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
00178         for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
00179         delete newAbund; // Was "new" in the beginning of the name space
00180       }
00181     }
00182     std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
00183     std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
00184     IsoProbInEl.push_back(SPI);
00185     std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
00186     ElIsoN.push_back(IsN);
00187     G4int nIs=cs->size();                   // A#Of Isotopes in the Element
00188 #ifdef debug
00189     G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
00190 #endif
00191     G4double susi=0.;                       // sum of CS over isotopes
00192     if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
00193     {
00194       std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
00195       G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
00196       IsN->push_back(N);                    // Remember Min N for the Element
00197 #ifdef debug
00198       G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl;
00199 #endif
00200       G4bool ccsf=true;                    // Extract inelastic Ion-Ion cross-section
00201 #ifdef debug
00202       G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl;
00203 #endif
00204       G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope
00205 #ifdef debug
00206       G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
00207             <<Momentum<<", XSec="<<CSI/millibarn<<G4endl;
00208 #endif
00209       curIs->second = CSI;
00210       susi+=CSI;                            // Make a sum per isotopes
00211       SPI->push_back(susi);                 // Remember summed cross-section
00212     } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
00213     sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
00214 #ifdef debug
00215     G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
00216           <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
00217 #endif
00218     ElProbInMat.push_back(sigma);
00219   } // End of LOOP over Elements
00220   // Check that cross section is not zero and return the mean free path
00221 #ifdef debug
00222   G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
00223 #endif
00224   if(sigma > 0.000000001) return 1./sigma;                 // Mean path [distance] 
00225   return DBL_MAX;
00226 }
00227 
00228 G4bool G4QLowEnergy::IsApplicable(const G4ParticleDefinition& particle) 
00229 {
00230   G4int Z=static_cast<G4int>(particle.GetPDGCharge());
00231   G4int A=particle.GetBaryonNumber();
00232   if      (particle == *(     G4Proton::Proton()     )) return true;
00233   else if (particle == *(    G4Neutron::Neutron()    )) return true;
00234   else if (particle == *(   G4Deuteron::Deuteron()   )) return true;
00235   else if (particle == *(      G4Alpha::Alpha()      )) return true;
00236   else if (particle == *(     G4Triton::Triton()     )) return true;
00237   else if (particle == *(        G4He3::He3()        )) return true;
00238   else if (particle == *( G4GenericIon::GenericIon() )) return true;
00239   else if (Z > 0 && A > 0)                              return true;
00240 #ifdef debug
00241   G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<", A="
00242         <<A<<", Z="<<Z<<G4endl;
00243 #endif
00244   return false;
00245 }
00246 
00247 G4VParticleChange* G4QLowEnergy::PostStepDoIt(const G4Track& track, const G4Step& step)
00248 {
00249   static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV
00250   static const G4double mPro2= mProt*mProt;                    // CHIPS sq proton Mass
00251   static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV
00252   static const G4double mNeu2= mNeut*mNeut;                    // CHIPS sq neutron Mass
00253   static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV
00254   static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV;
00255   static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV;
00256   static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV;
00257   static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV;
00258   static const G4double mFm= 250*MeV;
00259   static const G4double third= 1./3.;
00260   static const G4ThreeVector zeroMom(0.,0.,0.);
00261   static G4ParticleDefinition* aGamma    = G4Gamma::Gamma();
00262   static G4ParticleDefinition* aPiZero   = G4PionZero::PionZero();
00263   static G4ParticleDefinition* aPiPlus   = G4PionPlus::PionPlus();
00264   static G4ParticleDefinition* aPiMinus  = G4PionMinus::PionMinus();
00265   static G4ParticleDefinition* aProton   = G4Proton::Proton();
00266   static G4ParticleDefinition* aNeutron  = G4Neutron::Neutron();
00267   static G4ParticleDefinition* aLambda   = G4Lambda::Lambda();
00268   static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
00269   static G4ParticleDefinition* aTriton   = G4Triton::Triton();
00270   static G4ParticleDefinition* aHe3      = G4He3::He3();
00271   static G4ParticleDefinition* anAlpha   = G4Alpha::Alpha();
00272   static const G4int nCh=26;                         // #of combinations
00273   static G4QNucleus Nuc;                             // A fake nucleus to call Evaporation
00274   //
00275   //-------------------------------------------------------------------------------------
00276   static G4bool CWinit = true;                       // CHIPS Warld needs to be initted
00277   if(CWinit)
00278   {
00279     CWinit=false;
00280     G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
00281   }
00282   //-------------------------------------------------------------------------------------
00283   const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00284   const G4ParticleDefinition* particle=projHadron->GetDefinition();
00285 #ifdef pdebug
00286   G4cout<<"G4QLowEnergy::PostStepDoIt: *** Called *** In4M="
00287         <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
00288         <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
00289 #endif
00290   //G4ForceCondition cond=NotForced;
00291   //GetMeanFreePath(track, -27., &cond);              // @@ ?? jus to update parameters?
00292 #ifdef debug
00293   G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
00294 #endif
00295   std::vector<G4Track*> result;
00296   G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
00297   G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
00298   G4double Momentum = proj4M.rho();                   // @@ Just for the test purposes
00299   if(std::fabs(Momentum-momentum)>.000001)
00300     G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
00301 #ifdef debug
00302   G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M,m="
00303         <<proj4M<<proj4M.m()<<G4endl;
00304 #endif
00305   if (!IsApplicable(*particle))  // Check applicability
00306   {
00307     G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl;
00308     return 0;
00309   }
00310   const G4Material* material = track.GetMaterial();      // Get the current material
00311   const G4ElementVector* theElementVector = material->GetElementVector();
00312   G4int nE=material->GetNumberOfElements();
00313 #ifdef debug
00314   G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00315 #endif
00316   G4int projPDG=0;                           // PDG Code prototype for the captured hadron
00317   // Not all these particles are implemented yet (see Is Applicable)
00318   G4int Z=static_cast<G4int>(particle->GetPDGCharge());
00319   G4int A=particle->GetBaryonNumber();
00320   if      (particle ==      G4Proton::Proton()     ) projPDG= 2212;
00321   else if (particle ==     G4Neutron::Neutron()    ) projPDG= 2112;
00322   else if (particle ==    G4Deuteron::Deuteron()   ) projPDG= 1000010020;
00323   else if (particle ==       G4Alpha::Alpha()      ) projPDG= 1000020040;
00324   else if (particle ==      G4Triton::Triton()     ) projPDG= 1000010030;
00325   else if (particle ==         G4He3::He3()        ) projPDG= 1000020030;
00326   else if (particle ==  G4GenericIon::GenericIon() || (Z > 0 && A > 0))
00327   {
00328     projPDG=particle->GetPDGEncoding();
00329 #ifdef debug
00330     G4int prPDG=1000000000+10000*Z+10*A; // just for the testing purposes
00331     G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl;
00332 #endif
00333   }
00334   else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl;
00335 #ifdef pdebug
00336   G4int prPDG=particle->GetPDGEncoding();
00337   G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
00338 #endif
00339   if(!projPDG)
00340   {
00341     G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
00342     return 0;
00343   }
00344   // Element treatment
00345   G4int EPIM=ElProbInMat.size();
00346 #ifdef debug
00347   G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00348 #endif
00349   G4int i=0;
00350   if(EPIM>1)
00351   {
00352     G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00353     for(i=0; i<nE; ++i)
00354     {
00355 #ifdef debug
00356       G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl;
00357 #endif
00358       if (rnd<ElProbInMat[i]) break;
00359     }
00360     if(i>=nE) i=nE-1;                        // Top limit for the Element
00361   }
00362   G4Element* pElement=(*theElementVector)[i];
00363   G4int tZ=static_cast<G4int>(pElement->GetZ());
00364 #ifdef debug
00365     G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl;
00366 #endif
00367   if(tZ<=0)
00368   {
00369     G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl;
00370     if(tZ<0) return 0;
00371   }
00372   std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
00373   std::vector<G4int>* IsN = ElIsoN[i];     // Vector of "#of neutrons" in the isotope El[i]
00374   G4int nofIsot=SPI->size();               // #of isotopes in the element i
00375 #ifdef debug
00376   G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
00377 #endif
00378   G4int j=0;
00379   if(nofIsot>1)
00380   {
00381     G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
00382     for(j=0; j<nofIsot; ++j)
00383     {
00384 #ifdef debug
00385       G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
00386 #endif
00387       if(rndI < (*SPI)[j]) break;
00388     }
00389     if(j>=nofIsot) j=nofIsot-1;            // Top limit for the isotope
00390   }
00391   G4int tN =(*IsN)[j]; ;                    // Randomized number of neutrons
00392 #ifdef debug
00393   G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl;
00394 #endif
00395   if(tN<0)
00396   {
00397     G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl;
00398     return 0;
00399   }
00400   nOfNeutrons=tN;                           // Remember it for the energy-momentum check
00401 #ifdef debug
00402   G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl;
00403 #endif
00404   if(tN<0)
00405   {
00406     G4cerr<<"***Warning***G4QLowEnergy::PostStepDoIt: Element with N="<<tN<< G4endl;
00407     return 0;
00408   }
00409   aParticleChange.Initialize(track);
00410 #ifdef debug
00411   G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl;
00412 #endif
00413   G4double weight        = track.GetWeight();
00414   G4double localtime     = track.GetGlobalTime();
00415   G4ThreeVector position = track.GetPosition();
00416 #ifdef debug
00417   G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl;
00418 #endif
00419   G4TouchableHandle trTouchable = track.GetTouchableHandle();
00420 #ifdef debug
00421   G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl;
00422 #endif
00423   G4int      targPDG=90000000+tZ*1000+tN;    // Target PDG code
00424   G4QPDGCode targQPDG(targPDG);              // To get CHIPS properties
00425   G4double tM=targQPDG.GetMass();            // CHIPS target nucleus mass in MeV
00426   G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness
00427   G4int pZ=static_cast<G4int>(particle->GetPDGCharge());  // Charge of the projectile
00428   G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile
00429   G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV
00430   G4double cosp=-14*Momentum*(tM-pM)/tM/pM;  // Asymmetry power for angular distribution
00431 #ifdef debug
00432   G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
00433         <<","<<tN<<"), cosp="<<cosp<<G4endl;
00434 #endif
00435   G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
00436   G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
00437   G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM); // Target's 4-mom
00438   G4LorentzVector tot4M =proj4M+targ4M;      // Total 4-mom of the reaction
00439 #ifdef pdebug
00440   G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
00441 #endif
00442   EnMomConservation=tot4M;                   // Total 4mom of reaction for E/M conservation
00443   // --- Start of Coulomb barrier check ---
00444   G4double dER = tot4M.e() - tM - pM;        // Energy of the reaction
00445   G4double pA = pZ+pN;                       // Atomic weight of the projectile (chged blw)
00446   G4double tA = tZ+tN;                       // Atomic weight of the target (changed below)
00447   G4double CBE = Nuc.CoulombBarGen(tZ, tA, pZ, pA); // Coulomb Barrier of the reaction
00448 #ifdef debug
00449   G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
00450         <<","<<tN<<"), dE="<<dER<<", CB="<<CBE<<G4endl;
00451 #endif
00452   if(dER<CBE) // The cross-section iz 0 -> Do Nothing
00453   {
00454 #ifdef debug
00455     G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Below Coulomb barrier* PDG="<<projPDG
00456           <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
00457 #endif
00458     //Do Nothing Action insead of the reaction
00459     aParticleChange.ProposeEnergy(kinEnergy);
00460     aParticleChange.ProposeLocalEnergyDeposit(0.);
00461     aParticleChange.ProposeMomentumDirection(dir) ;
00462     return G4VDiscreteProcess::PostStepDoIt(track,step);
00463   }
00464   // --- End of Coulomb barrier check ---
00465   // @@ Probably this is not necessary any more
00466 #ifdef debug
00467   G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl;
00468 #endif
00469   G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection
00470 #ifdef pdebug
00471   G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",tZ="<<tZ<<",N="<<tN<<",XS="
00472         <<xSec/millibarn<<G4endl;
00473 #endif
00474 #ifdef nandebug
00475   if(xSec>0. || xSec<0. || xSec==0);
00476   else  G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
00477 #endif
00478   // @@ check a possibility to separate p, n, or alpha (!)
00479   if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
00480   {
00481 #ifdef debug
00482     G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Zero cross-section* PDG="<<projPDG
00483           <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
00484 #endif
00485     //Do Nothing Action insead of the reaction
00486     aParticleChange.ProposeEnergy(kinEnergy);
00487     aParticleChange.ProposeLocalEnergyDeposit(0.);
00488     aParticleChange.ProposeMomentumDirection(dir) ;
00489     return G4VDiscreteProcess::PostStepDoIt(track,step);
00490   }
00491   // Kill interacting hadron
00492   aParticleChange.ProposeEnergy(0.);
00493   aParticleChange.ProposeTrackStatus(fStopAndKill);
00494   G4int tB=tZ+tN;
00495   G4int pB=pZ+pN;
00496 #ifdef pdebug
00497   G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<", tA="<<tB<<", pA="<<pB<<G4endl;
00498 #endif
00499   // algorithm implementation STARTS HERE --- All calculations are in IU --------
00500   tA=tB;
00501   pA=pB;
00502   G4double tR=1.1;                    // target nucleus R in fm
00503   if(tB > 1) tR*=std::pow(tA,third);  // in fm
00504   G4double pR=1.1;                    // projectile nucleus R in fm
00505   if(pB > 1) pR*=std::pow(pA,third);  // in fm
00506   G4double R=tR+pR;                   // total radius
00507   G4double R2=R*R;
00508   G4int tD=0;
00509   G4int pD=0;
00510   G4int nAt=0;
00511   G4int nAtM=27;
00512   G4int nSec = 0;
00513   G4double tcM=0.;
00514   G4double tnM=1.;
00515 #ifdef edebug
00516   G4int totChg=0;
00517   G4int totBaN=0;
00518   G4LorentzVector tch4M =tot4M;       // Total 4-mom of the reaction
00519 #endif
00520   while((!tD && !pD && ++nAt<nAtM) || tcM<tnM)
00521   {
00522 #ifdef edebug
00523     totChg=tZ+pZ;
00524     totBaN=tB+pB;
00525     tch4M =tot4M;
00526     G4cout<<">G4QLEn::PSDI:#"<<nAt<<",tC="<<totChg<<",tA="<<totBaN<<",t4M="<<tch4M<<G4endl;
00527 #endif
00528     G4LorentzVector tt4M=tot4M;
00529     G4int resultSize=result.size();
00530     if(resultSize)
00531     {
00532       for(i=0; i<resultSize; ++i) delete result[i];
00533       result.clear();
00534     }
00535     G4double D=std::sqrt(R2*G4UniformRand());
00536 #ifdef pdebug
00537     G4cout<<"G4QLowEn::PSDI: D="<<D<<", tR="<<tR<<", pR="<<pR<<G4endl;
00538 #endif
00539     if(D > std::fabs(tR-pR))          // leading parts can be separated
00540     {
00541       nSec = 0;
00542       G4double DtR=D-tR;
00543       G4double DpR=D-pR;
00544       G4double DtR2=DtR*DtR;
00545       G4double DpR2=DpR*DpR;
00546       //G4double DtR3=DtR2*DtR;
00547       G4double DpR3=DpR2*DpR;
00548       //G4double DtR4=DtR3*DtR;
00549       G4double DpR4=DpR3*DpR;
00550       G4double tR2=tR*tR;
00551       G4double pR2=pR*pR;
00552       G4double tR3=tR2*tR;
00553       //G4double pR3=pR2*pR;
00554       G4double tR4=tR3*tR;
00555       //G4double pR4=pR3*pR;
00556       G4double HD=16.*D;
00557       G4double tF=tA*(6*tR2*(pR2-DtR2)-4*D*(tR3-DpR3)+3*(tR4-DpR4))/HD/tR3;
00558       G4double pF=tF;
00559       tD=static_cast<G4int>(tF);
00560       pD=static_cast<G4int>(pF);
00561       //if(G4UniformRand() < tF-tD) ++tD; // Simple solution
00562       //if(G4UniformRand() < pF-pD) ++pD;
00563       // Enhance alphas solution
00564       if(std::fabs(tF-4.) < 1.) tD=4;       // @@ 1. is the enhancement parameter
00565       else if(G4UniformRand() < tF-tD) ++tD;
00566       if(std::fabs(pF-4.) < 1.) pD=4;
00567       else if(G4UniformRand() < pF-pD) ++pD;
00568       if(tD > tB) tD=tB;
00569       if(pD > pB) pD=tB;
00570       // @@ Quasi Free is not debugged @@ The following close it
00571       if(tD < 1) tD=1;
00572       if(pD < 1) pD=1;
00573 #ifdef pdebug
00574       G4cout<<"G4QLE::PSDI:#"<<nAt<<",pF="<<pF<<",tF="<<tF<<",pD="<<pD<<",tD="<<tD<<G4endl;
00575 #endif
00576       G4int pC=0;                     // charge of the projectile fraction
00577       G4int tC=0;                     // charge of the target fraction
00578       if((tD || pD) && tD<tB && pD<pB)// Periferal interaction
00579       {
00580         if(!tD || !pD)                // Quasi-Elastic: B+A->(B-1)+N+A || ->B+N+(A-1)
00581         {
00582           G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
00583           G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
00584           G4int pPDG=2112;            // Proto of the nucleon PDG (proton)
00585           G4double prM =mNeut;        // Proto of the nucleon mass
00586           G4double prM2=mNeu2;        // Proto of the nucleon sq mass
00587           if     (!tD)                // Quasi-elastic scattering of the proj QF nucleon
00588           {
00589             if(pD > 1) pD=1;
00590             if(!pN || (pZ && pA*G4UniformRand() < pZ) ) // proj QF proton
00591             {
00592               pPDG=2212;
00593               prM=mProt;
00594               prM2=mPro2;
00595               pC=1;
00596             }
00597             G4double Mom=0.;
00598             G4LorentzVector com4M = targ4M; // Proto of 4mom of compound
00599             G4double tgM=targ4M.e();
00600             G4double rNM=0.;
00601             G4LorentzVector rNuc4M(0.,0.,0.,0.);
00602             if(pD==pB)
00603             {
00604               Mom=proj4M.rho();
00605               com4M += proj4M;        // Total 4-mom for scattering
00606               rNM=prM;
00607             }
00608             else                      // It is necessary to split the nucleon (with fermiM)
00609             {
00610               G4ThreeVector fm=mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
00611               rNM=G4QPDGCode(2112).GetNuclMass(pZ-pC, pN, 0);
00612               G4double rNE=std::sqrt(fm*fm+rNM*rNM);
00613               rNuc4M=G4LorentzVector(fm,rNE);
00614               G4ThreeVector boostV=proj4M.vect()/proj4M.e();
00615               rNuc4M.boost(boostV);
00616               G4LorentzVector qfN4M=proj4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
00617               com4M += qfN4M;         // Calculate Total 4Mom for NA scattering
00618               G4double pNE = qfN4M.e(); // Energy of the QF nucleon in LS
00619               if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2); // Mom(s) fake value
00620               else break;             // Break the while loop
00621             }
00622             xSec=0.;
00623             if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
00624             else           xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
00625             if( xSec <= 0. ) break;   // Break the while loop
00626             G4double mint=0.;                 // Prototype of functional randomized -t
00627             if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
00628             else           mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
00629             G4double maxt=0.;                           // Prototype of maximum -t
00630             if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
00631             else           maxt=NCSmanager->GetHMaxT(); // maximum -t
00632             G4double cost=1.-mint/maxt;       // cos(theta) in CMS
00633             if(cost>1. || cost<-1.) break; // Break the while loop
00634             G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM); // 4mom of recoil target
00635             G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
00636             G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
00637             if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
00638             {
00639               G4cout<<"G4QLE::Pt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
00640               break;                  // Break the while loop
00641             }
00642             G4Track* projSpect = 0;
00643             G4Track* aNucTrack = 0;
00644             if(pB > pD)               // Fill the proj residual nucleus
00645             {
00646               G4int rZ=pZ-pC;
00647               G4int rA=pB-1;
00648               G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
00649               if(rA==1)
00650               {
00651                 if(!rZ) theDefinition = aNeutron;
00652                 else    theDefinition = aProton;
00653               }
00654               else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
00655               G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
00656               projSpect = new G4Track(resN, localtime, position );
00657               projSpect->SetWeight(weight);                    //    weighted
00658               projSpect->SetTouchableHandle(trTouchable);
00659 #ifdef pdebug
00660              G4cout<<"G4QLowEn::PSDI:-->ProjQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
00661 #endif
00662               ++nSec;
00663             }
00664             G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
00665             if(pPDG==2212) theDefinition = G4Proton::Proton();
00666             G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
00667             aNucTrack = new G4Track(scatN, localtime, position );
00668             aNucTrack->SetWeight(weight);                    //    weighted
00669             aNucTrack->SetTouchableHandle(trTouchable);
00670 #ifdef pdebug
00671              G4cout<<"G4QLowEn::PSDI:-->TgQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
00672 #endif
00673             ++nSec;
00674             G4Track* aFraTrack=0;
00675             if(tA==1)
00676             {
00677               if(!tZ) theDefinition = aNeutron;
00678               else    theDefinition = aProton;
00679             }
00680             else if(tA==8 && tC==4)             // Be8 decay
00681             {
00682               theDefinition = anAlpha;
00683               G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
00684               G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
00685               if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
00686               {
00687                 G4cout<<"*G4QLE::TS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
00688               }
00689               G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
00690               aFraTrack = new G4Track(pAl, localtime, position );
00691               aFraTrack->SetWeight(weight);                    //    weighted
00692               aFraTrack->SetTouchableHandle(trTouchable);
00693 #ifdef pdebug
00694               G4cout<<"G4QLowEn::PSDI:-->TgRQFA4M="<<f4M<<G4endl;
00695 #endif
00696               ++nSec;
00697               reco4M=s4M;
00698             }
00699             else if(tA==5 && (tC==2 || tC==3))   // He5/Li5 decay
00700             {
00701               theDefinition = aProton;
00702               G4double mNuc = mProt;
00703               if(tC==2)
00704               {
00705                 theDefinition = aNeutron;
00706                 mNuc = mNeut;
00707               }
00708               G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc);  // 4mom of the nucleon
00709               G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
00710               if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
00711               {
00712                 G4cout<<"*G4QLE::TS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
00713               }
00714               G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
00715               aFraTrack = new G4Track(pAl, localtime, position );
00716               aFraTrack->SetWeight(weight);                    //    weighted
00717               aFraTrack->SetTouchableHandle(trTouchable);
00718 #ifdef pdebug
00719               G4cout<<"G4QLowEn::PSDI:-->TgQFRN4M="<<f4M<<G4endl;
00720 #endif
00721               ++nSec;
00722               theDefinition = anAlpha;
00723               reco4M=s4M;
00724             }
00725             else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(tZ,tB,0,tZ);
00726             ++nSec;
00727 #ifdef pdebug
00728             G4cout<<"G4QLowEn::PSDI:-->PQF_nSec="<<nSec<<G4endl;
00729 #endif
00730             aParticleChange.SetNumberOfSecondaries(nSec); 
00731             if(projSpect) aParticleChange.AddSecondary(projSpect);
00732             if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
00733             if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
00734             G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
00735             G4Track* aResTrack = new G4Track(resA, localtime, position );
00736             aResTrack->SetWeight(weight);                    //    weighted
00737             aResTrack->SetTouchableHandle(trTouchable);
00738 #ifdef pdebug
00739             G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
00740 #endif
00741             aParticleChange.AddSecondary(aResTrack);
00742           }
00743           else                         // !pD : QF target Nucleon on the whole Projectile
00744           {
00745             if(tD > 1) tD=1;
00746             if(!tN || (tZ && tA*G4UniformRand() < tZ) ) // target QF proton
00747             {
00748               pPDG=2212;
00749               prM=mProt;
00750               prM2=mPro2;
00751               tC=1;
00752             }
00753             G4double Mom=0.;
00754             G4LorentzVector com4M=proj4M; // Proto of 4mom of compound
00755             prM=proj4M.m();
00756             G4double rNM=0.;
00757             G4LorentzVector rNuc4M(0.,0.,0.,0.);
00758             if(tD==tB)
00759             {
00760               Mom=prM*proj4M.rho()/proj4M.m();
00761               com4M += targ4M;        // Total 4-mom for scattering
00762               rNM=prM;
00763             }
00764             else                      // It is necessary to split the nucleon (with fermiM)
00765             {
00766               G4ThreeVector fm=250.*std::pow(G4UniformRand(),third)*G4RandomDirection();
00767               rNM=G4QPDGCode(2112).GetNuclMass(tZ-tC, tN, 0)/MeV;
00768               G4double rNE=std::sqrt(fm*fm+rNM*rNM);
00769               rNuc4M=G4LorentzVector(fm,rNE);
00770               G4LorentzVector qfN4M=targ4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
00771               com4M += qfN4M;                 // Calculate Total 4Mom for NA scattering
00772               G4ThreeVector boostV=proj4M.vect()/proj4M.e();
00773               qfN4M.boost(-boostV);
00774               G4double tNE = qfN4M.e();       // Energy of the QF nucleon in LS
00775               if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value
00776               else break;                     // Break the while loop
00777             }
00778             xSec=0.;
00779             if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
00780             else           xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
00781             if( xSec <= 0. ) break;           // Break the while loop
00782             G4double mint=0.;                 // Prototype of functional randomized -t
00783             if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
00784             else           mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
00785             G4double maxt=0.;                           // Prototype of maximum -t
00786             if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
00787             else           maxt=NCSmanager->GetHMaxT(); // maximum -t
00788             G4double cost=1.-mint/maxt;                 // cos(theta) in CMS
00789             if(cost>1. || cost<-1.) break;    // Break the while loop
00790             G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM); // 4mom of recoil target
00791             G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
00792             G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
00793             if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
00794             {
00795               G4cout<<"G4QLE::Tt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
00796               break;                          // Break the while loop
00797             }
00798             G4Track* targSpect = 0;
00799             G4Track* aNucTrack = 0;
00800             if(tB > tD)                       // Fill the residual nucleus
00801             {
00802               G4int rZ=tZ-tC;
00803               G4int rA=tB-1;
00804               G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
00805               if(rA==1)
00806               {
00807                 if(!rZ) theDefinition = aNeutron;
00808                 else    theDefinition = aProton;
00809               }
00810               else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
00811               G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
00812               targSpect = new G4Track(resN, localtime, position );
00813               targSpect->SetWeight(weight);                    //    weighted
00814               targSpect->SetTouchableHandle(trTouchable);
00815 #ifdef pdebug
00816              G4cout<<"G4QLowEn::PSDI:-->TargQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
00817 #endif
00818               ++nSec;
00819             }
00820             G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
00821             if(pPDG==2212) theDefinition = G4Proton::Proton();
00822             G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
00823             aNucTrack = new G4Track(scatN, localtime, position );
00824             aNucTrack->SetWeight(weight);                    //    weighted
00825             aNucTrack->SetTouchableHandle(trTouchable);
00826 #ifdef pdebug
00827              G4cout<<"G4QLowEn::PSDI:-->PrQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
00828 #endif
00829             ++nSec;
00830             G4Track* aFraTrack=0;
00831             if(pA==1)
00832             {
00833               if(!pZ) theDefinition = aNeutron;
00834               else    theDefinition = aProton;
00835             }
00836             else if(pA==8 && pC==4)             // Be8 decay
00837             {
00838               theDefinition = anAlpha;
00839               G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
00840               G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
00841               if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
00842               {
00843                 G4cout<<"*G4QLE::PS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
00844               }
00845               G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
00846               aFraTrack = new G4Track(pAl, localtime, position );
00847               aFraTrack->SetWeight(weight);                    //    weighted
00848               aFraTrack->SetTouchableHandle(trTouchable);
00849 #ifdef pdebug
00850               G4cout<<"G4QLowEn::PSDI:-->PrRQFA4M="<<f4M<<G4endl;
00851 #endif
00852               ++nSec;
00853               reco4M=s4M;
00854             }
00855             else if(pA==5 && (pC==2 || pC==3))   // He5/Li5 decay
00856             {
00857               theDefinition = aProton;
00858               G4double mNuc = mProt;
00859               if(pC==2)
00860               {
00861                 theDefinition = aNeutron;
00862                 mNuc = mNeut;
00863               }
00864               G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc);  // 4mom of the nucleon
00865               G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
00866               if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
00867               {
00868                 G4cout<<"*G4QLE::PS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
00869               }
00870               G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
00871               aFraTrack = new G4Track(pAl, localtime, position );
00872               aFraTrack->SetWeight(weight);                    //    weighted
00873               aFraTrack->SetTouchableHandle(trTouchable);
00874 #ifdef pdebug
00875               G4cout<<"G4QLowEn::PSDI:-->PrQFRN4M="<<f4M<<G4endl;
00876 #endif
00877               ++nSec;
00878               theDefinition = anAlpha;
00879               reco4M=s4M;
00880             }
00881             else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(pZ,pB,0,pZ);
00882             ++nSec;
00883 #ifdef pdebug
00884             G4cout<<"G4QLowEn::PSDI:-->TQF_nSec="<<nSec<<G4endl;
00885 #endif
00886             aParticleChange.SetNumberOfSecondaries(nSec); 
00887             if(targSpect) aParticleChange.AddSecondary(targSpect);
00888             if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
00889             if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
00890             G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
00891             G4Track* aResTrack = new G4Track(resA, localtime, position );
00892             aResTrack->SetWeight(weight);                    //    weighted
00893             aResTrack->SetTouchableHandle(trTouchable);
00894 #ifdef pdebug
00895             G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
00896 #endif
00897             aParticleChange.AddSecondary( aResTrack );
00898           }
00899 #ifdef debug
00900           G4cout<<"G4QLowEnergy::PostStepDoIt:***PostStepDoIt is done:Quasi-El***"<<G4endl;
00901 #endif
00902           return G4VDiscreteProcess::PostStepDoIt(track, step); // ===> Reaction is DONE
00903         }
00904         else                          // The cental region compound can be created
00905         {
00906           // First calculate the isotopic state of the parts of the compound
00907           if     (!pZ) pC =  0;
00908           else if(!pN) pC = pD;
00909           else
00910           {
00911 #ifdef pdebug
00912             G4cout<<"G4QLowEn::PSDI: pD="<<pD<<", pZ="<<pZ<<", pA="<<pA<<G4endl;
00913 #endif
00914             G4double C=pD*pZ/pA;
00915             pC=static_cast<G4int>(C); 
00916             if(G4UniformRand() < C-pC) ++pC;
00917           }
00918           if(!tZ) tC=0;
00919           else if(!tN) tC=tD;
00920           else
00921           {
00922 #ifdef pdebug
00923             G4cout<<"G4QLowEn::PSDI: tD="<<tD<<", tZ="<<tZ<<", tA="<<tA<<G4endl;
00924 #endif
00925             G4double C=tD*tZ/tA;
00926             tC=static_cast<G4int>(C); 
00927             if(G4UniformRand() < C-tC) ++tC;
00928           }
00929           // calculate the transferred momentum
00930           G4ThreeVector pFM(0.,0.,0.);
00931           if(pD < pB)                    // The projectile nucleus must be splitted
00932           {
00933             G4int nc=pD;
00934             if(pD+pD>pB) nc=pB-pD;
00935             pFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
00936             for(i=1; i < nc; ++i)
00937                              pFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
00938           }
00939           G4ThreeVector tFM(0.,0.,0.);
00940           if(tD<tB)                    // The projectile nucleus must be splitted
00941           {
00942             G4int nc=pD;
00943             if(tD+tD>tB) nc=tB-tD;
00944             tFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
00945             for(i=1; i < nc; ++i)
00946                              tFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
00947           }
00948 #ifdef pdebug
00949           G4cout<<"G4QLE::PSDI:pC="<<pC<<", tC="<<tC<<", pFM="<<pFM<<", tFM="<<tFM<<G4endl;
00950 #endif
00951           // Split the projectile spectator
00952           G4int rpZ=pZ-pC;            // number of protons in the projectile spectator
00953           G4int pF_value=pD-pC;       // number of neutrons in the projectile part of comp
00954           G4int rpN=pN-pF_value;      // number of neutrons in the projectile spectator
00955           G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0); // Mass of the spectator
00956           G4ThreeVector boostV=proj4M.vect()/proj4M.e(); // Antilab Boost Vector
00957           G4double rpE=std::sqrt(rpNM*rpNM+pFM.mag2());
00958           G4LorentzVector rp4M(pFM,rpE);
00959 #ifdef pdebug
00960             G4cout<<"G4QLE::PSDI: boostV="<<boostV<<",rp4M="<<rp4M<<",pr4M="<<proj4M<<G4endl;
00961 #endif
00962           rp4M.boost(boostV);
00963 #ifdef pdebug
00964             G4cout<<"G4QLE::PSDI: After boost, rp4M="<<rp4M<<G4endl;
00965 #endif
00966           G4ParticleDefinition* theDefinition; // Prototype of projSpectatorNuclDefinition
00967           G4int rpA=rpZ+rpN;
00968           G4Track* aFraPTrack = 0;
00969           theDefinition = 0;
00970           if(rpA==1)
00971           {
00972             if(!rpZ) theDefinition = G4Neutron::Neutron();
00973             else     theDefinition = G4Proton::Proton();
00974 #ifdef pdebug
00975             G4cout<<"G4QLE::PSDI: rpA=1, rpZ"<<rpZ<<G4endl;
00976 #endif
00977           }
00978           else if(rpA==2 && rpZ==0)            // nn decay
00979           {
00980             theDefinition = aNeutron;
00981             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
00982             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
00983 #ifdef pdebug
00984             G4cout<<"G4QLE::CPS->n+n,nn="<<rp4M.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
00985 #endif
00986             if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
00987             {
00988               G4cout<<"*W*G4QLE::CPS->n+n,t="<<rp4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
00989             }
00990             G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
00991             aFraPTrack = new G4Track(pNeu, localtime, position );
00992             aFraPTrack->SetWeight(weight);                    //    weighted
00993             aFraPTrack->SetTouchableHandle(trTouchable);
00994             tt4M-=f4M;
00995 #ifdef edebug
00996             totBaN-=2;
00997             tch4M -=f4M;
00998             G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
00999 #endif
01000 #ifdef pdebug
01001             G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
01002 #endif
01003             rp4M=s4M;
01004           }
01005           else if(rpA>2 && rpZ==0)            // Z=0 decay
01006           {
01007             theDefinition = aNeutron;
01008             G4LorentzVector f4M=rp4M/rpA;     // 4mom of 1st neutron
01009 #ifdef pdebug
01010             G4cout<<"G4QLE::CPS->Nn,M="<<rp4M.m()<<" >? N*MNeutron="<<rpA*mNeutron<<G4endl;
01011 #endif
01012             for(G4int it=1; it<rpA; ++it)     // Fill (N-1) neutrons to output
01013             {
01014               G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
01015               G4Track* aNTrack = new G4Track(pNeu, localtime, position );
01016               aNTrack->SetWeight(weight);                    //    weighted
01017               aNTrack->SetTouchableHandle(trTouchable);
01018               result.push_back(aNTrack);
01019             }
01020             G4int nesc = rpA-1;
01021             tt4M-=f4M*nesc;
01022 #ifdef edebug
01023             totBaN-=nesc;
01024             tch4M -=f4M*nesc;
01025             G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01026 #endif
01027 #ifdef pdebug
01028             G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
01029 #endif
01030             rp4M=f4M;
01031           }
01032           else if(rpA==8 && rpZ==4)            // Be8 decay
01033           {
01034             theDefinition = anAlpha;
01035             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
01036             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
01037 #ifdef pdebug
01038             G4cout<<"G4QLE::CPS->A+A,mBe8="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
01039 #endif
01040             if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
01041             {
01042               G4cout<<"*W*G4QLE::CPS->A+A,t="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
01043             }
01044             G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
01045             aFraPTrack = new G4Track(pAl, localtime, position );
01046             aFraPTrack->SetWeight(weight);                    //    weighted
01047             aFraPTrack->SetTouchableHandle(trTouchable);
01048             tt4M-=f4M;
01049 #ifdef edebug
01050             totChg-=2;
01051             totBaN-=4;
01052             tch4M -=f4M;
01053             G4cout<<">>G4QLEn::PSDI:1,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01054 #endif
01055 #ifdef pdebug
01056             G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
01057 #endif
01058             rp4M=s4M;
01059           }
01060           else if(rpA==5 && (rpZ==2 || rpZ==3)) // He5/Li5 decay
01061           {
01062             theDefinition = aProton;
01063             G4double mNuc = mProt;
01064             if(rpZ==2)
01065             {
01066               theDefinition = aNeutron;
01067               mNuc = mNeut;
01068             }
01069             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc);  // 4mom of the nucleon
01070             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
01071 #ifdef pdebug
01072             G4cout<<"G4QLowE::CPS->N+A, tM5="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
01073 #endif
01074             if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
01075             {
01076               G4cout<<"*W*G4QLE::CPS->N+A,t="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
01077             }
01078             G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
01079             aFraPTrack = new G4Track(pAl, localtime, position );
01080             aFraPTrack->SetWeight(weight);                    //    weighted
01081             aFraPTrack->SetTouchableHandle(trTouchable);
01082             tt4M-=f4M;
01083 #ifdef edebug
01084             if(theDefinition == aProton) totChg-=1;
01085             totBaN-=1;
01086             tch4M -=f4M;
01087             G4cout<<">>G4QLEn::PSDI:2,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01088 #endif
01089 #ifdef pdebug
01090             G4cout<<"G4QLowEn::PSDI:-->ProjSpectN4M="<<f4M<<G4endl;
01091 #endif
01092             theDefinition = anAlpha;
01093             rp4M=s4M;
01094           }
01095           else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rpZ,rpA,0,rpZ);
01096           if(!theDefinition)
01097             {
01098               // G4cout<<"-Warning-G4QLowEn::PSDI: pDef=0, Z="<<rpZ<<", A="<<rpA<<G4endl;
01099               // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
01100               G4ExceptionDescription ed;
01101               ed << "Particle definition is a null pointer: pDef=0, Z= " << rpZ
01102                  << ", A=" << rpA << G4endl;
01103               G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01104                           JustWarning, ed);
01105             } 
01106 #ifdef edebug
01107           if(theDefinition == anAlpha)
01108           {
01109             totChg-=2;
01110             totBaN-=4;
01111           }
01112           else
01113           {
01114             totChg-=rpZ;
01115             totBaN-=rpA;
01116           }
01117           tch4M -=rp4M;
01118           G4cout<<">>G4QLEn::PSDI:3, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
01119 #endif
01120           G4DynamicParticle* rpD = new G4DynamicParticle(theDefinition, rp4M);
01121           G4Track* aNewPTrack = new G4Track(rpD, localtime, position);
01122           aNewPTrack->SetWeight(weight);//    weighted
01123           aNewPTrack->SetTouchableHandle(trTouchable);
01124           tt4M-=rp4M;
01125 #ifdef pdebug
01126           G4cout<<"G4QLowEn::PSDI:-->ProjSpectR4M="<<rp4M<<",Z="<<rpZ<<",A="<<rpA<<G4endl;
01127 #endif
01128           //
01129           // Split the target spectator
01130           G4int rtZ=tZ-tC;            // number of protons in the target spectator
01131           G4int tF_value=tD-tC;       // number of neutrons in the target part of comp
01132           G4int rtN=tN-tF_value;      // number of neutrons in the target spectator
01133           G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0); // Mass of the spectator
01134           G4double rtE=std::sqrt(rtNM*rtNM+tFM.mag2());
01135           G4LorentzVector rt4M(tFM,rtE);
01136           G4int rtA=rtZ+rtN;
01137           G4Track* aFraTTrack = 0;
01138           theDefinition = 0;
01139           if(rtA==1)
01140           {
01141             if(!rtZ) theDefinition = G4Neutron::Neutron();
01142             else     theDefinition = G4Proton::Proton();
01143 #ifdef pdebug
01144             G4cout<<"G4QLE::PSDI: rtA=1, rtZ"<<rtZ<<G4endl;
01145 #endif
01146           }
01147           else if(rtA==2 && rtZ==0)            // nn decay
01148           {
01149             theDefinition = aNeutron;
01150             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
01151             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
01152 #ifdef pdebug
01153             G4cout<<"G4QLE::CPS->n+n,nn="<<rptM.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
01154 #endif
01155             if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
01156               G4cout<<"*W*G4QLE::CPS->n+n,t="<<rt4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
01157             G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
01158             aFraPTrack = new G4Track(pNeu, localtime, position );
01159             aFraPTrack->SetWeight(weight);                    //    weighted
01160             aFraPTrack->SetTouchableHandle(trTouchable);
01161             tt4M-=f4M;
01162 #ifdef edebug
01163             totBaN-=2;
01164             tch4M -=f4M;
01165             G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01166 #endif
01167 #ifdef pdebug
01168             G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
01169 #endif
01170             rt4M=s4M;
01171           }
01172           else if(rtA>2 && rtZ==0)            // Z=0 decay
01173           {
01174             theDefinition = aNeutron;
01175             G4LorentzVector f4M=rt4M/rtA;     // 4mom of 1st neutron
01176 #ifdef pdebug
01177             G4cout<<"G4QLE::CPS->Nn,M="<<rt4M.m()<<" >? N*MNeutron="<<rtA*mNeutron<<G4endl;
01178 #endif
01179             for(G4int it=1; it<rtA; ++it)     // Fill (N-1) neutrons to output
01180             {
01181               G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
01182               G4Track* aNTrack = new G4Track(pNeu, localtime, position );
01183               aNTrack->SetWeight(weight);                    //    weighted
01184               aNTrack->SetTouchableHandle(trTouchable);
01185               result.push_back(aNTrack);
01186             }
01187             G4int nesc = rtA-1;
01188             tt4M-=f4M*nesc;
01189 #ifdef edebug
01190             totBaN-=nesc;
01191             tch4M -=f4M*nesc;
01192             G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01193 #endif
01194 #ifdef pdebug
01195             G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
01196 #endif
01197             rt4M=f4M;
01198           }
01199           else if(rtA==8 && rtZ==4)            // Be8 decay
01200           {
01201             theDefinition = anAlpha;
01202             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
01203             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
01204 #ifdef pdebug
01205             G4cout<<"G4QLE::CPS->A+A,mBe8="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
01206 #endif
01207             if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
01208             {
01209               G4cout<<"*W*G4QLE::CTS->A+A,t="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
01210             }
01211             G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
01212             aFraTTrack = new G4Track(pAl, localtime, position );
01213             aFraTTrack->SetWeight(weight);                        // weighted
01214             aFraTTrack->SetTouchableHandle(trTouchable);
01215             tt4M-=f4M;
01216 #ifdef edebug
01217             totChg-=2;
01218             totBaN-=4;
01219             tch4M -=f4M;
01220             G4cout<<">>G4QLEn::PSDI:4,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01221 #endif
01222 #ifdef pdebug
01223             G4cout<<"G4QLowEn::PSDI:-->TargSpectA4M="<<f4M<<G4endl;
01224 #endif
01225             rt4M=s4M;
01226           }
01227           else if(rtA==5 && (rtZ==2 || rtZ==3)) // He5/Li5 decay
01228           {
01229             theDefinition = aProton;
01230             G4double mNuc = mProt;
01231             if(rpZ==2)
01232             {
01233               theDefinition = aNeutron;
01234               mNuc = mNeut;
01235             }
01236             G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc);  // 4mom of the nucleon
01237             G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
01238 #ifdef pdebug
01239             G4cout<<"G4QLowE::CPS->N+A, tM5="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
01240 #endif
01241             if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
01242             {
01243               G4cout<<"*W*G4QLE::CTS->N+A,t="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
01244             }
01245             G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
01246             aFraTTrack = new G4Track(pAl, localtime, position );
01247             aFraTTrack->SetWeight(weight);                        // weighted
01248             aFraTTrack->SetTouchableHandle(trTouchable);
01249             tt4M-=f4M;
01250 #ifdef edebug
01251             if(theDefinition == aProton) totChg-=1;
01252             totBaN-=1;
01253             tch4M -=f4M;
01254             G4cout<<">>G4QLEn::PSDI:5,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
01255 #endif
01256 #ifdef pdebug
01257             G4cout<<"G4QLowEn::PSDI:-->TargSpectN4M="<<f4M<<G4endl;
01258 #endif
01259             theDefinition = anAlpha;
01260             rt4M=s4M;
01261           }
01262           else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rtZ,rtA,0,rtZ);
01263           if(!theDefinition)
01264           {
01265             // G4cout<<"-Warning-G4QLowEn::PSDI: tDef=0, Z="<<rtZ<<", A="<<rtA<<G4endl;
01266             // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
01267             G4ExceptionDescription ed;
01268             ed << "Particle definition is a null pointer: tDef=0, Z= " << rtZ
01269                << ", A=" << rtA << G4endl;
01270             G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01271                         JustWarning, ed);
01272           }
01273 #ifdef edebug
01274           if(theDefinition == anAlpha)
01275           {
01276             totChg-=2;
01277             totBaN-=4;
01278           }
01279           else
01280           {
01281             totChg-=rtZ;
01282             totBaN-=rtA;
01283           }
01284           tch4M -=rt4M;
01285           G4cout<<">>G4QLEn::PSDI:6, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
01286 #endif
01287           G4DynamicParticle* rtD = new G4DynamicParticle(theDefinition, rt4M);
01288           G4Track* aNewTTrack = new G4Track(rtD, localtime, position);
01289           aNewTTrack->SetWeight(weight);                         // weighted
01290           aNewTTrack->SetTouchableHandle(trTouchable);
01291           tt4M-=rt4M;
01292 #ifdef pdebug
01293           G4cout<<"G4QLowEn::PSDI:-->TargSpectR4M="<<rt4M<<",Z="<<rtZ<<",A="<<rtA<<G4endl;
01294 #endif
01295           if(aFraPTrack) result.push_back(aFraPTrack);
01296           if(aNewPTrack) result.push_back(aNewPTrack);
01297           if(aFraTTrack) result.push_back(aFraTTrack);
01298           if(aNewTTrack) result.push_back(aNewTTrack);
01299           tcM = tt4M.m();            // total CMS mass of the compound (after reduction!)
01300           G4int sN=tF_value+pF_value;
01301           G4int sZ=tC+pC;
01302           tnM = targQPDG.GetNuclMass(sZ,sN,0); // GSM
01303 #ifdef pdebug
01304           G4cout<<"G4QLEn::PSDI:At#"<<nAt<<",totM="<<tcM<<",gsM="<<tnM<<",Z="<<sZ<<",N="
01305                 <<sN<<G4endl;
01306 #endif
01307           if(tcM > tnM)              // !! The totalresidual  reaction is modified !
01308           {
01309             pZ=pC;
01310             pN=pF_value;
01311             tZ=tC;
01312             tN=tF_value;
01313             tot4M=tt4M;
01314           }
01315           //else
01316           //{
01317           //  G4cout<<"***G4QLowEnergy::PostStepDoIt: totM="<<tcM<<" <= GSM="<<tnM<<G4endl;
01318           //  throw G4QException("G4QLowEnergy::PostStepDoIt: ResibualNucl below GSM shell");
01319           //}
01320         }
01321       }                               // At least one is splitted
01322       else if(tD==tB || pD==pB)       // Total absorption
01323       {
01324         tD=tZ+tN;
01325         pD=pZ+pN;
01326         tcM=tnM+1.;
01327       }
01328     }
01329     else                              // Total compound (define tD to go out of the while)
01330     {
01331       tD=tZ+tN;
01332       pD=pZ+pN;
01333       tcM=tnM+1.;
01334     }
01335   } // End of the interaction WHILE
01336   G4double totM=tot4M.m();            // total CMS mass of the reaction
01337   G4int totN=tN+pN;
01338   G4int totZ=tZ+pZ;
01339 #ifdef edebug
01340   G4cout<<">>>G4QLEn::PSDI: dZ="<<totChg-totZ<<", dB="<<totBaN-totN-totZ<<", d4M="
01341         <<tch4M-tot4M<<",N="<<totN<<",Z="<<totZ<<G4endl;
01342 #endif
01343   // @@ Here mass[i] can be calculated if mass=0
01344   G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
01345                       0.,0.,0.,0.,0.,0.};
01346   mass[0] = tM = targQPDG.GetNuclMass(totZ, totN, 0); // gamma+gamma and GSM update
01347 #ifdef pdebug
01348   G4cout<<"G4QLEn::PSDI:TotM="<<totM<<",NucM="<<tM<<",totZ="<<totZ<<",totN="<<totN<<G4endl;
01349 #endif
01350   if (totN>0 && totZ>0)
01351   {
01352     mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron
01353     mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton
01354   }
01355   if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d
01356   if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t
01357   if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3
01358   if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a
01359   if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ  ,totN-2,0); //n+n
01360     mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron)
01361   if ( totZ > 2 )             mass[9] = targQPDG.GetNuclMass(totZ-2,totN  ,0); //p+p
01362     mass[10] = mass[5]; // proton+deuteron (the same as He3)
01363     mass[11] = mass[4]; // neutron+deuteron (the same as t)
01364     mass[12] = mass[6]; // deuteron+deuteron (the same as alpha)
01365     mass[13] = mass[6]; // proton+tritium (the same as alpha)
01366   if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);//n+t
01367   if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);//He3+p
01368     mass[16] = mass[6]; // neutron+He3 (the same as alpha)
01369   if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);//pa
01370   if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);//na
01371   if(pL>0) // @@ Not debugged @@
01372   {
01373     G4int pL1=pL-1;
01374     if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ  ,totN  ,pL1);// Lambda+gamma
01375     if( (totN > 0 && totZ > 0) ||        totZ > 1 ) 
01376       mass[20]=targQPDG.GetNuclMass(totZ-1,totN  ,pL1);//Lp
01377     if( (totN > 0 && totZ > 0) || totN > 0        ) 
01378       mass[21]=targQPDG.GetNuclMass(totZ  ,totN-1,pL1);//Ln
01379     if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) ) 
01380       mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld
01381     if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) ) 
01382       mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt
01383     if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) ) 
01384       mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3
01385     if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) ) 
01386       mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La
01387   }
01388 #ifdef debug
01389   G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl;
01390 #endif
01391   tA=tZ+tN; 
01392   pA=pZ+pN; 
01393   G4double prZ=pZ/pA+tZ/tA;
01394   G4double prN=pN/pA+tN/tA;
01395   G4double prD=prN*prZ;
01396   G4double prA=prD*prD;
01397   G4double prH=prD*prZ;
01398   G4double prT=prD*prN;
01399   G4double fhe3=6.*std::sqrt(tA);
01400   G4double prL=0.;
01401   if(pL>0) prL=pL/pA;
01402   G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
01403                       0.,0.,0.,0.,0.,0.};
01404   qval[ 0] = (totM - mass[ 0])/137./137.;
01405   qval[ 1] = (totM - mass[ 1] - mNeut)*prN/137.;
01406   qval[ 2] = (totM - mass[ 2] - mProt)*prZ/137.;
01407   qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3./137.;
01408   qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6./137.;
01409   qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3/137.;
01410   qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9./137.;
01411   qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN;
01412   qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD;
01413   qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ;
01414   qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.;
01415   qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.;
01416   qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.;
01417   qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.;
01418   qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.;
01419   qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3;
01420   qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3;
01421   qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.;
01422   qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.;
01423   if(pZ>0)
01424   {
01425     qval[19] = (totM - mass[19] - mLamb)*prL;
01426     qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ;
01427     qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN;
01428     qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.;
01429     qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.;
01430     qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3;
01431     qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4;
01432   }
01433 #ifdef debug
01434   G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<", prA="<<pA<<G4endl;
01435 #endif
01436   
01437   G4double qv = 0.0;                        // Total sum of probabilities (q-values)
01438   for(i=0; i<nCh; ++i )
01439   {
01440 #ifdef sdebug
01441     G4cout<<"G4QLowEn::PSDI: i="<<i<<", q="<<qval[i]<<",m="<<mass[i]<<G4endl;
01442 #endif
01443     if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels (mesons)
01444     if( qval[i] < 0.0 )      qval[i] = 0.0; // Close the splitting impossible channels
01445     qv += qval[i];
01446   }
01447   // Select the channel
01448   G4double qv1 = 0.0;
01449   G4double ran = qv*G4UniformRand();
01450   G4int index  = 0;
01451   for( index=0; index<nCh; ++index ) if( qval[index] > 0.0 )
01452   {
01453     qv1 += qval[index];
01454     if( ran <= qv1 ) break;
01455   }
01456 #ifdef debug
01457   G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl;
01458 #endif
01459   if(index == nCh)
01460   {
01461     G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl;
01462     G4cout<<"G4QLowEn::PoStDI:Reaction "<<projPDG<<"+"<<targPDG<<", P="<<Momentum<<G4endl;
01463     G4int nRes=result.size();
01464     if(nRes)G4cout<<"G4QLowEnergy::PoStDI: result exists with N="<<nRes<<" tracks"<<G4endl;
01465 //    else throw G4QException("***G4QLowEnergy::PostStepDoIt: Can't decay ResidualCompound");
01466     else {
01467       G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01468                   FatalException, "Can't decay ResidualCompound"); 
01469     }
01470     G4double minEx=1000000.;                       // Prototype of minimal excess
01471     G4bool found = false;                          // The solution is not found yet
01472     G4int cInd = 0;                                // Prototype of the best index
01473     G4int ctN  = 0;                                // To
01474     G4int ctZ  = 0;                                //    avoid
01475     G4LorentzVector c4M=tot4M;                     //          recalculation
01476     G4double ctM=0.;                               //                        of found
01477     for(G4int it=0; it<nRes; ++it)
01478     {
01479       G4Track* iTrack = result[it];
01480       const G4DynamicParticle* iHadron = iTrack->GetDynamicParticle();
01481       G4ParticleDefinition* iParticle=iHadron->GetDefinition();
01482       G4int iPDG = iParticle->GetPDGEncoding();
01483       G4LorentzVector ih4M = projHadron->Get4Momentum();
01484       G4cout<<"     G4QLowEn::PoStDI: H["<<it<<"] = [ "<<iPDG<<", "<<ih4M<<" ]"<<G4endl;
01485       G4int iB = iParticle->GetBaryonNumber();     // A of the secondary
01486       G4int iC = G4int(iParticle->GetPDGCharge()); // Z of the secondary
01487       G4LorentzVector new4M = tot4M + ih4M;        // Make temporary tot4M
01488       G4int ntN=totN + iB-iC;                      // Make temporary totN
01489       G4int ntZ=totZ + iC;                         // Make temporary totZ
01490       G4double ntM = targQPDG.GetNuclMass(ntZ, ntN, 0); // Make temporary GSM
01491       G4double ttM = new4M.m();
01492       if(ttM > ntM) // >>> This is at least one possible solution ! <<<
01493       {
01494         G4double nEx = ttM - ntM;
01495         if(found && nEx < minEx) // This one is better than the previous found
01496         {
01497           cInd = it;
01498           minEx= nEx;
01499           ctN  = ntN;
01500           ctZ  = ntZ;
01501           c4M  = new4M;
01502           ctM  = ttM;
01503             mass[0] = ntM;
01504         }
01505         found = true;
01506       }
01507     }
01508     if(found)
01509     {
01510       tot4M = c4M;
01511       totM  = ctM;
01512       totN  = ctN;
01513       totZ  = ctZ;
01514       delete result[cInd];
01515       G4int nR1 = nRes -1;
01516       if(cInd < nR1) result[cInd] = result[nR1];
01517       result.pop_back();
01518       //nRes=nR1;                            // @@ can be used for the two pt correction
01519       index = 0;                             // @@ emergency gamma+gamma decay is selected
01520     }
01521     else
01522     {
01523       G4cout<<"***G4QLowEnergy::PoStDI: One hadron coddection did not succeed***"<<G4endl;
01524       if(nRes>1) G4cout<<"***G4QLowEnergy::PoStDI:nRes's big enough to use 2PtCor"<<G4endl;
01525       // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't correct by one particle.");
01526       G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01527                   FatalException, "Can't correct by one particle");
01528     }
01529   }
01530   // @@ Convert it to G4QHadrons    
01531   G4DynamicParticle* ResSec = new G4DynamicParticle;
01532   G4DynamicParticle* FstSec = new G4DynamicParticle;
01533   G4DynamicParticle* SecSec = new G4DynamicParticle;
01534 #ifdef debug
01535   G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl;
01536 #endif
01537 
01538   G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype
01539   G4double mF=0.;
01540   G4double mS=0.;
01541   G4int    rA=totZ+totN;
01542   G4int    rZ=totZ;
01543   G4int    rL=pL;
01544   G4int complete=3;
01545   G4ParticleDefinition* theDefinition;  // Prototype for qfNucleon
01546   switch( index )
01547   {
01548    case 0:
01549      if(!evaporate || rA<2)
01550      {
01551        if(!rZ) theDefinition=aNeutron;
01552        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01553        if(!theDefinition)
01554        {
01555          // G4cerr<<"-Warning-G4LE::PSDI: notDef(0), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01556          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01557          G4ExceptionDescription ed;
01558          ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
01559             << ", A=" << rA << ", L=" << rL << G4endl;
01560          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01561                         JustWarning, ed);
01562        }
01563        ResSec->SetDefinition( theDefinition );
01564        FstSec->SetDefinition( aGamma );
01565        SecSec->SetDefinition( aGamma );
01566      }
01567      else
01568      {
01569        delete ResSec;
01570        delete FstSec;
01571        delete SecSec;
01572        complete=0;
01573      }
01574      break;
01575    case 1:
01576      rA-=1;                                              // gamma+n
01577      if(!evaporate || rA<2)
01578      {
01579        if(!rZ) theDefinition=aNeutron;
01580        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01581        if(!theDefinition)
01582        {
01583          // G4cerr<<"-Warning-G4LE::PSDI: notDef(1), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01584          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01585          G4ExceptionDescription ed;
01586          ed << "Particle definition is a null pointer: notDef(1), Z= " << rZ
01587             << ", A=" << rA << ", L=" << rL << G4endl;
01588          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0001",
01589                         JustWarning, ed);
01590        }
01591        ResSec->SetDefinition( theDefinition );
01592        SecSec->SetDefinition( aGamma );
01593      }
01594      else
01595      {
01596        delete ResSec;
01597        delete SecSec;
01598        complete=1;
01599      }
01600      FstSec->SetDefinition( aNeutron );
01601      mF=mNeut; // First hadron 4-momentum
01602      break;
01603    case 2:
01604      rA-=1;
01605      rZ-=1;                                             // gamma+p
01606      if(!evaporate || rA<2)
01607      {
01608        if(!rZ) theDefinition=aNeutron;
01609        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01610        if(!theDefinition)
01611        {
01612          // G4cerr<<"-Warning-G4LE::PSDI: notDef(2), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01613          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01614           G4ExceptionDescription ed;
01615           ed << "Particle definition is a null pointer: notDef(2), Z= " << rZ
01616              << ", A=" << rA << ", L="<< rL << G4endl;
01617           G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0002",
01618                       JustWarning, ed);
01619        }
01620        ResSec->SetDefinition( theDefinition );
01621        SecSec->SetDefinition( aGamma );
01622      }
01623      else
01624      {
01625        delete ResSec;
01626        delete SecSec;
01627        complete=1;
01628      }
01629      FstSec->SetDefinition( aProton );
01630      mF=mProt; // First hadron 4-momentum
01631      break;
01632    case 3:
01633      rA-=2;
01634      rZ-=1;                                             // gamma+d
01635      if(!evaporate || rA<2)
01636      {
01637        if(!rZ) theDefinition=aNeutron;
01638        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01639        if(!theDefinition)
01640        {
01641          // G4cerr<<"-Warning-G4LE::PSDI: notDef(3), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01642          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01643          G4ExceptionDescription ed;
01644          ed << "Particle definition is a null pointer: notDef(3), Z= " << rZ
01645             << ", A=" << rA << ", L=" << rL << G4endl;
01646          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0003",
01647                         JustWarning, ed);
01648        }
01649        ResSec->SetDefinition( theDefinition );
01650        SecSec->SetDefinition( aGamma );
01651      }
01652      else
01653      {
01654        delete ResSec;
01655        delete SecSec;
01656        complete=1;
01657      }
01658      FstSec->SetDefinition( aDeuteron );
01659      mF=mDeut; // First hadron 4-momentum
01660      break;
01661    case 4:
01662      rA-=3;                                             // gamma+t
01663      rZ-=1;
01664      if(!evaporate || rA<2)
01665      {
01666        if(!rZ) theDefinition=aNeutron;
01667        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01668        if(!theDefinition)
01669        {
01670          // G4cerr<<"-Warning-G4LE::PSDI: notDef(4), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01671          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01672          G4ExceptionDescription ed;
01673          ed << "Particle definition is a null pointer: notDef(4), Z= " << rZ
01674             << ", A=" << rA << ", L=" << rL << G4endl;
01675          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0004",
01676                      JustWarning, ed);
01677        }
01678        ResSec->SetDefinition( theDefinition );
01679        SecSec->SetDefinition( aGamma );
01680      }
01681      else
01682      {
01683        delete ResSec;
01684        delete SecSec;
01685        complete=1;
01686      }
01687      FstSec->SetDefinition( aTriton );
01688      mF=mTrit; // First hadron 4-momentum
01689      break;
01690   case 5:                                            // gamma+He3
01691      rA-=3;
01692      rZ-=2;
01693      if(!evaporate || rA<2)
01694      {
01695        if(!rZ) theDefinition=aNeutron;
01696        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01697        if(!theDefinition)
01698        {
01699          // G4cerr<<"-Warning-G4LE::PSDI: notDef(5), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01700          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01701          G4ExceptionDescription ed;
01702          ed << "Particle definition is a null pointer: notDef(5), Z= " << rZ
01703             << ", A=" << rA << ", L=" << rL << G4endl;
01704          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0005",
01705                      JustWarning, ed);
01706        }
01707        ResSec->SetDefinition( theDefinition );
01708        SecSec->SetDefinition( aGamma );
01709      }
01710      else
01711      {
01712        delete ResSec;
01713        delete SecSec;
01714        complete=1;
01715      }
01716      FstSec->SetDefinition( aHe3);
01717      mF=mHel3; // First hadron 4-momentum
01718      break;
01719    case 6:
01720      rA-=4;
01721      rZ-=2;                                         // gamma+He4
01722      if(!evaporate || rA<2)
01723      {
01724        if(!rZ) theDefinition=aNeutron;
01725        else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01726        if(!theDefinition)
01727        {
01728          // G4cerr<<"-Warning-G4LE::PSDI: notDef(6), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01729          // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01730          G4ExceptionDescription ed;
01731          ed << "Particle definition is a null pointer: notDef(6), Z= " << rZ
01732             << ", A=" << rA << ", L=" << rL << G4endl;
01733          G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0006",
01734                      JustWarning, ed);
01735        }
01736        ResSec->SetDefinition( theDefinition );
01737        SecSec->SetDefinition( aGamma );
01738      }
01739      else
01740      {
01741        delete ResSec;
01742        delete SecSec;
01743        complete=1;
01744      }
01745      FstSec->SetDefinition( anAlpha );
01746      mF=mAlph; // First hadron 4-momentum
01747      break;
01748    case 7:
01749      rA-=2;                                          // n+n
01750      if(rA==1 && !rZ) theDefinition=aNeutron;
01751      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01752      if(!theDefinition)
01753      {
01754        // G4cerr<<"-Warning-G4LE::PSDI: notDef(7), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01755        // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
01756        G4ExceptionDescription ed;
01757        ed << "Particle definition is a null pointer: notDef(7), Z= " << rZ
01758           << ", A=" << rA << ", L=" << rL << G4endl;
01759        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0007",
01760                    JustWarning, ed);
01761      }
01762      ResSec->SetDefinition( theDefinition );
01763      FstSec->SetDefinition( aNeutron );
01764      SecSec->SetDefinition( aNeutron );
01765      mF=mNeut; // First hadron 4-momentum
01766      mS=mNeut; // Second hadron 4-momentum
01767      break;
01768    case 8:
01769      rZ-=1;
01770      rA-=2;                                           // n+p
01771      if(rA==1 && !rZ) theDefinition=aNeutron;
01772      else if(rA==2 && !rZ)
01773      {
01774        index=7;
01775        ResSec->SetDefinition( aDeuteron);
01776        FstSec->SetDefinition( aNeutron );
01777        SecSec->SetDefinition( aNeutron );
01778        mF=mNeut; // First hadron 4-momentum
01779        mS=mNeut; // Second hadron 4-momentum
01780        break;
01781      }
01782      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01783      if(!theDefinition)
01784      {
01785        // G4cerr<<"-Warning-G4LE::PSDI: notDef(8), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01786        // throw G4QException("G4QLowEn::PostStepDoIt: Darticle Definition is a null pointer");
01787        G4ExceptionDescription ed;
01788        ed << "Particle definition is a null pointer: notDef(8), Z= " << rZ
01789           << ", A=" << rA << ", L=" << rL << G4endl;
01790        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0008",
01791                       JustWarning, ed);
01792      }
01793      ResSec->SetDefinition( theDefinition );
01794      FstSec->SetDefinition( aNeutron );
01795      SecSec->SetDefinition( aProton );
01796      mF=mNeut; // First hadron 4-momentum
01797      mS=mProt; // Second hadron 4-momentum
01798      break;
01799    case 9:
01800      rZ-=2;
01801      rA-=2;                                           // p+p
01802      if(rA==1 && !rZ) theDefinition=aNeutron;
01803      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01804      if(!theDefinition)
01805      {
01806        // G4cerr<<"-Warning-G4LE::PSDI: notDef(9), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01807        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01808        G4ExceptionDescription ed;
01809        ed << "Particle definition is a null pointer: notDef(9), Z= " << rZ
01810           << ", A=" << rA << ", L=" << rL << G4endl;
01811        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0009",
01812                    JustWarning, ed);
01813      }
01814      ResSec->SetDefinition( theDefinition );
01815      FstSec->SetDefinition( aProton );
01816      SecSec->SetDefinition( aProton );
01817      mF=mProt; // First hadron 4-momentum
01818      mS=mProt; // Second hadron 4-momentum
01819      break;
01820    case 10:
01821      rZ-=2;
01822      rA-=3;                                            // p+d
01823      if(rA==1 && !rZ) theDefinition=aNeutron;
01824      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01825      if(!theDefinition)
01826      {
01827        // G4cerr<<"-Warning-G4LE::PSDI: notDef(10), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01828        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01829        G4ExceptionDescription ed;
01830        ed << "Particle definition is a null pointer: notDef(10), Z= " << rZ
01831           << ", A=" << rA << ", L=" << rL << G4endl;
01832        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0010",
01833                    JustWarning, ed);
01834      }
01835      ResSec->SetDefinition( theDefinition );
01836      FstSec->SetDefinition( aProton );
01837      SecSec->SetDefinition( aDeuteron );
01838      mF=mProt; // First hadron 4-momentum
01839      mS=mDeut; // Second hadron 4-momentum
01840      break;
01841    case 11:
01842      rZ-=1;
01843      rA-=3;                                            // n+d
01844      if(rA==1 && !rZ) theDefinition=aNeutron;
01845      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01846      if(!theDefinition)
01847      {
01848        // G4cerr<<"-Warning-G4LE::PSDI: notDef(11), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01849        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01850        G4ExceptionDescription ed;
01851        ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
01852           << ", A=" << rA << ", L=" << rL << G4endl;
01853        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0011",
01854                       JustWarning, ed);
01855      }
01856      ResSec->SetDefinition( theDefinition );
01857      FstSec->SetDefinition( aNeutron );
01858      SecSec->SetDefinition( aDeuteron );
01859      mF=mNeut; // First hadron 4-momentum
01860      mS=mDeut; // Second hadron 4-momentum
01861      break;
01862    case 12:
01863      rZ-=2;
01864      rA-=4;                                            // d+d
01865      if(rA==1 && !rZ) theDefinition=aNeutron;
01866      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01867      if(!theDefinition)
01868      {
01869        // G4cerr<<"-Warning-G4LE::PSDI: notDef(12), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01870        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01871        G4ExceptionDescription ed;
01872        ed << "Particle definition is a null pointer: notDef(12), Z= " << rZ
01873           << ", A=" << rA << ", L=" << rL << G4endl;
01874        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0012",
01875                    JustWarning, ed);
01876      }
01877      ResSec->SetDefinition( theDefinition );
01878      FstSec->SetDefinition( aDeuteron );
01879      SecSec->SetDefinition( aDeuteron );
01880      mF=mDeut; // First hadron 4-momentum
01881      mS=mDeut; // Second hadron 4-momentum
01882      break;
01883    case 13:
01884      rZ-=2;
01885      rA-=4;                                            // p+t
01886      if(rA==1 && !rZ) theDefinition=aNeutron;
01887      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01888      if(!theDefinition)
01889      {
01890        // G4cerr<<"-Warning-G4LE::PSDI: notDef(13), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01891        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01892        G4ExceptionDescription ed;
01893        ed << "Particle definition is a null pointer: notDef(13), Z= " << rZ
01894           << ", A=" << rA << ", L=" << rL << G4endl;
01895        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0013",
01896                    JustWarning, ed);
01897      }
01898      ResSec->SetDefinition( theDefinition );
01899      FstSec->SetDefinition( aProton );
01900      SecSec->SetDefinition( aTriton );
01901      mF=mProt; // First hadron 4-momentum
01902      mS=mTrit; // Second hadron 4-momentum
01903      break;
01904    case 14:
01905      rZ-=1;
01906      rA-=4;                                             // n+t
01907      if(rA==1 && !rZ) theDefinition=aNeutron;
01908      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01909      if(!theDefinition)
01910      {
01911        // G4cerr<<"-Warning-G4LE::PSDI: notDef(14), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01912        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01913        G4ExceptionDescription ed;
01914        ed << "Particle definition is a null pointer: notDef(14), Z= " << rZ
01915           << ", A=" << rA << ", L=" << rL << G4endl;
01916        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0014",
01917                    JustWarning, ed);
01918      }
01919      ResSec->SetDefinition( theDefinition );
01920      FstSec->SetDefinition( aNeutron );
01921      SecSec->SetDefinition( aTriton );
01922      mF=mNeut; // First hadron 4-momentum
01923      mS=mTrit; // Second hadron 4-momentum
01924      break;
01925    case 15:
01926      rZ-=3;
01927      rA-=4;                                             // p+He3
01928      if(rA==1 && !rZ) theDefinition=aNeutron;
01929      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01930      if(!theDefinition)
01931      {
01932        // G4cerr<<"-Warning-G4LE::PSDI: notDef(15), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01933        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01934        G4ExceptionDescription ed;
01935        ed << "Particle definition is a null pointer: notDef(15), Z= " << rZ
01936           << ", A=" << rA << ", L=" << rL << G4endl;
01937        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0015",
01938                    JustWarning, ed);
01939      }
01940      ResSec->SetDefinition( theDefinition );
01941      FstSec->SetDefinition( aProton);
01942      SecSec->SetDefinition( aHe3 );
01943      mF=mProt; // First hadron 4-momentum
01944      mS=mHel3; // Second hadron 4-momentum
01945      break;
01946    case 16:
01947      rZ-=2;
01948      rA-=4;                                             // n+He3
01949      if(rA==1 && !rZ) theDefinition=aNeutron;
01950      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01951      if(!theDefinition)
01952      {
01953        // G4cerr<<"-Warning-G4LE::PSDI: notDef(16), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01954        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01955        G4ExceptionDescription ed;
01956        ed << "Particle definition is a null pointer: notDef(16), Z= " << rZ
01957           << ", A=" << rA << ", L=" << rL << G4endl;
01958        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0016",
01959                       JustWarning, ed);
01960      }
01961      ResSec->SetDefinition( theDefinition );
01962      FstSec->SetDefinition( aNeutron );
01963      SecSec->SetDefinition( aHe3 );
01964      mF=mNeut; // First hadron 4-momentum
01965      mS=mHel3; // Second hadron 4-momentum
01966      break;
01967    case 17:
01968      rZ-=3;
01969      rA-=5;                                            // p+alph
01970      if(rA==1 && !rZ) theDefinition=aNeutron;
01971      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01972      if(!theDefinition)
01973      {
01974        // G4cerr<<"-Warning-G4LE::PSDI: notDef(17), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01975        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01976        G4ExceptionDescription ed;
01977        ed << "Particle definition is a null pointer: notDef(17), Z= " << rZ
01978           << ", A=" << rA << ", L=" << rL << G4endl;
01979        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0017",
01980                    JustWarning, ed);
01981      }
01982      ResSec->SetDefinition( theDefinition );
01983      FstSec->SetDefinition( aProton );
01984      SecSec->SetDefinition( anAlpha );
01985      mF=mProt; // First hadron 4-momentum
01986      mS=mAlph; // Second hadron 4-momentum
01987      break;
01988    case 18:
01989      rZ-=2;
01990      rA-=5;                                              // n+alph
01991      if(rA==1 && !rZ) theDefinition=aNeutron;
01992      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01993      if(!theDefinition)
01994      {
01995        // G4cerr<<"-Warning-G4LE::PSDI: notDef(18), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
01996        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
01997        G4ExceptionDescription ed;
01998        ed << "Particle definition is a null pointer: notDef(18), Z= " << rZ
01999           << ", A=" << rA << ", L=" << rL << G4endl;
02000        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0018",
02001                    JustWarning, ed);
02002      }
02003      ResSec->SetDefinition( theDefinition );
02004      FstSec->SetDefinition( aNeutron );
02005      SecSec->SetDefinition( anAlpha );
02006      mF=mNeut; // First hadron 4-momentum
02007      mS=mAlph; // Second hadron 4-momentum
02008      break;
02009    case 19:
02010      rL-=1;                                              // L+gamma (@@ rA inludes rL?)
02011      rA-=1;
02012      if(rA==1 && !rZ) theDefinition=aNeutron;
02013      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02014      if(!theDefinition)
02015      {
02016        // G4cerr<<"-Warning-G4LE::PSDI: notDef(19), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02017        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
02018        G4ExceptionDescription ed;
02019        ed << "Particle definition is a null pointer: notDef(19), Z= " << rZ
02020           << ", A=" << rA << ", L=" << rL << G4endl;
02021        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0019",
02022                       JustWarning, ed);
02023      }
02024      ResSec->SetDefinition( theDefinition );
02025      FstSec->SetDefinition( aLambda );
02026      SecSec->SetDefinition( aGamma );
02027      mF=mLamb; // First hadron 4-momentum
02028      break;
02029    case 20:
02030      rL-=1;                                              // L+p (@@ rA inludes rL?)
02031      rZ-=1;
02032      rA-=2;
02033      if(rA==1 && !rZ) theDefinition=aNeutron;
02034      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02035      if(!theDefinition)
02036      {
02037        // G4cerr<<"-Warning-G4LE::PSDI: notDef(20), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02038        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
02039        G4ExceptionDescription ed;
02040        ed << "Particle definition is a null pointer: notDef(20), Z= " << rZ
02041           << ", A=" << rA << ", L=" << rL << G4endl;
02042        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0020",
02043                       JustWarning, ed);
02044      }
02045      ResSec->SetDefinition( theDefinition );
02046      FstSec->SetDefinition( aProton );
02047      SecSec->SetDefinition( aLambda );
02048      mF=mProt; // First hadron 4-momentum
02049      mS=mLamb; // Second hadron 4-momentum
02050      break;
02051    case 21:
02052      rL-=1;                                              // L+n (@@ rA inludes rL?)
02053      rA-=2;
02054      if(rA==1 && !rZ) theDefinition=aNeutron;
02055      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02056      if(!theDefinition)
02057      {
02058        // G4cerr<<"-Warning-G4LE::PSDI: notDef(21), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02059        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Pefinition is a null pointer");
02060        G4ExceptionDescription ed;
02061        ed << "Particle definition is a null pointer: notDef(21), Z= " << rZ
02062           << ", A=" << rA << ", L=" << rL << G4endl;
02063        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0021",
02064                       JustWarning, ed);
02065      }
02066      ResSec->SetDefinition( theDefinition );
02067      FstSec->SetDefinition( aNeutron );
02068      SecSec->SetDefinition( aLambda );
02069      mF=mNeut; // First hadron 4-momentum
02070      mS=mLamb; // Second hadron 4-momentum
02071      break;
02072    case 22:
02073      rL-=1;                                              // L+d (@@ rA inludes rL?)
02074      rZ-=1;
02075      rA-=3;
02076      if(rA==1 && !rZ) theDefinition=aNeutron;
02077      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02078      if(!theDefinition)
02079      {
02080        // G4cerr<<"-Warning-G4LE::PSDI: notDef(22), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02081        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
02082        G4ExceptionDescription ed;
02083        ed << "Particle definition is a null pointer: notDef(22), Z= " << rZ
02084           << ", A=" << rA << ", L=" << rL << G4endl;
02085        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0022",
02086                       JustWarning, ed);
02087      }
02088     ResSec->SetDefinition( theDefinition );
02089      FstSec->SetDefinition( aDeuteron );
02090      SecSec->SetDefinition( aLambda );
02091      mF=mDeut; // First hadron 4-momentum
02092      mS=mLamb; // Second hadron 4-momentum
02093      break;
02094    case 23:
02095      rL-=1;                                              // L+t (@@ rA inludes rL?)
02096      rZ-=1;
02097      rA-=4;
02098      if(rA==1 && !rZ) theDefinition=aNeutron;
02099      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02100      if(!theDefinition)
02101      {
02102        // G4cerr<<"-Warning-G4LE::PSDI: notDef(23), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02103        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
02104        G4ExceptionDescription ed;
02105        ed << "Particle definition is a null pointer: notDef(23), Z= " << rZ
02106           << ", A=" << rA << ", L=" << rL << G4endl;
02107        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0023",
02108                       JustWarning, ed);
02109      }
02110     ResSec->SetDefinition( theDefinition );
02111      FstSec->SetDefinition( aTriton );
02112      SecSec->SetDefinition( aLambda );
02113      mF=mTrit; // First hadron 4-momentum
02114      mS=mLamb; // Second hadron 4-momentum
02115      break;
02116    case 24:
02117      rL-=1;                                              // L+He3 (@@ rA inludes rL?)
02118      rZ-=2;
02119      rA-=4;
02120      if(rA==1 && !rZ) theDefinition=aNeutron;
02121      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02122      if(!theDefinition) 
02123      {
02124        // G4cerr<<"-Warning-G4LE::PSDI: notDef(24), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02125        // throw G4QException("G4QLowEn::PostStepDoIt: particle definition is a null pointer");
02126        G4ExceptionDescription ed;
02127        ed << "Particle definition is a null pointer: notDef(24), Z= " << rZ
02128           << ", A=" << rA << ", L=" << rL << G4endl;
02129        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0024",
02130                       JustWarning, ed);
02131      }
02132      ResSec->SetDefinition( theDefinition );
02133      FstSec->SetDefinition( aHe3 );
02134      SecSec->SetDefinition( aLambda );
02135      mF=mHel3; // First hadron 4-momentum
02136      mS=mLamb; // Second hadron 4-momentum
02137      break;
02138    case 25:
02139      rL-=1;                                              // L+alph (@@ rA inludes rL?)
02140      rZ-=2;
02141      rA-=5;
02142      if(rA==1 && !rZ) theDefinition=aNeutron;
02143      else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
02144      if(!theDefinition)
02145      {
02146        // G4cerr<<"-Warning-G4LE::PSDI: notDef(25), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
02147        // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
02148        G4ExceptionDescription ed;
02149        ed << "Particle definition is a null pointer: notDef(25), Z= " << rZ
02150           << ", A=" << rA << ", L=" << rL << G4endl;
02151        G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0025",
02152                       JustWarning, ed);
02153      }
02154      ResSec->SetDefinition( theDefinition );
02155      FstSec->SetDefinition( anAlpha );
02156      SecSec->SetDefinition( aLambda );
02157      mF=mAlph; // First hadron 4-momentum
02158      mS=mLamb; // Second hadron 4-momentum
02159      break;
02160   }
02161 #ifdef debug
02162   G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl;
02163 #endif
02164   G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum
02165   G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum
02166   G4LorentzVector dir4Mom=tot4M;       // Prototype of the resN decay direction 4-momentum
02167   dir4Mom.setE(tot4M.e()/2.);          // Get half energy and total 3-momentum
02168   // @@ Can be repeated to take into account the Coulomb Barrier
02169   if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp))
02170   {
02171     // G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2="
02172     //  <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl;
02173     // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound");
02174     G4ExceptionDescription ed;
02175     ed << "Can't decay the Compound: i=" << index << ",tM=" << totM << "->M1="
02176        << res4Mom.m() << "+M2=" << fst4Mom.m() << "+M3=" << snd4Mom.m() << "=="
02177        << res4Mom.m()+fst4Mom.m()+snd4Mom.m() << G4endl;
02178     G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
02179                 FatalException, ed);
02180   }
02181 #ifdef debug
02182   G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl;
02183 #endif
02184   G4Track* aNewTrack = 0;
02185   if(complete)
02186   {
02187     FstSec->Set4Momentum(fst4Mom);
02188     aNewTrack = new G4Track(FstSec, localtime, position );
02189     aNewTrack->SetWeight(weight);                                   //    weighted
02190     aNewTrack->SetTouchableHandle(trTouchable);
02191     result.push_back( aNewTrack );
02192     EnMomConservation-=fst4Mom;
02193 #ifdef debug
02194     G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom
02195           <<", PDG="<<FstSec->GetDefinition()->GetPDGEncoding()<<G4endl;
02196 #endif
02197     if(complete>2)                     // Final solution
02198     {
02199       ResSec->Set4Momentum(res4Mom);
02200       aNewTrack = new G4Track(ResSec, localtime, position );
02201       aNewTrack->SetWeight(weight);                                   //    weighted
02202       aNewTrack->SetTouchableHandle(trTouchable);
02203       result.push_back( aNewTrack );
02204       EnMomConservation-=res4Mom;
02205 #ifdef debug
02206       G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL="
02207             <<rL<<G4endl;
02208 #endif
02209       SecSec->Set4Momentum(snd4Mom);
02210       aNewTrack = new G4Track(SecSec, localtime, position );
02211       aNewTrack->SetWeight(weight);                                   //    weighted
02212       aNewTrack->SetTouchableHandle(trTouchable);
02213       result.push_back( aNewTrack );
02214       EnMomConservation-=snd4Mom;
02215 #ifdef debug
02216       G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom
02217           <<", PDG="<<SecSec->GetDefinition()->GetPDGEncoding()<<G4endl;
02218 #endif
02219     }
02220     else res4Mom+=snd4Mom;
02221   }
02222   else res4Mom=tot4M;
02223   if(complete<3)                        // Evaporation of the residual must be done
02224   {
02225     G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus
02226     G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!)
02227     Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear !
02228     G4int nOut=evaHV->size();
02229     for(i=0; i<nOut; i++)
02230     {
02231       G4QHadron* curH = (*evaHV)[i];
02232       G4int hPDG=curH->GetPDGCode();
02233       G4LorentzVector h4Mom=curH->Get4Momentum();
02234       EnMomConservation-=h4Mom;
02235 #ifdef debug
02236       G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl;
02237 #endif
02238       if     (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron;
02239       else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton;
02240       else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda;
02241       else if(hPDG==     22 )               theDefinition = aGamma;
02242       else if(hPDG==     111)               theDefinition = aPiZero;
02243       else if(hPDG==     211)               theDefinition = aPiPlus;
02244       else if(hPDG==    -211)               theDefinition = aPiMinus;
02245       else
02246       {
02247         G4int hZ=curH->GetCharge();
02248         G4int hA=curH->GetBaryonNumber();
02249         G4int hS=curH->GetStrangeness();
02250         theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion
02251       }
02252       if(theDefinition)
02253       {
02254         G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom);
02255         G4Track* evaQH = new G4Track(theEQH, localtime, position );
02256         evaQH->SetWeight(weight);                                   //    weighted
02257         evaQH->SetTouchableHandle(trTouchable);
02258         result.push_back( evaQH );         
02259       }
02260       else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl;
02261     }
02262   }
02263   // algorithm implementation --- STOPS HERE
02264   G4int nres=result.size();
02265   aParticleChange.SetNumberOfSecondaries(nres);
02266   for(i=0; i<nres; ++i) aParticleChange.AddSecondary(result[i]);
02267 #ifdef debug
02268   G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
02269 #endif
02270   return G4VDiscreteProcess::PostStepDoIt(track, step);
02271 }
02272 
02273 G4double G4QLowEnergy::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG) 
02274 {
02275   static G4bool first=true;
02276   static G4VQCrossSection* CSmanager;
02277   if(first)                              // Connection with a singletone
02278   {
02279     CSmanager=G4QIonIonCrossSection::GetPointer();
02280     if(PDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer();
02281     first=false;
02282   }
02283 #ifdef debug
02284   G4cout<<"G4QLowE::CXS: *DONE* p="<<p<<",Z="<<Z<<",N="<<N<<",PDG="<<PDG<<G4endl;
02285 #endif
02286   return CSmanager->GetCrossSection(true, p, Z, N, PDG);
02287 }

Generated on Mon May 27 17:49:36 2013 for Geant4 by  doxygen 1.4.7