Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtBsquark.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: EvtTauScalarnu.cc
12 //
13 // Description: The leptonic decay of the tau meson.
14 //              E.g., tau- -> e- nueb nut
15 //
16 // Modification history:
17 //
18 //    RYD       January 17, 1997       Module created
19 //
20 //------------------------------------------------------------------------
21 //
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtPatches.hh"
24 #include <iostream>
25 #include <string>
26 #include "EvtGenBase/EvtParticle.hh"
27 #include "EvtGenBase/EvtDiracParticle.hh"
28 #include "EvtGenBase/EvtPDL.hh"
29 #include "EvtGenBase/EvtIdSet.hh"
30 #include "EvtGenBase/EvtGenKine.hh"
31 #include "EvtGenModels/EvtBsquark.hh"
32 #include "EvtGenBase/EvtDiracSpinor.hh"
33 #include "EvtGenBase/EvtGammaMatrix.hh"
34 #include "EvtGenBase/EvtReport.hh"
35
36 EvtBsquark::~EvtBsquark() {}
37
38 std::string EvtBsquark::getName(){
39
40   return "BSQUARK";     
41
42 }
43
44
45 EvtDecayBase* EvtBsquark::clone(){
46
47   return new EvtBsquark;
48
49 }
50
51 void EvtBsquark::init(){
52
53   // check that there are 5 arguments
54   checkNArg(5);
55 }
56
57 void EvtBsquark::initProbMax(){
58
59   //For now do not set a maximum.
60
61   //SetProbMax(0.000000000005);
62
63 }
64
65 void EvtBsquark::decay(EvtParticle *p){
66
67   static EvtId cquark=EvtPDL::getId("c");
68   static EvtId anticquark=EvtPDL::getId("anti-c");
69
70   static EvtIdSet leptons("e-","mu-","tau-");
71
72   p->initializePhaseSpace(getNDaug(),getDaugs());
73
74
75   int charge=1;
76
77   EvtParticle* lepton;
78   lepton = p->getDaug(1);
79   if (leptons.contains(lepton->getId())){
80     charge=-1;
81   }
82
83   
84   EvtDiracParticle charmquark;
85
86   //this is a very crude approximation...
87   if (charge==-1){
88     charmquark.init(cquark,p->getDaug(0)->getP4());
89   }
90   else{
91     charmquark.init(anticquark,p->getDaug(0)->getP4());
92   }
93     
94   EvtVector4R p4c = p->getDaug(0)->getP4();
95
96   EvtVector4R p4sn = p->getDaug(2)->getP4();
97
98   EvtVector4R p4b(p->mass(),0.0,0.0,0.0);
99
100   EvtComplex M[2][2];
101
102   int il,ic;
103
104   //project out the right handed current
105   EvtGammaMatrix PR=0.5*(EvtGammaMatrix::id()+EvtGammaMatrix::g5());
106
107   double tanbeta=getArg(1);
108   double cosbeta=cos(atan(tanbeta));
109   double sinbeta=sin(atan(tanbeta));
110
111   double mb=4.9;
112   double mc=1.3;
113   double mw=80.4;
114
115   double Mass=getArg(2);
116   double mu=getArg(3);
117   double mchargino=getArg(4);
118
119
120   double tan2phim=2*sqrt(2.0)*mw*(mu*cosbeta+Mass*sinbeta)/
121     (Mass*Mass-mu*mu+2*mw*mw*cos(2*atan(tanbeta)));
122
123   double phim=0.5*atan(tan2phim);
124   
125   EvtComplex U11=cos(phim);
126   EvtComplex U12=sin(phim);
127   EvtComplex U21=-sin(phim);
128   EvtComplex U22=cos(phim);
129
130   double tan2phip=2*sqrt(2.0)*mw*(mu*cosbeta+Mass*sinbeta)/
131     (Mass*Mass-mu*mu-2*mw*mw*cos(2*atan(tanbeta)));
132
133   double phip=0.5*atan(tan2phip);
134
135   EvtComplex V11=cos(phip);
136   EvtComplex V12=sin(phip);
137   EvtComplex V21=-sin(phip);
138   EvtComplex V22=cos(phip);
139
140
141   double theta=getArg(0);
142   double ctheta=cos(theta);
143   double stheta=sin(theta);
144
145   double vcsb=0.08;
146   double mchi1=mchargino;  
147   double mchi2=mchargino;
148
149   //overall scale factor
150   double g=1.0;
151
152   EvtComplex a1=mchi1*(U11*ctheta-mb*U12*stheta/(sqrt(2.0)*mw*cosbeta));
153   EvtComplex a2=mchi2*(U21*ctheta-mb*U22*stheta/(sqrt(2.0)*mw*cosbeta));
154   
155   EvtComplex b1=mc*conj(V12)*ctheta/(sqrt(2.0)*mw*sinbeta);
156   EvtComplex b2=mc*conj(V22)*ctheta/(sqrt(2.0)*mw*sinbeta);
157
158   EvtComplex f1=-(g*g*V11*vcsb)/((p4b-p4c).mass2()-mchi1*mchi1);
159   EvtComplex f2=-(g*g*V21*vcsb)/((p4b-p4c).mass2()-mchi1*mchi2);
160
161   //report(INFO,"EvtGen") <<g<<" "<<V11<<" "<<FL<<" "<<vcsb<<" "<<mchi1<<endl;
162   //report(INFO,"EvtGen") << "f1:"<<f1<<" "<<(p4b-p4c).mass2()<<endl;
163   //report(INFO,"EvtGen") << "f2:"<<f2<<" "<<(p4b-p4c).mass2()<<endl;
164
165   //report(INFO,"EvtGen") << "p4sn:"<<p4sn<<endl;
166
167   EvtGammaMatrix pslash=p4sn.get(0)*EvtGammaMatrix::g0()
168                        -p4sn.get(1)*EvtGammaMatrix::g1()
169                        -p4sn.get(2)*EvtGammaMatrix::g2()
170                        -p4sn.get(3)*EvtGammaMatrix::g3();
171
172   //report(INFO,"EvtGen") << "pslash:"<<pslash<<endl;
173
174
175
176   for(il=0;il<2;il++){
177     for(ic=0;ic<2;ic++){
178
179       EvtComplex a=0.0;
180       EvtComplex b=0.0;
181
182       if (charge==-1){
183         a=charmquark.spParent(ic)*(PR*lepton->spParent(il));
184         b=charmquark.spParent(ic)*((pslash*PR)*lepton->spParent(il));
185       }
186       else{
187         a=lepton->spParent(il)*(PR*charmquark.spParent(ic));
188         b=lepton->spParent(il)*((pslash*PR)*charmquark.spParent(ic));
189       }
190
191       //report(INFO,"EvtGen") <<"pslash*PR:"<<pslash*PR<<endl;
192       //report(INFO,"EvtGen") <<"sp charm:"<<charmquark.spParent(ic)<<endl;
193       //report(INFO,"EvtGen") <<"sp lepton:"<<lepton->spParent(il)<<endl;
194
195       M[ic][il]=f1*(a1*a+b1*b)+f2*(a2*a+b2*b);
196
197       //report(INFO,"EvtGen") << "Contr1:"<<a1<<" "<<a<<" "<<b1<<" "<<b<<endl;
198       //report(INFO,"EvtGen") << "Contr2:"<<a2<<" "<<a<<" "<<b2<<" "<<b<<endl;
199
200       //report(INFO,"EvtGen") <<"case1:"<<f1<<" "<<a1<<" "<<b1<<" "<<a<<" "<<b<<endl;
201       //report(INFO,"EvtGen") <<"case2:"<<f2<<" "<<a2<<" "<<b2<<" "<<a<<" "<<b<<endl;
202
203     }
204   }
205  
206   double prob=real(M[0][0]*conj(M[0][0])+
207                    M[1][0]*conj(M[1][0])+
208                    M[0][1]*conj(M[0][1])+
209                    M[1][1]*conj(M[1][1]));
210
211   //report(INFO,"EvtGen") <<"prob:"<<prob<<endl;
212
213   setProb(prob);
214
215   return;
216
217 }
218