Line data Source code
1 : // $Id: EvtSSD_DirectCP.cpp,v 1.2 2009-03-16 16:24:05 robbep Exp $
2 : // Generation of direct CP violation in hadronic environment
3 : // Patrick Robbe, LHCb, 08 Nov 2006
4 : //
5 : #include <stdlib.h>
6 : #include "EvtGenBase/EvtParticle.hh"
7 : #include "EvtGenBase/EvtRandom.hh"
8 : #include "EvtGenBase/EvtGenKine.hh"
9 : #include "EvtGenBase/EvtCPUtil.hh"
10 : #include "EvtGenBase/EvtPDL.hh"
11 : #include "EvtGenBase/EvtReport.hh"
12 : #include "EvtGenBase/EvtVector4C.hh"
13 : #include "EvtGenBase/EvtTensor4C.hh"
14 : #include "EvtGenModels/EvtSSD_DirectCP.hh"
15 : #include <string>
16 : #include "EvtGenBase/EvtConst.hh"
17 :
18 0 : EvtSSD_DirectCP::~EvtSSD_DirectCP() {}
19 :
20 : std::string EvtSSD_DirectCP::getName( ){
21 :
22 0 : return "SSD_DirectCP" ;
23 :
24 : }
25 :
26 :
27 : EvtDecayBase* EvtSSD_DirectCP::clone(){
28 :
29 0 : return new EvtSSD_DirectCP;
30 :
31 0 : }
32 :
33 : void EvtSSD_DirectCP::init(){
34 :
35 : // check that there is 1 argument and 2-body decay
36 :
37 0 : checkNArg(1);
38 0 : checkNDaug(2);
39 :
40 0 : EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
41 0 : EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
42 :
43 0 : if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
44 0 : (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||
45 0 : (d2type==EvtSpinType::TENSOR)))||
46 0 : (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||
47 0 : (d1type==EvtSpinType::TENSOR)))
48 : ) {
49 0 : report(Severity::Error,"EvtGen") << "EvtSSD_DirectCP generator expected "
50 0 : << "one of the daugters to be a scalar, "
51 0 : << "the other either scalar, vector, or tensor, "
52 0 : << "found:"
53 0 : << EvtPDL::name(getDaug(0)).c_str()
54 0 : <<" and "
55 0 : <<EvtPDL::name(getDaug(1)).c_str()<<std::endl;
56 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<std::endl;
57 0 : ::abort();
58 : }
59 :
60 0 : _acp = getArg( 0 ) ; // A_CP defined as A_CP = (BR(fbar)-BR(f))/(BR(fbar)+BR(f))
61 :
62 0 : }
63 :
64 : void EvtSSD_DirectCP::initProbMax() {
65 : double theProbMax = 1. ;
66 :
67 0 : EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
68 0 : EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
69 0 : if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
70 :
71 0 : setProbMax(theProbMax);
72 0 : }
73 :
74 : void EvtSSD_DirectCP::decay( EvtParticle *p) {
75 :
76 : bool flip = false ;
77 0 : EvtId daugs[2];
78 :
79 : // decide it is B or Bbar:
80 0 : if ( EvtRandom::Flat(0.,1.) < ( ( 1. - _acp ) / 2. ) ) {
81 : // it is a B
82 0 : if ( EvtPDL::getStdHep( getParentId() ) < 0 ) flip = true ;
83 : } else {
84 : // it is a Bbar
85 0 : if ( EvtPDL::getStdHep( getParentId() ) > 0 ) flip = true ;
86 : }
87 :
88 0 : if ( flip ) {
89 0 : if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
90 0 : p->getParent()
91 0 : ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
92 0 : p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
93 0 : }
94 : else {
95 0 : p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
96 : }
97 : }
98 :
99 0 : if (!flip) {
100 0 : daugs[0]=getDaug(0);
101 0 : daugs[1]=getDaug(1);
102 0 : }
103 : else{
104 0 : daugs[0]=EvtPDL::chargeConj(getDaug(0));
105 0 : daugs[1]=EvtPDL::chargeConj(getDaug(1));
106 : }
107 :
108 : EvtParticle *d;
109 0 : p->initializePhaseSpace(2, daugs);
110 :
111 0 : EvtVector4R p4_parent=p->getP4Restframe();
112 0 : double m_parent=p4_parent.mass();
113 :
114 0 : EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
115 :
116 0 : EvtVector4R momv;
117 0 : EvtVector4R moms;
118 :
119 0 : if (d2type==EvtSpinType::SCALAR){
120 0 : d2type=EvtPDL::getSpinType(getDaug(0));
121 0 : d= p->getDaug(0);
122 0 : momv = d->getP4();
123 0 : moms = p->getDaug(1)->getP4();
124 0 : }
125 : else{
126 0 : d= p->getDaug(1);
127 0 : momv = d->getP4();
128 0 : moms = p->getDaug(0)->getP4();
129 : }
130 :
131 0 : if (d2type==EvtSpinType::SCALAR) {
132 0 : vertex(1.);
133 0 : }
134 :
135 0 : if (d2type==EvtSpinType::VECTOR) {
136 :
137 0 : double norm=momv.mass()/(momv.d3mag()*p->mass());
138 :
139 0 : vertex(0,norm*p4_parent*(d->epsParent(0)));
140 0 : vertex(1,norm*p4_parent*(d->epsParent(1)));
141 0 : vertex(2,norm*p4_parent*(d->epsParent(2)));
142 :
143 0 : }
144 :
145 0 : if (d2type==EvtSpinType::TENSOR) {
146 :
147 : double norm=
148 0 : d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
149 :
150 :
151 0 : vertex(0,norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
152 0 : vertex(1,norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
153 0 : vertex(2,norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
154 0 : vertex(3,norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
155 0 : vertex(4,norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);
156 0 : }
157 0 : }
158 :
159 : bool EvtSSD_DirectCP::isB0Mixed ( EvtParticle * p ) {
160 0 : if ( ! ( p->getParent() ) ) return false ;
161 :
162 0 : static EvtId B0 =EvtPDL::getId("B0");
163 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
164 :
165 0 : if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
166 :
167 0 : if ( ( p->getParent()->getId() == B0 ) ||
168 0 : ( p->getParent()->getId() == B0B ) ) return true ;
169 :
170 0 : return false ;
171 0 : }
172 :
173 : bool EvtSSD_DirectCP::isBsMixed ( EvtParticle * p ) {
174 0 : if ( ! ( p->getParent() ) ) return false ;
175 :
176 0 : static EvtId BS0=EvtPDL::getId("B_s0");
177 0 : static EvtId BSB=EvtPDL::getId("anti-B_s0");
178 :
179 0 : if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
180 :
181 0 : if ( ( p->getParent()->getId() == BS0 ) ||
182 0 : ( p->getParent()->getId() == BSB ) ) return true ;
183 :
184 0 : return false ;
185 0 : }
186 :
187 : std::string EvtSSD_DirectCP::getParamName(int i) {
188 0 : switch(i) {
189 : case 0:
190 0 : return "ACP";
191 : default:
192 0 : return "";
193 : }
194 0 : }
|
