Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
53 G4int G4QLowEnergy::nPartCWorld=85; // #of particles initialized in CHIPS Reduced
54 std::vector<G4int> G4QLowEnergy::ElementZ; // Z of element(i) in theLastCalc
55 std::vector<G4double> G4QLowEnergy::ElProbInMat; // SumProbOfElem in Material
56 std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;// N of isotope(j), E(i)
57 std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;//SumProbIsotE(i)
58 
59 // Constructor
60 G4QLowEnergy::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
73 G4QLowEnergy::~G4QLowEnergy() {}
74 
75 
76 G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;}
77 
78 G4int 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 ************
83 G4double 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 
228 G4bool 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 
247 G4VParticleChange* 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 
2273 G4double 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 }
Generated on Sat May 25 2013 14:33:51 for Geant4 by   doxygen 1.8.4