[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 | // |
---|
[1055] | 27 | // $Id: G4QuasiFreeRatios.cc,v 1.21 2009/04/09 08:25:46 mkossov Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
[819] | 29 | // |
---|
| 30 | // |
---|
| 31 | // G4 Physics class: G4QuasiFreeRatios for N+A elastic cross sections |
---|
| 32 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01 |
---|
| 33 | // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06 |
---|
| 34 | // |
---|
[1055] | 35 | //======================================================================= |
---|
| 36 | // Short description: Provides percentage of quasi-free and quasi-elastic |
---|
| 37 | // reactions in the inelastic reactions. |
---|
| 38 | // ---------------------------------------------------------------------- |
---|
[819] | 39 | |
---|
| 40 | //#define debug |
---|
| 41 | //#define pdebug |
---|
| 42 | //#define ppdebug |
---|
| 43 | //#define nandebug |
---|
| 44 | |
---|
| 45 | #include "G4QuasiFreeRatios.hh" |
---|
| 46 | |
---|
| 47 | // initialisation of statics |
---|
| 48 | std::vector<G4double*> G4QuasiFreeRatios::vT; // Vector of pointers to LinTable in C++ heap |
---|
| 49 | std::vector<G4double*> G4QuasiFreeRatios::vL; // Vector of pointers to LogTable in C++ heap |
---|
| 50 | std::vector<std::pair<G4double,G4double>*> G4QuasiFreeRatios::vX; // ETPointers to LogTable |
---|
| 51 | |
---|
| 52 | G4QuasiFreeRatios::G4QuasiFreeRatios() |
---|
| 53 | { |
---|
| 54 | #ifdef pdebug |
---|
[1055] | 55 | G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl; |
---|
[819] | 56 | #endif |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | G4QuasiFreeRatios::~G4QuasiFreeRatios() |
---|
| 60 | { |
---|
| 61 | std::vector<G4double*>::iterator pos; |
---|
| 62 | for(pos=vT.begin(); pos<vT.end(); pos++) |
---|
| 63 | { delete [] *pos; } |
---|
| 64 | vT.clear(); |
---|
| 65 | for(pos=vL.begin(); pos<vL.end(); pos++) |
---|
| 66 | { delete [] *pos; } |
---|
| 67 | vL.clear(); |
---|
| 68 | |
---|
| 69 | std::vector<std::pair<G4double,G4double>*>::iterator pos2; |
---|
| 70 | for(pos2=vX.begin(); pos2<vX.end(); pos2++) |
---|
| 71 | { delete [] *pos2; } |
---|
| 72 | vX.clear(); |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | // Returns Pointer to the G4VQCrossSection class |
---|
| 76 | G4QuasiFreeRatios* G4QuasiFreeRatios::GetPointer() |
---|
| 77 | { |
---|
| 78 | static G4QuasiFreeRatios theRatios; // *** Static body of the QEl Cross Section *** |
---|
| 79 | return &theRatios; |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | // Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree) |
---|
| 83 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetRatios(G4double pIU, G4int pPDG, |
---|
| 84 | G4int tgZ, G4int tgN) |
---|
| 85 | { |
---|
[1055] | 86 | #ifdef pdebug |
---|
| 87 | G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl; |
---|
| 88 | #endif |
---|
[819] | 89 | G4double R=0.; |
---|
[1055] | 90 | G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot |
---|
| 91 | G4int tgA=tgZ+tgN; |
---|
| 92 | if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon |
---|
| 93 | std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU) |
---|
| 94 | //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE |
---|
| 95 | if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE |
---|
| 96 | else if(ElTot.second>0.) |
---|
[819] | 97 | { |
---|
[1055] | 98 | R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units |
---|
| 99 | QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio |
---|
[819] | 100 | } |
---|
[1055] | 101 | #ifdef pdebug |
---|
| 102 | G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl; |
---|
| 103 | #endif |
---|
[819] | 104 | return std::make_pair(QF2In,R); |
---|
| 105 | } |
---|
| 106 | |
---|
| 107 | // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A |
---|
| 108 | G4double G4QuasiFreeRatios::GetQF2IN_Ratio(G4double s, G4int A) |
---|
| 109 | { |
---|
| 110 | static const G4int nps=150; // Number of steps in the R(s) LinTable |
---|
| 111 | static const G4int mps=nps+1; // Number of elements in the R(s) LinTable |
---|
| 112 | static const G4double sma=150.; // The first LinTabEl(s=0)=1., s>sma -> logTab |
---|
| 113 | static const G4double ds=sma/nps; // Step of the linear Table |
---|
| 114 | static const G4int nls=100; // Number of steps in the R(lns) logTable |
---|
| 115 | static const G4int mls=nls+1; // Number of elements in the R(lns) logTable |
---|
| 116 | static const G4double lsi=5.; // The min ln(s) logTabEl(s=148.4 < sma=150.) |
---|
| 117 | static const G4double lsa=9.; // The max ln(s) logTabEl(s=148.4 - 8103. mb) |
---|
| 118 | static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 148.4 mb) |
---|
| 119 | static const G4double ms=std::exp(lsa);// The max s of logTabEl(~ 8103. mb) |
---|
| 120 | static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table |
---|
| 121 | static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table |
---|
| 122 | static const G4double toler=.01; // The tolarence mb defining the same cross-section |
---|
| 123 | static G4double lastS=0.; // The last sigma value for which R was calculated |
---|
| 124 | static G4double lastR=0.; // The last ratio R which was calculated |
---|
| 125 | // Local Associative Data Base: |
---|
| 126 | static std::vector<G4int> vA; // Vector of calculated A |
---|
| 127 | static std::vector<G4double> vH; // Vector of max s initialized in the LinTable |
---|
| 128 | static std::vector<G4int> vN; // Vector of topBin number initialized in LinTable |
---|
| 129 | static std::vector<G4double> vM; // Vector of rel max ln(s) initialized in LogTable |
---|
| 130 | static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable |
---|
| 131 | // Last values of the Associative Data Base: |
---|
| 132 | static G4int lastA=0; // theLast of calculated A |
---|
| 133 | static G4double lastH=0.; // theLast of max s initialized in the LinTable |
---|
| 134 | static G4int lastN=0; // theLast of topBin number initialized in LinTable |
---|
| 135 | static G4double lastM=0.; // theLast of rel max ln(s) initialized in LogTable |
---|
| 136 | static G4int lastK=0; // theLast of topBin number initialized in LogTable |
---|
| 137 | static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap |
---|
| 138 | static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap |
---|
| 139 | // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei |
---|
[1055] | 140 | #ifdef pdebug |
---|
| 141 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s<<G4endl; |
---|
| 142 | #endif |
---|
[819] | 143 | if(s<toler || A<2) return 1.; |
---|
| 144 | if(s>ms) return 0.; |
---|
| 145 | if(A>238) |
---|
| 146 | { |
---|
| 147 | G4cout<<"-Warning-G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl; |
---|
| 148 | return 0.; |
---|
| 149 | } |
---|
| 150 | G4int nDB=vA.size(); // A number of nuclei already initialized in AMDB |
---|
| 151 | // if(nDB && lastA==A && std::fabs(s-lastS)<toler) return lastR; |
---|
| 152 | if(nDB && lastA==A && s==lastS) return lastR; // VI do not use tolerance |
---|
| 153 | G4bool found=false; |
---|
| 154 | G4int i=-1; |
---|
[1055] | 155 | if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB |
---|
[819] | 156 | { |
---|
| 157 | found=true; // The A value is found |
---|
| 158 | break; |
---|
| 159 | } |
---|
[1055] | 160 | #ifdef pdebug |
---|
| 161 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl; |
---|
| 162 | #endif |
---|
[819] | 163 | if(!nDB || !found) // Create new line in the AMDB |
---|
[1055] | 164 | { |
---|
[819] | 165 | lastA = A; |
---|
| 166 | #ifdef pdebug |
---|
[1055] | 167 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
[819] | 168 | #endif |
---|
| 169 | lastT = new G4double[mps]; // Create the linear Table |
---|
| 170 | lastN = static_cast<int>(s/ds)+1; // MaxBin to be initialized |
---|
| 171 | if(lastN>nps) |
---|
| 172 | { |
---|
| 173 | lastN=nps; |
---|
| 174 | lastH=sma; |
---|
| 175 | } |
---|
| 176 | else lastH = lastN*ds; // Calculate max initialized s for LinTab |
---|
| 177 | G4double sv=0; |
---|
| 178 | lastT[0]=1.; |
---|
| 179 | for(G4int j=1; j<=lastN; j++) // Calculate LogTab values |
---|
| 180 | { |
---|
| 181 | sv+=ds; |
---|
| 182 | lastT[j]=CalcQF2IN_Ratio(sv,A); |
---|
| 183 | } |
---|
[1055] | 184 | lastL=new G4double[mls]; // Create the logarithmic Table |
---|
[819] | 185 | if(s>sma) // Initialize the logarithmic Table |
---|
| 186 | { |
---|
| 187 | #ifdef pdebug |
---|
[1055] | 188 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
[819] | 189 | #endif |
---|
| 190 | G4double ls=std::log(s); |
---|
| 191 | lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
| 192 | if(lastK>nls) |
---|
| 193 | { |
---|
| 194 | lastK=nls; |
---|
| 195 | lastM=lsa-lsi; |
---|
| 196 | } |
---|
| 197 | else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab |
---|
| 198 | sv=mi; |
---|
| 199 | for(G4int j=0; j<=lastK; j++) // Calculate LogTab values |
---|
| 200 | { |
---|
| 201 | lastL[j]=CalcQF2IN_Ratio(sv,A); |
---|
[1055] | 202 | if(j!=lastK) sv*=edl; |
---|
[819] | 203 | } |
---|
| 204 | } |
---|
| 205 | else // LogTab is not initialized |
---|
| 206 | { |
---|
| 207 | lastK = 0; |
---|
| 208 | lastM = 0.; |
---|
| 209 | } |
---|
| 210 | i++; // Make a new record to AMDB and position on it |
---|
| 211 | vA.push_back(lastA); |
---|
| 212 | vH.push_back(lastH); |
---|
| 213 | vN.push_back(lastN); |
---|
| 214 | vM.push_back(lastM); |
---|
| 215 | vK.push_back(lastK); |
---|
| 216 | vT.push_back(lastT); |
---|
| 217 | vL.push_back(lastL); |
---|
| 218 | } |
---|
| 219 | else // The A value was found in AMDB |
---|
[1055] | 220 | { |
---|
[819] | 221 | lastA=vA[i]; |
---|
| 222 | lastH=vH[i]; |
---|
| 223 | lastN=vN[i]; |
---|
| 224 | lastM=vM[i]; |
---|
| 225 | lastK=vK[i]; |
---|
| 226 | lastT=vT[i]; |
---|
| 227 | lastL=vL[i]; |
---|
[1055] | 228 | #ifdef pdebug |
---|
| 229 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s<<", lastM="<<lastM<<G4endl; |
---|
| 230 | #endif |
---|
[819] | 231 | if(s>lastM) // At least LinTab must be updated |
---|
| 232 | { |
---|
| 233 | G4int nextN=lastN+1; // The next bin to be initialized |
---|
[1055] | 234 | #ifdef pdebug |
---|
| 235 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl; |
---|
| 236 | #endif |
---|
[819] | 237 | if(lastN<nps) |
---|
| 238 | { |
---|
| 239 | lastN = static_cast<int>(s/ds)+1;// MaxBin to be initialized |
---|
| 240 | if(lastN>nps) |
---|
| 241 | { |
---|
| 242 | lastN=nps; |
---|
| 243 | lastH=sma; |
---|
| 244 | } |
---|
| 245 | else lastH = lastN*ds; // Calculate max initialized s for LinTab |
---|
| 246 | G4double sv=lastM; |
---|
| 247 | for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values |
---|
| 248 | { |
---|
| 249 | sv+=ds; |
---|
| 250 | lastT[j]=CalcQF2IN_Ratio(sv,A); |
---|
| 251 | } |
---|
[1055] | 252 | #ifdef pdebug |
---|
| 253 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
| 254 | #endif |
---|
[819] | 255 | } // End of LinTab update |
---|
[1055] | 256 | #ifdef pdebug |
---|
| 257 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl; |
---|
| 258 | #endif |
---|
[819] | 259 | if(lastN>=nextN) |
---|
[1055] | 260 | { |
---|
[819] | 261 | vH[i]=lastH; |
---|
| 262 | vN[i]=lastN; |
---|
| 263 | } |
---|
| 264 | G4int nextK=lastK+1; |
---|
[1055] | 265 | if(!lastK) nextK=0; |
---|
| 266 | #ifdef pdebug |
---|
| 267 | G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl; |
---|
| 268 | #endif |
---|
[819] | 269 | if(s>sma && lastK<nls) // LogTab must be updated |
---|
[1055] | 270 | { |
---|
[819] | 271 | G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed) |
---|
| 272 | G4double ls=std::log(s); |
---|
| 273 | lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
| 274 | if(lastK>nls) |
---|
| 275 | { |
---|
| 276 | lastK=nls; |
---|
| 277 | lastM=lsa-lsi; |
---|
| 278 | } |
---|
| 279 | else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab |
---|
[1055] | 280 | #ifdef pdebug |
---|
| 281 | G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl; |
---|
| 282 | #endif |
---|
[819] | 283 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
| 284 | { |
---|
[1055] | 285 | sv*=edl; |
---|
| 286 | #ifdef pdebug |
---|
| 287 | G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl; |
---|
| 288 | #endif |
---|
[819] | 289 | lastL[j]=CalcQF2IN_Ratio(sv,A); |
---|
| 290 | } |
---|
[1055] | 291 | #ifdef pdebug |
---|
| 292 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
| 293 | #endif |
---|
[819] | 294 | } // End of LogTab update |
---|
[1055] | 295 | #ifdef pdebug |
---|
| 296 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl; |
---|
| 297 | #endif |
---|
[819] | 298 | if(lastK>=nextK) |
---|
[1055] | 299 | { |
---|
[819] | 300 | vM[i]=lastM; |
---|
| 301 | vK[i]=lastK; |
---|
| 302 | } |
---|
| 303 | } |
---|
| 304 | } |
---|
[1055] | 305 | #ifdef pdebug |
---|
| 306 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s<<", sma="<<sma<<G4endl; |
---|
| 307 | #endif |
---|
[819] | 308 | // Now one can use tabeles to calculate the value |
---|
| 309 | if(s<sma) // Use linear table |
---|
[1055] | 310 | { |
---|
[819] | 311 | G4int n=static_cast<int>(s/ds); // Low edge number of the bin |
---|
| 312 | G4double d=s-n*ds; // Linear shift |
---|
| 313 | G4double v=lastT[n]; // Base |
---|
| 314 | lastR=v+d*(lastT[n+1]-v)/ds; // Result |
---|
| 315 | } |
---|
[1055] | 316 | else // Use log table |
---|
| 317 | { |
---|
[819] | 318 | G4double ls=std::log(s)-lsi; // ln(s)-l_min |
---|
| 319 | G4int n=static_cast<int>(ls/dl); // Low edge number of the bin |
---|
| 320 | G4double d=ls-n*dl; // Log shift |
---|
| 321 | G4double v=lastL[n]; // Base |
---|
| 322 | lastR=v+d*(lastL[n+1]-v)/dl; // Result |
---|
| 323 | } |
---|
| 324 | if(lastR<0.) lastR=0.; |
---|
| 325 | if(lastR>1.) lastR=1.; |
---|
[1055] | 326 | #ifdef pdebug |
---|
| 327 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl; |
---|
| 328 | #endif |
---|
[819] | 329 | return lastR; |
---|
| 330 | } // End of CalcQF2IN_Ratio |
---|
| 331 | |
---|
| 332 | // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A |
---|
| 333 | G4double G4QuasiFreeRatios::CalcQF2IN_Ratio(G4double s, G4int A) |
---|
| 334 | { |
---|
| 335 | static const G4double C=1.246; |
---|
[1055] | 336 | G4double s2=s*s; |
---|
[819] | 337 | G4double s4=s2*s2; |
---|
[1055] | 338 | G4double ss=std::sqrt(std::sqrt(s)); |
---|
[819] | 339 | G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2); |
---|
| 340 | G4double E=.2644+.016/(1.+std::exp((29.54-s)/2.49)); |
---|
| 341 | G4double F=ss*.1526*std::exp(-s2*ss*.0000859); |
---|
[1055] | 342 | return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P); |
---|
[819] | 343 | } // End of CalcQF2IN_Ratio |
---|
| 344 | |
---|
| 345 | // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot) |
---|
| 346 | std::pair<G4double,G4double> G4QuasiFreeRatios::CalcElTot(G4double p, G4int I) |
---|
| 347 | { |
---|
| 348 | // ---------> Each parameter set can have not more than nPoints=128 parameters |
---|
| 349 | static const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola |
---|
| 350 | static const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient |
---|
| 351 | static const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient |
---|
| 352 | static const G4double pmi=.1; // Below that fast LE calculation is made |
---|
| 353 | static const G4double pma=1000.; // Above that fast HE calculation is made |
---|
| 354 | // ====================================================================================== |
---|
| 355 | G4double El=0.; // prototype of the elastic hN cross-section |
---|
| 356 | G4double To=0.; // prototype of the total hN cross-section |
---|
| 357 | if(p<=0.) |
---|
| 358 | { |
---|
| 359 | G4cout<<"-Warning-G4QuasiFreeRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl; |
---|
| 360 | return std::make_pair(El,To); |
---|
| 361 | } |
---|
| 362 | if (!I) // pp/nn |
---|
[1055] | 363 | { |
---|
[819] | 364 | #ifdef debug |
---|
[1055] | 365 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl; |
---|
[819] | 366 | #endif |
---|
| 367 | if(p<pmi) |
---|
| 368 | { |
---|
| 369 | G4double p2=p*p; |
---|
| 370 | El=1./(.00012+p2*.2); |
---|
| 371 | To=El; |
---|
| 372 | #ifdef debug |
---|
[1055] | 373 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl; |
---|
[819] | 374 | #endif |
---|
| 375 | } |
---|
| 376 | else if(p>pma) |
---|
| 377 | { |
---|
| 378 | G4double lp=std::log(p)-lmi; |
---|
| 379 | G4double lp2=lp*lp; |
---|
| 380 | El=pbe*lp2+6.72; |
---|
| 381 | To=pbt*lp2+38.2; |
---|
| 382 | #ifdef debug |
---|
[1055] | 383 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl; |
---|
[819] | 384 | #endif |
---|
| 385 | } |
---|
| 386 | else |
---|
| 387 | { |
---|
| 388 | G4double p2=p*p; |
---|
| 389 | G4double LE=1./(.00012+p2*.2); |
---|
| 390 | G4double lp=std::log(p)-lmi; |
---|
| 391 | G4double lp2=lp*lp; |
---|
| 392 | G4double rp2=1./p2; |
---|
| 393 | El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p); |
---|
| 394 | To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2); |
---|
| 395 | #ifdef debug |
---|
[1055] | 396 | G4cout<<"G4QuasiFreeR::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl; |
---|
[819] | 397 | #endif |
---|
| 398 | } |
---|
| 399 | } |
---|
| 400 | else if(I==1) // np/pn |
---|
[1055] | 401 | { |
---|
[819] | 402 | if(p<pmi) |
---|
| 403 | { |
---|
| 404 | G4double p2=p*p; |
---|
| 405 | El=1./(.00012+p2*(.051+.1*p2)); |
---|
| 406 | To=El; |
---|
| 407 | } |
---|
| 408 | else if(p>pma) |
---|
| 409 | { |
---|
| 410 | G4double lp=std::log(p)-lmi; |
---|
| 411 | G4double lp2=lp*lp; |
---|
| 412 | El=pbe*lp2+6.72; |
---|
| 413 | To=pbt*lp2+38.2; |
---|
| 414 | } |
---|
| 415 | else |
---|
| 416 | { |
---|
| 417 | G4double p2=p*p; |
---|
| 418 | G4double LE=1./(.00012+p2*(.051+.1*p2)); |
---|
| 419 | G4double lp=std::log(p)-lmi; |
---|
| 420 | G4double lp2=lp*lp; |
---|
| 421 | G4double rp2=1./p2; |
---|
| 422 | El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p); |
---|
| 423 | To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2); |
---|
| 424 | } |
---|
| 425 | } |
---|
| 426 | else if(I==2) // pimp/pipn |
---|
[1055] | 427 | { |
---|
[819] | 428 | G4double lp=std::log(p); |
---|
| 429 | if(p<pmi) |
---|
| 430 | { |
---|
| 431 | G4double lr=lp+1.27; |
---|
| 432 | El=1.53/(lr*lr+.0676); |
---|
| 433 | To=El*3; |
---|
| 434 | } |
---|
| 435 | else if(p>pma) |
---|
| 436 | { |
---|
| 437 | G4double ld=lp-lmi; |
---|
| 438 | G4double ld2=ld*ld; |
---|
| 439 | G4double sp=std::sqrt(p); |
---|
| 440 | El=pbe*ld2+2.4+7./sp; |
---|
| 441 | To=pbt*ld2+22.3+12./sp; |
---|
| 442 | } |
---|
| 443 | else |
---|
| 444 | { |
---|
| 445 | G4double lr=lp+1.27; |
---|
| 446 | G4double LE=1.53/(lr*lr+.0676); |
---|
| 447 | G4double ld=lp-lmi; |
---|
| 448 | G4double ld2=ld*ld; |
---|
| 449 | G4double p2=p*p; |
---|
| 450 | G4double p4=p2*p2; |
---|
| 451 | G4double sp=std::sqrt(p); |
---|
| 452 | G4double lm=lp+.36; |
---|
| 453 | G4double md=lm*lm+.04; |
---|
| 454 | G4double lh=lp-.017; |
---|
| 455 | G4double hd=lh*lh+.0025; |
---|
| 456 | El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd; |
---|
| 457 | To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd; |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | else if(I==3) // pipp/pimn |
---|
[1055] | 461 | { |
---|
[819] | 462 | G4double lp=std::log(p); |
---|
| 463 | if(p<pmi) |
---|
| 464 | { |
---|
| 465 | G4double lr=lp+1.27; |
---|
| 466 | G4double lr2=lr*lr; |
---|
| 467 | El=13./(lr2+lr2*lr2+.0676); |
---|
| 468 | To=El; |
---|
| 469 | } |
---|
| 470 | else if(p>pma) |
---|
| 471 | { |
---|
| 472 | G4double ld=lp-lmi; |
---|
| 473 | G4double ld2=ld*ld; |
---|
| 474 | G4double sp=std::sqrt(p); |
---|
| 475 | El=pbe*ld2+2.4+6./sp; |
---|
| 476 | To=pbt*ld2+22.3+5./sp; |
---|
| 477 | } |
---|
| 478 | else |
---|
| 479 | { |
---|
| 480 | G4double lr=lp+1.27; |
---|
| 481 | G4double lr2=lr*lr; |
---|
| 482 | G4double LE=13./(lr2+lr2*lr2+.0676); |
---|
| 483 | G4double ld=lp-lmi; |
---|
| 484 | G4double ld2=ld*ld; |
---|
| 485 | G4double p2=p*p; |
---|
| 486 | G4double p4=p2*p2; |
---|
| 487 | G4double sp=std::sqrt(p); |
---|
| 488 | G4double lm=lp-.32; |
---|
| 489 | G4double md=lm*lm+.0576; |
---|
| 490 | El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; |
---|
| 491 | To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md; |
---|
| 492 | } |
---|
| 493 | } |
---|
[1055] | 494 | else if(I==4) // Kmp/Kmn/K0p/K0n |
---|
| 495 | { |
---|
[819] | 496 | |
---|
| 497 | if(p<pmi) |
---|
| 498 | { |
---|
| 499 | G4double psp=p*std::sqrt(p); |
---|
| 500 | El=5.2/psp; |
---|
| 501 | To=14./psp; |
---|
| 502 | } |
---|
| 503 | else if(p>pma) |
---|
| 504 | { |
---|
| 505 | G4double ld=std::log(p)-lmi; |
---|
| 506 | G4double ld2=ld*ld; |
---|
| 507 | El=pbe*ld2+2.23; |
---|
| 508 | To=pbt*ld2+19.5; |
---|
| 509 | } |
---|
| 510 | else |
---|
| 511 | { |
---|
| 512 | G4double ld=std::log(p)-lmi; |
---|
| 513 | G4double ld2=ld*ld; |
---|
| 514 | G4double sp=std::sqrt(p); |
---|
| 515 | G4double psp=p*sp; |
---|
| 516 | G4double p2=p*p; |
---|
| 517 | G4double p4=p2*p2; |
---|
| 518 | G4double lm=p-.39; |
---|
| 519 | G4double md=lm*lm+.000156; |
---|
| 520 | G4double lh=p-1.; |
---|
| 521 | G4double hd=lh*lh+.0156; |
---|
| 522 | El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd; |
---|
| 523 | To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd; |
---|
| 524 | } |
---|
| 525 | } |
---|
| 526 | else if(I==5) // Kpp/Kpn/aKp/aKn |
---|
[1055] | 527 | { |
---|
[819] | 528 | if(p<pmi) |
---|
| 529 | { |
---|
| 530 | G4double lr=p-.38; |
---|
| 531 | G4double lm=p-1.; |
---|
| 532 | G4double md=lm*lm+.372; |
---|
| 533 | El=.7/(lr*lr+.0676)+2./md; |
---|
| 534 | To=El+.6/md; |
---|
| 535 | } |
---|
| 536 | else if(p>pma) |
---|
| 537 | { |
---|
| 538 | G4double ld=std::log(p)-lmi; |
---|
| 539 | G4double ld2=ld*ld; |
---|
| 540 | El=pbe*ld2+2.23; |
---|
| 541 | To=pbt*ld2+19.5; |
---|
| 542 | } |
---|
| 543 | else |
---|
| 544 | { |
---|
| 545 | G4double ld=std::log(p)-lmi; |
---|
| 546 | G4double ld2=ld*ld; |
---|
| 547 | G4double lr=p-.38; |
---|
| 548 | G4double LE=.7/(lr*lr+.0676); |
---|
| 549 | G4double sp=std::sqrt(p); |
---|
| 550 | G4double p2=p*p; |
---|
| 551 | G4double p4=p2*p2; |
---|
| 552 | G4double lm=p-1.; |
---|
| 553 | G4double md=lm*lm+.372; |
---|
| 554 | El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md; |
---|
| 555 | To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md; |
---|
| 556 | } |
---|
| 557 | } |
---|
| 558 | else if(I==6) // hyperon-N |
---|
[1055] | 559 | { |
---|
[819] | 560 | if(p<pmi) |
---|
| 561 | { |
---|
| 562 | G4double p2=p*p; |
---|
| 563 | El=1./(.002+p2*(.12+p2)); |
---|
| 564 | To=El; |
---|
| 565 | } |
---|
| 566 | else if(p>pma) |
---|
| 567 | { |
---|
| 568 | G4double lp=std::log(p)-lmi; |
---|
| 569 | G4double lp2=lp*lp; |
---|
| 570 | G4double sp=std::sqrt(p); |
---|
| 571 | El=(pbe*lp2+6.72)/(1.+2./sp); |
---|
| 572 | To=(pbt*lp2+38.2+900./sp)/(1.+27./sp); |
---|
| 573 | } |
---|
| 574 | else |
---|
| 575 | { |
---|
| 576 | G4double p2=p*p; |
---|
| 577 | G4double LE=1./(.002+p2*(.12+p2)); |
---|
| 578 | G4double lp=std::log(p)-lmi; |
---|
| 579 | G4double lp2=lp*lp; |
---|
| 580 | G4double p4=p2*p2; |
---|
| 581 | G4double sp=std::sqrt(p); |
---|
| 582 | El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4); |
---|
| 583 | To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4); |
---|
| 584 | } |
---|
| 585 | } |
---|
| 586 | else if(I==7) // antibaryon-N |
---|
[1055] | 587 | { |
---|
[819] | 588 | if(p>pma) |
---|
| 589 | { |
---|
| 590 | G4double lp=std::log(p)-lmi; |
---|
| 591 | G4double lp2=lp*lp; |
---|
| 592 | El=pbe*lp2+6.72; |
---|
| 593 | To=pbt*lp2+38.2; |
---|
| 594 | } |
---|
| 595 | else |
---|
| 596 | { |
---|
| 597 | G4double ye=std::pow(p,1.25); |
---|
| 598 | G4double yt=std::pow(p,.35); |
---|
| 599 | G4double lp=std::log(p)-lmi; |
---|
| 600 | G4double lp2=lp*lp; |
---|
| 601 | El=80./(ye+1.)+pbe*lp2+6.72; |
---|
| 602 | To=(80./yt+.3)/yt+pbt*lp2+38.2; |
---|
| 603 | } |
---|
| 604 | } |
---|
| 605 | else |
---|
| 606 | { |
---|
| 607 | G4cout<<"*Error*G4QuasiFreeRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl; |
---|
| 608 | G4Exception("G4QuasiFreeRatios::CalcElTot:","23",FatalException,"CHIPScrash"); |
---|
| 609 | } |
---|
| 610 | if(El>To) El=To; |
---|
| 611 | return std::make_pair(El,To); |
---|
| 612 | } // End of CalcElTot |
---|
| 613 | |
---|
| 614 | // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron |
---|
| 615 | std::pair<G4double,G4double> G4QuasiFreeRatios::FetchElTot(G4double p, G4int PDG, G4bool F) |
---|
| 616 | { |
---|
| 617 | static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step) |
---|
| 618 | static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable |
---|
| 619 | static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
| 620 | static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
| 621 | static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c) |
---|
| 622 | static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV) |
---|
| 623 | static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table |
---|
| 624 | static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table |
---|
[1055] | 625 | //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum" |
---|
[819] | 626 | static G4double lastP=0.; // The last momentum for which XS was calculated |
---|
| 627 | static G4int lastH=0; // The last projPDG for which XS was calculated |
---|
| 628 | static G4bool lastF=true; // The last nucleon for which XS was calculated |
---|
| 629 | static std::pair<G4double,G4double> lastR=std::make_pair(0.,0.); // The last result |
---|
| 630 | // Local Associative Data Base: |
---|
| 631 | static std::vector<G4int> vI; // Vector of index for which XS was calculated |
---|
| 632 | static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable |
---|
| 633 | static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable |
---|
| 634 | // Last values of the Associative Data Base: |
---|
| 635 | static G4int lastI=0; // The Last index for which XS was calculated |
---|
| 636 | static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable |
---|
| 637 | static G4int lastK=0; // The Last topBin number initialized in LogTable |
---|
| 638 | static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap |
---|
| 639 | // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei |
---|
| 640 | G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB |
---|
| 641 | #ifdef pdebug |
---|
| 642 | G4cout<<"G4QuasiFreeR::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl; |
---|
| 643 | #endif |
---|
[1055] | 644 | if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler. |
---|
[819] | 645 | // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
| 646 | lastH=PDG; |
---|
| 647 | lastF=F; |
---|
| 648 | G4int ind=-1; // Prototipe of the index of the PDG/F combination |
---|
| 649 | // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p), |
---|
[1055] | 650 | // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) |
---|
[819] | 651 | G4bool kfl=true; // Flag of K0/aK0 oscillation |
---|
| 652 | G4bool kf=false; |
---|
| 653 | if(PDG==130||PDG==310) |
---|
| 654 | { |
---|
| 655 | kf=true; |
---|
| 656 | if(G4UniformRand()>.5) kfl=false; |
---|
| 657 | } |
---|
| 658 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) { |
---|
| 659 | ind=0; // pp/nn |
---|
| 660 | } else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) { |
---|
| 661 | ind=1; // np/pn |
---|
| 662 | } else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) { |
---|
| 663 | ind=2; // pimp/pipn |
---|
| 664 | } else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) { |
---|
| 665 | ind=3; // pipp/pimn |
---|
| 666 | } else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) { |
---|
| 667 | ind=4; // KmN/K0N |
---|
| 668 | } else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) { |
---|
| 669 | ind=5; // KpN/aK0N |
---|
| 670 | } else if ( PDG > 3000 && PDG < 3335) { |
---|
| 671 | ind=6; // @@ for all hyperons - take Lambda |
---|
| 672 | } else if ( PDG < -2000 && PDG > -3335 ) { |
---|
| 673 | ind=7; // @@ for all anti-baryons - anti-p/anti-n |
---|
| 674 | } else { |
---|
| 675 | G4cout<<"*Error*G4QuasiFreeRatios::FetchElTot: PDG="<<PDG |
---|
| 676 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
| 677 | G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash"); |
---|
| 678 | } |
---|
| 679 | if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler |
---|
| 680 | // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
| 681 | if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?) |
---|
| 682 | G4bool found=false; |
---|
| 683 | G4int i=-1; |
---|
[1055] | 684 | if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i]) // Sirch for this index in AMDB |
---|
[819] | 685 | { |
---|
| 686 | found=true; // The index is found |
---|
| 687 | break; |
---|
| 688 | } |
---|
| 689 | G4double lp=std::log(p); |
---|
| 690 | #ifdef pdebug |
---|
[1055] | 691 | G4cout<<"G4QuasiFreeR::FetchElTot:I="<<ind<<",i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl; |
---|
[819] | 692 | #endif |
---|
| 693 | if(!nDB || !found) // Create new line in the AMDB |
---|
[1055] | 694 | { |
---|
[819] | 695 | #ifdef pdebug |
---|
[1055] | 696 | G4cout<<"G4QuasiFreeRatios::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl; |
---|
[819] | 697 | #endif |
---|
| 698 | lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot |
---|
| 699 | lastI = ind; // Remember the initialized inex |
---|
| 700 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
| 701 | if(lastK>nlp) |
---|
| 702 | { |
---|
| 703 | lastK=nlp; |
---|
| 704 | lastM=lpa-lpi; |
---|
| 705 | } |
---|
| 706 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
| 707 | G4double pv=mi; |
---|
| 708 | for(G4int j=0; j<=lastK; j++) // Calculate LogTab values |
---|
| 709 | { |
---|
| 710 | lastX[j]=CalcElTot(pv,ind); |
---|
| 711 | #ifdef pdebug |
---|
[1055] | 712 | G4cout<<"G4QuasiFreeR::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
[819] | 713 | <<lastX[j].second<<G4endl; |
---|
| 714 | #endif |
---|
| 715 | if(j!=lastK) pv*=edl; |
---|
| 716 | } |
---|
| 717 | i++; // Make a new record to AMDB and position on it |
---|
| 718 | vI.push_back(lastI); |
---|
| 719 | vM.push_back(lastM); |
---|
| 720 | vK.push_back(lastK); |
---|
| 721 | vX.push_back(lastX); |
---|
| 722 | } |
---|
| 723 | else // The A value was found in AMDB |
---|
[1055] | 724 | { |
---|
[819] | 725 | lastI=vI[i]; |
---|
| 726 | lastM=vM[i]; |
---|
| 727 | lastK=vK[i]; |
---|
| 728 | lastX=vX[i]; |
---|
| 729 | G4int nextK=lastK+1; |
---|
| 730 | G4double lpM=lastM+lpi; |
---|
| 731 | #ifdef pdebug |
---|
[1055] | 732 | G4cout<<"G4QuasiFreeR::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl; |
---|
[819] | 733 | #endif |
---|
| 734 | if(lp>lpM && lastK<nlp) // LogTab must be updated |
---|
[1055] | 735 | { |
---|
[819] | 736 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab |
---|
| 737 | #ifdef pdebug |
---|
[1055] | 738 | G4cout<<"G4QuasiFreeR::FetET:K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl; |
---|
[819] | 739 | #endif |
---|
| 740 | if(lastK>nlp) |
---|
| 741 | { |
---|
| 742 | lastK=nlp; |
---|
| 743 | lastM=lpa-lpi; |
---|
| 744 | } |
---|
| 745 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
| 746 | G4double pv=std::exp(lpM); // momentum of the last calculated beam |
---|
| 747 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
| 748 | { |
---|
| 749 | pv*=edl; |
---|
| 750 | lastX[j]=CalcElTot(pv,ind); |
---|
| 751 | #ifdef pdebug |
---|
[1055] | 752 | G4cout<<"G4QuasiFreeR::FetchElTot:U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
[819] | 753 | <<lastX[j].second<<G4endl; |
---|
| 754 | #endif |
---|
| 755 | } |
---|
| 756 | } // End of LogTab update |
---|
| 757 | if(lastK>=nextK) // The AMDB was apdated |
---|
[1055] | 758 | { |
---|
[819] | 759 | vM[i]=lastM; |
---|
| 760 | vK[i]=lastK; |
---|
| 761 | } |
---|
| 762 | } |
---|
| 763 | // Now one can use tabeles to calculate the value |
---|
| 764 | G4double dlp=lp-lpi; // Shifted log(p) value |
---|
| 765 | G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin |
---|
| 766 | G4double d=dlp-n*dl; // Log shift |
---|
| 767 | G4double e=lastX[n].first; // E-Base |
---|
| 768 | lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result |
---|
| 769 | if(lastR.first<0.) lastR.first = 0.; |
---|
| 770 | G4double t=lastX[n].second; // T-Base |
---|
| 771 | lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result |
---|
| 772 | if(lastR.second<0.) lastR.second= 0.; |
---|
[1055] | 773 | #ifdef pdebug |
---|
| 774 | G4cout<<"=O=>G4QuasiFreeR::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl; |
---|
| 775 | #endif |
---|
[819] | 776 | if(lastR.first>lastR.second) lastR.first = lastR.second; |
---|
| 777 | return lastR; |
---|
| 778 | } // End of FetchElTot |
---|
| 779 | |
---|
| 780 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
| 781 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG, |
---|
| 782 | G4int Z, G4int N) |
---|
| 783 | { |
---|
| 784 | G4double pGeV=pIU/gigaelectronvolt; |
---|
[1055] | 785 | #ifdef pdebug |
---|
| 786 | G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl; |
---|
| 787 | #endif |
---|
[819] | 788 | if(Z<1 && N<1) |
---|
| 789 | { |
---|
| 790 | G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
| 791 | return std::make_pair(0.,0.); |
---|
| 792 | } |
---|
| 793 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
| 794 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
[1055] | 795 | #ifdef pdebug |
---|
| 796 | G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<"," |
---|
| 797 | <<hn.second<<")"<<G4endl; |
---|
| 798 | #endif |
---|
[819] | 799 | G4double A=(Z+N)/millibarn; // To make the result in independent units(IU) |
---|
| 800 | return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A); |
---|
| 801 | } // End of GetElTot |
---|
| 802 | |
---|
| 803 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
| 804 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG, |
---|
| 805 | G4int Z, G4int N) |
---|
| 806 | { |
---|
| 807 | G4double pGeV=pIU/gigaelectronvolt; |
---|
| 808 | G4double resP=0.; |
---|
| 809 | G4double resN=0.; |
---|
| 810 | if(Z<1 && N<1) |
---|
| 811 | { |
---|
| 812 | G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
| 813 | return std::make_pair(resP,resN); |
---|
| 814 | } |
---|
| 815 | G4double A=Z+N; |
---|
| 816 | G4double pf=0.; // Possibility to interact with a proton |
---|
| 817 | G4double nf=0.; // Possibility to interact with a neutron |
---|
| 818 | if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N); |
---|
| 819 | else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z); |
---|
| 820 | else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310) |
---|
| 821 | { |
---|
| 822 | G4double dA=A+A; |
---|
| 823 | pf=Z/(dA+N+N); |
---|
| 824 | nf=N/(dA+Z+Z); |
---|
| 825 | } |
---|
| 826 | G4double mult=1.; // Factor of increasing multiplicity ( ? @@) |
---|
| 827 | if(pGeV>.5) |
---|
| 828 | { |
---|
| 829 | mult=1./(1.+std::log(pGeV+pGeV))/pGeV; |
---|
| 830 | if(mult>1.) mult=1.; |
---|
| 831 | } |
---|
| 832 | if(pf) |
---|
| 833 | { |
---|
| 834 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
| 835 | resP=pf*(hp.second/hp.first-1.)*mult; |
---|
| 836 | } |
---|
| 837 | if(nf) |
---|
| 838 | { |
---|
| 839 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
| 840 | resN=nf*(hn.second/hn.first-1.)*mult; |
---|
| 841 | } |
---|
| 842 | return std::make_pair(resP,resN); |
---|
| 843 | } // End of GetChExFactor |
---|
| 844 | |
---|
| 845 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
| 846 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
| 847 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG, |
---|
| 848 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
| 849 | { |
---|
| 850 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
| 851 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
| 852 | static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron |
---|
| 853 | static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium |
---|
| 854 | static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3 |
---|
| 855 | static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha |
---|
| 856 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) |
---|
| 857 | N4M/=megaelectronvolt; |
---|
| 858 | G4LorentzVector tot4M=N4M+p4M; |
---|
[1055] | 859 | #ifdef ppdebug |
---|
| 860 | G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl; |
---|
| 861 | #endif |
---|
[819] | 862 | G4double mT=mNeut; |
---|
| 863 | G4int Z=0; |
---|
| 864 | G4int N=1; |
---|
| 865 | if(NPDG==2212||NPDG==90001000) |
---|
| 866 | { |
---|
| 867 | mT=mProt; |
---|
| 868 | Z=1; |
---|
| 869 | N=0; |
---|
| 870 | } |
---|
| 871 | else if(NPDG==90001001) |
---|
| 872 | { |
---|
| 873 | mT=mDeut; |
---|
| 874 | Z=1; |
---|
| 875 | N=1; |
---|
| 876 | } |
---|
| 877 | else if(NPDG==90002001) |
---|
| 878 | { |
---|
| 879 | mT=mHel3; |
---|
| 880 | Z=2; |
---|
| 881 | N=1; |
---|
| 882 | } |
---|
| 883 | else if(NPDG==90001002) |
---|
| 884 | { |
---|
| 885 | mT=mTrit; |
---|
| 886 | Z=1; |
---|
| 887 | N=2; |
---|
| 888 | } |
---|
| 889 | else if(NPDG==90002002) |
---|
| 890 | { |
---|
| 891 | mT=mAlph; |
---|
| 892 | Z=2; |
---|
| 893 | N=2; |
---|
| 894 | } |
---|
| 895 | else if(NPDG!=2112&&NPDG!=90000001) |
---|
| 896 | { |
---|
| 897 | G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
| 898 | G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain"); |
---|
| 899 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 900 | } |
---|
| 901 | G4double mT2=mT*mT; |
---|
| 902 | G4double mP2=pr4M.m2(); |
---|
| 903 | G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); |
---|
| 904 | #ifdef pdebug |
---|
| 905 | G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl; |
---|
| 906 | #endif |
---|
| 907 | G4double E2=E*E; |
---|
| 908 | if(E<0. || E2<mP2) |
---|
| 909 | { |
---|
[1055] | 910 | #ifdef ppdebug |
---|
[819] | 911 | G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; |
---|
| 912 | #endif |
---|
| 913 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 914 | } |
---|
[1055] | 915 | G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system |
---|
[819] | 916 | G4VQCrossSection* CSmanager=G4QElasticCrossSection::GetPointer(); |
---|
[1055] | 917 | #ifdef ppdebug |
---|
[819] | 918 | G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
| 919 | #endif |
---|
| 920 | // @@ Temporary NN t-dependence for all hadrons |
---|
| 921 | if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QElast::Scatter: pPDG="<<pPDG<<G4endl; |
---|
| 922 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
| 923 | if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG |
---|
| 924 | G4double xSec=CSmanager->GetCrossSection(false, P, Z, N, PDG); // Rec.CrossSect *TMP* |
---|
| 925 | //G4double xSec=CSmanager->GetCrossSection(false, P, Z, N, pPDG); // Rec.CrossSect |
---|
[1055] | 926 | #ifdef ppdebug |
---|
[819] | 927 | G4cout<<"G4QElast::Scatter:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
| 928 | #endif |
---|
| 929 | #ifdef nandebug |
---|
| 930 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 931 | else G4cout<<"*Warning*G4QElast::Scatter: xSec="<<xSec/millibarn<<G4endl; |
---|
| 932 | #endif |
---|
| 933 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
| 934 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
| 935 | { |
---|
[1055] | 936 | #ifdef ppdebug |
---|
[819] | 937 | G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
| 938 | #endif |
---|
| 939 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
| 940 | } |
---|
| 941 | G4double mint=CSmanager->GetExchangeT(Z,N,PDG); // functional randomized -t (MeV^2) *TMP* |
---|
| 942 | //G4double mint=CSmanager->GetExchangeT(Z,N,pPDG); // functional randomized -t in MeV^2 |
---|
| 943 | G4double maxt=CSmanager->GetHMaxT(); // max possible -t |
---|
| 944 | #ifdef ppdebug |
---|
[1055] | 945 | G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z=" |
---|
| 946 | <<Z<<",N="<<N<<G4endl; |
---|
[819] | 947 | #endif |
---|
| 948 | #ifdef nandebug |
---|
| 949 | if(mint>-.0000001); |
---|
| 950 | else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl; |
---|
| 951 | #endif |
---|
| 952 | G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS |
---|
| 953 | #ifdef ppdebug |
---|
| 954 | G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl; |
---|
| 955 | #endif |
---|
| 956 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
| 957 | { |
---|
| 958 | if (cost>1.) cost=1.; |
---|
| 959 | else if(cost<-1.) cost=-1.; |
---|
| 960 | else |
---|
[1055] | 961 | { |
---|
[819] | 962 | G4cerr<<"G4QFR::S:*NAN*c="<<cost<<",t="<<mint<<",tm="<<CSmanager->GetHMaxT()<<G4endl; |
---|
| 963 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 964 | } |
---|
| 965 | } |
---|
| 966 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
| 967 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
| 968 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
| 969 | { |
---|
| 970 | G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl; |
---|
| 971 | //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp"); |
---|
| 972 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 973 | } |
---|
| 974 | #ifdef ppdebug |
---|
[1055] | 975 | G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl; |
---|
[819] | 976 | #endif |
---|
[1055] | 977 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
[819] | 978 | } // End of Scatter |
---|
| 979 | |
---|
| 980 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
| 981 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
| 982 | // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.) |
---|
| 983 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG, |
---|
| 984 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
| 985 | { |
---|
| 986 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
| 987 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
| 988 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M) |
---|
| 989 | N4M/=megaelectronvolt; |
---|
| 990 | G4LorentzVector tot4M=N4M+p4M; |
---|
| 991 | G4int Z=0; |
---|
| 992 | G4int N=1; |
---|
| 993 | G4int RPDG=2212; // PDG of the recoil nucleon |
---|
| 994 | G4int sPDG=0; // PDG code of the scattered hadron |
---|
| 995 | G4double mS=0.; // proto of mass of scattered hadron |
---|
| 996 | G4double mT=mProt; // mass of the recoil nucleon |
---|
| 997 | if(NPDG==2212) |
---|
| 998 | { |
---|
| 999 | mT=mNeut; |
---|
| 1000 | Z=1; |
---|
| 1001 | N=0; |
---|
| 1002 | RPDG=2112; // PDG of the recoil nucleon |
---|
[1055] | 1003 | if(pPDG==-211) sPDG=111; // pi+ -> pi0 |
---|
| 1004 | else if(pPDG==-321) |
---|
[819] | 1005 | { |
---|
| 1006 | sPDG=310; // K+ -> K0S |
---|
| 1007 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
| 1008 | } |
---|
| 1009 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?) |
---|
| 1010 | else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0 |
---|
| 1011 | else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+ |
---|
| 1012 | else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0 |
---|
| 1013 | } |
---|
| 1014 | else if(NPDG==2112) // Default |
---|
| 1015 | { |
---|
[1055] | 1016 | if(pPDG==211) sPDG=111; // pi+ -> pi0 |
---|
| 1017 | else if(pPDG==321) |
---|
[819] | 1018 | { |
---|
| 1019 | sPDG=310; // K+ -> K0S |
---|
| 1020 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
| 1021 | } |
---|
| 1022 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?) |
---|
| 1023 | else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0 |
---|
| 1024 | else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma- |
---|
| 1025 | else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi- |
---|
| 1026 | } |
---|
| 1027 | else |
---|
| 1028 | { |
---|
| 1029 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
| 1030 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
| 1031 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 1032 | } |
---|
| 1033 | if(sPDG) mS=G4QPDGCode(2112).GetMass(); |
---|
| 1034 | else |
---|
| 1035 | { |
---|
| 1036 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl; |
---|
| 1037 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
| 1038 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 1039 | } |
---|
| 1040 | G4double mT2=mT*mT; |
---|
| 1041 | G4double mS2=mS*mS; |
---|
| 1042 | G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT); |
---|
| 1043 | G4double E2=E*E; |
---|
| 1044 | if(E<0. || E2<mS2) |
---|
| 1045 | { |
---|
| 1046 | #ifdef pdebug |
---|
| 1047 | G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl; |
---|
| 1048 | #endif |
---|
| 1049 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1050 | } |
---|
[1055] | 1051 | G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system |
---|
[819] | 1052 | G4VQCrossSection* CSmanager=G4QElasticCrossSection::GetPointer(); |
---|
| 1053 | #ifdef debug |
---|
| 1054 | G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
| 1055 | #endif |
---|
| 1056 | // @@ Temporary NN t-dependence for all hadrons |
---|
| 1057 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
| 1058 | G4double xSec=CSmanager->GetCrossSection(false, P, Z, N, PDG); // Rec.CrossSect *TMP* |
---|
| 1059 | //G4double xSec=CSmanager->GetCrossSection(false, P, Z, N, pPDG); // Rec.CrossSect |
---|
| 1060 | #ifdef debug |
---|
| 1061 | G4cout<<"G4QElast::ChExer:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
| 1062 | #endif |
---|
| 1063 | #ifdef nandebug |
---|
| 1064 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 1065 | else G4cout<<"*Warning*G4QElast::ChExer: xSec="<<xSec/millibarn<<G4endl; |
---|
| 1066 | #endif |
---|
| 1067 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
| 1068 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
| 1069 | { |
---|
| 1070 | #ifdef pdebug |
---|
| 1071 | G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
| 1072 | #endif |
---|
| 1073 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
| 1074 | } |
---|
| 1075 | G4double mint=CSmanager->GetExchangeT(Z,N,PDG); // functional randomized -t (MeV^2) *TMP* |
---|
| 1076 | //G4double mint=CSmanager->GetExchangeT(Z,N,pPDG); // functional randomized -t in MeV^2 |
---|
| 1077 | #ifdef pdebug |
---|
| 1078 | G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl; |
---|
| 1079 | #endif |
---|
| 1080 | #ifdef nandebug |
---|
| 1081 | if(mint>-.0000001); |
---|
| 1082 | else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl; |
---|
| 1083 | #endif |
---|
| 1084 | G4double cost=1.-mint/CSmanager->GetHMaxT(); // cos(theta) in CMS |
---|
| 1085 | #ifdef pdebug |
---|
| 1086 | G4cout<<"G4QFR::ChEx:-t="<<mint<<",dpc2="<<CSmanager->GetHMaxT()<<",cost="<<cost<<G4endl; |
---|
| 1087 | #endif |
---|
| 1088 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
| 1089 | { |
---|
| 1090 | if (cost>1.) cost=1.; |
---|
| 1091 | else if(cost<-1.) cost=-1.; |
---|
| 1092 | else |
---|
[1055] | 1093 | { |
---|
[819] | 1094 | G4cerr<<"G4QFR::C:*NAN*c="<<cost<<",t="<<mint<<",tm="<<CSmanager->GetHMaxT()<<G4endl; |
---|
| 1095 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1096 | } |
---|
| 1097 | } |
---|
| 1098 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
| 1099 | pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron |
---|
| 1100 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
| 1101 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
| 1102 | { |
---|
| 1103 | G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl; |
---|
| 1104 | //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp"); |
---|
| 1105 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1106 | } |
---|
| 1107 | #ifdef debug |
---|
[1055] | 1108 | G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl; |
---|
[819] | 1109 | #endif |
---|
[1055] | 1110 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
[819] | 1111 | } // End of ChExer |
---|
| 1112 | |
---|
| 1113 | // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile) |
---|
| 1114 | G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) |
---|
| 1115 | { |
---|
| 1116 | p/=MeV; // Converted from independent units |
---|
| 1117 | G4double A=Z+N; |
---|
| 1118 | if(A<1.5) return 0.; |
---|
| 1119 | G4double C=0.; |
---|
| 1120 | if (pPDG==2212) C=N/(A+Z); |
---|
| 1121 | else if(pPDG==2112) C=Z/(A+N); |
---|
| 1122 | else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl; |
---|
| 1123 | C*=C; // Coherent processes squares the amplitude |
---|
| 1124 | // @@ This is true only for nucleons: other projectiles must be treated differently |
---|
| 1125 | G4double sp=std::sqrt(p); |
---|
| 1126 | G4double p2=p*p; |
---|
| 1127 | G4double p4=p2*p2; |
---|
[1055] | 1128 | G4double dl1=std::log(p)-5.; |
---|
[819] | 1129 | G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013); |
---|
| 1130 | G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; |
---|
| 1131 | G4double R=U/T; |
---|
| 1132 | return C*R*R; |
---|
| 1133 | } |
---|