[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 | // hpw: done, but low quality at present. |
---|
| 27 | |
---|
| 28 | #include "globals.hh" |
---|
| 29 | #include "G4AngularDistribution.hh" |
---|
| 30 | #include "Randomize.hh" |
---|
| 31 | |
---|
| 32 | G4AngularDistribution::G4AngularDistribution(G4bool symmetrize) |
---|
| 33 | : sym(symmetrize) |
---|
| 34 | { |
---|
| 35 | // The following are parameters of the model - not to be confused with the PDG values! |
---|
| 36 | |
---|
| 37 | mSigma = 0.55; |
---|
| 38 | cmSigma = 1.20; |
---|
| 39 | gSigma = 9.4; |
---|
| 40 | |
---|
| 41 | mOmega = 0.783; |
---|
| 42 | cmOmega = 0.808; |
---|
| 43 | gOmega = 10.95; |
---|
| 44 | |
---|
| 45 | mPion = 0.138; |
---|
| 46 | cmPion = 0.51; |
---|
| 47 | gPion = 7.27; |
---|
| 48 | |
---|
| 49 | mNucleon = 0.938; |
---|
| 50 | |
---|
| 51 | // Definition of constants for pion-Term (no s-dependence) |
---|
| 52 | |
---|
| 53 | m42 = 4. * mNucleon * mNucleon; |
---|
| 54 | mPion2 = mPion * mPion; |
---|
| 55 | cmPion2 = cmPion * cmPion; |
---|
| 56 | dPion1 = cmPion2-mPion2; |
---|
| 57 | dPion2 = dPion1 * dPion1; |
---|
| 58 | cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2; |
---|
| 59 | |
---|
| 60 | cPion_3 = -(cm6gp/3.); |
---|
| 61 | cPion_2 = -(cm6gp * mPion2/dPion1); |
---|
| 62 | cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2); |
---|
| 63 | cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2); |
---|
| 64 | cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1); |
---|
| 65 | cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m); |
---|
| 66 | |
---|
| 67 | // Definition of constants for sigma-Term (no s-dependence) |
---|
| 68 | |
---|
| 69 | G4double gSigmaSq = gSigma * gSigma; |
---|
| 70 | |
---|
| 71 | mSigma2 = mSigma * mSigma; |
---|
| 72 | cmSigma2 = cmSigma * cmSigma; |
---|
| 73 | cmSigma4 = cmSigma2 * cmSigma2; |
---|
| 74 | cmSigma6 = cmSigma2 * cmSigma4; |
---|
| 75 | dSigma1 = m42 - cmSigma2; |
---|
| 76 | dSigma2 = m42 - mSigma2; |
---|
| 77 | dSigma3 = cmSigma2 - mSigma2; |
---|
| 78 | |
---|
| 79 | G4double dSigma1Sq = dSigma1 * dSigma1; |
---|
| 80 | G4double dSigma2Sq = dSigma2 * dSigma2; |
---|
| 81 | G4double dSigma3Sq = dSigma3 * dSigma3; |
---|
| 82 | |
---|
| 83 | cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq; |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | cSigma_3 = -(cm2gs * dSigma1Sq / 3.); |
---|
| 87 | cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3); |
---|
| 88 | cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq); |
---|
| 89 | cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq); |
---|
| 90 | cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq)); |
---|
| 91 | cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m); |
---|
| 92 | |
---|
| 93 | // Definition of constants for omega-Term |
---|
| 94 | |
---|
| 95 | G4double gOmegaSq = gOmega * gOmega; |
---|
| 96 | |
---|
| 97 | mOmega2 = mOmega * mOmega; |
---|
| 98 | cmOmega2 = cmOmega * cmOmega; |
---|
| 99 | cmOmega4 = cmOmega2 * cmOmega2; |
---|
| 100 | cmOmega6 = cmOmega2 * cmOmega4; |
---|
| 101 | dOmega1 = m42 - cmOmega2; |
---|
| 102 | dOmega2 = m42 - mOmega2; |
---|
| 103 | dOmega3 = cmOmega2 - mOmega2; |
---|
| 104 | sOmega1 = cmOmega2 + mOmega2; |
---|
| 105 | |
---|
| 106 | G4double dOmega3Sq = dOmega3 * dOmega3; |
---|
| 107 | |
---|
| 108 | cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq; |
---|
| 109 | |
---|
| 110 | cOmega_3 = cm2go / 3.; |
---|
| 111 | cOmega_2 = -(cm2go * cmOmega2 / dOmega3); |
---|
| 112 | cOmega_1 = cm2go * cmOmega4 / dOmega3Sq; |
---|
| 113 | cOmega_m = cm2go * cmOmega6 / (dOmega3Sq * mOmega2); |
---|
| 114 | cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq)); |
---|
| 115 | |
---|
| 116 | // Definition of constants for mix-Term |
---|
| 117 | |
---|
| 118 | G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2); |
---|
| 119 | fac1 = -(fac1Tmp * fac1Tmp * m42); |
---|
| 120 | dMix1 = cmOmega2 - cmSigma2; |
---|
| 121 | dMix2 = cmOmega2 - mSigma2; |
---|
| 122 | dMix3 = cmSigma2 - mOmega2; |
---|
| 123 | |
---|
| 124 | G4double dMix1Sq = dMix1 * dMix1; |
---|
| 125 | G4double dMix2Sq = dMix2 * dMix2; |
---|
| 126 | G4double dMix3Sq = dMix3 * dMix3; |
---|
| 127 | |
---|
| 128 | cMix_o1 = fac1 / (cmOmega2 * dMix1Sq * dMix2 * dOmega3); |
---|
| 129 | cMix_s1 = fac1 / (cmSigma2 * dMix1Sq * dMix3 * dSigma3); |
---|
| 130 | cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2)); |
---|
| 131 | cMix_sm = fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2)); |
---|
| 132 | fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq); |
---|
| 133 | fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq); |
---|
| 134 | |
---|
| 135 | cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4 - cmOmega4 * cmSigma2 |
---|
| 136 | - 2. * cmOmega4 * mOmega2 - 2. * cmOmega4 * mSigma2 |
---|
| 137 | + cmOmega2 * mOmega2 * mSigma2 + cmSigma2 * mOmega2 * mSigma2 |
---|
| 138 | - 4. * cmOmega4 * m42 + 2. * cmOmega2 * cmSigma2 * m42 |
---|
| 139 | + 3. * cmOmega2 * mOmega2 * m42 - cmSigma2 * mOmega2 * m42 |
---|
| 140 | + 3. * cmOmega2 * mSigma2 * m42 - cmSigma2 * mSigma2 * m42 |
---|
| 141 | - 2. * mOmega2 * mSigma2 * m42); |
---|
| 142 | |
---|
| 143 | cMix_oLs = fac2 * (8. * cmOmega4 - 4. * cmOmega2 * cmSigma2 |
---|
| 144 | - 6. * cmOmega2 * mOmega2 + 2. * cmSigma2 * mOmega2 |
---|
| 145 | - 6. * cmOmega2 * mSigma2 + 2. * cmSigma2 * mSigma2 |
---|
| 146 | + 4. * mOmega2 * mSigma2); |
---|
| 147 | |
---|
| 148 | cMix_sLc = fac3 * (cmOmega2 * cmSigma4 - 3. * cmSigma6 |
---|
| 149 | + 2. * cmSigma4 * mOmega2 + 2. * cmSigma4 * mSigma2 |
---|
| 150 | - cmOmega2 * mOmega2 * mSigma2 - cmSigma2 * mOmega2 * mSigma2 |
---|
| 151 | - 2. * cmOmega2 * cmSigma2 * m42 + 4. * cmSigma4 * m42 |
---|
| 152 | + cmOmega2 * mOmega2 * m42 - 3. * cmSigma2 * mOmega2 * m42 |
---|
| 153 | + cmOmega2 * mSigma2 * m42 - 3. * cmSigma2 * mSigma2 * m42 |
---|
| 154 | + 2. * mOmega2 * mSigma2 * m42); |
---|
| 155 | |
---|
| 156 | cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2 - 8. * cmSigma4 |
---|
| 157 | - 2. * cmOmega2 * mOmega2 + 6. * cmSigma2 * mOmega2 |
---|
| 158 | - 2. * cmOmega2 * mSigma2 + 6. * cmSigma2 * mSigma2 |
---|
| 159 | - 4. * mOmega2 * mSigma2); |
---|
| 160 | } |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | G4AngularDistribution::~ |
---|
| 164 | G4AngularDistribution() |
---|
| 165 | { } |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | G4double G4AngularDistribution::CosTheta(G4double s, G4double m1, G4double m2) const |
---|
| 169 | { |
---|
| 170 | G4double random = G4UniformRand(); |
---|
| 171 | G4double dCosTheta = 2.; |
---|
| 172 | G4double cosTheta = -1.; |
---|
| 173 | |
---|
| 174 | // For jmax=12 the accuracy is better than 0.1 degree |
---|
| 175 | G4int jMax = 12; |
---|
| 176 | |
---|
| 177 | for (G4int j = 1; j <= jMax; ++j) |
---|
| 178 | { |
---|
| 179 | // Accuracy is 2^-jmax |
---|
| 180 | dCosTheta *= 0.5; |
---|
| 181 | G4double cosTh = cosTheta + dCosTheta; |
---|
| 182 | if(DifferentialCrossSection(s, m1, m2, cosTh) <= random) cosTheta = cosTh; |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | // Randomize in final interval in order to avoid discrete angles |
---|
| 186 | cosTheta += G4UniformRand() * dCosTheta; |
---|
| 187 | |
---|
| 188 | |
---|
| 189 | if (cosTheta > 1. || cosTheta < -1.) |
---|
| 190 | throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range"); |
---|
| 191 | |
---|
| 192 | return cosTheta; |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m1, G4double m2, |
---|
| 197 | G4double cosTheta) const |
---|
| 198 | { |
---|
| 199 | // local calculus is in GeV, ie. normalize input |
---|
| 200 | sIn = sIn/sqr(GeV)+m42/2.; |
---|
| 201 | m1 = m1/GeV; |
---|
| 202 | m2 = m2/GeV; |
---|
| 203 | // G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl; |
---|
| 204 | // scaling from masses other than p,p. |
---|
| 205 | G4double s = sIn - (m1+m2) * (m1+m2) + m42; |
---|
| 206 | G4double tMax = s - m42; |
---|
| 207 | G4double tp = 0.5 * (cosTheta + 1.) * tMax; |
---|
| 208 | G4double twoS = 2. * s; |
---|
| 209 | |
---|
| 210 | // Define s-dependent stuff for omega-Term |
---|
| 211 | G4double brak1 = (twoS-m42) * (twoS-m42); |
---|
| 212 | G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1); |
---|
| 213 | G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1); |
---|
| 214 | G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2 |
---|
| 215 | - 2. * mOmega2*mOmega2 |
---|
| 216 | - 2. * (cmOmega2 + 2 * mOmega2) * twoS |
---|
| 217 | - 3. * brak1); |
---|
| 218 | G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1); |
---|
| 219 | G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * s + brak1); |
---|
| 220 | G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m); |
---|
| 221 | |
---|
| 222 | // Define s-dependent stuff for mix-Term |
---|
| 223 | G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS); |
---|
| 224 | G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS); |
---|
| 225 | G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS); |
---|
| 226 | G4double bMix_sm = cMix_sm * (dSigma2 - twoS); |
---|
| 227 | G4double bMix_oL = cMix_oLc + cMix_oLs * s; |
---|
| 228 | G4double bMix_sL = cMix_sLc + cMix_sLs * s; |
---|
| 229 | |
---|
| 230 | G4double t1_Pion = 1. / (1. + tMax / cmPion2); |
---|
| 231 | G4double t2_Pion = 1. + tMax / mPion2; |
---|
| 232 | G4double t1_Sigma = 1. / (1. + tMax / cmSigma2); |
---|
| 233 | G4double t2_Sigma = 1. + tMax / mSigma2; |
---|
| 234 | G4double t1_Omega = 1. / (1. + tMax / cmOmega2); |
---|
| 235 | G4double t2_Omega = 1. + tMax / mOmega2; |
---|
| 236 | |
---|
| 237 | G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
| 238 | t2_Pion, t2_Sigma, t2_Omega, |
---|
| 239 | bMix_o1, bMix_s1, bMix_Omega, |
---|
| 240 | bMix_sm, bMix_oL, bMix_sL, |
---|
| 241 | bOmega_0, bOmega_1, bOmega_2, |
---|
| 242 | bOmega_3, bOmega_m, bOmega_L); |
---|
| 243 | |
---|
| 244 | t1_Pion = 1. / (1. + tp / cmPion2); |
---|
| 245 | t2_Pion = 1. + tp / mPion2; |
---|
| 246 | t1_Sigma = 1. / (1. + tp / cmSigma2); |
---|
| 247 | t2_Sigma = 1. + tp / mSigma2; |
---|
| 248 | t1_Omega = 1. / (1. + tp / cmOmega2); |
---|
| 249 | t2_Omega = 1. + tp / mOmega2; |
---|
| 250 | |
---|
| 251 | G4double dSigma; |
---|
| 252 | if (sym) |
---|
| 253 | { |
---|
| 254 | G4double to; |
---|
| 255 | norm = 2. * norm; |
---|
| 256 | to = tMax - tp; |
---|
| 257 | G4double t3_Pion = 1. / (1. + to / cmPion2); |
---|
| 258 | G4double t4_Pion = 1. + to / mPion2; |
---|
| 259 | G4double t3_Sigma = 1. / (1. + to / cmSigma2); |
---|
| 260 | G4double t4_Sigma = 1. + to / mSigma2; |
---|
| 261 | G4double t3_Omega = 1. / (1. + to / cmOmega2); |
---|
| 262 | G4double t4_Omega = 1. + to / mOmega2; |
---|
| 263 | |
---|
| 264 | dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
| 265 | t2_Pion,t2_Sigma, t2_Omega, |
---|
| 266 | bMix_o1, bMix_s1, bMix_Omega, |
---|
| 267 | bMix_sm, bMix_oL, bMix_sL, |
---|
| 268 | bOmega_0, bOmega_1, bOmega_2, |
---|
| 269 | bOmega_3, bOmega_m, bOmega_L) - |
---|
| 270 | Cross(t3_Pion,t3_Sigma, t3_Omega, |
---|
| 271 | t4_Pion, t4_Sigma, t4_Omega, |
---|
| 272 | bMix_o1, bMix_s1, bMix_Omega, |
---|
| 273 | bMix_sm, bMix_oL, bMix_sL, |
---|
| 274 | bOmega_0, bOmega_1, bOmega_2, |
---|
| 275 | bOmega_3, bOmega_m, bOmega_L) ) |
---|
| 276 | / norm + 0.5; |
---|
| 277 | } |
---|
| 278 | else |
---|
| 279 | { |
---|
| 280 | dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
| 281 | t2_Pion, t2_Sigma, t2_Omega, |
---|
| 282 | bMix_o1, bMix_s1, bMix_Omega, |
---|
| 283 | bMix_sm, bMix_oL, bMix_sL, |
---|
| 284 | bOmega_0, bOmega_1, bOmega_2, |
---|
| 285 | bOmega_3, bOmega_m, bOmega_L) |
---|
| 286 | / norm; |
---|
| 287 | } |
---|
| 288 | |
---|
| 289 | return dSigma; |
---|
| 290 | } |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | G4double G4AngularDistribution::Cross(G4double tpPion, |
---|
| 294 | G4double tpSigma, |
---|
| 295 | G4double tpOmega, |
---|
| 296 | G4double tmPion, |
---|
| 297 | G4double tmSigma, |
---|
| 298 | G4double tmOmega, |
---|
| 299 | G4double bMix_o1, |
---|
| 300 | G4double bMix_s1, |
---|
| 301 | G4double bMix_Omega, |
---|
| 302 | G4double bMix_sm, |
---|
| 303 | G4double bMix_oL, |
---|
| 304 | G4double bMix_sL, |
---|
| 305 | G4double bOmega_0, |
---|
| 306 | G4double bOmega_1, |
---|
| 307 | G4double bOmega_2, |
---|
| 308 | G4double bOmega_3, |
---|
| 309 | G4double bOmega_m, |
---|
| 310 | G4double bOmega_L) const |
---|
| 311 | { |
---|
| 312 | G4double cross = 0; |
---|
| 313 | // Pion |
---|
| 314 | cross += ((cPion_3 * tpPion + cPion_2) * tpPion + cPion_1) * tpPion + cPion_m/tmPion + cPion_0 + cPion_L * std::log(tpPion*tmPion); |
---|
| 315 | // G4cout << "cross1 "<< cross<<G4endl; |
---|
| 316 | // Sigma |
---|
| 317 | cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * std::log(tpSigma*tmSigma); |
---|
| 318 | // G4cout << "cross2 "<< cross<<G4endl; |
---|
| 319 | // Omega |
---|
| 320 | cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * std::log(tpOmega*tmOmega) |
---|
| 321 | // Mix |
---|
| 322 | + bMix_o1 * (tpOmega - 1.) |
---|
| 323 | + bMix_s1 * (tpSigma - 1.) |
---|
| 324 | + bMix_Omega * std::log(tmOmega) |
---|
| 325 | + bMix_sm * std::log(tmSigma) |
---|
| 326 | + bMix_oL * std::log(tpOmega) |
---|
| 327 | + bMix_sL * std::log(tpSigma); |
---|
| 328 | /* G4cout << "cross3 "<< cross<<" " |
---|
| 329 | <<bMix_o1<<" " |
---|
| 330 | <<bMix_s1<<" " |
---|
| 331 | <<bMix_Omega<<" " |
---|
| 332 | <<bMix_sm<<" " |
---|
| 333 | <<bMix_oL<<" " |
---|
| 334 | <<bMix_sL<<" " |
---|
| 335 | <<tpOmega<<" " |
---|
| 336 | <<tpSigma<<" " |
---|
| 337 | <<tmOmega<<" " |
---|
| 338 | <<tmSigma<<" " |
---|
| 339 | <<tpOmega<<" " |
---|
| 340 | <<tpSigma |
---|
| 341 | <<G4endl; |
---|
| 342 | */ |
---|
| 343 | return cross; |
---|
| 344 | |
---|
| 345 | } |
---|