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