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: EvtGen/EvtYmSToYnSpipiCLEO.hh
12 : //
13 : // Description: This model is based on matrix element method used by
14 : // CLEO in Phys.Rev.D76:072001,2007 to model the dipion mass
15 : // and helicity angle distribution in the decays Y(mS) -> pi pi Y(nS),
16 : // where m,n are integers and m>n and m<4.
17 : // This model has two parameters, Re(B/A) and Im(B/A), which
18 : // are coefficients of the matrix element's terms determined by
19 : // the CLEO fits.
20 : //
21 : // Example:
22 : //
23 : // Decay Upsilon(3S)
24 : // 1.0000 Upsilon pi+ pi- YMSTOYNSPIPICLEO -2.523 1.189;
25 : // Enddecay
26 : // Decay Upsilon(3S)
27 : // 1.0000 Upsilon(2S) pi+ pi- YMSTOYNSPIPICLEO -0.395 0.001;
28 : // Enddecay
29 : // Decay Upsilon(2S)
30 : // 1.0000 Upsilon pi+ pi- YMSTOYNSPIPICLEO -0.753 0.000;
31 : // Enddecay
32 : //
33 : // --> the order of parameters is: Re(B/A) Im(B/A)
34 : //
35 : // Modification history:
36 : //
37 : // SEKULA Jan. 28, 2008 Module created
38 : //
39 : //------------------------------------------------------------------------
40 :
41 :
42 : #include "EvtGenBase/EvtPatches.hh"
43 : #include <stdlib.h>
44 : #include "EvtGenBase/EvtParticle.hh"
45 : #include "EvtGenBase/EvtGenKine.hh"
46 : #include "EvtGenBase/EvtPDL.hh"
47 : #include "EvtGenBase/EvtVector4C.hh"
48 : #include "EvtGenBase/EvtTensor4C.hh"
49 : #include "EvtGenBase/EvtReport.hh"
50 : #include "EvtGenBase/EvtRandom.hh"
51 : #include "EvtGenModels/EvtYmSToYnSpipiCLEO.hh"
52 : #include <string>
53 : using std::endl;
54 :
55 0 : EvtYmSToYnSpipiCLEO::~EvtYmSToYnSpipiCLEO() {}
56 :
57 : std::string EvtYmSToYnSpipiCLEO::getName(){
58 :
59 0 : return "YMSTOYNSPIPICLEO";
60 :
61 : }
62 :
63 :
64 : EvtDecayBase* EvtYmSToYnSpipiCLEO::clone(){
65 :
66 0 : return new EvtYmSToYnSpipiCLEO;
67 :
68 0 : }
69 :
70 : void EvtYmSToYnSpipiCLEO::init(){
71 :
72 0 : static EvtId PIP=EvtPDL::getId("pi+");
73 0 : static EvtId PIM=EvtPDL::getId("pi-");
74 0 : static EvtId PI0=EvtPDL::getId("pi0");
75 :
76 : // check that there are 2 arguments
77 0 : checkNArg(2);
78 0 : checkNDaug(3);
79 :
80 0 : checkSpinParent(EvtSpinType::VECTOR);
81 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
82 :
83 :
84 :
85 0 : if ((!(getDaug(1)==PIP&&getDaug(2)==PIM))&&
86 0 : (!(getDaug(1)==PI0&&getDaug(2)==PI0))) {
87 0 : report(Severity::Error,"EvtGen") << "EvtYmSToYnSpipiCLEO generator expected "
88 0 : << " pi+ and pi- (or pi0 and pi0) "
89 0 : << "as 2nd and 3rd daughter. "<<endl;
90 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
91 0 : ::abort();
92 : }
93 :
94 0 : }
95 :
96 : void EvtYmSToYnSpipiCLEO::initProbMax() {
97 0 : setProbMax(2.0);
98 0 : }
99 :
100 : void EvtYmSToYnSpipiCLEO::decay( EvtParticle *p){
101 :
102 :
103 : // We want to simulate the following process:
104 : //
105 : // Y(mS) -> Y(nS) X, X -> pi+ pi- (pi0 pi0)
106 : //
107 : // The CLEO analysis assumed such an intermediate process
108 : // were occurring, and wrote down the matrix element
109 : // and its components according to this assumption.
110 : //
111 : //
112 :
113 :
114 0 : double ReB_over_A = getArg(0);
115 0 : double ImB_over_A = getArg(1);
116 :
117 0 : p->makeDaughters(getNDaug(),getDaugs());
118 : EvtParticle *v,*s1,*s2;
119 0 : v=p->getDaug(0);
120 0 : s1=p->getDaug(1);
121 0 : s2=p->getDaug(2);
122 :
123 0 : double m_pi = s1->getP4().mass();
124 0 : double M_mS = p->getP4().mass();
125 0 : double M_nS = v->getP4().mass();
126 :
127 : // // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "M_nS = " << v->getP4().mass() << endl;
128 :
129 0 : EvtVector4R P_nS;
130 0 : EvtVector4R P_pi1;
131 0 : EvtVector4R P_pi2;
132 :
133 : // Perform a simple accept/reject until we get a configuration of the
134 : // dipion system that passes
135 : bool acceptX = false;
136 :
137 0 : while( false == acceptX )
138 : {
139 :
140 : // Begin by generating a random X mass between the kinematic
141 : // boundaries, 2*m_pi and M(mS) - M(nS)
142 :
143 0 : double mX = EvtRandom::Flat(2.0 * m_pi, M_mS-M_nS);
144 :
145 : // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "m_X = " << mX << endl;
146 :
147 : // Now create a two-body decay from the Y(mS) in its rest frame
148 : // of Y(mS) -> Y(nS) + X
149 :
150 0 : double masses[2];
151 0 : masses[0] = M_nS;
152 0 : masses[1] = mX;
153 :
154 0 : EvtVector4R p4[2];
155 :
156 0 : EvtGenKine::PhaseSpace( 2, masses, p4, M_mS );
157 :
158 0 : P_nS = p4[0];
159 0 : EvtVector4R P_X = p4[1];
160 :
161 : // Now create the four-vectors for the two pions in the X
162 : // rest frame, X -> pi pi
163 :
164 0 : masses[0] = s1->mass();
165 0 : masses[1] = s2->mass();
166 :
167 0 : EvtGenKine::PhaseSpace( 2, masses, p4, P_X.mass() );
168 :
169 : // compute cos(theta), the polar helicity angle between a pi+ and
170 : // the direction opposite the Y(mS) in the X rest frame. If the pions are pi0s, then
171 : // choose the one where cos(theta) = [0:1].
172 :
173 0 : EvtVector4R P_YmS_X = boostTo(p->getP4(), P_X);
174 0 : double costheta = - p4[0].dot(P_YmS_X)/(p4[0].d3mag()*P_YmS_X.d3mag());
175 0 : if (EvtPDL::name(s1->getId()) == "pi0") {
176 0 : if (costheta < 0) {
177 0 : costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
178 0 : }
179 : }
180 0 : if (EvtPDL::name(s1->getId()) == "pi-") {
181 0 : costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
182 0 : }
183 :
184 : // // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "cos(theta) = " << costheta << endl;
185 :
186 :
187 :
188 : // Now boost the pion four vectors into the Y(mS) rest frame
189 0 : P_pi1 = boostTo(p4[0],P_YmS_X);
190 0 : P_pi2 = boostTo(p4[1],P_YmS_X);
191 :
192 : // Use a simple accept-reject to test this dipion system
193 :
194 : // Now compute the components of the matrix-element squared
195 : //
196 : // M(x,y)^2 = Q(x,y)^2 + |B/A|^2 * E1E2(x,y)^2 + 2*Re(B/A)*Q(x,y)*E1E2(x,y)
197 : //
198 : // x=m_pipi^2 and y = cos(theta), and where
199 : //
200 : // Q(x,y) = (x^2 + 2*m_pi^2)
201 : //
202 : // E1E2(x,y) = (1/4) * ( (E1 + E2)^2 - (E2 - E1)^2_max * cos(theta)^2 )
203 : //
204 : // and E1 + E2 = M_mS - M_nS and (E2 - E1)_max is the maximal difference
205 : // in the energy of the two pions allowed for a given mX value.
206 : //
207 :
208 0 : double Q = (mX*mX - 2.0 * m_pi * m_pi);
209 :
210 : double deltaEmax =
211 0 : - 2.0 *
212 0 : sqrt( P_nS.get(0)*P_nS.get(0) - M_nS*M_nS ) *
213 0 : sqrt( 0.25 - pow(m_pi/mX,2.0));
214 :
215 0 : double sumE = (M_mS*M_mS - M_nS*M_nS + mX*mX)/(2.0 * M_mS);
216 :
217 0 : double E1E2 = 0.25 * ( pow(sumE, 2.0) - pow( deltaEmax * costheta, 2.0) );
218 :
219 0 : double M2 = Q*Q + (pow(ReB_over_A,2.0) + pow(ImB_over_A,2.0)) * E1E2*E1E2 + 2.0 * ReB_over_A * Q * E1E2;
220 :
221 : // phase space factor
222 : //
223 : // this is given as d(PS) = C * p(*)_X * p(X)_{pi+} * d(cosTheta) * d(m_X)
224 : //
225 : // where C is a normalization constant, p(*)_X is the X momentum magnitude in the
226 : // Y(mS) rest frame, and p(X)_{pi+} is the pi+/pi0 momentum in the X rest frame
227 : //
228 :
229 : double dPS =
230 0 : sqrt( (M_mS*M_mS - pow(M_nS + mX,2.0)) * (M_mS*M_mS - pow(M_nS - mX,2.0)) ) * // p(*)_X
231 0 : sqrt(mX*mX - 4*m_pi*m_pi); // p(X)_{pi}
232 :
233 : // the double-differential decay rate dG/(dcostheta dmX)
234 0 : double dG = M2 * dPS;
235 :
236 : // Throw a uniform random number from 0 --> probMax and do accept/reject on this
237 :
238 0 : double rnd = EvtRandom::Flat(0.0,getProbMax(0.0));
239 :
240 0 : if (rnd < dG)
241 0 : acceptX = true;
242 :
243 0 : }
244 :
245 :
246 : // initialize the daughters
247 0 : v->init( getDaugs()[0], P_nS);
248 0 : s1->init( getDaugs()[1], P_pi1);
249 0 : s2->init( getDaugs()[2], P_pi2);
250 :
251 : // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "M_nS = " << v->getP4().mass() << endl;
252 : // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "m_pi = " << s1->getP4().mass() << endl;
253 : // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "m_pi = " << s2->getP4().mass() << endl;
254 : // report(Severity::Info,"EvtYmSToYnSpipiCLEO") << "M2 = " << M2 << endl;
255 :
256 : // Pass the polarization of the parent Upsilon
257 0 : EvtVector4C ep0,ep1,ep2;
258 :
259 0 : ep0=p->eps(0);
260 0 : ep1=p->eps(1);
261 0 : ep2=p->eps(2);
262 :
263 :
264 0 : vertex(0,0,(ep0*v->epsParent(0).conj()));
265 0 : vertex(0,1,(ep0*v->epsParent(1).conj()));
266 0 : vertex(0,2,(ep0*v->epsParent(2).conj()));
267 :
268 0 : vertex(1,0,(ep1*v->epsParent(0).conj()));
269 0 : vertex(1,1,(ep1*v->epsParent(1).conj()));
270 0 : vertex(1,2,(ep1*v->epsParent(2).conj()));
271 :
272 0 : vertex(2,0,(ep2*v->epsParent(0).conj()));
273 0 : vertex(2,1,(ep2*v->epsParent(1).conj()));
274 0 : vertex(2,2,(ep2*v->epsParent(2).conj()));
275 :
276 :
277 : return ;
278 :
279 0 : }
280 :
281 :
282 :
283 :
|
