GooFit  v2.1.3
Classes | Enumerations | Functions | Variables
pipipi0DPFit.cpp File Reference
#include <TCanvas.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TLegend.h>
#include <TLine.h>
#include <TRandom.h>
#include <TRandom3.h>
#include <TStyle.h>
#include <TText.h>
#include <cassert>
#include <climits>
#include <fstream>
#include <goofit/Application.h>
#include <goofit/PDFs/GooPdf.h>
#include <goofit/PDFs/basic/BinTransformPdf.h>
#include <goofit/PDFs/basic/CrystalBallPdf.h>
#include <goofit/PDFs/basic/ExpGausPdf.h>
#include <goofit/PDFs/basic/ExpPdf.h>
#include <goofit/PDFs/basic/GaussianPdf.h>
#include <goofit/PDFs/basic/JohnsonSUPdf.h>
#include <goofit/PDFs/basic/PolynomialPdf.h>
#include <goofit/PDFs/basic/SmoothHistogramPdf.h>
#include <goofit/PDFs/basic/StepPdf.h>
#include <goofit/PDFs/basic/TrigThresholdPdf.h>
#include <goofit/PDFs/basic/VoigtianPdf.h>
#include <goofit/PDFs/combine/AddPdf.h>
#include <goofit/PDFs/combine/ConvolutionPdf.h>
#include <goofit/PDFs/combine/EventWeightedAddPdf.h>
#include <goofit/PDFs/combine/MappedPdf.h>
#include <goofit/PDFs/combine/ProdPdf.h>
#include <goofit/PDFs/physics/DalitzVetoPdf.h>
#include <goofit/PDFs/physics/IncoherentSumPdf.h>
#include <goofit/PDFs/physics/ResonancePdf.h>
#include <goofit/PDFs/physics/TddpPdf.h>
#include <goofit/PDFs/physics/ThreeGaussResolution_Aux.h>
#include <goofit/PDFs/physics/TruthResolution_Aux.h>
#include <goofit/Variable.h>
#include <goofit/utilities/Uncertain.h>
#include <goofit/fitting/FitManagerMinuit1.h>
#include <goofit/fitting/FitManagerMinuit2.h>
#include <goofit/FunctorWriter.h>
#include <goofit/PDFs/combine/CompositePdf.h>
#include <goofit/UnbinnedDataSet.h>

Go to the source code of this file.

Classes

struct  BigBin
 
struct  ChisqInfo
 

Enumerations

enum  Bkg2Model { Histogram, Parameter, Sideband }
 
enum  Resolutions { DplotRes = 1, TimeRes = 2, Efficiency = 4 }
 

Functions

SmoothHistogramPdfmakeBackgroundHistogram (int bkgnum, std::string overridename="")
 
void makeToyDalitzPlots (GooPdf *overallSignal, std::string plotdir="./plots_from_toy_mixfit/")
 
void getBackgroundFile (int bkgType)
 
double calcGauInteg (double x1, double x2)
 
double calcToyWeight (double sigratio, double m)
 
void printMemoryStatus (std::string file, int line)
 
void loadDataFile (std::string fname, UnbinnedDataSet **setToFill=0, int effSkip=3)
 
int runBackgroundDalitzFit (int bkgType, bool plots=false)
 
void normalize (TH1F *dat)
 
fptype cpuGetM23 (fptype massPZ, fptype massPM)
 
