CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
EvtPsi3Sdecay Class Reference

#include <EvtPsi3Sdecay.hh>

Public Member Functions

 EvtPsi3Sdecay (double ecms, EvtParticle *parent)
 
 EvtPsi3Sdecay ()
 
virtual ~EvtPsi3Sdecay ()
 
int findPoints ()
 
double polint (std::vector< double > points)
 
bool choseDecay ()
 
EvtParticle * choseDecay (EvtParticle *par)
 
int getDecay (double ecms)
 
double theProb (std::vector< double > myxs, int ich)
 
int findMode ()
 
int getMode ()
 
std::vector< EvtId > getVId (int mode)
 
void PHSPDecay (EvtParticle *par)
 
std::vector< EvtId > getDaugId ()
 
std::vector< EvtVector4R > getDaugP4 ()
 
bool AngSam (EvtVector4R parent_p4cm, EvtVector4R son_p4cm, double alpha)
 
void setMode (int m)
 

Detailed Description

Definition at line 33 of file EvtPsi3Sdecay.hh.

Constructor & Destructor Documentation

◆ EvtPsi3Sdecay() [1/2]

EvtPsi3Sdecay::EvtPsi3Sdecay ( double  ecms,
EvtParticle *  parent 
)
inline

Definition at line 36 of file EvtPsi3Sdecay.hh.

