diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx index 01b866d..405312a 100644 --- a/src/MINERvA/MINERvA_SignalDef.cxx +++ b/src/MINERvA/MINERvA_SignalDef.cxx @@ -1,512 +1,521 @@ // Copyright 2016-2021 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 "FitUtils.h" #include "SignalDef.h" #include "MINERvA_SignalDef.h" #include "MINERvAUtils.h" namespace SignalDef { // ********************************* // MINERvA CC1pi+/- signal definition (2015 release) // Note: There is a 2016 release which is different to this (listed below), but // it is CCNpi+ and has a different W cut -// Note2: The W cut is implemented in the class implementation in MINERvA/ -// rather than here so we can draw events that don't pass the W cut (W cut is -// 1.4 GeV) -// Could possibly be changed for slight speed increase since less events -// would be used // // MINERvA signal is slightly different to MiniBooNE // // Exactly one negative muon // Exactly one charged pion (both + and -); however, there is a Michel e- // requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be // applied) // Exactly 1 charged pion exits (so + and - charge), however, has Michel // electron requirement, so look for + only here? // No restriction on neutral pions or other mesons // MINERvA has unfolded and not unfolded muon phase space for 2015 +// W_TRUE < 1.4 GeV // // Possible issues with the data: // 1) pi- is allowed in signal even when Michel cut included; most pi- is // efficiency corrected in GENIE 2) There is a T_pi < 350 MeV cut coming from // requiring a stopping pion; this is efficiency corrected in GENIE 3) There is // a 1.5 < Enu < 10.0 cut in signal definition 4) There is an angular muon cut // which is sometimes efficiency corrected (why we have bool isRestricted below) // // Nice things: // Much data given: with and without muon angle cuts and with and without shape // only data + covariance // bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted) { // ********************************* // Signal is both pi+ and pi- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Allow pi+/pi- int piPDG[] = {211, -211}; int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // MINERvA released another CC1pi+ xsec without muon unfolding! // here the muon angle is < 20 degrees (seen in MINOS) TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; if (isRestricted) { double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20) return false; } - // Extract Hadronic Mass - double hadMass = FitUtils::Wrec(pnu, pmu); + double hadMass = 9999.99; // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition. #ifdef __GENIE_ENABLED__ if (event->fType == kGENIE) { EventRecord *gevent = static_cast(event->genie_event->event); const Interaction *interaction = gevent->Summary(); const Kinematics &kine = interaction->Kine(); double Ws = kine.W(true); - // std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " << - // kine.W(false) << std::endl; hadMass = Ws * 1000.0; } +#else + // Extract Hadronic Mass + // Cut is *INDEED* on Wtrue, not Wrec, so need to pass initial state nucleon too + // Either a proton or a nucleon + int nucPDG[] = {2212, 2112}; + // This won't work perfectly for 2p2h though + FitParticle *part = event->GetHMISParticle(nucPDG); + // There may not be an initial state nucleon if it's a coherent event + if (part == NULL) { + return 9999.999; + } + TLorentzVector pnuc = part->P4(); + hadMass = FitUtils::Wtrue(pnu, pmu, pnuc); + #endif if (hadMass > 1400.0) return false; return true; }; // Updated MINERvA 2017 Signal using Wexp and no restriction on angle bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax) { // Signal is both pi+ and pi- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Allow pi+/pi- int piPDG[] = {211, -211}; int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // Get Muon and Lepton Kinematics TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; // Extract Hadronic Mass + // This time it's Wrec, not Wtrue double hadMass = FitUtils::Wrec(pnu, pmu); // Cut on 2017 data is still 1.4 GeV if (hadMass > 1400.0) return false; return true; }; // ********************************* // MINERvA CCNpi+/- signal definition from 2016 publication // Different to CC1pi+/- listed above; additional has W < 1.8 GeV // // For notes on strangeness of signal definition, see CC1pip_MINERvA // // One negative muon // At least one charged pion // 1.5 < Enu < 10 // No restrictions on pi0 or other mesons or baryons // W_reconstructed (ignoring initial state motion) cut at 1.8 GeV // // Also writes number of pions (nPions) if studies on this want to be done... bool isCCNpip_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted, bool isWtrue) { // ********************************* // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; int nLeptons = event->NumFSLeptons(); // Write the number of pions to the measurement class... // Maybe better to just do that inside the class? int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions); // Check that there is a pion! if (nPions == 0) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // Need the muon and the neutrino to check angles and W TLorentzVector pnu = event->GetNeutrinoIn()->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; // MINERvA released some data with restricted muon angle // Here the muon angle is < 20 degrees (seen in MINOS) if (isRestricted) { double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.) return false; } // Lastly check the W cut (always at 1.8 GeV) double Wrec = FitUtils::Wrec(pnu, pmu) + 0.; // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition. if (isWtrue) { #ifdef __GENIE_ENABLED__ if (event->fType == kGENIE) { GHepRecord *ghep = static_cast(event->genie_event->event); const Interaction *interaction = ghep->Summary(); const Kinematics &kine = interaction->Kine(); double Ws = kine.W(true); Wrec = Ws * 1000.0; // Say Wrec is Ws } #endif } if (Wrec > 1800. || Wrec < 0.0) return false; return true; }; //******************************************************************** bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool fullphasespace) { //******************************************************************** if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; double ThetaMu = pnu.Vect().Angle(pmu.Vect()); double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true); // If Restricted phase space if (!fullphasespace && ThetaMu > 0.34906585) return false; // restrict energy range if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false; return true; }; //******************************************************************** bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool fullphasespace) { //******************************************************************** if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false; TLorentzVector pnu = event->GetHMISParticle(-14)->fP; TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; double ThetaMu = pnu.Vect().Angle(pmu.Vect()); double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true); // If Restricted phase space if (!fullphasespace && ThetaMu > 0.34906585) return false; // restrict energy range if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false; return true; } //******************************************************************** bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) { //******************************************************************** if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Need at least one muon if (event->NumFSParticle(13) < 1) return false; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pnu = event->GetHMISParticle(14)->fP; // Cut on muon angle greated than 20deg if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078) return false; // Cut on muon energy < 1.5 GeV if (pmu.E() / 1000.0 < 1.5) return false; return true; } // Used in 2014 muon+proton analysis // Events with muon angles up to 70 degrees // One right sign muon, at least one proton, no pions // proton kinetic energies greater than 100 MeV bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) { bool signal = (isCC0pi(event, 14, enumin, enumax) && // Require numu CC0pi event HasProtonKEAboveThreshold(event, 110.0) && // With proton above threshold (event->GetHMFSMuon())->P3().Angle((event->GetNeutrinoIn())->P3()) * 180. / M_PI < 70 // And muon within production angle ); return signal; } // 2015 analysis just asks for 1pi0 and no pi+/pi- bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax); return CC1pi0_anu; } // 2016 analysis just asks for 1pi0 and no other charged tracks. Half-open to // interpretation: we go with "charged tracks" meaning pions. You'll be forgiven // for thinking proton tracks should be included here too but we checked with // MINERvA bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax); /* // Additionally look for charged proton track // Comment: This is _NOT_ in the signal definition but was tested bool HasProton = event->HasFSParticle(2212); if (CC1pi0_anu) { if (!HasProton) { return true; } else { return false; } } else { return false; } */ return CC1pi0_anu; } // 2016 analysis just asks for 1pi0 and no other charged tracks bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_nu = SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax); return CC1pi0_nu; } //******************************************************************** bool isCC0pi_MINERvAPTPZ(FitEvent *event, int nuPDG, double emin, double emax) { //******************************************************************** // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false; // Make Angle Cut > 20.0 TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.0) return false; int genie_n_muons = 0; int genie_n_mesons = 0; int genie_n_heavy_baryons_plus_pi0s = 0; int genie_n_photons = 0; for (unsigned int i = 0; i < event->NParticles(); ++i) { FitParticle *p = event->GetParticle(i); if (p->Status() != kFinalState) continue; int pdg = p->fPID; double energy = p->fP.E(); if (pdg == 13) { genie_n_muons++; } else if (pdg == 22 && energy > 10.0) { genie_n_photons++; } else if (abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331) { genie_n_mesons++; } else if (pdg == 3112 || pdg == 3122 || pdg == 3212 || pdg == 3222 || pdg == 4112 || pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 411 || pdg == 421 || pdg == 111) { genie_n_heavy_baryons_plus_pi0s++; } } if (genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0) return true; + return false; } // ************************************************** // Section VI Event Selection of // https://journals.aps.org/prd/pdf/10.1103/PhysRevD.97.052002 Anti-neutrino // charged-current Post-FSI final states without // mesons, // prompt photons above nuclear deexcitation // energies heavy baryons protons above kinetic // energy of 120 MeV // Muon-neutrino angle of 20 degrees // Parallel muon momentum: 1.5 GeV < P|| < 15 GeV --- N.B. APPARENTLY NOT // INCLUDED, see below Transverse muon momentum: pT < 1.5 GeV --- N.B. // APPARENTLY NOT INCLUDED, see below bool isCC0pi_anti_MINERvAPTPZ(FitEvent *event, int nuPDG, double emin, double emax) { // ************************************************** // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false; TLorentzVector pnu = event->GetNeutrinoIn()->fP; TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; // Make Angle Cut > 20.0 double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.0) return false; // Heidi Schellman (schellmh@science.oregonstate.edu) assured me that the p_t // and p_z (or p_||) cuts aren't actually implemented as a signal definition: // they're only implemented in the binning for p_t and p_z (but not Q2QE and // EnuQE) /* // Cut on pT and pZ Double_t px = pmu.X()/1.E3; Double_t py = pmu.Y()/1.E3; Double_t pt = sqrt(px*px+py*py); // Don't want to assume the event generators all have neutrino coming along z // pz is muon momentum projected onto the neutrino direction Double_t pz = pmu.Vect().Dot(pnu.Vect()*(1.0/pnu.Vect().Mag()))/1.E3; if (pz > 15 || pz < 1.5) return false; if (pt > 1.5) return false; */ // Find if there are any protons above 120 MeV kinetic energy if (HasProtonKEAboveThreshold(event, 120.0)) return false; // Particle counters int genie_n_muons = 0; int genie_n_mesons = 0; int genie_n_heavy_baryons_plus_pi0s = 0; int genie_n_photons = 0; // Loop over the particles in the event and count them up for (unsigned int i = 0; i < event->NParticles(); ++i) { FitParticle *p = event->GetParticle(i); if (p->Status() != kFinalState) continue; int pdg = p->fPID; double energy = p->fP.E(); // Any charged muons if (abs(pdg) == 13) { genie_n_muons++; // De-excitation photons } else if (pdg == 22 && energy > 10.0) { genie_n_photons++; // Mesons } else if (abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331) { genie_n_mesons++; // Heavy baryons and pi0s } else if (abs(pdg) == 3112 || abs(pdg) == 3122 || abs(pdg) == 3212 || abs(pdg) == 3222 || abs(pdg) == 4112 || abs(pdg) == 4122 || abs(pdg) == 4212 || abs(pdg) == 4222 || abs(pdg) == 411 || abs(pdg) == 421 || abs(pdg) == 111) { genie_n_heavy_baryons_plus_pi0s++; } } // Look for one muon with no mesons, heavy baryons or deexcitation photons if (genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0) return true; return false; } // MINERvA CC0pi transverse variables signal defintion bool isCC0piNp_MINERvA_STV(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; // Muon momentum cuts if (pmu.Vect().Mag() < 1500 || pmu.Vect().Mag() > 10000) return false; // Muon angle cuts if (pmu.Vect().Angle(pnu.Vect()) > (M_PI / 180.0) * 20.0) return false; const double ctheta_cut = cos((M_PI / 180.0) * 70.0); //Did you find a proton in the PS? if (MINERvAUtils::GetProtonInRange(event, 450, 1200, ctheta_cut).E() == 0) return false; return true; }; } // namespace SignalDef diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx index c08c342..912470e 100644 --- a/src/Utils/PlotUtils.cxx +++ b/src/Utils/PlotUtils.cxx @@ -1,1183 +1,1183 @@ // Copyright 2016-2021 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 "PlotUtils.h" #include "FitEvent.h" #include "StatUtils.h" // MOVE TO GENERAL UTILS? bool PlotUtils::CheckObjectWithName(TFile *inFile, std::string substring) { TIter nextkey(inFile->GetListOfKeys()); TKey *key; while ((key = (TKey *)nextkey())) { std::string test(key->GetName()); if (test.find(substring) != std::string::npos) return true; } return false; }; // MOVE TO GENERAL UTILS? std::string PlotUtils::GetObjectWithName(TFile *inFile, std::string substring) { TIter nextkey(inFile->GetListOfKeys()); TKey *key; std::string output = ""; while ((key = (TKey *)nextkey())) { std::string test(key->GetName()); if (test.find(substring) != std::string::npos) output = test; } return output; }; void PlotUtils::CreateNeutModeArray(TH1 *hist, TH1 *neutarray[]) { for (int i = 0; i < 60; i++) { neutarray[i] = (TH1 *)hist->Clone(Form("%s_NMODE_%i", hist->GetName(), i)); } return; }; void PlotUtils::DeleteNeutModeArray(TH1 *neutarray[]) { for (int i = 0; i < 60; i++) { delete neutarray[i]; } return; }; void PlotUtils::FillNeutModeArray(TH1D *hist[], int mode, double xval, double weight) { if (abs(mode) > 60) return; hist[abs(mode)]->Fill(xval, weight); return; }; void PlotUtils::FillNeutModeArray(TH2D *hist[], int mode, double xval, double yval, double weight) { if (abs(mode) > 60) return; hist[abs(mode)]->Fill(xval, yval, weight); return; }; THStack PlotUtils::GetNeutModeStack(std::string title, TH1 *ModeStack[], int option) { (void)option; THStack allmodes = THStack(title.c_str(), title.c_str()); for (int i = 0; i < 60; i++) { allmodes.Add(ModeStack[i]); } // Credit to Clarence for copying all this out. // CC ModeStack[1]->SetTitle("CCQE"); ModeStack[1]->SetFillColor(kBlue); // ModeStack[1]->SetFillStyle(3444); ModeStack[1]->SetLineColor(kBlue); ModeStack[2]->SetTitle("2p/2h Nieves"); ModeStack[2]->SetFillColor(kRed); // ModeStack[2]->SetFillStyle(3344); ModeStack[2]->SetLineColor(kRed); // ModeStack[11]->SetTitle("#it{#nu + p #rightarrow l^{-} + p + #pi^{+}}"); ModeStack[11]->SetTitle("CC1#pi^{+} on p"); ModeStack[11]->SetFillColor(kGreen); // ModeStack[11]->SetFillStyle(3004); ModeStack[11]->SetLineColor(kGreen); // ModeStack[12]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #pi^{0}}"); ModeStack[12]->SetTitle("CC1#pi^{0} on n"); ModeStack[12]->SetFillColor(kGreen + 3); // ModeStack[12]->SetFillStyle(3005); ModeStack[12]->SetLineColor(kGreen); // ModeStack[13]->SetTitle("#it{#nu + n #rightarrow l^{-} + n + #pi^{+}}"); ModeStack[13]->SetTitle("CC1#pi^{+} on n"); ModeStack[13]->SetFillColor(kGreen - 2); // ModeStack[13]->SetFillStyle(3004); ModeStack[13]->SetLineColor(kGreen); ModeStack[16]->SetTitle("CC coherent"); - ModeStack[16]->SetFillColor(kBlue); + ModeStack[16]->SetFillColor(kBlue-6); // ModeStack[16]->SetFillStyle(3644); - ModeStack[16]->SetLineColor(kBlue); + ModeStack[16]->SetLineColor(kBlue-6); // ModeStack[17]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #gamma}"); ModeStack[17]->SetTitle("CC1#gamma"); ModeStack[17]->SetFillColor(kMagenta); // ModeStack[17]->SetFillStyle(3001); ModeStack[17]->SetLineColor(kMagenta); ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)"); ModeStack[21]->SetFillColor(kYellow); // ModeStack[21]->SetFillStyle(3005); ModeStack[21]->SetLineColor(kYellow); // ModeStack[22]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #eta^{0}}"); ModeStack[22]->SetTitle("CC1#eta^{0} on n"); ModeStack[22]->SetFillColor(kYellow - 2); // ModeStack[22]->SetFillStyle(3013); ModeStack[22]->SetLineColor(kYellow - 2); // ModeStack[23]->SetTitle("#it{#nu + n #rightarrow l^{-} + #Lambda + // K^{+}}"); ModeStack[23]->SetTitle("CC1#Labda1K^{+}"); ModeStack[23]->SetFillColor(kYellow - 6); // ModeStack[23]->SetFillStyle(3013); ModeStack[23]->SetLineColor(kYellow - 6); ModeStack[26]->SetTitle("DIS (W > 2.0)"); - ModeStack[26]->SetFillColor(kRed); + ModeStack[26]->SetFillColor(kRed-6); // ModeStack[26]->SetFillStyle(3006); - ModeStack[26]->SetLineColor(kRed); + ModeStack[26]->SetLineColor(kRed-6); // NC // ModeStack[31]->SetTitle("#it{#nu + n #rightarrow #nu + n + #pi^{0}}"); ModeStack[31]->SetTitle("NC1#pi^{0} on n"); - ModeStack[31]->SetFillColor(kBlue); + ModeStack[31]->SetFillColor(kBlue-4); // ModeStack[31]->SetFillStyle(3004); - ModeStack[31]->SetLineColor(kBlue); + ModeStack[31]->SetLineColor(kBlue-4); // ModeStack[32]->SetTitle("#it{#nu + p #rightarrow #nu + p + #pi^{0}}"); ModeStack[32]->SetTitle("NC1#pi^{0} on p"); ModeStack[32]->SetFillColor(kBlue + 3); // ModeStack[32]->SetFillStyle(3004); ModeStack[32]->SetLineColor(kBlue + 3); // ModeStack[33]->SetTitle("#it{#nu + n #rightarrow #nu + p + #pi^{-}}"); ModeStack[33]->SetTitle("NC1#pi^{-} on n"); ModeStack[33]->SetFillColor(kBlue - 2); // ModeStack[33]->SetFillStyle(3005); ModeStack[33]->SetLineColor(kBlue - 2); // ModeStack[34]->SetTitle("#it{#nu + p #rightarrow #nu + n + #pi^{+}}"); ModeStack[34]->SetTitle("NC1#pi^{+} on p"); ModeStack[34]->SetFillColor(kBlue - 8); // ModeStack[34]->SetFillStyle(3005); ModeStack[34]->SetLineColor(kBlue - 8); ModeStack[36]->SetTitle("NC Coherent"); ModeStack[36]->SetFillColor(kBlue + 8); // ModeStack[36]->SetFillStyle(3644); ModeStack[36]->SetLineColor(kBlue + 8); // ModeStack[38]->SetTitle("#it{#nu + n #rightarrow #nu + n + #gamma}"); ModeStack[38]->SetTitle("NC1#gamma on n"); ModeStack[38]->SetFillColor(kMagenta); // ModeStack[38]->SetFillStyle(3001); ModeStack[38]->SetLineColor(kMagenta); // ModeStack[39]->SetTitle("#it{#nu + p #rightarrow #nu + p + #gamma}"); ModeStack[39]->SetTitle("NC1#gamma on p"); ModeStack[39]->SetFillColor(kMagenta - 10); // ModeStack[39]->SetFillStyle(3001); ModeStack[39]->SetLineColor(kMagenta - 10); ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)"); ModeStack[41]->SetFillColor(kBlue - 10); // ModeStack[41]->SetFillStyle(3005); ModeStack[41]->SetLineColor(kBlue - 10); // ModeStack[42]->SetTitle("#it{#nu + n #rightarrow #nu + n + #eta^{0}}"); ModeStack[42]->SetTitle("NC1#eta^{0} on n"); ModeStack[42]->SetFillColor(kYellow - 2); // ModeStack[42]->SetFillStyle(3013); ModeStack[42]->SetLineColor(kYellow - 2); // ModeStack[43]->SetTitle("#it{#nu + p #rightarrow #nu + p + #eta^{0}}"); ModeStack[43]->SetTitle("NC1#eta^{0} on p"); ModeStack[43]->SetFillColor(kYellow - 4); // ModeStack[43]->SetFillStyle(3013); ModeStack[43]->SetLineColor(kYellow - 4); // ModeStack[44]->SetTitle("#it{#nu + n #rightarrow #nu + #Lambda + K^{0}}"); ModeStack[44]->SetTitle("NC1#Lambda1K^{0} on n"); ModeStack[44]->SetFillColor(kYellow - 6); // ModeStack[44]->SetFillStyle(3014); ModeStack[44]->SetLineColor(kYellow - 6); // ModeStack[45]->SetTitle("#it{#nu + p #rightarrow #nu + #Lambda + K^{+}}"); ModeStack[45]->SetTitle("NC1#Lambda1K^{+}"); ModeStack[45]->SetFillColor(kYellow - 10); // ModeStack[45]->SetFillStyle(3014); ModeStack[45]->SetLineColor(kYellow - 10); ModeStack[46]->SetTitle("DIS (W > 2.0)"); - ModeStack[46]->SetFillColor(kRed); + ModeStack[46]->SetFillColor(kRed-6); // ModeStack[46]->SetFillStyle(3006); - ModeStack[46]->SetLineColor(kRed); + ModeStack[46]->SetLineColor(kRed-6); // ModeStack[51]->SetTitle("#it{#nu + p #rightarrow #nu + p}"); ModeStack[51]->SetTitle("NC on p"); ModeStack[51]->SetFillColor(kBlack); // ModeStack[51]->SetFillStyle(3444); ModeStack[51]->SetLineColor(kBlack); // ModeStack[52]->SetTitle("#it{#nu + n #rightarrow #nu + n}"); ModeStack[52]->SetTitle("NC on n"); ModeStack[52]->SetFillColor(kGray); // ModeStack[52]->SetFillStyle(3444); ModeStack[52]->SetLineColor(kGray); return allmodes; }; TLegend PlotUtils::GenerateStackLegend(THStack stack, int xlow, int ylow, int xhigh, int yhigh) { TLegend leg = TLegend(xlow, ylow, xhigh, yhigh); TObjArray *histarray = stack.GetStack(); int nhist = histarray->GetEntries(); for (int i = 0; i < nhist; i++) { TH1 *hist = (TH1 *)(histarray->At(i)); leg.AddEntry((hist), ((TH1 *)histarray->At(i))->GetTitle(), "fl"); } leg.SetName(Form("%s_LEG", stack.GetName())); return leg; }; void PlotUtils::ScaleNeutModeArray(TH1 *hist[], double factor, std::string option) { for (int i = 0; i < 60; i++) { if (hist[i]) hist[i]->Scale(factor, option.c_str()); } return; }; void PlotUtils::ResetNeutModeArray(TH1 *hist[]) { for (int i = 0; i < 60; i++) { if (hist[i]) hist[i]->Reset(); } return; }; //******************************************************************** // This assumes the Enu axis is the x axis, as is the case for MiniBooNE 2D // distributions void PlotUtils::FluxUnfoldedScaling(TH2D *fMCHist, TH1D *fhist, TH1D *ehist, double scalefactor) { //******************************************************************** // Make clones to avoid changing stuff TH1D *eventhist = (TH1D *)ehist->Clone(); TH1D *fFluxHist = (TH1D *)fhist->Clone(); // Undo width integral in SF fMCHist->Scale(scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width")); // Standardise The Flux eventhist->Scale(1.0 / fFluxHist->Integral()); fFluxHist->Scale(1.0 / fFluxHist->Integral()); // Do interpolation for 2D plots? // fFluxHist = PlotUtils::InterpolateFineHistogram(fFluxHist,100,"width"); // eventhist = PlotUtils::InterpolateFineHistogram(eventhist,100,"width"); // eventhist->Scale(1.0/fFluxHist->Integral()); // fFluxHist->Scale(1.0/fFluxHist->Integral()); // Scale fMCHist by eventhist integral fMCHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1)); // Find which axis is the Enu axis bool EnuOnXaxis = false; std::string xaxis = fMCHist->GetXaxis()->GetTitle(); if (xaxis.find("E") != std::string::npos && xaxis.find("nu") != std::string::npos) EnuOnXaxis = true; std::string yaxis = fMCHist->GetYaxis()->GetTitle(); if (yaxis.find("E") != std::string::npos && xaxis.find("nu") != std::string::npos) { // First check that xaxis didn't also find Enu if (EnuOnXaxis) { NUIS_ERR(FTL, fMCHist->GetTitle() << " error:"); NUIS_ERR(FTL, "Found Enu in xaxis title: " << xaxis); NUIS_ERR(FTL, "AND"); NUIS_ERR(FTL, "Found Enu in yaxis title: " << yaxis); NUIS_ABORT("Enu on x and Enu on y flux unfolded scaling isn't " "implemented, please modify " << __FILE__ << ":" << __LINE__); } EnuOnXaxis = false; } // Now Get a flux PDF assuming X axis is Enu TH1D *pdfflux = NULL; // If xaxis is Enu if (EnuOnXaxis) pdfflux = (TH1D *)fMCHist->ProjectionX()->Clone(); // If yaxis is Enu else pdfflux = (TH1D *)fMCHist->ProjectionY()->Clone(); // pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str()); pdfflux->Reset(); // Awful MiniBooNE Check for the time being // Needed because the flux is in GeV whereas the measurement is in MeV bool ismb = std::string(fMCHist->GetName()).find("MiniBooNE") != std::string::npos; for (int i = 0; i < pdfflux->GetNbinsX(); i++) { double Ml = pdfflux->GetXaxis()->GetBinLowEdge(i + 1); double Mh = pdfflux->GetXaxis()->GetBinLowEdge(i + 2); // double Mc = pdfflux->GetXaxis()->GetBinCenter(i+1); // double Mw = pdfflux->GetBinWidth(i+1); double fluxint = 0.0; // Scaling to match flux for MB if (ismb) { Ml /= 1.E3; Mh /= 1.E3; // Mc /= 1.E3; // Mw /= 1.E3; } for (int j = 0; j < fFluxHist->GetNbinsX(); j++) { // double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1); double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1); double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2); double Fe = fFluxHist->GetBinContent(j + 1); double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1); if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; } else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; } else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; } else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; } else { continue; } } pdfflux->SetBinContent(i + 1, fluxint); } // Then finally divide by the bin-width in for (int i = 0; i < fMCHist->GetNbinsX(); i++) { for (int j = 0; j < fMCHist->GetNbinsY(); j++) { if (pdfflux->GetBinContent(i + 1) == 0.0) continue; // Different scaling depending on if Enu is on x or y axis double scaling = 1.0; // If Enu is on the x-axis, we want the ith entry of the flux // And to divide by the bin width of the jth bin if (EnuOnXaxis) { double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) - fMCHist->GetYaxis()->GetBinLowEdge(j + 1); scaling = pdfflux->GetBinContent(i + 1) * binWidth; } else { double binWidth = fMCHist->GetXaxis()->GetBinLowEdge(i + 2) - fMCHist->GetXaxis()->GetBinLowEdge(i + 1); scaling = pdfflux->GetBinContent(j + 1) * binWidth; } // fMCHist->SetBinContent(i + 1, j + 1, // fMCHist->GetBinContent(i + 1, j + 1) / // pdfflux->GetBinContent(i + 1) / binWidth); // fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) / // pdfflux->GetBinContent(i + 1) / // binWidth); fMCHist->SetBinContent(i + 1, j + 1, fMCHist->GetBinContent(i + 1, j + 1) / scaling); fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) / scaling); } } delete eventhist; delete fFluxHist; }; TH1D *PlotUtils::InterpolateFineHistogram(TH1D *hist, int res, std::string opt) { int nbins = hist->GetNbinsX(); double elow = hist->GetXaxis()->GetBinLowEdge(1); double ehigh = hist->GetXaxis()->GetBinLowEdge(nbins + 1); bool width = true; // opt.find("width") != std::string::npos; TH1D *fine = new TH1D("fine", "fine", nbins * res, elow, ehigh); TGraph *temp = new TGraph(); for (int i = 0; i < nbins; i++) { double E = hist->GetXaxis()->GetBinCenter(i + 1); double C = hist->GetBinContent(i + 1); double W = hist->GetXaxis()->GetBinWidth(i + 1); if (!width) W = 1.0; if (W != 0.0) temp->SetPoint(temp->GetN(), E, C / W); } for (int i = 0; i < fine->GetNbinsX(); i++) { double E = fine->GetXaxis()->GetBinCenter(i + 1); double W = fine->GetBinWidth(i + 1); if (!width) W = 1.0; fine->SetBinContent(i + 1, temp->Eval(E, 0, "S") * W); } fine->Scale(hist->Integral(1, hist->GetNbinsX() + 1) / fine->Integral(1, fine->GetNbinsX() + 1)); // std::cout << "Interpolation Difference = " //<< fine->Integral(1, fine->GetNbinsX() + 1) << "/" //<< hist->Integral(1, hist->GetNbinsX() + 1) << std::endl; return fine; } //******************************************************************** // This interpolates the flux by a TGraph instead of requiring the flux and MC // flux to have the same binning void PlotUtils::FluxUnfoldedScaling(TH1D *mcHist, TH1D *fhist, TH1D *ehist, double scalefactor, int nevents) { //******************************************************************** TH1D *eventhist = (TH1D *)ehist->Clone(); TH1D *fFluxHist = (TH1D *)fhist->Clone(); std::string name = std::string(mcHist->GetName()); if (FitPar::Config().GetParB("save_flux_debug")) { mcHist->Write((name + "_UNF_MC").c_str()); fFluxHist->Write((name + "_UNF_FLUX").c_str()); eventhist->Write((name + "_UNF_EVT").c_str()); TH1D *scalehist = new TH1D("scalehist", "scalehist", 1, 0.0, 1.0); scalehist->SetBinContent(1, scalefactor); scalehist->SetBinContent(2, nevents); scalehist->Write((name + "_UNF_SCALE").c_str()); } // Undo width integral in SF mcHist->Scale(scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width")); // Standardise The Flux eventhist->Scale(1.0 / fFluxHist->Integral()); fFluxHist->Scale(1.0 / fFluxHist->Integral()); // Scale mcHist by eventhist integral mcHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1)); // Now Get a flux PDF TH1D *pdfflux = (TH1D *)mcHist->Clone(); pdfflux->Reset(); for (int i = 0; i < mcHist->GetNbinsX(); i++) { double Ml = mcHist->GetXaxis()->GetBinLowEdge(i + 1); double Mh = mcHist->GetXaxis()->GetBinLowEdge(i + 2); // double Mc = mcHist->GetXaxis()->GetBinCenter(i+1); // double Me = mcHist->GetBinContent(i+1); // double Mw = mcHist->GetBinWidth(i+1); double fluxint = 0.0; for (int j = 0; j < fFluxHist->GetNbinsX(); j++) { // double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1); double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1); double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2); double Fe = fFluxHist->GetBinContent(j + 1); double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1); if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; } else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; } else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; } else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; } else { continue; } } pdfflux->SetBinContent(i + 1, fluxint); } if (FitPar::Config().GetParB("save_flux_debug")) { pdfflux->Write((name + "_UNF_SCALEHIST").c_str()); } // Scale MC hist by pdfflux for (int i = 0; i < mcHist->GetNbinsX(); i++) { if (pdfflux->GetBinContent(i + 1) == 0.0) continue; mcHist->SetBinContent(i + 1, mcHist->GetBinContent(i + 1) / pdfflux->GetBinContent(i + 1)); mcHist->SetBinError(i + 1, mcHist->GetBinError(i + 1) / pdfflux->GetBinContent(i + 1)); } delete eventhist; delete fFluxHist; }; // MOVE TO GENERAL UTILS //******************************************************************** void PlotUtils::Set2DHistFromText(std::string dataFile, TH2 *hist, double norm, bool skipbins) { //******************************************************************** std::string line; std::ifstream data(dataFile.c_str(), std::ifstream::in); int yBin = 0; while (std::getline(data >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToDbl(line, " "); // Loop over entries and insert them into the histogram for (uint xBin = 0; xBin < entries.size(); xBin++) { if (!skipbins || entries[xBin] != -1.0) hist->SetBinContent(xBin + 1, yBin + 1, entries[xBin] * norm); } yBin++; } return; } // MOVE TO GENERAL UTILS TH1D *PlotUtils::GetTH1DFromFile(std::string dataFile, std::string title, std::string fPlotTitles, std::string alt_name) { TH1D *tempPlot; // If format is a root file if (dataFile.find(".root") != std::string::npos) { TFile *temp_infile = new TFile(dataFile.c_str(), "READ"); tempPlot = (TH1D *)temp_infile->Get(title.c_str()); tempPlot->SetDirectory(0); temp_infile->Close(); delete temp_infile; // Else its a space separated txt file } else { // Make a TGraph Errors TGraphErrors *gr = new TGraphErrors(dataFile.c_str(), "%lg %lg %lg"); if (gr->IsZombie()) { NUIS_ABORT( dataFile << " is a zombie and could not be read. Are you sure it exists?" << std::endl); } double *bins = gr->GetX(); double *values = gr->GetY(); double *errors = gr->GetEY(); int npoints = gr->GetN(); // Fill the histogram from it tempPlot = new TH1D(title.c_str(), title.c_str(), npoints - 1, bins); for (int i = 0; i < npoints; ++i) { tempPlot->SetBinContent(i + 1, values[i]); // If only two columns are present in the input file, use the sqrt(values) // as the error equivalent to assuming that the error is statistical. Also // check that we're looking at an event rate rather than a cross section if (!errors[i] && values[i] > 1E-30) { tempPlot->SetBinError(i + 1, sqrt(values[i])); } else { tempPlot->SetBinError(i + 1, errors[i]); } } delete gr; } // Allow alternate naming for root files if (!alt_name.empty()) { tempPlot->SetNameTitle(alt_name.c_str(), alt_name.c_str()); } // Allow alternate axis titles if (!fPlotTitles.empty()) { tempPlot->SetNameTitle( tempPlot->GetName(), (std::string(tempPlot->GetTitle()) + fPlotTitles).c_str()); } return tempPlot; }; TH1D *PlotUtils::GetRatioPlot(TH1D *hist1, TH1D *hist2, TH1D *new_hist) { // If the hist to save into doesn't exist, make copy of first hist if (!new_hist) new_hist = (TH1D *)hist1->Clone(); // Do bins and errors ourselves as scales can go awkward for (int i = 0; i < new_hist->GetNbinsX(); i++) { double binVal = 0; double binErr = 0; if (hist2->GetBinContent(i+1) && hist1->GetBinContent(i+1)) { binVal = hist1->GetBinContent(i+1)/hist2->GetBinContent(i+1); double fractErr1 = hist1->GetBinError(i+1)/hist1->GetBinContent(i+1); double fractErr2 = hist2->GetBinError(i+1)/hist2->GetBinContent(i+1); binErr = binVal * sqrt(fractErr1*fractErr1 + fractErr2*fractErr2); } new_hist->SetBinContent(i+1, binVal); new_hist->SetBinError(i+1, binErr); } return new_hist; }; TH1D *PlotUtils::GetRenormalisedPlot(TH1D *hist1, TH1D *hist2) { // make copy of first hist TH1D *new_hist = (TH1D *)hist1->Clone(); if (hist1->Integral("width") == 0 or hist2->Integral("width") == 0) { new_hist->Reset(); return new_hist; } Double_t scaleF = hist2->Integral("width") / hist1->Integral("width"); new_hist->Scale(scaleF); return new_hist; }; TH1D *PlotUtils::GetShapePlot(TH1D *hist1) { // make copy of first hist TH1D *new_hist = (TH1D *)hist1->Clone(); if (hist1->Integral("width") == 0) { new_hist->Reset(); return new_hist; } Double_t scaleF1 = 1.0 / hist1->Integral("width"); new_hist->Scale(scaleF1); return new_hist; }; TH1D *PlotUtils::GetShapeRatio(TH1D *hist1, TH1D *hist2) { TH1D *new_hist1 = GetShapePlot(hist1); TH1D *new_hist2 = GetShapePlot(hist2); // Do bins and errors ourselves as scales can go awkward for (int i = 0; i < new_hist1->GetNbinsX(); i++) { if (hist2->GetBinContent(i + 1) == 0) { new_hist1->SetBinContent(i + 1, 0.0); } new_hist1->SetBinContent(i + 1, new_hist1->GetBinContent(i + 1) / new_hist2->GetBinContent(i + 1)); new_hist1->SetBinError(i + 1, new_hist1->GetBinError(i + 1) / new_hist2->GetBinContent(i + 1)); } delete new_hist2; return new_hist1; }; TH2D *PlotUtils::GetCovarPlot(TMatrixDSym *cov, std::string name, std::string title) { TH2D *CovarPlot; if (cov) CovarPlot = new TH2D((*cov)); else CovarPlot = new TH2D(name.c_str(), title.c_str(), 1, 0, 1, 1, 0, 1); CovarPlot->SetName(name.c_str()); CovarPlot->SetTitle(title.c_str()); return CovarPlot; } TH2D *PlotUtils::GetFullCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_COV", name + "_COV;Bins;Bins;Covariance (#times10^{-76})"); } TH2D *PlotUtils::GetInvCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_INVCOV", name + "_INVCOV;Bins;Bins;Inv. Covariance (#times10^{-76})"); } TH2D *PlotUtils::GetDecompCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_DECCOV", name + "_DECCOV;Bins;Bins;Decomp Covariance (#times10^{-76})"); } TH1D *PlotUtils::GetTH1DFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); TH1D *tempHist = (TH1D *)rootHistFile->Get(name.c_str())->Clone(); if (tempHist == NULL) { NUIS_ABORT("Could not find distribution " << name << " in file " << file); } tempHist->SetDirectory(0); rootHistFile->Close(); return tempHist; } TH2D *PlotUtils::GetTH2DFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); TH2D *tempHist = (TH2D *)rootHistFile->Get(name.c_str())->Clone(); tempHist->SetDirectory(0); rootHistFile->Close(); delete rootHistFile; return tempHist; } TH1 *PlotUtils::GetTH1FromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << file << "\"."); } TH1 *tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!tempHist) { NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" << file << "\"."); } tempHist->SetDirectory(0); rootHistFile->Close(); delete rootHistFile; return tempHist; } TGraph *PlotUtils::GetTGraphFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TDirectory *olddir = gDirectory; TFile *rootHistFile = new TFile(file.c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << file << "\"."); } TDirectory *newdir = gDirectory; TGraph *temp = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!temp) { NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" << file << "\"."); } newdir->Remove(temp); olddir->Append(temp); rootHistFile->Close(); olddir->cd(); return temp; } /// Returns a vector of named TH1*s found in a single input file. /// /// Expects a descriptor like: file.root[hist1|hist2|...] std::vector PlotUtils::GetTH1sFromRootFile(std::string const &descriptor) { std::vector descriptors = GeneralUtils::ParseToStr(descriptor, ","); std::vector hists; for (size_t d_it = 0; d_it < descriptors.size(); ++d_it) { std::string &d = descriptors[d_it]; std::vector fname = GeneralUtils::ParseToStr(d, "["); if (!fname.size() || !fname[0].length()) { NUIS_ABORT("Couldn't find input file when attempting to parse : \"" << d << "\". Expected input.root[hist1|hist2|...]."); } if (fname[1][fname[1].length() - 1] == ']') { fname[1] = fname[1].substr(0, fname[1].length() - 1); } std::vector histnames = GeneralUtils::ParseToStr(fname[1], "|"); if (!histnames.size()) { NUIS_ABORT( "Couldn't find any histogram name specifiers when attempting to " "parse " ": \"" << fname[1] << "\". Expected hist1|hist2|..."); } TFile *rootHistFile = new TFile(fname[0].c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << fname[0] << "\"."); } for (size_t i = 0; i < histnames.size(); ++i) { TH1 *tempHist = dynamic_cast(rootHistFile->Get(histnames[i].c_str())->Clone()); if (!tempHist) { NUIS_ABORT("Couldn't retrieve: \"" << histnames[i] << "\" from root file: \"" << fname[0] << "\"."); } tempHist->SetDirectory(0); hists.push_back(tempHist); } rootHistFile->Close(); } return hists; } // Create an array from an input file std::vector PlotUtils::GetArrayFromTextFile(std::string DataFile) { std::string line; std::ifstream data(DataFile.c_str(), std::ifstream::in); // Get first line std::getline(data >> std::ws, line, '\n'); // Convert from a string into a vector of double std::vector entries = GeneralUtils::ParseToDbl(line, " "); return entries; } // Get a 2D array from a text file std::vector > PlotUtils::Get2DArrayFromTextFile(std::string DataFile) { std::string line; std::vector > DataArray; std::ifstream data(DataFile.c_str(), std::ifstream::in); while (std::getline(data >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToDbl(line, " "); DataArray.push_back(entries); } return DataArray; } TH2D *PlotUtils::GetTH2DFromTextFile(std::string data, std::string binx, std::string biny) { // First read in the binning // Array of x binning std::vector xbins = GetArrayFromTextFile(binx); // Array of y binning std::vector ybins = GetArrayFromTextFile(biny); // Read in the data std::vector > Data = Get2DArrayFromTextFile(data); // And finally fill the data TH2D *DataPlot = new TH2D("TempHist", "TempHist", xbins.size() - 1, &xbins[0], ybins.size() - 1, &ybins[0]); int nBinsX = 0; int nBinsY = 0; for (std::vector >::iterator it = Data.begin(); it != Data.end(); ++it) { nBinsX++; // Get the inner vector std::vector temp = *it; // Save the previous number[of bins to make sure it's uniform binning int oldBinsY = nBinsY; // Reset the counter nBinsY = 0; for (std::vector::iterator jt = temp.begin(); jt != temp.end(); ++jt) { nBinsY++; DataPlot->SetBinContent(nBinsX, nBinsY, *jt); DataPlot->SetBinError(nBinsX, nBinsY, 0.0); } if (oldBinsY > 0 && oldBinsY != nBinsY) { NUIS_ERR(FTL, "Found non-uniform y-binning in " << data); NUIS_ERR(FTL, "Previous slice: " << oldBinsY); NUIS_ERR(FTL, "Current slice: " << nBinsY); NUIS_ABORT("Non-uniform binning is not supported in " "PlotUtils::GetTH2DFromTextFile"); } } // Check x bins if (size_t(nBinsX + 1) != xbins.size()) { NUIS_ERR(FTL, "Number of x bins in data histogram does not match the binning " "histogram!"); NUIS_ERR( FTL, "Are they the wrong way around (i.e. xbinning should be ybinning)?"); NUIS_ERR(FTL, "Data: " << nBinsX); NUIS_ABORT("From " << binx << " binning: " << xbins.size()); } // Check y bins if (size_t(nBinsY + 1) != ybins.size()) { NUIS_ERR(FTL, "Number of y bins in data histogram does not match the binning " "histogram!"); NUIS_ERR( FTL, "Are they the wrong way around (i.e. xbinning should be ybinning)?"); NUIS_ERR(FTL, "Data: " << nBinsY); NUIS_ABORT("From " << biny << " binning: " << ybins.size()); } return DataPlot; } TH1D *PlotUtils::GetSliceY(TH2D *Hist, int SliceNo) { TH1D *Slice = Hist->ProjectionX(Form("%s_SLICEY%i", Hist->GetName(), SliceNo), SliceNo, SliceNo, "e"); Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetYaxis()->GetTitle(), Hist->GetYaxis()->GetBinLowEdge(SliceNo), Hist->GetYaxis()->GetBinLowEdge(SliceNo + 1))); Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle()); return Slice; } TH1D *PlotUtils::GetSliceX(TH2D *Hist, int SliceNo) { TH1D *Slice = Hist->ProjectionY(Form("%s_SLICEX%i", Hist->GetName(), SliceNo), SliceNo, SliceNo, "e"); Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetXaxis()->GetTitle(), Hist->GetXaxis()->GetBinLowEdge(SliceNo), Hist->GetXaxis()->GetBinLowEdge(SliceNo + 1))); Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle()); return Slice; } void PlotUtils::AddNeutModeArray(TH1D *hist1[], TH1D *hist2[], double scaling) { for (int i = 0; i < 60; i++) { if (!hist2[i]) continue; if (!hist1[i]) continue; hist1[i]->Add(hist2[i], scaling); } return; } void PlotUtils::ScaleToData(TH1D *data, TH1D *mc, TH1I *mask) { double scaleF = GetDataMCRatio(data, mc, mask); mc->Scale(scaleF); return; } void PlotUtils::MaskBins(TH1D *hist, TH1I *mask) { for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1) <= 0.5) continue; hist->SetBinContent(i + 1, 0.0); hist->SetBinError(i + 1, 0.0); NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " to 0.0 +- 0.0"); } return; } void PlotUtils::MaskBins(TH2D *hist, TH2I *mask) { for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1) <= 0.5) continue; hist->SetBinContent(i + 1, j + 1, 0.0); hist->SetBinError(i + 1, j + 1, 0.0); NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " " << j + 1 << " to 0.0 +- 0.0"); } } return; } double PlotUtils::GetDataMCRatio(TH1D *data, TH1D *mc, TH1I *mask) { double rat = 1.0; TH1D *newmc = (TH1D *)mc->Clone(); TH1D *newdt = (TH1D *)data->Clone(); if (mask) { MaskBins(newmc, mask); MaskBins(newdt, mask); } rat = newdt->Integral() / newmc->Integral(); return rat; } TH2D *PlotUtils::GetCorrelationPlot(TH2D *cov, std::string name) { TH2D *cor = (TH2D *)cov->Clone(); cor->Reset(); for (int i = 0; i < cov->GetNbinsX(); i++) { for (int j = 0; j < cov->GetNbinsY(); j++) { if (cov->GetBinContent(i + 1, i + 1) != 0.0 and cov->GetBinContent(j + 1, j + 1) != 0.0) cor->SetBinContent(i + 1, j + 1, cov->GetBinContent(i + 1, j + 1) / (sqrt(cov->GetBinContent(i + 1, i + 1) * cov->GetBinContent(j + 1, j + 1)))); } } if (!name.empty()) { cor->SetNameTitle(name.c_str(), (name + ";;correlation").c_str()); } cor->SetMinimum(-1); cor->SetMaximum(1); return cor; } TH2D *PlotUtils::GetDecompPlot(TH2D *cov, std::string name) { TMatrixDSym *covarmat = new TMatrixDSym(cov->GetNbinsX()); for (int i = 0; i < cov->GetNbinsX(); i++) for (int j = 0; j < cov->GetNbinsY(); j++) (*covarmat)(i, j) = cov->GetBinContent(i + 1, j + 1); TMatrixDSym *decompmat = StatUtils::GetDecomp(covarmat); TH2D *dec = (TH2D *)cov->Clone(); for (int i = 0; i < cov->GetNbinsX(); i++) for (int j = 0; j < cov->GetNbinsY(); j++) dec->SetBinContent(i + 1, j + 1, (*decompmat)(i, j)); delete covarmat; delete decompmat; dec->SetNameTitle(name.c_str(), (name + ";;;decomposition").c_str()); return dec; } TH2D *PlotUtils::MergeIntoTH2D(TH1D *xhist, TH1D *yhist, std::string zname) { std::vector xedges, yedges; for (int i = 0; i < xhist->GetNbinsX() + 2; i++) { xedges.push_back(xhist->GetXaxis()->GetBinLowEdge(i + 1)); } for (int i = 0; i < yhist->GetNbinsX() + 2; i++) { yedges.push_back(yhist->GetXaxis()->GetBinLowEdge(i + 1)); } int nbinsx = xhist->GetNbinsX(); int nbinsy = yhist->GetNbinsX(); std::string name = std::string(xhist->GetName()) + "_vs_" + std::string(yhist->GetName()); std::string titles = ";" + std::string(xhist->GetXaxis()->GetTitle()) + ";" + std::string(yhist->GetXaxis()->GetTitle()) + ";" + zname; TH2D *newplot = new TH2D(name.c_str(), (name + titles).c_str(), nbinsx, &xedges[0], nbinsy, &yedges[0]); return newplot; } //*************************************************** void PlotUtils::MatchEmptyBins(TH1D *data, TH1D *mc) { //************************************************** for (int i = 0; i < data->GetNbinsX(); i++) { if (data->GetBinContent(i + 1) == 0.0 or data->GetBinError(i + 1) == 0.0) mc->SetBinContent(i + 1, 0.0); } return; } //*************************************************** void PlotUtils::MatchEmptyBins(TH2D *data, TH2D *mc) { //************************************************** for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsY(); j++) { if (data->GetBinContent(i + 1, j + 1) == 0.0 or data->GetBinError(i + 1, j + 1) == 0.0) mc->SetBinContent(i + 1, j + 1, 0.0); } } return; } //*************************************************** TH1D *PlotUtils::GetProjectionX(TH2D *hist, TH2I *mask) { //*************************************************** TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask); // This includes the underflow/overflow TH1D *hist_X = maskedhist->ProjectionX("_px", 1, maskedhist->GetXaxis()->GetNbins()); hist_X->SetTitle(Form("%s x no under/overflow", hist_X->GetTitle())); delete maskedhist; return hist_X; } //*************************************************** TH1D *PlotUtils::GetProjectionY(TH2D *hist, TH2I *mask) { //*************************************************** TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask); // This includes the underflow/overflow TH1D *hist_Y = maskedhist->ProjectionY("_py", 1, maskedhist->GetYaxis()->GetNbins()); hist_Y->SetTitle(Form("%s y no under/overflow", hist_Y->GetTitle())); delete maskedhist; return hist_Y; }