Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtPto3PAmp.cpp
1 #include "EvtGenBase/EvtPatches.hh"
2 /*******************************************************************************
3  * Project: BaBar detector at the SLAC PEP-II B-factory
4  * Package: EvtGenBase
5  *    File: $Id: EvtPto3PAmp.cpp,v 1.3 2009-03-16 15:44:04 robbep Exp $
6  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
7  *
8  * Copyright (C) 2002 Caltech
9  *
10  * Modified: May 30, 2005, Denis Dujmic (ddujmic@slac.stanford.edu): flip sign 
11  *           for vector resonances (-i -> i at resonance mass).
12  *           Introduce Flatte lineshape.
13  *******************************************************************************/
14
15 #include <assert.h>
16 #include <math.h>
17 #include <iostream>
18 #include "EvtGenBase/EvtComplex.hh"
19 #include "EvtGenBase/EvtPto3PAmp.hh"
20 #include "EvtGenBase/EvtDalitzCoord.hh"
21 #include "EvtGenBase/EvtdFunction.hh"
22 #include "EvtGenBase/EvtCyclic3.hh"
23 #include "EvtGenBase/EvtReport.hh"
24
25 using std::endl;
26 using EvtCyclic3::Index;
27 using EvtCyclic3::Pair;
28
29
30 EvtPto3PAmp::EvtPto3PAmp(EvtDalitzPlot dp, Pair pairAng, Pair pairRes, 
31                          EvtSpinType::spintype spin, 
32                          const EvtPropagator& prop, NumType typeN) 
33   : EvtAmplitude<EvtDalitzPoint>(),
34   _pairAng(pairAng), _pairRes(pairRes),
35   _spin(spin),
36   _typeN(typeN),
37   _prop((EvtPropagator*) prop.clone()), 
38   _g0(prop.g0()),
39   _min(0),
40   _max(0),
41   _vb(prop.m0(),dp.m(EvtCyclic3::other(pairRes)),dp.bigM(),spin), 
42   _vd(dp.m(EvtCyclic3::first(pairRes)),dp.m(EvtCyclic3::second(pairRes)),prop.m0(),spin)
43 {}
44
45
46
47 EvtPto3PAmp::EvtPto3PAmp(const EvtPto3PAmp& other) 
48   : EvtAmplitude<EvtDalitzPoint>(other),
49   _pairAng(other._pairAng),
50   _pairRes(other._pairRes),
51   _spin(other._spin),
52   _typeN(other._typeN),
53   _prop( (other._prop) ? (EvtPropagator*) other._prop->clone() : 0),
54   _g0(other._g0),
55   _min(other._min),
56   _max(other._max),
57   _vb(other._vb), _vd(other._vd)
58 {}
59
60
61 EvtPto3PAmp::~EvtPto3PAmp()
62 {
63   if(_prop) delete _prop;
64 }
65
66
67 void EvtPto3PAmp::set_fd(double R)
68 {
69   _vd.set_f(R);
70 }
71
72 void EvtPto3PAmp::set_fb(double R)
73 {
74   _vb.set_f(R);
75 }
76
77  
78 EvtComplex EvtPto3PAmp::amplitude(const EvtDalitzPoint& x) const
79 {
80   EvtComplex amp(1.0,0.0);
81
82   double m = sqrt(x.q(_pairRes));
83
84   if ((_max>0 && m > _max) || (_min>0 && m < _min))  return EvtComplex(0.0,0.0);
85
86   EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
87                     x.m(EvtCyclic3::second(_pairRes)),m);
88   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
89
90
91   // Compute mass-dependent width for relativistic propagators
92
93   if(_typeN!=NBW && _typeN!=FLATTE) {
94     
95     _prop->set_g0(_g0*_vd.widthFactor(vd));
96   }
97   
98
99   // Compute propagator
100
101   amp *= evalPropagator(m);
102
103
104
105   // Compute form-factors
106
107   amp *= _vd.formFactor(vd);
108   amp *= _vb.formFactor(vb);
109
110   amp *= numerator(x);
111
112   return amp;
113 }
114
115
116 EvtComplex EvtPto3PAmp::numerator(const EvtDalitzPoint& x) const
117 {
118   EvtComplex ret(0.,0.);
119   double m = sqrt(x.q(_pairRes)); 
120   EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
121                     x.m(EvtCyclic3::second(_pairRes)),m);
122   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
123
124   // Non-relativistic Breit-Wigner
125
126   if(NBW == _typeN) {
127
128     ret = angDep(x);
129   }
130
131   // Standard relativistic Zemach propagator
132
133   else if(RBW_ZEMACH == _typeN) {
134
135     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
136   }
137     
138   // Kuehn-Santamaria normalization:
139   
140   else if(RBW_KUEHN == _typeN) {
141
142     ret = _prop->m0()*_prop->m0() * angDep(x);
143   }
144
145
146   // CLEO amplitude is not factorizable
147   //
148   // The CLEO amplitude numerator is proportional to:
149   //
150   // m2_AC - m2_BC + (m2_D - m2_C)(m2_B - m2_A)/m2_0
151   //
152   // m2_AC = (eA + eC)^2 + (P - P_C cosTh(BC))^2
153   // m2_BC = (eB + eC)^2 + (P + P_C cosTh(BC))^2
154   //
155   // The first term m2_AB-m2_BC is therefore a p-wave term
156   // - 4PP_C cosTh(BC) 
157   // The second term is an s-wave, the amplitude
158   // does not factorize!
159   //
160   // The first term is just Zemach. However, the sign is flipped! 
161   // Let's consistently use the convention in which the amplitude
162   // is proportional to +cosTh(BC). In the CLEO expressions, I will
163   // therefore exchange AB to get rid of the sign flip.
164   
165
166   if(RBW_CLEO == _typeN || FLATTE == _typeN || GS == _typeN) {
167
168     Index iA = EvtCyclic3::other(_pairAng);           // A = other(BC)
169     Index iB = EvtCyclic3::common(_pairRes,_pairAng); // B = common(AB,BC)
170     Index iC = EvtCyclic3::other(_pairRes);           // C = other(AB)
171     
172     double M = x.bigM();
173     double mA = x.m(iA);
174     double mB = x.m(iB);
175     double mC = x.m(iC);
176     double qAB = x.q(EvtCyclic3::combine(iA,iB));
177     double qBC = x.q(EvtCyclic3::combine(iB,iC));
178     double qCA = x.q(EvtCyclic3::combine(iC,iA));
179     
180     //double m0 = _prop->m0();
181
182     if(_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
183     else
184       if(_spin == EvtSpinType::VECTOR) {
185
186         //ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
187         ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
188       }
189       else
190         if(_spin == EvtSpinType::TENSOR) {
191       
192           //double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;       
193           double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;       
194           double x2 = M*M - mC*mC;      
195           //double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/m0/m0;      
196           double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/qAB;      
197           double x4 = mB*mB - mA*mA;
198           //double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/m0/m0;
199           double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/qAB;
200           ret = (x1*x1 - 1./3.*x3*x5);
201         }
202         else assert(0);
203   }
204
205   return ret;
206 }
207
208
209 double EvtPto3PAmp::angDep(const EvtDalitzPoint& x) const 
210
211   // Angular dependece for factorizable amplitudes  
212   // unphysical cosines indicate we are in big trouble
213
214   double cosTh = x.cosTh(_pairAng,_pairRes);  
215   if(fabs(cosTh) > 1.) {
216     
217     report(INFO,"EvtGen") << "cosTh " << cosTh << endl; 
218     assert(0);
219   }
220   
221   // in units of half-spin
222   
223   return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
224 }
225