source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QHadronInelasticDataSet.cc @ 1340

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

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