Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
4 : // This software is part of the EvtGen package. If you use all or part
5 : // of it, please give an appropriate acknowledgement.
6 : //
7 : // Copyright Information: See EvtGen/COPYRIGHT
8 : //
9 : // Module: EvtGenModels/EvtBcToNPi.hh
10 : //
11 : // Description: General decay model for Bc -> V + npi and Bc -> P + npi
12 : //
13 : // Modification history:
14 : //
15 : // A.Berezhnoy, A.Likhoded, A.Luchinsky April 2011 Module created
16 : //
17 : //------------------------------------------------------------------------
18 :
19 : #include "EvtGenBase/EvtPatches.hh"
20 :
21 : #include "EvtGenBase/EvtPDL.hh"
22 : #include "EvtGenBase/EvtReport.hh"
23 : #include "EvtGenBase/EvtVector4C.hh"
24 : #include "EvtGenBase/EvtTensor4C.hh"
25 : #include "EvtGenBase/EvtIdSet.hh"
26 : #include "EvtGenBase/EvtSpinType.hh"
27 : #include "EvtGenBase/EvtScalarParticle.hh"
28 :
29 : #include "EvtGenModels/EvtBcToNPi.hh"
30 :
31 : #include <iostream>
32 : using std::endl;
33 :
34 0 : EvtBcToNPi::EvtBcToNPi(bool printAuthorInfo) {
35 0 : nCall=0; maxAmp2=0;
36 0 : if (printAuthorInfo == true) {this->printAuthorInfo();}
37 0 : }
38 :
39 0 : EvtBcToNPi::~EvtBcToNPi() {
40 0 : }
41 :
42 : std::string EvtBcToNPi::getName(){
43 :
44 0 : return "EvtBcToNPi";
45 :
46 : }
47 :
48 : EvtDecayBase* EvtBcToNPi::clone(){
49 :
50 0 : return new EvtBcToNPi;
51 :
52 0 : }
53 :
54 :
55 : void EvtBcToNPi::init(){
56 :
57 : // check spins
58 0 : checkSpinParent(EvtSpinType::SCALAR);
59 :
60 :
61 : // the others are scalar
62 0 : for (int i=1; i<=(getNDaug()-1);i++) {
63 0 : checkSpinDaughter(i,EvtSpinType::SCALAR);
64 : };
65 :
66 0 : _beta=-0.108; _mRho=0.775; _gammaRho=0.149;
67 0 : _mRhopr=1.364; _gammaRhopr=0.400; _mA1=1.23; _gammaA1=0.4;
68 :
69 : // read arguments
70 0 : if( EvtPDL::getSpinType(getDaug(0)) == EvtSpinType::VECTOR) {
71 0 : checkNArg(10);
72 : int n=0;
73 0 : _maxProb=getArg(n++);
74 0 : FA0_N=getArg(n++);
75 0 : FA0_c1=getArg(n++);
76 0 : FA0_c2=getArg(n++);
77 0 : FAp_N=getArg(n++);
78 0 : FAp_c1=getArg(n++);
79 0 : FAp_c2=getArg(n++);
80 0 : FV_N=getArg(n++);
81 0 : FV_c1=getArg(n++);
82 0 : FV_c2=getArg(n++);
83 0 : FAm_N=0;
84 0 : FAm_c1=0;
85 0 : FAm_c2=0;
86 0 : }
87 0 : else if( EvtPDL::getSpinType(getDaug(0)) == EvtSpinType::SCALAR) {
88 0 : checkNArg(4);
89 : int n=0;
90 0 : _maxProb=getArg(n++);
91 0 : Fp_N=getArg(n++);
92 0 : Fp_c1=getArg(n++);
93 0 : Fp_c2=getArg(n++);
94 0 : Fm_N=0;
95 0 : Fm_c1=0;
96 0 : Fm_c2=0;
97 : }
98 : else {
99 0 : report(Severity::Error,"EvtGen") << "Have not yet implemented this final state in BCPSINPI model" << endl;
100 0 : report(Severity::Error,"EvtGen") << "Ndaug="<<getNDaug() << endl;
101 0 : for ( int id=0; id<(getNDaug()-1); id++ )
102 0 : report(Severity::Error,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
103 0 : return;
104 :
105 : };
106 :
107 0 : if(getNDaug()<2 || getNDaug()>4) {
108 0 : report(Severity::Error,"EvtGen") << "Have not yet implemented this final state in BCPSINPI model" << endl;
109 0 : report(Severity::Error,"EvtGen") << "Ndaug="<<getNDaug() << endl;
110 0 : for ( int id=0; id<(getNDaug()-1); id++ )
111 0 : report(Severity::Error,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
112 0 : return;
113 : }
114 :
115 0 : }
116 :
117 : double EvtBcToNPi::_ee(double M, double m1, double m2) {
118 0 : return (M*M+m1*m1-m2*m2)/(2*M);
119 : };
120 :
121 : double EvtBcToNPi::_pp(double M, double m1, double m2) {
122 0 : double __ee=_ee(M,m1,m2);
123 0 : return sqrt(__ee*__ee-m1*m1);
124 : };
125 :
126 :
127 : void EvtBcToNPi::initProbMax(){
128 0 : if(_maxProb>0.) setProbMax(_maxProb);
129 : else {
130 0 : EvtId id=getParentId();
131 0 : EvtScalarParticle *p=new EvtScalarParticle();
132 0 : p->init(id, EvtPDL::getMass(id),0., 0., 0.);
133 0 : p->setDiagonalSpinDensity();
134 : // add daughters
135 0 : p->makeDaughters(getNDaug(), getDaugs() );
136 :
137 : // fill the momenta
138 0 : if(getNDaug()==2) {
139 0 : double M=EvtPDL::getMass(id), m1=EvtPDL::getMass(getDaug(0)), m2=EvtPDL::getMass(getDaug(1));
140 0 : double __pp=_pp(M,m1,m2);
141 0 : p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,m2), 0., 0., __pp) );
142 0 : p->getDaug(1)->setP4( EvtVector4R( _ee(M,m2,m1), 0., 0., -__pp) );
143 0 : }
144 0 : else if( getNDaug()==3) {
145 0 : double M=EvtPDL::getMass(id),
146 0 : m1=EvtPDL::getMass(getDaug(0)),
147 0 : m2=EvtPDL::getMass(getDaug(1)),
148 0 : m3=EvtPDL::getMass(getDaug(2));
149 0 : double __ppRho=_pp(M,m1,_mRho), __ppPi=_pp(_mRho,m2,m3);
150 0 : p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,_mRho), 0., 0., __ppRho) );
151 0 : EvtVector4R _pRho( _ee(M,_mRho,m1), 0., 0., -__ppRho);
152 0 : EvtVector4R _p2( _ee(_mRho, m2, m3), 0., 0., __ppPi); _p2.applyBoostTo(_pRho);
153 0 : EvtVector4R _p3( _ee(_mRho, m2, m3), 0., 0., -__ppPi); _p3.applyBoostTo(_pRho);
154 0 : p->getDaug(1)->setP4(_p2);
155 0 : p->getDaug(2)->setP4(_p3);
156 :
157 0 : }
158 0 : else if( getNDaug()==4) {
159 0 : double M=EvtPDL::getMass(id),
160 0 : m1=EvtPDL::getMass(getDaug(0)),
161 0 : m2=EvtPDL::getMass(getDaug(1)),
162 0 : m3=EvtPDL::getMass(getDaug(2)),
163 0 : m4=EvtPDL::getMass(getDaug(3));
164 0 : if(M<m1+_mA1) return;
165 0 : double __ppA1=_pp(M,m1,_mA1),
166 0 : __ppRho=_pp(_mA1,_mRho,m4),
167 0 : __ppPi=_pp(_mRho, m2, m3);
168 0 : p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,_mRho), 0., 0., __ppA1) );
169 0 : EvtVector4R _pA1( _ee(M,_mA1,m1), 0., 0., -__ppA1);
170 0 : EvtVector4R _pRho(_ee(_mA1, _mRho, m4), 0, 0, __ppRho);
171 0 : _pRho.applyBoostTo(_pA1);
172 0 : EvtVector4R _p4( _ee(_mA1, m4, _mRho), 0, 0, -__ppRho); _p4.applyBoostTo(_pA1);
173 0 : p->getDaug(3)->setP4(_p4);
174 0 : EvtVector4R _p2( _ee(_mRho, m2, m3), 0, 0, __ppPi); _p2.applyBoostTo(_pRho);
175 0 : p->getDaug(1)->setP4(_p2);
176 0 : EvtVector4R _p3( _ee(_mRho, m2, m3), 0, 0, -__ppPi); _p2.applyBoostTo(_pRho);
177 0 : p->getDaug(2)->setP4(_p3);
178 0 : };
179 :
180 :
181 0 : _amp2.init(p->getId(),getNDaug(),getDaugs());
182 :
183 0 : decay(p);
184 :
185 0 : EvtSpinDensity rho=_amp2.getSpinDensity();
186 :
187 0 : double prob=p->getSpinDensityForward().normalizedProb(rho);
188 :
189 0 : if(prob>0) setProbMax(0.9*prob);
190 :
191 :
192 0 : };
193 0 : }
194 :
195 : void EvtBcToNPi::decay( EvtParticle *root_particle ){
196 0 : ++nCall;
197 :
198 0 : EvtIdSet thePis("pi+","pi-","pi0");
199 0 : EvtComplex I=EvtComplex(0.0, 1.0);
200 :
201 :
202 0 : root_particle->initializePhaseSpace(getNDaug(),getDaugs());
203 :
204 0 : EvtVector4R
205 0 : p(root_particle->mass(), 0., 0., 0.), // Bc momentum
206 0 : k=root_particle->getDaug(0)->getP4(), // J/psi momenta
207 0 : Q=p-k;
208 :
209 0 : double Q2=Q.mass2();
210 :
211 :
212 : // check pi-mesons and calculate hadronic current
213 0 : EvtVector4C hardCur;
214 : bool foundHadCurr=false;
215 :
216 0 : if ( getNDaug() == 2 ) // Bc -> psi pi+
217 : {
218 0 : hardCur=Q;
219 : foundHadCurr=true;
220 0 : }
221 0 : else if ( getNDaug() == 3 ) // Bc -> psi pi+ pi0
222 : {
223 0 : EvtVector4R p1,p2;
224 0 : p1=root_particle->getDaug(1)->getP4(), // pi+ momenta
225 0 : p2=root_particle->getDaug(2)->getP4(), // pi0 momentum
226 0 : hardCur=Fpi(p1,p2)*(p1-p2);
227 : foundHadCurr=true;
228 0 : }
229 0 : else if( getNDaug()==4 ) // Bc -> psi pi+ pi pi
230 : {
231 : int diffPi(0),samePi1(0),samePi2(0);
232 0 : if ( getDaug( 1) == getDaug( 2) ) {diffPi= 3; samePi1= 1; samePi2= 2;}
233 0 : if ( getDaug( 1) == getDaug( 3) ) {diffPi= 2; samePi1= 1; samePi2= 3;}
234 0 : if ( getDaug( 2) == getDaug( 3) ) {diffPi= 1; samePi1= 2; samePi2= 3;}
235 :
236 0 : EvtVector4R p1=root_particle->getDaug(samePi1)->getP4();
237 0 : EvtVector4R p2=root_particle->getDaug(samePi2)->getP4();
238 0 : EvtVector4R p3=root_particle->getDaug(diffPi)->getP4();
239 :
240 0 : EvtComplex BA1;
241 0 : double GA1=_gammaA1*pi3G(Q2,samePi1)/pi3G(_mA1*_mA1,samePi1);
242 0 : EvtComplex denBA1(_mA1*_mA1 - Q.mass2(),-1.*_mA1*GA1);
243 0 : BA1 = _mA1*_mA1 / denBA1;
244 :
245 0 : hardCur = BA1*( (p1-p3) - (Q*(Q*(p1-p3))/Q2)*Fpi(p2,p3) +
246 0 : (p2-p3) - (Q*(Q*(p2-p3))/Q2)*Fpi(p1,p3) );
247 : foundHadCurr=true;
248 0 : }
249 :
250 0 : if ( !foundHadCurr ) {
251 0 : report(Severity::Error,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
252 0 : report(Severity::Error,"EvtGen") << "Ndaug="<<getNDaug() << endl;
253 : int id;
254 0 : for ( id=0; id<(getNDaug()-1); id++ )
255 0 : report(Severity::Error,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
256 0 : ::abort();
257 : };
258 :
259 0 : EvtTensor4C H;
260 : double amp2=0.;
261 0 : if( root_particle->getDaug(0)->getSpinType() == EvtSpinType::VECTOR) {
262 0 : double FA0=FA0_N*exp(FA0_c1*Q2 + FA0_c2*Q2*Q2);
263 0 : double FAp=FAp_N*exp(FAp_c1*Q2 + FAp_c2*Q2*Q2);
264 0 : double FAm=FAm_N*exp(FAm_c1*Q2 + FAm_c2*Q2*Q2);
265 0 : double FV= FV_N* exp( FV_c1*Q2 + FV_c2*Q2*Q2 );
266 0 : H=-FA0*EvtTensor4C::g()
267 0 : -FAp*EvtGenFunctions::directProd(p,p+k)
268 0 : +FAm*EvtGenFunctions::directProd(p,p-k)
269 0 : +2*I*FV*dual(EvtGenFunctions::directProd(p,k));
270 0 : EvtVector4C Heps=H.cont2(hardCur);
271 :
272 0 : for(int i=0; i<4; i++) {
273 0 : EvtVector4C eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
274 0 : EvtComplex amp=eps*Heps;
275 0 : vertex(i,amp);
276 0 : amp2+=pow( abs(amp),2);
277 0 : }
278 0 : }
279 0 : else if( root_particle->getDaug(0)->getSpinType() == EvtSpinType::SCALAR) {
280 0 : double Fp=Fp_N*exp(Fp_c1*Q2 + Fp_c2*Q2*Q2);
281 0 : double Fm=Fm_N*exp(Fm_c1*Q2 + Fm_c2*Q2*Q2);
282 0 : EvtVector4C H=Fp*(p+k)+Fm*(p-k);
283 0 : EvtComplex amp=H*hardCur;
284 0 : vertex(amp);
285 0 : amp2+=pow( abs(amp),2);
286 0 : };
287 0 : if(amp2>maxAmp2) maxAmp2=amp2;
288 :
289 : return ;
290 0 : }
291 :
292 : EvtComplex EvtBcToNPi::Fpi( EvtVector4R q1, EvtVector4R q2) {
293 0 : double m1=q1.mass();
294 0 : double m2=q2.mass();
295 :
296 0 : EvtVector4R Q = q1 + q2;
297 0 : double mQ2= Q*Q;
298 :
299 : // momenta in the rho->pipi decay
300 0 : double dRho= _mRho*_mRho - m1*m1 - m2*m2;
301 0 : double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
302 :
303 0 : double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
304 0 : double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
305 :
306 0 : double dQ= mQ2 - m1*m1 - m2*m2;
307 0 : double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
308 :
309 :
310 0 : double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
311 0 : EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
312 0 : EvtComplex BRho= _mRho*_mRho / BRhoDem;
313 :
314 0 : double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
315 0 : EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
316 0 : EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
317 :
318 0 : return (BRho + _beta*BRhopr)/(1+_beta);
319 0 : }
320 :
321 : double EvtBcToNPi::pi3G(double m2,int dupD) {
322 0 : double mPi= EvtPDL::getMeanMass(getDaug(dupD));
323 0 : if ( m2 > (_mRho+mPi) ) {
324 0 : return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
325 : }
326 : else {
327 0 : double t1=m2-9.0*mPi*mPi;
328 0 : return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
329 : }
330 0 : }
331 :
332 : void EvtBcToNPi::printAuthorInfo() {
333 :
334 0 : report(Severity::Info,"EvtGen")<<"Defining EvtBcToNPi model: Bc -> V + npi and Bc -> P + npi decays\n"
335 0 : <<"from A.V. Berezhnoy, A.K. Likhoded, A.V. Luchinsky: "
336 0 : <<"Phys.Rev.D 82, 014012 (2010) and arXiV:1104.0808."<<endl;
337 :
338 0 : }
339 :
|
