Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QDiffraction.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 //
28 // ---------------- G4QDiffraction class -----------------
29 // by Mikhail Kossov, Aug 2007.
30 // G4QDiffraction class of the CHIPS Simulation Branch in GEANT4
31 // ---------------------------------------------------------------
32 // Short description: This is a process, which describes the diffraction
33 // excitation of the projectile and the nucleus. On nuclei in addition there
34 // can be a coherent diffraction process for the projectile, but it is
35 // comparably small. The most important part of the diffraction is the
36 // progectile diffraction excitation, as in this interaction proton can lose
37 // only a small part of its energy and make the shower longer. This is because
38 // only 1-2 (n) pions are produce in the diffraction escitation, and the mean
39 // kept energy of the nucleon is (1-n/7)=80%. For kaons the kept energy is much
40 // smaller (1-n/3.5)=60%, and for pions it is less important (about 40%).
41 // ----------------------------------------------------------------------------
42 
43 //#define debug
44 //#define pdebug
45 //#define tdebug
46 //#define nandebug
47 //#define ppdebug
48 
49 #include "G4QDiffraction.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4HadronicDeprecate.hh"
52 
53 
54 // Initialization of static vectors
55 //G4int G4QDiffraction::nPartCWorld=152; // #of particles initialized in CHIPS
56 //G4int G4QDiffraction::nPartCWorld=122; // #of particles initialized in CHIPS
57 G4int G4QDiffraction::nPartCWorld=85; // #of particles initialized in CHIPS Reduced
58 std::vector<G4int> G4QDiffraction::ElementZ; // Z of element(i) in theLastCalc
59 std::vector<G4double> G4QDiffraction::ElProbInMat; // SumProbOfElem in Material
60 std::vector<std::vector<G4int>*> G4QDiffraction::ElIsoN;// N of isotope(j), E(i)
61 std::vector<std::vector<G4double>*>G4QDiffraction::IsoProbInEl;//SumProbIsotE(i)
62 
63 // Constructor
64 G4QDiffraction::G4QDiffraction(const G4String& processName):
65  G4VDiscreteProcess(processName, fHadronic)
66 {
67  G4HadronicDeprecate("G4QDiffraction");
68 
69 #ifdef debug
70  G4cout<<"G4QDiffraction::Constructor is called processName="<<processName<<G4endl;
71 #endif
72  if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
73  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
74 }
75 
76 // Destructor
77 G4QDiffraction::~G4QDiffraction() {}
78 
79 
80 G4LorentzVector G4QDiffraction::GetEnegryMomentumConservation(){return EnMomConservation;}
81 
82 G4int G4QDiffraction::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
83 
84 // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
85 // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
86 // ********** All CHIPS cross sections are calculated in the surface units ************
87 G4double G4QDiffraction::GetMeanFreePath(const G4Track& Track,G4double,G4ForceCondition *F)
88 {
89  *F = NotForced;
90  const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
91  G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
92  if( !IsApplicable(*incidentParticleDefinition))
93  G4cout<<"-Warning-G4QDiffraction::GetMeanFreePath for notImplemented Particle"<<G4endl;
94  // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
95  G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
96 #ifdef debug
97  G4double KinEn = incidentParticle->GetKineticEnergy();
98  G4cout<<"G4QDiffraction::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
99 #endif
100  const G4Material* material = Track.GetMaterial(); // Get the current material
101  const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
102  const G4ElementVector* theElementVector = material->GetElementVector();
103  G4int nE=material->GetNumberOfElements();
104 #ifdef debug
105  G4cout<<"G4QDiffraction::GetMeanFreePath:"<<nE<<" Elems in Material="<<*material<<G4endl;
106 #endif
107  G4int pPDG=0;
108  // @@ At present it is made only for n & p, but can be extended if inXS are available
109  if (incidentParticleDefinition == G4Proton::Proton() ) pPDG=2212;
110  else if(incidentParticleDefinition == G4Neutron::Neutron()) pPDG=2112;
111  else G4cout<<"G4QDiffraction::GetMeanFreePath: only nA & pA are implemented"<<G4endl;
112 
113  G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
114  G4double sigma=0.; // Sums over elements for the material
115  G4int IPIE=IsoProbInEl.size(); // How many old elements?
116  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
117  {
118  std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
119  SPI->clear();
120  delete SPI;
121  std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
122  IsN->clear();
123  delete IsN;
124  }
125  ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
126  ElementZ.clear(); // Clear the body vector for Z of Elements
127  IsoProbInEl.clear(); // Clear the body vector for SPI
128  ElIsoN.clear(); // Clear the body vector for N of Isotopes
129  for(G4int i=0; i<nE; ++i)
130  {
131  G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
132  G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
133  ElementZ.push_back(Z); // Remember Z of the Element
134  G4int isoSize=0; // The default for the isoVectorLength is 0
135  G4int indEl=0; // Index of non-natural element or 0(default)
136  G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
137  if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
138 #ifdef debug
139  G4cout<<"G4QDiffraction::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
140 #endif
141  if(isoSize) // The Element has non-trivial abundance set
142  {
143  indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order
144 #ifdef debug
145  G4cout<<"G4QDiffr::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
146 #endif
147  if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
148  {
149  std::vector<std::pair<G4int,G4double>*>* newAbund =
150  new std::vector<std::pair<G4int,G4double>*>;
151  G4double* abuVector=pElement->GetRelativeAbundanceVector();
152  for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
153  {
154  G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
155  if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
156  <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
157  G4double abund=abuVector[j];
158  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
159 #ifdef debug
160  G4cout<<"G4QDiffract::GetMeanFreePath:pair#"<<j<<",N="<<N<<",ab="<<abund<<G4endl;
161 #endif
162  newAbund->push_back(pr);
163  }
164 #ifdef debug
165  G4cout<<"G4QDiffract::GetMeanFreePath:pairVectorLength="<<newAbund->size()<<G4endl;
166 #endif
167  indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
168  for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
169  delete newAbund; // Was "new" in the beginning of the name space
170  }
171  }
172  std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
173  std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
174  IsoProbInEl.push_back(SPI);
175  std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
176  ElIsoN.push_back(IsN);
177  G4int nIs=cs->size(); // A#Of Isotopes in the Element
178 #ifdef debug
179  G4cout<<"G4QDiffract::GetMFP:=***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
180 #endif
181  G4double susi=0.; // sum of CS over isotopes
182  if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
183  {
184  std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
185  G4int N=curIs->first; // #of Neuterons in the isotope j of El i
186  IsN->push_back(N); // Remember Min N for the Element
187 #ifdef debug
188  G4cout<<"G4Q::GMFP:j="<<j<<",P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PD="<<pPDG<<G4endl;
189 #endif
190  G4double CSI=CalculateXS(Momentum, Z, N, pPDG); // XS(j,i) for theIsotope
191 
192 #ifdef debug
193  G4cout<<"G4QDiffraction::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
194  <<Momentu<<", XSec="<<CSI/millibarn<<G4endl;
195 #endif
196  curIs->second = CSI;
197  susi+=CSI; // Make a sum per isotopes
198  SPI->push_back(susi); // Remember summed cross-section
199  } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
200  sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
201 #ifdef debug
202  G4cout<<"G4QDiffraction::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
203  <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
204 #endif
205  ElProbInMat.push_back(sigma);
206  } // End of LOOP over Elements
207  // Check that cross section is not zero and return the mean free path
208 #ifdef debug
209  G4cout<<"G4QDiffraction::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
210 #endif
211  if(sigma > 0.) return 1./sigma; // Mean path [distance]
212  return DBL_MAX;
213 }
214 
215 G4bool G4QDiffraction::IsApplicable(const G4ParticleDefinition& particle)
216 {
217  if (particle == *( G4Proton::Proton() )) return true;
218  else if (particle == *( G4Neutron::Neutron() )) return true;
219  //else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
220  //else if (particle == *( G4TauPlus::TauPlus() )) return true;
221  //else if (particle == *( G4TauMinus::TauMinus() )) return true;
222  //else if (particle == *( G4Electron::Electron() )) return true;
223  //else if (particle == *( G4Positron::Positron() )) return true;
224  //else if (particle == *( G4Gamma::Gamma() )) return true;
225  //else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
226  //else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
227  //else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
228  //else if (particle == *( G4PionMinus::PionMinus() )) return true;
229  //else if (particle == *( G4PionPlus::PionPlus() )) return true;
230  //else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
231  //else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
232  //else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
233  //else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
234  //else if (particle == *( G4Lambda::Lambda() )) return true;
235  //else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
236  //else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
237  //else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
238  //else if (particle == *( G4XiMinus::XiMinus() )) return true;
239  //else if (particle == *( G4XiZero::XiZero() )) return true;
240  //else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
241  //else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
242  //else if (particle == *( G4AntiProton::AntiProton() )) return true;
243 #ifdef debug
244  G4cout<<"***>>G4QDiffraction::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl;
245 #endif
246  return false;
247 }
248 
249 G4VParticleChange* G4QDiffraction::PostStepDoIt(const G4Track& track, const G4Step& step)
250 {
251  static const G4double mProt= G4QPDGCode(2212).GetMass(); // CHIPS proton Mass in MeV
252  static const G4double mNeut= G4QPDGCode(2112).GetMass(); // CHIPS neutron Mass in MeV
253  static const G4double mPion= G4QPDGCode(111).GetMass(); // CHIPS Pi0 Mass in MeV
254  static G4QDiffractionRatio* diffRatio;
255  //
256  //-------------------------------------------------------------------------------------
257  static G4bool CWinit = true; // CHIPS Warld needs to be initted
258  if(CWinit)
259  {
260  CWinit=false;
261  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
262  diffRatio=G4QDiffractionRatio::GetPointer();
263  }
264  //-------------------------------------------------------------------------------------
265  const G4DynamicParticle* projHadron = track.GetDynamicParticle();
266  const G4ParticleDefinition* particle=projHadron->GetDefinition();
267 #ifdef debug
268  G4cout<<"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M="
269  <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
270  <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
271 #endif
272  G4ForceCondition cond=NotForced;
273  GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters?
274 #ifdef debug
275  G4cout<<"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
276 #endif
277  G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
278  G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
279  G4double Momentum = proj4M.rho(); // @@ Just for the test purposes
280  if(std::fabs(Momentum-momentum)>.000001)
281  G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
282 #ifdef pdebug
283  G4cout<<"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum
284  <<",proj4M="<<proj4M<<", projM="<<proj4M.m()<<G4endl;
285 #endif
286  if (!IsApplicable(*particle)) // Check applicability
287  {
288  G4cerr<<"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<G4endl;
289  return 0;
290  }
291  const G4Material* material = track.GetMaterial(); // Get the current material
292  G4int Z=0;
293  const G4ElementVector* theElementVector = material->GetElementVector();
294  G4int nE=material->GetNumberOfElements();
295 #ifdef debug
296  G4cout<<"G4QDiffraction::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 == G4PionMinus::PionMinus() ) projPDG= -211;
303  //else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
304  //else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321;
305  //else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
306  //else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
307  //else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
308  //else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
309  //else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
310  //else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
311  //else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
312  //else if (particle == G4Electron::Electron() ) projPDG= 11;
313  //else if (particle == G4Positron::Positron() ) projPDG= -11;
314  //else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
315  //else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
316  //else if (particle == G4Gamma::Gamma() ) projPDG= 22;
317  //else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
318  //else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
319  //else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
320  //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
321  //else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
322  //else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
323  //else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
324  //else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
325  //else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
326  //else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
327  //else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
328  //else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
329  //else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
330 #ifdef debug
331  G4int prPDG=particle->GetPDGEncoding();
332  G4cout<<"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
333 #endif
334  if(!projPDG)
335  {
336  G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
337  return 0;
338  }
339  //G4double pM2=proj4M.m2(); // in MeV^2
340  //G4double pM=std::sqrt(pM2); // in MeV
341  G4double pM=mNeut;
342  if(projPDG==2112) pM=mProt;
343  // Element treatment
344  G4int EPIM=ElProbInMat.size();
345 #ifdef debug
346  G4cout<<"G4QDiffra::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
347 #endif
348  G4int i=0;
349  if(EPIM>1)
350  {
351  G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
352  for(i=0; i<nE; ++i)
353  {
354 #ifdef debug
355  G4cout<<"G4QDiffra::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
356 #endif
357  if (rnd<ElProbInMat[i]) break;
358  }
359  if(i>=nE) i=nE-1; // Top limit for the Element
360  }
361  G4Element* pElement=(*theElementVector)[i];
362  Z=static_cast<G4int>(pElement->GetZ());
363 #ifdef debug
364  G4cout<<"G4QDiffraction::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
365 #endif
366  if(Z<=0)
367  {
368  G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<G4endl;
369  if(Z<0) return 0;
370  }
371  std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
372  std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
373  G4int nofIsot=SPI->size(); // #of isotopes in the element i
374 #ifdef debug
375  G4cout<<"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
376 #endif
377  G4int j=0;
378  if(nofIsot>1)
379  {
380  G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
381  for(j=0; j<nofIsot; ++j)
382  {
383 #ifdef debug
384  G4cout<<"G4QDiffraction::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
385 #endif
386  if(rndI < (*SPI)[j]) break;
387  }
388  if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
389  }
390  G4int N =(*IsN)[j]; ; // Randomized number of neutrons
391 #ifdef debug
392  G4cout<<"G4QDiffraction::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl;
393 #endif
394  if(N<0)
395  {
396  G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<" has 0>N="<<N<<G4endl;
397  return 0;
398  }
399  nOfNeutrons=N; // Remember it for the energy-momentum check
400 #ifdef debug
401  G4cout<<"G4QDiffraction::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
402 #endif
403  if(N<0)
404  {
405  G4cerr<<"*Warning*G4QDiffraction::PostStepDoIt:Element with N="<<N<< G4endl;
406  return 0;
407  }
408  aParticleChange.Initialize(track);
409 #ifdef debug
410  G4cout<<"G4QDiffraction::PostStepDoIt: track is initialized"<<G4endl;
411 #endif
412  G4double weight = track.GetWeight();
413  G4double localtime = track.GetGlobalTime();
414  G4ThreeVector position = track.GetPosition();
415 #ifdef debug
416  G4cout<<"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<G4endl;
417 #endif
418  G4TouchableHandle trTouchable = track.GetTouchableHandle();
419 #ifdef debug
420  G4cout<<"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<G4endl;
421 #endif
422  G4int targPDG=90000000+Z*1000+N; // CHIPS PDG Code of the target nucleus
423  G4QPDGCode targQPDG(targPDG); // @@ use G4Ion and get rid of CHIPS World
424  G4double tM=targQPDG.GetMass(); // CHIPS final nucleus mass in MeV
425  G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
426  G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
427  G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
428 #ifdef debug
429  G4cout<<"G4QDiffraction::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
430 #endif
431  EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation
432  // @@ Probably this is not necessary any more
433 #ifdef debug
434  G4cout<<"G4QDiff::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl;
435 #endif
436  G4double xSec=CalculateXS(Momentum, Z, N, projPDG); // Recalculate CrossSection
437 #ifdef debug
438  G4cout<<"G4QDiffra::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl;
439 #endif
440 #ifdef nandebug
441  if(xSec>0. || xSec<0. || xSec==0);
442  else G4cout<<"-Warning-G4QDiffraction::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
443 #endif
444  // @@ check a possibility to separate p, n, or alpha (!)
445  if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
446  {
447 #ifdef pdebug
448  G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Zero cross-section* PDG="<<projPDG
449  <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl;
450 #endif
451  //Do Nothing Action insead of the reaction
452  aParticleChange.ProposeEnergy(kinEnergy);
453  aParticleChange.ProposeLocalEnergyDeposit(0.);
454  aParticleChange.ProposeMomentumDirection(dir) ;
455  return G4VDiscreteProcess::PostStepDoIt(track,step);
456  }
457  G4double totCMMass=tot4M.m(); // Total CM mass, pM=projectileMass, tM=targetMass
458  if(totCMMass < mPion+pM+tM) // The diffraction reaction is impossible -> Do Nothing
459  {
460 #ifdef pdebug
461  G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Below Diffraction Threshold* cmM="<<totCMMass
462  <<">pM="<<pM<<"+tM="<<tM<<"+pi0="<<mPion<<"=="<<pM+tM+mPion<<G4endl;
463 #endif
464  //Do Nothing Action insead of the reaction
465  aParticleChange.ProposeEnergy(kinEnergy);
466  aParticleChange.ProposeLocalEnergyDeposit(0.);
467  aParticleChange.ProposeMomentumDirection(dir) ;
468  return G4VDiscreteProcess::PostStepDoIt(track,step);
469  }
470  // Kill interacting hadron
471  aParticleChange.ProposeTrackStatus(fStopAndKill);
472  G4QHadronVector* out=diffRatio->TargFragment(projPDG, proj4M, Z, N);
473  G4int nSec=out->size(); // #of secondaries in the diffraction reaction
474  G4DynamicParticle* theSec=0; // A prototype for secondary for the secondary
475  G4LorentzVector dif4M(0.,0.,0.,0.); // Prototype for the secondary 4-momentum
476  G4int difPDG=0; // PDG code of the secondary
477  G4QHadron* difQH=0; // Prototype for a Q-secondary
478 #ifdef pdebug
479  G4cout<<"G4QDiffraction::PostStepDoIt: =---=found=---= nSecondaries="<<nSec<<G4endl;
480 #endif
481  for(i=0; i<nSec; i++)
482  {
483  difQH = (*out)[i];
484  difPDG= difQH->GetPDGCode();
485  G4ParticleDefinition* theDefinition=0;
486  if (difPDG==2212 || difPDG==90001000) theDefinition=G4Proton::Proton();
487  else if(difPDG==2112 || difPDG==90000001) theDefinition=G4Neutron::Neutron();
488  else if(difPDG== 22) theDefinition=G4Gamma::Gamma();
489  else if(difPDG== 111) theDefinition=G4PionZero::PionZero();
490  else if(difPDG==-211 || difPDG==89999001) theDefinition=G4PionMinus::PionMinus();
491  else if(difPDG== 211 || difPDG==90000999) theDefinition=G4PionPlus::PionPlus();
492  else if(difPDG== 321 || difPDG==89001000) theDefinition=G4KaonPlus::KaonPlus();
493  else if(difPDG==-321 || difPDG==90999000) theDefinition=G4KaonMinus::KaonMinus();
494  else if(difPDG== 130 || difPDG==-311 || difPDG==89000001)
495  theDefinition=G4KaonZeroLong::KaonZeroLong();
496  else if(difPDG== 310 || difPDG== 311 || difPDG==90999999)
497  theDefinition=G4KaonZeroShort::KaonZeroShort();
498  else if(difPDG==3122 || difPDG==91000000) theDefinition=G4Lambda::Lambda();
499  else if(difPDG== 3222) theDefinition=G4SigmaPlus::SigmaPlus();
500  else if(difPDG== 3112) theDefinition=G4SigmaMinus::SigmaMinus();
501  else if(difPDG== 3212) theDefinition=G4SigmaZero::SigmaZero();
502  else if(difPDG== 3312) theDefinition=G4XiMinus::XiMinus();
503  else if(difPDG== 3322) theDefinition=G4XiZero::XiZero();
504  else if(difPDG== 3334) theDefinition=G4OmegaMinus::OmegaMinus();
505  else if(difPDG==-2112) theDefinition=G4AntiNeutron::AntiNeutron();
506  else if(difPDG==-2212) theDefinition=G4AntiProton::AntiProton();
507  else if(difPDG==-3122) theDefinition=G4AntiLambda::AntiLambda();
508  else if(difPDG==-3222) theDefinition=G4AntiSigmaPlus::AntiSigmaPlus();
509  else if(difPDG==-3112) theDefinition=G4AntiSigmaMinus::AntiSigmaMinus();
510  else if(difPDG==-3212) theDefinition=G4AntiSigmaZero::AntiSigmaZero();
511  else if(difPDG==-3312) theDefinition=G4AntiXiMinus::AntiXiMinus();
512  else if(difPDG==-3322) theDefinition=G4AntiXiZero::AntiXiZero();
513  else if(difPDG==-3334) theDefinition=G4AntiOmegaMinus::AntiOmegaMinus();
514  else if(difPDG== -11) theDefinition=G4Electron::Electron();
515  else if(difPDG== -13) theDefinition=G4MuonMinus::MuonMinus();
516  else if(difPDG== 11) theDefinition=G4Positron::Positron();
517  else if(difPDG== 13) theDefinition=G4MuonPlus::MuonPlus();
518  else
519  {
520  Z = difQH->GetCharge();
521  G4int B = difQH->GetBaryonNumber();
522  G4int S = difQH->GetStrangeness();
523  if(S||Z>B||Z<0)G4cout<<"-Warning-G4QDif::PoStDoIt:Z="<<Z<<",A="<<B<<",S="<<S<<G4endl;
524  theDefinition = G4ParticleTable::GetParticleTable()->FindIon(Z,B,0,0);
525 #ifdef pdebug
526  G4cout<<"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<",A="<<B<<",D="<<theDefinition<<G4endl;
527 #endif
528  }
529  if(theDefinition)
530  {
531  theSec = new G4DynamicParticle; // A secondary for the recoil hadron
532  theSec->SetDefinition(theDefinition);
533  dif4M = difQH->Get4Momentum();
534  EnMomConservation-=dif4M;
535  theSec->Set4Momentum(dif4M);
536  G4Track* aNewTrack = new G4Track(theSec, localtime, position );
537  aNewTrack->SetWeight(weight); // weighted
538  aNewTrack->SetTouchableHandle(trTouchable);
539  aParticleChange.AddSecondary( aNewTrack );
540 #ifdef pdebug
541  G4cout<<"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<", PDG="<<difPDG<<G4endl;
542 #endif
543  }
544  else G4cout<<"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<", Z="<<difQH->GetCharge()
545  <<", A="<<difQH->GetBaryonNumber()<<",S ="<<difQH->GetStrangeness()<<G4endl;
546  delete difQH; // Clean up the output QHadrons
547  }
548  delete out; // Delete the output QHadron-vector
549 #ifdef debug
550  G4cout<<"G4QDiffraction::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
551 #endif
552  return G4VDiscreteProcess::PostStepDoIt(track, step);
553 }
554 
555 G4double G4QDiffraction::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG)
556 {
557  static G4bool first=true;
558  //static G4VQCrossSection* CSmanager;
559  static G4QDiffractionRatio* diffRatio;
560  if(first) // Connection with a singletone
561  {
562  //CSmanager=G4QProtonNuclearCrossSection::GetPointer();
563  diffRatio=G4QDiffractionRatio::GetPointer();
564  first=false;
565  }
566  //G4double x=CSmanager->GetCrossSection(true, p, Z, N, PDG); // inelastic XS
567  //G4double pIU=p*GeV; // IndependentUnistMomentum
568  //G4double r=diffRatio->GetRatio(pIU, PDG, Z, N); // Proj. Diffraction Part
569  //G4double xs_value=x*r; // XS for proj. diffraction
570  G4double xs_value=diffRatio->GetTargSingDiffXS(p, PDG, Z, N); // XS for target diffraction
571 #ifdef debug
572  G4cout<<"G4QDiff::CXS:p="<<p<<",Z="<<Z<<",N="<<N<<",C="<<PDG<<",XS="<<xs_value<<G4endl;
573 #endif
574  return xs_value;
575 }
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4