Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtbTosllVectorAmp.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) 2000      Caltech, UCSB
10 //
11 // Module: EvtbTosllVectorAmp.cc
12 //
13 // Description: Routine to implement bTosll decays to vector
14 //              mesons. 
15 //
16 // Modification history:
17 //
18 //    Ryd       January 5,2000       Module created
19 //
20 //------------------------------------------------------------------------
21 //
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtPatches.hh"
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtGenKine.hh"
26 #include "EvtGenBase/EvtPDL.hh"
27 #include "EvtGenBase/EvtReport.hh"
28 #include "EvtGenBase/EvtVector4C.hh"
29 #include "EvtGenBase/EvtTensor4C.hh"
30 #include "EvtGenBase/EvtDiracSpinor.hh"
31 #include "EvtGenModels/EvtbTosllVectorAmp.hh"
32 #include "EvtGenBase/EvtId.hh"
33 #include "EvtGenBase/EvtIdSet.hh"
34 #include "EvtGenBase/EvtAmp.hh"
35 #include "EvtGenModels/EvtbTosllAmp.hh"
36 #include "EvtGenModels/EvtbTosllFF.hh"
37
38 void EvtbTosllVectorAmp::CalcAmp( EvtParticle *parent,
39                                   EvtAmp& amp,
40                                   EvtbTosllFF *formFactors ) {
41
42   //Add the lepton and neutrino 4 momenta to find q2
43
44   EvtVector4R q = parent->getDaug(1)->getP4() 
45                     + parent->getDaug(2)->getP4();
46   double q2 = (q.mass2());
47
48
49   double a1,a2,a0,v,t1,t2,t3;
50   double mesonmass = parent->getDaug(0)->mass();
51   double parentmass = parent->mass();
52
53   formFactors->getVectorFF(parent->getId(),
54                            parent->getDaug(0)->getId(),
55                            q2,
56                            mesonmass,
57                            a1,a2,a0,v,t1,t2,t3);
58
59
60   EvtId daught = parent->getDaug(0)->getId();
61   bool btod = false;
62   bool nnlo = true;
63  if 
64   (
65    daught == EvtPDL::getId(std::string("rho+")) ||
66    daught == EvtPDL::getId(std::string("rho-")) ||
67    daught == EvtPDL::getId(std::string("rho0")) ||
68    daught == EvtPDL::getId(std::string("omega"))
69   ) 
70    btod = true;
71
72   EvtVector4R p4b;
73   p4b.set(parent->mass(),0.0,0.0,0.0);
74   EvtVector4R p4meson = parent->getDaug(0)->getP4();
75  
76   EvtVector4C l11,l12;
77   EvtVector4C l21,l22;
78
79   EvtVector4C a11,a12;
80   EvtVector4C a21,a22;
81
82   EvtId parentID = parent->getId();
83
84       //EvtId l_num = parent->getDaug(1)->getId();
85
86   EvtVector4R pbhat=p4b/parentmass;
87   EvtVector4R qhat=q/parentmass;
88   EvtVector4R pkstarhat=p4meson/parentmass;
89   EvtVector4R phat=pbhat+pkstarhat;
90
91   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(q2,nnlo);
92   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(q2,nnlo,btod);
93   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(q2,nnlo);
94   EvtComplex uniti(0.0,1.0);
95
96   double mhatb=4.4/(parentmass); 
97   double mhatkstar=mesonmass/(parentmass);
98   double shat=q2/(parentmass*parentmass);
99
100
101   EvtComplex a;
102   a=c9eff*v*2/(1+mhatkstar)+4*mhatb*c7eff*t1/shat;
103   EvtComplex b;
104   b=(1+mhatkstar)*(c9eff*a1+2*mhatb*(1-mhatkstar)*c7eff*t2/shat);
105   EvtComplex c;
106   c=((1-mhatkstar)*c9eff*a2+
107             2*mhatb*c7eff*(t3+(1-mhatkstar*mhatkstar)*t2/shat))/
108     (1-mhatkstar*mhatkstar);
109   EvtComplex d;
110   d=(c9eff*((1+mhatkstar)*a1-(1-mhatkstar)*a2-2*mhatkstar*a0)
111             -2*mhatb*c7eff*t3)/shat;
112   EvtComplex e;
113   e=2*c10eff*v/(1+mhatkstar);
114   EvtComplex f;
115   f=(1+mhatkstar)*c10eff*a1;
116   EvtComplex g;
117   g=c10eff*a2/(1+mhatkstar);
118   EvtComplex h;
119   h=c10eff*((1+mhatkstar)*a1-(1-mhatkstar)*a2-2*mhatkstar*a0)/shat;
120   
121   EvtTensor4C T1,T2;
122   
123   static EvtIdSet bmesons("B-","anti-B0","anti-B_s0");
124   static EvtIdSet bbarmesons("B+","B0","B_s0");
125   
126   EvtParticle* lepPlus=0;
127   EvtParticle* lepMinus=0;
128   
129   int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
130   int charge2 = EvtPDL::chg3(parent->getDaug(2)->getId());
131   
132   lepPlus = (charge1 > charge2) ? parent->getDaug(1) : parent->getDaug(2);
133   lepMinus = (charge1 < charge2) ? parent->getDaug(1) : parent->getDaug(2);
134
135   if (bmesons.contains(parentID)) {
136
137     T1=a*dual(EvtGenFunctions::directProd(pbhat,pkstarhat))
138        -b*uniti*EvtTensor4C::g()
139        +c*uniti*EvtGenFunctions::directProd(pbhat,phat)
140        +d*uniti*EvtGenFunctions::directProd(pbhat,qhat);
141     
142     T2=e*dual(EvtGenFunctions::directProd(pbhat,pkstarhat))
143        -f*uniti*EvtTensor4C::g()
144        +g*uniti*EvtGenFunctions::directProd(pbhat,phat)
145        +h*uniti*EvtGenFunctions::directProd(pbhat,qhat);
146     
147     l11=EvtLeptonVCurrent(lepPlus->spParent(0),
148                           lepMinus->spParent(0));
149     l21=EvtLeptonVCurrent(lepPlus->spParent(1),
150                           lepMinus->spParent(0));
151     l12=EvtLeptonVCurrent(lepPlus->spParent(0), 
152                           lepMinus->spParent(1));
153     l22=EvtLeptonVCurrent(lepPlus->spParent(1),
154                           lepMinus->spParent(1));
155
156     a11=EvtLeptonACurrent(lepPlus->spParent(0),
157                           lepMinus->spParent(0));
158     a21=EvtLeptonACurrent(lepPlus->spParent(1),
159                           lepMinus->spParent(0));
160     a12=EvtLeptonACurrent(lepPlus->spParent(0),
161                           lepMinus->spParent(1));
162     a22=EvtLeptonACurrent(lepPlus->spParent(1),
163                           lepMinus->spParent(1));
164
165   } else {
166     
167     if (bbarmesons.contains(parentID)) {
168
169       T1=-a*dual(EvtGenFunctions::directProd(pbhat,pkstarhat))
170          -b*uniti*EvtTensor4C::g()
171          +c*uniti*EvtGenFunctions::directProd(pbhat,phat)
172          +d*uniti*EvtGenFunctions::directProd(pbhat,qhat);
173       
174       T2=-e*dual(EvtGenFunctions::directProd(pbhat,pkstarhat))
175          -f*uniti*EvtTensor4C::g()
176          +g*uniti*EvtGenFunctions::directProd(pbhat,phat)
177          +h*uniti*EvtGenFunctions::directProd(pbhat,qhat);
178       
179       l11=EvtLeptonVCurrent(lepPlus->spParent(1),
180                             lepMinus->spParent(1));
181       l21=EvtLeptonVCurrent(lepPlus->spParent(0),
182                             lepMinus->spParent(1));
183       l12=EvtLeptonVCurrent(lepPlus->spParent(1),
184                             lepMinus->spParent(0));
185       l22=EvtLeptonVCurrent(lepPlus->spParent(0),
186                             lepMinus->spParent(0));
187       
188       a11=EvtLeptonACurrent(lepPlus->spParent(1),
189                             lepMinus->spParent(1));
190       a21=EvtLeptonACurrent(lepPlus->spParent(0),
191                             lepMinus->spParent(1));
192       a12=EvtLeptonACurrent(lepPlus->spParent(1),
193                             lepMinus->spParent(0));
194       a22=EvtLeptonACurrent(lepPlus->spParent(0),
195                             lepMinus->spParent(0));
196       
197     }
198     else{
199       report(ERROR,"EvtGen") << "Wrong lepton number\n";
200       T1.zero(); T2.zero(); // Set all tensor terms to zero.
201     }    
202   }
203
204
205   int i;
206
207   for(i=0;i<3;i++){
208     EvtVector4C eps=parent->getDaug(0)->epsParent(i).conj();
209
210     EvtVector4C E1=T1.cont1(eps);
211     EvtVector4C E2=T2.cont1(eps);
212
213     amp.vertex(i,0,0,l11*E1+a11*E2);
214     amp.vertex(i,0,1,l12*E1+a12*E2);
215     amp.vertex(i,1,0,l21*E1+a21*E2);
216     amp.vertex(i,1,1,l22*E1+a22*E2);
217   } 
218 }
219
220
221
222
223
224
225