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