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