source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QHadronElasticDataSet.cc @ 1316

Last change on this file since 1316 was 1316, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 14.2 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: G4QHadronElasticDataSet.cc,v 1.3 2010/05/26 12:19:06 mkossov Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// GEANT4 physics class: G4QHadronElasticDataSet -- header file
30// Created by M. Kosov (Mikhail.Kossov@cern.ch) 21.01.10
31//
32// ----------------------------------------------------------------------
33// Short description: G4hadr wrapper for CHIPS elastic hA cross-sections.
34// ----------------------------------------------------------------------
35//
36
37#include "G4QHadronElasticDataSet.hh"
38
39// Initialization of static vectors
40std::vector<G4int> G4QHadronElasticDataSet::ElementZ; // Z of the element(i) in LastCalc
41std::vector<std::vector<G4int>*> G4QHadronElasticDataSet::ElIsoN; // N of iso(j) of El(i)
42std::vector<std::vector<G4double>*> G4QHadronElasticDataSet::IsoProbInEl;//SumProbIsoInEl
43
44G4QHadronElasticDataSet::G4QHadronElasticDataSet()
45{
46  Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
47}
48
49G4bool G4QHadronElasticDataSet::IsApplicable(const G4DynamicParticle* P,const G4Element*)
50{
51  G4ParticleDefinition* particle = P->GetDefinition();
52  if      (particle ==         G4Neutron::Neutron()        ) return true; 
53  else if (particle ==          G4Proton::Proton()         ) return true;
54  else if (particle ==       G4PionMinus::PionMinus()      ) return true;
55  else if (particle ==        G4PionPlus::PionPlus()       ) return true;
56  else if (particle ==        G4KaonPlus::KaonPlus()       ) return true;
57  else if (particle ==       G4KaonMinus::KaonMinus()      ) return true;
58  else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) return true;
59  else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) return true;
60  else if (particle ==          G4Lambda::Lambda()         ) return true;
61  else if (particle ==       G4SigmaPlus::SigmaPlus()      ) return true;
62  else if (particle ==      G4SigmaMinus::SigmaMinus()     ) return true;
63  else if (particle ==       G4SigmaZero::SigmaZero()      ) return true;
64  else if (particle ==         G4XiMinus::XiMinus()        ) return true;
65  else if (particle ==          G4XiZero::XiZero()         ) return true;
66  else if (particle ==      G4OmegaMinus::OmegaMinus()     ) return true;
67  else if (particle ==     G4AntiNeutron::AntiNeutron()    ) return true;
68  else if (particle ==      G4AntiProton::AntiProton()     ) return true;
69  else if (particle ==      G4AntiLambda::AntiLambda()     ) return true;
70  else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) return true;
71  else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
72  else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) return true;
73  else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) return true;
74  else if (particle ==      G4AntiXiZero::AntiXiZero()     ) return true;
75  else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
76  return false;
77}
78
79G4bool G4QHadronElasticDataSet::IsZAApplicable(const G4DynamicParticle* Pt,
80                                                 G4double, G4double)
81{
82  G4ParticleDefinition* particle = Pt->GetDefinition();
83  if      (particle ==         G4Neutron::Neutron()        ) return true; // @@ isotopes?
84  else if (particle ==          G4Proton::Proton()         ) return true;
85  else if (particle ==       G4PionMinus::PionMinus()      ) return true;
86  else if (particle ==        G4PionPlus::PionPlus()       ) return true;
87  else if (particle ==        G4KaonPlus::KaonPlus()       ) return true;
88  else if (particle ==       G4KaonMinus::KaonMinus()      ) return true;
89  else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) return true;
90  else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) return true;
91  else if (particle ==          G4Lambda::Lambda()         ) return true;
92  else if (particle ==       G4SigmaPlus::SigmaPlus()      ) return true;
93  else if (particle ==      G4SigmaMinus::SigmaMinus()     ) return true;
94  else if (particle ==       G4SigmaZero::SigmaZero()      ) return true;
95  else if (particle ==         G4XiMinus::XiMinus()        ) return true;
96  else if (particle ==          G4XiZero::XiZero()         ) return true;
97  else if (particle ==      G4OmegaMinus::OmegaMinus()     ) return true;
98  else if (particle ==     G4AntiNeutron::AntiNeutron()    ) return true;
99  else if (particle ==      G4AntiProton::AntiProton()     ) return true;
100  else if (particle ==      G4AntiLambda::AntiLambda()     ) return true;
101  else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) return true;
102  else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
103  else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) return true;
104  else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) return true;
105  else if (particle ==      G4AntiXiZero::AntiXiZero()     ) return true;
106  else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
107  return false;
108}
109
110G4double G4QHadronElasticDataSet::GetCrossSection(const G4DynamicParticle* Pt,
111                                                    const G4Element* pElement,
112                                                    G4double)
113{
114  G4int IPIE=IsoProbInEl.size();          // How many old elements?
115  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
116  {
117    std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
118    SPI->clear();
119    delete SPI;
120    std::vector<G4int>* IsN=ElIsoN[ip];   // Pointer to the N vector
121    IsN->clear();
122    delete IsN;
123  }
124  ElementZ.clear();                       // Clear the body vector for Z of Elements
125  IsoProbInEl.clear();                    // Clear the body vector for SPI
126  ElIsoN.clear();                         // Clear the body vector for N of Isotopes
127  G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
128  ElementZ.push_back(Z);                  // Remember Z of the Element
129  G4int isoSize=0;                        // The default for the isoVectorLength is 0
130  G4int indEl=0;                          // Index of non-trivial element or 0(default)
131  G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
132  if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
133  if(isoSize)                             // The Element has non-trivial abundance set
134  {
135    indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
136    if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
137    {
138      std::vector<std::pair<G4int,G4double>*>* newAbund =
139                                             new std::vector<std::pair<G4int,G4double>*>;
140      G4double* abuVector=pElement->GetRelativeAbundanceVector();
141      for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
142      {
143        G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
144        if(pElement->GetIsotope(j)->GetZ()!=Z)
145          G4cerr<<"G4QHadronElasticDataSet::GetCrossSection"<<": Z="
146                <<pElement->GetIsotope(j)->GetZ()<<" # "<<Z<<G4endl;
147        G4double abund=abuVector[j];
148        std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
149        newAbund->push_back(pr);
150      }
151      indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
152      for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
153      delete newAbund; // Was "new" in the beginning of the name space
154    }
155  }
156  std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
157  std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
158  IsoProbInEl.push_back(SPI);
159  std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
160  ElIsoN.push_back(IsN);
161  G4int nIs=cs->size();                   // A#Of Isotopes in the Element
162  G4double susi=0.;                       // sum of CS over isotopes
163  if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
164  {
165    std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
166    G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
167    IsN->push_back(N);                    // Remember Min N for the Element
168    G4double CSI=GetIsoZACrossSection(Pt,Z,Z+N,0.);//CrossSection(j,i) for the isotope
169    curIs->second = CSI;
170    susi+=CSI;                            // Make a sum per isotopes
171    SPI->push_back(susi);                 // Remember summed cross-section
172  } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
173  return Isotopes->GetMeanCrossSection(Z,indEl); // MeanCS over isotopes
174}
175
176G4double G4QHadronElasticDataSet::GetIsoZACrossSection(const G4DynamicParticle* Pt,
177                                                         G4double Z, G4double A, G4double)
178{
179  G4ParticleDefinition* particle = Pt->GetDefinition();
180  G4double Momentum=Pt->GetTotalMomentum();
181  G4VQCrossSection* CSmanager=0;
182  G4VQCrossSection* CSmanager2=0;
183  //G4VQCrossSection* CSmanager2=0;
184  G4int pPDG=0;
185  if(particle == G4Neutron::Neutron())
186  {
187    CSmanager=G4QNeutronElasticCrossSection::GetPointer();
188    pPDG=2112;
189  }
190  else if(particle == G4Proton::Proton())
191  {
192    CSmanager=G4QProtonElasticCrossSection::GetPointer();
193    pPDG=2212;
194  }
195  else if(particle == G4PionMinus::PionMinus())
196  {
197    CSmanager=G4QPionMinusElasticCrossSection::GetPointer();
198    pPDG=-211;
199  }
200  else if(particle == G4PionPlus::PionPlus())
201  {
202    CSmanager=G4QPionPlusElasticCrossSection::GetPointer();
203    pPDG=211;
204  }
205  else if(particle == G4KaonMinus::KaonMinus())
206  {
207    CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
208    pPDG=-321;
209  }
210  else if(particle == G4KaonPlus::KaonPlus())
211  {
212    CSmanager=G4QKaonPlusElasticCrossSection::GetPointer();
213    pPDG=321;
214  }
215  else if(particle == G4KaonZeroLong::KaonZeroLong()   ||
216          particle == G4KaonZeroShort::KaonZeroShort() ||
217          particle == G4KaonZero::KaonZero()           ||
218          particle == G4AntiKaonZero::AntiKaonZero()   )
219  {
220    CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
221    CSmanager2=G4QKaonPlusElasticCrossSection::GetPointer();
222    if(G4UniformRand() > 0.5) pPDG= 321;
223    else                      pPDG=-321;
224  }
225  else if(particle == G4Lambda::Lambda())
226  {
227    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
228    pPDG=3122;
229  }
230  else if(particle == G4SigmaPlus::SigmaPlus())
231  {
232    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
233    pPDG=3222;
234  }
235  else if(particle == G4SigmaMinus::SigmaMinus())
236  {
237    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
238    pPDG=3112;
239  }
240  else if(particle == G4SigmaZero::SigmaZero())
241  {
242    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
243    pPDG=3212;
244  }
245  else if(particle == G4XiMinus::XiMinus())
246  {
247    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
248    pPDG=3312;
249  }
250  else if(particle == G4XiZero::XiZero())
251  {
252    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
253    pPDG=3322;
254  }
255  else if(particle == G4OmegaMinus::OmegaMinus())
256  {
257    CSmanager=G4QHyperonElasticCrossSection::GetPointer();
258    pPDG=3334;
259  }
260  else if(particle == G4AntiNeutron::AntiNeutron())
261  {
262    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
263    pPDG=-2112;
264  }
265  else if(particle == G4AntiProton::AntiProton())
266  {
267    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
268    pPDG=-2212;
269  }
270  else if(particle == G4AntiLambda::AntiLambda())
271  {
272    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
273    pPDG=-3122;
274  }
275  else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
276  {
277    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
278    pPDG=-3222;
279  }
280  else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
281  {
282    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
283    pPDG=-3112;
284  }
285  else if(particle == G4AntiSigmaZero::AntiSigmaZero())
286  {
287    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
288    pPDG=-3212;
289  }
290  else if(particle == G4AntiXiMinus::AntiXiMinus())
291  {
292    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
293    pPDG=-3312;
294  }
295  else if(particle == G4AntiXiZero::AntiXiZero())
296  {
297    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
298    pPDG=-3322;
299  }
300  else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
301  {
302    CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
303    pPDG=-3334;
304  }
305  else G4cout<<"-Warning-G4QHadronElasticDataSet::GetIsoZACrossSection: PDG="
306             <<particle->GetPDGEncoding()<<" isn't supported by CHIPS"<<G4endl;
307  G4int tZ=(G4int)(Z);
308  G4int tN=(G4int)(A-Z);
309  G4double CSI=CSmanager->GetCrossSection(true, Momentum, tZ, tN, pPDG); // CS(j,i) basic
310  if(CSmanager2) CSI= (CSI  +CSmanager2->GetCrossSection(true, Momentum, tZ, tN, pPDG))/2.;
311  return CSI;
312}
Note: See TracBrowser for help on using the repository browser.