GooFit  v2.1.3
dalitz.cpp
Go to the documentation of this file.
1 // ROOT stuff
2 #include <TCanvas.h>
3 #include <TFile.h>
4 #include <TH1F.h>
5 #include <TH2F.h>
6 #include <TLegend.h>
7 #include <TLine.h>
8 #include <TRandom.h>
9 #include <TRandom3.h>
10 
11 // System stuff
12 #include <fstream>
13 #include <sys/time.h>
14 #include <sys/times.h>
15 
16 // GooFit stuff
17 #include <goofit/Application.h>
18 #include <goofit/FitManager.h>
19 #include <goofit/PDFs/GooPdf.h>
27 #include <goofit/UnbinnedDataSet.h>
28 #include <goofit/Variable.h>
29 #include <goofit/utilities/Style.h>
30 
31 using namespace std;
32 using namespace GooFit;
33 
34 Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8);
35 Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2);
36 
37 const fptype _mD0 = 1.86484;
38 const fptype _mD02 = _mD0 * _mD0;
39 const fptype _mD02inv = 1. / _mD02;
40 const fptype piPlusMass = 0.13957018;
41 const fptype piZeroMass = 0.1349766;
42 
43 // Constants used in more than one PDF component.
44 Variable motherM("motherM", _mD0);
45 Variable chargeM("chargeM", piPlusMass);
46 Variable neutrlM("neutrlM", piZeroMass);
47 Variable massSum("massSum", _mD0 *_mD0 + 2 * piPlusMass * piPlusMass + piZeroMass * piZeroMass); // = 3.53481
48 Variable constantOne("constantOne", 1);
49 Variable constantZero("constantZero", 0);
50 
51 fptype cpuGetM23(fptype massPZ, fptype massPM) {
52  return (_mD02 + piZeroMass * piZeroMass + piPlusMass * piPlusMass + piPlusMass * piPlusMass - massPZ - massPM);
53 }
54 
56  toyFileName = app.get_filename(toyFileName, "examples/dalitz");
57 
58  auto obs = data.getObservables();
59  Observable m12 = obs.at(0);
60  Observable m13 = obs.at(1);
61  Observable eventNumber = obs.at(2);
62 
63  TH2F dalitzplot("dalitzplot",
64  "Original Data",
65  m12.getNumBins(),
66  m12.getLowerLimit(),
67  m12.getUpperLimit(),
68  m13.getNumBins(),
69  m13.getLowerLimit(),
70  m13.getUpperLimit());
71  std::vector<Observable> vars;
72  vars.push_back(m12);
73  vars.push_back(m13);
74  vars.push_back(eventNumber);
75 
76  std::ifstream reader(toyFileName);
77  std::string buffer;
78 
79  while(reader >> buffer) {
80  if(buffer == "====")
81  break;
82  std::cout << buffer;
83  }
84 
85  double dummy = 0;
86 
87  while(reader >> dummy
88  // m23, m(pi+ pi-), called m12 in processToyRoot convention.
89  // Already swapped according to D* charge. m12 = m(pi+pi0)
90  >> dummy >> m12 >> m13
91 
92  // Errors on Dalitz variables
93  >> dummy >> dummy >> dummy
94 
95  // Decay time, sigma_t
96  >> dummy >> dummy
97 
98  // Md0, deltaM, Prob, Sig, Dst charge, Run, Event, Signal and four bkg fractions
99  >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy) {
100  // EXERCISE 1 (preliminary): Impose an artificial reconstruction efficiency
101  // by throwing out events with a probability linear in m12.
102  // NB! This exercise continues below.
103 
104  // EXERCISE 2: Instead of the above efficiency, impose a
105  // K0 veto, by throwing out events with 0.475 < m23 < 0.505.
106 
107  // EXERCISE 3: Use both the above.
108 
109  eventNumber.setValue(data.getNumEvents());
110  data.addEvent();
111 
112  dalitzplot.Fill(m12.getValue(), m13.getValue());
113  }
114 
115  GOOFIT_INFO("Read in {} events", data.getNumEvents());
116 
117  TCanvas foo;
118  dalitzplot.SetStats(false);
119  dalitzplot.Draw("colz");
120  foo.SaveAs("dalitzplot.png");
121 }
122 
124 
127  dtop0pp.motherMass = _mD0;
128  dtop0pp.daug1Mass = piZeroMass;
129  dtop0pp.daug2Mass = piPlusMass;
130  dtop0pp.daug3Mass = piPlusMass;
131  dtop0pp.meson_radius = 1.5;
132 
133  ResonancePdf *rhop = new Resonances::RBW(
134  "rhop", Variable("rhop_amp_real", 1), Variable("rhop_amp_imag", 0), fixedRhoMass, fixedRhoWidth, 1, PAIR_12);
135 
136  bool fixAmps = false; // Takes ~400x longer
137 
138  ResonancePdf *rhom = new Resonances::RBW(
139  "rhom",
140  fixAmps ? Variable("rhom_amp_real", 0.714) : Variable("rhom_amp_real", 0.714, 0.001, 0, 0),
141  fixAmps ? Variable("rhom_amp_imag", -0.025) : Variable("rhom_amp_imag", -0.025, 0.1, 0, 0),
142  fixedRhoMass,
144  1,
145  PAIR_13);
146 
147  ResonancePdf *rho0 = new Resonances::RBW(
148  "rho0",
149  fixAmps ? Variable("rho0_amp_real", 0.565) : Variable("rho0_amp_real", 0.565, 0.001, 0, 0),
150  fixAmps ? Variable("rho0_amp_imag", 0.164) : Variable("rho0_amp_imag", 0.164, 0.1, 0, 0),
151  fixedRhoMass,
153  1,
154  PAIR_23);
155 
156  Variable sharedMass("rhop_1450_mass", 1.465, 0.01, 1.0, 2.0);
157  Variable shareWidth("rhop_1450_width", 0.400, 0.01, 0.01, 5.0);
158 
159  ResonancePdf *rhop_1450 = new Resonances::RBW(
160  "rhop_1450",
161  fixAmps ? Variable("rhop_1450_amp_real", -0.174) : Variable("rhop_1450_amp_real", -0.174, 0.001, 0, 0),
162  fixAmps ? Variable("rhop_1450_amp_imag", -0.117) : Variable("rhop_1450_amp_imag", -0.117, 0.1, 0, 0),
163  sharedMass,
164  shareWidth,
165  1,
166  PAIR_12);
167 
168  ResonancePdf *rho0_1450 = new Resonances::RBW(
169  "rho0_1450",
170  fixAmps ? Variable("rho0_1450_amp_real", 0.325) : Variable("rho0_1450_amp_real", 0.325, 0.001, 0, 0),
171  fixAmps ? Variable("rho0_1450_amp_imag", 0.057) : Variable("rho0_1450_amp_imag", 0.057, 0.1, 0, 0),
172  sharedMass,
173  shareWidth,
174  1,
175  PAIR_23);
176 
177  ResonancePdf *rhom_1450 = new Resonances::RBW(
178  "rhom_1450",
179  fixAmps ? Variable("rhom_1450_amp_real", 0.788) : Variable("rhom_1450_amp_real", 0.788, 0.001, 0, 0),
180  fixAmps ? Variable("rhom_1450_amp_imag", 0.226) : Variable("rhom_1450_amp_imag", 0.226, 0.1, 0, 0),
181  sharedMass,
182  shareWidth,
183  1,
184  PAIR_13);
185 
186  Variable sharedMass2("rhop_1700_mass", 1.720, 0.01, 1.6, 1.9);
187  Variable shareWidth2("rhop_1700_width", 0.250, 0.01, 0.1, 1.0);
188 
189  ResonancePdf *rhop_1700 = new Resonances::RBW(
190  "rhop_1700",
191  fixAmps ? Variable("rhop_1700_amp_real", 2.151) : Variable("rhop_1700_amp_real", 2.151, 0.001, 0, 0),
192  fixAmps ? Variable("rhop_1700_amp_imag", -0.658) : Variable("rhop_1700_amp_imag", -0.658, 0.1, 0, 0),
193  sharedMass2,
194  shareWidth2,
195  1,
196  PAIR_12);
197 
198  ResonancePdf *rho0_1700 = new Resonances::RBW(
199  "rho0_1700",
200  fixAmps ? Variable("rho0_1700_amp_real", 2.400) : Variable("rho0_1700_amp_real", 2.400, 0.001, 0, 0),
201  fixAmps ? Variable("rho0_1700_amp_imag", -0.734) : Variable("rho0_1700_amp_imag", -0.734, 0.1, 0, 0),
202  sharedMass2,
203  shareWidth2,
204  1,
205  PAIR_23);
206 
207  ResonancePdf *rhom_1700 = new Resonances::RBW(
208  "rhom_1700",
209  fixAmps ? Variable("rhom_1700_amp_real", 1.286) : Variable("rhom_1700_amp_real", 1.286, 0.001, 0, 0),
210  fixAmps ? Variable("rhom_1700_amp_imag", -1.532) : Variable("rhom_1700_amp_imag", -1.532, 0.1, 0, 0),
211  sharedMass2,
212  shareWidth2,
213  1,
214  PAIR_13);
215 
216  ResonancePdf *f0_980 = new Resonances::RBW("f0_980",
217  fixAmps ? Variable("f0_980_amp_real", 0.008 * (-_mD02))
218  : Variable("f0_980_amp_real", 0.008 * (-_mD02), 0.001, 0, 0),
219  fixAmps ? Variable("f0_980_amp_imag", -0.013 * (-_mD02))
220  : Variable("f0_980_amp_imag", -0.013 * (-_mD02), 0.1, 0, 0),
221  Variable("f0_980_mass", 0.980, 0.01, 0.8, 1.2),
222  Variable("f0_980_width", 0.044, 0.001, 0.001, 0.08),
223  0,
224  PAIR_23);
225 
226  ResonancePdf *f0_1370 = new Resonances::RBW("f0_1370",
227  fixAmps ? Variable("f0_1370_amp_real", -0.058 * (-_mD02))
228  : Variable("f0_1370_amp_real", -0.058 * (-_mD02), 0.001, 0, 0),
229  fixAmps ? Variable("f0_1370_amp_imag", 0.026 * (-_mD02))
230  : Variable("f0_1370_amp_imag", 0.026 * (-_mD02), 0.1, 0, 0),
231  Variable("f0_1370_mass", 1.434, 0.01, 1.2, 1.6),
232  Variable("f0_1370_width", 0.173, 0.01, 0.01, 0.4),
233  0,
234  PAIR_23);
235 
236  ResonancePdf *f0_1500 = new Resonances::RBW("f0_1500",
237  fixAmps ? Variable("f0_1500_amp_real", 0.057 * (-_mD02))
238  : Variable("f0_1500_amp_real", 0.057 * (-_mD02), 0.001, 0, 0),
239  fixAmps ? Variable("f0_1500_amp_imag", 0.012 * (-_mD02))
240  : Variable("f0_1500_amp_imag", 0.012 * (-_mD02), 0.1, 0, 0),
241  Variable("f0_1500_mass", 1.507, 0.01, 1.3, 1.7),
242  Variable("f0_1500_width", 0.109, 0.01, 0.01, 0.3),
243  0,
244  PAIR_23);
245 
246  ResonancePdf *f0_1710 = new Resonances::RBW("f0_1710",
247  fixAmps ? Variable("f0_1710_amp_real", 0.070 * (-_mD02))
248  : Variable("f0_1710_amp_real", 0.070 * (-_mD02), 0.001, 0, 0),
249  fixAmps ? Variable("f0_1710_amp_imag", 0.087 * (-_mD02))
250  : Variable("f0_1710_amp_imag", 0.087 * (-_mD02), 0.1, 0, 0),
251  Variable("f0_1710_mass", 1.714, 0.01, 1.5, 2.9),
252  Variable("f0_1710_width", 0.140, 0.01, 0.01, 0.5),
253  0,
254  PAIR_23);
255 
256  ResonancePdf *f2_1270
257  = new Resonances::RBW("f2_1270",
258  fixAmps ? Variable("f2_1270_amp_real", -1.027 * (-_mD02inv))
259  : Variable("f2_1270_amp_real", -1.027 * (-_mD02inv), 0.001, 0, 0),
260  fixAmps ? Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv))
261  : Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv), 0.1, 0, 0),
262  Variable("f2_1270_mass", 1.2754, 0.01, 1.0, 1.5),
263  Variable("f2_1270_width", 0.1851, 0.01, 0.01, 0.4),
264  2,
265  PAIR_23);
266 
267  ResonancePdf *f0_600 = new Resonances::RBW("f0_600",
268  fixAmps ? Variable("f0_600_amp_real", 0.068 * (-_mD02))
269  : Variable("f0_600_amp_real", 0.068 * (-_mD02), 0.001, 0, 0),
270  fixAmps ? Variable("f0_600_amp_imag", 0.010 * (-_mD02))
271  : Variable("f0_600_amp_imag", 0.010 * (-_mD02), 0.1, 0, 0),
272  Variable("f0_600_mass", 0.500, 0.01, 0.3, 0.7),
273  Variable("f0_600_width", 0.400, 0.01, 0.2, 0.6),
274  0,
275  PAIR_23);
276 
277  ResonancePdf *nonr = new Resonances::NonRes(
278  "nonr",
279  fixAmps ? Variable("nonr_amp_real", 0.5595 * (-1)) : Variable("nonr_amp_real", 0.5595 * (-1), 0.001, 0, 0),
280  fixAmps ? Variable("nonr_amp_imag", -0.108761 * (-1)) : Variable("nonr_amp_imag", -0.108761 * (-1), 0.1, 0, 0));
281 
282  dtop0pp.resonances.push_back(nonr);
283  dtop0pp.resonances.push_back(rhop);
284  dtop0pp.resonances.push_back(rho0);
285  dtop0pp.resonances.push_back(rhom);
286  dtop0pp.resonances.push_back(rhop_1450);
287  dtop0pp.resonances.push_back(rho0_1450);
288  dtop0pp.resonances.push_back(rhom_1450);
289  dtop0pp.resonances.push_back(rhop_1700);
290  dtop0pp.resonances.push_back(rho0_1700);
291  dtop0pp.resonances.push_back(rhom_1700);
292  dtop0pp.resonances.push_back(f0_980);
293  dtop0pp.resonances.push_back(f0_1370);
294  dtop0pp.resonances.push_back(f0_1500);
295  dtop0pp.resonances.push_back(f0_1710);
296  dtop0pp.resonances.push_back(f2_1270);
297  dtop0pp.resonances.push_back(f0_600);
298 
299  bool fitMasses = false;
300 
301  if(!fitMasses) {
302  for(vector<ResonancePdf *>::iterator res = dtop0pp.resonances.begin(); res != dtop0pp.resonances.end(); ++res) {
303  (*res)->setParameterConstantness(true);
304  }
305  }
306 
307  if(!eff) {
308  // By default create a constant efficiency.
309  vector<Variable> offsets = {constantZero, constantZero};
310  vector<Observable> observables = {m12, m13};
311  vector<Variable> coefficients = {constantOne};
312 
313  eff = new PolynomialPdf("constantEff", observables, coefficients, offsets, 0);
314  }
315 
316  return new DalitzPlotPdf("signalPDF", m12, m13, eventNumber, dtop0pp, eff);
317 }
318 
320  // EXERCISE 1 (real part): Create a PolynomialPdf which models
321  // the efficiency you imposed in the preliminary, and use it in constructing
322  // the signal PDF.
323 
324  // EXERCISE 2: Create a K0 veto function and use it as the efficiency.
325 
326  // EXERCISE 3: Make the efficiency a product of the two functions
327  // from the previous exercises.
328 
329  signal->setData(data);
330  signal->setDataSize(data->getNumEvents());
331  FitManager datapdf(signal);
332 
333  datapdf.fit();
334 
335  ProdPdf prodpdf{"prodpdf", {signal}};
336 
337  DalitzPlotter plotter(&prodpdf, signal);
338 
339  TCanvas foo;
340  TH2F *dalitzplot = plotter.make2D();
341  dalitzplot->Draw("colz");
342 
343  foo.SaveAs("dalitzpdf.png");
344 
345  return datapdf;
346 }
347 
348 int main(int argc, char **argv) {
349  GooFit::Application app("Dalitz example", argc, argv);
350 
351  std::string filename = "dalitz_toyMC_000.txt";
352  app.add_option("-f,--filename,filename", filename, "File to read in", true)->check(GooFit::ExistingFile);
353 
354  bool make_toy;
355  app.add_flag("-m,--make-toy", make_toy, "Make a toy instead of reading a file in");
356 
357  GOOFIT_PARSE(app);
358 
359  GooFit::setROOTStyle();
360 
361  // Observables setup
362  Observable m12("m12", 0, 3);
363  Observable m13("m13", 0, 3);
364  EventNumber eventNumber("eventNumber");
365  m12.setNumBins(240);
366  m13.setNumBins(240);
367 
368  // Prepare the data
369  UnbinnedDataSet data({m12, m13, eventNumber});
370 
371  // Set up the model
372  DalitzPlotPdf *signal = makeSignalPdf(m12, m13, eventNumber);
373 
374  // A wrapper for plotting without complex number segfault
375  ProdPdf prodpdf{"prodpdf", {signal}};
376 
377  // Add nice tool for making data or plotting
378  DalitzPlotter dplotter{&prodpdf, signal};
379 
380  // Read in data
381  if(make_toy) {
382  dplotter.fillDataSetMC(data, 1000000);
383  } else {
384  getToyData(filename, app, data);
385  }
386 
387  try {
388  return runToyFit(signal, &data);
389  } catch(const std::runtime_error &e) {
390  std::cerr << e.what() << std::endl;
391  return 7;
392  }
393 }
size_t getNumEvents() const
Definition: DataSet.h:47
double fptype
const fptype piZeroMass
Definition: dalitz.cpp:41
std::vector< ResonancePdf * > resonances
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
#define GOOFIT_INFO(...)
Definition: Log.h:10
Nonresonant constructor.
Definition: ResonancePdf.h:117
Variable constantZero("constantZero", 0)
void fillDataSetMC(UnbinnedDataSet &dataset, size_t nTotal)
Fill a dataset with MC events.
Definition: DalitzPlotter.h:64
Variable neutrlM("neutrlM", piZeroMass)
Variable motherM("motherM", _mD0)
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
Observable * m12
Special class for observables. Used in DataSets.
Definition: Variable.h:109
const fptype _mD02inv
Definition: dalitz.cpp:39
void getToyData(std::string toyFileName, GooFit::Application &app, DataSet &data)
Definition: dalitz.cpp:55
std::string toyFileName
DecayInfo3t dtop0pp
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
fptype cpuGetM23(fptype massPZ, fptype massPM)
Definition: dalitz.cpp:51
UnbinnedDataSet * data
const fptype piPlusMass
Definition: dalitz.cpp:40
const fptype _mD02
Definition: dalitz.cpp:38
__host__ void setData(DataSet *data)
Observable * m13
const fptype _mD0
Definition: dalitz.cpp:37
Variable constantOne("constantOne", 1)
Variable chargeM("chargeM", piPlusMass)
TCanvas foo
Definition: chisquare.cpp:22
Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2)
EventNumber * eventNumber
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
void makeToyData(DalitzPlotter &dplotter, UnbinnedDataSet &data)
Definition: dalitz.cpp:123
const std::vector< Observable > & getObservables() const
Definition: DataSet.cpp:49
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
int main(int argc, char **argv)
Definition: dalitz.cpp:348
__host__ void setDataSize(unsigned int dataSize, unsigned int evtSize=3)
ROOT::Minuit2::FunctionMinimum fit()
This runs the fit.
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
virtual void addEvent()=0
Relativistic Breit-Wigner.
Definition: ResonancePdf.h:68
DalitzPlotPdf * makeSignalPdf(Observable m12, Observable m13, EventNumber eventNumber, GooPdf *eff=0)
Definition: dalitz.cpp:125
This class makes it easy to make plots over 3 body Dalitz PDFs. You can use ROOT style value access o...
Definition: DalitzPlotter.h:19
int runToyFit(DalitzPlotPdf *signal, UnbinnedDataSet *data)
Definition: dalitz.cpp:319
Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8)
Variable massSum("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
std::string get_filename(const std::string &input_str, std::string base="") const
bool fitMasses