GooFit  v2.1.3
TDDP4.cu
Go to the documentation of this file.
1 // ROOT
2 #include <TFile.h>
3 #include <TTree.h>
4 
5 // GooFit stuff
6 #include <fstream>
7 #include <goofit/Application.h>
8 #include <goofit/PDFs/basic/PolynomialPdf.h>
9 #include <goofit/PDFs/combine/AddPdf.h>
10 #include <goofit/PDFs/physics/DP4Pdf.h>
11 #include <goofit/PDFs/physics/Tddp4Pdf.h>
12 #include <goofit/PDFs/physics/TruthResolution_Aux.h>
13 #include <goofit/UnbinnedDataSet.h>
14 #include <goofit/Variable.h>
15 #include <thrust/count.h>
16 
17 using namespace std;
18 using namespace GooFit;
19 
20 // Constants used in more than one PDF component.
21 const fptype _mD0 = 1.8645;
22 const fptype piPlusMass = 0.13957018;
23 const fptype KmMass = .493677;
24 
25 int main(int argc, char **argv) {
26  GooFit::Application app("Time dependent Dalitz plot, 4 particles", argc, argv);
27 
28  TString output = "test_10_15.output";
29  app.add_option("-o,--output,output", output, "File to output", true)->check(GooFit::NonexistentPath);
30 
31  int trials = 100;
32  app.add_option("-t,--trials,output", trials, "Number of trials", true);
33 
34  GOOFIT_PARSE(app);
35 
36  DecayInfo4t DK3P_DI{Variable("tau", 0.4101, 0.001, 0.300, 0.500),
37  Variable("xmixing", 0.005, 0.001, 0, 0),
38  Variable("ymixing", 0.01, 0.001, 0, 0),
39  Variable("SqWStoRSrate", 1.0 / sqrt(300.0))};
40 
41  DK3P_DI.meson_radius = 1.5;
42  DK3P_DI.particle_masses.push_back(_mD0);
43  DK3P_DI.particle_masses.push_back(piPlusMass);
44  DK3P_DI.particle_masses.push_back(piPlusMass);
45  DK3P_DI.particle_masses.push_back(KmMass);
46  DK3P_DI.particle_masses.push_back(piPlusMass);
47 
48  Variable RhoMass{"rho_mass", 0.77526, 0.01, 0.7, 0.8};
49  Variable RhoWidth{"rho_width", 0.1478, 0.01, 0.1, 0.2};
50  Variable KstarM{"KstarM", 0.89581, 0.01, 0.9, 0.1};
51  Variable KstarW{"KstarW", 0.0474, 0.01, 0.1, 0.2};
52 
53  // Variable* f600M = new Variable("f600M", 0.519, 0.01, 0.75, 0.85);
54  // Variable* f600W = new Variable("f600W", 0.454, 0.01, 0.75, 0.85);
55  // Variable* a1M = new Variable("a1M", 1.23, 0.01, 1.2, 1.3);
56  // Variable* a1W = new Variable("a1W", 0.42, 0.01, 0.37, 0.47);
57  // Variable* K1M = new Variable("K1M", 1.272, 0.01, 1.2, 1.3);
58  // Variable* K1W = new Variable("K1W", 0.09, 0.01, 0.08, 0.1);
59  // Variable* K1430M = new Variable("K1430M", 1.414, 0.01, 1.4, 1.5);
60  // Variable* K1430W = new Variable("K1430W", .29, 0.01, 0.25, 0.35);
61 
62  // Spin factors: we have two due to the bose symmetrization of the two pi+
63  std::vector<SpinFactor *> SFKRS = {new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_S, 0, 1, 2, 3),
64  new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_S, 3, 1, 2, 0)};
65 
66  std::vector<SpinFactor *> SFKRP;
67  SFKRP.push_back(new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_P, 0, 1, 2, 3));
68  SFKRP.push_back(new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_P, 3, 1, 2, 0));
69 
70  std::vector<SpinFactor *> SFKRD;
71  SFKRD.push_back(new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_D, 0, 1, 2, 3));
72  SFKRD.push_back(new SpinFactor("SF", SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_D, 3, 1, 2, 0));
73 
74  std::vector<SpinFactor *> SFKF;
75  SFKF.push_back(new SpinFactor("SF", SF_4Body::DtoVS_VtoP1P2_StoP3P4, 2, 3, 0, 1));
76  SFKF.push_back(new SpinFactor("SF", SF_4Body::DtoVS_VtoP1P2_StoP3P4, 2, 0, 3, 1));
77 
78  std::vector<SpinFactor *> SFKK;
79  SFKK.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoSP2_StoP3P4, 0, 1, 3, 2));
80  SFKK.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoSP2_StoP3P4, 3, 1, 0, 2));
81 
82  std::vector<SpinFactor *> SFK1R;
83  SFK1R.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 3, 2, 0, 1));
84  SFK1R.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 0, 2, 3, 1));
85 
86  std::vector<SpinFactor *> SFA1R;
87  SFA1R.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 2, 3, 0, 1));
88  SFA1R.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2_VtoP3P4, 2, 0, 3, 1));
89 
90  std::vector<SpinFactor *> SFA1RD;
91  SFA1RD.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2Dwave_VtoP3P4, 2, 3, 0, 1));
92  SFA1RD.push_back(new SpinFactor("SF", SF_4Body::DtoAP1_AtoVP2Dwave_VtoP3P4, 2, 0, 3, 1));
93 
94  // Lineshapes, also for both pi+ configurations
95  std::vector<Lineshape *> LSKRS = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12, FF::BL2),
96  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34, FF::BL2),
97  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24, FF::BL2),
98  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13, FF::BL2)};
99 
100  std::vector<Lineshape *> LSKRP = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12, FF::BL2),
101  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34, FF::BL2),
102  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24, FF::BL2),
103  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13, FF::BL2)};
104 
105  std::vector<Lineshape *> LSKRD = {new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_12, FF::BL2),
106  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_34, FF::BL2),
107  new Lineshapes::RBW("rho(770)", RhoMass, RhoWidth, 1, M_24, FF::BL2),
108  new Lineshapes::RBW("K*(892)bar", KstarM, KstarW, 1, M_13, FF::BL2)};
109 
110  // the very last parameter means that we have two permutations. so the first half of the Lineshapes
111  // and the first half of the spinfactors are amplitude 1, rest is amplitude 2
112  // This means that it is important for symmetrized amplitudes that the spinfactors and lineshapes are in the "right"
113  // order
114 
115  Amplitude Bose_symmetrized_AMP_S{
116  "K*(892)rho(770)_S", Variable("amp_real1", 1.0), Variable("amp_imag1", 0.0), LSKRS, SFKRS, 2};
117  Amplitude Bose_symmetrized_AMP_P{
118  "K*(892)rho(770)_P", Variable("amp_real2", 0.526), Variable("amp_imag2", -0.626), LSKRP, SFKRP, 2};
119  Amplitude Bose_symmetrized_AMP_D{
120  "K*(892)rho(770)_D", Variable("amp_real3", 26.537), Variable("amp_imag3", 12.284), LSKRD, SFKRD, 2};
121 
122  Amplitude Bose_symmetrized_AMP_S_B{
123  "B_K*(892)rho(770)_S", Variable("amp_real1", 1.0), Variable("amp_imag1", 0), LSKRS, SFKRS, 2};
124  Amplitude Bose_symmetrized_AMP_P_B{
125  "B_K*(892)rho(770)_P", Variable("amp_real2", -0.145), Variable("amp_imag2", 0.86), LSKRP, SFKRP, 2};
126  Amplitude Bose_symmetrized_AMP_D_B{
127  "B_K*(892)rho(770)_D", Variable("amp_real3", 24.343), Variable("amp_imag3", 5.329), LSKRD, SFKRD, 2};
128 
129  DK3P_DI.amplitudes_B.push_back(&Bose_symmetrized_AMP_S);
130  DK3P_DI.amplitudes_B.push_back(&Bose_symmetrized_AMP_P);
131  DK3P_DI.amplitudes_B.push_back(&Bose_symmetrized_AMP_D);
132 
133  DK3P_DI.amplitudes.push_back(&Bose_symmetrized_AMP_S_B);
134  DK3P_DI.amplitudes.push_back(&Bose_symmetrized_AMP_P_B);
135  DK3P_DI.amplitudes.push_back(&Bose_symmetrized_AMP_D_B);
136 
137  Observable m12{"m12", 0, 3};
138  Observable m34{"m34", 0, 3};
139  Observable cos12{"cos12", -1, 1};
140  Observable cos34{"m12", -1, 1};
141  Observable phi{"phi", -3.5, 3.5};
142  EventNumber eventNumber{"eventNumber"};
143  Observable dtime{"dtime", 0, 10};
144  Observable sigmat{"sigmat", -3, 3};
145  Variable constantOne{"constantOne", 1};
146  Variable constantZero{"constantZero", 0};
147 
148  vector<Observable> observables{m12, m34, cos12, cos34, phi, eventNumber, dtime, sigmat};
149  vector<Variable> offsets{constantZero, constantZero};
150  vector<Variable> coefficients{constantOne};
151 
152  TruthResolution dat;
153  PolynomialPdf eff{"constantEff", observables, coefficients, offsets, 0};
154  TDDP4 dp{"test", observables, DK3P_DI, &dat, &eff, 0, 1};
155 
156  TFile *file = new TFile(output, "RECREATE");
157  TTree *tree = new TTree("events", "events");
158 
159  double tm12, tm34, tc12, tc34, tphi, tdtime, D0_E, D0_Px, D0_Py, D0_Pz, Kplus_E, Kplus_Px, Kplus_Py, Kplus_Pz,
160  Piminus1_E, Piminus1_Px, Piminus1_Py, Piminus1_Pz, Piminus2_E, Piminus2_Px, Piminus2_Py, Piminus2_Pz, Piplus_E,
161  Piplus_Px, Piplus_Py, Piplus_Pz;
162  int D0_pdg, Kplus_pdg, Piminus1_pdg, Piminus2_pdg, Piplus_pdg;
163 
164  tree->Branch("m12", &tm12, "m12/D");
165  tree->Branch("m34", &tm34, "m34/D");
166  tree->Branch("c12", &tc12, "c12/D");
167  tree->Branch("c34", &tc34, "c34/D");
168  tree->Branch("phi", &tphi, "phi/D");
169  tree->Branch("dtime", &tdtime, "dtime/D");
170  tree->Branch("D0_E", &D0_E, "D0_E/D");
171  tree->Branch("D0_Px", &D0_Px, "D0_Px/D");
172  tree->Branch("D0_Py", &D0_Py, "D0_Py/D");
173  tree->Branch("D0_Pz", &D0_Pz, "D0_Pz/D");
174  tree->Branch("D0_pdg", &D0_pdg, "D0_pdg/I");
175  tree->Branch("Kplus_E", &Kplus_E, "Kplus_E/D");
176  tree->Branch("Kplus_Px", &Kplus_Px, "Kplus_Px/D");
177  tree->Branch("Kplus_Py", &Kplus_Py, "Kplus_Py/D");
178  tree->Branch("Kplus_Pz", &Kplus_Pz, "Kplus_Pz/D");
179  tree->Branch("Kplus_pdg", &Kplus_pdg, "Kplus_pdg/I");
180  tree->Branch("Piminus1_E", &Piminus1_E, "Piminus1_E/D");
181  tree->Branch("Piminus1_Px", &Piminus1_Px, "Piminus1_Px/D");
182  tree->Branch("Piminus1_Py", &Piminus1_Py, "Piminus1_Py/D");
183  tree->Branch("Piminus1_Pz", &Piminus1_Pz, "Piminus1_Pz/D");
184  tree->Branch("Piminus1_pdg", &Piminus1_pdg, "Piminus1_pdg/I");
185  tree->Branch("Piminus2_E", &Piminus2_E, "Piminus2_E/D");
186  tree->Branch("Piminus2_Px", &Piminus2_Px, "Piminus2_Px/D");
187  tree->Branch("Piminus2_Py", &Piminus2_Py, "Piminus2_Py/D");
188  tree->Branch("Piminus2_Pz", &Piminus2_Pz, "Piminus2_Pz/D");
189  tree->Branch("Piminus2_pdg", &Piminus2_pdg, "Piminus2_pdg/I");
190  tree->Branch("Piplus_E", &Piplus_E, "Piplus_E/D");
191  tree->Branch("Piplus_Px", &Piplus_Px, "Piplus_Px/D");
192  tree->Branch("Piplus_Py", &Piplus_Py, "Piplus_Py/D");
193  tree->Branch("Piplus_Pz", &Piplus_Pz, "Piplus_Pz/D");
194  tree->Branch("Piplus_pdg", &Piplus_pdg, "Piplus_pdg/I");
195 
196  int total_accepted = 0;
197 
198  for(int k = 0; k < trials; ++k) {
199  int numEvents = 800000;
200  dp.setGenerationOffset(k * numEvents);
201 
202  mcbooster::ParticlesSet_h particles; // typedef for std::vector<Particles_h *>
203  mcbooster::VariableSet_h variables;
204  mcbooster::RealVector_h weights;
205  mcbooster::BoolVector_h flags;
206 
207  std::tie(particles, variables, weights, flags) = dp.GenerateSig(numEvents);
208 
209  int accepted = thrust::count_if(flags.begin(), flags.end(), thrust::identity<bool>());
210  total_accepted += accepted;
211 
212  GOOFIT_INFO(
213  "Run #{}: Using accept-reject method would leave you with {} out of {} events", k, accepted, numEvents);
214 
215  for(int i = 0; i < weights.size(); ++i) {
216  if(flags[i] == 1) {
217  // printf("%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n", (*(variables[0]))[i], (*(variables[1]))[i],
218  // (*(variables[2]))[i], (*(variables[3]))[i], (*(variables[4]))[i], weights[i], flags[i]);
219  tm12 = (*(variables[0]))[i];
220  tm34 = (*(variables[1]))[i];
221  tc12 = (*(variables[2]))[i];
222  tc34 = (*(variables[3]))[i];
223  tphi = (*(variables[4]))[i];
224  tdtime = (*(variables[5]))[i];
225  D0_E = 1864;
226  D0_Px = 0.0;
227  D0_Py = 0.0;
228  D0_Pz = 0.0;
229  D0_pdg = 421;
230  Kplus_E = 1000 * (*(particles[2]))[i].get(0);
231  Kplus_Px = 1000 * (*(particles[2]))[i].get(1);
232  Kplus_Py = 1000 * (*(particles[2]))[i].get(2);
233  Kplus_Pz = 1000 * (*(particles[2]))[i].get(3);
234  Kplus_pdg = -321;
235  Piminus1_E = 1000 * (*(particles[3]))[i].get(0);
236  Piminus1_Px = 1000 * (*(particles[3]))[i].get(1);
237  Piminus1_Py = 1000 * (*(particles[3]))[i].get(2);
238  Piminus1_Pz = 1000 * (*(particles[3]))[i].get(3);
239  Piminus1_pdg = 211;
240  Piminus2_E = 1000 * (*(particles[0]))[i].get(0);
241  Piminus2_Px = 1000 * (*(particles[0]))[i].get(1);
242  Piminus2_Py = 1000 * (*(particles[0]))[i].get(2);
243  Piminus2_Pz = 1000 * (*(particles[0]))[i].get(3);
244  Piminus2_pdg = 211;
245  Piplus_E = 1000 * (*(particles[1]))[i].get(0);
246  Piplus_Px = 1000 * (*(particles[1]))[i].get(1);
247  Piplus_Py = 1000 * (*(particles[1]))[i].get(2);
248  Piplus_Pz = 1000 * (*(particles[1]))[i].get(3);
249  Piplus_pdg = -211;
250 
251  tree->Fill();
252  }
253  }
254 
255  delete variables[0];
256  delete variables[1];
257  delete variables[2];
258  delete variables[3];
259  delete variables[4];
260  delete variables[5];
261 
262  delete particles[0];
263  delete particles[1];
264  delete particles[2];
265  delete particles[3];
266  }
267 
268  tree->Write();
269  file->Close();
270 
271  if(total_accepted > 0)
272  return 0;
273  else {
274  GOOFIT_ERROR("Total accepted was 0! Something is wrong.");
275  return 1;
276  }
277 }
Generated on Sat Apr 21 2018 20:13:17 for GooFit by   doxygen 1.8.13