Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
4 : // This software is part of the EvtGen package developed jointly
5 : // for the BaBar and CLEO collaborations. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtSVVHelCPMix.cc
12 : //
13 : // Description: Routine to decay scalar -> 2 vectors
14 : // by specifying the helicity amplitudes, taking appropriate
15 : // weak phases into account to get mixing and CP violation through
16 : // interference. Based on EvtSVVHelAmp. Particularly appropriate for
17 : // Bs->J/Psi+Phi
18 : //
19 : // Modification history:
20 : //
21 : // RYD November 24, 1996 EvtSVVHelAmp Module
22 : // CATMORE March 2004 Modified to EvtSVVHelCPMix
23 : //
24 : // Model takes the following as user-specified arguments:
25 : // deltaM, averageM - mass diference and average of light and heavy mass eigenstates (real scalars)
26 : // gamma, deltagamma - average width and width difference of the l and h eigenstates (real scalars)
27 : // delta1, delta2 - strong phases (real scalars)
28 : // direct weak phase (real scalar) (for Bs->JPsiPhi this will be zero)
29 : // weak mixing phase (real scalar) (this is equal to 2*arg(Vts Vtb) for Bs->JPsiPhi)
30 : // Magnitudes of helicity amplitudes as in SVV_HELAMP
31 : // See Phys Rev D 34 p1404 - p1417 and chapters 5 and 7 of Physics Reports 370 p537-680 for more details
32 : //------------------------------------------------------------------------
33 : //
34 : #ifdef WIN32
35 : #pragma warning( disable : 4786 )
36 : // Disable anoying warning about symbol size
37 : #endif
38 : #include <iostream>
39 : #include <fstream>
40 : #include <stdlib.h>
41 : #include <ctype.h>
42 : #include "EvtGenBase/EvtParticle.hh"
43 : #include "EvtGenBase/EvtGenKine.hh"
44 : #include "EvtGenBase/EvtPDL.hh"
45 : #include "EvtGenBase/EvtVector4C.hh"
46 : #include "EvtGenBase/EvtTensor4C.hh"
47 : #include "EvtGenBase/EvtVector3C.hh"
48 : #include "EvtGenBase/EvtVector3R.hh"
49 : #include "EvtGenBase/EvtTensor3C.hh"
50 : #include "EvtGenBase/EvtReport.hh"
51 : #include "EvtGenModels/EvtSVVHelCPMix.hh"
52 : #include "EvtGenBase/EvtId.hh"
53 : #include <string>
54 :
55 0 : EvtSVVHelCPMix::~EvtSVVHelCPMix() {}
56 :
57 : std::string EvtSVVHelCPMix::getName(){
58 :
59 0 : return "SVVHELCPMIX";
60 :
61 : }
62 :
63 :
64 : EvtDecayBase* EvtSVVHelCPMix::clone(){
65 :
66 0 : return new EvtSVVHelCPMix;
67 :
68 0 : }
69 :
70 : void EvtSVVHelCPMix::init(){
71 :
72 : // check that there are 12 arguments
73 0 : checkNArg(12);
74 0 : checkNDaug(2);
75 :
76 0 : checkSpinParent(EvtSpinType::SCALAR);
77 :
78 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
79 0 : checkSpinDaughter(1,EvtSpinType::VECTOR);
80 :
81 0 : hp = EvtComplex(getArg(0)*cos(getArg(1)),getArg(0)*sin(getArg(1)));
82 0 : h0 = EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
83 0 : hm = EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
84 0 : averageM = getArg(6);
85 0 : deltaM = getArg(7);
86 0 : gamma = getArg(8);
87 0 : deltagamma = getArg(9);
88 0 : weakmixingphase = EvtComplex(cos(getArg(10)),sin(getArg(10)));
89 0 : weakdirectphase = EvtComplex(cos(getArg(11)),sin(getArg(11)));
90 0 : }
91 :
92 :
93 : void EvtSVVHelCPMix::initProbMax(){
94 :
95 0 : setProbMax(getArg(0)*getArg(0)+getArg(2)*getArg(2)+getArg(4)*getArg(4));
96 :
97 0 : }
98 :
99 :
100 : void EvtSVVHelCPMix::decay( EvtParticle *p){
101 :
102 : EvtParticle* parent = p;
103 0 : EvtAmp& amp = _amp2;
104 0 : EvtId n_v1 = getDaug(0);
105 0 : EvtId n_v2 = getDaug(1);
106 :
107 : // Routine to decay a vector into a vector and scalar. Started
108 : // by ryd on Oct 17, 1996.
109 : // Modified by J.Catmore to take account of CP-violation and mixing
110 :
111 : int tndaug = 2;
112 0 : EvtId tdaug[2];
113 0 : EvtId Bs=EvtPDL::getId("B_s0");
114 0 : EvtId antiBs=EvtPDL::getId("anti-B_s0");
115 0 : tdaug[0] = n_v1;
116 0 : tdaug[1] = n_v2;
117 :
118 : // Phase space and kinematics
119 :
120 0 : parent->initializePhaseSpace(tndaug,tdaug);
121 :
122 : EvtParticle *v1,*v2;
123 0 : v1 = parent->getDaug(0);
124 0 : v2 = parent->getDaug(1);
125 :
126 0 : EvtVector4R momv1 = v1->getP4();
127 :
128 0 : EvtVector3R v1dir(momv1.get(1),momv1.get(2),momv1.get(3));
129 0 : v1dir=v1dir/v1dir.d3mag();
130 :
131 : // Definition of quantities used in construction of complex amplitudes:
132 :
133 0 : EvtTensor3C M; // Tensor as defined in EvtGen manual, equ 117
134 0 : EvtComplex a,b,c; // Helicity amplitudes; EvtGen manual eqns 126-128, also see Phys Lett B 369 p144-150 eqn 15
135 : //EvtComplex deltamu = EvtComplex(deltaM, -0.5*deltagamma); // See Phys Rev D 34 p1404
136 :
137 : // conversion from times in mm/c to natural units [GeV]^-1
138 0 : double t = ((parent->getLifetime())/2.998e11)*6.58e-25;
139 :
140 : // The following two quantities defined in Phys Rev D 34 p1404
141 0 : EvtComplex fplus = EvtComplex(cos(averageM*t),-1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*
142 0 : (cos(0.5*deltaM*t)*cosh(0.25*deltagamma*t)+EvtComplex(0.0,sin(0.5*deltaM*t)*sinh(0.25*deltagamma*t)));
143 0 : EvtComplex fminus = EvtComplex(cos(averageM*t), -1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*EvtComplex(0.0,1.0)*
144 0 : (sin(0.5*deltaM*t)*cosh(0.25*deltagamma*t)-EvtComplex(0.0,1.0)*sinh(0.25*deltagamma*t)*cos(0.5*deltaM*t));
145 :
146 : // See EvtGen manual pp 106-107
147 :
148 0 : a=-0.5*(hp+hm);
149 0 : b=EvtComplex(0.0,0.5)*(hp-hm);
150 0 : c=(h0+0.5*(hp+hm));
151 :
152 0 : M=a*EvtTensor3C::id()+
153 0 : b*EvtGenFunctions::eps(v1dir)+
154 0 : c*EvtGenFunctions::directProd(v1dir,v1dir);
155 :
156 0 : EvtVector3C t0=M.cont1(v1->eps(0).vec().conj());
157 0 : EvtVector3C t1=M.cont1(v1->eps(1).vec().conj());
158 0 : EvtVector3C t2=M.cont1(v1->eps(2).vec().conj());
159 :
160 0 : EvtVector3C eps0=v2->eps(0).vec().conj();
161 0 : EvtVector3C eps1=v2->eps(1).vec().conj();
162 0 : EvtVector3C eps2=v2->eps(2).vec().conj();
163 :
164 : // We need two sets of equations, one for mesons which were in the Bs state at t=0, and another
165 : // for those which were in the antiBs state. Each equation consists of a sum of amplitudes - mod-squaring gives the interference terms.
166 :
167 0 : EvtComplex amplSum00, amplSum01, amplSum02;
168 0 : EvtComplex amplSum10, amplSum11, amplSum12;
169 0 : EvtComplex amplSum20, amplSum21, amplSum22;
170 :
171 : // First the Bs state:
172 :
173 0 : if (parent->getId()==Bs) {
174 :
175 0 : amplSum00 = (fplus*weakdirectphase*t0*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps0);
176 0 : amplSum01 = (fplus*weakdirectphase*t0*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps1);
177 0 : amplSum02 = (fplus*weakdirectphase*t0*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps2);
178 :
179 0 : amplSum10 = (fplus*weakdirectphase*t1*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps0);
180 0 : amplSum11 = (fplus*weakdirectphase*t1*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps1);
181 0 : amplSum12 = (fplus*weakdirectphase*t1*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps2);
182 :
183 0 : amplSum20 = (fplus*weakdirectphase*t2*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps0);
184 0 : amplSum21 = (fplus*weakdirectphase*t2*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps1);
185 0 : amplSum22 = (fplus*weakdirectphase*t2*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps2);
186 0 : }
187 :
188 : // Now the anti-Bs state:
189 :
190 0 : if (parent->getId()==antiBs) {
191 :
192 0 : amplSum00 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps0) + (fplus*(1.0/weakdirectphase)*t0*eps0);
193 0 : amplSum01 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps1) + (fplus*(1.0/weakdirectphase)*t0*eps1);
194 0 : amplSum02 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps2) + (fplus*(1.0/weakdirectphase)*t0*eps2);
195 :
196 0 : amplSum10 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps0) + (fplus*(1.0/weakdirectphase)*t1*eps0);
197 0 : amplSum11 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps1) + (fplus*(1.0/weakdirectphase)*t1*eps1);
198 0 : amplSum12 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps2) + (fplus*(1.0/weakdirectphase)*t1*eps2);
199 :
200 0 : amplSum20 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps0) + (fplus*(1.0/weakdirectphase)*t2*eps0);
201 0 : amplSum21 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps1) + (fplus*(1.0/weakdirectphase)*t2*eps1);
202 0 : amplSum22 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps2) + (fplus*(1.0/weakdirectphase)*t2*eps2);
203 :
204 0 : }
205 :
206 : // Now set the amplitude
207 :
208 0 : amp.vertex(0,0,amplSum00);
209 0 : report(Severity::Info,"EvtGen") << "00: " << amplSum00 << std::endl;
210 0 : amp.vertex(0,1,amplSum01);
211 0 : report(Severity::Info,"EvtGen") << "01: " << amplSum01 << std::endl;
212 0 : amp.vertex(0,2,amplSum02);
213 0 : report(Severity::Info,"EvtGen") << "02: " << amplSum02 << std::endl;
214 :
215 0 : amp.vertex(1,0,amplSum10);
216 0 : report(Severity::Info,"EvtGen") << "10: " << amplSum10 << std::endl;
217 0 : amp.vertex(1,1,amplSum11);
218 0 : report(Severity::Info,"EvtGen") << "11: " << amplSum11 << std::endl;
219 0 : amp.vertex(1,2,amplSum12);
220 0 : report(Severity::Info,"EvtGen") << "12: " << amplSum12 << std::endl;
221 :
222 0 : amp.vertex(2,0,amplSum20);
223 0 : report(Severity::Info,"EvtGen") << "20: " << amplSum20 << std::endl;
224 0 : amp.vertex(2,1,amplSum21);
225 0 : report(Severity::Info,"EvtGen") << "21: " << amplSum21 << std::endl;
226 0 : amp.vertex(2,2,amplSum22);
227 0 : report(Severity::Info,"EvtGen") << "22: " << amplSum22 << std::endl;
228 :
229 : return ;
230 :
231 0 : }
232 :
233 : std::string EvtSVVHelCPMix::getParamName(int i) {
234 0 : switch(i) {
235 : case 0:
236 0 : return "plusHelAmp";
237 : case 1:
238 0 : return "plusHelAmpPhase";
239 : case 2:
240 0 : return "zeroHelAmp";
241 : case 3:
242 0 : return "zeroHelAmpPhase";
243 : case 4:
244 0 : return "minusHelAmp";
245 : case 5:
246 0 : return "minusHelAmpPhase";
247 : case 6:
248 0 : return "averageM";
249 : case 7:
250 0 : return "deltaM";
251 : case 8:
252 0 : return "gamma";
253 : case 9:
254 0 : return "deltaGamma";
255 : case 10:
256 0 : return "weakMixPhase";
257 : case 11:
258 0 : return "weakDirectPhase";
259 : default:
260 0 : return "";
261 : }
262 0 : }
263 :
264 : std::string EvtSVVHelCPMix::getParamDefault(int i) {
265 0 : switch(i) {
266 : case 0:
267 0 : return "1.0";
268 : case 1:
269 0 : return "0.0";
270 : case 2:
271 0 : return "1.0";
272 : case 3:
273 0 : return "0.0";
274 : case 4:
275 0 : return "1.0";
276 : case 5:
277 0 : return "0.0";
278 : default:
279 0 : return "";
280 : }
281 0 : }
|
