Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtResonance.cpp
CommitLineData
da0e9ce3 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: EvtResonance.cc
12//
13// Description: resonance-defining class
14//
15// Modification history:
16//
17// NK September 4, 1997 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtGenBase/EvtPatches.hh"
22#include <math.h>
23#include "EvtGenBase/EvtVector4R.hh"
24#include "EvtGenBase/EvtKine.hh"
25#include "EvtGenBase/EvtComplex.hh"
26#include "EvtGenBase/EvtResonance.hh"
27#include "EvtGenBase/EvtReport.hh"
28#include "EvtGenBase/EvtConst.hh"
29using std::endl;
30
31EvtResonance::~EvtResonance(){}
32
33
34EvtResonance& EvtResonance::operator = ( const EvtResonance &n)
35{
36 if ( &n == this ) return *this;
37 _p4_p = n._p4_p;
38 _p4_d1 = n._p4_d1;
39 _p4_d2 = n._p4_d2;
40 _ampl = n._ampl;
41 _theta = n._theta;
42 _gamma = n._gamma;
43 _spin = n._spin;
44 _bwm = n._bwm;
45 return *this;
46}
47
48
49EvtResonance::EvtResonance(const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
50 const EvtVector4R& p4_d2, double ampl,
51 double theta, double gamma, double bwm, int spin):
52 _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2),_ampl(ampl), _theta(theta),
53 _gamma(gamma), _bwm(bwm), _spin(spin) {}
54
55EvtComplex EvtResonance::resAmpl() {
56
57 double pi180inv = 1.0/EvtConst::radToDegrees;
58
59 EvtComplex ampl;
60 //EvtVector4R _p4_d3 = _p4_p-_p4_d1-_p4_d2;
61
62 //get cos of the angle between the daughters from their 4-momenta
63 //and the 4-momentum of the parent
64
65 //in general, EvtDecayAngle(parent, part1+part2, part1) gives the angle
66 //the missing particle (not listed in the arguments) makes
67 //with part2 in the rest frame of both
68 //listed particles (12)
69
70 //angle 3 makes with 2 in rest frame of 12 (CS3)
71 double cos_phi_0 = EvtDecayAngle(_p4_p, _p4_d1+_p4_d2, _p4_d1);
72 //angle 3 makes with 1 in 12 is, of course, -cos_phi_0
73
74 switch (_spin) {
75
76 case 0 :
77 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
78 sqrt(_gamma/EvtConst::twoPi)*
79 (1.0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
80 break;
81
82 case 1 :
83 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
84 sqrt(_gamma/EvtConst::twoPi)*
85 (cos_phi_0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
86 break;
87
88 case 2:
89 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
90 sqrt(_gamma/EvtConst::twoPi)*
91 ((1.5*cos_phi_0*cos_phi_0-0.5)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0, 0.5*_gamma))));
92 break;
93
94 case 3:
95 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
96 sqrt(_gamma/EvtConst::twoPi)*
97 ((2.5*cos_phi_0*cos_phi_0*cos_phi_0-1.5*cos_phi_0)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0, 0.5*_gamma))));
98 break;
99
100 default:
101 report(DEBUG,"EvtGen") << "EvtGen: wrong spin in EvtResonance" << endl;
102 ampl = EvtComplex(0.0);
103 break;
104
105 }
106
107 return ampl;
108}
109
110EvtComplex EvtResonance::relBrWig(int i) {
111
112//this function returns relativistic Breit-Wigner amplitude
113//for a given resonance (for P-wave decays of scalars only at the moment!)
114
115 EvtComplex BW;
116 EvtVector4R _p4_d3 = _p4_p-_p4_d1-_p4_d2;
117 EvtVector4R _p4_12 = _p4_d1 + _p4_d2;
118
119 double msq13 = (_p4_d1 + _p4_d3).mass2();
120 double msq23 = (_p4_d2 + _p4_d3).mass2();
121 double msqParent = _p4_p.mass2();
122 double msq1 = _p4_d1.mass2();
123 double msq2 = _p4_d2.mass2();
124 double msq3 = _p4_d3.mass2();
125
126 double M;
127
128 double p2 = sqrt((_p4_12.mass2() - (_p4_d1.mass() + _p4_d2.mass())*(_p4_d1.mass() + _p4_d2.mass()))*(_p4_12.mass2() - (_p4_d1.mass() - _p4_d2.mass())*(_p4_d1.mass() - _p4_d2.mass())))/(2.0*_p4_12.mass());
129
130 double p2R = sqrt((_bwm*_bwm - (_p4_d1.mass() + _p4_d2.mass())*(_p4_d1.mass() + _p4_d2.mass()))*(_bwm*_bwm - (_p4_d1.mass() - _p4_d2.mass())*(_p4_d1.mass() - _p4_d2.mass())))/(2.0*_bwm);
131
132 double gam, R;
133
134 if (i == 1) {
135
136 R = 2.0/(0.197);
137
138 }
139 else R = 5.0/(0.197);
140
141 gam = _gamma*(_bwm/_p4_12.mass())*(p2/p2R)*(p2/p2R)*(p2/p2R)*((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
142 M = (msq13 - msq23 - (msqParent - msq3)*(msq1 - msq2)/(_bwm*_bwm))*sqrt((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
143
144 BW = sqrt(_gamma)*M/((_bwm*_bwm - _p4_12.mass2()) - EvtComplex(0.0,1.0)*gam*_bwm);
145
146 return BW;
147
148}
149
150