Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtSSDCP.cpp
CommitLineData
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"
35using std::endl;
36
37EvtSSDCP::~EvtSSDCP() {}
38
39std::string EvtSSDCP::getName(){
40
41 return "SSD_CP";
42
43}
44
45
46EvtDecayBase* EvtSSDCP::clone(){
47
48 return new EvtSSDCP;
49
50}
51
52void 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
161void 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
177void 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 316std::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
351std::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}