Merge pull request #68 from lcorcodilos/dev

Merge development for Beta 2.0

TIMBER/CollectionOrganizer.py

@@ -38,20 +38,23 @@ def _parsetype(self, t):
         Returns:
             str: TIMBER-friendly version of the type.
         '''
-        if not t.startswith('ROOT::VecOps::RVec<'):
-            collType = False
+        nRVecs = len(re.findall('ROOT::VecOps::RVec<',t))
+        if nRVecs == 0:
+            collType = t
+            isVect = False
         else:
-            collType = str(t).replace('ROOT::VecOps::RVec<','')
-            if collType.endswith('>'):
-                collType = collType[:-1]
-            collType += '&'
+            isVect = True
+            collType = t.strip()
+            collType = re.sub('ROOT::VecOps::RVec<','',collType,count=1)
+            collType = re.sub('>','',collType,count=1)
+            collType += ' &'
             if 'Bool_t' in collType:
-                collType = collType.replace('Bool_t&','std::_Bit_reference')
+                collType = collType.replace('Bool_t &','Bool_t&').replace('Bool_t&','std::_Bit_reference')
 
-        if collType == '&':
+        if collType == ' &':
             collType = ''
 
-        return collType
+        return collType, isVect
 
     def AddCollection(self, c):
         '''Add a collection to tracking.
