source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QLowEnergy.cc @ 1196

Last change on this file since 1196 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 44.8 KB
Line 
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: G4QLowEnergy.cc,v 1.9 2009/03/09 15:41:17 mkossov Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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
39//#define debug
40//#define pdebug
41//#define tdebug
42//#define nandebug
43//#define ppdebug
44
45#include "G4QLowEnergy.hh"
46
47// Initialization of static vectors
48G4int    G4QLowEnergy::nPartCWorld=152;  // #of particles initialized in CHIPS
49std::vector<G4int> G4QLowEnergy::ElementZ; // Z of element(i) in theLastCalc
50std::vector<G4double> G4QLowEnergy::ElProbInMat; // SumProbOfElem in Material
51std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;// N of isotope(j), E(i)
52std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;//SumProbIsotE(i)
53
54// Constructor
55G4QLowEnergy::G4QLowEnergy(const G4String& processName):
56  G4VDiscreteProcess(processName, fHadronic), evaporate(true)
57{
58#ifdef debug
59  G4cout<<"G4QLowEnergy::Constructor is called processName="<<processName<<G4endl;
60#endif
61  if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created "<<G4endl;
62  SetProcessSubType(fHadronInelastic);
63  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
64}
65
66// Destructor
67G4QLowEnergy::~G4QLowEnergy() {}
68
69
70G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;}
71
72G4int G4QLowEnergy::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
73
74// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
75// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
76// ********** All CHIPS cross sections are calculated in the surface units ************
77G4double G4QLowEnergy::GetMeanFreePath(const G4Track&Track, G4double, G4ForceCondition*F)
78{
79  *F = NotForced;
80  const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
81  G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
82  if( !IsApplicable(*incidentParticleDefinition))
83    G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl;
84  // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
85  G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
86#ifdef debug
87  G4double KinEn = incidentParticle->GetKineticEnergy();
88  G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
89#endif
90  const G4Material* material = Track.GetMaterial();        // Get the current material
91  const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
92  const G4ElementVector* theElementVector = material->GetElementVector();
93  G4int nE=material->GetNumberOfElements();
94#ifdef debug
95  G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl;
96#endif
97  G4int pPDG=0;
98  if      ( incidentParticleDefinition ==  G4Proton::Proton()         ) pPDG = 2212;
99  else if ( incidentParticleDefinition ==  G4Deuteron::Deuteron()     ) pPDG = 100001002;
100  else if ( incidentParticleDefinition ==  G4Alpha::Alpha()           ) pPDG = 100002004;
101  else if ( incidentParticleDefinition ==  G4Triton::Triton()         ) pPDG = 100001003;
102  else if ( incidentParticleDefinition ==  G4He3::He3()               ) pPDG = 100002003;
103  else if ( incidentParticleDefinition ==  G4GenericIon::GenericIon() )
104  {
105    pPDG=incidentParticleDefinition->GetPDGEncoding();
106#ifdef debug
107    G4int B=incidentParticleDefinition->GetBaryonNumber();
108    G4int C=incidentParticleDefinition->GetPDGCharge();
109    prPDG=100000000+1000*C+B;
110    G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl;
111#endif
112  }
113  else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl;
114  G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
115  if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer();
116  Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A
117  G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
118  G4double sigma=0.;                        // Sums over elements for the material
119  G4int IPIE=IsoProbInEl.size();            // How many old elements?
120  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
121  {
122    std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
123    SPI->clear();
124    delete SPI;
125    std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
126    IsN->clear();
127    delete IsN;
128  }
129  ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
130  ElementZ.clear();                         // Clear the body vector for Z of Elements
131  IsoProbInEl.clear();                      // Clear the body vector for SPI
132  ElIsoN.clear();                           // Clear the body vector for N of Isotopes
133  for(G4int i=0; i<nE; ++i)
134  {
135    G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
136    G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
137    ElementZ.push_back(Z);                  // Remember Z of the Element
138    G4int isoSize=0;                        // The default for the isoVectorLength is 0
139    G4int indEl=0;                          // Index of non-natural element or 0(default)
140    G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
141    if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
142#ifdef debug
143    G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
144#endif
145    if(isoSize)                             // The Element has non-trivial abundance set
146    {
147      indEl=pElement->GetIndex()+1;         // Index of the non-trivial element is an order
148#ifdef debug
149      G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
150#endif
151      if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
152      {
153        std::vector<std::pair<G4int,G4double>*>* newAbund =
154                                               new std::vector<std::pair<G4int,G4double>*>;
155        G4double* abuVector=pElement->GetRelativeAbundanceVector();
156        for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
157        {
158          G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
159          if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
160                                         <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
161          G4double abund=abuVector[j];
162          std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
163#ifdef debug
164          G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl;
165#endif
166          newAbund->push_back(pr);
167        }
168#ifdef debug
169        G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
170#endif
171        indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
172        for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
173        delete newAbund; // Was "new" in the beginning of the name space
174      }
175    }
176    std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
177    std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
178    IsoProbInEl.push_back(SPI);
179    std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
180    ElIsoN.push_back(IsN);
181    G4int nIs=cs->size();                   // A#Of Isotopes in the Element
182#ifdef debug
183    G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
184#endif
185    G4double susi=0.;                       // sum of CS over isotopes
186    if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
187    {
188      std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
189      G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
190      IsN->push_back(N);                    // Remember Min N for the Element
191#ifdef debug
192      G4cout<<"G4QLowE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
193#endif
194      G4bool ccsf=true;                    // Extract inelastic Ion-Ion cross-section
195#ifdef debug
196      G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl;
197#endif
198      G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope
199#ifdef debug
200      G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
201            <<Momentum<<", XSec="<<CSI/millibarn<<G4endl;
202#endif
203      curIs->second = CSI;
204      susi+=CSI;                            // Make a sum per isotopes
205      SPI->push_back(susi);                 // Remember summed cross-section
206    } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
207    sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
208#ifdef debug
209    G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
210          <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
211#endif
212    ElProbInMat.push_back(sigma);
213  } // End of LOOP over Elements
214  // Check that cross section is not zero and return the mean free path
215#ifdef debug
216  G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
217#endif
218  if(sigma > 0.) return 1./sigma;                 // Mean path [distance]
219  return DBL_MAX;
220}
221
222G4bool G4QLowEnergy::IsApplicable(const G4ParticleDefinition& particle) 
223{
224  if      (particle == *(     G4Proton::Proton()     )) return true;
225  else if (particle == *(    G4Neutron::Neutron()    )) return true;
226  else if (particle == *(   G4Deuteron::Deuteron()   )) return true;
227  else if (particle == *(      G4Alpha::Alpha()      )) return true;
228  else if (particle == *(     G4Triton::Triton()     )) return true;
229  else if (particle == *(        G4He3::He3()        )) return true;
230  else if (particle == *( G4GenericIon::GenericIon() )) return true;
231#ifdef debug
232  G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl;
233#endif
234  return false;
235}
236
237G4VParticleChange* G4QLowEnergy::PostStepDoIt(const G4Track& track, const G4Step& step)
238{
239  static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV
240  static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV
241  static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV
242  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV;
243  static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV;
244  static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV;
245  static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV;
246  static const G4ThreeVector zeroMom(0.,0.,0.);
247  static G4ParticleDefinition* aGamma    = G4Gamma::Gamma();
248  static G4ParticleDefinition* aProton   = G4Proton::Proton();
249  static G4ParticleDefinition* aNeutron  = G4Neutron::Neutron();
250  static G4ParticleDefinition* aLambda   = G4Lambda::Lambda();
251  static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
252  static G4ParticleDefinition* aTriton   = G4Triton::Triton();
253  static G4ParticleDefinition* aHe3      = G4He3::He3();
254  static G4ParticleDefinition* anAlpha   = G4Alpha::Alpha();
255  static const G4int nCh=26;                         // #of combinations
256  static G4QNucleus Nuc;                             // A fake nucleus to call Evaporation
257  //
258  //-------------------------------------------------------------------------------------
259  static G4bool CWinit = true;                       // CHIPS Warld needs to be initted
260  if(CWinit)
261  {
262    CWinit=false;
263    G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
264  }
265  //-------------------------------------------------------------------------------------
266  const G4DynamicParticle* projHadron = track.GetDynamicParticle();
267  const G4ParticleDefinition* particle=projHadron->GetDefinition();
268#ifdef debug
269  G4cout<<"G4QLowEnergy::PostStepDoIt: Before the GetMeanFreePath is called In4M="
270        <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
271        <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
272#endif
273  //G4ForceCondition cond=NotForced;
274  //GetMeanFreePath(track, -27., &cond);              // @@ ?? jus to update parameters?
275#ifdef debug
276  G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
277#endif
278  G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
279  G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
280  G4double Momentum = proj4M.rho();                   // @@ Just for the test purposes
281  if(std::fabs(Momentum-momentum)>.000001)
282    G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
283#ifdef pdebug
284  G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M/m="
285        <<proj4M<<proj4M.m()<<G4endl;
286#endif
287  if (!IsApplicable(*particle))  // Check applicability
288  {
289    G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl;
290    return 0;
291  }
292  const G4Material* material = track.GetMaterial();      // Get the current material
293  const G4ElementVector* theElementVector = material->GetElementVector();
294  G4int nE=material->GetNumberOfElements();
295#ifdef debug
296  G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
297#endif
298  G4int projPDG=0;                           // PDG Code prototype for the captured hadron
299  // Not all these particles are implemented yet (see Is Applicable)
300  if      (particle ==      G4Proton::Proton()     ) projPDG= 2212;
301  else if (particle ==     G4Neutron::Neutron()    ) projPDG= 2112;
302  else if (particle ==    G4Deuteron::Deuteron()   ) projPDG= 100001002;
303  else if (particle ==       G4Alpha::Alpha()      ) projPDG= 100002004;
304  else if (particle ==      G4Triton::Triton()     ) projPDG= 100001003;
305  else if (particle ==         G4He3::He3()        ) projPDG= 100002003;
306  else if (particle ==  G4GenericIon::GenericIon() )
307  {
308    projPDG=particle->GetPDGEncoding();
309#ifdef debug
310    G4int B=particle->GetBaryonNumber();
311    G4int C=particle->GetPDGCharge();
312    prPDG=100000000+1000*C+B;
313    G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl;
314#endif
315  }
316  else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl;
317#ifdef debug
318  G4int prPDG=particle->GetPDGEncoding();
319  G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
320#endif
321  if(!projPDG)
322  {
323    G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
324    return 0;
325  }
326  // Element treatment
327  G4int EPIM=ElProbInMat.size();
328#ifdef debug
329  G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
330#endif
331  G4int i=0;
332  if(EPIM>1)
333  {
334    G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
335    for(i=0; i<nE; ++i)
336    {
337#ifdef debug
338      G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl;
339#endif
340      if (rnd<ElProbInMat[i]) break;
341    }
342    if(i>=nE) i=nE-1;                        // Top limit for the Element
343  }
344  G4Element* pElement=(*theElementVector)[i];
345  G4int tZ=static_cast<G4int>(pElement->GetZ());
346#ifdef debug
347    G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl;
348#endif
349  if(tZ<=0)
350  {
351    G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl;
352    if(tZ<0) return 0;
353  }
354  std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
355  std::vector<G4int>* IsN = ElIsoN[i];     // Vector of "#of neutrons" in the isotope El[i]
356  G4int nofIsot=SPI->size();               // #of isotopes in the element i
357#ifdef debug
358  G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
359#endif
360  G4int j=0;
361  if(nofIsot>1)
362  {
363    G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
364    for(j=0; j<nofIsot; ++j)
365    {
366#ifdef debug
367      G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
368#endif
369      if(rndI < (*SPI)[j]) break;
370    }
371    if(j>=nofIsot) j=nofIsot-1;            // Top limit for the isotope
372  }
373  G4int tN =(*IsN)[j]; ;                    // Randomized number of neutrons
374#ifdef debug
375  G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl;
376#endif
377  if(tN<0)
378  {
379    G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl;
380    return 0;
381  }
382  nOfNeutrons=tN;                           // Remember it for the energy-momentum check
383#ifdef debug
384  G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl;
385#endif
386  if(tN<0)
387  {
388    G4cerr<<"*Warning*G4QLowEnergy::PostStepDoIt:Element with N="<<tN<< G4endl;
389    return 0;
390  }
391  aParticleChange.Initialize(track);
392#ifdef debug
393  G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl;
394#endif
395  G4double weight        = track.GetWeight();
396  G4double localtime     = track.GetGlobalTime();
397  G4ThreeVector position = track.GetPosition();
398#ifdef debug
399  G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl;
400#endif
401  G4TouchableHandle trTouchable = track.GetTouchableHandle();
402#ifdef debug
403  G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl;
404#endif
405  G4QPDGCode targQPDG(90000000+tZ*1000+tN);  // @@ use G4Ion and get rid of CHIPS World
406  G4double tM=targQPDG.GetMass();            // CHIPS target nucleus mass in MeV
407  G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness
408  G4int pZ=static_cast<G4int>(particle->GetPDGCharge());  // Charge of the projectile
409  G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile
410  G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV
411  G4double cosp=-14*Momentum*(tM-pM)/tM/pM;  // Asymmetry power for angular distribution
412#ifdef debug
413  G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
414        <<","<<tN<<"), cosp="<<cosp<<G4endl;
415#endif
416  G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
417  G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
418  G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
419#ifdef debug
420  G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
421#endif
422  EnMomConservation=tot4M;                 // Total 4-mom of reaction for E/M conservation
423  // @@ Probably this is not necessary any more
424#ifdef debug
425  G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl;
426#endif
427  G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection
428#ifdef debug
429  G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl;
430#endif
431#ifdef nandebug
432  if(xSec>0. || xSec<0. || xSec==0);
433  else  G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
434#endif
435  // @@ check a possibility to separate p, n, or alpha (!)
436  if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
437  {
438#ifdef pdebug
439    G4cerr<<"-Warning-G4QLowEnergy::PSDoIt:*Zero cross-section* PDG="<<projPDG
440          <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
441#endif
442    //Do Nothing Action insead of the reaction
443    aParticleChange.ProposeEnergy(kinEnergy);
444    aParticleChange.ProposeLocalEnergyDeposit(0.);
445    aParticleChange.ProposeMomentumDirection(dir) ;
446    return G4VDiscreteProcess::PostStepDoIt(track,step);
447  }
448  // Kill interacting hadron
449  aParticleChange.ProposeEnergy(0.);
450  aParticleChange.ProposeTrackStatus(fStopAndKill);
451#ifdef debug
452  G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<G4endl;
453#endif
454  // algorithm implementation --- STARTS HERE --- All calculations are in IU --------
455  G4double totM=tot4M.m(); // total CMS mass of the reaction
456  G4int totN=tN+pN;
457  G4int totZ=tZ+pZ;
458  // @@ Here mass[i] can be calculated if mass=0
459  G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
460                      0.,0.,0.,0.,0.,0.};
461  mass[0] = targQPDG.GetNuclMass(totZ,totN,0); // gamma+gamma
462#ifdef debug
463  G4cout<<"G4QLowEn::PSDI: Nuclear Mass="<<mass[0]<<G4endl;
464#endif
465  if (totN>0 && totZ>0)
466  {
467    mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron
468    mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton
469  }
470  if ( (totZ > 0 && totN > 1) || (totN > 0 && totZ > 1) ) 
471     mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d
472  if ( (totZ > 0 && totN > 2) || (totN > 1 && totZ > 1) ) 
473     mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t
474  if ( (totZ > 2 && totN > 0) || (totN > 1 && totZ > 1) ) 
475     mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3
476  if ( (totZ > 2 && totN > 1) || (totN > 2 && totZ > 1) ) 
477     mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a
478  if ( (totZ > 0 && totN > 1) || totN > 2 ) 
479     mass[7] = targQPDG.GetNuclMass(totZ,totN-2,0); // n+n
480  mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron)
481  if ( (totZ > 1 && totN > 0) || totZ > 2 ) 
482     mass[9] = targQPDG.GetNuclMass(totZ-2,totN,0); // p+p
483  mass[10] = mass[5]; // proton+deuteron (the same as He3)
484  mass[11] = mass[4]; // neutron+deuteron (the same as t)
485  mass[12] = mass[6]; // deuteron+deuteron (the same as alpha)
486  mass[13] = mass[6]; // proton+tritium (the same as alpha)
487  if ( totN > 3 || (totN > 2 && totZ > 0) ) 
488     mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0); // n+t
489  if ( totZ > 3 || (totZ > 2 && totN > 0) ) 
490     mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0); // He3+p
491  mass[16] = mass[6]; // neutron+He3 (the same as alpha)
492  if ( (totZ > 3 && totN > 1) || (totZ > 2 && totN > 2) ) 
493    mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0); // pa
494  if ( (totZ > 1 && totN > 3) || (totZ > 2 && totN > 2) ) 
495    mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0); // na
496  if(pL>0)
497  {
498    G4int pL1=pL-1;
499    if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ  ,totN  ,pL1);// Lambda+gamma
500    if( (totN > 0 && totZ > 0) ||        totZ > 1 ) 
501      mass[20]=targQPDG.GetNuclMass(totZ-1,totN  ,pL1);//Lp
502    if( (totN > 0 && totZ > 0) || totN > 0        ) 
503      mass[21]=targQPDG.GetNuclMass(totZ  ,totN-1,pL1);//Ln
504    if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) ) 
505      mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld
506    if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) ) 
507      mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt
508    if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) ) 
509      mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3
510    if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) ) 
511      mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La
512  }
513#ifdef debug
514  G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl;
515#endif
516  G4double tA=tZ+tN; 
517  G4double pA=pZ+pN; 
518  G4double prZ=pZ/pA+tZ/tA;
519  G4double prN=pN/pA+tN/tA;
520  G4double prD=prN*prZ;
521  G4double prA=prD*prD;
522  G4double prH=prD*prZ;
523  G4double prT=prD*prN;
524  G4double fhe3=6.*std::sqrt(tA);
525  G4double prL=0.;
526  if(pL>0) prL=pL/pA;
527  G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
528                      0.,0.,0.,0.,0.,0.};
529  qval[ 0] =  totM - mass[ 0];
530  qval[ 1] = (totM - mass[ 1] - mNeut)*prN;
531  qval[ 2] = (totM - mass[ 2] - mProt)*prZ;
532  qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3.;
533  qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6.;
534  qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3;
535  qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9.;
536  qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN;
537  qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD;
538  qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ;
539  qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.;
540  qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.;
541  qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.;
542  qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.;
543  qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.;
544  qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3;
545  qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3;
546  qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.;
547  qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.;
548  if(pZ>0)
549  {
550    qval[19] = (totM - mass[19] - mLamb)*prL;
551    qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ;
552    qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN;
553    qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.;
554    qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.;
555    qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3;
556    qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4;
557  }
558#ifdef debug
559  G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<"prA="<<pA<<G4endl;
560#endif
561
562  if( !pZ && pN==1)
563  {
564    if(G4UniformRand()>(tA-1.)*(tA-1.)/52900.) qval[0] = 0.0;
565    if(G4UniformRand()>kinEnergy/7.925*tA) qval[2]=qval[3]=qval[4]=qval[5]=qval[9]=0.;
566  }
567  else qval[0] = 0.0;
568 
569  G4double qv = 0.0;                        // Total sum of probabilities (q-values)
570  for(G4int i=0; i<nCh; ++i )
571  {
572    if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels
573    if( qval[i] < 0.0 )      qval[i] = 0.0; // Close the splitting impossible channels
574    qv += qval[i];
575  }
576  // Select the channel
577  G4double qv1 = 0.0;
578  G4double ran = G4UniformRand();
579  G4int index  = 0;
580  for( index=0; index<nCh; ++index )
581  {
582    if( qval[index] > 0.0 )
583    {
584      qv1 += qval[index]/qv;
585      if( ran <= qv1 ) break;
586    }
587  }
588#ifdef debug
589  G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl;
590#endif
591  if(index == nCh)
592  {
593    G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl;
594    throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound");
595  }
596  // @@ Convert it to G4QHadrons   
597  G4DynamicParticle* ResSec = new G4DynamicParticle;
598  G4DynamicParticle* FstSec = new G4DynamicParticle;
599  G4DynamicParticle* SecSec = new G4DynamicParticle;
600#ifdef debug
601  G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl;
602#endif
603
604  G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype
605  G4double mF=0.;
606  G4double mS=0.;
607  G4int    rA=totZ+totN;
608  G4int    rZ=totZ;
609  G4int    rL=pL;
610  G4int complete=3;
611  G4ParticleDefinition* theDefinition;  // Prototype for qfNucleon
612  switch( index )
613  {
614   case 0:
615     if(!evaporate || rA<2)
616     {
617       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
618       if(!theDefinition)
619         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
620       ResSec->SetDefinition( theDefinition );
621       FstSec->SetDefinition( aGamma );
622       SecSec->SetDefinition( aGamma );
623     }
624     else
625     {
626       delete ResSec;
627       delete FstSec;
628       delete SecSec;
629       complete=0;
630     }
631     break;
632   case 1:
633     rA-=1;
634     if(!evaporate || rA<2)
635     {
636       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
637       if(!theDefinition)
638         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
639       ResSec->SetDefinition( theDefinition );
640       SecSec->SetDefinition( aGamma );
641     }
642     else
643     {
644       delete ResSec;
645       delete SecSec;
646       complete=1;
647     }
648     FstSec->SetDefinition( aNeutron );
649     mF=mNeut; // First hadron 4-momentum
650     break;
651   case 2:
652     rA-=1;
653     rZ-=1;
654     if(!evaporate || rA<2)
655     {
656       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
657       if(!theDefinition)
658         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
659       ResSec->SetDefinition( theDefinition );
660       SecSec->SetDefinition( aGamma );
661     }
662     else
663     {
664       delete ResSec;
665       delete SecSec;
666       complete=1;
667     }
668     FstSec->SetDefinition( aProton );
669     mF=mProt; // First hadron 4-momentum
670     break;
671   case 3:
672     rA-=2;
673     rZ-=1;
674     if(!evaporate || rA<2)
675     {
676       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
677       if(!theDefinition)
678         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
679       ResSec->SetDefinition( theDefinition );
680       SecSec->SetDefinition( aGamma );
681     }
682     else
683     {
684       delete ResSec;
685       delete SecSec;
686       complete=1;
687     }
688     FstSec->SetDefinition( aDeuteron );
689     mF=mDeut; // First hadron 4-momentum
690     break;
691   case 4:
692     rA-=3;
693     rZ-=1;
694     if(!evaporate || rA<2)
695     {
696       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
697       if(!theDefinition)
698         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
699       ResSec->SetDefinition( theDefinition );
700       SecSec->SetDefinition( aGamma );
701     }
702     else
703     {
704       delete ResSec;
705       delete SecSec;
706       complete=1;
707     }
708     FstSec->SetDefinition( aTriton );
709     mF=mTrit; // First hadron 4-momentum
710     break;
711   case 5:
712     rA-=3;
713     rZ-=2;
714     if(!evaporate || rA<2)
715     {
716       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
717       if(!theDefinition)
718         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
719       ResSec->SetDefinition( theDefinition );
720       SecSec->SetDefinition( aGamma );
721     }
722     else
723     {
724       delete ResSec;
725       delete SecSec;
726       complete=1;
727     }
728     FstSec->SetDefinition( aHe3);
729     mF=mHel3; // First hadron 4-momentum
730     break;
731   case 6:
732     rA-=4;
733     rZ-=2;
734     if(!evaporate || rA<2)
735     {
736       theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
737       if(!theDefinition)
738         G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
739       ResSec->SetDefinition( theDefinition );
740       SecSec->SetDefinition( aGamma );
741     }
742     else
743     {
744       delete ResSec;
745       delete SecSec;
746       complete=1;
747     }
748     FstSec->SetDefinition( anAlpha );
749     mF=mAlph; // First hadron 4-momentum
750     break;
751   case 7:
752     rA-=2;
753     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
754     if(!theDefinition)
755       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
756     ResSec->SetDefinition( theDefinition );
757     FstSec->SetDefinition( aNeutron );
758     SecSec->SetDefinition( aNeutron );
759     mF=mNeut; // First hadron 4-momentum
760     mS=mNeut; // Second hadron 4-momentum
761     break;
762   case 8:
763     rZ-=1;
764     rA-=2;
765     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
766     if(!theDefinition)
767       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
768     ResSec->SetDefinition( theDefinition );
769     FstSec->SetDefinition( aNeutron );
770     SecSec->SetDefinition( aProton );
771     mF=mNeut; // First hadron 4-momentum
772     mS=mProt; // Second hadron 4-momentum
773     break;
774   case 9:
775     rZ-=2;
776     rA-=2;
777     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
778     if(!theDefinition)
779       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
780     ResSec->SetDefinition( theDefinition );
781     FstSec->SetDefinition( aProton );
782     SecSec->SetDefinition( aProton );
783     mF=mProt; // First hadron 4-momentum
784     mS=mProt; // Second hadron 4-momentum
785     break;
786   case 10:
787     rZ-=2;
788     rA-=3;
789     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
790     if(!theDefinition)
791       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
792     ResSec->SetDefinition( theDefinition );
793     FstSec->SetDefinition( aProton );
794     SecSec->SetDefinition( aDeuteron );
795     mF=mProt; // First hadron 4-momentum
796     mS=mDeut; // Second hadron 4-momentum
797     break;
798   case 11:
799     rZ-=1;
800     rA-=3;
801     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
802     if(!theDefinition)
803       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
804     ResSec->SetDefinition( theDefinition );
805     FstSec->SetDefinition( aNeutron );
806     SecSec->SetDefinition( aDeuteron );
807     mF=mNeut; // First hadron 4-momentum
808     mS=mDeut; // Second hadron 4-momentum
809     break;
810   case 12:
811     rZ-=2;
812     rA-=4;
813     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
814     if(!theDefinition)
815       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
816     ResSec->SetDefinition( theDefinition );
817     FstSec->SetDefinition( aDeuteron );
818     SecSec->SetDefinition( aDeuteron );
819     mF=mDeut; // First hadron 4-momentum
820     mS=mDeut; // Second hadron 4-momentum
821     break;
822   case 13:
823     rZ-=2;
824     rA-=4;
825     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
826     if(!theDefinition)
827       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
828     ResSec->SetDefinition( theDefinition );
829     FstSec->SetDefinition( aProton );
830     SecSec->SetDefinition( aTriton );
831     mF=mProt; // First hadron 4-momentum
832     mS=mTrit; // Second hadron 4-momentum
833     break;
834   case 14:
835     rZ-=1;
836     rA-=4;
837     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
838     if(!theDefinition)
839       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
840     ResSec->SetDefinition( theDefinition );
841     FstSec->SetDefinition( aNeutron );
842     SecSec->SetDefinition( aTriton );
843     mF=mNeut; // First hadron 4-momentum
844     mS=mTrit; // Second hadron 4-momentum
845     break;
846   case 15:
847     rZ-=3;
848     rA-=4;
849     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
850     if(!theDefinition)
851       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
852     ResSec->SetDefinition( theDefinition );
853     FstSec->SetDefinition( aProton);
854     SecSec->SetDefinition( aHe3 );
855     mF=mProt; // First hadron 4-momentum
856     mS=mHel3; // Second hadron 4-momentum
857     break;
858   case 16:
859     rZ-=2;
860     rA-=4;
861     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
862     if(!theDefinition)
863       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
864     ResSec->SetDefinition( theDefinition );
865     FstSec->SetDefinition( aNeutron );
866     SecSec->SetDefinition( aHe3 );
867     mF=mNeut; // First hadron 4-momentum
868     mS=mHel3; // Second hadron 4-momentum
869     break;
870   case 17:
871     rZ-=3;
872     rA-=5;
873     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
874     if(!theDefinition)
875       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
876     ResSec->SetDefinition( theDefinition );
877     FstSec->SetDefinition( aProton );
878     SecSec->SetDefinition( anAlpha );
879     mF=mProt; // First hadron 4-momentum
880     mS=mAlph; // Second hadron 4-momentum
881     break;
882   case 18:
883     rZ-=2;
884     rA-=5;
885     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
886     if(!theDefinition)
887       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
888     ResSec->SetDefinition( theDefinition );
889     FstSec->SetDefinition( aNeutron );
890     SecSec->SetDefinition( anAlpha );
891     mF=mNeut; // First hadron 4-momentum
892     mS=mAlph; // Second hadron 4-momentum
893     break;
894   case 19:
895     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
896     if(!theDefinition)
897       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
898     ResSec->SetDefinition( theDefinition );
899     FstSec->SetDefinition( aLambda );
900     SecSec->SetDefinition( aGamma );
901     mF=mLamb; // First hadron 4-momentum
902     break;
903   case 20:
904     rZ-=1;
905     rA-=1;
906     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
907     if(!theDefinition)
908       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
909     ResSec->SetDefinition( theDefinition );
910     FstSec->SetDefinition( aProton );
911     SecSec->SetDefinition( aLambda );
912     mF=mProt; // First hadron 4-momentum
913     mS=mLamb; // Second hadron 4-momentum
914     break;
915   case 21:
916     rA-=1;
917     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
918     if(!theDefinition)
919       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
920     ResSec->SetDefinition( theDefinition );
921     FstSec->SetDefinition( aNeutron );
922     SecSec->SetDefinition( aLambda );
923     mF=mNeut; // First hadron 4-momentum
924     mS=mLamb; // Second hadron 4-momentum
925     break;
926   case 22:
927     rZ-=1;
928     rA-=2;
929     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
930     if(!theDefinition)
931       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
932     ResSec->SetDefinition( theDefinition );
933     FstSec->SetDefinition( aDeuteron );
934     SecSec->SetDefinition( aLambda );
935     mF=mDeut; // First hadron 4-momentum
936     mS=mLamb; // Second hadron 4-momentum
937     break;
938   case 23:
939     rZ-=1;
940     rA-=3;
941     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
942     if(!theDefinition)
943       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
944     ResSec->SetDefinition( theDefinition );
945     FstSec->SetDefinition( aTriton );
946     SecSec->SetDefinition( aLambda );
947     mF=mTrit; // First hadron 4-momentum
948     mS=mLamb; // Second hadron 4-momentum
949     break;
950   case 24:
951     rZ-=2;
952     rA-=3;
953     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
954     if(!theDefinition)
955       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
956     ResSec->SetDefinition( theDefinition );
957     FstSec->SetDefinition( aHe3 );
958     SecSec->SetDefinition( aLambda );
959     mF=mHel3; // First hadron 4-momentum
960     mS=mLamb; // Second hadron 4-momentum
961     break;
962   case 25:
963     rZ-=2;
964     rA-=4;
965     theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
966     if(!theDefinition)
967       G4cerr<<"-Warning-G4LE::PSDI: notDefined, Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
968     ResSec->SetDefinition( theDefinition );
969     FstSec->SetDefinition( anAlpha );
970     SecSec->SetDefinition( aLambda );
971     mF=mAlph; // First hadron 4-momentum
972     mS=mLamb; // Second hadron 4-momentum
973     break;
974  }
975#ifdef debug
976  G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl;
977#endif
978  G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum
979  G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum
980  G4LorentzVector dir4Mom=tot4M;       // Prototype of the resN decay direction 4-momentum
981  dir4Mom.setE(tot4M.e()/2.);          // Get half energy and total 3-momentum
982  // @@ Can be repeated to take into account the Coulomb Barrier
983  if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp))
984  {                   //                                         
985    G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2="
986     <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl;
987    throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound");
988  }                   //                                         
989#ifdef debug
990  G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl;
991#endif
992  G4Track* aNewTrack = 0;
993  if(complete)
994  {
995    FstSec->Set4Momentum(fst4Mom);
996    aNewTrack = new G4Track(FstSec, localtime, position );
997    aNewTrack->SetWeight(weight);                                   //    weighted
998    aNewTrack->SetTouchableHandle(trTouchable);
999    aParticleChange.AddSecondary( aNewTrack );
1000    EnMomConservation-=fst4Mom;
1001#ifdef debug
1002    G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom<<G4endl;
1003#endif
1004    if(complete>2)                     // Final solution
1005    {
1006      ResSec->Set4Momentum(res4Mom);
1007      aNewTrack = new G4Track(ResSec, localtime, position );
1008      aNewTrack->SetWeight(weight);                                   //    weighted
1009      aNewTrack->SetTouchableHandle(trTouchable);
1010      aParticleChange.AddSecondary( aNewTrack );
1011      EnMomConservation-=res4Mom;
1012#ifdef debug
1013      G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL="
1014            <<rL<<G4endl;
1015#endif
1016      SecSec->Set4Momentum(snd4Mom);
1017      aNewTrack = new G4Track(SecSec, localtime, position );
1018      aNewTrack->SetWeight(weight);                                   //    weighted
1019      aNewTrack->SetTouchableHandle(trTouchable);
1020      aParticleChange.AddSecondary( aNewTrack );
1021      EnMomConservation-=snd4Mom;
1022#ifdef debug
1023      G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom<<G4endl;
1024#endif
1025    }
1026    else res4Mom+=snd4Mom;
1027  }
1028  else res4Mom=tot4M;
1029  if(complete<3)                        // Evaporation of the residual must be done
1030  {
1031    G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus
1032    G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!)
1033    Nuc.EvaporateNucleus(rHadron, evaHV);
1034    G4int nOut=evaHV->size();
1035    for(G4int i=0; i<nOut; i++)
1036    {
1037      G4QHadron* curH = (*evaHV)[i];
1038      G4int hPDG=curH->GetPDGCode();
1039      G4LorentzVector h4Mom=curH->Get4Momentum();
1040      EnMomConservation-=h4Mom;
1041#ifdef debug
1042      G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl;
1043#endif
1044      if     (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron;
1045      else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton;
1046      else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda;
1047      else if(hPDG==     22 )               theDefinition = aGamma;
1048      else
1049      {
1050        G4int hZ=curH->GetCharge();
1051        G4int hA=curH->GetBaryonNumber();
1052        G4int hS=curH->GetStrangeness();
1053        theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion
1054      }
1055      if(theDefinition)
1056      {
1057        G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom);
1058        G4Track* evaQH = new G4Track(theEQH, localtime, position );
1059        evaQH->SetWeight(weight);                                   //    weighted
1060        evaQH->SetTouchableHandle(trTouchable);
1061        aParticleChange.AddSecondary( evaQH );
1062      }
1063      else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl;
1064    }
1065  }
1066  // algorithm implementation --- STOPS HERE
1067#ifdef debug
1068  G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
1069#endif
1070  return G4VDiscreteProcess::PostStepDoIt(track, step);
1071}
1072
1073G4double G4QLowEnergy::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG) 
1074{
1075  static G4bool first=true;
1076  static G4VQCrossSection* CSmanager;
1077  if(first)                              // Connection with a singletone
1078  {
1079    CSmanager=G4QIonIonCrossSection::GetPointer();
1080    if(PDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer();
1081    first=false;
1082  }
1083#ifdef debug
1084  G4cout<<"G4QLowE::CXS: *DONE* p="<<p<<",Z="<<Z<<",N="<<N<<",PDG="<<PDG<<G4endl;
1085#endif
1086  return CSmanager->GetCrossSection(true, p, Z, N, PDG);
1087}
Note: See TracBrowser for help on using the repository browser.