| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // hpw: done, but low quality at present. 26 // hpw: done, but low quality at present.
27 27
28 #include "globals.hh" 28 #include "globals.hh"
29 #include "G4Log.hh" <<
30 #include "G4SystemOfUnits.hh" 29 #include "G4SystemOfUnits.hh"
31 #include "G4AngularDistribution.hh" 30 #include "G4AngularDistribution.hh"
32 #include "Randomize.hh" 31 #include "Randomize.hh"
33 32
34 G4AngularDistribution::G4AngularDistribution(G 33 G4AngularDistribution::G4AngularDistribution(G4bool symmetrize)
35 : sym(symmetrize) 34 : sym(symmetrize)
36 { 35 {
37 // The following are parameters of the model 36 // The following are parameters of the model - not to be confused with the PDG values!
38 37
39 mSigma = 0.55; 38 mSigma = 0.55;
40 cmSigma = 1.20; 39 cmSigma = 1.20;
41 gSigma = 9.4; 40 gSigma = 9.4;
42 41
43 mOmega = 0.783; 42 mOmega = 0.783;
44 cmOmega = 0.808; 43 cmOmega = 0.808;
45 gOmega = 10.95; 44 gOmega = 10.95;
46 45
47 mPion = 0.138; 46 mPion = 0.138;
48 cmPion = 0.51; 47 cmPion = 0.51;
49 gPion = 7.27; 48 gPion = 7.27;
50 49
51 mNucleon = 0.938; 50 mNucleon = 0.938;
52 51
53 // Definition of constants for pion-Term 52 // Definition of constants for pion-Term (no s-dependence)
54 53
55 m42 = 4. * mNucleon * mNucleon; 54 m42 = 4. * mNucleon * mNucleon;
56 mPion2 = mPion * mPion; 55 mPion2 = mPion * mPion;
57 cmPion2 = cmPion * cmPion; 56 cmPion2 = cmPion * cmPion;
58 dPion1 = cmPion2-mPion2; 57 dPion1 = cmPion2-mPion2;
59 dPion2 = dPion1 * dPion1; 58 dPion2 = dPion1 * dPion1;
60 cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) 59 cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2;
61 60
62 cPion_3 = -(cm6gp/3.); 61 cPion_3 = -(cm6gp/3.);
63 cPion_2 = -(cm6gp * mPion2/dPion1); 62 cPion_2 = -(cm6gp * mPion2/dPion1);
64 cPion_1 = -(cm6gp * mPion2 * (2. * cmPio 63 cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2);
65 cPion_m = -(cm6gp * cmPion2 * mPion2 / d 64 cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2);
66 cPion_L = -(cm6gp * 2. * cmPion2 * mPion 65 cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1);
67 cPion_0 = -(cPion_3 + cPion_2 + cPion_1 66 cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m);
68 67
69 // Definition of constants for sigma-Term 68 // Definition of constants for sigma-Term (no s-dependence)
70 69
71 G4double gSigmaSq = gSigma * gSigma; 70 G4double gSigmaSq = gSigma * gSigma;
72 71
73 mSigma2 = mSigma * mSigma; 72 mSigma2 = mSigma * mSigma;
74 cmSigma2 = cmSigma * cmSigma; 73 cmSigma2 = cmSigma * cmSigma;
75 cmSigma4 = cmSigma2 * cmSigma2; 74 cmSigma4 = cmSigma2 * cmSigma2;
76 cmSigma6 = cmSigma2 * cmSigma4; 75 cmSigma6 = cmSigma2 * cmSigma4;
77 dSigma1 = m42 - cmSigma2; 76 dSigma1 = m42 - cmSigma2;
78 dSigma2 = m42 - mSigma2; 77 dSigma2 = m42 - mSigma2;
79 dSigma3 = cmSigma2 - mSigma2; 78 dSigma3 = cmSigma2 - mSigma2;
80 79
81 G4double dSigma1Sq = dSigma1 * dSigma1; 80 G4double dSigma1Sq = dSigma1 * dSigma1;
82 G4double dSigma2Sq = dSigma2 * dSigma2; 81 G4double dSigma2Sq = dSigma2 * dSigma2;
83 G4double dSigma3Sq = dSigma3 * dSigma3; 82 G4double dSigma3Sq = dSigma3 * dSigma3;
84 83
85 cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigma 84 cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq;
86 85
87 86
88 cSigma_3 = -(cm2gs * dSigma1Sq / 3.); 87 cSigma_3 = -(cm2gs * dSigma1Sq / 3.);
89 cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 88 cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3);
90 cSigma_1 = -(cm2gs * cmSigma4 * (2. * dS 89 cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq);
91 cSigma_m = -(cm2gs * cmSigma6 * dSigma2S 90 cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq);
92 cSigma_L = -(cm2gs * cmSigma6 * dSigma2 91 cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq));
93 cSigma_0 = -(cSigma_3 + cSigma_2 + cSigm 92 cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m);
94 93
95 // Definition of constants for omega-Ter 94 // Definition of constants for omega-Term
96 95
97 G4double gOmegaSq = gOmega * gOmega; 96 G4double gOmegaSq = gOmega * gOmega;
98 97
99 mOmega2 = mOmega * mOmega; 98 mOmega2 = mOmega * mOmega;
100 cmOmega2 = cmOmega * cmOmega; 99 cmOmega2 = cmOmega * cmOmega;
101 cmOmega4 = cmOmega2 * cmOmega2; 100 cmOmega4 = cmOmega2 * cmOmega2;
102 cmOmega6 = cmOmega2 * cmOmega4; 101 cmOmega6 = cmOmega2 * cmOmega4;
103 dOmega1 = m42 - cmOmega2; 102 dOmega1 = m42 - cmOmega2;
104 dOmega2 = m42 - mOmega2; 103 dOmega2 = m42 - mOmega2;
105 dOmega3 = cmOmega2 - mOmega2; 104 dOmega3 = cmOmega2 - mOmega2;
106 sOmega1 = cmOmega2 + mOmega2; 105 sOmega1 = cmOmega2 + mOmega2;
107 106
108 G4double dOmega3Sq = dOmega3 * dOmega3; 107 G4double dOmega3Sq = dOmega3 * dOmega3;
109 108
110 cm2go = 0.5 * cmOmega2 * gOmegaSq * gOme 109 cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq;
111 110
112 cOmega_3 = cm2go / 3.; 111 cOmega_3 = cm2go / 3.;
113 cOmega_2 = -(cm2go * cmOmega2 / dOmega3) 112 cOmega_2 = -(cm2go * cmOmega2 / dOmega3);
114 cOmega_1 = cm2go * cmOmega4 / dOmega3Sq 113 cOmega_1 = cm2go * cmOmega4 / dOmega3Sq;
115 cOmega_m = cm2go * cmOmega6 / (dOmega3S 114 cOmega_m = cm2go * cmOmega6 / (dOmega3Sq * mOmega2);
116 cOmega_L = -(cm2go * cmOmega6 * 4. / (dO 115 cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq));
117 116
118 // Definition of constants for mix-Term 117 // Definition of constants for mix-Term
119 118
120 G4double fac1Tmp = (gSigma * gOmega * cm 119 G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2);
121 fac1 = -(fac1Tmp * fac1Tmp * m42); 120 fac1 = -(fac1Tmp * fac1Tmp * m42);
122 dMix1 = cmOmega2 - cmSigma2; 121 dMix1 = cmOmega2 - cmSigma2;
123 dMix2 = cmOmega2 - mSigma2; 122 dMix2 = cmOmega2 - mSigma2;
124 dMix3 = cmSigma2 - mOmega2; 123 dMix3 = cmSigma2 - mOmega2;
125 124
126 G4double dMix1Sq = dMix1 * dMix1; 125 G4double dMix1Sq = dMix1 * dMix1;
127 G4double dMix2Sq = dMix2 * dMix2; 126 G4double dMix2Sq = dMix2 * dMix2;
128 G4double dMix3Sq = dMix3 * dMix3; 127 G4double dMix3Sq = dMix3 * dMix3;
129 128
130 cMix_o1 = fac1 / (cmOmega2 * dMix1Sq 129 cMix_o1 = fac1 / (cmOmega2 * dMix1Sq * dMix2 * dOmega3);
131 cMix_s1 = fac1 / (cmSigma2 * dMix1Sq 130 cMix_s1 = fac1 / (cmSigma2 * dMix1Sq * dMix3 * dSigma3);
132 cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq 131 cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2));
133 cMix_sm = fac1 / (dSigma3Sq * dMix2Sq 132 cMix_sm = fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2));
134 fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega 133 fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq);
135 fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma 134 fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq);
136 135
137 cMix_oLc = fac2 * (3. * cmOmega2*cmOmega 136 cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4 - cmOmega4 * cmSigma2
138 - 2. * cmOmega4 * mOmega2 137 - 2. * cmOmega4 * mOmega2 - 2. * cmOmega4 * mSigma2
139 + cmOmega2 * mOmega2 * mSigma2 138 + cmOmega2 * mOmega2 * mSigma2 + cmSigma2 * mOmega2 * mSigma2
140 - 4. * cmOmega4 * m42 139 - 4. * cmOmega4 * m42 + 2. * cmOmega2 * cmSigma2 * m42
141 + 3. * cmOmega2 * mOmega2 * m42 140 + 3. * cmOmega2 * mOmega2 * m42 - cmSigma2 * mOmega2 * m42
142 + 3. * cmOmega2 * mSigma2 * m42 141 + 3. * cmOmega2 * mSigma2 * m42 - cmSigma2 * mSigma2 * m42
143 - 2. * mOmega2 * mSigma2 * m42); 142 - 2. * mOmega2 * mSigma2 * m42);
144 143
145 cMix_oLs = fac2 * (8. * cmOmega4 144 cMix_oLs = fac2 * (8. * cmOmega4 - 4. * cmOmega2 * cmSigma2
146 - 6. * cmOmega2 * mOmega2 + 2. 145 - 6. * cmOmega2 * mOmega2 + 2. * cmSigma2 * mOmega2
147 - 6. * cmOmega2 * mSigma2 + 2. 146 - 6. * cmOmega2 * mSigma2 + 2. * cmSigma2 * mSigma2
148 + 4. * mOmega2 * mSigma2); 147 + 4. * mOmega2 * mSigma2);
149 148
150 cMix_sLc = fac3 * (cmOmega2 * cmSigma4 149 cMix_sLc = fac3 * (cmOmega2 * cmSigma4 - 3. * cmSigma6
151 + 2. * cmSigma4 * mOmega2 + 2. 150 + 2. * cmSigma4 * mOmega2 + 2. * cmSigma4 * mSigma2
152 - cmOmega2 * mOmega2 * mSigma2 - cmS 151 - cmOmega2 * mOmega2 * mSigma2 - cmSigma2 * mOmega2 * mSigma2
153 - 2. * cmOmega2 * cmSigma2 * m42 + 4. 152 - 2. * cmOmega2 * cmSigma2 * m42 + 4. * cmSigma4 * m42
154 + cmOmega2 * mOmega2 * m42 - 3. 153 + cmOmega2 * mOmega2 * m42 - 3. * cmSigma2 * mOmega2 * m42
155 + cmOmega2 * mSigma2 * m42 - 3. 154 + cmOmega2 * mSigma2 * m42 - 3. * cmSigma2 * mSigma2 * m42
156 + 2. * mOmega2 * mSigma2 * m42); 155 + 2. * mOmega2 * mSigma2 * m42);
157 156
158 cMix_sLs = fac3 * (4. * cmOmega2 * cmSig 157 cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2 - 8. * cmSigma4
159 - 2. * cmOmega2 * mOmega2 158 - 2. * cmOmega2 * mOmega2 + 6. * cmSigma2 * mOmega2
160 - 2. * cmOmega2 * mSigma2 159 - 2. * cmOmega2 * mSigma2 + 6. * cmSigma2 * mSigma2
161 - 4. * mOmega2 * mSigma2); 160 - 4. * mOmega2 * mSigma2);
162 } 161 }
163 162
164 163
165 G4AngularDistribution::~ 164 G4AngularDistribution::~
166 G4AngularDistribution() 165 G4AngularDistribution()
167 { } 166 { }
168 167
169 168
170 G4double G4AngularDistribution::CosTheta(G4dou 169 G4double G4AngularDistribution::CosTheta(G4double S, G4double m_1, G4double m_2) const
171 { 170 {
172 G4double random = G4UniformRand(); 171 G4double random = G4UniformRand();
173 G4double dCosTheta = 2.; 172 G4double dCosTheta = 2.;
174 G4double cosTheta = -1.; 173 G4double cosTheta = -1.;
175 174
176 // For jmax=12 the accuracy is better than 175 // For jmax=12 the accuracy is better than 0.1 degree
177 G4int jMax = 12; 176 G4int jMax = 12;
178 177
179 for (G4int j = 1; j <= jMax; ++j) 178 for (G4int j = 1; j <= jMax; ++j)
180 { 179 {
181 // Accuracy is 2^-jmax 180 // Accuracy is 2^-jmax
182 dCosTheta *= 0.5; 181 dCosTheta *= 0.5;
183 G4double cosTh = cosTheta + dCosTheta; 182 G4double cosTh = cosTheta + dCosTheta;
184 if(DifferentialCrossSection(S, m_1, m_2, 183 if(DifferentialCrossSection(S, m_1, m_2, cosTh) <= random) cosTheta = cosTh;
185 } 184 }
186 185
187 // Randomize in final interval in order to 186 // Randomize in final interval in order to avoid discrete angles
188 cosTheta += G4UniformRand() * dCosTheta; 187 cosTheta += G4UniformRand() * dCosTheta;
189 188
190 189
191 if (cosTheta > 1. || cosTheta < -1.) 190 if (cosTheta > 1. || cosTheta < -1.)
192 throw G4HadronicException(__FILE__, __LIN 191 throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range");
193 192
194 return cosTheta; 193 return cosTheta;
195 } 194 }
196 195
197 196
198 G4double G4AngularDistribution::DifferentialCr 197 G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m_1, G4double m_2,
199 G4double cosTheta) const 198 G4double cosTheta) const
200 { 199 {
201 // local calculus is in GeV, ie. normalize inp 200 // local calculus is in GeV, ie. normalize input
202 sIn = sIn/sqr(GeV)+m42/2.; 201 sIn = sIn/sqr(GeV)+m42/2.;
203 m_1 = m_1/GeV; 202 m_1 = m_1/GeV;
204 m_2 = m_2/GeV; 203 m_2 = m_2/GeV;
205 // G4cout << "Here we go"<<sIn << " "<<m1 << 204 // G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl;
206 // scaling from masses other than p,p. 205 // scaling from masses other than p,p.
207 G4double S = sIn - (m_1+m_2) * (m_1+m_2) + m 206 G4double S = sIn - (m_1+m_2) * (m_1+m_2) + m42;
208 G4double tMax = S - m42; 207 G4double tMax = S - m42;
209 G4double tp = 0.5 * (cosTheta + 1.) * tMax; 208 G4double tp = 0.5 * (cosTheta + 1.) * tMax;
210 G4double twoS = 2. * S; 209 G4double twoS = 2. * S;
211 210
212 // Define s-dependent stuff for omega-Term 211 // Define s-dependent stuff for omega-Term
213 G4double brak1 = (twoS-m42) * (twoS-m42); 212 G4double brak1 = (twoS-m42) * (twoS-m42);
214 G4double bOmega_3 = cOmega_3 * (-2. * cmOmeg 213 G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1);
215 G4double bOmega_2 = cOmega_2 * ( 2. * cmOmeg 214 G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1);
216 G4double bOmega_1 = cOmega_1 * (-4. * cmOmeg 215 G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2
217 - 2. * mOmega2*mOmega2 216 - 2. * mOmega2*mOmega2
218 - 2. * (cmOmega2 + 2 * mOmega2) * twoS 217 - 2. * (cmOmega2 + 2 * mOmega2) * twoS
219 - 3. * brak1); 218 - 3. * brak1);
220 G4double bOmega_m = cOmega_m * (-2. * mOmega 219 G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1);
221 G4double bOmega_L = cOmega_L * (sOmega1 * mO 220 G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * S + brak1);
222 G4double bOmega_0 = -(bOmega_3 + bOmega_2 + 221 G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m);
223 222
224 // Define s-dependent stuff for mix-Term 223 // Define s-dependent stuff for mix-Term
225 G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS 224 G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS);
226 G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS 225 G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS);
227 G4double bMix_Omega = cMix_Omega * (dOmega2 226 G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS);
228 G4double bMix_sm = cMix_sm * (dSigma2 - twoS 227 G4double bMix_sm = cMix_sm * (dSigma2 - twoS);
229 G4double bMix_oL = cMix_oLc + cMix_oLs * S; 228 G4double bMix_oL = cMix_oLc + cMix_oLs * S;
230 G4double bMix_sL = cMix_sLc + cMix_sLs * S; 229 G4double bMix_sL = cMix_sLc + cMix_sLs * S;
231 230
232 G4double t1_Pion = 1. / (1. + tMax / cmPion2 231 G4double t1_Pion = 1. / (1. + tMax / cmPion2);
233 G4double t2_Pion = 1. + tMax / mPion2; 232 G4double t2_Pion = 1. + tMax / mPion2;
234 G4double t1_Sigma = 1. / (1. + tMax / cmSigm 233 G4double t1_Sigma = 1. / (1. + tMax / cmSigma2);
235 G4double t2_Sigma = 1. + tMax / mSigma2; 234 G4double t2_Sigma = 1. + tMax / mSigma2;
236 G4double t1_Omega = 1. / (1. + tMax / cmOmeg 235 G4double t1_Omega = 1. / (1. + tMax / cmOmega2);
237 G4double t2_Omega = 1. + tMax / mOmega2; 236 G4double t2_Omega = 1. + tMax / mOmega2;
238 237
239 G4double norm = Cross(t1_Pion, t1_Sigma, t1_ 238 G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega,
240 t2_Pion, t2_Sigma, t2_Omega, 239 t2_Pion, t2_Sigma, t2_Omega,
241 bMix_o1, bMix_s1, bMix_Omega, 240 bMix_o1, bMix_s1, bMix_Omega,
242 bMix_sm, bMix_oL, bMix_sL, 241 bMix_sm, bMix_oL, bMix_sL,
243 bOmega_0, bOmega_1, bOmega_2, 242 bOmega_0, bOmega_1, bOmega_2,
244 bOmega_3, bOmega_m, bOmega_L); 243 bOmega_3, bOmega_m, bOmega_L);
245 244
246 t1_Pion = 1. / (1. + tp / cmPion2); 245 t1_Pion = 1. / (1. + tp / cmPion2);
247 t2_Pion = 1. + tp / mPion2; 246 t2_Pion = 1. + tp / mPion2;
248 t1_Sigma = 1. / (1. + tp / cmSigma2); 247 t1_Sigma = 1. / (1. + tp / cmSigma2);
249 t2_Sigma = 1. + tp / mSigma2; 248 t2_Sigma = 1. + tp / mSigma2;
250 t1_Omega = 1. / (1. + tp / cmOmega2); 249 t1_Omega = 1. / (1. + tp / cmOmega2);
251 t2_Omega = 1. + tp / mOmega2; 250 t2_Omega = 1. + tp / mOmega2;
252 251
253 G4double dSigma; 252 G4double dSigma;
254 if (sym) 253 if (sym)
255 { 254 {
256 G4double to; 255 G4double to;
257 norm = 2. * norm; 256 norm = 2. * norm;
258 to = tMax - tp; 257 to = tMax - tp;
259 G4double t3_Pion = 1. / (1. + to / cmPio 258 G4double t3_Pion = 1. / (1. + to / cmPion2);
260 G4double t4_Pion = 1. + to / mPion2; 259 G4double t4_Pion = 1. + to / mPion2;
261 G4double t3_Sigma = 1. / (1. + to / cmSi 260 G4double t3_Sigma = 1. / (1. + to / cmSigma2);
262 G4double t4_Sigma = 1. + to / mSigma2; 261 G4double t4_Sigma = 1. + to / mSigma2;
263 G4double t3_Omega = 1. / (1. + to / cmOm 262 G4double t3_Omega = 1. / (1. + to / cmOmega2);
264 G4double t4_Omega = 1. + to / mOmega2; 263 G4double t4_Omega = 1. + to / mOmega2;
265 264
266 dSigma = ( Cross(t1_Pion, t1_Sigma, t1_O 265 dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega,
267 t2_Pion,t2_Sigma, t2_Omega, 266 t2_Pion,t2_Sigma, t2_Omega,
268 bMix_o1, bMix_s1, bMix_Omega, 267 bMix_o1, bMix_s1, bMix_Omega,
269 bMix_sm, bMix_oL, bMix_sL, 268 bMix_sm, bMix_oL, bMix_sL,
270 bOmega_0, bOmega_1, bOmega_2, 269 bOmega_0, bOmega_1, bOmega_2,
271 bOmega_3, bOmega_m, bOmega_L) - 270 bOmega_3, bOmega_m, bOmega_L) -
272 Cross(t3_Pion,t3_Sigma, t3_Omega, 271 Cross(t3_Pion,t3_Sigma, t3_Omega,
273 t4_Pion, t4_Sigma, t4_Omega, 272 t4_Pion, t4_Sigma, t4_Omega,
274 bMix_o1, bMix_s1, bMix_Omega, 273 bMix_o1, bMix_s1, bMix_Omega,
275 bMix_sm, bMix_oL, bMix_sL, 274 bMix_sm, bMix_oL, bMix_sL,
276 bOmega_0, bOmega_1, bOmega_2, 275 bOmega_0, bOmega_1, bOmega_2,
277 bOmega_3, bOmega_m, bOmega_L) ) 276 bOmega_3, bOmega_m, bOmega_L) )
278 / norm + 0.5; 277 / norm + 0.5;
279 } 278 }
280 else 279 else
281 { 280 {
282 dSigma = Cross(t1_Pion, t1_Sigma, t1_Ome 281 dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega,
283 t2_Pion, t2_Sigma, t2_Omega, 282 t2_Pion, t2_Sigma, t2_Omega,
284 bMix_o1, bMix_s1, bMix_Omega, 283 bMix_o1, bMix_s1, bMix_Omega,
285 bMix_sm, bMix_oL, bMix_sL, 284 bMix_sm, bMix_oL, bMix_sL,
286 bOmega_0, bOmega_1, bOmega_2, 285 bOmega_0, bOmega_1, bOmega_2,
287 bOmega_3, bOmega_m, bOmega_L) 286 bOmega_3, bOmega_m, bOmega_L)
288 / norm; 287 / norm;
289 } 288 }
290 289
291 return dSigma; 290 return dSigma;
292 } 291 }
293 292
294 293
295 G4double G4AngularDistribution::Cross(G4double 294 G4double G4AngularDistribution::Cross(G4double tpPion,
296 G4double tpSigma, 295 G4double tpSigma,
297 G4double tpOmega, 296 G4double tpOmega,
298 G4double tmPion, 297 G4double tmPion,
299 G4double tmSigma, 298 G4double tmSigma,
300 G4double tmOmega, 299 G4double tmOmega,
301 G4double bMix_o1, 300 G4double bMix_o1,
302 G4double bMix_s1, 301 G4double bMix_s1,
303 G4double bMix_Omega, 302 G4double bMix_Omega,
304 G4double bMix_sm, 303 G4double bMix_sm,
305 G4double bMix_oL, 304 G4double bMix_oL,
306 G4double bMix_sL, 305 G4double bMix_sL,
307 G4double bOmega_0, 306 G4double bOmega_0,
308 G4double bOmega_1, 307 G4double bOmega_1,
309 G4double bOmega_2, 308 G4double bOmega_2,
310 G4double bOmega_3, 309 G4double bOmega_3,
311 G4double bOmega_m, 310 G4double bOmega_m,
312 G4double bOmega_L) const 311 G4double bOmega_L) const
313 { 312 {
314 G4double cross = 0; 313 G4double cross = 0;
315 // Pion 314 // Pion
316 cross += ((cPion_3 * tpPion + cPion_2) * << 315 cross += ((cPion_3 * tpPion + cPion_2) * tpPion + cPion_1) * tpPion + cPion_m/tmPion + cPion_0 + cPion_L * std::log(tpPion*tmPion);
317 // G4cout << "cross1 "<< cross<<G4endl; 316 // G4cout << "cross1 "<< cross<<G4endl;
318 // Sigma 317 // Sigma
319 cross += ((cSigma_3 * tpSigma + cSigma_2) << 318 cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * std::log(tpSigma*tmSigma);
320 // G4cout << "cross2 "<< cross<<G4endl; 319 // G4cout << "cross2 "<< cross<<G4endl;
321 // Omega 320 // Omega
322 cross += ((bOmega_3 * tpOmega + bOmega_2) << 321 cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * std::log(tpOmega*tmOmega)
323 // Mix 322 // Mix
324 + bMix_o1 * (tpOmega - 1.) 323 + bMix_o1 * (tpOmega - 1.)
325 + bMix_s1 * (tpSigma - 1.) 324 + bMix_s1 * (tpSigma - 1.)
326 + bMix_Omega * G4Log(tmOmega) << 325 + bMix_Omega * std::log(tmOmega)
327 + bMix_sm * G4Log(tmSigma) << 326 + bMix_sm * std::log(tmSigma)
328 + bMix_oL * G4Log(tpOmega) << 327 + bMix_oL * std::log(tpOmega)
329 + bMix_sL * G4Log(tpSigma); << 328 + bMix_sL * std::log(tpSigma);
330 /* G4cout << "cross3 "<< cross<<" " 329 /* G4cout << "cross3 "<< cross<<" "
331 <<bMix_o1<<" " 330 <<bMix_o1<<" "
332 <<bMix_s1<<" " 331 <<bMix_s1<<" "
333 <<bMix_Omega<<" " 332 <<bMix_Omega<<" "
334 <<bMix_sm<<" " 333 <<bMix_sm<<" "
335 <<bMix_oL<<" " 334 <<bMix_oL<<" "
336 <<bMix_sL<<" " 335 <<bMix_sL<<" "
337 <<tpOmega<<" " 336 <<tpOmega<<" "
338 <<tpSigma<<" " 337 <<tpSigma<<" "
339 <<tmOmega<<" " 338 <<tmOmega<<" "
340 <<tmSigma<<" " 339 <<tmSigma<<" "
341 <<tpOmega<<" " 340 <<tpOmega<<" "
342 <<tpSigma 341 <<tpSigma
343 <<G4endl; 342 <<G4endl;
344 */ 343 */
345 return cross; 344 return cross;
346 345
347 } 346 }
348 347