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