bool cpuDalitz (fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
 
void plotLoHiSigma ()
 
void plotFit (Observable var, UnbinnedDataSet *dat, GooPdf *fit)
 
bool readWrapper (std::ifstream &reader, std::string fname=strbuffer)
 
void getToyData (float sigweight=0.9)
 
GooPdfmakeEfficiencyPdf ()
 
GooPdfmakeKzeroVeto ()
 
GooPdfmakeEfficiencyPoly ()
 
TddpPdfmakeSignalPdf (MixingTimeResolution *resolution=0, GooPdf *eff=0)
 
GooPdfmakeFlatBkgDalitzPdf (bool fixem=true)
 
int runToyFit (int ifile, int nfile, bool noPlots=true)
 
void makeFullFitVariables ()
 
GooPdfmakeSignalJSU_gg (int idx, bool fixem=true)
 
GooPdfmakeMikhailJSU_gg (bool fixem=true)
 
GooPdfmakeSigmaMap ()
 
GooPdfmake1BinSigmaMap ()
 
GooPdfmake4BinSigmaMap ()
 
void coarseBin (TH2F &dalitzHist, int grain)
 
ChisqInfogetAdaptiveChisquare (TH2F *datPlot, TH2F *pdfPlot)
 
void makeDalitzPlots (GooPdf *overallSignal, std::string plotdir="./plots_from_mixfit/")
 
GooPdfmake_m23_transform ()
 
GooPdfmakeSigmaHists ()
 
GooPdfmakeBkg_sigma_strips (int bkgnum)
 
void createWeightHistogram ()
 
GooPdfmakeOverallSignal ()
 
int runTruthMCFit (std::string fname, bool noPlots=true)
 
int runGeneratedMCFit (std::string fname, int genResolutions, double dplotres)
 
GooPdfmakeBkg2_sigma ()
 
GooPdfmakeBkg4_sigma ()
 
GooPdfmakeBkg3_sigma ()
 
GooPdfmakeGaussianTimePdf (int bkg)
 
GooPdfmakeExpGausTimePdf (int bkg)
 
GooPdfmakeBkg2DalitzPdf (bool fixem=true)
 
GooPdfmakeBkg3Eff ()
 
GooPdfmakeBackground3DalitzParam ()
 
GooPdfmakeBackground4DalitzParam ()
 
GooPdfmakeBkg3DalitzPdf (bool fixem=true)
 
GooPdfmakeBkg4DalitzPdf (bool fixem=true)
 
int runCanonicalFit (std::string fname, bool noPlots=true)
 
int runSigmaFit (const char *fname)
 
int runEfficiencyFit (int which)
 
void makeTimePlots (std::string fname)
 
int runBackgroundSigmaFit (int bkgType)
 
void writeBackgroundHistograms (int bkg)
 
void set_bkg_model_from_string ()
 
void parseArg (GooFit::App *app)
 
int main (int argc, char **argv)
 

Variables

TCanvas * foo
 
TCanvas * foodal
 
UnbinnedDataSetdata = nullptr
 
UnbinnedDataSeteffdata = nullptr
 
BinnedDataSetbinEffData = nullptr
 
TH2F * weightHistogram = nullptr
 
TH2F * underlyingBins = nullptr
 
GooFit::Applicationapp_ptr
 
bool minuit1
 
Observablem12 = nullptr
 
Observablem13 = nullptr
 
EventNumbereventNumber = nullptr
 
Observablemassd0 = nullptr
 
Observabledeltam = nullptr
 
Observabledtime = nullptr
 
Observablesigma = nullptr
 
ObservablewSig0 = nullptr
 
ObservablewBkg1 = nullptr
 
ObservablewBkg2 = nullptr
 
ObservablewBkg3 = nullptr
 
ObservablewBkg4 = nullptr
 
bool fitMasses = false
 
Variable fixedRhoMass ("rho_mass", 0.7758, 0.01, 0.7, 0.8)
 
Variable fixedRhoWidth ("rho_width", 0.1503, 0.01, 0.1, 0.2)
 
double luckyFrac = 0.5
 
double mesonRad = 1.5
 
int normBinning = 240
 
int blindSeed = 4
 
int mdslices = 1
 
double md0offset = 0
 
int md0_lower_window = -2
 
int md0_upper_window = 2
 
double deltam_lower_window = -2
 
double deltam_upper_window = 2
 
double lowerTime = -2
 
double upperTime = 3
 
double maxSigma = 0.8
 
bool polyEff = false
 
bool saveEffPlot = true
 
bool useHistogramSigma = false
 
Bkg2Model bkg2Model = Sideband
 
std::string bkg2Model_str = "sideband"
 
bool notUseBackground3Hist = false
 
bool notUseBackground4Hist = false
 
bool makePlots = false
 
int m23Slices = 6
 
bool drop_rho_1450 = false
 
bool drop_rho_1700 = false
 
bool drop_f0_980 = false
 
bool drop_f0_1370 = false
 
bool drop_f0_1500 = false
 
bool drop_f0_1710 = false
 
bool drop_f2_1270 = false
 
bool drop_f0_600 = false
 
bool gaussBkgTime = false
 
bool mikhailSetup = false
 
int bkgHistBins = 80
 
string paramUp = ""
 
string paramDn = ""
 
int bkgHistRandSeed = -1
 
const fptype _mD0 = 1.86484
 
const fptype _mD02 = _mD0 * _mD0
 
const fptype _mD02inv = 1. / _mD02
 
const fptype piPlusMass = 0.13957018
 
const fptype piZeroMass = 0.1349766
 
size_t maxEvents = 100000
 
Variable motherM ("motherM", 1.86484)
 
Variable chargeM ("chargeM", 0.13957018)
 
Variable neutrlM ("neutrlM", 0.1349766)
 
Variable constantBigM ("constantBigM", _mD02+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
 
Variable constantOne ("constantOne", 1)
 
Variable constantTwo ("constantTwo", 2)
 
Variable constantZero ("constantZero", 0)
 
Variable constantMinusOne ("constantMinusOne", -1)
 
Variable minDalitzX ("minDalitzX", pow(piPlusMass+piZeroMass, 2))
 
Variable maxDalitzX ("maxDalitzX", pow(_mD0 - piPlusMass, 2))
 
Variable minDalitzY ("minDalitzY", pow(piPlusMass+piZeroMass, 2))
 
Variable maxDalitzY ("maxDalitzY", pow(_mD0 - piPlusMass, 2))
 
Variable minDalitzZ ("minDalitzZ", pow(piPlusMass+piPlusMass, 2))
 
Variable maxDalitzZ ("maxDalitzZ", pow(_mD0 - piZeroMass, 2))
 
std::vector< Variableweights
 
std::vector< Observableobsweights
 
std::vector< PdfBase * > comps
 
TH1F * dataTimePlot = nullptr
 
TH1F * loM23Sigma = nullptr
 
TH1F * hiM23Sigma = nullptr
 
TddpPdfsignalDalitz = nullptr
 
IncoherentSumPdfincsum1 = nullptr
 
IncoherentSumPdfincsum2 = nullptr
 
IncoherentSumPdfincsum3 = nullptr
 
IncoherentSumPdfincsum4 = nullptr
 
IncoherentSumPdfincsum5 = nullptr
 
IncoherentSumPdfincsum6 = nullptr
 
GooPdfsig0_jsugg = nullptr
 
GooPdfbkg2_jsugg = nullptr
 
GooPdfbkg3_jsugg = nullptr
 
GooPdfbkg4_jsugg = nullptr
 
Variable massSum ("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
 
GooPdfkzero_veto = nullptr
 
bool doToyStudy = false
 
float md0_toy_width = 0.0075
 
float md0_toy_mean = 1.8654
 
const float toySigFraction = 0.6
 
const float toyBkgTimeMean = 0.0
 
const float toyBkgTimeRMS = 0.7
 
std::string toyFileName
 
char strbuffer [1000]
 
double intGaus = -1
 
DecayInfo3t dtop0pp
 
const int numSigmaBins = 36
 
TH1F ** sigma_dat_hists = 0
 
TH1F ** sigma_pdf_hists = 0
 
UnbinnedDataSet ** sigma_data = 0
 
vector< GooPdf * > jsuList
 

Enumeration Type Documentation

◆ Bkg2Model

enum Bkg2Model
Enumerator
Histogram 
Parameter 
Sideband 

Definition at line 103 of file pipipi0DPFit.cpp.

◆ Resolutions

Enumerator
DplotRes 
TimeRes 
Efficiency 

Definition at line 130 of file pipipi0DPFit.cpp.

Function Documentation

◆ calcGauInteg()

double calcGauInteg ( double  x1,
double  x2 
)

Definition at line 186 of file pipipi0DPFit.cpp.

Referenced by calcToyWeight().

186  {
187  if(x1 > x2) {
188  double swp = x2;
189  x2 = x1;
190  x1 = swp;
191  }
192 
193  double sum1 = x1, sum2 = x2;
194  double value1 = x1, value2 = x2;
195 
196  for(int i = 1; i < 100; i++) {
197  value1 = value1 * x1 * x1 / (2 * i + 1);
198  sum1 += value1;
199  value2 = value2 * x2 * x2 / (2 * i + 1);
200  sum2 += value2;
201  }
202 
203  return sum2 * exp(-(x2 * x2) / 2) - sum1 * exp(-(x1 * x1) / 2);
204 }

◆ calcToyWeight()

double calcToyWeight ( double  sigratio,
double  m 
)

Definition at line 206 of file pipipi0DPFit.cpp.

References calcGauInteg(), intGaus, loadDataFile(), md0_lower_window, md0_toy_mean, md0_toy_width, md0_upper_window, printMemoryStatus(), and runBackgroundDalitzFit().

Referenced by getToyData().

206  {
207  if(intGaus < 0)
209 
210  double t = (m - md0_toy_mean) / md0_toy_width;
211  double fsig = sigratio * exp(-t * t / 2.) / intGaus;
212  double fbkg = (1 - sigratio) / ((md0_upper_window - md0_lower_window) * 0.0075 / md0_toy_width);
213  return fsig / (fsig + fbkg);
214 }
double intGaus
int md0_lower_window
int md0_upper_window
double calcGauInteg(double x1, double x2)
float md0_toy_mean
float md0_toy_width

◆ coarseBin()

void coarseBin ( TH2F &  dalitzHist,
int  grain 
)

Definition at line 1501 of file pipipi0DPFit.cpp.

References GooFit::Observable::getNumBins().

Referenced by makeDalitzPlots().

1501  {
1502  // Move from too-fine-to-see binning down to reasonable binning in Dalitz plots.
1503  for(int i = 1; i < m12->getNumBins(); i += grain) {
1504  for(int j = 1; j < m13->getNumBins(); j += grain) {
1505  double total = 0;
1506 
1507  for(int k = 0; k < grain; ++k) {
1508  for(int l = 0; l < grain; ++l) {
1509  total += dalitzHist.GetBinContent(i + k, j + l);
1510  }
1511  }
1512 
1513  for(int k = 0; k < grain; ++k) {
1514  for(int l = 0; l < grain; ++l) {
1515  dalitzHist.SetBinContent(i + k, j + l, total);
1516  }
1517  }
1518  }
1519  }
1520 }
Observable * m12
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
Observable * m13

◆ cpuDalitz()

bool cpuDalitz ( fptype  m12,
fptype  m13,
fptype  bigM,
fptype  dm1,
fptype  dm2,
fptype  dm3 
)

Definition at line 239 of file pipipi0DPFit.cpp.

Referenced by createWeightHistogram(), getAdaptiveChisquare(), getToyData(), makeDalitzPlots(), makeToyDalitzPlots(), runEfficiencyFit(), and runSigmaFit().

239  {
240  if(m12 < pow(dm1 + dm2, 2))
241  return false; // This m12 cannot exist, it's less than the square of the (1,2) particle mass.
242 
243  if(m12 > pow(bigM - dm3, 2))
244  return false; // This doesn't work either, there's no room for an at-rest 3 daughter.
245 
246  // Calculate energies of 1 and 3 particles in m12 rest frame.
247  fptype e1star = 0.5 * (m12 - dm2 * dm2 + dm1 * dm1) / sqrt(m12);
248  fptype e3star = 0.5 * (bigM * bigM - m12 - dm3 * dm3) / sqrt(m12);
249 
250  // Bounds for m13 at this value of m12.
251  fptype minimum
252  = pow(e1star + e3star, 2) - pow(sqrt(e1star * e1star - dm1 * dm1) + sqrt(e3star * e3star - dm3 * dm3), 2);
253 
254  if(m13 < minimum)
255  return false;
256 
257  fptype maximum
258  = pow(e1star + e3star, 2) - pow(sqrt(e1star * e1star - dm1 * dm1) - sqrt(e3star * e3star - dm3 * dm3), 2);
259 
260  if(m13 > maximum)
261  return false;
262 
263  return true;
264 }
double fptype
Observable * m12
Observable * m13

◆ cpuGetM23()

fptype cpuGetM23 ( fptype  massPZ,
fptype  massPM 
)

Definition at line 235 of file pipipi0DPFit.cpp.

References _mD02, piPlusMass, and piZeroMass.

Referenced by loadDataFile(), makeDalitzPlots(), and runSigmaFit().

235  {
236  return (_mD02 + piZeroMass * piZeroMass + piPlusMass * piPlusMass + piPlusMass * piPlusMass - massPZ - massPM);
237 }
const fptype piPlusMass
const fptype piZeroMass
const fptype _mD02

◆ createWeightHistogram()

void createWeightHistogram ( )

Definition at line 2600 of file pipipi0DPFit.cpp.

References _mD0, cpuDalitz(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::Indexable::getUpperLimit(), piPlusMass, piZeroMass, underlyingBins, and weightHistogram.

Referenced by makeOverallSignal(), and runGeneratedMCFit().

2600  {
2601  weightHistogram = new TH2F("weightHistogram",
2602  "",
2603  m12->getNumBins(),
2604  m12->getLowerLimit(),
2605  m12->getUpperLimit(),
2606  m13->getNumBins(),
2607  m13->getLowerLimit(),
2608  m13->getUpperLimit());
2609  weightHistogram->SetStats(false);
2610  double step12 = (m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins();
2611  double step13 = (m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins();
2612 
2613  for(int i = 1; i <= m12->getNumBins(); ++i) {
2614  for(int j = 1; j < m13->getNumBins(); ++j) {
2615  double maxCount = 0;
2616  double count = 0;
2617 
2618  for(double currM12 = m12->getLowerLimit() + step12 * (i - 1) + 0.05 * step12;
2619  currM12 < m12->getLowerLimit() + step12 * i;
2620  currM12 += 0.1 * step12) {
2621  for(double currM13 = m13->getLowerLimit() + step13 * (j - 1) + 0.05 * step13;
2622  currM13 < m13->getLowerLimit() + step13 * j;
2623  currM13 += 0.1 * step13) {
2624  maxCount++;
2625 
2626  if(!cpuDalitz(currM12, currM13, _mD0, piZeroMass, piPlusMass, piPlusMass))
2627  continue;
2628 
2629  count++;
2630  }
2631  }
2632 
2633  if(0.1 > maxCount)
2634  continue;
2635 
2636  count /= maxCount;
2637  weightHistogram->SetBinContent(i, j, count);
2638  }
2639  }
2640 
2641  underlyingBins = new TH2F("underlyingBins",
2642  "",
2643  m12->getNumBins(),
2644  m12->getLowerLimit(),
2645  m12->getUpperLimit(),
2646  m13->getNumBins(),
2647  m13->getLowerLimit(),
2648  m13->getUpperLimit());
2649  underlyingBins->SetStats(false);
2650 }
TH2F * underlyingBins
const fptype piPlusMass
bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
Observable * m12
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
const fptype _mD0
Observable * m13
const fptype piZeroMass
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
TH2F * weightHistogram

◆ getAdaptiveChisquare()

ChisqInfo* getAdaptiveChisquare ( TH2F *  datPlot,
TH2F *  pdfPlot 
)

Definition at line 1563 of file pipipi0DPFit.cpp.

References _mD0, ChisqInfo::chisq, ChisqInfo::ChisqInfo(), ChisqInfo::contribPlot, cpuDalitz(), ChisqInfo::dof, BigBin::getContent(), BigBin::height, piPlusMass, piZeroMass, BigBin::width, BigBin::xbin, and BigBin::ybin.

Referenced by makeDalitzPlots().

1563  {
1564  bool acceptable = false;
1565  int binSize = 1;
1566  vector<BigBin> binlist;
1567  double limit = 26;
1568 
1569  while(!acceptable) {
1570  binlist.clear();
1571  std::cout << "Attempting bin generation with size " << binSize << std::endl;
1572 
1573  for(int xbin = 1; xbin <= datPlot->GetNbinsX(); xbin += binSize) {
1574  for(int ybin = 1; ybin <= datPlot->GetNbinsY(); ybin += binSize) {
1575  double lox = datPlot->GetXaxis()->GetBinLowEdge(xbin + 0);
1576  double hix = datPlot->GetXaxis()->GetBinLowEdge(xbin + 1 + binSize);
1577  double loy = datPlot->GetYaxis()->GetBinLowEdge(ybin + 0);
1578  double hiy = datPlot->GetYaxis()->GetBinLowEdge(ybin + 1 + binSize);
1579  bool corner = false;
1580 
1581  if(cpuDalitz(lox, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1582  corner = true;
1583  else if(cpuDalitz(lox, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1584  corner = true;
1585  else if(cpuDalitz(hix, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1586  corner = true;
1587  else if(cpuDalitz(hix, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1588  corner = true;
1589 
1590  if(!corner)
1591  continue;
1592 
1593  BigBin curr;
1594  curr.xbin = xbin;
1595  curr.ybin = ybin;
1596  curr.width = curr.height = binSize;
1597  binlist.push_back(curr);
1598  }
1599  }
1600 
1601  acceptable = true;
1602 
1603  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1604  if((*bin).getContent(datPlot) >= limit)
1605  continue;
1606 
1607  acceptable = false;
1608  binSize *= 2;
1609  std::cout << "Couldn't get good bins, retry.\n";
1610  break;
1611  }
1612  }
1613 
1614  std::cout << "Good bins at size " << binSize << ", beginning splits.\n";
1615 
1616  // Now attempt to split bins.
1617  int numSplits = 1;
1618 
1619  while(0 < numSplits) {
1620  numSplits = 0;
1621  vector<BigBin> newbins;
1622 
1623  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1624  if(1 == (*bin).width * (*bin).height) {
1625  newbins.push_back(*bin);
1626  continue;
1627  }
1628 
1629  BigBin lolef;
1630  BigBin lorig;
1631  BigBin hilef;
1632  BigBin hirig;
1633  lolef.xbin = (*bin).xbin;
1634  hilef.xbin = (*bin).xbin;
1635  lorig.xbin = (*bin).xbin + (*bin).width / 2;
1636  hirig.xbin = (*bin).xbin + (*bin).width / 2;
1637  lolef.ybin = (*bin).ybin;
1638  hilef.ybin = (*bin).ybin + (*bin).height / 2;
1639  lorig.ybin = (*bin).ybin;
1640  hirig.ybin = (*bin).ybin + (*bin).height / 2;
1641 
1642  lolef.width = (*bin).width / 2;
1643  lorig.width = (*bin).width / 2;
1644  hilef.width = (*bin).width / 2;
1645  hirig.width = (*bin).width / 2;
1646  lolef.height = (*bin).height / 2;
1647  lorig.height = (*bin).height / 2;
1648  hilef.height = (*bin).height / 2;
1649  hirig.height = (*bin).height / 2;
1650 
1651  int mask = 0;
1652 
1653  if(limit < lolef.getContent(datPlot))
1654  mask += 1;
1655 
1656  if(limit < lorig.getContent(datPlot))
1657  mask += 2;
1658 
1659  if(limit < hilef.getContent(datPlot))
1660  mask += 4;
1661 
1662  if(limit < hirig.getContent(datPlot))
1663  mask += 8;
1664 
1665  if(mask != 15) {
1666  newbins.push_back(*bin);
1667  } else {
1668  newbins.push_back(lolef);
1669  newbins.push_back(lorig);
1670  newbins.push_back(hilef);
1671  newbins.push_back(hirig);
1672  numSplits++;
1673  }
1674  }
1675 
1676  binlist.clear();
1677 
1678  for(vector<BigBin>::iterator i = newbins.begin(); i != newbins.end(); ++i)
1679  binlist.push_back(*i);
1680 
1681  std::cout << "Split " << numSplits << " bins.\n";
1682  }
1683 
1684  ChisqInfo *ret = new ChisqInfo();
1685  ret->dof = binlist.size();
1686  ret->contribPlot = new TH2F("contribPlot",
1687  "",
1688  datPlot->GetNbinsX(),
1689  datPlot->GetXaxis()->GetBinLowEdge(1),
1690  datPlot->GetXaxis()->GetBinLowEdge(datPlot->GetNbinsX() + 1),
1691  datPlot->GetNbinsY(),
1692  datPlot->GetYaxis()->GetBinLowEdge(1),
1693  datPlot->GetYaxis()->GetBinLowEdge(datPlot->GetNbinsY() + 1));
1694 
1695  double totalDat = 0;
1696  double totalPdf = 0;
1697 
1698  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1699  double dat = (*bin).getContent(datPlot);
1700  double pdf = (*bin).getContent(pdfPlot);
1701  double term = (dat - pdf) / sqrt(dat);
1702  ret->chisq += term * term;
1703 
1704  /*
1705  std::cout << "Bin (" << (*bin).xbin << ", " << (*bin).ybin << ") "
1706  << (*bin).width << " " << (*bin).height << " : "
1707  << dat << " " << pdf << " "
1708  << term << std::endl;
1709  */
1710  for(int i = 0; i < (*bin).width; ++i) {
1711  for(int j = 0; j < (*bin).height; ++j) {
1712  bool corner = false;
1713  double lox = datPlot->GetXaxis()->GetBinLowEdge((*bin).xbin + i);
1714  double hix = datPlot->GetXaxis()->GetBinLowEdge((*bin).xbin + i + 1);
1715  double loy = datPlot->GetYaxis()->GetBinLowEdge((*bin).ybin + j);
1716  double hiy = datPlot->GetYaxis()->GetBinLowEdge((*bin).ybin + j + 1);
1717 
1718  if(cpuDalitz(lox, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1719  corner = true;
1720  else if(cpuDalitz(lox, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1721  corner = true;
1722  else if(cpuDalitz(hix, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1723  corner = true;
1724  else if(cpuDalitz(hix, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1725  corner = true;
1726 
1727  if(!corner)
1728  continue;
1729 
1730  ret->contribPlot->SetBinContent((*bin).xbin + i, (*bin).ybin + j, term);
1731  }
1732  }
1733 
1734  totalPdf += pdf;
1735  totalDat += dat;
1736  }
1737 
1738  return ret;
1739 }
const fptype piPlusMass
bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
double getContent(TH2F *plot)
const fptype _mD0
const fptype piZeroMass
TH2F * contribPlot

◆ getBackgroundFile()

void getBackgroundFile ( int  bkgType)

Definition at line 4667 of file pipipi0DPFit.cpp.

References bkg2Model, m23Slices, mikhailSetup, notUseBackground3Hist, notUseBackground4Hist, and Sideband.

Referenced by runBackgroundDalitzFit(), and runCanonicalFit().

4667  {
4668  if(mikhailSetup)
4669  sprintf(strbuffer, "./bkg_%i_mikhail.txt", bkgType);
4670  else {
4671  if(2 == bkgType) {
4672  if(Sideband == bkg2Model)
4673  sprintf(strbuffer, "./bkg_2_pdf_sideband_%islices.txt", m23Slices);
4674  else
4675  sprintf(strbuffer, "./bkg_2_pdf_%islices.txt", m23Slices);
4676  } else {
4677  std::string pdftype;
4678 
4679  if(((3 == bkgType) && (notUseBackground3Hist)) || ((4 == bkgType) && (notUseBackground4Hist)))
4680  pdftype = "_param";
4681 
4682  sprintf(strbuffer, "./bkg_%i_pdf%s.txt", bkgType, pdftype.c_str());
4683  }
4684  }
4685 }
int m23Slices
bool notUseBackground3Hist
char strbuffer[1000]
Bkg2Model bkg2Model
bool mikhailSetup
bool notUseBackground4Hist

◆ getToyData()

void getToyData ( float  sigweight = 0.9)

Definition at line 338 of file pipipi0DPFit.cpp.

References _mD0, GooFit::UnbinnedDataSet::addEvent(), calcToyWeight(), cpuDalitz(), GooFit::Indexable::getLowerLimit(), GooFit::DataSet::getNumEvents(), GooFit::Indexable::getUpperLimit(), GooFit::Indexable::getValue(), GOOFIT_INFO, maxEvents, md0_lower_window, md0_toy_mean, md0_upper_window, piPlusMass, piZeroMass, readWrapper(), and GooFit::Indexable::setValue().

Referenced by runToyFit().

338  {
339  if(!data) {
340  std::vector<Observable> vars;
341  vars.push_back(*m12);
342  vars.push_back(*m13);
343  vars.push_back(*dtime);
344  vars.push_back(*sigma);
345  vars.push_back(*eventNumber);
346  vars.push_back(*wSig0);
347  // vars.push_back(wBkg1);
348  // vars.push_back(wBkg2);
349  data = new UnbinnedDataSet(vars);
350  }
351 
352  std::ifstream reader;
353  readWrapper(reader, toyFileName);
354  std::string buffer;
355 
356  while(!reader.eof()) {
357  reader >> buffer;
358 
359  if(buffer == "====")
360  break;
361 
362  std::cout << buffer;
363  }
364 
365  TRandom3 donram(42);
366 
367  int nsig = 0;
368  double sigprob = 0;
369  double dummy = 0;
370  double md0 = md0_toy_mean;
371 
372  while(reader >> dummy
373 
374  // Ignoring m23, m(pi+ pi-), called m12 in processToyRoot convention.
375  // Already swapped m12 m13 according to D* charge. m12 = m(pi+pi0)
376  >> dummy >> *m12 >> *m13
377 
378  // Errors on Dalitz variables
379  >> dummy >> dummy >> dummy
380 
381  // Remaining two interesting values
382  >> *dtime >> *sigma
383 
384  // Md0, deltaM, ProbSig, Dst charge, Run, Event, Signal and four bkg fractions
385  >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy) {
386  double resolution = donram.Gaus(0, 1);
387  dtime->setValue(dtime->getValue() + resolution * sigma->getValue());
388 
390 
391  do {
392  md0 = donram.Gaus(md0_toy_mean, md0_toy_width);
393  } while(md0 <= 1.8654 + 0.0075 * md0_lower_window || md0 >= 1.8654 + 0.0075 * md0_upper_window);
394 
395  // wSig0->getValue() = sigweight;
396  wSig0->setValue(calcToyWeight(sigweight, md0));
397  sigprob += wSig0->getValue();
398  data->addEvent();
399  nsig++;
400 
401  if(data->getNumEvents() < 10) {
402  std::cout << data->getNumEvents() << " : " << m12->getValue() << " " << m13->getValue() << " "
403  << dtime->getValue() << " " << sigma->getValue() << " " << std::endl;
404  }
405 
406  if(data->getNumEvents() >= maxEvents)
407  break;
408  }
409 
410  GOOFIT_INFO("Read in {} events", data->getNumEvents())
411 
412  for(int ib = 0; ib < nsig * (1 - sigweight) / sigweight; ib++) {
413  do {
414  m12->setValue(donram.Uniform() * (m12->getUpperLimit() - m12->getLowerLimit()) + m12->getLowerLimit());
415  m13->setValue(donram.Uniform() * (m13->getUpperLimit() - m13->getLowerLimit()) + m13->getLowerLimit());
417 
418  do {
419  dtime->setValue(donram.Gaus(toyBkgTimeMean, toyBkgTimeRMS));
420  } while(!(dtime->getValue() > dtime->getLowerLimit() && dtime->getValue() < dtime->getUpperLimit()));
421 
423  md0 = donram.Uniform(1.8654 + 0.0075 * md0_lower_window, 1.8654 + 0.0075 * md0_upper_window);
424  // wSig0->getValue() = sigweight;
425  wSig0->setValue(calcToyWeight(sigweight, md0));
426  sigprob += wSig0->getValue();
427  data->addEvent();
428  }
429 
430  reader.close();
431 }
size_t getNumEvents() const
Definition: DataSet.h:47
#define GOOFIT_INFO(...)
Definition: Log.h:10
const fptype piPlusMass
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
Observable * m12
std::string toyFileName
UnbinnedDataSet * data
size_t maxEvents
Observable * wSig0
const fptype _mD0
Observable * m13
Observable * sigma
const float toyBkgTimeRMS
const fptype piZeroMass
const float toyBkgTimeMean
EventNumber * eventNumber
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
int md0_lower_window
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
int md0_upper_window
float md0_toy_mean
Observable * dtime
bool readWrapper(std::ifstream &reader, std::string fname=strbuffer)
float md0_toy_width
double calcToyWeight(double sigratio, double m)

◆ loadDataFile()

void loadDataFile ( std::string  fname,
UnbinnedDataSet **  setToFill = 0,
int  effSkip = 3 
)

Definition at line 921 of file pipipi0DPFit.cpp.

References GooFit::BinnedDataSet::addWeightedEvent(), cpuGetM23(), data, deltam, deltam_lower_window, deltam_upper_window, dtime, GooFit::Indexable::getLowerLimit(), GooFit::Indexable::getUpperLimit(), GooFit::Indexable::getValue(), luckyFrac, m12, m13, massd0, md0_lower_window, md0_upper_window, md0offset, readWrapper(), GooFit::Indexable::setValue(), sigma, underlyingBins, wBkg1, wBkg2, wBkg3, wBkg4, weightHistogram, and wSig0.

Referenced by calcToyWeight(), makeOverallSignal(), makeTimePlots(), runBackgroundDalitzFit(), runBackgroundSigmaFit(), runCanonicalFit(), runEfficiencyFit(), runGeneratedMCFit(), runSigmaFit(), and runTruthMCFit().

921  {
922  if(!setToFill)
923  setToFill = &data;
924 
925  std::vector<Observable> vars;
926  vars.push_back(*m12);
927  vars.push_back(*m13);
928  vars.push_back(*dtime);
929  vars.push_back(*sigma);
930  vars.push_back(*eventNumber);
931  vars.push_back(*wSig0);
932  vars.push_back(*wBkg1);
933  vars.push_back(*wBkg2);
934  vars.push_back(*wBkg3);
935  vars.push_back(*wBkg4);
936 
937  if(massd0)
938  vars.push_back(*massd0);
939 
940  (*setToFill) = new UnbinnedDataSet(vars);
941  std::ifstream reader;
942  readWrapper(reader, fname.c_str());
943  std::string buffer;
944 
945  while(!reader.eof()) {
946  reader >> buffer;
947 
948  if(buffer == "====")
949  break;
950 
951  std::cout << buffer;
952  }
953 
954  double integralWeights[5] = {0, 0, 0, 0, 0};
955 
956  double dummy = 0;
957  double mass = 0;
958  double delm = 0;
959  int events = 0;
960 
961  while(!reader.eof()) {
962  reader >> dummy;
963 
964  if(reader.eof())
965  break;
966 
967  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
968  reader >> *m12; // Already swapped according to D* charge
969  reader >> *m13;
970 
971  // Errors on Dalitz variables
972  reader >> dummy;
973  reader >> dummy;
974  reader >> dummy;
975 
976  reader >> *dtime;
977  reader >> *sigma;
978 
979  if(massd0) {
980  reader >> *massd0; // Md0
981  } else
982  reader >> mass;
983 
984  if(deltam) {
985  reader >> *deltam;
986  } else
987  reader >> delm;
988 
989  reader >> dummy; // ProbSig
990  reader >> dummy; // Dst charge
991  reader >> dummy; // Run
992  reader >> dummy; // Event
993  reader >> *wSig0;
994  reader >> *wBkg1;
995  reader >> *wBkg2;
996  reader >> *wBkg3;
997  reader >> *wBkg4;
998 
999  if(massd0) {
1000  if(massd0->getValue() <= massd0->getLowerLimit())
1001  continue;
1002 
1003  if(massd0->getValue() >= massd0->getUpperLimit())
1004  continue;
1005  } else {
1006  // Enforce signal box on all data sets!
1007  if(mass <= 1.8654 + 0.0075 * md0_lower_window + md0offset)
1008  continue;
1009 
1010  if(mass >= 1.8654 + 0.0075 * md0_upper_window + md0offset)
1011  continue;
1012  }
1013 
1014  if(deltam) {
1015  if(deltam->getValue() >= deltam->getUpperLimit())
1016  continue;
1017 
1018  if(deltam->getValue() <= deltam->getLowerLimit())
1019  continue;
1020  } else {
1021  if(delm >= 0.1454 + 0.0003 * deltam_upper_window)
1022  continue;
1023 
1024  if(delm <= 0.1454 + 0.0003 * deltam_lower_window)
1025  continue;
1026  }
1027 
1028  if(dtime->getValue() < dtime->getLowerLimit())
1029  continue;
1030 
1031  if(dtime->getValue() > dtime->getUpperLimit())
1032  continue;
1033 
1034  if(sigma->getValue() > sigma->getUpperLimit())
1035  continue; // Lower limit is 0, and it can't be lower than that, so whatever.
1036 
1037  integralWeights[0] += wSig0->getValue();
1038  integralWeights[1] += wBkg1->getValue();
1039  integralWeights[2] += wBkg2->getValue();
1040  integralWeights[3] += wBkg3->getValue();
1041  integralWeights[4] += wBkg4->getValue();
1042  eventNumber->setValue((*setToFill)->getNumEvents());
1043 
1044  // See comments in TddpPdf.hh for explanation of this.
1045  double mistag = wSig0->getValue() + wBkg1->getValue() * luckyFrac;
1046  wSig0->setValue(wSig0->getValue() + wBkg1->getValue());
1047  wBkg1->setValue(mistag / wSig0->getValue());
1048 
1049  if((binEffData) && (0 == events % effSkip)) {
1050  double weight = weightHistogram->GetBinContent(weightHistogram->FindBin(m12->getValue(), m13->getValue()));
1051  // weight = 1;
1052  binEffData->addWeightedEvent(weight);
1053 
1054  // binEffData->addEvent();
1055  if(underlyingBins)
1056  underlyingBins->Fill(m12->getValue(), m13->getValue(), weight);
1057  } else
1058  (*setToFill)->addEvent();
1059 
1060  events++;
1061 
1062  if(loM23Sigma) {
1063  double currM23 = cpuGetM23(m12->getValue(), m13->getValue());
1064 
1065  if(currM23 < 1.5)
1066  loM23Sigma->Fill(sigma->getValue());
1067  else
1068  hiM23Sigma->Fill(sigma->getValue());
1069 
1070  // std::cout << "Filled sigma with " << sigma->getValue() << " " << currM23 << std::endl;
1071  }
1072 
1073  if(wSig0->getValue() > 1.0)
1074  std::cout << "Problem with event " << (*setToFill)->getNumEvents() << ", too-large signal weight "
1075  << wSig0->getValue() << std::endl;
1076 
1077  /*
1078  if ((*setToFill)->getNumEvents() < 10) {
1079  std::cout << (*setToFill)->getNumEvents() << " : "
1080  << m12->getValue() << " "
1081  << m13->getValue() << " "
1082  << dtime->getValue() << " "
1083  << sigma->getValue() << " "
1084  << wSig0->getValue() << " "
1085  << wBkg1->getValue() << " "
1086  << wBkg2->getValue() << " "
1087  << wBkg3->getValue() << " "
1088  << wBkg4->getValue() << " "
1089  << std::endl;
1090  }
1091  */
1092  }
1093 
1094  std::cout << "Loaded " << (*setToFill)->getNumEvents() << " events.\n";
1095  std::cout << "Integrals: " << integralWeights[0] << " " << integralWeights[1] << " " << integralWeights[2] << " "
1096  << integralWeights[3] << " " << integralWeights[4] << "\n";
1097 }
Observable * wBkg3
TH2F * underlyingBins
Observable * massd0
double deltam_lower_window
TH1F * loM23Sigma
fptype cpuGetM23(fptype massPZ, fptype massPM)
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
Observable * m12
Observable * deltam
TH1F * hiM23Sigma
UnbinnedDataSet * data
Observable * wSig0
BinnedDataSet * binEffData
Observable * m13
Observable * sigma
double deltam_upper_window
EventNumber * eventNumber
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
Observable * wBkg4
int md0_lower_window
double luckyFrac
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
void addWeightedEvent(double weight) override
int md0_upper_window
double md0offset
Observable * wBkg2
Observable * wBkg1
Observable * dtime
bool readWrapper(std::ifstream &reader, std::string fname=strbuffer)
TH2F * weightHistogram

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 4920 of file pipipi0DPFit.cpp.

References data, DplotRes, foo, foodal, GOOFIT_PARSE, m23Slices, makePlots, makeTimePlots(), maxEvents, minuit1, parseArg(), runBackgroundDalitzFit(), runBackgroundSigmaFit(), runCanonicalFit(), runEfficiencyFit(), runGeneratedMCFit(), runSigmaFit(), runToyFit(), runTruthMCFit(), set_bkg_model_from_string(), and writeBackgroundHistograms().

4920  {
4921  gStyle->SetCanvasBorderMode(0);
4922  gStyle->SetCanvasColor(10);
4923  gStyle->SetFrameFillColor(10);
4924  gStyle->SetFrameBorderMode(0);
4925  gStyle->SetPadColor(0);
4926  gStyle->SetTitleColor(1);
4927  gStyle->SetStatColor(0);
4928  gStyle->SetFillColor(0);
4929  gStyle->SetFuncWidth(1);
4930  gStyle->SetLineWidth(1);
4931  gStyle->SetLineColor(1);
4932  gStyle->SetPalette(1, 0);
4933  foo = new TCanvas();
4934  foodal = new TCanvas();
4935  foodal->Size(10, 10);
4936 
4937  int retval = 0;
4938 
4939  GooFit::Application app("pipipi0 Dalitz fit example", argc, argv);
4940  app_ptr = &app;
4941  app.require_subcommand();
4942 
4943  app.add_flag("--minuit1", minuit1, "Use Minuit 1 instead of Minuit 2");
4944 
4945  std::string data;
4946  int sample = 0;
4947  int load = 1;
4948  bool plots;
4949  int genResolutions = 0;
4950  double dplotres = 0;
4951 
4952  auto toy = app.add_subcommand("toy", "Toy MC Performance evaluation");
4953  toy->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4954  toy->add_option("-l,--load,load", load, "Number of times to load", true);
4955  toy->add_option("-m,--max", maxEvents, "Maximum number of events to read", true);
4956  toy->add_flag("-p,--plot", plots, "Also make plots");
4957  toy->set_callback([&]() { retval = runToyFit(sample, load, plots); });
4958 
4959  auto truth_fit = app.add_subcommand("truth", "Truth Monte Carlo fit");
4960  truth_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4961  truth_fit->set_callback([&]() { retval = runTruthMCFit(data, false); });
4962 
4963  auto sigma_fit = app.add_subcommand("sigma", "Run sigma fit");
4964  sigma_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4965  sigma_fit->add_option("-s,--slices,slices", m23Slices, "m23 slices")->required();
4966  sigma_fit->set_callback([&]() { retval = runSigmaFit(data.c_str()); });
4967 
4968  auto efficiency_fit = app.add_subcommand("efficiency", "Run efficiency fit");
4969  efficiency_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4970  efficiency_fit->set_callback([&]() { retval = runEfficiencyFit(sample); });
4971 
4972  auto canonical_fit = app.add_subcommand("canonical", "Run the canonical fit");
4973  canonical_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4974  parseArg(canonical_fit);
4975  canonical_fit->set_callback([&]() {
4977  retval = runCanonicalFit(data, !makePlots);
4978  });
4979 
4980  auto background_dalitz_fit = app.add_subcommand("background_dalitz", "Run the background Dalitz fit");
4981  background_dalitz_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4982  parseArg(background_dalitz_fit);
4983  background_dalitz_fit->set_callback([&]() {
4985  retval = runBackgroundDalitzFit(sample, true);
4986  });
4987 
4988  auto background_sigma_fit = app.add_subcommand("background_sigma", "Run background sigma fit");
4989  background_sigma_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4990  background_sigma_fit->set_callback([&]() { retval = runBackgroundSigmaFit(sample); });
4991 
4992  auto write_background_histograms = app.add_subcommand("background_histograms", "Write background histograms");
4993  write_background_histograms->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4994  write_background_histograms->set_callback([&]() { writeBackgroundHistograms(sample); });
4995 
4996  auto run_gen_mc_fit = app.add_subcommand("run_gen_mc", "Run generated Monte Carlo fit");
4997  run_gen_mc_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4998  run_gen_mc_fit->add_option("-g,--genres,gen-resolutions", genResolutions)->required();
4999  run_gen_mc_fit->add_option("-p,--dplotres,dplotres", dplotres);
5000  run_gen_mc_fit->set_callback([&]() {
5001  if(!(DplotRes & genResolutions))
5002  dplotres = 0;
5003  retval = runGeneratedMCFit(data, genResolutions, dplotres);
5004  });
5005 
5006  auto make_time_plots = app.add_subcommand("make_time_plots", "Make time plots");
5007  make_time_plots->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
5008  make_time_plots->set_callback([&]() { makeTimePlots(data); });
5009 
5010  GOOFIT_PARSE(app);
5011 
5012  return retval;
5013 }
int runGeneratedMCFit(std::string fname, int genResolutions, double dplotres)
int m23Slices
int runSigmaFit(const char *fname)
void parseArg(GooFit::App *app)
int runToyFit(int ifile, int nfile, bool noPlots=true)
UnbinnedDataSet * data
size_t maxEvents
TCanvas * foodal
int runCanonicalFit(std::string fname, bool noPlots=true)
int runBackgroundDalitzFit(int bkgType, bool plots=false)
int runBackgroundSigmaFit(int bkgType)
bool minuit1
bool makePlots
void set_bkg_model_from_string()
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
void makeTimePlots(std::string fname)
int runTruthMCFit(std::string fname, bool noPlots=true)
TCanvas * foo
int runEfficiencyFit(int which)
GooFit::Application * app_ptr
void writeBackgroundHistograms(int bkg)

◆ make1BinSigmaMap()

GooPdf* make1BinSigmaMap ( )

Definition at line 1345 of file pipipi0DPFit.cpp.

References GooFit::UnbinnedDataSet::addEvent(), GooFit::FitManagerMinuit2::fit(), GooFit::FitManagerMinuit1::fit(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::DataSet::getNumEvents(), GooFit::Indexable::getUpperLimit(), GooFit::UnbinnedDataSet::getValue(), GooFit::Indexable::getValue(), makeSignalJSU_gg(), minuit1, GooFit::PdfBase::setData(), GooFit::GooPdf::setParameterConstantness(), and GooFit::Indexable::setValue().

1345  {
1346  sigma_dat_hists = new TH1F *[1];
1347  sigma_pdf_hists = new TH1F *[1];
1348  sigma_data = new UnbinnedDataSet *[1];
1349 
1350  std::vector<Observable> vars;
1351  vars.push_back(*sigma);
1352 
1353  for(int i = 0; i < 1; ++i) {
1354  sprintf(strbuffer, "sigma_data_hist_%i", i);
1355  sigma_dat_hists[i]
1356  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1357  sigma_dat_hists[i]->SetStats(false);
1358  sigma_dat_hists[i]->SetMarkerStyle(8);
1359  sigma_dat_hists[i]->SetMarkerSize(0.7);
1360  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1361  sigma_pdf_hists[i]
1362  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1363  sigma_pdf_hists[i]->SetStats(false);
1364  sigma_pdf_hists[i]->SetLineWidth(3);
1365  sigma_pdf_hists[i]->SetLineColor(kBlue);
1366 
1367  sigma_data[i] = new UnbinnedDataSet(vars);
1368  }
1369 
1370  int numEvents = data->getNumEvents();
1371 
1372  for(int i = 0; i < numEvents; ++i) {
1373  m12->setValue(data->getValue(*m12, i));
1374  m13->setValue(data->getValue(*m13, i));
1375  sigma->setValue(data->getValue(*sigma, i));
1376 
1377  int overallbin = 0;
1378  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1379  sigma_data[overallbin]->addEvent();
1380  }
1381 
1382  // vector<GooPdf*> jsuList;
1383  for(int i = 0; i < 1; ++i) {
1384  GooPdf *js = makeSignalJSU_gg(i, false);
1385  jsuList.push_back(js);
1386 
1387  if(0 == sigma_data[i]->getNumEvents())
1388  js->setParameterConstantness(true);
1389  else {
1390  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1391  js->setData(sigma_data[i]);
1392  if(minuit1) {
1393  GooFit::FitManagerMinuit1 currpdf(js);
1394  currpdf.fit();
1395  } else {
1396  GooFit::FitManagerMinuit2 currpdf(js);
1397  currpdf.fit();
1398  }
1399  js->setParameterConstantness(true);
1400  // js->clearCurrentFit();
1401  std::cout << "Done with sigma box " << i << "\n";
1402  }
1403  }
1404 
1405  vector<Observable> obses;
1406  obses.push_back(*m12);
1407  obses.push_back(*m13);
1408  vector<double> limits;
1409  limits.push_back(0);
1410  limits.push_back(0);
1411  vector<double> binSizes;
1412  binSizes.push_back(3);
1413  binSizes.push_back(3);
1414  vector<int> numBins;
1415  numBins.push_back(1);
1416  numBins.push_back(1);
1417  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1418 
1419  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1420 }
size_t getNumEvents() const
Definition: DataSet.h:47
vector< GooPdf * > jsuList
__host__ void setParameterConstantness(bool constant=true)
char strbuffer[1000]
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
TH1F ** sigma_pdf_hists
Observable * m12
fptype getValue(const Observable &var, size_t idx) const
Get the value at a specific variable and event number.
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * data
__host__ void setData(DataSet *data)
Observable * m13
Observable * sigma
bool minuit1
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
TH1F ** sigma_dat_hists
GooPdf * makeSignalJSU_gg(int idx, bool fixem=true)
UnbinnedDataSet ** sigma_data

◆ make4BinSigmaMap()

GooPdf* make4BinSigmaMap ( )

Definition at line 1422 of file pipipi0DPFit.cpp.

References GooFit::UnbinnedDataSet::addEvent(), GooFit::FitManagerMinuit2::fit(), GooFit::FitManagerMinuit1::fit(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::DataSet::getNumEvents(), GooFit::Indexable::getUpperLimit(), GooFit::UnbinnedDataSet::getValue(), GooFit::Indexable::getValue(), makeSignalJSU_gg(), minuit1, GooFit::PdfBase::setData(), GooFit::GooPdf::setParameterConstantness(), and GooFit::Indexable::setValue().

1422  {
1423  sigma_dat_hists = new TH1F *[4];
1424  sigma_pdf_hists = new TH1F *[4];
1425  sigma_data = new UnbinnedDataSet *[4];
1426 
1427  std::vector<Observable> vars;
1428  vars.push_back(*sigma);
1429 
1430  for(int i = 0; i < 4; ++i) {
1431  sprintf(strbuffer, "sigma_data_hist_%i", i);
1432  sigma_dat_hists[i]
1433  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1434  sigma_dat_hists[i]->SetStats(false);
1435  sigma_dat_hists[i]->SetMarkerStyle(8);
1436  sigma_dat_hists[i]->SetMarkerSize(0.7);
1437  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1438  sigma_pdf_hists[i]
1439  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1440  sigma_pdf_hists[i]->SetStats(false);
1441  sigma_pdf_hists[i]->SetLineWidth(3);
1442  sigma_pdf_hists[i]->SetLineColor(kBlue);
1443 
1444  sigma_data[i] = new UnbinnedDataSet(vars);
1445  }
1446 
1447  int numEvents = data->getNumEvents();
1448 
1449  for(int i = 0; i < numEvents; ++i) {
1450  m12->setValue(data->getValue(*m12, i));
1451  m13->setValue(data->getValue(*m13, i));
1452  sigma->setValue(data->getValue(*sigma, i));
1453 
1454  int xbin = (int)floor(m12->getValue() / 1.5);
1455  int ybin = (int)floor(m13->getValue() / 1.5);
1456  int overallbin = ybin * 2 + xbin;
1457  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1458  sigma_data[overallbin]->addEvent();
1459  }
1460 
1461  // vector<GooPdf*> jsuList;
1462  for(int i = 0; i < 4; ++i) {
1463  GooPdf *js = makeSignalJSU_gg(i, false);
1464  jsuList.push_back(js);
1465 
1466  if(0 == sigma_data[i]->getNumEvents())
1467  js->setParameterConstantness(true);
1468  else {
1469  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1470  js->setData(sigma_data[i]);
1471  if(minuit1) {
1472  GooFit::FitManagerMinuit1 currpdf(js);
1473  currpdf.fit();
1474  } else {
1475  GooFit::FitManagerMinuit2 currpdf(js);
1476  currpdf.fit();
1477  }
1478  js->setParameterConstantness(true);
1479  // js->clearCurrentFit();
1480  std::cout << "Done with sigma box " << i << "\n";
1481  }
1482  }
1483 
1484  vector<Observable> obses;
1485  obses.push_back(*m12);
1486  obses.push_back(*m13);
1487  vector<double> limits;
1488  limits.push_back(0);
1489  limits.push_back(0);
1490  vector<double> binSizes;
1491  binSizes.push_back(1.5);
1492  binSizes.push_back(1.5);
1493  vector<int> numBins;
1494  numBins.push_back(2);
1495  numBins.push_back(2);
1496  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1497 
1498  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1499 }
size_t getNumEvents() const
Definition: DataSet.h:47
vector< GooPdf * > jsuList
__host__ void setParameterConstantness(bool constant=true)
char strbuffer[1000]
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
TH1F ** sigma_pdf_hists
Observable * m12
fptype getValue(const Observable &var, size_t idx) const
Get the value at a specific variable and event number.
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * data
__host__ void setData(DataSet *data)
Observable * m13
Observable * sigma
bool minuit1
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
TH1F ** sigma_dat_hists
GooPdf * makeSignalJSU_gg(int idx, bool fixem=true)
UnbinnedDataSet ** sigma_data

◆ make_m23_transform()

GooPdf* make_m23_transform ( )

Definition at line 2459 of file pipipi0DPFit.cpp.

References constantBigM, constantMinusOne, constantZero, and m23Slices.

Referenced by makeBkg_sigma_strips(), and makeSigmaHists().

2459  {
2460  // This is complicated. I want to make a function such that the parameters
2461  // depend on position in the Dalitz plot. I want to use stripes of m23, because
2462  // that seems to be the real shape of the dependence, and I haven't got so much
2463  // data for bkg3; so stripes work better than boxes of m12, m13. So I need a
2464  // transform from m12, m13 to m23. Then I need a transform from m23 to a bin number.
2465  // Finally I need a transform from bin number to function. Keep the tongue straight
2466  // in the mouth, now!
2467 
2468  vector<Observable> obses;
2469  vector<Variable> offsets;
2470  vector<Variable> coefficients;
2471  vector<PdfBase *> components;
2472  vector<double> limits;
2473  vector<double> binSizes;
2474  vector<int> numBins;
2475 
2476  // First the transform to m23.
2477  // m23 = _mD02 + piZeroMass*piZeroMass + piPlusMass*piPlusMass + piPlusMass*piPlusMass - massPZ - massPM
2478 
2479  obses.push_back(*m12);
2480  obses.push_back(*m13);
2481  coefficients.push_back(constantBigM);
2482  coefficients.push_back(constantMinusOne);
2483  coefficients.push_back(constantMinusOne);
2484  offsets.push_back(constantZero);
2485  offsets.push_back(constantZero);
2486  PolynomialPdf *m23_transform = new PolynomialPdf("m23_transform", obses, coefficients, offsets, 1);
2487 
2488  // Now create the BinTransform which will make bins out of m23 values.
2489  obses.clear();
2490  coefficients.clear();
2491  obses.push_back(*m12); // Fake dependence; CompositePdf will create a fake event using this index.
2492  limits.push_back(0); // Bins of m23 start at 0.
2493  binSizes.push_back(3.0 / m23Slices);
2494  numBins.push_back(m23Slices);
2495  BinTransformPdf *m23_binMap = new BinTransformPdf("m23_binMap", obses, limits, binSizes, numBins);
2496 
2497  // Now make a composite, so that the BinTransform takes the Polynomial result as input.
2498  CompositePdf *m23_composite = new CompositePdf("m23_composite", m23_transform, m23_binMap);
2499  return m23_composite;
2500 }
Variable constantZero("constantZero", 0)
int m23Slices
Observable * m12
Variable constantBigM("constantBigM", _mD02+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
Observable * m13
Variable constantMinusOne("constantMinusOne", -1)

◆ makeBackground3DalitzParam()

GooPdf* makeBackground3DalitzParam ( )

Definition at line 3604 of file pipipi0DPFit.cpp.

References _mD0, GooFit::PdfBase::addSpecialMask(), constantMinusOne, constantOne, constantZero, GooFit::DecayInfo3::daug1Mass, GooFit::DecayInfo3::daug2Mass, GooFit::DecayInfo3::daug3Mass, GooFit::DecayInfo3::meson_radius, GooFit::DecayInfo3::motherMass, GooFit::PAIR_12, GooFit::PAIR_13, GooFit::PAIR_23, piPlusMass, piZeroMass, and GooFit::DecayInfo3::resonances.

Referenced by makeBkg3DalitzPdf().

3604  {
3605  // I can't make this thing describe the background 3 data.
3606 
3607  // GooPdf* bkg3_eff = makeBkg3Eff();
3608  weights.clear();
3609 
3610  vector<Variable> offsets;
3611  vector<Observable> observables;
3612  vector<Variable> coefficients;
3613  offsets.push_back(constantOne);
3614  offsets.push_back(constantOne);
3615 
3616  observables.push_back(*m12);
3617  observables.push_back(*m13);
3618  // Recurring factor 3 offsets division by total weight in AddPdf.
3619  double weightOffset = 1;
3620 
3621  coefficients.push_back(Variable("bkg3_x0y0", 1.00 * weightOffset));
3622  // coefficients.push_back(new Variable("bkg3_x0y0", 1.10 * weightOffset, 0.01, 0.01 * weightOffset, 1.50 *
3623  // weightOffset));
3624  coefficients.push_back(
3625  Variable("bkg3_x1y0", -0.36937 * weightOffset, 0.01, -1.50 * weightOffset, 0.00 * weightOffset));
3626  coefficients.push_back(
3627  Variable("bkg3_x2y0", 1.36184 * weightOffset, 0.01, -0.10 * weightOffset, 1.60 * weightOffset));
3628  // coefficients.push_back(new Variable("bkg3_x3y0", -0.43177 * weightOffset, 0.01,-1.60*weightOffset,
3629  // 0.60*weightOffset));
3630  coefficients.push_back(
3631  Variable("bkg3_x0y1", -0.27691 * weightOffset, 0.01, -1.50 * weightOffset, 0.00 * weightOffset));
3632  coefficients.push_back(
3633  Variable("bkg3_x1y1", 2.16029 * weightOffset, 0.01, 0.30 * weightOffset, 4.50 * weightOffset));
3634  // coefficients.push_back(new Variable("bkg3_x2y1", -2.04133 * weightOffset, 0.01,-2.50*weightOffset,
3635  // 1.50*weightOffset));
3636  coefficients.push_back(
3637  Variable("bkg3_x0y2", 1.33100 * weightOffset, 0.01, 1.00 * weightOffset, 2.00 * weightOffset));
3638  // coefficients.push_back(new Variable("bkg3_x1y2", -1.88226 * weightOffset, 0.01,-2.20*weightOffset,
3639  // 1.00*weightOffset));
3640  // coefficients.push_back(new Variable("bkg3_x0y3", -0.58920 * weightOffset, 0.01,-1.00*weightOffset,
3641  // 2.00*weightOffset));
3642  // PolynomialPdf* poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 3);
3643  PolynomialPdf *poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 2);
3644 
3645  // coefficients.push_back(new Variable("bkg3_x0y0", 0.10 * weightOffset, 0.01, 0.01 * weightOffset, 0.20 *
3646  // weightOffset));
3647  // PolynomialPdf* poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 0);
3648 
3649  // Background 3 does not have a trig threshold in Mikhail's fit
3650  // - that is, it has one, but the dec variable is set above the
3651  // threshold, so it doesn't do anything.
3652  // That's Mikhail's fit; I'm putting one in to try to deal with
3653  // the asymmetry in the rho+.
3654  // Didn't work.
3655  // Variable* bkg3_decZmin = new Variable("bkg3_decZmin", 3.30411, 0.1, 1, 5);
3656  // Variable* bkg3_conZmin = new Variable("bkg3_conZmin", 0.29909, 0.01, 0.1, 0.9);
3657  // TrigThresholdPdf* bkg3_loZ = new TrigThresholdPdf("bkg3_loZ", m12, m13, minDalitzZ, bkg3_decZmin, bkg3_conZmin,
3658  // massSum, false);
3659 
3660  comps.clear();
3661  comps.push_back(poly);
3662  comps.push_back(kzero_veto);
3663  // comps.push_back(bkg3_eff);
3664  // comps.push_back(bkg3_loZ);
3665 
3666  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3667 
3668  // One misIDpi0.
3669  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3670  DecayInfo3 special_pi0_decay;
3671  special_pi0_decay.motherMass = _mD0;
3672  special_pi0_decay.daug1Mass = piZeroMass;
3673  special_pi0_decay.daug2Mass = piPlusMass;
3674  special_pi0_decay.daug3Mass = piPlusMass;
3675  special_pi0_decay.meson_radius = 1.5;
3676 
3677  ResonancePdf *bkg3_pi0_ref = new Resonances::Gauss(
3678  "bkg3_pi0_ref",
3679  Variable("bkg3_pi0_ref_amp", 0.01189 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3680  constantZero,
3681  Variable("bkg3_pi0_ref_mass", 1.65766, 0.01, 1.4, 1.8),
3682  Variable("bkg3_pi0_ref_width", 0.05018, 0.01, 0.02, 0.20),
3683  PAIR_23);
3684  special_pi0_decay.resonances.push_back(bkg3_pi0_ref);
3685 
3686  // Mikhail defines 'transverse Z' as y - x - (parameter).
3687  // Variable* bkg3_pi0_transZ_offset = new Variable("bkg3_pi0_transZ_offset", -0.04381, 0.001, -0.5, 0.5);
3688  Variable bkg3_pi0_transZ_offset("bkg3_pi0_transZ_offset", -0.04381);
3689  offsets.clear();
3690  observables.clear();
3691  coefficients.clear();
3692  observables.push_back(*m12);
3693  observables.push_back(*m13);
3694  coefficients.push_back(bkg3_pi0_transZ_offset);
3695  coefficients.push_back(constantMinusOne);
3696  coefficients.push_back(constantOne);
3697  offsets.push_back(constantZero);
3698  offsets.push_back(constantZero);
3699  PolynomialPdf *bkg3_pi0_transZ = new PolynomialPdf("bkg3_pi0_transZ", observables, coefficients, offsets, 1);
3700 
3701  // Now we're going to take (1 - tz^2 * (parameter)) and multiply that into the misID pi0.
3702  // Variable* bkg3_pi0_transZ_quad = new Variable("bkg3_pi0_transZ_quad", 2.12277, 0.01, -1.5, 6.0);
3703  Variable bkg3_pi0_transZ_quad("bkg3_pi0_transZ_quad", 2.12277);
3704  coefficients.clear();
3705  coefficients.push_back(constantOne);
3706  coefficients.push_back(constantZero);
3707  coefficients.push_back(bkg3_pi0_transZ_quad);
3708  // Notice the fake dependence of the polynomial on m12; in fact CompositePdf
3709  // will send it a fake event, we just have to supply a reasonable index.
3710  PolynomialPdf *bkg3_pi0_shell = new PolynomialPdf("bkg3_pi0_shell", *m12, coefficients);
3711  CompositePdf *bkg3_pi0_transZ_total = new CompositePdf("bkg3_pi0_transZ_total", bkg3_pi0_transZ, bkg3_pi0_shell);
3712 
3713  comps.clear();
3714  comps.push_back(kzero_veto);
3715  comps.push_back(bkg3_pi0_transZ_total);
3716  // comps.push_back(bkg3_eff);
3717  // comps.push_back(bkg3_loZ);
3718  // ProdPdf* bkg3_pi0_mods = new ProdPdf("bkg3_pi0_mods", comps);
3719  // incsum3 = new IncoherentSumPdf("incsum3", m12, m13, eventNumber, special_pi0_decay, bkg3_pi0_mods);
3720 
3721  // Three spin-1 rho resonances to be added incoherently.
3722  DecayInfo3 incoherent_rhos;
3723  incoherent_rhos.motherMass = _mD0;
3724  incoherent_rhos.daug1Mass = piZeroMass;
3725  incoherent_rhos.daug2Mass = piPlusMass;
3726  incoherent_rhos.daug3Mass = piPlusMass;
3727  incoherent_rhos.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3728 
3729  ResonancePdf *bkg3_incRho0 = new Resonances::RBW(
3730  "bkg3_incRho0",
3731  Variable("bkg3_incRho0_amp", 0.00807 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3732  constantZero,
3733  Variable("bkg3_incRho0_mass", 0.800, 0.01, 0.6, 1.0),
3734  Variable("bkg3_incRho0_width", 0.15, 0.01, 0.10, 0.40),
3735  1, // These rhos are spin 1, being bad signal.
3736  PAIR_23);
3737  incoherent_rhos.resonances.push_back(bkg3_incRho0);
3738 
3739  ResonancePdf *bkg3_incRhoP = new Resonances::RBW(
3740  "bkg3_incRhoP",
3741  Variable("bkg3_incRhoP_amp", 0.01683 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3742  constantZero,
3743  Variable("bkg3_incRhoP_mass", 0.800, 0.01, 0.6, 1.0),
3744  Variable("bkg3_incRhoP_width", 0.15, 0.01, 0.10, 0.40),
3745  1,
3746  PAIR_12);
3747  incoherent_rhos.resonances.push_back(bkg3_incRhoP);
3748 
3749  ResonancePdf *bkg3_incRhoM = new Resonances::RBW(
3750  "bkg3_incRhoM",
3751  Variable("bkg3_incRhoM_amp", 0.01645 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3752  constantZero,
3753  Variable("bkg3_incRhoM_mass", 0.900, 0.01, 0.6, 1.0),
3754  Variable("bkg3_incRhoM_width", 0.35, 0.01, 0.10, 0.60),
3755  1,
3756  PAIR_13);
3757  incoherent_rhos.resonances.push_back(bkg3_incRhoM);
3758 
3759  comps.clear();
3760  comps.push_back(kzero_veto);
3761  // comps.push_back(bkg3_loZ);
3762  // comps.push_back(bkg3_eff);
3763  // ProdPdf* bkg3_rho_mods = new ProdPdf("bkg3_rho_mods", comps);
3764 
3765  // incsum4 = new IncoherentSumPdf("incsum4", m12, m13, eventNumber, incoherent_rhos, bkg3_rho_mods);
3766  // incsum4 = new IncoherentSumPdf("incsum4", m12, m13, eventNumber, incoherent_rhos, kzero_veto);
3767 
3768  weights.clear();
3769  weights.push_back(constantOne);
3770  // weights.push_back(constantOne);
3771  // weights.push_back(constantOne);
3772  comps.clear();
3773  comps.push_back(poly_x_veto);
3774  // comps.push_back(incsum3);
3775  // comps.push_back(incsum4);
3776 
3777  AddPdf *bkg3_dalitz = new AddPdf("bkg3_dalitz", weights, comps);
3778  bkg3_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3779  return bkg3_dalitz;
3780 }
std::vector< ResonancePdf * > resonances
Variable constantZero("constantZero", 0)
std::vector< Variable > weights
const fptype piPlusMass
Gaussian constructor.
Definition: ResonancePdf.h:110
Observable * m12
Variable constantOne("constantOne", 1)
const fptype _mD0
Observable * m13
std::vector< PdfBase * > comps
const fptype piZeroMass
Variable constantMinusOne("constantMinusOne", -1)
Relativistic Breit-Wigner.
Definition: ResonancePdf.h:68
GooPdf * kzero_veto
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88

◆ makeBackground4DalitzParam()

GooPdf* makeBackground4DalitzParam ( )

Definition at line 3782 of file pipipi0DPFit.cpp.

References _mD0, GooFit::PdfBase::addSpecialMask(), constantMinusOne, constantOne, constantZero, GooFit::DecayInfo3::daug1Mass, GooFit::DecayInfo3::daug2Mass, GooFit::DecayInfo3::daug3Mass, fixedRhoMass, fixedRhoWidth, massSum, GooFit::DecayInfo3::meson_radius, minDalitzZ, GooFit::DecayInfo3::motherMass, GooFit::PAIR_12, GooFit::PAIR_13, GooFit::PAIR_23, piPlusMass, piZeroMass, and GooFit::DecayInfo3::resonances.

Referenced by makeBkg4DalitzPdf().

3782  {
3783  vector<Variable> offsets;
3784  vector<Observable> observables;
3785  vector<Variable> coefficients;
3786  offsets.push_back(constantOne);
3787  offsets.push_back(constantOne);
3788 
3789  observables.push_back(*m12);
3790  observables.push_back(*m13);
3791  // Recurring factor 3 offsets division by total weight in AddPdf.
3792 
3793  double weightOffset = 3;
3794  // coefficients.push_back(new Variable("bkg4_x0y0", 1.0 * weightOffset, 0.01, 0.50*weightOffset,
3795  // 1.50*weightOffset));
3796  coefficients.push_back(Variable("bkg4_x0y0", 1.0 * weightOffset));
3797  coefficients.push_back(
3798  Variable("bkg4_x1y0", -0.18594 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3799  coefficients.push_back(
3800  Variable("bkg4_x2y0", 0.45459 * weightOffset, 0.01, 0.25 * weightOffset, 0.75 * weightOffset));
3801  coefficients.push_back(
3802  Variable("bkg4_x3y0", -0.20869 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3803  coefficients.push_back(
3804  Variable("bkg4_x0y1", -0.65061 * weightOffset, 0.01, -1.50 * weightOffset, 0.50 * weightOffset));
3805  coefficients.push_back(
3806  Variable("bkg4_x1y1", 0.11000 * weightOffset, 0.01, 0.00 * weightOffset, 0.50 * weightOffset));
3807  coefficients.push_back(
3808  Variable("bkg4_x2y1", 0.42009 * weightOffset, 0.01, 0.25 * weightOffset, 1.00 * weightOffset));
3809  coefficients.push_back(
3810  Variable("bkg4_x0y2", -0.06151 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3811  coefficients.push_back(
3812  Variable("bkg4_x1y2", 0.58508 * weightOffset, 0.01, 0.20 * weightOffset, 1.50 * weightOffset));
3813  coefficients.push_back(
3814  Variable("bkg4_x0y3", 0.54740 * weightOffset, 0.01, 0.20 * weightOffset, 1.50 * weightOffset));
3815 
3816  PolynomialPdf *poly = new PolynomialPdf("bkg4Pdf", observables, coefficients, offsets, 3);
3817 
3818  Variable bkg4_decZmin("bkg4_decZmin", 2.77576);
3819  Variable bkg4_conZmin("bkg4_conZmin", 0.23328);
3820  TrigThresholdPdf *bkg4_loZ
3821  = new TrigThresholdPdf("bkg4_loZ", *m12, *m13, minDalitzZ, bkg4_decZmin, bkg4_conZmin, massSum, false);
3822 
3823  comps.clear();
3824  comps.push_back(poly);
3825  comps.push_back(bkg4_loZ);
3826  comps.push_back(kzero_veto);
3827  // Separate PDF to avoid triggering numerical normalisation over all four observables.
3828  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3829 
3830  // One pipi bump.
3831  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3832  DecayInfo3 special_pipi_decay;
3833  special_pipi_decay.motherMass = _mD0;
3834  special_pipi_decay.daug1Mass = piZeroMass;
3835  special_pipi_decay.daug2Mass = piPlusMass;
3836  special_pipi_decay.daug3Mass = piPlusMass;
3837  special_pipi_decay.meson_radius = 1.5;
3838 
3839  ResonancePdf *bkg4_pipi_ref = new Resonances::Gauss("bkg4_pipi_ref",
3840  Variable("bkg4_pipi_ref_amp", 0.00147 * weightOffset),
3841  constantZero,
3842  Variable("bkg4_pipi_ref_mass", 1.32447),
3843  Variable("bkg4_pipi_ref_width", 0.04675),
3844  PAIR_23);
3845  special_pipi_decay.resonances.push_back(bkg4_pipi_ref);
3846 
3847  // Mikhail defines 'transverse Z' as y - x - (parameter).
3848  Variable bkg4_pipi_transZ_offset("bkg4_pipi_transZ_offset", -0.39877);
3849 
3850  offsets.clear();
3851  observables.clear();
3852  coefficients.clear();
3853  observables.push_back(*m12);
3854  observables.push_back(*m13);
3855  coefficients.push_back(bkg4_pipi_transZ_offset);
3856  coefficients.push_back(constantMinusOne);
3857  coefficients.push_back(constantOne);
3858  offsets.push_back(constantZero);
3859  offsets.push_back(constantZero);
3860  PolynomialPdf *bkg4_pipi_transZ = new PolynomialPdf("bkg4_pipi_transZ", observables, coefficients, offsets, 1);
3861 
3862  // Now we're going to take (1 - tz^2 * (parameter)) and multiply that into the pipi bump.
3863  Variable bkg4_pipi_transZ_quad("bkg4_pipi_transZ_quad", -0.25640);
3864  coefficients.clear();
3865  coefficients.push_back(constantOne);
3866  coefficients.push_back(constantZero);
3867  coefficients.push_back(bkg4_pipi_transZ_quad);
3868  // Notice the fake dependence of the polynomial on m12; in fact CompositePdf
3869  // will send it a fake event, we just have to supply a reasonable index.
3870  PolynomialPdf *bkg4_pipi_shell = new PolynomialPdf("bkg4_pipi_shell", *m12, coefficients);
3871  CompositePdf *bkg4_pipi_transZ_total
3872  = new CompositePdf("bkg4_pipi_transZ_total", bkg4_pipi_transZ, bkg4_pipi_shell);
3873 
3874  comps.clear();
3875  comps.push_back(kzero_veto);
3876  comps.push_back(bkg4_loZ);
3877  comps.push_back(bkg4_pipi_transZ_total);
3878 
3879  ProdPdf *bkg4_pipi_mods = new ProdPdf("bkg4_pipi_mods", comps);
3880  incsum5 = new IncoherentSumPdf("incsum5", *m12, *m13, *eventNumber, special_pipi_decay, bkg4_pipi_mods);
3881 
3882  // Three spin-0 rho resonances to be added incoherently.
3883  DecayInfo3 incoherent_rho0s;
3884  incoherent_rho0s.motherMass = _mD0;
3885  incoherent_rho0s.daug1Mass = piZeroMass;
3886  incoherent_rho0s.daug2Mass = piPlusMass;
3887  incoherent_rho0s.daug3Mass = piPlusMass;
3888  incoherent_rho0s.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3889 
3890  ResonancePdf *bkg4_incRho0 = new Resonances::RBW("bkg4_incRho0",
3891  Variable("bkg4_incRho0_amp", 0.00429 * weightOffset),
3892  constantZero,
3893  fixedRhoMass,
3894  fixedRhoWidth,
3895  0, // These rhos are spin 0.
3896  PAIR_23);
3897  incoherent_rho0s.resonances.push_back(bkg4_incRho0);
3898 
3899  ResonancePdf *bkg4_incRhoP = new Resonances::RBW("bkg4_incRhoP",
3900  Variable("bkg4_incRhoP_amp", 0.00705 * weightOffset),
3901  constantZero,
3902  fixedRhoMass,
3903  fixedRhoWidth,
3904  0,
3905  PAIR_12);
3906  incoherent_rho0s.resonances.push_back(bkg4_incRhoP);
3907 
3908  ResonancePdf *bkg4_incRhoM = new Resonances::RBW("bkg4_incRhoM",
3909  Variable("bkg4_incRhoM_amp", -0.00043 * weightOffset),
3910  constantZero,
3911  fixedRhoMass,
3912  fixedRhoWidth,
3913  0,
3914  PAIR_13);
3915  incoherent_rho0s.resonances.push_back(bkg4_incRhoM);
3916 
3917  comps.clear();
3918  comps.push_back(kzero_veto);
3919  comps.push_back(bkg4_loZ);
3920  ProdPdf *bkg4_incrho_mods = new ProdPdf("bkg4_incrho_mods", comps);
3921  incsum6 = new IncoherentSumPdf("incsum6", *m12, *m13, *eventNumber, incoherent_rho0s, bkg4_incrho_mods);
3922  // incsum6 = new IncoherentSumPdf("incsum6", m12, m13, eventNumber, incoherent_rho0s, kzero_veto);
3923 
3924  weights.clear();
3925  weights.push_back(constantOne);
3926  weights.push_back(constantOne);
3927  weights.push_back(constantOne);
3928  comps.clear();
3929  comps.push_back(poly_x_veto);
3930  comps.push_back(incsum5);
3931  comps.push_back(incsum6);
3932 
3933  AddPdf *bkg4_dalitz = new AddPdf("bkg4_dalitz", weights, comps);
3934  bkg4_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3935  return bkg4_dalitz;
3936 }
std::vector< ResonancePdf * > resonances
Variable constantZero("constantZero", 0)
Variable massSum("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
std::vector< Variable > weights
const fptype piPlusMass
Gaussian constructor.
Definition: ResonancePdf.h:110
Observable * m12
Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2)
Variable minDalitzZ("minDalitzZ", pow(piPlusMass+piPlusMass, 2))
Variable constantOne("constantOne", 1)
const fptype _mD0
Observable * m13
std::vector< PdfBase * > comps
const fptype piZeroMass
EventNumber * eventNumber
Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8)
Variable constantMinusOne("constantMinusOne", -1)
Relativistic Breit-Wigner.
Definition: ResonancePdf.h:68
IncoherentSumPdf * incsum5
GooPdf * kzero_veto
IncoherentSumPdf * incsum6
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88

◆ makeBackgroundHistogram()

SmoothHistogramPdf * makeBackgroundHistogram ( int  bkgnum,
std::string  overridename = "" 
)

Definition at line 3532 of file pipipi0DPFit.cpp.

References GooFit::BinnedDataSet::addEvent(), bkgHistRandSeed, GooFit::Application::get_filename(), GooFit::BinnedDataSet::getBinContent(), GooFit::BinnedDataSet::getNumBins(), GooFit::DataSet::getNumEvents(), m12, m13, readWrapper(), and GooFit::BinnedDataSet::setBinContent().

Referenced by makeBkg2DalitzPdf(), makeBkg3DalitzPdf(), makeBkg4DalitzPdf(), and writeBackgroundHistograms().

3532  {
3533  std::ifstream reader;
3534  sprintf(strbuffer, "./dataFiles/bkgDalitz_%i.txt", bkgnum);
3535 
3536  if(overridename != "")
3537  sprintf(strbuffer, "%s", overridename.c_str());
3538 
3539  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
3540  readWrapper(reader, fname);
3541  std::string buffer;
3542 
3543  while(!reader.eof()) {
3544  reader >> buffer;
3545 
3546  if(buffer == "====")
3547  break;
3548 
3549  std::cout << buffer;
3550  }
3551 
3552  obsweights.clear();
3553  obsweights.push_back(*m12);
3554  obsweights.push_back(*m13);
3555  BinnedDataSet *bkg_binned_data = new BinnedDataSet(obsweights);
3556 
3557  double dummy = 0;
3558 
3559  while(!reader.eof()) {
3560  reader >> dummy;
3561  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
3562  reader >> *m12; // Already swapped according to D* charge. m12 = m(pi+pi0)
3563  reader >> *m13;
3564 
3565  // Don't need the rest.
3566  for(int i = 0; i < 16; ++i)
3567  reader >> dummy;
3568 
3569  bkg_binned_data->addEvent();
3570 
3571  // std::cout << m12->getValue() << " " << m13->getValue() << std::endl;
3572  }
3573 
3574  std::cout << "Read " << bkg_binned_data->getNumEvents() << " events for background " << bkgnum << std::endl;
3575  sprintf(strbuffer, "bkg%i_dalitz_smoothing", bkgnum);
3576  Variable smoothing(strbuffer, 1);
3577 
3578  if((-1 != bkgHistRandSeed) && ((3 == bkgnum) || (4 == bkgnum))) {
3579  std::cout << "Shuffling background " << bkgnum << " histogram with random seed " << bkgHistRandSeed
3580  << std::endl;
3581  TRandom donram(bkgHistRandSeed);
3582 
3583  for(unsigned int bin = 0; bin < bkg_binned_data->getNumBins(); ++bin) {
3584  double events = bkg_binned_data->getBinContent(bin);
3585 
3586  if(1 > events)
3587  continue;
3588 
3589  double newEvents = -1;
3590 
3591  while(0 > newEvents)
3592  newEvents = donram.Gaus(events, sqrt(events));
3593 
3594  bkg_binned_data->setBinContent(bin, newEvents);
3595  }
3596  }
3597 
3598  sprintf(strbuffer, "bkg%i_dalitz", bkgnum);
3599  SmoothHistogramPdf *bkg_dalitz = new SmoothHistogramPdf(strbuffer, bkg_binned_data, smoothing);
3600  // bkg_dalitz->setDebugMask(1);
3601  return bkg_dalitz;
3602 }
size_t getNumEvents() const
Definition: DataSet.h:47
fptype getBinContent(size_t bin) const
Get the content of a bin.
Definition: BinnedDataSet.h:26
void addEvent() override
char strbuffer[1000]
Observable * m12
void setBinContent(unsigned int bin, fptype value)
Definition: BinnedDataSet.h:45
Observable * m13
size_t getNumBins() const
int bkgHistRandSeed
bool readWrapper(std::ifstream &reader, std::string fname=strbuffer)
GooFit::Application * app_ptr
std::string get_filename(const std::string &input_str, std::string base="") const
std::vector< Observable > obsweights

◆ makeBkg2_sigma()

GooPdf* makeBkg2_sigma ( )

Definition at line 3021 of file pipipi0DPFit.cpp.

Referenced by runBackgroundSigmaFit().

3021  {
3022  // Variable* bkg2_sigma_js_meana = new Variable("bkg2_sigma_js_meana", 0.01, 0.01, -0.30, 0.30);
3023  // Variable* bkg2_sigma_js_sigma = new Variable("bkg2_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3024  // Variable* bkg2_sigma_js_gamma = new Variable("bkg2_sigma_js_gamma",-5.00, 0.10, -30, 0);
3025  // Variable* bkg2_sigma_js_delta = new Variable("bkg2_sigma_js_delta", 1.49, 0.01, 0.50, 5.00);
3026  // Variable* bkg2_sigma_frac_jsu = new Variable("bkg2_sigma_frac_jsu", 0.85, 0.01, 0.01, 1.00);
3027  // Variable* bkg2_sigma_frac_ga1 = new Variable("bkg2_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3028 
3029  Variable bkg2_sigma_num_jsu("bkg2_sigma_num_jsu", 9100, 200, 1000, 15000);
3030  Variable bkg2_sigma_num_ga1("bkg2_sigma_num_ga1", 2400, 200, 500, 7000);
3031  Variable bkg2_sigma_num_ga2("bkg2_sigma_num_ga2", 2900, 200, 500, 7000);
3032 
3033  Variable bkg2_sigma_g1_meana("bkg2_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3034  Variable bkg2_sigma_g1_sigma("bkg2_sigma_g1_sigma", 0.30, 0.01, 0.05, 0.55);
3035  Variable bkg2_sigma_g2_meana("bkg2_sigma_g2_meana", 0.80, 0.01, 0.01, 1.50);
3036  Variable bkg2_sigma_g2_sigma("bkg2_sigma_g2_sigma", 0.90, 0.01, 0.01, 2.75);
3037  // JohnsonSUPdf* bkg2_sigma_js = new JohnsonSUPdf("bkg2_sigma_js", sigma, bkg2_sigma_js_meana, bkg2_sigma_js_sigma,
3038  // bkg2_sigma_js_gamma, bkg2_sigma_js_delta);
3039 
3040  Variable bkg2_sigma_js_meana("bkg2_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3041  Variable bkg2_sigma_js_sigma("bkg2_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3042  Variable bkg2_sigma_js_gamma("bkg2_sigma_js_gamma", 2.00, 0.10, 0, 10);
3043  Variable bkg2_sigma_js_delta("bkg2_sigma_js_delta", 2);
3044  CrystalBallPdf *bkg2_sigma_js = new CrystalBallPdf(
3045  "bkg2_sigma_js", *sigma, bkg2_sigma_js_meana, bkg2_sigma_js_sigma, bkg2_sigma_js_gamma, bkg2_sigma_js_delta);
3046 
3047  GaussianPdf *bkg2_sigma_g1 = new GaussianPdf("bkg2_sigma_g1", *sigma, bkg2_sigma_g1_meana, bkg2_sigma_g1_sigma);
3048  GaussianPdf *bkg2_sigma_g2 = new GaussianPdf("bkg2_sigma_g2", *sigma, bkg2_sigma_g2_meana, bkg2_sigma_g2_sigma);
3049 
3050  weights.clear();
3051  weights.push_back(bkg2_sigma_num_jsu);
3052  weights.push_back(bkg2_sigma_num_ga1);
3053  weights.push_back(bkg2_sigma_num_ga2);
3054  comps.clear();
3055  comps.push_back(bkg2_sigma_js);
3056  comps.push_back(bkg2_sigma_g1);
3057  comps.push_back(bkg2_sigma_g2);
3058 
3059  GooPdf *ret = new AddPdf("bkg2_jsugg", weights, comps);
3060  return ret;
3061  // return bkg2_sigma_js;
3062 }
std::vector< Variable > weights
Observable * sigma
std::vector< PdfBase * > comps

◆ makeBkg2DalitzPdf()

GooPdf* makeBkg2DalitzPdf ( bool  fixem = true)

Definition at line 3276 of file pipipi0DPFit.cpp.

References _mD0, GooFit::PdfBase::addSpecialMask(), bkg2Model, constantOne, constantZero, GooFit::DecayInfo3::daug1Mass, GooFit::DecayInfo3::daug2Mass, GooFit::DecayInfo3::daug3Mass, fixedRhoMass, fixedRhoWidth, gaussBkgTime, GooFit::Application::get_filename(), Histogram, makeBackgroundHistogram(), makeBkg_sigma_strips(), makeExpGausTimePdf(), makeGaussianTimePdf(), makeKzeroVeto(), massSum, GooFit::DecayInfo3::meson_radius, minDalitzZ, GooFit::DecayInfo3::motherMass, GooFit::PAIR_12, GooFit::PAIR_13, GooFit::PAIR_23, Parameter, piPlusMass, piZeroMass, GooFit::DecayInfo3::resonances, GooFit::GooPdf::setParameterConstantness(), and Sideband.

Referenced by runBackgroundDalitzFit(), and runCanonicalFit().

3276  {
3277  if(!kzero_veto)
3278  makeKzeroVeto();
3279 
3280  GooPdf *bkg2_dalitz = nullptr;
3281 
3282  if(Parameter == bkg2Model) {
3283  comps.clear();
3284 
3285  vector<Variable> offsets;
3286  vector<Observable> observables;
3287  vector<Variable> coefficients;
3288  offsets.push_back(constantOne);
3289  offsets.push_back(constantOne);
3290  observables.push_back(*m12);
3291  observables.push_back(*m13);
3292  double weightOffset = 3;
3293  // Recurring factor 3 offsets division by total weight in AddPdf.
3294  coefficients.push_back(Variable("bkg2_x0y0", 1.0 * weightOffset));
3295  coefficients.push_back(
3296  Variable("bkg2_x1y0", 0.13184 * weightOffset, 0.01, 0.01 * weightOffset, 0.18 * weightOffset));
3297  coefficients.push_back(
3298  Variable("bkg2_x2y0", 0.02062 * weightOffset, 0.01, 0.00 * weightOffset, 0.17 * weightOffset));
3299  coefficients.push_back(
3300  Variable("bkg2_x3y0", 0.04688 * weightOffset, 0.01, 0.00 * weightOffset, 0.08 * weightOffset));
3301  coefficients.push_back(
3302  Variable("bkg2_x0y1", -0.02568 * weightOffset, 0.01, -0.15 * weightOffset, 0.04 * weightOffset));
3303  coefficients.push_back(
3304  Variable("bkg2_x1y1", 0.06805 * weightOffset, 0.01, 0.02 * weightOffset, 0.10 * weightOffset));
3305  coefficients.push_back(
3306  Variable("bkg2_x2y1", 0.38557 * weightOffset, 0.01, 0.30 * weightOffset, 0.50 * weightOffset));
3307  coefficients.push_back(
3308  Variable("bkg2_x0y2", 0.11252 * weightOffset, 0.01, 0.05 * weightOffset, 0.20 * weightOffset));
3309  coefficients.push_back(
3310  Variable("bkg2_x1y2", 0.24896 * weightOffset, 0.01, 0.20 * weightOffset, 0.50 * weightOffset));
3311  coefficients.push_back(
3312  Variable("bkg2_x0y3", 0.05605 * weightOffset, 0.01, -0.05 * weightOffset, 0.15 * weightOffset));
3313 
3314  PolynomialPdf *poly = new PolynomialPdf("bkg2Pdf", observables, coefficients, offsets, 3);
3315 
3316  Variable bkg2_decZmin("bkg2_decZmin", 3.30411);
3317  Variable bkg2_conZmin("bkg2_conZmin", 0.29909);
3318  TrigThresholdPdf *bkg2_loZ
3319  = new TrigThresholdPdf("bkg2_loZ", *m12, *m13, minDalitzZ, bkg2_decZmin, bkg2_conZmin, massSum, false);
3320  // bkg2_loZ->setDebugMask(1);
3321 
3322  comps.clear();
3323  comps.push_back(poly);
3324  comps.push_back(bkg2_loZ);
3325  comps.push_back(kzero_veto);
3326  // Separate PDF to avoid triggering numerical normalisation over all four observables.
3327  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3328  // poly_x_veto->setDebugMask(1);
3329 
3330  // One omega->pipipi0 reflection.
3331  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3332  DecayInfo3 special_rho_decay;
3333  special_rho_decay.motherMass = _mD0;
3334  special_rho_decay.daug1Mass = piZeroMass;
3335  special_rho_decay.daug2Mass = piPlusMass;
3336  special_rho_decay.daug3Mass = piPlusMass;
3337  special_rho_decay.meson_radius = 1.5;
3338 
3339  ResonancePdf *bkg2_rho_ref = new Resonances::Gauss(
3340  "bkg2_rho_ref",
3341  Variable("bkg2_rho_ref_amp", 0.00896 * weightOffset, 0.001, 0, 0.015 * weightOffset),
3342  constantZero,
3343  Variable("bkg2_rho_ref_mass", 0.53172),
3344  Variable("bkg2_rho_ref_width", 0.06426),
3345  PAIR_13);
3346  special_rho_decay.resonances.push_back(bkg2_rho_ref);
3347 
3348  Variable bkg2_rho_poly_offset("bkg2_rho_poly_offset", 1.64254);
3349  Variable bkg2_rho_poly_linear("bkg2_rho_poly_linear", 0);
3350  Variable bkg2_rho_poly_second("bkg2_rho_poly_second", -0.48166);
3351  weights.clear();
3352  weights.push_back(constantOne);
3353  weights.push_back(bkg2_rho_poly_linear);
3354  weights.push_back(bkg2_rho_poly_second);
3355  PolynomialPdf *bkg2_rho_poly = new PolynomialPdf("bkg2_rho_poly", *m12, weights, bkg2_rho_poly_offset);
3356  comps.clear();
3357  comps.push_back(kzero_veto);
3358  comps.push_back(bkg2_rho_poly);
3359  comps.push_back(bkg2_loZ);
3360  ProdPdf *bkg2_rho_mods = new ProdPdf("bkg2_rho_mods", comps);
3361  incsum1 = new IncoherentSumPdf("incsum1", *m12, *m13, *eventNumber, special_rho_decay, bkg2_rho_mods);
3362 
3363  // Three spin-0 rho resonances to be added incoherently.
3364  DecayInfo3 incoherent_rho0s;
3365  incoherent_rho0s.motherMass = _mD0;
3366  incoherent_rho0s.daug1Mass = piZeroMass;
3367  incoherent_rho0s.daug2Mass = piPlusMass;
3368  incoherent_rho0s.daug3Mass = piPlusMass;
3369  incoherent_rho0s.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3370 
3371  ResonancePdf *bkg2_incRho0 = new Resonances::RBW(
3372  "bkg2_incRho0",
3373  Variable("bkg2_incRho0_amp", 0.00304 * weightOffset, 0.001, 0.0, 0.006 * weightOffset),
3374  constantZero,
3375  fixedRhoMass,
3376  fixedRhoWidth,
3377  0, // Incoherent rho has effective spin 0.
3378  PAIR_23);
3379  incoherent_rho0s.resonances.push_back(bkg2_incRho0);
3380 
3381  ResonancePdf *bkg2_incRhoP = new Resonances::RBW(
3382  "bkg2_incRhoP",
3383  Variable("bkg2_incRhoP_amp", 0.00586 * weightOffset, 0.001, 0.0, 0.012 * weightOffset),
3384  constantZero,
3385  fixedRhoMass,
3386  fixedRhoWidth,
3387  0,
3388  PAIR_12);
3389  incoherent_rho0s.resonances.push_back(bkg2_incRhoP);
3390 
3391  ResonancePdf *bkg2_incRhoM = new Resonances::RBW(
3392  "bkg2_incRhoM",
3393  Variable("bkg2_incRhoM_amp", 0.00635 * weightOffset, 0.001, 0.0, 0.015 * weightOffset),
3394  constantZero,
3395  fixedRhoMass,
3396  fixedRhoWidth,
3397  0,
3398  PAIR_13);
3399  incoherent_rho0s.resonances.push_back(bkg2_incRhoM);
3400 
3401  comps.clear();
3402  comps.push_back(kzero_veto);
3403  comps.push_back(bkg2_loZ);
3404  ProdPdf *bkg2_rho_mods2 = new ProdPdf("bkg2_rho_mods2", comps);
3405 
3406  incsum2 = new IncoherentSumPdf("incsum2", *m12, *m13, *eventNumber, incoherent_rho0s, bkg2_rho_mods2);
3407 
3408  weights.clear();
3409  weights.push_back(constantOne);
3410  weights.push_back(constantOne);
3411  weights.push_back(constantOne);
3412  comps.clear();
3413  comps.push_back(poly_x_veto);
3414  comps.push_back(incsum1);
3415  comps.push_back(incsum2);
3416 
3417  bkg2_dalitz = new AddPdf("bkg2_dalitz", weights, comps);
3418  bkg2_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3419  // bkg2_dalitz->setDebugMask(1);
3420 
3421  } else if(Histogram == bkg2Model) {
3422  bkg2_dalitz = makeBackgroundHistogram(2);
3423  } else if(Sideband == bkg2Model) {
3424  comps.clear();
3425  // Originally had massD0 sidebands separated into deltaM high and low,
3426  // but these distributions were extremely similar - collapsed them
3427  // into just massD0 sidebands.
3428  std::string fname = app_ptr->get_filename("dataFiles/sideband1.txt", "examples/pipipi0DPFit");
3429  comps.push_back(makeBackgroundHistogram(101, fname));
3430  fname = app_ptr->get_filename("dataFiles/sideband2.txt", "examples/pipipi0DPFit");
3431  comps.push_back(makeBackgroundHistogram(102, fname));
3432  // comps.push_back(makeBackgroundHistogram(103, "./dataFiles/sideband3.txt"));
3433  // comps.push_back(makeBackgroundHistogram(104, "./dataFiles/sideband4.txt"));
3434  weights.clear();
3435  weights.push_back(Variable("sband1Weight", 300000, 1000, 100, 750000));
3436  weights.push_back(Variable("sband2Weight", 100000, 1000, 100, 500000));
3437  // weights.push_back(new Variable("sband3Weight", 150000, 1000, 100, 500000));
3438  // weights.push_back(new Variable("sband4Weight", 50000, 1000, 100, 500000));
3439  bkg2_dalitz = new AddPdf("bkg2_dalitz", weights, comps);
3440  } else {
3441  // This cannot possibly happen, and if it does something is wrong. Panic!
3442  assert(Sideband == bkg2Model);
3443  }
3444 
3445  GooPdf *bkg2_dtime = nullptr;
3446 
3447  if(gaussBkgTime)
3448  bkg2_dtime = makeGaussianTimePdf(2);
3449  else
3450  bkg2_dtime = makeExpGausTimePdf(2);
3451 
3452  // bkg2_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
3453  // Separate sigma_t distribution
3454  // bkg2_jsugg = makeBkg2_sigma();
3456  bkg2_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
3457 
3458  // Finally create overall product.
3459  comps.clear();
3460  // poly->setDebugMask(1);
3461  // bkg2_dalitz->setDebugMask(1);
3462  // incsum2->setDebugMask(1);
3463  comps.push_back(bkg2_dalitz);
3464  comps.push_back(bkg2_dtime);
3465  comps.push_back(bkg2_jsugg);
3466 
3467  GooPdf *ret = new ProdPdf("bkg2_total", comps);
3468 
3469  if(fixem)
3470  ret->setParameterConstantness(true);
3471 
3472  return ret;
3473 }
IncoherentSumPdf * incsum2
std::vector< ResonancePdf * > resonances
__host__ void setParameterConstantness(bool constant=true)
Variable constantZero("constantZero", 0)
Variable massSum("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
std::vector< Variable > weights
Bkg2Model bkg2Model
GooPdf * makeBkg_sigma_strips(int bkgnum)
const fptype piPlusMass
Gaussian constructor.
Definition: ResonancePdf.h:110
Observable * m12
Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2)
Variable minDalitzZ("minDalitzZ", pow(piPlusMass+piPlusMass, 2))
Variable constantOne("constantOne", 1)
const fptype _mD0
Observable * m13
std::vector< PdfBase * > comps
const fptype piZeroMass
bool gaussBkgTime
EventNumber * eventNumber
Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8)
GooPdf * makeKzeroVeto()
GooPdf * makeGaussianTimePdf(int bkg)
GooPdf * makeExpGausTimePdf(int bkg)
Relativistic Breit-Wigner.
Definition: ResonancePdf.h:68
IncoherentSumPdf * incsum1
GooPdf * bkg2_jsugg
GooPdf * kzero_veto
GooFit::Application * app_ptr
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88
std::string get_filename(const std::string &input_str, std::string base="") const
SmoothHistogramPdf * makeBackgroundHistogram(int bkgnum, std::string overridename="")

◆ makeBkg3_sigma()

GooPdf* makeBkg3_sigma ( )

Definition at line 3106 of file pipipi0DPFit.cpp.

3106  {
3107  // Variable* bkg3_sigma_js_meana = new Variable("bkg3_sigma_js_meana", 0.05, 0.01, -0.30, 0.30);
3108  // Variable* bkg3_sigma_js_sigma = new Variable("bkg3_sigma_js_sigma", 0.013, 0.01, 0, 0.2);
3109  // Variable* bkg3_sigma_js_gamma = new Variable("bkg3_sigma_js_gamma",-6.00, 1.00, -30, 0);
3110  // Variable* bkg3_sigma_js_delta = new Variable("bkg3_sigma_js_delta", 1.99, 0.10, 0.50, 5.00);
3111  Variable bkg3_sigma_frac_jsu("bkg3_sigma_frac_jsu", 0.50, 0.01, 0.01, 1.00);
3112  Variable bkg3_sigma_frac_ga1("bkg3_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3113 
3114  // Variable* bkg3_sigma_num_jsu = new Variable("bkg3_sigma_num_jsu", 11000, 200, 1000, 35000);
3115  // Variable* bkg3_sigma_num_ga1 = new Variable("bkg3_sigma_num_ga1", 9400, 200, 500, 10000);
3116  // Variable* bkg3_sigma_num_ga2 = new Variable("bkg3_sigma_num_ga2", 3900, 200, 500, 17000);
3117  // Variable* bkg3_sigma_num_ga3 = new Variable("bkg3_sigma_num_ga3", 3900, 200, 500, 7000);
3118 
3119  Variable bkg3_sigma_g1_meana("bkg3_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3120  Variable bkg3_sigma_g1_sigma("bkg3_sigma_g1_sigma", 0.10, 0.01, 0.01, 0.55);
3121  Variable bkg3_sigma_g2_meana("bkg3_sigma_g2_meana", 0.20, 0.01, 0.01, 1.50);
3122  Variable bkg3_sigma_g2_sigma("bkg3_sigma_g2_sigma", 0.10, 0.01, 0.001, 0.15);
3123  Variable bkg3_sigma_g2_gamma("bkg3_sigma_g2_gamma", -2.00, 1.00, -10, 10);
3124  Variable bkg3_sigma_g2_delta("bkg3_sigma_g2_delta", 2, 0.10, 0.50, 5.00);
3125  // Variable* bkg3_sigma_g3_meana = new Variable("bkg3_sigma_g3_meana", 0.20, 0.01, 0.01, 1.50);
3126  // Variable* bkg3_sigma_g3_sigma = new Variable("bkg3_sigma_g3_sigma", 0.20, 0.01, 0.01, 0.75);
3127  // Variable* bkg3_sigma_g2_width = new Variable("bkg3_sigma_g2_width", 0.10, 0.01, 0.01, 0.75);
3128  // JohnsonSUPdf* bkg3_sigma_js = new JohnsonSUPdf("bkg3_sigma_js", sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma,
3129  // bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3130 
3131  Variable bkg3_sigma_js_meana("bkg3_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3132  Variable bkg3_sigma_js_sigma("bkg3_sigma_js_sigma", 0.09, 0.01, 0, 0.40);
3133  Variable bkg3_sigma_js_gamma("bkg3_sigma_js_gamma", 2.00, 0.10, 0, 10);
3134  Variable bkg3_sigma_js_delta("bkg3_sigma_js_delta", 2);
3135  CrystalBallPdf *bkg3_sigma_js = new CrystalBallPdf(
3136  "bkg3_sigma_js", *sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma, bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3137  // JohnsonSUPdf* bkg3_sigma_js = new JohnsonSUPdf("bkg3_sigma_js", sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma,
3138  // bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3139 
3140  GaussianPdf *bkg3_sigma_g1 = new GaussianPdf("bkg3_sigma_g1", *sigma, bkg3_sigma_g1_meana, bkg3_sigma_g1_sigma);
3141  // GaussianPdf* bkg3_sigma_g2 = new GaussianPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma);
3142  // CrystalBallPdf* bkg3_sigma_g2 = new CrystalBallPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana,
3143  // bkg3_sigma_g2_sigma, bkg3_sigma_g2_gamma, bkg3_sigma_g2_delta);
3144  JohnsonSUPdf *bkg3_sigma_g2 = new JohnsonSUPdf(
3145  "bkg3_sigma_g2", *sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma, bkg3_sigma_g2_gamma, bkg3_sigma_g2_delta);
3146  // GaussianPdf* bkg3_sigma_g3 = new GaussianPdf("bkg3_sigma_g3", sigma, bkg3_sigma_g3_meana, bkg3_sigma_g3_sigma);
3147  // VoigtianPdf* bkg3_sigma_g2 = new VoigtianPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma,
3148  // bkg3_sigma_g2_width);
3149 
3150  weights.clear();
3151  // weights.push_back(bkg3_sigma_num_jsu);
3152  // weights.push_back(bkg3_sigma_num_ga1);
3153  // weights.push_back(bkg3_sigma_num_ga2);
3154  // weights.push_back(bkg3_sigma_num_ga3);
3155  weights.push_back(bkg3_sigma_frac_jsu);
3156  weights.push_back(bkg3_sigma_frac_ga1);
3157  comps.clear();
3158  comps.push_back(bkg3_sigma_js);
3159  comps.push_back(bkg3_sigma_g1);
3160  comps.push_back(bkg3_sigma_g2);
3161  // comps.push_back(bkg3_sigma_g3);
3162 
3163  GooPdf *ret = new AddPdf("bkg3_jsugg", weights, comps);
3164  return ret;
3165 }
std::vector< Variable > weights
Observable * sigma
std::vector< PdfBase * > comps

◆ makeBkg3DalitzPdf()

GooPdf* makeBkg3DalitzPdf ( bool  fixem = true)

Definition at line 3938 of file pipipi0DPFit.cpp.

References GooFit::PdfBase::addSpecialMask(), gaussBkgTime, makeBackground3DalitzParam(), makeBackgroundHistogram(), makeBkg_sigma_strips(), makeExpGausTimePdf(), makeGaussianTimePdf(), makeKzeroVeto(), notUseBackground3Hist, and GooFit::GooPdf::setParameterConstantness().

Referenced by runBackgroundDalitzFit(), and runCanonicalFit().

3938  {
3939  if(!kzero_veto)
3940  makeKzeroVeto();
3941 
3942  comps.clear();
3943  weights.clear();
3944 
3945  // Using a histogram for Dalitz description. Notice: The purpose of this
3946  // fit is to get a description for use in fitting actual data. When fitting
3947  // data I can just use the background 3 Monte Carlo histogram. When fitting
3948  // the MC, it doesn't matter what I do, because I'm going to be using the
3949  // histogram. So I load up all the MC data and use it in the histogram, either
3950  // way.
3951 
3952  GooPdf *bkg3_dalitz = nullptr;
3953 
3955  bkg3_dalitz = makeBackgroundHistogram(3);
3956  else
3957  bkg3_dalitz = makeBackground3DalitzParam();
3958 
3959  // bkg3_dalitz->setDebugMask(1);
3960 
3961  GooPdf *bkg3_dtime = nullptr;
3962 
3963  if(gaussBkgTime)
3964  bkg3_dtime = makeGaussianTimePdf(3);
3965  else
3966  bkg3_dtime = makeExpGausTimePdf(3);
3967 
3968  // bkg3_jsugg = makeBkg3_sigma();
3969  // bkg3_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
3971  bkg3_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
3972  // Otherwise ProdPdf tries to use the default overall integration,
3973  // because bkg3_jsugg depends on m12, m13 due to the striping, and that has
3974  // disastrous results for bkg3_dalitz. Note that this doesn't, actually,
3975  // contradict the ForceCommonNorm above, even though it looks like it should,
3976  // because CommonNorm applies to the AddPdf while SeparateNorm
3977  // applies to the ProdPdf.
3978 
3979  comps.clear();
3980  comps.push_back(bkg3_dalitz);
3981  // bkg3_dalitz->setDebugMask(1);
3982  // incsum3->setDebugMask(1);
3983  comps.push_back(bkg3_dtime);
3984  // bkg3_dtime->setDebugMask(1);
3985  // comps.push_back(bkg3_jsugg);
3986  // sig0_jsugg->setDebugMask(1);
3987 
3988  GooPdf *ret = new ProdPdf("bkg3_total", comps);
3989 
3990  if(fixem)
3991  ret->setParameterConstantness(true);
3992 
3993  return ret;
3994 }
__host__ void setParameterConstantness(bool constant=true)
bool notUseBackground3Hist
std::vector< Variable > weights
GooPdf * makeBkg_sigma_strips(int bkgnum)
GooPdf * bkg3_jsugg
GooPdf * makeBackground3DalitzParam()
std::vector< PdfBase * > comps
bool gaussBkgTime
GooPdf * makeKzeroVeto()
GooPdf * makeGaussianTimePdf(int bkg)
GooPdf * makeExpGausTimePdf(int bkg)
GooPdf * kzero_veto
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88
SmoothHistogramPdf * makeBackgroundHistogram(int bkgnum, std::string overridename="")

◆ makeBkg3Eff()

GooPdf* makeBkg3Eff ( )

Definition at line 3475 of file pipipi0DPFit.cpp.

References GooFit::BinnedDataSet::addEvent(), GooFit::Application::get_filename(), GooFit::Observable::getNumBins(), m12, m13, readWrapper(), and GooFit::Observable::setNumBins().

3475  {
3476  // Smoothed histogram from flat-file data.
3477  // Only 4500 events, so use large bins.
3478 
3479  obsweights.clear();
3480  obsweights.push_back(*m12);
3481  obsweights.push_back(*m13);
3482 
3483  int m12bins = m12->getNumBins();
3484  int m13bins = m13->getNumBins();
3485 
3486  m12->setNumBins(30);
3487  m13->setNumBins(30);
3488  BinnedDataSet *bkg3_eff_data = new BinnedDataSet(obsweights);
3489  std::ifstream reader;
3490  std::string fname = app_ptr->get_filename("dataFiles/efficiency_bkg3_flat.txt", "examples/pipipi0DPFit");
3491  readWrapper(reader, fname);
3492  std::string buffer;
3493 
3494  while(!reader.eof()) {
3495  reader >> buffer;
3496 
3497  if(buffer == "====")
3498  break;
3499 
3500  std::cout << buffer;
3501  }
3502 
3503  double dummy = 0;
3504 
3505  while(!reader.eof()) {
3506  reader >> dummy;
3507 
3508  if(reader.eof())
3509  break;
3510 
3511  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
3512  reader >> *m12; // Already swapped according to D* charge
3513  reader >> *m13;
3514 
3515  // Everything else is irrelevant for this purpose!
3516  for(int i = 0; i < 16; ++i)
3517  reader >> dummy;
3518 
3519  bkg3_eff_data->addEvent();
3520  }
3521 
3522  Variable bkg3_eff_smoothing("bkg3_eff_smoothing", 1.0, 0.1, 0, 1.25);
3523  // Variable* bkg3_eff_smoothing = new Variable("bkg3_eff_smoothing", 1.0);
3524  SmoothHistogramPdf *ret = new SmoothHistogramPdf("bkg3_efficiency", bkg3_eff_data, bkg3_eff_smoothing);
3525 
3526  m12->setNumBins(m12bins);
3527  m13->setNumBins(m13bins);
3528 
3529  return ret;
3530 }
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
void addEvent() override
Observable * m12
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
Observable * m13
bool readWrapper(std::ifstream &reader, std::string fname=strbuffer)
GooFit::Application * app_ptr
std::string get_filename(const std::string &input_str, std::string base="") const
std::vector< Observable > obsweights

◆ makeBkg4_sigma()

GooPdf* makeBkg4_sigma ( )

Definition at line 3064 of file pipipi0DPFit.cpp.

3064  {
3065  // Variable* bkg4_sigma_js_meana = new Variable("bkg4_sigma_js_meana", 0.01, 0.01, -0.30, 0.30);
3066  // Variable* bkg4_sigma_js_sigma = new Variable("bkg4_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3067  // Variable* bkg4_sigma_js_gamma = new Variable("bkg4_sigma_js_gamma",-5.00, 0.10, -30, 0);
3068  // Variable* bkg4_sigma_js_delta = new Variable("bkg4_sigma_js_delta", 1.49, 0.01, 0.50, 5.00);
3069  // Variable* bkg4_sigma_frac_jsu = new Variable("bkg4_sigma_frac_jsu", 0.85, 0.01, 0.01, 1.00);
3070  // Variable* bkg4_sigma_frac_ga1 = new Variable("bkg4_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3071 
3072  Variable bkg4_sigma_num_jsu("bkg4_sigma_num_jsu", 9100, 200, 1000, 15000);
3073  Variable bkg4_sigma_num_ga1("bkg4_sigma_num_ga1", 2400, 200, 500, 7000);
3074  Variable bkg4_sigma_num_ga2("bkg4_sigma_num_ga2", 2900, 200, 500, 7000);
3075 
3076  Variable bkg4_sigma_g1_meana("bkg4_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3077  Variable bkg4_sigma_g1_sigma("bkg4_sigma_g1_sigma", 0.30, 0.01, 0.05, 0.55);
3078  Variable bkg4_sigma_g2_meana("bkg4_sigma_g2_meana", 0.80, 0.01, 0.01, 1.50);
3079  Variable bkg4_sigma_g2_sigma("bkg4_sigma_g2_sigma", 0.90, 0.01, 0.01, 2.75);
3080  // JohnsonSUPdf* bkg4_sigma_js = new JohnsonSUPdf("bkg4_sigma_js", sigma, bkg4_sigma_js_meana, bkg4_sigma_js_sigma,
3081  // bkg4_sigma_js_gamma, bkg4_sigma_js_delta);
3082 
3083  Variable bkg4_sigma_js_meana("bkg4_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3084  Variable bkg4_sigma_js_sigma("bkg4_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3085  Variable bkg4_sigma_js_gamma("bkg4_sigma_js_gamma", 2.00, 0.10, 0, 10);
3086  Variable bkg4_sigma_js_delta("bkg4_sigma_js_delta", 2);
3087  CrystalBallPdf *bkg4_sigma_js = new CrystalBallPdf(
3088  "bkg4_sigma_js", *sigma, bkg4_sigma_js_meana, bkg4_sigma_js_sigma, bkg4_sigma_js_gamma, bkg4_sigma_js_delta);
3089 
3090  GaussianPdf *bkg4_sigma_g1 = new GaussianPdf("bkg4_sigma_g1", *sigma, bkg4_sigma_g1_meana, bkg4_sigma_g1_sigma);
3091  GaussianPdf *bkg4_sigma_g2 = new GaussianPdf("bkg4_sigma_g2", *sigma, bkg4_sigma_g2_meana, bkg4_sigma_g2_sigma);
3092 
3093  weights.clear();
3094  weights.push_back(bkg4_sigma_num_jsu);
3095  weights.push_back(bkg4_sigma_num_ga1);
3096  weights.push_back(bkg4_sigma_num_ga2);
3097  comps.clear();
3098  comps.push_back(bkg4_sigma_js);
3099  comps.push_back(bkg4_sigma_g1);
3100  comps.push_back(bkg4_sigma_g2);
3101 
3102  GooPdf *ret = new AddPdf("bkg4_jsugg", weights, comps);
3103  return ret;
3104 }
std::vector< Variable > weights
Observable * sigma
std::vector< PdfBase * > comps

◆ makeBkg4DalitzPdf()

GooPdf* makeBkg4DalitzPdf ( bool  fixem = true)

Definition at line 3996 of file pipipi0DPFit.cpp.

References GooFit::PdfBase::addSpecialMask(), gaussBkgTime, makeBackground4DalitzParam(), makeBackgroundHistogram(), makeBkg_sigma_strips(), makeExpGausTimePdf(), makeGaussianTimePdf(), makeKzeroVeto(), notUseBackground4Hist, and GooFit::GooPdf::setParameterConstantness().

Referenced by runBackgroundDalitzFit(), and runCanonicalFit().

3996  {
3997  if(!kzero_veto)
3998  makeKzeroVeto();
3999 
4000  comps.clear();
4001  weights.clear();
4002 
4003  GooPdf *bkg4_dalitz = nullptr;
4004 
4006  bkg4_dalitz = makeBackgroundHistogram(4);
4007  else
4008  bkg4_dalitz = makeBackground4DalitzParam();
4009 
4010  GooPdf *bkg4_dtime = nullptr;
4011 
4012  if(gaussBkgTime)
4013  bkg4_dtime = makeGaussianTimePdf(4);
4014  else
4015  bkg4_dtime = makeExpGausTimePdf(4);
4016 
4017  // Separate sigma_t distribution
4018  // bkg4_jsugg = makeBkg4_sigma();
4019  // bkg4_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
4021  bkg4_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm); // See comments to bkg3_jsugg.
4022 
4023  comps.clear();
4024  comps.push_back(bkg4_dalitz);
4025  comps.push_back(bkg4_dtime);
4026  // comps.push_back(bkg4_jsugg);
4027 
4028  ProdPdf *ret = new ProdPdf("bkg4_total", comps);
4029 
4030  if(fixem)
4031  ret->setParameterConstantness(true);
4032 
4033  return ret;
4034 }
__host__ void setParameterConstantness(bool constant=true)
std::vector< Variable > weights
GooPdf * makeBkg_sigma_strips(int bkgnum)
bool notUseBackground4Hist
GooPdf * bkg4_jsugg
std::vector< PdfBase * > comps
bool gaussBkgTime
GooPdf * makeKzeroVeto()
GooPdf * makeGaussianTimePdf(int bkg)
GooPdf * makeExpGausTimePdf(int bkg)
GooPdf * kzero_veto
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88
GooPdf * makeBackground4DalitzParam()
SmoothHistogramPdf * makeBackgroundHistogram(int bkgnum, std::string overridename="")

◆ makeBkg_sigma_strips()

GooPdf* makeBkg_sigma_strips ( int  bkgnum)

Definition at line 2570 of file pipipi0DPFit.cpp.

References jsuList, m23Slices, and make_m23_transform().

Referenced by makeBkg2DalitzPdf(), makeBkg3DalitzPdf(), makeBkg4DalitzPdf(), makeOverallSignal(), runBackgroundSigmaFit(), and runSigmaFit().

2570  {
2571  GooPdf *m23_composite = make_m23_transform();
2572 
2573  vector<GooPdf *> jsuList;
2574  vector<ConvolutionPdf *> convList;
2575  bool useShare = false;
2576 
2577  for(int i = 0; i < m23Slices; ++i) {
2578  sprintf(strbuffer, "bkg%i_sigma_slice%i_expalpha", bkgnum, i);
2579  Variable exp_alpha(strbuffer, 7.50, 0.10, 0, 10.00);
2580  sprintf(strbuffer, "bkg%i_sigma_slice%i_gauss_meana", bkgnum, i);
2581  Variable g_meana(strbuffer, 0.20, 0.01, 0.00, 0.80);
2582  sprintf(strbuffer, "bkg%i_sigma_slice%i_gauss_sigma", bkgnum, i);
2583  Variable g_sigma(strbuffer, 0.03, 0.01, 0.01, 0.20);
2584 
2585  sprintf(strbuffer, "bkg%i_sigma_slice%i_conv", bkgnum, i);
2586  ExpGausPdf *expfunc = new ExpGausPdf(strbuffer, *sigma, g_meana, g_sigma, exp_alpha);
2587  jsuList.push_back(expfunc);
2588  }
2589 
2590  if(useShare) {
2591  for(vector<ConvolutionPdf *>::iterator conv = convList.begin(); conv != convList.end(); ++conv) {
2592  (*conv)->registerOthers(convList);
2593  }
2594  }
2595 
2596  sprintf(strbuffer, "bkg%i_sigma_map", bkgnum);
2597  return new MappedPdf(strbuffer, m23_composite, jsuList);
2598 }
vector< GooPdf * > jsuList
int m23Slices
char strbuffer[1000]
Observable * sigma
GooPdf * make_m23_transform()

◆ makeDalitzPlots()

void makeDalitzPlots ( GooPdf overallSignal,
std::string  plotdir = "./plots_from_mixfit/" 
)

Definition at line 1850 of file pipipi0DPFit.cpp.

References _mD0, GooFit::UnbinnedDataSet::addEvent(), ChisqInfo::chisq, coarseBin(), ChisqInfo::contribPlot, cpuDalitz(), cpuGetM23(), ChisqInfo::dof, foo, foodal, getAdaptiveChisquare(), GooFit::GooPdf::getCompProbsAtDataPoints(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::DataSet::getNumEvents(), GooFit::PdfBase::getObservables(), GooFit::Indexable::getUpperLimit(), GooFit::UnbinnedDataSet::getValue(), GooFit::Indexable::getValue(), piPlusMass, piZeroMass, GooFit::PdfBase::setData(), GooFit::IncoherentSumPdf::setDataSize(), GooFit::TddpPdf::setDataSize(), GooFit::Indexable::setValue(), GooFit::UnbinnedDataSet::setValueForAllEvents(), and sigma.

Referenced by runCanonicalFit(), and runTruthMCFit().

1850  {
1851  std::string mkplotdir{"mkdir " + plotdir};
1852  system(mkplotdir.c_str());
1853  foo->cd();
1854 
1855  TH1F dtime_dat_hist("dtime_dat_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1856  dtime_dat_hist.SetStats(false);
1857  dtime_dat_hist.SetMarkerStyle(8);
1858  dtime_dat_hist.SetMarkerSize(1.2);
1859  dtime_dat_hist.GetXaxis()->SetTitle("Decay time [ps]");
1860  dtime_dat_hist.GetYaxis()->SetTitle("Events / 50 fs");
1861  TH1F dtime_pdf_hist("dtime_pdf_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1862  dtime_pdf_hist.SetStats(false);
1863  dtime_pdf_hist.SetLineColor(kBlue);
1864  dtime_pdf_hist.SetLineWidth(3);
1865 
1866  TH1F sigma_dat_hist("sigma_dat_hist", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1867  sigma_dat_hist.SetStats(false);
1868  sigma_dat_hist.SetMarkerStyle(8);
1869  sigma_dat_hist.SetMarkerSize(1.2);
1870  sigma_dat_hist.GetXaxis()->SetTitle("Decay time error [ps]");
1871  sigma_dat_hist.GetYaxis()->SetTitle("Events / 8 fs");
1872  TH1F sigma_pdf_hist("sigma_pdf_hist", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1873  sigma_pdf_hist.SetStats(false);
1874  sigma_pdf_hist.SetLineColor(kBlue);
1875  sigma_pdf_hist.SetLineWidth(3);
1876 
1877  TH1F m12_dat_hist("m12_dat_hist", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit());
1878  m12_dat_hist.SetStats(false);
1879  m12_dat_hist.SetMarkerStyle(8);
1880  m12_dat_hist.SetMarkerSize(1.2);
1881  m12_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1882  m12_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1883  TH1F m12_pdf_hist("m12_pdf_hist", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit());
1884  m12_pdf_hist.SetStats(false);
1885  m12_pdf_hist.SetLineColor(kBlue);
1886  m12_pdf_hist.SetLineWidth(3);
1887 
1888  TH1F m13_dat_hist("m13_dat_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1889  m13_dat_hist.SetStats(false);
1890  m13_dat_hist.SetMarkerStyle(8);
1891  m13_dat_hist.SetMarkerSize(1.2);
1892  m13_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1893  m13_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1894  TH1F m13_pdf_hist("m13_pdf_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1895  m13_pdf_hist.SetStats(false);
1896  m13_pdf_hist.SetLineColor(kBlue);
1897  m13_pdf_hist.SetLineWidth(3);
1898 
1899  TH1F m23_dat_hist("m23_dat_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1900  m23_dat_hist.SetStats(false);
1901  m23_dat_hist.SetMarkerStyle(8);
1902  m23_dat_hist.SetMarkerSize(1.2);
1903  m23_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{+}) [GeV]");
1904  m23_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1905  TH1F m23_pdf_hist("m23_pdf_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1906  m23_pdf_hist.SetStats(false);
1907  m23_pdf_hist.SetLineColor(kBlue);
1908  m23_pdf_hist.SetLineWidth(3);
1909 
1910  TH2F dalitzpm_dat_hist("dalitzpm_dat_hist",
1911  "",
1912  m12->getNumBins(),
1913  m12->getLowerLimit(),
1914  m12->getUpperLimit(),
1915  m13->getNumBins(),
1916  m13->getLowerLimit(),
1917  m13->getUpperLimit());
1918  dalitzpm_dat_hist.SetStats(false);
1919  dalitzpm_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1920  dalitzpm_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1921  TH2F dalitzpm_pdf_hist("dalitzpm_pdf_hist",
1922  "",
1923  m12->getNumBins(),
1924  m12->getLowerLimit(),
1925  m12->getUpperLimit(),
1926  m13->getNumBins(),
1927  m13->getLowerLimit(),
1928  m13->getUpperLimit());
1929  dalitzpm_pdf_hist.SetStats(false);
1930 
1931  TH2F dalitzp0_dat_hist("dalitzp0_dat_hist",
1932  "",
1933  m12->getNumBins(),
1934  m12->getLowerLimit(),
1935  m12->getUpperLimit(),
1936  m13->getNumBins(),
1937  m13->getLowerLimit(),
1938  m13->getUpperLimit());
1939  dalitzp0_dat_hist.SetStats(false);
1940  dalitzp0_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1941  dalitzp0_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{-} #pi^{+}) [GeV]");
1942  TH2F dalitzp0_pdf_hist("dalitzp0_pdf_hist",
1943  "",
1944  m12->getNumBins(),
1945  m12->getLowerLimit(),
1946  m12->getUpperLimit(),
1947  m13->getNumBins(),
1948  m13->getLowerLimit(),
1949  m13->getUpperLimit());
1950  dalitzp0_pdf_hist.SetStats(false);
1951 
1952  TH2F dalitzm0_dat_hist("dalitzm0_dat_hist",
1953  "",
1954  m12->getNumBins(),
1955  m12->getLowerLimit(),
1956  m12->getUpperLimit(),
1957  m13->getNumBins(),
1958  m13->getLowerLimit(),
1959  m13->getUpperLimit());
1960  dalitzm0_dat_hist.SetStats(false);
1961  dalitzm0_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1962  dalitzm0_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{+} #pi^{-}) [GeV]");
1963  TH2F dalitzm0_pdf_hist("dalitzm0_pdf_hist",
1964  "",
1965  m12->getNumBins(),
1966  m12->getLowerLimit(),
1967  m12->getUpperLimit(),
1968  m13->getNumBins(),
1969  m13->getLowerLimit(),
1970  m13->getUpperLimit());
1971  dalitzm0_pdf_hist.SetStats(false);
1972 
1973  TH1F *bkg3_pdfs[6];
1974  TH1F *bkg3_data[6];
1975  double num_sigma_dat[6];
1976  double num_sigma_pdf[6];
1977 
1978  for(int i = 0; i < 6; ++i) {
1979  sprintf(strbuffer, "bkg_sigma_pdf_%i", i);
1980  bkg3_pdfs[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1981  sprintf(strbuffer, "bkg_sigma_dat_%i", i);
1982  bkg3_data[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1983 
1984  num_sigma_dat[i] = 0;
1985  num_sigma_pdf[i] = 0;
1986 
1987  bkg3_data[i]->SetStats(false);
1988  bkg3_data[i]->SetMarkerStyle(8);
1989  bkg3_data[i]->SetMarkerSize(1.2);
1990  bkg3_data[i]->GetXaxis()->SetTitle("Decay time error [ps]");
1991  bkg3_data[i]->GetYaxis()->SetTitle("Events / 8 fs");
1992 
1993  bkg3_pdfs[i]->SetStats(false);
1994  bkg3_pdfs[i]->SetLineColor(kBlue);
1995  bkg3_pdfs[i]->SetLineWidth(3);
1996  }
1997 
1998  double totalPdf = 0;
1999  double totalDat = 0;
2000 
2001  for(unsigned int evt = 0; evt < data->getNumEvents(); ++evt) {
2002  double currTime = data->getValue(*dtime, evt);
2003  dtime_dat_hist.Fill(currTime);
2004 
2005  double currSigma = data->getValue(*sigma, evt);
2006  sigma_dat_hist.Fill(currSigma);
2007 
2008  double currm12 = data->getValue(*m12, evt);
2009  m12_dat_hist.Fill(currm12);
2010 
2011  double currm13 = data->getValue(*m13, evt);
2012  m13_dat_hist.Fill(currm13);
2013 
2014  dalitzpm_dat_hist.Fill(currm12, currm13);
2015 
2016  double currm23 = cpuGetM23(currm12, currm13);
2017  m23_dat_hist.Fill(currm23);
2018 
2019  dalitzp0_dat_hist.Fill(currm12, currm23);
2020  dalitzm0_dat_hist.Fill(currm13, currm23);
2021 
2022  totalDat++;
2023 
2024  int m23bin = (int)floor(currm23 / 0.5);
2025  bkg3_data[m23bin]->Fill(currSigma);
2026  num_sigma_dat[m23bin]++;
2027  }
2028 
2029  double maxBinContent = 0;
2030  int bestI = 0;
2031  int bestJ = 0;
2032 
2033  for(int i = 1; i <= m12->getNumBins(); ++i) {
2034  for(int j = 1; j <= m13->getNumBins(); ++j) {
2035  double curr = dalitzpm_dat_hist.GetBinContent(i, j);
2036 
2037  if(curr < maxBinContent)
2038  continue;
2039 
2040  maxBinContent = curr;
2041  bestI = i;
2042  bestJ = j;
2043  }
2044  }
2045 
2046  std::cout << "Max bin content: " << maxBinContent << " (" << bestI << ", " << bestJ << ")\n";
2047 
2048  bool dependsOnSigma = true;
2049  std::vector<Observable> obses = overallSignal->getObservables();
2050 
2051  if(std::find(obses.begin(), obses.end(), *sigma) == obses.end())
2052  dependsOnSigma = false;
2053 
2054  // overallSignal->setDebugMask(1);
2055 
2056  wBkg1->setValue(0);
2057  const int division = 2;
2058 
2059  for(int half = 0; half < division; ++half) {
2060  std::vector<Observable> vars;
2061  vars.push_back(*m12);
2062  vars.push_back(*m13);
2063  vars.push_back(*dtime);
2064  vars.push_back(*sigma);
2065  vars.push_back(*eventNumber);
2066  vars.push_back(*wBkg1);
2067 
2068  UnbinnedDataSet currData(vars);
2069 
2070  sigma->setValue(0.5);
2071  int evtCounter = 0;
2072 
2073  for(int i = 0; i < m12->getNumBins(); ++i) {
2075  + (m12->getUpperLimit() - m12->getLowerLimit()) * (i + 0.5) / m12->getNumBins());
2076 
2077  for(int j = 0; j < m13->getNumBins(); ++j) {
2079  + (m13->getUpperLimit() - m13->getLowerLimit()) * (j + 0.5) / m13->getNumBins());
2080 
2082  continue;
2083 
2084  for(int l = half; l < dtime->getNumBins(); l += division) {
2086  + (dtime->getUpperLimit() - dtime->getLowerLimit()) * (l + 0.5)
2087  / dtime->getNumBins());
2088  eventNumber->setValue(evtCounter);
2089  evtCounter++;
2090  currData.addEvent();
2091  }
2092  }
2093  }
2094 
2095  for(int k = 0; k < sigma->getNumBins(); ++k) {
2097  + (sigma->getUpperLimit() - sigma->getLowerLimit()) * (k + 0.5) / sigma->getNumBins());
2098 
2099  if(0 == k % 10)
2100  std::cout << "sigma iteration " << half << " " << k << std::endl;
2101 
2102  currData.setValueForAllEvents(*sigma);
2103  overallSignal->setData(&currData);
2104 
2105  if(0 == k) {
2106  if(signalDalitz)
2107  signalDalitz->setDataSize(currData.getNumEvents(), 6);
2108 
2109  if(incsum1)
2110  incsum1->setDataSize(currData.getNumEvents(), 6);
2111 
2112  if(incsum2)
2113  incsum2->setDataSize(currData.getNumEvents(), 6);
2114 
2115  if(incsum3)
2116  incsum3->setDataSize(currData.getNumEvents(), 6);
2117 
2118  if(incsum4)
2119  incsum4->setDataSize(currData.getNumEvents(), 6);
2120 
2121  if(incsum5)
2122  incsum5->setDataSize(currData.getNumEvents(), 6);
2123 
2124  if(incsum6)
2125  incsum6->setDataSize(currData.getNumEvents(), 6);
2126  }
2127 
2128  std::vector<std::vector<double>> pdfValues = overallSignal->getCompProbsAtDataPoints();
2129 
2130  for(unsigned int j = 0; j < pdfValues[0].size(); ++j) {
2131  double currTime = currData.getValue(*dtime, j);
2132  dtime_pdf_hist.Fill(currTime, pdfValues[0][j]);
2133 
2134  double currSigma = currData.getValue(*sigma, j);
2135  sigma_pdf_hist.Fill(currSigma, pdfValues[0][j]);
2136 
2137  // Um... these two are switched? Weirdness...
2138  double currm12 = currData.getValue(*m13, j);
2139  m12_pdf_hist.Fill(currm12, pdfValues[0][j]);
2140 
2141  double currm13 = currData.getValue(*m12, j);
2142  m13_pdf_hist.Fill(currm13, pdfValues[0][j]);
2143  dalitzpm_pdf_hist.Fill(currm12, currm13, pdfValues[0][j]);
2144 
2145  double currm23 = cpuGetM23(currm12, currm13);
2146  m23_pdf_hist.Fill(currm23, pdfValues[0][j]);
2147  dalitzp0_pdf_hist.Fill(currm12, currm23, pdfValues[0][j]);
2148  dalitzm0_pdf_hist.Fill(currm13, currm23, pdfValues[0][j]);
2149 
2150  int currM23Bin = (int)(currm23 / 0.5);
2151  bkg3_pdfs[currM23Bin]->Fill(currSigma, pdfValues[0][j]);
2152  num_sigma_pdf[currM23Bin] += pdfValues[0][j];
2153 
2154  totalPdf += pdfValues[0][j];
2155 
2156  if(std::isnan(pdfValues[0][j])) {
2157  std::cout << "Major problem: " << k << " " << j << std::endl;
2158  assert(false);
2159  }
2160 
2161  if(std::isinf(pdfValues[0][j])) {
2162  std::cout << "Infinity " << k << " " << j << std::endl;
2163  assert(false);
2164  }
2165  }
2166 
2167  // If PDF doesn't depend on sigma, don't project from that dimension.
2168  if(!dependsOnSigma)
2169  break;
2170  }
2171  }
2172 
2173  // std::cout << "Final values: " << totalDat << " " << totalPdf << std::endl;
2174 
2175  for(int i = 1; i <= dtime->getNumBins(); ++i) {
2176  dtime_pdf_hist.SetBinContent(i, dtime_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2177  }
2178 
2179  for(int i = 1; i <= sigma->getNumBins(); ++i) {
2180  sigma_pdf_hist.SetBinContent(i, sigma_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2181 
2182  for(int j = 0; j < 6; ++j) {
2183  bkg3_pdfs[j]->SetBinContent(i, bkg3_pdfs[j]->GetBinContent(i) * num_sigma_dat[j] / num_sigma_pdf[j]);
2184  }
2185  }
2186 
2187  for(int i = 1; i <= m12->getNumBins(); ++i) {
2188  m12_pdf_hist.SetBinContent(i, m12_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2189  }
2190 
2191  for(int i = 1; i <= m13->getNumBins(); ++i) {
2192  m13_pdf_hist.SetBinContent(i, m13_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2193  m23_pdf_hist.SetBinContent(i, m23_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2194  }
2195 
2196  for(int i = 1; i <= m12->getNumBins(); ++i) {
2197  for(int j = 1; j <= m13->getNumBins(); ++j) {
2198  dalitzpm_pdf_hist.SetBinContent(i, j, dalitzpm_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2199  dalitzp0_pdf_hist.SetBinContent(i, j, dalitzp0_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2200  dalitzm0_pdf_hist.SetBinContent(i, j, dalitzm0_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2201  }
2202  }
2203 
2204  ChisqInfo *chisq = getAdaptiveChisquare(&dalitzpm_dat_hist, &dalitzpm_pdf_hist);
2205 
2206  std::cout << "Chisquare: " << chisq->chisq << " / " << chisq->dof << std::endl;
2207  foodal->cd();
2208  foodal->SetLogz(false);
2209  chisq->contribPlot->SetStats(false);
2210  chisq->contribPlot->Draw("colz");
2211  foodal->SaveAs((plotdir + "/chisq.png").c_str());
2212  foo->cd();
2213 
2214  coarseBin(dalitzpm_pdf_hist, 2);
2215  coarseBin(dalitzpm_dat_hist, 2);
2216  coarseBin(dalitzp0_pdf_hist, 2);
2217  coarseBin(dalitzp0_dat_hist, 2);
2218  coarseBin(dalitzm0_pdf_hist, 2);
2219  coarseBin(dalitzm0_dat_hist, 2);
2220 
2221  dtime_dat_hist.Draw("p");
2222  dtime_pdf_hist.Draw("lsame");
2223  foo->SaveAs((plotdir + "/dtime_fit.png").c_str());
2224  foo->SetLogy(true);
2225  foo->SaveAs((plotdir + "/dtime_fit_log.png").c_str());
2226  foo->SetLogy(false);
2227 
2228  for(int i = 0; i < 6; ++i) {
2229  if(!dependsOnSigma) {
2230  bkg3_data[i]->Draw("p");
2231  } else if(bkg3_data[i]->GetMaximum() > bkg3_pdfs[i]->GetMaximum()) {
2232  bkg3_data[i]->Draw("p");
2233  bkg3_pdfs[i]->Draw("lsame");
2234  } else {
2235  bkg3_pdfs[i]->Draw("l");
2236  bkg3_data[i]->Draw("psame");
2237  }
2238 
2239  sprintf(strbuffer, "%i", i);
2240  TText slicenum;
2241  slicenum.DrawTextNDC(0.2, 0.8, strbuffer);
2242 
2243  foo->SaveAs((plotdir + bkg3_pdfs[i]->GetName() + ".png").c_str());
2244  foo->SetLogy(true);
2245  foo->SaveAs((plotdir + bkg3_pdfs[i]->GetName() + "_log.png").c_str());
2246  foo->SetLogy(false);
2247  }
2248 
2249  m13_dat_hist.Draw("p");
2250  m13_pdf_hist.Draw("lsame");
2251  foo->SaveAs((plotdir + "/m13_fit.png").c_str());
2252  foo->SetLogy(true);
2253  foo->SaveAs((plotdir + "/m13_fit_log.png").c_str());
2254  foo->SetLogy(false);
2255 
2256  sigma_dat_hist.Draw("p");
2257  sigma_pdf_hist.Draw("lsame");
2258  foo->SaveAs((plotdir + "/sigma_fit.png").c_str());
2259  foo->SetLogy(true);
2260  foo->SaveAs((plotdir + "/sigma_fit_log.png").c_str());
2261  foo->SetLogy(false);
2262 
2263  m12_dat_hist.Draw("p");
2264  m12_pdf_hist.Draw("lsame");
2265  foo->SaveAs((plotdir + "/m12_fit.png").c_str());
2266  foo->SetLogy(true);
2267  foo->SaveAs((plotdir + "/m12_fit_log.png").c_str());
2268  foo->SetLogy(false);
2269 
2270  m23_dat_hist.Draw("p");
2271  m23_pdf_hist.Draw("lsame");
2272  foo->SaveAs((plotdir + "/m23_fit.png").c_str());
2273  foo->SetLogy(true);
2274  foo->SaveAs((plotdir + "/m23_fit_log.png").c_str());
2275  foo->SetLogy(false);
2276 
2277  foodal->cd();
2278  dalitzpm_dat_hist.Draw("colz");
2279 
2280  for(int slice = 0; slice < 6; ++slice) {
2281  double line_m12 = cpuGetM23(0, 0.5 * (slice + 1));
2282  TLine sliceLine;
2283  sliceLine.SetLineWidth(2);
2284  sliceLine.DrawLine(0, line_m12, line_m12, 0);
2285  sprintf(strbuffer, "%i", slice);
2286  TText sliceNum;
2287  sliceNum.DrawText(0.25, line_m12 - 0.25, strbuffer);
2288  }
2289 
2290  foodal->SaveAs((plotdir + "/dalitzpm_dat.png").c_str());
2291  foodal->SetLogz(true);
2292  foodal->SaveAs((plotdir + "/dalitzpm_dat_log.png").c_str());
2293  foodal->SetLogz(false);
2294  dalitzp0_dat_hist.Draw("colz");
2295  foodal->SaveAs((plotdir + "/dalitzp0_dat.png").c_str());
2296  foodal->SetLogz(true);
2297  foodal->SaveAs((plotdir + "/dalitzp0_dat_log.png").c_str());
2298  foodal->SetLogz(false);
2299  dalitzm0_dat_hist.Draw("colz");
2300  foodal->SaveAs((plotdir + "/dalitzm0_dat.png").c_str());
2301  foodal->SetLogz(true);
2302  foodal->SaveAs((plotdir + "/dalitzm0_dat_log.png").c_str());
2303  foodal->SetLogz(false);
2304 
2305  dalitzpm_pdf_hist.Draw("colz");
2306  foodal->SaveAs((plotdir + "/dalitzpm_pdf.png").c_str());
2307  foodal->SetLogz(true);
2308  foodal->SaveAs((plotdir + "/dalitzpm_pdf_log.png").c_str());
2309  foodal->SetLogz(false);
2310  dalitzpm_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzpm_dat_hist.GetMaximum());
2311  dalitzpm_pdf_hist.Draw("colz");
2312  foodal->SaveAs((plotdir + "/dalitzpm_pdf_matched.png").c_str());
2313  foodal->SetLogz(true);
2314  foodal->SaveAs((plotdir + "/dalitzpm_pdf_matched_log.png").c_str());
2315  foodal->SetLogz(false);
2316  dalitzp0_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzp0_dat_hist.GetMaximum());
2317  dalitzp0_pdf_hist.Draw("colz");
2318  foodal->SaveAs((plotdir + "/dalitzp0_pdf.png").c_str());
2319  foodal->SetLogz(true);
2320  foodal->SaveAs((plotdir + "/dalitzp0_pdf_log.png").c_str());
2321  foodal->SetLogz(false);
2322  dalitzm0_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzm0_dat_hist.GetMaximum());
2323  dalitzm0_pdf_hist.Draw("colz");
2324  foodal->SaveAs((plotdir + "/dalitzm0_pdf.png").c_str());
2325  foodal->SetLogz(true);
2326  foodal->SaveAs((plotdir + "/dalitzm0_pdf_log.png").c_str());
2327  foodal->SetLogz(false);
2328 
2329  TH1F pull_pm_hist("pull_pm_hist", "", 100, -5, 5);
2330  pull_pm_hist.GetXaxis()->SetTitle("(Data - PDF) / sqrt(Data)");
2331  pull_pm_hist.GetYaxis()->SetTitle("Bins / 0.1");
2332 
2333  for(int i = 1; i <= m12->getNumBins(); ++i) {
2334  double m12loedge
2335  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i - 1);
2336  double m12hiedge
2337  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i);
2338 
2339  for(int j = 1; j <= m13->getNumBins(); ++j) {
2340  double m13loedge
2341  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2342 
2343  if(!cpuDalitz(m12loedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2344  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2345  continue;
2346  }
2347 
2348  if(!cpuDalitz(m12hiedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2349  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2350  continue;
2351  }
2352 
2353  double m13hiedge
2354  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2355 
2356  if(!cpuDalitz(m12loedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2357  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2358  continue;
2359  }
2360 
2361  if(!cpuDalitz(m12hiedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2362  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2363  continue;
2364  }
2365 
2366  double dat = dalitzpm_dat_hist.GetBinContent(i, j);
2367  double pdf = dalitzpm_pdf_hist.GetBinContent(i, j);
2368 
2369  double pullval = (dat - pdf) / sqrt(max(1.0, dat));
2370  dalitzpm_dat_hist.SetBinContent(i, j, pullval);
2371  pull_pm_hist.Fill(pullval);
2372  }
2373 
2374  for(int j = 1; j <= m13->getNumBins(); ++j) {
2375  double m23loedge
2376  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2377 
2378  // To get 12, 23 instead of 12, 13, just exchange 1 and 2.
2379  if(!cpuDalitz(m12loedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2380  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2381  continue;
2382  }
2383 
2384  if(!cpuDalitz(m12hiedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2385  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2386  continue;
2387  }
2388 
2389  double m23hiedge
2390  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2391 
2392  if(!cpuDalitz(m12loedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2393  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2394  continue;
2395  }
2396 
2397  if(!cpuDalitz(m12hiedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2398  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2399  continue;
2400  }
2401 
2402  double dat = dalitzp0_dat_hist.GetBinContent(i, j);
2403  double pdf = dalitzp0_pdf_hist.GetBinContent(i, j);
2404 
2405  dalitzp0_dat_hist.SetBinContent(i, j, (dat - pdf) / sqrt(max(1.0, dat)));
2406  }
2407 
2408  // NB, this exploits symmetry 12 and 13 by treating the outer loop as 13.
2409  for(int j = 1; j <= m13->getNumBins(); ++j) {
2410  double m23loedge
2411  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2412 
2413  if(!cpuDalitz(m12loedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2414  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2415  continue;
2416  }
2417 
2418  if(!cpuDalitz(m12hiedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2419  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2420  continue;
2421  }
2422 
2423  double m23hiedge
2424  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2425 
2426  if(!cpuDalitz(m12loedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2427  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2428  continue;
2429  }
2430 
2431  if(!cpuDalitz(m12hiedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2432  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2433  continue;
2434  }
2435 
2436  double dat = dalitzm0_dat_hist.GetBinContent(i, j);
2437  double pdf = dalitzm0_pdf_hist.GetBinContent(i, j);
2438 
2439  dalitzm0_dat_hist.SetBinContent(i, j, (dat - pdf) / sqrt(max(1.0, dat)));
2440  }
2441  }
2442 
2443  dalitzpm_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2444  dalitzpm_dat_hist.Draw("colz");
2445  foodal->SaveAs((plotdir + "/dalitzpm_pull.png").c_str());
2446  dalitzp0_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2447  dalitzp0_dat_hist.Draw("colz");
2448  foodal->SaveAs((plotdir + "/dalitzp0_pull.png").c_str());
2449  dalitzm0_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2450  dalitzm0_dat_hist.Draw("colz");
2451  foodal->SaveAs((plotdir + "/dalitzm0_pull.png").c_str());
2452 
2453  foo->cd();
2454 
2455  pull_pm_hist.Draw("");
2456  foo->SaveAs((plotdir + "/pull_pm_hist.png").c_str());
2457 }
IncoherentSumPdf * incsum2
size_t getNumEvents() const
Definition: DataSet.h:47
fptype cpuGetM23(fptype massPZ, fptype massPM)
char strbuffer[1000]
void coarseBin(TH2F &dalitzHist, int grain)
__host__ std::vector< std::vector< fptype > > getCompProbsAtDataPoints()
Produce a list of probabilies at points.
const fptype piPlusMass
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
Observable * m12
ChisqInfo * getAdaptiveChisquare(TH2F *datPlot, TH2F *pdfPlot)
fptype getValue(const Observable &var, size_t idx) const
Get the value at a specific variable and event number.
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * data
TddpPdf * signalDalitz
IncoherentSumPdf * incsum3
TCanvas * foodal
const fptype _mD0
__host__ void setDataSize(unsigned int dataSize, unsigned int evtSize=5)
__host__ void setData(DataSet *data)
Observable * m13
Observable * sigma
const fptype piZeroMass
EventNumber * eventNumber
virtual __host__ std::vector< Observable > getObservables() const
Definition: PdfBase.cpp:111
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
Observable * wBkg1
IncoherentSumPdf * incsum5
IncoherentSumPdf * incsum1
Observable * dtime
TCanvas * foo
__host__ void setDataSize(unsigned int dataSize, unsigned int evtSize=3)
TH2F * contribPlot
IncoherentSumPdf * incsum6
IncoherentSumPdf * incsum4

◆ makeEfficiencyPdf()

GooPdf* makeEfficiencyPdf ( )

Definition at line 433 of file pipipi0DPFit.cpp.

Referenced by makeOverallSignal(), runEfficiencyFit(), and runGeneratedMCFit().

433  {
434  Variable effSmoothing("effSmoothing", 1.0, 0.1, 0, 1.25);
435  // Variable* effSmoothing = new Variable("effSmoothing", 0);
436  SmoothHistogramPdf *ret = new SmoothHistogramPdf("efficiency", binEffData, effSmoothing);
437  return ret;
438 }
BinnedDataSet * binEffData

◆ makeEfficiencyPoly()

GooPdf* makeEfficiencyPoly ( )

Definition at line 452 of file pipipi0DPFit.cpp.

References constantOne, makeKzeroVeto(), massSum, maxDalitzX, maxDalitzY, maxDalitzZ, minDalitzX, minDalitzY, minDalitzZ, and GooFit::GooPdf::setParameterConstantness().

Referenced by makeOverallSignal().

452  {
453  if(!kzero_veto)
454  makeKzeroVeto();
455 
456  vector<Variable> offsets;
457  vector<Observable> observables;
458  vector<Variable> coefficients;
459  offsets.push_back(constantOne);
460  offsets.push_back(constantOne);
461 
462  observables.push_back(*m12);
463  observables.push_back(*m13);
464 
465  coefficients.push_back(Variable("x0y0", 1.0));
466  coefficients.push_back(Variable("x1y0", 0.07999, 0.01, -0.5, 0.5));
467  coefficients.push_back(Variable("x2y0", -0.23732, 0.01, -0.5, 0.5));
468  coefficients.push_back(Variable("x3y0", 0.10369, 0.01, -1.5, 0.5));
469  coefficients.push_back(Variable("x0y1", 0.10248, 0.01, -1.5, 0.5));
470  coefficients.push_back(Variable("x1y1", -0.28543, 0.01, -0.5, 0.5));
471  coefficients.push_back(Variable("x2y1", 0.15058, 0.01, -1.5, 0.5));
472  coefficients.push_back(Variable("x0y2", -0.20648, 0.01, -0.5, 0.5));
473  coefficients.push_back(Variable("x1y2", 0.14567, 0.01, -1.5, 0.5));
474  coefficients.push_back(Variable("x0y3", 0.06231, 0.01, -0.5, 0.5));
475 
476  PolynomialPdf *poly = new PolynomialPdf("efficiency", observables, coefficients, offsets, 3);
477  poly->setParameterConstantness(true);
478 
479  Variable decXmin("decXmin", 6.22596);
480  Variable decYmin("decYmin", 6.30722);
481  Variable decZmin("decZmin", 10.82390);
482  Variable conXmin("conXmin", 0.65621);
483  Variable conYmin("conYmin", 0.69527);
484  Variable conZmin("conZmin", 0.31764);
485  Variable decXmax("decXmax", 0.79181);
486  Variable decYmax("decYmax", 5.91211);
487  Variable decZmax("decZmax", 1.52031);
488  Variable conXmax("conXmax", 0.80918);
489  Variable conYmax("conYmax", 0.22158);
490  Variable conZmax("conZmax", 0.41866);
491 
492  TrigThresholdPdf *loX = new TrigThresholdPdf("loX", *m12, minDalitzX, decXmin, conXmin, false);
493  TrigThresholdPdf *hiX = new TrigThresholdPdf("hiX", *m12, maxDalitzX, decXmax, conXmax, true);
494 
495  TrigThresholdPdf *loY = new TrigThresholdPdf("loY", *m13, minDalitzY, decYmin, conYmin, false);
496  TrigThresholdPdf *hiY = new TrigThresholdPdf("hiY", *m13, maxDalitzY, decYmax, conYmax, true);
497 
498  TrigThresholdPdf *loZ = new TrigThresholdPdf("loZ", *m12, *m13, minDalitzZ, decZmin, conZmin, massSum, false);
499  TrigThresholdPdf *hiZ = new TrigThresholdPdf("hiZ", *m12, *m13, maxDalitzZ, decZmax, conZmax, massSum, true);
500 
501  comps.clear();
502  comps.push_back(poly);
503  comps.push_back(loX);
504  comps.push_back(hiX);
505  comps.push_back(loY);
506  comps.push_back(hiY);
507  comps.push_back(loZ);
508  comps.push_back(hiZ);
509  comps.push_back(kzero_veto);
510  ProdPdf *ret = new ProdPdf("efficiency_total", comps);
511 
512  // return poly;
513  return ret;
514 }
__host__ void setParameterConstantness(bool constant=true)
Variable maxDalitzY("maxDalitzY", pow(_mD0 - piPlusMass, 2))
Variable massSum("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
Variable minDalitzX("minDalitzX", pow(piPlusMass+piZeroMass, 2))
Observable * m12
Variable minDalitzY("minDalitzY", pow(piPlusMass+piZeroMass, 2))
Variable minDalitzZ("minDalitzZ", pow(piPlusMass+piPlusMass, 2))
Variable constantOne("constantOne", 1)
Observable * m13
std::vector< PdfBase * > comps
Variable maxDalitzX("maxDalitzX", pow(_mD0 - piPlusMass, 2))
Variable maxDalitzZ("maxDalitzZ", pow(_mD0 - piZeroMass, 2))
GooPdf * makeKzeroVeto()
GooPdf * kzero_veto

◆ makeExpGausTimePdf()

GooPdf* makeExpGausTimePdf ( int  bkg)

Definition at line 3240 of file pipipi0DPFit.cpp.

Referenced by makeBkg2DalitzPdf(), makeBkg3DalitzPdf(), and makeBkg4DalitzPdf().

3240  {
3241  std::string bkgname = "";
3242 
3243  switch(bkg) {
3244  case 4:
3245  bkgname = "bkg4";
3246  break;
3247 
3248  case 3:
3249  bkgname = "bkg3";
3250  break;
3251 
3252  case 2:
3253  default:
3254  bkgname = "bkg2";
3255  break;
3256  }
3257 
3258  Variable g1_mean((bkgname + "_dtime_gmean1"), 0, 0.01, -0.5, 0.5);
3259  Variable g1_sigm((bkgname + "_dtime_gsigm1"), 0.2, 0.01, 0.01, 0.8);
3260  Variable e1_alph((bkgname + "_dtime_alpha1"), 2.5, 0.01, 0.01, 7.5);
3261 
3262  Variable g2_mean((bkgname + "_dtime_gmean2"), -0.3, 0.01, -0.85, 0.85);
3263  Variable g2_sigm((bkgname + "_dtime_gsigm2"), 0.2, 0.01, 0.01, 0.8);
3264  Variable e2_alph((bkgname + "_dtime_alpha2"), 0.5, 0.01, 0.01, 10.0);
3265 
3266  ExpGausPdf *exp1 = new ExpGausPdf((bkgname + "_exp1"), *dtime, g1_mean, g1_sigm, e1_alph);
3267  ExpGausPdf *exp2 = new ExpGausPdf((bkgname + "_exp2"), *dtime, g2_mean, g2_sigm, e2_alph);
3268 
3269  Variable frac1((bkgname + "_dtime_frac1"), 0.1, 0.01, 0, 0.8);
3270 
3271  GooPdf *ret = new AddPdf((bkgname + "_dtime"), frac1, exp1, exp2);
3272 
3273  return ret;
3274 }
Observable * dtime

◆ makeFlatBkgDalitzPdf()

GooPdf* makeFlatBkgDalitzPdf ( bool  fixem = true)

Definition at line 814 of file pipipi0DPFit.cpp.

References constantOne, constantZero, and GooFit::GooPdf::setParameterConstantness().

Referenced by runToyFit().

814  {
815  vector<Variable> offsets;
816  vector<Observable> observables;
817  vector<Variable> coefficients;
818  offsets.push_back(constantZero);
819  offsets.push_back(constantZero);
820  observables.push_back(*m12);
821  observables.push_back(*m13);
822  coefficients.push_back(constantOne);
823 
824  PolynomialPdf *poly = new PolynomialPdf("flatbkgPdf", observables, coefficients, offsets, 0);
825  Variable g_mean("g_mean", toyBkgTimeMean, 0.01, -0.2, 0.5);
826  Variable g_sigma("g_sigma", toyBkgTimeRMS, 0.01, 0.15, 1.5);
827  GooPdf *gt = new GaussianPdf("flatbkg_timepdf", *dtime, g_mean, g_sigma);
828  comps.clear();
829  comps.push_back(poly);
830  comps.push_back(gt);
831  GooPdf *ret = new ProdPdf("flatbkg_total", comps);
832 
833  if(fixem)
834  ret->setParameterConstantness(true);
835 
836  return ret;
837 }
__host__ void setParameterConstantness(bool constant=true)
Variable constantZero("constantZero", 0)
Observable * m12
Variable constantOne("constantOne", 1)
Observable * m13
const float toyBkgTimeRMS
std::vector< PdfBase * > comps
const float toyBkgTimeMean
Observable * dtime

◆ makeFullFitVariables()

void makeFullFitVariables ( )

Definition at line 1099 of file pipipi0DPFit.cpp.

References lowerTime, maxSigma, normBinning, GooFit::Observable::setNumBins(), and upperTime.

Referenced by makeTimePlots(), runBackgroundDalitzFit(), runBackgroundSigmaFit(), runCanonicalFit(), runEfficiencyFit(), runGeneratedMCFit(), runSigmaFit(), runTruthMCFit(), and writeBackgroundHistograms().

1099  {
1100  static bool exists = false;
1101 
1102  if(exists)
1103  return;
1104 
1105  exists = true;
1106 
1107  dtime = new Observable("dtime", lowerTime, upperTime);
1108  dtime->setNumBins(floor((upperTime - lowerTime) / 0.05 + 0.5));
1109  sigma = new Observable("sigma", 0, maxSigma);
1110  sigma->setNumBins(floor((maxSigma / 0.01) + 0.5));
1111  m12 = new Observable("m12", 0, 3);
1112  m13 = new Observable("m13", 0, 3);
1115  eventNumber = new EventNumber("eventNumber", 0, INT_MAX);
1116  wSig0 = new Observable("wSig0", 0, 1);
1117  wBkg1 = new Observable("wBkg1", 0, 1);
1118  wBkg2 = new Observable("wBkg2", 0, 1);
1119  wBkg3 = new Observable("wBkg3", 0, 1);
1120  wBkg4 = new Observable("wBkg4", 0, 1);
1121 }
Observable * wBkg3
double maxSigma
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
double lowerTime
Observable * m12
Special class for observables. Used in DataSets.
Definition: Variable.h:109
Observable * wSig0
Observable * m13
Observable * sigma
double upperTime
EventNumber * eventNumber
Observable * wBkg4
int normBinning
Observable * wBkg2
Observable * wBkg1
Observable * dtime

◆ makeGaussianTimePdf()

GooPdf* makeGaussianTimePdf ( int  bkg)

Definition at line 3167 of file pipipi0DPFit.cpp.

Referenced by makeBkg2DalitzPdf(), makeBkg3DalitzPdf(), and makeBkg4DalitzPdf().

3167  {
3168  // Gaussians for decay time.
3169  Variable *frac_ga2;
3170  Variable *frac_ga3;
3171  Variable *g1_meana;
3172  Variable *g1_sigma;
3173  Variable *g2_meana;
3174  Variable *g2_sigma;
3175  Variable *g3_meana;
3176  Variable *g3_sigma;
3177 
3178  std::string bkgname = "";
3179 
3180  switch(bkg) {
3181  case 4:
3182  frac_ga2 = new Variable("bkg4_frac_ga2", 0.60780, 0.01, 0.40, 0.80);
3183  frac_ga3 = new Variable("bkg4_frac_ga3", 0.04776, 0.01, 0.00, 0.20);
3184  g1_meana = new Variable("bkg4_g1_meana", 0.10164, 0.01, 0.00, 0.80);
3185  g1_sigma = new Variable("bkg4_g1_sigma", 0.27504, 0.01, 0.10, 0.80);
3186  g2_meana = new Variable("bkg4_g2_meana", 0.19974, 0.01, 0.10, 0.80);
3187  g2_sigma = new Variable("bkg4_g2_sigma", 0.63765, 0.01, 0.40, 0.80);
3188  g3_meana = new Variable("bkg4_g3_meana", 0.45817, 0.01, 0.20, 0.80);
3189  g3_sigma = new Variable("bkg4_g3_sigma", 1.52905, 0.01, 1.40, 1.80);
3190  bkgname = "bkg4";
3191  break;
3192 
3193  case 3:
3194  frac_ga2 = new Variable("bkg3_frac_ga2", 0.51448, 0.01, 0.25, 0.75);
3195  frac_ga3 = new Variable("bkg3_frac_ga3", 0.04169, 0.01, 0.00, 0.40);
3196  g1_meana = new Variable("bkg3_g1_meana", 0.25101, 0.01, 0.01, 0.50);
3197  g1_sigma = new Variable("bkg3_g1_sigma", 0.31953, 0.01, 0.02, 0.50);
3198  g2_meana = new Variable("bkg3_g2_meana", 0.49139, 0.01, 0.25, 0.75);
3199  g2_sigma = new Variable("bkg3_g2_sigma", 0.65443, 0.01, 0.10, 1.00);
3200  g3_meana = new Variable("bkg3_g3_meana", 0.83600, 0.01, 0.50, 1.00);
3201  g3_sigma = new Variable("bkg3_g3_sigma", 1.51839, 0.01, 0.10, 2.00);
3202  bkgname = "bkg3";
3203  break;
3204 
3205  case 2:
3206  default:
3207  frac_ga2 = new Variable("frac_ga2", 0.48994, 0.01, 0.1, 0.6);
3208  frac_ga3 = new Variable("frac_ga3", 0.04721, 0.01, 0.0, 0.1);
3209  g1_meana = new Variable("g1_meana", 0.01216, 0.01, -0.1, 0.1);
3210  g1_sigma = new Variable("g1_sigma", 0.25813, 0.01, 0.15, 0.35);
3211  g2_meana = new Variable("g2_meana", 0.05335, 0.01, -0.1, 0.1);
3212  g2_sigma = new Variable("g2_sigma", 0.58651, 0.01, 0.5, 1.2);
3213  g3_meana = new Variable("g3_meana", 0.17451, 0.01, 0.1, 1.85);
3214  g3_sigma = new Variable("g3_sigma", 1.15125, 0.01, 0.5, 1.3);
3215  bkgname = "bkg2";
3216  break;
3217  }
3218 
3219  GaussianPdf *g1 = new GaussianPdf((bkgname + "_g1").c_str(), *dtime, *g1_meana, *g1_sigma);
3220  GaussianPdf *g2 = new GaussianPdf((bkgname + "_g2").c_str(), *dtime, *g2_meana, *g2_sigma);
3221  GaussianPdf *g3 = new GaussianPdf((bkgname + "_g3").c_str(), *dtime, *g3_meana, *g3_sigma);
3222 
3223  weights.clear();
3224  weights.push_back(*frac_ga2);
3225 
3226  if(3 == bkg)
3227  weights.push_back(*frac_ga3);
3228 
3229  comps.clear();
3230  comps.push_back(g2);
3231 
3232  if(3 == bkg)
3233  comps.push_back(g3);
3234 
3235  comps.push_back(g1);
3236  AddPdf *bkg_dtime = new AddPdf((bkgname + "_dtime").c_str(), weights, comps);
3237  return bkg_dtime;
3238 }
std::vector< Variable > weights
std::vector< PdfBase * > comps
Observable * dtime

◆ makeKzeroVeto()

GooPdf* makeKzeroVeto ( )

Definition at line 440 of file pipipi0DPFit.cpp.

References chargeM, kzero_veto, motherM, neutrlM, and GooFit::PAIR_23.

Referenced by makeBkg2DalitzPdf(), makeBkg3DalitzPdf(), makeBkg4DalitzPdf(), makeEfficiencyPoly(), makeOverallSignal(), runBackgroundDalitzFit(), and runEfficiencyFit().

440  {
441  if(kzero_veto)
442  return kzero_veto;
443 
444  VetoInfo kVetoInfo(Variable("veto_min", 0.475 * 0.475), Variable("veto_max", 0.505 * 0.505), PAIR_23);
445 
446  vector<VetoInfo> vetos;
447  vetos.push_back(kVetoInfo);
448  kzero_veto = new DalitzVetoPdf("kzero_veto", *m12, *m13, motherM, neutrlM, chargeM, chargeM, vetos);
449  return kzero_veto;
450 }
Variable neutrlM("neutrlM", 0.1349766)
Variable chargeM("chargeM", 0.13957018)
Observable * m12
Observable * m13
Variable motherM("motherM", 1.86484)
GooPdf * kzero_veto

◆ makeMikhailJSU_gg()

GooPdf* makeMikhailJSU_gg ( bool  fixem = true)

Definition at line 1207 of file pipipi0DPFit.cpp.

References GooFit::Indexable::getValue(), GooFit::Variable::setFixed(), and GooFit::Indexable::setValue().

1207  {
1208  // Values from TSigma, 'Mikhail default'.
1209 
1210  Variable *js_meana = new Variable("js_meana", 0.0593279, 0.01, 0, 0.2);
1211  js_meana->setFixed(fixem);
1212  Variable *js_sigma = new Variable("js_sigma", 0.000474171, 0.0001, 0, 0.001);
1213  js_sigma->setFixed(fixem);
1214  Variable *js_gamma = new Variable("js_gamma", -10.1942, 1, -30, 0);
1215  js_gamma->setFixed(fixem);
1216  Variable *js_delta = new Variable("js_delta", 1.4907, 0.1, 0.5, 5);
1217  js_delta->setFixed(fixem);
1218  Variable *frac_jsu = new Variable("frac_jsu", 0.951638, 0.01, 0.5, 1.0);
1219  frac_jsu->setFixed(fixem);
1220  Variable *frac_ga1 = new Variable("frac_ga1", 0.0184522, 0.00001, 0.001, 0.3);
1221  frac_ga1->setFixed(fixem);
1222  Variable *g1_meana = new Variable("g1_meana", 0.257796, 0.003, 0.1, 0.5);
1223  g1_meana->setFixed(fixem);
1224  Variable *g1_sigma = new Variable("g1_sigma", 0.0308619, 0.01, 0.005, 0.25);
1225  g1_sigma->setFixed(fixem);
1226  Variable *g2_meana = new Variable("g2_meana", 0.319952, 0.01, 0.1, 0.5);
1227  g2_meana->setFixed(fixem);
1228  Variable *g2_sigma = new Variable("g2_sigma", 0.0582524, 0.01, 0.01, 0.1);
1229  g2_sigma->setFixed(fixem);
1230 
1231  // Variable* g2_sigma = new Variable("g2_sigma", 0.5825);
1232 
1233  frac_ga1->setValue(frac_ga1->getValue() * (1 - frac_jsu->getValue()));
1234 
1235  JohnsonSUPdf *js = new JohnsonSUPdf("js", *sigma, *js_meana, *js_sigma, *js_gamma, *js_delta);
1236  GaussianPdf *g1 = new GaussianPdf("g1", *sigma, *g1_meana, *g1_sigma);
1237  GaussianPdf *g2 = new GaussianPdf("g2", *sigma, *g2_meana, *g2_sigma);
1238 
1239  weights.clear();
1240  weights.push_back(*frac_jsu);
1241  weights.push_back(*frac_ga1);
1242  comps.clear();
1243  comps.push_back(js);
1244  comps.push_back(g1);
1245  comps.push_back(g2);
1246 
1247  GooPdf *ret = new AddPdf("signal_sigma", weights, comps);
1248  return ret;
1249 }
void setFixed(bool fix)
Set the fixedness of a variable.
Definition: Variable.h:205
std::vector< Variable > weights
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
Observable * sigma
std::vector< PdfBase * > comps
fptype getValue() const
Get the value.
Definition: Variable.h:68

◆ makeOverallSignal()

GooPdf* makeOverallSignal ( )

Definition at line 2652 of file pipipi0DPFit.cpp.

References GooFit::PdfBase::addSpecialMask(), binEffData, createWeightHistogram(), effdata, GooFit::FitManagerMinuit2::fit(), GooFit::FitManagerMinuit1::fit(), foo, foodal, GooFit::Application::get_filename(), GooFit::Observable::getNumBins(), loadDataFile(), m23Slices, makeBkg_sigma_strips(), makeEfficiencyPdf(), makeEfficiencyPoly(), makeKzeroVeto(), makeSigmaHists(), makeSignalPdf(), minuit1, polyEff, GooFit::readFromFile(), saveEffPlot, GooFit::PdfBase::setData(), GooFit::Observable::setNumBins(), GooFit::GooPdf::setParameterConstantness(), underlyingBins, and useHistogramSigma.

Referenced by runCanonicalFit(), and runTruthMCFit().

2652  {
2653  makeKzeroVeto();
2654 
2655  int oldBins1 = m12->getNumBins();
2656  int oldBins2 = m13->getNumBins();
2657  // Too fine a binning here leads to bad results due to fluctuations.
2658  m12->setNumBins(120);
2659  m13->setNumBins(120);
2660  vector<Observable> lvars;
2661  lvars.push_back(*m12);
2662  lvars.push_back(*m13);
2663  binEffData = new BinnedDataSet(lvars);
2665  std::cout << "Loading efficiency data" << std::endl;
2666  std::string fname = app_ptr->get_filename("./dataFiles/efficiency_flat.txt", "examples/pipipi0DPFit");
2667  loadDataFile(fname, &effdata);
2668 
2669  if(saveEffPlot) {
2670  system("mkdir plots_from_mixfit");
2671  foodal->cd();
2672  underlyingBins->Draw("colz");
2673  foodal->SaveAs("plots_from_mixfit/efficiency_bins.png");
2674  foodal->SetLogz(true);
2675  foodal->SaveAs("plots_from_mixfit/efficiency_bins_log.png");
2676  foo->cd();
2677  }
2678 
2679  GooPdf *eff = 0;
2680 
2681  // Polynomial version
2682  if(polyEff)
2683  eff = makeEfficiencyPoly();
2684  // SmoothHistogram version
2685  else
2686  eff = makeEfficiencyPdf();
2687 
2688  eff->setData(effdata);
2689 
2690  if(minuit1) {
2691  GooFit::FitManagerMinuit1 effpdf(eff);
2692  effpdf.fit();
2693  } else {
2694  GooFit::FitManagerMinuit2 effpdf(eff);
2695  effpdf.fit();
2696  }
2697 
2698  eff->setParameterConstantness(true);
2699  delete binEffData;
2700  binEffData = nullptr;
2701  delete effdata;
2702  effdata = nullptr;
2703 
2704  m12->setNumBins(oldBins1);
2705  m13->setNumBins(oldBins2);
2706 
2707  comps.clear();
2708  comps.push_back(eff);
2709  comps.push_back(kzero_veto);
2710  ProdPdf *effWithVeto = new ProdPdf("effWithVeto", comps);
2711 
2712  std::cout << "Creating signal PDF\n";
2713  signalDalitz = makeSignalPdf(0, effWithVeto);
2714 
2715  std::cout << "Creating sigma PDF\n";
2716 
2717  // sig0_jsugg = makeSigmaMap();
2718  // sig0_jsugg = make1BinSigmaMap();
2719  // sig0_jsugg = make4BinSigmaMap();
2720  // sig0_jsugg = makeMikhailJSU_gg();
2721  if(useHistogramSigma)
2723  else
2725 
2726  sig0_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
2727  // sig0_jsugg = makeSignalJSU_gg(-1, false);
2728 
2729  /*
2730  sig0_jsugg->setData(data);
2731  FitManager jsupdf(sig0_jsugg);
2732  gettimeofday(&startTime, nullptr);
2733  jsupdf.fit();
2734  gettimeofday(&stopTime, nullptr);
2735  timersub(&stopTime, &startTime, &totalTime);
2736  std::cout << "Time for sigma fit : " << totalTime.tv_sec + totalTime.tv_usec/1000000.0 << " seconds." << std::endl;
2737  */
2738  sprintf(strbuffer, "signal_sigma_%islices_pdf.txt", m23Slices);
2739  fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
2742 
2743  comps.clear();
2744  comps.push_back(signalDalitz);
2745  comps.push_back(sig0_jsugg);
2746  std::cout << "Creating overall PDF\n";
2747  ProdPdf *overallSignal = new ProdPdf("overallSignal", comps);
2748  return overallSignal;
2749 }
TH2F * underlyingBins
__host__ void setParameterConstantness(bool constant=true)
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
int m23Slices
TddpPdf * makeSignalPdf(MixingTimeResolution *resolution=0, GooPdf *eff=0)
char strbuffer[1000]
GooPdf * makeBkg_sigma_strips(int bkgnum)
void loadDataFile(std::string fname, UnbinnedDataSet **setToFill=0, int effSkip=3)
GooPdf * makeEfficiencyPdf()
Observable * m12
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * effdata
TddpPdf * signalDalitz
BinnedDataSet * binEffData
TCanvas * foodal
bool polyEff
__host__ void setData(DataSet *data)
Observable * m13
std::vector< PdfBase * > comps
void createWeightHistogram()
bool minuit1
bool useHistogramSigma
GooPdf * makeKzeroVeto()
void readFromFile(PdfBase *pdf, const char *fname)
bool saveEffPlot
GooPdf * kzero_veto
TCanvas * foo
GooFit::Application * app_ptr
GooPdf * makeSigmaHists()
GooPdf * sig0_jsugg
__host__ void addSpecialMask(int m)
Definition: PdfBase.h:88
std::string get_filename(const std::string &input_str, std::string base="") const
GooPdf * makeEfficiencyPoly()

◆ makeSigmaHists()

GooPdf* makeSigmaHists ( )

Definition at line 2502 of file pipipi0DPFit.cpp.

References constantZero, dtime, GooFit::Application::get_filename(), jsuList, m23Slices, make_m23_transform(), readWrapper(), and sigma.

Referenced by makeOverallSignal().

2502  {
2503  GooPdf *m23_composite = make_m23_transform();
2504  std::vector<std::unique_ptr<BinnedDataSet>> sigmaHists;
2505 
2506  for(int i = 0; i < m23Slices; ++i)
2507  sigmaHists.emplace_back(new BinnedDataSet(*sigma));
2508 
2509  std::ifstream reader;
2510  std::string fname = app_ptr->get_filename("./dataFiles/signalMC_truth_mm_0.txt", "examples/pipipi0DPFit");
2511  readWrapper(reader, fname);
2512  std::string buffer;
2513 
2514  while(!reader.eof()) {
2515  reader >> buffer;
2516 
2517  if(buffer == "====")
2518  break;
2519 
2520  std::cout << buffer;
2521  }
2522 
2523  double dummy = 0;
2524  double m23 = 0;
2525 
2526  while(!reader.eof()) {
2527  reader >> dummy;
2528  reader >> m23;
2529  reader >> dummy;
2530  reader >> dummy;
2531  reader >> dummy;
2532  reader >> dummy;
2533  reader >> dummy;
2534  reader >> *dtime;
2535  reader >> *sigma;
2536 
2537  reader >> dummy;
2538  reader >> dummy;
2539  reader >> dummy;
2540  reader >> dummy;
2541  reader >> dummy;
2542  reader >> dummy;
2543  reader >> dummy;
2544  reader >> dummy;
2545  reader >> dummy;
2546  reader >> dummy;
2547  reader >> dummy;
2548 
2549  if(dtime->getValue() < dtime->getLowerLimit())
2550  continue;
2551 
2552  if(dtime->getValue() > dtime->getUpperLimit())
2553  continue;
2554 
2555  int bin = (int)floor((m23 / 3.0) * m23Slices);
2556  sigmaHists[bin]->addEvent();
2557  }
2558 
2559  vector<GooPdf *> jsuList;
2560 
2561  for(int i = 0; i < m23Slices; ++i) {
2562  sprintf(strbuffer, "signal_sigma_hist_%i", i);
2563  SmoothHistogramPdf *hist = new SmoothHistogramPdf(strbuffer, sigmaHists[i].get(), constantZero);
2564  jsuList.push_back(hist);
2565  }
2566 
2567  return new MappedPdf("signalSigmaHist", m23_composite, jsuList);
2568 }
vector< GooPdf * > jsuList
Variable constantZero("constantZero", 0)
int m23Slices
char strbuffer[1000]
Observable * sigma
GooPdf * make_m23_transform()
Observable * dtime
bool readWrapper(std::ifstream &reader, std::string fname=strbuffer)
GooFit::Application * app_ptr
std::string get_filename(const std::string &input_str, std::string base="") const

◆ makeSigmaMap()

GooPdf* makeSigmaMap ( )

Definition at line 1257 of file pipipi0DPFit.cpp.

References GooFit::UnbinnedDataSet::addEvent(), GooFit::FitManagerMinuit2::fit(), GooFit::FitManagerMinuit1::fit(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::DataSet::getNumEvents(), GooFit::Indexable::getUpperLimit(), GooFit::UnbinnedDataSet::getValue(), GooFit::Indexable::getValue(), makeSignalJSU_gg(), minuit1, numSigmaBins, GooFit::PdfBase::setData(), GooFit::GooPdf::setParameterConstantness(), and GooFit::Indexable::setValue().

1257  {
1258  sigma_dat_hists = new TH1F *[numSigmaBins];
1259  sigma_pdf_hists = new TH1F *[numSigmaBins];
1261 
1262  std::vector<Observable> vars;
1263  vars.push_back(*sigma);
1264 
1265  for(int i = 0; i < numSigmaBins; ++i) {
1266  sprintf(strbuffer, "sigma_data_hist_%i", i);
1267  sigma_dat_hists[i]
1268  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1269  sigma_dat_hists[i]->SetStats(false);
1270  sigma_dat_hists[i]->SetMarkerStyle(8);
1271  sigma_dat_hists[i]->SetMarkerSize(0.7);
1272  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1273  sigma_pdf_hists[i]
1274  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1275  sigma_pdf_hists[i]->SetStats(false);
1276  sigma_pdf_hists[i]->SetLineWidth(3);
1277  sigma_pdf_hists[i]->SetLineColor(kBlue);
1278 
1279  sigma_data[i] = new UnbinnedDataSet(vars);
1280  }
1281 
1282  int numEvents = data->getNumEvents();
1283 
1284  for(int i = 0; i < numEvents; ++i) {
1285  m12->setValue(data->getValue(*m12, i));
1286  m13->setValue(data->getValue(*m13, i));
1287  sigma->setValue(data->getValue(*sigma, i));
1288 
1289  int xbin = (int)floor(m12->getValue() / 0.5);
1290  int ybin = (int)floor(m13->getValue() / 0.5);
1291  int overallbin = ybin * 6 + xbin;
1292  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1293  sigma_data[overallbin]->addEvent();
1294  }
1295 
1296  // vector<GooPdf*> jsuList;
1297  for(int i = 0; i < numSigmaBins; ++i) {
1298  GooPdf *js = makeSignalJSU_gg(i, false);
1299  jsuList.push_back(js);
1300 
1301  // int xbin = i % 6;
1302  // int ybin = i / 6;
1303 
1304  /*int inDalitzPlot = 0;
1305  if (cpuDalitz(0.5*(xbin+0), 0.5*(ybin+0), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1306  if (cpuDalitz(0.5*(xbin+1), 0.5*(ybin+0), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1307  if (cpuDalitz(0.5*(xbin+0), 0.5*(ybin+1), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1308  if (cpuDalitz(0.5*(xbin+1), 0.5*(ybin+1), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1309  */
1310  if(0 == sigma_data[i]->getNumEvents())
1311  js->setParameterConstantness(true);
1312  else {
1313  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1314  js->setData(sigma_data[i]);
1315  if(minuit1) {
1316  GooFit::FitManagerMinuit1 currpdf(js);
1317  currpdf.fit();
1318  } else {
1319  GooFit::FitManagerMinuit2 currpdf(js);
1320  currpdf.fit();
1321  }
1322  js->setParameterConstantness(true);
1323  // js->clearCurrentFit();
1324  std::cout << "Done with sigma box " << i << "\n";
1325  }
1326  }
1327 
1328  vector<Observable> obses;
1329  obses.push_back(*m12);
1330  obses.push_back(*m13);
1331  vector<double> limits;
1332  limits.push_back(0);
1333  limits.push_back(0);
1334  vector<double> binSizes;
1335  binSizes.push_back(0.5);
1336  binSizes.push_back(0.5);
1337  vector<int> numBins;
1338  numBins.push_back(6);
1339  numBins.push_back(6);
1340  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1341 
1342  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1343 }
size_t getNumEvents() const
Definition: DataSet.h:47
vector< GooPdf * > jsuList
__host__ void setParameterConstantness(bool constant=true)
char strbuffer[1000]
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
TH1F ** sigma_pdf_hists
const int numSigmaBins
Observable * m12
fptype getValue(const Observable &var, size_t idx) const
Get the value at a specific variable and event number.
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * data
__host__ void setData(DataSet *data)
Observable * m13
Observable * sigma
bool minuit1
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
TH1F ** sigma_dat_hists
GooPdf * makeSignalJSU_gg(int idx, bool fixem=true)
UnbinnedDataSet ** sigma_data

◆ makeSignalJSU_gg()

GooPdf* makeSignalJSU_gg ( int  idx,
bool  fixem = true 
)

Definition at line 1123 of file pipipi0DPFit.cpp.

References GooFit::Variable::setFixed(), GooFit::Indexable::setUpperLimit(), and GooFit::Indexable::setValue().

Referenced by make1BinSigmaMap(), make4BinSigmaMap(), and makeSigmaMap().

1123  {
1124  // Values from TSigma, 'Mikhail default'.
1125  static int jsugg_num = -1;
1126  jsugg_num++;
1127 
1128  sprintf(strbuffer, "js_meana_%i", jsugg_num);
1129  Variable js_meana(strbuffer, 0.0593, 0.01, 0, 0.2);
1130  js_meana.setFixed(fixem);
1131  sprintf(strbuffer, "js_sigma_%i", jsugg_num);
1132  Variable js_sigma(strbuffer, 0.000474, 0.0001, 0, 0.001);
1133  js_sigma.setFixed(fixem);
1134  sprintf(strbuffer, "js_gamma_%i", jsugg_num);
1135  Variable js_gamma(strbuffer, -10.1942, 1, -30, 0);
1136  js_gamma.setFixed(fixem);
1137  sprintf(strbuffer, "js_delta_%i", jsugg_num);
1138  Variable js_delta(strbuffer, 1.4907, 0.1, 0.5, 5);
1139  js_delta.setFixed(fixem);
1140  sprintf(strbuffer, "frac_jsu_%i", jsugg_num);
1141  Variable frac_jsu(strbuffer, 0.9516, 0.01, 0.5, 1.0);
1142  frac_jsu.setFixed(fixem);
1143  sprintf(strbuffer, "frac_ga1_%i", jsugg_num);
1144  Variable frac_ga1(strbuffer, 0.001845, 0.0005, 0.0001, 0.3);
1145  frac_ga1.setFixed(fixem);
1146  sprintf(strbuffer, "g1_meana_%i", jsugg_num);
1147  Variable g1_meana(strbuffer, 0.2578, 0.003, 0.1, 0.5);
1148  g1_meana.setFixed(fixem);
1149  sprintf(strbuffer, "g1_sigma_%i", jsugg_num);
1150  Variable g1_sigma(strbuffer, 0.03086, 0.01, 0.005, 0.25);
1151  g1_sigma.setFixed(fixem);
1152  sprintf(strbuffer, "g2_meana_%i", jsugg_num);
1153  Variable g2_meana(strbuffer, 0.32, 0.01, 0.1, 0.5);
1154  g2_meana.setFixed(fixem);
1155  sprintf(strbuffer, "g2_sigma_%i", jsugg_num);
1156  Variable g2_sigma(strbuffer, 0.05825, 0.01, 0.01, 0.1);
1157  g2_sigma.setFixed(fixem);
1158  // Variable* g2_sigma = new Variable("g2_sigma", 0.5825);
1159 
1160  sprintf(strbuffer, "js_%i", jsugg_num);
1161  JohnsonSUPdf *js = new JohnsonSUPdf(strbuffer, *sigma, js_meana, js_sigma, js_gamma, js_delta);
1162  sprintf(strbuffer, "g1_%i", jsugg_num);
1163  GaussianPdf *g1 = new GaussianPdf(strbuffer, *sigma, g1_meana, g1_sigma);
1164  sprintf(strbuffer, "g2_%i", jsugg_num);
1165  // GaussianPdf* g2 = new GaussianPdf(strbuffer, sigma, g2_meana, g2_sigma);
1166 
1167  weights.clear();
1168  weights.push_back(frac_jsu);
1169  // weights.push_back(frac_ga1);
1170  comps.clear();
1171  comps.push_back(js);
1172  comps.push_back(g1);
1173  // comps.push_back(g2);
1174 
1175  // Deal with special indices to get different starting points
1176  switch(idx) {
1177  case 1:
1178  g2_sigma.setUpperLimit(0.15);
1179  break;
1180 
1181  case 2:
1182  frac_jsu.setValue(0.80);
1183  break;
1184 
1185  // case 5:
1186  // g1_sigma->getLowerLimit() = 0.005;
1187  // break;
1188  case 7:
1189  frac_jsu.setValue(0.80);
1190  // frac_ga1->getValue() = 0.05;
1191  break;
1192 
1193  case 11:
1194  frac_ga1.setUpperLimit(0.4);
1195  frac_jsu.setValue(0.80);
1196  break;
1197 
1198  default:
1199  break;
1200  }
1201 
1202  sprintf(strbuffer, "signal_sigma_%i", jsugg_num);
1203  AddPdf *signal_sigma = new AddPdf(strbuffer, weights, comps);
1204  return signal_sigma;
1205 }
char strbuffer[1000]
std::vector< Variable > weights
Observable * sigma
std::vector< PdfBase * > comps

◆ makeSignalPdf()

TddpPdf* makeSignalPdf ( MixingTimeResolution resolution = 0,
GooPdf eff = 0 
)

Definition at line 520 of file pipipi0DPFit.cpp.

References _mD0, _mD02, _mD02inv, constantOne, constantZero, GooFit::DecayInfo3::daug1Mass, GooFit::DecayInfo3::daug2Mass, GooFit::DecayInfo3::daug3Mass, drop_f0_1370, drop_f0_1500, drop_f0_1710, drop_f0_600, drop_f0_980, drop_f2_1270, drop_rho_1450, drop_rho_1700, dtop0pp, fitMasses, fixedRhoMass, fixedRhoWidth, mdslices, GooFit::DecayInfo3::meson_radius, mesonRad, GooFit::DecayInfo3::motherMass, GooFit::PAIR_12, GooFit::PAIR_13, GooFit::PAIR_23, piPlusMass, piZeroMass, GooFit::DecayInfo3::resonances, GooFit::Variable::setFixed(), and GooFit::Indexable::setValue().

Referenced by makeOverallSignal(), runGeneratedMCFit(), and runToyFit().

520  {
526 
527  // dtop0pp._tau = new Variable("tau", 0.4116, 0.001, 0.300, 0.500);
528  // Setting limits causes trouble with MNHESSE - why is this?
529  // dtop0pp._xmixing = new Variable("xmixing", 0.01, 0.001, -0.05, 0.05);
530  // dtop0pp._ymixing = new Variable("ymixing", 0.01, 0.001, -0.05, 0.05);
531 
532  // dtop0pp._tau.fixed = true;
533  // dtop0pp._xmixing.fixed = true;
534  // dtop0pp._ymixing.fixed = true;
535 
536  ResonancePdf *rhop = new Resonances::RBW(
537  "rhop", Variable("rhop_amp_real", 1), Variable("rhop_amp_imag", 0), fixedRhoMass, fixedRhoWidth, 1, PAIR_12);
538 
539  bool fixAmps = false;
540 
541  ResonancePdf *rhom = new Resonances::RBW(
542  "rhom",
543  fixAmps ? Variable("rhom_amp_real", 0.714) : Variable("rhom_amp_real", 0.714, 0.001, 0, 0),
544  fixAmps ? Variable("rhom_amp_imag", -0.025) : Variable("rhom_amp_imag", -0.025, 0.1, 0, 0),
545  fixedRhoMass,
547  1,
548  PAIR_13);
549 
550  ResonancePdf *rho0
551  = new Resonances::GS("rho0",
552  fixAmps ? Variable("rho0_amp_real", 0.565) : Variable("rho0_amp_real", 0.565, 0.001, 0, 0),
553  fixAmps ? Variable("rho0_amp_imag", 0.164) : Variable("rho0_amp_imag", 0.164, 0.1, 0, 0),
554  fixedRhoMass,
556  1,
557  PAIR_23);
558 
559  Variable rho1450Mass("rhop_1450_mass", 1.465, 0.01, 1.0, 2.0);
560  Variable rho1450Width("rhop_1450_width", 0.400, 0.01, 0.01, 5.0);
561 
562  ResonancePdf *rhop_1450 = new Resonances::RBW(
563  "rhop_1450",
564  fixAmps ? Variable("rhop_1450_amp_real", -0.174) : Variable("rhop_1450_amp_real", -0.174, 0.001, 0, 0),
565  fixAmps ? Variable("rhop_1450_amp_imag", -0.117) : Variable("rhop_1450_amp_imag", -0.117, 0.1, 0, 0),
566  rho1450Mass,
567  rho1450Width,
568  1,
569  PAIR_12);
570 
571  ResonancePdf *rho0_1450 = new Resonances::RBW(
572  "rho0_1450",
573  fixAmps ? Variable("rho0_1450_amp_real", 0.325) : Variable("rho0_1450_amp_real", 0.325, 0.001, 0, 0),
574  fixAmps ? Variable("rho0_1450_amp_imag", 0.057) : Variable("rho0_1450_amp_imag", 0.057, 0.1, 0, 0),
575  rho1450Mass,
576  rho1450Width,
577  1,
578  PAIR_23);
579 
580  ResonancePdf *rhom_1450 = new Resonances::RBW(
581  "rhom_1450",
582  fixAmps ? Variable("rhom_1450_amp_real", 0.788) : Variable("rhom_1450_amp_real", 0.788, 0.001, 0, 0),
583  fixAmps ? Variable("rhom_1450_amp_imag", 0.226) : Variable("rhom_1450_amp_imag", 0.226, 0.1, 0, 0),
584  rho1450Mass,
585  rho1450Width,
586  1,
587  PAIR_13);
588 
589  Variable rho1700Mass("rhop_1700_mass", 1.720, 0.01, 1.6, 1.9);
590  Variable rho1700Width("rhop_1700_width", 0.250, 0.01, 0.1, 1.0);
591 
592  ResonancePdf *rhop_1700 = new Resonances::RBW(
593  "rhop_1700",
594  fixAmps ? Variable("rhop_1700_amp_real", 2.151) : Variable("rhop_1700_amp_real", 2.151, 0.001, 0, 0),
595  fixAmps ? Variable("rhop_1700_amp_imag", -0.658) : Variable("rhop_1700_amp_imag", -0.658, 0.1, 0, 0),
596  rho1700Mass,
597  rho1700Width,
598  1,
599  PAIR_12);
600 
601  ResonancePdf *rho0_1700 = new Resonances::RBW(
602  "rho0_1700",
603  fixAmps ? Variable("rho0_1700_amp_real", 2.400) : Variable("rho0_1700_amp_real", 2.400, 0.001, 0, 0),
604  fixAmps ? Variable("rho0_1700_amp_imag", -0.734) : Variable("rho0_1700_amp_imag", -0.734, 0.1, 0, 0),
605  rho1700Mass,
606  rho1700Width,
607  1,
608  PAIR_23);
609 
610  ResonancePdf *rhom_1700 = new Resonances::RBW(
611  "rhom_1700",
612  fixAmps ? Variable("rhom_1700_amp_real", 1.286) : Variable("rhom_1700_amp_real", 1.286, 0.001, 0, 0),
613  fixAmps ? Variable("rhom_1700_amp_imag", -1.532) : Variable("rhom_1700_amp_imag", -1.532, 0.1, 0, 0),
614  rho1700Mass,
615  rho1700Width,
616  1,
617  PAIR_13);
618 
619  ResonancePdf *f0_980 = new Resonances::RBW("f0_980",
620  fixAmps ? Variable("f0_980_amp_real", 0.008 * (-_mD02))
621  : Variable("f0_980_amp_real", 0.008 * (-_mD02), 0.001, 0, 0),
622  fixAmps ? Variable("f0_980_amp_imag", -0.013 * (-_mD02))
623  : Variable("f0_980_amp_imag", -0.013 * (-_mD02), 0.1, 0, 0),
624  Variable("f0_980_mass", 0.980, 0.01, 0.8, 1.2),
625  Variable("f0_980_width", 0.044, 0.001, 0.001, 0.08),
626  0,
627  PAIR_23);
628 
629  ResonancePdf *f0_1370 = new Resonances::RBW("f0_1370",
630  fixAmps ? Variable("f0_1370_amp_real", -0.058 * (-_mD02))
631  : Variable("f0_1370_amp_real", -0.058 * (-_mD02), 0.001, 0, 0),
632  fixAmps ? Variable("f0_1370_amp_imag", 0.026 * (-_mD02))
633  : Variable("f0_1370_amp_imag", 0.026 * (-_mD02), 0.1, 0, 0),
634  Variable("f0_1370_mass", 1.434, 0.01, 1.2, 1.6),
635  Variable("f0_1370_width", 0.173, 0.01, 0.01, 0.4),
636  0,
637  PAIR_23);
638 
639  ResonancePdf *f0_1500 = new Resonances::RBW("f0_1500",
640  fixAmps ? Variable("f0_1500_amp_real", 0.057 * (-_mD02))
641  : Variable("f0_1500_amp_real", 0.057 * (-_mD02), 0.001, 0, 0),
642  fixAmps ? Variable("f0_1500_amp_imag", 0.012 * (-_mD02))
643  : Variable("f0_1500_amp_imag", 0.012 * (-_mD02), 0.1, 0, 0),
644  Variable("f0_1500_mass", 1.507, 0.01, 1.3, 1.7),
645  Variable("f0_1500_width", 0.109, 0.01, 0.01, 0.3),
646  0,
647  PAIR_23);
648 
649  ResonancePdf *f0_1710 = new Resonances::RBW("f0_1710",
650  fixAmps ? Variable("f0_1710_amp_real", 0.070 * (-_mD02))
651  : Variable("f0_1710_amp_real", 0.070 * (-_mD02), 0.001, 0, 0),
652  fixAmps ? Variable("f0_1710_amp_imag", 0.087 * (-_mD02))
653  : Variable("f0_1710_amp_imag", 0.087 * (-_mD02), 0.1, 0, 0),
654  Variable("f0_1710_mass", 1.714, 0.01, 1.5, 2.9),
655  Variable("f0_1710_width", 0.140, 0.01, 0.01, 0.5),
656  0,
657  PAIR_23);
658 
659  ResonancePdf *f2_1270
660  = new Resonances::RBW("f2_1270",
661  fixAmps ? Variable("f2_1270_amp_real", -1.027 * (-_mD02inv))
662  : Variable("f2_1270_amp_real", -1.027 * (-_mD02inv), 0.001, 0, 0),
663  fixAmps ? Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv))
664  : Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv), 0.1, 0, 0),
665  Variable("f2_1270_mass", 1.2754, 0.01, 1.0, 1.5),
666  Variable("f2_1270_width", 0.1851, 0.01, 0.01, 0.4),
667  2,
668  PAIR_23);
669 
670  ResonancePdf *f0_600 = new Resonances::RBW("f0_600",
671  fixAmps ? Variable("f0_600_amp_real", 0.068 * (-_mD02))
672  : Variable("f0_600_amp_real", 0.068 * (-_mD02), 0.001, 0, 0),
673  fixAmps ? Variable("f0_600_amp_imag", 0.010 * (-_mD02))
674  : Variable("f0_600_amp_imag", 0.010 * (-_mD02), 0.1, 0, 0),
675  Variable("f0_600_mass", 0.500, 0.01, 0.3, 0.7),
676  Variable("f0_600_width", 0.400, 0.01, 0.2, 0.6),
677  0,
678  PAIR_23);
679 
680  ResonancePdf *nonr = new Resonances::NonRes(
681  "nonr",
682  fixAmps ? Variable("nonr_amp_real", 0.5595 * (-1)) : Variable("nonr_amp_real", 0.5595 * (-1), 0.001, 0, 0),
683  fixAmps ? Variable("nonr_amp_imag", -0.108761 * (-1)) : Variable("nonr_amp_imag", -0.108761 * (-1), 0.1, 0, 0));
684 
685  dtop0pp.resonances.push_back(nonr);