Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtBcVNpi.cpp
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 EvtBcVNpi::~EvtBcVNpi() {
45 //   cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
46   
47 }
48
49 std::string EvtBcVNpi::getName(){ return "BC_VNPI";}
50
51 EvtDecayBase* EvtBcVNpi::clone() {  return new EvtBcVNpi;}
52
53 //======================================================
54 void 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 //======================================================
87 void 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 //======================================================
101 void 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