@@ -90,17 +93,32 @@ def AddBranch(self, b, btype=''):
         '''
         collname = b.split('_')[0]
         varname = '_'.join(b.split('_')[1:])
-        typeStr = self._parsetype(btype)
+        typeStr, isVect = self._parsetype(btype)
 
         if typeStr == False or varname == '' or 'n'+collname not in self._baseBranches:
-            self._otherBranches[b] = {
-                'type': typeStr,
-                'alias': False
-            }
+            matches = [m for m in self._otherBranches.keys() if (m.startswith(collname) and '_'.join(m.split('_')[1:]) != '')]
+            if len(matches) == 0:
+                self._otherBranches[b] = {
+                    'type': typeStr,
+                    'isVect': isVect,
+                    'alias': False
+                }
+            else:
+                if varname != '':
+                    self.AddCollection(collname)
+                    self._collectionDict[collname][varname] = {
+                        'type': typeStr,
+                        'isVect': isVect,
+                        'alias': False
+                    }
+                    for match in matches:
+                        self._collectionDict[collname]['_'.join(match.split('_')[1:])] = self._otherBranches[match]
+                        del self._otherBranches[match]
         elif varname != '':
             self.AddCollection(collname)
             self._collectionDict[collname][varname] = {
                 'type': typeStr,
+                'isVect': isVect,
                 'alias': False
             }
 
@@ -147,7 +165,8 @@ def BuildCppCollection(self,collection,node,silent=True):
         newNode = node
         attributes = []
         for aname in self.GetCollectionAttributes(collection):
-            attributes.append('%s %s'%(self._collectionDict[collection][aname]['type'], aname))
+            if self._collectionDict[collection][aname]['isVect']:
+                attributes.append('%s %s'%(self._collectionDict[collection][aname]['type'], aname))
 
         if collection+'s' not in self._builtCollections:
             self._builtCollections.append(collection+'s')

TIMBER/Framework/include/EffLoader.h

@@ -1,3 +1,5 @@
+#ifndef _TIMBER_EFFLOADER
+#define _TIMBER_EFFLOADER
 #include <string>
 #include "TFile.h"
 #include "TEfficiency.h"
@@ -53,5 +55,5 @@ class EffLoader {
          * @return std::vector<float> {nominal value, up error+nominal, down error+nominal}
          */
         std::vector<float> eval(float xval, float yval = 0, float zval = 0);
-            
-};
+};
+#endif

TIMBER/Framework/include/JER_weight.h

@@ -34,7 +34,7 @@ class JER_weight {
          *  variations for that jet.
          */
         template <class Tjet, class TgenJet>
-        RVec< RVec<float> > eval(std::vector<Tjet> jets, std::vector<TgenJet> genJets, float fixedGridRhoFastjetAll){
+        RVec<RVec<float>> eval(std::vector<Tjet> jets, std::vector<TgenJet> genJets, float fixedGridRhoFastjetAll){
             RVec< RVec<float> > out (jets.size());
             for (size_t ijet = 0; ijet < jets.size(); ijet++) {
                 out[ijet] = (RVec<float>)_smearer.GetSmearValsPt(hardware::TLvector(jets[ijet]), hardware::TLvector(genJets), fixedGridRhoFastjetAll);

TIMBER/Framework/include/JES_weight.h

@@ -54,8 +54,8 @@ class JES_weight {
          *  variations for that jet.
          */
         template <class T>
-        RVec< RVec<float> > eval(std::vector<T> jets, float fixedGridRhoFastjetAll){
-            RVec< RVec<float> > out (jets.size());
+        RVec<RVec<float>> eval(std::vector<T> jets, float fixedGridRhoFastjetAll){
+            RVec<RVec<float>> out (jets.size());
 
             for (size_t ijet = 0; ijet < jets.size(); ijet++) {
                 RVec<float> ijet_out {1.0, 1.0, 1.0};

TIMBER/Framework/include/JMR_weight.h

@@ -48,8 +48,8 @@ class JMR_weight {
          *  variations for that jet.
          */
         template <class Tjet, class TgenJet>
-        RVec< RVec<float> > eval(std::vector<Tjet> jets, std::vector<TgenJet> genJets){
-            RVec< RVec<float> > out (jets.size());
+        RVec<RVec<float>> eval(std::vector<Tjet> jets, std::vector<TgenJet> genJets){
+            RVec<RVec<float>> out (jets.size());
 
             for (size_t ijet = 0; ijet < jets.size(); ijet++) {
                 RVec<LorentzV> genJetsV = hardware::TLvector(genJets);

TIMBER/Framework/include/JMS_weight.h

@@ -41,6 +41,6 @@ class JMS_weight {
          *  is the jet index and the inner vector is the nominal (0), up (1), and down (2
          *  variations for that jet.
          */
-        RVec< RVec<float>> eval(size_t nJets);
+        RVec<RVec<float>> eval(size_t nJets);
 };
 #endif

TIMBER/Framework/include/LumiFilter.h

@@ -27,7 +27,7 @@ class LumiFilter {
         LumiFilter(std::string filename);
         /**
          * @brief Construct a new Lumi Filter object for a given year.
-         * Supports formatting as `2017` or `17`. 2016 not needed currently.
+         * Supports formatting as (for example) `2017` or `17`.
          * 
          * @param year 
          */

TIMBER/Framework/include/Prefire_weight.h

@@ -21,47 +21,14 @@ class Prefire_weight {
         std::vector<std::pair<int, std::string>> _variations = {{0,"PrefireWeight"},{1,"PrefireWeight_Up"},{-1,"PrefireWeight_Down"}};
         int _variation;
 
-        template <class Photon, class Electron>
-        float EGvalue(int jetidx, std::vector<Photon> PhoColl, std::vector<Electron> EleColl) {
-            float phopf = 1.0;
-            float phopf_temp, elepf_temp;
-            std::vector<int> PhotonInJet;
-            Photon *pho;
-            for (size_t pid = 0; pid < PhoColl.size(); pid++) {
-                pho = &PhoColl[pid];
-                if (pho->jetIdx == jetidx) {
-                    if (PhotonInRange(pho->pt,pho->eta)) {
-                        phopf_temp = 1 - GetPrefireProbability(_photonmap, pho->eta, pho->pt, _photonPtRange.second);
-                        elepf_temp = 1.0;
-                        if (pho->electronIdx > -1) {
-                            Electron *matchEle = &EleColl[pho->electronIdx];
-                            if (PhotonInRange(matchEle->pt,matchEle->eta)) {
-                                elepf_temp = 1 - GetPrefireProbability(_photonmap,matchEle->eta,matchEle->pt,_photonPtRange.second);
-                            }
-                        }
-                        phopf *= std::min(phopf_temp,elepf_temp);
-                        PhotonInJet.push_back(pid);
-                    }
-                }
-            }
-            Electron *ele;
-            for (size_t eid = 0; eid < EleColl.size(); eid++){
-                ele = &EleColl[eid];
-                if ((ele->jetIdx == jetidx) && (Pythonic::InList(ele->photonIdx, PhotonInJet) == -1)) {
-                    if (PhotonInRange(ele->pt, ele->eta)) {
-                        phopf *= 1 - GetPrefireProbability(_photonmap, ele->eta, ele->pt, _photonPtRange.second);
-                    }
-                }
-            }
-            return phopf;
-        }
+        float EGvalue(int jetidx, RVec<float> Photon_pt, RVec<float> Photon_eta, RVec<int> Photon_jetIdx, RVec<int> Photon_electronIdx,
+                                  RVec<float> Electron_pt, RVec<float> Electron_eta, RVec<int> Electron_jetIdx, RVec<int> Electron_photonIdx);
 
         float GetPrefireProbability(TH2F* map, float eta, float pt, float maxpt);
         bool PhotonInRange(float pt, float eta);
         bool JetInRange(float pt, float eta);
         bool ObjInRange(float pt, float eta, std::pair<float,float> ptRange, std::pair<float,float> etaRange);
 
-
     public:
         /**
          * @brief Construct a new Prefire_weight object
@@ -74,38 +41,20 @@ class Prefire_weight {
         /**
          * @brief Calculate the value of the weight.
          * 
-         * @tparam Jet 
-         * @tparam Photon 
-         * @tparam Electron 
-         * @param Jets JetStruct from TIMBER
-         * @param Photons PhotonStruct from TIMBER
-         * @param Electrons ElectronStruct from TIMBER
+         * @param Jet_pt 
+         * @param Jet_eta 
+         * @param Photon_pt 
+         * @param Photon_eta 
+         * @param Photon_jetIdx 
+         * @param Photon_electronIdx 
+         * @param Electron_pt 
+         * @param Electron_eta 
+         * @param Electron_jetIdx 
+         * @param Electron_photonIdx 
          * @return ROOT::VecOps::RVec<float> 
          */
-        template <class Jet, class Photon, class Electron>
-        ROOT::VecOps::RVec<float> eval(std::vector<Jet> Jets, std::vector<Photon> Photons, std::vector<Electron> Electrons) {
-            Jet *jet;
-            float weight, jetpf, phopf;
-            std::string branchname;
-            ROOT::VecOps::RVec<float> out(3);
-            for (size_t i = 0; i<_variations.size(); i++){
-                _variation = _variations[i].first;
-                branchname = _variations[i].second;
-                weight = 1.0;
-                for (size_t ijet = 0; ijet<Jets.size(); ijet++) {
-                    jetpf = 1.0;
-                    jet = &Jets[ijet];
-
-                    if (JetInRange(jet->pt,jet->eta)) {
-                        jetpf *= 1 - GetPrefireProbability(_jetmap, jet->eta, jet->pt, _jetPtRange.second);
-                    }
-                    phopf = EGvalue(ijet, Photons, Electrons);
-                    weight *= std::min(jetpf, phopf);
-                }
-                weight *= EGvalue(-1, Photons, Electrons);
-                out[i] = weight;
-            }
-            return out;
-        }
+        ROOT::VecOps::RVec<float> eval(RVec<float> Jet_pt, RVec<float> Jet_eta,
+                                       RVec<float> Photon_pt, RVec<float> Photon_eta, RVec<int> Photon_jetIdx, RVec<int> Photon_electronIdx,
+                                       RVec<float> Electron_pt, RVec<float> Electron_eta, RVec<int> Electron_jetIdx, RVec<int> Electron_photonIdx);
 };
 #endif

TIMBER/Framework/include/common.h

@@ -32,24 +32,27 @@ namespace hardware {
      * just does the char + string formatting for you.
      * 
      * @param file Ex. TIMBER/data/example.root
+     * @param inTIMBER If file is in TIMBERPATH. Defaults to true.
      * @param option Defaults to "READ" for read-only.
      * @return TFile* 
      */
-    TFile *Open(std::string file, const char *option = "READ");
+    TFile *Open(std::string file, bool inTIMBER = true, const char *option = "READ");
     /**
      * @brief Generically open a histogram from a file into memory (closing
      * the file in the process).
      * 
      * @param filename 
      * @param histname 
+     * @param inTIMBER If file is in TIMBERPATH. Defaults to true.
      * @return TH1* 
      */
-    TH1 *LoadHist(std::string filename, std::string histname);
+    TH1 *LoadHist(std::string filename, std::string histname, bool inTIMBER = true);
     /**
      * @brief Hadamard product of two vectors (`v3[i] = v1[i]*v2[i]`)
      * 
      * @param v1 
      * @param v2 
+     * @param inTIMBER If file is in TIMBERPATH. Defaults to true.
      * @return RVec<float> 
      */
     RVec<float> HadamardProduct(RVec<float> v1, RVec<float> v2);
@@ -64,6 +67,14 @@ namespace hardware {
      * @return RVec<float> 
      */
     RVec<float> HadamardProduct(RVec<float> v1, RVec<RVec<float>> v2, int v2subindex);
+    /**
+     * @brief Hadamard product of a base vector and a list of N more vectors (`vout[i] = v1[i]*v2[i]*v3[i]...`).
+     * 
+     * @param v1 
+     * @param Multiv2 
+     * @return RVec<float> 
+     */
+    RVec<float> MultiHadamardProduct(RVec<float> v1, RVec<RVec<float>> Multiv2);
     /**
      * @brief Hadamard product of a base vector and a list of N more vectors (`vout[i] = v1[i]*v2[i][v2subindex]*v3[i][v2subindex]...`)
      * where v<N> has multiple sub-elements, only one of which should be accessed 

TIMBER/Framework/src/JMS_weight.cc

@@ -2,8 +2,8 @@
 
 JMS_weight::JMS_weight(int year) : _jmsVals(_jmsTable[year]) {};
 
-RVec< RVec<float>> JMS_weight::eval(size_t nFatJet){
-	RVec< RVec<float>> out;
+RVec<RVec<float>> JMS_weight::eval(size_t nFatJet){
+	RVec<RVec<float>> out;
 	out.reserve(nFatJet);
 	for (size_t i=0; i<nFatJet; i++){
 		out.push_back(this->_jmsVals);

TIMBER/Framework/src/LumiFilter.cc

@@ -5,10 +5,12 @@ LumiFilter::LumiFilter(std::string filename) {
 }
 
 LumiFilter::LumiFilter(int year){
-    if (year == 17 or year == 2017) {
-        boost::property_tree::read_json(std::getenv("TIMBERPATH")+std::string("TIMBER/data/LumiJSON/Cert_294927-306462_13TeV_EOY2017ReReco_Collisions17_JSON.txt"), ptree);
+    if (year == 16 or year == 2016) {
+        boost::property_tree::read_json(std::getenv("TIMBERPATH")+std::string("TIMBER/data/LumiJSON/Cert_271036-284044_13TeV_Legacy2016_Collisions16_JSON.txt"), ptree);
+    } else if (year == 17 or year == 2017) {
+        boost::property_tree::read_json(std::getenv("TIMBERPATH")+std::string("TIMBER/data/LumiJSON/Cert_294927-306462_13TeV_UL2017_Collisions17_GoldenJSON.txt"), ptree);
     } else if (year == 18 or year == 2018) {
-        boost::property_tree::read_json(std::getenv("TIMBERPATH")+std::string("TIMBER/data/LumiJSON/Cert_314472-325175_13TeV_17SeptEarlyReReco2018ABC_PromptEraD_Collisions18_JSON.txt"), ptree);
+        boost::property_tree::read_json(std::getenv("TIMBERPATH")+std::string("TIMBER/data/LumiJSON/Cert_314472-325175_13TeV_Legacy2018_Collisions18_JSON.txt"), ptree);
     } else {
         std::cout << "Year " << year << " not supported by LumiFilter" << std::endl;
         throw;

TIMBER/Framework/src/Pileup_weight.cc

@@ -17,7 +17,7 @@ void Pileup_weight::init(std::string filename_mc, std::string filename_data,
 
     if (filename_mc != "auto") {
         _autoPU = false;
-        _mcHist = hardware::LoadHist(filename_mc,histname_mc);
+        _mcHist = hardware::LoadHist(filename_mc,histname_mc,false);
     } else {
         _autoPU = true;
         gROOT->cd();

TIMBER/Framework/src/Prefire_weight.cc

@@ -55,4 +55,62 @@ bool Prefire_weight::ObjInRange(float pt, float eta, std::pair<float,float> ptRa
             out = true;
     }
     return out;
+}
+
+float Prefire_weight::EGvalue(int jetidx, RVec<float> Photon_pt, RVec<float> Photon_eta, RVec<int> Photon_jetIdx, RVec<int> Photon_electronIdx,
+                          RVec<float> Electron_pt, RVec<float> Electron_eta, RVec<int> Electron_jetIdx, RVec<int> Electron_photonIdx) {
+    float phopf = 1.0;
+    float phopf_temp, elepf_temp;
+    std::vector<int> PhotonInJet;
+    for (size_t pid = 0; pid < Photon_pt.size(); pid++) {
+        if (Photon_jetIdx[pid] == jetidx) {
+            if (PhotonInRange(Photon_pt[pid],Photon_eta[pid])) {
+                phopf_temp = 1 - GetPrefireProbability(_photonmap, Photon_eta[pid], Photon_pt[pid], _photonPtRange.second);
+                elepf_temp = 1.0;
+                if (Photon_electronIdx[pid] > -1) {
+                    int eid = Photon_electronIdx[pid];
+                    if (PhotonInRange(Electron_pt[eid],Electron_eta[eid])) {
+                        elepf_temp = 1 - GetPrefireProbability(_photonmap,Electron_eta[eid],Electron_pt[eid],_photonPtRange.second);
+                    }
+                }
+                phopf *= std::min(phopf_temp,elepf_temp);
+                PhotonInJet.push_back(pid);
+            }
+        }
+    }
+    for (size_t eid = 0; eid < Electron_pt.size(); eid++){
+        if ((Electron_jetIdx[eid] == jetidx) && (Pythonic::InList(Electron_photonIdx[eid], PhotonInJet) == -1)) {
+            if (PhotonInRange(Electron_pt[eid], Electron_eta[eid])) {
+                phopf *= 1 - GetPrefireProbability(_photonmap, Electron_eta[eid], Electron_pt[eid], _photonPtRange.second);
+            }
+        }
+    }
+    return phopf;
+}
+
+ROOT::VecOps::RVec<float> Prefire_weight::eval(RVec<float> Jet_pt, RVec<float> Jet_eta,
+                               RVec<float> Photon_pt, RVec<float> Photon_eta, RVec<int> Photon_jetIdx, RVec<int> Photon_electronIdx,
+                               RVec<float> Electron_pt, RVec<float> Electron_eta, RVec<int> Electron_jetIdx, RVec<int> Electron_photonIdx) {
+    float weight, jetpf, phopf;
+    std::string branchname;
+    ROOT::VecOps::RVec<float> out(3);
+    for (size_t i = 0; i<_variations.size(); i++){
+        _variation = _variations[i].first;
+        branchname = _variations[i].second;
+        weight = 1.0;
+        for (size_t ijet = 0; ijet<Jet_pt.size(); ijet++) {
+            jetpf = 1.0;
+
+            if (JetInRange(Jet_pt[ijet],Jet_eta[ijet])) {
+                jetpf *= 1 - GetPrefireProbability(_jetmap, Jet_eta[ijet], Jet_pt[ijet], _jetPtRange.second);
+            }
+            phopf = EGvalue(ijet, Photon_pt, Photon_eta, Photon_jetIdx, Photon_electronIdx,
+                                  Electron_pt, Electron_eta, Electron_jetIdx, Electron_photonIdx);
+            weight *= std::min(jetpf, phopf);
+        }
+        weight *= EGvalue(-1, Photon_pt, Photon_eta, Photon_jetIdx, Photon_electronIdx,
+                              Electron_pt, Electron_eta, Electron_jetIdx, Electron_photonIdx);
+        out[i] = weight;
+    }
+    return out;
 }