BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
D0ToKSpipi.cxx
Go to the documentation of this file.
1#include "QCMCFilterAlg/D0ToKSpipi.h"
2#include <stdlib.h>
3#include <iostream>
4#include <string>
5#include <complex>
6#include <vector>
7#include <math.h>
8#include "TMath.h"
9
10#include "CLHEP/Random/RandFlat.h"
11#include "CLHEP/Matrix/Vector.h"
12#include "CLHEP/Matrix/Matrix.h"
13#include "CLHEP/Matrix/SymMatrix.h"
14#include "CLHEP/Vector/ThreeVector.h"
15#include "CLHEP/Vector/LorentzVector.h"
16#include "CLHEP/Vector/TwoVector.h"
17using CLHEP::HepVector;
18using CLHEP::Hep3Vector;
19using CLHEP::Hep2Vector;
20using CLHEP::HepLorentzVector;
21using namespace std;
22
23D0ToKSpipi::~D0ToKSpipi() {}
24
25void D0ToKSpipi::init(){
26 //std::cout << "D0ToKSpipi ==> Initialization !" << std::endl;
27
28 _nd = 3;
29 tan2thetaC = (0.22650*0.22650)/(1.-(0.22650*0.22650)) ; //sin(theta_C) = 0.22650 +/- 0.00048
30 pi180inv = 1.0*3.1415926/180;
31
32 mass_R[0]= 0.77526; width_R[0]= 0.14740; spin_R[0]= 1; ar[0]= 1; phir[0]= 0;
33 mass_R[1]= 0.78266; width_R[1]= 0.00868; spin_R[1]= 1; ar[1]= 0.037606; phir[1]= 109.677;
34 mass_R[2]= 1.27550; width_R[2]= 0.18670; spin_R[2]= 2; ar[2]= 1.54909; phir[2]= -42.7425;
35 mass_R[3]= 1.46500; width_R[3]= 0.40000; spin_R[3]= 1; ar[3]= 3.70735; phir[3]= 103.644;
36 mass_R[4]= 0.89167; width_R[4]= 0.0514; spin_R[4]= 1; ar[4]= 1.86093; phir[4]= 136.529;
37 mass_R[5]= 1.42730; width_R[5]= 0.10000; spin_R[5]= 2; ar[5]= 1.74288; phir[5]= -48.0968;
38 mass_R[6]= 1.71800; width_R[6]= 0.3220; spin_R[6]= 1; ar[6]= 3.31; phir[6]= -118.2;
39 mass_R[7]= 1.41400; width_R[7]= 0.2320; spin_R[7]= 1; ar[7]= 0.171672; phir[7]= -68.41;
40 mass_R[8]= 0.89167; width_R[8]= 0.0514; spin_R[8]= 1; ar[8]= 0.164; phir[8]= -42.2;
41 mass_R[9]= 1.42730; width_R[9]= 0.1000; spin_R[9]= 2; ar[9]= 0.1; phir[9]= -89.6;
42 mass_R[10]= 1.41400; width_R[10]= 0.2320; spin_R[10]= 1; ar[10]= 0.21; phir[10]= 150.2;
43 mass_R[11]= 1.42500; width_R[11]= 0.2700; spin_R[11]= 1; ar[11]= 2.78276; phir[11]= 97.9608;
44 mass_R[12]= 1.42500; width_R[12]= 0.2700; spin_R[12]= 1; ar[12]= 0.11; phir[12]= 162.3;
45
46 beta[0] = complex<double>( 0.255303*cos( 47.8861 *pi180inv), 0.255303*sin( 47.8861 *pi180inv));
47 beta[1] = complex<double>(13.4446 *cos( -5.11127*pi180inv), 13.4446 *sin( -5.11127*pi180inv));
48 beta[2] = complex<double>(38.8496 *cos(-30.06 *pi180inv), 38.8496 *sin(-30.06 *pi180inv));
49 beta[3] = complex<double>(13.1086 *cos(-81.4148 *pi180inv), 13.1086 *sin(-81.4148 *pi180inv));
50 beta[4] = complex<double>(0., 0.);
51
52 fprod[0] = complex<double>(5.08049*cos(-182.312*pi180inv), 5.08049*sin(-182.312*pi180inv));
53 fprod[1] = complex<double>(17.2388*cos(-219.209*pi180inv), 17.2388*sin(-219.209*pi180inv));
54 fprod[2] = complex<double>(19.0145*cos(-76.9884*pi180inv), 19.0145*sin(-76.9884*pi180inv));
55 fprod[3] = complex<double>(11.9875*cos(-190.502*pi180inv), 11.9875*sin(-190.502*pi180inv));
56 fprod[4] = complex<double>(0., 0.);
57 //beta.push_back( complex<double>( 0.255303*cos( 47.8861 *pi180inv), 0.255303*sin( 47.8861 *pi180inv)) );
58 //beta.push_back( complex<double>(13.4446 *cos( -5.11127*pi180inv), 13.4446 *sin( -5.11127*pi180inv)) );
59 //beta.push_back( complex<double>(38.8496 *cos(-30.06 *pi180inv), 38.8496 *sin(-30.06 *pi180inv)) );
60 //beta.push_back( complex<double>(13.1086 *cos(-81.4148 *pi180inv), 13.1086 *sin(-81.4148 *pi180inv)) );
61 //beta.push_back( complex<double>(0., 0.) );
62
63 //fprod.push_back( complex<double>(5.08049*cos(-182.312*pi180inv), 5.08049*sin(-182.312*pi180inv)));
64 //fprod.push_back( complex<double>(17.2388*cos(-219.209*pi180inv), 17.2388*sin(-219.209*pi180inv)));
65 //fprod.push_back( complex<double>(19.0145*cos(-76.9884*pi180inv), 19.0145*sin(-76.9884*pi180inv)));
66 //fprod.push_back( complex<double>(11.9875*cos(-190.502*pi180inv), 11.9875*sin(-190.502*pi180inv)));
67 //fprod.push_back( complex<double>(0., 0.));
68
69 ma[0]= 0.651; g[0][0]= 0.22889; g[0][1]= -0.55377; g[0][2]= 0; g[0][3]= -0.39899; g[0][4]= -0.34639;
70 ma[1]= 1.20360; g[1][0]= 0.94128; g[1][1]= 0.55095; g[1][2]= 0; g[1][3]= 0.39065; g[1][4]= 0.31503;
71 ma[2]= 1.55817; g[2][0]= 0.36856; g[2][1]= 0.23888; g[2][2]= 0.55639; g[2][3]= 0.18340; g[2][4]= 0.18681;
72 ma[3]= 1.21000; g[3][0]= 0.33650; g[3][1]= 0.40907; g[3][2]= 0.85679; g[3][3]= 0.19906; g[3][4]= -0.00984;
73 ma[4]= 1.82206; g[4][0]= 0.18171; g[4][1]= -0.17558; g[4][2]= -0.79658; g[4][3]= -0.00355; g[4][4]= 0.22358;
74
75 // Hadronic parameters for tag modes: 0=no-specified, 1=Kpi, 2=Kpipi0, 3=K3pi
76 rd[0] = 0.0;
77 rd[1] = 0.0586;
78 rd[2] = 0.0440;
79 rd[3] = 0.0546;
80 deltad[0] = 0.0;
81 deltad[1] = 194.7*pi180inv;
82 deltad[2] = 196.0*pi180inv;
83 deltad[3] = 167.0*pi180inv;
84 Rf[0] = 0.0;
85 Rf[1] = 1.0;
86 Rf[2] = 0.78;
87 Rf[3] = 0.52;
88
89 return;
90}
91
92complex<double> D0ToKSpipi::Amp_PFT(vector<double> k0l, vector<double> pip, vector<double> pim) {
93 // Breit-Wigner lineshapes
94 vector<double> pD;pD.clear();
95 if(k0l.size()!=4||pip.size()!=4||pim.size()!=4)cout<<"ERROR in KSPIPI daughter 4 momentum"<<endl;
96 for(int i=0;i<k0l.size();i++){
97 pD.push_back(k0l[i] + pip[i] + pim[i]);
98 }
99
100 complex<double> DK2piRes0 = Resonance2(pD, pip, pim, ar[0], phir[0], width_R[0], mass_R[0], spin_R[0]); //ar, phir, width, mass, spin Rho770
101 complex<double> DK2piRes1 = Resonance2(pD, pip, pim, ar[1], phir[1], width_R[1], mass_R[1], spin_R[1]); //ar, phir, width, mass, spin Omega782
102 complex<double> DK2piRes2 = Resonance2(pD, pip, pim, ar[2], phir[2], width_R[2], mass_R[2], spin_R[2]); //ar, phir, width, mass, spin ftwo1270
103 complex<double> DK2piRes3 = Resonance2(pD, pip, pim, ar[3], phir[3], width_R[3], mass_R[3], spin_R[3]); //ar, phir, width, mass, spin Rho1450
104 complex<double> DK2piRes4 = Resonance2(pD, k0l, pim, ar[4], phir[4], width_R[4], mass_R[4], spin_R[4]); //ar, phir, width, mass, spin Kstar892-
105 complex<double> DK2piRes5 = Resonance2(pD, k0l, pim, ar[5], phir[5], width_R[5], mass_R[5], spin_R[5]); //ar, phir, width, mass, spin K2star1430-
106 complex<double> DK2piRes6 = Resonance2(pD, k0l, pim, ar[6], phir[6], width_R[6], mass_R[6], spin_R[6]); //ar, phir, width, mass, spin Kstar1680-
107 complex<double> DK2piRes7 = Resonance2(pD, k0l, pim, ar[7], phir[7], width_R[7], mass_R[7], spin_R[7]); //ar, phir, width, mass, spin Kstar1410-
108 complex<double> DK2piRes8 = Resonance2(pD, k0l, pip, ar[8], phir[8], width_R[8], mass_R[8], spin_R[8]); //ar, phir, width, mass, spin Kstar892+
109 complex<double> DK2piRes9 = Resonance2(pD, k0l, pip, ar[9], phir[9], width_R[9], mass_R[9], spin_R[9]); //ar, phir, width, mass, spin K2star1430+
110 complex<double> DK2piRes10 = Resonance2(pD, k0l, pip, ar[10], phir[10], width_R[10], mass_R[10], spin_R[10]); //ar, phir, width, mass, spin Kstar1410+
111 // K-matrix for pipi S-wave
112 complex<double> pipi_s_wave = K_matrix(pip, pim);
113 if(pipi_s_wave == complex<double>(9999., 9999.)) return 1e-20;
114 // LASS parametrization for Kpi S-wave
115 complex<double> kpi_s_wave = amplitude_LASS(k0l, pip, pim, "k0spim", ar[11], phir[11]*pi180inv);
116 complex<double> dcs_kpi_s_wave = amplitude_LASS(k0l, pip, pim, "k0spip", ar[12], phir[12]*pi180inv);
117
118 complex<double> _tmpAmp = DK2piRes0 + DK2piRes1 + DK2piRes2 + DK2piRes3 + pipi_s_wave;
119 //complex<double> TOT_PFT_AMP = DK2piRes0+ DK2piRes1+ DK2piRes2+ DK2piRes3+ DK2piRes4+ DK2piRes5+ DK2piRes6+ DK2piRes7+ DK2piRes8+ DK2piRes9+ DK2piRes10+ pipi_s_wave + kpi_s_wave+ dcs_kpi_s_wave ;
120 complex<double> TOT_PFT_AMP = _tmpAmp + DK2piRes4+ DK2piRes5+ DK2piRes6+ DK2piRes7+ DK2piRes8+ DK2piRes9+ DK2piRes10 + kpi_s_wave+ dcs_kpi_s_wave ;
121 // Coherent sum for pure-flavor-tagged amplitudes (PFT)
122 return TOT_PFT_AMP;
123}
124
125complex<double> D0ToKSpipi::Resonance2(vector<double> p4_p, vector<double> p4_d1, vector<double> p4_d2, double mag, double theta, double gamma, double bwm, int spin) {
126
127 complex<double> ampl;
128
129 //EvtVector4R p4_d3 = p4_p - p4_d1 - p4_d2;
130 HepLorentzVector _p4_p;_p4_p.setX(p4_p[0]);_p4_p.setY(p4_p[1]);_p4_p.setZ(p4_p[2]);_p4_p.setT(p4_p[3]);
131 HepLorentzVector _p4_d1;_p4_d1.setX(p4_d1[0]);_p4_d1.setY(p4_d1[1]);_p4_d1.setZ(p4_d1[2]);_p4_d1.setT(p4_d1[3]);
132 HepLorentzVector _p4_d2;_p4_d2.setX(p4_d2[0]);_p4_d2.setY(p4_d2[1]);_p4_d2.setZ(p4_d2[2]);_p4_d2.setT(p4_d2[3]);
133 HepLorentzVector _p4_d3=_p4_p-_p4_d1-_p4_d2;
134
135
136 double mAB= (_p4_d1 + _p4_d2).invariantMass();
137 double mBC= (_p4_d2 + _p4_d3).invariantMass();
138 double mAC= (_p4_d1 + _p4_d3).invariantMass();
139 double mA = _p4_d1.invariantMass();
140 double mB = _p4_d2.invariantMass();
141 double mD = _p4_p.invariantMass();
142 double mC = _p4_d3.invariantMass();
143
144
145 double mR = bwm;
146 double gammaR = gamma;
147 double pAB = sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
148 double pR = sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
149
150 double pD= (((mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0) - mR*mR*mC*mC)/(mD*mD);
151 if ( pD>0 ) { pD = sqrt(pD); }
152 else { pD = 0;}
153 double pDAB=sqrt( (((mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0) - mAB*mAB*mC*mC)/(mD*mD));
154 double fR = 1;
155 double fD = 1;
156 int power = 0;
157 switch (spin) {
158 case 0:
159 fR = 1.0;
160 fD = 1.0;
161 power = 1;
162 break;
163 case 1:
164 fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
165 fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
166 power = 3;
167 break;
168 case 2:
169 fR = sqrt( (9+3*pow((1.5*pR),2)+pow((1.5*pR),4))/(9+3*pow((1.5*pAB),2) +pow((1.5*pAB) ,4)) );
170 fD = sqrt( (9+3*pow((5.0*pD),2)+pow((5.0*pD),4))/(9+3*pow((5.0*pDAB),2)+pow((5.0*pDAB),4)) );
171 power = 5;
172 break;
173 default:
174 cout << "Incorrect spin in D0ToKSpipi::EvtResonance2.cc\n" <<endl;
175 }
176
177 double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
178 switch (spin) {
179 case 0:
180 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*fR*fD/(mR*mR-mAB*mAB-complex<double>(0.0,mR*gammaAB));
181 break;
182 case 1:
183 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*
184 (fR*fD*(mAC*mAC-mBC*mBC+((mD*mD-mC*mC)*(mB*mB-mA*mA)/(mAB*mAB)))/(mR*mR-mAB*mAB-complex<double>(0.0,mR*gammaAB)));
185 break;
186 case 2:
187 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*
188 (fR*fD/(mR*mR-mAB*mAB-complex<double>(0.0,mR*gammaAB)))*
189 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mAB*mAB)),2)-
190 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mAB, 2))*
191 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mAB,2)));
192 break;
193 default:
194 cout << "Incorrect spin in D0ToKSpipi::Resonance2.cc\n" <<endl;
195 }
196
197 return ampl;
198}
199
200complex<double> D0ToKSpipi::K_matrix(vector<double> p_pip, vector<double> p_pim) {
201 const double mD0 = 1.86483;
202 const double mKl = 0.49761;
203 const double mPi = 0.13957;
204 bool reject = false;
205
206 HepLorentzVector _p_pip(p_pip[0],p_pip[1],p_pip[2],p_pip[3]);
207 HepLorentzVector _p_pim(p_pim[0],p_pim[1],p_pim[2],p_pim[3]);
208
209 double mAB = (_p_pip + _p_pim).m() ;
210 double s = mAB*mAB;
211
212 complex<double> n11,n12,n13,n14,n15,n21,n22,n23,n24,n25,n31,n32,n33,n34,n35,n41,n42,n43,n44,n45,n51,n52,n53,n54,n55;
213 double rho1sq,rho2sq, rho4sq,rho5sq;
214 complex<double> rho1, rho2,rho3,rho4, rho5;
215 vector< complex<double> > rho;rho.clear();
216 complex<double> pole,SVT,Adler;
217 complex<double> det;
218 //vector< vector< complex<double> > > i;
219 double f[5][5];
220
221 const double mpi = 0.13957;
222 const double mK = 0.493677;
223 const double meta = 0.54775;
224 const double metap = 0.95778;
225
226 // Init matrices and vectors with zeros
227 //vector< vector< complex<double> > > K;
228 complex<double>K[5][5];
229 complex<double>i[5][5];
230 for(int k=0;k<5;k++) {
231 //vector< complex<double> > _itemp;
232 //vector< complex<double> > _Ktemp;
233 for(int l=0;l<5;l++) {
234 //_itemp.push_back(complex<double>(0.,0.));
235 //_Ktemp.push_back(complex<double>(0.,0.));
236 i[k][l]=complex<double>(0.,0.);
237 K[k][l]=complex<double>(0.,0.);
238 f[k][l]=0.;
239 }
240 //i.push_back(_itemp);
241 //K.push_back(_Ktemp);
242 //rho.pus_back(0.);
243 }
244
245 // Fill scattering data values
246 double s_scatt = -3.92637;
247 double sa = 1.0;
248 double sa_0 = -0.15;
249
250 // f_scattering (At least one of the two channels must be pi+pi-)
251 f[0][0] = 0.23399;
252 f[0][1] = 0.15044;
253 f[0][2] = -0.20545;
254 f[0][3] = 0.32825;
255 f[0][4] = 0.35412;
256
257 f[1][0] = f[0][1];
258 f[2][0] = f[0][2];
259 f[3][0] = f[0][3];
260 f[4][0] = f[0][4];
261
262 // Compute phase space factors
263 // rho_0
264 rho1sq=(1.0-(pow((mpi+mpi),2)/s));
265 if(rho1sq >=0.) rho1=complex<double>(sqrt(rho1sq),0.);
266 else rho1=complex<double>(0.,sqrt(-rho1sq));
267 rho.push_back(rho1);
268
269 // rho_1
270 rho2sq=(1.0-(pow((mK+mK),2)/s));
271 if(rho2sq >=0.) rho2=complex<double>(sqrt(rho2sq),0.);
272 else rho2=complex<double>(0.,sqrt(-rho2sq));
273 rho.push_back(rho2);
274
275 // rho_2
276 rho3=complex<double>(0.,0.);
277 if(s<=1) {
278 double real = 1.2274+0.00370909/(s*s) - (0.111203)/(s) - 6.39017*s +16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
279 double cont32=sqrt(1.0-(16.0*mpi*mpi));
280 rho3=complex<double>(cont32*real,0.);
281 }
282 else rho3=complex<double>(sqrt(1.0-(16.0*mpi*mpi/s)),0.);
283 rho.push_back(rho3);
284
285 // rho_3
286 rho4sq=(1.0-(pow((meta+meta),2)/s));
287 if(rho4sq>=0.) rho4=complex<double>(sqrt(rho4sq),0.);
288 else rho4=complex<double>(0.,sqrt(-rho4sq));
289 rho.push_back(rho4);
290
291 // rho_4
292 rho5sq=(1.0-(pow((meta+metap),2)/s));
293 if(rho5sq >=0.) rho5=complex<double>(sqrt(rho5sq),0.);
294 else rho5=complex<double>(0.,sqrt(-rho5sq));
295 rho.push_back(rho5);
296
297 // Sum over the poles [Intermediate channel(k) -> pole(pole_index) -> final channel(l)]
298 for(int k=0;k<5;k++) {
299 for(int l=0;l<5;l++) {
300 for (int pole_index=0;pole_index<5;pole_index++) {
301 double A=g[pole_index][k]*g[pole_index][l];
302 double B=ma[pole_index]*ma[pole_index]-s;
303 K[k][l]=K[k][l]+complex<double>(A/B,0.);
304 }
305 }
306 }
307
308 // Direct scattering term [k -> l]
309 for(int k=0;k<5;k++) {
310 for(int l=0;l<5;l++) {
311 double C=f[k][l]*(1.0-s_scatt);
312 double D=(s-s_scatt);
313 K[k][l]=K[k][l]+complex<double>(C/D,0.);
314 }
315 }
316
317 // Multiplying the "Adler zero" term
318 for(int k=0;k<5;k++) {
319 for(int l=0;l<5;l++) {
320 double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
321 double F=(s-sa_0);
322 K[k][l]=K[k][l]*complex<double>(E/F,0.);
323 }
324 }
325
326 // (1 - i rho K)_ij
327 n11=complex<double>(1.,0.)-complex<double>(0.,1.)*K[0][0]*rho[0];
328 n12=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][1]*rho[1];
329 n13=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][2]*rho[2];
330 n14=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][3]*rho[3];
331 n15=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][4]*rho[4];
332
333 n21=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][0]*rho[0];
334 n22=complex<double>(1.,0.)-complex<double>(0.,1.)*K[1][1]*rho[1];
335 n23=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][2]*rho[2];
336 n24=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][3]*rho[3];
337 n25=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][4]*rho[4];
338
339 n31=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][0]*rho[0];
340 n32=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][1]*rho[1];
341 n33=complex<double>(1.,0.)-complex<double>(0.,1.)*K[2][2]*rho[2];
342 n34=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][3]*rho[3];
343 n35=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][4]*rho[4];
344
345 n41=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][0]*rho[0];
346 n42=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][1]*rho[1];
347 n43=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][2]*rho[2];
348 n44=complex<double>(1.,0.)-complex<double>(0.,1.)*K[3][3]*rho[3];
349 n45=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][4]*rho[4];
350
351 n51=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][0]*rho[0];
352 n52=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][1]*rho[1];
353 n53=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][2]*rho[2];
354 n54=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][3]*rho[3];
355 n55=complex<double>(1.,0.)-complex<double>(0.,1.)*K[4][4]*rho[4];
356
357 // Compute the determinant for inverse [Looks horrible but TMatrixT does not support complex quantities; python bindings may help, working on it.]
358 det = (n15*n24*n33*n42*n51 - n14*n25*n33*n42*n51 - n15*n23*n34*n42*n51 +
359 n13*n25*n34*n42*n51 + n14*n23*n35*n42*n51 - n13*n24*n35*n42*n51 -
360 n15*n24*n32*n43*n51 + n14*n25*n32*n43*n51 + n15*n22*n34*n43*n51 -
361 n12*n25*n34*n43*n51 - n14*n22*n35*n43*n51 + n12*n24*n35*n43*n51 +
362 n15*n23*n32*n44*n51 - n13*n25*n32*n44*n51 - n15*n22*n33*n44*n51 +
363 n12*n25*n33*n44*n51 + n13*n22*n35*n44*n51 - n12*n23*n35*n44*n51 -
364 n14*n23*n32*n45*n51 + n13*n24*n32*n45*n51 + n14*n22*n33*n45*n51 -
365 n12*n24*n33*n45*n51 - n13*n22*n34*n45*n51 + n12*n23*n34*n45*n51 -
366 n15*n24*n33*n41*n52 + n14*n25*n33*n41*n52 + n15*n23*n34*n41*n52 -
367 n13*n25*n34*n41*n52 - n14*n23*n35*n41*n52 + n13*n24*n35*n41*n52 +
368 n15*n24*n31*n43*n52 - n14*n25*n31*n43*n52 - n15*n21*n34*n43*n52 +
369 n11*n25*n34*n43*n52 + n14*n21*n35*n43*n52 - n11*n24*n35*n43*n52 -
370 n15*n23*n31*n44*n52 + n13*n25*n31*n44*n52 + n15*n21*n33*n44*n52 -
371 n11*n25*n33*n44*n52 - n13*n21*n35*n44*n52 + n11*n23*n35*n44*n52 +
372 n14*n23*n31*n45*n52 - n13*n24*n31*n45*n52 - n14*n21*n33*n45*n52 +
373 n11*n24*n33*n45*n52 + n13*n21*n34*n45*n52 - n11*n23*n34*n45*n52 +
374 n15*n24*n32*n41*n53 - n14*n25*n32*n41*n53 - n15*n22*n34*n41*n53 +
375 n12*n25*n34*n41*n53 + n14*n22*n35*n41*n53 - n12*n24*n35*n41*n53 -
376 n15*n24*n31*n42*n53 + n14*n25*n31*n42*n53 + n15*n21*n34*n42*n53 -
377 n11*n25*n34*n42*n53 - n14*n21*n35*n42*n53 + n11*n24*n35*n42*n53 +
378 n15*n22*n31*n44*n53 - n12*n25*n31*n44*n53 - n15*n21*n32*n44*n53 +
379 n11*n25*n32*n44*n53 + n12*n21*n35*n44*n53 - n11*n22*n35*n44*n53 -
380 n14*n22*n31*n45*n53 + n12*n24*n31*n45*n53 + n14*n21*n32*n45*n53 -
381 n11*n24*n32*n45*n53 - n12*n21*n34*n45*n53 + n11*n22*n34*n45*n53 -
382 n15*n23*n32*n41*n54 + n13*n25*n32*n41*n54 + n15*n22*n33*n41*n54 -
383 n12*n25*n33*n41*n54 - n13*n22*n35*n41*n54 + n12*n23*n35*n41*n54 +
384 n15*n23*n31*n42*n54 - n13*n25*n31*n42*n54 - n15*n21*n33*n42*n54 +
385 n11*n25*n33*n42*n54 + n13*n21*n35*n42*n54 - n11*n23*n35*n42*n54 -
386 n15*n22*n31*n43*n54 + n12*n25*n31*n43*n54 + n15*n21*n32*n43*n54 -
387 n11*n25*n32*n43*n54 - n12*n21*n35*n43*n54 + n11*n22*n35*n43*n54 +
388 n13*n22*n31*n45*n54 - n12*n23*n31*n45*n54 - n13*n21*n32*n45*n54 +
389 n11*n23*n32*n45*n54 + n12*n21*n33*n45*n54 - n11*n22*n33*n45*n54 +
390 n14*n23*n32*n41*n55 - n13*n24*n32*n41*n55 - n14*n22*n33*n41*n55 +
391 n12*n24*n33*n41*n55 + n13*n22*n34*n41*n55 - n12*n23*n34*n41*n55 -
392 n14*n23*n31*n42*n55 + n13*n24*n31*n42*n55 + n14*n21*n33*n42*n55 -
393 n11*n24*n33*n42*n55 - n13*n21*n34*n42*n55 + n11*n23*n34*n42*n55 +
394 n14*n22*n31*n43*n55 - n12*n24*n31*n43*n55 - n14*n21*n32*n43*n55 +
395 n11*n24*n32*n43*n55 + n12*n21*n34*n43*n55 - n11*n22*n34*n43*n55 -
396 n13*n22*n31*n44*n55 + n12*n23*n31*n44*n55 + n13*n21*n32*n44*n55 -
397 n11*n23*n32*n44*n55 - n12*n21*n33*n44*n55 + n11*n22*n33*n44*n55);
398
399 if(det == complex<double>(0., 0.)) reject=true;
400
401 // The 1st row of the inverse matrix [(1-i\rhoK)^-1]_0j
402 i[0][0] = ( n25*n34*n43*n52 -
403 n24*n35*n43*n52 - n25*n33*n44*n52 + n23*n35*n44*n52 +
404 n24*n33*n45*n52 - n23*n34*n45*n52 - n25*n34*n42*n53 +
405 n24*n35*n42*n53 + n25*n32*n44*n53 - n22*n35*n44*n53 -
406 n24*n32*n45*n53 + n22*n34*n45*n53 + n25*n33*n42*n54 -
407 n23*n35*n42*n54 - n25*n32*n43*n54 + n22*n35*n43*n54 +
408 n23*n32*n45*n54 - n22*n33*n45*n54 - n24*n33*n42*n55 +
409 n23*n34*n42*n55 + n24*n32*n43*n55 - n22*n34*n43*n55 -
410 n23*n32*n44*n55 + n22*n33*n44*n55)/det;
411
412 i[0][1] = ( -n15*n34*n43*n52 +
413 n14*n35*n43*n52 + n15*n33*n44*n52 - n13*n35*n44*n52 -
414 n14*n33*n45*n52 + n13*n34*n45*n52 + n15*n34*n42*n53 -
415 n14*n35*n42*n53 - n15*n32*n44*n53 + n12*n35*n44*n53 +
416 n14*n32*n45*n53 - n12*n34*n45*n53 - n15*n33*n42*n54 +
417 n13*n35*n42*n54 + n15*n32*n43*n54 - n12*n35*n43*n54 -
418 n13*n32*n45*n54 + n12*n33*n45*n54 + n14*n33*n42*n55 -
419 n13*n34*n42*n55 - n14*n32*n43*n55 + n12*n34*n43*n55 +
420 n13*n32*n44*n55 - n12*n33*n44*n55)/det;
421
422 i[0][2] = ( n15*n24*n43*n52 -
423 n14*n25*n43*n52 - n15*n23*n44*n52 + n13*n25*n44*n52 +
424 n14*n23*n45*n52 - n13*n24*n45*n52 - n15*n24*n42*n53 +
425 n14*n25*n42*n53 + n15*n22*n44*n53 - n12*n25*n44*n53 -
426 n14*n22*n45*n53 + n12*n24*n45*n53 + n15*n23*n42*n54 -
427 n13*n25*n42*n54 - n15*n22*n43*n54 + n12*n25*n43*n54 +
428 n13*n22*n45*n54 - n12*n23*n45*n54 - n14*n23*n42*n55 +
429 n13*n24*n42*n55 + n14*n22*n43*n55 - n12*n24*n43*n55 -
430 n13*n22*n44*n55 + n12*n23*n44*n55)/det;
431
432 i[0][3] = ( -n15*n24*n33*n52 +
433 n14*n25*n33*n52 + n15*n23*n34*n52 - n13*n25*n34*n52 -
434 n14*n23*n35*n52 + n13*n24*n35*n52 + n15*n24*n32*n53 -
435 n14*n25*n32*n53 - n15*n22*n34*n53 + n12*n25*n34*n53 +
436 n14*n22*n35*n53 - n12*n24*n35*n53 - n15*n23*n32*n54 +
437 n13*n25*n32*n54 + n15*n22*n33*n54 - n12*n25*n33*n54 -
438 n13*n22*n35*n54 + n12*n23*n35*n54 + n14*n23*n32*n55 -
439 n13*n24*n32*n55 - n14*n22*n33*n55 + n12*n24*n33*n55 +
440 n13*n22*n34*n55 - n12*n23*n34*n55)/det;
441
442 i[0][4] = ( n15*n24*n33*n42 -
443 n14*n25*n33*n42 - n15*n23*n34*n42 + n13*n25*n34*n42 +
444 n14*n23*n35*n42 - n13*n24*n35*n42 - n15*n24*n32*n43 +
445 n14*n25*n32*n43 + n15*n22*n34*n43 - n12*n25*n34*n43 -
446 n14*n22*n35*n43 + n12*n24*n35*n43 + n15*n23*n32*n44 -
447 n13*n25*n32*n44 - n15*n22*n33*n44 + n12*n25*n33*n44 +
448 n13*n22*n35*n44 - n12*n23*n35*n44 - n14*n23*n32*n45 +
449 n13*n24*n32*n45 + n14*n22*n33*n45 - n12*n24*n33*n45 -
450 n13*n22*n34*n45 + n12*n23*n34*n45)/det;
451
452 double s0_prod = -0.07;
453
454 complex<double> value0(0., 0.);
455 complex<double> value1(0., 0.);
456
457 // [(1-i\rhoK)^-1]_0j*P_j {P_j: Production vector}
458 for(int k=0;k<5;k++) {
459 double u1j_re = real(i[0][k]);
460 double u1j_im = imag(i[0][k]);
461 if(u1j_re==0. || u1j_im==0.) reject=true;
462
463 // Initial state to K-matrix pole couplings * Pole to intermediate channels coupling
464 for(int pole_index=0;pole_index<5;pole_index++) {
465 complex<double> A = beta[pole_index]*g[pole_index][k];
466 value0 += (i[0][k]*A)/(ma[pole_index]*ma[pole_index]-s);
467 }
468
469 // Direct initial state to intermediate channels couplings
470 value1 += i[0][k]*fprod[k];
471
472 }
473
474 // Slowly varying polynomial term for the direct coupling
475 value1 *= (1.-s0_prod)/(s-s0_prod) ;
476
477 if(reject==true) return complex<double>(9999., 9999.);
478 else return (value0+value1);
479
480}
481
482complex<double> D0ToKSpipi::amplitude_LASS(vector<double> p_k0l, vector<double> p_pip, vector<double> p_pim, string reso, double A_r, double Phi_r) {
483 double mR = 1.425 ;
484 double gammaR = 0.27 ;
485 double mab2 = 0.0;
486
487 HepLorentzVector _p_k0l(p_k0l[0],p_k0l[1],p_k0l[2],p_k0l[3]);
488 HepLorentzVector _p_pip(p_pip[0],p_pip[1],p_pip[2],p_pip[3]);
489 HepLorentzVector _p_pim(p_pim[0],p_pim[1],p_pim[2],p_pim[3]);
490 if (reso == "k0spim") mab2 = pow((_p_k0l + _p_pim).m(),2);
491 else if(reso == "k0spip") mab2 = pow((_p_k0l + _p_pip).m(),2);
492 double s = mab2;
493
494 const double mD0 = 1.86483;
495 const double mKl = 0.49761;
496 const double mPi = 0.13957;
497
498 double _a = 0.113;
499 double _r = -33.8;
500 double _R = 1.0;
501 double _F = 0.96;
502 double _phiR = -1.9146;
503 double _phiF = 0.0017;
504 double fR = 1.0; // K*0(1430) has spin zero
505 int power = 1; // Power is 1 for spin zero
506
507 double mAB = sqrt(mab2);
508 double mA = mKl;
509 double mB = mPi;
510 double mC = mPi;
511 double mD = mD0;
512
513 double pAB=sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
514 double q=pAB;
515
516 double pR=sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
517 double q0=pR;
518
519 // Running width.
520 double g = gammaR*pow(q/q0,power)*(mR/mAB)*fR*fR;
521 complex<double> propagator_relativistic_BreitWigner = 1./(mR*mR - mAB*mAB - complex<double>(0.,mR*g));
522
523 // Non-resonant phase shift
524 double cot_deltaF = 1.0/(_a*q) + 0.5*_r*q;
525 double qcot_deltaF = 1.0/_a + 0.5*_r*q*q;
526
527 // Compute resonant part
528 complex<double> expi2deltaF = complex<double>(qcot_deltaF, q)/ complex<double>(qcot_deltaF, -q);
529 complex<double> resonant_term_T = _R * complex<double>(cos(_phiR + 2 * _phiF), sin(_phiR + 2 * _phiF)) * propagator_relativistic_BreitWigner * mR * gammaR * mR / q0 * expi2deltaF;
530
531 // Compute non-resonant part
532 complex<double> non_resonant_term_F = _F * complex<double>(cos(_phiF), sin(_phiF)) * (cos(_phiF) + cot_deltaF * sin(_phiF)) * sqrt(s) / complex<double>(qcot_deltaF, -q);
533
534 // Add non-resonant and resonant terms
535 complex<double> LASS_contribution = non_resonant_term_F + resonant_term_T;
536
537 return complex<double>(A_r*cos(Phi_r), A_r*sin(Phi_r)) * LASS_contribution;
538}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:246
meta
double meta
Definition EvtPhokharaDef.hh:21
mPi
double mPi
Definition FarmMonitorAlg.cxx:30
mpi
const double mpi
Definition Gam4pikp.cxx:47
s
XmlRpcServer s
Definition HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
mD0
const double mD0
Definition MyConst.h:5
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
D0ToKSpipi::Amp_PFT
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition D0ToKSpipi.cxx:92
D0ToKSpipi::init
void init()
Definition D0ToKSpipi.cxx:25
D0ToKSpipi::~D0ToKSpipi
virtual ~D0ToKSpipi()
Definition D0ToKSpipi.cxx:23