GooFit  v2.1.3
DP4.cu
Go to the documentation of this file.
1 #include <fstream>
2 
3 // GooFit stuff
4 #include <goofit/Application.h>
5 #include <goofit/FitManager.h>
6 #include <goofit/PDFs/basic/PolynomialPdf.h>
7 #include <goofit/PDFs/combine/AddPdf.h>
8 #include <goofit/PDFs/physics/DP4Pdf.h>
9 #include <goofit/UnbinnedDataSet.h>
10 #include <goofit/Variable.h>
11 
12 using namespace std;
13 using namespace GooFit;
14 
15 const fptype _mD0 = 1.8645;
16 const fptype piPlusMass = 0.13957018;
17 const fptype KmMass = .493677;
18 // Constants used in more than one PDF component.
19 
20 int main(int argc, char **argv) {
21  GooFit::Application app("Dalitz 4 daughter example", argc, argv);
22 
23  // Set all host normalisation just to make sure no errors
24  for(int i = 0; i < maxParams; i++)
25  host_normalisation[i] = -7;
26 
27  GOOFIT_PARSE(app, argc, argv);
28 
29  Observable m12{"m12", 0, 3};
30  Observable m34{"m34", 0, 3};
31  Observable cos12{"cos12", -1, 1};
32  Observable cos34{"m12", -1, 1};
33  Observable phi{"phi", -3.5, 3.5};
34  EventNumber eventNumber{"eventNumber", 0, INT_MAX};
35 
36  UnbinnedDataSet currData{m12, m34, cos12, cos34, phi, eventNumber};
37 
38  unsigned int MCevents = 0;
39 
40  std::string input_str = app.get_filename("ToyMC.txt", "examples/DP4");
41 
42  { // FStream block
43  fstream input(input_str, std::ios_base::in);
44 
45  while(input >> m12 >> m34 >> cos12 >> cos34 >> phi) {
46  eventNumber.setValue(MCevents++);
47  currData.addEvent();
48  }
49  }
50 
51  GOOFIT_INFO("Read in {} events", MCevents);
52 
53  DecayInfo4 DK3P_DI;
54  DK3P_DI.meson_radius = 1.5;
55  DK3P_DI.particle_masses = {_mD0, piPlusMass, piPlusMass, KmMass, piPlusMass};
56 
57  Variable RhoMass{"rho_mass", 0.77526};
58  Variable RhoWidth{"rho_width", 0.1478};
59  Variable KstarM{"KstarM", 0.89581};
60  Variable KstarW{"KstarW", 0.0474};
61  Variable f600M{"f600M", 0.519};
62  Variable f600W{"f600W", 0.454};
63  Variable a1M{"a1M", 1.23};
64  Variable a1W{"a1W", 0.42};
65  Variable K1M{"K1M", 1.272};
66  Variable K1W{"K1W", 0.09};
67  Variable K1430M{"K1430M", 1.414};
68  Variable K1430W{"K1430W", 0.29};
69 
70  // Spin factors: we have two due to the bose symmetrization of the two pi+
71  std::vector<SpinFactor *> SFKRS = {new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_S, 0, 1, 2, 3),
72  new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_S, 3, 1, 2, 0)};
73 
74  std::vector<SpinFactor *> SFKRP = {new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_P, 0, 1, 2, 3),
75  new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_P, 3, 1, 2, 0)};
76 
77  std::vector<SpinFactor *> SFKRD = {new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_D, 0, 1, 2, 3),
78  new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_D, 3, 1, 2, 0)};
79 
80  std::vector<SpinFactor *> SFKF = {new SpinFactor("SF", SF_4Body::DtoVS_VtoP1P2_StoP3P4, 2, 3, 0, 1),
81  new SpinFactor("SF", SF_4Body::DtoVS_VtoP1P2_StoP3P4, 2, 0, 3, 1)};
82 
83  std::vector<SpinFactor *> SFKK = {new SpinFactor("SF", SF_4Body::DtoAP1_AtoSP2_StoP3P4, 0, 1, 3, 2),
84  new SpinFactor("SF", SF_4Body::DtoAP1_AtoSP2_StoP3P4, 3, 1, 0, 2)};
85 
86  std::vector<SpinFactor *> SFK1R = {new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 3, 2, 0, 1),
87  new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 0, 2, 3, 1)};
88 
89  std::vector<SpinFactor *> SFA1R = {new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 2, 3, 0, 1),
90  new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 2, 0, 3, 1)};
91 
92  std::vector<SpinFactor *> SFA1RD = {new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2Dwave_VtoP3P4, 2, 3, 0, 1),
93  new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2Dwave_VtoP3P4, 2, 0, 3, 1)};
94 
95  // Lineshapes, also for both pi+ configurations
96  std::vector<Lineshape *> LSKRS = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
97  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34),
98  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24),
99  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13)};
100 
101  std::vector<Lineshape *> LSKRP = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
102  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34),
103  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24),
104  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13)};
105 
106  std::vector<Lineshape *> LSKRD = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
107  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34),
108  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24),
109  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13)};
110 
111  std::vector<Lineshape *> LSKF = {new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34),
112  new Lineshapes::Bugg("f600", f600M, f600W, 0, M_12),
113  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13),
114  new Lineshapes::Bugg("f600", f600M, f600W, 0, M_24)};
115 
116  std::vector<Lineshape *> LSKK = {new Lineshapes::SBW("K(1)(1270)bar", K1M, K1W, 1, M_34_2),
117  new Lineshapes::LASS("K(0)*(1430)bar", K1430M, K1430W, 0, M_34),
118  new Lineshapes::SBW("K(1)(1270)bar2", K1M, K1W, 1, M_13_2),
119  new Lineshapes::LASS("K(0)*(1430)bar2", K1430M, K1430W, 0, M_13)};
120 
121  std::vector<Lineshape *> LSK1R = {new Lineshapes::SBW("K(1)(1270)bar", K1M, K1W, 0, M_12_3),
122  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
123  new Lineshapes::SBW("K(1)(1270)bar", K1M, K1W, 0, M_24_3),
124  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24)};
125 
126  std::vector<Lineshape *> LSA1R = {new Lineshapes::SBW("a(1)(1260)+", a1M, a1W, 0, M_12_4),
127  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
128  new Lineshapes::SBW("a(1)(1260)+", a1M, a1W, 0, M_24_1),
129  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24)};
130 
131  std::vector<Lineshape *> LSA1RD = {new Lineshapes::SBW("a(1)(1260)+", a1M, a1W, 2, M_12_4),
132  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12),
133  new Lineshapes::SBW("a(1)(1260)+", a1M, a1W, 2, M_24_1),
134  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24)};
135 
136  // the very last parameter means that we have two permutations. so the first half of the Lineshapes
137  // and the first half of the spinfactors are amplitude 1, rest is amplitude
138  // This means that it is important for symmetrized amplitueds that the spinfactors and lineshapes are in the "right"
139  // order
140 
141  // Amplitudes with the correct fixed values.
142  // Amplitude* Bose_symmetrized_AMP_S = new Amplitude( "K*(892)rho(770)_S", new Variable("amp_real1", -0.115177), new
143  // Variable("amp_imag1", 0.153976), LSKRS, SFKRS, 2);
144  // Amplitude* Bose_symmetrized_AMP_P = new Amplitude( "K*(892)rho(770)_P", new Variable("amp_real2", -0.0298697),
145  // new Variable("amp_imag2", -0.0722874), LSKRP, SFKRP, 2);
146  // Amplitude* Bose_symmetrized_AMP_D = new Amplitude( "K*(892)rho(770)_D", new Variable("amp_real3", -0.452212), new
147  // Variable("amp_imag3", 0.426521), LSKRD, SFKRD, 2);
148 
149  // Amplitudes with floating slightly different values to be fitted.
150 
151  Amplitude Bose_symmetrized_AMP_S{"K*(892)rho(770)_S",
152  Variable("amp_real1", -0.1, 0.001, 0, 0),
153  Variable("amp_imag1", 0.1, 0.001, 0, 0),
154  LSKRS,
155  SFKRS,
156  2};
157 
158  Amplitude Bose_symmetrized_AMP_P{"K*(892)rho(770)_P",
159  Variable("amp_real2", -0.02, 0.001, 0, 0),
160  Variable("amp_imag2", -0.07, 0.001, 0, 0),
161  LSKRP,
162  SFKRP,
163  2};
164  Amplitude Bose_symmetrized_AMP_D{"K*(892)rho(770)_D",
165  Variable("amp_real3", -0.4, 0.001, 0, 0),
166  Variable("amp_imag3", 0.4, 0.001, 0, 0),
167  LSKRD,
168  SFKRD,
169  2};
170 
171  Amplitude Bose_symmetrized_KF{
172  "KF", Variable("amp_real4", 0.0120787), Variable("amp_imag4", -0.0332525), LSKF, SFKF, 2};
173  Amplitude Bose_symmetrized_KK{
174  "LSKK", Variable("amp_real5", 0.0109033), Variable("amp_imag5", -0.00186219), LSKK, SFKK, 2};
175  Amplitude Bose_symmetrized_K1R{
176  "LSK1R", Variable("amp_real6", -0.10728), Variable("amp_imag6", -0.130213), LSK1R, SFK1R, 2};
177  Amplitude Bose_symmetrized_A1R{"LSA1R", Variable("amp_real7", 1.0), Variable("amp_imag7", 0.0), LSA1R, SFA1R, 2};
178  Amplitude Bose_symmetrized_A1RD{
179  "LSA1RD", Variable("amp_real8", -0.94921), Variable("amp_imag8", -1.73407), LSA1RD, SFA1RD, 2};
180 
181  DK3P_DI.amplitudes = {&Bose_symmetrized_KF,
182  &Bose_symmetrized_AMP_S,
183  &Bose_symmetrized_AMP_P,
184  &Bose_symmetrized_AMP_D,
185  &Bose_symmetrized_KK,
186  &Bose_symmetrized_K1R,
187  &Bose_symmetrized_A1R,
188  &Bose_symmetrized_A1RD};
189 
190  for(auto res : LSKF)
191  res->setParameterConstantness(true);
192 
193  for(auto res : LSKRS)
194  res->setParameterConstantness(true);
195 
196  for(auto res : LSKRP)
197  res->setParameterConstantness(true);
198 
199  for(auto res : LSKRD)
200  res->setParameterConstantness(true);
201 
202  for(auto res : LSKK)
203  res->setParameterConstantness(true);
204 
205  for(auto res : LSK1R)
206  res->setParameterConstantness(true);
207 
208  for(auto res : LSA1R)
209  res->setParameterConstantness(true);
210 
211  for(auto res : LSA1RD)
212  res->setParameterConstantness(true);
213 
214  Variable constantOne{"constantOne", 1};
215  Variable constantZero{"constantZero", 0};
216 
217  vector<Observable> observables = {m12, m34, cos12, cos34, phi, eventNumber};
218  vector<Variable> coefficients = {constantOne};
219  vector<Variable> offsets = {constantZero, constantZero};
220 
221  PolynomialPdf eff{"constantEff", observables, coefficients, offsets, 0};
222  DPPdf dp{"test", observables, DK3P_DI, &eff, 1000000};
223 
224  Variable constant{"constant", 0.1};
225  Variable constant2{"constant2", 1.0};
226  PolynomialPdf backgr{"backgr", m12, {constant}};
227  AddPdf signal{"signal", constant2, &dp, &backgr};
228 
229  signal.setData(&currData);
230 
231  dp.setDataSize(currData.getNumEvents(), 6);
232 
233  FitManager datapdf{&signal};
234  datapdf.fit();
235 
236  return datapdf;
237 }
Generated on Sat Apr 21 2018 20:13:17 for GooFit by   doxygen 1.8.13