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