Geant4: G4QHadronElasticDataSet.cc Source 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: G4QHadronElasticDataSet -- header file
00029 // Created by M. Kosov (Mikhail.Kossov@cern.ch) 21.01.10
00030 //
00031 // ----------------------------------------------------------------------
00032 // Short description: G4hadr wrapper for CHIPS elastic hA cross-sections.
00033 // ----------------------------------------------------------------------
00034 // 
00035 
00036 #include "G4QHadronElasticDataSet.hh"
00037 
00038 // Initialization of static vectors
00039 //std::vector<G4int> G4QHadronElasticDataSet::ElementZ;
00040 // Z of the element(i) in LastCalc
00041 
00042 //std::vector<std::vector<G4int>*> G4QHadronElasticDataSet::ElIsoN;
00043 // N of iso(j) of El(i)
00044 
00045 //std::vector<std::vector<G4double>*> G4QHadronElasticDataSet::IsoProbInEl;
00046 //SumProbIsoInEl
00047 
00048 
00049 G4QHadronElasticDataSet::G4QHadronElasticDataSet(const G4String& dataSetName)
00050  : G4VCrossSectionDataSet(dataSetName)
00051 {
00052   //Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
00053   Description();
00054 }
00055 
00056 void G4QHadronElasticDataSet::Description() const
00057 {
00058   char* dirName = getenv("G4PhysListDocDir");
00059   if (dirName) {
00060     std::ofstream outFile;
00061     G4String outFileName = GetName() + ".html";
00062     G4String pathName = G4String(dirName) + "/" + outFileName;
00063 
00064     outFile.open(pathName);
00065     outFile << "<html>\n";
00066     outFile << "<head>\n";
00067 
00068     outFile << "<title>Description of CHIPSElasticXS</title>\n";
00069     outFile << "</head>\n";
00070     outFile << "<body>\n";
00071 
00072     outFile << "CHIPSElasticXS provides hadron-nuclear elastic cross\n"
00073             << "sections for all hadrons at all energies.  These cross\n"
00074             << "sections represent parameterizations developed by M. Kossov.\n";
00075 
00076     outFile << "</body>\n";
00077     outFile << "</html>\n";
00078     outFile.close();
00079   }
00080 }
00081 
00082 G4bool G4QHadronElasticDataSet::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00083                                                 const G4Element*, const G4Material*)
00084 {
00085   G4ParticleDefinition* particle = Pt->GetDefinition();
00086   if      (particle ==         G4Neutron::Neutron()        ) return true; // @@ isotopes?
00087   else if (particle ==          G4Proton::Proton()         ) return true;
00088   else if (particle ==       G4PionMinus::PionMinus()      ) return true;
00089   else if (particle ==        G4PionPlus::PionPlus()       ) return true;
00090   else if (particle ==        G4KaonPlus::KaonPlus()       ) return true;
00091   else if (particle ==       G4KaonMinus::KaonMinus()      ) return true;
00092   else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) return true;
00093   else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) return true;
00094   else if (particle ==          G4Lambda::Lambda()         ) return true;
00095   else if (particle ==       G4SigmaPlus::SigmaPlus()      ) return true;
00096   else if (particle ==      G4SigmaMinus::SigmaMinus()     ) return true;
00097   else if (particle ==       G4SigmaZero::SigmaZero()      ) return true;
00098   else if (particle ==         G4XiMinus::XiMinus()        ) return true;
00099   else if (particle ==          G4XiZero::XiZero()         ) return true;
00100   else if (particle ==      G4OmegaMinus::OmegaMinus()     ) return true;
00101   else if (particle ==     G4AntiNeutron::AntiNeutron()    ) return true;
00102   else if (particle ==      G4AntiProton::AntiProton()     ) return true;
00103   else if (particle ==      G4AntiLambda::AntiLambda()     ) return true;
00104   else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) return true;
00105   else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) return true;
00106   else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) return true;
00107   else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) return true;
00108   else if (particle ==      G4AntiXiZero::AntiXiZero()     ) return true;
00109   else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) return true;
00110   return false;
00111 }
00112 
00113 // Now it is in G4VCrossSectionDataSet...
00114 /*
00115 G4double G4QHadronElasticDataSet::GetCrossSection(const G4DynamicParticle* Pt,
00116                                                     const G4Element* pElement,
00117                                                     G4double)
00118 {
00119   G4int IPIE=IsoProbInEl.size();          // How many old elements?
00120   if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
00121   {
00122     std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
00123     SPI->clear();
00124     delete SPI;
00125     std::vector<G4int>* IsN=ElIsoN[ip];   // Pointer to the N vector
00126     IsN->clear();
00127     delete IsN;
00128   }
00129   ElementZ.clear();                       // Clear the body vector for Z of Elements
00130   IsoProbInEl.clear();                    // Clear the body vector for SPI
00131   ElIsoN.clear();                         // Clear the body vector for N of Isotopes
00132   G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
00133   ElementZ.push_back(Z);                  // Remember Z of the Element
00134   G4int isoSize=0;                        // The default for the isoVectorLength is 0
00135   G4int indEl=0;                          // Index of non-trivial element or 0(default)
00136   G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
00137   if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
00138   if(isoSize)                             // The Element has non-trivial abundance set
00139   {
00140     indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
00141     if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
00142     {
00143       std::vector<std::pair<G4int,G4double>*>* newAbund =
00144                                              new std::vector<std::pair<G4int,G4double>*>;
00145       G4double* abuVector=pElement->GetRelativeAbundanceVector();
00146       for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
00147       {
00148         G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
00149         if(pElement->GetIsotope(j)->GetZ()!=Z)
00150           G4cerr<<"G4QHadronElasticDataSet::GetCrossSection"<<": Z="
00151                 <<pElement->GetIsotope(j)->GetZ()<<" # "<<Z<<G4endl;
00152         G4double abund=abuVector[j];
00153         std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00154         newAbund->push_back(pr);
00155       }
00156       indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
00157       for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
00158       delete newAbund; // Was "new" in the beginning of the name space
00159     }
00160   }
00161   std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
00162   std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
00163   IsoProbInEl.push_back(SPI);
00164   std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
00165   ElIsoN.push_back(IsN);
00166   G4int nIs=cs->size();                   // A#Of Isotopes in the Element
00167   G4double susi=0.;                       // sum of CS over isotopes
00168   if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
00169   {
00170     std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
00171     G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
00172     IsN->push_back(N);                    // Remember Min N for the Element
00173     G4double CSI=GetIsoZACrossSection(Pt,Z,Z+N,0.);//CrossSection(j,i) for the isotope
00174     curIs->second = CSI;
00175     susi+=CSI;                            // Make a sum per isotopes
00176     SPI->push_back(susi);                 // Remember summed cross-section
00177   } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
00178   return Isotopes->GetMeanCrossSection(Z,indEl); // MeanCS over isotopes
00179 }
00180 */
00181 
00182 G4double G4QHadronElasticDataSet::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int Z,
00183                                                      G4int A, const G4Isotope*,
00184                                                      const G4Element*, const G4Material*)
00185 {
00186   G4ParticleDefinition* particle = Pt->GetDefinition();
00187   G4double Momentum=Pt->GetTotalMomentum();
00188   G4VQCrossSection* CSmanager=0;
00189   G4VQCrossSection* CSmanager2=0;
00190 
00191   G4int pPDG=0;
00192   if(particle == G4Neutron::Neutron())
00193   {
00194     CSmanager=G4QNeutronElasticCrossSection::GetPointer();
00195     pPDG=2112;
00196   }
00197   else if(particle == G4Proton::Proton())
00198   {
00199     CSmanager=G4QProtonElasticCrossSection::GetPointer();
00200     pPDG=2212;
00201   }
00202   else if(particle == G4PionMinus::PionMinus())
00203   {
00204     CSmanager=G4QPionMinusElasticCrossSection::GetPointer();
00205     pPDG=-211;
00206   }
00207   else if(particle == G4PionPlus::PionPlus())
00208   {
00209     CSmanager=G4QPionPlusElasticCrossSection::GetPointer();
00210     pPDG=211;
00211   }
00212   else if(particle == G4KaonMinus::KaonMinus())
00213   {
00214     CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00215     pPDG=-321;
00216   }
00217   else if(particle == G4KaonPlus::KaonPlus())
00218   {
00219     CSmanager=G4QKaonPlusElasticCrossSection::GetPointer();
00220     pPDG=321;
00221   }
00222   else if(particle == G4KaonZeroLong::KaonZeroLong()   ||
00223           particle == G4KaonZeroShort::KaonZeroShort() ||
00224           particle == G4KaonZero::KaonZero()           ||
00225           particle == G4AntiKaonZero::AntiKaonZero()   )
00226   {
00227     CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00228     CSmanager2=G4QKaonPlusElasticCrossSection::GetPointer();
00229     if(G4UniformRand() > 0.5) pPDG= 321;
00230     else                      pPDG=-321;
00231   }
00232   else if(particle == G4Lambda::Lambda())
00233   {
00234     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00235     pPDG=3122;
00236   }
00237   else if(particle == G4SigmaPlus::SigmaPlus())
00238   {
00239     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00240     pPDG=3222;
00241   }
00242   else if(particle == G4SigmaMinus::SigmaMinus())
00243   {
00244     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00245     pPDG=3112;
00246   }
00247   else if(particle == G4SigmaZero::SigmaZero())
00248   {
00249     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00250     pPDG=3212;
00251   }
00252   else if(particle == G4XiMinus::XiMinus())
00253   {
00254     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00255     pPDG=3312;
00256   }
00257   else if(particle == G4XiZero::XiZero())
00258   {
00259     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00260     pPDG=3322;
00261   }
00262   else if(particle == G4OmegaMinus::OmegaMinus())
00263   {
00264     CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00265     pPDG=3334;
00266   }
00267   else if(particle == G4AntiNeutron::AntiNeutron())
00268   {
00269     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00270     pPDG=-2112;
00271   }
00272   else if(particle == G4AntiProton::AntiProton())
00273   {
00274     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00275     pPDG=-2212;
00276   }
00277   else if(particle == G4AntiLambda::AntiLambda())
00278   {
00279     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00280     pPDG=-3122;
00281   }
00282   else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
00283   {
00284     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00285     pPDG=-3222;
00286   }
00287   else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
00288   {
00289     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00290     pPDG=-3112;
00291   }
00292   else if(particle == G4AntiSigmaZero::AntiSigmaZero())
00293   {
00294     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00295     pPDG=-3212;
00296   }
00297   else if(particle == G4AntiXiMinus::AntiXiMinus())
00298   {
00299     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00300     pPDG=-3312;
00301   }
00302   else if(particle == G4AntiXiZero::AntiXiZero())
00303   {
00304     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00305     pPDG=-3322;
00306   }
00307   else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
00308   {
00309     CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00310     pPDG=-3334;
00311   }
00312   else
00313   {
00314     G4cerr << "-ERROR-G4QHadronElasticDataSet::GetIsoZACrossSection: PDG="
00315            << particle->GetPDGEncoding() << " isn't supported by CHIPS" << G4endl;
00316     //throw G4HadronicException(__FILE__, __LINE__,
00317     // "G4QHadronElasticDataSet::GetIsoZACrossSection: Particle isn't supported by CHIPS");
00318     return 0; 
00319   }
00320   G4int N=A-Z;
00321   G4double CSI=CSmanager->GetCrossSection(true, Momentum, Z, N, pPDG); // CS(j,i) basic
00322   if(CSmanager2) CSI = (CSI + CSmanager2->GetCrossSection(true, Momentum, Z, N, pPDG))/2;
00323   return CSI;
00324 }