G4GlauberGribovCrossSection.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 // author: V. Grichine
00027 // 
00028 // 17.07.06 V. Grichine - first implementation
00029 // 22.01.07 V.Ivanchenko - add interface with Z and A
00030 // 05.03.07 V.Ivanchenko - add IfZAApplicable
00031 // 11.06.10 V. Grichine - update for antiprotons
00032 // 10.11.11 V. Grichine - update for kaons
00033 
00034 #include "G4GlauberGribovCrossSection.hh"
00035 
00036 #include "G4PhysicalConstants.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4ParticleTable.hh"
00039 #include "G4IonTable.hh"
00040 #include "G4ParticleDefinition.hh"
00041 #include "G4HadronNucleonXsc.hh"
00042 
00043 // factory
00044 #include "G4CrossSectionFactory.hh"
00045 //
00046 G4_DECLARE_XS_FACTORY(G4GlauberGribovCrossSection);
00047 
00049 
00050 const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionTot[93] = {
00051 
00052 1.0, 1.0,     1.118517e+00, 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00, 
00053 1.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00, 
00054 1.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00,
00055 1.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00, 
00056 1.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00, 
00057 1.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00, 
00058 1.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00, 
00059 1.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00, 
00060 1.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00, 
00061 1.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00, 
00062 1.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00, 
00063 1.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00, 
00064 1.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00, 
00065 1.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00, 
00066 1.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00, 
00067 1.037075e+00, 1.034721e+00
00068 
00069 };
00070 
00071 const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionIn[93] = {
00072 
00073 1.0, 1.0,     1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00, 
00074 1.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00, 
00075 1.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00,
00076 1.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00, 
00077 1.131515e+00, 1.130338e+00, 1.134171e+00, 1.139206e+00, 1.141474e+00, 1.142189e+00, 
00078 1.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00, 
00079 1.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00, 
00080 1.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00, 
00081 1.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00, 
00082 1.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00, 
00083 1.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00, 
00084 1.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00, 
00085 1.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00, 
00086 1.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00, 
00087 1.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00, 
00088 1.046435e+00, 1.042614e+00
00089 
00090 };
00091 
00092 const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionTot[93] = {
00093 
00094 1.0, 1.0,     
00095 1.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00, 
00096 1.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00, 
00097 1.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00, 
00098 1.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00, 
00099 1.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00, 
00100 1.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00, 
00101 1.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00, 
00102 1.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00, 
00103 1.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00, 
00104 1.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00, 
00105 1.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00, 
00106 1.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00, 
00107 1.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00, 
00108 1.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00, 
00109 1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00, 
00110 1.034720e+00 
00111 
00112 };
00113 
00114 const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionIn[93] = {
00115 
00116 1.0, 1.0,     
00117 1.167419e+00, 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00, 
00118 1.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00, 
00119 1.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00, 
00120 1.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00, 
00121 1.130337e+00, 1.134170e+00, 1.139205e+00, 1.141472e+00, 1.142188e+00, 1.140724e+00, 
00122 1.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00, 
00123 1.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00, 
00124 1.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00, 
00125 1.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00, 
00126 1.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00, 
00127 1.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00, 
00128 1.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00, 
00129 1.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00, 
00130 1.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00, 
00131 1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00, 
00132 1.042613e+00 
00133 
00134 };
00135 
00136 
00137 const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionTot[93] = {
00138 
00139 1.0, 1.0,     
00140 1.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00, 
00141 1.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00, 
00142 1.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00, 
00143 1.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00, 
00144 1.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00, 
00145 1.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00, 
00146 1.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00, 
00147 1.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00, 
00148 1.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00, 
00149 1.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00, 
00150 1.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00, 
00151 1.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00, 
00152 1.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00, 
00153 1.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00, 
00154 1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00, 
00155 1.152974e+00
00156  
00157 };
00158 
00159 const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionIn[93] = {
00160 
00161 1.0, 1.0,    
00162 1.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.044495e+00, 1.062622e+00, 
00163 1.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.065100e+00, 
00164 1.070480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00, 
00165 1.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00, 
00166 1.137312e+00, 1.126263e+00, 1.126459e+00, 1.115191e+00, 1.116986e+00, 1.117184e+00, 
00167 1.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00, 
00168 1.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00, 
00169 1.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00, 
00170 1.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00, 
00171 1.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00, 
00172 1.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00, 
00173 1.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00, 
00174 1.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00, 
00175 1.071107e+00, 1.069955e+00, 1.064856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00, 
00176 1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00, 
00177 1.059658e+00
00178  
00179 };
00180 
00181 
00182 const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionTot[93] = {
00183 
00184 1.0, 1.0,     
00185 1.075927e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00, 
00186 1.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00, 
00187 1.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00, 
00188 1.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00, 
00189 1.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00, 
00190 1.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00, 
00191 1.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00, 
00192 1.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00, 
00193 1.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00, 
00194 1.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00, 
00195 1.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00, 
00196 1.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00, 
00197 1.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00, 
00198 1.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00, 
00199 1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00, 
00200 1.157267e+00 
00201 };
00202 
00203 
00204 const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionIn[93] = {
00205 
00206 1.0, 1.0,    
00207 1.140246e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.044514e+00, 1.062628e+00, 
00208 1.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00, 
00209 1.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00, 
00210 1.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00, 
00211 1.138502e+00, 1.127678e+00, 1.127244e+00, 1.116634e+00, 1.118347e+00, 1.118988e+00, 
00212 1.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00, 
00213 1.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00, 
00214 1.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00, 
00215 1.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00, 
00216 1.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00, 
00217 1.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00, 
00218 1.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00, 
00219 1.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00, 
00220 1.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00, 
00221 1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00, 
00222 1.062699e+00 
00223 
00224 };
00225 
00226 
00228 //
00229 
00230 G4GlauberGribovCrossSection::G4GlauberGribovCrossSection() 
00231  : G4VCrossSectionDataSet(Default_Name()),
00232    fUpperLimit(100000*GeV), fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
00233    fRadiusConst(1.08*fermi),  // 1.1, 1.3 ?
00234    fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
00235    fDiffractionXsc(0.0), fHadronNucleonXsc(0.0)
00236 {
00237   theGamma    = G4Gamma::Gamma();
00238   theProton   = G4Proton::Proton();
00239   theNeutron  = G4Neutron::Neutron();
00240   theAProton  = G4AntiProton::AntiProton();
00241   theANeutron = G4AntiNeutron::AntiNeutron();
00242   thePiPlus   = G4PionPlus::PionPlus();
00243   thePiMinus  = G4PionMinus::PionMinus();
00244   thePiZero   = G4PionZero::PionZero();
00245   theKPlus    = G4KaonPlus::KaonPlus();
00246   theKMinus   = G4KaonMinus::KaonMinus();
00247   theK0S      = G4KaonZeroShort::KaonZeroShort();
00248   theK0L      = G4KaonZeroLong::KaonZeroLong();
00249   theL        = G4Lambda::Lambda();
00250   theAntiL    = G4AntiLambda::AntiLambda();
00251   theSPlus    = G4SigmaPlus::SigmaPlus();
00252   theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
00253   theSMinus   = G4SigmaMinus::SigmaMinus();
00254   theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
00255   theS0       = G4SigmaZero::SigmaZero();
00256   theAS0      = G4AntiSigmaZero::AntiSigmaZero();
00257   theXiMinus  = G4XiMinus::XiMinus();
00258   theXi0      = G4XiZero::XiZero();
00259   theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00260   theAXi0     = G4AntiXiZero::AntiXiZero();
00261   theOmega    = G4OmegaMinus::OmegaMinus();
00262   theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
00263   theD        = G4Deuteron::Deuteron();
00264   theT        = G4Triton::Triton();
00265   theA        = G4Alpha::Alpha();
00266   theHe3      = G4He3::He3();
00267 
00268   hnXsc = new G4HadronNucleonXsc();
00269 }
00270 
00272 //
00273 //
00274 
00275 G4GlauberGribovCrossSection::~G4GlauberGribovCrossSection()
00276 {
00277   if (hnXsc) delete hnXsc;
00278 }
00279 
00281 
00282 G4bool 
00283 G4GlauberGribovCrossSection::IsIsoApplicable(const G4DynamicParticle* aDP, 
00284                                              G4int Z, G4int /*A*/, 
00285                                              const G4Element*,
00286                                              const G4Material*)
00287 {
00288   G4bool applicable      = false;
00289   // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
00290   G4double kineticEnergy = aDP->GetKineticEnergy();
00291 
00292   const G4ParticleDefinition* theParticle = aDP->GetDefinition();
00293  
00294   if ( ( kineticEnergy  >= fLowerLimit &&
00295          Z > 1 &&      // >=  He
00296        ( theParticle == theAProton   ||
00297          theParticle == theGamma     ||
00298          theParticle == theKPlus     ||
00299          theParticle == theKMinus    || 
00300          theParticle == theK0L     ||
00301          theParticle == theK0S    || 
00302          theParticle == theSMinus    ||  
00303          theParticle == theProton    ||
00304          theParticle == theNeutron   ||   
00305          theParticle == thePiPlus    ||
00306          theParticle == thePiMinus       ) )    ) applicable = true;
00307 
00308   return applicable;
00309 }
00310 
00312 //
00313 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
00314 // Glauber model with Gribov correction calculated in the dipole approximation on
00315 // light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
00316 // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
00317 
00318 G4double 
00319 G4GlauberGribovCrossSection::GetIsoCrossSection(const G4DynamicParticle* aParticle, 
00320                                                 G4int Z, G4int A,  
00321                                                 const G4Isotope*,
00322                                                 const G4Element*,
00323                                                 const G4Material*)
00324 {
00325   G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00326   G4double hpInXsc(0.), hnInXsc(0.);
00327   G4double R             = GetNucleusRadius(A); 
00328 
00329   G4int N = A - Z;              // number of neutrons
00330   if (N < 0) N = 0;
00331 
00332   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00333 
00334   if( theParticle == theProton  || 
00335       theParticle == theNeutron ||
00336       theParticle == thePiPlus  || 
00337       theParticle == thePiMinus      )
00338   {
00339     // sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00340 
00341     sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
00342 
00343     hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00344 
00345     sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
00346 
00347     hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00348 
00349     cofInelastic = 2.4;
00350     cofTotal     = 2.0;
00351   }
00352   else if( theParticle == theKPlus   || 
00353            theParticle == theKMinus  || 
00354            theParticle == theK0S     || 
00355            theParticle == theK0L        ) 
00356   {
00357     sigma        = GetKaonNucleonXscVector(aParticle, A, Z);
00358     cofInelastic = 2.2;
00359     cofTotal     = 2.0;
00360     R = 1.3*fermi;
00361     R *= std::pow(G4double(A), 0.3333);
00362   }
00363   else
00364   {
00365     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00366     cofInelastic = 2.2;
00367     cofTotal     = 2.0;
00368   }
00369   // cofInelastic = 2.0;
00370 
00371   if( A > 1 )
00372   { 
00373     nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00374     ratio = sigma/nucleusSquare;
00375 
00376     xsection =  nucleusSquare*std::log( 1. + ratio );
00377 
00378     xsection *= GetParticleBarCorTot(theParticle, Z);
00379 
00380     fTotalXsc = xsection;
00381 
00382   
00383 
00384     fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00385 
00386     fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
00387 
00388     fElasticXsc   = fTotalXsc - fInelasticXsc;
00389 
00390     if(fElasticXsc < 0.) fElasticXsc = 0.;
00391     
00392     G4double difratio = ratio/(1.+ratio);
00393 
00394     fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00395 
00396 
00397     // sigma = GetHNinelasticXsc(aParticle, A, Z);
00398 
00399     sigma = Z*hpInXsc + N*hnInXsc;
00400 
00401     ratio = sigma/nucleusSquare;
00402 
00403     fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00404 
00405     if (fElasticXsc < 0.) fElasticXsc = 0.;
00406   }
00407   else // H
00408   {
00409     fTotalXsc = sigma;
00410     xsection  = sigma;
00411     
00412     if ( theParticle != theAProton ) 
00413     {
00414       sigma         = GetHNinelasticXsc(aParticle, A, Z);
00415       fInelasticXsc = sigma;
00416       fElasticXsc   = fTotalXsc - fInelasticXsc;      
00417     }
00418     else
00419     {
00420       fElasticXsc   = fTotalXsc - fInelasticXsc;
00421     }
00422     if (fElasticXsc < 0.) fElasticXsc = 0.;
00423       
00424   }
00425   return xsection; 
00426 }
00427 
00429 //
00430 // Return single-diffraction/inelastic cross-section ratio
00431 
00432 G4double G4GlauberGribovCrossSection::
00433 GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
00434 {
00435   G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00436   G4double R             = GetNucleusRadius(A); 
00437 
00438   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00439 
00440   if( theParticle == theProton  || 
00441       theParticle == theNeutron ||
00442       theParticle == thePiPlus  || 
00443       theParticle == thePiMinus      )
00444   {
00445     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00446     cofInelastic = 2.4;
00447     cofTotal     = 2.0;
00448   }
00449   else
00450   {
00451     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00452     cofInelastic = 2.2;
00453     cofTotal     = 2.0;
00454   }
00455   nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00456   ratio = sigma/nucleusSquare;
00457 
00458   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00459    
00460   G4double difratio = ratio/(1.+ratio);
00461 
00462   fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00463 
00464   if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
00465   else                    ratio = 0.;
00466 
00467   return ratio; 
00468 }
00469 
00471 //
00472 // Return suasi-elastic/inelastic cross-section ratio
00473 
00474 G4double G4GlauberGribovCrossSection::
00475 GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
00476 {
00477   G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00478   G4double R             = GetNucleusRadius(A); 
00479 
00480   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00481 
00482   if( theParticle == theProton  || 
00483       theParticle == theNeutron ||
00484       theParticle == thePiPlus  || 
00485       theParticle == thePiMinus      )
00486   {
00487     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00488     cofInelastic = 2.4;
00489     cofTotal     = 2.0;
00490   }
00491   else
00492   {
00493     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00494     cofInelastic = 2.2;
00495     cofTotal     = 2.0;
00496   }
00497   nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00498   ratio = sigma/nucleusSquare;
00499 
00500   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00501 
00502   sigma = GetHNinelasticXsc(aParticle, A, Z);
00503   ratio = sigma/nucleusSquare;
00504 
00505   fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00506 
00507   if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
00508   else                                ratio = 0.;
00509   if ( ratio < 0. )                   ratio = 0.;
00510 
00511   return ratio; 
00512 }
00513 
00515 //
00516 // Returns hadron-nucleon Xsc according to differnt parametrisations:
00517 // [2] E. Levin, hep-ph/9710546
00518 // [3] U. Dersch, et al, hep-ex/9910052
00519 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 
00520 
00521 G4double 
00522 G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, 
00523                                                  const G4Element* anElement)
00524 {
00525   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00526   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00527 
00528   return GetHadronNucleonXsc(aParticle, At, Zt);
00529 }
00530 
00532 //
00533 // Returns hadron-nucleon Xsc according to differnt parametrisations:
00534 // [2] E. Levin, hep-ph/9710546
00535 // [3] U. Dersch, et al, hep-ex/9910052
00536 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 
00537 
00538 G4double 
00539 G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, 
00540                                                  G4int At, G4int /*Zt*/)
00541 {
00542   G4double xsection;
00543 
00544   //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
00545 
00546   G4double targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00547 
00548   G4double proj_mass     = aParticle->GetMass();
00549   G4double proj_momentum = aParticle->GetMomentum().mag();
00550   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00551 
00552   sMand /= GeV*GeV;  // in GeV for parametrisation
00553   proj_momentum /= GeV;
00554 
00555   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00556   
00557   G4double aa = At;
00558 
00559   if(theParticle == theGamma) 
00560   {
00561     xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
00562   } 
00563   else if(theParticle == theNeutron) // as proton ??? 
00564   {
00565     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00566   } 
00567   else if(theParticle == theProton) 
00568   {
00569     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00570     // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
00571     // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
00572   } 
00573   else if(theParticle == theAProton) 
00574   {
00575     xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
00576   } 
00577   else if(theParticle == thePiPlus) 
00578   {
00579     xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
00580   } 
00581   else if(theParticle == thePiMinus) 
00582   {
00583     // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
00584     xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
00585   } 
00586   else if(theParticle == theKPlus) 
00587   {
00588     xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
00589   } 
00590   else if(theParticle == theKMinus) 
00591   {
00592     xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
00593   }
00594   else  // as proton ??? 
00595   {
00596     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00597   } 
00598   xsection *= millibarn;
00599   return xsection;
00600 }
00601 
00602 
00604 //
00605 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
00606 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
00607 
00608 G4double 
00609 G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, 
00610                                                     const G4Element* anElement)
00611 {
00612   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00613   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00614 
00615   return GetHadronNucleonXscPDG(aParticle, At, Zt);
00616 }
00617 
00618 
00619 
00620 
00622 //
00623 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
00624 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
00625 //  At = number of nucleons,  Zt = number of protons 
00626 
00627 G4double 
00628 G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, 
00629                                                     G4int At, G4int Zt)
00630 {
00631   G4double xsection;
00632 
00633   G4int Nt = At-Zt;              // number of neutrons
00634   if (Nt < 0) Nt = 0;
00635   
00636   G4double zz = Zt;
00637   G4double aa = At;
00638   G4double nn = Nt;
00639 
00640   G4double targ_mass = G4ParticleTable::GetParticleTable()->
00641     GetIonTable()->GetIonMass(Zt, At);
00642 
00643   targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00644 
00645   G4double proj_mass     = aParticle->GetMass(); 
00646   G4double proj_momentum = aParticle->GetMomentum().mag();
00647 
00648   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00649 
00650   sMand         /= GeV*GeV;  // in GeV for parametrisation
00651 
00652   // General PDG fit constants
00653 
00654   G4double s0   = 5.38*5.38; // in Gev^2
00655   G4double eta1 = 0.458;
00656   G4double eta2 = 0.458;
00657   G4double B    = 0.308;
00658 
00659 
00660   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00661   
00662 
00663   if(theParticle == theNeutron) // proton-neutron fit 
00664   {
00665     xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00666                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00667     xsection  += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00668                       + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
00669   } 
00670   else if(theParticle == theProton) 
00671   {
00672       
00673       xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00674                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00675 
00676       xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00677                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00678   } 
00679   else if(theParticle == theAProton) 
00680   {
00681     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00682                           + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00683 
00684     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00685                           + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00686   } 
00687   else if(theParticle == thePiPlus) 
00688   {
00689     xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
00690                           + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
00691   } 
00692   else if(theParticle == thePiMinus) 
00693   {
00694     xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
00695                           + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
00696   } 
00697   else if(theParticle == theKPlus || theParticle == theK0L ) 
00698   {
00699     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
00700                           + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
00701 
00702     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
00703                           + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
00704   } 
00705   else if(theParticle == theKMinus || theParticle == theK0S ) 
00706   {
00707     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
00708                           + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
00709 
00710     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
00711                           + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
00712   }
00713   else if(theParticle == theSMinus) 
00714   {
00715     xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
00716                           - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
00717   } 
00718   else if(theParticle == theGamma) // modify later on
00719   {
00720     xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
00721                           + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
00722    
00723   } 
00724   else  // as proton ??? 
00725   {
00726     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00727                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00728 
00729     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00730                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00731   } 
00732   xsection *= millibarn; // parametrised in mb
00733   return xsection;
00734 }
00735 
00736 
00738 //
00739 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
00740 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
00741 
00742 G4double 
00743 G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle, 
00744                                                    const G4Element* anElement)
00745 {
00746   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00747   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00748 
00749   return GetHadronNucleonXscNS(aParticle, At, Zt);
00750 }
00751 
00752 
00753 
00754 
00756 //
00757 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
00758 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
00759 
00760 G4double 
00761 G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle, 
00762                                                    G4int At, G4int Zt)
00763 {
00764   G4double xsection(0);
00765   // G4double Delta;   DHW 19 May 2011: variable set but not used
00766   G4double A0, B0;
00767   G4double hpXscv(0);
00768   G4double hnXscv(0);
00769 
00770   G4int Nt = At-Zt;              // number of neutrons
00771   if (Nt < 0) Nt = 0;  
00772 
00773   G4double aa = At;
00774   G4double zz = Zt;
00775   G4double nn = Nt;
00776 
00777   G4double targ_mass = G4ParticleTable::GetParticleTable()->
00778   GetIonTable()->GetIonMass(Zt, At);
00779 
00780   targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00781 
00782   G4double proj_mass     = aParticle->GetMass();
00783   G4double proj_energy   = aParticle->GetTotalEnergy(); 
00784   G4double proj_momentum = aParticle->GetMomentum().mag();
00785 
00786   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00787 
00788   sMand         /= GeV*GeV;  // in GeV for parametrisation
00789   proj_momentum /= GeV;
00790   proj_energy   /= GeV;
00791   proj_mass     /= GeV;
00792 
00793   // General PDG fit constants
00794 
00795   G4double s0   = 5.38*5.38; // in Gev^2
00796   G4double eta1 = 0.458;
00797   G4double eta2 = 0.458;
00798   G4double B    = 0.308;
00799 
00800 
00801   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00802   
00803 
00804   if(theParticle == theNeutron) 
00805   {
00806     if( proj_momentum >= 373.)
00807     {
00808       return GetHadronNucleonXscPDG(aParticle,At,Zt);
00809     }
00810     else if( proj_momentum >= 10.)
00811     // if( proj_momentum >= 2.)
00812     {
00813       //  Delta = 1.;  // DHW 19 May 2011: variable set but not used
00814       // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
00815 
00816       if(proj_momentum >= 10.)
00817       {
00818         B0 = 7.5;
00819         A0 = 100. - B0*std::log(3.0e7);
00820 
00821         xsection = A0 + B0*std::log(proj_energy) - 11
00822                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00823                      0.93827*0.93827,-0.165);        //  mb
00824       }
00825       xsection *= zz + nn;
00826     }
00827     else
00828     {
00829       // nn to be pp
00830 
00831       if( proj_momentum < 0.73 )
00832       {
00833         hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00834       }
00835       else if( proj_momentum < 1.05  )
00836       {
00837        hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00838                          (std::log(proj_momentum/0.73));
00839       }
00840       else  // if( proj_momentum < 10.  )
00841       {
00842          hnXscv = 39.0+
00843               75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00844       }
00845       // pn to be np
00846 
00847       if( proj_momentum < 0.8 )
00848       {
00849         hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00850       }      
00851       else if( proj_momentum < 1.4 )
00852       {
00853         hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00854       }
00855       else    // if( proj_momentum < 10.  )
00856       {
00857         hpXscv = 33.3+
00858               20.8*(std::pow(proj_momentum,2.0)-1.35)/
00859                  (std::pow(proj_momentum,2.50)+0.95);
00860       }
00861       xsection = hpXscv*zz + hnXscv*nn;
00862     }
00863   } 
00864   else if(theParticle == theProton) 
00865   {
00866     if( proj_momentum >= 373.)
00867     {
00868       return GetHadronNucleonXscPDG(aParticle,At,Zt);
00869     }
00870     else if( proj_momentum >= 10.)
00871     // if( proj_momentum >= 2.)
00872     {
00873       // Delta = 1.;  DHW 19 May 2011: variable set but not used
00874       // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
00875 
00876       if(proj_momentum >= 10.)
00877       {
00878         B0 = 7.5;
00879         A0 = 100. - B0*std::log(3.0e7);
00880 
00881         xsection = A0 + B0*std::log(proj_energy) - 11
00882                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00883                      0.93827*0.93827,-0.165);        //  mb
00884       }
00885       xsection *= zz + nn;
00886     }
00887     else
00888     {
00889       // pp
00890 
00891       if( proj_momentum < 0.73 )
00892       {
00893         hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00894       }
00895       else if( proj_momentum < 1.05  )
00896       {
00897        hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00898                          (std::log(proj_momentum/0.73));
00899       }
00900       else    // if( proj_momentum < 10.  )
00901       {
00902          hpXscv = 39.0+
00903               75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00904       }
00905       // pn to be np
00906 
00907       if( proj_momentum < 0.8 )
00908       {
00909         hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00910       }      
00911       else if( proj_momentum < 1.4 )
00912       {
00913         hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00914       }
00915       else   // if( proj_momentum < 10.  )
00916       {
00917         hnXscv = 33.3+
00918               20.8*(std::pow(proj_momentum,2.0)-1.35)/
00919                  (std::pow(proj_momentum,2.50)+0.95);
00920       }
00921       xsection = hpXscv*zz + hnXscv*nn;
00922       // xsection = hpXscv*(Zt + Nt);
00923       // xsection = hnXscv*(Zt + Nt);
00924     }    
00925     // xsection *= 0.95;
00926   } 
00927   else if( theParticle == theAProton ) 
00928   {
00929     // xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00930     //                       + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00931 
00932     // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00933     //                    + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00934 
00935     G4double logP = std::log(proj_momentum);
00936 
00937     if( proj_momentum <= 1.0 )
00938     {
00939       xsection  = zz*(65.55 + 53.84/(proj_momentum+1.e-6)  );
00940     }
00941     else
00942     {
00943       xsection  = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68) 
00944                        + 0.293*logP*logP - 1.82*logP );
00945     }
00946     if ( nn > 0.)  
00947     {
00948       xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
00949     }
00950     else // H
00951     {
00952       fInelasticXsc =   38.0 + 38.0*std::pow( proj_momentum, -0.96) 
00953                         - 0.169*logP*logP;
00954       fInelasticXsc *=  millibarn;
00955     }    
00956   } 
00957   else if( theParticle == thePiPlus ) 
00958   {
00959     if(proj_momentum < 0.4)
00960     {
00961       G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
00962       hpXscv      = Ex3+20.0;
00963     }
00964     else if( proj_momentum < 1.15 )
00965     {
00966       G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
00967       hpXscv = Ex4+14.0;
00968     }
00969     else if(proj_momentum < 3.5)
00970     {
00971       G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
00972       G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
00973       hpXscv = Ex1+Ex2+27.5;
00974     }
00975     else //  if(proj_momentum > 3.5) // mb
00976     {
00977       hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
00978     }
00979     // pi+n = pi-p??
00980 
00981     if(proj_momentum < 0.37)
00982     {
00983       hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
00984     }
00985     else if(proj_momentum<0.65)
00986     {
00987        hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
00988     }
00989     else if(proj_momentum<1.3)
00990     {
00991       hnXscv = 36.1+
00992                 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
00993                 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
00994     }
00995     else if(proj_momentum<3.0)
00996     {
00997       hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
00998                 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
00999                 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01000     }
01001     else   // mb
01002     {
01003       hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
01004     }
01005     xsection = hpXscv*zz + hnXscv*nn;
01006   } 
01007   else if(theParticle == thePiMinus) 
01008   {
01009     // pi-n = pi+p??
01010 
01011     if(proj_momentum < 0.4)
01012     {
01013       G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
01014       hnXscv      = Ex3+20.0;
01015     }
01016     else if(proj_momentum < 1.15)
01017     {
01018       G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
01019       hnXscv = Ex4+14.0;
01020     }
01021     else if(proj_momentum < 3.5)
01022     {
01023       G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
01024       G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
01025       hnXscv = Ex1+Ex2+27.5;
01026     }
01027     else //  if(proj_momentum > 3.5) // mb
01028     {
01029       hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
01030     }
01031     // pi-p
01032 
01033     if(proj_momentum < 0.37)
01034     {
01035       hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
01036     }
01037     else if(proj_momentum<0.65)
01038     {
01039        hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
01040     }
01041     else if(proj_momentum<1.3)
01042     {
01043       hpXscv = 36.1+
01044                 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
01045                 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
01046     }
01047     else if(proj_momentum<3.0)
01048     {
01049       hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
01050                 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
01051                 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01052     }
01053     else   // mb
01054     {
01055       hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
01056     }
01057     xsection = hpXscv*zz + hnXscv*nn;
01058   } 
01059   else if(theParticle == theKPlus) 
01060   {
01061     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
01062                           + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
01063 
01064     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
01065                           + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
01066   } 
01067   else if(theParticle == theKMinus) 
01068   {
01069     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
01070                           + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
01071 
01072     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
01073                           + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
01074   }
01075   else if(theParticle == theSMinus) 
01076   {
01077     xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
01078                           - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
01079   } 
01080   else if(theParticle == theGamma) // modify later on
01081   {
01082     xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
01083                           + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
01084    
01085   } 
01086   else  // as proton ??? 
01087   {
01088     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
01089                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
01090 
01091     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
01092                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
01093   } 
01094   xsection *= millibarn; // parametrised in mb
01095   return xsection;
01096 }
01097 
01098 G4double 
01099 G4GlauberGribovCrossSection::GetKaonNucleonXscVector(const G4DynamicParticle* aParticle, 
01100                                                    G4int At, G4int Zt)
01101 {
01102   G4double Tkin, logTkin, xsc, xscP, xscN;
01103   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01104 
01105   G4int Nt = At-Zt;              // number of neutrons
01106   if (Nt < 0) Nt = 0;  
01107 
01108   Tkin = aParticle->GetKineticEnergy(); // Tkin in MeV
01109 
01110   if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);
01111 
01112   logTkin = std::log(Tkin); // Tkin in MeV!!!
01113 
01114  if( theParticle == theKPlus )
01115  {
01116    xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
01117    xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
01118  }
01119  else if( theParticle == theKMinus )
01120  {
01121    xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
01122    xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
01123  }
01124  else // K-zero as half of K+ and K-
01125  {
01126    xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
01127    xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
01128  }
01129  xsc = xscP*Zt + xscN*Nt;
01130   return xsc;
01131 }
01133 //
01134 // Returns hadron-nucleon inelastic cross-section based on proper parametrisation 
01135 
01136 G4double 
01137 G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle, 
01138                                                const G4Element* anElement)
01139 {
01140   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
01141   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
01142 
01143   return GetHNinelasticXsc(aParticle, At, Zt);
01144 }
01145 
01147 //
01148 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation 
01149 
01150 G4double 
01151 G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle, 
01152                                                      G4int At,  G4int Zt)
01153 {
01154   G4ParticleDefinition* hadron = aParticle->GetDefinition();
01155   G4double sumInelastic;
01156   G4int Nt = At - Zt;
01157   if(Nt < 0) Nt = 0;
01158   
01159   if( hadron == theKPlus )
01160   {
01161     sumInelastic =  GetHNinelasticXscVU(aParticle, At, Zt);
01162   }
01163   else
01164   {
01165     //sumInelastic  = Zt*GetHadronNucleonXscMK(aParticle, theProton);
01166     // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);    
01167     sumInelastic  = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
01168     sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);    
01169   } 
01170   return sumInelastic;
01171 }
01172 
01173 
01175 //
01176 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation 
01177 
01178 G4double 
01179 G4GlauberGribovCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle, 
01180                                                  G4int At, G4int Zt)
01181 {
01182   G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
01183   G4int absPDGcode = std::abs(PDGcode);
01184 
01185   G4double Elab = aParticle->GetTotalEnergy();              
01186                           // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
01187   G4double Plab = aParticle->GetMomentum().mag();            
01188                           // std::sqrt(Elab * Elab - 0.88);
01189 
01190   Elab /= GeV;
01191   Plab /= GeV;
01192 
01193   G4double LogPlab    = std::log( Plab );
01194   G4double sqrLogPlab = LogPlab * LogPlab;
01195 
01196   //G4cout<<"Plab = "<<Plab<<G4endl;
01197 
01198   G4double NumberOfTargetProtons = G4double(Zt); 
01199   G4double NumberOfTargetNucleons = G4double(At);
01200   G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
01201 
01202   if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
01203 
01204   G4double Xtotal, Xelastic, Xinelastic;
01205 
01206   if( absPDGcode > 1000 )  //------Projectile is baryon --------
01207   {
01208        G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
01209                          0.522*sqrLogPlab - 4.51*LogPlab;
01210 
01211        G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
01212                          0.513*sqrLogPlab - 4.27*LogPlab;
01213 
01214        G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01215                          0.169*sqrLogPlab - 1.85*LogPlab;
01216 
01217        G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01218                          0.169*sqrLogPlab - 1.85*LogPlab;
01219 
01220        Xtotal          = (NumberOfTargetProtons * XtotPP +
01221                           NumberOfTargetNeutrons * XtotPN);
01222 
01223        Xelastic        = (NumberOfTargetProtons * XelPP +
01224                           NumberOfTargetNeutrons * XelPN);
01225   }
01226   else if( PDGcode ==  211 ) //------Projectile is PionPlus -------
01227   {
01228        G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01229                           0.19 *sqrLogPlab - 0.0 *LogPlab;
01230 
01231        G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01232                           0.456*sqrLogPlab - 4.03*LogPlab;
01233 
01234        G4double XelPiP  =  0.0 + 11.4*std::pow(Plab,-0.40) +
01235                            0.079*sqrLogPlab - 0.0 *LogPlab;
01236 
01237        G4double XelPiN  = 1.76 + 11.2*std::pow(Plab,-0.64) +
01238                           0.043*sqrLogPlab - 0.0 *LogPlab;
01239 
01240        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01241                             NumberOfTargetNeutrons * XtotPiN  );
01242 
01243        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01244                             NumberOfTargetNeutrons * XelPiN   );
01245   }
01246   else if( PDGcode == -211 ) //------Projectile is PionMinus -------
01247   {
01248        G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01249                           0.456*sqrLogPlab - 4.03*LogPlab;
01250 
01251        G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01252                           0.19 *sqrLogPlab - 0.0 *LogPlab;
01253 
01254        G4double XelPiP  = 1.76 + 11.2*std::pow(Plab,-0.64) +
01255                           0.043*sqrLogPlab - 0.0 *LogPlab;
01256 
01257        G4double XelPiN  =  0.0 + 11.4*std::pow(Plab,-0.40) +
01258                            0.079*sqrLogPlab - 0.0 *LogPlab;
01259 
01260        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01261                             NumberOfTargetNeutrons * XtotPiN  );
01262 
01263        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01264                             NumberOfTargetNeutrons * XelPiN   );
01265   }
01266   else if( PDGcode ==  111 )  //------Projectile is PionZero  -------
01267   {
01268        G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
01269                           0.19 *sqrLogPlab - 0.0 *LogPlab +   //Pi+
01270                           33.0 + 14.0 *std::pow(Plab,-1.36) +
01271                           0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
01272 
01273        G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
01274                           0.456*sqrLogPlab - 4.03*LogPlab +   //Pi+
01275                           16.4 + 19.3 *std::pow(Plab,-0.42) +
01276                           0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01277 
01278        G4double XelPiP  =( 0.0 + 11.4*std::pow(Plab,-0.40) +
01279                            0.079*sqrLogPlab - 0.0 *LogPlab +    //Pi+
01280                            1.76 + 11.2*std::pow(Plab,-0.64) +
01281                            0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01282 
01283        G4double XelPiN  =( 1.76 + 11.2*std::pow(Plab,-0.64) +
01284                            0.043*sqrLogPlab - 0.0 *LogPlab +   //Pi+
01285                            0.0  + 11.4*std::pow(Plab,-0.40) +
01286                            0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01287 
01288        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01289                             NumberOfTargetNeutrons * XtotPiN  );
01290 
01291        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01292                             NumberOfTargetNeutrons * XelPiN   );
01293   }
01294   else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
01295   {
01296        G4double XtotKP = 18.1 +  0. *std::pow(Plab, 0.  ) +
01297                          0.26 *sqrLogPlab - 1.0 *LogPlab;
01298        G4double XtotKN = 18.7 +  0. *std::pow(Plab, 0.  ) +
01299                          0.21 *sqrLogPlab - 0.89*LogPlab;
01300 
01301        G4double XelKP  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
01302                           0.16 *sqrLogPlab - 1.3 *LogPlab;
01303 
01304        G4double XelKN  =  7.3 +  0. *std::pow(Plab,-0.  ) +
01305                           0.29 *sqrLogPlab - 2.4 *LogPlab;
01306 
01307        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01308                            NumberOfTargetNeutrons * XtotKN  );
01309 
01310        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01311                            NumberOfTargetNeutrons * XelKN   );
01312   }
01313   else if( PDGcode ==-321 )  //------Projectile is KaonMinus ------
01314   {
01315        G4double XtotKP = 32.1 +  0. *std::pow(Plab, 0.  ) +
01316                          0.66 *sqrLogPlab - 5.6 *LogPlab;
01317        G4double XtotKN = 25.2 +  0. *std::pow(Plab, 0.  ) +
01318                          0.38 *sqrLogPlab - 2.9 *LogPlab;
01319 
01320        G4double XelKP  =  7.3 +  0. *std::pow(Plab,-0.  ) +
01321                           0.29 *sqrLogPlab - 2.4 *LogPlab;
01322 
01323        G4double XelKN  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
01324                           0.16 *sqrLogPlab - 1.3 *LogPlab;
01325 
01326        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01327                            NumberOfTargetNeutrons * XtotKN  );
01328 
01329        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01330                            NumberOfTargetNeutrons * XelKN   );
01331   }
01332   else if( PDGcode == 311 ) //------Projectile is KaonZero ------
01333   {
01334        G4double XtotKP = ( 18.1 +  0. *std::pow(Plab, 0.  ) +
01335                           0.26 *sqrLogPlab - 1.0 *LogPlab +   //K+
01336                           32.1 +  0. *std::pow(Plab, 0.  ) +
01337                           0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
01338 
01339        G4double XtotKN = ( 18.7 +  0. *std::pow(Plab, 0.  ) +
01340                           0.21 *sqrLogPlab - 0.89*LogPlab +   //K+
01341                           25.2 +  0. *std::pow(Plab, 0.  ) +
01342                           0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
01343 
01344        G4double XelKP  = (  5.0 +  8.1*std::pow(Plab,-1.8 )
01345                            + 0.16 *sqrLogPlab - 1.3 *LogPlab +   //K+
01346                            7.3 +  0. *std::pow(Plab,-0.  ) +
01347                            0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
01348 
01349        G4double XelKN  = (  7.3 +  0. *std::pow(Plab,-0.  ) +
01350                            0.29 *sqrLogPlab - 2.4 *LogPlab +   //K+
01351                            5.0 +  8.1*std::pow(Plab,-1.8 ) +
01352                            0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
01353 
01354        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01355                            NumberOfTargetNeutrons * XtotKN  );
01356 
01357        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01358                            NumberOfTargetNeutrons * XelKN   );
01359   }
01360   else  //------Projectile is undefined, Nucleon assumed
01361   {
01362        G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
01363                          0.522*sqrLogPlab - 4.51*LogPlab;
01364 
01365        G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
01366                          0.513*sqrLogPlab - 4.27*LogPlab;
01367 
01368        G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01369                          0.169*sqrLogPlab - 1.85*LogPlab;
01370        G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01371                          0.169*sqrLogPlab - 1.85*LogPlab;
01372 
01373        Xtotal          = ( NumberOfTargetProtons  * XtotPP +
01374                            NumberOfTargetNeutrons * XtotPN  );
01375 
01376        Xelastic        = ( NumberOfTargetProtons  * XelPP  +
01377                            NumberOfTargetNeutrons * XelPN   );
01378   }
01379   Xinelastic = Xtotal - Xelastic;
01380 
01381   if( Xinelastic < 0.) Xinelastic = 0.;
01382 
01383   return Xinelastic*= millibarn;
01384 }
01385 
01387 //
01388 //
01389 
01390 G4double 
01391 G4GlauberGribovCrossSection::GetNucleusRadius(const G4DynamicParticle* , 
01392                                               const G4Element* anElement)
01393 {
01394   G4int At = G4lrint(anElement->GetN());
01395   G4double oneThird = 1.0/3.0;
01396   G4double cubicrAt = std::pow(G4double(At), oneThird); 
01397 
01398   G4double R;  // = fRadiusConst*cubicrAt;
01399   /*  
01400   G4double tmp = std::pow( cubicrAt-1., 3.);
01401   tmp         += At;
01402   tmp         *= 0.5;
01403 
01404   if (At > 20.)   // 20.
01405   {
01406     R = fRadiusConst*std::pow (tmp, oneThird); 
01407   }
01408   else
01409   {
01410     R = fRadiusConst*cubicrAt; 
01411   }
01412   */
01413   
01414   R = fRadiusConst*cubicrAt;
01415 
01416   G4double meanA  = 21.;
01417 
01418   G4double tauA1  = 40.; 
01419   G4double tauA2  = 10.; 
01420   G4double tauA3  = 5.; 
01421 
01422   G4double a1 = 0.85;
01423   G4double b1 = 1. - a1;
01424 
01425   G4double b2 = 0.3;
01426   G4double b3 = 4.;
01427 
01428   if (At > 20)   // 20.
01429   {
01430     R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) ); 
01431   }
01432   else if (At > 3)
01433   {
01434     R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) ); 
01435   }
01436   else 
01437   {
01438     R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) ); 
01439   }  
01440   return R;
01441  
01442 }
01444 //
01445 //
01446 
01447 G4double 
01448 G4GlauberGribovCrossSection::GetNucleusRadius(G4int At)
01449 {
01450   G4double oneThird = 1.0/3.0;
01451   G4double cubicrAt = std::pow(G4double(At), oneThird); 
01452 
01453   G4double R;  // = fRadiusConst*cubicrAt;
01454 
01455   /*
01456   G4double tmp = std::pow( cubicrAt-1., 3.);
01457   tmp         += At;
01458   tmp         *= 0.5;
01459 
01460   if (At > 20.)
01461   {
01462     R = fRadiusConst*std::pow (tmp, oneThird); 
01463   }
01464   else
01465   {
01466     R = fRadiusConst*cubicrAt; 
01467   }
01468   */
01469 
01470   R = fRadiusConst*cubicrAt;
01471 
01472   G4double meanA = 20.;
01473   G4double tauA  = 20.; 
01474 
01475   if (At > 20)   // 20.
01476   {
01477     R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) ); 
01478   }
01479   else
01480   {
01481     R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) ); 
01482   }
01483 
01484   return R;
01485 }
01486 
01488 //
01489 //
01490 
01491 G4double G4GlauberGribovCrossSection::CalculateEcmValue( const G4double mp , 
01492                                                          const G4double mt , 
01493                                                          const G4double Plab )
01494 {
01495   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01496   G4double Ecm  = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
01497   // G4double Pcm  = Plab * mt / Ecm;
01498   // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
01499 
01500   return Ecm ; // KEcm;
01501 }
01502 
01504 //
01505 //
01506 
01507 G4double G4GlauberGribovCrossSection::CalcMandelstamS( const G4double mp , 
01508                                                        const G4double mt , 
01509                                                        const G4double Plab )
01510 {
01511   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01512   G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
01513 
01514   return sMand;
01515 }
01516 
01518 //
01519 //
01520 
01521 void G4GlauberGribovCrossSection::CrossSectionDescription(std::ostream& outFile) const
01522 {
01523   outFile << "G4GlauberGribovCrossSection calculates total, inelastic and\n"
01524           << "elastic cross sections for hadron-nucleus cross sections using\n"
01525           << "the Glauber model with Gribov corrections.  It is valid for all\n"
01526           << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
01527           << "pi+, pi-, K+, K- and gammas with energies above 3 GeV.  This is\n"
01528           << "a cross section component which is to be used to build a cross\n"
01529           << "data set.\n";
01530 }
01531 
01532 //
01533 //

Generated on Mon May 27 17:48:23 2013 for Geant4 by  doxygen 1.4.7