BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0toKpiomegaPlot.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: EvtD0toKpiomegaPlot.cc
12//
13// Modification history:
14//
15// Liaoyuan Dong August, 2022 Module created
16//
17//------------------------------------------------------------------------
18#include "EvtD0toKpiomegaPlot.hh"
19#include "../EvtGenBase/EvtGenKine.hh"
20#include "../EvtGenBase/EvtPDL.hh"
21#include "../EvtGenBase/EvtParticle.hh"
22#include "../EvtGenBase/EvtPatches.hh"
23#include "../EvtGenBase/EvtRandom.hh"
24#include "../EvtGenBase/EvtReport.hh"
25#include <stdlib.h>
26#include <string>
27
28#include "TApplication.h"
29#include "TAxis.h"
30#include "TFile.h"
31#include "TH1.h"
32#include "TH2.h"
33#include "TROOT.h"
34using std::endl;
35
36EvtD0toKpiomegaPlot::~EvtD0toKpiomegaPlot() {}
37
38void EvtD0toKpiomegaPlot::getName( std::string& model_name ) {
39 model_name = "D0toKpiomegaPlot";
40}
41
42EvtDecayBase* EvtD0toKpiomegaPlot::clone() { return new EvtD0toKpiomegaPlot; }
43
44void EvtD0toKpiomegaPlot::init() {
45
46 checkNArg( 0 );
47
48 bool idN = getDaugs()[0] == EvtPDL::getId( std::string( "K+" ) ) ||
49 getDaugs()[0] == EvtPDL::getId( std::string( "K-" ) );
50 bool idKs = getDaugs()[1] == EvtPDL::getId( std::string( "pi+" ) ) ||
51 getDaugs()[1] == EvtPDL::getId( std::string( "pi-" ) );
52 bool idPi = getDaugs()[2] == EvtPDL::getId( std::string( "omega" ) );
53 if ( !( idN && idKs && idPi ) )
54 {
55 std::cout << "EvtD0toKpiomegaPlot: the daughter sequence should be K+ pi- omega"
56 << std::endl;
57 abort();
58 }
59
60 EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getParentId() );
61 /*
62 m_inputFileName=getenv("BESEVTGENROOT");
63 f1=m_inputFileName+"/src/EvtGen/EvtGenModels/"+"Body3root/EvtD0toKpiomegaPlot.root";
64
65 TFile f_1(f1.c_str()); //TFile f_4(f4.c_str());
66 TH2D *h_1 = (TH2D*)f_1.Get("h_13"); //m_12:m_13
67 TAxis* dx=h_1->GetXaxis();
68 TAxis* dy=h_1->GetYaxis();
69 int BINS1 =dx->GetLast();
70 int BINS2 =dy->GetLast();
71
72 cout << "X LowEdge " << dx->GetXmin() << endl;
73 cout << "X UpEdge " << dx->GetXmax() << endl;
74 cout << "X width " << dx->GetBinWidth(1) << " Nbin = " << dx->GetNbins() <<
75 endl; cout << "Y LowEdge " << dy->GetXmin() << endl; cout << "Y UpEdge " << dy->GetXmax()
76 << endl; cout << "Y width " << dy->GetBinWidth(1) << " Nbin = " << dy->GetNbins() <<
77 endl;
78
79 cout << "INIT EvtD0toKpiomegaPlot.cc BINS1, BINS2 = " << BINS1 << ", " <<
80 BINS2 << endl; double av1; avm1=0.0; for(int i=0;i<BINS1+2;i++){ for (int j=0; j<BINS2+2;
81 j++) { av1=h_1->GetBinContent(i,j); HisPDF[i][j] = av1; if(av1>avm1) avm1=av1;
82 }
83 }
84 std::cout<<"INIT EvtD0toKpiomegaPlot.cc avm1= " << avm1 << endl;
85
86 for(int i=0;i<BINS1+2;i++){
87 cout << "{";
88 for (int j=0; j<BINS2+2; j++) {
89 if (j==(BINS2+1)) cout << HisPDF[i][j] << "}," ;
90 else cout << HisPDF[i][j] << ", ";
91 }
92 cout << endl;
93 }
94 cout << endl;
95 */
96 Xmin = 0.400689;
97 Xmax = 1.17072;
98 Xwid = 0.0385018;
99 Ymin = 0.850084;
100 Ymax = 1.87964;
101 Ywid = 0.0514778;
102 avm1 = 3.266;
103 double HisPDFtmp[22][22] = {
104 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
105 { 0, 1e-08, 1e-08, 1e-08, 0.153529, 0.125029, 0.0999474, 0.0921067,
106 0.235419, 0.0732537, 0.244, 0.222612, 0.0362105, 1e-08, 1e-08, 1e-08,
107 1e-08, 1e-08, 1e-08, 1e-08, 1e-08, 0 },
108 { 0, 1e-08, 0.076, 0.175484, 0.0854667, 0.122373, 0.201829, 0.0794211,
109 0.14544, 0.0757183, 0.136829, 0.1594, 0.164526, 0.140187, 0.073, 0.149333,
110 1e-08, 1e-08, 1e-08, 1e-08, 1e-08, 0 },
111 { 0, 0.0752, 0.413333, 0.13014, 0.163396, 0.288067, 0.0853143, 0.104611,
112 0.145151, 0.103304, 0.216986, 0.308468, 0.307484, 0.391385, 0.260129, 0.247548,
113 0.0598261, 0.1055, 1e-08, 1e-08, 1e-08, 0 },
114 { 0, 1e-08, 0.137059, 0.132308, 0.176, 0.131097, 0.0684762, 0.124595,
115 0.291517, 0.222442, 0.0900267, 0.161487, 0.321813, 0.354053, 0.320571, 0.360848,
116 0.1766, 0.0409231, 0.171059, 1e-08, 1e-08, 0 },
117 { -0.156, 0.11619, 0.123043, 0.126873, 0.135803, 0.113026, 0.147288, 0.472,
118 0.132958, 0.260787, 0.276351, 0.394613, 0.351871, 0.279692, 0.239781, 0.202687,
119 0.132475, 0.0924118, 0.0709268, 0.224, 0.5, 0 },
120 { 0.172, 0.146514, 0.10814, 0.0553333, 0.0663133, 0.123048, 0.317266, 0.271873,
121 0.345944, 0.245123, 0.140451, 0.313574, 0.376778, 0.336216, 0.392276, 0.242783,
122 0.176113, 0.0955942, 0.224, 0.0511515, 1e-08, 0 },
123 { 3, 0.1064, 0.151429, 0.0543125, 0.171186, 0.0552727, 0.346757, 0.274806,
124 0.382867, 0.322328, 0.332562, 0.38503, 0.2928, 0.408762, 0.333543, 0.139677,
125 0.312542, 0.231683, 0.354439, 0.147167, 0.0844, 0 },
126 { -0.156, 0.278545, 0.270817, 0.392865, 0.223314, 0.230448, 0.319944, 0.196468,
127 0.138933, 0.437086, 0.544945, 0.407806, 0.508919, 0.54419, 0.62568, 0.408491,
128 0.437585, 0.402298, 0.264622, 0.277217, 0.666667, 0 },
129 { 1.454, 0.715314, 0.408417, 0.365697, 0.396971, 0.745923, 0.287053, 0.408438,
130 0.546725, 0.503155, 0.524258, 0.485722, 0.503, 0.443273, 0.591803, 0.43719,
131 0.717714, 0.517739, 0.428121, 0.23632, 0.3165, 0 },
132 { 1.266, 1.11545, 1.18407, 1.18716, 0.805808, 1.03558, 0.756932, 0.822611,
133 0.833455, 1.24262, 1.1553, 1.0626, 1.24459, 0.9324, 0.877433, 1.0099,
134 0.808188, 1.02157, 0.732585, 0.832774, 0.187636, 0 },
135 { 0, 1.855, 1.38925, 1.45507, 1.53681, 1.42868, 1.20651, 1.32688,
136 1.59969, 1.46761, 1.34576, 1.44618, 1.37882, 1.62752, 1.30176, 1.883,
137 1.93709, 2.41573, 1.74571, 1.95983, 0.6908, 0 },
138 { 0, 0.542, 0.66806, 0.684552, 0.733136, 1.05503, 0.780238, 0.581659,
139 0.87226, 1.03668, 0.894293, 1.79318, 1.16124, 1.48548, 1.37559, 1.42393,
140 1.23634, 1.22988, 1.09569, 1.69173, 1.384, 0 },
141 { 0, 0.46, 0.243667, 0.489558, 0.54961, 0.416471, 0.906316, 0.823714,
142 0.593829, 0.671209, 1.11564, 0.850765, 1.26434, 0.957657, 0.963943, 0.5821,
143 1.29176, 0.750035, 0.859412, 1.45, 3.266, 0 },
144 { 0, 1e-08, 0.44512, 0.358027, 0.254211, 0.648062, 0.382085, 0.453908,
145 0.660164, 0.6399, 0.669714, 0.780053, 0.840615, 0.611846, 1.03192, 0.755437,
146 0.538031, 0.585724, 0.371333, 1.25855, 0.422, 0 },
147 { 0, 1e-08, 0.294857, 0.122467, 0.258182, 0.211697, 0.269333, 0.487704,
148 0.37463, 0.380565, 0.698588, 0.760491, 0.831701, 0.902189, 0.690623, 0.936963,
149 0.603042, 0.9766, 0.926182, 1e-08, 1e-08, 0 },
150 { 0, 1e-08, 0.188, 0.0591111, 0.198588, 0.252526, 0.338974, 0.391897,
151 0.639178, 0.59129, 0.650435, 0.83071, 0.627353, 0.617973, 0.627358, 0.867344,
152 0.865091, 0.781273, 1.14067, 1e-08, 1e-08, 0 },
153 { 0, 1e-08, 1e-08, 0.3816, 0.09008, 0.317097, 0.303, 0.25525,
154 0.478684, 0.522815, 0.484421, 0.56464, 0.658195, 0.533506, 0.738067, 0.736197,
155 0.60128, 0.5, 1e-08, 1e-08, 1e-08, 0 },
156 { 0, 1e-08, 1e-08, 1e-08, 0.433143, 0.199186, 0.227155, 0.311676,
157 0.388952, 0.351394, 0.444396, 0.547948, 0.482377, 0.51403, 0.663027, 0.812541,
158 0.477714, 0.383, 1e-08, 1e-08, 1e-08, 0 },
159 { 0, 1e-08, 1e-08, 1e-08, 0.422, 0.127667, 0.257302, 0.387368,
160 0.392984, 0.396156, 0.588373, 0.709227, 0.633825, 0.877182, 1.02372, 0.461,
161 1, 1e-08, 1e-08, 1e-08, 1e-08, 0 },
162 { 0, 1e-08, 1e-08, 1e-08, 1e-08, 1e-08, 0.24025, 0.145867,
163 0.307442, 0.390074, 0.406828, 0.474902, 0.513756, 0.406286, 0.75, 1e-08,
164 1e-08, 1e-08, 1e-08, 1e-08, 1e-08, 0 },
165 { 0, 0, 0, 0, 0, 0, -0.039, 0.807333, 0.409778, 0.140667, 0.139273,
166 0.422, -0.039, 0.688, 0, 0, 0, 0, 0, 0, 0, 0 } };
167
168 for ( int i = 0; i < 22; i++ )
169 {
170 for ( int j = 0; j < 22; j++ ) { HisPDF[i][j] = HisPDFtmp[i][j]; }
171 }
172}
173
174void EvtD0toKpiomegaPlot::initProbMax() { noProbMax(); }
175
176void EvtD0toKpiomegaPlot::decay( EvtParticle* p ) {
177
178loop:
179 p->initializePhaseSpace( getNDaug(), getDaugs() );
180
181 EvtParticle *id1, *id2, *id3;
182 EvtVector4R pd1, pd2, pd3;
183 double xmass13, xmass12, xmass23;
184
185 id1 = p->getDaug( 0 );
186 id2 = p->getDaug( 1 );
187 id3 = p->getDaug( 2 );
188
189 pd1 = id1->getP4Lab();
190 pd2 = id2->getP4Lab();
191 pd3 = id3->getP4Lab();
192
193 xmass12 = ( pd1 + pd2 ).mass() * ( pd1 + pd2 ).mass(); // M_kpi
194 // xmass13=(pd1+pd3).mass()*(pd1+pd3).mass(); // M_kw
195 xmass23 = ( pd2 + pd3 ).mass() * ( pd2 + pd3 ).mass(); // M_piw
196
197 int xbin = FindXBin( xmass12 );
198 int ybin = FindYBin( xmass23 );
199 double xratio12 = HisPDF[xbin][ybin] / avm1;
200
201 if ( xratio12 == 0 ) goto loop;
202
203 double rd12 = EvtRandom::Flat( 0.0, 1.0 );
204
205 if ( rd12 > xratio12 ) goto loop;
206
207 return;
208}
209
210int EvtD0toKpiomegaPlot::FindXBin( double mass2 ) {
211 if ( mass2 < Xmin ) { return 0; }
212 else if ( mass2 >= Xmax ) { return 21; }
213 else { return int( ( mass2 - Xmin ) / Xwid ) + 1; }
214}
215
216int EvtD0toKpiomegaPlot::FindYBin( double mass2 ) {
217 if ( mass2 < Ymin ) { return 0; }
218 else if ( mass2 >= Ymax ) { return 21; }
219 else { return int( ( mass2 - Ymin ) / Ywid ) + 1; }
220}
mass
double mass
Definition CosmicGenerator.cxx:128
EvtD0toKpiomegaPlot::FindXBin
int FindXBin(double mass2)
Definition EvtD0toKpiomegaPlot.cc:210
EvtD0toKpiomegaPlot::getName
void getName(std::string &name)
Definition EvtD0toKpiomegaPlot.cc:38
EvtD0toKpiomegaPlot::FindYBin
int FindYBin(double mass2)
Definition EvtD0toKpiomegaPlot.cc:216
EvtD0toKpiomegaPlot::decay
void decay(EvtParticle *p)
Definition EvtD0toKpiomegaPlot.cc:176
EvtDecayBase::getParentId
EvtId getParentId()
Definition EvtDecayBase.hh:63
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtPDL::getSpinType
static EvtSpinType::spintype getSpinType(EvtId i)
Definition EvtPDL.hh:66
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
EvtParticle::getP4Lab
EvtVector4R getP4Lab()
Definition EvtParticle.cc:694
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cc:69