diff --git a/src/FitBase/StackBase.cxx b/src/FitBase/StackBase.cxx index 1a721a1..bd6aee0 100644 --- a/src/FitBase/StackBase.cxx +++ b/src/FitBase/StackBase.cxx @@ -1,216 +1,250 @@ #include "StackBase.h" void StackBase::AddMode(std::string name, std::string title, int linecolor, int linewidth, int fillstyle) { // int ncur = fAllLabels.size(); fAllLabels.push_back(name); fAllTitles.push_back(title); std::vector temp; temp.push_back(linecolor); temp.push_back(linewidth); temp.push_back(fillstyle); fAllStyles.push_back(temp); } void StackBase::FluxUnfold(TH1D *flux, TH1D *events, double scalefactor, int nevents) { for (size_t i = 0; i < fAllLabels.size(); i++) { if (fNDim == 1) { PlotUtils::FluxUnfoldedScaling((TH1D *)fAllHists[i], flux, events, scalefactor, nevents); } else if (fNDim == 2) { PlotUtils::FluxUnfoldedScaling((TH2D *)fAllHists[i], flux, events, scalefactor); } } } void StackBase::AddMode(int index, std::string name, std::string title, int linecolor, int linewidth, int fillstyle) { while (fAllLabels.size() <= (UInt_t)index) { fAllLabels.push_back(""); fAllTitles.push_back(""); fAllStyles.push_back(std::vector(1, 1)); } fAllLabels[index] = (name); fAllTitles[index] = (title); std::vector temp; temp.push_back(linecolor); temp.push_back(linewidth); temp.push_back(fillstyle); fAllStyles[index] = temp; } bool StackBase::IncludeInStack(TH1 *hist) { if (!FitPar::Config().GetParB("includeemptystackhists") and hist->Integral() == 0.0) return false; return true; } bool StackBase::IncludeInStack(int index) { return true; } void StackBase::SetupStack(TH1 *hist) { fTemplate = (TH1 *)hist->Clone(fName.c_str()); fTemplate->Reset(); // Determine template dim fNDim = fTemplate->GetDimension(); for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists.push_back( (TH1 *)fTemplate->Clone((fName + "_" + fAllLabels[i]).c_str())); } }; void StackBase::Scale(double sf, std::string opt) { for (size_t i = 0; i < fAllLabels.size(); i++) { // std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl; fAllHists[i]->Scale(sf, opt.c_str()); } }; void StackBase::Reset() { for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists[i]->Reset(); } }; void StackBase::FillStack(int index, double x, double y, double z, double weight) { if (index < 0 or (UInt_t) index >= fAllLabels.size()) { NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " " - << fAllLabels.size()); + << fAllLabels.size()); return; } if (fNDim == 1) fAllHists[index]->Fill(x, y); else if (fNDim == 2) { // std::cout << "Filling 2D Stack " << index << " " << x << " " << y << " " // << z << std::endl; ((TH2 *)fAllHists[index])->Fill(x, y, z); } else if (fNDim == 3) ((TH3 *)fAllHists[index])->Fill(x, y, z, weight); } +void StackBase::SetBinContentStack(int index, int binx, int biny, int binz, + double content) { + if (index < 0 or (UInt_t) index >= fAllLabels.size()) { + NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " " + << fAllLabels.size()); + return; + } + + if (fNDim == 1) { + fAllHists[index]->SetBinContent(binx, content); + } else if (fNDim == 2) { + ((TH2 *)fAllHists[index])->SetBinContent(binx, biny, content); + } else if (fNDim == 3) { + ((TH3 *)fAllHists[index])->SetBinContent(binx, biny, binz, content); + } +} + +void StackBase::SetBinErrorStack(int index, int binx, int biny, int binz, + double error) { + if (index < 0 or (UInt_t) index >= fAllLabels.size()) { + NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " " + << fAllLabels.size()); + return; + } + + if (fNDim == 1) { + fAllHists[index]->SetBinError(binx, error); + } else if (fNDim == 2) { + ((TH2 *)fAllHists[index])->SetBinError(binx, biny, error); + } else if (fNDim == 3) { + ((TH3 *)fAllHists[index])->SetBinError(binx, biny, binz, error); + } +} + void StackBase::Write() { THStack *st = new THStack(); // Loop and add all histograms bool saveseperate = FitPar::Config().GetParB("WriteSeperateStacks"); for (size_t i = 0; i < fAllLabels.size(); i++) { if (!IncludeInStack(fAllHists[i])) continue; if (!IncludeInStack(i)) continue; fAllHists[i]->SetTitle(fAllTitles[i].c_str()); fAllHists[i]->GetXaxis()->SetTitle(fXTitle.c_str()); fAllHists[i]->GetYaxis()->SetTitle(fYTitle.c_str()); fAllHists[i]->GetZaxis()->SetTitle(fZTitle.c_str()); fAllHists[i]->SetLineColor(fAllStyles[i][0]); fAllHists[i]->SetLineWidth(fAllStyles[i][1]); fAllHists[i]->SetFillStyle(fAllStyles[i][2]); fAllHists[i]->SetFillColor(fAllStyles[i][0]); if (saveseperate) fAllHists[i]->Write(); st->Add(fAllHists[i]); } st->SetTitle(fTitle.c_str()); st->SetName(fName.c_str()); st->Write(); delete st; }; void StackBase::Multiply(TH1 *hist) { for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists[i]->Multiply(hist); } } void StackBase::Divide(TH1 *hist) { for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists[i]->Divide(hist); } } void StackBase::Add(TH1 *hist, double scale) { for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists[i]->Add(hist, scale); } } void StackBase::Add(StackBase *hist, double scale) { if (hist->GetType() != fType) { NUIS_ERR(WRN, "Trying to add two StackBases of different types!"); NUIS_ERR(WRN, fType << " + " << hist->GetType() << " = Undefined."); NUIS_ERR(WRN, "Doing nothing..."); return; } for (size_t i = 0; i < fAllLabels.size(); i++) { fAllHists[i]->Add(hist->GetHist(i)); } } TH1 *StackBase::GetHist(int entry) { return fAllHists[entry]; } TH1 *StackBase::GetHist(std::string label) { TH1 *hist = NULL; std::vector splitlabels = GeneralUtils::ParseToStr(label, "+"); for (size_t j = 0; j < splitlabels.size(); j++) { std::string newlabel = splitlabels[j]; for (size_t i = 0; i < fAllLabels.size(); i++) { if (newlabel == fAllLabels[i]) { if (!hist) hist = (TH1 *)fAllHists[i]->Clone(); else hist->Add(fAllHists[i]); } } } return hist; } THStack StackBase::GetStack() { THStack st = THStack(); for (size_t i = 0; i < fAllLabels.size(); i++) { st.Add(fAllHists[i]); } return st; } void StackBase::AddNewHist(std::string name, TH1 *hist) { AddMode(fAllLabels.size(), name, hist->GetTitle(), hist->GetLineColor()); fAllHists.push_back((TH1 *)hist->Clone()); } void StackBase::AddToCategory(std::string name, TH1 *hist) { for (size_t i = 0; i < fAllLabels.size(); i++) { if (name == fAllLabels[i]) { fAllHists[i]->Add(hist); break; } } } void StackBase::AddToCategory(int index, TH1 *hist) { fAllHists[index]->Add(hist); } diff --git a/src/FitBase/StackBase.h b/src/FitBase/StackBase.h index c227813..a2bfe1e 100644 --- a/src/FitBase/StackBase.h +++ b/src/FitBase/StackBase.h @@ -1,126 +1,129 @@ #ifndef STACK_BASE_H #define STACK_BASE_H -#include "MeasurementVariableBox.h" #include "FitLogger.h" #include "GeneralUtils.h" +#include "MeasurementVariableBox.h" #include "TH1.h" #include "TH1D.h" -#include "TH2D.h" -#include "THStack.h" #include "TH2.h" +#include "TH2D.h" #include "TH3.h" +#include "THStack.h" #include "PlotUtils.h" class StackBase { public: - StackBase() {}; - ~StackBase() {}; - - virtual void AddMode(std::string name, std::string title, - int linecolor = 1, int linewidth = 1, int fillstyle = 1001); - virtual void AddMode(int index, std::string name, std::string title, - int linecolor = 1, int linewidth = 1, int fillstyle = 1001); - - virtual bool IncludeInStack(TH1* hist); - virtual bool IncludeInStack(int index); - - virtual void SetupStack(TH1* hist); - virtual void Scale(double sf, std::string opt = ""); - virtual void FluxUnfold(TH1D* flux, TH1D* events, double scalefactor, int nevents); - virtual void Reset(); - virtual void FillStack(int index, double x, double y = 1.0, double z = 1.0, double weight = 1.0); - virtual void Write(); - - virtual void Add(StackBase* hist, double scale); - virtual void Add(TH1* hist, double scale); - - virtual void AddNewHist(std::string name, TH1* hist); - virtual void AddToCategory(std::string name, TH1* hist); - virtual void AddToCategory(int index, TH1* hist); - - virtual void Divide(TH1* hist); - virtual void Multiply(TH1* hist); - virtual TH1* GetHist(int entry); - virtual TH1* GetHist(std::string label); - virtual THStack GetStack(); - - std::string GetType(){return fType;}; - - std::string fName; - std::string fTitle; - std::string fXTitle; - std::string fYTitle; - std::string fZTitle; - std::string fType; - - TH1* fTemplate; - int fNDim; - - // Maps incase we want to be selective - std::vector< std::vector > fAllStyles; - std::vector fAllTitles; - std::vector fAllLabels; - std::vector fAllHists; + StackBase(){}; + ~StackBase(){}; + + virtual void AddMode(std::string name, std::string title, int linecolor = 1, + int linewidth = 1, int fillstyle = 1001); + virtual void AddMode(int index, std::string name, std::string title, + int linecolor = 1, int linewidth = 1, + int fillstyle = 1001); + + virtual bool IncludeInStack(TH1 *hist); + virtual bool IncludeInStack(int index); + + virtual void SetupStack(TH1 *hist); + virtual void Scale(double sf, std::string opt = ""); + virtual void FluxUnfold(TH1D *flux, TH1D *events, double scalefactor, + int nevents); + virtual void Reset(); + virtual void FillStack(int index, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); + virtual void SetBinContentStack(int index, int binx, int biny, int binz, + double content); + virtual void SetBinErrorStack(int index, int binx, int biny, int binz, + double error); + virtual void Write(); + + virtual void Add(StackBase *hist, double scale); + virtual void Add(TH1 *hist, double scale); + + virtual void AddNewHist(std::string name, TH1 *hist); + virtual void AddToCategory(std::string name, TH1 *hist); + virtual void AddToCategory(int index, TH1 *hist); + + virtual void Divide(TH1 *hist); + virtual void Multiply(TH1 *hist); + virtual TH1 *GetHist(int entry); + virtual TH1 *GetHist(std::string label); + virtual THStack GetStack(); + + std::string GetType() { return fType; }; + + std::string fName; + std::string fTitle; + std::string fXTitle; + std::string fYTitle; + std::string fZTitle; + std::string fType; + + TH1 *fTemplate; + int fNDim; + + // Maps incase we want to be selective + std::vector > fAllStyles; + std::vector fAllTitles; + std::vector fAllLabels; + std::vector fAllHists; }; - /* class NuSpeciesStack : public StackBase { public: - SetupStack(TH1* hist) { - AddMode("numu", "numu", kBlue, 2, 3004); - AddMode("numubar", "numubar", kRed, 2, 3004 ); - AddMode("nue", "nue", kGreen, 2, 3004 ); - StackBase::SetupStack(hist); - }; - - void NuSpeciesStack::FillStack(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0) { - Stackbase::FillStack(ConvertSpeciesToIndex(species), x, y, z, weight); - } - - int ConvertSpeciesToIndex(int species) { - switch (species) { - case 14: return 0; - case -14: return 1; - case 11: return 2; - default: return -1; - } - }; + SetupStack(TH1* hist) { + AddMode("numu", "numu", kBlue, 2, 3004); + AddMode("numubar", "numubar", kRed, 2, 3004 ); + AddMode("nue", "nue", kGreen, 2, 3004 ); + StackBase::SetupStack(hist); + }; + + void NuSpeciesStack::FillStack(int species, double x, double y = 1.0, +double z = 1.0, double weight = 1.0) { + Stackbase::FillStack(ConvertSpeciesToIndex(species), x, y, z, +weight); + } + + int ConvertSpeciesToIndex(int species) { + switch (species) { + case 14: return 0; + case -14: return 1; + case 11: return 2; + default: return -1; + } + }; }; class TargetStack : public StackBase { public: - SetupStack(TH1* hist) { - AddMode("C", "C", kBlue, 2, 3004); - AddMode("H", "H", kRed, 2, 3004 ); - AddMode("O", "O", kGreen, 2, 3004 ); - StackBase::SetupStack(hist); - }; - - void NuSpeciesStack::FillStack(int species, double x, - double y = 1.0, double z = 1.0, - double weight = 1.0) { - Stackbase::FillStack(ConvertTargetToIndex(species), x, y, z, weight); - } - - int ConvertTargetToIndex(int target) { - switch (species) { - case 1000010010: return 0; - case 1000: return 1; - case 1000: return 2; - default: return -1; - } - }; + SetupStack(TH1* hist) { + AddMode("C", "C", kBlue, 2, 3004); + AddMode("H", "H", kRed, 2, 3004 ); + AddMode("O", "O", kGreen, 2, 3004 ); + StackBase::SetupStack(hist); + }; + + void NuSpeciesStack::FillStack(int species, double x, + double y = 1.0, double z = 1.0, + double weight = 1.0) { + Stackbase::FillStack(ConvertTargetToIndex(species), x, y, z, +weight); + } + + int ConvertTargetToIndex(int target) { + switch (species) { + case 1000010010: return 0; + case 1000: return 1; + case 1000: return 2; + default: return -1; + } + }; } */ - #endif - - - - - diff --git a/src/FitBase/StandardStacks.cxx b/src/FitBase/StandardStacks.cxx index da2c562..6088382 100644 --- a/src/FitBase/StandardStacks.cxx +++ b/src/FitBase/StandardStacks.cxx @@ -1,428 +1,527 @@ #include "StandardStacks.h" /// Fake stack functions -FakeStack::FakeStack(TH1* hist) { - fTemplate = (TH1*)hist; +FakeStack::FakeStack(TH1 *hist) { + fTemplate = (TH1 *)hist; fNDim = fTemplate->GetDimension(); fFakeType = "TH1"; fTGraphObject = NULL; fTF1Object = NULL; } - -FakeStack::FakeStack(TF1* f){ +FakeStack::FakeStack(TF1 *f) { fTemplate = NULL; fNDim = 1; fTF1Object = f; fFakeType = "TF1"; fTGraphObject = NULL; } -FakeStack::FakeStack(TGraph* gr){ +FakeStack::FakeStack(TGraph *gr) { fTemplate = NULL; fNDim = 1; fTGraphObject = gr; fFakeType = "TGRAPH"; fTF1Object = NULL; } FakeStack::~FakeStack() { fTemplate = NULL; fNDim = 0; fFakeType = "NULL"; fTF1Object = NULL; fTGraphObject = NULL; } void FakeStack::Fill(double x, double y, double z, double weight) { - if (fTemplate){ - if (fNDim == 1) fTemplate->Fill(x, y); - else if (fNDim == 2) ((TH2*)fTemplate)->Fill(x, y, z); - else if (fNDim == 3) ((TH3*)fTemplate)->Fill(x, y, z, weight); + if (fTemplate) { + if (fNDim == 1) + fTemplate->Fill(x, y); + else if (fNDim == 2) + ((TH2 *)fTemplate)->Fill(x, y, z); + else if (fNDim == 3) + ((TH3 *)fTemplate)->Fill(x, y, z, weight); } - if (fTGraphObject){ - fTGraphObject->SetPoint( fTGraphObject->GetN(), x, y); + if (fTGraphObject) { + fTGraphObject->SetPoint(fTGraphObject->GetN(), x, y); } } void FakeStack::Scale(double norm, std::string opt) { - if (fTemplate){ fTemplate->Scale(norm, opt.c_str()); } - if (fTGraphObject){ - for (int i = 0; i < fTGraphObject->GetN(); i++){ - const double* x = fTGraphObject->GetX(); - const double* y = fTGraphObject->GetY(); + if (fTemplate) { + fTemplate->Scale(norm, opt.c_str()); + } + if (fTGraphObject) { + for (int i = 0; i < fTGraphObject->GetN(); i++) { + const double *x = fTGraphObject->GetX(); + const double *y = fTGraphObject->GetY(); fTGraphObject->SetPoint(i, x[i], y[i] * norm); } } } void FakeStack::Reset() { - if (fTemplate){ fTemplate->Reset(); } + if (fTemplate) { + fTemplate->Reset(); + } } void FakeStack::Write() { - if (fTemplate){ fTemplate->Write(); } - if (fTF1Object){ fTF1Object->Write(); } - if (fTGraphObject){ fTGraphObject->Write(); } + if (fTemplate) { + fTemplate->Write(); + } + if (fTF1Object) { + fTF1Object->Write(); + } + if (fTGraphObject) { + fTGraphObject->Write(); + } } - - /// TrueModeStack Functions -TrueModeStack::TrueModeStack(std::string name, std::string title, TH1* hist) { +TrueModeStack::TrueModeStack(std::string name, std::string title, TH1 *hist) { fName = name; fTitle = title; // CC - AddMode(0, "CCQE", "CCQE", kBlue, 2, 1001); - AddMode(1, "CC2p2h", "2p2h", kRed, 2, 1001); - AddMode(2, "CC1piponp", "CC1#pi^{+} on p", kGreen, 2, 1001); - AddMode(3, "CC1pi0onn", "CC1#pi^{0} on n", kGreen + 3, 2, 1001); - AddMode(4, "CC1piponn", "CC1#pi^{+} on n", kGreen - 2, 2, 1001); - AddMode(5, "CCcoh", "CC coherent", kBlue, 2, 1001); - AddMode(6, "CC1gamma", "CC1#gamma", kMagenta, 2, 1001); + AddMode(0, "CCQE", "CCQE", kBlue, 2, 1001); + AddMode(1, "CC2p2h", "2p2h", kRed, 2, 1001); + AddMode(2, "CC1piponp", "CC1#pi^{+} on p", kGreen, 2, 1001); + AddMode(3, "CC1pi0onn", "CC1#pi^{0} on n", kGreen + 3, 2, 1001); + AddMode(4, "CC1piponn", "CC1#pi^{+} on n", kGreen - 2, 2, 1001); + AddMode(5, "CCcoh", "CC coherent", kBlue, 2, 1001); + AddMode(6, "CC1gamma", "CC1#gamma", kMagenta, 2, 1001); AddMode(7, "CCMultipi", "Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001); - AddMode(8, "CC1eta", "CC1#eta^{0} on n", kYellow - 2, 2, 1001); - AddMode(9, "CC1lamkp", "CC1#Lambda1K^{+}", kYellow - 6, 2, 1001); - AddMode(10, "CCDIS", "DIS (W > 2.0)", kRed, 2, 1001); + AddMode(8, "CC1eta", "CC1#eta^{0} on n", kYellow - 2, 2, 1001); + AddMode(9, "CC1lamkp", "CC1#Lambda1K^{+}", kYellow - 6, 2, 1001); + AddMode(10, "CCDIS", "DIS (W > 2.0)", kRed, 2, 1001); // NC - AddMode(11, "NC1pi0onn", "NC1#pi^{0} on n", kBlue, 2, 3004); + AddMode(11, "NC1pi0onn", "NC1#pi^{0} on n", kBlue, 2, 3004); AddMode(12, "NC1pi0onp", "NC1#pi^{0} on p", kBlue + 3, 2, 3004); - AddMode(13, "NC1pimonn", "NC1#pi^{-} on n", kBlue - 2, 2, 3004); + AddMode(13, "NC1pimonn", "NC1#pi^{-} on n", kBlue - 2, 2, 3004); AddMode(14, "NC1piponp", "NC1#pi^{+} on p", kBlue - 8, 2, 3004); - AddMode(15, "NCcoh", "NC Coherent", kBlue + 8, 2, 3004); - AddMode(16, "NC1gammaonn", "NC1#gamma on n", kMagenta, 2, 3004); - AddMode(17, "NC1gammaonp", "NC1#gamma on p", kMagenta - 10, 2, 3004); + AddMode(15, "NCcoh", "NC Coherent", kBlue + 8, 2, 3004); + AddMode(16, "NC1gammaonn", "NC1#gamma on n", kMagenta, 2, 3004); + AddMode(17, "NC1gammaonp", "NC1#gamma on p", kMagenta - 10, 2, 3004); AddMode(18, "NCMultipi", "Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004); AddMode(19, "NC1etaonn", "NC1#eta^{0} on n", kYellow - 2, 2, 3004); AddMode(20, "NC1etaonp", "NC1#eta^{0} on p", kYellow - 4, 2, 3004); - AddMode(21, "NC1kamk0", "NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004); - AddMode(22, "NC1lamkp", "NC1#Lambda1K^{+}", kYellow - 10, 2, 3004); + AddMode(21, "NC1kamk0", "NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004); + AddMode(22, "NC1lamkp", "NC1#Lambda1K^{+}", kYellow - 10, 2, 3004); AddMode(23, "NCDIS", "DIS (W > 2.0)", kRed, 2, 3004); AddMode(24, "NCELonp", "NCEL on p", kBlack, 2, 3004); AddMode(25, "NCELonn", "NCEL on n", kGray, 2, 3004); - AddMode(26, "NC2p2h", "NC 2p2h", kRed+1, 2, 3004); + AddMode(26, "NC2p2h", "NC 2p2h", kRed + 1, 2, 3004); // Undefined AddMode(27, "UNDEFINED", "Undefined", kRed + 1, 2, 3000); StackBase::SetupStack(hist); }; int TrueModeStack::ConvertModeToIndex(int mode) { // std::cout << "Converting Mode " << (mode) << std::endl; switch (abs(mode)) { - case 1: return 0; // CCQE - case 2: return 1; // CC2p2h - case 11: return 2; // CC1piponp - case 12: return 3; // CC1pi0onn - case 13: return 4; // CC1piponn - case 16: return 5; // CCcoh - case 17: return 6; // CC1gamma - case 21: return 7; // CCMultipi - case 22: return 8; // CC1eta - case 23: return 9; // CC1lamkp - case 26: return 10; // CCDIS - - case 31: return 11; // NC1pi0onn - case 32: return 12; // NC1pi0onp - case 33: return 13; // NC1pimonn - case 34: return 14; // NC1piponp - case 36: return 15; // NCcoh - case 38: return 16; // NC1gammaonn - case 39: return 17; // NC1gammaonp - case 41: return 18; // NC Multipi - case 42: return 19; // NC1etaonn - case 43: return 20; // NC1etaonp - case 44: return 21; // NC1kamk0 - case 45: return 22; // NC1lamkp - case 46: return 23; // NCDIS - case 51: return 24; // NCEL on p - case 52: return 25; // NCEL on n - case 53: return 26; // NC 2p2h - default: return 27; // Undefined + case 1: + return 0; // CCQE + case 2: + return 1; // CC2p2h + case 11: + return 2; // CC1piponp + case 12: + return 3; // CC1pi0onn + case 13: + return 4; // CC1piponn + case 16: + return 5; // CCcoh + case 17: + return 6; // CC1gamma + case 21: + return 7; // CCMultipi + case 22: + return 8; // CC1eta + case 23: + return 9; // CC1lamkp + case 26: + return 10; // CCDIS + + case 31: + return 11; // NC1pi0onn + case 32: + return 12; // NC1pi0onp + case 33: + return 13; // NC1pimonn + case 34: + return 14; // NC1piponp + case 36: + return 15; // NCcoh + case 38: + return 16; // NC1gammaonn + case 39: + return 17; // NC1gammaonp + case 41: + return 18; // NC Multipi + case 42: + return 19; // NC1etaonn + case 43: + return 20; // NC1etaonp + case 44: + return 21; // NC1kamk0 + case 45: + return 22; // NC1lamkp + case 46: + return 23; // NCDIS + case 51: + return 24; // NCEL on p + case 52: + return 25; // NCEL on n + case 53: + return 26; // NC 2p2h + default: + return 27; // Undefined } }; -void TrueModeStack::Fill(int mode, double x, double y, double z, double weight) { +void TrueModeStack::Fill(int mode, double x, double y, double z, + double weight) { StackBase::FillStack(ConvertModeToIndex(mode), x, y, z, weight); }; -void TrueModeStack::Fill(FitEvent* evt, double x, double y, double z, double weight) { +void TrueModeStack::SetBinContent(int mode, int binx, int biny, int binz, double content) { + StackBase::SetBinContentStack(ConvertModeToIndex(mode), binx, biny, binz, + content); +} +void TrueModeStack::SetBinError(int mode, int binx, int biny, int binz, double error) { + StackBase::SetBinErrorStack(ConvertModeToIndex(mode), binx, biny, binz, + error); +} + +void TrueModeStack::Fill(FitEvent *evt, double x, double y, double z, + double weight) { StackBase::FillStack(ConvertModeToIndex(evt->Mode), x, y, z, weight); }; -void TrueModeStack::Fill(BaseFitEvt* evt, double x, double y, double z, double weight) { +void TrueModeStack::Fill(BaseFitEvt *evt, double x, double y, double z, + double weight) { StackBase::FillStack(ConvertModeToIndex(evt->Mode), x, y, z, weight); }; - - /// TrueModeStack Functions -NuNuBarTrueModeStack::NuNuBarTrueModeStack(std::string name, std::string title, TH1* hist) { +NuNuBarTrueModeStack::NuNuBarTrueModeStack(std::string name, std::string title, + TH1 *hist) { fName = name; fTitle = title; // Neutrino // CC - AddMode(0, "NU_CCQE", "#nu CCQE", kBlue, 2, 1001); - AddMode(1, "NU_CC2p2h", "#nu 2p2h", kRed, 2, 1001); - AddMode(2, "NU_CC1piponp", "#nu CC1#pi^{+} on p", kGreen, 2, 1001); - AddMode(3, "NU_CC1pi0onn", "#nu CC1#pi^{0} on n", kGreen + 3, 2, 1001); - AddMode(4, "NU_CC1piponn", "#nu CC1#pi^{+} on n", kGreen - 2, 2, 1001); - AddMode(5, "NU_CCcoh", "#nu CC coherent", kBlue, 2, 1001); - AddMode(6, "NU_CC1gamma", "#nu CC1#gamma", kMagenta, 2, 1001); - AddMode(7, "NU_CCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001); - AddMode(8, "NU_CC1eta", "#nu CC1#eta^{0} on n", kYellow - 2, 2, 1001); - AddMode(9, "NU_CC1lamkp", "#nu CC1#Lambda1K^{+}", kYellow - 6, 2, 1001); - AddMode(10, "NU_CCDIS", "#nu DIS (W > 2.0)", kRed, 2, 1001); + AddMode(0, "NU_CCQE", "#nu CCQE", kBlue, 2, 1001); + AddMode(1, "NU_CC2p2h", "#nu 2p2h", kRed, 2, 1001); + AddMode(2, "NU_CC1piponp", "#nu CC1#pi^{+} on p", kGreen, 2, 1001); + AddMode(3, "NU_CC1pi0onn", "#nu CC1#pi^{0} on n", kGreen + 3, 2, 1001); + AddMode(4, "NU_CC1piponn", "#nu CC1#pi^{+} on n", kGreen - 2, 2, 1001); + AddMode(5, "NU_CCcoh", "#nu CC coherent", kBlue, 2, 1001); + AddMode(6, "NU_CC1gamma", "#nu CC1#gamma", kMagenta, 2, 1001); + AddMode(7, "NU_CCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001); + AddMode(8, "NU_CC1eta", "#nu CC1#eta^{0} on n", kYellow - 2, 2, 1001); + AddMode(9, "NU_CC1lamkp", "#nu CC1#Lambda1K^{+}", kYellow - 6, 2, 1001); + AddMode(10, "NU_CCDIS", "#nu DIS (W > 2.0)", kRed, 2, 1001); // NC - AddMode(11, "NU_NC1pi0onn", "#nu NC1#pi^{0} on n", kBlue, 2, 3004); - AddMode(12, "NU_NC1pi0onp", "#nu NC1#pi^{0} on p", kBlue + 3, 2, 3004); - AddMode(13, "NU_NC1pimonn", "#nu NC1#pi^{-} on n", kBlue - 2, 2, 3004); - AddMode(14, "NU_NC1piponp", "#nu NC1#pi^{+} on p", kBlue - 8, 2, 3004); - AddMode(15, "NU_NCcoh", "#nu NC Coherent", kBlue + 8, 2, 3004); - AddMode(16, "NU_NC1gammaonn", "#nu NC1#gamma on n", kMagenta, 2, 3004); - AddMode(17, "NU_NC1gammaonp", "#nu NC1#gamma on p", kMagenta - 10, 2, 3004); - AddMode(18, "NU_NCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004); - AddMode(19, "NU_NC1etaonn", "#nu NC1#eta^{0} on n", kYellow - 2, 2, 3004); - AddMode(20, "NU_NC1etaonp", "#nu NC1#eta^{0} on p", kYellow - 4, 2, 3004); - AddMode(21, "NU_NC1kamk0", "#nu NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004); - AddMode(22, "NU_NC1lamkp", "#nu NC1#Lambda1K^{+}", kYellow - 10, 2, 3004); - AddMode(23, "NU_NCDIS", "#nu DIS (W > 2.0)", kRed, 2, 3004); - AddMode(24, "NU_NCELonp", "#nu NCEL on p", kBlack, 2, 3004); - AddMode(25, "NU_NCELonn", "#nu NCEL on n", kGray, 2, 3004); + AddMode(11, "NU_NC1pi0onn", "#nu NC1#pi^{0} on n", kBlue, 2, 3004); + AddMode(12, "NU_NC1pi0onp", "#nu NC1#pi^{0} on p", kBlue + 3, 2, 3004); + AddMode(13, "NU_NC1pimonn", "#nu NC1#pi^{-} on n", kBlue - 2, 2, 3004); + AddMode(14, "NU_NC1piponp", "#nu NC1#pi^{+} on p", kBlue - 8, 2, 3004); + AddMode(15, "NU_NCcoh", "#nu NC Coherent", kBlue + 8, 2, 3004); + AddMode(16, "NU_NC1gammaonn", "#nu NC1#gamma on n", kMagenta, 2, 3004); + AddMode(17, "NU_NC1gammaonp", "#nu NC1#gamma on p", kMagenta - 10, 2, 3004); + AddMode(18, "NU_NCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, + 3004); + AddMode(19, "NU_NC1etaonn", "#nu NC1#eta^{0} on n", kYellow - 2, 2, 3004); + AddMode(20, "NU_NC1etaonp", "#nu NC1#eta^{0} on p", kYellow - 4, 2, 3004); + AddMode(21, "NU_NC1kamk0", "#nu NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004); + AddMode(22, "NU_NC1lamkp", "#nu NC1#Lambda1K^{+}", kYellow - 10, 2, 3004); + AddMode(23, "NU_NCDIS", "#nu DIS (W > 2.0)", kRed, 2, 3004); + AddMode(24, "NU_NCELonp", "#nu NCEL on p", kBlack, 2, 3004); + AddMode(25, "NU_NCELonn", "#nu NCEL on n", kGray, 2, 3004); // Undefined - AddMode(26, "NU_UNDEFINED", "#nu Undefined", kRed + 2, 2, 3000); + AddMode(26, "NU_UNDEFINED", "#nu Undefined", kRed + 2, 2, 3000); // CC - AddMode(27, "ANTINU_CCQE", "#bar{#nu} CCQE", kBlue, 2, 1001); - AddMode(28, "ANTINU_CC2p2h", "#bar{#nu} 2p2h", kRed, 2, 1001); - AddMode(29, "ANTINU_CC1piponp", "#bar{#nu} CC1#pi^{+} on p", kGreen, 2, 1001); - AddMode(30, "ANTINU_CC1pi0onn", "#bar{#nu} CC1#pi^{0} on n", kGreen + 3, 2, 1001); - AddMode(31, "ANTINU_CC1piponn", "#bar{#nu} CC1#pi^{+} on n", kGreen - 2, 2, 1001); - AddMode(32, "ANTINU_CCcoh", "#bar{#nu} CC coherent", kBlue, 2, 1001); - AddMode(33, "ANTINU_CC1gamma", "#bar{#nu} CC1#gamma", kMagenta, 2, 1001); - AddMode(34, "ANTINU_CCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001); - AddMode(35, "ANTINU_CC1eta", "#bar{#nu} CC1#eta^{0} on n", kYellow - 2, 2, 1001); - AddMode(36, "ANTINU_CC1lamkp", "#bar{#nu} CC1#Lambda1K^{+}", kYellow - 6, 2, 1001); - AddMode(37, "ANTINU_CCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 1001); + AddMode(27, "ANTINU_CCQE", "#bar{#nu} CCQE", kBlue, 2, 1001); + AddMode(28, "ANTINU_CC2p2h", "#bar{#nu} 2p2h", kRed, 2, 1001); + AddMode(29, "ANTINU_CC1piponp", "#bar{#nu} CC1#pi^{+} on p", kGreen, 2, 1001); + AddMode(30, "ANTINU_CC1pi0onn", "#bar{#nu} CC1#pi^{0} on n", kGreen + 3, 2, + 1001); + AddMode(31, "ANTINU_CC1piponn", "#bar{#nu} CC1#pi^{+} on n", kGreen - 2, 2, + 1001); + AddMode(32, "ANTINU_CCcoh", "#bar{#nu} CC coherent", kBlue, 2, 1001); + AddMode(33, "ANTINU_CC1gamma", "#bar{#nu} CC1#gamma", kMagenta, 2, 1001); + AddMode(34, "ANTINU_CCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", + kYellow, 2, 1001); + AddMode(35, "ANTINU_CC1eta", "#bar{#nu} CC1#eta^{0} on n", kYellow - 2, 2, + 1001); + AddMode(36, "ANTINU_CC1lamkp", "#bar{#nu} CC1#Lambda1K^{+}", kYellow - 6, 2, + 1001); + AddMode(37, "ANTINU_CCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 1001); // NC - AddMode(38, "ANTINU_NC1pi0onn", "#bar{#nu} NC1#pi^{0} on n", kBlue, 2, 3004); - AddMode(39, "ANTINU_NC1pi0onp", "#bar{#nu} NC1#pi^{0} on p", kBlue + 3, 2, 3004); - AddMode(40, "ANTINU_NC1pimonn", "#bar{#nu} NC1#pi^{-} on n", kBlue - 2, 2, 3004); - AddMode(41, "ANTINU_NC1piponp", "#bar{#nu} NC1#pi^{+} on p", kBlue - 8, 2, 3004); - AddMode(42, "ANTINU_NCcoh", "#bar{#nu} NC Coherent", kBlue + 8, 2, 3004); - AddMode(43, "ANTINU_NC1gammaonn", "#bar{#nu} NC1#gamma on n", kMagenta, 2, 3004); - AddMode(44, "ANTINU_NC1gammaonp", "#bar{#nu} NC1#gamma on p", kMagenta - 10, 2, 3004); - AddMode(45, "ANTINU_NCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004); - AddMode(46, "ANTINU_NC1etaonn", "#bar{#nu} NC1#eta^{0} on n", kYellow - 2, 2, 3004); - AddMode(47, "ANTINU_NC1etaonp", "#bar{#nu} NC1#eta^{0} on p", kYellow - 4, 2, 3004); - AddMode(48, "ANTINU_NC1kamk0", "#bar{#nu} NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004); - AddMode(49, "ANTINU_NC1lamkp", "#bar{#nu} NC1#Lambda1K^{+}", kYellow - 10, 2, 3004); - AddMode(50, "ANTINU_NCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 3004); - AddMode(51, "ANTINU_NCELonp", "#bar{#nu} NCEL on p", kBlack, 2, 3004); - AddMode(52, "ANTINU_NCELonn", "#bar{#nu} NCEL on n", kGray, 2, 3004); + AddMode(38, "ANTINU_NC1pi0onn", "#bar{#nu} NC1#pi^{0} on n", kBlue, 2, 3004); + AddMode(39, "ANTINU_NC1pi0onp", "#bar{#nu} NC1#pi^{0} on p", kBlue + 3, 2, + 3004); + AddMode(40, "ANTINU_NC1pimonn", "#bar{#nu} NC1#pi^{-} on n", kBlue - 2, 2, + 3004); + AddMode(41, "ANTINU_NC1piponp", "#bar{#nu} NC1#pi^{+} on p", kBlue - 8, 2, + 3004); + AddMode(42, "ANTINU_NCcoh", "#bar{#nu} NC Coherent", kBlue + 8, 2, 3004); + AddMode(43, "ANTINU_NC1gammaonn", "#bar{#nu} NC1#gamma on n", kMagenta, 2, + 3004); + AddMode(44, "ANTINU_NC1gammaonp", "#bar{#nu} NC1#gamma on p", kMagenta - 10, + 2, 3004); + AddMode(45, "ANTINU_NCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", + kBlue - 10, 2, 3004); + AddMode(46, "ANTINU_NC1etaonn", "#bar{#nu} NC1#eta^{0} on n", kYellow - 2, 2, + 3004); + AddMode(47, "ANTINU_NC1etaonp", "#bar{#nu} NC1#eta^{0} on p", kYellow - 4, 2, + 3004); + AddMode(48, "ANTINU_NC1kamk0", "#bar{#nu} NC1#Lambda1K^{0} on n", kYellow - 6, + 2, 3004); + AddMode(49, "ANTINU_NC1lamkp", "#bar{#nu} NC1#Lambda1K^{+}", kYellow - 10, 2, + 3004); + AddMode(50, "ANTINU_NCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 3004); + AddMode(51, "ANTINU_NCELonp", "#bar{#nu} NCEL on p", kBlack, 2, 3004); + AddMode(52, "ANTINU_NCELonn", "#bar{#nu} NCEL on n", kGray, 2, 3004); // Undefined - AddMode(53, "NU_UNDEFINED", "#bar{#nu} Undefined", kRed + 2, 2, 3000); + AddMode(53, "NU_UNDEFINED", "#bar{#nu} Undefined", kRed + 2, 2, 3000); // Non Neutrino AddMode(54, "UNDEFINED", "Non-#nu Undefined", kBlack, 2, 3000); StackBase::SetupStack(hist); }; int NuNuBarTrueModeStack::ConvertModeToIndex(int mode) { switch (abs(mode)) { - case 1: return 0; // CCQE - case 2: return 1; // CC2p2h - case 11: return 2; // CC1piponp - case 12: return 3; // CC1pi0onn - case 13: return 4; // CC1piponn - case 16: return 5; // CCcoh - case 17: return 6; // CC1gamma - case 21: return 7; // CCMultipi - case 22: return 8; // CC1eta - case 23: return 9; // CC1lamkp - case 26: return 10; // CCDIS - - case 31: return 11; // NC1pi0onn - case 32: return 12; // NC1pi0onp - case 33: return 13; // NC1pimonn - case 34: return 14; // NC1piponp - case 36: return 15; // NCcoh - case 38: return 16; // NC1gammaonn - case 39: return 17; // NC1gammaonp - case 41: return 18; // NC Multipi - case 42: return 19; // NC1etaonn - case 43: return 20; // NC1etaonp - case 44: return 21; // NC1kamk0 - case 45: return 22; // NC1lamkp - case 46: return 23; // NCDIS - case 51: return 24; // NCEL on p - case 52: return 25; // NCEL on n - default: return 26; // Undefined + case 1: + return 0; // CCQE + case 2: + return 1; // CC2p2h + case 11: + return 2; // CC1piponp + case 12: + return 3; // CC1pi0onn + case 13: + return 4; // CC1piponn + case 16: + return 5; // CCcoh + case 17: + return 6; // CC1gamma + case 21: + return 7; // CCMultipi + case 22: + return 8; // CC1eta + case 23: + return 9; // CC1lamkp + case 26: + return 10; // CCDIS + + case 31: + return 11; // NC1pi0onn + case 32: + return 12; // NC1pi0onp + case 33: + return 13; // NC1pimonn + case 34: + return 14; // NC1piponp + case 36: + return 15; // NCcoh + case 38: + return 16; // NC1gammaonn + case 39: + return 17; // NC1gammaonp + case 41: + return 18; // NC Multipi + case 42: + return 19; // NC1etaonn + case 43: + return 20; // NC1etaonp + case 44: + return 21; // NC1kamk0 + case 45: + return 22; // NC1lamkp + case 46: + return 23; // NCDIS + case 51: + return 24; // NCEL on p + case 52: + return 25; // NCEL on n + default: + return 26; // Undefined } }; -void NuNuBarTrueModeStack::Fill(int species, int mode, double x, double y, double z, double weight) { +void NuNuBarTrueModeStack::Fill(int species, int mode, double x, double y, + double z, double weight) { int modeindex = ConvertModeToIndex(mode); int index = 54; // undefined - if (species == 12 or species == 14 or species == 16) index = modeindex; - else if (species == -12 or species == -14 or species == -16) index = modeindex + 27; + if (species == 12 or species == 14 or species == 16) + index = modeindex; + else if (species == -12 or species == -14 or species == -16) + index = modeindex + 27; StackBase::FillStack(index, x, y, z, weight); }; - - - // Species Stack Functions -BeamSpeciesStack::BeamSpeciesStack(std::string name, std::string title, TH1* hist) { +BeamSpeciesStack::BeamSpeciesStack(std::string name, std::string title, + TH1 *hist) { fName = name; fTitle = title; // charged eptons AddMode(0, "electron", "e^{-}", kBlue, 2, 1001); AddMode(1, "positron", "e^{+}", kBlue - 2, 2, 3004); AddMode(2, "muon", "#mu^{-}", kRed, 2, 1001); AddMode(3, "antimuon", "#mu^{+}", kRed - 2, 2, 3004); AddMode(4, "tau", "#tau^{-}", kGreen, 2, 1001); AddMode(5, "antitau", "#tau^{+}", kGreen - 2, 2, 3004); // neutrinos AddMode(6, "nue", "#nu_e", kBlue, 2, 1001); AddMode(7, "antinue", "#bar{#nu}_e", kBlue - 2, 2, 3004); AddMode(8, "numu", "#nu_#mu", kRed, 2, 1001); AddMode(9, "antinumu", "#bar{#nu}_#mu", kRed - 2, 2, 3004); AddMode(10, "nutau", "#nu_#tau", kGreen, 2, 1001); AddMode(11, "antinutau", "#bar{#nu}_#tau", kGreen - 2, 2, 3004); StackBase::SetupStack(hist); }; int BeamSpeciesStack::ConvertSpeciesToIndex(int species) { switch (species) { - case 11: return 0; // e- - case -11: return 1; // e+ - case 13: return 2; // mu- - case -13: return 3; // mu+ - case 15: return 4; // tau- - case -15: return 5; // tau+ - case 12: return 6; // nue - case -12: return 7; // nuebar - case 14: return 8; // numu - case -14: return 9; // numubar - case 16: return 10; // nutau - case -16: return 11; //nutaubar - default: return 12; + case 11: + return 0; // e- + case -11: + return 1; // e+ + case 13: + return 2; // mu- + case -13: + return 3; // mu+ + case 15: + return 4; // tau- + case -15: + return 5; // tau+ + case 12: + return 6; // nue + case -12: + return 7; // nuebar + case 14: + return 8; // numu + case -14: + return 9; // numubar + case 16: + return 10; // nutau + case -16: + return 11; // nutaubar + default: + return 12; } }; -void BeamSpeciesStack::Fill(int species, double x, double y, double z, double weight) { +void BeamSpeciesStack::Fill(int species, double x, double y, double z, + double weight) { StackBase::FillStack(ConvertSpeciesToIndex(species), x, y, z, weight); } - - // Target Stack Functions -TargetTypeStack::TargetTypeStack(std::string name, std::string title, TH1* hist){ +TargetTypeStack::TargetTypeStack(std::string name, std::string title, + TH1 *hist) { fName = name; fTitle = title; AddMode(0, "H", "Hydrogen", kBlue, 2, 1001); AddMode(1, "C", "Carbon", kRed, 2, 1001); AddMode(2, "O", "Oxygen", kViolet, 2, 1001); - AddMode(3, "UNDEFINED", "Undefined", kBlack, 2, 1001 ); + AddMode(3, "UNDEFINED", "Undefined", kBlack, 2, 1001); StackBase::SetupStack(hist); } -void TargetTypeStack::Fill(int pdg, double x, double y, double z, double weight){ +void TargetTypeStack::Fill(int pdg, double x, double y, double z, + double weight) { int index = ConvertPDGToIndex(pdg); StackBase::FillStack(index, x, y, z, weight); } -int TargetTypeStack::ConvertPDGToIndex(int pdg){ - switch(pdg){ - case 1000010010: return 0; // H - case 1000060120: return 1; // C - case 1000080160: return 2; // O - default: return 3; // Undef +int TargetTypeStack::ConvertPDGToIndex(int pdg) { + switch (pdg) { + case 1000010010: + return 0; // H + case 1000060120: + return 1; // C + case 1000080160: + return 2; // O + default: + return 3; // Undef } } -bool TargetTypeStack::IncludeInStack(TH1* hist){ +bool TargetTypeStack::IncludeInStack(TH1 *hist) { return (hist->Integral() > 0.0); } - - - // CC Topology Stack Functions -CCTopologyStack::CCTopologyStack(std::string name, std::string title, TH1* hist) { +CCTopologyStack::CCTopologyStack(std::string name, std::string title, + TH1 *hist) { fName = name; fTitle = title; AddMode(0, "CC0pi", "CC-0#pi", kBlue, 2, 1001); AddMode(1, "CC1pip", "CC-1#pi^{+}", kRed, 2, 1001); AddMode(2, "CC1pim", "CC-1#pi^{-}", kGreen, 2, 1001); AddMode(3, "CC1pi0", "CC-1#pi^{0}", kYellow, 2, 1001); AddMode(4, "CCNpi", "CC-N#pi", kGray, 2, 1001); AddMode(5, "CCOther", "CC-Other", kViolet, 2, 1001); - AddMode(6, "NC", "NC", kMagenta, 2, 1001); + AddMode(6, "NC", "NC", kMagenta, 2, 1001); AddMode(7, "UNDEFINED", "Undefined", kBlack, 2, 1001); - - } -void CCTopologyStack::Fill(FitEvent* evt, double x, double y, double z, double weight) { +void CCTopologyStack::Fill(FitEvent *evt, double x, double y, double z, + double weight) { int index = GetIndexFromEventParticles(evt); StackBase::FillStack(index, x, y, z, weight); } -int CCTopologyStack::GetIndexFromEventParticles(FitEvent* evt) { +int CCTopologyStack::GetIndexFromEventParticles(FitEvent *evt) { int nleptons = evt->NumFSLeptons(); - int npiplus = evt->NumFSParticle(211); - int npineg = evt->NumFSParticle(-211); - int npi0 = evt->NumFSParticle(111); + int npiplus = evt->NumFSParticle(211); + int npineg = evt->NumFSParticle(-211); + int npi0 = evt->NumFSParticle(111); int npions = npiplus + npineg + npi0; if (nleptons == 1) { if (npions == 0) { return 0; // CC0pi } else if (npions == 1) { - if (npiplus == 1) return 1; //CC1pi+ - else if (npineg == 1) return 2; //CC1pi- - else if (npi0 == 1) return 3; //CC1pi0 + if (npiplus == 1) + return 1; // CC1pi+ + else if (npineg == 1) + return 2; // CC1pi- + else if (npi0 == 1) + return 3; // CC1pi0 } else if (npions > 1) { return 4; // CCNpi } } else if (nleptons > 1) { return 5; // CCOther } else if (nleptons < 1) { return 6; } return 7; // Undefined? - } - - - - - - - - - - - diff --git a/src/FitBase/StandardStacks.h b/src/FitBase/StandardStacks.h index 4603608..322832b 100644 --- a/src/FitBase/StandardStacks.h +++ b/src/FitBase/StandardStacks.h @@ -1,133 +1,129 @@ #ifndef STANDARD_STACKS_H #define STANDARD_STACKS_H -#include "StackBase.h" -#include "FitEvent.h" #include "BaseFitEvt.h" +#include "FitEvent.h" +#include "StackBase.h" -/// Single stack class, for consistent handling of TH1* and StackBase* histograms. +/// Single stack class, for consistent handling of TH1* and StackBase* +/// histograms. class FakeStack : public StackBase { public: - /// Sets Template to exact pointer to hist without cloning. - FakeStack(TH1* hist) ; + /// Sets Template to exact pointer to hist without cloning. + FakeStack(TH1 *hist); - /// Sets NULL template and saves TF1 object instead - FakeStack(TF1* f); + /// Sets NULL template and saves TF1 object instead + FakeStack(TF1 *f); - /// Sets NULL template and saves TGraph object instead - FakeStack(TGraph* gr); + /// Sets NULL template and saves TGraph object instead + FakeStack(TGraph *gr); - /// Unlinks pointer to original histogram - ~FakeStack(); + /// Unlinks pointer to original histogram + ~FakeStack(); - /// Fills the normal fTemplate histogram - void Fill(double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// Fills the normal fTemplate histogram + void Fill(double x, double y = 1.0, double z = 1.0, double weight = 1.0); - /// Scales the normal fTemplate histogram - void Scale(double norm, std::string opt); + /// Scales the normal fTemplate histogram + void Scale(double norm, std::string opt); - /// Resets the normal fTemplate histogram - void Reset(); + /// Resets the normal fTemplate histogram + void Reset(); - /// Writes the normal fTemplate histogram - void Write(); + /// Writes the normal fTemplate histogram + void Write(); - std::string fFakeType; - TGraph* fTGraphObject; - TF1* fTF1Object; + std::string fFakeType; + TGraph *fTGraphObject; + TF1 *fTF1Object; }; - - - - /// True Mode stack, an array of neut true interaction channels class TrueModeStack : public StackBase { public: + /// Main constructor listing true mode categories. + TrueModeStack(std::string name, std::string title, TH1 *hist); - /// Main constructor listing true mode categories. - TrueModeStack(std::string name, std::string title, TH1* hist); - - /// List to convert Modes to Index. - /// Should be kept in sync with constructor. - int ConvertModeToIndex(int mode); + /// List to convert Modes to Index. + /// Should be kept in sync with constructor. + static int ConvertModeToIndex(int mode); - /// Fill fromgiven mode integer - void Fill(int mode, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// Fill fromgiven mode integer + void Fill(int mode, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); - /// Extracts Mode from FitEvent and fills - void Fill(FitEvent* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + void SetBinContent(int mode, int binx, int biny, int binz, double content); + void SetBinError(int mode, int binx, int biny, int binz, double error); - /// Extracts Mode from BaseFitEvt - void Fill(BaseFitEvt* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// Extracts Mode from FitEvent and fills + void Fill(FitEvent *evt, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); + /// Extracts Mode from BaseFitEvt + void Fill(BaseFitEvt *evt, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); }; /// True Mode NuNuBar stack, array of true channels split by nu/nubar class NuNuBarTrueModeStack : public StackBase { public: + /// Main constructor listing true mode categories. + NuNuBarTrueModeStack(std::string name, std::string title, TH1 *hist); - /// Main constructor listing true mode categories. - NuNuBarTrueModeStack(std::string name, std::string title, TH1* hist); + /// List to convert Modes to Index. + /// Should be kept in sync with constructor. + int ConvertModeToIndex(int mode); - /// List to convert Modes to Index. - /// Should be kept in sync with constructor. - int ConvertModeToIndex(int mode); - - /// Fill fromgiven mode integer - void Fill(int species, int mode, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// Fill fromgiven mode integer + void Fill(int species, int mode, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); }; /// Species stack to look at contributions from multiple beam leptons class BeamSpeciesStack : public StackBase { public: + /// main constructor listing beam categories + BeamSpeciesStack(std::string name, std::string title, TH1 *hist); - /// main constructor listing beam categories - BeamSpeciesStack(std::string name, std::string title, TH1* hist); - - /// Fills stack using neutrino species - void Fill(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// Fills stack using neutrino species + void Fill(int species, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); - /// List converts PDG Beam to index. - /// Should be kept in sync with constructor. - int ConvertSpeciesToIndex(int species); + /// List converts PDG Beam to index. + /// Should be kept in sync with constructor. + int ConvertSpeciesToIndex(int species); }; - /// Species stack to look at contributions from multiple beam leptons class TargetTypeStack : public StackBase { public: + /// main constructor listing beam categories + TargetTypeStack(std::string name, std::string title, TH1 *hist); - /// main constructor listing beam categories - TargetTypeStack(std::string name, std::string title, TH1* hist); + /// Fills stack using target pdg + void Fill(int pdg, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); - /// Fills stack using target pdg - void Fill(int pdg, double x, double y = 1.0, double z = 1.0, double weight = 1.0); + /// List converts PDG Beam to index. + /// Should be kept in sync with constructor. + int ConvertPDGToIndex(int pdg); - /// List converts PDG Beam to index. - /// Should be kept in sync with constructor. - int ConvertPDGToIndex(int pdg); - - /// Specie empty histograms to not be included in final stack - bool IncludeInStack(TH1* hist); + /// Specie empty histograms to not be included in final stack + bool IncludeInStack(TH1 *hist); }; - - /// CC Topology Stack, categories defined by final state particle counts. class CCTopologyStack : public StackBase { public: + /// main constructor listing beam categories + CCTopologyStack(std::string name, std::string title, TH1 *hist); - /// main constructor listing beam categories - CCTopologyStack(std::string name, std::string title, TH1* hist); - - /// Fills stack using FitEvent - void Fill(FitEvent* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0); - - /// Extracts index from evt particle counts - int GetIndexFromEventParticles(FitEvent* evt); - + /// Fills stack using FitEvent + void Fill(FitEvent *evt, double x, double y = 1.0, double z = 1.0, + double weight = 1.0); + /// Extracts index from evt particle counts + int GetIndexFromEventParticles(FitEvent *evt); }; #endif diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx index 0875981..be6e097 100644 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx @@ -1,330 +1,373 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "T2K_SignalDef.h" #include "T2K_CC0pi_XSec_2DPcos_nu_I.h" static size_t nangbins = 9; static double angular_binning_costheta[] = {-1, 0, 0.6, 0.7, 0.8, 0.85, 0.9, 0.94, 0.98, 1}; //******************************************************************** T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n" "Target: CH \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" "Signal: CC0pi\n" "https://journals.aps.org/prd/abstract/10.1103/" "PhysRevD.93.112012 Analysis I"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("P_{#mu} (GeV)"); fSettings.SetYTitle("cos#theta_{#mu}"); fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV", "FIX/FULL"); fSettings.DefineAllowedTargets("C,H"); fSettings.SetEnuRangeFromFlux(fFluxHist); // CCQELike plot information fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); fMaskMomOverflow = false; if (samplekey.Has("mask_mom_overflow")) { fMaskMomOverflow = samplekey.GetB("mask_mom_overflow"); } // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * 1E-38 / (TotalIntegratedFlux())); assert(std::isnormal(fScaleFactor)); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) { return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I); }; void T2K_CC0pi_XSec_2DPcos_nu_I::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double pmu = Pmu.Vect().Mag() / 1000.; double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); fXVar = pmu; fYVar = CosThetaMu; + fMode = event->Mode; + return; }; void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() { Measurement1D::FillHistograms(); if (Signal) { fMCHist_Fine2D->Fill(fXVar, fYVar, Weight); FillMCSlice(fXVar, fYVar, Weight); } } // Modification is needed after the full reconfigure to move bins around // Otherwise this would need to be replaced by a TH2Poly which is too awkward. void T2K_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() { // Do standard conversion. Measurement1D::ConvertEventRates(); // Scale MC slices by their bin area for (size_t i = 0; i < nangbins; ++i) { fMCHist_Slices[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]), "width"); } // Now Convert into 1D list fMCHist->Reset(); int bincount = 0; for (size_t i = 0; i < nangbins; i++) { for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fMCHist->SetBinContent(bincount + 1, fMCHist_Slices[i]->GetBinContent(j + 1)); fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1)); + bincount++; } } + if (fMCHist_Modes) { + fMCHist_Modes->Reset(); + + for (std::map >::iterator it = + fMCModeHists_Slices.begin(); + it != fMCModeHists_Slices.end(); ++it) { + + bincount = 0; + for (size_t i = 0; i < nangbins; i++) { + it->second[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] - + angular_binning_costheta[i]), + "width"); + + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { + fMCHist_Modes->SetBinContent(it->first, bincount + 1, 0, 0, + it->second[i]->GetBinContent(j + 1)); + fMCHist_Modes->SetBinError(it->first, bincount + 1, 0, 0, + it->second[i]->GetBinError(j + 1)); + + bincount++; + + } + } + } + } + return; } void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) { + if (fMCHist_Modes && !fMCModeHists_Slices.count(fMode)) { + for (size_t i = 0; i < nangbins; ++i) { + std::stringstream ss(""); + ss << "T2K_CC0pi_XSec_2DPcos_nu_I_MODE_ " << fMode << "_slice" << i + << std::endl; + fMCModeHists_Slices[fMode].push_back( + static_cast(fMCHist_Slices[i]->Clone(ss.str().c_str()))); + fMCModeHists_Slices[fMode].back()->Reset(); + } + } + for (size_t i = 0; i < nangbins; ++i) { if ((y >= angular_binning_costheta[i]) && (y < angular_binning_costheta[i + 1])) { fMCHist_Slices[i]->Fill(x, w); + if (fMCHist_Modes) { + fMCModeHists_Slices[fMode][i]->Fill(x, w); + } } } } void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() { // Read in 1D Data Histograms TFile input( (FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root") .c_str(), "READ"); fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", "T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0, 100, -1.0, 1.0); fMCHist_Fine2D->SetDirectory(NULL); SetAutoProcessTH1(fMCHist_Fine2D); TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov"); fFullCovar = new TMatrixDSym(tempcov->GetNbinsX()); for (int i = 0; i < tempcov->GetNbinsX(); i++) { for (int j = 0; j < tempcov->GetNbinsX(); j++) { (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); } } // Read in 2D Data Slices and Make MC Slices int bincount = 0; for (size_t i = 0; i < nangbins; i++) { fDataHist_Slices.push_back( (TH1D *)input.Get(Form("dataslice_%lu", i))->Clone()); fDataHist_Slices[i]->SetNameTitle( Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu", i), (Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu, cos(#theta) [%f,%f] ", i, angular_binning_costheta[i], angular_binning_costheta[i + 1]))); fDataHist_Slices.back()->SetDirectory(NULL); // Loop over nbins and set errors from covar for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fDataHist_Slices[i]->SetBinError( j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38); bincount++; } } assert(bincount == tempcov->GetNbinsX()); if (fMaskMomOverflow) { MaskMomOverflow(); bincount = fFullCovar->GetNcols(); } std::vector > data_slice_bcbes; for (size_t i = 0; i < nangbins; i++) { for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { data_slice_bcbes.push_back( std::make_pair(fDataHist_Slices[i]->GetBinContent(j + 1), fDataHist_Slices[i]->GetBinError(j + 1))); } } for (size_t i = 0; i < nangbins; i++) { fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone()); fMCHist_Slices.back()->SetDirectory(NULL); fMCHist_Slices.back()->Reset(); fMCHist_Slices.back()->SetLineColor(kRed); fMCHist_Slices[i]->SetNameTitle( Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu", i), (Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu, cos(#theta) [%f,%f] ", i, angular_binning_costheta[i], angular_binning_costheta[i + 1]))); } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); fDataHist = new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", "T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount); fDataHist->SetDirectory(NULL); for (size_t i = 0; i < data_slice_bcbes.size(); ++i) { fDataHist->SetBinContent(i + 1, data_slice_bcbes[i].first); fDataHist->SetBinError(i + 1, data_slice_bcbes[i].second); } SetShapeCovar(); fMCHist = (TH1D *)fDataHist->Clone("T2K_CC0pi_XSec_2DPcos_nu_I_MC_1D"); fMCHist->SetDirectory(NULL); return; } void T2K_CC0pi_XSec_2DPcos_nu_I::MaskMomOverflow() { std::vector bins_to_cut; size_t nallbins = 0; for (size_t i = 0; i < nangbins; i++) { std::vector slice_bin_edges; slice_bin_edges.push_back( fDataHist_Slices[i]->GetXaxis()->GetBinLowEdge(1)); for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) { slice_bin_edges.push_back( fDataHist_Slices[i]->GetXaxis()->GetBinUpEdge(j + 1)); nallbins++; } bins_to_cut.push_back(nallbins++); TH1D *tmp = new TH1D(fDataHist_Slices[i]->GetName(), fDataHist_Slices[i]->GetTitle(), slice_bin_edges.size() - 1, slice_bin_edges.data()); tmp->SetDirectory(NULL); for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) { tmp->SetBinContent(j + 1, fDataHist_Slices[i]->GetBinContent(j + 1)); tmp->SetBinError(j + 1, fDataHist_Slices[i]->GetBinError(j + 1)); } delete fDataHist_Slices[i]; fDataHist_Slices[i] = tmp; } TMatrixDSym *tmpcovar = new TMatrixDSym(nallbins - bins_to_cut.size()); int icut = 0; for (int ifull = 0; ifull < fFullCovar->GetNcols(); ifull++) { if (std::find(bins_to_cut.begin(), bins_to_cut.end(), ifull) != bins_to_cut.end()) { continue; } int jcut = 0; for (int jfull = 0; jfull < fFullCovar->GetNcols(); jfull++) { if (std::find(bins_to_cut.begin(), bins_to_cut.end(), jfull) != bins_to_cut.end()) { continue; } (*tmpcovar)[icut][jcut] = (*fFullCovar)[ifull][jfull]; jcut++; } icut++; } delete fFullCovar; fFullCovar = tmpcovar; } void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) { this->Measurement1D::Write(drawopt); for (size_t i = 0; i < nangbins; i++) { fMCHist_Slices[i]->Write(); fDataHist_Slices[i]->Write(); } if (fResidualHist) { std::vector MCResidual_Slices; size_t tb_it = 0; for (size_t i = 0; i < fMCHist_Slices.size(); ++i) { std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%lu", i); MCResidual_Slices.push_back( static_cast(fMCHist_Slices[i]->Clone(name.c_str()))); MCResidual_Slices.back()->Reset(); for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { double bc = fResidualHist->GetBinContent(tb_it + 1); MCResidual_Slices.back()->SetBinContent(j + 1, bc); tb_it++; } MCResidual_Slices.back()->Write(); } } if (fChi2LessBinHist) { std::vector MCChi2LessBin_Slices; size_t tb_it = 0; for (size_t i = 0; i < fMCHist_Slices.size(); ++i) { std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%lu", i); MCChi2LessBin_Slices.push_back( static_cast(fMCHist_Slices[i]->Clone(name.c_str()))); MCChi2LessBin_Slices.back()->Reset(); for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { double bc = fChi2LessBinHist->GetBinContent(tb_it + 1); MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc); tb_it++; } MCChi2LessBin_Slices.back()->Write(); } } } diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h index 2570f1e..63f7a26 100644 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h @@ -1,72 +1,74 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* -* This file is part of NUISANCE. -* -* NUISANCE is free software: you can redistribute it and/or modify -* it under the terms of the GNU General Public License as published by -* the Free Software Foundation, either version 3 of the License, or -* (at your option) any later version. -* -* NUISANCE is distributed in the hope that it will be useful, -* but WITHOUT ANY WARRANTY; without even the implied warranty of -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -* GNU General Public License for more details. -* -* You should have received a copy of the GNU General Public License -* along with NUISANCE. If not, see . -*******************************************************************************/ + * This file is part of NUISANCE. + * + * NUISANCE is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * NUISANCE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NUISANCE. If not, see . + *******************************************************************************/ #ifndef T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN #define T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN #include "Measurement1D.h" -#include "TH2Poly.h" #include "MeasurementVariableBox2D.h" +#include "TH2Poly.h" class T2K_CC0pi_XSec_2DPcos_nu_I : public Measurement1D { public: - T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey); /// Virtual Destructor - ~T2K_CC0pi_XSec_2DPcos_nu_I() {}; + ~T2K_CC0pi_XSec_2DPcos_nu_I(){}; /// Numu CC0PI Signal Definition /// /// /item bool isSignal(FitEvent *nvect); /// Read histograms in a special way because format is different. - /// Read from FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root" + /// Read from + /// FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root" void SetHistograms(); /// Bin Tmu CosThetaMu - void FillEventVariables(FitEvent* customEvent); + void FillEventVariables(FitEvent *customEvent); // Fill Histograms void FillHistograms(); /// Have to do a weird event scaling for analysis 1 void ConvertEventRates(); void Write(std::string drawopt); /// \brief Create Q2 Box to save correction info - inline MeasurementVariableBox* CreateBox(){ return new MeasurementVariableBox2D(); }; - - private: + inline MeasurementVariableBox *CreateBox() { + return new MeasurementVariableBox2D(); + }; +private: bool fIsSystCov, fIsStatCov, fIsNormCov; bool fMaskMomOverflow; + int fMode; - TH2D* fMCHist_Fine2D; + TH2D *fMCHist_Fine2D; - std::vector fMCHist_Slices; - std::vector fDataHist_Slices; + std::vector fMCHist_Slices; + std::map > fMCModeHists_Slices; + std::vector fDataHist_Slices; void FillMCSlice(double x, double y, double w); void MaskMomOverflow(); - }; #endif