[1199] | 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 | // |
---|
[1315] | 27 | // $Id: G4QuasiFreeRatios.cc,v 1.3 2010/01/22 17:02:48 mkossov Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ |
---|
[1199] | 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 | // |
---|
| 35 | //======================================================================= |
---|
| 36 | // Short description: Provides percentage of quasi-free and quasi-elastic |
---|
| 37 | // reactions in the inelastic reactions. |
---|
| 38 | // ---------------------------------------------------------------------- |
---|
| 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 |
---|
| 55 | G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl; |
---|
| 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 | { |
---|
| 86 | #ifdef pdebug |
---|
| 87 | G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl; |
---|
| 88 | #endif |
---|
| 89 | G4double R=0.; |
---|
| 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.) |
---|
| 97 | { |
---|
| 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 |
---|
| 100 | } |
---|
| 101 | #ifdef pdebug |
---|
| 102 | G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl; |
---|
| 103 | #endif |
---|
| 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 |
---|
| 140 | #ifdef pdebug |
---|
| 141 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s<<G4endl; |
---|
| 142 | #endif |
---|
| 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; |
---|
| 155 | if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB |
---|
| 156 | { |
---|
| 157 | found=true; // The A value is found |
---|
| 158 | break; |
---|
| 159 | } |
---|
| 160 | #ifdef pdebug |
---|
| 161 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl; |
---|
| 162 | #endif |
---|
| 163 | if(!nDB || !found) // Create new line in the AMDB |
---|
| 164 | { |
---|
| 165 | lastA = A; |
---|
| 166 | #ifdef pdebug |
---|
| 167 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
| 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 | } |
---|
| 184 | lastL=new G4double[mls]; // Create the logarithmic Table |
---|
| 185 | if(s>sma) // Initialize the logarithmic Table |
---|
| 186 | { |
---|
| 187 | #ifdef pdebug |
---|
| 188 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
| 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); |
---|
| 202 | if(j!=lastK) sv*=edl; |
---|
| 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 |
---|
| 220 | { |
---|
| 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]; |
---|
| 228 | #ifdef pdebug |
---|
| 229 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s<<", lastM="<<lastM<<G4endl; |
---|
| 230 | #endif |
---|
| 231 | if(s>lastM) // At least LinTab must be updated |
---|
| 232 | { |
---|
| 233 | G4int nextN=lastN+1; // The next bin to be initialized |
---|
| 234 | #ifdef pdebug |
---|
| 235 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl; |
---|
| 236 | #endif |
---|
| 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 | } |
---|
| 252 | #ifdef pdebug |
---|
| 253 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
| 254 | #endif |
---|
| 255 | } // End of LinTab update |
---|
| 256 | #ifdef pdebug |
---|
| 257 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl; |
---|
| 258 | #endif |
---|
| 259 | if(lastN>=nextN) |
---|
| 260 | { |
---|
| 261 | vH[i]=lastH; |
---|
| 262 | vN[i]=lastN; |
---|
| 263 | } |
---|
| 264 | G4int nextK=lastK+1; |
---|
| 265 | if(!lastK) nextK=0; |
---|
| 266 | #ifdef pdebug |
---|
| 267 | G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl; |
---|
| 268 | #endif |
---|
| 269 | if(s>sma && lastK<nls) // LogTab must be updated |
---|
| 270 | { |
---|
| 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 |
---|
| 280 | #ifdef pdebug |
---|
| 281 | G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl; |
---|
| 282 | #endif |
---|
| 283 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
| 284 | { |
---|
| 285 | sv*=edl; |
---|
| 286 | #ifdef pdebug |
---|
| 287 | G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl; |
---|
| 288 | #endif |
---|
| 289 | lastL[j]=CalcQF2IN_Ratio(sv,A); |
---|
| 290 | } |
---|
| 291 | #ifdef pdebug |
---|
| 292 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
| 293 | #endif |
---|
| 294 | } // End of LogTab update |
---|
| 295 | #ifdef pdebug |
---|
| 296 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl; |
---|
| 297 | #endif |
---|
| 298 | if(lastK>=nextK) |
---|
| 299 | { |
---|
| 300 | vM[i]=lastM; |
---|
| 301 | vK[i]=lastK; |
---|
| 302 | } |
---|
| 303 | } |
---|
| 304 | } |
---|
| 305 | #ifdef pdebug |
---|
| 306 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s<<", sma="<<sma<<G4endl; |
---|
| 307 | #endif |
---|
| 308 | // Now one can use tabeles to calculate the value |
---|
| 309 | if(s<sma) // Use linear table |
---|
| 310 | { |
---|
| 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 | } |
---|
| 316 | else // Use log table |
---|
| 317 | { |
---|
| 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.; |
---|
| 326 | #ifdef pdebug |
---|
| 327 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl; |
---|
| 328 | #endif |
---|
| 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; |
---|
| 336 | G4double s2=s*s; |
---|
| 337 | G4double s4=s2*s2; |
---|
| 338 | G4double ss=std::sqrt(std::sqrt(s)); |
---|
| 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); |
---|
| 342 | return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P); |
---|
| 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 |
---|
| 363 | { |
---|
| 364 | #ifdef debug |
---|
| 365 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl; |
---|
| 366 | #endif |
---|
| 367 | if(p<pmi) |
---|
| 368 | { |
---|
| 369 | G4double p2=p*p; |
---|
| 370 | El=1./(.00012+p2*.2); |
---|
| 371 | To=El; |
---|
| 372 | #ifdef debug |
---|
| 373 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl; |
---|
| 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 |
---|
| 383 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl; |
---|
| 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 |
---|
| 396 | G4cout<<"G4QuasiFreeR::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl; |
---|
| 397 | #endif |
---|
| 398 | } |
---|
| 399 | } |
---|
| 400 | else if(I==1) // np/pn |
---|
| 401 | { |
---|
| 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 |
---|
| 427 | { |
---|
| 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 | { |
---|
[1315] | 445 | G4double lr=lp+1.27; // p1 |
---|
| 446 | G4double LE=1.53/(lr*lr+.0676); // p2, p3 |
---|
| 447 | G4double ld=lp-lmi; // p4 (lmi=3.5) |
---|
[1199] | 448 | G4double ld2=ld*ld; |
---|
| 449 | G4double p2=p*p; |
---|
| 450 | G4double p4=p2*p2; |
---|
| 451 | G4double sp=std::sqrt(p); |
---|
[1315] | 452 | G4double lm=lp+.36; // p5 |
---|
| 453 | G4double md=lm*lm+.04; // p6 |
---|
| 454 | G4double lh=lp-.017; // p7 |
---|
| 455 | G4double hd=lh*lh+.0025; // p8 |
---|
| 456 | El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14 |
---|
[1199] | 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 |
---|
| 461 | { |
---|
| 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 | { |
---|
[1315] | 480 | G4double lr=lp+1.27; // p1 |
---|
[1199] | 481 | G4double lr2=lr*lr; |
---|
[1315] | 482 | G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3 |
---|
| 483 | G4double ld=lp-lmi; // p4 (lmi=3.5) |
---|
[1199] | 484 | G4double ld2=ld*ld; |
---|
| 485 | G4double p2=p*p; |
---|
| 486 | G4double p4=p2*p2; |
---|
| 487 | G4double sp=std::sqrt(p); |
---|
[1315] | 488 | G4double lm=lp-.32; // p5 |
---|
| 489 | G4double md=lm*lm+.0576; // p6 |
---|
| 490 | El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11 |
---|
[1199] | 491 | To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md; |
---|
| 492 | } |
---|
| 493 | } |
---|
| 494 | else if(I==4) // Kmp/Kmn/K0p/K0n |
---|
| 495 | { |
---|
| 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 |
---|
| 527 | { |
---|
| 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 |
---|
| 559 | { |
---|
| 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 |
---|
| 587 | { |
---|
| 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 | // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb) |
---|
| 615 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTotXS(G4double p, G4int PDG, G4bool F) |
---|
| 616 | { |
---|
| 617 | G4int ind=0; // Prototype of the reaction index |
---|
| 618 | G4bool kfl=true; // Flag of K0/aK0 oscillation |
---|
| 619 | G4bool kf=false; |
---|
| 620 | if(PDG==130||PDG==310) |
---|
| 621 | { |
---|
| 622 | kf=true; |
---|
| 623 | if(G4UniformRand()>.5) kfl=false; |
---|
| 624 | } |
---|
| 625 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn |
---|
| 626 | else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn |
---|
| 627 | else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn |
---|
| 628 | else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn |
---|
| 629 | else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N |
---|
| 630 | else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N |
---|
| 631 | else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda |
---|
| 632 | else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) |
---|
| 633 | else { |
---|
| 634 | G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG |
---|
| 635 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
| 636 | G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash"); |
---|
| 637 | } |
---|
| 638 | return CalcElTot(p,ind); |
---|
| 639 | } |
---|
| 640 | |
---|
| 641 | // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron |
---|
| 642 | std::pair<G4double,G4double> G4QuasiFreeRatios::FetchElTot(G4double p, G4int PDG, G4bool F) |
---|
| 643 | { |
---|
| 644 | static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step) |
---|
| 645 | static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable |
---|
| 646 | static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
| 647 | static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
| 648 | static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c) |
---|
| 649 | static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV) |
---|
| 650 | static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table |
---|
| 651 | static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table |
---|
| 652 | //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum" |
---|
| 653 | static G4double lastP=0.; // The last momentum for which XS was calculated |
---|
| 654 | static G4int lastH=0; // The last projPDG for which XS was calculated |
---|
| 655 | static G4bool lastF=true; // The last nucleon for which XS was calculated |
---|
| 656 | static std::pair<G4double,G4double> lastR=std::make_pair(0.,0.); // The last result |
---|
| 657 | // Local Associative Data Base: |
---|
| 658 | static std::vector<G4int> vI; // Vector of index for which XS was calculated |
---|
| 659 | static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable |
---|
| 660 | static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable |
---|
| 661 | // Last values of the Associative Data Base: |
---|
| 662 | static G4int lastI=0; // The Last index for which XS was calculated |
---|
| 663 | static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable |
---|
| 664 | static G4int lastK=0; // The Last topBin number initialized in LogTable |
---|
| 665 | static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap |
---|
| 666 | // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei |
---|
| 667 | G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB |
---|
| 668 | #ifdef pdebug |
---|
| 669 | G4cout<<"G4QuasiFreeR::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl; |
---|
| 670 | #endif |
---|
| 671 | if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler. |
---|
| 672 | // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
| 673 | lastH=PDG; |
---|
| 674 | lastF=F; |
---|
| 675 | G4int ind=-1; // Prototipe of the index of the PDG/F combination |
---|
| 676 | // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p), |
---|
| 677 | // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) |
---|
| 678 | G4bool kfl=true; // Flag of K0/aK0 oscillation |
---|
| 679 | G4bool kf=false; |
---|
| 680 | if(PDG==130||PDG==310) |
---|
| 681 | { |
---|
| 682 | kf=true; |
---|
| 683 | if(G4UniformRand()>.5) kfl=false; |
---|
| 684 | } |
---|
| 685 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn |
---|
| 686 | else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn |
---|
| 687 | else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn |
---|
| 688 | else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn |
---|
| 689 | else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N |
---|
| 690 | else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N |
---|
| 691 | else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda |
---|
| 692 | else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) |
---|
| 693 | else { |
---|
| 694 | G4cout<<"*Error*G4QuasiFreeRatios::FetchElTot: PDG="<<PDG |
---|
| 695 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
| 696 | G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash"); |
---|
| 697 | } |
---|
| 698 | if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler |
---|
| 699 | // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
| 700 | if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?) |
---|
| 701 | G4bool found=false; |
---|
| 702 | G4int i=-1; |
---|
| 703 | if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i]) // Sirch for this index in AMDB |
---|
| 704 | { |
---|
| 705 | found=true; // The index is found |
---|
| 706 | break; |
---|
| 707 | } |
---|
| 708 | G4double lp=std::log(p); |
---|
| 709 | #ifdef pdebug |
---|
| 710 | G4cout<<"G4QuasiFreeR::FetchElTot:I="<<ind<<",i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl; |
---|
| 711 | #endif |
---|
| 712 | if(!nDB || !found) // Create new line in the AMDB |
---|
| 713 | { |
---|
| 714 | #ifdef pdebug |
---|
| 715 | G4cout<<"G4QuasiFreeRatios::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl; |
---|
| 716 | #endif |
---|
| 717 | lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot |
---|
| 718 | lastI = ind; // Remember the initialized inex |
---|
| 719 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
| 720 | if(lastK>nlp) |
---|
| 721 | { |
---|
| 722 | lastK=nlp; |
---|
| 723 | lastM=lpa-lpi; |
---|
| 724 | } |
---|
| 725 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
| 726 | G4double pv=mi; |
---|
| 727 | for(G4int j=0; j<=lastK; j++) // Calculate LogTab values |
---|
| 728 | { |
---|
| 729 | lastX[j]=CalcElTot(pv,ind); |
---|
| 730 | #ifdef pdebug |
---|
| 731 | G4cout<<"G4QuasiFreeR::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
| 732 | <<lastX[j].second<<G4endl; |
---|
| 733 | #endif |
---|
| 734 | if(j!=lastK) pv*=edl; |
---|
| 735 | } |
---|
| 736 | i++; // Make a new record to AMDB and position on it |
---|
| 737 | vI.push_back(lastI); |
---|
| 738 | vM.push_back(lastM); |
---|
| 739 | vK.push_back(lastK); |
---|
| 740 | vX.push_back(lastX); |
---|
| 741 | } |
---|
| 742 | else // The A value was found in AMDB |
---|
| 743 | { |
---|
| 744 | lastI=vI[i]; |
---|
| 745 | lastM=vM[i]; |
---|
| 746 | lastK=vK[i]; |
---|
| 747 | lastX=vX[i]; |
---|
| 748 | G4int nextK=lastK+1; |
---|
| 749 | G4double lpM=lastM+lpi; |
---|
| 750 | #ifdef pdebug |
---|
| 751 | G4cout<<"G4QuasiFreeR::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl; |
---|
| 752 | #endif |
---|
| 753 | if(lp>lpM && lastK<nlp) // LogTab must be updated |
---|
| 754 | { |
---|
| 755 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab |
---|
| 756 | #ifdef pdebug |
---|
| 757 | G4cout<<"G4QuasiFreeR::FetET:K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl; |
---|
| 758 | #endif |
---|
| 759 | if(lastK>nlp) |
---|
| 760 | { |
---|
| 761 | lastK=nlp; |
---|
| 762 | lastM=lpa-lpi; |
---|
| 763 | } |
---|
| 764 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
| 765 | G4double pv=std::exp(lpM); // momentum of the last calculated beam |
---|
| 766 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
| 767 | { |
---|
| 768 | pv*=edl; |
---|
| 769 | lastX[j]=CalcElTot(pv,ind); |
---|
| 770 | #ifdef pdebug |
---|
| 771 | G4cout<<"G4QuasiFreeR::FetchElTot:U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
| 772 | <<lastX[j].second<<G4endl; |
---|
| 773 | #endif |
---|
| 774 | } |
---|
| 775 | } // End of LogTab update |
---|
| 776 | if(lastK>=nextK) // The AMDB was apdated |
---|
| 777 | { |
---|
| 778 | vM[i]=lastM; |
---|
| 779 | vK[i]=lastK; |
---|
| 780 | } |
---|
| 781 | } |
---|
| 782 | // Now one can use tabeles to calculate the value |
---|
| 783 | G4double dlp=lp-lpi; // Shifted log(p) value |
---|
| 784 | G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin |
---|
| 785 | G4double d=dlp-n*dl; // Log shift |
---|
| 786 | G4double e=lastX[n].first; // E-Base |
---|
| 787 | lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result |
---|
| 788 | if(lastR.first<0.) lastR.first = 0.; |
---|
| 789 | G4double t=lastX[n].second; // T-Base |
---|
| 790 | lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result |
---|
| 791 | if(lastR.second<0.) lastR.second= 0.; |
---|
| 792 | #ifdef pdebug |
---|
| 793 | G4cout<<"=O=>G4QuasiFreeR::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl; |
---|
| 794 | #endif |
---|
| 795 | if(lastR.first>lastR.second) lastR.first = lastR.second; |
---|
| 796 | return lastR; |
---|
| 797 | } // End of FetchElTot |
---|
| 798 | |
---|
| 799 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
| 800 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG, |
---|
| 801 | G4int Z, G4int N) |
---|
| 802 | { |
---|
| 803 | G4double pGeV=pIU/gigaelectronvolt; |
---|
| 804 | #ifdef pdebug |
---|
| 805 | G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl; |
---|
| 806 | #endif |
---|
| 807 | if(Z<1 && N<1) |
---|
| 808 | { |
---|
| 809 | G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
| 810 | return std::make_pair(0.,0.); |
---|
| 811 | } |
---|
| 812 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
| 813 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
| 814 | #ifdef pdebug |
---|
| 815 | G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<"," |
---|
| 816 | <<hn.second<<")"<<G4endl; |
---|
| 817 | #endif |
---|
| 818 | G4double A=(Z+N)/millibarn; // To make the result in independent units(IU) |
---|
| 819 | return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A); |
---|
| 820 | } // End of GetElTot |
---|
| 821 | |
---|
| 822 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
| 823 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG, |
---|
| 824 | G4int Z, G4int N) |
---|
| 825 | { |
---|
| 826 | G4double pGeV=pIU/gigaelectronvolt; |
---|
| 827 | G4double resP=0.; |
---|
| 828 | G4double resN=0.; |
---|
| 829 | if(Z<1 && N<1) |
---|
| 830 | { |
---|
| 831 | G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
| 832 | return std::make_pair(resP,resN); |
---|
| 833 | } |
---|
| 834 | G4double A=Z+N; |
---|
| 835 | G4double pf=0.; // Possibility to interact with a proton |
---|
| 836 | G4double nf=0.; // Possibility to interact with a neutron |
---|
| 837 | if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N); |
---|
| 838 | else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z); |
---|
| 839 | else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310) |
---|
| 840 | { |
---|
| 841 | G4double dA=A+A; |
---|
| 842 | pf=Z/(dA+N+N); |
---|
| 843 | nf=N/(dA+Z+Z); |
---|
| 844 | } |
---|
| 845 | G4double mult=1.; // Factor of increasing multiplicity ( ? @@) |
---|
| 846 | if(pGeV>.5) |
---|
| 847 | { |
---|
| 848 | mult=1./(1.+std::log(pGeV+pGeV))/pGeV; |
---|
| 849 | if(mult>1.) mult=1.; |
---|
| 850 | } |
---|
| 851 | if(pf) |
---|
| 852 | { |
---|
| 853 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
| 854 | resP=pf*(hp.second/hp.first-1.)*mult; |
---|
| 855 | } |
---|
| 856 | if(nf) |
---|
| 857 | { |
---|
| 858 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
| 859 | resN=nf*(hn.second/hn.first-1.)*mult; |
---|
| 860 | } |
---|
| 861 | return std::make_pair(resP,resN); |
---|
| 862 | } // End of GetChExFactor |
---|
| 863 | |
---|
| 864 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
| 865 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
| 866 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG, |
---|
| 867 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
| 868 | { |
---|
| 869 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
| 870 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
| 871 | static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron |
---|
| 872 | static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium |
---|
| 873 | static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3 |
---|
| 874 | static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha |
---|
| 875 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) |
---|
| 876 | N4M/=megaelectronvolt; |
---|
| 877 | G4LorentzVector tot4M=N4M+p4M; |
---|
| 878 | #ifdef ppdebug |
---|
| 879 | G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl; |
---|
| 880 | #endif |
---|
| 881 | G4double mT=mNeut; |
---|
| 882 | G4int Z=0; |
---|
| 883 | G4int N=1; |
---|
| 884 | if(NPDG==2212||NPDG==90001000) |
---|
| 885 | { |
---|
| 886 | mT=mProt; |
---|
| 887 | Z=1; |
---|
| 888 | N=0; |
---|
| 889 | } |
---|
| 890 | else if(NPDG==90001001) |
---|
| 891 | { |
---|
| 892 | mT=mDeut; |
---|
| 893 | Z=1; |
---|
| 894 | N=1; |
---|
| 895 | } |
---|
| 896 | else if(NPDG==90002001) |
---|
| 897 | { |
---|
| 898 | mT=mHel3; |
---|
| 899 | Z=2; |
---|
| 900 | N=1; |
---|
| 901 | } |
---|
| 902 | else if(NPDG==90001002) |
---|
| 903 | { |
---|
| 904 | mT=mTrit; |
---|
| 905 | Z=1; |
---|
| 906 | N=2; |
---|
| 907 | } |
---|
| 908 | else if(NPDG==90002002) |
---|
| 909 | { |
---|
| 910 | mT=mAlph; |
---|
| 911 | Z=2; |
---|
| 912 | N=2; |
---|
| 913 | } |
---|
| 914 | else if(NPDG!=2112&&NPDG!=90000001) |
---|
| 915 | { |
---|
| 916 | G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
| 917 | G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain"); |
---|
| 918 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 919 | } |
---|
| 920 | G4double mT2=mT*mT; |
---|
| 921 | G4double mP2=pr4M.m2(); |
---|
| 922 | G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); |
---|
| 923 | #ifdef pdebug |
---|
| 924 | G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl; |
---|
| 925 | #endif |
---|
| 926 | G4double E2=E*E; |
---|
| 927 | if(E<0. || E2<mP2) |
---|
| 928 | { |
---|
| 929 | #ifdef ppdebug |
---|
| 930 | G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; |
---|
| 931 | #endif |
---|
| 932 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 933 | } |
---|
| 934 | G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system |
---|
[1315] | 935 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
| 936 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
[1199] | 937 | #ifdef ppdebug |
---|
| 938 | G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
| 939 | #endif |
---|
| 940 | // @@ Temporary NN t-dependence for all hadrons |
---|
| 941 | if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QElast::Scatter: pPDG="<<pPDG<<G4endl; |
---|
| 942 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
| 943 | if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG |
---|
[1315] | 944 | if(!Z && N==1) // Change for Quasi-Elastic on neutron |
---|
| 945 | { |
---|
| 946 | Z=1; |
---|
| 947 | N=0; |
---|
| 948 | if (PDG==2212) PDG=2112; |
---|
| 949 | else if(PDG==2112) PDG=2212; |
---|
| 950 | } |
---|
| 951 | G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* |
---|
| 952 | if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP* |
---|
| 953 | else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP* |
---|
[1199] | 954 | #ifdef ppdebug |
---|
| 955 | G4cout<<"G4QElast::Scatter:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
| 956 | #endif |
---|
| 957 | #ifdef nandebug |
---|
| 958 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 959 | else G4cout<<"*Warning*G4QElast::Scatter: xSec="<<xSec/millibarn<<G4endl; |
---|
| 960 | #endif |
---|
| 961 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
[1315] | 962 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
[1199] | 963 | { |
---|
| 964 | #ifdef ppdebug |
---|
| 965 | G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
| 966 | #endif |
---|
| 967 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
| 968 | } |
---|
[1315] | 969 | G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* |
---|
| 970 | if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* |
---|
| 971 | else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* |
---|
| 972 | G4double maxt=0.; // Prototype of max possible -t |
---|
| 973 | if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t |
---|
| 974 | else maxt=NCSmanager->GetHMaxT(); // max possible -t |
---|
[1199] | 975 | #ifdef ppdebug |
---|
| 976 | G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z=" |
---|
| 977 | <<Z<<",N="<<N<<G4endl; |
---|
| 978 | #endif |
---|
| 979 | #ifdef nandebug |
---|
| 980 | if(mint>-.0000001); |
---|
| 981 | else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl; |
---|
| 982 | #endif |
---|
| 983 | G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS |
---|
| 984 | #ifdef ppdebug |
---|
| 985 | G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl; |
---|
| 986 | #endif |
---|
| 987 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
| 988 | { |
---|
| 989 | if (cost>1.) cost=1.; |
---|
| 990 | else if(cost<-1.) cost=-1.; |
---|
| 991 | else |
---|
| 992 | { |
---|
[1315] | 993 | G4double tm=0.; |
---|
| 994 | if(PDG==2212) tm=PCSmanager->GetHMaxT(); |
---|
| 995 | else tm=NCSmanager->GetHMaxT(); |
---|
| 996 | G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl; |
---|
[1199] | 997 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 998 | } |
---|
| 999 | } |
---|
| 1000 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
| 1001 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
| 1002 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
| 1003 | { |
---|
| 1004 | G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl; |
---|
| 1005 | //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp"); |
---|
| 1006 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1007 | } |
---|
| 1008 | #ifdef ppdebug |
---|
| 1009 | G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl; |
---|
| 1010 | #endif |
---|
| 1011 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
| 1012 | } // End of Scatter |
---|
| 1013 | |
---|
| 1014 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
| 1015 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
| 1016 | // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.) |
---|
| 1017 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG, |
---|
| 1018 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
| 1019 | { |
---|
| 1020 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
| 1021 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
| 1022 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M) |
---|
| 1023 | N4M/=megaelectronvolt; |
---|
| 1024 | G4LorentzVector tot4M=N4M+p4M; |
---|
| 1025 | G4int Z=0; |
---|
| 1026 | G4int N=1; |
---|
| 1027 | G4int RPDG=2212; // PDG of the recoil nucleon |
---|
| 1028 | G4int sPDG=0; // PDG code of the scattered hadron |
---|
| 1029 | G4double mS=0.; // proto of mass of scattered hadron |
---|
| 1030 | G4double mT=mProt; // mass of the recoil nucleon |
---|
| 1031 | if(NPDG==2212) |
---|
| 1032 | { |
---|
| 1033 | mT=mNeut; |
---|
| 1034 | Z=1; |
---|
| 1035 | N=0; |
---|
| 1036 | RPDG=2112; // PDG of the recoil nucleon |
---|
| 1037 | if(pPDG==-211) sPDG=111; // pi+ -> pi0 |
---|
| 1038 | else if(pPDG==-321) |
---|
| 1039 | { |
---|
| 1040 | sPDG=310; // K+ -> K0S |
---|
| 1041 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
| 1042 | } |
---|
| 1043 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?) |
---|
| 1044 | else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0 |
---|
| 1045 | else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+ |
---|
| 1046 | else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0 |
---|
| 1047 | } |
---|
| 1048 | else if(NPDG==2112) // Default |
---|
| 1049 | { |
---|
| 1050 | if(pPDG==211) sPDG=111; // pi+ -> pi0 |
---|
| 1051 | else if(pPDG==321) |
---|
| 1052 | { |
---|
| 1053 | sPDG=310; // K+ -> K0S |
---|
| 1054 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
| 1055 | } |
---|
| 1056 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?) |
---|
| 1057 | else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0 |
---|
| 1058 | else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma- |
---|
| 1059 | else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi- |
---|
| 1060 | } |
---|
| 1061 | else |
---|
| 1062 | { |
---|
| 1063 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
| 1064 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
| 1065 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 1066 | } |
---|
| 1067 | if(sPDG) mS=mNeut; |
---|
| 1068 | else |
---|
| 1069 | { |
---|
| 1070 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl; |
---|
| 1071 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
| 1072 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
| 1073 | } |
---|
| 1074 | G4double mT2=mT*mT; |
---|
| 1075 | G4double mS2=mS*mS; |
---|
| 1076 | G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT); |
---|
| 1077 | G4double E2=E*E; |
---|
| 1078 | if(E<0. || E2<mS2) |
---|
| 1079 | { |
---|
| 1080 | #ifdef pdebug |
---|
| 1081 | G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl; |
---|
| 1082 | #endif |
---|
| 1083 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1084 | } |
---|
| 1085 | G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system |
---|
[1315] | 1086 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
| 1087 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
[1199] | 1088 | #ifdef debug |
---|
| 1089 | G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
| 1090 | #endif |
---|
| 1091 | // @@ Temporary NN t-dependence for all hadrons |
---|
| 1092 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
[1315] | 1093 | if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG |
---|
| 1094 | if(!Z && N==1) // Change for Quasi-Elastic on neutron |
---|
| 1095 | { |
---|
| 1096 | Z=1; |
---|
| 1097 | N=0; |
---|
| 1098 | if (PDG==2212) PDG=2112; |
---|
| 1099 | else if(PDG==2112) PDG=2212; |
---|
| 1100 | } |
---|
| 1101 | G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* |
---|
| 1102 | if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP* |
---|
| 1103 | else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP* |
---|
[1199] | 1104 | #ifdef debug |
---|
| 1105 | G4cout<<"G4QElast::ChExer:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
| 1106 | #endif |
---|
| 1107 | #ifdef nandebug |
---|
| 1108 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 1109 | else G4cout<<"*Warning*G4QElast::ChExer: xSec="<<xSec/millibarn<<G4endl; |
---|
| 1110 | #endif |
---|
| 1111 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
| 1112 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
| 1113 | { |
---|
| 1114 | #ifdef pdebug |
---|
| 1115 | G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
| 1116 | #endif |
---|
| 1117 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
| 1118 | } |
---|
[1315] | 1119 | G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* |
---|
| 1120 | if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* |
---|
| 1121 | else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* |
---|
[1199] | 1122 | #ifdef pdebug |
---|
| 1123 | G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl; |
---|
| 1124 | #endif |
---|
| 1125 | #ifdef nandebug |
---|
| 1126 | if(mint>-.0000001); |
---|
| 1127 | else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl; |
---|
| 1128 | #endif |
---|
[1315] | 1129 | G4double maxt=0.; // Prototype of max possible -t |
---|
| 1130 | if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t |
---|
| 1131 | else maxt=NCSmanager->GetHMaxT(); // max possible -t |
---|
| 1132 | G4double cost=1.-mint/maxt; // cos(theta) in CMS |
---|
[1199] | 1133 | #ifdef pdebug |
---|
[1315] | 1134 | G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl; |
---|
[1199] | 1135 | #endif |
---|
| 1136 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
| 1137 | { |
---|
| 1138 | if (cost>1.) cost=1.; |
---|
| 1139 | else if(cost<-1.) cost=-1.; |
---|
| 1140 | else |
---|
| 1141 | { |
---|
[1315] | 1142 | G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl; |
---|
[1199] | 1143 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1144 | } |
---|
| 1145 | } |
---|
| 1146 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
| 1147 | pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron |
---|
| 1148 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
| 1149 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
| 1150 | { |
---|
| 1151 | G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl; |
---|
| 1152 | //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp"); |
---|
| 1153 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
| 1154 | } |
---|
| 1155 | #ifdef debug |
---|
| 1156 | G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl; |
---|
| 1157 | #endif |
---|
| 1158 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
| 1159 | } // End of ChExer |
---|
| 1160 | |
---|
| 1161 | // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile) |
---|
| 1162 | G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) |
---|
| 1163 | { |
---|
| 1164 | p/=MeV; // Converted from independent units |
---|
| 1165 | G4double A=Z+N; |
---|
| 1166 | if(A<1.5) return 0.; |
---|
| 1167 | G4double C=0.; |
---|
| 1168 | if (pPDG==2212) C=N/(A+Z); |
---|
| 1169 | else if(pPDG==2112) C=Z/(A+N); |
---|
| 1170 | else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl; |
---|
| 1171 | C*=C; // Coherent processes squares the amplitude |
---|
| 1172 | // @@ This is true only for nucleons: other projectiles must be treated differently |
---|
| 1173 | G4double sp=std::sqrt(p); |
---|
| 1174 | G4double p2=p*p; |
---|
| 1175 | G4double p4=p2*p2; |
---|
| 1176 | G4double dl1=std::log(p)-5.; |
---|
| 1177 | G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013); |
---|
| 1178 | G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; |
---|
| 1179 | G4double R=U/T; |
---|
| 1180 | return C*R*R; |
---|
| 1181 | } |
---|