Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QHadronInelasticDataSet.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 // GEANT4 physics class: G4QHadronInelasticDataSet -- header file
29 // Created by M. Kosov (Mikhail.Kossov@cern.ch) 11.11.09
30 //
31 // ------------------------------------------------------------------------
32 // Short description: G4hadr wrapper for CHIPS inelastic hA cross-sections.
33 // ------------------------------------------------------------------------
34 //
35 
36 #include "G4QHadronInelasticDataSet.hh"
37 #include <iostream>
38 
40 //std::vector<G4int> G4QHadronInelasticDataSet::ElementZ;
42 //
43 //std::vector<std::vector<G4int>*> G4QHadronInelasticDataSet::ElIsoN;
45 //
46 //std::vector<std::vector<G4double>*> G4QHadronInelasticDataSet::IsoProbInEl;
48 
49 G4QHadronInelasticDataSet::G4QHadronInelasticDataSet(const G4String& input_name)
50  : G4VCrossSectionDataSet(input_name)
51 {
52  //CHIPSpAin = G4QProtonNuclearCrossSection::GetPointer();
53  //CHIPSnAin = G4QNeutronNuclearCrossSection::GetPointer();
54  //CHIPSpimAin = G4QPionMinusNuclearCrossSection::GetPointer();
55  //CHIPSpipAin = G4QPionPlusNuclearCrossSection::GetPointer();
56  //CHIPSkpAin = G4QKaonPlusNuclearCrossSection::GetPointer();
57  //CHIPSkmAin = G4QKaonMinusNuclearCrossSection::GetPointer();
58  //CHIPSk0Ain = G4QKaonZeroNuclearCrossSection::GetPointer();
59  //CHIPShAin = G4QHyperonNuclearCrossSection::GetPointer();
60  //CHIPShpAin = G4QHyperonPlusNuclearCrossSection::GetPointer();
61  //CHIPSabpAin = G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
62  //CHIPSabAin = G4QAntiBaryonNuclearCrossSection::GetPointer();
73 
74  //Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
75  Description();
76 }
77 
78 void G4QHadronInelasticDataSet::Description() const
79 {
80  char* dirName = getenv("G4PhysListDocDir");
81  if (dirName) {
82  std::ofstream outFile;
83  G4String outFileName = GetName() + ".html";
84  G4String pathName = G4String(dirName) + "/" + outFileName;
85 
86  outFile.open(pathName);
87  outFile << "<html>\n";
88  outFile << "<head>\n";
89 
90  outFile << "<title>Description of CHIPSInelasticXS</title>\n";
91  outFile << "</head>\n";
92  outFile << "<body>\n";
93 
94  outFile << "CHIPSInelasticXS provides hadron-nuclear inelastic cross\n"
95  << "sections for all hadrons at all energies. These cross\n"
96  << "sections represent parameterizations developed by M. Kossov.\n";
97 
98  outFile << "</body>\n";
99  outFile << "</html>\n";
100  outFile.close();
101  }
102 }
103 
104 G4bool G4QHadronInelasticDataSet::IsIsoApplicable(const G4DynamicParticle* Pt, G4int,G4int,
105  const G4Element*, const G4Material*)
106 {
107  G4ParticleDefinition* particle = Pt->GetDefinition();
108  if (particle == G4Neutron::Neutron() ) return true; // @@ isotopes?
109  else if (particle == G4Proton::Proton() ) return true;
110  else if (particle == G4PionMinus::PionMinus() ) return true;
111  else if (particle == G4PionPlus::PionPlus() ) return true;
112  else if (particle == G4KaonPlus::KaonPlus() ) return true;
113  else if (particle == G4KaonMinus::KaonMinus() ) return true;
114  else if (particle == G4KaonZeroLong::KaonZeroLong() ) return true;
115  else if (particle == G4KaonZeroShort::KaonZeroShort() ) return true;
116  else if (particle == G4Lambda::Lambda() ) return true;
117  else if (particle == G4SigmaPlus::SigmaPlus() ) return true;
118  else if (particle == G4SigmaMinus::SigmaMinus() ) return true;
119  else if (particle == G4SigmaZero::SigmaZero() ) return true;
120  else if (particle == G4XiMinus::XiMinus() ) return true;
121  else if (particle == G4XiZero::XiZero() ) return true;
122  else if (particle == G4OmegaMinus::OmegaMinus() ) return true;
123  else if (particle == G4AntiNeutron::AntiNeutron() ) return true;
124  else if (particle == G4AntiProton::AntiProton() ) return true;
125  else if (particle == G4AntiLambda::AntiLambda() ) return true;
126  else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) return true;
127  else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
128  else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) return true;
129  else if (particle == G4AntiXiMinus::AntiXiMinus() ) return true;
130  else if (particle == G4AntiXiZero::AntiXiZero() ) return true;
131  else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
132  //else if (particle == G4MuonPlus::MuonPlus() ) return true;
133  //else if (particle == G4MuonMinus::MuonMinus() ) return true;
134  //else if (particle == G4Gamma::Gamma() ) return true;
135  //else if (particle == G4Electron::Electron() ) return true;
136  //else if (particle == G4Positron::Positron() ) return true;
137  //else if (particle == G4TauPlus::TauPlus() ) return true;
138  //else if (particle == G4TauMinus::TauMinus() ) return true;
139  //else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) return true;
140  //else if (particle == G4NeutrinoE::NeutrinoE() ) return true;
141  //else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) return true;
142  //else if (particle == G4NeutrinoMu::NeutrinoMu() ) return true;
143  //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) return true;
144  //else if (particle == G4NeutrinoTau::NeutrinoTau() ) return true;
145  return false;
146 }
147 
148 // Now it is in G4VCrossSectionDataSet...
149 /*
150 G4double G4QHadronInelasticDataSet::GetCrossSection(const G4DynamicParticle* Pt,
151  const G4Element* pElement, G4double)
152 {
153  G4int IPIE=IsoProbInEl.size(); // How many old elements?
154  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
155  {
156  std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
157  SPI->clear();
158  delete SPI;
159  std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
160  IsN->clear();
161  delete IsN;
162  }
163  ElementZ.clear(); // Clear the body vector for Z of Elements
164  IsoProbInEl.clear(); // Clear the body vector for SPI
165  ElIsoN.clear(); // Clear the body vector for N of Isotopes
166  G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
167  ElementZ.push_back(Z); // Remember Z of the Element
168  G4int isoSize=0; // The default for the isoVectorLength is 0
169  G4int indEl=0; // Index of non-trivial element or 0(default)
170  G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
171  if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
172  if(isoSize) // The Element has non-trivial abundance set
173  {
174  indEl=pElement->GetIndex()+1; // Index of the non-trivial element
175  if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
176  {
177  std::vector<std::pair<G4int,G4double>*>* newAbund =
178  new std::vector<std::pair<G4int,G4double>*>;
179  G4double* abuVector=pElement->GetRelativeAbundanceVector();
180  for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
181  {
182  G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
183  if(pElement->GetIsotope(j)->GetZ()!=Z)
184  G4cerr<<"G4QHadronInelasticDataSet::GetCrossSection"<<": Z="
185  <<pElement->GetIsotope(j)->GetZ()<<" # "<<Z<<G4endl;
186  G4double abund=abuVector[j];
187  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
188  newAbund->push_back(pr);
189  }
190  indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
191  for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
192  delete newAbund; // Was "new" in the beginning of the name space
193  }
194  }
195  std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
196  std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
197  IsoProbInEl.push_back(SPI);
198  std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
199  ElIsoN.push_back(IsN);
200  G4int nIs=cs->size(); // A#Of Isotopes in the Element
201  G4double susi=0.; // sum of CS over isotopes
202  if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
203  {
204  std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
205  G4int N=curIs->first; // #of Neuterons in the isotope j of El i
206  IsN->push_back(N); // Remember Min N for the Element
207  G4double CSI=GetIsoZACrossSection(Pt,Z,Z+N,0.);//CrossSection(j,i) for the isotope
208  curIs->second = CSI;
209  susi+=CSI; // Make a sum per isotopes
210  SPI->push_back(susi); // Remember summed cross-section
211  } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
212  return Isotopes->GetMeanCrossSection(Z,indEl); // MeanCS over isotopes
213 }
214 */
215 
216 G4double G4QHadronInelasticDataSet::GetIsoCrossSection(const G4DynamicParticle* Pt,G4int Z,
217  G4int A, const G4Isotope*,
218  const G4Element*, const G4Material*)
219 {
220  G4ParticleDefinition* particle = Pt->GetDefinition();
221  G4double Momentum=Pt->GetTotalMomentum();
222  G4VQCrossSection* CSmanager=0;
223 
224  G4int pPDG=0;
225  if(particle == G4Neutron::Neutron())
226  {
227  CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
228  pPDG=2112;
229  }
230  else if(particle == G4Proton::Proton())
231  {
232  CSmanager=G4QProtonNuclearCrossSection::GetPointer();
233  pPDG=2212;
234  }
235  else if(particle == G4PionMinus::PionMinus())
236  {
237  CSmanager=G4QPionMinusNuclearCrossSection::GetPointer();
238  pPDG=-211;
239  }
240  else if(particle == G4PionPlus::PionPlus())
241  {
242  CSmanager=G4QPionPlusNuclearCrossSection::GetPointer();
243  pPDG=211;
244  }
245  else if(particle == G4KaonMinus::KaonMinus())
246  {
247  CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer();
248  pPDG=-321;
249  }
250  else if(particle == G4KaonPlus::KaonPlus())
251  {
252  CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer();
253  pPDG=321;
254  }
255  else if(particle == G4KaonZeroLong::KaonZeroLong() ||
256  particle == G4KaonZeroShort::KaonZeroShort() ||
257  particle == G4KaonZero::KaonZero() ||
258  particle == G4AntiKaonZero::AntiKaonZero() )
259  {
260  CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer();
261  if(G4UniformRand() > 0.5) pPDG= 311;
262  else pPDG=-311;
263  }
264  else if(particle == G4Lambda::Lambda())
265  {
266  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
267  pPDG=3122;
268  }
269  else if(particle == G4SigmaPlus::SigmaPlus())
270  {
271  CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer();
272  pPDG=3222;
273  }
274  else if(particle == G4SigmaMinus::SigmaMinus())
275  {
276  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
277  pPDG=3112;
278  }
279  else if(particle == G4SigmaZero::SigmaZero())
280  {
281  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
282  pPDG=3212;
283  }
284  else if(particle == G4XiMinus::XiMinus())
285  {
286  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
287  pPDG=3312;
288  }
289  else if(particle == G4XiZero::XiZero())
290  {
291  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
292  pPDG=3322;
293  }
294  else if(particle == G4OmegaMinus::OmegaMinus())
295  {
296  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
297  pPDG=3334;
298  }
299  else if(particle == G4AntiNeutron::AntiNeutron())
300  {
301  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
302  pPDG=-2112;
303  }
304  else if(particle == G4AntiProton::AntiProton())
305  {
306  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
307  pPDG=-2212;
308  }
309  else if(particle == G4AntiLambda::AntiLambda())
310  {
311  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
312  pPDG=-3122;
313  }
314  else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
315  {
316  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
317  pPDG=-3222;
318  }
319  else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
320  {
321  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
322  pPDG=-3112;
323  }
324  else if(particle == G4AntiSigmaZero::AntiSigmaZero())
325  {
326  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
327  pPDG=-3212;
328  }
329  else if(particle == G4AntiXiMinus::AntiXiMinus())
330  {
331  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
332  pPDG=-3312;
333  }
334  else if(particle == G4AntiXiZero::AntiXiZero())
335  {
336  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
337  pPDG=-3322;
338  }
339  else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
340  {
341  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
342  pPDG=-3334;
343  }
344  //else if(particle == G4Gamma::Gamma())
345  //{
346  // CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
347  // pPDG=22;
348  //}
349  //else if(particle == G4Electron::Electron() ||
350  // particle == G4Positron::Positron())
351  //{
352  // CSmanager=G4QElectronNuclearCrossSection::GetPointer();
353  // pPDG=11;
354  //}
355  //else if(particle == G4MuonPlus::MuonPlus() ||
356  // particle == G4MuonMinus::MuonMinus())
357  //{
358  // CSmanager=G4QMuonNuclearCrossSection::GetPointer();
359  // pPDG=13;
360  //}
361  //else if(particle == G4TauPlus::TauPlus() ||
362  // particle == G4TauMinus::TauMinus())
363  //{
364  // CSmanager=G4QTauNuclearCrossSection::GetPointer();
365  // pPDG=15;
366  //}
367  //else if(particle == G4NeutrinoMu::NeutrinoMu() )
368  //{
369  // CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
370  // CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
371  // pPDG=14;
372  //}
373  //else if(particle == G4AntiNeutrinoMu::AntiNeutrinoMu() )
374  //{
375  // CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
376  // CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
377  // pPDG=-14;
378  //}
379  //else if(particle == G4NeutrinoE::NeutrinoE() )
380  //{
381  // CSmanager=G4QNuENuclearCrossSection::GetPointer();
382  // CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
383  // pPDG=12;
384  //}
385  //else if(particle == G4AntiNeutrinoE::AntiNeutrinoE() )
386  //{
387  // CSmanager=G4QANuENuclearCrossSection::GetPointer();
388  // CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
389  // pPDG=-12;
390  //}
391  else
392  {
393  G4cerr << "-ERROR-G4QHadronInelasticDataSet::GetIsoZACrossSection: PDG="
394  << particle->GetPDGEncoding() << " isn't supported by CHIPS" << G4endl;
395  //throw G4HadronicException(__FILE__, __LINE__,
396  //"G4QHadronInelasticDataSet::GetIsoZACrossSection: Particle not supported by CHIPS");
397  return 0;
398  }
399  G4int N=A-Z;
400  G4double CSI=CSmanager->GetCrossSection(true, Momentum, Z, N, pPDG); // CS(j,i) basic
401  return CSI;
402 }
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4