diff --git a/Documentation/FinalPresentation/Presentation.pdf b/Documentation/FinalPresentation/Presentation.pdf index 3d3b65c3..dd003911 100644 Binary files a/Documentation/FinalPresentation/Presentation.pdf and b/Documentation/FinalPresentation/Presentation.pdf differ diff --git a/Documentation/FinalPresentation/Presentation.tex b/Documentation/FinalPresentation/Presentation.tex index ecf5f86d..8f8917fd 100644 --- a/Documentation/FinalPresentation/Presentation.tex +++ b/Documentation/FinalPresentation/Presentation.tex @@ -1,14 +1,11 @@ \documentclass{beamer} \usepackage[utf8]{inputenc} -\setbeamertemplate{navigation symbols}{} - -% 8/10 -% \usetheme{Antibes} - -% 8.1 -% \usetheme{Copenhagen} +\usepackage{graphicx} +\usepackage{textcomp} +\newcommand{\textapprox}{\raisebox{0.5ex}{\texttildelow}} +\setbeamertemplate{navigation symbols}{} \usetheme{Luebeck} \newcommand\pro{\item[$+$]} @@ -19,10 +16,10 @@ \title[GEANT4-GPU (McMaster University)]{Running GEANT4 Functions on a GPU} \subtitle{Discussion of Results} \institute{McMaster University} -\author[S. Douglas, M .Pagnan, R. Gorrie, V. Reginato]{ +\author[S. Douglas, R. Gorrie, M .Pagnan, V. Reginato]{ Stuart Douglas -- dougls2 -\\Matthew Pagnan -- pagnanmm \\Rob Gorrie -- gorrierw +\\Matthew Pagnan -- pagnanmm \\Victor Reginato -- reginavp } @@ -40,6 +37,10 @@ \section{Introduction} \subsection{Brief Project Overview} \begin{frame} \frametitle{Brief Project Overview} +Take an existing particle simulation toolkit - GEANT4 - and have some functions run on a GPU device to improve performance. +\begin{block}{Definition: GEANT4} +GEANT4 is +\end{block} \end{frame} \begin{frame} @@ -47,12 +48,23 @@ \subsection{Brief Project Overview} \end{frame} \subsection{Explanation of Terms} +%Victor \begin{frame} \frametitle{What is GEANT4} +\begin{itemize} +\item Geant4 is a toolkit that is meant to simulate the passage of particles through matter. +\item It has been developed over the years through collaborative effort of many different institutions and individuals. +\item Geant4 has many different applications, including applications in high energy physics, space and radiation, medical. +\end{itemize} \end{frame} \begin{frame} +\begin{itemize} \frametitle{What is GP-GPU} +\item General purpose graphic processing unit computing is a re-purposing of graphics hardware +\item Allows GPUs to perform computations that would typically be computed on the CPU +\item If problems are suitable to mass parallelization +\end{itemize} \end{frame} \subsection{Scope} @@ -80,7 +92,7 @@ \section{Discussion} \end{itemize} \end{frame} -\subsection{Completely on GPU} +\subsection{Entire G4ParticleHPVector Object on GPU} \begin{frame} \frametitle{Completely on GPU} \begin{itemize} @@ -127,14 +139,29 @@ \subsection{Implementation 2} \end{itemize} \end{frame} -\subsubsection{Intensive Functions on GPU} +\subsubsection{Performance} \begin{frame} -\frametitle{Intensive Functions on GPU} +\frametitle{Performance Summary} +\begin{itemize} +\item Most functions slower on GPU until \textapprox 10,000 entries +\item Most \emph{commonly-used} functions significantly slower on GPU +\begin{itemize} +\item Lots of data accesses +\end{itemize} +\item Many problems in vector class not well-suited to parallelism +\end{itemize} \end{frame} -\subsubsection{Performance} \begin{frame} -\frametitle{Performance Results} +\frametitle{Performance Results - \texttt{Times}} +\begin{itemize} +\item Multiplies each point in vector by factor +\end{itemize} +\begin{figure} +\centering +\includegraphics[width=0.9\textwidth]{images/times_line.png} +\caption{Runtime vs. Number of Data Points -- \texttt{Times}} +\end{figure} \end{frame} \begin{frame} diff --git a/Documentation/FinalPresentation/images/get15_bar.png b/Documentation/FinalPresentation/images/get15_bar.png new file mode 100644 index 00000000..5e9965a0 Binary files /dev/null and b/Documentation/FinalPresentation/images/get15_bar.png differ diff --git a/Documentation/FinalPresentation/images/get15_line.png b/Documentation/FinalPresentation/images/get15_line.png new file mode 100644 index 00000000..054865cb Binary files /dev/null and b/Documentation/FinalPresentation/images/get15_line.png differ diff --git a/Documentation/FinalPresentation/images/get50_bar.png b/Documentation/FinalPresentation/images/get50_bar.png new file mode 100644 index 00000000..22ec604a Binary files /dev/null and b/Documentation/FinalPresentation/images/get50_bar.png differ diff --git a/Documentation/FinalPresentation/images/get50_line.png b/Documentation/FinalPresentation/images/get50_line.png new file mode 100644 index 00000000..4bef5f43 Binary files /dev/null and b/Documentation/FinalPresentation/images/get50_line.png differ diff --git a/Documentation/FinalPresentation/images/getxsec_e_bar.png b/Documentation/FinalPresentation/images/getxsec_e_bar.png new file mode 100644 index 00000000..288a9061 Binary files /dev/null and b/Documentation/FinalPresentation/images/getxsec_e_bar.png differ diff --git a/Documentation/FinalPresentation/images/getxsec_e_line.png b/Documentation/FinalPresentation/images/getxsec_e_line.png new file mode 100644 index 00000000..a0a9972f Binary files /dev/null and b/Documentation/FinalPresentation/images/getxsec_e_line.png differ diff --git a/Documentation/FinalPresentation/images/getxsec_e_min_bar.png b/Documentation/FinalPresentation/images/getxsec_e_min_bar.png new file mode 100644 index 00000000..b6067acd Binary files /dev/null and b/Documentation/FinalPresentation/images/getxsec_e_min_bar.png differ diff --git a/Documentation/FinalPresentation/images/getxsec_e_min_line.png b/Documentation/FinalPresentation/images/getxsec_e_min_line.png new file mode 100644 index 00000000..6ae5772b Binary files /dev/null and b/Documentation/FinalPresentation/images/getxsec_e_min_line.png differ diff --git a/Documentation/FinalPresentation/images/init_bar.png b/Documentation/FinalPresentation/images/init_bar.png new file mode 100644 index 00000000..7b04c554 Binary files /dev/null and b/Documentation/FinalPresentation/images/init_bar.png differ diff --git a/Documentation/FinalPresentation/images/init_line.png b/Documentation/FinalPresentation/images/init_line.png new file mode 100644 index 00000000..8a88c242 Binary files /dev/null and b/Documentation/FinalPresentation/images/init_line.png differ diff --git a/Documentation/FinalPresentation/images/samplelin_bar.png b/Documentation/FinalPresentation/images/samplelin_bar.png new file mode 100644 index 00000000..dcf6221f Binary files /dev/null and b/Documentation/FinalPresentation/images/samplelin_bar.png differ diff --git a/Documentation/FinalPresentation/images/samplelin_line.png b/Documentation/FinalPresentation/images/samplelin_line.png new file mode 100644 index 00000000..b94c5be6 Binary files /dev/null and b/Documentation/FinalPresentation/images/samplelin_line.png differ diff --git a/Documentation/FinalPresentation/images/times_bar.png b/Documentation/FinalPresentation/images/times_bar.png new file mode 100644 index 00000000..5f16b66c Binary files /dev/null and b/Documentation/FinalPresentation/images/times_bar.png differ diff --git a/Documentation/FinalPresentation/images/times_line.png b/Documentation/FinalPresentation/images/times_line.png new file mode 100644 index 00000000..c47eafc1 Binary files /dev/null and b/Documentation/FinalPresentation/images/times_line.png differ diff --git a/Examples/Hadr04-Uranium/src/src/DetectorConstruction.cc b/Examples/Hadr04-Uranium/src/src/DetectorConstruction.cc index 7478414c..4a33edeb 100644 --- a/Examples/Hadr04-Uranium/src/src/DetectorConstruction.cc +++ b/Examples/Hadr04-Uranium/src/src/DetectorConstruction.cc @@ -86,24 +86,24 @@ G4NistManager* man = G4NistManager::Instance(); // specific element name for thermal neutronHP // (see G4NeutronHPThermalScatteringNames.cc) - G4Element* H = new G4Element("TS_H_of_Water" ,"H" , 1., 1.0079*g/mole); - G4Element* O = new G4Element("Oxygen" ,"O" , 8., 16.00*g/mole); - - G4Material* H2O = - new G4Material("Water_ts", 1.000*g/cm3, ncomponents=2, - kStateLiquid, 593*kelvin, 150*bar); - H2O->AddElement(H, natoms=2); - H2O->AddElement(O, natoms=1); - H2O->GetIonisation()->SetMeanExcitationEnergy(78.0*eV); - - // graphite - G4Isotope* C12 = new G4Isotope("C12", 6, 12); - G4Element* elC = new G4Element("TS_C_of_Graphite","C", ncomponents=1); - elC->AddIsotope(C12, 100.*perCent); - G4Material* graphite = - new G4Material("graphite", 2.27*g/cm3, ncomponents=1, - kStateSolid, 293*kelvin, 1*atmosphere); - graphite->AddElement(elC, natoms=1); + // G4Element* H = new G4Element("TS_H_of_Water" ,"H" , 1., 1.0079*g/mole); + // G4Element* O = new G4Element("Oxygen" ,"O" , 8., 16.00*g/mole); + + // G4Material* H2O = + // new G4Material("Water_ts", 1.000*g/cm3, ncomponents=2, + // kStateLiquid, 593*kelvin, 150*bar); + // H2O->AddElement(H, natoms=2); + // H2O->AddElement(O, natoms=1); + // H2O->GetIonisation()->SetMeanExcitationEnergy(78.0*eV); + + // // graphite + // G4Isotope* C12 = new G4Isotope("C12", 6, 12); + // G4Element* elC = new G4Element("TS_C_of_Graphite","C", ncomponents=1); + // elC->AddIsotope(C12, 100.*perCent); + // G4Material* graphite = + // new G4Material("graphite", 2.27*g/cm3, ncomponents=1, + // kStateSolid, 293*kelvin, 1*atmosphere); + // graphite->AddElement(elC, natoms=1); ///G4cout << *(G4Material::GetMaterialTable()) << G4endl; diff --git a/geant4.10.02/source/externals/cuda/include/G4ParticleHPVector.hh b/geant4.10.02/source/externals/cuda/include/G4ParticleHPVector.hh index ddca5d31..709d1aac 100644 --- a/geant4.10.02/source/externals/cuda/include/G4ParticleHPVector.hh +++ b/geant4.10.02/source/externals/cuda/include/G4ParticleHPVector.hh @@ -73,7 +73,7 @@ class G4ParticleHPVector G4ParticleHPVector & operator = (const G4ParticleHPVector & right); G4double GetXsec(G4double e); - void GetXsecBuffer(G4double * queryList, G4int length); + void GetXsecList(G4double * queryList, G4int length); void Dump(); void ThinOut(G4double precision); void Merge(G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector * active, G4ParticleHPVector * passive); diff --git a/geant4.10.02/source/processes/hadronic/models/particle_hp/include/G4ParticleHPVector.hh b/geant4.10.02/source/processes/hadronic/models/particle_hp/include/G4ParticleHPVector.hh index ba5edee3..41da8528 100644 --- a/geant4.10.02/source/processes/hadronic/models/particle_hp/include/G4ParticleHPVector.hh +++ b/geant4.10.02/source/processes/hadronic/models/particle_hp/include/G4ParticleHPVector.hh @@ -49,8 +49,7 @@ #include #include -#define GEANT4_ENABLE_CUDA 0 -#if GEANT4_ENABLE_CUDA +#if 1 #include "/Users/stuart/Documents/4th_Year/CS_4ZP6/GEANT4-GPU/geant4.10.02/source/externals/cuda/include/G4ParticleHPVector_CUDA.hh" #endif @@ -62,17 +61,19 @@ class G4ParticleHPVector { friend G4ParticleHPVector & operator + (G4ParticleHPVector & left, G4ParticleHPVector & right); + /****************************************** * PUBLIC ******************************************/ public: + G4ParticleHPVector(); G4ParticleHPVector(G4int n); ~G4ParticleHPVector(); G4ParticleHPVector & operator = (const G4ParticleHPVector & right); G4double GetXsec(G4double e); - void GetXsecList(G4double * energiesIn_xSecsOut, G4int length); + void GetXsecList(G4double * queryList, G4int length); void Dump(); void ThinOut(G4double precision); void Merge(G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector * active, G4ParticleHPVector * passive); @@ -85,171 +86,228 @@ class G4ParticleHPVector ******************************************/ inline const G4ParticleHPDataPoint & GetPoint(G4int i) const { - if (i < 0) { - i = 0; - } else if (i >= GetVectorLength()) { - i = GetVectorLength()-1; - } - return theData[i]; + #if GEANT4_ENABLE_CUDA + return cudaVector->GetPoint(i); + #else + if (i < 0) { + i=0; + } else if(i >= GetVectorLength()) { + i=GetVectorLength()-1; + } + return theData[i]; + #endif } inline G4int GetVectorLength() const { - return nEntries; + #if GEANT4_ENABLE_CUDA + return cudaVector->GetVectorLength(); + #else + return nEntries; + #endif } inline G4double GetEnergy(G4int i) const { - if (i < 0) { - i = 0; - } - if(i >= GetVectorLength()) { - i = GetVectorLength()-1; - } - return theData[i].GetX(); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetX(i); + #else + if (i<0) { + i=0; + } + if(i>=GetVectorLength()) { + i=GetVectorLength()-1; + } + return theData[i].GetX(); + #endif } inline G4double GetX(G4int i) const { - if (i < 0) { - i = 0; - } - if (i >= GetVectorLength()) { - i = GetVectorLength()-1; - } - return theData[i].GetX(); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetX(i); + #else + if (i<0) { + i=0; + } + if(i>=GetVectorLength()) { + i=GetVectorLength()-1; + } + return theData[i].GetX(); + #endif } inline G4double GetXsec(G4int i) { - if (i < 0) { - i = 0; - } - if (i >= GetVectorLength()) { - i = GetVectorLength()-1; - } - return theData[i].GetY(); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetY(i); + #else + if (i<0) { + i=0; + } + if(i>=GetVectorLength()) { + i=GetVectorLength()-1; + } + return theData[i].GetY(); + #endif } + // ---- Can put on GPU G4double GetXsec(G4double e, G4int min) - { - if (GetVectorLength() <= 0) { - return 0.0; - } else if (min >= nEntries) { - return theData[0].GetY(); - } - - min = (min >= 0) ? min : 0; - G4int i; - for(i = min; i < nEntries; i++) - { - if (theData[i].GetX() > e) { - break; - } - } - G4int low = i-1; - G4int high = i; - if (i == 0) - { - low = 0; - high = 1; - } - else if (i == nEntries) - { - low = nEntries-2; - high = nEntries-1; - } - G4double y; - if (e < theData[nEntries-1].GetX()) - { - // Protect against doubled-up x values - // Note: added std::abs check and check that theData[high] has non-zero energy to match - // GetXSec(e), as defined in G4ParticleHPVector.cc - Geant4-GPU team - if (theData[high].GetX() != 0 - && (std::abs((theData[high].GetX()-theData[low].GetX())/theData[high].GetX()) < 0.000001)) + { + #if GEANT4_ENABLE_CUDA + return cudaVector->GetXsec(e, min); + #else + if (GetVectorLength() <= 0) { + return 0.0; + } else if (min >= nEntries) { + return theData[0].GetY(); + } + + min = (min >= 0) ? min : 0; + G4int i; + for(i=min ; ie) { + break; + } + } + G4int low = i-1; + G4int high = i; + if(i==0) { - y = theData[low].GetY(); + low = 0; + high = 1; + } + else if(i==nEntries) + { + low = nEntries-2; + high = nEntries-1; + } + G4double y; + if(eGetXsec(x); + #else + return GetXsec(x); + #endif } inline G4double GetY(G4int i) { - if (i < 0) { - i = 0; - } - if(i >= GetVectorLength()){ - i = GetVectorLength()-1; - } - return theData[i].GetY(); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetY(i); + #else + if (i<0) { + i=0; + } + if(i>=GetVectorLength()){ + i=GetVectorLength()-1; + } + return theData[i].GetY(); + #endif } inline G4double GetY(G4int i) const { - if (i < 0) { - i = 0; - } - if(i >= GetVectorLength()){ - i = GetVectorLength()-1; - } - return theData[i].GetY(); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetY(i); + #else + if (i<0) { + i=0; + } + if(i>=GetVectorLength()){ + i=GetVectorLength()-1; + } + return theData[i].GetY(); + #endif } - + // Should be able to run this on GPU G4double GetMeanX() { - G4double result; - G4double running = 0; - G4double weighted = 0; - for(G4int i = 1; i < nEntries; i++) - { - running += theInt.GetBinIntegral(theManager.GetScheme(i-1), - theData[i-1].GetX(), theData[i].GetX(), - theData[i-1].GetY(), theData[i].GetY()); - weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1), - theData[i-1].GetX(), theData[i].GetX(), - theData[i-1].GetY(), theData[i].GetY()); - } - result = weighted / running; - return result; + #if GEANT4_ENABLE_CUDA + return cudaVector->GetMeanX(); + #else + G4double result; + G4double running = 0; + G4double weighted = 0; + for(G4int i=1; iGetLabel(); + #else + return label; + #endif } inline G4double GetIntegral() // linear interpolation; use with care { - if (totalIntegral < -0.5) { - Integrate(); - } - return totalIntegral; + #if GEANT4_ENABLE_CUDA + return cudaVector->GetIntegral(); + #else + if(totalIntegral<-0.5) { + Integrate(); + } + return totalIntegral; + #endif } inline const G4InterpolationManager & GetInterpolationManager() const { - return theManager; + #if GEANT4_ENABLE_CUDA + return cudaVector->GetInterpolationManager(); + #else + return theManager; + #endif } inline G4InterpolationScheme GetScheme(G4int anIndex) { - return theManager.GetScheme(anIndex); + #if GEANT4_ENABLE_CUDA + return cudaVector->GetScheme(anIndex); + #else + return theManager.GetScheme(anIndex); + #endif } @@ -258,79 +316,115 @@ class G4ParticleHPVector ******************************************/ inline void SetVerbose(G4int ff) { - Verbose = ff; + #if GEANT4_ENABLE_CUDA + cudaVector->SetVerbose(ff); + #else + Verbose = ff; + #endif } inline void SetPoint(G4int i, const G4ParticleHPDataPoint & it) { - G4double x = it.GetX(); - G4double y = it.GetY(); - SetData(i, x, y); + #if GEANT4_ENABLE_CUDA + cudaVector->SetPoint(i, it); + #else + G4double x = it.GetX(); + G4double y = it.GetY(); + SetData(i, x, y); + #endif } inline void SetData(G4int i, G4double x, G4double y) { - Check(i); - if (y > maxValue) { - maxValue = y; - } - theData[i].SetData(x, y); + #if GEANT4_ENABLE_CUDA + cudaVector->SetData(i,x,y); + #else + Check(i); + if(y>maxValue) { + maxValue=y; + } + theData[i].SetData(x, y); + #endif } inline void SetX(G4int i, G4double e) { - Check(i); - theData[i].SetX(e); + #if GEANT4_ENABLE_CUDA + cudaVector->SetX(i,e); + #else + Check(i); + theData[i].SetX(e); + #endif } inline void SetEnergy(G4int i, G4double e) { - Check(i); - theData[i].SetX(e); + #if GEANT4_ENABLE_CUDA + cudaVector->SetEnergy(i,e); + #else + Check(i); + theData[i].SetX(e); + #endif } inline void SetY(G4int i, G4double x) { - Check(i); - if (x > maxValue) { - maxValue = x; - } - theData[i].SetY(x); + #if GEANT4_ENABLE_CUDA + cudaVector->SetY(i,x); + #else + Check(i); + if(x>maxValue) maxValue=x; + theData[i].SetY(x); + #endif } inline void SetXsec(G4int i, G4double x) { - Check(i); - if (x > maxValue) { - maxValue = x; - } - theData[i].SetY(x); + #if GEANT4_ENABLE_CUDA + cudaVector->SetXsec(i,x); + #else + Check(i); + if(x>maxValue) { + maxValue=x; + } + theData[i].SetY(x); + #endif } inline void SetLabel(G4double aLabel) { - label = aLabel; + #if GEANT4_ENABLE_CUDA + cudaVector->SetLabel(aLabel); + #else + label = aLabel; + #endif } inline void SetInterpolationManager(const G4InterpolationManager & aManager) { #if GEANT4_ENABLE_CUDA cudaVector->SetInterpolationManager(aManager); - #endif - theManager = aManager; + #else + theManager = aManager; + #endif } inline void SetInterpolationManager(G4InterpolationManager & aMan) { #if GEANT4_ENABLE_CUDA cudaVector->SetInterpolationManager(aMan); + #else + theManager = aMan; #endif - theManager = aMan; } inline void SetScheme(G4int aPoint, const G4InterpolationScheme & aScheme) { - theManager.AppendScheme(aPoint, aScheme); + #if GEANT4_ENABLE_CUDA + cudaVector->SetScheme(aPoint, aScheme); + #else + theManager.AppendScheme(aPoint, aScheme); + #endif } @@ -339,287 +433,334 @@ class G4ParticleHPVector ******************************************/ void Init(std::istream & aDataFile, G4int total, G4double ux=1., G4double uy=1.) { - G4double x,y; - for (G4int i = 0; i < total; i++) - { - aDataFile >> x >> y; - x *= ux; - y *= uy; - SetData(i,x,y); - if (0 == nEntries % 10) + #if GEANT4_ENABLE_CUDA + cudaVector->Init(aDataFile, total, ux, uy); + #else + G4double x,y; + for (G4int i=0;i> x >> y; + x*=ux; + y*=uy; + SetData(i,x,y); + if(0 == nEntries%10) + { + theHash.SetData(nEntries-1, x, y); + } } - } + #endif } void Init(std::istream & aDataFile,G4double ux=1., G4double uy=1.) { - G4int total; - aDataFile >> total; - if (theData != 0) { - delete [] theData; - } - theData = new G4ParticleHPDataPoint[total]; - nPoints = total; - nEntries = 0; - theManager.Init(aDataFile); - Init(aDataFile, total, ux, uy); + #if GEANT4_ENABLE_CUDA + cudaVector->Init(aDataFile, ux, uy); + #else + G4int total; + aDataFile >> total; + if(theData!=0) { + delete [] theData; + } + theData = new G4ParticleHPDataPoint[total]; + nPoints=total; + nEntries=0; + theManager.Init(aDataFile); + Init(aDataFile, total, ux, uy); + #endif } inline void InitInterpolation(std::istream & aDataFile) { - theManager.Init(aDataFile); + #if GEANT4_ENABLE_CUDA + cudaVector->InitInterpolation(aDataFile); + #else + theManager.Init(aDataFile); + #endif } void Hash() { - G4int i; - G4double x, y; - for(i = 0; i < nEntries; i++) - { - if (0 == (i+1) % 10) + #if GEANT4_ENABLE_CUDA + // no hashing on cuda + #else + G4int i; + G4double x, y; + for(i=0 ; iCleanUp(); + #else + nEntries=0; + theManager.CleanUp(); + maxValue = -DBL_MAX; + theHash.Clear(); + //080811 TK DB + delete[] theIntegral; + theIntegral = NULL; + #endif } // merges the vectors active and passive into *this inline void Merge(G4ParticleHPVector * active, G4ParticleHPVector * passive) { - CleanUp(); - G4int s_tmp = 0, n=0, m_tmp=0; - G4ParticleHPVector * tmp; - G4int a = s_tmp, p = n, t; - while (a < active->GetVectorLength() && p < passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi - { - if (active->GetEnergy(a) <= passive->GetEnergy(p)) + #if GEANT4_ENABLE_CUDA + cudaVector->Merge(&(*(active)->cudaVector), &(*(passive)->cudaVector)); + #else + CleanUp(); + G4int s_tmp = 0, n=0, m_tmp=0; + G4ParticleHPVector * tmp; + G4int a = s_tmp, p = n, t; + while (aGetVectorLength()&&pGetVectorLength()) // Loop checking, 11.05.2015, T. Koi { - G4double xa = active->GetEnergy(a); - G4double yy = active->GetXsec(a); - SetData(m_tmp, xa, yy); - theManager.AppendScheme(m_tmp, active->GetScheme(a)); - m_tmp++; - a++; - G4double xp = passive->GetEnergy(p); - - //080409 TKDB - //if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++; - if (!(xa == 0) && std::abs(std::abs(xp-xa)/xa) < 0.001) { - p++; + if(active->GetEnergy(a) <= passive->GetEnergy(p)) + { + G4double xa = active->GetEnergy(a); + G4double yy = active->GetXsec(a); + SetData(m_tmp, xa, yy); + theManager.AppendScheme(m_tmp, active->GetScheme(a)); + m_tmp++; + a++; + G4double xp = passive->GetEnergy(p); + + //080409 TKDB + //if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++; + if ( !( xa == 0 ) && std::abs(std::abs(xp-xa)/xa)<0.001 ) p++; + } else { + tmp = active; + t=a; + active = passive; + a=p; + passive = tmp; + p=t; } - } else { - tmp = active; - t=a; - active = passive; - a=p; - passive = tmp; - p=t; - } - } - while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi - { - SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a)); - theManager.AppendScheme(m_tmp++, active->GetScheme(a)); - a++; - } - while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi - { - if(std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001) - //if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001) + } + while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi { - SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p)); - theManager.AppendScheme(m_tmp++, active->GetScheme(p)); + SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a)); + theManager.AppendScheme(m_tmp++, active->GetScheme(a)); + a++; } - p++; - } + while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi + { + if(std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001) + //if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001) + { + SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p)); + theManager.AppendScheme(m_tmp++, active->GetScheme(p)); + } + p++; + } + #endif } G4double SampleLin() // Samples X according to distribution Y, linear int { - G4double result; - if(theIntegral==0) - { - IntegrateAndNormalise(); - } - - if (GetVectorLength() == 0) - { - result = 0; - } - else if (GetVectorLength()==1) - { - result = theData[0].GetX(); - } - else - { - G4int i; - G4double rand = G4UniformRand(); - - // this was replaced - // for(i=1;i= 0 ; i--) + #if GEANT4_ENABLE_CUDA + return cudaVector->SampleLin(); + #else + G4double result; + if(theIntegral==0) { - if (rand > theIntegral[i]/theIntegral[GetVectorLength()-1]) { - break; - } + IntegrateAndNormalise(); } - if (i != GetVectorLength()-1) { - i++; + + if (GetVectorLength() == 0) + { + result = 0; } - // until this (end) - - G4double x1, x2, y1, y2; - y1 = theData[i-1].GetX(); - x1 = theIntegral[i-1]; - y2 = theData[i].GetX(); - x2 = theIntegral[i]; - if (std::abs((y2-y1)/y2)<0.0000001) // not really necessary, since the case is excluded by construction + else if (GetVectorLength()==1) { - y1 = theData[i-2].GetX(); - x1 = theIntegral[i-2]; + result = theData[0].GetX(); + } + else + { + G4int i; + G4double rand = G4UniformRand(); + + // this was replaced + // for(i=1;i= 0 ; i--) + { + if(rand > theIntegral[i]/theIntegral[GetVectorLength()-1]) + break; + } + if(i!=GetVectorLength()-1) { + i++; + } + // until this (end) + + G4double x1, x2, y1, y2; + y1 = theData[i-1].GetX(); + x1 = theIntegral[i-1]; + y2 = theData[i].GetX(); + x2 = theIntegral[i]; + if(std::abs((y2-y1)/y2)<0.0000001) // not really necessary, since the case is excluded by construction + { + y1 = theData[i-2].GetX(); + x1 = theIntegral[i-2]; + } + result = theLin.Lin(rand, x1, x2, y1, y2); } - result = theLin.Lin(rand, x1, x2, y1, y2); - } - return result; + return result; + #endif } G4double * Debug() { - return theIntegral; + #if GEANT4_ENABLE_CUDA + return cudaVector->Debug(); + #else + return theIntegral; + #endif } inline void Integrate() { - G4int i; - if (nEntries == 1) - { - totalIntegral = 0; - return; - } - G4double sum = 0; - for (i = 1; i < GetVectorLength(); i++) - { - if (std::abs((theData[i].GetX() - theData[i-1].GetX()) / theData[i].GetX()) > 0.0000001) + #if GEANT4_ENABLE_CUDA + cudaVector->Integrate(); + #else + G4int i; + if(nEntries == 1) { - G4double x1 = theData[i-1].GetX(); - G4double x2 = theData[i].GetX(); - G4double y1 = theData[i-1].GetY(); - G4double y2 = theData[i].GetY(); - G4InterpolationScheme aScheme = theManager.GetScheme(i); - if (aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN) - { - sum+= 0.5*(y2+y1)*(x2-x1); - } - else if (aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG) - { - G4double a = y1; - G4double b = (y2-y1)/(G4Log(x2)-G4Log(x1)); - sum+= (a-b)*(x2-x1) + b*(x2*G4Log(x2)-x1*G4Log(x1)); - } - else if (aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN) - { - G4double a = G4Log(y1); - G4double b = (G4Log(y2)-G4Log(y1))/(x2-x1); - sum += (G4Exp(a)/b)*(G4Exp(b*x2)-G4Exp(b*x1)); - } - else if (aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO) - { - sum+= y1*(x2-x1); - } - else if (aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG) - { - G4double a = G4Log(y1); - G4double b = (G4Log(y2)-G4Log(y1))/(G4Log(x2)-G4Log(x1)); - sum += (G4Exp(a)/(b+1))*(G4Pow::GetInstance()->powA(x2,b+1)-G4Pow::GetInstance()->powA(x1,b+1)); - } - else + totalIntegral = 0; + return; + } + G4double sum = 0; + // G4InterpolationScheme scheme1 = theManager.GetScheme(1); + for(i=1;i0.0000001) { - throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate"); - } + G4double x1 = theData[i-1].GetX(); + G4double x2 = theData[i].GetX(); + G4double y1 = theData[i-1].GetY(); + G4double y2 = theData[i].GetY(); + G4InterpolationScheme aScheme = theManager.GetScheme(i); + if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN) + { + sum+= 0.5*(y2+y1)*(x2-x1); + } + else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG) + { + G4double a = y1; + G4double b = (y2-y1)/(G4Log(x2)-G4Log(x1)); + sum+= (a-b)*(x2-x1) + b*(x2*G4Log(x2)-x1*G4Log(x1)); + } + else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN) + { + G4double a = G4Log(y1); + G4double b = (G4Log(y2)-G4Log(y1))/(x2-x1); + sum += (G4Exp(a)/b)*(G4Exp(b*x2)-G4Exp(b*x1)); + } + else if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO) + { + sum+= y1*(x2-x1); + } + else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG) + { + G4double a = G4Log(y1); + G4double b = (G4Log(y2)-G4Log(y1))/(G4Log(x2)-G4Log(x1)); + sum += (G4Exp(a)/(b+1))*(G4Pow::GetInstance()->powA(x2,b+1)-G4Pow::GetInstance()->powA(x1,b+1)); + } + else + { + throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate"); + } + } } - } - totalIntegral = sum; + totalIntegral = sum; + #endif } inline void IntegrateAndNormalise() { - G4int i; - if (theIntegral != 0) { return; } - - theIntegral = new G4double[nEntries]; - if (nEntries == 1) - { - theIntegral[0] = 1; - return; - } - theIntegral[0] = 0; - G4double sum = 0; - G4double x1 = 0; - G4double x0 = 0; - for(i = 1; i < GetVectorLength(); i++) - { - x1 = theData[i].GetX(); - x0 = theData[i-1].GetX(); - if (std::abs(x1-x0) > std::abs(x1*0.0000001)) + #if GEANT4_ENABLE_CUDA + cudaVector->IntegrateAndNormalise(); + #else + G4int i; + if(theIntegral!=0) return; + theIntegral = new G4double[nEntries]; + if(nEntries == 1) + { + theIntegral[0] = 1; + return; + } + theIntegral[0] = 0; + G4double sum = 0; + G4double x1 = 0; + G4double x0 = 0; + for(i=1;i std::abs(x1*0.0000001) ) + { + G4InterpolationScheme aScheme = theManager.GetScheme(i); + G4double y0 = theData[i-1].GetY(); + G4double y1 = theData[i].GetY(); + G4double integ=theInt.GetBinIntegral(aScheme,x0,x1,y0,y1); + #if defined WIN32-VC + if(!_finite(integ)){integ=0;} + #elif defined __IBMCPP__ + if(isinf(integ)||isnan(integ)){integ=0;} + #else + if(std::isinf(integ)||std::isnan(integ)){integ=0;} + #endif + sum+=integ; + //******************************************************************** + } + theIntegral[i] = sum; + } + G4double total = theIntegral[GetVectorLength()-1]; + for(i=1;iTimes(factor); + #else + G4int i; + for(i=0; i theBlocked; - std::vector theBuffered; - G4double the15percentBorderCash; - G4double the50percentBorderCash; + //======================================= + #if GEANT4_ENABLE_CUDA + public: + G4ParticleHPVector_CUDA * cudaVector; + #else + G4ParticleHPInterpolator theLin; + G4double totalIntegral; + + G4ParticleHPDataPoint * theData; + G4InterpolationManager theManager; // knows how to interpolate the data. + G4double * theIntegral; + G4int nEntries; + G4int nPoints; + G4double label; + + G4ParticleHPInterpolator theInt; + G4int Verbose; + G4int isFreed; // debug only + + G4ParticleHPHash theHash; + G4double maxValue; + std::vector theBlocked; + // std::vector theBuffered; + G4double the15percentBorderCash; + G4double the50percentBorderCash; + #endif + //======================================= + }; -#endif \ No newline at end of file +#endif