diff --git a/README.md b/README.md index 76556166..a1905a12 100644 --- a/README.md +++ b/README.md @@ -244,7 +244,7 @@ You can plot kinematic distribution of your sample using ```bash ./bin/KinemDistributionTutorial TutorialConfigs/FitterConfig.yaml ``` -Notice same config being used. In other words you can add or disable sample, modify cuts in same way as was discussed. +Notice same config being used. In other words you can add or disable sample, modify cuts in same way as was discussed. At the bottom of the config, you can specify any individual variables you would like to plot with this executable, along with any selection cuts for each plot. Example of plot you can see here: Kinematic example diff --git a/SamplesTutorial/SampleHandlerTutorial.cpp b/SamplesTutorial/SampleHandlerTutorial.cpp index 5201a717..03e624ae 100755 --- a/SamplesTutorial/SampleHandlerTutorial.cpp +++ b/SamplesTutorial/SampleHandlerTutorial.cpp @@ -8,6 +8,11 @@ SampleHandlerTutorial::SampleHandlerTutorial(const std::string& config_name, Par KinematicParameters = &KinematicParametersTutorial; ReversedKinematicParameters = &ReversedKinematicParametersTutorial; + // === JM assign kinematic vector maps === + KinematicVectors = &KinematicVectorsTutorial; + ReversedKinematicVectors = &ReversedKinematicVectorsTutorial; + // ======================================= + isATM = false; Initialise(); } @@ -103,6 +108,10 @@ void SampleHandlerTutorial::SetupWeightPointers() { } } +void SampleHandlerTutorial::CleanMemoryBeforeFit() { + CleanVector(TutorialPlottingSamples); +} + // ************************************************ int SampleHandlerTutorial::SetupExperimentMC() { // ************************************************ @@ -115,6 +124,7 @@ int SampleHandlerTutorial::SetupExperimentMC() { delete _Chain; TutorialSamples.resize(nEntries); + TutorialPlottingSamples.resize(nEntries); int TotalEventCounter = 0.; for(int iChannel = 0 ; iChannel < static_cast(OscChannels.size()); iChannel++) { @@ -182,8 +192,25 @@ int SampleHandlerTutorial::SetupExperimentMC() { */ _data->GetEntry(0); + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<> unif(0,3); + std::normal_distribution<> mu_angle(0,M_PI/8); + std::normal_distribution<> pi_angle(0,M_PI/2); + std::uniform_real_distribution<> nucl_angle(-M_PI,M_PI); + std::exponential_distribution<> pi_energy(1./0.5); + std::exponential_distribution<> nucl_energy(1./2); + for (int i = 0; i < nEntries; ++i) { // Loop through tree _data->GetEntry(i); + + // === JM resize particle-level vectors === + // JM: We don't have particle-level info in the tutorial sample, so will fake it + int nParticles = 5; //fake number of particles in event + TutorialPlottingSamples[TotalEventCounter].particle_energy.resize(nParticles); + TutorialPlottingSamples[TotalEventCounter].particle_pdg.resize(nParticles); + TutorialPlottingSamples[TotalEventCounter].particle_beamangle.resize(nParticles); + // ======================================== TutorialSamples[TotalEventCounter].TrueEnu = Enu_true; // HH: We don't have Erec in the tutorial sample, so we set it to the true energy @@ -206,6 +233,51 @@ int SampleHandlerTutorial::SetupExperimentMC() { if(isATM) { TutorialSamples[TotalEventCounter].TrueCosZenith = trueCZ; } + + // === JM loop through particles in event === + for (int iParticle = 0; iParticle < nParticles; ++iParticle) { + //JM: No particle-level data in sample, so fake it + if (iParticle==0) { + TutorialPlottingSamples[TotalEventCounter].particle_pdg[iParticle] = PDGLep; + TutorialPlottingSamples[TotalEventCounter].particle_energy[iParticle] = ELep; + TutorialPlottingSamples[TotalEventCounter].particle_beamangle[iParticle] = mu_angle(gen); + } + else { + int particle_seed = unif(gen); + double angle, energy; + int pdg; + switch (particle_seed) { + case 0: + pdg = 211; + energy = pi_energy(gen); + angle = pi_angle(gen); + while (angle>M_PI || angle<-M_PI) angle = pi_angle(gen); + break; + case 1: + pdg = -211; + energy = pi_energy(gen); + angle = pi_angle(gen); + while (angle>M_PI || angle<-M_PI) angle = pi_angle(gen); + break; + case 2: + pdg = 2212; + energy = nucl_energy(gen); + angle = nucl_angle(gen); + break; + case 3: + pdg = 2112; + energy = nucl_energy(gen); + angle = nucl_angle(gen); + break; + default: + break; + } + TutorialPlottingSamples[TotalEventCounter].particle_energy[iParticle] = energy; + TutorialPlottingSamples[TotalEventCounter].particle_beamangle[iParticle] = angle; + TutorialPlottingSamples[TotalEventCounter].particle_pdg[iParticle] = pdg; + } + } + // ========================================== TotalEventCounter++; } _sampleFile->Close(); @@ -230,6 +302,32 @@ double SampleHandlerTutorial::ReturnKinematicParameter(std::string KinematicPara return ReturnKinematicParameter(KinPar, iEvent); } +// === JM Define ReturnKinematicVector functions === +std::vector SampleHandlerTutorial::ReturnKinematicVector(KinematicParticleVecs KinVec, int iEvent) { + switch (KinVec) { + case kParticleEnergy: + return TutorialPlottingSamples[iEvent].particle_energy; + case kParticlePDG: + return TutorialPlottingSamples[iEvent].particle_pdg; + case kParticleBeamAngle: + return TutorialPlottingSamples[iEvent].particle_beamangle; + default: + MACH3LOG_ERROR("Unrecognized Kinematic Vector: {}", static_cast(KinVec)); + throw MaCh3Exception(__FILE__, __LINE__); + } +} + +std::vector SampleHandlerTutorial::ReturnKinematicVector(int KinematicVector, int iEvent) { + KinematicParticleVecs KinVec = static_cast(std::round(KinematicVector)); + return ReturnKinematicVector(KinVec, iEvent); +} + +std::vector SampleHandlerTutorial::ReturnKinematicVector(std::string KinematicVector, int iEvent) { + KinematicParticleVecs KinVec = static_cast(ReturnKinematicVectorFromString(KinematicVector)); + return ReturnKinematicVector(KinVec, iEvent); +} +// ================================================= + const double* SampleHandlerTutorial::GetPointerToKinematicParameter(KinematicTypes KinPar, int iEvent) { switch (KinPar) { case kTrueNeutrinoEnergy: @@ -270,6 +368,8 @@ void SampleHandlerTutorial::SetupFDMC() { } std::vector SampleHandlerTutorial::ReturnKinematicParameterBinning(std::string KinematicParameterStr) { + if (IsSubEventVarString(KinematicParameterStr)) return ReturnKinematicVectorBinning(KinematicParameterStr); + std::vector binningVector; KinematicTypes KinematicParameter = static_cast(ReturnKinematicParameterFromString(KinematicParameterStr)); @@ -300,3 +400,40 @@ std::vector SampleHandlerTutorial::ReturnKinematicParameterBinning(std:: return binningVector; } +std::vector SampleHandlerTutorial::ReturnKinematicVectorBinning(std::string KinematicVectorStr) { + std::vector binningVector; + KinematicParticleVecs KinematicVector = static_cast(ReturnKinematicVectorFromString(KinematicVectorStr)); + + int nBins = 0; + double bin_start = 0; + double bin_end = 0; + + switch(KinematicVector){ + case(kParticleEnergy): + nBins = 20; + bin_start = 0; + bin_end = 6; + break; + case(kParticlePDG): + nBins = 100; + bin_start = -2000; + bin_end = 2000; + break; + case(kParticleBeamAngle): + nBins = 20; + bin_start = -M_PI; + bin_end = M_PI; + break; + default: + MACH3LOG_ERROR("Binning for {} not found", KinematicVectorStr); + throw MaCh3Exception(__FILE__,__LINE__); + break; + } + double step = (bin_end - bin_start)/nBins; + + for(double bin_i = bin_start; bin_i <= bin_end; bin_i+=step){ + binningVector.push_back(bin_i); + } + + return binningVector; +} diff --git a/SamplesTutorial/SampleHandlerTutorial.h b/SamplesTutorial/SampleHandlerTutorial.h index 7bd0271e..a4a19540 100755 --- a/SamplesTutorial/SampleHandlerTutorial.h +++ b/SamplesTutorial/SampleHandlerTutorial.h @@ -3,6 +3,7 @@ #include "Samples/SampleHandlerFD.h" #include "StructsTutorial.h" #include "SplinesTutorial/BinnedSplinesTutorial.h" +#include class SampleHandlerTutorial : public SampleHandlerFD { @@ -10,7 +11,12 @@ class SampleHandlerTutorial : public SampleHandlerFD SampleHandlerTutorial(const std::string& config_name, ParameterHandlerGeneric* parameter_handler, const std::shared_ptr& Oscillator_ = nullptr); virtual ~SampleHandlerTutorial(); + enum KinematicTypes {kTrueNeutrinoEnergy, kTrueQ2, kM3Mode, kRecoNeutrinoEnergy}; + + // === JM enum for particle-level parameters === + enum KinematicParticleVecs {kParticleEnergy, kParticlePDG, kParticleBeamAngle}; + // ============================================= std::vector ReturnKinematicParameterBinning(std::string KinematicParameter) override; @@ -24,11 +30,18 @@ class SampleHandlerTutorial : public SampleHandlerFD int SetupExperimentMC() override; - void CleanMemoryBeforeFit() override {}; + void CleanMemoryBeforeFit() override; double ReturnKinematicParameter(KinematicTypes KinPar, int iEvent); double ReturnKinematicParameter(int KinematicVariable, int iEvent); double ReturnKinematicParameter(std::string KinematicParameter, int iEvent); + + // === JM ReturnKinematicVector declarations for particle-level parameters === + std::vector ReturnKinematicVector(KinematicParticleVecs KinVec, int iEvent); + std::vector ReturnKinematicVector(int KinematicVector, int iEvent); + std::vector ReturnKinematicVector(std::string KinematicVector, int iEvent); + std::vector ReturnKinematicVectorBinning(std::string KinematicParameter); + // =========================================================================== const double* GetPointerToKinematicParameter(KinematicTypes KinPar, int iEvent); const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iEvent) override; @@ -38,6 +51,7 @@ class SampleHandlerTutorial : public SampleHandlerFD void CalcWeightFunc(int iEvent) override {return; (void)iEvent;} std::vector TutorialSamples; + std::vector TutorialPlottingSamples; const std::unordered_map KinematicParametersTutorial = { {"TrueNeutrinoEnergy", kTrueNeutrinoEnergy}, @@ -52,6 +66,20 @@ class SampleHandlerTutorial : public SampleHandlerFD {kM3Mode,"Mode"}, {kRecoNeutrinoEnergy, "RecoNeutrinoEnergy"}, }; + + // === JM maps for particle-level parameters === + const std::unordered_map KinematicVectorsTutorial = { + {"ParticleEnergy", kParticleEnergy}, + {"ParticlePDG", kParticlePDG}, + {"ParticleBeamAngle", kParticleBeamAngle}, + }; + + const std::unordered_map ReversedKinematicVectorsTutorial = { + {kParticleEnergy, "ParticleEnergy"}, + {kParticlePDG, "ParticlePDG"}, + {kParticleBeamAngle, "kParticleBeamAngle"}, + }; + // ============================================= double pot; bool isATM; diff --git a/SamplesTutorial/StructsTutorial.h b/SamplesTutorial/StructsTutorial.h index a4e217d5..0f9fe41e 100755 --- a/SamplesTutorial/StructsTutorial.h +++ b/SamplesTutorial/StructsTutorial.h @@ -31,3 +31,14 @@ struct TutorialMCInfo { /// Lepton energy double ELep = M3::_BAD_DOUBLE_; }; + +struct TutorialMCPlottingInfo { + // === JM Particle-level kinematic parameters === + /// particle energy + std::vector particle_energy = {}; + /// particle pdg + std::vector particle_pdg = {}; + /// particle angle to beam + std::vector particle_beamangle = {}; + // ============================================== +}; diff --git a/Tutorial/KinemDistributionTutorial.cpp b/Tutorial/KinemDistributionTutorial.cpp index f190647b..dc52d966 100755 --- a/Tutorial/KinemDistributionTutorial.cpp +++ b/Tutorial/KinemDistributionTutorial.cpp @@ -2,7 +2,21 @@ #include "Fitters/MaCh3Factory.h" #include "SamplesTutorial/SampleHandlerTutorial.h" +struct PlotKinematicCut { + std::string ParamToCutOn = ""; + double LowerBound = M3::_BAD_DOUBLE_; + double UpperBound = M3::_BAD_DOUBLE_; +}; + +struct Plot { + std::string Name; + std::vector VarStrings; + std::vector> BinEdges; + std::vector SelectionCuts = {}; +}; + int main(int argc, char **argv) { + //JM OscillationChannel kinematic variable does not exist, so only supports looping modes int Selection = 0; //0 to loop modes, 1 to loop osc channels // Initialise manger responsible for config handling @@ -23,6 +37,59 @@ int main(int argc, char **argv) { TString OutputName = TString(OutName.str()) + "_KinemPlot" + ".pdf"; std::vector vecParams = {"TrueNeutrinoEnergy", "TrueQ2", "RecoNeutrinoEnergy"}; + + //JM Read in kinematic distribution plots from config + std::vector PlotsToDraw = {}; + auto ConfigPlots = FitManager->raw()["KinematicDistributionPlots"]; + + for (const auto& ConfigPlot : ConfigPlots) { + Plot PlotToDraw; + PlotToDraw.Name = Get(ConfigPlot["Name"], __FILE__, __LINE__); + PlotToDraw.VarStrings = Get>(ConfigPlot["VarStrings"], __FILE__, __LINE__); + if (PlotToDraw.VarStrings.size() != 1 && PlotToDraw.VarStrings.size() != 2) { + MACH3LOG_ERROR("Error in yaml file: In KinemDistribtuion Plot {}, VarStrings is of size {}. VarString should be of size 1 or 2 (higher dimensional histogram plotting is not yet supported)"); + throw MaCh3Exception(__FILE__,__LINE__); + } + PlotToDraw.BinEdges = Get>>(ConfigPlot["VarBins"], __FILE__,__LINE__); + if (PlotToDraw.BinEdges.size() != 1 && PlotToDraw.BinEdges.size() != 2) { + MACH3LOG_ERROR("Error in yaml file: In KinemDistribtuion Plot {}, BinEdges is of size {}. VarString should be of size 1 or 2 (higher dimensional histogram plotting is not yet supported)"); + throw MaCh3Exception(__FILE__,__LINE__); + } + + //If binning vector is of size 3, treat as [nbins, xmin, xmax] (otherwise treat as bin edges) + for (unsigned int iBinning=0; iBinning(); + std::vector range = Cut["Range"].as>(); + + if (range.size() != 2) { + MACH3LOG_ERROR("Error in yaml file: In KinemDistribution Plot {}, KinematicCut {} has range of size {}. Range should be of size 2.", PlotToDraw.Name, SelectionCut.ParamToCutOn, range.size()); + throw MaCh3Exception(__FILE__,__LINE__); + } + SelectionCut.LowerBound = range[0]; + SelectionCut.UpperBound = range[1]; + PlotToDraw.SelectionCuts.push_back(SelectionCut); + } + PlotsToDraw.push_back(PlotToDraw); + } TCanvas* Canv = new TCanvas("Canv",""); Canv->Divide(1,2); @@ -33,7 +100,7 @@ int main(int argc, char **argv) { for (size_t iPDF=0; iPDF < mySamples.size(); iPDF++) { MACH3LOG_INFO("Number of samples: {}", mySamples[iPDF]->GetNsamples()); - + THStack* Stack = new THStack(*mySamples[iPDF]->ReturnStackedHistBySelection1D(vecParams[iParam], Selection)); TLegend* Legend = new TLegend(*mySamples[iPDF]->ReturnStackHistLegend()); @@ -54,13 +121,47 @@ int main(int argc, char **argv) { delete Canv; Canv = new TCanvas("Canv",""); - if(vecParams.size() == 2) - { + TH1* Hist; + int WeightStyle = 1; + for (size_t iHist=0; iHist PlotVar_Str = PlotsToDraw[iHist].VarStrings; + int histdim = PlotVar_Str.size(); + TAxis AxisX = TAxis(PlotsToDraw[iHist].BinEdges[0].size()-1,PlotsToDraw[iHist].BinEdges[0].data()); + TAxis AxisY; + if (histdim == 2) AxisY = TAxis(PlotsToDraw[iHist].BinEdges[1].size()-1,PlotsToDraw[iHist].BinEdges[1].data()); + for (size_t iPDF = 0;iPDF < mySamples.size(); iPDF++) { - TH2D* Hist = static_cast(mySamples[iPDF]->Get2DVarHist(vecParams[0], vecParams[1])); + std::vector EventSelectionVector = {}; + std::vector SubEventSelectionVector = {}; + for (size_t iCut=0; iCutIsSubEventVarString(PlotsToDraw[iHist].SelectionCuts[iCut].ParamToCutOn)) { + Selection.ParamToCutOnIt = mySamples[iPDF]->ReturnKinematicVectorFromString(PlotsToDraw[iHist].SelectionCuts[iCut].ParamToCutOn); + SubEventSelectionVector.push_back(Selection); + } + else { + Selection.ParamToCutOnIt = mySamples[iPDF]->ReturnKinematicParameterFromString(PlotsToDraw[iHist].SelectionCuts[iCut].ParamToCutOn); + EventSelectionVector.push_back(Selection); + } + } + + if (histdim == 1) { + Hist = (TH1*)mySamples[iPDF]->Get1DVarHist(PlotVar_Str[0], EventSelectionVector, WeightStyle, &AxisX, SubEventSelectionVector); + Hist->GetYaxis()->SetTitle("Events"); + } + else { + Hist = (TH1*)mySamples[iPDF]->Get2DVarHist(PlotVar_Str[0], PlotVar_Str[1], EventSelectionVector, WeightStyle, &AxisX, &AxisY, SubEventSelectionVector); + Hist->GetYaxis()->SetTitle(PlotVar_Str[1].c_str()); + } Canv->cd(1); - Hist->SetTitle(mySamples[iPDF]->GetTitle().c_str()); + Hist->SetTitle(PlotsToDraw[iHist].Name.c_str()); + Hist->GetXaxis()->SetTitle(PlotVar_Str[0].c_str()); Hist->SetStats(false); Hist->Draw("COLZ"); Canv->Print(OutputName); @@ -69,6 +170,5 @@ int main(int argc, char **argv) { } Canv->Print(OutputName+"]"); delete Canv; - return 0; } diff --git a/TutorialConfigs/FitterConfig.yaml b/TutorialConfigs/FitterConfig.yaml index b88a6561..95ab1996 100755 --- a/TutorialConfigs/FitterConfig.yaml +++ b/TutorialConfigs/FitterConfig.yaml @@ -69,3 +69,14 @@ LLHScan: ScanRanges: #delta_cp: [-3.14,3.14] +#KinemDistribution extra plots to draw +KinematicDistributionPlots: +- Name: "True Neutrino Energy (0.5 GeV < Enu < 3 GeV)" + VarStrings: ["TrueNeutrinoEnergy"] + VarBins: [[40,0,4]] + KinematicCuts: + - VarString: "TrueNeutrinoEnergy" + Range: [0.5,3] +- Name: "Q2 vs True Enu" + VarStrings: ["TrueNeutrinoEnergy", "TrueQ2"] + VarBins: [[50,0,4],[50,0,6]]