Generalized folding code

.gitignore

@@ -17,3 +17,7 @@ Pipfile.lock
 *.run
 *.DS_Store
 .vscode/settings.json
+*.save
+*.tar.gz
+cpptools/src/TGlauberMC/
+cpptools/src/TennGen/

cpptools/src/rutil/rutil.hh

@@ -19,7 +19,7 @@ namespace RUtil
         virtual ~Test() {;}
         void setMember(Double_t v) {fMember = v;}
         Double_t getMember() {return fMember;}
-        
+
     private:
         Double_t fMember;
 
@@ -41,13 +41,15 @@ namespace RUtil
                               const double & prior_variation_parameter);
     double prior_scale_func_4(const double & obs_true, const double & content,
                               const double & prior_variation_parameter);
+    double prior_scale_func_5(const double & obs_true, const double & content,
+                              const double & prior_variation_parameter);
     double prior_scale_func_def(const double & obs_true, const double & content,
                                 const double & prior_variation_parameter);
 
-	void delete_h(TH2* h);
-	void delete_h(THn* h);
+    void delete_h(TH2* h);
+    void delete_h(THn* h);
 
-	bool* sorted_match(const int* a, const int a_len, const int* b, const int b_len);
+    bool* sorted_match(const int* a, const int a_len, const int* b, const int b_len);
 
     //------------------------------------------------------
     // Rebinning utilities
@@ -57,7 +59,7 @@ namespace RUtil
         HistUtils() : TObject()
         {;}
         virtual ~HistUtils() {;}
-        
+
         // Rebin 2D histogram h with name hname using axes given by x_bins and y_bins
         TH2F* rebin_th2(TH2F & h_to_rebin, char* hname, double* x_bins, int n_x_bins,
                         double* y_bins, int n_y_bins, bool move_y_underflow = false);
@@ -88,19 +90,71 @@ namespace RUtil
                         const bool use_miss_fake=false,
                         const bool do_roounfoldresponse=true);
 
-		//------------------------------------------------------
-		// Convolution of nonperturbative shape functions
-
-		// Create and return 2D histogram, convolving h with shape function
-		TH2D* convolve_F_np(const double & Omega, const double & R, const double & beta,
-							const double* ob_bins, const int & n_ob_bins, const double* obs,
-							const double* ob_bin_width,
-							const double* pT_bins, const int & n_pT_bins, const double* pTs,
-							const TH2D & h, const std::string & name, const bool groomed = false,
-							const double & sd_beta = 0, const double & sd_zcut = 0.2,
-							const std::string & option = "");
-
-		double find_cell(double val, const double * cell, const int range, bool phi);
+        // Rebin N-dimensional THn to a new histogram with name name_thn_rebinned using provided axes
+        // WARNING: currently requires n_dim = 4
+        // Alternate definition to take TH2 as prior variation
+        THnF* rebin_thn_th2prior(
+            const std::string & response_file_name,
+            const THnF* thn,
+            const std::string & name_thn_rebinned,
+            const std::string & name_roounfold,
+            const int & n_dim,
+            const int & n_pt_bins_det,
+            const double* det_pt_bin_array,
+            const int & n_obs_bins_det,
+            const double* det_bin_array,
+            const int & n_pt_bins_truth,
+            const double* truth_pt_bin_array,
+            const int & n_obs_bins_truth,
+            const double* truth_bin_array,
+            const std::string & label="",
+            const double & prior_variation_parameter=0.,
+            const TH2* prior_variation=nullptr,
+            const bool move_underflow=false,
+            const bool use_miss_fake=false,
+            const bool do_roounfoldresponse=true);
+
+        //---------------------------------------------------------------
+        // Remove outliers from a TH1 via "simple" method:
+        //   delete any bin contents with N counts < limit
+        // Modifies histogram in-place and returns its pointer
+        //---------------------------------------------------------------
+        TH1* simpleRemoveOutliers(TH1* hist, bool verbose, int limit);
+
+        //---------------------------------------------------------------
+        // Remove outliers from a TH1 via pT-hat method:
+        //   delete any bin contents with pT > limit
+        // Modifies histogram in-place and returns its pointer
+        //---------------------------------------------------------------
+        TH1* pThatRemoveOutliers(TH1* hist, bool verbose, const double & limit);
+
+        //---------------------------------------------------------------
+        // Remove outliers from a THn via "simple" method:
+        //   delete any bin contents with N counts < limit
+        // Modifies histogram in-place and returns its pointer
+        //---------------------------------------------------------------
+        THn* simpleRemoveOutliersTHn(THn* hist, bool verbose, int limit, int dim);
+
+        //---------------------------------------------------------------
+        // Remove outliers from a THn via pT-hat method:
+        //   delete any bin contents with pT_truth > limit
+        // Modifies histogram in-place and returns its pointer
+        //---------------------------------------------------------------
+        THn* pThatRemoveOutliers(THn* hist, bool verbose, const double & limit, int dim, int pTdim);
+
+        //------------------------------------------------------
+        // Convolution of nonperturbative shape functions
+
+        // Create and return 2D histogram, convolving h with shape function
+        TH2D* convolve_F_np(const double & Omega, const double & R, const double & beta,
+                            const double* ob_bins, const int & n_ob_bins, const double* obs,
+                            const double* ob_bin_width,
+                            const double* pT_bins, const int & n_pT_bins, const double* pTs,
+                            const TH2D & h, const std::string & name, const bool groomed = false,
+                            const double & sd_beta = 0, const double & sd_zcut = 0.2,
+                            const std::string & option = "");
+
+        double find_cell(double val, const double * cell, const int range, bool phi);
 
     private:
         // Create empty THn using provided axes
