Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBcVNpi.cpp
CommitLineData
0ca57c2f 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
44EvtBcVNpi::~EvtBcVNpi() {
45// cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
46
47}
48
49std::string EvtBcVNpi::getName(){ return "BC_VNPI";}
50
51EvtDecayBase* EvtBcVNpi::clone() { return new EvtBcVNpi;}
52
53//======================================================
54void EvtBcVNpi::init(){
55 //cout<<"BcVNpi::init()"<<endl;
56
57 checkNArg(1);
58 checkSpinParent(EvtSpinType::SCALAR);
59 checkSpinDaughter(0,EvtSpinType::VECTOR);
60 for (int i=1; i<=(getNDaug()-1);i++) {
61 checkSpinDaughter(i,EvtSpinType::SCALAR);
62 };
63
64 if(getNDaug()<2 || getNDaug()>6) {
65 report(ERROR,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl;
66 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
67 for ( int id=0; id<(getNDaug()-1); id++ )
68 report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
69 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 idVector = getDaug(0).getId();
77 whichfit = int(getArg(0)+0.1);
78// cout<<"BcVNpi: whichfit ="<<whichfit<<" idVector="<<idVector<<endl;
79 ffmodel = new EvtBCVFF(idVector,whichfit);
80
81 wcurr = new EvtWnPi();
82
83 nCall = 0;
84}
85
86//======================================================
87void EvtBcVNpi::initProbMax() {
88// cout<<"BcVNpi::initProbMax()"<<endl;
89 if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==6) setProbMax(720000.);
90 else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==6) setProbMax(471817.);
91 else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==4) setProbMax(42000.);
92 else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==4) setProbMax(16000.);
93
94 else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==4) setProbMax(1200.);
95 else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==4) setProbMax(2600.);
96 else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==6) setProbMax(40000.);
97 else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==6) setProbMax(30000.);
98}
99
100//======================================================
101void EvtBcVNpi::decay( EvtParticle *root_particle ) {
102 ++nCall;
103// cout<<"BcVNpi::decay()"<<endl;
104 root_particle->initializePhaseSpace(getNDaug(),getDaugs());
105
106 EvtVector4R
107 p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum
108 p4meson=root_particle->getDaug(0)->getP4(), // J/psi momenta
109 Q=p4b-p4meson;
110 double Q2=Q.mass2();
111
112
113// check pi-mesons and calculate hadronic current
114 EvtVector4C hardCur;
115// bool foundHadCurr=false;
116 if( getNDaug() == 2) {
117 hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() );
118// foundHadCurr=true;
119 }
120 else if( getNDaug() == 3) {
121 hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
122 root_particle->getDaug(2)->getP4()
123 );
124// foundHadCurr=true;
125 }
126 else if( getNDaug() == 4) {
127 hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
128 root_particle->getDaug(2)->getP4(),
129 root_particle->getDaug(3)->getP4()
130 );
131// foundHadCurr=true;
132 }
133 else if( getNDaug() == 6) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162
134 {
135
136 hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(),
137 root_particle->getDaug(2)->getP4(),
138 root_particle->getDaug(3)->getP4(),
139 root_particle->getDaug(4)->getP4(),
140 root_particle->getDaug(5)->getP4()
141 );
142// foundHadCurr=true;
143 }
144 else {
145 report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
146 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
147 int id;
148 for ( id=0; id<(getNDaug()-1); id++ )
149 report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
150 ::abort();
151 };
152
153// calculate Bc -> V W form-factors
154 double a1f, a2f, vf, a0f;
155 double m_meson = root_particle->getDaug(0)->mass();
156 double m_b = root_particle->mass();
157 ffmodel->getvectorff(root_particle->getId(),
158 root_particle->getDaug(0)->getId(),
159 Q2,
160 m_meson,
161 &a1f,
162 &a2f,
163 &vf,
164 &a0f);
165 double a3f = ((m_b+m_meson)/(2.0*m_meson))*a1f -
166 ((m_b-m_meson)/(2.0*m_meson))*a2f;
167
168// calculate Bc -> V W current
169 EvtTensor4C H;
170 H = a1f*(m_b+m_meson)*EvtTensor4C::g();
171 H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson);
172 H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson));
173 H.addDirProd((a0f-a3f)*2.0*(m_meson/Q2)*p4b,p4b-p4meson);
174 EvtVector4C Heps=H.cont2(hardCur);
175
176 for(int i=0; i<4; i++) {
177 EvtVector4C eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
178 EvtComplex amp=eps*Heps;
179 vertex(i,amp);
180 };
181
182}
183