Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QLowEnergy.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id$
27// GEANT4 tag $Name: geant4-09-04 $
28//
29// ---------------- G4QLowEnergy class -----------------
30// by Mikhail Kossov, Aug 2007.
31// G4QLowEnergy class of the CHIPS Simulation Branch in GEANT4
32// ---------------------------------------------------------------
33// Short description: This is a fast low energy algorithm for the
34// inelastic interactions of nucleons and nuclei (ions) with nuclei.
35// This is a fase-space algorithm, but not quark level. Provides
36// nuclear fragments upto alpha only. Never was tumed (but can be).
37// ---------------------------------------------------------------
38//#define debug
39//#define pdebug
40//#define edebug
41//#define tdebug
42//#define nandebug
43//#define ppdebug
44
45#include "G4QLowEnergy.hh"
46#include "G4SystemOfUnits.hh"
47#include "G4HadronicDeprecate.hh"
48
49
50// Initialization of static vectors
51//G4int G4QLowEnergy::nPartCWorld=152; // #of particles initialized in CHIPS
52//G4int G4QLowEnergy::nPartCWorld=122; // #of particles initialized in CHIPS
53G4int G4QLowEnergy::nPartCWorld=85; // #of particles initialized in CHIPS Reduced
54std::vector<G4int> G4QLowEnergy::ElementZ; // Z of element(i) in theLastCalc
55std::vector<G4double> G4QLowEnergy::ElProbInMat; // SumProbOfElem in Material
56std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;// N of isotope(j), E(i)
57std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;//SumProbIsotE(i)
58
59// Constructor
60G4QLowEnergy::G4QLowEnergy(const G4String& processName):
61 G4VDiscreteProcess(processName, fHadronic), evaporate(true)
62{
63 G4HadronicDeprecate("G4QLowEnergy");
64
65#ifdef debug
66 G4cout<<"G4QLowEnergy::Constructor is called processName="<<processName<<G4endl;
67#endif
68 if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created "<<G4endl;
69 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
70}
71
72// Destructor
73G4QLowEnergy::~G4QLowEnergy() {}
74
75
76G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;}
77
78G4int G4QLowEnergy::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
79
80// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
81// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
82// ********** All CHIPS cross sections are calculated in the surface units ************
83G4double G4QLowEnergy::GetMeanFreePath(const G4Track&Track, G4double, G4ForceCondition*F)
84{
85 *F = NotForced;
86 const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
87 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
88 if( !IsApplicable(*incidentParticleDefinition))
89 G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl;
90 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
91 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
92#ifdef debug
93 G4double KinEn = incidentParticle->GetKineticEnergy();
94 G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
95#endif
96 const G4Material* material = Track.GetMaterial(); // Get the current material
97 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
98 const G4ElementVector* theElementVector = material->GetElementVector();
99 G4int nE=material->GetNumberOfElements();
100#ifdef debug
101 G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl;
102#endif
103 G4int pPDG=0;
104 G4int Z=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge());
105 G4int A=incidentParticleDefinition->GetBaryonNumber();
106 if ( incidentParticleDefinition == G4Proton::Proton() ) pPDG = 2212;
107 else if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020;
108 else if ( incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040;
109 else if ( incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030;
110 else if ( incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030;
111 else if ( incidentParticleDefinition == G4GenericIon::GenericIon() || (Z > 0 && A > 0))
112 {
113 pPDG=incidentParticleDefinition->GetPDGEncoding();
114#ifdef debug
115 G4int prPDG=1000000000+10000*A+10*Z;
116 G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl;
117#endif
118 }
119 else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl;
120 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
121 if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); // just to test
122 Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A
123 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
124 G4double sigma=0.; // Sums over elements for the material
125 G4int IPIE=IsoProbInEl.size(); // How many old elements?
126 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
127 {
128 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
129 SPI->clear();
130 delete SPI;
131 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
132 IsN->clear();
133 delete IsN;
134 }
135 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
136 ElementZ.clear(); // Clear the body vector for Z of Elements
137 IsoProbInEl.clear(); // Clear the body vector for SPI
138 ElIsoN.clear(); // Clear the body vector for N of Isotopes
139 for(G4int i=0; i<nE; ++i)
140 {
141 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
142 Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
143 ElementZ.push_back(Z); // Remember Z of the Element
144 G4int isoSize=0; // The default for the isoVectorLength is 0
145 G4int indEl=0; // Index of non-natural element or 0(default)
146 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
147 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
148#ifdef debug
149 G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
150#endif
151 if(isoSize) // The Element has non-trivial abundance set
152 {
153 indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order
154#ifdef debug
155 G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
156#endif
157 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
158 {
159 std::vector<std::pair<G4int,G4double>*>* newAbund =
160 new std::vector<std::pair<G4int,G4double>*>;
161 G4double* abuVector=pElement->GetRelativeAbundanceVector();
162 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
163 {
164 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
165 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
166 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
167 G4double abund=abuVector[j];
168 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
169#ifdef debug
170 G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl;
171#endif
172 newAbund->push_back(pr);
173 }
174#ifdef debug
175 G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
176#endif
177 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
178 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
179 delete newAbund; // Was "new" in the beginning of the name space
180 }
181 }
182 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
183 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
184 IsoProbInEl.push_back(SPI);
185 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
186 ElIsoN.push_back(IsN);
187 G4int nIs=cs->size(); // A#Of Isotopes in the Element
188#ifdef debug
189 G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
190#endif
191 G4double susi=0.; // sum of CS over isotopes
192 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
193 {
194 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
195 G4int N=curIs->first; // #of Neuterons in the isotope j of El i
196 IsN->push_back(N); // Remember Min N for the Element
197#ifdef debug
198 G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl;
199#endif
200 G4bool ccsf=true; // Extract inelastic Ion-Ion cross-section
201#ifdef debug
202 G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl;
203#endif
204 G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope
205#ifdef debug
206 G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
207 <<Momentum<<", XSec="<<CSI/millibarn<<G4endl;
208#endif
209 curIs->second = CSI;
210 susi+=CSI; // Make a sum per isotopes
211 SPI->push_back(susi); // Remember summed cross-section
212 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
213 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
214#ifdef debug
215 G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
216 <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
217#endif
218 ElProbInMat.push_back(sigma);
219 } // End of LOOP over Elements
220 // Check that cross section is not zero and return the mean free path
221#ifdef debug
222 G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
223#endif
224 if(sigma > 0.000000001) return 1./sigma; // Mean path [distance]
225 return DBL_MAX;
226}
227
228G4bool G4QLowEnergy::IsApplicable(const G4ParticleDefinition& particle)
229{
230 G4int Z=static_cast<G4int>(particle.GetPDGCharge());
231 G4int A=particle.GetBaryonNumber();
232 if (particle == *( G4Proton::Proton() )) return true;
233 else if (particle == *( G4Neutron::Neutron() )) return true;
234 else if (particle == *( G4Deuteron::Deuteron() )) return true;
235 else if (particle == *( G4Alpha::Alpha() )) return true;
236 else if (particle == *( G4Triton::Triton() )) return true;
237 else if (particle == *( G4He3::He3() )) return true;
238 else if (particle == *( G4GenericIon::GenericIon() )) return true;
239 else if (Z > 0 && A > 0) return true;
240#ifdef debug
241 G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<", A="
242 <<A<<", Z="<<Z<<G4endl;
243#endif
244 return false;
245}
246
247G4VParticleChange* G4QLowEnergy::PostStepDoIt(const G4Track& track, const G4Step& step)
248{
249 static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV
250 static const G4double mPro2= mProt*mProt; // CHIPS sq proton Mass
251 static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV
252 static const G4double mNeu2= mNeut*mNeut; // CHIPS sq neutron Mass
253 static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV
254 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV;
255 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV;
256 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV;
257 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV;
258 static const G4double mFm= 250*MeV;
259 static const G4double third= 1./3.;
260 static const G4ThreeVector zeroMom(0.,0.,0.);
261 static G4ParticleDefinition* aGamma = G4Gamma::Gamma();
262 static G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
263 static G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
264 static G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
265 static G4ParticleDefinition* aProton = G4Proton::Proton();
266 static G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
267 static G4ParticleDefinition* aLambda = G4Lambda::Lambda();
268 static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
269 static G4ParticleDefinition* aTriton = G4Triton::Triton();
270 static G4ParticleDefinition* aHe3 = G4He3::He3();
271 static G4ParticleDefinition* anAlpha = G4Alpha::Alpha();
272 static const G4int nCh=26; // #of combinations
273 static G4QNucleus Nuc; // A fake nucleus to call Evaporation
274 //
275 //-------------------------------------------------------------------------------------
276 static G4bool CWinit = true; // CHIPS Warld needs to be initted
277 if(CWinit)
278 {
279 CWinit=false;
280 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
281 }
282 //-------------------------------------------------------------------------------------
283 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
284 const G4ParticleDefinition* particle=projHadron->GetDefinition();
285#ifdef pdebug
286 G4cout<<"G4QLowEnergy::PostStepDoIt: *** Called *** In4M="
287 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
288 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
289#endif
290 //G4ForceCondition cond=NotForced;
291 //GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters?
292#ifdef debug
293 G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
294#endif
295 std::vector<G4Track*> result;
296 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
297 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
298 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes
299 if(std::fabs(Momentum-momentum)>.000001)
300 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
301#ifdef debug
302 G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M,m="
303 <<proj4M<<proj4M.m()<<G4endl;
304#endif
305 if (!IsApplicable(*particle)) // Check applicability
306 {
307 G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl;
308 return 0;
309 }
310 const G4Material* material = track.GetMaterial(); // Get the current material
311 const G4ElementVector* theElementVector = material->GetElementVector();
312 G4int nE=material->GetNumberOfElements();
313#ifdef debug
314 G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
315#endif
316 G4int projPDG=0; // PDG Code prototype for the captured hadron
317 // Not all these particles are implemented yet (see Is Applicable)
318 G4int Z=static_cast<G4int>(particle->GetPDGCharge());
319 G4int A=particle->GetBaryonNumber();
320 if (particle == G4Proton::Proton() ) projPDG= 2212;
321 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
322 else if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020;
323 else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040;
324 else if (particle == G4Triton::Triton() ) projPDG= 1000010030;
325 else if (particle == G4He3::He3() ) projPDG= 1000020030;
326 else if (particle == G4GenericIon::GenericIon() || (Z > 0 && A > 0))
327 {
328 projPDG=particle->GetPDGEncoding();
329#ifdef debug
330 G4int prPDG=1000000000+10000*Z+10*A; // just for the testing purposes
331 G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl;
332#endif
333 }
334 else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl;
335#ifdef pdebug
336 G4int prPDG=particle->GetPDGEncoding();
337 G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
338#endif
339 if(!projPDG)
340 {
341 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
342 return 0;
343 }
344 // Element treatment
345 G4int EPIM=ElProbInMat.size();
346#ifdef debug
347 G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
348#endif
349 G4int i=0;
350 if(EPIM>1)
351 {
352 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
353 for(i=0; i<nE; ++i)
354 {
355#ifdef debug
356 G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl;
357#endif
358 if (rnd<ElProbInMat[i]) break;
359 }
360 if(i>=nE) i=nE-1; // Top limit for the Element
361 }
362 G4Element* pElement=(*theElementVector)[i];
363 G4int tZ=static_cast<G4int>(pElement->GetZ());
364#ifdef debug
365 G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl;
366#endif
367 if(tZ<=0)
368 {
369 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl;
370 if(tZ<0) return 0;
371 }
372 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
373 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
374 G4int nofIsot=SPI->size(); // #of isotopes in the element i
375#ifdef debug
376 G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
377#endif
378 G4int j=0;
379 if(nofIsot>1)
380 {
381 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
382 for(j=0; j<nofIsot; ++j)
383 {
384#ifdef debug
385 G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
386#endif
387 if(rndI < (*SPI)[j]) break;
388 }
389 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
390 }
391 G4int tN =(*IsN)[j]; ; // Randomized number of neutrons
392#ifdef debug
393 G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl;
394#endif
395 if(tN<0)
396 {
397 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl;
398 return 0;
399 }
400 nOfNeutrons=tN; // Remember it for the energy-momentum check
401#ifdef debug
402 G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl;
403#endif
404 if(tN<0)
405 {
406 G4cerr<<"***Warning***G4QLowEnergy::PostStepDoIt: Element with N="<<tN<< G4endl;
407 return 0;
408 }
409 aParticleChange.Initialize(track);
410#ifdef debug
411 G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl;
412#endif
413 G4double weight = track.GetWeight();
414 G4double localtime = track.GetGlobalTime();
415 G4ThreeVector position = track.GetPosition();
416#ifdef debug
417 G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl;
418#endif
419 G4TouchableHandle trTouchable = track.GetTouchableHandle();
420#ifdef debug
421 G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl;
422#endif
423 G4int targPDG=90000000+tZ*1000+tN; // Target PDG code
424 G4QPDGCode targQPDG(targPDG); // To get CHIPS properties
425 G4double tM=targQPDG.GetMass(); // CHIPS target nucleus mass in MeV
426 G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness
427 G4int pZ=static_cast<G4int>(particle->GetPDGCharge()); // Charge of the projectile
428 G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile
429 G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV
430 G4double cosp=-14*Momentum*(tM-pM)/tM/pM; // Asymmetry power for angular distribution
431#ifdef debug
432 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
433 <<","<<tN<<"), cosp="<<cosp<<G4endl;
434#endif
435 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
436 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
437 G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM); // Target's 4-mom
438 G4LorentzVector tot4M =proj4M+targ4M; // Total 4-mom of the reaction
439#ifdef pdebug
440 G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
441#endif
442 EnMomConservation=tot4M; // Total 4mom of reaction for E/M conservation
443 // --- Start of Coulomb barrier check ---
444 G4double dER = tot4M.e() - tM - pM; // Energy of the reaction
445 G4double pA = pZ+pN; // Atomic weight of the projectile (chged blw)
446 G4double tA = tZ+tN; // Atomic weight of the target (changed below)
447 G4double CBE = Nuc.CoulombBarGen(tZ, tA, pZ, pA); // Coulomb Barrier of the reaction
448#ifdef debug
449 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
450 <<","<<tN<<"), dE="<<dER<<", CB="<<CBE<<G4endl;
451#endif
452 if(dER<CBE) // The cross-section iz 0 -> Do Nothing
453 {
454#ifdef debug
455 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Below Coulomb barrier* PDG="<<projPDG
456 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
457#endif
458 //Do Nothing Action insead of the reaction
459 aParticleChange.ProposeEnergy(kinEnergy);
460 aParticleChange.ProposeLocalEnergyDeposit(0.);
461 aParticleChange.ProposeMomentumDirection(dir) ;
462 return G4VDiscreteProcess::PostStepDoIt(track,step);
463 }
464 // --- End of Coulomb barrier check ---
465 // @@ Probably this is not necessary any more
466#ifdef debug
467 G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl;
468#endif
469 G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection
470#ifdef pdebug
471 G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",tZ="<<tZ<<",N="<<tN<<",XS="
472 <<xSec/millibarn<<G4endl;
473#endif
474#ifdef nandebug
475 if(xSec>0. || xSec<0. || xSec==0);
476 else G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
477#endif
478 // @@ check a possibility to separate p, n, or alpha (!)
479 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
480 {
481#ifdef debug
482 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Zero cross-section* PDG="<<projPDG
483 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
484#endif
485 //Do Nothing Action insead of the reaction
486 aParticleChange.ProposeEnergy(kinEnergy);
487 aParticleChange.ProposeLocalEnergyDeposit(0.);
488 aParticleChange.ProposeMomentumDirection(dir) ;
489 return G4VDiscreteProcess::PostStepDoIt(track,step);
490 }
491 // Kill interacting hadron
492 aParticleChange.ProposeEnergy(0.);
493 aParticleChange.ProposeTrackStatus(fStopAndKill);
494 G4int tB=tZ+tN;
495 G4int pB=pZ+pN;
496#ifdef pdebug
497 G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<", tA="<<tB<<", pA="<<pB<<G4endl;
498#endif
499 // algorithm implementation STARTS HERE --- All calculations are in IU --------
500 tA=tB;
501 pA=pB;
502 G4double tR=1.1; // target nucleus R in fm
503 if(tB > 1) tR*=std::pow(tA,third); // in fm
504 G4double pR=1.1; // projectile nucleus R in fm
505 if(pB > 1) pR*=std::pow(pA,third); // in fm
506 G4double R=tR+pR; // total radius
507 G4double R2=R*R;
508 G4int tD=0;
509 G4int pD=0;
510 G4int nAt=0;
511 G4int nAtM=27;
512 G4int nSec = 0;
513 G4double tcM=0.;
514 G4double tnM=1.;
515#ifdef edebug
516 G4int totChg=0;
517 G4int totBaN=0;
518 G4LorentzVector tch4M =tot4M; // Total 4-mom of the reaction
519#endif
520 while((!tD && !pD && ++nAt<nAtM) || tcM<tnM)
521 {
522#ifdef edebug
523 totChg=tZ+pZ;
524 totBaN=tB+pB;
525 tch4M =tot4M;
526 G4cout<<">G4QLEn::PSDI:#"<<nAt<<",tC="<<totChg<<",tA="<<totBaN<<",t4M="<<tch4M<<G4endl;
527#endif
528 G4LorentzVector tt4M=tot4M;
529 G4int resultSize=result.size();
530 if(resultSize)
531 {
532 for(i=0; i<resultSize; ++i) delete result[i];
533 result.clear();
534 }
535 G4double D=std::sqrt(R2*G4UniformRand());
536#ifdef pdebug
537 G4cout<<"G4QLowEn::PSDI: D="<<D<<", tR="<<tR<<", pR="<<pR<<G4endl;
538#endif
539 if(D > std::fabs(tR-pR)) // leading parts can be separated
540 {
541 nSec = 0;
542 G4double DtR=D-tR;
543 G4double DpR=D-pR;
544 G4double DtR2=DtR*DtR;
545 G4double DpR2=DpR*DpR;
546 //G4double DtR3=DtR2*DtR;
547 G4double DpR3=DpR2*DpR;
548 //G4double DtR4=DtR3*DtR;
549 G4double DpR4=DpR3*DpR;
550 G4double tR2=tR*tR;
551 G4double pR2=pR*pR;
552 G4double tR3=tR2*tR;
553 //G4double pR3=pR2*pR;
554 G4double tR4=tR3*tR;
555 //G4double pR4=pR3*pR;
556 G4double HD=16.*D;
557 G4double tF=tA*(6*tR2*(pR2-DtR2)-4*D*(tR3-DpR3)+3*(tR4-DpR4))/HD/tR3;
558 G4double pF=tF;
559 tD=static_cast<G4int>(tF);
560 pD=static_cast<G4int>(pF);
561 //if(G4UniformRand() < tF-tD) ++tD; // Simple solution
562 //if(G4UniformRand() < pF-pD) ++pD;
563 // Enhance alphas solution
564 if(std::fabs(tF-4.) < 1.) tD=4; // @@ 1. is the enhancement parameter
565 else if(G4UniformRand() < tF-tD) ++tD;
566 if(std::fabs(pF-4.) < 1.) pD=4;
567 else if(G4UniformRand() < pF-pD) ++pD;
568 if(tD > tB) tD=tB;
569 if(pD > pB) pD=tB;
570 // @@ Quasi Free is not debugged @@ The following close it
571 if(tD < 1) tD=1;
572 if(pD < 1) pD=1;
573#ifdef pdebug
574 G4cout<<"G4QLE::PSDI:#"<<nAt<<",pF="<<pF<<",tF="<<tF<<",pD="<<pD<<",tD="<<tD<<G4endl;
575#endif
576 G4int pC=0; // charge of the projectile fraction
577 G4int tC=0; // charge of the target fraction
578 if((tD || pD) && tD<tB && pD<pB)// Periferal interaction
579 {
580 if(!tD || !pD) // Quasi-Elastic: B+A->(B-1)+N+A || ->B+N+(A-1)
581 {
582 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
583 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
584 G4int pPDG=2112; // Proto of the nucleon PDG (proton)
585 G4double prM =mNeut; // Proto of the nucleon mass
586 G4double prM2=mNeu2; // Proto of the nucleon sq mass
587 if (!tD) // Quasi-elastic scattering of the proj QF nucleon
588 {
589 if(pD > 1) pD=1;
590 if(!pN || (pZ && pA*G4UniformRand() < pZ) ) // proj QF proton
591 {
592 pPDG=2212;
593 prM=mProt;
594 prM2=mPro2;
595 pC=1;
596 }
597 G4double Mom=0.;
598 G4LorentzVector com4M = targ4M; // Proto of 4mom of compound
599 G4double tgM=targ4M.e();
600 G4double rNM=0.;
601 G4LorentzVector rNuc4M(0.,0.,0.,0.);
602 if(pD==pB)
603 {
604 Mom=proj4M.rho();
605 com4M += proj4M; // Total 4-mom for scattering
606 rNM=prM;
607 }
608 else // It is necessary to split the nucleon (with fermiM)
609 {
610 G4ThreeVector fm=mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
611 rNM=G4QPDGCode(2112).GetNuclMass(pZ-pC, pN, 0);
612 G4double rNE=std::sqrt(fm*fm+rNM*rNM);
613 rNuc4M=G4LorentzVector(fm,rNE);
614 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
615 rNuc4M.boost(boostV);
616 G4LorentzVector qfN4M=proj4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
617 com4M += qfN4M; // Calculate Total 4Mom for NA scattering
618 G4double pNE = qfN4M.e(); // Energy of the QF nucleon in LS
619 if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2); // Mom(s) fake value
620 else break; // Break the while loop
621 }
622 xSec=0.;
623 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
624 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
625 if( xSec <= 0. ) break; // Break the while loop
626 G4double mint=0.; // Prototype of functional randomized -t
627 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
628 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
629 G4double maxt=0.; // Prototype of maximum -t
630 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
631 else maxt=NCSmanager->GetHMaxT(); // maximum -t
632 G4double cost=1.-mint/maxt; // cos(theta) in CMS
633 if(cost>1. || cost<-1.) break; // Break the while loop
634 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM); // 4mom of recoil target
635 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
636 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
637 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
638 {
639 G4cout<<"G4QLE::Pt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
640 break; // Break the while loop
641 }
642 G4Track* projSpect = 0;
643 G4Track* aNucTrack = 0;
644 if(pB > pD) // Fill the proj residual nucleus
645 {
646 G4int rZ=pZ-pC;
647 G4int rA=pB-1;
648 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
649 if(rA==1)
650 {
651 if(!rZ) theDefinition = aNeutron;
652 else theDefinition = aProton;
653 }
654 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
655 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
656 projSpect = new G4Track(resN, localtime, position );
657 projSpect->SetWeight(weight); // weighted
658 projSpect->SetTouchableHandle(trTouchable);
659#ifdef pdebug
660 G4cout<<"G4QLowEn::PSDI:-->ProjQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
661#endif
662 ++nSec;
663 }
664 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
665 if(pPDG==2212) theDefinition = G4Proton::Proton();
666 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
667 aNucTrack = new G4Track(scatN, localtime, position );
668 aNucTrack->SetWeight(weight); // weighted
669 aNucTrack->SetTouchableHandle(trTouchable);
670#ifdef pdebug
671 G4cout<<"G4QLowEn::PSDI:-->TgQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
672#endif
673 ++nSec;
674 G4Track* aFraTrack=0;
675 if(tA==1)
676 {
677 if(!tZ) theDefinition = aNeutron;
678 else theDefinition = aProton;
679 }
680 else if(tA==8 && tC==4) // Be8 decay
681 {
682 theDefinition = anAlpha;
683 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
684 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
685 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
686 {
687 G4cout<<"*G4QLE::TS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
688 }
689 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
690 aFraTrack = new G4Track(pAl, localtime, position );
691 aFraTrack->SetWeight(weight); // weighted
692 aFraTrack->SetTouchableHandle(trTouchable);
693#ifdef pdebug
694 G4cout<<"G4QLowEn::PSDI:-->TgRQFA4M="<<f4M<<G4endl;
695#endif
696 ++nSec;
697 reco4M=s4M;
698 }
699 else if(tA==5 && (tC==2 || tC==3)) // He5/Li5 decay
700 {
701 theDefinition = aProton;
702 G4double mNuc = mProt;
703 if(tC==2)
704 {
705 theDefinition = aNeutron;
706 mNuc = mNeut;
707 }
708 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
709 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
710 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
711 {
712 G4cout<<"*G4QLE::TS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
713 }
714 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
715 aFraTrack = new G4Track(pAl, localtime, position );
716 aFraTrack->SetWeight(weight); // weighted
717 aFraTrack->SetTouchableHandle(trTouchable);
718#ifdef pdebug
719 G4cout<<"G4QLowEn::PSDI:-->TgQFRN4M="<<f4M<<G4endl;
720#endif
721 ++nSec;
722 theDefinition = anAlpha;
723 reco4M=s4M;
724 }
725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(tZ,tB,0,tZ);
726 ++nSec;
727#ifdef pdebug
728 G4cout<<"G4QLowEn::PSDI:-->PQF_nSec="<<nSec<<G4endl;
729#endif
730 aParticleChange.SetNumberOfSecondaries(nSec);
731 if(projSpect) aParticleChange.AddSecondary(projSpect);
732 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
733 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
734 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
735 G4Track* aResTrack = new G4Track(resA, localtime, position );
736 aResTrack->SetWeight(weight); // weighted
737 aResTrack->SetTouchableHandle(trTouchable);
738#ifdef pdebug
739 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
740#endif
741 aParticleChange.AddSecondary(aResTrack);
742 }
743 else // !pD : QF target Nucleon on the whole Projectile
744 {
745 if(tD > 1) tD=1;
746 if(!tN || (tZ && tA*G4UniformRand() < tZ) ) // target QF proton
747 {
748 pPDG=2212;
749 prM=mProt;
750 prM2=mPro2;
751 tC=1;
752 }
753 G4double Mom=0.;
754 G4LorentzVector com4M=proj4M; // Proto of 4mom of compound
755 prM=proj4M.m();
756 G4double rNM=0.;
757 G4LorentzVector rNuc4M(0.,0.,0.,0.);
758 if(tD==tB)
759 {
760 Mom=prM*proj4M.rho()/proj4M.m();
761 com4M += targ4M; // Total 4-mom for scattering
762 rNM=prM;
763 }
764 else // It is necessary to split the nucleon (with fermiM)
765 {
766 G4ThreeVector fm=250.*std::pow(G4UniformRand(),third)*G4RandomDirection();
767 rNM=G4QPDGCode(2112).GetNuclMass(tZ-tC, tN, 0)/MeV;
768 G4double rNE=std::sqrt(fm*fm+rNM*rNM);
769 rNuc4M=G4LorentzVector(fm,rNE);
770 G4LorentzVector qfN4M=targ4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
771 com4M += qfN4M; // Calculate Total 4Mom for NA scattering
772 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
773 qfN4M.boost(-boostV);
774 G4double tNE = qfN4M.e(); // Energy of the QF nucleon in LS
775 if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value
776 else break; // Break the while loop
777 }
778 xSec=0.;
779 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
780 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
781 if( xSec <= 0. ) break; // Break the while loop
782 G4double mint=0.; // Prototype of functional randomized -t
783 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
784 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
785 G4double maxt=0.; // Prototype of maximum -t
786 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
787 else maxt=NCSmanager->GetHMaxT(); // maximum -t
788 G4double cost=1.-mint/maxt; // cos(theta) in CMS
789 if(cost>1. || cost<-1.) break; // Break the while loop
790 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM); // 4mom of recoil target
791 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
792 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
793 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
794 {
795 G4cout<<"G4QLE::Tt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
796 break; // Break the while loop
797 }
798 G4Track* targSpect = 0;
799 G4Track* aNucTrack = 0;
800 if(tB > tD) // Fill the residual nucleus
801 {
802 G4int rZ=tZ-tC;
803 G4int rA=tB-1;
804 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
805 if(rA==1)
806 {
807 if(!rZ) theDefinition = aNeutron;
808 else theDefinition = aProton;
809 }
810 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
811 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
812 targSpect = new G4Track(resN, localtime, position );
813 targSpect->SetWeight(weight); // weighted
814 targSpect->SetTouchableHandle(trTouchable);
815#ifdef pdebug
816 G4cout<<"G4QLowEn::PSDI:-->TargQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
817#endif
818 ++nSec;
819 }
820 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
821 if(pPDG==2212) theDefinition = G4Proton::Proton();
822 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
823 aNucTrack = new G4Track(scatN, localtime, position );
824 aNucTrack->SetWeight(weight); // weighted
825 aNucTrack->SetTouchableHandle(trTouchable);
826#ifdef pdebug
827 G4cout<<"G4QLowEn::PSDI:-->PrQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
828#endif
829 ++nSec;
830 G4Track* aFraTrack=0;
831 if(pA==1)
832 {
833 if(!pZ) theDefinition = aNeutron;
834 else theDefinition = aProton;
835 }
836 else if(pA==8 && pC==4) // Be8 decay
837 {
838 theDefinition = anAlpha;
839 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
840 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
841 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
842 {
843 G4cout<<"*G4QLE::PS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
844 }
845 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
846 aFraTrack = new G4Track(pAl, localtime, position );
847 aFraTrack->SetWeight(weight); // weighted
848 aFraTrack->SetTouchableHandle(trTouchable);
849#ifdef pdebug
850 G4cout<<"G4QLowEn::PSDI:-->PrRQFA4M="<<f4M<<G4endl;
851#endif
852 ++nSec;
853 reco4M=s4M;
854 }
855 else if(pA==5 && (pC==2 || pC==3)) // He5/Li5 decay
856 {
857 theDefinition = aProton;
858 G4double mNuc = mProt;
859 if(pC==2)
860 {
861 theDefinition = aNeutron;
862 mNuc = mNeut;
863 }
864 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
865 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
866 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
867 {
868 G4cout<<"*G4QLE::PS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
869 }
870 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
871 aFraTrack = new G4Track(pAl, localtime, position );
872 aFraTrack->SetWeight(weight); // weighted
873 aFraTrack->SetTouchableHandle(trTouchable);
874#ifdef pdebug
875 G4cout<<"G4QLowEn::PSDI:-->PrQFRN4M="<<f4M<<G4endl;
876#endif
877 ++nSec;
878 theDefinition = anAlpha;
879 reco4M=s4M;
880 }
881 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(pZ,pB,0,pZ);
882 ++nSec;
883#ifdef pdebug
884 G4cout<<"G4QLowEn::PSDI:-->TQF_nSec="<<nSec<<G4endl;
885#endif
886 aParticleChange.SetNumberOfSecondaries(nSec);
887 if(targSpect) aParticleChange.AddSecondary(targSpect);
888 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
889 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
890 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
891 G4Track* aResTrack = new G4Track(resA, localtime, position );
892 aResTrack->SetWeight(weight); // weighted
893 aResTrack->SetTouchableHandle(trTouchable);
894#ifdef pdebug
895 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
896#endif
897 aParticleChange.AddSecondary( aResTrack );
898 }
899#ifdef debug
900 G4cout<<"G4QLowEnergy::PostStepDoIt:***PostStepDoIt is done:Quasi-El***"<<G4endl;
901#endif
902 return G4VDiscreteProcess::PostStepDoIt(track, step); // ===> Reaction is DONE
903 }
904 else // The cental region compound can be created
905 {
906 // First calculate the isotopic state of the parts of the compound
907 if (!pZ) pC = 0;
908 else if(!pN) pC = pD;
909 else
910 {
911#ifdef pdebug
912 G4cout<<"G4QLowEn::PSDI: pD="<<pD<<", pZ="<<pZ<<", pA="<<pA<<G4endl;
913#endif
914 G4double C=pD*pZ/pA;
915 pC=static_cast<G4int>(C);
916 if(G4UniformRand() < C-pC) ++pC;
917 }
918 if(!tZ) tC=0;
919 else if(!tN) tC=tD;
920 else
921 {
922#ifdef pdebug
923 G4cout<<"G4QLowEn::PSDI: tD="<<tD<<", tZ="<<tZ<<", tA="<<tA<<G4endl;
924#endif
925 G4double C=tD*tZ/tA;
926 tC=static_cast<G4int>(C);
927 if(G4UniformRand() < C-tC) ++tC;
928 }
929 // calculate the transferred momentum
930 G4ThreeVector pFM(0.,0.,0.);
931 if(pD < pB) // The projectile nucleus must be splitted
932 {
933 G4int nc=pD;
934 if(pD+pD>pB) nc=pB-pD;
935 pFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
936 for(i=1; i < nc; ++i)
937 pFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
938 }
939 G4ThreeVector tFM(0.,0.,0.);
940 if(tD<tB) // The projectile nucleus must be splitted
941 {
942 G4int nc=pD;
943 if(tD+tD>tB) nc=tB-tD;
944 tFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
945 for(i=1; i < nc; ++i)
946 tFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
947 }
948#ifdef pdebug
949 G4cout<<"G4QLE::PSDI:pC="<<pC<<", tC="<<tC<<", pFM="<<pFM<<", tFM="<<tFM<<G4endl;
950#endif
951 // Split the projectile spectator
952 G4int rpZ=pZ-pC; // number of protons in the projectile spectator
953 G4int pF_value=pD-pC; // number of neutrons in the projectile part of comp
954 G4int rpN=pN-pF_value; // number of neutrons in the projectile spectator
955 G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0); // Mass of the spectator
956 G4ThreeVector boostV=proj4M.vect()/proj4M.e(); // Antilab Boost Vector
957 G4double rpE=std::sqrt(rpNM*rpNM+pFM.mag2());
958 G4LorentzVector rp4M(pFM,rpE);
959#ifdef pdebug
960 G4cout<<"G4QLE::PSDI: boostV="<<boostV<<",rp4M="<<rp4M<<",pr4M="<<proj4M<<G4endl;
961#endif
962 rp4M.boost(boostV);
963#ifdef pdebug
964 G4cout<<"G4QLE::PSDI: After boost, rp4M="<<rp4M<<G4endl;
965#endif
966 G4ParticleDefinition* theDefinition; // Prototype of projSpectatorNuclDefinition
967 G4int rpA=rpZ+rpN;
968 G4Track* aFraPTrack = 0;
969 theDefinition = 0;
970 if(rpA==1)
971 {
972 if(!rpZ) theDefinition = G4Neutron::Neutron();
973 else theDefinition = G4Proton::Proton();
974#ifdef pdebug
975 G4cout<<"G4QLE::PSDI: rpA=1, rpZ"<<rpZ<<G4endl;
976#endif
977 }
978 else if(rpA==2 && rpZ==0) // nn decay
979 {
980 theDefinition = aNeutron;
981 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
982 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
983#ifdef pdebug
984 G4cout<<"G4QLE::CPS->n+n,nn="<<rp4M.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
985#endif
986 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
987 {
988 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rp4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
989 }
990 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
991 aFraPTrack = new G4Track(pNeu, localtime, position );
992 aFraPTrack->SetWeight(weight); // weighted
993 aFraPTrack->SetTouchableHandle(trTouchable);
994 tt4M-=f4M;
995#ifdef edebug
996 totBaN-=2;
997 tch4M -=f4M;
998 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
999#endif
1000#ifdef pdebug
1001 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1002#endif
1003 rp4M=s4M;
1004 }
1005 else if(rpA>2 && rpZ==0) // Z=0 decay
1006 {
1007 theDefinition = aNeutron;
1008 G4LorentzVector f4M=rp4M/rpA; // 4mom of 1st neutron
1009#ifdef pdebug
1010 G4cout<<"G4QLE::CPS->Nn,M="<<rp4M.m()<<" >? N*MNeutron="<<rpA*mNeutron<<G4endl;
1011#endif
1012 for(G4int it=1; it<rpA; ++it) // Fill (N-1) neutrons to output
1013 {
1014 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1015 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
1016 aNTrack->SetWeight(weight); // weighted
1017 aNTrack->SetTouchableHandle(trTouchable);
1018 result.push_back(aNTrack);
1019 }
1020 G4int nesc = rpA-1;
1021 tt4M-=f4M*nesc;
1022#ifdef edebug
1023 totBaN-=nesc;
1024 tch4M -=f4M*nesc;
1025 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1026#endif
1027#ifdef pdebug
1028 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1029#endif
1030 rp4M=f4M;
1031 }
1032 else if(rpA==8 && rpZ==4) // Be8 decay
1033 {
1034 theDefinition = anAlpha;
1035 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
1036 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
1037#ifdef pdebug
1038 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1039#endif
1040 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
1041 {
1042 G4cout<<"*W*G4QLE::CPS->A+A,t="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1043 }
1044 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1045 aFraPTrack = new G4Track(pAl, localtime, position );
1046 aFraPTrack->SetWeight(weight); // weighted
1047 aFraPTrack->SetTouchableHandle(trTouchable);
1048 tt4M-=f4M;
1049#ifdef edebug
1050 totChg-=2;
1051 totBaN-=4;
1052 tch4M -=f4M;
1053 G4cout<<">>G4QLEn::PSDI:1,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1054#endif
1055#ifdef pdebug
1056 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1057#endif
1058 rp4M=s4M;
1059 }
1060 else if(rpA==5 && (rpZ==2 || rpZ==3)) // He5/Li5 decay
1061 {
1062 theDefinition = aProton;
1063 G4double mNuc = mProt;
1064 if(rpZ==2)
1065 {
1066 theDefinition = aNeutron;
1067 mNuc = mNeut;
1068 }
1069 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
1070 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
1071#ifdef pdebug
1072 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1073#endif
1074 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
1075 {
1076 G4cout<<"*W*G4QLE::CPS->N+A,t="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1077 }
1078 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1079 aFraPTrack = new G4Track(pAl, localtime, position );
1080 aFraPTrack->SetWeight(weight); // weighted
1081 aFraPTrack->SetTouchableHandle(trTouchable);
1082 tt4M-=f4M;
1083#ifdef edebug
1084 if(theDefinition == aProton) totChg-=1;
1085 totBaN-=1;
1086 tch4M -=f4M;
1087 G4cout<<">>G4QLEn::PSDI:2,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1088#endif
1089#ifdef pdebug
1090 G4cout<<"G4QLowEn::PSDI:-->ProjSpectN4M="<<f4M<<G4endl;
1091#endif
1092 theDefinition = anAlpha;
1093 rp4M=s4M;
1094 }
1095 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rpZ,rpA,0,rpZ);
1096 if(!theDefinition)
1097 {
1098 // G4cout<<"-Warning-G4QLowEn::PSDI: pDef=0, Z="<<rpZ<<", A="<<rpA<<G4endl;
1099 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
1100 G4ExceptionDescription ed;
1101 ed << "Particle definition is a null pointer: pDef=0, Z= " << rpZ
1102 << ", A=" << rpA << G4endl;
1103 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1104 JustWarning, ed);
1105 }
1106#ifdef edebug
1107 if(theDefinition == anAlpha)
1108 {
1109 totChg-=2;
1110 totBaN-=4;
1111 }
1112 else
1113 {
1114 totChg-=rpZ;
1115 totBaN-=rpA;
1116 }
1117 tch4M -=rp4M;
1118 G4cout<<">>G4QLEn::PSDI:3, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
1119#endif
1120 G4DynamicParticle* rpD = new G4DynamicParticle(theDefinition, rp4M);
1121 G4Track* aNewPTrack = new G4Track(rpD, localtime, position);
1122 aNewPTrack->SetWeight(weight);// weighted
1123 aNewPTrack->SetTouchableHandle(trTouchable);
1124 tt4M-=rp4M;
1125#ifdef pdebug
1126 G4cout<<"G4QLowEn::PSDI:-->ProjSpectR4M="<<rp4M<<",Z="<<rpZ<<",A="<<rpA<<G4endl;
1127#endif
1128 //
1129 // Split the target spectator
1130 G4int rtZ=tZ-tC; // number of protons in the target spectator
1131 G4int tF_value=tD-tC; // number of neutrons in the target part of comp
1132 G4int rtN=tN-tF_value; // number of neutrons in the target spectator
1133 G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0); // Mass of the spectator
1134 G4double rtE=std::sqrt(rtNM*rtNM+tFM.mag2());
1135 G4LorentzVector rt4M(tFM,rtE);
1136 G4int rtA=rtZ+rtN;
1137 G4Track* aFraTTrack = 0;
1138 theDefinition = 0;
1139 if(rtA==1)
1140 {
1141 if(!rtZ) theDefinition = G4Neutron::Neutron();
1142 else theDefinition = G4Proton::Proton();
1143#ifdef pdebug
1144 G4cout<<"G4QLE::PSDI: rtA=1, rtZ"<<rtZ<<G4endl;
1145#endif
1146 }
1147 else if(rtA==2 && rtZ==0) // nn decay
1148 {
1149 theDefinition = aNeutron;
1150 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
1151 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
1152#ifdef pdebug
1153 G4cout<<"G4QLE::CPS->n+n,nn="<<rptM.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
1154#endif
1155 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1156 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rt4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
1157 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1158 aFraPTrack = new G4Track(pNeu, localtime, position );
1159 aFraPTrack->SetWeight(weight); // weighted
1160 aFraPTrack->SetTouchableHandle(trTouchable);
1161 tt4M-=f4M;
1162#ifdef edebug
1163 totBaN-=2;
1164 tch4M -=f4M;
1165 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1166#endif
1167#ifdef pdebug
1168 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1169#endif
1170 rt4M=s4M;
1171 }
1172 else if(rtA>2 && rtZ==0) // Z=0 decay
1173 {
1174 theDefinition = aNeutron;
1175 G4LorentzVector f4M=rt4M/rtA; // 4mom of 1st neutron
1176#ifdef pdebug
1177 G4cout<<"G4QLE::CPS->Nn,M="<<rt4M.m()<<" >? N*MNeutron="<<rtA*mNeutron<<G4endl;
1178#endif
1179 for(G4int it=1; it<rtA; ++it) // Fill (N-1) neutrons to output
1180 {
1181 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1182 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
1183 aNTrack->SetWeight(weight); // weighted
1184 aNTrack->SetTouchableHandle(trTouchable);
1185 result.push_back(aNTrack);
1186 }
1187 G4int nesc = rtA-1;
1188 tt4M-=f4M*nesc;
1189#ifdef edebug
1190 totBaN-=nesc;
1191 tch4M -=f4M*nesc;
1192 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1193#endif
1194#ifdef pdebug
1195 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1196#endif
1197 rt4M=f4M;
1198 }
1199 else if(rtA==8 && rtZ==4) // Be8 decay
1200 {
1201 theDefinition = anAlpha;
1202 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
1203 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
1204#ifdef pdebug
1205 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1206#endif
1207 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1208 {
1209 G4cout<<"*W*G4QLE::CTS->A+A,t="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1210 }
1211 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1212 aFraTTrack = new G4Track(pAl, localtime, position );
1213 aFraTTrack->SetWeight(weight); // weighted
1214 aFraTTrack->SetTouchableHandle(trTouchable);
1215 tt4M-=f4M;
1216#ifdef edebug
1217 totChg-=2;
1218 totBaN-=4;
1219 tch4M -=f4M;
1220 G4cout<<">>G4QLEn::PSDI:4,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1221#endif
1222#ifdef pdebug
1223 G4cout<<"G4QLowEn::PSDI:-->TargSpectA4M="<<f4M<<G4endl;
1224#endif
1225 rt4M=s4M;
1226 }
1227 else if(rtA==5 && (rtZ==2 || rtZ==3)) // He5/Li5 decay
1228 {
1229 theDefinition = aProton;
1230 G4double mNuc = mProt;
1231 if(rpZ==2)
1232 {
1233 theDefinition = aNeutron;
1234 mNuc = mNeut;
1235 }
1236 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
1237 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
1238#ifdef pdebug
1239 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1240#endif
1241 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1242 {
1243 G4cout<<"*W*G4QLE::CTS->N+A,t="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1244 }
1245 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1246 aFraTTrack = new G4Track(pAl, localtime, position );
1247 aFraTTrack->SetWeight(weight); // weighted
1248 aFraTTrack->SetTouchableHandle(trTouchable);
1249 tt4M-=f4M;
1250#ifdef edebug
1251 if(theDefinition == aProton) totChg-=1;
1252 totBaN-=1;
1253 tch4M -=f4M;
1254 G4cout<<">>G4QLEn::PSDI:5,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1255#endif
1256#ifdef pdebug
1257 G4cout<<"G4QLowEn::PSDI:-->TargSpectN4M="<<f4M<<G4endl;
1258#endif
1259 theDefinition = anAlpha;
1260 rt4M=s4M;
1261 }
1262 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rtZ,rtA,0,rtZ);
1263 if(!theDefinition)
1264 {
1265 // G4cout<<"-Warning-G4QLowEn::PSDI: tDef=0, Z="<<rtZ<<", A="<<rtA<<G4endl;
1266 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
1267 G4ExceptionDescription ed;
1268 ed << "Particle definition is a null pointer: tDef=0, Z= " << rtZ
1269 << ", A=" << rtA << G4endl;
1270 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1271 JustWarning, ed);
1272 }
1273#ifdef edebug
1274 if(theDefinition == anAlpha)
1275 {
1276 totChg-=2;
1277 totBaN-=4;
1278 }
1279 else
1280 {
1281 totChg-=rtZ;
1282 totBaN-=rtA;
1283 }
1284 tch4M -=rt4M;
1285 G4cout<<">>G4QLEn::PSDI:6, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
1286#endif
1287 G4DynamicParticle* rtD = new G4DynamicParticle(theDefinition, rt4M);
1288 G4Track* aNewTTrack = new G4Track(rtD, localtime, position);
1289 aNewTTrack->SetWeight(weight); // weighted
1290 aNewTTrack->SetTouchableHandle(trTouchable);
1291 tt4M-=rt4M;
1292#ifdef pdebug
1293 G4cout<<"G4QLowEn::PSDI:-->TargSpectR4M="<<rt4M<<",Z="<<rtZ<<",A="<<rtA<<G4endl;
1294#endif
1295 if(aFraPTrack) result.push_back(aFraPTrack);
1296 if(aNewPTrack) result.push_back(aNewPTrack);
1297 if(aFraTTrack) result.push_back(aFraTTrack);
1298 if(aNewTTrack) result.push_back(aNewTTrack);
1299 tcM = tt4M.m(); // total CMS mass of the compound (after reduction!)
1300 G4int sN=tF_value+pF_value;
1301 G4int sZ=tC+pC;
1302 tnM = targQPDG.GetNuclMass(sZ,sN,0); // GSM
1303#ifdef pdebug
1304 G4cout<<"G4QLEn::PSDI:At#"<<nAt<<",totM="<<tcM<<",gsM="<<tnM<<",Z="<<sZ<<",N="
1305 <<sN<<G4endl;
1306#endif
1307 if(tcM > tnM) // !! The totalresidual reaction is modified !
1308 {
1309 pZ=pC;
1310 pN=pF_value;
1311 tZ=tC;
1312 tN=tF_value;
1313 tot4M=tt4M;
1314 }
1315 //else
1316 //{
1317 // G4cout<<"***G4QLowEnergy::PostStepDoIt: totM="<<tcM<<" <= GSM="<<tnM<<G4endl;
1318 // throw G4QException("G4QLowEnergy::PostStepDoIt: ResibualNucl below GSM shell");
1319 //}
1320 }
1321 } // At least one is splitted
1322 else if(tD==tB || pD==pB) // Total absorption
1323 {
1324 tD=tZ+tN;
1325 pD=pZ+pN;
1326 tcM=tnM+1.;
1327 }
1328 }
1329 else // Total compound (define tD to go out of the while)
1330 {
1331 tD=tZ+tN;
1332 pD=pZ+pN;
1333 tcM=tnM+1.;
1334 }
1335 } // End of the interaction WHILE
1336 G4double totM=tot4M.m(); // total CMS mass of the reaction
1337 G4int totN=tN+pN;
1338 G4int totZ=tZ+pZ;
1339#ifdef edebug
1340 G4cout<<">>>G4QLEn::PSDI: dZ="<<totChg-totZ<<", dB="<<totBaN-totN-totZ<<", d4M="
1341 <<tch4M-tot4M<<",N="<<totN<<",Z="<<totZ<<G4endl;
1342#endif
1343 // @@ Here mass[i] can be calculated if mass=0
1344 G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1345 0.,0.,0.,0.,0.,0.};
1346 mass[0] = tM = targQPDG.GetNuclMass(totZ, totN, 0); // gamma+gamma and GSM update
1347#ifdef pdebug
1348 G4cout<<"G4QLEn::PSDI:TotM="<<totM<<",NucM="<<tM<<",totZ="<<totZ<<",totN="<<totN<<G4endl;
1349#endif
1350 if (totN>0 && totZ>0)
1351 {
1352 mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron
1353 mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton
1354 }
1355 if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d
1356 if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t
1357 if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3
1358 if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a
1359 if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ ,totN-2,0); //n+n
1360 mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron)
1361 if ( totZ > 2 ) mass[9] = targQPDG.GetNuclMass(totZ-2,totN ,0); //p+p
1362 mass[10] = mass[5]; // proton+deuteron (the same as He3)
1363 mass[11] = mass[4]; // neutron+deuteron (the same as t)
1364 mass[12] = mass[6]; // deuteron+deuteron (the same as alpha)
1365 mass[13] = mass[6]; // proton+tritium (the same as alpha)
1366 if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);//n+t
1367 if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);//He3+p
1368 mass[16] = mass[6]; // neutron+He3 (the same as alpha)
1369 if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);//pa
1370 if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);//na
1371 if(pL>0) // @@ Not debugged @@
1372 {
1373 G4int pL1=pL-1;
1374 if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ ,totN ,pL1);// Lambda+gamma
1375 if( (totN > 0 && totZ > 0) || totZ > 1 )
1376 mass[20]=targQPDG.GetNuclMass(totZ-1,totN ,pL1);//Lp
1377 if( (totN > 0 && totZ > 0) || totN > 0 )
1378 mass[21]=targQPDG.GetNuclMass(totZ ,totN-1,pL1);//Ln
1379 if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) )
1380 mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld
1381 if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) )
1382 mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt
1383 if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) )
1384 mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3
1385 if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) )
1386 mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La
1387 }
1388#ifdef debug
1389 G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl;
1390#endif
1391 tA=tZ+tN;
1392 pA=pZ+pN;
1393 G4double prZ=pZ/pA+tZ/tA;
1394 G4double prN=pN/pA+tN/tA;
1395 G4double prD=prN*prZ;
1396 G4double prA=prD*prD;
1397 G4double prH=prD*prZ;
1398 G4double prT=prD*prN;
1399 G4double fhe3=6.*std::sqrt(tA);
1400 G4double prL=0.;
1401 if(pL>0) prL=pL/pA;
1402 G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1403 0.,0.,0.,0.,0.,0.};
1404 qval[ 0] = (totM - mass[ 0])/137./137.;
1405 qval[ 1] = (totM - mass[ 1] - mNeut)*prN/137.;
1406 qval[ 2] = (totM - mass[ 2] - mProt)*prZ/137.;
1407 qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3./137.;
1408 qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6./137.;
1409 qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3/137.;
1410 qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9./137.;
1411 qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN;
1412 qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD;
1413 qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ;
1414 qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.;
1415 qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.;
1416 qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.;
1417 qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.;
1418 qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.;
1419 qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3;
1420 qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3;
1421 qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.;
1422 qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.;
1423 if(pZ>0)
1424 {
1425 qval[19] = (totM - mass[19] - mLamb)*prL;
1426 qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ;
1427 qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN;
1428 qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.;
1429 qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.;
1430 qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3;
1431 qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4;
1432 }
1433#ifdef debug
1434 G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<", prA="<<pA<<G4endl;
1435#endif
1436
1437 G4double qv = 0.0; // Total sum of probabilities (q-values)
1438 for(i=0; i<nCh; ++i )
1439 {
1440#ifdef sdebug
1441 G4cout<<"G4QLowEn::PSDI: i="<<i<<", q="<<qval[i]<<",m="<<mass[i]<<G4endl;
1442#endif
1443 if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels (mesons)
1444 if( qval[i] < 0.0 ) qval[i] = 0.0; // Close the splitting impossible channels
1445 qv += qval[i];
1446 }
1447 // Select the channel
1448 G4double qv1 = 0.0;
1449 G4double ran = qv*G4UniformRand();
1450 G4int index = 0;
1451 for( index=0; index<nCh; ++index ) if( qval[index] > 0.0 )
1452 {
1453 qv1 += qval[index];
1454 if( ran <= qv1 ) break;
1455 }
1456#ifdef debug
1457 G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl;
1458#endif
1459 if(index == nCh)
1460 {
1461 G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl;
1462 G4cout<<"G4QLowEn::PoStDI:Reaction "<<projPDG<<"+"<<targPDG<<", P="<<Momentum<<G4endl;
1463 G4int nRes=result.size();
1464 if(nRes)G4cout<<"G4QLowEnergy::PoStDI: result exists with N="<<nRes<<" tracks"<<G4endl;
1465// else throw G4QException("***G4QLowEnergy::PostStepDoIt: Can't decay ResidualCompound");
1466 else {
1467 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1468 FatalException, "Can't decay ResidualCompound");
1469 }
1470 G4double minEx=1000000.; // Prototype of minimal excess
1471 G4bool found = false; // The solution is not found yet
1472 G4int cInd = 0; // Prototype of the best index
1473 G4int ctN = 0; // To
1474 G4int ctZ = 0; // avoid
1475 G4LorentzVector c4M=tot4M; // recalculation
1476 G4double ctM=0.; // of found
1477 for(G4int it=0; it<nRes; ++it)
1478 {
1479 G4Track* iTrack = result[it];
1480 const G4DynamicParticle* iHadron = iTrack->GetDynamicParticle();
1481 G4ParticleDefinition* iParticle=iHadron->GetDefinition();
1482 G4int iPDG = iParticle->GetPDGEncoding();
1483 G4LorentzVector ih4M = projHadron->Get4Momentum();
1484 G4cout<<" G4QLowEn::PoStDI: H["<<it<<"] = [ "<<iPDG<<", "<<ih4M<<" ]"<<G4endl;
1485 G4int iB = iParticle->GetBaryonNumber(); // A of the secondary
1486 G4int iC = G4int(iParticle->GetPDGCharge()); // Z of the secondary
1487 G4LorentzVector new4M = tot4M + ih4M; // Make temporary tot4M
1488 G4int ntN=totN + iB-iC; // Make temporary totN
1489 G4int ntZ=totZ + iC; // Make temporary totZ
1490 G4double ntM = targQPDG.GetNuclMass(ntZ, ntN, 0); // Make temporary GSM
1491 G4double ttM = new4M.m();
1492 if(ttM > ntM) // >>> This is at least one possible solution ! <<<
1493 {
1494 G4double nEx = ttM - ntM;
1495 if(found && nEx < minEx) // This one is better than the previous found
1496 {
1497 cInd = it;
1498 minEx= nEx;
1499 ctN = ntN;
1500 ctZ = ntZ;
1501 c4M = new4M;
1502 ctM = ttM;
1503 mass[0] = ntM;
1504 }
1505 found = true;
1506 }
1507 }
1508 if(found)
1509 {
1510 tot4M = c4M;
1511 totM = ctM;
1512 totN = ctN;
1513 totZ = ctZ;
1514 delete result[cInd];
1515 G4int nR1 = nRes -1;
1516 if(cInd < nR1) result[cInd] = result[nR1];
1517 result.pop_back();
1518 //nRes=nR1; // @@ can be used for the two pt correction
1519 index = 0; // @@ emergency gamma+gamma decay is selected
1520 }
1521 else
1522 {
1523 G4cout<<"***G4QLowEnergy::PoStDI: One hadron coddection did not succeed***"<<G4endl;
1524 if(nRes>1) G4cout<<"***G4QLowEnergy::PoStDI:nRes's big enough to use 2PtCor"<<G4endl;
1525 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't correct by one particle.");
1526 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1527 FatalException, "Can't correct by one particle");
1528 }
1529 }
1530 // @@ Convert it to G4QHadrons
1531 G4DynamicParticle* ResSec = new G4DynamicParticle;
1532 G4DynamicParticle* FstSec = new G4DynamicParticle;
1533 G4DynamicParticle* SecSec = new G4DynamicParticle;
1534#ifdef debug
1535 G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl;
1536#endif
1537
1538 G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype
1539 G4double mF=0.;
1540 G4double mS=0.;
1541 G4int rA=totZ+totN;
1542 G4int rZ=totZ;
1543 G4int rL=pL;
1544 G4int complete=3;
1545 G4ParticleDefinition* theDefinition; // Prototype for qfNucleon
1546 switch( index )
1547 {
1548 case 0:
1549 if(!evaporate || rA<2)
1550 {
1551 if(!rZ) theDefinition=aNeutron;
1552 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1553 if(!theDefinition)
1554 {
1555 // G4cerr<<"-Warning-G4LE::PSDI: notDef(0), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1556 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1557 G4ExceptionDescription ed;
1558 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
1559 << ", A=" << rA << ", L=" << rL << G4endl;
1560 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1561 JustWarning, ed);
1562 }
1563 ResSec->SetDefinition( theDefinition );
1564 FstSec->SetDefinition( aGamma );
1565 SecSec->SetDefinition( aGamma );
1566 }
1567 else
1568 {
1569 delete ResSec;
1570 delete FstSec;
1571 delete SecSec;
1572 complete=0;
1573 }
1574 break;
1575 case 1:
1576 rA-=1; // gamma+n
1577 if(!evaporate || rA<2)
1578 {
1579 if(!rZ) theDefinition=aNeutron;
1580 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1581 if(!theDefinition)
1582 {
1583 // G4cerr<<"-Warning-G4LE::PSDI: notDef(1), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1584 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1585 G4ExceptionDescription ed;
1586 ed << "Particle definition is a null pointer: notDef(1), Z= " << rZ
1587 << ", A=" << rA << ", L=" << rL << G4endl;
1588 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0001",
1589 JustWarning, ed);
1590 }
1591 ResSec->SetDefinition( theDefinition );
1592 SecSec->SetDefinition( aGamma );
1593 }
1594 else
1595 {
1596 delete ResSec;
1597 delete SecSec;
1598 complete=1;
1599 }
1600 FstSec->SetDefinition( aNeutron );
1601 mF=mNeut; // First hadron 4-momentum
1602 break;
1603 case 2:
1604 rA-=1;
1605 rZ-=1; // gamma+p
1606 if(!evaporate || rA<2)
1607 {
1608 if(!rZ) theDefinition=aNeutron;
1609 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1610 if(!theDefinition)
1611 {
1612 // G4cerr<<"-Warning-G4LE::PSDI: notDef(2), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1613 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1614 G4ExceptionDescription ed;
1615 ed << "Particle definition is a null pointer: notDef(2), Z= " << rZ
1616 << ", A=" << rA << ", L="<< rL << G4endl;
1617 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0002",
1618 JustWarning, ed);
1619 }
1620 ResSec->SetDefinition( theDefinition );
1621 SecSec->SetDefinition( aGamma );
1622 }
1623 else
1624 {
1625 delete ResSec;
1626 delete SecSec;
1627 complete=1;
1628 }
1629 FstSec->SetDefinition( aProton );
1630 mF=mProt; // First hadron 4-momentum
1631 break;
1632 case 3:
1633 rA-=2;
1634 rZ-=1; // gamma+d
1635 if(!evaporate || rA<2)
1636 {
1637 if(!rZ) theDefinition=aNeutron;
1638 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1639 if(!theDefinition)
1640 {
1641 // G4cerr<<"-Warning-G4LE::PSDI: notDef(3), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1642 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1643 G4ExceptionDescription ed;
1644 ed << "Particle definition is a null pointer: notDef(3), Z= " << rZ
1645 << ", A=" << rA << ", L=" << rL << G4endl;
1646 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0003",
1647 JustWarning, ed);
1648 }
1649 ResSec->SetDefinition( theDefinition );
1650 SecSec->SetDefinition( aGamma );
1651 }
1652 else
1653 {
1654 delete ResSec;
1655 delete SecSec;
1656 complete=1;
1657 }
1658 FstSec->SetDefinition( aDeuteron );
1659 mF=mDeut; // First hadron 4-momentum
1660 break;
1661 case 4:
1662 rA-=3; // gamma+t
1663 rZ-=1;
1664 if(!evaporate || rA<2)
1665 {
1666 if(!rZ) theDefinition=aNeutron;
1667 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1668 if(!theDefinition)
1669 {
1670 // G4cerr<<"-Warning-G4LE::PSDI: notDef(4), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1671 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1672 G4ExceptionDescription ed;
1673 ed << "Particle definition is a null pointer: notDef(4), Z= " << rZ
1674 << ", A=" << rA << ", L=" << rL << G4endl;
1675 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0004",
1676 JustWarning, ed);
1677 }
1678 ResSec->SetDefinition( theDefinition );
1679 SecSec->SetDefinition( aGamma );
1680 }
1681 else
1682 {
1683 delete ResSec;
1684 delete SecSec;
1685 complete=1;
1686 }
1687 FstSec->SetDefinition( aTriton );
1688 mF=mTrit; // First hadron 4-momentum
1689 break;
1690 case 5: // gamma+He3
1691 rA-=3;
1692 rZ-=2;
1693 if(!evaporate || rA<2)
1694 {
1695 if(!rZ) theDefinition=aNeutron;
1696 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1697 if(!theDefinition)
1698 {
1699 // G4cerr<<"-Warning-G4LE::PSDI: notDef(5), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1700 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1701 G4ExceptionDescription ed;
1702 ed << "Particle definition is a null pointer: notDef(5), Z= " << rZ
1703 << ", A=" << rA << ", L=" << rL << G4endl;
1704 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0005",
1705 JustWarning, ed);
1706 }
1707 ResSec->SetDefinition( theDefinition );
1708 SecSec->SetDefinition( aGamma );
1709 }
1710 else
1711 {
1712 delete ResSec;
1713 delete SecSec;
1714 complete=1;
1715 }
1716 FstSec->SetDefinition( aHe3);
1717 mF=mHel3; // First hadron 4-momentum
1718 break;
1719 case 6:
1720 rA-=4;
1721 rZ-=2; // gamma+He4
1722 if(!evaporate || rA<2)
1723 {
1724 if(!rZ) theDefinition=aNeutron;
1725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1726 if(!theDefinition)
1727 {
1728 // G4cerr<<"-Warning-G4LE::PSDI: notDef(6), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1729 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1730 G4ExceptionDescription ed;
1731 ed << "Particle definition is a null pointer: notDef(6), Z= " << rZ
1732 << ", A=" << rA << ", L=" << rL << G4endl;
1733 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0006",
1734 JustWarning, ed);
1735 }
1736 ResSec->SetDefinition( theDefinition );
1737 SecSec->SetDefinition( aGamma );
1738 }
1739 else
1740 {
1741 delete ResSec;
1742 delete SecSec;
1743 complete=1;
1744 }
1745 FstSec->SetDefinition( anAlpha );
1746 mF=mAlph; // First hadron 4-momentum
1747 break;
1748 case 7:
1749 rA-=2; // n+n
1750 if(rA==1 && !rZ) theDefinition=aNeutron;
1751 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1752 if(!theDefinition)
1753 {
1754 // G4cerr<<"-Warning-G4LE::PSDI: notDef(7), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1755 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1756 G4ExceptionDescription ed;
1757 ed << "Particle definition is a null pointer: notDef(7), Z= " << rZ
1758 << ", A=" << rA << ", L=" << rL << G4endl;
1759 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0007",
1760 JustWarning, ed);
1761 }
1762 ResSec->SetDefinition( theDefinition );
1763 FstSec->SetDefinition( aNeutron );
1764 SecSec->SetDefinition( aNeutron );
1765 mF=mNeut; // First hadron 4-momentum
1766 mS=mNeut; // Second hadron 4-momentum
1767 break;
1768 case 8:
1769 rZ-=1;
1770 rA-=2; // n+p
1771 if(rA==1 && !rZ) theDefinition=aNeutron;
1772 else if(rA==2 && !rZ)
1773 {
1774 index=7;
1775 ResSec->SetDefinition( aDeuteron);
1776 FstSec->SetDefinition( aNeutron );
1777 SecSec->SetDefinition( aNeutron );
1778 mF=mNeut; // First hadron 4-momentum
1779 mS=mNeut; // Second hadron 4-momentum
1780 break;
1781 }
1782 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1783 if(!theDefinition)
1784 {
1785 // G4cerr<<"-Warning-G4LE::PSDI: notDef(8), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1786 // throw G4QException("G4QLowEn::PostStepDoIt: Darticle Definition is a null pointer");
1787 G4ExceptionDescription ed;
1788 ed << "Particle definition is a null pointer: notDef(8), Z= " << rZ
1789 << ", A=" << rA << ", L=" << rL << G4endl;
1790 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0008",
1791 JustWarning, ed);
1792 }
1793 ResSec->SetDefinition( theDefinition );
1794 FstSec->SetDefinition( aNeutron );
1795 SecSec->SetDefinition( aProton );
1796 mF=mNeut; // First hadron 4-momentum
1797 mS=mProt; // Second hadron 4-momentum
1798 break;
1799 case 9:
1800 rZ-=2;
1801 rA-=2; // p+p
1802 if(rA==1 && !rZ) theDefinition=aNeutron;
1803 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1804 if(!theDefinition)
1805 {
1806 // G4cerr<<"-Warning-G4LE::PSDI: notDef(9), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1807 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1808 G4ExceptionDescription ed;
1809 ed << "Particle definition is a null pointer: notDef(9), Z= " << rZ
1810 << ", A=" << rA << ", L=" << rL << G4endl;
1811 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0009",
1812 JustWarning, ed);
1813 }
1814 ResSec->SetDefinition( theDefinition );
1815 FstSec->SetDefinition( aProton );
1816 SecSec->SetDefinition( aProton );
1817 mF=mProt; // First hadron 4-momentum
1818 mS=mProt; // Second hadron 4-momentum
1819 break;
1820 case 10:
1821 rZ-=2;
1822 rA-=3; // p+d
1823 if(rA==1 && !rZ) theDefinition=aNeutron;
1824 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1825 if(!theDefinition)
1826 {
1827 // G4cerr<<"-Warning-G4LE::PSDI: notDef(10), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1828 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1829 G4ExceptionDescription ed;
1830 ed << "Particle definition is a null pointer: notDef(10), Z= " << rZ
1831 << ", A=" << rA << ", L=" << rL << G4endl;
1832 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0010",
1833 JustWarning, ed);
1834 }
1835 ResSec->SetDefinition( theDefinition );
1836 FstSec->SetDefinition( aProton );
1837 SecSec->SetDefinition( aDeuteron );
1838 mF=mProt; // First hadron 4-momentum
1839 mS=mDeut; // Second hadron 4-momentum
1840 break;
1841 case 11:
1842 rZ-=1;
1843 rA-=3; // n+d
1844 if(rA==1 && !rZ) theDefinition=aNeutron;
1845 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1846 if(!theDefinition)
1847 {
1848 // G4cerr<<"-Warning-G4LE::PSDI: notDef(11), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1849 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1850 G4ExceptionDescription ed;
1851 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
1852 << ", A=" << rA << ", L=" << rL << G4endl;
1853 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0011",
1854 JustWarning, ed);
1855 }
1856 ResSec->SetDefinition( theDefinition );
1857 FstSec->SetDefinition( aNeutron );
1858 SecSec->SetDefinition( aDeuteron );
1859 mF=mNeut; // First hadron 4-momentum
1860 mS=mDeut; // Second hadron 4-momentum
1861 break;
1862 case 12:
1863 rZ-=2;
1864 rA-=4; // d+d
1865 if(rA==1 && !rZ) theDefinition=aNeutron;
1866 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1867 if(!theDefinition)
1868 {
1869 // G4cerr<<"-Warning-G4LE::PSDI: notDef(12), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1870 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1871 G4ExceptionDescription ed;
1872 ed << "Particle definition is a null pointer: notDef(12), Z= " << rZ
1873 << ", A=" << rA << ", L=" << rL << G4endl;
1874 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0012",
1875 JustWarning, ed);
1876 }
1877 ResSec->SetDefinition( theDefinition );
1878 FstSec->SetDefinition( aDeuteron );
1879 SecSec->SetDefinition( aDeuteron );
1880 mF=mDeut; // First hadron 4-momentum
1881 mS=mDeut; // Second hadron 4-momentum
1882 break;
1883 case 13:
1884 rZ-=2;
1885 rA-=4; // p+t
1886 if(rA==1 && !rZ) theDefinition=aNeutron;
1887 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1888 if(!theDefinition)
1889 {
1890 // G4cerr<<"-Warning-G4LE::PSDI: notDef(13), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1891 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1892 G4ExceptionDescription ed;
1893 ed << "Particle definition is a null pointer: notDef(13), Z= " << rZ
1894 << ", A=" << rA << ", L=" << rL << G4endl;
1895 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0013",
1896 JustWarning, ed);
1897 }
1898 ResSec->SetDefinition( theDefinition );
1899 FstSec->SetDefinition( aProton );
1900 SecSec->SetDefinition( aTriton );
1901 mF=mProt; // First hadron 4-momentum
1902 mS=mTrit; // Second hadron 4-momentum
1903 break;
1904 case 14:
1905 rZ-=1;
1906 rA-=4; // n+t
1907 if(rA==1 && !rZ) theDefinition=aNeutron;
1908 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1909 if(!theDefinition)
1910 {
1911 // G4cerr<<"-Warning-G4LE::PSDI: notDef(14), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1912 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1913 G4ExceptionDescription ed;
1914 ed << "Particle definition is a null pointer: notDef(14), Z= " << rZ
1915 << ", A=" << rA << ", L=" << rL << G4endl;
1916 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0014",
1917 JustWarning, ed);
1918 }
1919 ResSec->SetDefinition( theDefinition );
1920 FstSec->SetDefinition( aNeutron );
1921 SecSec->SetDefinition( aTriton );
1922 mF=mNeut; // First hadron 4-momentum
1923 mS=mTrit; // Second hadron 4-momentum
1924 break;
1925 case 15:
1926 rZ-=3;
1927 rA-=4; // p+He3
1928 if(rA==1 && !rZ) theDefinition=aNeutron;
1929 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1930 if(!theDefinition)
1931 {
1932 // G4cerr<<"-Warning-G4LE::PSDI: notDef(15), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1933 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1934 G4ExceptionDescription ed;
1935 ed << "Particle definition is a null pointer: notDef(15), Z= " << rZ
1936 << ", A=" << rA << ", L=" << rL << G4endl;
1937 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0015",
1938 JustWarning, ed);
1939 }
1940 ResSec->SetDefinition( theDefinition );
1941 FstSec->SetDefinition( aProton);
1942 SecSec->SetDefinition( aHe3 );
1943 mF=mProt; // First hadron 4-momentum
1944 mS=mHel3; // Second hadron 4-momentum
1945 break;
1946 case 16:
1947 rZ-=2;
1948 rA-=4; // n+He3
1949 if(rA==1 && !rZ) theDefinition=aNeutron;
1950 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1951 if(!theDefinition)
1952 {
1953 // G4cerr<<"-Warning-G4LE::PSDI: notDef(16), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1954 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1955 G4ExceptionDescription ed;
1956 ed << "Particle definition is a null pointer: notDef(16), Z= " << rZ
1957 << ", A=" << rA << ", L=" << rL << G4endl;
1958 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0016",
1959 JustWarning, ed);
1960 }
1961 ResSec->SetDefinition( theDefinition );
1962 FstSec->SetDefinition( aNeutron );
1963 SecSec->SetDefinition( aHe3 );
1964 mF=mNeut; // First hadron 4-momentum
1965 mS=mHel3; // Second hadron 4-momentum
1966 break;
1967 case 17:
1968 rZ-=3;
1969 rA-=5; // p+alph
1970 if(rA==1 && !rZ) theDefinition=aNeutron;
1971 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1972 if(!theDefinition)
1973 {
1974 // G4cerr<<"-Warning-G4LE::PSDI: notDef(17), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1975 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1976 G4ExceptionDescription ed;
1977 ed << "Particle definition is a null pointer: notDef(17), Z= " << rZ
1978 << ", A=" << rA << ", L=" << rL << G4endl;
1979 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0017",
1980 JustWarning, ed);
1981 }
1982 ResSec->SetDefinition( theDefinition );
1983 FstSec->SetDefinition( aProton );
1984 SecSec->SetDefinition( anAlpha );
1985 mF=mProt; // First hadron 4-momentum
1986 mS=mAlph; // Second hadron 4-momentum
1987 break;
1988 case 18:
1989 rZ-=2;
1990 rA-=5; // n+alph
1991 if(rA==1 && !rZ) theDefinition=aNeutron;
1992 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1993 if(!theDefinition)
1994 {
1995 // G4cerr<<"-Warning-G4LE::PSDI: notDef(18), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1996 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1997 G4ExceptionDescription ed;
1998 ed << "Particle definition is a null pointer: notDef(18), Z= " << rZ
1999 << ", A=" << rA << ", L=" << rL << G4endl;
2000 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0018",
2001 JustWarning, ed);
2002 }
2003 ResSec->SetDefinition( theDefinition );
2004 FstSec->SetDefinition( aNeutron );
2005 SecSec->SetDefinition( anAlpha );
2006 mF=mNeut; // First hadron 4-momentum
2007 mS=mAlph; // Second hadron 4-momentum
2008 break;
2009 case 19:
2010 rL-=1; // L+gamma (@@ rA inludes rL?)
2011 rA-=1;
2012 if(rA==1 && !rZ) theDefinition=aNeutron;
2013 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2014 if(!theDefinition)
2015 {
2016 // G4cerr<<"-Warning-G4LE::PSDI: notDef(19), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2017 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2018 G4ExceptionDescription ed;
2019 ed << "Particle definition is a null pointer: notDef(19), Z= " << rZ
2020 << ", A=" << rA << ", L=" << rL << G4endl;
2021 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0019",
2022 JustWarning, ed);
2023 }
2024 ResSec->SetDefinition( theDefinition );
2025 FstSec->SetDefinition( aLambda );
2026 SecSec->SetDefinition( aGamma );
2027 mF=mLamb; // First hadron 4-momentum
2028 break;
2029 case 20:
2030 rL-=1; // L+p (@@ rA inludes rL?)
2031 rZ-=1;
2032 rA-=2;
2033 if(rA==1 && !rZ) theDefinition=aNeutron;
2034 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2035 if(!theDefinition)
2036 {
2037 // G4cerr<<"-Warning-G4LE::PSDI: notDef(20), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2038 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2039 G4ExceptionDescription ed;
2040 ed << "Particle definition is a null pointer: notDef(20), Z= " << rZ
2041 << ", A=" << rA << ", L=" << rL << G4endl;
2042 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0020",
2043 JustWarning, ed);
2044 }
2045 ResSec->SetDefinition( theDefinition );
2046 FstSec->SetDefinition( aProton );
2047 SecSec->SetDefinition( aLambda );
2048 mF=mProt; // First hadron 4-momentum
2049 mS=mLamb; // Second hadron 4-momentum
2050 break;
2051 case 21:
2052 rL-=1; // L+n (@@ rA inludes rL?)
2053 rA-=2;
2054 if(rA==1 && !rZ) theDefinition=aNeutron;
2055 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2056 if(!theDefinition)
2057 {
2058 // G4cerr<<"-Warning-G4LE::PSDI: notDef(21), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2059 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Pefinition is a null pointer");
2060 G4ExceptionDescription ed;
2061 ed << "Particle definition is a null pointer: notDef(21), Z= " << rZ
2062 << ", A=" << rA << ", L=" << rL << G4endl;
2063 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0021",
2064 JustWarning, ed);
2065 }
2066 ResSec->SetDefinition( theDefinition );
2067 FstSec->SetDefinition( aNeutron );
2068 SecSec->SetDefinition( aLambda );
2069 mF=mNeut; // First hadron 4-momentum
2070 mS=mLamb; // Second hadron 4-momentum
2071 break;
2072 case 22:
2073 rL-=1; // L+d (@@ rA inludes rL?)
2074 rZ-=1;
2075 rA-=3;
2076 if(rA==1 && !rZ) theDefinition=aNeutron;
2077 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2078 if(!theDefinition)
2079 {
2080 // G4cerr<<"-Warning-G4LE::PSDI: notDef(22), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2081 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2082 G4ExceptionDescription ed;
2083 ed << "Particle definition is a null pointer: notDef(22), Z= " << rZ
2084 << ", A=" << rA << ", L=" << rL << G4endl;
2085 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0022",
2086 JustWarning, ed);
2087 }
2088 ResSec->SetDefinition( theDefinition );
2089 FstSec->SetDefinition( aDeuteron );
2090 SecSec->SetDefinition( aLambda );
2091 mF=mDeut; // First hadron 4-momentum
2092 mS=mLamb; // Second hadron 4-momentum
2093 break;
2094 case 23:
2095 rL-=1; // L+t (@@ rA inludes rL?)
2096 rZ-=1;
2097 rA-=4;
2098 if(rA==1 && !rZ) theDefinition=aNeutron;
2099 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2100 if(!theDefinition)
2101 {
2102 // G4cerr<<"-Warning-G4LE::PSDI: notDef(23), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2103 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2104 G4ExceptionDescription ed;
2105 ed << "Particle definition is a null pointer: notDef(23), Z= " << rZ
2106 << ", A=" << rA << ", L=" << rL << G4endl;
2107 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0023",
2108 JustWarning, ed);
2109 }
2110 ResSec->SetDefinition( theDefinition );
2111 FstSec->SetDefinition( aTriton );
2112 SecSec->SetDefinition( aLambda );
2113 mF=mTrit; // First hadron 4-momentum
2114 mS=mLamb; // Second hadron 4-momentum
2115 break;
2116 case 24:
2117 rL-=1; // L+He3 (@@ rA inludes rL?)
2118 rZ-=2;
2119 rA-=4;
2120 if(rA==1 && !rZ) theDefinition=aNeutron;
2121 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2122 if(!theDefinition)
2123 {
2124 // G4cerr<<"-Warning-G4LE::PSDI: notDef(24), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2125 // throw G4QException("G4QLowEn::PostStepDoIt: particle definition is a null pointer");
2126 G4ExceptionDescription ed;
2127 ed << "Particle definition is a null pointer: notDef(24), Z= " << rZ
2128 << ", A=" << rA << ", L=" << rL << G4endl;
2129 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0024",
2130 JustWarning, ed);
2131 }
2132 ResSec->SetDefinition( theDefinition );
2133 FstSec->SetDefinition( aHe3 );
2134 SecSec->SetDefinition( aLambda );
2135 mF=mHel3; // First hadron 4-momentum
2136 mS=mLamb; // Second hadron 4-momentum
2137 break;
2138 case 25:
2139 rL-=1; // L+alph (@@ rA inludes rL?)
2140 rZ-=2;
2141 rA-=5;
2142 if(rA==1 && !rZ) theDefinition=aNeutron;
2143 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2144 if(!theDefinition)
2145 {
2146 // G4cerr<<"-Warning-G4LE::PSDI: notDef(25), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2147 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2148 G4ExceptionDescription ed;
2149 ed << "Particle definition is a null pointer: notDef(25), Z= " << rZ
2150 << ", A=" << rA << ", L=" << rL << G4endl;
2151 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0025",
2152 JustWarning, ed);
2153 }
2154 ResSec->SetDefinition( theDefinition );
2155 FstSec->SetDefinition( anAlpha );
2156 SecSec->SetDefinition( aLambda );
2157 mF=mAlph; // First hadron 4-momentum
2158 mS=mLamb; // Second hadron 4-momentum
2159 break;
2160 }
2161#ifdef debug
2162 G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl;
2163#endif
2164 G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum
2165 G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum
2166 G4LorentzVector dir4Mom=tot4M; // Prototype of the resN decay direction 4-momentum
2167 dir4Mom.setE(tot4M.e()/2.); // Get half energy and total 3-momentum
2168 // @@ Can be repeated to take into account the Coulomb Barrier
2169 if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp))
2170 {
2171 // G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2="
2172 // <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl;
2173 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound");
2174 G4ExceptionDescription ed;
2175 ed << "Can't decay the Compound: i=" << index << ",tM=" << totM << "->M1="
2176 << res4Mom.m() << "+M2=" << fst4Mom.m() << "+M3=" << snd4Mom.m() << "=="
2177 << res4Mom.m()+fst4Mom.m()+snd4Mom.m() << G4endl;
2178 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
2179 FatalException, ed);
2180 }
2181#ifdef debug
2182 G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl;
2183#endif
2184 G4Track* aNewTrack = 0;
2185 if(complete)
2186 {
2187 FstSec->Set4Momentum(fst4Mom);
2188 aNewTrack = new G4Track(FstSec, localtime, position );
2189 aNewTrack->SetWeight(weight); // weighted
2190 aNewTrack->SetTouchableHandle(trTouchable);
2191 result.push_back( aNewTrack );
2192 EnMomConservation-=fst4Mom;
2193#ifdef debug
2194 G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom
2195 <<", PDG="<<FstSec->GetDefinition()->GetPDGEncoding()<<G4endl;
2196#endif
2197 if(complete>2) // Final solution
2198 {
2199 ResSec->Set4Momentum(res4Mom);
2200 aNewTrack = new G4Track(ResSec, localtime, position );
2201 aNewTrack->SetWeight(weight); // weighted
2202 aNewTrack->SetTouchableHandle(trTouchable);
2203 result.push_back( aNewTrack );
2204 EnMomConservation-=res4Mom;
2205#ifdef debug
2206 G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL="
2207 <<rL<<G4endl;
2208#endif
2209 SecSec->Set4Momentum(snd4Mom);
2210 aNewTrack = new G4Track(SecSec, localtime, position );
2211 aNewTrack->SetWeight(weight); // weighted
2212 aNewTrack->SetTouchableHandle(trTouchable);
2213 result.push_back( aNewTrack );
2214 EnMomConservation-=snd4Mom;
2215#ifdef debug
2216 G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom
2217 <<", PDG="<<SecSec->GetDefinition()->GetPDGEncoding()<<G4endl;
2218#endif
2219 }
2220 else res4Mom+=snd4Mom;
2221 }
2222 else res4Mom=tot4M;
2223 if(complete<3) // Evaporation of the residual must be done
2224 {
2225 G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus
2226 G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!)
2227 Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear !
2228 G4int nOut=evaHV->size();
2229 for(i=0; i<nOut; i++)
2230 {
2231 G4QHadron* curH = (*evaHV)[i];
2232 G4int hPDG=curH->GetPDGCode();
2233 G4LorentzVector h4Mom=curH->Get4Momentum();
2234 EnMomConservation-=h4Mom;
2235#ifdef debug
2236 G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl;
2237#endif
2238 if (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron;
2239 else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton;
2240 else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda;
2241 else if(hPDG== 22 ) theDefinition = aGamma;
2242 else if(hPDG== 111) theDefinition = aPiZero;
2243 else if(hPDG== 211) theDefinition = aPiPlus;
2244 else if(hPDG== -211) theDefinition = aPiMinus;
2245 else
2246 {
2247 G4int hZ=curH->GetCharge();
2248 G4int hA=curH->GetBaryonNumber();
2249 G4int hS=curH->GetStrangeness();
2250 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion
2251 }
2252 if(theDefinition)
2253 {
2254 G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom);
2255 G4Track* evaQH = new G4Track(theEQH, localtime, position );
2256 evaQH->SetWeight(weight); // weighted
2257 evaQH->SetTouchableHandle(trTouchable);
2258 result.push_back( evaQH );
2259 }
2260 else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl;
2261 }
2262 }
2263 // algorithm implementation --- STOPS HERE
2264 G4int nres=result.size();
2265 aParticleChange.SetNumberOfSecondaries(nres);
2266 for(i=0; i<nres; ++i) aParticleChange.AddSecondary(result[i]);
2267#ifdef debug
2268 G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
2269#endif
2270 return G4VDiscreteProcess::PostStepDoIt(track, step);
2271}
2272
2273G4double G4QLowEnergy::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG)
2274{
2275 static G4bool first=true;
2276 static G4VQCrossSection* CSmanager;
2277 if(first) // Connection with a singletone
2278 {
2279 CSmanager=G4QIonIonCrossSection::GetPointer();
2280 if(PDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer();
2281 first=false;
2282 }
2283#ifdef debug
2284 G4cout<<"G4QLowE::CXS: *DONE* p="<<p<<",Z="<<Z<<",N="<<N<<",PDG="<<PDG<<G4endl;
2285#endif
2286 return CSmanager->GetCrossSection(true, p, Z, N, PDG);
2287}
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4ForceCondition
G4ForceCondition
Definition: G4ForceCondition.hh:51
NotForced
@ NotForced
Definition: G4ForceCondition.hh:56
G4HadronicDeprecate
#define G4HadronicDeprecate(name)
Definition: G4HadronicDeprecate.hh:33
G4IsotopeVector
std::vector< G4Isotope * > G4IsotopeVector
Definition: G4IsotopeVector.hh:44
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
fHadronic
@ fHadronic
Definition: G4ProcessType.hh:49
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4RandomDirection
G4ThreeVector G4RandomDirection()
Definition: G4RandomDirection.hh:55
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:56
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:271
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:351
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4Element::GetRelativeAbundanceVector
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Element::GetIsotope
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:182
G4Element::GetIsotopeVector
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4GenericIon::GenericIon
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:92
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4Isotope::GetZ
G4int GetZ() const
Definition: G4Isotope.hh:91
G4Isotope::GetN
G4int GetN() const
Definition: G4Isotope.hh:94
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4ParticleChange::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChange.cc:228
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:136
G4ParticleDefinition::GetQuarkContent
G4int GetQuarkContent(G4int flavor) const
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:137
G4ParticleDefinition::GetAntiQuarkContent
G4int GetAntiQuarkContent(G4int flavor) const
G4ParticleTable::FindIon
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
Definition: G4ParticleTable.cc:399
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:63
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4QCHIPSWorld::Get
static G4QCHIPSWorld * Get()
Definition: G4QCHIPSWorld.cc:79
G4QCHIPSWorld::GetParticles
G4QParticleVector * GetParticles(G4int nOfParts=0)
Definition: G4QCHIPSWorld.cc:95
G4QHadron::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QHadron.hh:79
G4QHadron::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QHadron.hh:181
G4QHadron::GetCharge
G4int GetCharge() const
Definition: G4QHadron.hh:179
G4QHadron::GetPDGCode
G4int GetPDGCode() const
Definition: G4QHadron.hh:170
G4QHadron::GetStrangeness
G4int GetStrangeness() const
Definition: G4QHadron.hh:180
G4QIonIonCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QIonIonCrossSection.cc:65
G4QIsotope::GetCSVector
std::vector< std::pair< G4int, G4double > * > * GetCSVector(G4int Z, G4int index=0)
Definition: G4QIsotope.cc:1737
G4QIsotope::IsDefined
G4bool IsDefined(G4int Z, G4int Ind)
Definition: G4QIsotope.cc:1660
G4QIsotope::GetMeanCrossSection
G4double GetMeanCrossSection(G4int Z, G4int index=0)
Definition: G4QIsotope.cc:1842
G4QIsotope::Get
static G4QIsotope * Get()
Definition: G4QIsotope.cc:720
G4QIsotope::InitElement
G4int InitElement(G4int Z, G4int index, std::vector< std::pair< G4int, G4double > * > *abund)
Definition: G4QIsotope.cc:1557
G4QLowEnergy::GetEnegryMomentumConservation
G4LorentzVector GetEnegryMomentumConservation()
Definition: G4QLowEnergy.cc:76
G4QLowEnergy::GetMeanFreePath
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
Definition: G4QLowEnergy.cc:83
G4QLowEnergy::PostStepDoIt
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4QLowEnergy.cc:247
G4QLowEnergy::IsApplicable
G4bool IsApplicable(const G4ParticleDefinition &particle)
Definition: G4QLowEnergy.cc:228
G4QLowEnergy::GetNumberOfNeutronsInTarget
G4int GetNumberOfNeutronsInTarget()
Definition: G4QLowEnergy.cc:78
G4QLowEnergy::G4QLowEnergy
G4QLowEnergy(const G4String &processName="CHIPS_LowEnergyIonIonInelastic")
Definition: G4QLowEnergy.cc:60
G4QNucleus::CoulombBarGen
G4double CoulombBarGen(const G4double &rZ, const G4double &rA, const G4double &cZ, const G4double &cA)
Definition: G4QNucleus.cc:3358
G4QNucleus::EvaporateNucleus
void EvaporateNucleus(G4QHadron *hA, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:4171
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QPDGCode::GetNuclMass
G4double GetNuclMass(G4int Z, G4int N, G4int S)
Definition: G4QPDGCode.cc:766
G4QProtonElasticCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonElasticCrossSection.cc:143
G4QProtonNuclearCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonNuclearCrossSection.cc:63
G4Step
Definition: G4Step.hh:78
G4Track::GetWeight
G4double GetWeight() const
G4Track::SetWeight
void SetWeight(G4double aValue)
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Track::GetGlobalTime
G4double GetGlobalTime() const
G4Track::GetMaterial
G4Material * GetMaterial() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4VDiscreteProcess::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
Definition: G4VDiscreteProcess.cc:112
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)
G4VProcess::aParticleChange
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4VProcess::verboseLevel
G4int verboseLevel
Definition: G4VProcess.hh:368
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379
G4VQCrossSection::GetExchangeT
virtual G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4VQCrossSection.cc:70
G4VQCrossSection::GetCrossSection
virtual G4double GetCrossSection(G4bool, G4double, G4int, G4int, G4int pPDG=0)
Definition: G4VQCrossSection.hh:99
G4VQCrossSection::GetHMaxT
virtual G4double GetHMaxT()
Definition: G4VQCrossSection.cc:72
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83