BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDDalitz.cc
Go to the documentation of this file.
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: EvtDDalitz.cc
12//
13// Description: Routine to handle three-body decays of D0/D0_bar or D(s)+/D(s)-
14//
15// Modification history:
16//
17// NK September 3, 1997 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtGenBase/EvtPatches.hh"
22#include <stdlib.h>
23#include "EvtGenBase/EvtParticle.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtPDL.hh"
26#include "EvtGenBase/EvtReport.hh"
27#include "EvtGenBase/EvtResonance.hh"
28#include "EvtGenBase/EvtResonance2.hh"
29#include "EvtGenBase/EvtCubicSpline.hh"
30#include "EvtGenModels/EvtDDalitz.hh"
31#include <string>
32#include "EvtGenBase/EvtConst.hh"
33//#include "EvtGenBase/EvtFlatte.hh"
34#include "EvtGenBase/EvtDecayTable.hh"
35using std::endl;
36
37EvtDDalitz::~EvtDDalitz() {}
38
39void EvtDDalitz::getName(std::string& model_name){
40
41 model_name="D_DALITZ";
42
43}
44
45
46EvtDecayBase* EvtDDalitz::clone(){
47
48 return new EvtDDalitz;
49
50}
51
52void EvtDDalitz::init(){
53
54 // check that there are 0 arguments
55
56 static EvtId DM=EvtPDL::getId("D-");
57 static EvtId DP=EvtPDL::getId("D+");
58 static EvtId D0=EvtPDL::getId("D0");
59 static EvtId D0B=EvtPDL::getId("anti-D0");
60 static EvtId DSP=EvtPDL::getId("D_s+");
61 static EvtId DSM=EvtPDL::getId("D_s-");
62 static EvtId KM=EvtPDL::getId("K-");
63 static EvtId KP=EvtPDL::getId("K+");
64 static EvtId K0=EvtPDL::getId("K0");
65 static EvtId KB=EvtPDL::getId("anti-K0");
66 static EvtId KL=EvtPDL::getId("K_L0");
67 static EvtId KS=EvtPDL::getId("K_S0");
68 static EvtId PIM=EvtPDL::getId("pi-");
69 static EvtId PIP=EvtPDL::getId("pi+");
70 static EvtId PI0=EvtPDL::getId("pi0");
71
72 static double MPI = EvtPDL::getMeanMass(PI0);
73 static double MKP = EvtPDL::getMeanMass(KP);
74
75
76 checkNArg(0);
77 checkNDaug(3);
78
79 checkSpinParent(EvtSpinType::SCALAR);
80
81 checkSpinDaughter(0,EvtSpinType::SCALAR);
82 checkSpinDaughter(1,EvtSpinType::SCALAR);
83 checkSpinDaughter(2,EvtSpinType::SCALAR);
84
85
86 EvtId parnum=getParentId();
87 EvtId d1=getDaug(0);
88 EvtId d2=getDaug(1);
89 EvtId d3=getDaug(2);
90 _flag=0;
91 if ( parnum == D0 ) {
92 //look for either a K- pi+ pi0 or K0bar pi+ pi-
93 if ( d1==KM && d2==PIP && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;}
94 if ( d1==KM && d3==PIP && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
95 if ( d2==KM && d1==PIP && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
96 if ( d2==KM && d3==PIP && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
97 if ( d3==KM && d1==PIP && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
98 if ( d3==KM && d2==PIP && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
99
100 if ( d1==KB && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
101 if ( d1==KB && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
102 if ( d2==KB && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
103 if ( d2==KB && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
104 if ( d3==KB && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
105 if ( d3==KB && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
106
107 if ( d1==KL && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
108 if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
109 if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
110 if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
111 if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
112 if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
113
114 if ( d1==KS && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
115 if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
116 if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
117 if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
118 if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
119 if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
120
121 if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
122 if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
123 if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
124 if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
125 if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
126 if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
127
128 if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
129 if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
130 if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
131 if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
132 if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
133 if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
134
135 if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
136 if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
137 if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
138 if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
139 if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
140 if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
141
142 if ( d1==PIM && d2==PIP && d3==PI0 ) { _flag=11; _d1=0; _d2=1; _d3=2;} // DLY: D0 -> pi- pi+ pi0
143 if ( d1==PIM && d3==PIP && d2==PI0 ) { _flag=11; _d1=0; _d2=2; _d3=1;}
144 if ( d2==PIM && d1==PIP && d3==PI0 ) { _flag=11; _d1=1; _d2=0; _d3=2;}
145 if ( d2==PIM && d3==PIP && d1==PI0 ) { _flag=11; _d1=1; _d2=2; _d3=0;}
146 if ( d3==PIM && d1==PIP && d2==PI0 ) { _flag=11; _d1=2; _d2=0; _d3=1;}
147 if ( d3==PIM && d2==PIP && d1==PI0 ) { _flag=11; _d1=2; _d2=1; _d3=0;}
148
149 }
150 if ( parnum == D0B ) {
151 //look for either a K+ pi- pi0 or K0 pi+ pi-
152 if ( d1==KP && d2==PIM && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;}
153 if ( d1==KP && d3==PIM && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
154 if ( d2==KP && d1==PIM && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
155 if ( d2==KP && d3==PIM && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
156 if ( d3==KP && d1==PIM && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
157 if ( d3==KP && d2==PIM && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
158
159 if ( d1==K0 && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
160 if ( d1==K0 && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
161 if ( d2==K0 && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
162 if ( d2==K0 && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
163 if ( d3==K0 && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
164 if ( d3==K0 && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
165
166 if ( d1==KL && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
167 if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
168 if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
169 if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
170 if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
171 if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
172
173 if ( d1==KS && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
174 if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
175 if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
176 if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
177 if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
178 if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
179
180 if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
181 if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
182 if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
183 if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
184 if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
185 if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
186
187 if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
188 if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
189 if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
190 if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
191 if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
192 if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
193
194 if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
195 if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
196 if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
197 if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
198 if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
199 if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
200
201 if ( d1==PIP && d2==PIM && d3==PI0 ) { _flag=11; _d1=0; _d2=1; _d3=2;} // DLY: D0B -> pi+ pi- pi0
202 if ( d1==PIP && d3==PIM && d2==PI0 ) { _flag=11; _d1=0; _d2=2; _d3=1;}
203 if ( d2==PIP && d1==PIM && d3==PI0 ) { _flag=11; _d1=1; _d2=0; _d3=2;}
204 if ( d2==PIP && d3==PIM && d1==PI0 ) { _flag=11; _d1=1; _d2=2; _d3=0;}
205 if ( d3==PIP && d1==PIM && d2==PI0 ) { _flag=11; _d1=2; _d2=0; _d3=1;}
206 if ( d3==PIP && d2==PIM && d1==PI0 ) { _flag=11; _d1=2; _d2=1; _d3=0;}
207
208 }
209
210 if ( parnum == DP ) {
211 //look for K- pi+ pi+
212 if ( d1==KB && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
213 if ( d1==KB && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
214 if ( d2==KB && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
215 if ( d2==KB && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
216 if ( d3==KB && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
217 if ( d3==KB && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
218
219 if ( d1==KL && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
220 if ( d1==KL && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
221 if ( d2==KL && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
222 if ( d2==KL && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
223 if ( d3==KL && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
224 if ( d3==KL && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
225
226 if ( d1==KS && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
227 if ( d1==KS && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
228 if ( d2==KS && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
229 if ( d2==KS && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
230 if ( d3==KS && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
231 if ( d3==KS && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
232
233 if ( d1==KM && d2==PIP && d3==PIP ) { _flag=1; _d1=0; _d2=1; _d3=2;}
234 if ( d2==KM && d1==PIP && d3==PIP ) { _flag=1; _d1=1; _d2=0; _d3=2;}
235 if ( d3==KM && d1==PIP && d2==PIP ) { _flag=1; _d1=2; _d2=0; _d3=1;}
236
237 if ( d1==KM && d2==KP && d3==PIP ) { _flag=8; _d1=0; _d2=1; _d3=2;} // DLY: D+ -> K- K+ pi+
238 if ( d1==KM && d3==KP && d2==PIP ) { _flag=8; _d1=0; _d2=2; _d3=1;}
239 if ( d2==KM && d1==KP && d3==PIP ) { _flag=8; _d1=1; _d2=0; _d3=2;}
240 if ( d2==KM && d3==KP && d1==PIP ) { _flag=8; _d1=1; _d2=2; _d3=0;}
241 if ( d3==KM && d1==KP && d2==PIP ) { _flag=8; _d1=2; _d2=0; _d3=1;}
242 if ( d3==KM && d2==KP && d1==PIP ) { _flag=8; _d1=2; _d2=1; _d3=0;}
243
244 if ( d1==PIM && d2==PIP && d3==KP ) { _flag=9; _d1=0; _d2=1; _d3=2;} // DLY: D+ -> pi- pi+ K+
245 if ( d1==PIM && d3==PIP && d2==KP ) { _flag=9; _d1=0; _d2=2; _d3=1;}
246 if ( d2==PIM && d1==PIP && d3==KP ) { _flag=9; _d1=1; _d2=0; _d3=2;}
247 if ( d2==PIM && d3==PIP && d1==KP ) { _flag=9; _d1=1; _d2=2; _d3=0;}
248 if ( d3==PIM && d1==PIP && d2==KP ) { _flag=9; _d1=2; _d2=0; _d3=1;}
249 if ( d3==PIM && d2==PIP && d1==KP ) { _flag=9; _d1=2; _d2=1; _d3=0;}
250
251 if ( d1==PIM && d2==PIP && d3==PIP ) { _flag=10; _d1=0; _d2=1; _d3=2;} // DLY: D+ -> pi- pi+ pi+
252 if ( d2==PIM && d1==PIP && d3==PIP ) { _flag=10; _d1=1; _d2=0; _d3=2;}
253 if ( d3==PIM && d1==PIP && d2==PIP ) { _flag=10; _d1=2; _d2=0; _d3=1;}
254
255 }
256
257 if ( parnum == DM ) {
258 //look for K- pi+ pi+
259 if ( d1==K0 && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
260 if ( d1==K0 && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
261 if ( d2==K0 && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
262 if ( d2==K0 && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
263 if ( d3==K0 && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
264 if ( d3==K0 && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
265
266 if ( d1==KL && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
267 if ( d1==KL && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
268 if ( d2==KL && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
269 if ( d2==KL && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
270 if ( d3==KL && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
271 if ( d3==KL && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
272
273 if ( d1==KS && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;}
274 if ( d1==KS && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
275 if ( d2==KS && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
276 if ( d2==KS && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
277 if ( d3==KS && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
278 if ( d3==KS && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
279
280 if ( d1==KP && d2==PIM && d3==PIM ) { _flag=1; _d1=0; _d2=1; _d3=2;}
281 if ( d2==KP && d1==PIM && d3==PIM ) { _flag=1; _d1=1; _d2=0; _d3=2;}
282 if ( d3==KP && d1==PIM && d2==PIM ) { _flag=1; _d1=2; _d2=0; _d3=1;}
283
284 if ( d1==KP && d2==KM && d3==PIM ) { _flag=8; _d1=0; _d2=1; _d3=2;} // DLY: D- -> K+ K- pi-
285 if ( d1==KP && d3==KM && d2==PIM ) { _flag=8; _d1=0; _d2=2; _d3=1;}
286 if ( d2==KP && d1==KM && d3==PIM ) { _flag=8; _d1=1; _d2=0; _d3=2;}
287 if ( d2==KP && d3==KM && d1==PIM ) { _flag=8; _d1=1; _d2=2; _d3=0;}
288 if ( d3==KP && d1==KM && d2==PIM ) { _flag=8; _d1=2; _d2=0; _d3=1;}
289 if ( d3==KP && d2==KM && d1==PIM ) { _flag=8; _d1=2; _d2=1; _d3=0;}
290
291 if ( d1==PIP && d2==PIM && d3==KM ) { _flag=9; _d1=0; _d2=1; _d3=2;} // DLY: D- -> pi+ pi- K-
292 if ( d1==PIP && d3==PIM && d2==KM ) { _flag=9; _d1=0; _d2=2; _d3=1;}
293 if ( d2==PIP && d1==PIM && d3==KM ) { _flag=9; _d1=1; _d2=0; _d3=2;}
294 if ( d2==PIP && d3==PIM && d1==KM ) { _flag=9; _d1=1; _d2=2; _d3=0;}
295 if ( d3==PIP && d1==PIM && d2==KM ) { _flag=9; _d1=2; _d2=0; _d3=1;}
296 if ( d3==PIP && d2==PIM && d1==KM ) { _flag=9; _d1=2; _d2=1; _d3=0;}
297
298 if ( d1==PIP && d2==PIM && d3==PIM ) { _flag=10; _d1=0; _d2=1; _d3=2;} // DLY: D- -> pi+ pi- pi-
299 if ( d2==PIP && d1==PIM && d3==PIM ) { _flag=10; _d1=1; _d2=0; _d3=2;}
300 if ( d3==PIP && d1==PIM && d2==PIM ) { _flag=10; _d1=2; _d2=0; _d3=1;}
301
302 }
303
304 if ( parnum == DSP ) {
305 if ( d1==KM && d2==KP && d3==PIP ) { _flag=6; _d1=0; _d2=1; _d3=2; }
306 if ( d1==KM && d3==KP && d2==PIP ) { _flag=6; _d1=0; _d2=2; _d3=1; }
307 if ( d2==KM && d1==KP && d3==PIP ) { _flag=6; _d1=1; _d2=0; _d3=2; }
308 if ( d2==KM && d3==KP && d1==PIP ) { _flag=6; _d1=1; _d2=2; _d3=0; }
309 if ( d3==KM && d1==KP && d2==PIP ) { _flag=6; _d1=2; _d2=0; _d3=1; }
310 if ( d3==KM && d2==KP && d1==PIP ) { _flag=6; _d1=2; _d2=1; _d3=0; }
311
312 if ( d1==PIM && d2==PIP && d3==PIP ) { _flag=7; _d1=0; _d2=1; _d3=2;}
313 if ( d2==PIM && d1==PIP && d3==PIP ) { _flag=7; _d1=1; _d2=0; _d3=2;}
314 if ( d3==PIM && d1==PIP && d2==PIP ) { _flag=7; _d1=2; _d2=0; _d3=1;}
315 }
316
317 if ( parnum == DSM ) {
318 if ( d1==KP && d2==KM && d3==PIM ) { _flag=6; _d1=0; _d2=1; _d3=2; }
319 if ( d1==KP && d3==KM && d2==PIM ) { _flag=6; _d1=0; _d2=2; _d3=1; }
320 if ( d2==KP && d1==KM && d3==PIM ) { _flag=6; _d1=1; _d2=0; _d3=2; }
321 if ( d2==KP && d3==KM && d1==PIM ) { _flag=6; _d1=1; _d2=2; _d3=0; }
322 if ( d3==KP && d1==KM && d2==PIM ) { _flag=6; _d1=2; _d2=0; _d3=1; }
323 if ( d3==KP && d2==KM && d1==PIM ) { _flag=6; _d1=2; _d2=1; _d3=0; }
324
325 if ( d1==PIP && d2==PIM && d3==PIM ) { _flag=7; _d1=0; _d2=1; _d3=2;}
326 if ( d2==PIP && d1==PIM && d3==PIM ) { _flag=7; _d1=1; _d2=0; _d3=2;}
327 if ( d3==PIP && d1==PIM && d2==PIM ) { _flag=7; _d1=2; _d2=0; _d3=1;}
328 }
329
330/*
331 cout << "DDalitz : flag = " << _flag;
332 if (_flag==1) {
333 cout << " = D+ -> K- pi+ pi+"<<endl;
334 } else if (_flag==2) {
335 cout << " = D+ -> K0bar pi+ pi0"<<endl;
336 } else if (_flag==3) {
337 cout << " = D0 -> K0bar pi+ pi-"<<endl;
338 } else if (_flag==4) {
339 cout << " = D0 -> K- pi+ pi0"<<endl;
340 } else if (_flag==5) {
341 cout << " = D0 -> K0 K+ K-"<<endl;
342 } else if (_flag==6) {
343 cout << " = Ds+ -> K- K+ pi+"<<endl;
344 } else if (_flag==7) {
345 cout << " = Ds+ -> pi+ pi- pi+"<<endl;
346 } else if (_flag==8) {
347 cout << " = D+ -> K- K+ pi+"<<endl;
348 } else if (_flag==9) {
349 cout << " = D+ -> pi- pi+ K+ WS (DCS)"<<endl;
350 } else if (_flag==10) {
351 cout << " = D+ -> pi- pi+ pi+"<<endl;
352 } else if (_flag==11) {
353 cout << " = D0 -> pi- pi+ pi0"<<endl;
354 } else if (_flag==0) {
355 cout << " = Invalid mode."<<endl;
356 }
357*/
358
359 if ( _flag==6) {
360 _kkpi_params.push_back(EvtFlatteParam(MPI, MPI, 0.406));
361 _kkpi_params.push_back(EvtFlatteParam(MKP, MKP, 0.800));
362 }
363
364 if ( _flag==0) {
365 report(ERROR,"EvtGen") << "EvtDDaltiz: Invalid mode."<<endl;
366 assert(0);
367 }
368}
369
370void EvtDDalitz::initProbMax() {
371 if ( _flag==1 ) {setProbMax(2500.0);}
372 if ( _flag==2 ) {setProbMax(10000);} //old if ( _flag==2 ) {setProbMax(147.9);}
373 if ( _flag==3 ) {setProbMax(5000.0);}
374 if ( _flag==4 ) {setProbMax(3000.0);}
375 if ( _flag==5 ) {setProbMax(10000000.0);}
376 if ( _flag==6 ) {setProbMax(50000.0);}
377 if ( _flag==7 ) {setProbMax(50000.0);}
378 if ( _flag==8 ) {setProbMax(35000.0);}
379 if ( _flag==9 ) {setProbMax(2500.0);}
380 if ( _flag==10 ){setProbMax(1300.0);}
381 if ( _flag==11 ){setProbMax(1000.0);}
382}
383
384void EvtDDalitz::decay( EvtParticle *p){
385
386 //added by Lange Jan4,2000
387 static EvtId BP = EvtPDL::getId("B+");
388 static EvtId BM = EvtPDL::getId("B-");
389 static EvtId B0 = EvtPDL::getId("B0");
390 static EvtId B0B = EvtPDL::getId("anti-B0");
391 static EvtId DM=EvtPDL::getId("D-");
392 static EvtId DP=EvtPDL::getId("D+");
393 static EvtId D0=EvtPDL::getId("D0");
394 static EvtId D0B=EvtPDL::getId("anti-D0");
395 static EvtId KM=EvtPDL::getId("K-");
396 static EvtId KP=EvtPDL::getId("K+");
397 static EvtId K0=EvtPDL::getId("K0");
398 static EvtId KB=EvtPDL::getId("anti-K0");
399 static EvtId PIM=EvtPDL::getId("pi-");
400 static EvtId PIP=EvtPDL::getId("pi+");
401 static EvtId PI0=EvtPDL::getId("pi0");
402
403 double oneby2 = 0.707106782;
404
405 bool isBToDK=false;
406 if ( p -> getParent () ) {
407 std::string name1;
408 EvtDecayTable::getDecayFunc(p->getParent())->getName(name1);
409
410 EvtId parId = p -> getParent()->getId ();
411 if ( ( BP == parId ) || ( BM == parId ) || ( B0 == parId ) ||
412 ( B0B == parId ) )
413 if (name1 == "BTODDALITZCPK") isBToDK=true;
414 }
415
416
417//same structure for all of these decays
418
419 p->initializePhaseSpace(getNDaug(),getDaugs());
420 EvtVector4R moms1 = p->getDaug(_d1)->getP4();
421 EvtVector4R moms2 = p->getDaug(_d2)->getP4();
422 EvtVector4R moms3 = p->getDaug(_d3)->getP4();
423
424 EvtVector4R p4_p;
425 p4_p.set(p->mass(),0.0,0.0,0.0);
426
427 EvtComplex amp(1.0,0.0);
428
429//now determine which D and which decay
430
431//data from Anjos et al, Phys.Rev.D 1993, v.48,num.1,p.56 (E691 resuls)
432//for D+ -> K- pi+ pi+, and from Adler et al, Phys.Lett. B196 (1987), 107
433//(Mark III results) for D+ -> K0bar pi+ pi0.
434 //CLEO results for D0->k-pi+pi0
435
436 if ( _flag==1) {
437
438//have a D+ -> K- pi+ pi+ decay, or charge conjugate
439//Anjos etal e691 - Phys Rev D48, 56 (1993)
440// EvtResonance DplusRes11(p4_p,moms1,moms2,0.78,-60.0,0.0498,0.89610,1);
441// EvtResonance DplusRes12(p4_p,moms3,moms1,0.78,-60.0,0.0498,0.89610,1);//K*(892)
442//
443// EvtResonance DplusRes21(p4_p,moms1,moms2,0.53,132.0,0.287,1.429,0);
444// EvtResonance DplusRes22(p4_p,moms3,moms1,0.53,132.0,0.287,1.429,0);//K*(1430)
445//
446// EvtResonance DplusRes31(p4_p,moms1,moms2,0.47,-51.0,0.323,1.714,1);
447// EvtResonance DplusRes32(p4_p,moms3,moms1,0.47,-51.0,0.323,1.714,1);//K*(1680)
448//
449// amp = amp + oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl()) + oneby2*(DplusRes21.resAmpl() + DplusRes22.resAmpl()) + oneby2*(-DplusRes31.resAmpl()+ DplusRes32.resAmpl());
450 //CLEO-c p15,arxiv:0802.4214v2
451 EvtResonance2 DplusRes11( p4_p, moms1, moms2, 1.0, 0.0, 0.0503, 0.896,
452 1, true );
453 EvtResonance2 DplusRes12( p4_p, moms3, moms1, 1.0, 0.0, 0.0503, 0.896,
454 1, true ); //K*(892)
455 EvtResonance2 DplusRes21( p4_p, moms1, moms2, 3.0, 49.7 - 180.0, 0.164,
456 1.463, 0 );
457 EvtResonance2 DplusRes22( p4_p, moms3, moms1, 3.0, 49.7 - 180.0, 0.164,
458 1.463, 0 ); //K*(1430)
459 EvtResonance2 DplusRes31( p4_p, moms1, moms2, 0.96, -29.9 + 180.0,
460 0.109, 1.4324, 2, true );
461 EvtResonance2 DplusRes32( p4_p, moms3, moms1, 0.96, -29.9 + 180.0,
462 0.109, 1.4324, 2, true ); // K*_2(1430)
463 EvtResonance2 DplusRes41( p4_p, moms1, moms2, 6.5, 29.0, 0.323, 1.717,
464 1, true );
465 EvtResonance2 DplusRes42( p4_p, moms3, moms1, 6.5, 29.0, 0.323, 1.717,
466 1, true ); //K*(1680)
467 EvtResonance2 DplusRes51( p4_p, moms1, moms2, 5.01, -163.7 + 180.0,
468 0.470, 0.809, 0 );
469 EvtResonance2 DplusRes52( p4_p, moms3, moms1, 5.01, -163.7 + 180.0,
470 0.470, 0.809, 0 ); //kappa(800)
471 double pi180inv = 1.0 / EvtConst::radToDegrees;
472 amp = EvtComplex( 7.4 * cos( ( -18.4 + 180.0 ) * pi180inv ),
473 7.4 * sin( ( -18.4 + 180.0 ) * pi180inv ) ) +
474 oneby2 * ( -DplusRes11.resAmpl() + DplusRes12.resAmpl() ) +
475 oneby2 * ( DplusRes21.resAmpl() + DplusRes22.resAmpl() ) +
476 oneby2 * ( DplusRes31.resAmpl() + DplusRes32.resAmpl() ) +
477 oneby2 * ( -DplusRes41.resAmpl() + DplusRes42.resAmpl() ) +
478 oneby2 * ( DplusRes51.resAmpl() + DplusRes52.resAmpl() );
479
480 }
481
482 if ( _flag==2) {
483
484//have a D+ -> K0bar pi+ pi0 decay (daughters must be listed in this order!)
485//adler etal MarkIII - Phys Lett B196, 107 (1987)
486// Results in this paper:
487// Kbar rho+ FitFraction = 68+/-8+/-12 Phase 0
488// Kbar* pi+ 19+/-6+/-6 43+/-23
489// nonres 13+/-7+/-8 250+/-19
490// These numbers below seem not to be exactly the same
491// the phases are equiv to -106=254 and 41
492//
493// old EvtResonance DplusKpipi0Res1(p4_p,moms2,moms3,1.00,0.00,0.1512,0.7699,1); //rho+
494// old EvtResonance DplusKpipi0Res2(p4_p,moms3,moms1,0.8695,0.7191,0.0498,0.89159,1); //K*0
495// old amp = 0.9522*EvtComplex(cos(-1.8565),sin(-1.8565)) + 1.00*DplusKpipi0Res1.relBrWig(0) + 0.8695*EvtComplex(cos(0.7191),sin(0.7191))*DplusKpipi0Res2.relBrWig(1);
496 double pi180inv = 1.0/EvtConst::radToDegrees;
497 // NR rho K*(892) K0(1430) K2(1430) K(1680) rho(1450) kappa
498 //double amp1[8]={0.920322, 1, 0.198943, 0.604074, 0.225741, 1.257776, 1.989557, 1.332243};
499 //double phi[8]={278.810706, 0, 288.035519, 324.018021, 263.172850, 276.218704, 169.42471, 88.091929};
500 //double amp1[8]={0.797929, 1, 0.194377, 0.332583, 0.286707, 1.440519, 2.278415, 0.611784};
501 //double phi[8]={272.941122, 0, 288.439024, 329.600012, 258.197983,273.091907, 172.456067, 85.072272};
502 //double amp1[8]={0.8, 1, 0.2, 0.3, 0, 0, 2.5, 0};
503 //double phi[8] ={273, 0, 288, 330, 258, 0, 172, 85};
504 // Final result
505 double amp1[8]={0.797929, 1, 0.194377, 0.332583, 0, 1.440519, 2.278415, 0.611784};
506 double phi[8]={279, 0, 293, 334, 0, 252, 187, 93};
507 EvtResonance2 DplusKpipi0Res1(p4_p,moms2,moms3,amp1[1],phi[1],0.1494,0.77549,1);
508 EvtResonance2 DplusKpipi0Res2(p4_p,moms1,moms3,amp1[2],phi[2],0.0503,0.896 ,1);
509 EvtResonance2 DplusKpipi0Res3(p4_p,moms1,moms3,amp1[3],phi[3],0.175 ,1.465 ,0);
510 EvtResonance2 DplusKpipi0Res4(p4_p,moms1,moms3,amp1[4],phi[4],0.109 ,1.4324 ,2);
511 EvtResonance2 DplusKpipi0Res5(p4_p,moms1,moms3,amp1[5],phi[5],0.322 ,1.717 ,1);
512 // Fitted
513 EvtResonance2 DplusKpipi0Res6(p4_p,moms2,moms3,amp1[6],phi[6],0.19 ,1.464 ,1);
514 EvtResonance2 DplusKpipi0Res7(p4_p,moms1,moms3,1,0, 0.175, 1.465, 0);
515 amp = amp1[0]*EvtComplex(cos(phi[0]*pi180inv),sin(phi[0]*pi180inv)) + DplusKpipi0Res1.resAmpl()+DplusKpipi0Res2.resAmpl()
516 + DplusKpipi0Res3.resAmpl()+DplusKpipi0Res4.resAmpl()+ DplusKpipi0Res5.resAmpl() +DplusKpipi0Res6.resAmpl();
517 //EvtComplex kappa=1./(EvtComplex(0.731086,-0.244799)*EvtComplex(0.731086,-0.244799)-(moms1+moms3).mass2());
518 //EvtComplex kappa=1./(EvtComplex(0.710843,-0.251841)*EvtComplex(0.710843,-0.251841)-(moms1+moms3).mass2());
519 // Final result
520 EvtComplex kappa=1./(EvtComplex(0.752,-0.229)*EvtComplex(0.752,-0.229)-(moms1+moms3).mass2());
521 amp += amp1[7]*EvtComplex(cos(phi[7]*pi180inv),sin(phi[7]*pi180inv))*kappa;
522 }
523
524 if(_flag==3) {
525
526//know it's a D0 -> K0_bar pi+ pi- or charge conjugate
527
528// EvtResonance DK2piRes1(p4_p,moms1,moms2,2.31,109.0,0.0498,0.89610,1);//K*(892)
529// EvtResonance DK2piRes2(p4_p,moms3,moms2,1.59,-123.0,0.1491,0.7683,1);//RHO(770)
530
531// amp = amp + DK2piRes1.resAmpl() + DK2piRes2.resAmpl();
532
533
534 // D0 -> K0 pi+ pi- + CC
535 // If it does not come from a B->DK, decay it as D0 or D0bar separatly
536 // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
537 // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
538
539 if ( isBToDK ) {
540 // Gamma angle in rad.
541 double gamma = EvtDecayTable::getDecayFunc( p->getParent() )
542 -> getArg( 0 ) ;
543 // Strong phase in rad.
544 double delta = EvtDecayTable::getDecayFunc( p->getParent() )
545 -> getArg( 1 ) ;
546 // Ratio between B->D0K and B->D0barK
547 double A = EvtDecayTable::getDecayFunc( p->getParent() )
548 -> getArg( 2 ) ;
549
550 EvtComplex Factor( fabs( A ) * cos ( delta ) ,
551 fabs( A ) * sin ( delta ) ) ;
552
553 if ( ( p->getParent()->getId() == BP ) ||
554 ( p->getParent()->getId() == B0 ) ) {
555 // the ratio D/Dbar
556 Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
557 if ( p->getId() == D0 ) {
558 // the flavor of the particle has no meaning. But we need
559 // it to know which daughter is pi+ or pi-
560 // M( B+ or B0 ) = f(Dbar) + factor * f(D)
561 // f(Dbar) = amplDtoK0PiPi(pD, K0, pi+, pi-)
562 // f(D) = amplDtoK0PiPi(pD, K0, pi-, pi+)
563 // Then ...
564 amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
565 Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
566 }
567 else {
568 amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
569 Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
570 }
571 }
572 else if ( ( p->getParent() -> getId() == BM ) ||
573 ( p->getParent() -> getId() == B0B ) ) {
574 Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
575 // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
576 if ( p->getId() == D0 ) {
577 amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
578 Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
579 }
580 else {
581 amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
582 Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
583 }
584 }
585 }
586 else {
587 amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
588 }
589 }
590
591
592 if(_flag==4) {
593
594 EvtResonance2 DKpipi0Res1(p4_p,moms2,moms3,1.0 ,0.0 ,0.1507,0.770 ,1); //rho
595 EvtResonance2 DKpipi0Res2(p4_p,moms1,moms2,0.39, -0.2 ,0.0505,0.8961,1); //k*0
596 EvtResonance2 DKpipi0Res3(p4_p,moms1,moms3,0.44, 163.0 ,0.050 ,0.8915,1); //k*-
597
598 EvtResonance2 DKpipi0Res4(p4_p,moms1,moms3,0.77 ,55.5 ,0.294 ,1.412 ,0); //k01430-
599 EvtResonance2 DKpipi0Res5(p4_p,moms1,moms2,0.85 ,166.0 ,0.294 ,1.412 ,0); //k01430bar
600 EvtResonance2 DKpipi0Res6(p4_p,moms2,moms3,2.5 ,171.0 ,0.240 ,1.700 ,1); //rho1700
601 EvtResonance2 DKpipi0Res7(p4_p,moms1,moms3,2.5 ,103.0 ,0.322 ,1.717 ,1); //K*1680-
602
603
604
605 double pi180inv = 1.0/EvtConst::radToDegrees;
606
607 amp = EvtComplex(1.75*cos(31.2*pi180inv),1.75*sin(31.2*pi180inv))
608 + DKpipi0Res1.resAmpl() + DKpipi0Res2.resAmpl() + DKpipi0Res3.resAmpl()
609 + DKpipi0Res4.resAmpl() + DKpipi0Res5.resAmpl()
610 + DKpipi0Res6.resAmpl()
611 + DKpipi0Res7.resAmpl();
612
613 }
614
615
616 if(_flag==5) {
617
618 // D0 -> K0 K+ K- + CC
619 // If it does not come from a B->DK, decay it as D0 or D0bar separatly
620 // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
621 // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
622
623 if ( isBToDK ){
624 // Gamma angle in rad.
625 double gamma = EvtDecayTable::getDecayFunc( p->getParent() )
626 -> getArg( 0 ) ;
627 // Strong phase in rad.
628 double delta = EvtDecayTable::getDecayFunc( p->getParent() )
629 -> getArg( 1 ) ;
630 // Ratio between B->D0K and B->D0barK
631 double A = EvtDecayTable::getDecayFunc( p->getParent() )
632 -> getArg( 2 ) ;
633
634 EvtComplex Factor( fabs( A ) * cos ( delta ) ,
635 fabs( A ) * sin ( delta ) ) ;
636
637 if ( ( p->getParent()->getId() == BP ) ||
638 ( p->getParent()->getId() == B0 ) ) {
639 // the ratio D/Dbar
640 Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
641 if ( p->getId() == D0 ) {
642 // the flavor of the particle has no meaning. But we need
643 // it to know which daughter is pi+ or pi-
644 // M( B+ or B0 ) = f(Dbar) + factor * f(D)
645 // f(Dbar) = amplDtoK0PiPi(pD, K0, K+, K-)
646 // f(D) = amplDtoK0PiPi(pD, K0, K-, K+)
647 // Then ...
648 amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
649 Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
650 }
651 else {
652 amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
653 Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
654 }
655 }
656 else if ( ( p->getParent() -> getId() == BM ) ||
657 ( p->getParent() -> getId() == B0B ) ) {
658 Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
659 // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
660 if ( p->getId() == D0 ) {
661 amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
662 Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
663 }
664 else {
665 amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
666 Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
667 }
668 }
669 }
670 else {
671 amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
672 }
673 }
674
675 // Ds -> K K pi
676 if(_flag==6) {
677 EvtResonance2 DsKKpiRes1(p4_p, moms3, moms1, 1.0, 0.0, 0.0455, 0.8944, 1, true); // K*(892)
678 EvtResonance2 DsKKpiRes2(p4_p, moms3, moms1, 1.48, 138., 0.290, 1.414, 0); // K*_0(1430)
679 EvtFlatte DsKKpiRes3(p4_p, moms1, moms2, 5.07, 156., 0.965, _kkpi_params); // f_0(980)
680 EvtResonance2 DsKKpiRes4(p4_p, moms1, moms2, 1.15, -10., 0.00426, 1.019455, 1, true); // phi(1020)
681 EvtResonance2 DsKKpiRes5(p4_p, moms1, moms2, 1.28, 53., 0.265, 1.350, 0); // f_0(1370)
682 EvtResonance2 DsKKpiRes6(p4_p, moms1, moms2, 1.19, 87., 0.137, 1.724, 0); // f_0(1710)
683 amp = DsKKpiRes1.resAmpl() + DsKKpiRes2.resAmpl() + DsKKpiRes3.resAmpl()
684 + DsKKpiRes4.resAmpl() + DsKKpiRes5.resAmpl() + DsKKpiRes6.resAmpl();
685 }
686
687 // D -> pi pi pi
688 if(_flag==7) {
689// BaBar decay amplitudes for Ds+ -> pi- pi+ pi+, or charge conjugate
690// p4_p is Ds+
691// moms1 is pi-
692// moms2 is pi+
693// moms3 is pi+
694// Amplitudes and phases are taken from BaBar hep-ex/0808.0971 with MIPWA s-wave component
695
696 EvtResonance2 DspipipiRes1(p4_p, moms1, moms2, 1.0, 0.0, 0.1867, 1.2755, 2, true,false); // f2(1270), 12
697 EvtResonance2 DspipipiRes2(p4_p, moms1, moms3, 1.0, 0.0, 0.1867, 1.2755, 2, true,false); // f2(1270), 13
698
699 EvtResonance2 DspipipiRes3(p4_p,moms1,moms2, 1.25706e-01,5.44277*EvtConst::radToDegrees ,0.1491,0.77526 ,1, true,false); //rho, 12
700 EvtResonance2 DspipipiRes4(p4_p,moms1,moms3, 1.25706e-01,5.44277*EvtConst::radToDegrees ,0.1491,0.77526 ,1, true,false); //rho, 13
701 EvtResonance2 DspipipiRes5(p4_p,moms1,moms2, 9.06766e-01,1.03370*EvtConst::radToDegrees ,0.400,1.465,1, true,false); //rho(1450), 12
702 EvtResonance2 DspipipiRes6(p4_p,moms1,moms3, 9.06766e-01,1.03370*EvtConst::radToDegrees ,0.400,1.465,1, true,false); //rho(1450), 13
703
704 const double xl[] = {0.28, 0.448, 0.55, 0.647, 0.736, 0.803, 0.873, 0.921, 0.951, 0.968, 0.981, 0.993, 1.024, 1.078, 1.135, 1.193, 1.235, 1.267, 1.297, 1.323, 1.35, 1.376, 1.402, 1.427, 1.455, 1.492, 1.557, 1.64, 1.735};
705 const double ym[] = {1.23487, 2.80323, 3.42035, 3.31677, 5.45084, 6.21976, 7.88104, 11.8548, 16.8433, 21.7424, 26.455, 18.6398, 11.1707, 8.00181, 6.74089, 6.09654, 6.62815, 6.26746, 6.50165, 7.49954, 7.2659, 7.53306, 8.49471, 8.07989, 8.28178, 5.8168, 1.63515, 1.37752, 2.08794};
706 const double yp[] = {-3.59059, -3.82322, -3.87059, -3.74264, -3.37687, -3.09568, -2.59601, -2.15795, -1.76578, -1.21117, -0.577655, -0.246062, 0.17292, 0.550796, 0.981734, 1.27572, 1.32191, 1.56281, 1.46746, 1.60172, 1.75168, 1.79924, 1.94462, 2.09212, 2.53836, 3.06809, 3.05283, 7.0558, 7.31884};
707 EvtCubicSpline::setParams(29, xl, ym,yp);
708 EvtCubicSpline Dspipipiswave1(p4_p,moms1,moms2); //S-wave, 12
709 EvtCubicSpline Dspipipiswave2(p4_p,moms1,moms3); //S-wave, 13
710
711 amp = DspipipiRes1.resAmpl() + DspipipiRes2.resAmpl() + DspipipiRes3.resAmpl()
712 + DspipipiRes4.resAmpl() + DspipipiRes5.resAmpl() + DspipipiRes6.resAmpl()
713 + Dspipipiswave1.resAmpl() + Dspipipiswave2.resAmpl();
714
715 }
716
717 //D+ -> K- K+ pi+
718 //CLEO PRD 78, 072003 (2008) Fit A
719 if(_flag==8) {
720
721 EvtResonance2 DpKKpiRes1( p4_p, moms3, moms1, 1.0, 0.0, 0.0503, 0.8960, 1, true ); // K*(892)
722 EvtResonance2 DpKKpiRes2( p4_p, moms3, moms1, 3.7, 73.0, 0.290, 1.414, 0 ); // K*_0(1430)
723 EvtResonance2 DpKKpiRes3( p4_p, moms1, moms2, 1.189, -179.0 + 180.0, 0.00426, 1.019455, 1, true ); // phi(1020)
724 EvtResonance2 DpKKpiRes4( p4_p, moms1, moms2, 1.72, 123., 0.265, 1.474, 0 ); // a_0(1450)
725 EvtResonance2 DpKKpiRes5( p4_p, moms1, moms2, 1.9, -52.0 + 180.0, 0.15, 1.68, 1, true ); // phi(1680)
726 EvtResonance2 DpKKpiRes6( p4_p, moms3, moms1, 6.4, 150., 0.109, 1.4324, 2, true ); // K*_2(1430)
727 double pi180inv = 1.0 / EvtConst::radToDegrees;
728 amp = EvtComplex( 5.1 * cos( ( 53.0 ) * pi180inv ), 5.1 * sin( ( 53.0 ) * pi180inv ) ) +
729 DpKKpiRes1.resAmpl() + DpKKpiRes2.resAmpl() + DpKKpiRes3.resAmpl() +
730 DpKKpiRes4.resAmpl() + DpKKpiRes5.resAmpl() + DpKKpiRes6.resAmpl();
731
732 }
733
734 //D+ -> pi- pi+ K+ WS (DCS)
735 //FOCUS PLB 601 10 (2004) ; amplitudes there are individually normalized (although not explicit in the paper)
736 // thus the magnitudes appearing below come from dividing the ones appearing in the paper by the sqrt of the
737 // integral over the DP of the corresponding squared amplitude. Writing as pi- pi+ K+ so pipi resonances are (12)
738 // and Kpi resonances are (31); masses and widths corresponds to PDG 2010
739 if(_flag==9) {
740 EvtResonance2 DpKpipiDCSRes1( p4_p, moms1, moms2, 1.0, 0.0, 0.149, 0.775, 1, true ); // rho(770)
741 EvtResonance2 DpKpipiDCSRes2( p4_p, moms3, moms1, 1.0971, -167.1, 0.0487, 0.896, 1, true ); // K*(890)
742 EvtResonance2 DpKpipiDCSRes3( p4_p, moms1, moms2, 0.4738, -134.5, 0.059, 0.972, 0 ); // f0(980) as simple BW
743 EvtResonance2 DpKpipiDCSRes4( p4_p, moms3, moms1, 2.2688, 54.4, 0.109, 1.432, 2, true ); // K*2(1430)
744 amp = DpKpipiDCSRes1.resAmpl() + DpKpipiDCSRes2.resAmpl() +
745 DpKpipiDCSRes3.resAmpl() + DpKpipiDCSRes4.resAmpl();
746 }
747
748 // D+ -> pi- pi+ pi+ from E791 [PRL 86 770 (2001)]
749 // masses and widths below correspond to what they used; there, the amplitudes were individually normalized
750 // (although not explicit) so magnitudes here are obtained after correcting for that
751 // Breit-Wigner has a factor of (-1) there which changes the relative phase of the NR wrt to the resonances
752 // thus the NR magnitude is set as negative
753 if(_flag==10) {
754 EvtResonance2 DppipipiRes11( p4_p, moms1, moms2, 1.0, 0.0, 0.150, 0.769, 1, true ); // rho(770)
755 EvtResonance2 DppipipiRes12( p4_p, moms3, moms1, 1.0, 0.0, 0.150, 0.769, 1, true ); // rho(770)
756 EvtResonance2 DppipipiRes21( p4_p, moms1, moms2, 2.2811, 205.7, 0.324, 0.478, 0 ); // sigma(500)
757 EvtResonance2 DppipipiRes22( p4_p, moms3, moms1, 2.2811, 205.7, 0.324, 0.478, 0 ); // sigma(500)
758 EvtResonance2 DppipipiRes31( p4_p, moms1, moms2, 0.4265, 165.0, 0.044, 0.977, 0 ); // f0(980) simple BW
759 EvtResonance2 DppipipiRes32( p4_p, moms3, moms1, 0.4265, 165.0, 0.044, 0.977, 0 ); // f0(980) simple BW
760 EvtResonance2 DppipipiRes41( p4_p, moms1, moms2, 2.0321, 57.3, 0.185, 1.275, 2, true ); // f2(1270)
761 EvtResonance2 DppipipiRes42( p4_p, moms3, moms1, 2.0321, 57.3, 0.185, 1.275, 2, true ); // f2(1270)
762 EvtResonance2 DppipipiRes51( p4_p, moms1, moms2, 0.7888, 105.4, 0.173, 1.434, 0 ); // f0(1370)
763 EvtResonance2 DppipipiRes52( p4_p, moms3, moms1, 0.7888, 105.4, 0.173, 1.434, 0 ); // f0(1370)
764 EvtResonance2 DppipipiRes61( p4_p, moms1, moms2, 0.7363, 319.1, 0.310, 1.465, 1, true ); // rho(1450)
765 EvtResonance2 DppipipiRes62( p4_p, moms3, moms1, 0.7363, 319.1, 0.310, 1.465, 1, true ); // rho(1450)
766 double pi180inv = 1.0 / EvtConst::radToDegrees;
767 amp = EvtComplex( -3.98 * cos( 57.3 * pi180inv ), -3.98 * sin( 57.3 * pi180inv ) ) +
768 ( DppipipiRes11.resAmpl() - DppipipiRes12.resAmpl() ) //spin1
769 + ( DppipipiRes21.resAmpl() + DppipipiRes22.resAmpl() ) +
770 ( DppipipiRes31.resAmpl() + DppipipiRes32.resAmpl() ) +
771 ( DppipipiRes41.resAmpl() + DppipipiRes42.resAmpl() ) +
772 ( DppipipiRes51.resAmpl() + DppipipiRes52.resAmpl() ) +
773 ( DppipipiRes61.resAmpl() - DppipipiRes62.resAmpl() ); //spin1
774 }
775
776 //D0 -> pi- pi+ pi0
777 //PRL 99, 251801 (2007) arXiv:hep-ex/0703037
778 if(_flag==11) {
779 EvtResonance2 DpipipiRes1p( p4_p, moms2, moms3, 1.0, 0.0, 0.149, 0.775, 1, true ); //rho+(770)
780 EvtResonance2 DpipipiRes1( p4_p, moms1, moms2, 0.588, 16.2, 0.149, 0.775, 1, true ); //rho0(770)
781 EvtResonance2 DpipipiRes1m( p4_p, moms3, moms1, 0.714, -2.0, 0.149, 0.775, 1, true ); //rho-(770)
782 EvtResonance2 DpipipiRes2p( p4_p, moms2, moms3, 0.21, -146.0, 0.400, 1.465, 1, true ); //rho+(1450)
783 EvtResonance2 DpipipiRes2( p4_p, moms1, moms2, 0.33, 10.0, 0.400, 1.465, 1, true ); //rho0(1450)
784 EvtResonance2 DpipipiRes2m( p4_p, moms3, moms1, 0.82, 16.0, 0.400, 1.465, 1, true ); //rho-(1450)
785 EvtResonance2 DpipipiRes3p( p4_p, moms2, moms3, 2.25, -17.0, 0.250, 1.720, 1, true ); //rho+(1700)
786 EvtResonance2 DpipipiRes3( p4_p, moms1, moms2, 2.51, -17.0, 0.250, 1.720, 1, true ); //rho0(1700)
787 EvtResonance2 DpipipiRes3m( p4_p, moms3, moms1, 2.00, -50.0, 0.250, 1.720, 1, true ); //rho-(1700)
788 EvtResonance2 DpipipiRes4( p4_p, moms1, moms2, 0.015, -59.0, 0.07, 0.980, 0 ); //f0(980)
789 EvtResonance2 DpipipiRes5( p4_p, moms1, moms2, 0.063, 156.0, 0.350, 1.370, 0 ); //f0(1370)
790 EvtResonance2 DpipipiRes6( p4_p, moms1, moms2, 0.058, 12.0, 0.109, 1.505, 0 ); //f0(1500)
791 EvtResonance2 DpipipiRes7( p4_p, moms1, moms2, 0.112, 51.0, 0.135, 1.720, 0 ); //f0(1720)
792 EvtResonance2 DpipipiRes8( p4_p, moms1, moms2, 1.04, -171.0, 0.185, 1.275, 2, true ); //f2(1270)
793 EvtResonance2 DpipipiRes9( p4_p, moms1, moms2, 0.069, 8.0, 0.600, 0.400, 0 ); //sigma(400)
794
795 double pi180inv = 1.0 / EvtConst::radToDegrees;
796 amp = EvtComplex( 0.57 * cos( -11.0 * pi180inv ), 0.57 * sin( -11.0 * pi180inv ) ) +
797 DpipipiRes1p.resAmpl() + DpipipiRes1.resAmpl() +
798 DpipipiRes1m.resAmpl() + DpipipiRes2p.resAmpl() +
799 DpipipiRes2.resAmpl() + DpipipiRes2m.resAmpl() +
800 DpipipiRes3p.resAmpl() + DpipipiRes3.resAmpl() +
801 DpipipiRes3m.resAmpl() + DpipipiRes4.resAmpl() +
802 DpipipiRes5.resAmpl() + DpipipiRes6.resAmpl() +
803 DpipipiRes7.resAmpl() + DpipipiRes8.resAmpl() +
804 DpipipiRes9.resAmpl();
805 }
806
807 vertex(amp);
808
809 return ;
810}
811
812EvtComplex EvtDDalitz::amplDtoK0PiPi(EvtVector4R p4_p, EvtVector4R moms1,
813 EvtVector4R moms2, EvtVector4R moms3) {
814
815 //K*(892)-
816 EvtResonance2 DK2piRes1(p4_p,moms1,moms2,1.418,-190.0,0.0508,0.89166,1);
817 //K0*(1430)
818 EvtResonance2 DK2piRes2(p4_p,moms1,moms2,1.818,-337.0,0.294 ,1.412 ,0);
819 //K2*(1430)
820 EvtResonance2 DK2piRes3(p4_p,moms1,moms2,0.909, -5.0,0.0985,1.4256 ,2);
821 //K*(1680)
822 EvtResonance2 DK2piRes4(p4_p,moms1,moms2,5.091,-166.0,0.322 ,1.717 ,1);
823 //DCS K*(892)
824 EvtResonance2 DK2piRes5(p4_p,moms1,moms3,0.100, -19.0,0.0508,0.89166,1);
825
826 //Rho
827 EvtResonance2 DK2piRes6(p4_p,moms3,moms2,0.909,-340.0,0.1502,0.7693,1);
828 //Omega
829 EvtResonance2 DK2piRes7(p4_p,moms3,moms2,.0336,-226.0,0.00844,0.78257,1);
830 //f0(980)
831 EvtResonance2 DK2piRes8(p4_p,moms3,moms2,0.309,-152.0,0.05,0.977,0);
832 //f0(1370)
833 EvtResonance2 DK2piRes9(p4_p,moms3,moms2,1.636,-255.0,0.272,1.31,0);
834 //f2(1270)
835 EvtResonance2 DK2piRes10(p4_p,moms3,moms2,0.636,-32.0,0.1851,1.2754,2);
836
837 return EvtComplex(1.0,0.0) +
838 DK2piRes1.resAmpl() + DK2piRes2.resAmpl() +
839 DK2piRes3.resAmpl() + DK2piRes4.resAmpl() +
840 DK2piRes5.resAmpl() + DK2piRes6.resAmpl() +
841 DK2piRes7.resAmpl() + DK2piRes8.resAmpl() +
842 DK2piRes9.resAmpl() + DK2piRes10.resAmpl();
843}
844
845//
846// BaBar decay amplitudes for D0->Ks K+ K-
847//
848// p4_p is D0
849// moms1 is K0s
850// moms2 is K+
851// moms3 is K-
852// Amplitudes and phases are taken from BaBar hep-ex/0207089
853// with convention : Non Resonant = Amp 1. / Phase 0.
854
855EvtComplex EvtDDalitz::amplDtoK0KK(EvtVector4R p4_p, EvtVector4R moms1,
856 EvtVector4R moms2, EvtVector4R moms3) {
857
858 //phi
859 EvtResonance DK0KKRes1( p4_p, moms2, moms3, 113.75, -40.0, 0.0043,
860 1.019456, 1 ) ;
861 //a0(980)
862 EvtResonance DK0KKRes2( p4_p, moms2, moms3, 152.25, 69.0, 0.1196 , 0.9847,
863 0 ) ;
864 //f0(980)
865 EvtResonance DK0KKRes3( p4_p, moms2, moms3, 30.5, -201.0, 0.05, 0.980 ,
866 0 ) ;
867 //a0(980)+
868 EvtResonance DK0KKRes4( p4_p, moms1, moms2, 85.75, -93.0, 0.1196 , 0.9847,
869 0 ) ;
870 //a0(980)-
871 EvtResonance DK0KKRes5( p4_p, moms3, moms1, 8. , -53.0 ,0.1196, 0.9847,
872 0 ) ;
873
874
875 return EvtComplex(1.0,0.0) +
876
877 DK0KKRes1.resAmpl() + DK0KKRes2.resAmpl() +
878
879 DK0KKRes3.resAmpl() + DK0KKRes4.resAmpl() +
880
881 DK0KKRes5.resAmpl() ;
882
883
884}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
delta
const double delta
Definition: CheckRung/check/gainVsDoca.cxx:32
report
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
ERROR
@ ERROR
Definition: EvtReport.hh:49
EvtConst::radToDegrees
static const double radToDegrees
Definition: EvtConst.hh:30
EvtCubicSpline::resAmpl
EvtComplex resAmpl()
Definition: EvtCubicSpline.cc:195
EvtCubicSpline::setParams
static void setParams(const vector< double > x, const vector< double > ym, const vector< double > yp)
Definition: EvtCubicSpline.cc:71
EvtDDalitz::clone
EvtDecayBase * clone()
Definition: EvtDDalitz.cc:46
EvtDDalitz::initProbMax
void initProbMax()
Definition: EvtDDalitz.cc:370
EvtDDalitz::getName
void getName(std::string &name)
Definition: EvtDDalitz.cc:39
EvtDDalitz::decay
void decay(EvtParticle *p)
Definition: EvtDDalitz.cc:384
EvtDDalitz::~EvtDDalitz
virtual ~EvtDDalitz()
Definition: EvtDDalitz.cc:37
EvtDDalitz::init
void init()
Definition: EvtDDalitz.cc:52
EvtDecayAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtDecayAmp.hh:37
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:533
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::getName
virtual void getName(std::string &name)=0
EvtDecayBase::getParentId
EvtId getParentId()
Definition: EvtDecayBase.hh:60
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
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
EvtDecayTable::getDecayFunc
static EvtDecayBase * getDecayFunc(EvtParticle *)
Definition: EvtDecayTable.cc:73
EvtFlatte::resAmpl
EvtComplex resAmpl()
Definition: EvtFlatte.cc:71
EvtId
Definition: EvtId.hh:27
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:112
EvtParticle::getParent
EvtParticle * getParent()
Definition: EvtParticle.cc:86
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:126
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtResonance2::resAmpl
EvtComplex resAmpl()
Definition: EvtResonance2.cc:61
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:183