BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSpipi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of models developed at BES collaboration
5// based on the EvtGen framework. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/BesCopyright
9// Copyright (A) 2006 Ping Rong-Gang @IHEP
10//
11// Module: EvtD0ToKSpipi.cc
12//
13// Description: Model provided by user, see BAM-600 (arXiv:2212.09048)
14//
15// Modification history:
16//
17// Liaoyuan Dong Aug. 11, 2022 Module created
18//
19//------------------------------------------------------------------------
20#include "EvtGenBase/EvtPatches.hh"
21#include "EvtGenBase/EvtParticle.hh"
22#include "EvtGenBase/EvtGenKine.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtVector4R.hh"
25#include "EvtGenBase/EvtReport.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include "EvtGenModels/EvtD0ToKSpipi.hh"
29#include <stdlib.h>
30using namespace std;
31
32EvtD0ToKSpipi::~EvtD0ToKSpipi() {}
33
34void EvtD0ToKSpipi::getName(std::string& model_name){
35 model_name="D0ToKSpipi";
36}
37
38EvtDecayBase* EvtD0ToKSpipi::clone(){
39 return new EvtD0ToKSpipi;
40}
41
42void EvtD0ToKSpipi::init(){
43 checkNArg(1);
44 checkNDaug(3);
45 tagmode = getArg(0); // Specify the tag mode to be used. Key: {0:no-specified, 1:Kpi, 2:Kpipi0, 3:K3pi}
46 checkSpinParent(EvtSpinType::SCALAR);
47 //std::cout << "Initializing EvtD0ToKSpipi" << std::endl;
48
49 _nd = 3;
50 tan2thetaC = (0.22650*0.22650)/(1.-(0.22650*0.22650)) ; //sin(theta_C) = 0.22650 +/- 0.00048
51 pi180inv = 1.0/EvtConst::radToDegrees;
52
53 mass_R[0]= 0.77526; width_R[0]= 0.14740; spin_R[0]= 1; ar[0]= 1; phir[0]= 0;
54 mass_R[1]= 0.78266; width_R[1]= 0.00868; spin_R[1]= 1; ar[1]= 0.037606; phir[1]= 109.677;
55 mass_R[2]= 1.27550; width_R[2]= 0.18670; spin_R[2]= 2; ar[2]= 1.54909; phir[2]= -42.7425;
56 mass_R[3]= 1.46500; width_R[3]= 0.40000; spin_R[3]= 1; ar[3]= 3.70735; phir[3]= 103.644;
57 mass_R[4]= 0.89167; width_R[4]= 0.0514; spin_R[4]= 1; ar[4]= 1.86093; phir[4]= 136.529;
58 mass_R[5]= 1.42730; width_R[5]= 0.10000; spin_R[5]= 2; ar[5]= 1.74288; phir[5]= -48.0968;
59 mass_R[6]= 1.71800; width_R[6]= 0.3220; spin_R[6]= 1; ar[6]= 3.31; phir[6]= -118.2;
60 mass_R[7]= 1.41400; width_R[7]= 0.2320; spin_R[7]= 1; ar[7]= 0.171672; phir[7]= -68.41;
61 mass_R[8]= 0.89167; width_R[8]= 0.0514; spin_R[8]= 1; ar[8]= 0.164; phir[8]= -42.2;
62 mass_R[9]= 1.42730; width_R[9]= 0.1000; spin_R[9]= 2; ar[9]= 0.1; phir[9]= -89.6;
63 mass_R[10]= 1.41400; width_R[10]= 0.2320; spin_R[10]= 1; ar[10]= 0.21; phir[10]= 150.2;
64 mass_R[11]= 1.42500; width_R[11]= 0.2700; spin_R[11]= 1; ar[11]= 2.78276; phir[11]= 97.9608;
65 mass_R[12]= 1.42500; width_R[12]= 0.2700; spin_R[12]= 1; ar[12]= 0.11; phir[12]= 162.3;
66
67 beta[0] = EvtComplex( 0.255303*cos( 47.8861 *pi180inv), 0.255303*sin( 47.8861 *pi180inv));
68 beta[1] = EvtComplex(13.4446 *cos( -5.11127*pi180inv), 13.4446 *sin( -5.11127*pi180inv));
69 beta[2] = EvtComplex(38.8496 *cos(-30.06 *pi180inv), 38.8496 *sin(-30.06 *pi180inv));
70 beta[3] = EvtComplex(13.1086 *cos(-81.4148 *pi180inv), 13.1086 *sin(-81.4148 *pi180inv));
71 beta[4] = EvtComplex(0., 0.);
72
73 fprod[0] = EvtComplex(5.08049*cos(-182.312*pi180inv), 5.08049*sin(-182.312*pi180inv));
74 fprod[1] = EvtComplex(17.2388*cos(-219.209*pi180inv), 17.2388*sin(-219.209*pi180inv));
75 fprod[2] = EvtComplex(19.0145*cos(-76.9884*pi180inv), 19.0145*sin(-76.9884*pi180inv));
76 fprod[3] = EvtComplex(11.9875*cos(-190.502*pi180inv), 11.9875*sin(-190.502*pi180inv));
77 fprod[4] = EvtComplex(0., 0.);
78
79 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;
80 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;
81 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;
82 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;
83 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;
84
85 // Hadronic parameters for tag modes: 0=no-specified, 1=Kpi, 2=Kpipi0, 3=K3pi
86 rd[0] = 0.0;
87 rd[1] = 0.0586;
88 rd[2] = 0.0440;
89 rd[3] = 0.0546;
90 deltad[0] = 0.0;
91 deltad[1] = 194.7*pi180inv;
92 deltad[2] = 196.0*pi180inv;
93 deltad[3] = 167.0*pi180inv;
94 Rf[0] = 0.0;
95 Rf[1] = 1.0;
96 Rf[2] = 0.78;
97 Rf[3] = 0.52;
98
99}
100
101void EvtD0ToKSpipi::initProbMax(){
102 double ProbMax = 12500.0;
103 if (tagmode==1) ProbMax = 12500.0;
104 else if (tagmode==2) ProbMax = 12000.0;
105 else if (tagmode==3) ProbMax = 12100.0;
106 else ProbMax = 12000.0;
107 setProbMax(ProbMax);
108}
109
110void EvtD0ToKSpipi::decay( EvtParticle *p ){
111/*
112 double maxprob = 0.0;
113 for(int ir=0;ir<=60000000;ir++){
114 p->initializePhaseSpace(getNDaug(),getDaugs());
115 for(int i=0; i<_nd; i++){
116 _p4Lab[i]=p->getDaug(i)->getP4Lab();
117 _p4CM[i]=p->getDaug(i)->getP4();
118 }
119 double Prob = AmplitudeSquare();
120 if(Prob>maxprob) {
121 maxprob=Prob;
122 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
123 }
124 }
125 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
126*/
127 p->initializePhaseSpace(getNDaug(),getDaugs());
128 for(int i=0; i<_nd; i++){
129 _p4Lab[i]=p->getDaug(i)->getP4Lab();
130 _p4CM[i]=p->getDaug(i)->getP4();
131 }
132 double prob = AmplitudeSquare();
133 setProb(prob);
134 return;
135}
136
137double EvtD0ToKSpipi::AmplitudeSquare() {
138 EvtVector4R k0l = GetDaugMomLab(0);
139 EvtVector4R pip = GetDaugMomLab(1);
140 EvtVector4R pim = GetDaugMomLab(2);
141 EvtVector4R pD = k0l + pip + pim ;
142
143 double total_prob, Interference=0.;
144
145 // Breit-Wigner lineshapes
146 EvtComplex DK2piRes0 = Resonance2(pD, pip, pim, ar[0], phir[0], width_R[0], mass_R[0], spin_R[0]); //ar, phir, width, mass, spin Rho770
147 EvtComplex DK2piRes1 = Resonance2(pD, pip, pim, ar[1], phir[1], width_R[1], mass_R[1], spin_R[1]); //ar, phir, width, mass, spin Omega782
148 EvtComplex DK2piRes2 = Resonance2(pD, pip, pim, ar[2], phir[2], width_R[2], mass_R[2], spin_R[2]); //ar, phir, width, mass, spin ftwo1270
149 EvtComplex DK2piRes3 = Resonance2(pD, pip, pim, ar[3], phir[3], width_R[3], mass_R[3], spin_R[3]); //ar, phir, width, mass, spin Rho1450
150 EvtComplex DK2piRes4 = Resonance2(pD, k0l, pim, ar[4], phir[4], width_R[4], mass_R[4], spin_R[4]); //ar, phir, width, mass, spin Kstar892-
151 EvtComplex DK2piRes5 = Resonance2(pD, k0l, pim, ar[5], phir[5], width_R[5], mass_R[5], spin_R[5]); //ar, phir, width, mass, spin K2star1430-
152 EvtComplex DK2piRes6 = Resonance2(pD, k0l, pim, ar[6], phir[6], width_R[6], mass_R[6], spin_R[6]); //ar, phir, width, mass, spin Kstar1680-
153 EvtComplex DK2piRes7 = Resonance2(pD, k0l, pim, ar[7], phir[7], width_R[7], mass_R[7], spin_R[7]); //ar, phir, width, mass, spin Kstar1410-
154 EvtComplex DK2piRes8 = Resonance2(pD, k0l, pip, ar[8], phir[8], width_R[8], mass_R[8], spin_R[8]); //ar, phir, width, mass, spin Kstar892+
155 EvtComplex DK2piRes9 = Resonance2(pD, k0l, pip, ar[9], phir[9], width_R[9], mass_R[9], spin_R[9]); //ar, phir, width, mass, spin K2star1430+
156 EvtComplex DK2piRes10 = Resonance2(pD, k0l, pip, ar[10], phir[10], width_R[10], mass_R[10], spin_R[10]); //ar, phir, width, mass, spin Kstar1410+
157
158 // K-matrix for pipi S-wave
159 EvtComplex pipi_s_wave = K_matrix(pip, pim);
160 if(pipi_s_wave == EvtComplex(9999., 9999.)) return 1e-20;
161
162 // LASS parametrization for Kpi S-wave
163 EvtComplex kpi_s_wave = amplitude_LASS(k0l, pip, pim, "k0spim", ar[11], phir[11]*pi180inv);
164 EvtComplex dcs_kpi_s_wave = amplitude_LASS(k0l, pip, pim, "k0spip", ar[12], phir[12]*pi180inv);
165
166 // Coherent sum for pure-flavor-tagged amplitudes (PFT)
167 EvtComplex TOT_PFT_AMP = DK2piRes0+ DK2piRes1+ DK2piRes2+ DK2piRes3+ DK2piRes4+ DK2piRes5+ DK2piRes6+ DK2piRes7+ DK2piRes8+ DK2piRes9+ DK2piRes10+ pipi_s_wave + kpi_s_wave+ dcs_kpi_s_wave ;
168
169 if (tagmode==1||tagmode==2||tagmode==3) {
170 // Amplitudes for DCS tag decay modes (wherever required)
171 EvtComplex DK2piRes4_dcs = Resonance2(pD, k0l, pim, ar[8], phir[8], width_R[4], mass_R[4], spin_R[4]); //Kstar892minus
172 EvtComplex DK2piRes5_dcs = Resonance2(pD, k0l, pim, ar[9], phir[9], width_R[5], mass_R[5], spin_R[5]); //K2star1430minus
173 EvtComplex DK2piRes6_dcs = Resonance2(pD, k0l, pip, ar[6], phir[6], width_R[6], mass_R[6], spin_R[6]); //Kstar1680plus
174 EvtComplex DK2piRes7_dcs = Resonance2(pD, k0l, pim, ar[10], phir[10], width_R[7], mass_R[7], spin_R[7]); //Kstar1410minus
175 EvtComplex DK2piRes8_dcs = Resonance2(pD, k0l, pip, ar[4], phir[4], width_R[8], mass_R[8], spin_R[8]); //Kstar892plus
176 EvtComplex DK2piRes9_dcs = Resonance2(pD, k0l, pip, ar[5], phir[5], width_R[9], mass_R[9], spin_R[9]); //K2star1430plus
177 EvtComplex DK2piRes10_dcs = Resonance2(pD, k0l, pip, ar[7], phir[7], width_R[10], mass_R[10], spin_R[10]); //Kstar1410plus
178
179 EvtComplex kpi_s_wave_dcs = amplitude_LASS(k0l, pip, pim, "k0spip", ar[11], phir[11]*pi180inv) ;
180 EvtComplex dcs_kpi_s_wave_dcs = amplitude_LASS(k0l, pip, pim, "k0spim", ar[12], phir[12]*pi180inv) ;
181
182 // Coherent sum of amplitudes in case of tag side DCS decay modes
183 EvtComplex TOT_DCS_AMP = DK2piRes0+ DK2piRes1+ DK2piRes2+ DK2piRes3+ DK2piRes4_dcs+ DK2piRes5_dcs+ DK2piRes6_dcs+ DK2piRes7_dcs+ DK2piRes8_dcs+ DK2piRes9_dcs+ DK2piRes10_dcs+ pipi_s_wave+ kpi_s_wave_dcs + dcs_kpi_s_wave_dcs ;
184
185 // Interference between PFT and DCS
186 Interference = real( EvtComplex(1.0*cos(deltad[tagmode]), -1.0*sin(deltad[tagmode]))*TOT_DCS_AMP*conj(TOT_PFT_AMP) + EvtComplex(1.0*cos(deltad[tagmode]), 1.0*sin(deltad[tagmode]))*TOT_PFT_AMP*conj(TOT_DCS_AMP) );
187
188 // Coherently added PFT and DCS probabilities
189 total_prob = abs2(TOT_PFT_AMP) + rd[tagmode]*rd[tagmode]*abs2(TOT_DCS_AMP) - rd[tagmode]*Rf[tagmode]*Interference ;
190 } else {
191 total_prob = abs2(TOT_PFT_AMP);
192 }
193
194 return (total_prob <= 0) ? 1e-20 : total_prob;
195
196}
197
198EvtComplex EvtD0ToKSpipi::Resonance2(EvtVector4R p4_p, EvtVector4R p4_d1, const EvtVector4R p4_d2, double mag, double theta, double gamma, double bwm, int spin) {
199
200 EvtComplex ampl;
201 EvtVector4R p4_d3 = p4_p - p4_d1 - p4_d2;
202
203 double mAB= (p4_d1 + p4_d2).mass();
204 double mBC= (p4_d2 + p4_d3).mass();
205 double mAC= (p4_d1 + p4_d3).mass();
206 double mA = p4_d1.mass();
207 double mB = p4_d2.mass();
208 double mD = p4_p.mass();
209 double mC = p4_d3.mass();
210
211 double mR = bwm;
212 double gammaR = gamma;
213 double pAB = sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
214 double pR = sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
215
216 double pD= (((mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0) - mR*mR*mC*mC)/(mD*mD);
217 if ( pD>0 ) { pD = sqrt(pD); }
218 else { pD = 0;}
219 double pDAB=sqrt( (((mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0) - mAB*mAB*mC*mC)/(mD*mD));
220 double fR = 1;
221 double fD = 1;
222 int power = 0;
223 switch (spin) {
224 case 0:
225 fR = 1.0;
226 fD = 1.0;
227 power = 1;
228 break;
229 case 1:
230 fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
231 fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
232 power = 3;
233 break;
234 case 2:
235 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)) );
236 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)) );
237 power = 5;
238 break;
239 default:
240 report(INFO,"EvtGen") << "Incorrect spin in EvtD0ToKSpipi::EvtResonance2.cc\n";
241 }
242
243 double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
244 switch (spin) {
245 case 0:
246 ampl=mag*EvtComplex(cos(theta*pi180inv),sin(theta*pi180inv))*fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB));
247 break;
248 case 1:
249 ampl=mag*EvtComplex(cos(theta*pi180inv),sin(theta*pi180inv))*
250 (fR*fD*(mAC*mAC-mBC*mBC+((mD*mD-mC*mC)*(mB*mB-mA*mA)/(mAB*mAB)))/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB)));
251 break;
252 case 2:
253 ampl=mag*EvtComplex(cos(theta*pi180inv),sin(theta*pi180inv))*
254 (fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB)))*
255 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mAB*mAB)),2)-
256 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mAB, 2))*
257 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mAB,2)));
258 break;
259 default:
260 report(INFO,"EvtGen") << "Incorrect spin in EvtD0ToKSpipi::Resonance2.cc\n";
261 }
262
263 return ampl;
264}
265
266EvtComplex EvtD0ToKSpipi::K_matrix(EvtVector4R p_pip, EvtVector4R p_pim) {
267 const double mD0 = 1.86483;
268 const double mKl = 0.49761;
269 const double mPi = 0.13957;
270 bool reject = false;
271
272 double mAB = (p_pip + p_pim).mass() ;
273 double s = mAB*mAB;
274
275 EvtComplex 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;
276 double rho1sq,rho2sq, rho4sq,rho5sq;
277 EvtComplex rho1, rho2,rho3,rho4, rho5;
278 EvtComplex rho[5];
279 EvtComplex pole,SVT,Adler;
280 EvtComplex det;
281 EvtComplex i[5][5];
282 double f[5][5];
283
284 const double mpi = 0.13957;
285 const double mK = 0.493677;
286 const double meta = 0.54775;
287 const double metap = 0.95778;
288
289 // Init matrices and vectors with zeros
290 EvtComplex K[5][5];
291 for(int k=0;k<5;k++) {
292 for(int l=0;l<5;l++) {
293 i[k][l]=EvtComplex(0.,0.);
294 K[k][l]=EvtComplex(0.,0.);
295 f[k][l]=0.;
296 }
297 rho[k]=0.;
298 }
299
300 // Fill scattering data values
301 double s_scatt = -3.92637;
302 double sa = 1.0;
303 double sa_0 = -0.15;
304
305 // f_scattering (At least one of the two channels must be pi+pi-)
306 f[0][0] = 0.23399;
307 f[0][1] = 0.15044;
308 f[0][2] = -0.20545;
309 f[0][3] = 0.32825;
310 f[0][4] = 0.35412;
311
312 f[1][0] = f[0][1];
313 f[2][0] = f[0][2];
314 f[3][0] = f[0][3];
315 f[4][0] = f[0][4];
316
317 // Compute phase space factors
318 // rho_0
319 rho1sq=(1.0-(pow((mpi+mpi),2)/s));
320 if(rho1sq >=0.) rho1=EvtComplex(sqrt(rho1sq),0.);
321 else rho1=EvtComplex(0.,sqrt(-rho1sq));
322 rho[0]=rho1;
323
324 // rho_1
325 rho2sq=(1.0-(pow((mK+mK),2)/s));
326 if(rho2sq >=0.) rho2=EvtComplex(sqrt(rho2sq),0.);
327 else rho2=EvtComplex(0.,sqrt(-rho2sq));
328 rho[1]=rho2;
329
330 // rho_2
331 rho3=EvtComplex(0.,0.);
332 if(s<=1) {
333 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;
334 double cont32=sqrt(1.0-(16.0*mpi*mpi));
335 rho3=EvtComplex(cont32*real,0.);
336 }
337 else rho3=EvtComplex(sqrt(1.0-(16.0*mpi*mpi/s)),0.);
338 rho[2]=rho3;
339
340 // rho_3
341 rho4sq=(1.0-(pow((meta+meta),2)/s));
342 if(rho4sq>=0.) rho4=EvtComplex(sqrt(rho4sq),0.);
343 else rho4=EvtComplex(0.,sqrt(-rho4sq));
344 rho[3]=rho4;
345
346 // rho_4
347 rho5sq=(1.0-(pow((meta+metap),2)/s));
348 if(rho5sq >=0.) rho5=EvtComplex(sqrt(rho5sq),0.);
349 else rho5=EvtComplex(0.,sqrt(-rho5sq));
350 rho[4]=rho5;
351
352 // Sum over the poles [Intermediate channel(k) -> pole(pole_index) -> final channel(l)]
353 for(int k=0;k<5;k++) {
354 for(int l=0;l<5;l++) {
355 for (int pole_index=0;pole_index<5;pole_index++) {
356 double A=g[pole_index][k]*g[pole_index][l];
357 double B=ma[pole_index]*ma[pole_index]-s;
358 K[k][l]=K[k][l]+EvtComplex(A/B,0.);
359 }
360 }
361 }
362
363 // Direct scattering term [k -> l]
364 for(int k=0;k<5;k++) {
365 for(int l=0;l<5;l++) {
366 double C=f[k][l]*(1.0-s_scatt);
367 double D=(s-s_scatt);
368 K[k][l]=K[k][l]+EvtComplex(C/D,0.);
369 }
370 }
371
372 // Multiplying the "Adler zero" term
373 for(int k=0;k<5;k++) {
374 for(int l=0;l<5;l++) {
375 double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
376 double F=(s-sa_0);
377 K[k][l]=K[k][l]*EvtComplex(E/F,0.);
378 }
379 }
380
381 // (1 - i rho K)_ij
382 n11=EvtComplex(1.,0.)-EvtComplex(0.,1.)*K[0][0]*rho[0];
383 n12=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[0][1]*rho[1];
384 n13=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[0][2]*rho[2];
385 n14=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[0][3]*rho[3];
386 n15=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[0][4]*rho[4];
387
388 n21=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[1][0]*rho[0];
389 n22=EvtComplex(1.,0.)-EvtComplex(0.,1.)*K[1][1]*rho[1];
390 n23=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[1][2]*rho[2];
391 n24=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[1][3]*rho[3];
392 n25=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[1][4]*rho[4];
393
394 n31=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[2][0]*rho[0];
395 n32=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[2][1]*rho[1];
396 n33=EvtComplex(1.,0.)-EvtComplex(0.,1.)*K[2][2]*rho[2];
397 n34=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[2][3]*rho[3];
398 n35=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[2][4]*rho[4];
399
400 n41=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[3][0]*rho[0];
401 n42=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[3][1]*rho[1];
402 n43=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[3][2]*rho[2];
403 n44=EvtComplex(1.,0.)-EvtComplex(0.,1.)*K[3][3]*rho[3];
404 n45=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[3][4]*rho[4];
405
406 n51=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[4][0]*rho[0];
407 n52=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[4][1]*rho[1];
408 n53=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[4][2]*rho[2];
409 n54=EvtComplex(0.,0.)-EvtComplex(0.,1.)*K[4][3]*rho[3];
410 n55=EvtComplex(1.,0.)-EvtComplex(0.,1.)*K[4][4]*rho[4];
411
412 // Compute the determinant for inverse [Looks horrible but TMatrixT does not support complex quantities; python bindings may help, working on it.]
413 det = (n15*n24*n33*n42*n51 - n14*n25*n33*n42*n51 - n15*n23*n34*n42*n51 +
414 n13*n25*n34*n42*n51 + n14*n23*n35*n42*n51 - n13*n24*n35*n42*n51 -
415 n15*n24*n32*n43*n51 + n14*n25*n32*n43*n51 + n15*n22*n34*n43*n51 -
416 n12*n25*n34*n43*n51 - n14*n22*n35*n43*n51 + n12*n24*n35*n43*n51 +
417 n15*n23*n32*n44*n51 - n13*n25*n32*n44*n51 - n15*n22*n33*n44*n51 +
418 n12*n25*n33*n44*n51 + n13*n22*n35*n44*n51 - n12*n23*n35*n44*n51 -
419 n14*n23*n32*n45*n51 + n13*n24*n32*n45*n51 + n14*n22*n33*n45*n51 -
420 n12*n24*n33*n45*n51 - n13*n22*n34*n45*n51 + n12*n23*n34*n45*n51 -
421 n15*n24*n33*n41*n52 + n14*n25*n33*n41*n52 + n15*n23*n34*n41*n52 -
422 n13*n25*n34*n41*n52 - n14*n23*n35*n41*n52 + n13*n24*n35*n41*n52 +
423 n15*n24*n31*n43*n52 - n14*n25*n31*n43*n52 - n15*n21*n34*n43*n52 +
424 n11*n25*n34*n43*n52 + n14*n21*n35*n43*n52 - n11*n24*n35*n43*n52 -
425 n15*n23*n31*n44*n52 + n13*n25*n31*n44*n52 + n15*n21*n33*n44*n52 -
426 n11*n25*n33*n44*n52 - n13*n21*n35*n44*n52 + n11*n23*n35*n44*n52 +
427 n14*n23*n31*n45*n52 - n13*n24*n31*n45*n52 - n14*n21*n33*n45*n52 +
428 n11*n24*n33*n45*n52 + n13*n21*n34*n45*n52 - n11*n23*n34*n45*n52 +
429 n15*n24*n32*n41*n53 - n14*n25*n32*n41*n53 - n15*n22*n34*n41*n53 +
430 n12*n25*n34*n41*n53 + n14*n22*n35*n41*n53 - n12*n24*n35*n41*n53 -
431 n15*n24*n31*n42*n53 + n14*n25*n31*n42*n53 + n15*n21*n34*n42*n53 -
432 n11*n25*n34*n42*n53 - n14*n21*n35*n42*n53 + n11*n24*n35*n42*n53 +
433 n15*n22*n31*n44*n53 - n12*n25*n31*n44*n53 - n15*n21*n32*n44*n53 +
434 n11*n25*n32*n44*n53 + n12*n21*n35*n44*n53 - n11*n22*n35*n44*n53 -
435 n14*n22*n31*n45*n53 + n12*n24*n31*n45*n53 + n14*n21*n32*n45*n53 -
436 n11*n24*n32*n45*n53 - n12*n21*n34*n45*n53 + n11*n22*n34*n45*n53 -
437 n15*n23*n32*n41*n54 + n13*n25*n32*n41*n54 + n15*n22*n33*n41*n54 -
438 n12*n25*n33*n41*n54 - n13*n22*n35*n41*n54 + n12*n23*n35*n41*n54 +
439 n15*n23*n31*n42*n54 - n13*n25*n31*n42*n54 - n15*n21*n33*n42*n54 +
440 n11*n25*n33*n42*n54 + n13*n21*n35*n42*n54 - n11*n23*n35*n42*n54 -
441 n15*n22*n31*n43*n54 + n12*n25*n31*n43*n54 + n15*n21*n32*n43*n54 -
442 n11*n25*n32*n43*n54 - n12*n21*n35*n43*n54 + n11*n22*n35*n43*n54 +
443 n13*n22*n31*n45*n54 - n12*n23*n31*n45*n54 - n13*n21*n32*n45*n54 +
444 n11*n23*n32*n45*n54 + n12*n21*n33*n45*n54 - n11*n22*n33*n45*n54 +
445 n14*n23*n32*n41*n55 - n13*n24*n32*n41*n55 - n14*n22*n33*n41*n55 +
446 n12*n24*n33*n41*n55 + n13*n22*n34*n41*n55 - n12*n23*n34*n41*n55 -
447 n14*n23*n31*n42*n55 + n13*n24*n31*n42*n55 + n14*n21*n33*n42*n55 -
448 n11*n24*n33*n42*n55 - n13*n21*n34*n42*n55 + n11*n23*n34*n42*n55 +
449 n14*n22*n31*n43*n55 - n12*n24*n31*n43*n55 - n14*n21*n32*n43*n55 +
450 n11*n24*n32*n43*n55 + n12*n21*n34*n43*n55 - n11*n22*n34*n43*n55 -
451 n13*n22*n31*n44*n55 + n12*n23*n31*n44*n55 + n13*n21*n32*n44*n55 -
452 n11*n23*n32*n44*n55 - n12*n21*n33*n44*n55 + n11*n22*n33*n44*n55);
453
454 if(det == EvtComplex(0., 0.)) reject=true;
455
456 // The 1st row of the inverse matrix [(1-i\rhoK)^-1]_0j
457 i[0][0] = ( n25*n34*n43*n52 -
458 n24*n35*n43*n52 - n25*n33*n44*n52 + n23*n35*n44*n52 +
459 n24*n33*n45*n52 - n23*n34*n45*n52 - n25*n34*n42*n53 +
460 n24*n35*n42*n53 + n25*n32*n44*n53 - n22*n35*n44*n53 -
461 n24*n32*n45*n53 + n22*n34*n45*n53 + n25*n33*n42*n54 -
462 n23*n35*n42*n54 - n25*n32*n43*n54 + n22*n35*n43*n54 +
463 n23*n32*n45*n54 - n22*n33*n45*n54 - n24*n33*n42*n55 +
464 n23*n34*n42*n55 + n24*n32*n43*n55 - n22*n34*n43*n55 -
465 n23*n32*n44*n55 + n22*n33*n44*n55)/det;
466
467 i[0][1] = ( -n15*n34*n43*n52 +
468 n14*n35*n43*n52 + n15*n33*n44*n52 - n13*n35*n44*n52 -
469 n14*n33*n45*n52 + n13*n34*n45*n52 + n15*n34*n42*n53 -
470 n14*n35*n42*n53 - n15*n32*n44*n53 + n12*n35*n44*n53 +
471 n14*n32*n45*n53 - n12*n34*n45*n53 - n15*n33*n42*n54 +
472 n13*n35*n42*n54 + n15*n32*n43*n54 - n12*n35*n43*n54 -
473 n13*n32*n45*n54 + n12*n33*n45*n54 + n14*n33*n42*n55 -
474 n13*n34*n42*n55 - n14*n32*n43*n55 + n12*n34*n43*n55 +
475 n13*n32*n44*n55 - n12*n33*n44*n55)/det;
476
477 i[0][2] = ( n15*n24*n43*n52 -
478 n14*n25*n43*n52 - n15*n23*n44*n52 + n13*n25*n44*n52 +
479 n14*n23*n45*n52 - n13*n24*n45*n52 - n15*n24*n42*n53 +
480 n14*n25*n42*n53 + n15*n22*n44*n53 - n12*n25*n44*n53 -
481 n14*n22*n45*n53 + n12*n24*n45*n53 + n15*n23*n42*n54 -
482 n13*n25*n42*n54 - n15*n22*n43*n54 + n12*n25*n43*n54 +
483 n13*n22*n45*n54 - n12*n23*n45*n54 - n14*n23*n42*n55 +
484 n13*n24*n42*n55 + n14*n22*n43*n55 - n12*n24*n43*n55 -
485 n13*n22*n44*n55 + n12*n23*n44*n55)/det;
486
487 i[0][3] = ( -n15*n24*n33*n52 +
488 n14*n25*n33*n52 + n15*n23*n34*n52 - n13*n25*n34*n52 -
489 n14*n23*n35*n52 + n13*n24*n35*n52 + n15*n24*n32*n53 -
490 n14*n25*n32*n53 - n15*n22*n34*n53 + n12*n25*n34*n53 +
491 n14*n22*n35*n53 - n12*n24*n35*n53 - n15*n23*n32*n54 +
492 n13*n25*n32*n54 + n15*n22*n33*n54 - n12*n25*n33*n54 -
493 n13*n22*n35*n54 + n12*n23*n35*n54 + n14*n23*n32*n55 -
494 n13*n24*n32*n55 - n14*n22*n33*n55 + n12*n24*n33*n55 +
495 n13*n22*n34*n55 - n12*n23*n34*n55)/det;
496
497 i[0][4] = ( n15*n24*n33*n42 -
498 n14*n25*n33*n42 - n15*n23*n34*n42 + n13*n25*n34*n42 +
499 n14*n23*n35*n42 - n13*n24*n35*n42 - n15*n24*n32*n43 +
500 n14*n25*n32*n43 + n15*n22*n34*n43 - n12*n25*n34*n43 -
501 n14*n22*n35*n43 + n12*n24*n35*n43 + n15*n23*n32*n44 -
502 n13*n25*n32*n44 - n15*n22*n33*n44 + n12*n25*n33*n44 +
503 n13*n22*n35*n44 - n12*n23*n35*n44 - n14*n23*n32*n45 +
504 n13*n24*n32*n45 + n14*n22*n33*n45 - n12*n24*n33*n45 -
505 n13*n22*n34*n45 + n12*n23*n34*n45)/det;
506
507 double s0_prod = -0.07;
508
509 EvtComplex value0(0., 0.);
510 EvtComplex value1(0., 0.);
511
512 // [(1-i\rhoK)^-1]_0j*P_j {P_j: Production vector}
513 for(int k=0;k<5;k++) {
514 double u1j_re = real(i[0][k]);
515 double u1j_im = imag(i[0][k]);
516 if(u1j_re==0. || u1j_im==0.) reject=true;
517
518 // Initial state to K-matrix pole couplings * Pole to intermediate channels coupling
519 for(int pole_index=0;pole_index<5;pole_index++) {
520 EvtComplex A = beta[pole_index]*g[pole_index][k];
521 value0 += (i[0][k]*A)/(ma[pole_index]*ma[pole_index]-s);
522 }
523
524 // Direct initial state to intermediate channels couplings
525 value1 += i[0][k]*fprod[k];
526
527 }
528
529 // Slowly varying polynomial term for the direct coupling
530 value1 *= (1.-s0_prod)/(s-s0_prod) ;
531
532 if(reject==true) return EvtComplex(9999., 9999.);
533 else return (value0+value1);
534
535}
536
537EvtComplex EvtD0ToKSpipi::amplitude_LASS(EvtVector4R p_k0l, EvtVector4R p_pip, EvtVector4R p_pim, string reso, double A_r, double Phi_r) {
538 double mR = 1.425 ;
539 double gammaR = 0.27 ;
540 double mab2 = 0.0;
541 if (reso == "k0spim") mab2 = pow((p_k0l + p_pim).mass(),2);
542 else if(reso == "k0spip") mab2 = pow((p_k0l + p_pip).mass(),2);
543 double s = mab2;
544
545 const double mD0 = 1.86483;
546 const double mKl = 0.49761;
547 const double mPi = 0.13957;
548
549 double _a = 0.113;
550 double _r = -33.8;
551 double _R = 1.0;
552 double _F = 0.96;
553 double _phiR = -1.9146;
554 double _phiF = 0.0017;
555 double fR = 1.0; // K*0(1430) has spin zero
556 int power = 1; // Power is 1 for spin zero
557
558 double mAB = sqrt(mab2);
559 double mA = mKl;
560 double mB = mPi;
561 double mC = mPi;
562 double mD = mD0;
563
564 double pAB=sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mAB*mAB));
565 double q=pAB;
566
567 double pR=sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) - mA*mA*mB*mB)/(mR*mR));
568 double q0=pR;
569
570 // Running width.
571 double g = gammaR*pow(q/q0,power)*(mR/mAB)*fR*fR;
572 EvtComplex propagator_relativistic_BreitWigner = 1./(mR*mR - mAB*mAB - EvtComplex(0.,mR*g));
573
574 // Non-resonant phase shift
575 double cot_deltaF = 1.0/(_a*q) + 0.5*_r*q;
576 double qcot_deltaF = 1.0/_a + 0.5*_r*q*q;
577
578 // Compute resonant part
579 EvtComplex expi2deltaF = EvtComplex(qcot_deltaF, q)/ EvtComplex(qcot_deltaF, -q);
580 EvtComplex resonant_term_T = _R * EvtComplex(cos(_phiR + 2 * _phiF), sin(_phiR + 2 * _phiF)) * propagator_relativistic_BreitWigner * mR * gammaR * mR / q0 * expi2deltaF;
581
582 // Compute non-resonant part
583 EvtComplex non_resonant_term_F = _F * EvtComplex(cos(_phiF), sin(_phiF)) * (cos(_phiF) + cot_deltaF * sin(_phiF)) * sqrt(s) / EvtComplex(qcot_deltaF, -q);
584
585 // Add non-resonant and resonant terms
586 EvtComplex LASS_contribution = non_resonant_term_F + resonant_term_T;
587
588 return EvtComplex(A_r*cos(Phi_r), A_r*sin(Phi_r)) * LASS_contribution;
589}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cc:175
imag
double imag(const EvtComplex &c)
Definition: EvtComplex.hh:246
abs2
double abs2(const EvtComplex &c)
Definition: EvtComplex.hh:221
meta
double meta
Definition: EvtPhokharaDef.hh:21
report
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
INFO
@ INFO
Definition: EvtReport.hh:52
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
EvtConst::radToDegrees
static const double radToDegrees
Definition: EvtConst.hh:30
EvtD0ToKSpipi::decay
void decay(EvtParticle *p)
Definition: EvtD0ToKSpipi.cc:110
EvtD0ToKSpipi::clone
EvtDecayBase * clone()
Definition: EvtD0ToKSpipi.cc:38
EvtD0ToKSpipi::getName
void getName(std::string &name)
Definition: EvtD0ToKSpipi.cc:34
EvtD0ToKSpipi::~EvtD0ToKSpipi
virtual ~EvtD0ToKSpipi()
Definition: EvtD0ToKSpipi.cc:32
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:520
EvtDecayBase::getArg
double getArg(int j)
Definition: EvtDecayBase.cc:564
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:34
EvtParticle::getP4Lab
EvtVector4R getP4Lab()
Definition: EvtParticle.cc:684
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39