diff --git a/mcmc/PSO.cpp b/mcmc/PSO.cpp index afeaaf521..96e42b42f 100644 --- a/mcmc/PSO.cpp +++ b/mcmc/PSO.cpp @@ -2,8 +2,9 @@ #include +// *************** PSO::PSO(manager *man) : LikelihoodFit(man) { - +// *************** fConstriction = fitMan->raw()["General"]["PSO"]["Constriction"].as(); fInertia = fitMan->raw()["General"]["PSO"]["Inertia"].as()*fConstriction; fOne = fitMan->raw()["General"]["PSO"]["One"].as()*fConstriction; @@ -20,8 +21,9 @@ PSO::PSO(manager *man) : LikelihoodFit(man) { } } +// *************** void PSO::runMCMC(){ - +// *************** PrepareFit(); if(fTestLikelihood){ @@ -41,17 +43,16 @@ void PSO::runMCMC(){ return; } - +// ************************* void PSO::init(){ - +// ************************* fBestValue = 1234567890.0; //KS: For none PCA this will be eqaul to normal parameters //const int NparsPSOFull = NPars; //const int NparsPSO = NParsPCA; - std::cout << "Preparing PSO" << std::endl; - + MACH3LOG_INFO("Preparing PSO"); // Initialise bounds on parameters if(fTestLikelihood){ for (int i = 0; i < fDim; i++){ @@ -105,9 +106,9 @@ void PSO::init(){ } } - std::cout << "Printing Minimums and Maximums of Variables to be minimised" << std::endl; - for (int i =0; i > PSO::bisection(std::vectorposition,double minimum, double range, double precision){ +// ************************* std::vector> uncertainties_list; - for (unsigned int i = 0; i< position.size(); ++i){ - std::cout << i << std::endl; + for (unsigned int i = 0; i< position.size(); ++i) { + MACH3LOG_INFO("{}", i); //std::vector uncertainties; std::vector new_position = position; new_position[i] = position[i]-range; double val_1 = CalcChi(new_position)-minimum-1.0; @@ -201,7 +204,8 @@ std::vector > PSO::bisection(std::vectorposition,dou position_list_p[1] = new_bisect_position_p; res = std::abs(position[2]-position[0]); res_p = std::abs(position_list_p[1][i]-position_list_p[0][i]); - //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl; + MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]); + } else{ std::vector new_bisect_position_p = position_list_p[1];new_bisect_position_p[i] += (position_list_p[2][i]-position_list_p[1][i])/2; @@ -211,17 +215,20 @@ std::vector > PSO::bisection(std::vectorposition,dou value_list_p[1] = new_val_p; position_list_p[1] = new_bisect_position_p; res_p = std::abs(position_list_p[2][i]-position_list_p[1][i]); - //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl; + MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]); } } uncertainties_list.push_back({std::abs(position[i]-position_list[1][i]),std::abs(position[i]-position_list_p[1][i])}); - std::cout << "Uncertainty finished for d = "<< i << std::endl; - std::cout << std::setprecision(10)<< "LLR values for ± positive and negative uncertainties are " << CalcChi(position_list[1]) << " and " << CalcChi(position_list_p[1]) << std::endl; + MACH3LOG_INFO("Uncertainty finished for d = {}", i); + MACH3LOG_INFO("LLR values for ± positive and negative uncertainties are {:<10.2f} and {:<10.2f}", + CalcChi(position_list[1]), CalcChi(position_list_p[1])); } return uncertainties_list; } -std::vector> PSO::calc_uncertainty(std::vectorposition,double minimum) { +// ************************* +std::vector> PSO::calc_uncertainty(std::vectorposition, double minimum) { +// ************************* std::vector pos_uncertainty(position.size()); std::vector neg_uncertainty(position.size()); int num = 200; @@ -255,7 +262,7 @@ std::vector> PSO::calc_uncertainty(std::vectorpositi } } neg_uncertainty[i] = x[closest_index]; - std::cout << "Neg" << std::endl; + MACH3LOG_INFO("Neg"); x.assign(num, 0); y.assign(num, 0); StepPoint = (pos_stop-start) / (num - 1); @@ -283,11 +290,13 @@ std::vector> PSO::calc_uncertainty(std::vectorpositi return res; } +// ************************* void PSO::uncertainty_check(std::vector previous_pos){ +// ************************* std::vector> x_list; std::vector> y_list; std::vector position = previous_pos; - int num = 5000; + constexpr int num = 5000; for (unsigned int i = 0;i previous_pos){ std::vector y(num); double StepPoint = (stop - start) / (num - 1); double value = start; - // std::cout << "result for fDim " << 1 << std::endl; - for (int j =0;j< num; ++j){ + MACH3LOG_TRACE("result for fDim: {}", 1); + for (int j = 0;j < num; ++j) { position[i] = value; double LLR = CalcChi(position); x[j] = value; @@ -305,21 +314,24 @@ void PSO::uncertainty_check(std::vector previous_pos){ value += StepPoint; } position[i] = curr_ival; - std::cout << " " << std::endl; - std::cout << "For fDim" << i+1 << " x list is " ; - for (unsigned int k= 0;k total_pos(fDim,0.0); for (int i = 0; i < fParticles; ++i) { @@ -377,8 +389,9 @@ double PSO::swarmIterate(){ return mean_dist_sq; } +// ************************* void PSO::run() { - +// ************************* double mean_dist_sq = 0; int iter = 0; @@ -398,11 +411,11 @@ void PSO::run() { accCount++; if (i%100 == 0){ - std::cout << "Mean Dist Sq = " << mean_dist_sq <get_personal_best_position()[j]); } } if(fConvergence > 0.0){ @@ -411,18 +424,20 @@ void PSO::run() { } } } - std::cout << "Finished after " << iter <<" runs out of "<< fIterations << std::endl; - std::cout << "Mean Dist " << mean_dist_sq <get_personal_best_value() << std::endl; - + MACH3LOG_INFO("Finished after {} runs out of {}", iter, fIterations); + MACH3LOG_INFO("Mean Dist: {:.2f}", mean_dist_sq); + MACH3LOG_INFO("Best LLR: {:.2f}", get_best_particle()->get_personal_best_value()); uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005); - std::cout << "Position for Global Minimum = "<get_personal_best_position()[i], uncertainties[i][1], uncertainties[i][0]); } } +// ************************* void PSO::WriteOutput(){ +// ************************* + outputFile->cd(); TVectorD* PSOParValue = new TVectorD(fDim); diff --git a/mcmc/PSO.h b/mcmc/PSO.h index f5b3480d5..9baf0758e 100644 --- a/mcmc/PSO.h +++ b/mcmc/PSO.h @@ -1,9 +1,5 @@ #pragma once -// -// Created by Emily Ip on 24/2/2023. -// -// Created by Emily Ip on 26/1/2023. -// + // C++ includes #include #include @@ -15,6 +11,7 @@ /// @brief Class particle - stores the position, velocity and personal best /// With functions which move particle and update velocity +/// @note Created by Emily Ip on 24/2/2023. class particle{ public: particle(){}; @@ -72,6 +69,7 @@ class particle{ /// @brief Class PSO, consist of a vector of object Class Particle and global best /// Takes in the size (number of particle) and number of iteration /// functions includes: finding global best, updating velocity, actual minimisation function +/// @note Created by Emily Ip on 24/2/2023. class PSO : public LikelihoodFit { public: /// @brief constructor diff --git a/mcmc/SampleSummary.cpp b/mcmc/SampleSummary.cpp index 643dce61a..9988383f1 100644 --- a/mcmc/SampleSummary.cpp +++ b/mcmc/SampleSummary.cpp @@ -42,9 +42,7 @@ SampleSummary::SampleSummary(const int n_Samples, const std::string &Filename, s first_pass = true; Outputfile = nullptr; OutputTree = nullptr; - Dir = nullptr; - - rnd = new TRandom3(); + rnd = std::make_unique(); DataHist = new TH2Poly*[nSamples]; DataHist_ProjectX = new TH1D*[nSamples]; @@ -60,7 +58,7 @@ SampleSummary::SampleSummary(const int n_Samples, const std::string &Filename, s if(DoBetaParam) BetaHist = new TH1D**[nSamples]; else BetaHist = nullptr; - maxBins = new int[nSamples]; + maxBins.resize(nSamples); lnLHist_Mean = new TH2Poly*[nSamples]; lnLHist_Mode = new TH2Poly*[nSamples]; @@ -152,21 +150,6 @@ SampleSummary::SampleSummary(const int n_Samples, const std::string &Filename, s lnLHist_Sample_DrawflucDraw[i] = NULL; lnLHist_Sample_PredflucDraw[i] = NULL; }//end loop over samples - - llh_data_draw = NULL; - llh_drawfluc_draw = NULL; - llh_predfluc_draw = NULL; - llh_rate_data_draw = NULL; - llh_rate_predfluc_draw = NULL; - - llh_data_drawfluc = NULL; - llh_data_predfluc = NULL; - llh_draw_pred = NULL; - llh_drawfluc_pred = NULL; - - llh_predfluc_pred = NULL; - llh_drawfluc_predfluc = NULL; - llh_datafluc_draw = NULL; DoByModePlots = false; MeanHist_ByMode = NULL; @@ -186,8 +169,6 @@ SampleSummary::~SampleSummary() { //ROOT is weird and once you write TFile claim ownership of histograms. Best is to first delete histograms and then close file Outputfile->Close(); delete Outputfile; - - if(rnd != nullptr) delete rnd; if(lnLHist != NULL) delete lnLHist; if(lnLHist_drawdata != NULL) delete lnLHist_drawdata; @@ -271,28 +252,10 @@ SampleSummary::~SampleSummary() { delete[] lnLHist_Sample_DrawData; delete[] lnLHist_Sample_DrawflucDraw; delete[] lnLHist_Sample_PredflucDraw; - - delete[] maxBins; - + delete[] PosteriorHist; delete[] w2Hist; if(DoBetaParam) delete[] BetaHist; - - delete[] llh_data_draw; - delete[] llh_data_drawfluc; - delete[] llh_data_predfluc; - delete[] llh_rate_data_draw; - delete[] llh_rate_predfluc_draw; - delete[] llh_draw_pred; - delete[] llh_drawfluc_pred; - delete[] llh_drawfluc_predfluc; - delete[] llh_drawfluc_draw; - delete[] llh_predfluc_pred; - delete[] llh_predfluc_draw; - delete[] llh_datafluc_draw; - - delete[] llh_data_draw_ProjectX; - delete[] llh_drawfluc_draw_ProjectX; } // ******************* @@ -375,9 +338,9 @@ void SampleSummary::AddNominal(std::vector &Nominal, std::vector(NomW2[i]->Clone()); lnLHist_Mean[i] = static_cast(NominalHist[i]->Clone()); - lnLHist_Mean[i]->SetDirectory(0); + lnLHist_Mean[i]->SetDirectory(nullptr); lnLHist_Mode[i] = static_cast(NominalHist[i]->Clone()); - lnLHist_Mode[i]->SetDirectory(0); + lnLHist_Mode[i]->SetDirectory(nullptr); lnLHist_Mean_ProjectX[i] = static_cast(DataHist_ProjectX[i]->Clone()); MeanHist[i] = static_cast(NominalHist[i]->Clone()); if(DoBetaParam) MeanHistCorrected[i] = static_cast(NominalHist[i]->Clone()); @@ -440,12 +403,12 @@ void SampleSummary::AddNominal(std::vector &Nominal, std::vectorGetYaxis()->SetTitle("Events"); ViolinHists_ProjectX[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetXaxis()->GetTitle()).c_str() ); - ViolinHists_ProjectX[i]->SetDirectory(0); + ViolinHists_ProjectX[i]->SetDirectory(nullptr); ViolinHists_ProjectY[i] = new TH2D((name+"_Violin_ProjectY").c_str(), (name+"_Violin_ProjectY").c_str(), int(ybins.size()-1), &ybins[0] , 400, 0, MaxBinning); ViolinHists_ProjectY[i]->GetYaxis()->SetTitle("Events"); ViolinHists_ProjectY[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetYaxis()->GetTitle()).c_str()); - ViolinHists_ProjectY[i]->SetDirectory(0); + ViolinHists_ProjectY[i]->SetDirectory(nullptr); ModeHist[i]->SetNameTitle((name+"_mode").c_str(), (name+"_mode").c_str()); ModeHist[i]->Reset(""); @@ -529,9 +492,9 @@ void SampleSummary::AddThrow(std::vector &SampleVector, std::vectorGetYMin() << "-" << bin->GetYMax() << ")"; PosteriorHist[SampleNum][i] = new TH1D(ss2.str().c_str(), ss2.str().c_str(),nXBins, 1, -1); - PosteriorHist[SampleNum][i]->SetDirectory(0); + PosteriorHist[SampleNum][i]->SetDirectory(nullptr); w2Hist[SampleNum][i] = new TH1D(("w2_"+ss2.str()).c_str(), ("w2_"+ss2.str()).c_str(),nXBins, 1, -1); - w2Hist[SampleNum][i]->SetDirectory(0); + w2Hist[SampleNum][i]->SetDirectory(nullptr); if(DoBetaParam) { std::string betaName = "#beta_param_"; @@ -592,7 +555,7 @@ void SampleSummary::AddThrowByMode(std::vector> &SampleVec for (int j = 0; j < Modes->GetNModes()+1; j++) { PosteriorHist_ByMode[SampleNum][j] = new TH1D*[maxBins[SampleNum]+1]; - const int nXBins = 500; + constexpr int nXBins = 500; std::string name = std::string(NominalHist[SampleNum]->GetName()); name = name.substr(0, name.find("_nom")); @@ -653,36 +616,36 @@ void SampleSummary::PrepareOutput() { // The array of doubles we write to the TTree // Data vs Draw - llh_data_draw = new double[nSamples]; + llh_data_draw.resize(nSamples); // Fluctuated Draw vs Draw - llh_drawfluc_draw = new double[nSamples]; + llh_drawfluc_draw.resize(nSamples); // Fluctuated Predicitve vs Draw - llh_predfluc_draw = new double[nSamples]; + llh_predfluc_draw.resize(nSamples); // Data vs Draw using Rate - llh_rate_data_draw = new double[nSamples]; + llh_rate_data_draw.resize(nSamples); // Data vs Fluctuated Predictive using Rate - llh_rate_predfluc_draw = new double[nSamples]; + llh_rate_predfluc_draw.resize(nSamples); // Data vs Fluctuated Draw - llh_data_drawfluc = new double[nSamples]; + llh_data_drawfluc.resize(nSamples); // Data vs Fluctuated Predictive - llh_data_predfluc = new double[nSamples]; + llh_data_predfluc.resize(nSamples); // Draw vs Predictive - llh_draw_pred = new double[nSamples]; + llh_draw_pred.resize(nSamples); // Fluctuated Draw vs Predictive - llh_drawfluc_pred = new double[nSamples]; + llh_drawfluc_pred.resize(nSamples); // Fluctuated Draw vs Fluctuated Predictive - llh_drawfluc_predfluc = new double[nSamples]; + llh_drawfluc_predfluc.resize(nSamples); // Fluctuated Predictive vs Predictive - llh_predfluc_pred = new double[nSamples]; + llh_predfluc_pred.resize(nSamples); // Fluctuated Data vs Draw - llh_datafluc_draw = new double[nSamples]; + llh_datafluc_draw.resize(nSamples); // Data vs Draw for 1D projection - llh_data_draw_ProjectX = new double[nSamples]; - llh_drawfluc_draw_ProjectX = new double[nSamples]; + llh_data_draw_ProjectX.resize(nSamples); + llh_drawfluc_draw_ProjectX.resize(nSamples); // The output tree we're going to write to OutputTree = new TTree("LLH_draws", "LLH_draws"); @@ -746,7 +709,7 @@ void SampleSummary::PrepareOutput() { OutputTree->Branch("total_llh_drawfluc_draw_ProjectX", &total_llh_drawfluc_draw_ProjectX); Outputfile->cd(); - Dir = new TDirectory*[nSamples]; + Dir.resize(nSamples); for (int i = 0; i < nSamples; ++i) { // Make a new directory @@ -943,11 +906,13 @@ void SampleSummary::Write() { PosteriorHist[i][b]->Write(); std::string Title = PosteriorHist[i][b]->GetName(); - TLine *TempLine = new TLine(NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(), NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum()); + auto TempLine = std::make_unique(NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(), + NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum()); TempLine->SetLineColor(kRed); TempLine->SetLineWidth(2); - TLine *TempLineData = new TLine(DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(), DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum()); + auto TempLineData = std::make_unique(DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(), + DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum()); TempLineData->SetLineColor(kGreen); TempLineData->SetLineWidth(2); @@ -956,13 +921,13 @@ void SampleSummary::Write() { PosteriorHist[i][b]->Fit(Fitter, "RQ"); Fitter->SetLineColor(kRed-5); - TLegend *Legend = new TLegend(0.4, 0.75, 0.98, 0.90); + auto Legend = std::make_unique(0.4, 0.75, 0.98, 0.90); Legend->SetFillColor(0); Legend->SetFillStyle(0); Legend->SetLineWidth(0); Legend->SetLineColor(0); - Legend->AddEntry(TempLineData, Form("Data #mu=%.2f", DataHist[i]->GetBinContent(b)), "l"); - Legend->AddEntry(TempLine, Form("Prior #mu=%.2f", NominalHist[i]->GetBinContent(b)), "l"); + Legend->AddEntry(TempLineData.get(), Form("Data #mu=%.2f", DataHist[i]->GetBinContent(b)), "l"); + Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.2f", NominalHist[i]->GetBinContent(b)), "l"); Legend->AddEntry(PosteriorHist[i][b], Form("Post, #mu=%.2f#pm%.2f", PosteriorHist[i][b]->GetMean(), PosteriorHist[i][b]->GetRMS()), "l"); Legend->AddEntry(Fitter, Form("Gauss, #mu=%.2f#pm%.2f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l"); std::string TempTitle = std::string(PosteriorHist[i][b]->GetName()); @@ -983,11 +948,8 @@ void SampleSummary::Write() { Legend->Draw("same"); TempCanvas->Write(); - delete TempLine; - delete TempLineData; delete TempCanvas; delete Fitter; - delete Legend; //This isn't useful check only in desperation if(Debug > 1) w2Hist[i][b]->Write(); } @@ -1576,7 +1538,7 @@ void SampleSummary::MakeCutLLH1D(TH1D *Histogram, double llh_ref) { Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str()); // Write a TCanvas and make a line and a filled histogram - TLine *TempLine = new TLine(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum()); + auto TempLine = std::make_unique(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum()); TempLine->SetLineColor(kBlack); TempLine->SetLineWidth(2); @@ -1592,12 +1554,12 @@ void SampleSummary::MakeCutLLH1D(TH1D *Histogram, double llh_ref) { } } - TLegend *Legend = new TLegend(0.6, 0.6, 0.9, 0.9); + auto Legend = std::make_unique(0.6, 0.6, 0.9, 0.9); Legend->SetFillColor(0); Legend->SetFillStyle(0); Legend->SetLineWidth(0); Legend->SetLineColor(0); - Legend->AddEntry(TempLine, Form("Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue), "l"); + Legend->AddEntry(TempLine.get(), Form("Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue), "l"); Legend->AddEntry(Histogram, Form("LLH, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l"); std::string Title = Histogram->GetName(); Title += "_canv"; @@ -1611,17 +1573,14 @@ void SampleSummary::MakeCutLLH1D(TH1D *Histogram, double llh_ref) { TempCanvas->Write(); - delete TempLine; delete TempHistogram; delete TempCanvas; - delete Legend; } // **************** // Make the 2D cut distribution and give the 2D p-value void SampleSummary::MakeCutLLH2D(TH2D *Histogram) { // **************** - const double TotalIntegral = Histogram->Integral(); // Count how many fills are above y=x axis // This is the 2D p-value @@ -1653,7 +1612,7 @@ void SampleSummary::MakeCutLLH2D(TH2D *Histogram) { maximum += Histogram->GetYaxis()->GetBinWidth(Histogram->GetYaxis()->GetNbins()); } else maximum += Histogram->GetXaxis()->GetBinWidth(Histogram->GetXaxis()->GetNbins()); - TLine *TempLine = new TLine(minimum, minimum, maximum, maximum); + auto TempLine = std::make_unique(minimum, minimum, maximum, maximum); TempLine->SetLineColor(kRed); TempLine->SetLineWidth(2); @@ -1669,7 +1628,6 @@ void SampleSummary::MakeCutLLH2D(TH2D *Histogram) { TempLine->Draw("same"); TempCanvas->Write(); - delete TempLine; delete TempCanvas; } @@ -1678,7 +1636,7 @@ void SampleSummary::MakeCutLLH2D(TH2D *Histogram) { void SampleSummary::MakeCutEventRate(TH1D *Histogram, const double DataRate) { // **************** // For the event rate histogram add a TLine to the data rate - TLine *TempLine = new TLine(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum()); + auto TempLine = std::make_unique(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum()); TempLine->SetLineColor(kRed); TempLine->SetLineWidth(2); // Also fit a Gaussian because why not? @@ -1694,12 +1652,12 @@ void SampleSummary::MakeCutEventRate(TH1D *Histogram, const double DataRate) { } } const double pvalue = Above/Histogram->Integral(); - TLegend *Legend = new TLegend(0.4, 0.75, 0.98, 0.90); + auto Legend = std::make_unique(0.4, 0.75, 0.98, 0.90); Legend->SetFillColor(0); Legend->SetFillStyle(0); Legend->SetLineWidth(0); Legend->SetLineColor(0); - Legend->AddEntry(TempLine, Form("Data, %.0f, p-value=%.2f", DataRate, pvalue), "l"); + Legend->AddEntry(TempLine.get(), Form("Data, %.0f, p-value=%.2f", DataRate, pvalue), "l"); Legend->AddEntry(Histogram, Form("MC, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l"); Legend->AddEntry(Fitter, Form("Gauss, #mu=%.1f#pm%.1f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l"); std::string TempTitle = std::string(Histogram->GetName()); @@ -1719,10 +1677,8 @@ void SampleSummary::MakeCutEventRate(TH1D *Histogram, const double DataRate) { TempCanvas->Write(); Histogram->Write(); - delete TempLine; delete TempCanvas; delete Fitter; - delete Legend; } // **************** @@ -1854,7 +1810,7 @@ void SampleSummary::PlotBetaParameters() { BetaHist[i][j]->Fit(Fitter, "RQ"); Fitter->SetLineColor(kRed-5); - TLegend *Legend = new TLegend(0.4, 0.75, 0.98, 0.90); + auto Legend = std::make_unique(0.4, 0.75, 0.98, 0.90); Legend->SetFillColor(0); Legend->SetFillStyle(0); Legend->SetLineWidth(0); @@ -1883,7 +1839,6 @@ void SampleSummary::PlotBetaParameters() { delete TempLine; delete TempCanvas; delete Fitter; - delete Legend; } DirBeta[i]->Write(); delete DirBeta[i]; @@ -1904,7 +1859,7 @@ void SampleSummary::StudyKinematicCorrelations() { timer.Start(); // Data vs Draw for 1D projection - double* NEvents_Sample = new double[nSamples]; + std::vector NEvents_Sample(nSamples); double event_rate = 0.; // The output tree we're going to write to @@ -1930,7 +1885,6 @@ void SampleSummary::StudyKinematicCorrelations() { // Holds the event rate for the distribution TH1D **SumHist = new TH1D*[nSamples]; - for (int i = 0; i < nSamples; ++i) { std::string name = std::string(NominalHist[i]->GetName()); @@ -1985,14 +1939,13 @@ void SampleSummary::StudyKinematicCorrelations() { MakeCutEventRate(SumHist[SampleNum], NoOverflowIntegral(DataHist[SampleNum])); } - // Make a new direcotry + // Make a new directory TDirectory *CorrDir = Outputfile->mkdir("Correlations"); CorrDir->cd(); TMatrixDSym* SampleCorrelation = new TMatrixDSym(nSamples); TH2D*** SamCorr = new TH2D**[nSamples](); - for (int i = 0; i < nSamples; ++i) { SamCorr[i] = new TH2D*[nSamples](); @@ -2008,7 +1961,7 @@ void SampleSummary::StudyKinematicCorrelations() { // TH2D to hold the Correlation SamCorr[i][j] = new TH2D(Form("SamCorr_%i_%i",i,j), Form("SamCorr_%i_%i",i,j), 70, Min_i, Max_i, 70, Min_j, Max_j); - SamCorr[i][j]->SetDirectory(0); + SamCorr[i][j]->SetDirectory(nullptr); SamCorr[i][j]->SetMinimum(0); SamCorr[i][j]->GetXaxis()->SetTitle(SampleNames[i].c_str()); SamCorr[i][j]->GetYaxis()->SetTitle(SampleNames[j].c_str()); @@ -2041,7 +1994,7 @@ void SampleSummary::StudyKinematicCorrelations() { }// End i loop TH2D* hSamCorr = new TH2D("Sample Correlation", "Sample Correlation", nSamples, 0, nSamples, nSamples, 0, nSamples); - hSamCorr->SetDirectory(0); + hSamCorr->SetDirectory(nullptr); hSamCorr->GetZaxis()->SetTitle("Correlation"); hSamCorr->SetMinimum(-1); hSamCorr->SetMaximum(1); @@ -2121,7 +2074,7 @@ void SampleSummary::StudyKinematicCorrelations() { // TH2D to hold the Correlation KinCorr[i][j] = new TH2D(Form("Kin_%i_%i_%i",SampleNum,i,j), Form("Kin_%i_%i_%i",SampleNum,i,j), 70, Min_i, Max_i, 70, Min_j, Max_j); - KinCorr[i][j]->SetDirectory(0); + KinCorr[i][j]->SetDirectory(nullptr); KinCorr[i][j]->SetMinimum(0); KinCorr[i][j]->GetXaxis()->SetTitle(ss2.str().c_str()); @@ -2162,7 +2115,7 @@ void SampleSummary::StudyKinematicCorrelations() { TH2D* hKinCorr = new TH2D(SampleNames[SampleNum].c_str(), SampleNames[SampleNum].c_str(), maxBins[SampleNum], 0, maxBins[SampleNum], maxBins[SampleNum], 0, maxBins[SampleNum]); - hKinCorr->SetDirectory(0); + hKinCorr->SetDirectory(nullptr); hKinCorr->GetZaxis()->SetTitle("Correlation"); hKinCorr->SetMinimum(-1); hKinCorr->SetMaximum(1); @@ -2233,8 +2186,7 @@ void SampleSummary::StudyKinematicCorrelations() { if(DoByModePlots) { // Holds the total event rate by mode - TH1D ** EventHist_ByMode = new TH1D*[Modes->GetNModes()+1]; - + std::vector EventHist_ByMode(Modes->GetNModes()+1); for (int j = 0; j < Modes->GetNModes()+1; j++) { std::string ModeName = Modes->GetMaCh3ModeName(j); @@ -2284,7 +2236,7 @@ void SampleSummary::StudyKinematicCorrelations() { // TH2D to hold the Correlation ModeCorr[i][j] = new TH2D(Form("ModeCorr_%i_%i",i,j), Form("ModeCorr_%i_%i",i,j), 70, Min_i, Max_i, 70, Min_j, Max_j); - ModeCorr[i][j]->SetDirectory(0); + ModeCorr[i][j]->SetDirectory(nullptr); ModeCorr[i][j]->SetMinimum(0); ModeCorr[i][j]->GetXaxis()->SetTitle(Modes->GetMaCh3ModeName(i).c_str()); ModeCorr[i][j]->GetYaxis()->SetTitle(Modes->GetMaCh3ModeName(j).c_str()); @@ -2324,7 +2276,7 @@ void SampleSummary::StudyKinematicCorrelations() { }// End i loop TH2D* hModeCorr = new TH2D("Mode Correlation", "Mode Correlation", Modes->GetNModes()+1, 0, Modes->GetNModes()+1, Modes->GetNModes()+1, 0, Modes->GetNModes()+1); - hModeCorr->SetDirectory(0); + hModeCorr->SetDirectory(nullptr); hModeCorr->GetZaxis()->SetTitle("Correlation"); hModeCorr->SetMinimum(-1); hModeCorr->SetMaximum(1); @@ -2374,7 +2326,6 @@ void SampleSummary::StudyKinematicCorrelations() { { delete EventHist_ByMode[j]; } - delete[] EventHist_ByMode; } for (int i = 0; i < nSamples; ++i) @@ -2665,7 +2616,6 @@ void SampleSummary::StudyDIC() { MACH3LOG_INFO("Effective number of parameters following DIC formalism is equal to: {:.2f}", p_D); MACH3LOG_INFO("DIC test statistic = {:.2f}", DIC_stat); MACH3LOG_INFO("******************************"); - return; } // **************** @@ -2677,5 +2627,4 @@ void SampleSummary::FastViolinFill(TH2D* violin, TH1D* hist_1d){ const double y = violin->GetYaxis()->FindBin(hist_1d->GetBinContent(x+1)); violin->SetBinContent(x+1, y, violin->GetBinContent(x+1, y)+1); } - return; } diff --git a/mcmc/SampleSummary.h b/mcmc/SampleSummary.h index 0d070d064..dfa3e8705 100644 --- a/mcmc/SampleSummary.h +++ b/mcmc/SampleSummary.h @@ -120,7 +120,7 @@ class SampleSummary { inline void NormaliseTH2Poly(TH2Poly* Histogram); /// Random number generator - TRandom3* rnd; + std::unique_ptr rnd; /// KS: Hacky flag to let us know if this is first toy bool first_pass; @@ -236,7 +236,7 @@ class SampleSummary { unsigned int nThrows; /// Max Number of Bins per each sample - int* maxBins; + std::vector maxBins; /// Total LLH for the posterior predictive distribution double llh_total; @@ -246,41 +246,41 @@ class SampleSummary { /// Output filename TFile *Outputfile; /// Directory for each sample - TDirectory **Dir; + std::vector Dir; /// TTree which we save useful information to TTree *OutputTree; /// Data vs Draw - double *llh_data_draw; + std::vector llh_data_draw; /// Fluctuated Draw vs Draw - double *llh_drawfluc_draw; + std::vector llh_drawfluc_draw; /// Fluctuated Predictive vs Draw - double *llh_predfluc_draw; + std::vector llh_predfluc_draw; /// Data vs Draw using rate only - double *llh_rate_data_draw; + std::vector llh_rate_data_draw; /// Fluctuated Predictive vs Draw using rate only - double *llh_rate_predfluc_draw; + std::vector llh_rate_predfluc_draw; /// Data vs Fluctuated Draw - double *llh_data_drawfluc; + std::vector llh_data_drawfluc; /// Data vs Fluctuated Predictive - double *llh_data_predfluc; + std::vector llh_data_predfluc; /// Draw vs Predictive - double *llh_draw_pred; + std::vector llh_draw_pred; /// Fluctuated Draw vs Predictive - double *llh_drawfluc_pred; + std::vector llh_drawfluc_pred; /// Fluctuated Predictive vs Predictive - double *llh_predfluc_pred; + std::vector llh_predfluc_pred; /// Fluctuated Draw vs Fluctuated Predictive - double *llh_drawfluc_predfluc; + std::vector llh_drawfluc_predfluc; /// Fluctuated Data vs Draw - double *llh_datafluc_draw; + std::vector llh_datafluc_draw; /// Projection X (most likely muon momentum) of LLH - double *llh_data_draw_ProjectX; - double *llh_drawfluc_draw_ProjectX; + std::vector llh_data_draw_ProjectX; + std::vector llh_drawfluc_draw_ProjectX; /// LLH penalty for each throw double llh_penalty; diff --git a/plotting/CMakeLists.txt b/plotting/CMakeLists.txt index 60c73655d..88b1eaba8 100644 --- a/plotting/CMakeLists.txt +++ b/plotting/CMakeLists.txt @@ -7,7 +7,7 @@ foreach(app MatrixPlotter ) add_executable(${app} ${app}.cpp ) - target_link_libraries(${app} MaCh3::Plotting) + target_link_libraries(${app} MaCh3::Plotting MaCh3Warnings) install(TARGETS ${app} DESTINATION ${CMAKE_BINARY_DIR}/bin) endforeach(app) diff --git a/plotting/GetPostfitParamPlots.cpp b/plotting/GetPostfitParamPlots.cpp index b1f02b0bd..b816b37d2 100644 --- a/plotting/GetPostfitParamPlots.cpp +++ b/plotting/GetPostfitParamPlots.cpp @@ -36,6 +36,10 @@ /// /// @note Originally written by Clarence, with changes by Will, updates by Kamil, and converted to a generic plotter by Ewan. +//this file has lots of usage of the ROOT plotting interface that only takes floats, turn this warning off for this CU for now +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#pragma GCC diagnostic ignored "-Wconversion" + MaCh3Plotting::PlottingManager *man; TH1D *Prefit; @@ -114,7 +118,7 @@ void ReadSettings(std::shared_ptr File1) { MACH3LOG_DEBUG("Reading settings for file {}", File1->GetName()); File1->ls(); - TTree *Settings = (TTree*)(File1->Get("Settings")); + TTree *Settings = (File1->Get("Settings")); MACH3LOG_DEBUG("Got settings tree"); Settings->Print(); @@ -146,7 +150,7 @@ void ReadSettings(std::shared_ptr File1) inline TH1D* makeRatio(TH1D *PrefitCopy, TH1D *PostfitCopy, bool setAxes){ // set up the ratio hist - TH1D *Ratio = (TH1D*)PrefitCopy->Clone(); + TH1D* Ratio = static_cast(PrefitCopy->Clone()); Ratio->GetYaxis()->SetTitle("(x_{Post}-#mu_{Prior})/#sigma_{Prior}"); Ratio->SetMinimum(-3.7); Ratio->SetMaximum(3.7); @@ -224,7 +228,7 @@ inline void DrawPlots(TCanvas *plotCanv, TH1D* PrefitCopy, std::vectorPo PrefitCopy->GetYaxis()->SetTitleOffset(1.3); PrefitCopy->Draw("e2"); - for(int fileId=0; fileId < (int)PostfitVec.size(); fileId++){ + for (int fileId = 0; fileId < static_cast(PostfitVec.size()); fileId++) { TH1D *postFitHist = PostfitVec[fileId]; postFitHist->SetMarkerColor(TColor::GetColorPalette(fileId)); @@ -254,7 +258,7 @@ inline void DrawPlots(TCanvas *plotCanv, TH1D* PrefitCopy, std::vectorPo ratioHists.push_back(makeRatio(PrefitCopy, PostfitVec[0], true)); ratioHists[0]->Draw("p"); - for(int postFitIdx = 1; postFitIdx < (int)PostfitVec.size(); postFitIdx++){ + for(int postFitIdx = 1; postFitIdx < static_cast(PostfitVec.size()); postFitIdx++){ ratioHists.push_back(makeRatio(PrefitCopy, PostfitVec[postFitIdx], true)); ratioHists[postFitIdx]->SetMarkerColor(TColor::GetColorPalette(postFitIdx)); @@ -295,7 +299,7 @@ void MakeXsecPlots() { // get the names of the blocks of parameters to group together std::vector const blockNames = man->getOption>("paramGroups"); - const int XsecPlots = (int)blockNames.size(); + const int XsecPlots = static_cast(blockNames.size()); for (int i = 0; i < XsecPlots; i++) { @@ -307,10 +311,11 @@ void MakeXsecPlots() std::vector blockContents = paramBlock[2].as>(); // get num of params in the block - int nParams = (int)blockContents.size(); + const int nParams = static_cast(blockContents.size()); // set some plot things - TH1D *blockHist_prefit = new TH1D(blockName.c_str(), blockTitle.c_str(), nParams, 0.0, (double)nParams); + TH1D *blockHist_prefit = new TH1D(blockName.c_str(), + blockTitle.c_str(), nParams, 0.0, static_cast(nParams)); man->style().setTH1Style(blockHist_prefit, man->getOption("prefitHistStyle")); @@ -323,9 +328,10 @@ void MakeXsecPlots() // now set for the postfit blocks for all files std::vector blockHist_postfit_Vec; - for(int fileId = 0; fileId < man->getNFiles(); fileId++){ + for(unsigned int fileId = 0; fileId < man->getNFiles(); fileId++){ - TH1D *blockHist_postfit = new TH1D((blockName + man->getFileName(fileId)).c_str(), blockTitle.c_str(), nParams, 0.0, (double)nParams); + TH1D *blockHist_postfit = new TH1D((blockName + man->getFileName(fileId)).c_str(), + blockTitle.c_str(), nParams, 0.0, static_cast(nParams)); // loop throught all the parameters in this block and set the contents in the blocks TH1 for(int localBin=0; localBin < nParams; localBin ++){ @@ -357,7 +363,7 @@ void MakeFluxPlots() std::vector const fluxBlockNames = man->getOption>("fluxGroups"); auto const fluxBinningTable = man->getOption("FluxBinning"); - const int FluxPlots = (int)fluxBlockNames.size(); + const int FluxPlots = static_cast(fluxBlockNames.size()); for (int i = 0; i < FluxPlots; i++) { @@ -380,7 +386,7 @@ void MakeFluxPlots() MACH3LOG_CRITICAL(" Should have the form [, ]"); throw MaCh3Exception(__FILE__ , __LINE__ ); } - if(nParams != (int)binning.size() -1){ + if (nParams != static_cast(binning.size()) - 1) { MACH3LOG_CRITICAL("Binning provided for flux param block {} does not match the number of parameters specified for the block", fluxBlockName); MACH3LOG_CRITICAL(" Provided {} parameters but {} bins", nParams, binning.size() -1); throw MaCh3Exception(__FILE__ , __LINE__ ); @@ -400,7 +406,7 @@ void MakeFluxPlots() // now set for the postfit blocks for all files std::vector blockHist_postfit_Vec; - for(int fileId = 0; fileId < man->getNFiles(); fileId++){ + for(unsigned int fileId = 0; fileId < man->getNFiles(); fileId++){ TH1D *blockHist_postfit = new TH1D(fluxBlockName.c_str(), blockTitle.c_str(), nParams, binning.data()); for(int fluxParId = blockContents[0]; fluxParId <= blockContents[1]; fluxParId++){ @@ -440,7 +446,7 @@ void MakeNDDetPlots() NDbinCounter += NDSamplesBins[i]; std::vector PostfitNDDetHistVec(man->getNFiles()); - TH1D *PreFitNDDetHist = (TH1D*)man->input().getFile(0).file->Get(Form("param_%s_prefit", NDSamplesNames[i].c_str())); + TH1D *PreFitNDDetHist = man->input().getFile(0).file->Get(Form("param_%s_prefit", NDSamplesNames[i].c_str())); man->style().setTH1Style(PreFitNDDetHist, man->getOption("prefitHistStyle")); std::string temp = NDSamplesNames[i].c_str(); @@ -452,8 +458,8 @@ void MakeNDDetPlots() MACH3LOG_DEBUG(" Start bin: {} :: End bin: {}", Start, NDbinCounter); // set the x range for the postfits - for(int fileId = 0; fileId < man->getNFiles(); fileId++){ - PostfitNDDetHistVec[fileId] = (TH1D*)man->input().getFile(fileId).file->Get(Form("param_%s_%s", NDSamplesNames[i].c_str(), plotType.c_str())); + for(unsigned int fileId = 0; fileId < man->getNFiles(); fileId++){ + PostfitNDDetHistVec[fileId] = man->input().getFile(fileId).file->Get(Form("param_%s_%s", NDSamplesNames[i].c_str(), plotType.c_str())); } //KS: We dont' need name for every nd param @@ -526,8 +532,6 @@ void MakeOscPlots() void MakeXsecRidgePlots() { - - gStyle->SetPalette(51); TCanvas *blankCanv = new TCanvas ("blankCanv", "blankCanv", 2048, 2048); @@ -535,7 +539,7 @@ void MakeXsecRidgePlots() // get the names of the blocks of parameters to group together std::vector const blockNames = man->getOption>("paramGroups"); - const int XsecPlots = (int)blockNames.size(); + const int XsecPlots = static_cast(blockNames.size()); double padTopMargin = 0.9; double padBottomMargin = 0.1; @@ -544,47 +548,46 @@ void MakeXsecRidgePlots() for (int i = 0; i < XsecPlots; i++) { - // get the configuration for this parameter std::string blockName = blockNames[i]; auto const ¶mBlock = man->getOption(blockName); - std::string blockTitle = paramBlock[0].as(); - std::vector blockLimits = paramBlock[1].as>(); - std::vector blockContents = paramBlock[2].as>(); + auto blockTitle = paramBlock[0].as(); + auto blockLimits = paramBlock[1].as>(); + auto blockContents = paramBlock[2].as>(); // the directory of histograms - TDirectoryFile *posteriorDir = (TDirectoryFile *)man->input().getFile(0).file->Get("Post"); + TDirectoryFile *posteriorDir = man->input().getFile(0).file->Get("Post"); // get num of params in the block - int nParams = (int)blockContents.size(); + int nParams = static_cast(blockContents.size()); TCanvas *ridgeCanv = new TCanvas ("RidgePlotCanv", "RidgePlotCanv", 2048, 2048); ridgeCanv->Divide(1,1+nParams, 0.01, 0.0); - TLatex *title = new TLatex(); + auto title = std::make_unique(); title->SetTextAlign(21); title->SetTextSize(0.03); title->DrawLatex(0.5, 0.95, blockTitle.c_str()); - TLatex *label = new TLatex(); + auto label = std::make_unique(); label->SetTextAlign(31); label->SetTextSize(0.02); - TLine *line = new TLine(); + auto line = std::make_unique(); line->SetLineColor(kBlack); line->SetLineWidth(ridgeLineWidth); // use this to set the limits and also to plot the x axis and grid TH1D *axisPlot = new TH1D("axis plot", "", 1, blockLimits[0], blockLimits[1]); - for(int parId=0; parId < nParams; parId++){ + for(int parId = 0; parId < nParams; parId++){ std::string paramName = blockContents[parId]; - TCanvas *posteriorDistCanv = NULL; - TH1D *posteriorDist = NULL; + TCanvas *posteriorDistCanv = nullptr; + TH1D *posteriorDist = nullptr; // get the list of objects in the directory TIter next(posteriorDir->GetListOfKeys()); - while(TKey *key = (TKey*) next()){ + while (TKey* key = static_cast(next())) { // check if the end of the param name matches with the MaCh3 name, do this so we exclude things like nds_ at the start of the name std::string str(key->GetTitle()); std::string name = man->input().translateName(0, MaCh3Plotting::kPostFit, paramName); @@ -592,19 +595,20 @@ void MakeXsecRidgePlots() bool foundPar = (pos == str.length() - name.length()); if(foundPar){ - posteriorDistCanv = (TCanvas*)posteriorDir->Get(key->GetName()); - posteriorDist = (TH1D*)posteriorDistCanv->GetPrimitive(key->GetName()); + posteriorDistCanv = posteriorDir->Get(key->GetName()); + posteriorDist = static_cast(posteriorDistCanv->GetPrimitive(key->GetName())); } } - if(posteriorDist == NULL){ + if(posteriorDist == nullptr){ MACH3LOG_WARN("Couldnt find parameter {} when making ridgeline plots", paramName); continue; } // EM: do some funky scaling so that we always get evenly spaced pads in the range [bottomMargin, TopMargin] with the specified overlap - double padAnchor = padBottomMargin + ((float)(nParams - parId -1) / (float)(nParams -1)) * (padTopMargin - padBottomMargin); - double padWidth = ((padTopMargin - padBottomMargin) / (float)(nParams)); + double padAnchor = padBottomMargin + (static_cast(nParams - parId - 1) / + static_cast(nParams - 1)) * (padTopMargin - padBottomMargin); + double padWidth = (padTopMargin - padBottomMargin) / static_cast(nParams); double norm = (padTopMargin - padBottomMargin); double padTop = padWidth * (1.0 + padOverlap) * (padTopMargin - padAnchor) / norm + padAnchor; @@ -627,13 +631,14 @@ void MakeXsecRidgePlots() posteriorDist->SetTitle(""); posteriorDist->SetLineWidth(ridgeLineWidth); - TH1D *axisPlot_tmp = (TH1D*)axisPlot->Clone(Form("AxisPlot_%s", paramName.c_str())); + TH1D* axisPlot_tmp = static_cast(axisPlot->Clone(Form("AxisPlot_%s", paramName.c_str()))); axisPlot_tmp->Draw("A"); posteriorDist->Draw("H SAME"); axisPlot_tmp->GetYaxis()->SetRangeUser(0.0, 0.7 *posteriorDist->GetMaximum()); posteriorDist->SetLineColor(kWhite); - posteriorDist->SetFillColorAlpha(TColor::GetColorPalette(floor((float)parId * TColor::GetNumberOfColors()/ (float)nParams)), 0.85); + posteriorDist->SetFillColorAlpha(TColor::GetColorPalette(floor(static_cast(parId) * + TColor::GetNumberOfColors() / static_cast(nParams))), 0.85); posteriorDist->GetXaxis()->SetRangeUser(blockLimits[0], blockLimits[1]); posteriorDist->GetYaxis()->SetTitle(paramName.c_str()); @@ -651,7 +656,7 @@ void MakeXsecRidgePlots() ridgeCanv->cd(); ridgeCanv->SetGrid(1,1); - TPad *axisPad = new TPad("AxisPad", "", 0.3, 0.0, 0.9, 1.0, -1, 0, -1); + TPad *axisPad = new TPad("AxisPad", "", 0.3, 0.0, 0.9, 1.0, -1, 0, -1); axisPad->SetLeftMargin(0.0); axisPad->SetRightMargin(0.0); axisPad->Draw(); @@ -682,7 +687,6 @@ void MakeXsecRidgePlots() ridgeCanv->SaveAs("RidgePlots.pdf"); delete ridgeCanv; } - blankCanv->SaveAs("RidgePlots.pdf]"); } @@ -737,14 +741,13 @@ void GetPostfitParamPlots() p4->SetGrid(); // Make a Legend page - TLegend *leg = new TLegend(0.0, 0.0, 1.0, 1.0); - + auto leg = std::make_unique(0.0, 0.0, 1.0, 1.0); // make a dummy TH1 to set out legend Prefit = new TH1D(); man->style().setTH1Style(Prefit, man->getOption("prefitHistStyle")); leg->AddEntry(Prefit, "Prior", "lpf"); - for(int fileId = 0; fileId < man->getNFiles(); fileId++){ + for(unsigned int fileId = 0; fileId < man->getNFiles(); fileId++){ TH1D *postFitHist_tmp = new TH1D(); postFitHist_tmp->SetBit(kCanDelete); @@ -759,7 +762,6 @@ void GetPostfitParamPlots() canv->Clear(); leg->Draw(); canv->Print((SaveName).c_str()); - delete leg; MakeXsecPlots(); @@ -785,19 +787,18 @@ void GetPostfitParamPlots() inline TGraphAsymmErrors* MakeTGraphAsymmErrors(std::shared_ptr File) { - double* x = new double[nBins]; - double* y = new double[nBins]; - double* exl = new double[nBins]; - double* eyl = new double[nBins]; - double* exh = new double[nBins]; - double* eyh = new double[nBins]; - - TH1D* PostHist = (TH1D*)File->Get( ("param_xsec_"+plotType).c_str() )->Clone(); - - TVectorD* Errors_HPD_Positive = (TVectorD*)File->Get( "Errors_HPD_Positive" )->Clone(); - TVectorD* Errors_HPD_Negative = (TVectorD*)File->Get( "Errors_HPD_Negative" )->Clone(); - - //KS: I am tempted to multithred this... + std::vector x(nBins); + std::vector y(nBins); + std::vector exl(nBins); + std::vector eyl(nBins); + std::vector exh(nBins); + std::vector eyh(nBins); + + TH1D* PostHist = static_cast(File->Get( ("param_xsec_"+plotType).c_str() )->Clone()); + + TVectorD* Errors_HPD_Positive = static_cast(File->Get( "Errors_HPD_Positive" )->Clone()); + TVectorD* Errors_HPD_Negative = static_cast(File->Get( "Errors_HPD_Negative" )->Clone()); + //KS: I am tempted to multithread this... for(int i = 0; i < nBins; ++i) { //KS: We are extrcting value from three object each having different numbering scheme, I have checked carefully so this is correct please don't cahnge all these +1 +0.5 etc. it just work... @@ -810,18 +811,12 @@ inline TGraphAsymmErrors* MakeTGraphAsymmErrors(std::shared_ptr File) eyh[i] = (*Errors_HPD_Positive)(i); eyl[i] = (*Errors_HPD_Negative)(i); } - TGraphAsymmErrors* PostGraph = new TGraphAsymmErrors(nBins,x,y,exl,exh,eyl,eyh); + TGraphAsymmErrors* PostGraph = new TGraphAsymmErrors(nBins, x.data(), y.data(), exl.data(), exh.data(), eyl.data(), eyh.data()); PostGraph->SetTitle(""); delete PostHist; delete Errors_HPD_Positive; delete Errors_HPD_Negative; - delete[] x; - delete[] y; - delete[] exl; - delete[] eyl; - delete[] exh; - delete[] eyh; return PostGraph; } @@ -847,7 +842,7 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") if(PlotAssym) SaveName += "_Assym"; std::shared_ptr File1 = std::make_shared(FileName1.c_str()); - std::shared_ptr File2 = NULL; + std::shared_ptr File2 = nullptr; if(FileName2 != "") File2 = std::make_shared(FileName2.c_str()); canv = new TCanvas("canv", "canv", 1024, 1024); @@ -864,10 +859,10 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") canv->Print((SaveName+".pdf[").c_str()); canv->SetGrid(); - Violin = NULL; - ViolinPre = (TH2D*)File1->Get( "param_violin_prior" ); - Violin = (TH2D*)File1->Get( "param_violin" ); - if(Violin == NULL) + Violin = nullptr; + ViolinPre = File1->Get( "param_violin_prior" ); + Violin = File1->Get( "param_violin" ); + if(Violin == nullptr) { MACH3LOG_ERROR("Couldn't find violin plot, make sure method from MCMCProcessor is being called"); return; @@ -889,7 +884,7 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") Violin->SetMarkerStyle(20); Violin->SetMarkerSize(0.5); - TH1D* Postfit = (TH1D*)File1->Get( ("param_xsec_"+plotType).c_str() ); + TH1D* Postfit = File1->Get( ("param_xsec_"+plotType).c_str() ); Postfit->SetMarkerColor(kRed); Postfit->SetLineColor(kRed); Postfit->SetMarkerStyle(7); @@ -900,9 +895,9 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") PostGraph->SetMarkerStyle(7); PostGraph->SetLineWidth(2); PostGraph->SetLineStyle(kSolid); - if(File2 != NULL) + if(File2 != nullptr) { - Violin2 = (TH2D*)File2->Get( "param_violin" ); + Violin2 = File2->Get( "param_violin" ); Violin2->SetMarkerColor(kGreen); Violin2->SetLineColor(kGreen); Violin2->SetFillColor(kGreen); @@ -910,18 +905,17 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") } // Make a Legend page - TLegend *leg = new TLegend(0.0, 0.0, 1.0, 1.0); - if (ViolinPre != NULL) leg->AddEntry(ViolinPre, "Prior", "lpf"); - if (Violin != NULL) leg->AddEntry(Violin, "Posterior", "lpf"); - if (Violin2 != NULL) leg->AddEntry(Violin2, "Second Violin", "lpf"); - if(PlotAssym) leg->AddEntry(PostGraph, "HPD Assym", "lp"); - else leg->AddEntry(Postfit, "HPD", "lpf"); + auto leg = std::make_unique(0.0, 0.0, 1.0, 1.0); + if (ViolinPre != nullptr) leg->AddEntry(ViolinPre, "Prior", "lpf"); + if (Violin != nullptr) leg->AddEntry(Violin, "Posterior", "lpf"); + if (Violin2 != nullptr) leg->AddEntry(Violin2, "Second Violin", "lpf"); + if(PlotAssym) leg->AddEntry(PostGraph, "HPD Assym", "lp"); + else leg->AddEntry(Postfit, "HPD", "lpf"); canv->cd(); canv->Clear(); leg->Draw(); canv->Print((SaveName+".pdf").c_str()); - delete leg; // Do some fancy replacements PrettifyTitles(ViolinPre); @@ -939,47 +933,47 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") Violin->GetXaxis()->SetRangeUser(XsecOffset[i-1], XsecOffset[i]); Violin->GetYaxis()->SetRangeUser(-1.5, 2.5); - if(File2 != NULL) Violin2->GetYaxis()->SetRangeUser(-1.5, 2.5); + if(File2 != nullptr) Violin2->GetYaxis()->SetRangeUser(-1.5, 2.5); std::string name = ViolinPre->GetTitle(); if (name.find("SPP") != std::string::npos) { ViolinPre->GetYaxis()->SetRangeUser(-5., 25.); //For RES Eb we need huge range Violin->GetYaxis()->SetRangeUser(-5., 25.); - if(File2 != NULL) Violin2->GetYaxis()->SetRangeUser(-5., 25.); + if(File2 != nullptr) Violin2->GetYaxis()->SetRangeUser(-5., 25.); } else if (name.find("CCQE Binding Energy") != std::string::npos) { ViolinPre->GetYaxis()->SetRangeUser(-10., 16.); //For Eb we need bigger range Violin->GetYaxis()->SetRangeUser(-10., 16.); //For Eb we need bigger range - if(File2 != NULL) Violin2->GetYaxis()->SetRangeUser(-10., 16.); + if(File2 != nullptr) Violin2->GetYaxis()->SetRangeUser(-10., 16.); } else if (name.find("FSI") != std::string::npos) { ViolinPre->GetYaxis()->SetRangeUser(-0.5, 2.8); Violin->GetYaxis()->SetRangeUser(-0.5, 2.8); - if(File2 != NULL) Violin2->GetYaxis()->SetRangeUser(-0.5, 2.8); + if(File2 != nullptr) Violin2->GetYaxis()->SetRangeUser(-0.5, 2.8); } else if (name.find("CCQE") != std::string::npos) { ViolinPre->GetYaxis()->SetRangeUser(-2., 3.); Violin->GetYaxis()->SetRangeUser(-2., 3.); - if(File2 != NULL) Violin2->GetYaxis()->SetRangeUser(-2., 3.); + if(File2 != nullptr) Violin2->GetYaxis()->SetRangeUser(-2., 3.); } //KS: ROOT6 has some additional options, consider updating it. more https://root.cern/doc/master/classTHistPainter.html#HP140b ViolinPre->Draw("violinX(03100300)"); Violin->Draw("violinX(03100300) SAME"); - if(File2 != NULL) + if(File2 != nullptr) { Violin2->GetXaxis()->SetRangeUser(XsecOffset[i-1], XsecOffset[i]); Violin2->GetXaxis()->SetRangeUser(XsecOffset[i-1], XsecOffset[i]); Violin2->Draw("violinX(03100300) SAME"); } - if (Postfit != NULL) { + if (Postfit != nullptr) { Postfit->GetXaxis()->SetRangeUser(XsecOffset[i-1], XsecOffset[i]); if(!PlotAssym) Postfit->Draw("SAME"); } - if (PostGraph != NULL) { + if (PostGraph != nullptr) { PostGraph->GetXaxis()->SetRangeUser(XsecOffset[i-1], XsecOffset[i]); if(PlotAssym) PostGraph->Draw("P SAME"); } @@ -1008,7 +1002,7 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") //KS: ROOT6 has some additional options, consider updaiting it. more https://root.cern/doc/master/classTHistPainter.html#HP140b ViolinPre->Draw("violinX(03100300)"); Violin->Draw("violinX(03100300) SAME"); - if(File2 != NULL) + if(File2 != nullptr) { Violin2->GetYaxis()->SetRangeUser(0.7, 1.3); Violin2->GetXaxis()->SetRangeUser(nFlux, nFlux+FluxInterval); @@ -1033,7 +1027,7 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") //KS: ROOT6 has some additional options, consider updating it. more https://root.cern/doc/master/classTHistPainter.html#HP140b ViolinPre->Draw("violinX(03100300)"); Violin->Draw("violinX(03100300) SAME"); - if(File2 != NULL) + if(File2 != nullptr) { Violin2->GetXaxis()->SetRangeUser(CrossSectionParameters, nBins); Violin2->GetYaxis()->SetRangeUser(-3.4, 3.4); @@ -1047,11 +1041,11 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") delete canv; delete ViolinPre; delete Violin; - if(Violin2 != NULL) delete Violin2; + if(Violin2 != nullptr) delete Violin2; delete Postfit; delete PostGraph; File1->Close(); - if(File2 != NULL) + if(File2 != nullptr) { File2->Close(); } @@ -1059,45 +1053,45 @@ void GetViolinPlots(std::string FileName1 = "", std::string FileName2 = "") int main(int argc, char *argv[]) { - SetMaCh3LoggerFormat(); - - man = new MaCh3Plotting::PlottingManager(); - man->parseInputs(argc, argv); - std::cout << std::endl << std::endl << "====================" << std::endl; - man->input().getFile(0).file->ls(); - - man->setExec("GetPostfitParamPlots"); + SetMaCh3LoggerFormat(); - man->style().setPalette(man->getOption("colorPalette")); + man = new MaCh3Plotting::PlottingManager(); + man->parseInputs(argc, argv); + std::cout << std::endl << std::endl << "====================" << std::endl; + man->input().getFile(0).file->ls(); - GetPostfitParamPlots(); + man->setExec("GetPostfitParamPlots"); - if (argc != 2 && argc != 3 && argc !=4) - { - MACH3LOG_CRITICAL("Invalid command line options specified"); - MACH3LOG_CRITICAL("How to use: {} .root", argv[0]); - MACH3LOG_CRITICAL("You can add up to 3 different files"); - throw MaCh3Exception(__FILE__, __LINE__); - } - - if (argc == 2) - { - std::string filename = argv[1]; - GetViolinPlots(filename); - } - else if (argc == 3) - { - std::string filename1 = argv[1]; - std::string filename2 = argv[2]; - GetViolinPlots(filename1, filename2); - } - else if (argc == 4) - { - std::string filename1 = argv[1]; - std::string filename2 = argv[2]; - std::string filename3 = argv[3]; - //KS: Violin plot currently not supported by three file veriosn although it should be super easy to adapt - } + man->style().setPalette(man->getOption("colorPalette")); + + GetPostfitParamPlots(); + + if (argc != 2 && argc != 3 && argc !=4) + { + MACH3LOG_CRITICAL("Invalid command line options specified"); + MACH3LOG_CRITICAL("How to use: {} .root", argv[0]); + MACH3LOG_CRITICAL("You can add up to 3 different files"); + throw MaCh3Exception(__FILE__, __LINE__); + } + + if (argc == 2) + { + std::string filename = argv[1]; + GetViolinPlots(filename); + } + else if (argc == 3) + { + std::string filename1 = argv[1]; + std::string filename2 = argv[2]; + GetViolinPlots(filename1, filename2); + } + else if (argc == 4) + { + std::string filename1 = argv[1]; + std::string filename2 = argv[2]; + std::string filename3 = argv[3]; + //KS: Violin plot currently not supported by three file veriosn although it should be super easy to adapt + } return 0; } diff --git a/plotting/MatrixPlotter.cpp b/plotting/MatrixPlotter.cpp index 9556feefc..711757a6e 100755 --- a/plotting/MatrixPlotter.cpp +++ b/plotting/MatrixPlotter.cpp @@ -12,6 +12,10 @@ #include "plottingUtils/plottingUtils.h" #include "plottingUtils/plottingManager.h" +//this file has lots of usage of the ROOT plotting interface that only takes floats, turn this warning off for this CU for now +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#pragma GCC diagnostic ignored "-Wconversion" + TH2D* GetSubMatrix(TH2D *MatrixFull, const std::string& Title, const std::vector& Params) { std::vector ParamIndex(Params.size(), -999); diff --git a/plotting/PlotLLH.cpp b/plotting/PlotLLH.cpp index 82785a4dd..e1e1a90b3 100644 --- a/plotting/PlotLLH.cpp +++ b/plotting/PlotLLH.cpp @@ -5,6 +5,10 @@ // MACH3 PLOTTING #include "plottingUtils/plottingManager.h" +//this file has lots of usage of the ROOT plotting interface that only takes floats, turn this warning off for this CU for now +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#pragma GCC diagnostic ignored "-Wconversion" + // some options for the plots double ratioPlotSplit; double yTitleOffset; @@ -22,7 +26,6 @@ void getSplitSampleStack(int fileIdx, std::string parameterName, TH1D LLH_allSam THStack *sampleStack, TLegend *splitSamplesLegend, float baselineLLH_main = 0.00001) { - std::vector sampNames = man->input().getTaggedSamples(man->getOption>("sampleTags")); size_t nSamples = sampNames.size(); @@ -57,9 +60,9 @@ void getSplitSampleStack(int fileIdx, std::string parameterName, TH1D LLH_allSam LLH_indivSam->SetStats(0); LLH_indivSam->SetLineColor(TColor::GetColorPalette( - floor((float)i * TColor::GetNumberOfColors() / (float)nSamples))); + floor(static_cast(i) * TColor::GetNumberOfColors() / static_cast(nSamples)))); LLH_indivSam->SetFillColor(TColor::GetColorPalette( - floor((float)i * TColor::GetNumberOfColors() / (float)nSamples))); + floor(static_cast(i) * TColor::GetNumberOfColors() / static_cast(nSamples)))); sampleStack->Add(LLH_indivSam); splitSamplesLegend->AddEntry(LLH_indivSam, man->style().prettifySampleName(sampName).c_str(), "lf"); @@ -81,13 +84,9 @@ void getSplitSampleStack(int fileIdx, std::string parameterName, TH1D LLH_allSam { drawLabel[i] = false; } - MACH3LOG_DEBUG(" drawLabel = {}", drawLabel.back()); MACH3LOG_DEBUG(""); - } - - return; } // handy dandy helper function for drawing and nicely formatting a stack of ratio plots @@ -124,7 +123,7 @@ void makeLLHScanComparisons(std::string paramName, std::string LLHType, std::str // will use these to make comparisons of likelihoods THStack *compStack = new THStack("LLH_Stack", ""); THStack *ratioCompStack = new THStack("LLH_Ratio_Stack", ""); - TLegend *legend = new TLegend(0.3, 0.6, 0.7, 0.8); + auto legend = std::make_unique(0.3, 0.6, 0.7, 0.8); // get the sample reweight hist from the main file TH1D LLH_main = man->input().getLLHScan_TH1D(0, paramName, LLHType); @@ -137,9 +136,8 @@ void makeLLHScanComparisons(std::string paramName, std::string LLHType, std::str int nBins = LLH_main.GetNbinsX(); // go through the other files - for (int extraFileIdx = 1; extraFileIdx < man->input().getNInputFiles(); extraFileIdx++) + for (unsigned int extraFileIdx = 1; extraFileIdx < man->input().getNInputFiles(); extraFileIdx++) { - TH1D *compHist = new TH1D(man->input().getLLHScan_TH1D(extraFileIdx, paramName, LLHType)); compHist->SetBit(kCanDelete); // <- will allow this to be deleted by root once done plotting if (compHist->GetNbinsX() == 0) @@ -147,12 +145,12 @@ void makeLLHScanComparisons(std::string paramName, std::string LLHType, std::str // make them look different to each other compHist->SetLineColor( - TColor::GetColorPalette(floor((float)extraFileIdx * TColor::GetNumberOfColors() / - (float)man->input().getNInputFiles()))); + TColor::GetColorPalette(floor(static_cast(extraFileIdx) * TColor::GetNumberOfColors() / + static_cast(man->input().getNInputFiles())))); compHist->SetLineStyle(2 + extraFileIdx % 9); compHist->SetLineWidth(lineWidth); - TH1D *divHist = (TH1D *)compHist->Clone(Form("RatioHist_%i", extraFileIdx)); + TH1D *divHist = static_cast(compHist->Clone(Form("RatioHist_%i", extraFileIdx))); divHist->SetBit(kCanDelete); if (man->getPlotRatios()) @@ -206,12 +204,10 @@ void makeLLHScanComparisons(std::string paramName, std::string LLHType, std::str // will use these to make comparisons of likelihoods delete compStack; delete ratioCompStack; - delete legend; } void makeSplitSampleLLHScanComparisons(std::string paramName, std::string outputFileName, TCanvas *canv, TPad *LLHPad, TPad *ratioPad) { - MACH3LOG_DEBUG(" Making split sample LLH comparison"); canv->Clear(); canv->Draw(); @@ -276,7 +272,7 @@ void makeSplitSampleLLHScanComparisons(std::string paramName, std::string output } // now we plot the comparisson file plots - for (int extraFileIdx = 1; extraFileIdx < man->getNFiles(); extraFileIdx++) + for (unsigned int extraFileIdx = 1; extraFileIdx < man->getNFiles(); extraFileIdx++) { MACH3LOG_DEBUG(" - Adding plot for additional file {}", extraFileIdx); canv->cd(1 + extraFileIdx); @@ -342,7 +338,6 @@ void makeSplitSampleLLHScanComparisons(std::string paramName, std::string output if (man->getPlotRatios()) { - THStack *splitSamplesStackRatios = new THStack(paramName.c_str(), ""); TList *baselineHistList = baseSplitSamplesStack->GetHists(); @@ -356,10 +351,10 @@ void makeSplitSampleLLHScanComparisons(std::string paramName, std::string output Form("%s_%s_splitDiv", paramName.c_str(), man->getFileLabel(extraFileIdx).c_str()), compLLH_main.GetNbinsX(), compLLH_main.GetBinLowEdge(1), compLLH_main.GetBinLowEdge(compLLH_main.GetNbinsX() + 1)); - divHist->Divide((TH1D *)compHistList->At(sampleIdx), - (TH1D *)baselineHistList->At(sampleIdx)); + divHist->Divide(static_cast(compHistList->At(sampleIdx)), + static_cast(baselineHistList->At(sampleIdx))); splitSamplesStackRatios->Add(divHist); - divHist->SetLineColor(((TH1D *)compHistList->At(sampleIdx))->GetLineColor()); + divHist->SetLineColor((static_cast(compHistList->At(sampleIdx))->GetLineColor())); } canv->cd(2 + extraFileIdx);