@@ -131,18 +185,45 @@ namespace RUtil
                                const bool move_underflow=false,
                                const bool use_miss_fake=false);
 
+        // Fill empty thn_rebinned with data from thn
+        // Overloaded function to take TH2 as prior variation
+        void fill_rebinned_thn(const std::string & response_file_name, const THnF* thn,
+                               THnF* thn_rebinned, const unsigned int & n_dim,
+                               const bool do_roounfoldresponse=true,
+                               RooUnfoldResponse* roounfold_response=nullptr,
+                               const float min_det_pt=0.,
+                               const float min_truth_pt=0.,
+                               const float min_det=0.,
+                               const float min_truth=0.,
+                               const float max_det_pt=0.,
+                               const float max_truth_pt=0.,
+                               const float max_det=0.,
+                               const float max_truth=0.,
+                               const TH2* prior_variation=nullptr,
+                               const double & prior_variation_parameter=0.,
+                               const bool move_underflow=false,
+                               const bool use_miss_fake=false);
+
         // Set scaling of prior
         prior_scale_func prior_scale_factor_obs(const int & option);
 
-		//------------------------------------------------------
-		// Convolution of nonperturbative shape functions
+        // Recursive helper function for simpleRemoveOutliersTHn()
+        void simpleRemoveOutliersTHn_recurse(
+            THn* hist, int limit, int dim, int* n_bins, int* x, int dim_to_update);
+        // Recursive helper function for pThatRemoveOutliersTHn()
+        //void pThatRemoveOutliersTHn_recurse(
+        //   THn* hist, int limit, int dim, int pTdim, int max_bin_number,
+        //   int* n_bins, int* x, int dim_to_update, bool verbose);
+
+        //------------------------------------------------------
+        // Convolution of nonperturbative shape functions
 
         // Non-perturbative parameter with factored-out beta dependence
         // Omega is Omega_{a=0} == Omega_{beta=2}  (universal[?])
-		inline double Omega_beta(const double & Omega, const double & beta);
+        inline double Omega_beta(const double & Omega, const double & beta);
 
-		// Shape function for convolving nonperturbative effects
-		inline double F_np(const double & Omega, const double & k, const double & beta);
+        // Shape function for convolving nonperturbative effects
+        inline double F_np(const double & Omega, const double & k, const double & beta);
 
     ClassDef(HistUtils, 1)
     };

pyjetty/alice_analysis/analysis/user/ang/PbPb/pT_cut_comp.cxx

