]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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) 2001 Caltech | |
10 | // | |
11 | // Module: EvtSSDCP.cc | |
12 | // | |
13 | // Description: See EvtSSDCP.hh | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // RYD August 12, 2001 Module created | |
18 | // F. Sandrelli, Fernando M-V March 1, 2002 Debugged and added z parameter (CPT violation) | |
19 | //------------------------------------------------------------------------ | |
20 | // | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | #include <stdlib.h> | |
23 | #include "EvtGenBase/EvtParticle.hh" | |
24 | #include "EvtGenBase/EvtRandom.hh" | |
25 | #include "EvtGenBase/EvtGenKine.hh" | |
26 | #include "EvtGenBase/EvtCPUtil.hh" | |
27 | #include "EvtGenBase/EvtPDL.hh" | |
28 | #include "EvtGenBase/EvtReport.hh" | |
29 | #include "EvtGenBase/EvtVector4C.hh" | |
30 | #include "EvtGenBase/EvtTensor4C.hh" | |
31 | #include "EvtGenModels/EvtSSDCP.hh" | |
32 | #include <string> | |
33 | #include "EvtGenBase/EvtConst.hh" | |
34 | using std::endl; | |
35 | ||
36 | EvtSSDCP::~EvtSSDCP() {} | |
37 | ||
38 | std::string EvtSSDCP::getName(){ | |
39 | ||
40 | return "SSD_CP"; | |
41 | ||
42 | } | |
43 | ||
44 | ||
45 | EvtDecayBase* EvtSSDCP::clone(){ | |
46 | ||
47 | return new EvtSSDCP; | |
48 | ||
49 | } | |
50 | ||
51 | void EvtSSDCP::init(){ | |
52 | ||
53 | // check that there are 8 or 12 or 14 arguments | |
54 | ||
55 | checkNArg(14,12,8); | |
56 | checkNDaug(2); | |
57 | ||
58 | EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0)); | |
59 | EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1)); | |
60 | ||
61 | // Check it is a B0 or B0s | |
62 | if ( ( getParentId() != EvtPDL::getId( "B0" ) ) | |
63 | && ( getParentId() != EvtPDL::getId( "anti-B0" ) ) | |
64 | && ( getParentId() != EvtPDL::getId( "B_s0" ) ) | |
65 | && ( getParentId() != EvtPDL::getId( "anti-B_s0" ) ) ) { | |
66 | report( ERROR , "EvtGen" ) << "EvtSSDCP only decays B0 and B0s" | |
67 | << std::endl ; | |
68 | ::abort() ; | |
69 | } | |
70 | ||
71 | ||
72 | if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))|| | |
73 | (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||(d2type==EvtSpinType::TENSOR)))|| | |
74 | (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||(d1type==EvtSpinType::TENSOR))) | |
75 | ) { | |
76 | report(ERROR,"EvtGen") << "EvtSSDCP generator expected " | |
77 | << "one of the daugters to be a scalar, the other either scalar, vector, or tensor, found:" | |
78 | << EvtPDL::name(getDaug(0)).c_str()<<" and "<<EvtPDL::name(getDaug(1)).c_str()<<endl; | |
79 | report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; | |
80 | ::abort(); | |
81 | } | |
82 | ||
83 | _dm=getArg(0)/EvtConst::c; //units of 1/mm | |
84 | ||
85 | _dgog=getArg(1); | |
86 | ||
87 | ||
88 | _qoverp=getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3))); | |
89 | _poverq=1.0/_qoverp; | |
90 | ||
91 | _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5))); | |
92 | ||
93 | _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7))); | |
94 | ||
95 | if (getNArg()>=12){ | |
96 | _eigenstate=false; | |
97 | _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9))); | |
98 | _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11))); | |
99 | } | |
100 | else{ | |
101 | //I'm somewhat confused about this. For a CP eigenstate set the | |
102 | //amplitudes to the same. For a non CP eigenstate CPT invariance | |
103 | //is enforced. (ryd) | |
104 | if ( | |
105 | (getDaug(0)==EvtPDL::chargeConj(getDaug(0))&& | |
106 | getDaug(1)==EvtPDL::chargeConj(getDaug(1)))|| | |
107 | (getDaug(0)==EvtPDL::chargeConj(getDaug(1))&& | |
108 | getDaug(1)==EvtPDL::chargeConj(getDaug(0)))){ | |
109 | _eigenstate=true; | |
110 | }else{ | |
111 | _eigenstate=false; | |
112 | _A_fbar=conj(_Abar_f); | |
113 | _Abar_fbar=conj(_A_f); | |
114 | } | |
115 | } | |
116 | ||
117 | //FS: new check for z | |
118 | if (getNArg()==14){ //FS Set _z parameter if provided else set it 0 | |
119 | _z=EvtComplex(getArg(12),getArg(13)); | |
120 | } | |
121 | else{ | |
122 | _z=EvtComplex(0.0,0.0); | |
123 | } | |
124 | ||
125 | // FS substituted next 2 lines... | |
126 | ||
127 | ||
128 | // | |
129 | // _gamma=EvtPDL::getctau(EvtPDL::getId("B0")); //units of 1/mm | |
130 | //_dgamma=_gamma*0.5*_dgog; | |
131 | // | |
132 | // ...with: | |
133 | ||
134 | if ( ( getParentId() == EvtPDL::getId("B0") ) || | |
135 | ( getParentId() == EvtPDL::getId("anti-B0") ) ) { | |
136 | _gamma=1./EvtPDL::getctau(EvtPDL::getId("B0")); //gamma/c (1/mm) | |
137 | } | |
138 | else{ | |
139 | _gamma=1./EvtPDL::getctau(EvtPDL::getId("B_s0")) ; | |
140 | } | |
141 | ||
142 | _dgamma=_gamma*_dgog; //dgamma/c (1/mm) | |
143 | ||
144 | if (verbose()){ | |
145 | report(INFO,"EvtGen") << "SSD_CP will generate CP/CPT violation:" | |
146 | << endl << endl | |
147 | << " " << EvtPDL::name(getParentId()).c_str() << " --> " | |
148 | << EvtPDL::name(getDaug(0)).c_str() << " + " | |
149 | << EvtPDL::name(getDaug(1)).c_str() << endl << endl | |
150 | << "using parameters:" << endl << endl | |
151 | << " delta(m) = " << _dm << " hbar/ps" << endl | |
152 | << "dGamma = " << _dgamma <<" ps-1" <<endl | |
153 | << " q/p = " << _qoverp << endl | |
154 | << " z = " << _z << endl | |
155 | << " tau = " << 1./_gamma << " ps" << endl; | |
156 | } | |
157 | ||
158 | } | |
159 | ||
160 | void EvtSSDCP::initProbMax() { | |
161 | double theProbMax = | |
162 | abs(_A_f) * abs(_A_f) + | |
163 | abs(_Abar_f) * abs(_Abar_f) + | |
164 | abs(_A_fbar) * abs(_A_fbar) + | |
165 | abs(_Abar_fbar) * abs(_Abar_fbar); | |
166 | ||
167 | if (_eigenstate) theProbMax*=2; | |
168 | ||
169 | EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1)); | |
170 | EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0)); | |
171 | if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10; | |
172 | ||
173 | setProbMax(theProbMax); | |
174 | } | |
175 | ||
176 | void EvtSSDCP::decay( EvtParticle *p){ | |
177 | ||
178 | static EvtId B0=EvtPDL::getId("B0"); | |
179 | static EvtId B0B=EvtPDL::getId("anti-B0"); | |
180 | ||
181 | static EvtId B0s = EvtPDL::getId("B_s0"); | |
182 | static EvtId B0Bs = EvtPDL::getId("anti-B_s0"); | |
183 | ||
184 | double t; | |
185 | EvtId other_b; | |
186 | EvtId daugs[2]; | |
187 | ||
188 | int flip=0; | |
189 | if (!_eigenstate){ | |
190 | if (EvtRandom::Flat(0.0,1.0)<0.5) flip=1; | |
191 | } | |
192 | ||
193 | if (!flip) { | |
194 | daugs[0]=getDaug(0); | |
195 | daugs[1]=getDaug(1); | |
196 | } | |
197 | else{ | |
198 | daugs[0]=EvtPDL::chargeConj(getDaug(0)); | |
199 | daugs[1]=EvtPDL::chargeConj(getDaug(1)); | |
200 | } | |
201 | ||
202 | EvtParticle *d; | |
203 | p->initializePhaseSpace(2, daugs); | |
204 | ||
205 | EvtComplex amp; | |
206 | ||
207 | ||
208 | EvtCPUtil::OtherB(p,t,other_b,0.5); // t is c*Dt (mm) | |
209 | ||
210 | //if (flip) t=-t; | |
211 | ||
212 | //FS We assume DGamma=GammaLow-GammaHeavy and Dm=mHeavy-mLow | |
213 | EvtComplex expL=exp(-EvtComplex(-0.25*_dgamma*t,0.5*_dm*t)); | |
214 | EvtComplex expH=exp(EvtComplex(-0.25*_dgamma*t,0.5*_dm*t)); | |
215 | //FS Definition of gp and gm | |
216 | EvtComplex gp=0.5*(expL+expH); | |
217 | EvtComplex gm=0.5*(expL-expH); | |
218 | //FS Calculation os sqrt(1-z^2) | |
219 | EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2)); | |
220 | ||
221 | //EvtComplex BB=0.5*(expL+expH); // <B0|B0(t)> | |
222 | //EvtComplex barBB=_qoverp*0.5*(expL-expH); // <B0bar|B0(t)> | |
223 | //EvtComplex BbarB=_poverq*0.5*(expL-expH); // <B0|B0bar(t)> | |
224 | //EvtComplex barBbarB=BB; // <B0bar|B0bar(t)> | |
225 | // FS redefinition of these guys... (See BAD #188 eq.35 for ref.) | |
226 | // q/p is taken as in the BaBar Phys. Book (opposite sign wrt ref.) | |
227 | EvtComplex BB=gp+_z*gm; // <B0|B0(t)> | |
228 | EvtComplex barBB=sqz*_qoverp*gm; // <B0bar|B0(t)> | |
229 | EvtComplex BbarB=sqz*_poverq*gm; // <B0|B0bar(t)> | |
230 | EvtComplex barBbarB=gp-_z*gm; // <B0bar|B0bar(t)> | |
231 | ||
232 | if (!flip){ | |
233 | if (other_b==B0B||other_b==B0Bs){ | |
234 | //at t=0 we have a B0 | |
235 | //report(INFO,"EvtGen") << "B0B"<<endl; | |
236 | amp=BB*_A_f+BbarB*_Abar_f; | |
237 | //std::cout << "noflip B0B tag:"<<amp<<std::endl; | |
238 | //amp=0.0; | |
239 | } | |
240 | if (other_b==B0||other_b==B0s){ | |
241 | //report(INFO,"EvtGen") << "B0"<<endl; | |
242 | amp=barBB*_A_f+barBbarB*_Abar_f; | |
243 | } | |
244 | }else{ | |
245 | if (other_b==B0||other_b==B0s){ | |
246 | amp=barBB*_A_fbar+barBbarB*_Abar_fbar; | |
247 | //std::cout << "flip B0 tag:"<<amp<<std::endl; | |
248 | //amp=0.0; | |
249 | } | |
250 | if (other_b==B0B||other_b==B0Bs){ | |
251 | amp=BB*_A_fbar+BbarB*_Abar_fbar; | |
252 | } | |
253 | } | |
254 | ||
255 | ||
256 | EvtVector4R p4_parent=p->getP4Restframe(); | |
257 | double m_parent=p4_parent.mass(); | |
258 | ||
259 | EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1)); | |
260 | ||
261 | EvtVector4R momv; | |
262 | EvtVector4R moms; | |
263 | ||
264 | if (d2type==EvtSpinType::SCALAR){ | |
265 | d2type=EvtPDL::getSpinType(getDaug(0)); | |
266 | d= p->getDaug(0); | |
267 | momv = d->getP4(); | |
268 | moms = p->getDaug(1)->getP4(); | |
269 | } | |
270 | else{ | |
271 | d= p->getDaug(1); | |
272 | momv = d->getP4(); | |
273 | moms = p->getDaug(0)->getP4(); | |
274 | } | |
275 | ||
276 | ||
277 | ||
278 | if (d2type==EvtSpinType::SCALAR) { | |
279 | vertex(amp); | |
280 | } | |
281 | ||
282 | if (d2type==EvtSpinType::VECTOR) { | |
283 | ||
284 | double norm=momv.mass()/(momv.d3mag()*p->mass()); | |
285 | ||
286 | //std::cout << amp << " " << norm << " " << p4_parent << d->getP4()<< std::endl; | |
287 | // std::cout << EvtPDL::name(d->getId()) << " " << EvtPDL::name(p->getDaug(0)->getId()) << | |
288 | // " 1and2 " << EvtPDL::name(p->getDaug(1)->getId()) << std::endl; | |
289 | //std::cout << d->eps(0) << std::endl; | |
290 | //std::cout << d->epsParent(0) << std::endl; | |
291 | vertex(0,amp*norm*p4_parent*(d->epsParent(0))); | |
292 | vertex(1,amp*norm*p4_parent*(d->epsParent(1))); | |
293 | vertex(2,amp*norm*p4_parent*(d->epsParent(2))); | |
294 | ||
295 | } | |
296 | ||
297 | if (d2type==EvtSpinType::TENSOR) { | |
298 | ||
299 | double norm=d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag()); | |
300 | ||
301 | ||
302 | vertex(0,amp*norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent); | |
303 | vertex(1,amp*norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent); | |
304 | vertex(2,amp*norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent); | |
305 | vertex(3,amp*norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent); | |
306 | vertex(4,amp*norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent); | |
307 | ||
308 | } | |
309 | ||
310 | ||
311 | return ; | |
312 | } | |
313 |