[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 19-Aug-07 |
---|
[1055] | 28 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
[819] | 29 | // |
---|
| 30 | // |
---|
| 31 | // G4 Physics class: G4QIonIonCrossSection for gamma+A cross sections |
---|
| 32 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03 |
---|
| 33 | // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04 |
---|
| 34 | // -------------------------------------------------------------------------------- |
---|
| 35 | // **************************************************************************************** |
---|
| 36 | // ***** This HEADER is a property of the CHIPS hadronic package in Geant4 (M. Kosov) ***** |
---|
| 37 | // *********** DO NOT MAKE ANY CHANGE without approval of Mikhail.Kossov@cern.ch ********** |
---|
| 38 | // **************************************************************************************** |
---|
[1055] | 39 | // Short description: CHIPS cross-sectons for Ion-Ion interactions |
---|
| 40 | // --------------------------------------------------------------- |
---|
[819] | 41 | // |
---|
| 42 | //#define debug |
---|
| 43 | //#define pdebug |
---|
| 44 | //#define debug3 |
---|
| 45 | //#define debugn |
---|
| 46 | //#define debugs |
---|
| 47 | |
---|
| 48 | #include "G4QIonIonCrossSection.hh" |
---|
| 49 | |
---|
| 50 | // Initialization of the |
---|
| 51 | G4double* G4QIonIonCrossSection::lastLENI=0;// Pointer to the lastArray of LowEn Inelast CS |
---|
| 52 | G4double* G4QIonIonCrossSection::lastHENI=0;// Pointer to the lastArray of HighEn InelastCS |
---|
| 53 | G4double* G4QIonIonCrossSection::lastLENE=0;// Pointer to the lastArray of LowEn Elastic CS |
---|
| 54 | G4double* G4QIonIonCrossSection::lastHENE=0;// Pointer to the lastArray of HighEn ElasticCS |
---|
| 55 | G4int G4QIonIonCrossSection::lastPDG=0; // The last PDG code of the projectile |
---|
| 56 | G4int G4QIonIonCrossSection::lastN=0; // The last N of calculated nucleus |
---|
| 57 | G4int G4QIonIonCrossSection::lastZ=0; // The last Z of calculated nucleus |
---|
| 58 | G4double G4QIonIonCrossSection::lastP=0.; // Last used in cross section Momentum |
---|
| 59 | G4double G4QIonIonCrossSection::lastTH=0.; // Last threshold momentum |
---|
| 60 | G4double G4QIonIonCrossSection::lastICS=0.;// Last value of the Inelastic Cross Section |
---|
| 61 | G4double G4QIonIonCrossSection::lastECS=0.;// Last value of the Elastic Cross Section |
---|
| 62 | G4int G4QIonIonCrossSection::lastI=0; // The last position in the DAMDB |
---|
| 63 | |
---|
| 64 | // Returns Pointer to the G4VQCrossSection class |
---|
| 65 | G4VQCrossSection* G4QIonIonCrossSection::GetPointer() |
---|
| 66 | { |
---|
| 67 | static G4QIonIonCrossSection theCrossSection; //**Static body of Cross Section** |
---|
| 68 | return &theCrossSection; |
---|
| 69 | } |
---|
| 70 | |
---|
| 71 | // The main member function giving the collision cross section (P is in IU, CS is in mb) |
---|
| 72 | // Make pMom in independent units ! (Now it is MeV) |
---|
| 73 | G4double G4QIonIonCrossSection::GetCrossSection(G4bool fCS, G4double pMom, G4int tZ, |
---|
| 74 | G4int tN, G4int pPDG) |
---|
| 75 | { |
---|
| 76 | static G4int j; // A#0f records found in DB for this projectile |
---|
| 77 | static std::vector <G4int> colPDG;// Vector of the projectile PDG code |
---|
| 78 | static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) |
---|
| 79 | static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) |
---|
| 80 | static std::vector <G4double> colP; // Vector of last momenta for the reaction |
---|
| 81 | static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction |
---|
| 82 | static std::vector <G4double> colICS;// Vector of last inelastic cross-sections |
---|
| 83 | static std::vector <G4double> colECS;// Vector of last elastic cross-sections |
---|
| 84 | // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** |
---|
| 85 | #ifdef pdebug |
---|
| 86 | G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN |
---|
| 87 | <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz=" |
---|
| 88 | <<colN.size()<<G4endl; |
---|
[1055] | 89 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
[819] | 90 | #endif |
---|
| 91 | if(!pPDG) |
---|
| 92 | { |
---|
| 93 | #ifdef pdebug |
---|
| 94 | G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" ====> CS=0"<<G4endl; |
---|
| 95 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
| 96 | #endif |
---|
| 97 | return 0.; // projectile PDG=0 is a mistake (?!) @@ |
---|
| 98 | } |
---|
| 99 | G4bool in=false; // By default the isotope must be found in the AMDB |
---|
| 100 | if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope |
---|
| 101 | { |
---|
| 102 | in = false; // By default the isotope haven't be found in AMDB |
---|
| 103 | lastP = 0.; // New momentum history (nothing to compare with) |
---|
| 104 | lastPDG = pPDG; // The last PDG of the projectile |
---|
| 105 | lastN = tN; // The last N of the calculated nucleus |
---|
| 106 | lastZ = tZ; // The last Z of the calculated nucleus |
---|
| 107 | lastI = colN.size(); // Size of the Associative Memory DB in the heap |
---|
| 108 | j = 0; // A#0f records found in DB for this projectile |
---|
| 109 | if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found |
---|
[1055] | 110 | { // The nucleus with projPDG is found in AMDB |
---|
[819] | 111 | if(colN[i]==tN && colZ[i]==tZ) |
---|
[1055] | 112 | { |
---|
[819] | 113 | lastI=i; |
---|
| 114 | lastTH =colTH[i]; // Last THreshold (A-dependent) |
---|
| 115 | #ifdef pdebug |
---|
| 116 | G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; |
---|
| 117 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
| 118 | #endif |
---|
| 119 | if(pMom<=lastTH) |
---|
| 120 | { |
---|
| 121 | #ifdef pdebug |
---|
| 122 | G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl; |
---|
| 123 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
| 124 | #endif |
---|
| 125 | return 0.; // Energy is below the Threshold value |
---|
| 126 | } |
---|
| 127 | lastP =colP [i]; // Last Momentum (A-dependent) |
---|
| 128 | lastICS=colICS[i]; // Last Inelastic Cross-Section (A-dependent) |
---|
| 129 | lastECS=colECS[i]; // Last Elastic Cross-Section (A-dependent) |
---|
| 130 | if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
| 131 | { |
---|
| 132 | #ifdef pdebug |
---|
[1055] | 133 | G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS=" |
---|
| 134 | <<lastECS*millibarn<<G4endl; |
---|
[819] | 135 | #endif |
---|
| 136 | CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only |
---|
| 137 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
| 138 | return lastECS*millibarn; // Use theLastElasticCS |
---|
| 139 | } |
---|
| 140 | in = true; // This is the case when the isotop is found in DB |
---|
| 141 | // Momentum pMom is in IU ! @@ Units |
---|
| 142 | #ifdef pdebug |
---|
| 143 | G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; |
---|
| 144 | #endif |
---|
| 145 | lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update |
---|
| 146 | lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update |
---|
| 147 | #ifdef pdebug |
---|
| 148 | G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl; |
---|
| 149 | #endif |
---|
| 150 | if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold |
---|
| 151 | { |
---|
| 152 | #ifdef pdebug |
---|
| 153 | G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl; |
---|
| 154 | #endif |
---|
| 155 | lastTH=pMom; |
---|
| 156 | } |
---|
| 157 | break; // Go out of the LOOP |
---|
| 158 | } |
---|
| 159 | #ifdef pdebug |
---|
| 160 | G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] |
---|
| 161 | <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl; |
---|
| 162 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
| 163 | #endif |
---|
| 164 | j++; // Increment a#0f records found in DB for this pPDG |
---|
[1055] | 165 | } |
---|
| 166 | if(!in) // This nucleus has not been calculated previously |
---|
| 167 | { |
---|
[819] | 168 | #ifdef pdebug |
---|
| 169 | G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl; |
---|
| 170 | #endif |
---|
| 171 | //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU) |
---|
| 172 | lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create |
---|
| 173 | lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create |
---|
| 174 | if(lastICS<=0. || lastECS<=0.) |
---|
[1055] | 175 | { |
---|
[819] | 176 | lastTH = ThresholdEnergy(tZ, tN); // The Threshold Energy which is now the last |
---|
| 177 | #ifdef pdebug |
---|
| 178 | G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl; |
---|
| 179 | #endif |
---|
| 180 | if(pMom>lastTH) |
---|
| 181 | { |
---|
| 182 | #ifdef pdebug |
---|
| 183 | G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl; |
---|
| 184 | #endif |
---|
| 185 | lastTH=pMom; |
---|
| 186 | } |
---|
[1055] | 187 | } |
---|
[819] | 188 | #ifdef pdebug |
---|
| 189 | G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastICS<<",N="<<lastN<<",Z=" |
---|
| 190 | <<lastZ<<G4endl; |
---|
| 191 | #endif |
---|
| 192 | colN.push_back(tN); |
---|
| 193 | colZ.push_back(tZ); |
---|
| 194 | colPDG.push_back(pPDG); |
---|
| 195 | colP.push_back(pMom); |
---|
| 196 | colTH.push_back(lastTH); |
---|
| 197 | colICS.push_back(lastICS); |
---|
| 198 | colECS.push_back(lastECS); |
---|
| 199 | #ifdef pdebug |
---|
| 200 | G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn |
---|
| 201 | <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
| 202 | #endif |
---|
| 203 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
| 204 | return lastECS*millibarn; // Use theLastElasticCS |
---|
[1055] | 205 | } // End of creation of the new set of parameters |
---|
[819] | 206 | else |
---|
[1055] | 207 | { |
---|
[819] | 208 | #ifdef pdebug |
---|
| 209 | G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; |
---|
| 210 | #endif |
---|
| 211 | colP[lastI]=pMom; |
---|
| 212 | colPDG[lastI]=pPDG; |
---|
| 213 | colICS[lastI]=lastICS; |
---|
| 214 | colECS[lastI]=lastECS; |
---|
| 215 | } |
---|
| 216 | } // End of parameters udate |
---|
| 217 | else if(pMom<=lastTH) |
---|
| 218 | { |
---|
| 219 | #ifdef pdebug |
---|
| 220 | G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; |
---|
| 221 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
| 222 | #endif |
---|
| 223 | return 0.; // Momentum is below the Threshold Value -> CS=0 |
---|
| 224 | } |
---|
| 225 | else if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
| 226 | { |
---|
| 227 | #ifdef pdebug |
---|
| 228 | G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn |
---|
| 229 | <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
| 230 | #endif |
---|
| 231 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
| 232 | return lastECS*millibarn; // Use theLastElasticCS |
---|
| 233 | } |
---|
| 234 | else |
---|
| 235 | { |
---|
| 236 | #ifdef pdebug |
---|
| 237 | G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; |
---|
| 238 | #endif |
---|
| 239 | lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB |
---|
| 240 | lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB |
---|
| 241 | lastP=pMom; |
---|
| 242 | } |
---|
| 243 | #ifdef pdebug |
---|
| 244 | G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS=" |
---|
| 245 | <<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
| 246 | #endif |
---|
| 247 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
| 248 | return lastECS*millibarn; // Use theLastElasticCS |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | // The main member function giving the A-A cross section (Momentum in MeV, CS in mb) |
---|
| 252 | G4double G4QIonIonCrossSection::CalculateCrossSection(G4bool XS,G4int F,G4int I,G4int pPDG, |
---|
| 253 | G4int tZ,G4int tN, G4double TotMom) |
---|
| 254 | { |
---|
| 255 | //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)] |
---|
| 256 | //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21 |
---|
| 257 | // If change the following, please change in ::GetFunctions: |
---|
| 258 | static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold |
---|
| 259 | static const G4double dP=10.; // step for the LEN table |
---|
| 260 | static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106) |
---|
| 261 | static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part |
---|
| 262 | static const G4double Pmax=300000.; // maxE for the HighE part |
---|
| 263 | static const G4int nH=100; // A#of HResonance points in lnE |
---|
| 264 | static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part |
---|
| 265 | static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent) |
---|
| 266 | static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part |
---|
| 267 | // |
---|
| 268 | // Associative memory for acceleration |
---|
| 269 | static std::vector <G4double*> LENI; // Vector of pointers: LowEnIneIonIonCrossSection |
---|
| 270 | static std::vector <G4double*> HENI; // Vector of pointers: HighEnIneIonIonCrossSection |
---|
| 271 | static std::vector <G4double*> LENE; // Vector of pointers: LowEnElaIonIonCrossSection |
---|
| 272 | static std::vector <G4double*> HENE; // Vector of pointers: HighEnElaIonIonCrossSection |
---|
| 273 | #ifdef debug |
---|
| 274 | G4cout<<"G4QIonIonCrossSection::CalcCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<G4endl; |
---|
| 275 | #endif |
---|
| 276 | G4int dPDG=pPDG/10; // 10SZZZAAA |
---|
| 277 | G4int zPDG=dPDG/1000; // 10SZZZ (?) |
---|
| 278 | G4int zA=dPDG%1000; // proj A |
---|
| 279 | G4int pZ=zPDG%1000; // proj Z (?) |
---|
| 280 | G4int pN=zA-pZ; // proj N (?) |
---|
| 281 | G4double Momentum=TotMom/zA; // Momentum per nucleon |
---|
| 282 | if (Momentum<THmin) return 0.; // @@ This can be dangerouse for the heaviest nuc.! |
---|
| 283 | G4double sigma=0.; |
---|
| 284 | if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!! |
---|
| 285 | G4double tA=tN+tZ; // Target weight |
---|
| 286 | G4double pA=zA; // Projectile weight |
---|
| 287 | if(F<=0) // This isotope was not the last used isotop |
---|
| 288 | { |
---|
| 289 | if(F<0) // This isotope was found in DAMDB =========> RETRIEVE |
---|
[1055] | 290 | { |
---|
[819] | 291 | lastLENI=LENI[I]; // Pointer to Low Energy inelastic cross sections |
---|
| 292 | lastHENI=HENI[I]; // Pointer to High Energy inelastic cross sections |
---|
| 293 | lastLENE=LENE[I]; // Pointer to Low Energy inelastic cross sections |
---|
| 294 | lastHENE=HENE[I]; // Pointer to High Energy inelastic cross sections |
---|
| 295 | } |
---|
[1055] | 296 | else // This isotope wasn't calculated previously => CREATE |
---|
| 297 | { |
---|
[819] | 298 | lastLENI = new G4double[nL]; // Allocate memory for the new LEN cross sections |
---|
| 299 | lastHENI = new G4double[nH]; // Allocate memory for the new HEN cross sections |
---|
| 300 | lastLENE = new G4double[nL]; // Allocate memory for the new LEN cross sections |
---|
| 301 | lastHENE = new G4double[nH]; // Allocate memory for the new HEN cross sections |
---|
| 302 | G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE); |
---|
[1055] | 303 | if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl; |
---|
[819] | 304 | #ifdef debug |
---|
| 305 | G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA |
---|
| 306 | <<G4endl; |
---|
| 307 | #endif |
---|
| 308 | // *** The synchronization check *** |
---|
| 309 | G4int sync=LENI.size(); |
---|
| 310 | if(sync!=I) G4cerr<<"***G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl; |
---|
| 311 | LENI.push_back(lastLENI); // added LEN Inelastic |
---|
| 312 | HENI.push_back(lastHENI); // added HEN Inelastic |
---|
| 313 | LENE.push_back(lastLENE); // added LEN Elastic |
---|
| 314 | HENE.push_back(lastHENE); // added HEN Elastic |
---|
[1055] | 315 | } // End of creation of the new set of parameters |
---|
[819] | 316 | } // End of parameters udate |
---|
| 317 | // ============================== NOW the Magic Formula ================================= |
---|
| 318 | if (Momentum<lastTH) return 0.; // It must be already checked in the interface class |
---|
[1055] | 319 | else if (Momentum<Pmin) // LEN region (approximated in E, not in lnE) |
---|
[819] | 320 | { |
---|
| 321 | #ifdef debug |
---|
[1055] | 322 | G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl; |
---|
[819] | 323 | #endif |
---|
| 324 | if(tA<1. || pA<1.) |
---|
| 325 | { |
---|
[1055] | 326 | G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl; |
---|
[819] | 327 | sigma=0.; |
---|
| 328 | } |
---|
| 329 | else |
---|
| 330 | { |
---|
| 331 | if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLENI); |
---|
| 332 | else sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLENE); |
---|
| 333 | } |
---|
| 334 | #ifdef debugn |
---|
[1055] | 335 | if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P=" |
---|
[819] | 336 | <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl; |
---|
| 337 | #endif |
---|
| 338 | } |
---|
| 339 | else if (Momentum<Pmax) // High Energy region |
---|
| 340 | { |
---|
| 341 | G4double lP=std::log(Momentum); |
---|
| 342 | #ifdef debug |
---|
| 343 | G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl; |
---|
| 344 | #endif |
---|
| 345 | if(tA<=1. || pA<=1.) |
---|
| 346 | { |
---|
[1055] | 347 | G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<"or tA="<<tA<<" aren't composit"<<G4endl; |
---|
[819] | 348 | sigma=0.; |
---|
| 349 | } |
---|
| 350 | else |
---|
| 351 | { |
---|
| 352 | if(XS) sigma=EquLinearFit(lP,nH,milP,dlP,lastHENI); |
---|
| 353 | else sigma=EquLinearFit(lP,nH,milP,dlP,lastHENE); |
---|
| 354 | } |
---|
| 355 | } |
---|
| 356 | else // UltraHighE region (not frequent) |
---|
| 357 | { |
---|
[1055] | 358 | std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum); |
---|
[819] | 359 | if(XS) sigma=inelel.first; |
---|
| 360 | else sigma=inelel.second; |
---|
| 361 | } |
---|
| 362 | #ifdef debug |
---|
| 363 | G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl; |
---|
| 364 | #endif |
---|
| 365 | if(sigma<0.) return 0.; |
---|
| 366 | return sigma; |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | // Linear fit for YN[N] tabulated (from X0 with fixed step DX) function to X point |
---|
| 370 | |
---|
| 371 | // Calculate the functions for the log(A) |
---|
| 372 | G4int G4QIonIonCrossSection::GetFunctions(G4int pZ,G4int pN,G4int tZ,G4int tN,G4double* li, |
---|
| 373 | G4double* hi, G4double* le, G4double* he) |
---|
| 374 | { |
---|
| 375 | // If change the following, please change in ::CalculateCrossSection: |
---|
| 376 | static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold |
---|
| 377 | static const G4double dP=10.; // step for the LEN table |
---|
| 378 | static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106) |
---|
| 379 | static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part |
---|
| 380 | static const G4double Pmax=100000.; // maxE for the HighE part |
---|
| 381 | static const G4int nH=100; // A#of HResonance points in lnE |
---|
| 382 | static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part |
---|
| 383 | static const G4double malP=std::log(Pmax); // High logarithm energy |
---|
| 384 | static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part |
---|
| 385 | static const G4double lP=std::exp(dlP); // Multiplication factor in the HighE part |
---|
| 386 | // If the cross section approximation formula is changed - replace from file. |
---|
| 387 | if(pZ<1 || pN<0 || tZ<1 || pN<0) |
---|
| 388 | { |
---|
| 389 | G4cout<<"-W-G4QIonIonCS::GetFunct:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl; |
---|
| 390 | return -1; |
---|
| 391 | } |
---|
| 392 | G4double Mom=THmin; |
---|
| 393 | for(G4int k=0; k<nL; k++) |
---|
| 394 | { |
---|
| 395 | std::pair<G4double,G4double> len = CalculateXS(pZ, pN, tZ, tN, Mom); |
---|
| 396 | li[k]=len.first; |
---|
| 397 | le[k]=len.second; |
---|
| 398 | Mom+=dP; |
---|
| 399 | } |
---|
| 400 | G4double lMom=Pmin; |
---|
| 401 | for(G4int j=0; j<nH; j++) |
---|
| 402 | { |
---|
| 403 | std::pair<G4double,G4double> len = CalculateXS(pZ, pN, pZ, pN, Mom); |
---|
| 404 | hi[j]=len.first; |
---|
| 405 | he[j]=len.second; |
---|
| 406 | lMom*=lP; |
---|
| 407 | } |
---|
| 408 | #ifdef debug |
---|
| 409 | G4cout<<"G4QIonIonCS::GetFunctions:p="<<pA<<",t="<<tA<<" pair is calculated"<<G4endl; |
---|
| 410 | #endif |
---|
| 411 | return 1; |
---|
| 412 | } |
---|
| 413 | |
---|
| 414 | // Momentum (Mom=p/A) is in MeV/c, first=InelasticXS, second=ElasticXS (mb) |
---|
| 415 | std::pair<G4double,G4double> G4QIonIonCrossSection::CalculateXS(G4int pZ,G4int pN,G4int tZ, |
---|
| 416 | G4int tN, G4double Mom) |
---|
| 417 | { |
---|
| 418 | static G4VQCrossSection* ElCSman = G4QElasticCrossSection::GetPointer(); |
---|
| 419 | static G4VQCrossSection* InelPCSman = G4QProtonNuclearCrossSection::GetPointer(); |
---|
| 420 | static G4VQCrossSection* InelNCSman = G4QNeutronNuclearCrossSection::GetPointer(); |
---|
| 421 | G4double pA=pZ+pN; |
---|
| 422 | G4double tA=tZ+tN; |
---|
| 423 | if(pA<.9 || tA<.9 ||pA>239. || tA>239 || Mom < 0.) return std::make_pair(0.,0.); |
---|
| 424 | G4double inCS=0.; |
---|
| 425 | G4double elCS=0.; |
---|
| 426 | if(pA<1.1 || tA<1.1) // Ion-nucleon/nucleon-ion interaction use NA(in,el) |
---|
[1055] | 427 | { |
---|
[819] | 428 | if ( (pZ == 1 && !pN) || (tZ == 1 && !tN) ) // proton-nuclear |
---|
| 429 | { |
---|
| 430 | elCS=InelPCSman->GetCrossSection(true, Mom, tZ, tN, 2212); |
---|
| 431 | inCS=ElCSman->GetCrossSection(true, Mom, tZ, tN, 2212); |
---|
| 432 | } |
---|
| 433 | else if(pN==1 && !pZ) // neutron-nuclear |
---|
| 434 | { |
---|
| 435 | elCS=InelNCSman->GetCrossSection(true, Mom, tZ, tN, 2112); |
---|
| 436 | inCS=ElCSman->GetCrossSection(true, Mom, tZ, tN, 2112); |
---|
| 437 | } |
---|
[1055] | 438 | else G4cerr<<"-Warn-G4QIICS::CaCS:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl; |
---|
[819] | 439 | } |
---|
| 440 | else |
---|
[1055] | 441 | { |
---|
[819] | 442 | G4double P2=Mom*Mom; |
---|
| 443 | G4double P4=P2*P2; |
---|
| 444 | G4double P8=P4*P4; |
---|
| 445 | G4double T=ThresholdMomentum(pZ, pN, tZ, tN); // @@ Can be cashed as lastTH (?) |
---|
| 446 | G4double T2=T*T; |
---|
| 447 | G4double T4=T2*T2; |
---|
| 448 | G4double tot=CalculateTotal(pA, tA, Mom)*P8/(P8+T4*T4); // @@ convert to Indep. Units |
---|
| 449 | G4double rat=CalculateElTot(pA, tA, Mom); |
---|
| 450 | elCS=tot*rat; |
---|
| 451 | inCS=tot-elCS; |
---|
| 452 | } |
---|
| 453 | return std::make_pair(inCS,elCS); |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | // Total Ion-ion cross-section (mb), Momentum (Mom) here is p/A (MeV/c=IU) |
---|
| 457 | G4double G4QIonIonCrossSection::CalculateTotal(G4double pA, G4double tA, G4double Mom) |
---|
| 458 | { |
---|
| 459 | G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c |
---|
| 460 | G4double ab=pA+tA; |
---|
| 461 | G4double al=std::log(ab); |
---|
| 462 | G4double ap=std::log(pA*tA); |
---|
| 463 | G4double e=std::pow(pA,0.1)+std::pow(tA,0.1); |
---|
| 464 | G4double d=e-1.55/std::pow(al,0.2); |
---|
| 465 | G4double f=3.3+130./ab/ab+2.25/e; |
---|
| 466 | if(pA<4. || tA<4.) f=4.; |
---|
| 467 | G4double c=(385.+11.*ab/(1.+.02*ab*al)+(.5*ab+500./al/al)/al)*std::pow(d,f); |
---|
| 468 | G4double r=y-3.-4./ap/ap; |
---|
| 469 | return c+d*r*r; |
---|
| 470 | } |
---|
| 471 | |
---|
| 472 | // Ratio elastic/Total, Momentum (Mom) here is p/A (MeV/c=IU) |
---|
| 473 | G4double G4QIonIonCrossSection::CalculateElTot(G4double pA, G4double tA, G4double Mom) |
---|
| 474 | { |
---|
| 475 | G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c |
---|
| 476 | G4double ab=pA*tA; |
---|
| 477 | G4double ap=std::log(ab); |
---|
| 478 | G4double r=y-3.92-1.73/ap/ap; |
---|
| 479 | G4double d=.00166/(1.+.002*std::pow(ab,1.33333)); |
---|
| 480 | G4double al=std::log(pA+tA); |
---|
| 481 | G4double e=1.+.235*(std::fabs(ap-5.78)); |
---|
| 482 | if (std::fabs(pA-2.)<.1 && std::fabs(tA-2.)<.1) e=2.18; |
---|
| 483 | else if(std::fabs(pA-2.)<.1 && std::fabs(tA-3.)<.1) e=1.90; |
---|
| 484 | else if(std::fabs(pA-3.)<.1 && std::fabs(tA-2.)<.1) e=1.90; |
---|
| 485 | else if(std::fabs(pA-2.)<.1 && std::fabs(tA-4.)<.1) e=1.65; |
---|
| 486 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-2.)<.1) e=1.65; |
---|
| 487 | else if(std::fabs(pA-3.)<.1 && std::fabs(tA-4.)<.1) e=1.32; |
---|
| 488 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-3.)<.1) e=1.32; |
---|
| 489 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-4.)<.1) e=1.; |
---|
| 490 | G4double f=.37+.0188*al; |
---|
| 491 | G4double g=std::log(std::pow(pA,0.35)+std::pow(tA,0.35)); |
---|
| 492 | G4double h=g*g; |
---|
| 493 | G4double c=f/(1.+e/h/h); |
---|
| 494 | return c+d*r*r; |
---|
| 495 | } |
---|
| 496 | |
---|
| 497 | // Electromagnetic momentum/A-threshold (in MeV/c) |
---|
| 498 | G4double G4QIonIonCrossSection::ThresholdMomentum(G4int pZ, G4int pN, G4int tZ, G4int tN) |
---|
| 499 | { |
---|
| 500 | static const G4double third=1./3.; |
---|
| 501 | static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV |
---|
| 502 | static const G4double tpM= pM+pM; // Doubled proton mass (MeV) |
---|
| 503 | if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.; |
---|
| 504 | G4double tA=tZ+tN; |
---|
| 505 | G4double pA=pZ+pN; |
---|
| 506 | //G4double dE=1.263*tZ/(1.+std::pow(tA,third)); |
---|
| 507 | G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA |
---|
| 508 | return std::sqrt(dE*(tpM+dE)); |
---|
| 509 | } |
---|