Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtSSSCPpng.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: EvtSSSCPpng.cc
12 //
13 // Description: Routine to decay B -> 2 scalars taking into account penguin 
14 //              contributions (assuming single quark dominance for penguins)
15 //
16 // Modification history:
17 //
18 //    RYD/NK    December 3, 1997          Module created
19 //
20 //------------------------------------------------------------------------
21 // 
22 #include "EvtGenBase/EvtPatches.hh"
23 #include <stdlib.h>
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtRandom.hh"
26 #include "EvtGenBase/EvtGenKine.hh"
27 #include "EvtGenBase/EvtCPUtil.hh"
28 #include "EvtGenBase/EvtPDL.hh"
29 #include "EvtGenBase/EvtReport.hh"
30 #include "EvtGenModels/EvtSSSCPpng.hh"
31 #include "EvtGenBase/EvtId.hh"
32 #include <string>
33 #include "EvtGenBase/EvtConst.hh"
34
35 EvtSSSCPpng::~EvtSSSCPpng() {}
36
37 std::string EvtSSSCPpng::getName(){
38
39   return "SSS_CP_PNG";     
40
41 }
42
43
44 EvtDecayBase* EvtSSSCPpng::clone(){
45
46   return new EvtSSSCPpng;
47
48 }
49
50 void EvtSSSCPpng::init(){
51
52   // check that there are 7 arguments
53   checkNArg(7);
54   checkNDaug(2);
55
56   checkSpinParent(EvtSpinType::SCALAR);
57
58   checkSpinDaughter(0,EvtSpinType::SCALAR);
59   checkSpinDaughter(1,EvtSpinType::SCALAR);
60
61 }
62
63 void EvtSSSCPpng::initProbMax(){
64  
65    setProbMax(getArg(5)*getArg(5)*(1 + getArg(6)*getArg(6)));
66
67 }
68
69 void EvtSSSCPpng::decay( EvtParticle *p ){
70
71   //added by Lange Jan4,2000
72   static EvtId B0=EvtPDL::getId("B0");
73   static EvtId B0B=EvtPDL::getId("anti-B0");
74
75   double t;
76   EvtId other_b;
77  
78   p->initializePhaseSpace(getNDaug(),getDaugs());
79
80   EvtComplex amp;
81
82   EvtComplex A,Abar;
83   //EvtComplex ACC, AbarCC;
84
85 // assume single (top) quark dominance for the penguin.
86
87 // old: a0=alpha, a1=dm, a2=1, a3=1, a4=0, a5=1, a6=0
88 // new: a0=beta, a1=gamma, a2=delta, a3=dm, a4=1, a5=1=A_{T}, a6=A_{P}/A_{T}
89
90 // e.g., for B -> pi pi  
91 // A_{T} = |V_{ub} V_{ud}| T
92 // A_{P} = |V_{tb} V_{td}| P
93 // P and T are purely hadronic matrix elements       
94 // P/T = 0.055, A_{P}/A_{T} = 0.2 (see Marrocchesi and Paver, hep-ph/9702353)
95
96 // A = A_{T}( exp(i(beta+gamma)) + (A_{P}/A_{T}) exp(i(delta))
97 // A_bar = same, except for the sign of the weak phases 
98 // here, delta = delta_{p}-delta_{t} (rel. strong phase)
99
100    A=getArg(5)*(EvtComplex(cos(-getArg(0)-getArg(1)),sin(-getArg(0)-getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
101  
102    Abar=getArg(5)*(EvtComplex(cos(getArg(0)+getArg(1)),sin(getArg(0)+getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
103  
104 // get fraction of B0 tags with these amplitudes
105
106   //double xd = 0.65;
107   double ratio = 1/(1 + 0.65*0.65);
108   
109   EvtComplex rf, rbarf;
110
111   rf = EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*Abar/A;
112   rbarf = EvtComplex(1.0)/rf;
113
114   double A2 = real(A)*real(A) + imag(A)*imag(A);
115   double Abar2 = real(Abar)*real(Abar) + imag(Abar)*imag(Abar);
116   
117   double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);    
118   double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);    
119
120   //fraction of B0 _tags_
121   double fract =(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio)); 
122   
123   EvtCPUtil::getInstance()->OtherB(p,t,other_b,fract);
124
125 //this method works just as well -- NK
126 //randomly generate the tag (B0 or B0B) 
127
128 //  double tag = EvtRandom::Flat(0.0,1.0);
129 //  if (tag < 0.5) {
130 //
131 //   EvtCPUtil::OtherB(p,t,other_b,1.0);
132 //   other_b = B0;
133 //  }
134 //  else {
135 //   
136 //   EvtCPUtil::OtherB(p,t,other_b,0.0);
137 //   other_b = B0B;
138 //  }
139
140 //mixing angle = -beta
141
142    if (other_b==B0B){
143      amp=A*cos(getArg(3)*t/(2*EvtConst::c))+
144        EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
145        getArg(4)*EvtComplex(0.0,1.0)*Abar*sin(getArg(3)*t/(2*EvtConst::c));
146    }
147    if (other_b==B0){
148      amp=A*EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
149        EvtComplex(0.0,1.0)*sin(getArg(3)*t/(2*EvtConst::c))+       
150        getArg(4)*Abar*cos(getArg(3)*t/(2*EvtConst::c));
151    }
152
153    vertex(amp);
154
155   return ;
156 }
157
158