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