G4QHadronInelasticDataSet.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // GEANT4 physics class: G4QHadronInelasticDataSet -- header file
00029 // Created by M. Kosov (Mikhail.Kossov@cern.ch) 11.11.09
00030 //
00031 // ------------------------------------------------------------------------
00032 // Short description: G4hadr wrapper for CHIPS inelastic hA cross-sections.
00033 // ------------------------------------------------------------------------
00034 // 
00035 
00036 #include "G4QHadronInelasticDataSet.hh"
00037 #include <iostream>
00038 
00040 //std::vector<G4int> G4QHadronInelasticDataSet::ElementZ;
00042 //
00043 //std::vector<std::vector<G4int>*> G4QHadronInelasticDataSet::ElIsoN;
00045 //
00046 //std::vector<std::vector<G4double>*> G4QHadronInelasticDataSet::IsoProbInEl;
00048 
00049 G4QHadronInelasticDataSet::G4QHadronInelasticDataSet(const G4String& input_name)
00050  : G4VCrossSectionDataSet(input_name)
00051 {
00052   //CHIPSpAin    = G4QProtonNuclearCrossSection::GetPointer();
00053   //CHIPSnAin    = G4QNeutronNuclearCrossSection::GetPointer();
00054   //CHIPSpimAin  = G4QPionMinusNuclearCrossSection::GetPointer();
00055   //CHIPSpipAin  = G4QPionPlusNuclearCrossSection::GetPointer();
00056   //CHIPSkpAin   = G4QKaonPlusNuclearCrossSection::GetPointer();
00057   //CHIPSkmAin   = G4QKaonMinusNuclearCrossSection::GetPointer();
00058   //CHIPSk0Ain   = G4QKaonZeroNuclearCrossSection::GetPointer();
00059   //CHIPShAin    = G4QHyperonNuclearCrossSection::GetPointer();
00060   //CHIPShpAin   = G4QHyperonPlusNuclearCrossSection::GetPointer();
00061   //CHIPSabpAin  = G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00062   //CHIPSabAin   = G4QAntiBaryonNuclearCrossSection::GetPointer();
00073 
00074   //Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
00075   Description();
00076 }
00077 
00078 void G4QHadronInelasticDataSet::Description() const
00079 {
00080   char* dirName = getenv("G4PhysListDocDir");
00081   if (dirName) {
00082     std::ofstream outFile;
00083     G4String outFileName = GetName() + ".html";
00084     G4String pathName = G4String(dirName) + "/" + outFileName;
00085 
00086     outFile.open(pathName);
00087     outFile << "<html>\n";
00088     outFile << "<head>\n";
00089 
00090     outFile << "<title>Description of CHIPSInelasticXS</title>\n";
00091     outFile << "</head>\n";
00092     outFile << "<body>\n";
00093 
00094     outFile << "CHIPSInelasticXS provides hadron-nuclear inelastic cross\n"
00095             << "sections for all hadrons at all energies.  These cross\n"
00096             << "sections represent parameterizations developed by M. Kossov.\n";
00097 
00098     outFile << "</body>\n";
00099     outFile << "</html>\n";
00100     outFile.close();
00101   }
00102 }
00103 
00104 G4bool G4QHadronInelasticDataSet::IsIsoApplicable(const G4DynamicParticle* Pt, G4int,G4int,
00105                                                   const G4Element*, const G4Material*)
00106 {
00107   G4ParticleDefinition* particle = Pt->GetDefinition();
00108   if      (particle ==         G4Neutron::Neutron()        ) return true; // @@ isotopes?
00109   else if (particle ==          G4Proton::Proton()         ) return true;
00110   else if (particle ==       G4PionMinus::PionMinus()      ) return true;
00111   else if (particle ==        G4PionPlus::PionPlus()       ) return true;
00112   else if (particle ==        G4KaonPlus::KaonPlus()       ) return true;
00113   else if (particle ==       G4KaonMinus::KaonMinus()      ) return true;
00114   else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) return true;
00115   else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) return true;
00116   else if (particle ==          G4Lambda::Lambda()         ) return true;
00117   else if (particle ==       G4SigmaPlus::SigmaPlus()      ) return true;
00118   else if (particle ==      G4SigmaMinus::SigmaMinus()     ) return true;
00119   else if (particle ==       G4SigmaZero::SigmaZero()      ) return true;
00120   else if (particle ==         G4XiMinus::XiMinus()        ) return true;
00121   else if (particle ==          G4XiZero::XiZero()         ) return true;
00122   else if (particle ==      G4OmegaMinus::OmegaMinus()     ) return true;
00123   else if (particle ==     G4AntiNeutron::AntiNeutron()    ) return true;
00124   else if (particle ==      G4AntiProton::AntiProton()     ) return true;
00125   else if (particle ==      G4AntiLambda::AntiLambda()     ) return true;
00126   else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) return true;
00127   else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
00128   else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) return true;
00129   else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) return true;
00130   else if (particle ==      G4AntiXiZero::AntiXiZero()     ) return true;
00131   else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
00132   //else if (particle ==        G4MuonPlus::MuonPlus()       ) return true;
00133   //else if (particle ==       G4MuonMinus::MuonMinus()      ) return true;
00134   //else if (particle ==           G4Gamma::Gamma()          ) return true;
00135   //else if (particle ==        G4Electron::Electron()       ) return true;
00136   //else if (particle ==        G4Positron::Positron()       ) return true;
00137   //else if (particle ==         G4TauPlus::TauPlus()        ) return true;
00138   //else if (particle ==        G4TauMinus::TauMinus()       ) return true;
00139   //else if (particle ==   G4AntiNeutrinoE::AntiNeutrinoE()  ) return true;
00140   //else if (particle ==       G4NeutrinoE::NeutrinoE()      ) return true;
00141   //else if (particle ==  G4AntiNeutrinoMu::AntiNeutrinoMu() ) return true;
00142   //else if (particle ==      G4NeutrinoMu::NeutrinoMu()     ) return true;
00143   //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) return true;
00144   //else if (particle ==     G4NeutrinoTau::NeutrinoTau()    ) return true;
00145   return false;
00146 }
00147 
00148 // Now it is in G4VCrossSectionDataSet...
00149 /*
00150 G4double G4QHadronInelasticDataSet::GetCrossSection(const G4DynamicParticle* Pt,
00151                                                     const G4Element* pElement, G4double)
00152 {
00153   G4int IPIE=IsoProbInEl.size();          // How many old elements?
00154   if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
00155   {
00156     std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
00157     SPI->clear();
00158     delete SPI;
00159     std::vector<G4int>* IsN=ElIsoN[ip];   // Pointer to the N vector
00160     IsN->clear();
00161     delete IsN;
00162   }
00163   ElementZ.clear();                       // Clear the body vector for Z of Elements
00164   IsoProbInEl.clear();                    // Clear the body vector for SPI
00165   ElIsoN.clear();                         // Clear the body vector for N of Isotopes
00166   G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
00167   ElementZ.push_back(Z);                  // Remember Z of the Element
00168   G4int isoSize=0;                        // The default for the isoVectorLength is 0
00169   G4int indEl=0;                          // Index of non-trivial element or 0(default)
00170   G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
00171   if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
00172   if(isoSize)                             // The Element has non-trivial abundance set
00173   {
00174     indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
00175     if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
00176     {
00177       std::vector<std::pair<G4int,G4double>*>* newAbund =
00178                                              new std::vector<std::pair<G4int,G4double>*>;
00179       G4double* abuVector=pElement->GetRelativeAbundanceVector();
00180       for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
00181       {
00182         G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
00183         if(pElement->GetIsotope(j)->GetZ()!=Z)
00184           G4cerr<<"G4QHadronInelasticDataSet::GetCrossSection"<<": Z="
00185                 <<pElement->GetIsotope(j)->GetZ()<<" # "<<Z<<G4endl;
00186         G4double abund=abuVector[j];
00187         std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00188         newAbund->push_back(pr);
00189       }
00190       indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
00191       for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
00192       delete newAbund; // Was "new" in the beginning of the name space
00193     }
00194   }
00195   std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
00196   std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
00197   IsoProbInEl.push_back(SPI);
00198   std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
00199   ElIsoN.push_back(IsN);
00200   G4int nIs=cs->size();                   // A#Of Isotopes in the Element
00201   G4double susi=0.;                       // sum of CS over isotopes
00202   if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
00203   {
00204     std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
00205     G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
00206     IsN->push_back(N);                    // Remember Min N for the Element
00207     G4double CSI=GetIsoZACrossSection(Pt,Z,Z+N,0.);//CrossSection(j,i) for the isotope
00208     curIs->second = CSI;
00209     susi+=CSI;                            // Make a sum per isotopes
00210     SPI->push_back(susi);                 // Remember summed cross-section
00211   } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
00212   return Isotopes->GetMeanCrossSection(Z,indEl); // MeanCS over isotopes
00213 }
00214 */
00215 
00216 G4double G4QHadronInelasticDataSet::GetIsoCrossSection(const G4DynamicParticle* Pt,G4int Z,
00217                                                        G4int A, const G4Isotope*,
00218                                                        const G4Element*, const G4Material*)
00219 {
00220   G4ParticleDefinition* particle = Pt->GetDefinition();
00221   G4double Momentum=Pt->GetTotalMomentum();
00222   G4VQCrossSection* CSmanager=0;
00223 
00224   G4int pPDG=0;
00225   if(particle == G4Neutron::Neutron())
00226   {
00227     CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
00228     pPDG=2112;
00229   }
00230   else if(particle == G4Proton::Proton())
00231   {
00232     CSmanager=G4QProtonNuclearCrossSection::GetPointer();
00233     pPDG=2212;
00234   }
00235   else if(particle == G4PionMinus::PionMinus())
00236   {
00237     CSmanager=G4QPionMinusNuclearCrossSection::GetPointer();
00238     pPDG=-211;
00239   }
00240   else if(particle == G4PionPlus::PionPlus())
00241   {
00242     CSmanager=G4QPionPlusNuclearCrossSection::GetPointer();
00243     pPDG=211;
00244   }
00245   else if(particle == G4KaonMinus::KaonMinus())
00246   {
00247     CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer();
00248     pPDG=-321;
00249   }
00250   else if(particle == G4KaonPlus::KaonPlus())
00251   {
00252     CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer();
00253     pPDG=321;
00254   }
00255   else if(particle == G4KaonZeroLong::KaonZeroLong()   ||
00256           particle == G4KaonZeroShort::KaonZeroShort() ||
00257           particle == G4KaonZero::KaonZero()           ||
00258           particle == G4AntiKaonZero::AntiKaonZero()   )
00259   {
00260     CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer();
00261     if(G4UniformRand() > 0.5) pPDG= 311;
00262     else                      pPDG=-311;
00263   }
00264   else if(particle == G4Lambda::Lambda())
00265   {
00266     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00267     pPDG=3122;
00268   }
00269   else if(particle == G4SigmaPlus::SigmaPlus())
00270   {
00271     CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer();
00272     pPDG=3222;
00273   }
00274   else if(particle == G4SigmaMinus::SigmaMinus())
00275   {
00276     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00277     pPDG=3112;
00278   }
00279   else if(particle == G4SigmaZero::SigmaZero())
00280   {
00281     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00282     pPDG=3212;
00283   }
00284   else if(particle == G4XiMinus::XiMinus())
00285   {
00286     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00287     pPDG=3312;
00288   }
00289   else if(particle == G4XiZero::XiZero())
00290   {
00291     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00292     pPDG=3322;
00293   }
00294   else if(particle == G4OmegaMinus::OmegaMinus())
00295   {
00296     CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
00297     pPDG=3334;
00298   }
00299   else if(particle == G4AntiNeutron::AntiNeutron())
00300   {
00301     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00302     pPDG=-2112;
00303   }
00304   else if(particle == G4AntiProton::AntiProton())
00305   {
00306     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00307     pPDG=-2212;
00308   }
00309   else if(particle == G4AntiLambda::AntiLambda())
00310   {
00311     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00312     pPDG=-3122;
00313   }
00314   else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
00315   {
00316     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00317     pPDG=-3222;
00318   }
00319   else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
00320   {
00321     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00322     pPDG=-3112;
00323   }
00324   else if(particle == G4AntiSigmaZero::AntiSigmaZero())
00325   {
00326     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00327     pPDG=-3212;
00328   }
00329   else if(particle == G4AntiXiMinus::AntiXiMinus())
00330   {
00331     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00332     pPDG=-3312;
00333   }
00334   else if(particle == G4AntiXiZero::AntiXiZero())
00335   {
00336     CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
00337     pPDG=-3322;
00338   }
00339   else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
00340   {
00341     CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
00342     pPDG=-3334;
00343   }
00344   //else if(particle == G4Gamma::Gamma())
00345   //{
00346   //  CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
00347   //  pPDG=22;
00348   //}
00349   //else if(particle == G4Electron::Electron() ||
00350   //        particle == G4Positron::Positron())
00351   //{
00352   //  CSmanager=G4QElectronNuclearCrossSection::GetPointer();
00353   //  pPDG=11;
00354   //}
00355   //else if(particle == G4MuonPlus::MuonPlus() ||
00356   //        particle == G4MuonMinus::MuonMinus())
00357   //{
00358   //  CSmanager=G4QMuonNuclearCrossSection::GetPointer();
00359   //  pPDG=13;
00360   //}
00361   //else if(particle == G4TauPlus::TauPlus() ||
00362   //        particle == G4TauMinus::TauMinus())
00363   //{
00364   //  CSmanager=G4QTauNuclearCrossSection::GetPointer();
00365   //  pPDG=15;
00366   //}
00367   //else if(particle == G4NeutrinoMu::NeutrinoMu() )
00368   //{
00369   //  CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
00370   //  CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
00371   //  pPDG=14;
00372   //}
00373   //else if(particle == G4AntiNeutrinoMu::AntiNeutrinoMu() )
00374   //{
00375   //  CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
00376   //  CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
00377   //  pPDG=-14;
00378   //}
00379   //else if(particle == G4NeutrinoE::NeutrinoE() )
00380   //{
00381   //  CSmanager=G4QNuENuclearCrossSection::GetPointer();
00382   //  CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
00383   //  pPDG=12;
00384   //}
00385   //else if(particle == G4AntiNeutrinoE::AntiNeutrinoE() )
00386   //{
00387   //  CSmanager=G4QANuENuclearCrossSection::GetPointer();
00388   //  CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
00389   //  pPDG=-12;
00390   //}
00391   else
00392   {
00393     G4cerr << "-ERROR-G4QHadronInelasticDataSet::GetIsoZACrossSection: PDG="
00394            << particle->GetPDGEncoding() << " isn't supported by CHIPS" << G4endl;
00395     //throw G4HadronicException(__FILE__, __LINE__,
00396     //"G4QHadronInelasticDataSet::GetIsoZACrossSection: Particle not supported by CHIPS");
00397     return 0; 
00398   }
00399   G4int N=A-Z;
00400   G4double CSI=CSmanager->GetCrossSection(true, Momentum, Z, N, pPDG); // CS(j,i) basic
00401   return CSI;
00402 }

Generated on Mon May 27 17:49:35 2013 for Geant4 by  doxygen 1.4.7