BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
D0ToKLpipi.cxx
Go to the documentation of this file.
1#include "QCMCFilterAlg/D0ToKLpipi.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;
21
22using namespace std; //::endl;
23
24D0ToKLpipi::~D0ToKLpipi() {}
25
26void D0ToKLpipi::init(){
27 //std::cout << "D0ToKLpipi ==> Initialization !" << std::endl;
28
29 _nd = 3;
30 tan2thetaC = (0.22650*0.22650)/(1.-(0.22650*0.22650)) ; //sin(theta_C) = 0.22650 +/- 0.00048
31 pi180inv = 1.0*3.1415926/180;
32
33 //spin[0] = 1 ; //Rho770
34 //spin[1] = 1 ; //Omega782
35 //spin[2] = 2 ; //ftwo1270
36 //spin[3] = 1 ; //Rho1450
37 //spin[4] = 1 ; //Kstar892minus
38 //spin[5] = 2 ; //K2star1430minus
39 //spin[6] = 1 ; //Kstar1680minus
40 //spin[7] = 1 ; //Kstar1410minus
41 //spin[8] = 1 ; //Kstar892plus
42 //spin[9] = 2 ; //K2star1430plus
43 //spin[10]= 1 ; //Kstar1410plus
44 mass_R[0]= 0.77526; width_R[0]= 0.14740; spin_R[0]= 1; ar[0]= 1 ; phir[0]= 0 ;
45 mass_R[1]= 0.78266; width_R[1]= 0.00868; spin_R[1]= 1; ar[1]= 0.03928; phir[1]= 106.1 ;
46 mass_R[2]= 1.27550; width_R[2]= 0.18670; spin_R[2]= 2; ar[2]= 1.302 ; phir[2]= -39.97;
47 mass_R[3]= 1.46500; width_R[3]= 0.40000; spin_R[3]= 1; ar[3]= 1.656 ; phir[3]= 113 ;
48 mass_R[4]= 0.89167; width_R[4]= 0.0514; spin_R[4]= 1; ar[4]= 1.843 ; phir[4]= 138 ;
49 mass_R[5]= 1.42730; width_R[5]= 0.10000; spin_R[5]= 2; ar[5]= 1.537 ; phir[5]= -49.45;
50 mass_R[6]= 1.71800; width_R[6]= 0.3220; spin_R[6]= 1; ar[6]= 2.251 ; phir[6]= -196.3;
51 mass_R[7]= 1.41400; width_R[7]= 0.2320; spin_R[7]= 1; ar[7]= 0.4733 ; phir[7]= 1659 ;
52 mass_R[8]= 0.89167; width_R[8]= 0.0514; spin_R[8]= 1; ar[8]= 0.1617 ; phir[8]= -36.25;
53 mass_R[9]= 1.42730; width_R[9]= 0.1000; spin_R[9]= 2; ar[9]= 0.2267 ; phir[9]= -84.64;
54 mass_R[10]= 1.41400; width_R[10]= 0.2320; spin_R[10]= 1; ar[10]= 0.2185 ; phir[10]= 73.12 ;
55 mass_R[11]= 1.42500; width_R[11]= 0.2700; spin_R[11]= 1; ar[11]= 2.396 ; phir[11]= 98.14 ;
56 // beta[kb] = EvtComplex(mag*cos(phase*pi180inv), mag*sin(phase*pi180inv)) ;
57 // fprod[kf] = EvtComplex(mag*cos(phase*pi180inv), mag*sin(phase*pi180inv)) ;
58 beta[0] = complex<double>(8.5 *cos(68.5 *pi180inv), 8.5 *sin(68.5 *pi180inv));
59 beta[1] = complex<double>(12.2*cos(24 *pi180inv), 12.2 *sin(24 *pi180inv));
60 beta[2] = complex<double>(29.2*cos(-0.1 *pi180inv), 29.2 *sin(-0.1 *pi180inv));
61 beta[3] = complex<double>(10.8*cos(-51.9*pi180inv), 10.8 *sin(-51.9*pi180inv));
62 beta[4] = complex<double>(0., 0.);
63
64 fprod[0] = complex<double>(8 *cos(-126 *pi180inv), 8 *sin(-126 *pi180inv));
65 fprod[1] = complex<double>(26.3*cos(-152.3*pi180inv), 26.3*sin(-152.3*pi180inv));
66 fprod[2] = complex<double>(33 *cos(-93.2 *pi180inv), 33 *sin(-93.2 *pi180inv));
67 fprod[3] = complex<double>(26.2*cos(-121.4*pi180inv), 26.2*sin(-121.4*pi180inv));
68 fprod[4] = complex<double>(0., 0.);
69
70 //CP_mult[w] = (EvtComplex(1., 0.) - 2.*tan2thetaC*EvtComplex(r*cos(delta), r*sin(delta))) ;
71 CP_mult[0] = complex<double>(1.,0.)-2.*tan2thetaC*complex<double>(1.851 *cos(-94.07 *pi180inv),1.851 *sin(-94.07 *pi180inv));
72 CP_mult[1] = complex<double>(1.,0.)-2.*tan2thetaC*complex<double>(6.332 *cos(2.103 *pi180inv),6.332 *sin(2.103 *pi180inv));
73 CP_mult[2] = complex<double>(1.,0.)-2.*tan2thetaC*complex<double>(3.229 *cos(-60.05 *pi180inv),3.229 *sin(-60.05 *pi180inv));
74 CP_mult[3] = complex<double>(1.,0.)-2.*tan2thetaC*complex<double>(12.75 *cos(73.62 *pi180inv),12.75 *sin(73.62 *pi180inv));
75 CP_mult[4] = complex<double>(1.,0.)-2.*tan2thetaC*complex<double>(0.7116*cos(-177.149*pi180inv),0.7116*sin(-177.149*pi180inv));
76
77 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;
78 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;
79 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;
80 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;
81 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;
82
83 // Hadronic parameters for tag modes: 0=no-specified, 1=Kpi, 2=Kpipi0, 3=K3pi
84 rd[0] = 0.0;
85 rd[1] = 0.0586;
86 rd[2] = 0.0440;
87 rd[3] = 0.0546;
88 deltad[0] = 0.0;
89 deltad[1] = 194.7*pi180inv;
90 deltad[2] = 196.0*pi180inv;
91 deltad[3] = 167.0*pi180inv;
92 Rf[0] = 0.0;
93 Rf[1] = 1.0;
94 Rf[2] = 0.78;
95 Rf[3] = 0.52;
96
97 return;
98}
99
100complex<double> D0ToKLpipi::Amp_PFT(vector<double> k0l, vector<double> pip, vector<double> pim) {
101
102 vector<double> pD;pD.clear();
103 if(k0l.size()!=4||pip.size()!=4||pim.size()!=4)cout<<"ERROR in KSPIPI daughter 4 momentum"<<endl;
104 for(int i=0;i<k0l.size();i++){
105 pD.push_back(k0l[i] + pip[i] + pim[i]);
106 }
107
108 //EvtResonance2 DK2piRes0(pD, pip, pim, ar[0], phir[0], width_R[0], mass_R[0], spin[0]); //ar, phir, width, mass, spin //Rho770
109 //EvtResonance2 DK2piRes1(pD, pip, pim, ar[1], phir[1], width_R[1], mass_R[1], spin[1]); //ar, phir, width, mass, spin //Omega782
110 //EvtResonance2 DK2piRes2(pD, pip, pim, ar[2], phir[2], width_R[2], mass_R[2], spin[2]); //ar, phir, width, mass, spin //ftwo1270
111 //EvtResonance2 DK2piRes3(pD, pip, pim, ar[3], phir[3], width_R[3], mass_R[3], spin[3]); //ar, phir, width, mass, spin //Rho1450
112 //EvtResonance2 DK2piRes4(pD, k0l, pim, ar[4], phir[4], width_R[4], mass_R[4], spin[4]); //ar, phir, width, mass, spin //Kstar892minus
113 //EvtResonance2 DK2piRes5(pD, k0l, pim, ar[5], phir[5], width_R[5], mass_R[5], spin[5]); //ar, phir, width, mass, spin //K2star1430minus
114 //EvtResonance2 DK2piRes6(pD, k0l, pim, ar[6], phir[6], width_R[6], mass_R[6], spin[6]); //ar, phir, width, mass, spin //Kstar1680minus
115 //EvtResonance2 DK2piRes7(pD, k0l, pim, ar[7], phir[7], width_R[7], mass_R[7], spin[7]); //ar, phir, width, mass, spin //Kstar1410minus
116 //EvtResonance2 DK2piRes8(pD, k0l, pip, ar[8], phir[8], width_R[8], mass_R[8], spin[8]); //ar, phir, width, mass, spin //Kstar892plus
117 //EvtResonance2 DK2piRes9(pD, k0l, pip, ar[9], phir[9], width_R[9], mass_R[9], spin[9]); //ar, phir, width, mass, spin //K2star1430plus
118 //EvtResonance2 DK2piRes10(pD, k0l, pip, ar[10], phir[10], width_R[10], mass_R[10], spin[10]); //ar, phir, width, mass, spin //Kstar1410plus
119 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
120 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
121 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
122 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
123 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-
124 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-
125 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-
126 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-
127 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+
128 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+
129 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+
130
131 complex<double> pipi_s_wave = K_matrix(pip, pim) ;
132 if(pipi_s_wave == complex<double>(9999., 9999.)) return 1e-20 ;
133
134 complex<double> kpi_s_wave = amplitude_LASS(k0l, pip, pim, "k0lpim", ar[11], phir[11]*pi180inv) ;
135 //complex<double> kpi_s_wave_dcs = amplitude_LASS(k0l, pip, pim, "k0spip", ar[12], phir[12]*pi180inv) ; should be there but not observed yet GUESS
136
137 complex<double> TOT_PFT_AMP = DK2piRes0 * CP_mult[0]
138 + DK2piRes1 * CP_mult[1]
139 + DK2piRes2 * CP_mult[2]
140 + DK2piRes3 * CP_mult[3]
141 + DK2piRes4
142 + DK2piRes5
143 + DK2piRes6
144 + DK2piRes7
145 + DK2piRes8 * (-1.)
146 + DK2piRes9 * (-1.)
147 + DK2piRes10* (-1.)
148 + pipi_s_wave * CP_mult[4]
149 + kpi_s_wave ;
150
151
152 return TOT_PFT_AMP;
153
154}
155
156complex<double> D0ToKLpipi::Resonance2(vector<double> p4_p, vector<double> p4_d1, vector<double> p4_d2, double mag, double theta, double gamma, double bwm, int spin) {
157
158 complex<double> ampl;
159
160 //EvtVector4R p4_d3 = p4_p - p4_d1 - p4_d2;
161 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]);
162 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]);
163 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]);
164 HepLorentzVector _p4_d3=_p4_p-_p4_d1-_p4_d2;
165
166
167 double mAB= (_p4_d1 + _p4_d2).invariantMass();
168 double mBC= (_p4_d2 + _p4_d3).invariantMass();
169 double mAC= (_p4_d1 + _p4_d3).invariantMass();
170 double mA = _p4_d1.invariantMass();
171 double mB = _p4_d2.invariantMass();
172 double mD = _p4_p.invariantMass();
173 double mC = _p4_d3.invariantMass();
174
175
176 double mR = bwm;
177 double gammaR = gamma;
178 double pAB = sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
179 double pR = sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
180
181 double pD= (((mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0) - mR*mR*mC*mC)/(mD*mD);
182 if ( pD>0 ) { pD = sqrt(pD); }
183 else { pD = 0;}
184 double pDAB=sqrt( (((mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0) - mAB*mAB*mC*mC)/(mD*mD));
185 double fR = 1;
186 double fD = 1;
187 int power = 0;
188 switch (spin) {
189 case 0:
190 fR = 1.0;
191 fD = 1.0;
192 power = 1;
193 break;
194 case 1:
195 fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
196 fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
197 power = 3;
198 break;
199 case 2:
200 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)) );
201 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)) );
202 power = 5;
203 break;
204 default:
205 cout << "Incorrect spin in D0ToKLpipi::EvtResonance2.cc\n" <<endl;
206 }
207
208 double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
209 switch (spin) {
210 case 0:
211 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*fR*fD/(mR*mR-mAB*mAB-complex<double>(0.0,mR*gammaAB));
212 break;
213 case 1:
214 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*
215 (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)));
216 break;
217 case 2:
218 ampl=mag*complex<double>(cos(theta*pi180inv),sin(theta*pi180inv))*
219 (fR*fD/(mR*mR-mAB*mAB-complex<double>(0.0,mR*gammaAB)))*
220 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mAB*mAB)),2)-
221 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mAB, 2))*
222 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mAB,2)));
223 break;
224 default:
225 cout << "Incorrect spin in D0ToKSpipi::Resonance2.cc\n" <<endl;
226 }
227
228 return ampl;
229}
230
231complex<double> D0ToKLpipi::K_matrix(vector<double> p_pip, vector<double> p_pim) {
232
233 //double pi180inv = 1.0/EvtConst::radToDegrees;
234 bool reject=false;
235
236 const double mD0 = 1.86483;
237 const double mKl = 0.49761;
238 const double mPi = 0.13957;
239
240 HepLorentzVector _p_pip(p_pip[0],p_pip[1],p_pip[2],p_pip[3]);
241 HepLorentzVector _p_pim(p_pim[0],p_pim[1],p_pim[2],p_pim[3]);
242
243 double mAB = (_p_pip + _p_pim).m() ;
244
245 double s = mAB*mAB;
246
247 // Define the complex coupling constants
248 //Double_t g[5][5]; // g[Physical pole]Decay channel]
249
250 // pi+pi- channel
251 //g[0][0]=0.22889;
252 //g[1][0]=0.94128;
253 //g[2][0]=0.36856;
254 //g[3][0]=0.33650;
255 //g[4][0]=0.18171;
256
257 //// K+K- channel
258 //g[0][1]=-0.55377;
259 //g[1][1]=0.55095;
260 //g[2][1]=0.23888;
261 //g[3][1]=0.40907;
262 //g[4][1]=-0.17558;
263
264 //// 4pi channel
265 //g[0][2]=0;
266 //g[1][2]=0;
267 //g[2][2]=0.55639;
268 //g[3][2]=0.85679;
269 //g[4][2]=-0.79658;
270
271 //// eta eta channel
272 //g[0][3]=-0.39899;
273 //g[1][3]=0.39065;
274 //g[2][3]=0.18340;
275 //g[3][3]=0.19906;
276 //g[4][3]=-0.00355;
277
278 ////eta eta' channel
279 //g[0][4]=-0.34639;
280 //g[1][4]=0.31503;
281 //g[2][4]=0.18681;
282 //g[3][4]=-0.00984;
283 //g[4][4]=0.22358;
284
285 // Define masses of the physical poles (in GeV)
286 //Double_t ma[5];
287
288 //ma[0]=0.651;
289 //ma[1]=1.20360;
290 //ma[2]=1.55817;
291 //ma[3]=1.21000;
292 //ma[4]=1.82206;
293
294 // Define variables
295 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;
296 double rho1sq,rho2sq,rho4sq,rho5sq;
297 complex<double> rho1,rho2,rho3,rho4,rho5;
298 complex<double> rho[5];
299 complex<double> pole,SVT,Adler;
300 complex<double> det;
301 complex<double> i[5][5];
302 double f[5][5];
303
304 // pi+, K+, eta, and eta' PDG masses
305 double mpi=0.13957;
306 double mK=0.493677;
307 double meta=0.54775;
308 double metap=0.95778;
309
310 // Init matrices and vectors with zeros
311 complex<double> K[5][5];
312 for(Int_t k=0;k<5;k++) {
313 for(Int_t l=0;l<5;l++) {
314 i[k][l]=complex<double>(0.,0.);
315 K[k][l]=complex<double>(0.,0.);
316 f[k][l]=0.;
317 }
318 rho[k]=0.;
319 }
320
321 // Fill scattering data values
322 Double_t s_scatt=-3.92637;
323 Double_t sa=1.0;
324 Double_t sa_0=-0.15;
325
326 // f_scattering
327 f[0][0]=0.23399;
328 f[0][1]=0.15044;
329 f[0][2]=-0.20545;
330 f[0][3]=0.32825;
331 f[0][4]=0.35412;
332
333 f[1][0]=f[0][1];
334 f[2][0]=f[0][2];
335 f[3][0]=f[0][3];
336 f[4][0]=f[0][4];
337
338 // Compute phase space factors
339 rho1sq=(1.0-(pow((mpi+mpi),2)/s));
340 if(rho1sq >=0.) {
341 rho1=complex<double>(sqrt(rho1sq),0.);
342 }
343 else{
344 rho1=complex<double>(0.,sqrt(-rho1sq));
345 }
346 rho[0]=rho1;
347
348 rho2sq=(1.0-(pow((mK+mK),2)/s));
349 if(rho2sq >=0.) {
350 rho2=complex<double>(sqrt(rho2sq),0.);
351 }
352 else{
353 rho2=complex<double>(0.,sqrt(-rho2sq));
354 }
355
356 rho[1]=rho2;
357
358 rho3=complex<double>(0.,0.);
359
360 if(s<=1) {
361 Double_t 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;
362 Double_t cont32=sqrt(1.0-(16.0*mpi*mpi));
363 rho3=complex<double>(cont32*real,0.);
364 }
365 else{
366 rho3=complex<double>(sqrt(1.0-(16.0*mpi*mpi/s)),0.);
367 }
368 rho[2]=rho3;
369
370 rho4sq=(1.0-(pow((meta+meta),2)/s));
371 if(rho4sq>=0.) {
372 rho4=complex<double>(sqrt(rho4sq),0.);
373 }
374 else{
375 rho4=complex<double>(0.,sqrt(-rho4sq));
376 }
377 rho[3]=rho4;
378
379 rho5sq=(1.0-(pow((meta+metap),2)/s));
380 if(rho5sq >=0.) {
381 rho5=complex<double>(sqrt(rho5sq),0.);
382 }
383 else{
384 rho5=complex<double>(0.,sqrt(-rho5sq));
385 }
386 rho[4]=rho5;
387
388 // Sum over the poles
389 for(Int_t k=0;k<5;k++) {
390 for(Int_t l=0;l<5;l++) {
391 for (Int_t pole_index=0;pole_index<5;pole_index++) {
392 Double_t A=g[pole_index][k]*g[pole_index][l];
393 Double_t B=ma[pole_index]*ma[pole_index]-s;
394 K[k][l]=K[k][l]+complex<double>(A/B,0.);
395 }
396 }
397 }
398
399 for(Int_t k=0;k<5;k++) {
400 for(Int_t l=0;l<5;l++) {
401 Double_t C=f[k][l]*(1.0-s_scatt);
402 Double_t D=(s-s_scatt);
403 K[k][l]=K[k][l]+complex<double>(C/D,0.);
404 }
405 }
406
407 for(Int_t k=0;k<5;k++) {
408 for(Int_t l=0;l<5;l++) {
409 Double_t E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
410 Double_t F=(s-sa_0);
411 K[k][l]=K[k][l]*complex<double>(E/F,0.);
412 }
413 }
414
415 n11=complex<double>(1.,0.)-complex<double>(0.,1.)*K[0][0]*rho[0];
416 n12=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][1]*rho[1];
417 n13=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][2]*rho[2];
418 n14=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][3]*rho[3];
419 n15=complex<double>(0.,0.)-complex<double>(0.,1.)*K[0][4]*rho[4];
420
421 n21=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][0]*rho[0];
422 n22=complex<double>(1.,0.)-complex<double>(0.,1.)*K[1][1]*rho[1];
423 n23=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][2]*rho[2];
424 n24=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][3]*rho[3];
425 n25=complex<double>(0.,0.)-complex<double>(0.,1.)*K[1][4]*rho[4];
426
427 n31=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][0]*rho[0];
428 n32=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][1]*rho[1];
429 n33=complex<double>(1.,0.)-complex<double>(0.,1.)*K[2][2]*rho[2];
430 n34=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][3]*rho[3];
431 n35=complex<double>(0.,0.)-complex<double>(0.,1.)*K[2][4]*rho[4];
432
433 n41=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][0]*rho[0];
434 n42=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][1]*rho[1];
435 n43=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][2]*rho[2];
436 n44=complex<double>(1.,0.)-complex<double>(0.,1.)*K[3][3]*rho[3];
437 n45=complex<double>(0.,0.)-complex<double>(0.,1.)*K[3][4]*rho[4];
438
439 n51=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][0]*rho[0];
440 n52=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][1]*rho[1];
441 n53=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][2]*rho[2];
442 n54=complex<double>(0.,0.)-complex<double>(0.,1.)*K[4][3]*rho[3];
443 n55=complex<double>(1.,0.)-complex<double>(0.,1.)*K[4][4]*rho[4];
444
445 // Compute the determinant
446 det = (n15*n24*n33*n42*n51 - n14*n25*n33*n42*n51 - n15*n23*n34*n42*n51 +
447 n13*n25*n34*n42*n51 + n14*n23*n35*n42*n51 - n13*n24*n35*n42*n51 -
448 n15*n24*n32*n43*n51 + n14*n25*n32*n43*n51 + n15*n22*n34*n43*n51 -
449 n12*n25*n34*n43*n51 - n14*n22*n35*n43*n51 + n12*n24*n35*n43*n51 +
450 n15*n23*n32*n44*n51 - n13*n25*n32*n44*n51 - n15*n22*n33*n44*n51 +
451 n12*n25*n33*n44*n51 + n13*n22*n35*n44*n51 - n12*n23*n35*n44*n51 -
452 n14*n23*n32*n45*n51 + n13*n24*n32*n45*n51 + n14*n22*n33*n45*n51 -
453 n12*n24*n33*n45*n51 - n13*n22*n34*n45*n51 + n12*n23*n34*n45*n51 -
454 n15*n24*n33*n41*n52 + n14*n25*n33*n41*n52 + n15*n23*n34*n41*n52 -
455 n13*n25*n34*n41*n52 - n14*n23*n35*n41*n52 + n13*n24*n35*n41*n52 +
456 n15*n24*n31*n43*n52 - n14*n25*n31*n43*n52 - n15*n21*n34*n43*n52 +
457 n11*n25*n34*n43*n52 + n14*n21*n35*n43*n52 - n11*n24*n35*n43*n52 -
458 n15*n23*n31*n44*n52 + n13*n25*n31*n44*n52 + n15*n21*n33*n44*n52 -
459 n11*n25*n33*n44*n52 - n13*n21*n35*n44*n52 + n11*n23*n35*n44*n52 +
460 n14*n23*n31*n45*n52 - n13*n24*n31*n45*n52 - n14*n21*n33*n45*n52 +
461 n11*n24*n33*n45*n52 + n13*n21*n34*n45*n52 - n11*n23*n34*n45*n52 +
462 n15*n24*n32*n41*n53 - n14*n25*n32*n41*n53 - n15*n22*n34*n41*n53 +
463 n12*n25*n34*n41*n53 + n14*n22*n35*n41*n53 - n12*n24*n35*n41*n53 -
464 n15*n24*n31*n42*n53 + n14*n25*n31*n42*n53 + n15*n21*n34*n42*n53 -
465 n11*n25*n34*n42*n53 - n14*n21*n35*n42*n53 + n11*n24*n35*n42*n53 +
466 n15*n22*n31*n44*n53 - n12*n25*n31*n44*n53 - n15*n21*n32*n44*n53 +
467 n11*n25*n32*n44*n53 + n12*n21*n35*n44*n53 - n11*n22*n35*n44*n53 -
468 n14*n22*n31*n45*n53 + n12*n24*n31*n45*n53 + n14*n21*n32*n45*n53 -
469 n11*n24*n32*n45*n53 - n12*n21*n34*n45*n53 + n11*n22*n34*n45*n53 -
470 n15*n23*n32*n41*n54 + n13*n25*n32*n41*n54 + n15*n22*n33*n41*n54 -
471 n12*n25*n33*n41*n54 - n13*n22*n35*n41*n54 + n12*n23*n35*n41*n54 +
472 n15*n23*n31*n42*n54 - n13*n25*n31*n42*n54 - n15*n21*n33*n42*n54 +
473 n11*n25*n33*n42*n54 + n13*n21*n35*n42*n54 - n11*n23*n35*n42*n54 -
474 n15*n22*n31*n43*n54 + n12*n25*n31*n43*n54 + n15*n21*n32*n43*n54 -
475 n11*n25*n32*n43*n54 - n12*n21*n35*n43*n54 + n11*n22*n35*n43*n54 +
476 n13*n22*n31*n45*n54 - n12*n23*n31*n45*n54 - n13*n21*n32*n45*n54 +
477 n11*n23*n32*n45*n54 + n12*n21*n33*n45*n54 - n11*n22*n33*n45*n54 +
478 n14*n23*n32*n41*n55 - n13*n24*n32*n41*n55 - n14*n22*n33*n41*n55 +
479 n12*n24*n33*n41*n55 + n13*n22*n34*n41*n55 - n12*n23*n34*n41*n55 -
480 n14*n23*n31*n42*n55 + n13*n24*n31*n42*n55 + n14*n21*n33*n42*n55 -
481 n11*n24*n33*n42*n55 - n13*n21*n34*n42*n55 + n11*n23*n34*n42*n55 +
482 n14*n22*n31*n43*n55 - n12*n24*n31*n43*n55 - n14*n21*n32*n43*n55 +
483 n11*n24*n32*n43*n55 + n12*n21*n34*n43*n55 - n11*n22*n34*n43*n55 -
484 n13*n22*n31*n44*n55 + n12*n23*n31*n44*n55 + n13*n21*n32*n44*n55 -
485 n11*n23*n32*n44*n55 - n12*n21*n33*n44*n55 + n11*n22*n33*n44*n55);
486
487 if(det == complex<double>(0., 0.)) reject=true;
488
489 // The 1st row of the inverse matrix {(I-iKp)^-1}_0j
490 i[0][0] = (n25*n34*n43*n52 -
491 n24*n35*n43*n52 - n25*n33*n44*n52 + n23*n35*n44*n52 +
492 n24*n33*n45*n52 - n23*n34*n45*n52 - n25*n34*n42*n53 +
493 n24*n35*n42*n53 + n25*n32*n44*n53 - n22*n35*n44*n53 -
494 n24*n32*n45*n53 + n22*n34*n45*n53 + n25*n33*n42*n54 -
495 n23*n35*n42*n54 - n25*n32*n43*n54 + n22*n35*n43*n54 +
496 n23*n32*n45*n54 - n22*n33*n45*n54 - n24*n33*n42*n55 +
497 n23*n34*n42*n55 + n24*n32*n43*n55 - n22*n34*n43*n55 -
498 n23*n32*n44*n55 + n22*n33*n44*n55)/det;
499
500 i[0][1] = (-n15*n34*n43*n52 +
501 n14*n35*n43*n52 + n15*n33*n44*n52 - n13*n35*n44*n52 -
502 n14*n33*n45*n52 + n13*n34*n45*n52 + n15*n34*n42*n53 -
503 n14*n35*n42*n53 - n15*n32*n44*n53 + n12*n35*n44*n53 +
504 n14*n32*n45*n53 - n12*n34*n45*n53 - n15*n33*n42*n54 +
505 n13*n35*n42*n54 + n15*n32*n43*n54 - n12*n35*n43*n54 -
506 n13*n32*n45*n54 + n12*n33*n45*n54 + n14*n33*n42*n55 -
507 n13*n34*n42*n55 - n14*n32*n43*n55 + n12*n34*n43*n55 +
508 n13*n32*n44*n55 - n12*n33*n44*n55)/det;
509
510 i[0][2] = (n15*n24*n43*n52 -
511 n14*n25*n43*n52 - n15*n23*n44*n52 + n13*n25*n44*n52 +
512 n14*n23*n45*n52 - n13*n24*n45*n52 - n15*n24*n42*n53 +
513 n14*n25*n42*n53 + n15*n22*n44*n53 - n12*n25*n44*n53 -
514 n14*n22*n45*n53 + n12*n24*n45*n53 + n15*n23*n42*n54 -
515 n13*n25*n42*n54 - n15*n22*n43*n54 + n12*n25*n43*n54 +
516 n13*n22*n45*n54 - n12*n23*n45*n54 - n14*n23*n42*n55 +
517 n13*n24*n42*n55 + n14*n22*n43*n55 - n12*n24*n43*n55 -
518 n13*n22*n44*n55 + n12*n23*n44*n55)/det;
519
520 i[0][3] = (-n15*n24*n33*n52 +
521 n14*n25*n33*n52 + n15*n23*n34*n52 - n13*n25*n34*n52 -
522 n14*n23*n35*n52 + n13*n24*n35*n52 + n15*n24*n32*n53 -
523 n14*n25*n32*n53 - n15*n22*n34*n53 + n12*n25*n34*n53 +
524 n14*n22*n35*n53 - n12*n24*n35*n53 - n15*n23*n32*n54 +
525 n13*n25*n32*n54 + n15*n22*n33*n54 - n12*n25*n33*n54 -
526 n13*n22*n35*n54 + n12*n23*n35*n54 + n14*n23*n32*n55 -
527 n13*n24*n32*n55 - n14*n22*n33*n55 + n12*n24*n33*n55 +
528 n13*n22*n34*n55 - n12*n23*n34*n55)/det;
529
530 i[0][4] = (n15*n24*n33*n42 -
531 n14*n25*n33*n42 - n15*n23*n34*n42 + n13*n25*n34*n42 +
532 n14*n23*n35*n42 - n13*n24*n35*n42 - n15*n24*n32*n43 +
533 n14*n25*n32*n43 + n15*n22*n34*n43 - n12*n25*n34*n43 -
534 n14*n22*n35*n43 + n12*n24*n35*n43 + n15*n23*n32*n44 -
535 n13*n25*n32*n44 - n15*n22*n33*n44 + n12*n25*n33*n44 +
536 n13*n22*n35*n44 - n12*n23*n35*n44 - n14*n23*n32*n45 +
537 n13*n24*n32*n45 + n14*n22*n33*n45 - n12*n24*n33*n45 -
538 n13*n22*n34*n45 + n12*n23*n34*n45)/det;
539
540
541 //--------------------------------------------------------------------------------------------------------------------------------
542 double s0_prod = -0.07;
543
544 double u1j_re_limit[5][2], u1j_im_limit[5][2] ;
545 u1j_re_limit[0][0] = 0. ; u1j_re_limit[0][1] = 1. ;
546 u1j_re_limit[1][0] = -0.29 ; u1j_re_limit[1][1] = 0.12 ;
547 u1j_re_limit[2][0] = -0.17 ; u1j_re_limit[2][1] = 0.065 ;
548 u1j_re_limit[3][0] = -0.66 ; u1j_re_limit[3][1] = 0.1 ;
549 u1j_re_limit[4][0] = -1.36 ; u1j_re_limit[4][1] = 0.18 ;
550
551 u1j_im_limit[0][0] = -0.58 ; u1j_im_limit[0][1] = 0.58 ;
552 u1j_im_limit[1][0] = 0.00 ; u1j_im_limit[1][1] = 0.28 ;
553 u1j_im_limit[2][0] = -0.135 ; u1j_im_limit[2][1] = 0.10 ;
554 u1j_im_limit[3][0] = -0.13 ; u1j_im_limit[3][1] = 0.40 ;
555 u1j_im_limit[4][0] = -0.36 ; u1j_im_limit[4][1] = 0.80 ;
556
557 //for(int kk=0 ; kk<5 ; kk++)
558 //{
559 // cout<<"beta["<<kk<<"]: "<<abs(beta[kk])<<", \t"<<arg(beta[kk])*EvtConst::radToDegrees<<endl;
560 // cout<<"fprod["<<kk<<"]: "<<abs(fprod[kk])<<", \t"<<arg(fprod[kk])*EvtConst::radToDegrees<<endl;
561 //}
562
563 complex<double> value0(0., 0.) ;
564 complex<double> value1(0., 0.) ;
565
566 for(int l=0;l<5;l++)
567 {
568 double u1j_re = real(i[0][l]) ;
569 double u1j_im = imag(i[0][l]) ;
570 if(u1j_re==0. || u1j_im==0.) reject=true;
571 if(u1j_re < u1j_re_limit[l][0] || u1j_re > u1j_re_limit[l][1] || u1j_im < u1j_im_limit[l][0] || u1j_im > u1j_im_limit[l][1]) reject=true;
572
573 for(int pole_index=0;pole_index<5;pole_index++)
574 {
575 complex<double> A = beta[pole_index]*g[pole_index][l];
576 value0 += (i[0][l]*A)/(ma[pole_index]*ma[pole_index]-s) ;
577 //cout<<"value0["<<l<<"]["<<pole_index<<"] = "<<value0<<endl;
578 }
579 }
580
581 value1 += i[0][0]*fprod[0];
582 value1 += i[0][1]*fprod[1];
583 value1 += i[0][2]*fprod[2];
584 value1 += i[0][3]*fprod[3];
585 value1 += i[0][4]*fprod[4];
586 value1 *= (1.-s0_prod)/(s-s0_prod) ;
587 //cout<<"value1 = "<<value1<<endl;
588
589 if(reject==true) return complex<double>(9999., 9999.) ;
590 else return (value0+value1) ;
591 //return (value0+value1) ;
592
593}
594
595complex<double> D0ToKLpipi::amplitude_LASS(vector<double> p_k0l, vector<double> p_pip, vector<double> p_pim, string reso, double A_r, double Phi_r) {
596
597 double mR = 1.425 ;
598 double gammaR = 0.27 ;
599 double mab2=0.0;
600 HepLorentzVector _p_k0l(p_k0l[0],p_k0l[1],p_k0l[2],p_k0l[3]);
601 HepLorentzVector _p_pip(p_pip[0],p_pip[1],p_pip[2],p_pip[3]);
602 HepLorentzVector _p_pim(p_pim[0],p_pim[1],p_pim[2],p_pim[3]);
603 if (reso == "k0lpim") mab2 = pow((_p_k0l + _p_pim).m(),2);
604 else if(reso == "k0lpip") mab2 = pow((_p_k0l + _p_pip).m(),2);
605
606 double s = mab2;
607
608 const double mD0 = 1.86483;
609 const double mKl = 0.49761;
610 const double mPi = 0.13957;
611
612 double _a = 0.113 ;
613 double _r = -33.8 ;
614 double _R = 1.0 ;
615 double _F = 0.96 ;
616 double _phiR = -1.9146 ;
617 double _phiF = 0.0017 ;
618
619 double fR=1.0; // K*0(1430) has spin zero
620 int power=1; // Power is 1 for spin zero
621
622 double mAB = sqrt(mab2) ; // (_p4_d1+_p4_d2).mass();
623
624 double mA=mKl; // _p4_d1.mass();
625 double mB=mPi; // _p4_d2.mass();
626 double mC=mPi;
627 double mD = mD0;
628
629 double pAB=sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
630 double q=pAB;
631
632 double pR=sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
633 double q0=pR;
634
635 // Running width.
636 double g = gammaR*pow(q/q0,power)*(mR/mAB)*fR*fR;
637
638 complex<double> propagator_relativistic_BreitWigner = 1./(mR*mR - mAB*mAB - complex<double>(0.,mR*g));
639
640 // Non-resonant phase shift
641 double cot_deltaF = 1.0/(_a*q) + 0.5*_r*q;
642 double qcot_deltaF = 1.0/_a + 0.5*_r*q*q;
643
644 // Compute resonant part
645 complex<double> expi2deltaF = complex<double>(qcot_deltaF, q)/ complex<double>(qcot_deltaF, -q);
646
647 complex<double> resonant_term_T = _R * complex<double>(cos(_phiR + 2 * _phiF), sin(_phiR + 2 * _phiF)) * propagator_relativistic_BreitWigner * mR * gammaR * mR / q0 * expi2deltaF;
648
649 // Compute non-resonant part
650 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);
651
652 // Add non-resonant and resonant terms
653 complex<double> LASS_contribution = non_resonant_term_F + resonant_term_T;
654
655 return complex<double>(A_r*cos(Phi_r), A_r*sin(Phi_r)) * LASS_contribution;
656
657}
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
D0ToKLpipi::Amp_PFT
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition: D0ToKLpipi.cxx:100
D0ToKLpipi::~D0ToKLpipi
virtual ~D0ToKLpipi()
Definition: D0ToKLpipi.cxx:24
D0ToKLpipi::init
void init()
Definition: D0ToKLpipi.cxx:26