Merge pull request #10 from lcorcodilos/dev

Development of C++ module infrastructure, doc editing, and other re-organizing

.gitignore

@@ -6,4 +6,7 @@
 __pycache__
 build
 *.egg-info
-dist
+dist
+experimental/
+.vscode/
+examples/*.root

README.md

@@ -6,13 +6,14 @@ Default arguments assume the use of the NanoAOD format but any ROOT TTree can be
 
 ## Quick install
 Python 3 is recommended since 2.7 is now no longer supported. Remember to make sure your 
-ROOT version has been built with python3 compatibility.
+ROOT version has been built with Python 3 compatibility. For information on how to do this, see
+[this explanation](doxysetup/Python3.md).
 Working in a virtual environment is also recommended. Below are the commands for using virtualenv but
 you're obviously free to use your favorite tool for the job (you can install virtualenv for Python 3 with 
-`pip3 install virtualenv`).
+`pip install virtualenv` (`pip3` for Python 3)).
 
 ```
-virtualenv timber-env
+python -m virtualenv timber-env
 source timber-env/bin/activate
 git clone https://github.com/lcorcodilos/TIMBER.git
 cd TIMBER

TIMBER/Analyzer.py

@@ -119,10 +119,35 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs"):
                 if self.preV6: self.genEventCount+= RunChain.genEventCount
                 else: self.genEventCount+= RunChain.genEventCount_
 
+        # Get LHAID from LHEPdfWeights branch
+        self.lhaid = "-1"
+        if not self.isData:
+            pdfbranch = self.__eventsChain.GetBranch("LHEPdfWeight")
+            if pdfbranch != None:
+                branch_title = pdfbranch.GetTitle()
+                self.lhaid = ''
+                for c in branch_title:
+                    if c.isdigit(): 
+                        self.lhaid+=str(c)
+                    elif self.lhaid == '':
+                        continue
+                    else:
+                        break
+            self.lhaid = str(int(self.lhaid)-1)
+            print ('LHA ID: '+self.lhaid)
+
         # Cleanup
         del RunChain
         self.ActiveNode = self.BaseNode
 
+    @property
+    def DataFrame(self):
+        '''
+        Returns:
+            RDataFrame: Dataframe for the active node
+        '''        
+        return self.ActiveNode.DataFrame
+
     def SetActiveNode(self,node):
         """Sets the active node.
 
@@ -1223,6 +1248,8 @@ def LoadColumnNames(source=None):
     else:
         file = source
     f = open(file,'r')
-    cols = f.readlines()
+    cols = []
+    for c in f.readlines():
+        cols.append(c.strip('\n'))
     f.close()
     return cols

TIMBER/Framework/README.md

@@ -19,7 +19,7 @@ The name of the file is used to determine which type of output to expect for tem
 An output of the first type should have a file name that end in `_SF.cc` or `_weight.cc` and outputs of
 of the second type should have a file name that ends in `_uncert.cc`.
 
-Outputing a vector of vectors is not currently supported. One might want to do this
+Outputting a vector of vectors is not currently supported. One might want to do this
 if, for example, all jets in a vector need the same 
 correction in which case a vector of `[nominal,up,down]` vectors could be calculated where the index of each 
 `[nominal,up,down]` matches the index of the matching jet in the jet vector. While this would be 

TIMBER/Framework/include/Pythonic.h

@@ -0,0 +1,73 @@
+#include <string>
+#include <sstream>
+#include <algorithm>
+#include <iterator>
+#include <vector>
+#include <stdexcept>
+
+//https://stackoverflow.com/questions/13152252/is-there-a-compact-equivalent-to-python-range-in-c-stl
+template <typename IntType>
+std::vector<IntType> range(IntType start, IntType stop, IntType step)
+{
+  if (step == IntType(0))
+  {
+    throw std::invalid_argument("step for range must be non-zero");
+  }
+
+  std::vector<IntType> result;
+  IntType i = start;
+  while ((step > 0) ? (i < stop) : (i > stop))
+  {
+    result.push_back(i);
+    i += step;
+  }
+
+  return result;
+}
+
+template <typename IntType>
+std::vector<IntType> range(IntType start, IntType stop)
+{
+  return range(start, stop, IntType(1));
+}
+
+template <typename IntType>
+std::vector<IntType> range(IntType stop)
+{
+  return range(IntType(0), stop, IntType(1));
+}
+
+// Adapted from http://www.martinbroadhurst.com/how-to-split-a-string-in-c.html
+std::vector<std::string> split(const std::string& str, char delim = ' ') {
+    std::vector<std::string> out {};
+    std::stringstream ss(str);
+    std::string token;
+    while (std::getline(ss, token, delim)) {
+        out.push_back(token);
+    }
+
+    return out;
+}
+
+// Personal
+template<typename T>
+bool InList(T obj, std::vector<T> list) {
+    auto pos = std::find(std::begin(list), std::end(list), obj);
+    if (pos != std::end(list)){
+        return true;
+    } else {return false;}
+}
+
+bool InString(std::string sub, std::string main) {
+    auto found = main.find(sub);
+    if (found != std::string::npos){
+        return true;
+    } else {return false;}
+}
+
+template<typename T>
+void Extend(std::vector<T> base, std::vector<T> extension) {
+    for (int i = 0; i < extension.size(); i++) {
+        base.push_back(extension.at(i));
+    }
+}

TIMBER/Framework/src/Collection.cc

@@ -0,0 +1,18 @@
+#include <map>
+#include <string>
+#include <ROOT/RVec.hxx>
+
+using namespace ROOT::VecOps;
+using namespace std;
+
+/** Structure to store maps of the various types of objects
+ * in a collection. UChar not considered.
+ * Use by building each map as <branchName, branchValue> and
+ * then assigning to the correct struct member. */
+struct Collection {
+    map<string,int*> Int;
+    map<string,bool*> Bool;
+    map<string,RVec<float>*> RVecFloat;
+    map<string,RVec<int>*> RVecInt;
+};
+

