]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenModels/EvtBTo4piCP.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBTo4piCP.cxx
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) 1998 Caltech, UCSB
10//
11// Module: EvtBTo4piCP.cc
12//
13// Description: Routine to decay B->pi+ pi- pi+ pi-.
14//
15// Modification history:
16//
17// RYD March 2, 1997 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtGenBase/EvtPatches.hh"
22#include <stdlib.h>
23#include "EvtGenBase/EvtParticle.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtCPUtil.hh"
26#include "EvtGenBase/EvtPDL.hh"
27#include "EvtGenBase/EvtReport.hh"
28#include "EvtGenModels/EvtBTo4piCP.hh"
29#include "EvtGenBase/EvtId.hh"
30#include "EvtGenBase/EvtConst.hh"
31#include <string>
32
33EvtBTo4piCP::~EvtBTo4piCP() {}
34
35
36EvtComplex EvtAmpA2(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
37 const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
38
39 //added by Lange Jan4,2000
40 static EvtId A2M=EvtPDL::getId("a_2-");
41 static EvtId RHO0=EvtPDL::getId("rho0");
42
43 EvtVector4R p4a2,p4rho,p4b;
44
45 p4rho=p4pi1+p4pi2;
46
47 p4a2=p4rho+p4pi3;
48
49 p4b=p4a2+p4pi4;
50
51 EvtVector4R p4b_a2,p4rho_a2,p4pi1_a2,p4a2_a2;
52
53 p4b_a2=boostTo(p4b,p4a2);
54 p4rho_a2=boostTo(p4rho,p4a2);
55 p4pi1_a2=boostTo(p4pi1,p4a2);
56 p4a2_a2=boostTo(p4a2,p4a2);
57
58 EvtVector4R p4pi1_rho;
59
60 p4pi1_rho=boostTo(p4pi1_a2,p4rho_a2);
61
62 EvtVector4R vb,vrho,vpi,t;
63
64 vb=p4b_a2/p4b_a2.d3mag();
65 vrho=p4rho_a2/p4rho_a2.d3mag();
66 vpi=p4pi1_rho/p4pi1_rho.d3mag();
67
68 t.set(1.0,0.0,0.0,0.0);
69
70 // EvtComplex amp_a1,amp_a2;
71 EvtComplex amp_a2;
72
73 // double bwm_a1=EvtPDL::getMeanMass(A1M);
74 // double gamma_a1=EvtPDL::getWidth(A1M);
75 double bwm_a2=EvtPDL::getMeanMass(A2M);
76 double gamma_a2=EvtPDL::getWidth(A2M);
77 double bwm_rho=EvtPDL::getMeanMass(RHO0);
78 double gamma_rho=EvtPDL::getWidth(RHO0);
79
80 amp_a2=(sqrt(gamma_a2/EvtConst::twoPi)/
81 ((p4a2).mass()-bwm_a2-EvtComplex(0.0,0.5*gamma_a2)))*
82 (sqrt(gamma_rho/EvtConst::twoPi)/
83 ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
84
85 return amp_a2*
86 (vb.get(1)*vrho.get(1)+vb.get(2)*vrho.get(2)+vb.get(3)*vrho.get(3))*
87 (
88 vpi.get(1)*(vb.get(2)*vrho.get(3)-vb.get(3)*vrho.get(2))+
89 vpi.get(2)*(vb.get(3)*vrho.get(1)-vb.get(1)*vrho.get(3))+
90 vpi.get(3)*(vb.get(1)*vrho.get(2)-vb.get(2)*vrho.get(1))
91 );
92
93}
94
95EvtComplex EvtAmpA1(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
96 const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
97
98 //added by Lange Jan4,2000
99 static EvtId A1M=EvtPDL::getId("a_1-");
100 static EvtId RHO0=EvtPDL::getId("rho0");
101
102 EvtVector4R p4a1,p4rho,p4b;
103
104 p4rho=p4pi1+p4pi2;
105
106 p4a1=p4rho+p4pi3;
107
108 p4b=p4a1+p4pi4;
109
110 EvtVector4R p4b_a1,p4rho_a1,p4pi1_a1,p4a1_a1;
111
112 p4b_a1=boostTo(p4b,p4a1);
113 p4rho_a1=boostTo(p4rho,p4a1);
114 p4pi1_a1=boostTo(p4pi1,p4a1);
115 p4a1_a1=boostTo(p4a1,p4a1);
116
117 EvtVector4R p4pi1_rho;
118
119 p4pi1_rho=boostTo(p4pi1_a1,p4rho_a1);
120
121 EvtVector4R vb,vrho,vpi,t;
122
123 vb=p4b_a1/p4b_a1.d3mag();
124 vrho=p4rho_a1/p4rho_a1.d3mag();
125 vpi=p4pi1_rho/p4pi1_rho.d3mag();
126
127 t.set(1.0,0.0,0.0,0.0);
128
129 EvtComplex amp_a1;
130
131 double bwm_a1=EvtPDL::getMeanMass(A1M);
132 double gamma_a1=EvtPDL::getWidth(A1M);
133 // double bwm_a2=EvtPDL::getMeanMass(A2M);
134 // double gamma_a2=EvtPDL::getWidth(A2M);
135 double bwm_rho=EvtPDL::getMeanMass(RHO0);
136 double gamma_rho=EvtPDL::getWidth(RHO0);
137
138 amp_a1=(sqrt(gamma_a1/EvtConst::twoPi)/
139 ((p4a1).mass()-bwm_a1-EvtComplex(0.0,0.5*gamma_a1)))*
140 (sqrt(gamma_rho/EvtConst::twoPi)/
141 ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
142
143 return amp_a1*
144 (vb.get(1)*vpi.get(1)+vb.get(2)*vpi.get(2)+vb.get(3)*vpi.get(3));
145
146}
147
148
149std::string EvtBTo4piCP::getName(){
150
151 return "BTO4PI_CP";
152
153}
154
155
156EvtDecayBase* EvtBTo4piCP::clone(){
157
158 return new EvtBTo4piCP;
159
160}
161
162void EvtBTo4piCP::init(){
163
164 // check that there are 18 arguments
165 checkNArg(18);
166 checkNDaug(4);
167
168 checkSpinParent(EvtSpinType::SCALAR);
169
170 checkSpinDaughter(0,EvtSpinType::SCALAR);
171 checkSpinDaughter(1,EvtSpinType::SCALAR);
172 checkSpinDaughter(2,EvtSpinType::SCALAR);
173 checkSpinDaughter(3,EvtSpinType::SCALAR);
174}
175
176void EvtBTo4piCP::decay( EvtParticle *p){
177
178 //added by Lange Jan4,2000
179 static EvtId B0=EvtPDL::getId("B0");
180 static EvtId B0B=EvtPDL::getId("anti-B0");
181
182
183 double t;
184 EvtId other_b;
185
186 EvtCPUtil::OtherB(p,t,other_b);
187
188 p->initializePhaseSpace(getNDaug(),getDaugs());
189 EvtVector4R mom1 = p->getDaug(0)->getP4();
190 EvtVector4R mom2 = p->getDaug(1)->getP4();
191 EvtVector4R mom3 = p->getDaug(2)->getP4();
192 EvtVector4R mom4 = p->getDaug(3)->getP4();
193
194 // double alpha=getArg(0);
195 //double dm=getArg(1);
196
197 EvtComplex amp;
198
199 EvtComplex A,Abar;
200
201
202 EvtComplex A_a1p,Abar_a1p,A_a2p,Abar_a2p;
203 EvtComplex A_a1m,Abar_a1m,A_a2m,Abar_a2m;
204
205 A_a1p=EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
206 Abar_a1p=EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
207
208 A_a2p=EvtComplex(getArg(6)*cos(getArg(7)),getArg(6)*sin(getArg(7)));
209 Abar_a2p=EvtComplex(getArg(8)*cos(getArg(9)),getArg(8)*sin(getArg(9)));
210
211 A_a1m=EvtComplex(getArg(10)*cos(getArg(11)),getArg(10)*sin(getArg(11)));
212 Abar_a1m=EvtComplex(getArg(12)*cos(getArg(13)),getArg(12)*sin(getArg(13)));
213
214 A_a2m=EvtComplex(getArg(14)*cos(getArg(15)),getArg(14)*sin(getArg(15)));
215 Abar_a2m=EvtComplex(getArg(16)*cos(getArg(17)),getArg(16)*sin(getArg(17)));
216
217 EvtComplex a2p_amp=EvtAmpA2(mom1,mom2,mom3,mom4)+
218 EvtAmpA2(mom1,mom4,mom3,mom2)+
219 EvtAmpA2(mom3,mom2,mom1,mom4)+
220 EvtAmpA2(mom3,mom4,mom1,mom2);
221
222 EvtComplex a2m_amp=EvtAmpA2(mom2,mom3,mom4,mom1)+
223 EvtAmpA2(mom2,mom1,mom4,mom3)+
224 EvtAmpA2(mom4,mom3,mom2,mom1)+
225 EvtAmpA2(mom4,mom1,mom2,mom3);
226
227 EvtComplex a1p_amp=EvtAmpA1(mom1,mom2,mom3,mom4)+
228 EvtAmpA1(mom1,mom4,mom3,mom2)+
229 EvtAmpA1(mom3,mom2,mom1,mom4)+
230 EvtAmpA1(mom3,mom4,mom1,mom2);
231
232 EvtComplex a1m_amp=EvtAmpA1(mom2,mom3,mom4,mom1)+
233 EvtAmpA1(mom2,mom1,mom4,mom3)+
234 EvtAmpA1(mom4,mom3,mom2,mom1)+
235 EvtAmpA1(mom4,mom1,mom2,mom3);
236
237
238 A=A_a2p*a2p_amp+A_a1p*a1p_amp+
239 A_a2m*a2m_amp+A_a1m*a1m_amp;
240 Abar=Abar_a2p*a2p_amp+Abar_a1p*a1p_amp+
241 Abar_a2m*a2m_amp+Abar_a1m*a1m_amp;
242
243
244 if (other_b==B0B){
245 amp=A*cos(getArg(1)*t/(2*EvtConst::c))+
246 EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
247 getArg(2)*EvtComplex(0.0,1.0)*Abar*sin(getArg(1)*t/(2*EvtConst::c));
248 }
249 if (other_b==B0){
250 amp=A*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
251 EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
252 getArg(2)*Abar*cos(getArg(1)*t/(2*EvtConst::c));
253 }
254
255 vertex(amp);
256
257 return ;
258}
259