G4QInelastic.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 //
00028 //      ---------------- G4QInelastic class -----------------
00029 //                 by Mikhail Kossov, December 2003.
00030 // G4QInelastic class of the CHIPS Simulation Branch in GEANT4
00031 // ---------------------------------------------------------------
00032 // ****************************************************************************************
00033 // This Header is a part of the CHIPS physics package (author: M. Kosov)
00034 // ****************************************************************************************
00035 // Short description: This is a universal class for the incoherent (inelastic)
00036 // nuclear interactions in the CHIPS model.
00037 // ---------------------------------------------------------------------------
00038 //#define debug
00039 //#define pdebug
00040 //#define pickupd
00041 //#define ldebug
00042 //#define ppdebug
00043 //#define qedebug
00044 
00045 #include "G4QInelastic.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4HadronicDeprecate.hh"
00049 
00050 
00051 // Initialization of static vectors
00052 std::vector<G4int> G4QInelastic::ElementZ;            // Z of the element(i) in theLastCalc
00053 std::vector<G4double> G4QInelastic::ElProbInMat;      // SumProbabilityElements in Material
00054 std::vector<std::vector<G4int>*> G4QInelastic::ElIsoN;// N of isotope(j) of Element(i)
00055 std::vector<std::vector<G4double>*>G4QInelastic::IsoProbInEl;//SumProbabIsotopes inElementI
00056 
00057 G4QInelastic::G4QInelastic(const G4String& processName): 
00058  G4VDiscreteProcess(processName, fHadronic)
00059 {
00060   G4HadronicDeprecate("G4QInelastic");
00061 
00062   EnMomConservation = G4LorentzVector(0.,0.,0.,0.);
00063   nOfNeutrons       = 0;
00064 #ifdef debug
00065   G4cout<<"G4QInelastic::Constructor is called"<<G4endl;
00066 #endif
00067   if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
00068   //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max)
00069   G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
00070   G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
00071   G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
00072   //@@ Initialize here the G4QuasmonString parameters
00073 }
00074 
00075 G4bool   G4QInelastic::manualFlag=false; // If false then standard parameters are used
00076 G4double G4QInelastic::Temperature=140.; // Critical Temperature (sensitive at High En)
00077 G4double G4QInelastic::SSin2Gluons=0.3;  // Supression of s-quarks (in respect to u&d)
00078 G4double G4QInelastic::EtaEtaprime=0.3;  // Supression of eta mesons (gg->qq/3g->qq)
00079 G4double G4QInelastic::freeNuc=.04;       // Percentage of free nucleons on the surface
00080 G4double G4QInelastic::freeDib=.08;      // Percentage of free diBaryons on the surface
00081 G4double G4QInelastic::clustProb=4.;     // Nuclear clusterization parameter
00082 G4double G4QInelastic::mediRatio=1.;     // medium/vacuum hadronization ratio
00083 //G4int    G4QInelastic::nPartCWorld=152;// The#of particles initialized in CHIPS World
00084 //G4int    G4QInelastic::nPartCWorld=122;// The#of particles initialized in CHIPS World
00085 G4int    G4QInelastic::nPartCWorld=85;   // The#of particles initialized in CHIPS World Red
00086 G4double G4QInelastic::SolidAngle=0.5;   // Part of Solid Angle to capture (@@A-dep.)
00087 G4bool   G4QInelastic::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
00088 G4double G4QInelastic::PiPrThresh=141.4; // Pion Production Threshold for gammas
00089 G4double G4QInelastic::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
00090 G4double G4QInelastic::DiNuclMass=1870.; // DoubleNucleon Mass for VirtualNormalization
00091 G4double G4QInelastic::photNucBias=1.;   // BiasingParameter for photo(e,mu,tau)Nuclear
00092 G4double G4QInelastic::weakNucBias=1.;   // BiasingParameter for ChargedCurrents(nu,mu) 
00093 
00094 void G4QInelastic::SetManual()   {manualFlag=true;}
00095 void G4QInelastic::SetStandard() {manualFlag=false;}
00096 
00097 // Fill the private parameters
00098 void G4QInelastic::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
00099                                      G4double fN, G4double fD, G4double cP, G4double mR,
00100                                      G4int nParCW, G4double solAn, G4bool efFlag,
00101                                      G4double piThresh, G4double mpisq, G4double dinum)
00102 {
00103   Temperature=temper;
00104   SSin2Gluons=ssin2g;
00105   EtaEtaprime=etaetap;
00106   freeNuc=fN;
00107   freeDib=fD;
00108   clustProb=cP;
00109   mediRatio=mR;
00110   nPartCWorld = nParCW;
00111   EnergyFlux=efFlag;
00112   SolidAngle=solAn;
00113   PiPrThresh=piThresh;
00114   M2ShiftVir=mpisq;
00115   DiNuclMass=dinum;
00116   G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
00117   G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
00118   G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
00119   G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
00120 }
00121 
00122 void G4QInelastic::SetPhotNucBias(G4double phnB) {photNucBias=phnB;}
00123 void G4QInelastic::SetWeakNucBias(G4double ccnB) {weakNucBias=ccnB;}
00124 
00125 // Destructor
00126 
00127 G4QInelastic::~G4QInelastic()
00128 {
00129   // The following is just a copy of what is done in PostStepDoIt every interaction !
00130   // The correction is if(IPIE), so just for(...;ip<IPIE;...) does not work ! @@
00131   G4int IPIE=IsoProbInEl.size();            // How many old elements?
00132   if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
00133   {
00134     std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
00135     SPI->clear();
00136     delete SPI;
00137     std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
00138     IsN->clear();
00139     delete IsN;
00140   }
00141   ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
00142   ElementZ.clear();                         // Clear the body vector for Z of Elements
00143   IsoProbInEl.clear();                      // Clear the body vector for SPI
00144   ElIsoN.clear();                           // Clear the body vector for N of Isotopes
00145 }
00146 
00147 
00148 G4LorentzVector G4QInelastic::GetEnegryMomentumConservation()
00149 {
00150   return EnMomConservation;
00151 }
00152 
00153 G4int G4QInelastic::GetNumberOfNeutronsInTarget()
00154 {
00155   return nOfNeutrons;
00156 }
00157 
00158 // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
00159 // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
00160 // ********** All CHIPS cross sections are calculated in the surface units ************
00161 G4double G4QInelastic::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
00162 {
00163 #ifdef debug
00164   G4cout<<"G4QInelastic::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
00165 #endif
00166   *Fc = NotForced;
00167 #ifdef debug
00168   G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDynPart"<<G4endl;
00169 #endif
00170   const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
00171 #ifdef debug
00172   G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDef"<<G4endl;
00173 #endif
00174   G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00175   if( !IsApplicable(*incidentParticleDefinition))
00176   {
00177     G4cout<<"-W-G4QInelastic::GetMeanFreePath called for not implemented particle"<<G4endl;
00178     return DBL_MAX;
00179   }
00180   // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
00181   G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
00182 #ifdef debug
00183   G4cout<<"G4QInelastic::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
00184 #endif
00185   const G4Material* material = aTrack.GetMaterial();        // Get the current material
00186   const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00187   const G4ElementVector* theElementVector = material->GetElementVector();
00188   G4int nE=material->GetNumberOfElements();
00189 #ifdef debug
00190   G4cout<<"G4QInelastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
00191 #endif
00192   //G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear *Growing point*
00193   G4VQCrossSection* CSmanager=0;
00194   G4VQCrossSection* CSmanager2=0;
00195   G4QNeutronCaptureRatio* capMan=0;
00196   G4int pPDG=0;
00197   G4int pZ = incidentParticleDefinition->GetAtomicNumber();
00198   G4int pA = incidentParticleDefinition->GetAtomicMass();
00199   if(incidentParticleDefinition == G4Neutron::Neutron())
00200   {
00201     CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
00202     capMan=G4QNeutronCaptureRatio::GetPointer();
00203 #ifdef debug
00204     G4cout<<"G4QInelastic::GetMeanFreePath: CSmanager is defined for the neutron"<<G4endl;
00205 #endif
00206     pPDG=2112;
00207   }
00208   else if(incidentParticleDefinition == G4Proton::Proton())
00209   {
00210     CSmanager=G4QProtonNuclearCrossSection::GetPointer();
00211     pPDG=2212;
00212   }
00213   else if(incidentParticleDefinition == G4PionMinus::PionMinus())
00214   {
00215     CSmanager=G4QPionMinusNuclearCrossSection::GetPointer();
00216     pPDG=-211;
00217   }
00218   else if(incidentParticleDefinition == G4PionPlus::PionPlus())
00219   {
00220     CSmanager=G4QPionPlusNuclearCrossSection::GetPointer();
00221     pPDG=211;
00222   }
00223   else if(incidentParticleDefinition == G4KaonMinus::KaonMinus())
00224   {
00225     CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer();
00226     pPDG=-321;
00227   }
00228   else if(incidentParticleDefinition == G4KaonPlus::KaonPlus())
00229   {
00230     CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer();
00231     pPDG=321;
00232   }
00233   else if(incidentParticleDefinition == G4KaonZeroLong::KaonZeroLong()   ||
00234           incidentParticleDefinition == G4KaonZeroShort::KaonZeroShort() ||
00235           incidentParticleDefinition == G4KaonZero::KaonZero()           ||
00236           incidentParticleDefinition == G4AntiKaonZero::AntiKaonZero()   )
00237   {
00238     CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer();
00239     if(G4UniformRand() > 0.5) pPDG= 311;
00240     else                      pPDG=-311;
00241   }
00242   else if(incidentParticleDefinition == G4Lambda::Lambda())
00243   {
00244     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00245     pPDG=3122;
00246   }
00247   else if(pZ > 0 && pA > 1) // Ions (not implemented yet, should not be used)
00248   {
00249     G4cout<<"-Warning-G4QInelastic::GetMeanFreePath: G4QInelastic called for ions"<<G4endl;
00250     CSmanager=G4QProtonNuclearCrossSection::GetPointer();
00251     pPDG=90000000+999*pZ+pA;
00252   }
00253   else if(incidentParticleDefinition == G4SigmaPlus::SigmaPlus())
00254   {
00255     CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer();
00256     pPDG=3222;
00257   }
00258   else if(incidentParticleDefinition == G4SigmaMinus::SigmaMinus())
00259   {
00260     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00261     pPDG=3112;
00262   }
00263   else if(incidentParticleDefinition == G4SigmaZero::SigmaZero())
00264   {
00265     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00266     pPDG=3212;
00267   }
00268   else if(incidentParticleDefinition == G4XiMinus::XiMinus())
00269   {
00270     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00271     pPDG=3312;
00272   }
00273   else if(incidentParticleDefinition == G4XiZero::XiZero())
00274   {
00275     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00276     pPDG=3322;
00277   }
00278   else if(incidentParticleDefinition == G4OmegaMinus::OmegaMinus())
00279   {
00280     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00281     pPDG=3334;
00282   }
00283   else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() ||
00284           incidentParticleDefinition == G4MuonMinus::MuonMinus())
00285   {
00286     CSmanager=G4QMuonNuclearCrossSection::GetPointer();
00287     //leptoNuc=true;
00288     pPDG=13;
00289   }
00290   else if(incidentParticleDefinition == G4AntiNeutron::AntiNeutron())
00291   {
00292     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00293     pPDG=-2112;
00294   }
00295   else if(incidentParticleDefinition == G4AntiProton::AntiProton())
00296   {
00297     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00298     pPDG=-2212;
00299   }
00300   else if(incidentParticleDefinition == G4AntiLambda::AntiLambda())
00301   {
00302     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00303     pPDG=-3122;
00304   }
00305   else if(incidentParticleDefinition == G4AntiSigmaPlus::AntiSigmaPlus())
00306   {
00307     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00308     pPDG=-3222;
00309   }
00310   else if(incidentParticleDefinition == G4AntiSigmaMinus::AntiSigmaMinus())
00311   {
00312     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00313     pPDG=-3112;
00314   }
00315   else if(incidentParticleDefinition == G4AntiSigmaZero::AntiSigmaZero())
00316   {
00317     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00318     pPDG=-3212;
00319   }
00320   else if(incidentParticleDefinition == G4AntiXiMinus::AntiXiMinus())
00321   {
00322     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00323     pPDG=-3312;
00324   }
00325   else if(incidentParticleDefinition == G4AntiXiZero::AntiXiZero())
00326   {
00327     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00328     pPDG=-3322;
00329   }
00330   else if(incidentParticleDefinition == G4AntiOmegaMinus::AntiOmegaMinus())
00331   {
00332     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00333     pPDG=-3334;
00334   }
00335   else if(incidentParticleDefinition == G4Gamma::Gamma())
00336   {
00337     CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
00338     pPDG=22;
00339   }
00340   else if(incidentParticleDefinition == G4Electron::Electron() ||
00341           incidentParticleDefinition == G4Positron::Positron())
00342   {
00343     CSmanager=G4QElectronNuclearCrossSection::GetPointer();
00344     //leptoNuc=true;
00345     pPDG=11;
00346   }
00347   else if(incidentParticleDefinition == G4TauPlus::TauPlus() ||
00348           incidentParticleDefinition == G4TauMinus::TauMinus())
00349   {
00350     CSmanager=G4QTauNuclearCrossSection::GetPointer();
00351     //leptoNuc=true;
00352     pPDG=15;
00353   }
00354   else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() )
00355   {
00356     CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
00357     CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
00358     //leptoNuc=true;
00359     pPDG=14;
00360   }
00361   else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() )
00362   {
00363     CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
00364     CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
00365     //leptoNuc=true;
00366     pPDG=-14;
00367   }
00368   else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() )
00369   {
00370     CSmanager=G4QNuENuclearCrossSection::GetPointer();
00371     CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
00372     //leptoNuc=true;
00373     pPDG=12;
00374   }
00375   else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() )
00376   {
00377     CSmanager=G4QANuENuclearCrossSection::GetPointer();
00378     CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
00379     //leptoNuc=true;
00380     pPDG=-12;
00381   }
00382   else
00383   {
00384     G4cout<<"-Warning-G4QInelastic::GetMeanFreePath:Particle "
00385           <<incidentParticleDefinition->GetPDGEncoding()<<" isn't supported"<<G4endl;
00386     return DBL_MAX;
00387   }
00388   G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
00389   G4double sigma=0.;                        // Sums over elements for the material
00390   G4int IPIE=IsoProbInEl.size();            // How many old elements?
00391   if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
00392   {
00393     std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
00394     SPI->clear();
00395     delete SPI;
00396     std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
00397     IsN->clear();
00398     delete IsN;
00399   }
00400   ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
00401   ElementZ.clear();                         // Clear the body vector for Z of Elements
00402   IsoProbInEl.clear();                      // Clear the body vector for SPI
00403   ElIsoN.clear();                           // Clear the body vector for N of Isotopes
00404   for(G4int i=0; i<nE; ++i)
00405   {
00406     G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
00407     G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
00408     ElementZ.push_back(Z);                  // Remember Z of the Element
00409     G4int isoSize=0;                        // The default for the isoVectorLength is 0
00410     G4int indEl=0;                          // Index of non-trivial element or 0(default)
00411     G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
00412     if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
00413 #ifdef debug
00414     G4cout<<"G4QInelastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
00415 #endif
00416     if(isoSize)                             // The Element has non-trivial abundance set
00417     {
00418       indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
00419       if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
00420       {
00421         std::vector<std::pair<G4int,G4double>*>* newAbund =
00422                                                new std::vector<std::pair<G4int,G4double>*>;
00423         G4double* abuVector=pElement->GetRelativeAbundanceVector();
00424         for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
00425         {
00426           G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
00427           if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QInelastic::GetMeanFreePath"
00428                                  <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00429           G4double abund=abuVector[j];
00430           std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00431 #ifdef debug
00432           G4cout<<"G4QInelastic::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
00433 #endif
00434           newAbund->push_back(pr);
00435         }
00436 #ifdef debug
00437         G4cout<<"G4QInelastic::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
00438 #endif
00439         indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
00440         for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
00441         delete newAbund; // Was "new" in the beginning of the name space
00442       }
00443     }
00444     std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
00445     std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
00446     IsoProbInEl.push_back(SPI);
00447     std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
00448     ElIsoN.push_back(IsN);
00449     G4int nIs=cs->size();                   // A#Of Isotopes in the Element
00450     G4double susi=0.;                       // sum of CS over isotopes
00451 #ifdef debug
00452     G4cout<<"G4QInelastic::GetMeanFreePath: Before Loop nIs="<<nIs<<G4endl;
00453 #endif
00454     if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
00455     {
00456       std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
00457       G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
00458       IsN->push_back(N);                    // Remember Min N for the Element
00459 #ifdef debug
00460       G4cout<<"G4QInel::GetMeanFrP: Before CS, P="<<Momentum<<",Z="<<Z<<",N="<<N<<G4endl;
00461 #endif
00462       if(!pPDG) G4cout<<"-Warning-G4QInelastic::GetMeanFrP: (1) projectile PDG=0"<<G4endl;
00463       G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
00464       if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu
00465       if(capMan) CSI*=(1.-capMan->GetRatio(Momentum, Z, N));
00466 #ifdef debug
00467       G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
00468 #endif
00469       curIs->second = CSI;
00470       susi+=CSI;                            // Make a sum per isotopes
00471       SPI->push_back(susi);                 // Remember summed cross-section
00472     } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
00473     sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
00474     ElProbInMat.push_back(sigma);
00475   } // End of LOOP over Elements
00476 #ifdef debug
00477   G4cout<<"G4QInel::GetMeanFrPa:S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
00478 #endif
00479   // Check that cross section is not zero and return the mean free path
00480   if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma()         ||
00481                          incidentParticleDefinition == G4MuonPlus::MuonPlus()   ||
00482                          incidentParticleDefinition == G4MuonMinus::MuonMinus() ||
00483                          incidentParticleDefinition == G4Electron::Electron()   ||
00484                          incidentParticleDefinition == G4Positron::Positron()   ||  
00485                          incidentParticleDefinition == G4TauMinus::TauMinus()   ||
00486                          incidentParticleDefinition == G4TauPlus::TauPlus()       )
00487                                                                         sigma*=photNucBias;
00488   if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE()            ||
00489                          incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE()    ||
00490                          incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau()        ||
00491                          incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()||
00492                          incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu()          ||
00493                          incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu()   )
00494                                                                         sigma*=weakNucBias;
00495   if(sigma > 0.) return 1./sigma;                 // Mean path [distance] 
00496   return DBL_MAX;
00497 }
00498 
00499 
00500 G4bool G4QInelastic::IsApplicable(const G4ParticleDefinition& particle) 
00501 {
00502   G4int Z=particle.GetAtomicNumber();
00503   G4int A=particle.GetAtomicMass();
00504   if      (particle == *(       G4Neutron::Neutron()       )) return true; 
00505   else if (particle == *(        G4Proton::Proton()        )) return true;
00506   else if (particle == *(      G4MuonPlus::MuonPlus()      )) return true;
00507   else if (particle == *(     G4MuonMinus::MuonMinus()     )) return true;
00508   else if (particle == *(         G4Gamma::Gamma()         )) return true;
00509   else if (particle == *(      G4Electron::Electron()      )) return true;
00510   else if (particle == *(      G4Positron::Positron()      )) return true;
00511   else if (particle == *(     G4PionMinus::PionMinus()     )) return true;
00512   else if (particle == *(      G4PionPlus::PionPlus()      )) return true;
00513   else if (particle == *(      G4KaonPlus::KaonPlus()      )) return true;
00514   else if (particle == *(     G4KaonMinus::KaonMinus()     )) return true;
00515   else if (particle == *(  G4KaonZeroLong::KaonZeroLong()  )) return true;
00516   else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
00517   else if (particle == *(        G4Lambda::Lambda()        )) return true;
00518   else if (particle == *(     G4SigmaPlus::SigmaPlus()     )) return true;
00519   else if (particle == *(    G4SigmaMinus::SigmaMinus()    )) return true;
00520   else if (particle == *(     G4SigmaZero::SigmaZero()     )) return true;
00521   else if (particle == *(       G4XiMinus::XiMinus()       )) return true;
00522   else if (particle == *(        G4XiZero::XiZero()        )) return true;
00523   else if (particle == *(    G4OmegaMinus::OmegaMinus()    )) return true;
00524   else if (particle == *(G4GenericIon::GenericIon()) || (Z > 0 && A > 1)) return true;
00525   else if (particle == *(   G4AntiNeutron::AntiNeutron()   )) return true;
00526   else if (particle == *(    G4AntiProton::AntiProton()    )) return true;
00527   else if (particle == *(    G4AntiLambda::AntiLambda()    )) return true;
00528   else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true;
00529   else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true;
00530   else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true;
00531   else if (particle == *(   G4AntiXiMinus::AntiXiMinus()   )) return true;
00532   else if (particle == *(    G4AntiXiZero::AntiXiZero()    )) return true;
00533   else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true;
00534   else if (particle == *(       G4TauPlus::TauPlus()       )) return true;
00535   else if (particle == *(      G4TauMinus::TauMinus()      )) return true;
00536   else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true;
00537   else if (particle == *(     G4NeutrinoE::NeutrinoE()     )) return true;
00538   else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
00539   else if (particle == *(    G4NeutrinoMu::NeutrinoMu()    )) return true;
00540   //else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true;
00541   //else if (particle == *(    G4NeutrinoTau::NeutrinoTau()    )) return true;
00542 #ifdef debug
00543   G4cout<<"***G4QInelastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
00544 #endif
00545   return false;
00546 }
00547 
00548 G4VParticleChange* G4QInelastic::PostStepDoIt(const G4Track& track, const G4Step& step)
00549 {
00550   static const G4double third = 1./3.;
00551   static const G4double me=G4Electron::Electron()->GetPDGMass();   // electron mass
00552   static const G4double me2=me*me;                                 // squared electron mass
00553   static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
00554   static const G4double mu2=mu*mu;                                 // squared muon mass
00555   static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass();   // tau mass
00556   static const G4double mt2=mt*mt;                                 // squared tau mass
00557   //static const G4double dpi=M_PI+M_PI;   // 2*pi (for Phi distr.) ***changed to twopi***
00558   static const G4double mNeut= G4QPDGCode(2112).GetMass();
00559   static const G4double mNeut2= mNeut*mNeut;
00560   static const G4double mProt= G4QPDGCode(2212).GetMass();
00561   static const G4double mProt2= mProt*mProt;
00562   static const G4double dM=mProt+mNeut;                            // doubled nucleon mass
00563   static const G4double hdM=dM/2.;                                 // M of the "nucleon"
00564   static const G4double hdM2=hdM*hdM;                              // M2 of the "nucleon"
00565   static const G4double mLamb= G4QPDGCode(3122).GetMass();         // Mass of Lambda/antiL
00566   static const G4double mPi0 = G4QPDGCode(111).GetMass();
00567   static const G4double mPi0s= mPi0*mPi0;
00568   static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
00569   //static const G4double mDeut2=mDeut*mDeut;                      // SquaredMassOfDeuteron
00570   static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
00571   static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
00572   static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
00573   static const G4double mPi  = G4QPDGCode(211).GetMass();
00574   static const G4double tmPi = mPi+mPi;     // Doubled mass of the charged pion
00575   static const G4double stmPi= tmPi*tmPi;   // Squared Doubled mass of the charged pion
00576   static const G4double mPPi = mPi+mProt;   // Delta threshold
00577   static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
00578   //static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2)
00579   // Static definitions for electrons (nu,e) -----------------------------------------
00580   static const G4double meN = mNeut+me;
00581   static const G4double meN2= meN*meN;
00582   static const G4double fmeN= 4*mNeut2*me2;
00583   static const G4double mesN= mNeut2+me2;
00584   static const G4double meP = mProt+me;
00585   static const G4double meP2= meP*meP;
00586   static const G4double fmeP= 4*mProt2*me2;
00587   static const G4double mesP= mProt2+me2;
00588   static const G4double medM= me2/dM;       // for x limit
00589   static const G4double meD = mPPi+me;      // Multiperipheral threshold
00590   static const G4double meD2= meD*meD;
00591   // Static definitions for muons (nu,mu) -----------------------------------------
00592   static const G4double muN = mNeut+mu;
00593   static const G4double muN2= muN*muN;
00594   static const G4double fmuN= 4*mNeut2*mu2;
00595   static const G4double musN= mNeut2+mu2;
00596   static const G4double muP = mProt+mu;
00597   static const G4double muP2= muP*muP;      // +
00598   static const G4double fmuP= 4*mProt2*mu2; // +
00599   static const G4double musP= mProt2+mu2;
00600   static const G4double mudM= mu2/dM;       // for x limit
00601   static const G4double muD = mPPi+mu;      // Multiperipheral threshold
00602   static const G4double muD2= muD*muD;
00603   // Static definitions for muons (nu,nu) -----------------------------------------
00604   //static const G4double nuN = mNeut;
00605   //static const G4double nuN2= mNeut2;
00606   //static const G4double fnuN= 0.;
00607   //static const G4double nusN= mNeut2;
00608   //static const G4double nuP = mProt;
00609   //static const G4double nuP2= mProt2;
00610   //static const G4double fnuP= 0.;
00611   //static const G4double nusP= mProt2;
00612   //static const G4double nudM= 0.;           // for x limit
00613   //static const G4double nuD = mPPi;         // Multiperipheral threshold
00614   //static const G4double nuD2= mPPi2;
00615   static const G4LorentzVector vacuum4M(0.,0.,0.,0.);
00616   //-------------------------------------------------------------------------------------
00617   static G4bool CWinit = true;              // CHIPS Warld needs to be initted
00618   if(CWinit)
00619   {
00620     CWinit=false;
00621     G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
00622   }
00623   //-------------------------------------------------------------------------------------
00624   const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00625   const G4ParticleDefinition* particle=projHadron->GetDefinition();
00626 #ifdef debug
00627   G4cout<<"G4QInelastic::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
00628 #endif
00629   G4ForceCondition cond=NotForced;
00630   GetMeanFreePath(track, 1., &cond);
00631 #ifdef debug
00632   G4cout<<"G4QInelastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
00633 #endif
00634   G4bool scat=false;                                  // No CHEX in proj scattering
00635   G4int  scatPDG=0;                                   // Must be filled if true (CHEX)
00636   G4LorentzVector proj4M=projHadron->Get4Momentum();  // 4-momentum of the projectile (IU?)
00637   G4LorentzVector scat4M=proj4M;                      // Must be filled if true
00638   G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
00639   G4double Momentum=proj4M.rho();
00640   if(std::fabs(Momentum-momentum)>.001)
00641     G4cerr<<"*G4QInelastic::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
00642 #ifdef debug
00643   G4double mp=proj4M.m();
00644   G4cout<<"-->G4QInel::PostStDoIt:*called*, 4M="<<proj4M<<", P="<<Momentum<<"="<<momentum
00645         <<",m="<<mp<<G4endl;
00646 #endif
00647   if (!IsApplicable(*particle))  // Check applicability
00648   {
00649     G4cerr<<"G4QInelastic::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
00650           <<G4endl;
00651     return 0;
00652   }
00653   const G4Material* material = track.GetMaterial();      // Get the current material
00654   G4int Z=0;
00655   const G4ElementVector* theElementVector = material->GetElementVector();
00656   G4int nE=material->GetNumberOfElements();
00657 #ifdef debug
00658   G4cout<<"G4QInelastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00659 #endif
00660   G4int projPDG=0;                           // PDG Code prototype for the captured hadron
00661   G4int pZ=particle->GetAtomicNumber();
00662   G4int pA=particle->GetAtomicMass();
00663   if      (particle ==         G4Neutron::Neutron()        ) projPDG= 2112;
00664   else if (particle ==          G4Proton::Proton()         ) projPDG= 2212;
00665   else if (particle ==       G4PionMinus::PionMinus()      ) projPDG= -211;
00666   else if (particle ==        G4PionPlus::PionPlus()       ) projPDG=  211;
00667   else if (particle ==        G4KaonPlus::KaonPlus()       ) projPDG=  321;
00668   else if (particle ==       G4KaonMinus::KaonMinus()      ) projPDG= -321;
00669   else if (particle ==    G4KaonZeroLong::KaonZeroLong()  ||
00670            particle ==   G4KaonZeroShort::KaonZeroShort() ||
00671            particle ==        G4KaonZero::KaonZero()      ||
00672            particle ==    G4AntiKaonZero::AntiKaonZero()   )
00673   {
00674     if(G4UniformRand() > 0.5)                                projPDG=  311;
00675     else                                                     projPDG= -311;
00676   }
00677   else if (                      pZ > 0 && pA > 1          ) projPDG = 90000000+999*pZ+pA;
00678   else if (particle ==          G4Lambda::Lambda()         ) projPDG= 3122;
00679   else if (particle ==       G4SigmaPlus::SigmaPlus()      ) projPDG= 3222;
00680   else if (particle ==      G4SigmaMinus::SigmaMinus()     ) projPDG= 3112;
00681   else if (particle ==       G4SigmaZero::SigmaZero()      ) projPDG= 3212;
00682   else if (particle ==         G4XiMinus::XiMinus()        ) projPDG= 3312;
00683   else if (particle ==          G4XiZero::XiZero()         ) projPDG= 3322;
00684   else if (particle ==      G4OmegaMinus::OmegaMinus()     ) projPDG= 3334;
00685   else if (particle ==     G4AntiNeutron::AntiNeutron()    ) projPDG=-2112;
00686   else if (particle ==      G4AntiProton::AntiProton()     ) projPDG=-2212;
00687   else if (particle ==        G4MuonPlus::MuonPlus()       ) projPDG=  -13;
00688   else if (particle ==       G4MuonMinus::MuonMinus()      ) projPDG=   13;
00689   else if (particle ==           G4Gamma::Gamma()          ) projPDG=   22;
00690   else if (particle ==        G4Electron::Electron()       ) projPDG=   11;
00691   else if (particle ==        G4Positron::Positron()       ) projPDG=  -11;
00692   else if (particle ==      G4NeutrinoMu::NeutrinoMu()     ) projPDG=   14;
00693   else if (particle ==  G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG=  -14;
00694   else if (particle ==      G4AntiLambda::AntiLambda()     ) projPDG=-3122;
00695   else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) projPDG=-3222;
00696   else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
00697   else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) projPDG=-3212;
00698   else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) projPDG=-3312;
00699   else if (particle ==      G4AntiXiZero::AntiXiZero()     ) projPDG=-3322;
00700   else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
00701   else if (particle ==       G4NeutrinoE::NeutrinoE()      ) projPDG=   12;
00702   else if (particle ==   G4AntiNeutrinoE::AntiNeutrinoE()  ) projPDG=  -12;
00703   else if (particle ==         G4TauPlus::TauPlus()        ) projPDG=  -15;
00704   else if (particle ==        G4TauMinus::TauMinus()       ) projPDG=   15;
00705   else if (particle ==     G4NeutrinoTau::NeutrinoTau()    ) projPDG=   16;
00706   else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG=  -16;
00707   G4int aProjPDG=std::abs(projPDG);
00708 #ifdef debug
00709   G4int prPDG=particle->GetPDGEncoding();
00710   G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
00711 #endif
00712   if(!projPDG)
00713   {
00714     G4cerr<<"---Warning---G4QInelastic::PostStepDoIt:Undefined interacting hadron"<<G4endl;
00715     return 0;
00716   }
00717   G4int EPIM=ElProbInMat.size();
00718 #ifdef debug
00719   G4cout<<"G4QInel::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00720 #endif
00721   G4int i=0;
00722   if(EPIM>1)
00723   {
00724     G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00725     for(i=0; i<nE; ++i)
00726     {
00727 #ifdef debug
00728       G4cout<<"G4QInelastic::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
00729 #endif
00730       if (rnd<ElProbInMat[i]) break;
00731     }
00732     if(i>=nE) i=nE-1;                        // Top limit for the Element
00733   }
00734   G4Element* pElement=(*theElementVector)[i];
00735   Z=static_cast<G4int>(pElement->GetZ());
00736 #ifdef debug
00737     G4cout<<"G4QInelastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
00738 #endif
00739   if(Z<=0)
00740   {
00741     G4cerr<<"---Warning---G4QInelastic::PostStepDoIt: Element with Z="<<Z<<G4endl;
00742     if(Z<0) return 0;
00743   }
00744   std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
00745   std::vector<G4int>* IsN = ElIsoN[i];     // Vector of "#of neutrons" in the isotope El[i]
00746   G4int nofIsot=SPI->size();               // #of isotopes in the element i
00747 #ifdef debug
00748   G4cout<<"G4QInel::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
00749 #endif
00750   G4int j=0;
00751   if(nofIsot>1)
00752   {
00753     G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
00754     for(j=0; j<nofIsot; ++j)
00755     {
00756 #ifdef debug
00757       G4cout<<"G4QInelastic::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
00758 #endif
00759       if(rndI < (*SPI)[j]) break;
00760     }
00761     if(j>=nofIsot) j=nofIsot-1;            // Top limit for the isotope
00762   }
00763   G4int N =(*IsN)[j]; ;                    // Randomized number of neutrons
00764 #ifdef debug
00765   G4cout<<"G4QInelastic::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl;
00766 #endif
00767   G4double kinEnergy= projHadron->GetKineticEnergy();
00768   G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00769   if(projPDG==22 && Z==1 && !N && Momentum<150.)
00770   {
00771 #ifdef debug
00772     G4cerr<<"---LowPHOTON---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl;
00773 #endif
00774     //Do Nothing Action insead of the reaction
00775     aParticleChange.ProposeEnergy(kinEnergy);
00776     aParticleChange.ProposeLocalEnergyDeposit(0.);
00777     aParticleChange.ProposeMomentumDirection(dir);
00778     aParticleChange.ProposeTrackStatus(fAlive);
00779     return G4VDiscreteProcess::PostStepDoIt(track,step);
00780   }
00781   if(N<0)
00782   {
00783     G4cerr<<"-Warning-G4QInelastic::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
00784     return 0;
00785   }
00786   nOfNeutrons=N;                           // Remember it for the energy-momentum check
00787   //G4double am=Z+N;
00788   G4int am=Z+N;
00789   //G4double dd=0.025;
00790   //G4double sr=std::sqrt(am);
00791   //G4double dsr=0.01*(sr+sr);
00792   //if(dsr<dd)dsr=dd;
00793   //G4double medRA=mediRatio*G4QThd(am);
00794   G4double medRA=mediRatio;
00795   if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,medRA); // ManualP
00796   //else if(projPDG==-2212)G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,9.);//aP ClustPars
00797   //else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.); //Pi-ClustPars
00798   //else if(projPDG ==-2212 || projPDG ==-2112)
00799   //                       G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,1.);//aP ClustPars
00800   //else if(projPDG ==-211 ||projPDG == 211 )
00801   //                       G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,1.); //Pi-ClustPars
00802 #ifdef debug
00803   G4cout<<"G4QInelastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
00804 #endif
00805   aParticleChange.Initialize(track);
00806   G4double weight = track.GetWeight();
00807 #ifdef debug
00808   G4cout<<"G4QInelastic::PostStepDoIt: weight="<<weight<<G4endl;
00809 #endif
00810   if(photNucBias!=1.)      weight/=photNucBias;
00811   else if(weakNucBias!=1.) weight/=weakNucBias;
00812   G4double localtime = track.GetGlobalTime();
00813 #ifdef debug
00814   G4cout<<"G4QInelastic::PostStepDoIt: localtime="<<localtime<<G4endl;
00815 #endif
00816   G4ThreeVector position = track.GetPosition();
00817   G4TouchableHandle trTouchable = track.GetTouchableHandle();
00818 #ifdef debug
00819   G4cout<<"G4QInelastic::PostStepDoIt: position="<<position<<G4endl;
00820 #endif
00821   //
00822   G4int targPDG=90000000+Z*1000+N;            // PDG Code of the target nucleus
00823   G4QPDGCode targQPDG(targPDG);
00824   G4double tgM=targQPDG.GetMass();            // Target mass
00825   G4double tM=tgM;                            // Target mass (copy to be changed)
00826   G4QHadronVector* output=new G4QHadronVector;// Prototype of Fragm Output G4QHadronVector
00827   G4double absMom = 0.;                       // Prototype of absorbed by nucleus Moment
00828   G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVector
00829   G4LorentzVector lead4M(0.,0.,0.,0.);        // Prototype of LeadingQ 4-momentum
00830 #ifdef debug
00831   G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<aProjPDG<<", targPDG="<<targPDG<<G4endl;
00832 #endif
00833   //  
00834   // Leptons with photonuclear
00835   // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
00836   //
00837   if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15)
00838   {
00839 #ifdef debug
00840     G4cout<<"G4QInelastic::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
00841 #endif
00842     G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
00843     G4double ml=me;
00844     G4double ml2=me2;
00845     if(aProjPDG== 13)
00846     {
00847       CSmanager=G4QMuonNuclearCrossSection::GetPointer();
00848       ml=mu;
00849       ml2=mu2;
00850     }
00851     if(aProjPDG== 15)
00852     {
00853       CSmanager=G4QTauNuclearCrossSection::GetPointer();
00854       ml=mt;
00855       ml2=mt2;
00856     }
00857     // @@ Probably this is not necessary any more (?)
00858     if(!aProjPDG) G4cout<<"-Warning-G4QInelast::PostStepDoIt:(2) projectile PDG=0"<<G4endl;
00859     G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS
00860     // @@ check a possibility to separate p, n, or alpha (!)
00861     if(Z==1 && !N && Momentum<150.) xSec=0.;
00862 #ifdef debug
00863     G4cout<<"-Forse-G4QInel::PStDoIt:P="<<Momentum<<",PDG="<<projPDG<<",xS="<<xSec<<G4endl;
00864 #endif
00865     if(xSec <= 0.) // The cross-section is 0 -> Do Nothing
00866     {
00867 #ifdef debug
00868       G4cerr<<"---OUT---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl;
00869 #endif
00870       //Do Nothing Action insead of the reaction
00871       aParticleChange.ProposeEnergy(kinEnergy);
00872       aParticleChange.ProposeLocalEnergyDeposit(0.);
00873       aParticleChange.ProposeMomentumDirection(dir);
00874       aParticleChange.ProposeTrackStatus(fAlive);
00875       delete output;
00876       return G4VDiscreteProcess::PostStepDoIt(track,step);
00877     }
00878     G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
00879 #ifdef debug
00880     G4cout<<"G4QInel::PStDoIt:kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
00881 #endif
00882     if( kinEnergy < photonEnergy || photonEnergy < 0.)
00883     {
00884       //Do Nothing Action insead of the reaction
00885       G4cerr<<"--G4QInelastic::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
00886       aParticleChange.ProposeEnergy(kinEnergy);
00887       aParticleChange.ProposeLocalEnergyDeposit(0.);
00888       aParticleChange.ProposeMomentumDirection(dir);
00889       aParticleChange.ProposeTrackStatus(fAlive);
00890       delete output;
00891       return G4VDiscreteProcess::PostStepDoIt(track,step);
00892     }
00893     G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
00894     G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
00895     if(tM<999.) W-=mPi0+mPi0s/dM;       // Pion production threshold for a nucleon target
00896     if(W<0.) 
00897     {
00898       //Do Nothing Action insead of the reaction
00899 #ifdef debug
00900       G4cout<<"--G4QInelastic::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
00901 #endif
00902       aParticleChange.ProposeEnergy(kinEnergy);
00903       aParticleChange.ProposeLocalEnergyDeposit(0.);
00904       aParticleChange.ProposeMomentumDirection(dir);
00905       aParticleChange.ProposeTrackStatus(fAlive);
00906       delete output;
00907       return G4VDiscreteProcess::PostStepDoIt(track,step);
00908     }
00909     // Update G4VParticleChange for the scattered muon
00910     G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
00911     G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS
00912     G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22);       // Real XS
00913     G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
00914     if(sigNu*G4UniformRand()>sigK*rndFraction) 
00915     {
00916       //NothingToDo Action insead of the reaction
00917 #ifdef debug
00918       G4cout<<"-DoNoth-G4QInelastic::PostStepDoIt: probab. correction - DoNothing"<<G4endl;
00919 #endif
00920       aParticleChange.ProposeEnergy(kinEnergy);
00921       aParticleChange.ProposeLocalEnergyDeposit(0.);
00922       aParticleChange.ProposeMomentumDirection(dir);
00923       aParticleChange.ProposeTrackStatus(fAlive);
00924       delete output;
00925       return G4VDiscreteProcess::PostStepDoIt(track,step);
00926     }
00927     G4double iniE=kinEnergy+ml;          // Initial total energy of the lepton
00928     G4double finE=iniE-photonEnergy;     // Final total energy of the lepton
00929 #ifdef pdebug
00930     G4cout<<"G4QInelastic::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl;
00931 #endif
00932     aParticleChange.ProposeEnergy(finE-ml);
00933     if(finE<=ml)                         // Secondary lepton (e/mu/tau) at rest disappears
00934     {
00935       aParticleChange.ProposeEnergy(0.);
00936       if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill);
00937       else aParticleChange.ProposeTrackStatus(fStopButAlive);
00938       aParticleChange.ProposeMomentumDirection(dir);
00939     }
00940     else aParticleChange.ProposeTrackStatus(fAlive);
00941     G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron
00942     G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron
00943     G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
00944 #ifdef pdebug
00945     G4cout<<"G4QI::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
00946 #endif
00947     if(cost>1.) cost=1.;                 // To avoid the accuracy of calculation problem
00948     if(cost<-1.) cost=-1.;               // To avoid the accuracy of calculation problem
00949     //
00950     // Scatter the lepton ( @@ make the same thing for real photons)
00951     // At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart
00952     G4double absEn = G4QThd(am)*GeV;     // @@(b) Mean Energy Absorbed by a Nucleus
00953     //G4double absEn = std::pow(am,third)*GeV;  // @@(b) Mean Energy Absorbed by a Nucleus
00954     if(am>1 && absEn < photonEnergy)     // --> the absorption of energy can happen
00955     //if(absEn < photonEnergy)             // --> the absorption of energy can happen
00956     {
00957       G4double abtEn = absEn+hdM;        // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
00958       G4double abEn2 = abtEn*abtEn;      // Squared absorbed Energy + MN
00959       G4double abMo2 = abEn2-hdM2;       // Squared absorbed Momentum of compound system
00960       G4double phEn2 = photonEnergy*photonEnergy;
00961       G4double phMo2 = phEn2+photonQ2;   // Squared momentum of primary virtual photon
00962       G4double phMo  = std::sqrt(phMo2); // Momentum of the primary virtual photon
00963       absMom         = std::sqrt(abMo2); // Absorbed Momentum
00964       if(absMom < phMo)                  // --> the absorption of momentum can happen
00965       {
00966         G4double dEn = photonEnergy - absEn; // Leading energy
00967         G4double dMo = phMo - absMom;    // Leading momentum
00968         G4double sF  = dEn*dEn - dMo*dMo;// s of leading particle
00969 #ifdef ppdebug
00970         G4cout<<"-PhotoAbsorption-G4QIn::PStDoIt: sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
00971 #endif
00972         if(sF > stmPi)                   // --> Leading fragmentation is possible
00973         {
00974           photonEnergy = absEn;          // New value of the photon energy
00975           photonQ2=abMo2-absEn*absEn;    // New value of the photon Q2
00976           absEn = dEn;                   // Put energy of leading particle to absEn (!)
00977         }
00978         else absMom=0.;                  // Flag that nothing has happened
00979       }
00980       else absMom=0.;                    // Flag that nothing has happened
00981     }
00982     // ------------- End of ProjPart selection
00983     //
00984     // Scattering in respect to the derection of the incident muon is made impicitly:
00985     G4ThreeVector ort=dir.orthogonal();  // Not normed orthogonal vector (!) (to dir)
00986     G4ThreeVector ortx = ort.unit();     // First unit vector orthogonal to the direction
00987     G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
00988     G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
00989     G4double phi=twopi*G4UniformRand();  // phi of scattered electron
00990     G4double sinx=sint*std::sin(phi);    // x perpendicular component
00991     G4double siny=sint*std::cos(phi);    // y perpendicular component
00992     G4ThreeVector findir=cost*dir+sinx*ortx+siny*orty;
00993     aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
00994 #ifdef pdebug
00995     G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
00996           <<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
00997 #endif
00998     G4ThreeVector photon3M=iniP*dir-finP*findir;// 3D total momentum of photon
00999     if(absMom)                           // Photon must be reduced & LeadingSyst fragmented
01000     {
01001       G4double ptm=photon3M.mag();                    // 3M of the virtual photon
01002 #ifdef ppdebug
01003       G4cout<<"-Absorption-G4QInelastic::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
01004             <<iniP*dir<<", eFin3M="<<finP*findir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
01005 #endif
01006       G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q
01007       photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn)
01008       proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
01009 #ifdef ppdebug
01010       G4cout<<"-->G4QI::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
01011 #endif
01012       lead4M=proj4M;                        // Remember 4-mom for the total 4-momentum
01013       G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
01014       try                                                           //                    |
01015       {                                                             //                    |
01016         if(leadhs) delete leadhs;                                   //                    |
01017         G4QNucleus vac(90000000);                                   //                    |
01018         leadhs=pan->Fragment(vac,1);  // DELETED after it is copied to output vector      |
01019       }                                                             //                    |
01020       catch (G4QException& error)                                   //                    |
01021       {                                                             //                    |
01022         G4cerr<<"***G4QInelastic::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//|
01023         // G4Exception("G4QInelastic::PostStepDoIt:","72",FatalException,"QuasmonCrash");  //|
01024         G4Exception("G4QInelastic::PostStepDoIt()","HAD_CHPS_0072",
01025                     FatalException, "QuasmonCrash");
01026       }                                                             //                    |
01027       delete pan;                              // Delete the Nuclear Environment <----<---+
01028 #ifdef ppdebug
01029       G4cout<<"G4QInel::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
01030             <<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
01031 #endif
01032     }
01033     else
01034     {
01035       G4int qNH=0;
01036       if(leadhs) qNH=leadhs->size();
01037       if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
01038       if(leadhs) delete leadhs;
01039       leadhs=0;
01040     }
01041     projPDG=22;
01042     proj4M=G4LorentzVector(photon3M,photonEnergy);
01043 #ifdef debug
01044     G4cout<<"G4QInelastic::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
01045           <<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
01046 #endif
01047 
01048   //
01049   // neutrinoNuclear interactions (nu_e, nu_mu only)
01050   //
01051   }
01052   else if (aProjPDG == 12 || aProjPDG == 14)
01053   {
01054     kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
01055     G4double dKinE=kinEnergy+kinEnergy;  // doubled energy for s calculation
01056 #ifdef debug
01057     G4cout<<"G4QInelastic::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
01058 #endif
01059     dir = projHadron->GetMomentumDirection(); // unit vector
01060     G4double ml  = mu;
01061     G4double ml2 = mu2;
01062     //G4double mlN = muN;
01063     G4double mlN2= muN2;
01064     G4double fmlN= fmuN;
01065     G4double mlsN= musN;
01066     //G4double mlP = muP;
01067     G4double mlP2= muP2;
01068     G4double fmlP= fmuP;
01069     G4double mlsP= musP;
01070     G4double mldM= mudM;
01071     //G4double mlD = muD;
01072     G4double mlD2= muD2;
01073     if(aProjPDG==12)
01074     {
01075       ml  = me;
01076       ml2 = me2;
01077       //mlN = meN;
01078       mlN2= meN2;
01079       fmlN= fmeN;
01080       mlsN= mesN;
01081       //mlP = meP;
01082       mlP2= meP2;
01083       fmlP= fmeP;
01084       mlsP= mesP;
01085       mldM= medM;
01086       //mlD = meD;
01087       mlD2= meD2;
01088     }
01089     G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
01090     G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu)
01091     proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy);   // temporary
01092     G4bool nuanu=true;
01093     scatPDG=13;                          // Prototype = secondary scattered mu-
01094     if(projPDG==-14)
01095     {
01096       nuanu=false;                       // Anti-neutrino
01097       CSmanager=G4QANuMuNuclearCrossSection::GetPointer();  // (anu,mu+) CC @@ open
01098       CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
01099       scatPDG=-13;                       // secondary scattered mu+
01100     }
01101     else if(projPDG==12)
01102     {
01103       CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed)
01104       scatPDG=11;                        // secondary scattered e-
01105     }
01106     else if(projPDG==-12)
01107     {
01108       nuanu=false;                       // anti-neutrino
01109       CSmanager=G4QANuENuclearCrossSection::GetPointer();   // (anu,e+) CC @@ open
01110       CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
01111       scatPDG=-11;                       // secondary scattered e+
01112     }
01113     // @@ Probably this is not necessary any more
01114     if(!projPDG) G4cout<<"-Warning-G4QInelastic::PostStepDoIt:(3)projectile PDG=0"<<G4endl;
01115     G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS
01116     G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS
01117     G4double xSec=xSec1+xSec2;
01118     // @@ check a possibility to separate p, n, or alpha (!)
01119     if(xSec <= 0.) // The cross-section = 0 -> Do Nothing
01120     {
01121       G4cerr<<"G4QInelastic::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl;
01122       //Do Nothing Action insead of the reaction
01123       aParticleChange.ProposeEnergy(kinEnergy);
01124       aParticleChange.ProposeLocalEnergyDeposit(0.);
01125       aParticleChange.ProposeMomentumDirection(dir);
01126       aParticleChange.ProposeTrackStatus(fAlive);
01127       delete output;
01128       return G4VDiscreteProcess::PostStepDoIt(track,step);
01129     }
01130     G4bool secnu=false;
01131     if(xSec*G4UniformRand()>xSec1)               // recover neutrino/antineutrino
01132     {
01133       if(scatPDG>0) scatPDG++;
01134       else          scatPDG--;
01135       secnu=true;
01136     }
01137     scat=true;                                   // event with changed scattered projectile
01138     G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?)
01139     G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?)
01140     G4double totCS  = totCS1+totCS2;             // the last total cross section (isotope?)
01141     if(std::fabs(xSec-totCS*millibarn)/xSec>.0001)
01142           G4cout<<"-Warning-G4QInelastic::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl;
01143     G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1
01144     G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2
01145     G4double qelCS  = qelCS1+qelCS2;             // the last quasi-elastic cross section
01146     if(totCS - qelCS < 0.)                       // only at low energies
01147     {
01148       totCS  = qelCS;
01149       totCS1 = qelCS1;
01150       totCS2 = qelCS2;
01151     }
01152     // make different definitions for neutrino and antineutrino
01153     G4double mIN=mProt;                          // Just a prototype (for anu, Z=1, N=0)
01154     G4double mOT=mNeut;
01155     G4double OT=mlN2;
01156     //G4double mOT2=mNeut2;                        // Other formula: sd=s-mlsOT -Not used?-
01157     G4double mlOT=fmlN;
01158     G4double mlsOT=mlsN;
01159     if(secnu)
01160     {
01161       if(am*G4UniformRand()>Z)                   // Neutron target
01162       {
01163         targPDG-=1;                              // subtract neutron
01164         projPDG=2112;                            // neutron is going out
01165         mIN =mNeut;                     
01166         OT  =mNeut2;
01167         //mOT2=mNeut2;
01168         mlOT=0.;
01169         mlsOT=mNeut2;
01170       }
01171       else
01172       {
01173         targPDG-=1000;                           // subtract neutron
01174         projPDG=2212;                            // neutron is going out
01175         mOT  =mProt;
01176         OT   =mProt2;
01177         //mOT2 =mProt2;
01178         mlOT =0.;
01179         mlsOT=mProt2;
01180       }
01181       ml=0.;
01182       ml2=0.;
01183       mldM=0.;
01184       mlD2=mPPi2;
01185       G4QPDGCode temporary_targQPDG(targPDG); 
01186       targQPDG = temporary_targQPDG;
01187       G4double rM=targQPDG.GetMass();
01188       mIN=tM-rM;                                 // bounded in-mass of the neutron
01189       tM=rM;
01190     }
01191     else if(nuanu)
01192     {
01193       targPDG-=1;                                // Neutrino -> subtract neutron
01194       G4QPDGCode temporary_targQPDG(targPDG); 
01195       targQPDG = temporary_targQPDG;
01196       G4double rM=targQPDG.GetMass();
01197       mIN=tM-rM;                                 // bounded in-mass of the neutron
01198       tM=rM;
01199       mOT=mProt;
01200       OT=mlP2;
01201       //mOT2=mProt2;
01202       mlOT=fmlP;
01203       mlsOT=mlsP;
01204       projPDG=2212;                              // proton is going out
01205     }
01206     else
01207     {
01208       if(Z>1||N>0)                               // Calculate the splitted mass
01209       {
01210         targPDG-=1000;                           // Anti-Neutrino -> subtract proton
01211         G4QPDGCode temporary_targQPDG(targPDG); 
01212         targQPDG = temporary_targQPDG;
01213         G4double rM=targQPDG.GetMass();
01214         mIN=tM-rM;                               // bounded in-mass of the proton
01215         tM=rM;
01216       }
01217       else
01218       {
01219         targPDG=0;
01220         mIN=tM;
01221         tM=0.;
01222       }
01223       projPDG=2112;                              // neutron is going out
01224     }
01225     G4double s_value=mIN*(mIN+dKinE);            // s_value=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
01226 #ifdef debug
01227     G4cout<<"G4QInelastic::PostStDoIt: s="<<s_value<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
01228 #endif
01229     if(s_value<=OT)                                    // *** Do nothing solution ***
01230     {
01231       //Do NothingToDo Action insead of the reaction (@@ Can we make it common?)
01232       G4cout<<"G4QInelastic::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl;
01233       aParticleChange.ProposeEnergy(kinEnergy);
01234       aParticleChange.ProposeLocalEnergyDeposit(0.);
01235       aParticleChange.ProposeMomentumDirection(dir);
01236       aParticleChange.ProposeTrackStatus(fAlive);
01237       delete output;
01238       return G4VDiscreteProcess::PostStepDoIt(track,step);
01239     }
01240 #ifdef debug
01241     G4cout<<"G4QInelastic::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
01242 #endif
01243     aParticleChange.ProposeEnergy(0.);
01244     aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed
01245     // There is no way back from here !
01246     if ( ((secnu || !nuanu || N) && totCS*G4UniformRand() < qelCS) || s_value < mlD2 ) 
01247     {   // Quasi-Elastic interaction
01248       G4double Q2=0.;                           // Simulate transferred momentum, in MeV^2
01249       if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2();
01250       else      Q2=CSmanager->GetQEL_ExchangeQ2();
01251 #ifdef debug
01252       G4cout<<"G4QInelastic::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s_value<<",Q2="<<Q2<<G4endl;
01253 #endif
01254       //G4double ds=s_value+s_value;              // doubled s_value
01255       G4double sqs=std::sqrt(s_value);          // M_cm
01256       G4double dsqs=sqs+sqs;                    // 2*M_cm
01257       G4double p_init=(s_value-mIN*mIN)/dsqs;   // initial momentum in CMS (checked MK)
01258       G4double dpi=p_init+p_init;               // doubled initial momentum in CMS
01259       G4double sd=s_value-mlsOT;                // s_value-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept)
01260       G4double qo2=(sd*sd-mlOT)/dsqs;           // squared momentum of secondaries in CMS
01261       G4double qo=std::sqrt(qo2);               // momentum of secondaries in CMS
01262       G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK)
01263       G4LorentzVector t4M(0.,0.,0.,mIN);        // 4mom of the effective target
01264       G4LorentzVector c4M=t4M+proj4M;           // 4mom of the compound system
01265       t4M.setT(mOT);                            // now it is 4mom of the outgoing nucleon
01266       scat4M=G4LorentzVector(0.,0.,0.,ml);      // 4mom of the scattered muon
01267       if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost))
01268       {
01269         G4cerr<<"G4QIn::PStD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl;
01270         // throw G4QException("G4QInelastic::HadronizeQuasm: Can't dec QE nu,lept Compound");
01271         G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0000",
01272                     FatalException, "Hadronize quasmon: Can't dec QE nu,lept Compound");
01273       }
01274       proj4M=t4M;                               // 4mom of the new projectile nucleon
01275     }
01276     else                                        // ***** Non Quasi Elastic interaction
01277     {
01278       // Recover the target PDG Code if the quasi-elastic scattering did not happened
01279       if      ( (secnu && projPDG == 2212) || (!secnu && projPDG == 2112) ) targPDG+=1;
01280       else if ( (secnu && projPDG == 2112) || (!secnu && projPDG == 2212) ) targPDG+=1000;
01281       G4double Q2=0;                            // Simulate transferred momentum, in MeV^2
01282       if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2();
01283       else      Q2=CSmanager2->GetNQE_ExchangeQ2();
01284 #ifdef debug
01285       G4cout<<"G4QInel::PStDoIt: MultiPeriferal s="<<s_value<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
01286 #endif
01287       if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
01288       else      projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion
01289       //@@ Temporary made only for direct interaction and for N=3 (good for small Q2)
01290       //@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2...
01291       G4double r=G4UniformRand();
01292       G4double r1=0.5;                          // (1-x)
01293       if(r<0.5)      r1=std::sqrt(r+r)*(.5+.1579*(r-.5));
01294       else if(r>0.5) r1=1.-std::sqrt(2.-r-r)*(.5+.1579*(.5-r));
01295       G4double xn=1.-mldM/Momentum;             // Normalization of (1-x) [x>mldM/Mom]
01296       G4double x1=xn*r1;                        // (1-x)
01297       G4double x=1.-x1;                         // x=2k/M
01298       //G4double W2=(hdM2+Q2/x)*x1;               // W2 candidate
01299       G4double mx=hdM*x;                        // Part of the target to interact with
01300       G4double we=Q2/(mx+mx);                   // transfered energy
01301       if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg)
01302       G4double pot=kinEnergy-we;                // energy of the secondary lepton
01303       G4double mlQ2=ml2+Q2;
01304       G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta)
01305       if(std::fabs(cost)>1)
01306       {
01307 #ifdef debug
01308         G4cout<<"*G4QInelastic::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
01309               <<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
01310 #endif
01311         if(cost>1.) cost=1.;
01312         else        cost=-1.;
01313         pot=mlQ2/dKinE+dKinE*ml2/mlQ2;          // extreme output momentum
01314       }
01315       G4double lEn=std::sqrt(pot*pot+ml2);      // Lepton energy
01316       G4double lPl=pot*cost;                    // Lepton longitudinal momentum
01317       G4double lPt=pot*std::sqrt(1.-cost*cost); // Lepton transverse momentum
01318       std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi
01319       G4double lPx=lPt*d2d.first;
01320       G4double lPy=lPt*d2d.second;
01321       G4ThreeVector vdir=proj4M.vect();         // 3D momentum of the projectile
01322       G4ThreeVector vz= vdir.unit();            // Ort in the direction of the projectile
01323       G4ThreeVector vv= vz.orthogonal();        // Not normed orthogonal vector (!)
01324       G4ThreeVector vx= vv.unit();              // First ort orthogonal to the direction
01325       G4ThreeVector vy= vz.cross(vx);           // Second ort orthoganal to the direction
01326       G4ThreeVector lP= lPl*vz+lPx*vx+lPy*vy;   // 3D momentum of the scattered lepton
01327       scat4M=G4LorentzVector(lP,lEn);           // 4mom of the scattered lepton
01328       proj4M-=scat4M;                           // 4mom of the W/Z (effective pion/gamma)
01329 #ifdef debug
01330       G4cout<<"G4QInelastic::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
01331 #endif
01332       // Check that the en/mom transfer is possible, if not -> elastic
01333       G4int fintPDG=targPDG;                    // Prototype for the compound nucleus
01334       if(!secnu)
01335       {
01336         if(projPDG<0) fintPDG-= 999;
01337         else          fintPDG+= 999;
01338       }
01339       G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV)
01340       G4double fM2=fM*fM;
01341       G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM);
01342       G4LorentzVector c4M=tg4M+proj4M;
01343 #ifdef debug
01344       G4cout<<"G4QInelastic::PSDI:fM2="<<fM2<<"<? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl;
01345 #endif
01346       if(fM2>=c4M.m2())                         // Elastic scattering should be done
01347       {
01348         G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum
01349         s_value=tot4M.m2();
01350         G4double fs=s_value-fM2-ml2;
01351         G4double fMl=fM2*ml2;
01352         G4double hQ2max=(fs*fs/2-fMl-fMl)/s_value; // Maximum possible Q2/2
01353         cost=1.-Q2/hQ2max;                // cos(theta) in CMS (use MultProd Q2)
01354 #ifdef debug
01355         G4cout<<"G4QI::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl;
01356 #endif
01357         G4double acost=std::fabs(cost);
01358         if(acost>1.)
01359         {
01360           if(acost>1.001) G4cout<<"-Warning-G4QInelastic::PostStDoIt: cost="<<cost<<G4endl;
01361           if     (cost> 1.) cost= 1.;
01362           else if(cost<-1.) cost=-1.;
01363         }
01364         G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus
01365         scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton
01366         G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01);
01367         if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
01368         {
01369           G4cerr<<"G4QI::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl;
01370           //G4Exception("G4QInelastic::PostStepDoIt:","027",FatalException,"ElasticDecay");
01371         }
01372 #ifdef debug
01373         G4cout<<"G4QIn::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
01374 #endif
01375         // ----------------------------------------------------
01376         G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries
01377         // Fill scattered lepton
01378         if     (scatPDG==-11) theDefinition = G4Positron::Positron();
01379         else if(scatPDG== 11) theDefinition = G4Electron::Electron();
01380         else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus();
01381         else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus();
01382         //else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus();
01383         //else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus();
01384         if     (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE();
01385         else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE();
01386         else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu();
01387         else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu();
01388         //else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau();
01389         //else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau();
01390         else  G4cout<<"-Warning-G4QInelastic::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl;
01391         G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M);
01392         G4Track* scatLep = new G4Track(theScL, localtime, position ); //    scattered
01393         scatLep->SetWeight(weight);                                   //    weighted
01394         scatLep->SetTouchableHandle(trTouchable);                     //    residual
01395         aParticleChange.AddSecondary(scatLep);                        //    lepton
01396         // Fill residual nucleus
01397         if     (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron
01398         else if(fintPDG==90001000) theDefinition = G4Proton::Proton();   // proton
01399         else                                                             // ion
01400         {
01401           G4int fm=static_cast<G4int>(fintPDG/1000000);               // Strange part
01402           G4int ZN=fintPDG-1000000*fm;
01403           G4int rZ=static_cast<G4int>(ZN/1000);
01404           G4int rA=ZN-999*rZ;
01405           theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
01406         }
01407         G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M);
01408         G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
01409         scatReN->SetWeight(weight);                                   //    weighted
01410         scatReN->SetTouchableHandle(trTouchable);                     //    residual
01411         aParticleChange.AddSecondary(scatReN);                        //    nucleus
01412         delete output;
01413         return G4VDiscreteProcess::PostStepDoIt(track, step);
01414       }
01415     } // End of non-quasi-elastic interaction
01416   //
01417   // quasi-elastic (& pickup process) for p+A(Z,N)
01418   //
01419   }
01420   //else if ((projPDG == 2212 || projPDG == 2112) && Z > 0 && N > 0) // Never (in Fragm!)
01421   else if(2>3) // Possibility is closed as quasi-free scattering is made in fragmentation
01422   {
01423     if(momentum<500. && projPDG == 2112) // @@ It's reasonable to add proton capture too !
01424     {
01425       G4LorentzVector tot4M=G4LorentzVector(0.,0.,0.,tgM)+proj4M;
01426       G4double totM2=tot4M.m2();
01427       G4int tZ=Z;
01428       G4int tN=N;
01429       if(projPDG == 2212) tZ++;
01430       else                tN++; // @@ Now a neutron is the only alternative
01431       if(momentum<100.)         // Try the production threshold (@@ Why 5 MeV?)
01432       {
01433         G4QPDGCode FakeP(2112);
01434         //G4int m1p=tZ-1;
01435         G4int m2n=tN-2;
01436         // @@ For the secondary proton the Coulomb barrier should be taken into account
01437         //G4double mRP= FakeP.GetNuclMass(m1p,tN,0);    // Residual to the secondary proton
01438         G4double mR2N= FakeP.GetNuclMass(tZ,m2n,0);   // Residual to two secondary neutrons
01439         //G4double mRAl= FakeP.GetNuclMass(tZ-2,m2n,0); // Residual to the secondary alpha
01440         //G4double mR2N= FakeP.GetNuclMass(m1p,tN-1,0); // Residual to the p+n secondaries
01441         G4double minM=mR2N+mNeut+mNeut;
01442         // Compare with other possibilities (sciped for acceleration)
01443         if(totM2 < minM*minM) //  DoNothing (under reasonable threshold)
01444         {
01445           aParticleChange.ProposeEnergy(kinEnergy);
01446           aParticleChange.ProposeLocalEnergyDeposit(0.);
01447           aParticleChange.ProposeMomentumDirection(dir);
01448           aParticleChange.ProposeTrackStatus(fAlive);
01449           return G4VDiscreteProcess::PostStepDoIt(track,step);
01450         }
01451       }
01452     }
01453     // Now the quasi-free and PickUp processes should be done:
01454     G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
01455     std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, projPDG, Z, N);
01456     G4double qepart=fief.first*fief.second; // @@ charge exchange is not included @@
01457     // Keep QE part for the quasi-free nucleons
01458     const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
01459     G4double clProb[lCl]={.65,.85,.95};// N/P,D,t/He3,He4, integroProbab for .65,.2,.1,.05
01460     if(qepart>0.45) qepart=.45; // Quasielastic part is too large - shrink
01461     //else qepart/=clProb[0];     // Add corresponding number of 2N, 3N, & 4N clusters
01462     qepart=qepart/clProb[0]-qepart;// Add QE for 2N, 3N, & 4N clusters (N is made in G4QF)
01463     G4double pickup=1.-qepart;  // Estimate the rest of the cross-section
01464     G4double thresh=100.;
01465     //if(momentum >thresh) pickup*=50./momentum/std::pow(G4double(Z+N),third);// 50 is Par
01466     if(momentum > thresh) pickup*=50./momentum/G4QThd(Z+N);// 50 is Par
01467     // pickup = 0.;               // To exclude the pickup process
01468     if (N) pickup+=qepart;
01469     else   pickup =qepart;
01470     G4double rnd=G4UniformRand();
01471 #ifdef ldebug
01472     G4cout<<"-->G4QInelastic::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<")]="<<qepart
01473           <<", pickup="<<pickup<<G4endl;
01474 #endif
01475     if(rnd<pickup) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim,CoulBar,PauliBl
01476     {
01477       G4double dmom=91.;  // Fermi momentum (proto default for a deuteron)
01478       G4double fmm=287.;  // @@ Can be A-dependent
01479       if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
01480       // Values to be redefined for quasi-elastic
01481       G4LorentzVector r4M(0.,0.,0.,0.);      // Prototype of 4mom of the residual nucleus
01482       G4LorentzVector n4M(0.,0.,0.,0.);      // Prototype of 4mom of the quasi-cluster
01483       G4int nPDG=90000001;                   // Prototype for quasiCluster mass calculation
01484       G4int restPDG=targPDG;                 // Prototype should be reduced by quasiCluster
01485       G4int rA=Z+N-1;                        // Prototype for the residualNucl definition
01486       G4int rZ=Z;                            // residZ: OK for the quasi-free neutron
01487       G4int nA=1;                            // Prototype for the quasi-cluster definition
01488       G4int nZ=0;                            // nA=1,nZ=0: OK for the quasi-free neutron
01489       G4double qM=mNeut;                     // Free mass of the quasi-free cluster
01490       G4int A = Z + N;                       // Baryon number of the nucleus
01491       if(rnd<qepart)
01492       {
01493        // ===> First decay a nucleus in a nucleon and a residual (A-1) nucleus
01494        // Calculate clusterProbabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
01495        G4double base=1.;  // Base for randomization (can be reduced by totZ & totN)
01496        G4int max=lCl;     // Number of boundaries (can be reduced by totZ & totN)
01497        // Take into account that at least one nucleon must be left !
01498        if(Z<2 || N<2 || A < 6) base = clProb[--max]; // Alpha cluster is impossible
01499        if((Z > 1 && N < 2) || (Z < 2 && N > 1)) base=(clProb[max]+clProb[max-1])/2;// t/He3
01500        if ( (Z < 2 && N < 2) || A < 5) base=clProb[--max];// He3&t clusters are impossible
01501        if(A<3)           base=clProb[--max]; // Deuteron cluster is impossible
01502        //G4int cln=0;                        // Cluster#0 (Default for the selectedNucleon)
01503        G4int cln=1;                          // Cluster#1 (Default for the selectedDeutron)
01504        //if(max)                             // Not only nucleons are possible
01505        if(max>1)                             // Not only deuterons are possible
01506        {
01507          base-=clProb[0];                   // Exclude scattering on QF Nucleon
01508          G4double ran=+clProb[0]+base*G4UniformRand(); // Base can be reduced
01509          G4int ic=1;                        // Start from the smallest cluster boundary
01510          if(max>1) while(ic<max) if(ran>clProb[ic++]) cln=ic;
01511        }
01512        G4ParticleDefinition* theDefinition=0;// Prototype for qfNucleon
01513        G4bool cp1 = cln+2==A;               // A=ClusterBN+1 condition
01514        if(!cln || cp1)                      // Split in nucleon + (A-1) with Fermi momentum
01515        { 
01516         G4int nln=0;
01517         if(cln==2) nln=1;                   // @@ only for cp1: t/He3 choice from A=4
01518         // mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass)
01519         if ( ((!cln || cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) || 
01520              ((cln == 3 || cln == 1) && Z > N) )
01521         {
01522           nPDG=90001000;                          // Update quasi-free nucleon PDGCode to P
01523           nZ=1;                                   // Change charge of the quasiFree nucleon
01524           qM=mProt;                               // Update quasi-free nucleon mass
01525           rZ--;                                   // Reduce the residual Z
01526           restPDG-=1000;                          // Reduce the residual PDGCode
01527         }
01528         else restPDG--;
01529         G4LorentzVector t4M(0.,0.,0.,tM);         // 4m of the target nucleus to be decayed
01530         G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
01531         r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
01532         G4double rM2=rM*rM;
01533         G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-nucleon
01534         n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
01535 #ifdef qedebug
01536         G4cout<<"G4QInel::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
01537 #endif
01538         if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
01539         {
01540           //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+n="<<nM<<"="<<rM+nM<<G4endl;
01541           //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0002",
01542           //           FatalException,"Hadronize quasmon: Can't Dec totNuc->QENuc+ResNuc");
01543           G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl;
01544           aParticleChange.ProposeEnergy(kinEnergy);
01545           aParticleChange.ProposeLocalEnergyDeposit(0.);
01546           aParticleChange.ProposeMomentumDirection(dir);
01547           aParticleChange.ProposeTrackStatus(fAlive);
01548           return G4VDiscreteProcess::PostStepDoIt(track,step);
01549         }
01550 #ifdef qedebug
01551         G4cout<<"G4QIn::PStDoIt:QE-N, RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
01552 #endif
01553         if(cp1 && cln)                           // Quasi-cluster case: swap the output
01554         {
01555           qM=rM;                                 // Scattering will be made on a cluster
01556           nln=nPDG;
01557           nPDG=restPDG;
01558           restPDG=nln;
01559           t4M=n4M;
01560           n4M=r4M;
01561           r4M=t4M;
01562           nln=nZ;
01563           nZ=rZ;
01564           rZ=nln;
01565           nln=nA;
01566           nA=rA;
01567           rA=nln;
01568         }
01569        }
01570        else // Split the cluster (with or w/o "Fermi motion" and "Fermi decay")
01571        {
01572         if(cln==1)
01573         {
01574           nPDG=90001001;                  // Deuteron
01575           qM=mDeut;
01576           nA=2;
01577           nZ=1;
01578           restPDG-=1001;
01579           // Recalculate dmom
01580           G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
01581                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01582           dmom=sumv.mag();
01583         }
01584         else if(cln==2)
01585         {
01586           nA=3;
01587           if(G4UniformRand()*(A-2)>(N-1)) // He3
01588           {
01589             nPDG=90002001;
01590             qM=mHel3;
01591             nZ=2;
01592             restPDG-=2001;
01593           }
01594           else                            // tritium
01595           {
01596             nPDG=90001002;
01597             qM=mTrit;
01598             nZ=1;
01599             restPDG-=1002;
01600           }
01601           // Recalculate dmom
01602           G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
01603                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
01604                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01605           dmom=sumv.mag();
01606         }
01607         else
01608         {
01609           nPDG=90002002;                  // Alpha
01610           qM=mAlph;
01611           nA=4;
01612           nZ=2;
01613           restPDG-=2002;
01614           G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
01615                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
01616                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
01617                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01618           dmom=sumv.mag();
01619         }
01620         rA=A-nA;
01621         rZ=Z-nZ;
01622         // This is a simple case of cluster at rest
01623         //G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
01624         //r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
01625         //n4M=G4LorentzVector(0.,0.,0.,tM-rM);      // 4mom of the quasi-free cluster
01626         // --- End of the "simple case of cluster at rest"
01627         // Make a fake quasi-Fermi distribution for clusters (clusters are not at rest) 
01628         G4LorentzVector t4M(0.,0.,0.,tM);         // 4m of the target nucleus to be decayed
01629         G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
01630         r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
01631         G4double rM2=rM*rM;
01632         G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-cluster
01633         n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
01634 #ifdef qedebug
01635         G4cout<<"G4QInel::PStDoIt:QEC, p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
01636 #endif
01637         if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
01638         {
01639           //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+c="<<nM<<"="<<rM+nM<<G4endl;
01640           //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0003",
01641           //           FatalException,"Hadronize quasmon: Can't Dec totNuc->QEClu+ResNuc");
01642           G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl;
01643           aParticleChange.ProposeEnergy(kinEnergy);
01644           aParticleChange.ProposeLocalEnergyDeposit(0.);
01645           aParticleChange.ProposeMomentumDirection(dir);
01646           aParticleChange.ProposeTrackStatus(fAlive);
01647           return G4VDiscreteProcess::PostStepDoIt(track,step);
01648         }
01649         // --- End of the moving cluster implementation ---
01650 #ifdef qedebug
01651         G4cout<<"G4QIn::PStDoIt:QEC, RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
01652 #endif
01653        }
01654        G4LorentzVector s4M=n4M+proj4M;             // Tot 4-momentum for scattering
01655        G4double prjM2 = proj4M.m2();               // Squared mass of the projectile
01656        G4double prjM = std::sqrt(prjM2);           // @@ Get from pPDG (?)
01657        G4double minM = prjM+qM;                    // Min mass sum for the final products
01658        G4double cmM2 =s4M.m2();
01659        if(cmM2>minM*minM)
01660        {
01661 #ifdef qedebug
01662         G4cout<<"G4QInel::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl;
01663 #endif
01664         // Estimate and randomize charge-exchange with a quasi-free cluster
01665         G4bool chex=false;                        // Flag of the charge exchange scattering
01666         G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
01667         // This is reserved for the charge-exchange scattering (to help neutron spectras)
01668         //if(cln&&!cp1 &&(projPDG==2212&&rA>rZ || projPDG==2112&&rZ>1))// @@ Use proj chex
01669         if(2>3) // @@ charge exchange is not implemented yet @@
01670         {
01671 #ifdef qedebug
01672           G4cout<<"G4QIn::PStDoIt:-Enter, P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
01673 #endif
01674           G4double tprM=mProt;
01675           G4double tprM2=mProt2;
01676           G4int tprPDG=2212;
01677           G4int tresPDG=restPDG+999;
01678           if(projPDG==2212)
01679           {
01680             projpt=G4Neutron::Neutron();
01681             tprM=mNeut;
01682             tprM2=mNeut2;
01683             tprPDG=2112;
01684             tresPDG=restPDG-999;
01685           }
01686           minM=tprM+qM;
01687           G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM);
01688           G4double efP=std::sqrt(efE*efE-tprM2);
01689           G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
01690 #ifdef qedebug
01691           G4cout<<"G4QInel::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
01692 #endif
01693           if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl))     // minM is redefined
01694           {
01695             projPDG=tprPDG;
01696             prjM=tprM;
01697             G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus
01698             r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
01699             n4M=G4LorentzVector(0.,0.,0.,tM-rM);      // 4mom of the quasi-free cluster
01700             chex=true;                                // Confirm charge exchange scattering
01701           }
01702         }
01703         //
01704         std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M,
01705                                                                          projPDG, proj4M);
01706 #ifdef qedebug
01707         G4cout<<"G4QInelastic::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
01708               <<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
01709 #endif
01710         aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles
01711         // @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not
01712         if(chex)               // ==> Projectile is changed: fill everything to secondaries
01713         {
01714           aParticleChange.ProposeEnergy(0.);                // @@ ??
01715           aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed
01716           aParticleChange.SetNumberOfSecondaries(3); 
01717           G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second);
01718           G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex
01719           scatPrH->SetWeight(weight);                                   //    weighted
01720           scatPrH->SetTouchableHandle(trTouchable);                     //   projectile
01721           aParticleChange.AddSecondary(scatPrH);                        //     hadron
01722         }
01723         else               // ==> The leading particle is filled to the updated projectilee
01724         {
01725           aParticleChange.SetNumberOfSecondaries(2);        // @@ if proj=leading 
01726           G4double ldT=(sctout.second).e()-prjM;            // kin Energy of scat project.
01727           aParticleChange.ProposeEnergy(ldT);               // Change the kin Energy
01728           G4ThreeVector ldV=(sctout.second).vect();         // Change momentum direction
01729           aParticleChange.ProposeMomentumDirection(ldV/ldV.mag());
01730           aParticleChange.ProposeTrackStatus(fAlive);
01731         }
01732         // ---------------------------------------------------------
01733         // Fill scattered quasi-free nucleon/fragment
01734         if     (nPDG==90000001) theDefinition = G4Neutron::Neutron();
01735         else if(nPDG==90001000) theDefinition = G4Proton::Proton();
01736         else if(nZ>0 && nA>1)
01737                   theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);
01738 #ifdef debug
01739         else G4cout<<"-Warning_G4QIn::PSD:scatqfPDG="<<nPDG<<", Z="<<nZ<<",A="<<nA<<G4endl;
01740 #endif
01741         if(nZ>0 && nA>0)
01742         {
01743           G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first);
01744           G4Track* scatQFN = new G4Track(theQFN, localtime, position ); //   scattered
01745           scatQFN->SetWeight(weight);                                   //    weighted
01746           scatQFN->SetTouchableHandle(trTouchable);                     //   quasi-free
01747           aParticleChange.AddSecondary(scatQFN);                        //  nucleon/cluster
01748         }
01749         // ----------------------------------------------------
01750         // Fill residual nucleus
01751         if     (restPDG==90000001) theDefinition = G4Neutron::Neutron();
01752         else if(restPDG==90001000) theDefinition = G4Proton::Proton();
01753         else if(rZ>0 && rA>1)
01754                   theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
01755 #ifdef debug
01756         else G4cout<<"-Warning_G4QIn::PSD: resPDG="<<restPDG<<",Z="<<rZ<<",A="<<rA<<G4endl;
01757 #endif
01758         if(rZ>0 && rA>0)
01759         {
01760           G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
01761           G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
01762           scatReN->SetWeight(weight);                                   //    weighted
01763           scatReN->SetTouchableHandle(trTouchable);                     //    residual
01764           aParticleChange.AddSecondary(scatReN);                        //    nucleus
01765         }
01766         delete output;
01767         return G4VDiscreteProcess::PostStepDoIt(track, step);
01768        }
01769 #ifdef qedebug
01770        else G4cout<<"G4QInel::PSD:OUT,M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
01771 #endif
01772       } // end of proton quasi-elastic (including QE on NF)
01773       else // @@ make cost-condition @@ Pickup process (pickup 1 or 2 n and make d or t)
01774       {
01775        if(projPDG==2212) restPDG--; // Res. nucleus PDG is one neutron less than targ. PDG
01776        else
01777        {
01778          restPDG-=1000;          // Res. nucleus PDG is one proton less than targ. PDG
01779          rZ--;                   // Reduce the charge of the residual nucleus
01780        }
01781        G4double mPUF=mDeut;      // Prototype of the mass of the created pickup fragment
01782        G4ParticleDefinition* theDefinition = G4Deuteron::Deuteron(); // Default: make d
01783        //theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
01784        G4bool din=false;         // Flag of picking up 2 neutrons to create t (tritium)(p)
01785        G4bool dip=false;         // Flag of picking up 2 protons to create t (tritium) (n)
01786        G4bool pin=false;         // Flag of picking up 1 n + 1 p to create He3/t      (p/n)
01787        G4bool hin=false;         // Flag of PickUp creation of alpha (He4)            (p/n)
01788        G4double frn=G4UniformRand();
01789        if(N>2 && frn > 0.95)     // .95 is a parameter to pickup 2N (to cteate t or He3)
01790        {
01791          if(projPDG==2212)
01792          {
01793            if(N>1 && G4UniformRand()*(Z+.5*(N-1)) > Z)
01794            {
01795              theDefinition = G4Triton::Triton();
01796              mPUF=mTrit;         // The mass of the created pickup fragment is triton
01797              restPDG--;          // Res. nucleus PDG is two neutrons less than target PDG
01798              din=true;
01799            }
01800            else
01801            {
01802              theDefinition = G4He3::He3();
01803              mPUF=mHel3;         // The mass of the created pickup fragment is Helium3
01804              restPDG-=1000;      // Res. nucleus PDG is two neutrons less than target PDG
01805              rZ--;               // Reduce the charge of the residual nucleus
01806              pin=true;
01807            }
01808          }
01809          else // Neutron projectile (2112)
01810          {
01811            if(Z>1 && G4UniformRand()*(N+.5*(Z-1)) > N)
01812            {
01813              theDefinition = G4He3::He3();
01814              mPUF=mHel3;         // The mass of the created pickup fragment is Helium3
01815              restPDG-=1000;      // Res. nucleus PDG is two neutrons less than target PDG
01816              rZ--;               // Reduce the charge of the residual nucleus
01817              dip=true;
01818            }
01819            else
01820            {
01821              theDefinition = G4Triton::Triton();
01822              mPUF=mTrit;         // The mass of the created pickup fragment is triton
01823              restPDG--;          // Res. nucleus PDG is two neutrons less than target PDG
01824              pin=true;
01825            }
01826          }
01827          rA--;                   // Reduce the baryon number of the residual nucleus
01828          // Recalculate dmom
01829          G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
01830                        fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01831          dmom=sumv.mag();
01832          if(Z>1 && frn > 0.99)  // .99 is a parameter to pickup 2n1p || 1n2p & create alpha
01833          {
01834            theDefinition = G4Alpha::Alpha();
01835            if((din && projPDG==2212) || (pin && projPDG==2112))
01836            {
01837              restPDG-=1000;        // Residual nucleus PDG is reduced to create alpha
01838              rZ--;                 // Reduce the charge of the residual nucleus
01839            }
01840            else if((pin && projPDG==2212) || (dip && projPDG==2112)) restPDG--;
01841            else G4cout<<"-Warning-G4QIn::PSD: PickUp logic error, proj="<<projPDG<<G4endl;
01842            hin=true;
01843            mPUF=mAlph;           // The mass of the created pickup fragment (alpha)
01844            rA--;                 // Reduce the baryon number of the residual nucleus
01845            // Recalculate dmom
01846            sumv += fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01847            dmom=sumv.mag();
01848          }
01849        }
01850        G4double mPUF2=mPUF*mPUF;
01851        G4LorentzVector t4M(0.,0.,0.,tM);     // 4m of the target nucleus to be decayed
01852        G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
01853        G4double rM2=rM*rM;                   // Squared mass of the residual nucleus
01854        G4double nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom);// M2 of boundQF-nucleon(2N)
01855        if(nM2 < 0) nM2=0.;
01856        G4double den2=(dmom*dmom+nM2);        // squared energy of bQF-nucleon
01857        G4double den=std::sqrt(den2);         // energy of bQF-nucleon
01858 #ifdef qedebug
01859        G4double nM=std::sqrt(nM2);           // M of bQF-nucleon
01860        G4cout<<"G4QInel::PStDoIt:PiU, p="<<dmom<<",tM="<<tM<<", R="<<rM<<",N="<<nM<<G4endl;
01861 #endif
01862        G4double qp=momentum*dmom;
01863        G4double srm=proj4M.e()*den;
01864        G4double mNucl2=mProt2;
01865        if(projPDG == 2112) mNucl2=mNeut2;
01866        G4double cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
01867        G4int ict=0;
01868        while(std::fabs(cost)>1. && ict<3)
01869        {
01870          dmom=91.;                           // Fermi momentum (proDefault for a deuteron)
01871          if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
01872          if(din || pin || dip)               // Recalculate dmom for the final t/He3
01873          {
01874           G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
01875                         fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01876           if(hin)
01877                   sumv+=fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
01878           dmom=sumv.mag();
01879          }
01880          nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom); // M2 of bQF-nucleon
01881          if(nM2 < 0) nM2=0.;
01882          //nM=std::sqrt(nM2);                  // M of bQF-nucleon --Not used?--
01883          den2=(dmom*dmom+nM2);               // squared energy of bQF-nucleon
01884          den=std::sqrt(den2);                // energy of bQF-nucleon
01885          qp=momentum*dmom;
01886          srm=proj4M.e()*den;
01887          cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
01888          ict++;
01889 #ifdef ldebug
01890          if(ict>2)G4cout<<"G4QInelast::PStDoIt:i="<<ict<<",d="<<dmom<<",ct="<<cost<<G4endl;
01891 #endif
01892        }
01893        if(std::fabs(cost)<=1.)
01894        {// ---- @@ From here can be a MF for QF nucleon extraction (if used by others)
01895         G4ThreeVector vdir = proj4M.vect();  // 3-Vector in the projectile direction
01896         G4ThreeVector vx(0.,0.,1.);          // ProtoOrt in theDirection of the projectile
01897         G4ThreeVector vy(0.,1.,0.);          // First ProtoOrt orthogonal to the direction
01898         G4ThreeVector vz(1.,0.,0.);          // SecondProtoOrt orthoganal to the direction
01899         if(vdir.mag2() > 0.)                 // the projectile isn't at rest
01900         {
01901           vx = vdir.unit();                  // Ort in the direction of the projectile
01902           G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
01903           vy = vv.unit();                    // First ort orthogonal to the proj direction
01904           vz = vx.cross(vy);                 // Second ort orthoganal to the projDirection
01905         }
01906         // ---- @@ End of possible MF (Similar is in G4QEnvironment)
01907         G4double tdom=dmom*std::sqrt(1-cost*cost);// Transversal component of the Fermi-Mom
01908         G4double phi=twopi*G4UniformRand();  // Phi of the Fermi-Mom
01909         G4ThreeVector vedm=vx*dmom*cost+vy*tdom*std::sin(phi)+vz*tdom*std::cos(phi);// bQFN
01910         G4LorentzVector bqf(vedm,den);      // 4-mom of the bQF nucleon
01911         r4M=t4M-bqf;                         // 4mom of the residual nucleus
01912 #ifdef debug
01913         if(std::fabs(r4M.m()-rM)>.001) G4cout<<"G4QIn::PSD: rM="<<rM<<"#"<<r4M.m()<<G4endl;
01914 #endif
01915         G4LorentzVector f4M=proj4M+bqf;      // Prototype of 4-mom of Deuteron
01916 #ifdef debug
01917         if(std::fabs(f4M.m()-mPUF)>.001)G4cout<<"G4QI::PSD:m="<<mPUF<<"#"<<f4M.m()<<G4endl;
01918 #endif
01919         aParticleChange.ProposeEnergy(0.);   // @@ ??
01920         aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed
01921         aParticleChange.SetNumberOfSecondaries(2); 
01922         // Fill pickup hadron
01923         G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,f4M);
01924         G4Track* scatQFN = new G4Track(theQFN, localtime, position ); //    pickup
01925         scatQFN->SetWeight(weight);                                   //    weighted
01926         scatQFN->SetTouchableHandle(trTouchable);                     //    nuclear
01927         aParticleChange.AddSecondary(scatQFN);                        //    cluster
01928 #ifdef pickupd
01929         G4cout<<"G4QInelastic::PStDoIt: f="<<theDefinition<<",f4M="<<f4M.m()<<f4M<<G4endl;
01930 #endif
01931         // ----------------------------------------------------
01932         // Fill residual nucleus
01933         if     (restPDG==90000001) theDefinition = G4Neutron::Neutron();
01934         else if(restPDG==90001000) theDefinition = G4Proton::Proton();
01935         else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
01936         G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
01937         G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
01938         scatReN->SetWeight(weight);                                   //    weighted
01939         scatReN->SetTouchableHandle(trTouchable);                     //    residual
01940         aParticleChange.AddSecondary(scatReN);                        //    nucleus
01941 #ifdef pickupd
01942         G4cout<<"G4QInelastic::PStDoIt:rZ="<<rZ<<",rA="<<rA<<",r4M="<<r4M.m()<<r4M<<G4endl;
01943 #endif
01944         delete output;
01945         return G4VDiscreteProcess::PostStepDoIt(track, step);
01946          }
01947       }
01948     } // end of the quasi-elastic decoupled process for the neucleon projectile
01949   }  // end lepto-nuclear, neutrino-nuclear, proton/neutron decoupled processes
01950   EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM);    // Total 4-mom of the reaction
01951   if(absMom) EnMomConservation+=lead4M;         // Add E/M of leading System
01952 #ifdef debug
01953   G4cout<<"G4QInelastic::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl;
01954 #endif
01955 
01956   // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
01957   //G4bool elF=false; // Flag of the ellastic scattering is "false" by default
01958   //G4double eWei=1.;
01959   // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End  
01960 #ifdef debug
01961   G4cout<<"^^G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
01962 #endif
01963   const G4QNucleus targNuc(Z,N);                          // Define the target nucleus
01964   if(projPDG>9000000)
01965   {
01966     delete output;                                        // Before the output creation
01967     G4QNucleus projNuc(proj4M,projPDG);                   // Define the projectile nucleus
01968     G4QIonIonCollision IonIon(projNuc, targNuc);  // Define deep-inelastic ion-ion reaction
01969 #ifdef debug
01970     G4cout<<"G4QInel::PStDoIt: ProjNuc="<<projPDG<<proj4M<<", TargNuc="<<targPDG<<G4endl;
01971 #endif
01972     try {output = IonIon.Fragment();}                     // DINR AA interaction products
01973     catch (G4QException& error)
01974     {
01975       G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
01976       // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
01977       G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027",
01978                   FatalException, "CHIPS hA crash");
01979     }
01980   }
01981   else                                                    // --> A projectile hadron
01982   {
01983     // *** Let to produce elastic final states for acceleration ***
01984 #ifdef debug
01985     G4int maxCn=7;
01986 #endif
01987     G4int atCn=0;                                         // Attempts counter
01988     //G4bool inel=false;
01989     //while (inel==false && ++atCn <= maxCn)                // To evoid elastic final state
01990     //{
01991 #ifdef debug
01992       G4cout<<"G4QInelastic::PStDoIt: atCn="<<atCn<<" <= maxCn="<<maxCn<<G4endl;
01993 #endif
01994       G4int outN=output->size();
01995       if(outN)                                            // The output is not empty
01996       {
01997         std::for_each(output->begin(), output->end(), DeleteQHadron());
01998         output->clear();
01999       }
02000       delete output;                                      // Before the new output creation
02001       const G4QHadron projH(projPDG,proj4M);              // Define the projectile particle
02002 #ifdef debug
02003       G4cout<<"G4QInelastic::PStDoIt: Before Fragmentation"<<G4endl;
02004 #endif
02005       //if(aProjPDG==22 && projPDG==22 && proj4M.m2()<.01 && proj4M.e()<150.)
02006       //{
02007       //  G4QHadron* tgH= new G4QHadron(targNuc.GetPDGCode(),targNuc.Get4Momentum());
02008       //  G4QHadron* prH= new G4QHadron(projPDG,proj4M);
02009       //  output = new G4QHadronVector();
02010       //  output->push_back(prH);
02011       //  output->push_back(tgH);
02012       //}
02013       //else
02014       //{
02015         G4QFragmentation DINR(targNuc, projH);            // Define deep-inel. h-A reaction
02016 #ifdef debug
02017         G4cout<<"G4QInelastic::PStDoIt:Proj="<<projPDG<<proj4M<<",TgNuc="<<targPDG<<G4endl;
02018 #endif
02019         try {output = DINR.Fragment();}                   // DINR hA fragmentation
02020         catch (G4QException& error)
02021         {
02022           G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
02023           // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
02024           G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027",
02025                       FatalException, "CHIPS hA crash");
02026         }
02027         //}
02028       outN=output->size();
02029 #ifdef debug
02030       G4cout<<"G4QInelastic::PostStepDoIt:  At# "<<atCn<<", nSec="<<outN<<", fPDG="
02031             <<(*output)[0]->GetPDGCode()<<", pPDG="<< projPDG<<G4endl;
02032 #endif
02033       //inel=true;                                        // For while
02034       if(outN < 2)
02035       {
02036         G4cout<<"-Warning-G4QInelastic::PostStepDoIt: nSec="<<outN<<", At# "<<atCn<<G4endl;
02037         //inel=false;                                     // For while
02038       }
02039       //else if(outN==2 && (*output)[0]->GetPDGCode() == projPDG) inel=false; // For while
02040 #ifdef debug
02041       if(atCn==maxCn)G4cout<<"-Warning-G4QI::PostStDoIt:mAt="<<atCn<<" is reached"<<G4endl;
02042 #endif
02043     //}
02044   }
02045   //
02046   // --- the scattered hadron with changed nature can be added here ---
02047   if(scat)
02048   {
02049     G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M);
02050     output->push_back(scatHadron);
02051   }
02052   G4int qNH=0;
02053   if(leadhs) qNH=leadhs->size();
02054   if(absMom)
02055   {
02056     if(qNH) for(G4int iq=0; iq<qNH; iq++)
02057     {
02058       G4QHadron* loh=(*leadhs)[iq];   // Pointer to the output hadron
02059       output->push_back(loh);
02060     }
02061     if(leadhs) delete leadhs;
02062     leadhs=0;
02063   }
02064   else
02065   {
02066     if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
02067     if(leadhs) delete leadhs;
02068     leadhs=0;
02069   }
02070   // ------------- From here the secondaries are filled -------------------------
02071   G4int tNH = output->size();       // A#of hadrons in the output
02072   aParticleChange.SetNumberOfSecondaries(tNH); // @@ tNH can be changed (drop/decay!)
02073   // Now add nuclear fragments
02074 #ifdef debug
02075   G4cout<<"G4QInelastic::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
02076 #endif
02077 #ifdef ppdebug
02078   if(absMom)G4cout<<"G4QInelastic::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
02079 #endif
02080   G4int nOut=output->size();               // Real length of the output @@ Temporary
02081   if(tNH==1 && !scat)                      // @@ Temporary. Find out why it happened!
02082   {
02083     G4cout<<"-Warning-G4QInelastic::PostStepDoIt: 1 secondary! absMom="<<absMom;
02084     if(absMom) G4cout<<", qNH="<<qNH;
02085     G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
02086     G4cout<<G4endl;
02087     tNH=0;
02088     delete output->operator[](0);          // delete the creazy hadron
02089     output->pop_back();                    // clean up the output vector
02090   }
02091   if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QInelastic::PostStepDoIt: 2 # "<<nOut<<G4endl;
02092   // Deal with ParticleChange final state interface to GEANT4 output of the process
02093   //if(tNH==2) for(i=0; i<tNH; i++)          // @@ Temporary tNH==2 instead of just tNH
02094   if(tNH) for(i=0; i<tNH; i++)             // @@ Temporary tNH==2 instead of just tNH
02095   {
02096     // Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
02097     // Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
02098     // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
02099     G4QHadron* hadr=(*output)[i];          // Pointer to the output hadron    
02100     G4int PDGCode = hadr->GetPDGCode();
02101     G4int nFrag   = hadr->GetNFragments();
02102 #ifdef pdebug
02103     G4cout<<"G4QInelastic::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag
02104           <<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
02105 #endif
02106     if(nFrag)                              // Skip intermediate (decayed) hadrons
02107     {
02108 #ifdef debug
02109       G4cout<<"G4QInelastic::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
02110 #endif
02111       delete hadr;
02112       continue;
02113     }
02114     G4DynamicParticle* theSec = new G4DynamicParticle;
02115     G4ParticleDefinition* theDefinition;
02116     if     (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
02117     else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
02118     else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
02119     else if(PDGCode==311 || PDGCode==-311)
02120     {
02121       if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong();   // K_L
02122       else                   theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
02123     }
02124     else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
02125     else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
02126     else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
02127     else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
02128     else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
02129     else if(PDGCode==90000000)
02130     {
02131 #ifdef pdebug
02132       G4cout<<"-Warning-G4QInelastic::PostStepDoIt:Vacuum PDG="<<PDGCode
02133             <<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
02134 #endif
02135       theDefinition = G4Gamma::Gamma();
02136       G4double hM2=(hadr->Get4Momentum()).m2();
02137       if(std::fabs(hM2)>.01) G4cout<<"-Warning-G4QInel::PoStDoIt:90000000M2="<<hM2<<G4endl;
02138       else if(hadr->Get4Momentum()==vacuum4M)
02139       {
02140         delete hadr;
02141         continue;
02142       }
02143     }
02144     else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
02145     {
02146       G4int aZ = hadr->GetCharge();
02147       G4int aA = hadr->GetBaryonNumber();
02148 #ifdef pdebug
02149       G4cout<<"G4QInelastic::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
02150 #endif
02151       theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
02152     }
02153     //else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
02154     else if(PDGCode==89999998 || PDGCode==89998000 || PDGCode==88000000) // anti-di-baryon
02155     {
02156       G4double rM=mNeut;                                // Prototype of the baryon mass
02157       G4int  rPDG=-2112;                                // Prototype of the baryon PDG
02158       G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition
02159       if     (PDGCode==89998000)
02160       {
02161         rM=mProt;
02162         rPDG=-2212;
02163         BarDef=G4Proton::Proton();
02164       }
02165       else if(PDGCode==88000000)
02166       {
02167         rM=mLamb;
02168         rPDG=-3122;
02169         BarDef=G4Lambda::Lambda();
02170       }
02171       G4LorentzVector t4M=hadr->Get4Momentum();         // 4m of the di-baryon
02172       G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon
02173       G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon
02174 #ifdef qedebug
02175       G4cout<<"G4QInel::PStDoIt:AntiDiBar,t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl;
02176 #endif
02177       if(!G4QHadron(t4M).DecayIn2(f4M, s4M))
02178       {
02179         G4cerr<<"G4QIn::PostStDoIt: ADB, M="<<t4M.m()<<" < 2*rM="<<rM<<" = "<<2*rM<<G4endl;
02180         // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-dibaryon");
02181         G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0004",
02182                     FatalException, "Hadronize quasmon: Can't decay anti-dibaryon");
02183       }
02184       // --- End of the moving cluster implementation ---
02185 #ifdef qedebug
02186       G4cout<<"G4QInelastic::PStDoIt: ADB, H1="<<rM<<f4M<<", H2="<<rM<<s4M<<G4endl;
02187 #endif
02188       G4QHadron fH(rPDG,f4M);
02189       hadr->Set4Momentum(f4M);
02190       hadr->SetQPDG(fH.GetQPDG());
02191       theDefinition=BarDef;
02192 #ifdef debug
02193       G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h1="<<rPDG<<f4M<<G4endl;
02194 #endif
02195       G4QHadron* sH = new G4QHadron(rPDG,s4M);
02196       output->push_back(sH);
02197       ++tNH;
02198 #ifdef debug
02199       G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
02200 #endif
02201     }
02202     else if(PDGCode==90000997 || PDGCode==89997001) // anti-NDelta
02203     {
02204       G4double rM=mNeut;                                // Prototype of the baryon mass
02205       G4int  rPDG=-2112;                                // Prototype of the baryon PDG
02206       G4double iM=mPi;                                  // Prototype of the pion mass
02207       G4int  iPDG= 211;                                 // Prototype of the pion PDG
02208       G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition
02209       if(PDGCode==90000997)
02210       {
02211         rM=mProt;
02212         rPDG=-2212;
02213         iPDG=-211;
02214         BarDef=G4Proton::Proton();
02215       }
02216       G4LorentzVector t4M=hadr->Get4Momentum();         // 4m of the di-baryon
02217       G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon
02218       G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon
02219       G4LorentzVector u4M=G4LorentzVector(0.,0.,0.,iM); // 4mom of the pion
02220 #ifdef qedebug
02221       G4cout<<"G4QInel::PStDoIt:AntiNDelta, t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl;
02222 #endif
02223       if(!G4QHadron(t4M).DecayIn3(f4M, s4M, u4M))
02224       {
02225         G4cerr<<"G4QIn::PostStDoIt: AND, tM="<<t4M.m()<<" < 2*mB+mPi="<<2*rM+iM<<G4endl;
02226         // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-NDelta");
02227         G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0005",
02228                     FatalException, "Hadronize quasmon: Can't decay anti-NDelta");
02229       }
02230       // --- End of the moving cluster implementation ---
02231 #ifdef qedebug
02232       G4cout<<"G4QInel::PStDoI:AND,B1="<<rM<<f4M<<",B2="<<rM<<s4M<<",Pi="<<iM<<u4M<<G4endl;
02233 #endif
02234       G4QHadron fH(rPDG,f4M);
02235       hadr->Set4Momentum(f4M);
02236       hadr->SetQPDG(fH.GetQPDG());
02237       theDefinition=BarDef;
02238 #ifdef debug
02239       G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h1="<<rPDG<<f4M<<G4endl;
02240 #endif
02241       G4QHadron* sH = new G4QHadron(rPDG,s4M);
02242       output->push_back(sH);
02243       ++tNH;
02244 #ifdef debug
02245       G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
02246 #endif
02247       G4QHadron* uH = new G4QHadron(iPDG,u4M);
02248       output->push_back(uH);
02249       ++tNH;
02250 #ifdef debug
02251       G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
02252 #endif
02253     }
02254     else
02255     {
02256 #ifdef pdebug
02257       G4cout<<"G4QInelastic::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
02258 #endif
02259       theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
02260 #ifdef pdebug
02261       G4cout<<"G4QInelastic::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
02262 #endif
02263     }
02264     if(!theDefinition)
02265     { // !! Commenting of this print costs days of searching for E/mom nonconservation !!
02266       if(PDGCode!=90000000 || hadr->Get4Momentum()!=vacuum4M)
02267         G4cout<<"---Warning---G4QInelastic::PostStepDoIt: drop PDG="<<PDGCode<<", 4Mom="
02268               <<hadr->Get4Momentum()<<G4endl;
02269       delete hadr;
02270       continue;
02271     }
02272 #ifdef pdebug
02273     G4cout<<"G4QInelastic::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
02274 #endif
02275     theSec->SetDefinition(theDefinition);
02276     G4LorentzVector h4M=hadr->Get4Momentum();
02277     EnMomConservation-=h4M;
02278 #ifdef tdebug
02279     G4cout<<"G4QInelast::PSDI: "<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
02280 #endif
02281 #ifdef debug
02282     G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
02283 #endif
02284     theSec->Set4Momentum(h4M); //                                                         ^
02285     delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
02286 #ifdef debug
02287     G4ThreeVector curD=theSec->GetMomentumDirection();              //                    ^
02288     G4double curM=theSec->GetMass();                                //                    |
02289     G4double curE=theSec->GetKineticEnergy()+curM;                  //                    ^
02290     G4cout<<"G4QInelast::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;//|
02291 #endif
02292     G4Track* aNewTrack = new G4Track(theSec, localtime, position ); //                    ^
02293     aNewTrack->SetWeight(weight);                                   //    weighted        |
02294     aNewTrack->SetTouchableHandle(trTouchable);                     //                    |
02295     aParticleChange.AddSecondary( aNewTrack );                      //                    |
02296 #ifdef debug
02297     G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<" is done."<<G4endl; //                    |
02298 #endif
02299   } //                                                                                    |
02300   delete output; // instances of the G4QHadrons from the output are already deleted above +
02301   if(leadhs)     // To satisfy Valgrind ( How can that be?)
02302   {
02303     qNH=leadhs->size();
02304     if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
02305     delete leadhs;
02306     leadhs=0;
02307   }
02308 #ifdef debug
02309   G4cout<<"G4QInelastic::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
02310 #endif
02311   if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
02312     aParticleChange.ProposeTrackStatus(fStopAndKill);        // Kill the absorbed particle
02313 #ifdef pdebug
02314     G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()
02315           <<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
02316 #endif
02317 #ifdef ldebug
02318   G4cout<<"G4QInelastic::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St="
02319         <<aParticleChange.GetTrackStatus()<<G4endl;
02320 #endif
02321   return G4VDiscreteProcess::PostStepDoIt(track, step);
02322 }
02323 
02324 std::pair<G4double,G4double> G4QInelastic::Random2DDirection()
02325 {
02326   G4double sp=0;              // sin(phi)
02327   G4double cp=1.;             // cos(phi)
02328   G4double r2=2.;             // to enter the loop
02329   while(r2>1. || r2<.0001)    // pi/4 efficiency
02330   {
02331     G4double sine=G4UniformRand();
02332     G4double cosine=G4UniformRand();
02333     sp=1.-sine-sine;
02334     cp=1.-cosine-cosine;
02335     r2=sp*sp+cp*cp;
02336   }
02337   G4double norm=std::sqrt(r2);
02338   return std::make_pair(sp/norm,cp/norm);
02339 }

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