TIMBER/Framework/src/PDFweight_uncert.cc

@@ -0,0 +1,77 @@
+#include <cmath>
+#include <vector>
+#include <numeric>
+#include <fstream>
+#include "ROOT/RVec.hxx"
+#include </home/lucas/Projects/RDFanalyzer/TIMBER/TIMBER/Framework/include/Pythonic.h>
+
+using namespace ROOT::VecOps;
+
+class PDFweight_uncert
+{
+private:
+    bool hessian;
+    bool ignoreEmpty;
+    int lhaid;
+public:
+    PDFweight_uncert(int lhaID, bool ignoreEmptyBranch = false);
+    ~PDFweight_uncert();
+    std::vector<float> eval(RVec<float> LHEPdfWeight);
+};
+
+PDFweight_uncert::PDFweight_uncert(int lhaID, bool ignoreEmptyBranch){
+    ignoreEmpty = ignoreEmptyBranch;
+    lhaid = lhaID;
+    // auto determine whether these are replicas or hessian eigenvectors
+    // Search pdfsets.index, find matching lhaid, check if "hessian" is in the pdf set name
+    std::fstream lhaid_file;
+    lhaid_file.open(std::string(std::getenv("TIMBERPATH"))+"TIMBER/data/pdfsets.index",std::fstream::in);
+    std::string line;
+    std::vector<std::string> line_parts;
+    int idx = 0;
+    while (getline(lhaid_file, line)) {
+        line_parts = split(line,' ');
+        if (lhaid == (int)std::stoi(line_parts[0])) {
+            std::cout << "Found matching PDF set: " << line_parts[1] << " (" << lhaid << ")" << std::endl;
+            if (InString("hessian",line_parts[1])) {
+                hessian = true;
+            } else { hessian = false; }
+        };
+    }
+};
+
+PDFweight_uncert::~PDFweight_uncert(){};
+
+std::vector<float> PDFweight_uncert::eval(RVec<float> LHEPdfWeight) {
+    // [up,down]
+    std::vector<float> v;
+    float stddev;
+    float sumsquares = 0.0;
+    int size = LHEPdfWeight.size();
+    // check weights aren't empty (known bug that they all could)
+    if (size == 0 && !ignoreEmpty) {
+        throw "LHEPdfWeight vector empty. May be known bug in NanoAOD - see https://github.com/cms-nanoAOD/cmssw/issues/520. To ignore, set ignoreEmpty argument to true.";
+    }
+
+    if (hessian) { // Computes sqrt of sum of differences squared
+        float base_eigenv = LHEPdfWeight[0];
+        for (int ipdf = 1; ipdf < size; ipdf++) {
+            sumsquares = sumsquares + std::pow(LHEPdfWeight[ipdf] - base_eigenv,2);
+        }
+        stddev = sqrt(sumsquares);
+
+    } else { // Computes the std dev of the pdf MC replicas
+        float pdfavg = std::accumulate(LHEPdfWeight.begin(), LHEPdfWeight.end(), 0.0) / size;
+        for (int ipdf = 0; ipdf < size; ipdf++) {
+            sumsquares = sumsquares + std::pow(LHEPdfWeight[ipdf] - pdfavg,2);
+        }
+        stddev = sqrt(sumsquares/(size-1));
+    }
+
+    v = {
+        std::min((float)13.0,(float)1.0+stddev),
+        std::max((float)-13.0,(float)1.0-stddev)
+    };
+
+    return v;
+}

TIMBER/Tools/Common.py

@@ -44,7 +44,7 @@ def CompileCpp(blockcode,library=False):
         ROOT.gInterpreter.AddIncludePath(os.environ["TIMBERPATH"])
 
     if not library:
-        if '\n' in blockcode: # must be multiline string
+        if '\n' in blockcode or ';' in blockcode: # must be multiline string
             ROOT.gInterpreter.Declare(blockcode)
         else: # must be file name to compile
             if ('TIMBER/Framework/' in blockcode) and (os.environ['TIMBERPATH'] not in blockcode):