@@ -0,0 +1,66 @@
+// ROOT Macro for comparing ROOT THn before and after pT-hat based cut
+// Written by Ezra Lesser, August 2022
+#include <string>
+#include "THn.h"
+#include "TH2.h"
+
+int pT_cut_comp() {
+
+  // Open up trimmed and untrimmed ROOT files
+  std::string base_dir = "/rstorage/alice/AnalysisResults/ang/877553/";
+  std::string untrimmed_dir = base_dir + "Scaled_no_cut/";
+  std::string trimmed_dir = base_dir + "pt_trimmed_piecewise/";
+  std::string file_name = "AnalysisResults.root";
+
+  // Jet parameters
+  std::string obs = "ang";
+  std::string jetR = "0.2";  // "0.4";
+  std::string Rmax = "0.1";  // "0.25";
+  //std::vector<std::string> alphas = {"1", "1.5", "2", "3"};
+  std::vector<std::string> alphas = {"1"};
+
+  for (int pThatBin = 1; pThatBin <= 20; pThatBin++) { // loop over pT-hat bins
+    // Open ROOT files
+    TFile* file_untrimmed = new TFile((untrimmed_dir + std::to_string(pThatBin) + '/' + file_name).c_str(), "r");
+    TFile* file_trimmed = new TFile((trimmed_dir + std::to_string(pThatBin) + '/' + file_name).c_str(), "r");
+
+    for (const std::string & alpha : alphas) {
+      // Retrieve RM for making projections
+      std::string RM_name =
+        "hResponse_JetPt_" + obs + "_R" + jetR + '_' + alpha + "_Rmax" + Rmax + "Scaled";
+      THn* RM_untrimmed = (THn*) file_untrimmed->Get(RM_name.c_str());
+      RM_untrimmed->SetNameTitle("untrimmed", "untrimmed");
+      THn* RM_trimmed = (THn*) file_trimmed->Get(RM_name.c_str());
+
+      // Project out the x-y (pT) distributions
+      TH2D* pT_untrimmed = RM_untrimmed->Projection(1, 0);
+      TH2D* pT_trimmed = RM_trimmed->Projection(1, 0);
+
+      // Create ratio trimmed/untrimmed and write to file
+      TH2D* pT_ratio = (TH2D*) pT_trimmed->Clone((RM_name + "_ratio").c_str());
+      pT_ratio->Divide(pT_untrimmed);
+
+      const char canvas_name[] = "pT ratio trimmed/untrimmed";
+      TCanvas c(canvas_name, canvas_name, 700, 600);
+      c.Draw();
+
+      const char pad_name[] = "the pad";
+      TPad pad(pad_name, pad_name, 0, 0, 1, 1);
+      pad.SetTicks(1, 1);
+      pad.SetLeftMargin(0.1);
+      pad.SetRightMargin(0.05);
+      pad.Draw();
+
+      pT_ratio->Draw("colz");
+      gStyle->SetOptStat(0);
+      c.SaveAs((RM_name + "_ratio_pThat" + std::to_string(pThatBin) + ".pdf").c_str());
+      delete pT_untrimmed, pT_trimmed, pT_ratio;
+    }
+
+    file_untrimmed->Close();
+    file_trimmed->Close();
+    delete file_untrimmed, file_trimmed;
+  }
+
+  return 0;
+}

pyjetty/alice_analysis/analysis/user/ang/PbPb/plotting_utils_ang.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python3
+
+"""
+  Plotting utilities for jet substructure analysis with track dataframe.
+
+  Author: James Mulligan ([email protected])
+"""
+
+from __future__ import print_function
+
+# General
+import os
+import sys
+import math
+import yaml
+
+# Data analysis and plotting
+import numpy as np
+from array import *
+import ROOT
+
+# Base class
+from pyjetty.alice_analysis.analysis.user.james import plotting_utils_base
+
+################################################################
+class PlottingUtils(plotting_utils_base.PlottingUtilsBase):
+
+  #---------------------------------------------------------------
+  # Constructor
+  #---------------------------------------------------------------
+  def __init__(self, output_dir = '.', config_file = '', R_max = None, thermal = False, groomer_studies = False, **kwargs):
+    super(PlottingUtils, self).__init__(output_dir, config_file, R_max, thermal, groomer_studies, **kwargs)
+
+    print(self)
+