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