Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QHadronElasticDataSet.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: G4QHadronElasticDataSet -- header file
29 // Created by M. Kosov (Mikhail.Kossov@cern.ch) 21.01.10
30 //
31 // ----------------------------------------------------------------------
32 // Short description: G4hadr wrapper for CHIPS elastic hA cross-sections.
33 // ----------------------------------------------------------------------
34 //
35 
37 
38 // Initialization of static vectors
39 //std::vector<G4int> G4QHadronElasticDataSet::ElementZ;
40 // Z of the element(i) in LastCalc
41 
42 //std::vector<std::vector<G4int>*> G4QHadronElasticDataSet::ElIsoN;
43 // N of iso(j) of El(i)
44 
45 //std::vector<std::vector<G4double>*> G4QHadronElasticDataSet::IsoProbInEl;
46 //SumProbIsoInEl
47 
48 
50  : G4VCrossSectionDataSet(dataSetName)
51 {
52  //Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
53  Description();
54 }
55 
57 {
58  char* dirName = getenv("G4PhysListDocDir");
59  if (dirName) {
60  std::ofstream outFile;
61  G4String outFileName = GetName() + ".html";
62  G4String pathName = G4String(dirName) + "/" + outFileName;
63 
64  outFile.open(pathName);
65  outFile << "<html>\n";
66  outFile << "<head>\n";
67 
68  outFile << "<title>Description of CHIPSElasticXS</title>\n";
69  outFile << "</head>\n";
70  outFile << "<body>\n";
71 
72  outFile << "CHIPSElasticXS provides hadron-nuclear elastic cross\n"
73  << "sections for all hadrons at all energies. These cross\n"
74  << "sections represent parameterizations developed by M. Kossov.\n";
75 
76  outFile << "</body>\n";
77  outFile << "</html>\n";
78  outFile.close();
79  }
80 }
81 
83  const G4Element*, const G4Material*)
84 {
85  G4ParticleDefinition* particle = Pt->GetDefinition();
86  if (particle == G4Neutron::Neutron() ) return true; // @@ isotopes?
87  else if (particle == G4Proton::Proton() ) return true;
88  else if (particle == G4PionMinus::PionMinus() ) return true;
89  else if (particle == G4PionPlus::PionPlus() ) return true;
90  else if (particle == G4KaonPlus::KaonPlus() ) return true;
91  else if (particle == G4KaonMinus::KaonMinus() ) return true;
92  else if (particle == G4KaonZeroLong::KaonZeroLong() ) return true;
93  else if (particle == G4KaonZeroShort::KaonZeroShort() ) return true;
94  else if (particle == G4Lambda::Lambda() ) return true;
95  else if (particle == G4SigmaPlus::SigmaPlus() ) return true;
96  else if (particle == G4SigmaMinus::SigmaMinus() ) return true;
97  else if (particle == G4SigmaZero::SigmaZero() ) return true;
98  else if (particle == G4XiMinus::XiMinus() ) return true;
99  else if (particle == G4XiZero::XiZero() ) return true;
100  else if (particle == G4OmegaMinus::OmegaMinus() ) return true;
101  else if (particle == G4AntiNeutron::AntiNeutron() ) return true;
102  else if (particle == G4AntiProton::AntiProton() ) return true;
103  else if (particle == G4AntiLambda::AntiLambda() ) return true;
104  else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) return true;
105  else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
106  else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) return true;
107  else if (particle == G4AntiXiMinus::AntiXiMinus() ) return true;
108  else if (particle == G4AntiXiZero::AntiXiZero() ) return true;
109  else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
110  return false;
111 }
112 
113 // Now it is in G4VCrossSectionDataSet...
114 /*
115 G4double G4QHadronElasticDataSet::GetCrossSection(const G4DynamicParticle* Pt,
116  const G4Element* pElement,
117  G4double)
118 {
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  ElementZ.clear(); // Clear the body vector for Z of Elements
130  IsoProbInEl.clear(); // Clear the body vector for SPI
131  ElIsoN.clear(); // Clear the body vector for N of Isotopes
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-trivial 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  if(isoSize) // The Element has non-trivial abundance set
139  {
140  indEl=pElement->GetIndex()+1; // Index of the non-trivial element
141  if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
142  {
143  std::vector<std::pair<G4int,G4double>*>* newAbund =
144  new std::vector<std::pair<G4int,G4double>*>;
145  G4double* abuVector=pElement->GetRelativeAbundanceVector();
146  for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
147  {
148  G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
149  if(pElement->GetIsotope(j)->GetZ()!=Z)
150  G4cerr<<"G4QHadronElasticDataSet::GetCrossSection"<<": Z="
151  <<pElement->GetIsotope(j)->GetZ()<<" # "<<Z<<G4endl;
152  G4double abund=abuVector[j];
153  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
154  newAbund->push_back(pr);
155  }
156  indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
157  for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
158  delete newAbund; // Was "new" in the beginning of the name space
159  }
160  }
161  std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
162  std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
163  IsoProbInEl.push_back(SPI);
164  std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
165  ElIsoN.push_back(IsN);
166  G4int nIs=cs->size(); // A#Of Isotopes in the Element
167  G4double susi=0.; // sum of CS over isotopes
168  if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
169  {
170  std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
171  G4int N=curIs->first; // #of Neuterons in the isotope j of El i
172  IsN->push_back(N); // Remember Min N for the Element
173  G4double CSI=GetIsoZACrossSection(Pt,Z,Z+N,0.);//CrossSection(j,i) for the isotope
174  curIs->second = CSI;
175  susi+=CSI; // Make a sum per isotopes
176  SPI->push_back(susi); // Remember summed cross-section
177  } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
178  return Isotopes->GetMeanCrossSection(Z,indEl); // MeanCS over isotopes
179 }
180 */
181 
183  G4int A, const G4Isotope*,
184  const G4Element*, const G4Material*)
185 {
186  G4ParticleDefinition* particle = Pt->GetDefinition();
187  G4double Momentum=Pt->GetTotalMomentum();
188  G4VQCrossSection* CSmanager=0;
189  G4VQCrossSection* CSmanager2=0;
190 
191  G4int pPDG=0;
192  if(particle == G4Neutron::Neutron())
193  {
195  pPDG=2112;
196  }
197  else if(particle == G4Proton::Proton())
198  {
200  pPDG=2212;
201  }
202  else if(particle == G4PionMinus::PionMinus())
203  {
205  pPDG=-211;
206  }
207  else if(particle == G4PionPlus::PionPlus())
208  {
210  pPDG=211;
211  }
212  else if(particle == G4KaonMinus::KaonMinus())
213  {
215  pPDG=-321;
216  }
217  else if(particle == G4KaonPlus::KaonPlus())
218  {
220  pPDG=321;
221  }
222  else if(particle == G4KaonZeroLong::KaonZeroLong() ||
223  particle == G4KaonZeroShort::KaonZeroShort() ||
224  particle == G4KaonZero::KaonZero() ||
225  particle == G4AntiKaonZero::AntiKaonZero() )
226  {
229  if(G4UniformRand() > 0.5) pPDG= 321;
230  else pPDG=-321;
231  }
232  else if(particle == G4Lambda::Lambda())
233  {
235  pPDG=3122;
236  }
237  else if(particle == G4SigmaPlus::SigmaPlus())
238  {
240  pPDG=3222;
241  }
242  else if(particle == G4SigmaMinus::SigmaMinus())
243  {
245  pPDG=3112;
246  }
247  else if(particle == G4SigmaZero::SigmaZero())
248  {
250  pPDG=3212;
251  }
252  else if(particle == G4XiMinus::XiMinus())
253  {
255  pPDG=3312;
256  }
257  else if(particle == G4XiZero::XiZero())
258  {
260  pPDG=3322;
261  }
262  else if(particle == G4OmegaMinus::OmegaMinus())
263  {
265  pPDG=3334;
266  }
267  else if(particle == G4AntiNeutron::AntiNeutron())
268  {
270  pPDG=-2112;
271  }
272  else if(particle == G4AntiProton::AntiProton())
273  {
275  pPDG=-2212;
276  }
277  else if(particle == G4AntiLambda::AntiLambda())
278  {
280  pPDG=-3122;
281  }
282  else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
283  {
285  pPDG=-3222;
286  }
287  else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
288  {
290  pPDG=-3112;
291  }
292  else if(particle == G4AntiSigmaZero::AntiSigmaZero())
293  {
295  pPDG=-3212;
296  }
297  else if(particle == G4AntiXiMinus::AntiXiMinus())
298  {
300  pPDG=-3312;
301  }
302  else if(particle == G4AntiXiZero::AntiXiZero())
303  {
305  pPDG=-3322;
306  }
307  else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
308  {
310  pPDG=-3334;
311  }
312  else
313  {
314  G4cerr << "-ERROR-G4QHadronElasticDataSet::GetIsoZACrossSection: PDG="
315  << particle->GetPDGEncoding() << " isn't supported by CHIPS" << G4endl;
316  //throw G4HadronicException(__FILE__, __LINE__,
317  // "G4QHadronElasticDataSet::GetIsoZACrossSection: Particle isn't supported by CHIPS");
318  return 0;
319  }
320  G4int N=A-Z;
321  G4double CSI=CSmanager->GetCrossSection(true, Momentum, Z, N, pPDG); // CS(j,i) basic
322  if(CSmanager2) CSI = (CSI + CSmanager2->GetCrossSection(true, Momentum, Z, N, pPDG))/2;
323  return CSI;
324 }