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) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtSVSCPiso.cc
12 : //
13 : // Description: Routine to decay scalar -> vectors scalar
14 : // with CP violation and isospin amplitudes.
15 : // More specifically, it is indended to handle
16 : // decays like B->rho pi and B->a1 pi.
17 : //
18 : // Modification history:
19 : //
20 : // RYD/NK Febuary 16, 1998 Module created
21 : //
22 : //------------------------------------------------------------------------
23 : //
24 : #include "EvtGenBase/EvtPatches.hh"
25 : #include <stdlib.h>
26 : #include "EvtGenBase/EvtParticle.hh"
27 : #include "EvtGenBase/EvtRandom.hh"
28 : #include "EvtGenBase/EvtGenKine.hh"
29 : #include "EvtGenBase/EvtCPUtil.hh"
30 : #include "EvtGenBase/EvtPDL.hh"
31 : #include "EvtGenBase/EvtReport.hh"
32 : #include "EvtGenBase/EvtVector4C.hh"
33 : #include "EvtGenModels/EvtSVSCPiso.hh"
34 : #include "EvtGenBase/EvtId.hh"
35 : #include <string>
36 : #include "EvtGenBase/EvtConst.hh"
37 :
38 0 : EvtSVSCPiso::~EvtSVSCPiso() {}
39 :
40 : std::string EvtSVSCPiso::getName(){
41 :
42 0 : return "SVS_CP_ISO";
43 :
44 : }
45 :
46 :
47 : EvtDecayBase* EvtSVSCPiso::clone(){
48 :
49 0 : return new EvtSVSCPiso;
50 :
51 0 : }
52 :
53 : void EvtSVSCPiso::init(){
54 :
55 : // check that there are 27 arguments
56 0 : checkNArg(27);
57 0 : checkNDaug(2);
58 :
59 0 : checkSpinParent(EvtSpinType::SCALAR);
60 :
61 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
62 0 : checkSpinDaughter(1,EvtSpinType::SCALAR);
63 :
64 0 : }
65 :
66 :
67 : void EvtSVSCPiso::initProbMax(){
68 :
69 : //this might need some revision..
70 :
71 0 : if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
72 0 : setProbMax(2.0*(getArg(3)*getArg(3) + 4.0*getArg(23)*getArg(23)));
73 0 : }
74 :
75 0 : if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
76 0 : setProbMax(2.0*(getArg(5)*getArg(5) + 4.0*getArg(25)*getArg(25)));
77 0 : }
78 :
79 0 : if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
80 0 : setProbMax(2.0*(getArg(7)*getArg(7) + 4.0*getArg(23)*getArg(23)));
81 0 : }
82 :
83 0 : if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
84 0 : setProbMax(2.0*(getArg(9)*getArg(9) + 4.0*getArg(25)*getArg(25)));
85 0 : }
86 :
87 0 : if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
88 0 : setProbMax(2.0*(getArg(11)*getArg(11) + getArg(23)*getArg(23) +
89 0 : getArg(19)*getArg(19) + getArg(13)*getArg(13) +
90 0 : getArg(25)*getArg(25) + getArg(21)*getArg(21)));
91 0 : }
92 :
93 0 : if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
94 0 : setProbMax(2.0*(getArg(15)*getArg(15) + getArg(23)*getArg(23) +
95 0 : getArg(19)*getArg(19) + getArg(17)*getArg(17) +
96 0 : getArg(25)*getArg(25) + getArg(21)*getArg(21)));
97 0 : }
98 :
99 0 : if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
100 0 : setProbMax(2.0*(getArg(7)*getArg(7) + getArg(3)*getArg(3) + getArg(11)*getArg(11) +
101 0 : getArg(15)*getArg(15) + 4.0*getArg(19)*getArg(19) + getArg(9)*getArg(9)+
102 0 : getArg(5)*getArg(5) + getArg(13)*getArg(13) + getArg(17)*getArg(17) +
103 0 : 4.0*getArg(21)*getArg(21)));
104 0 : }
105 :
106 0 : }
107 :
108 :
109 : void EvtSVSCPiso::decay( EvtParticle *p){
110 :
111 : //added by Lange Jan4,2000
112 0 : static EvtId B0=EvtPDL::getId("B0");
113 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
114 :
115 0 : double t;
116 0 : EvtId other_b;
117 : int charged(0);
118 :
119 : int first_time=0;
120 : int flip=0;
121 0 : EvtId ds[2];
122 :
123 :
124 : //randomly generate the tag (B0 or B0B)
125 :
126 0 : double tag = EvtRandom::Flat(0.0,1.0);
127 0 : if (tag < 0.5) {
128 :
129 0 : EvtCPUtil::getInstance()->OtherB(p,t,other_b,1.0);
130 0 : other_b = B0;
131 0 : }
132 : else {
133 :
134 0 : EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.0);
135 0 : other_b = B0B;
136 : }
137 :
138 0 : if (p->getNDaug()==0) first_time=1;
139 :
140 0 : if (first_time){
141 0 : if (EvtRandom::Flat(0.0,1.0)<getArg(2)) flip=1;
142 : }
143 : else{
144 0 : if (getDaug(0)!=p->getDaug(0)->getId()) flip=1;
145 : }
146 :
147 0 : if (!flip) {
148 0 : ds[0]=getDaug(0);
149 0 : ds[1]=getDaug(1);
150 0 : }
151 : else{
152 0 : ds[0]=EvtPDL::chargeConj(getDaug(0));
153 0 : ds[1]=EvtPDL::chargeConj(getDaug(1));
154 : }
155 :
156 0 : p->initializePhaseSpace(getNDaug(),ds);
157 :
158 : EvtParticle *v,*s;
159 0 : v=p->getDaug(0);
160 0 : s=p->getDaug(1);
161 :
162 0 : EvtComplex amp;
163 :
164 0 : EvtComplex A_f,Abar_f;
165 0 : EvtComplex A_fbar,Abar_fbar;
166 0 : EvtComplex Apm, Apm_bar, Amp, Amp_bar;
167 :
168 0 : EvtComplex Tp0, Tp0_bar, T0p, T0p_bar,Tpm, Tpm_bar, Tmp, Tmp_bar;
169 0 : EvtComplex P1, P1_bar, P0, P0_bar;
170 :
171 0 : Tp0 = EvtComplex(getArg(3)*cos(getArg(4)),getArg(3)*sin(getArg(4)));
172 0 : Tp0_bar = EvtComplex(getArg(5)*cos(getArg(6)),getArg(5)*sin(getArg(6)));
173 0 : T0p = EvtComplex(getArg(7)*cos(getArg(8)),getArg(7)*sin(getArg(8)));
174 0 : T0p_bar = EvtComplex(getArg(9)*cos(getArg(10)),getArg(9)*sin(getArg(10)));
175 0 : Tpm = EvtComplex(getArg(11)*cos(getArg(12)),getArg(11)*sin(getArg(12)));
176 0 : Tpm_bar = EvtComplex(getArg(13)*cos(getArg(14)),getArg(13)*sin(getArg(14)));
177 0 : Tmp = EvtComplex(getArg(15)*cos(getArg(16)),getArg(15)*sin(getArg(16)));
178 0 : Tmp_bar = EvtComplex(getArg(17)*cos(getArg(18)),getArg(17)*sin(getArg(18)));
179 0 : P0 = EvtComplex(getArg(19)*cos(getArg(20)),getArg(19)*sin(getArg(20)));
180 0 : P0_bar = EvtComplex(getArg(21)*cos(getArg(22)),getArg(21)*sin(getArg(22)));
181 0 : P1 = EvtComplex(getArg(23)*cos(getArg(24)),getArg(23)*sin(getArg(24)));
182 0 : P1_bar = EvtComplex(getArg(25)*cos(getArg(26)),getArg(25)*sin(getArg(26)));
183 :
184 :
185 : //***********************charged modes****************************
186 :
187 0 : if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
188 :
189 : //V+ S0, so T+0 + 2 P1
190 :
191 : charged = 1;
192 0 : A_f = Tp0 + 2.0*P1;
193 0 : }
194 :
195 0 : if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
196 :
197 : //V- S0, so T+0_bar + 2P1_bar
198 :
199 : charged = 1;
200 0 : A_f = Tp0_bar + 2.0*P1_bar;
201 0 : }
202 :
203 0 : if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
204 :
205 : //V0 S+, so T0+ - 2 P1
206 :
207 : charged = 1;
208 0 : A_f = T0p - 2.0*P1;
209 0 : }
210 :
211 0 : if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
212 :
213 : //V0 S-, so T0+_bar - 2 P1_bar
214 :
215 : charged = 1;
216 0 : A_f = T0p_bar - 2.0*P1_bar;
217 0 : }
218 :
219 :
220 : //***********************neutral modes***************************
221 :
222 :
223 : //V+ S-, so Af = T+- + P1 + P0
224 0 : Apm = Tpm + P1 + P0;
225 0 : Apm_bar = Tpm_bar + P1_bar + P0_bar;
226 :
227 : //V- S+, so Af = T-+ - P1 + P0
228 0 : Amp = Tmp - P1 + P0;
229 0 : Amp_bar = Tmp_bar - P1_bar + P0;
230 :
231 :
232 0 : if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
233 :
234 : //V+ S-
235 : charged = 0;
236 0 : A_f = Apm;
237 0 : Abar_f = Apm_bar;
238 0 : A_fbar = Amp;
239 0 : Abar_fbar = Amp_bar;
240 :
241 0 : }
242 :
243 0 : if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
244 :
245 : //V- S+
246 : charged = 0;
247 0 : A_f = Amp;
248 0 : Abar_f = Amp_bar;
249 0 : A_fbar = Apm;
250 0 : Abar_fbar = Apm_bar;
251 :
252 0 : }
253 :
254 0 : if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
255 :
256 : //V0 S0
257 : charged = 0;
258 0 : A_f = T0p + Tp0 - Tpm - Tmp - 2.0*P0 ;
259 0 : Abar_f = T0p_bar + Tp0_bar - Tpm_bar - Tmp_bar - 2.0*P0_bar;
260 0 : A_fbar = A_f;
261 0 : Abar_fbar = Abar_f;
262 :
263 0 : }
264 :
265 0 : if (charged==0) {
266 :
267 0 : if (!flip) {
268 0 : if (other_b==B0B){
269 :
270 0 : amp=A_f*cos(getArg(1)*t/(2*EvtConst::c))+
271 0 : EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
272 0 : EvtComplex(0.0,1.0)*Abar_f*sin(getArg(1)*t/(2*EvtConst::c));
273 0 : }
274 0 : if (other_b==B0){
275 :
276 0 : amp=A_f*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
277 0 : EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
278 0 : Abar_f*cos(getArg(1)*t/(2*EvtConst::c));
279 0 : }
280 : }
281 : else{
282 0 : if (other_b==B0B){
283 :
284 0 : amp=A_fbar*cos(getArg(1)*t/(2*EvtConst::c))+
285 0 : EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
286 0 : EvtComplex(0.0,1.0)*Abar_fbar*sin(getArg(1)*t/(2*EvtConst::c));
287 0 : }
288 0 : if (other_b==B0){
289 :
290 0 : amp=A_fbar*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
291 0 : EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
292 0 : Abar_fbar*cos(getArg(1)*t/(2*EvtConst::c));
293 0 : }
294 : }
295 :
296 : }
297 0 : else amp = A_f;
298 :
299 0 : EvtVector4R p4_parent;
300 :
301 0 : p4_parent=v->getP4()+s->getP4();
302 :
303 0 : double norm=1.0/v->getP4().d3mag();
304 :
305 0 : vertex(0,amp*norm*p4_parent*(v->epsParent(0)));
306 0 : vertex(1,amp*norm*p4_parent*(v->epsParent(1)));
307 0 : vertex(2,amp*norm*p4_parent*(v->epsParent(2)));
308 :
309 : return ;
310 0 : }
311 :
|