36 { //for 2-body decays
37 //initializer
38 Ecms = ecms;
39 theParent = parent;
40 Ndaugs=parent->getNDaug();
41 nsize = 24;
42
43 _excflag=0;
44 x.clear();
45 // open charm cross section, see PRD 80, 072001, xs in unit pb.
46 double xx[24]={3.72968, 3.97, 3.99, 4.01, 4.015, 4.03, 4.06, 4.12, 4.14, 4.16, 4.17, 4.18, 4.2, 4.26, 4.30, 4.34, 4.38, 4.42, 4.46, 4.50, 4.54, 4.58, 4.62, 4.66 }; // 14 energy points
47 double y0[24]={0., 86, 133, 76, 10, 334, 410, 303, 177, 167, 177, 179, 180, 86, 31, 49, 65, 196, 52, 87, 166, 14, 33, 49};// 0) D0D0bar cross section
48 double y1[24]={0., 137, 90, 135, 38, 196, 480, 310, 200, 200, 182, 197, 181, 94, 108, 96, 154, 165, 171, 106, 27, 144, 36, 22}; // 1) D+D-
49 double y2[24]={0., 1140, 1370, 1660, 1920, 1600, 1115, 700, 675, 626, 636, 605, 515, 540, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 2)D0D*0bar
50 double y3[24]={0., 1140, 1370, 1660, 1920, 1600, 1115, 700, 675, 626, 636, 605, 515, 540, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 3)D0bar D*0
51 double y4[24]={0., 0, 0, 0, 213, 2000, 2290, 2550, 2443, 2566, 2363, 2173, 1830, 269, 822, 1045, 1020, 820, 398, 588, 690, 459, 360, 398};// 4)D*0 D*0bar
52 double y5[24]={0., 1115, 1375, 1650, 1851, 1650, 1085, 780, 688, 688, 642, 648, 535, 511, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 5)D*+D-
53 double y6[24]={0., 1115, 1375, 1650, 1851, 1650, 1085, 780, 688, 688, 642, 648, 535, 511, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 6)D*-D+
54 double y7[24]={0., 0, 0, 0, 0, 1400, 2390, 2280, 2556, 2479, 2357, 2145, 1564, 237, 822, 1045, 1020, 820, 398, 588, 690, 459, 360, 398};// 7)D*+D*-
55 double y8[24]={0., 102, 133, 269, 250, 174, 51, 26, 25, 15, 34, 7, 15, 47, 106, 70, 36, 10, 2, 28, 60, 60, 48, 36}; // 8)Ds+ Ds-
56 double y9[24]={0., 0, 0, 0, 0, 0, 0, 478, 684, 905, 916, 889, 812, 34, 314, 368, 318, 357, 292, 171, 66, 103, 190, 272};// 9)Ds*+ Ds-
57 double y10[24]={0., 0, 0, 0, 0, 0, 0, 478, 684, 905, 916, 889, 812, 34, 314, 368, 318, 357, 292, 171, 66, 103, 190, 272};//10)Ds*- Ds+
58 double y11[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 440, 398, 428, 310, 131, 0, 45, 126, 98, 39, 0}; //11)Ds*+ Ds*-
59 double y12[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6, 49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 };//12)D*+ D- pi0 //------ DD* pi----
60 double y13[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6, 49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //13)D*- D+ pi0
61 double y14[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //14)D*+ anti-D0 pi-
62 double y15[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //15)D*- D0 pi+
63 double y16[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //16)D+ anti-D*0 pi-
64 double y17[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //17)D- D*0 pi+
65 double y18[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //18) D*+ D*- pi0 //------D*D*pi, above 4.26Gev, assumed xs as D*D pi
66 double y19[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,107.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110 }; //19) anti-D*0 D*+ pi-
67 double y20[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,107.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110 }; //20) D*0 D*- pi+
68 double y21[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //21) D*0 D*0bar pi0
69 double y22[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //22)D*0 D0bar pi0 //------ DD* pi----
70 double y23[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //23)D*bar D0 pi0
71 d0d0bar.clear();
72 dpdm.clear();
73 d0dst0bar.clear();
74 dst0dst0bar.clear();
75 d0bardst0.clear();
76 dstpdm.clear();
77 dstmdp.clear();
78 dstpdstm.clear();
79 dspdsm.clear();
80 dsspdsm.clear();
81 dssmdsp.clear();
82 dsspdssm.clear();
83 xs12.clear();
84 xs13.clear();
85 xs14.clear();
86 xs15.clear();
87 xs16.clear();
88 xs17.clear();
89 xs18.clear();
90 xs19.clear();
91 xs20.clear();
92 xs21.clear();
93 xs22.clear();
94 xs23.clear();
95
96 for(int i=0;i<24;i++){
97 x.push_back(xx[i]);
98 d0d0bar.push_back(y0[i]);
99 dpdm.push_back(y1[i]);
100 d0dst0bar.push_back(y2[i]);
101 d0bardst0.push_back(y3[i]);
102 dst0dst0bar.push_back(y4[i]);
103 dstpdm.push_back( y5[i]);
104 dstmdp.push_back( y6[i]);
105 dstpdstm.push_back(y7[i]);
106 dspdsm.push_back( y8[i]);
107 dsspdsm.push_back( y9[i]);
108 dssmdsp.push_back( y10[i]);
109 dsspdssm.push_back( y11[i]);
110 xs12.push_back( y12[i] );
111 xs13.push_back( y13[i] );
112 xs14.push_back( y14[i] );
113 xs15.push_back( y15[i] );
114 xs16.push_back( y16[i] );
115 xs17.push_back( y17[i] );
116 xs18.push_back( y18[i] );
117 xs19.push_back( y19[i] );
118 xs20.push_back( y20[i] );
119 xs21.push_back( y21[i] );
120 xs22.push_back( y22[i] );
121 xs23.push_back( y23[i] );
122 }
123 }
EvtParticle::getNDaug
int getNDaug() const
Definition: EvtParticle.cc:125
ecms
const double ecms
Definition: inclkstar.cxx:42

◆ EvtPsi3Sdecay() [2/2]

EvtPsi3Sdecay::EvtPsi3Sdecay ( )
inline

Definition at line 126 of file EvtPsi3Sdecay.hh.

126 {//for 2-body and 3-body decays
127 //initializer
128 // Ecms = ecms;
129 nsize = 24;
130
131 x.clear();
132 // open charm cross section, see PRD 80, 072001, xross section in pb
133 double xx[24]={3.72968, 3.97, 3.99, 4.01, 4.015, 4.03, 4.06, 4.12, 4.14, 4.16, 4.17, 4.18, 4.2, 4.26, 4.30, 4.34, 4.38, 4.42, 4.46, 4.50, 4.54, 4.58, 4.62, 4.66 }; // 14 energy points
134 double y0[24]={0., 86, 133, 76, 10, 334, 410, 303, 177, 167, 177, 179, 180, 86, 31, 49, 65, 196, 52, 87, 166, 14, 33, 49};// 0) D0D0bar cross section
135 double y1[24]={0., 137, 90, 135, 38, 196, 480, 310, 200, 200, 182, 197, 181, 94, 108, 96, 154, 165, 171, 106, 27, 144, 36, 22}; // 1) D+D-
136 double y2[24]={0., 1140, 1370, 1660, 1920, 1600, 1115, 700, 675, 626, 636, 605, 515, 540, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 2)D0D*0bar
137 double y3[24]={0., 1140, 1370, 1660, 1920, 1600, 1115, 700, 675, 626, 636, 605, 515, 540, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 3)D0bar D*0
138 double y4[24]={0., 0, 0, 0, 213, 2000, 2290, 2550, 2443, 2566, 2363, 2173, 1830, 269, 822, 1045, 1020, 820, 398, 588, 690, 459, 360, 398};// 4)D*0 D*0bar
139 double y5[24]={0., 1115, 1375, 1650, 1851, 1650, 1085, 780, 688, 688, 642, 648, 535, 511, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 5)D*+D-
140 double y6[24]={0., 1115, 1375, 1650, 1851, 1650, 1085, 780, 688, 688, 642, 648, 535, 511, 748, 880, 556, 657, 477, 494, 320, 616, 575, 373};// 6)D*-D+
141 double y7[24]={0., 0, 0, 0, 0, 1400, 2390, 2280, 2556, 2479, 2357, 2145, 1564, 237, 822, 1045, 1020, 820, 398, 588, 690, 459, 360, 398};// 7)D*+D*-
142 double y8[24]={0., 102, 133, 269, 250, 174, 51, 26, 25, 15, 34, 7, 15, 47, 106, 70, 36, 10, 2, 28, 60, 60, 48, 36}; // 8)Ds+ Ds-
143 double y9[24]={0., 0, 0, 0, 0, 0, 0, 478, 684, 905, 916, 889, 812, 34, 314, 368, 318, 357, 292, 171, 66, 103, 190, 272};// 9)Ds*+ Ds-
144 double y10[24]={0., 0, 0, 0, 0, 0, 0, 478, 684, 905, 916, 889, 812, 34, 314, 368, 318, 357, 292, 171, 66, 103, 190, 272};//10)Ds*- Ds+
145 double y11[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 440, 398, 428, 310, 131, 0, 45, 126, 98, 39, 0}; //11)Ds*+ Ds*-
146 double y12[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6, 49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 };//12)D*+ D- pi0 //------ DD* pi----
147 double y13[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6, 49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //13)D*- D+ pi0
148 double y14[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //14)D*+ anti-D0 pi-
149 double y15[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //15)D*- D0 pi+
150 double y16[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //16)D+ anti-D*0 pi-
151 double y17[24]={0., 0, 0, 0, 0, 0, 24, 7.5, 68.6, 64.8, 73.3, 95.8, 122.5,106.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110}; //17)D- D*0 pi+
152 double y18[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //18) D*+ D*- pi0 //------D*D*pi, above 4.26Gev, assumed xs as D*D pi
153 double y19[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,107.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110 }; //19) anti-D*0 D*+ pi-
154 double y20[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,107.3,98, 378, 650, 564, 167, 304, 359, 381, 497, 110 }; //20) D*0 D*- pi+
155 double y21[24]={0., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //21) D*0 D*0bar pi0
156 double y22[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //22)D*0 D0bar pi0 //------ DD* pi----
157 double y23[24]={0., 0, 0, 0, 0, 0, 12, 3.8, 32.4, 32.4, 37, 48, 61.2, 53.6,49, 169, 325, 282, 83.5, 152, 179, 190, 248, 55 }; //23)D*bar D0 pi0
158
159 d0d0bar.clear();
160 dpdm.clear();
161 d0dst0bar.clear();
162 dst0dst0bar.clear();
163 d0bardst0.clear();
164 dstpdm.clear();
165 dstmdp.clear();
166 dstpdstm.clear();
167 dspdsm.clear();
168 dsspdsm.clear();
169 dssmdsp.clear();
170 dsspdssm.clear();
171 xs12.clear();
172 xs13.clear();
173 xs14.clear();
174 xs15.clear();
175 xs16.clear();
176 xs17.clear();
177 xs18.clear();
178 xs19.clear();
179 xs20.clear();
180 xs21.clear();
181 xs22.clear();
182 xs23.clear();
183 for(int i=0;i<24;i++){
184 x.push_back(xx[i]);
185 d0d0bar.push_back(y0[i]);
186 dpdm.push_back(y1[i]);
187 d0dst0bar.push_back(y2[i]);
188 d0bardst0.push_back(y3[i]);
189 dst0dst0bar.push_back(y4[i]);
190 dstpdm.push_back( y5[i]);
191 dstmdp.push_back( y6[i]);
192 dstpdstm.push_back(y7[i]);
193 dspdsm.push_back( y8[i]);
194 dsspdsm.push_back( y9[i]);
195 dssmdsp.push_back( y10[i]);
196 dsspdssm.push_back( y11[i]);
197 xs12.push_back( y12[i] );
198 xs13.push_back( y13[i] );
199 xs14.push_back( y14[i] );
200 xs15.push_back( y15[i] );
201 xs16.push_back( y16[i] );
202 xs17.push_back( y17[i] );
203 xs18.push_back( y18[i] );
204 xs19.push_back( y19[i] );
205 xs20.push_back( y20[i] );
206 xs21.push_back( y21[i] );
207 xs22.push_back( y22[i] );
208 xs23.push_back( y23[i] );
209 }
210
211//---- initilize Vmode
212
213 VmodeSon.clear();
214 //0: D0 anti-D0
215 Vson.clear();
216 Vson.push_back("D0"); Vson.push_back("anti-D0");
217 VmodeSon.push_back(Vson);
218
219 //1: D+ D-
220 Vson.clear();
221 Vson.push_back("D+"); Vson.push_back("D-");
222 VmodeSon.push_back(Vson);
223
224 //2: D0 anti-D*0
225 Vson.clear();
226 Vson.push_back("D0"); Vson.push_back("anti-D*0");
227 VmodeSon.push_back(Vson);
228
229 //3: anti-D0 D*0
230 Vson.clear();
231 Vson.push_back("anti-D0"); Vson.push_back("D*0");
232 VmodeSon.push_back(Vson);
233
234 //4: D*0 anti-D*0
235 Vson.clear();
236 Vson.push_back("D*0"); Vson.push_back("anti-D*0");
237 VmodeSon.push_back(Vson);
238
239 //5: D*+ D-
240 Vson.clear();
241 Vson.push_back("D*+"); Vson.push_back("D-");
242 VmodeSon.push_back(Vson);
243
244 //6: D*- D+
245 Vson.clear();
246 Vson.push_back("D*-"); Vson.push_back("D+");
247 VmodeSon.push_back(Vson);
248
249 //7: D*+ D*-
250 Vson.clear();
251 Vson.push_back("D*+"); Vson.push_back("D*-");
252 VmodeSon.push_back(Vson);
253
254 //8: D_s+ D_s-
255 Vson.clear();
256 Vson.push_back("D_s+"); Vson.push_back("D_s-");
257 VmodeSon.push_back(Vson);
258
259 //9: D_s*+ D_s-
260 Vson.clear();
261 Vson.push_back("D_s*+"); Vson.push_back("D_s-");
262 VmodeSon.push_back(Vson);
263
264 //10: D_s*- D_s+
265 Vson.clear();
266 Vson.push_back("D_s*-"); Vson.push_back("D_s+");
267 VmodeSon.push_back(Vson);
268
269 //11: D_s*+ D_s*-
270 Vson.clear();
271 Vson.push_back("D_s*+"); Vson.push_back("D_s*-");
272 VmodeSon.push_back(Vson);
273
274 //12: D*+ D- pi0
275 Vson.clear();
276 Vson.push_back("D*+"); Vson.push_back("D-");Vson.push_back("pi0");
277 VmodeSon.push_back(Vson);
278
279 //13: D*- D+ pi0
280 Vson.clear();
281 Vson.push_back("D*-"); Vson.push_back("D+");Vson.push_back("pi0");
282 VmodeSon.push_back(Vson);
283
284 //14: D*+ anti-D0 pi-
285 Vson.clear();
286 Vson.push_back("D*+"); Vson.push_back("anti-D0");Vson.push_back("pi-");
287 VmodeSon.push_back(Vson);
288
289 //15: D*- D0 pi+
290 Vson.clear();
291 Vson.push_back("D*-"); Vson.push_back("D0");Vson.push_back("pi+");
292 VmodeSon.push_back(Vson);
293
294 //16: D+ anti-D*0 pi-
295 Vson.clear();
296 Vson.push_back("D+"); Vson.push_back("anti-D*0");Vson.push_back("pi-");
297 VmodeSon.push_back(Vson);
298
299 //17: D- D*0 pi+
300 Vson.clear();
301 Vson.push_back("D-"); Vson.push_back("D*0");Vson.push_back("pi+");
302 VmodeSon.push_back(Vson);
303
304 //18: D*+ D*- pi0
305 Vson.clear();
306 Vson.push_back("D*+"); Vson.push_back("D*-");Vson.push_back("pi0");
307 VmodeSon.push_back(Vson);
308
309 //19: anti-D*0 D*+ pi-
310 Vson.clear();
311 Vson.push_back("anti-D*0"); Vson.push_back("D*+");Vson.push_back("pi-");
312 VmodeSon.push_back(Vson);
313
314 //20: D*0 D*- pi+
315 Vson.clear();
316 Vson.push_back("D*0"); Vson.push_back("D*-");Vson.push_back("pi+");
317 VmodeSon.push_back(Vson);
318
319 //21: D*0 D*0bar pi0
320 Vson.clear();
321 Vson.push_back("D*0"); Vson.push_back("anti-D*0");Vson.push_back("pi0");
322 VmodeSon.push_back(Vson);
323
324 //22: D0bar D*0 pi0
325 Vson.clear();
326 Vson.push_back("anti-D0"); Vson.push_back("D*0");Vson.push_back("pi0");
327 VmodeSon.push_back(Vson);
328
329 //23: D*0bar D0 pi0
330 Vson.clear();
331 Vson.push_back("anti-D*0"); Vson.push_back("D0");Vson.push_back("pi0");
332 VmodeSon.push_back(Vson);
333
334 }

◆ ~EvtPsi3Sdecay()

virtual EvtPsi3Sdecay::~EvtPsi3Sdecay ( )
inlinevirtual

Definition at line 337 of file EvtPsi3Sdecay.hh.

337{}

Member Function Documentation

◆ AngSam()

bool EvtPsi3Sdecay::AngSam ( EvtVector4R  parent_p4cm,
EvtVector4R  son_p4cm,
double  alpha 
)

Definition at line 412 of file EvtPsi3Sdecay.cc.

412 {
413 EvtHelSys angles(parent_p4cm,son_p4cm);
414 double theta=angles.getHelAng(1);
415 //double phi =angles.getHelAng(2);
416 //double gamma=0;
417 double costheta=cos(theta); //using helicity angles in parent system
418 double max_alpha;
419 if(alpha>=0) {max_alpha = 1+alpha;}else
420 {max_alpha=1;}
421 double ags = (1+alpha*costheta*costheta)/max_alpha;
422 double rand=EvtRandom::Flat(0.0, 1.0);
423 if(rand <=ags) {return true;}
424 else {return false;}
425}
alpha
const double alpha
Definition: FastVertexFit.cxx:4
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:73

Referenced by PHSPDecay().

◆ choseDecay() [1/2]

bool EvtPsi3Sdecay::choseDecay ( )

Definition at line 127 of file EvtPsi3Sdecay.cc.

127 { //determing accept or reject a generated decay
128
129 // findPoints();
130 double d0d0bar_xs=polint(d0d0bar);
131 double dpdm_xs = polint(dpdm);
132 double d0dst0bar_xs = polint(d0dst0bar);
133 double d0bardst0_xs = polint(d0bardst0);
134
135 double dst0dst0bar_xs = polint(dst0dst0bar);
136 double dstpdm_xs = polint(dstpdm);
137
138 double dstmdp_xs = polint(dstmdp);
139 double dstpdstm_xs = polint(dstpdstm);
140
141 double dspdsm_xs = polint(dspdsm);
142
143 double dsspdsm_xs = polint(dsspdsm);
144 double dssmdsp_xs = polint(dssmdsp);
145
146 double dsspdssm_xs = polint(dsspdssm);
147 //--- DDpi modes
148 double _xs12 = polint(xs12);
149 double _xs13 = polint(xs13);
150 double _xs14 = polint(xs14);
151 double _xs15 = polint(xs15);
152 double _xs16 = polint(xs16);
153 double _xs17 = polint(xs17);
154 double _xs18 = polint(xs18);
155 double _xs19 = polint(xs19);
156 double _xs20 = polint(xs20);
157 double _xs21 = polint(xs21);
158 double _xs22 = polint(xs22);
159 double _xs23 = polint(xs23);
160
161 int ich = findMode();
162 // std::cout<<"calculated XS "<< d0d0bar_xs<<" "<<dpdm_xs<<" "<<d0dst0bar_xs<<" "<<d0bardst0_xs<< " "<<dst0dst0bar_xs<<" "<<dstpdm_xs<< " "<<dstmdp_xs<<" "<<dstpdstm_xs<<" "<<dspdsm_xs<<std::endl;
163
164 double xmtotal=0;
165 for(int i=0;i<Vid.size();i++){
166 xmtotal += EvtPDL::getMinMass(Vid[i]);
167 }
168 double mparent= theParent->getP4().mass();
169 // std::cout<<"mparent= "<<mparent<<", xmtotal= "<<xmtotal<<endl;
170 if (mparent<xmtotal){return false;}
171
172
173 std::vector<double> myxs; myxs.clear();
174 myxs.push_back(d0d0bar_xs); //0
175 myxs.push_back(dpdm_xs);
176 myxs.push_back(d0dst0bar_xs); //2
177 myxs.push_back(d0bardst0_xs);
178 myxs.push_back(dst0dst0bar_xs); //4
179 myxs.push_back(dstpdm_xs);
180 myxs.push_back(dstmdp_xs); //6
181 myxs.push_back(dstpdstm_xs);
182 myxs.push_back(dspdsm_xs); //8
183 myxs.push_back(dsspdsm_xs);
184 myxs.push_back(dssmdsp_xs); //10
185 myxs.push_back(dsspdssm_xs); //11
186
187 myxs.push_back(_xs12); //12
188 myxs.push_back(_xs13); //13
189 myxs.push_back(_xs14); //14
190 myxs.push_back(_xs15); //15
191 myxs.push_back(_xs16); //16
192 myxs.push_back(_xs17); //17
193 myxs.push_back(_xs18); //18
194 myxs.push_back(_xs19); //19
195 myxs.push_back(_xs20); //20
196 myxs.push_back(_xs21); //18
197 myxs.push_back(_xs22); //19
198 myxs.push_back(_xs23); //20
199
200 double Prop0,Prop1;
201 if(ich==0){ Prop0=0;} else
202 {
203 Prop0 = theProb(myxs,ich-1);
204 }
205 Prop1 = theProb(myxs,ich);
206
207 double pm= EvtRandom::Flat(0.,1);
208 bool flag = false;
209 if( Prop0 < pm && pm<= Prop1 ) flag = true;
210
211 //--- debuging
212
213 if(flag) {
214 // std::cout<<"findMode= "<<ich<<std::endl;
215 // for(int i=0;i<myxs.size();i++){ std::cout<<"channel "<<i<<" myxs: "<<myxs[i]<<std::endl;}
216 //std::cout<<"prop0,prop1= "<<Prop0<<" "<<Prop1<<std::endl;
217 }
218
219 //-------------
220
221 return flag;
222}
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition: EvtPDL.hh:51
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
EvtPsi3Sdecay::theProb
double theProb(std::vector< double > myxs, int ich)
Definition: EvtPsi3Sdecay.cc:55
EvtPsi3Sdecay::polint
double polint(std::vector< double > points)
Definition: EvtPsi3Sdecay.cc:42
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39

◆ choseDecay() [2/2]

EvtParticle * EvtPsi3Sdecay::choseDecay ( EvtParticle *  par)

Definition at line 225 of file EvtPsi3Sdecay.cc.

225 { //decay par
226 double xm = par->mass();
227 int themode = getDecay(xm);
228 std::vector< EvtId > theid = getVId(themode);
229 int ndaugjs = theid.size();
230 EvtId myId[3];
231 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
232 par->makeDaughters(ndaugjs,myId);
233
234 for(int i=0;i<par->getNDaug();i++){
235 EvtParticle* di=par->getDaug(i);
236 double xmi=EvtPDL::getMinMass(di->getId());
237 di->setMass(xmi);
238 }
239 par->initializePhaseSpace(ndaugjs,myId);
240 _themode = themode;
241 return par;
242}
EvtId
Definition: EvtId.hh:27
EvtParticle::setMass
void setMass(double m)
Definition: EvtParticle.hh:358
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1175
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1069
EvtPsi3Sdecay::getDecay
int getDecay(double ecms)
Definition: EvtPsi3Sdecay.cc:307
EvtPsi3Sdecay::getVId
std::vector< EvtId > getVId(int mode)
Definition: EvtPsi3Sdecay.cc:401

◆ findMode()

int EvtPsi3Sdecay::findMode ( )

Definition at line 73 of file EvtPsi3Sdecay.cc.

73 {
74 // mode index: 0) D0D0bar, 1)D+D-; 2)D0D*0bar , 3)D0bar D*0, 4)D*0 D*0 bar, 5)D*+D- 6)D*-D+ 7)D*+ D*- 8) Ds+ Ds-
75 // more modes: 9) D_s^*+ D_s^- ;10) D_s^*- D_s^+; 11) D_s^*+ D_s^*-
76
77 std::string son0,son1,son2;
78 Vson.clear();
79 Vid.clear();
80 for(int i=0;i<Ndaugs;i++){
81 std::string nson=EvtPDL::name(theParent->getDaug(i)->getId());
82 if(nson!="gammaFSR" && nson!="gamma"){ Vson.push_back(nson);Vid.push_back(theParent->getDaug(i)->getId());}
83 }
84 int nh=Vson.size();
85 //debugging
86
87 if(nh == 2){
88 // std::cout<<"final states: "<<Vson[0]<<" "<<Vson[1]<<std::endl;
89 son0 = Vson[0];son1 = Vson[1];
90 if(son0 == "D0" && son1 == "anti-D0" || son1 == "D0" && son0 == "anti-D0") {return 0;} else
91 if(son0 == "D+" && son1 == "D-" || son1 == "D+" && son0 == "D-" ) {return 1;} else
92 if(son0 == "D0" && son1 == "anti-D*0" || son1 == "D0" && son0 == "anti-D*0") {return 2;} else
93 if(son0 == "anti-D0" && son1 == "D*0" || son1 == "anti-D0" && son0 == "D*0") {return 3;} else
94 if(son0 == "D*0" && son1 == "anti-D*0" || son1 == "D*0" && son0 == "anti-D*0") {return 4;} else
95 if(son0 == "D*+" && son1 == "D-" || son1 == "D*+" && son0 == "D-") {return 5;} else
96 if(son0 == "D*-" && son1 == "D+" || son1 == "D*-" && son0 == "D+") {return 6;} else
97 if(son0 == "D*+" && son1 == "D*-" || son1 == "D*+" && son0 == "D*-") {return 7;} else
98 if(son0 == "D_s+" && son1 == "D_s-" || son1 == "D_s+" && son0 == "D_s-") {return 8;} else
99 if(son0 == "D_s*+" && son1 == "D_s-" || son1 == "D_s*+" && son0 == "D_s-") {return 9;} else
100 if(son0 == "D_s*-" && son1 == "D_s+" || son1 == "D_s*-" && son0 == "D_s+") {return 10;}else
101 if(son0 == "D_s*+" && son1 == "D_s*-" || son1 == "D_s*+" && son0 == "D_s*-") {return 11;}
102 } else if(nh == 3){
103 son0 = Vson[0];son1 = Vson[1];son2 = Vson[2];
104 if(son0 == "D*+" && son1 == "D-" && son2 == "pi0" ) {return 12;} else
105 if(son0 == "D*-" && son1 == "D+" && son2 == "pi0" ) {return 13;} else
106 if(son0 == "D*+" && son1 == "anti-D0" && son2 == "pi-" ) {return 14;} else
107 if(son0 == "D*-" && son1 == "D0" && son2 == "pi+" ) {return 15;} else
108 if(son0 == "D+" && son1 == "anti-D*0" && son2 == "pi-" ) {return 16;} else
109 if(son0 == "D-" && son1 == "D*0" && son2 == "pi+" ) {return 17;} else
110 if(son0 == "D*+" && son1 == "D*-" && son2 == "pi0" ) {return 18;} else
111 if(son0 == "anti-D*0" &&son1 == "D*+" && son2 == "pi-" ) {return 19;} else
112 if(son0 == "D*0" && son1 == "D*-" && son2 == "pi+" ) {return 20;} else
113 if(son0 == "D*0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 21;} else
114 if(son0 == "D0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 22;} else
115 if(son0 == "anti-D0" && son1 == "D*0" && son2 == "pi0" ) {return 23;}
116 }else{
117 std::cout<<"Not open charm decay"<<std::endl;
118 for(int j=0;j<nh;j++){
119 std::cout<<"final states, J= "<<j<<" "<<Vson[j]<<std::endl;
120 }
121 ::abort();
122 }
123
124}
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.hh:64

Referenced by choseDecay().

◆ findPoints()

int EvtPsi3Sdecay::findPoints ( )

Definition at line 32 of file EvtPsi3Sdecay.cc.

32 {
33 if(Ecms < x[0] ){theLocation =0;} else
34 if(Ecms >=x[nsize-1]) {theLocation = nsize-1;} else{
35 for (int i=0;i<nsize-1;i++){
36 if( x[i] <= Ecms && x[i+1] > Ecms) {theLocation=i;}
37 }
38 }
39 return theLocation;
40}

Referenced by polint().

◆ getDaugId()

std::vector< EvtId > EvtPsi3Sdecay::getDaugId ( )
inline

Definition at line 351 of file EvtPsi3Sdecay.hh.

351{return Vid;}

Referenced by EvtOpenCharm::decay().

◆ getDaugP4()

std::vector< EvtVector4R > EvtPsi3Sdecay::getDaugP4 ( )
inline

Definition at line 352 of file EvtPsi3Sdecay.hh.

352{return v_p4;}

Referenced by EvtOpenCharm::decay().

◆ getDecay()

int EvtPsi3Sdecay::getDecay ( double  ecms)

Definition at line 307 of file EvtPsi3Sdecay.cc.

307 { //pick up a decay by the accept-reject sampling method
308 if(ecms<3.97 || ecms >4.66){std::cout<<"EvtPsi3Sdecay::getDecay: You need to set the CMS energy between 3.97~4.66 GeV, but you have ecms= "<<ecms<< " GeV.The lower end can be set at KKMC"<<std::endl; }
309 if(_excflag ==1) return _themode;
310 Ecms = ecms;
311 // findPoints();
312 double d0d0bar_xs=polint(d0d0bar);
313 double dpdm_xs = polint(dpdm);
314 double d0dst0bar_xs = polint(d0dst0bar);
315 double d0bardst0_xs = polint(d0bardst0);
316
317 double dst0dst0bar_xs = polint(dst0dst0bar);
318 double dstpdm_xs = polint(dstpdm);
319
320 double dstmdp_xs = polint(dstmdp);
321 double dstpdstm_xs = polint(dstpdstm);
322
323 double dspdsm_xs = polint(dspdsm);
324
325 double dsspdsm_xs = polint(dsspdsm);
326 double dssmdsp_xs = polint(dssmdsp);
327
328 double dsspdssm_xs = polint(dsspdssm);
329
330 double _xs12 = polint(xs12);
331 double _xs13 = polint(xs13);
332 double _xs14 = polint(xs14);
333 double _xs15 = polint(xs15);
334 double _xs16 = polint(xs16);
335 double _xs17 = polint(xs17);
336 double _xs18 = polint(xs18);
337 double _xs19 = polint(xs19);
338 double _xs20 = polint(xs20);
339 double _xs21 = polint(xs21);
340 double _xs22 = polint(xs22);
341 double _xs23 = polint(xs23);
342
343
344 std::vector<double> myxs; myxs.clear();
345 myxs.push_back(d0d0bar_xs); //0
346 myxs.push_back(dpdm_xs); //1
347 myxs.push_back(d0dst0bar_xs); //2
348 myxs.push_back(d0bardst0_xs); //3
349 myxs.push_back(dst0dst0bar_xs);//4
350 myxs.push_back(dstpdm_xs); //5
351 myxs.push_back(dstmdp_xs); //6
352 myxs.push_back(dstpdstm_xs); //7
353 myxs.push_back(dspdsm_xs); //8
354 myxs.push_back(dsspdsm_xs); //9
355 myxs.push_back(dssmdsp_xs); //10
356 myxs.push_back(dsspdssm_xs); //11
357 myxs.push_back(_xs12); //12
358 myxs.push_back(_xs13); //13
359 myxs.push_back(_xs14); //14
360 myxs.push_back(_xs15); //15
361 myxs.push_back(_xs16); //16
362 myxs.push_back(_xs17); //17
363 myxs.push_back(_xs18); //18
364 myxs.push_back(_xs19); //19
365 myxs.push_back(_xs20); //20
366 myxs.push_back(_xs21); //21
367 myxs.push_back(_xs22); //22
368 myxs.push_back(_xs23); //23
369
370 int niter = 0;
371 loop:
372 int ich = (int)24*EvtRandom::Flat(0.,1);//sampling modes over 24 channel
373
374 niter++;
375 if(niter>1000) {std::cout<<"EvtPsi3Sdecay:getDecay() do loops over 1000 times at Ecm= "<<ecms<<std::endl;abort();}
376
377 double xmtotal=0;
378 std::vector<EvtId> theVid;
379 theVid.clear();
380 theVid = getVId(ich);
381 for(int i=0;i<theVid.size();i++){
382 xmtotal += EvtPDL::getMinMass(theVid[i]);
383 }
384 if(Ecms < xmtotal ) {goto loop;}
385
386 double Prop0,Prop1;
387 if(ich==0){ Prop0=0;} else
388 {
389 Prop0 = theProb(myxs,ich-1);
390 }
391 Prop1 = theProb(myxs,ich);
392
393 double pm= EvtRandom::Flat(0.,1);
394
395 if( Prop0 < pm && pm<= Prop1 ) {return ich;}
396 else {goto loop;}
397
398}

Referenced by choseDecay(), and PHSPDecay().

◆ getMode()

int EvtPsi3Sdecay::getMode ( )
inline

Definition at line 346 of file EvtPsi3Sdecay.hh.

346{return _themode;};

Referenced by EvtOpenCharm::decay().

◆ getVId()

std::vector< EvtId > EvtPsi3Sdecay::getVId ( int  mode)

Definition at line 401 of file EvtPsi3Sdecay.cc.

401 {
402 std::vector<EvtId> theVid;
403 theVid.clear();
404 for(int i=0;i<VmodeSon[mode].size();i++){
405 EvtId theId = EvtPDL::getId(VmodeSon[mode][i]);
406 theVid.push_back(theId);
407 }
408 return theVid;
409}
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287

Referenced by choseDecay(), getDecay(), and PHSPDecay().

◆ PHSPDecay()

void EvtPsi3Sdecay::PHSPDecay ( EvtParticle *  par)

Definition at line 247 of file EvtPsi3Sdecay.cc.

247 {//decay par
248 v_p4.clear();Vid.clear();
249 double xm = par->mass();
250 EvtVector4R p_ini(xm,0,0,0);
251 EvtParticle* mypar= EvtParticleFactory::particleFactory(par->getId(),p_ini);
252
253 int themode = getDecay(xm);
254 std::vector< EvtId > theid = getVId(themode);
255 int ndaugjs = theid.size();
256 EvtId myId[3];
257 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
258 mypar->makeDaughters(ndaugjs,myId);
259
260 for(int i=0;i<mypar->getNDaug();i++){
261 EvtParticle* di=mypar->getDaug(i);
262 double xmi=EvtPDL::getMinMass(di->getId());
263 di->setMass(xmi);
264 }
265 loop:
266 mypar->initializePhaseSpace(ndaugjs,myId);
267
268 //-- here to do amgular distribution sampling
269 bool pp = (themode == 0||themode == 1 ||themode ==6); //V->PP mode, alpha = -1
270 bool vp = (themode >=2 && themode <=5 || themode >=7 && themode <=9 ); // V->VP mode, alpha = 1
271 bool flag1 = false;
272 double alpha;
273 if(pp){alpha=-1;}else if(vp){alpha=1;} else {alpha=0;}
274 EvtVector4R pp4=par->getP4();
275 EvtVector4R sp4=mypar->getDaug(1)->getP4();
276 flag1=AngSam(pp4,sp4,alpha);
277 if(alpha != 0 && !flag1 ) goto loop;
278 //--
279 _themode = themode;
280 for(int i=0;i<mypar->getNDaug();i++){
281 EvtParticle* di=mypar->getDaug(i);
282 v_p4.push_back( di->getP4() );
283 Vid.push_back(myId[i]);
284 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
285 }
286
287
288 /*********** same function can be realized
289 double mass[3];
290 EvtVector4R p4[30];
291 for(int i=0;i<ndaugjs;i++){mass[i]=mypar->getDaug(i)->mass();}
292 EvtGenKine::PhaseSpace(ndaugjs,mass,p4,xm);
293 _themode = themode;
294 for(int i=0;i<mypar->getNDaug();i++){
295 v_p4.push_back( p4[i] );
296 Vid.push_back(myId[i]);
297 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
298 }
299 *************/
300 mypar->deleteTree();
301 return;
302}
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition: EvtParticleFactory.cc:46
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:555
EvtPsi3Sdecay::AngSam
bool AngSam(EvtVector4R parent_p4cm, EvtVector4R son_p4cm, double alpha)
Definition: EvtPsi3Sdecay.cc:412

Referenced by EvtOpenCharm::decay().

◆ polint()

double EvtPsi3Sdecay::polint ( std::vector< double >  points)

Definition at line 42 of file EvtPsi3Sdecay.cc.

42 {
43
44 theLocation = findPoints();
45 double xs;
46 if(theLocation==nsize-1){xs = vy[nsize-1];} else{
47 xs = (vy[theLocation +1 ]- vy[theLocation]) / (x[theLocation+1]-x[theLocation])*(Ecms - x[theLocation])+vy[theLocation];
48 }
49 if(xs <0 ) xs = 0;
50 // for(int i=0;i<vy.size();i++){ std::cout<<"channel "<<i<<" xsection: "<<vy[i]<<std::endl;}
51 // std::cout<<"Ecms, theLocation= "<<Ecms<<" "<<theLocation<<" xs= "<<xs<<std::endl;
52 return xs;
53 }

Referenced by choseDecay(), and getDecay().

◆ setMode()

void EvtPsi3Sdecay::setMode ( int  m)
inline

Definition at line 354 of file EvtPsi3Sdecay.hh.

354 {
355 if(m<0 || m>24) {std::cout<<"EvtPsi3Decay::setMode: bad mode"<<std::endl;abort();}
356 _themode = m;_excflag=1;}

Referenced by EvtOpenCharm::decay().

◆ theProb()

double EvtPsi3Sdecay::theProb ( std::vector< double >  myxs,
int  ich 
)

Definition at line 55 of file EvtPsi3Sdecay.cc.

55 {
56 int Nchannels=myxs.size();
57 //---debuging
58 // std::cout<<"Nchannel= "<<Nchannels<<endl;
59
60 std::vector <double> thexs;
61 thexs.clear();
62 double sumxs=0;
63 for(int j=0;j<Nchannels;j++){
64 sumxs += myxs[j];
65 thexs.push_back(sumxs);
66 //----debugging
67 // std::cout<<"thexs["<<j<<"]= "<<myxs[j]<<std::endl;
68 }
69 if(sumxs == 0) {std::cout<<"EvtPsi3Sdecay::theProb, denominator is 0"<<std::endl;::abort();}
70 return thexs[ich] / sumxs ;
71}

Referenced by choseDecay(), and getDecay().

The documentation for this class was generated from the following files: