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: EvtBcVNpi.cc
12 : //
13 : // Description: Module to implement Bc -> psi + (n pi) decays
14 : //
15 : // Modification history:
16 : //
17 : // AVL July 6, 2012 Module created
18 : //
19 : //------------------------------------------------------------------------
20 : //
21 : #include "EvtGenBase/EvtPatches.hh"
22 : #include <iostream>
23 : #include <iomanip>
24 : #include <fstream>
25 : #include <ctype.h>
26 : #include <stdlib.h>
27 : #include <string.h>
28 : #include "EvtGenBase/EvtParticle.hh"
29 : #include "EvtGenBase/EvtPDL.hh"
30 : #include "EvtGenBase/EvtGenKine.hh"
31 : #include "EvtGenModels/EvtTauHadnu.hh"
32 : #include "EvtGenBase/EvtDiracSpinor.hh"
33 : #include "EvtGenBase/EvtReport.hh"
34 : #include "EvtGenBase/EvtVector4C.hh"
35 : #include "EvtGenBase/EvtTensor4C.hh"
36 : #include "EvtGenBase/EvtIdSet.hh"
37 : #include "EvtGenBase/EvtParser.hh"
38 :
39 : #include "EvtGenModels/EvtBcVNpi.hh"
40 : #include "EvtGenModels/EvtWnPi.hh"
41 :
42 :
43 :
44 0 : EvtBcVNpi::~EvtBcVNpi() {
45 : // cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
46 :
47 0 : }
48 :
49 0 : std::string EvtBcVNpi::getName(){ return "BC_VNPI";}
50 :
51 0 : EvtDecayBase* EvtBcVNpi::clone() { return new EvtBcVNpi;}
52 :
53 : //======================================================
54 : void EvtBcVNpi::init(){
55 : //cout<<"BcVNpi::init()"<<endl;
56 :
57 0 : checkNArg(1);
58 0 : checkSpinParent(EvtSpinType::SCALAR);
59 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
60 0 : for (int i=1; i<=(getNDaug()-1);i++) {
61 0 : checkSpinDaughter(i,EvtSpinType::SCALAR);
62 : };
63 :
64 0 : if(getNDaug()<2 || getNDaug()>6) {
65 0 : report(Severity::Error,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl;
66 0 : report(Severity::Error,"EvtGen") << "Ndaug="<<getNDaug() << endl;
67 0 : for ( int id=0; id<(getNDaug()-1); id++ )
68 0 : report(Severity::Error,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
69 0 : return;
70 : }
71 :
72 :
73 : // for(int i=0; i<getNDaug(); i++)
74 : // cout<<"BcVNpi::init \t\t daughter "<<i<<" : "<<getDaug(i).getId()<<" "<<EvtPDL::name(getDaug(i)).c_str()<<endl;
75 :
76 0 : idVector = getDaug(0).getId();
77 0 : whichfit = int(getArg(0)+0.1);
78 : // cout<<"BcVNpi: whichfit ="<<whichfit<<" idVector="<<idVector<<endl;
79 0 : ffmodel = new EvtBCVFF(idVector,whichfit);
80 :
81 0 : wcurr = new EvtWnPi();
82 :
83 0 : nCall = 0;
84 0 : }
85 :
86 : //======================================================
87 : void EvtBcVNpi::initProbMax() {
88 : // cout<<"BcVNpi::initProbMax()"<<endl;
89 0 : if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==6) setProbMax(720000.);
90 0 : else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==6) setProbMax(471817.);
91 0 : else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==4) setProbMax(42000.);
92 0 : else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==4) setProbMax(16000.);
93 :
94 0 : else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==4) setProbMax(1200.);
95 0 : else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==4) setProbMax(2600.);
96 0 : else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==6) setProbMax(40000.);
97 0 : else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==6) setProbMax(30000.);
98 0 : }
99 :
100 : //======================================================
101 : void EvtBcVNpi::decay( EvtParticle *root_particle ) {
102 0 : ++nCall;
103 : // cout<<"BcVNpi::decay()"<<endl;
104 0 : root_particle->initializePhaseSpace(getNDaug(),getDaugs());
105 :
106 0 : EvtVector4R
107 0 : p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum
108 0 : p4meson=root_particle->getDaug(0)->getP4(), // J/psi momenta
109 0 : Q=p4b-p4meson;
110 0 : double Q2=Q.mass2();
111 :
112 :
113 : // check pi-mesons and calculate hadronic current
114 0 : EvtVector4C hardCur;
115 : // bool foundHadCurr=false;
116 0 : if( getNDaug() == 2) {
117 0 : hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() );
118 : // foundHadCurr=true;
119 0 : }
120 0 : else if( getNDaug() == 3) {
121 0 : hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
122 0 : root_particle->getDaug(2)->getP4()
123 : );
124 : // foundHadCurr=true;
125 0 : }
126 0 : else if( getNDaug() == 4) {
127 0 : hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
128 0 : root_particle->getDaug(2)->getP4(),
129 0 : root_particle->getDaug(3)->getP4()
130 : );
131 : // foundHadCurr=true;
132 0 : }
133 0 : else if( getNDaug() == 6) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162
134 : {
135 :
136 0 : hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(),
137 0 : root_particle->getDaug(2)->getP4(),
138 0 : root_particle->getDaug(3)->getP4(),
139 0 : root_particle->getDaug(4)->getP4(),
140 0 : root_particle->getDaug(5)->getP4()
141 : );
142 : // foundHadCurr=true;
143 : }
144 : else {
145 0 : report(Severity::Error,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
146 0 : report(Severity::Error,"EvtGen") << "Ndaug="<<getNDaug() << endl;
147 : int id;
148 0 : for ( id=0; id<(getNDaug()-1); id++ )
149 0 : report(Severity::Error,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
150 0 : ::abort();
151 : };
152 :
153 : // calculate Bc -> V W form-factors
154 0 : double a1f, a2f, vf, a0f;
155 0 : double m_meson = root_particle->getDaug(0)->mass();
156 0 : double m_b = root_particle->mass();
157 0 : ffmodel->getvectorff(root_particle->getId(),
158 0 : root_particle->getDaug(0)->getId(),
159 : Q2,
160 : m_meson,
161 : &a1f,
162 : &a2f,
163 : &vf,
164 : &a0f);
165 0 : double a3f = ((m_b+m_meson)/(2.0*m_meson))*a1f -
166 0 : ((m_b-m_meson)/(2.0*m_meson))*a2f;
167 :
168 : // calculate Bc -> V W current
169 0 : EvtTensor4C H;
170 0 : H = a1f*(m_b+m_meson)*EvtTensor4C::g();
171 0 : H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson);
172 0 : H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson));
173 0 : H.addDirProd((a0f-a3f)*2.0*(m_meson/Q2)*p4b,p4b-p4meson);
174 0 : EvtVector4C Heps=H.cont2(hardCur);
175 :
176 0 : for(int i=0; i<4; i++) {
177 0 : EvtVector4C eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
178 0 : EvtComplex amp=eps*Heps;
179 0 : vertex(i,amp);
180 0 : };
181 :
182 0 : }
183 :
|
