27 {
28
29
30
31
32
33 double A0 = std::sqrt(params[0] * (1 - params[2]));
34 double Al = std::sqrt(params[1] * (1 - params[2]));
35 double As = std::sqrt(params[2]);
36 const double &GammaS =
params[3];
37 const double &DeltaGamma =
params[4];
38 const double &DeltaM =
params[5];
39 const double &phiS =
params[6];
40 const double &delta_p =
params[7];
41 const double &delta_l =
params[8];
43
45
46 double &tagprob =
x[4];
47
48
49 double Ap = 0;
50 if ( 1. - As*As - A0*A0 - Al*Al >= 0 ) Ap = std::sqrt(1. - As*As - A0*A0 - Al*Al);
51
52
53
54 double A[10], B1[10], B2[10],
C[10];
55
56
57 double sinphiS = std::sin(phiS);
58 double cosphiS = std::cos(phiS);
59
63 A[3] = 0.5 * A0*Al * std::cos(delta_l);
68 A[8] = 0.5 * As*Ap * std::sin(delta_p - delta_s);
70
71
72
73
74 double gammaH = GammaS - DeltaGamma/2.;
75 double gammaL = GammaS + DeltaGamma/2.;
76 double tauL = 1./gammaL;
77 double tauH = 1./gammaH;
78 double norm1 = (1./(tauL*(1.+cosphiS) + tauH*(1.-cosphiS)));
79 double norm2 = (1./(tauL*(1.-cosphiS) + tauH*(1.+cosphiS)));
80 double norm3 = (1./sqrt((tauL-tauH)*(tauL-tauH)*sinphiS*sinphiS+4.*tauL*tauH));
81
82 double ExpGSSinMT = 0.0;
83 double ExpGSCosMT = 0.0;
84
85
86
87
88 double ExpGL = std::exp(-time * gammaL);
89
90
91 double ExpGH = std::exp(-time * gammaH);
92
93 B1[0] = ( (1. + cosphiS) * ExpGL + (1. - cosphiS) * ExpGH ) * norm1;
94 B1[1] = B1[0];
95 B1[2] = ( (1. - cosphiS) * ExpGL + (1. + cosphiS) * ExpGH ) * norm2;
96 B1[3] = B1[0];
97 B1[4] = 0.5 * std::cos(delta_p - delta_l) * sinphiS * (ExpGL - ExpGH) * norm3;
98 B1[5] = 0.5 * std::cos(delta_p) * sinphiS * (ExpGL - ExpGH) * norm3;
99 B1[6] = B1[2];
100 B1[7] = 0.5 * sinphiS * std::sin(delta_l - delta_s) * (ExpGL - ExpGH) * norm3;
101 B1[8] = B1[2];
102 B1[9] = 0.5 * sinphiS * std::sin(delta_s) * (ExpGH - ExpGL) * norm3;
103
104
105
106
107
108
109 if( std::abs(tagprob - 0.5) > 1e-6 ){
110
111
112 ExpGSSinMT = std::exp(-time * GammaS) * std::sin(DeltaM * time);
113 ExpGSCosMT = std::exp(-time * GammaS) * std::cos(DeltaM * time);
114
115
116 B2[0] = 2. * ExpGSSinMT * sinphiS * norm1;
117 B2[1] = B2[0];
118 B2[2] = -2. * ExpGSSinMT * sinphiS * norm2;
119 B2[3] = B2[0];
120 B2[4] = ( ExpGSCosMT * std::sin(delta_p - delta_l) - cosphiS * std::cos(delta_p - delta_l) * ExpGSSinMT ) * norm3;
121 B2[5] = ( ExpGSCosMT * std::sin(delta_p) - cosphiS * std::cos(delta_p) * ExpGSSinMT ) * norm3;
122 B2[6] = B2[2];
123 B2[7] = ( ExpGSCosMT * std::cos(delta_l - delta_s) - cosphiS * std::sin(delta_l - delta_s) * ExpGSSinMT ) * norm3;
124 B2[8] = B2[2];
125 B2[9] = ( ExpGSCosMT * std::cos(delta_s) + cosphiS * std::sin(delta_s) * ExpGSSinMT ) * norm3;
126
127 }
128 else{
129
130 for (
int k = 0;
k<10;
k++) B2[k] = 0;
131
132 }
133
134
135 if (useHelicity) {
136
137
138 double &costhetal =
x[1];
139 double &costhetak =
x[2];
141
142 double sinsqthetal = 1. - (costhetal * costhetal);
143 double sinsqthetak = 1. - (costhetak * costhetak);
144 double cossqthetak = costhetak * costhetak;
145 double coschi = std::cos(chi);
146 double cossqchi = coschi * coschi;
147 double sinchi = std::sin(chi);
148 double sinsqchi = sinchi * sinchi;
149 double sin2thetak = 2. * std::sqrt(1. - costhetak * costhetak) * costhetak;
150 double sin2thetal = 2. * std::sqrt(1. - costhetal * costhetal) * costhetal;
151 double sin2chi = std::sin(2. * chi);
152 double sinthetak = std::sqrt(sinsqthetak);
153
154 C[0] = 2. * cossqthetak * sinsqthetal;
155 C[1] = sinsqthetak * (1 - sinsqthetal * cossqchi);
156 C[2] = sinsqthetak * (1 - sinsqthetal * sinsqchi);
157 C[3] = 1. / std::sqrt(2.) * sin2thetak * sin2thetal * coschi;
158 C[4] = -1. * sinsqthetak * sinsqthetal * sin2chi;
159 C[5] = 1. / std::sqrt(2.) * sin2thetak * sin2thetal * sinchi;
160 C[6] = 2. / 3. * sinsqthetal;
161 C[7] = 1. / 3. * std::sqrt(6.) * sinthetak * sin2thetal * coschi;
162 C[8] = 1. / 3. * std::sqrt(6.) * sinthetak * sin2thetal * sinchi;
163 C[9] = 4. / 3. * std::sqrt(3.) * costhetak * sinsqthetal;
164
165 }
166 else {
167
168
169 double &costheta =
x[1];
170 double &cospsi =
x[2];
172
173 double sinsqtheta = 1. - (costheta * costheta);
174 double sinsqpsi = 1. - (cospsi * cospsi);
175 double cossqpsi = cospsi * cospsi;
176 double sin2theta = 2. * std::sqrt(1. - costheta * costheta) * costheta;
177 double sin2psi = 2. * std::sqrt(1. - cospsi * cospsi) * cospsi;
178 double cosphi = std::cos(
phi);
179 double cossqphi = cosphi * cosphi;
180 double sinphi = std::sin(
phi);
181 double sinsqphi = sinphi * sinphi;
182 double sin2phi = std::sin(2. *
phi);
183 double sinpsi = std::sqrt(1. - cospsi*cospsi);
184
185 C[0] = 2. * cossqpsi * (1. - sinsqtheta * cossqphi);
186 C[1] = sinsqpsi * (1.-sinsqtheta * sinsqphi);
187 C[2] = sinsqpsi * sinsqtheta;
188 C[3] = 1. / sqrt(2.) * sin2psi * sinsqtheta * sin2phi;
189 C[4] = -1. * sinsqpsi * sin2theta * sinphi;
190 C[5] = 1. / sqrt(2.) * sin2psi * sin2theta * cosphi;
191 C[6] = 2. / 3. *(1. - sinsqtheta * cossqphi);
192 C[7] = 1. / 3. *sqrt(6.) *sinpsi * sinsqtheta * sin2phi;
193 C[8] = 1. / 3. *sqrt(6.) *sinpsi * sin2theta * cosphi;
194 C[9] = 4. / 3. *sqrt(3.) *cospsi * (1. - sinsqtheta * cossqphi);
195
196 }
197
198
199
200
201 double Wplus = 0;
202 double Wminus = 0;
203
204
205 for (
int i=0;
i<10; ++
i ) {
206
207 Wplus +=
A[
i] * ( B1[
i] + B2[
i] ) *
C[i];
208 Wminus +=
A[
i] * ( B1[
i] - B2[
i] ) *
C[i];
209
210 }
211
212 double value = tagprob * Wplus + (1. - tagprob) * Wminus;
213 if (value < 0)
216}
Scalar phi() const
phi method
time(flags, cells_name, *args, **kw)
hold the test vectors and ease the comparison