diff --git a/TIMBER/Analyzer.py b/TIMBER/Analyzer.py index 5b07534..7dd4eac 100755 --- a/TIMBER/Analyzer.py +++ b/TIMBER/Analyzer.py @@ -4,6 +4,7 @@ """ +from TIMBER.CollectionOrganizer import CollectionOrganizer from TIMBER.Utilities.CollectionGen import BuildCollectionDict, GetKeyValForBranch, StructDef, StructObj from TIMBER.Tools.Common import GetHistBinningTuple, CompileCpp, ConcatCols, GetStandardFlags, ExecuteCmd from clang import cindex @@ -94,20 +95,20 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs", createAll super(analyzer, self).__init__() self.fileName = fileName - self.__eventsTreeName = eventsTreeName + self._eventsTreeName = eventsTreeName self.silent = False # Setup TChains for multiple or single file - self.__eventsChain = ROOT.TChain(self.__eventsTreeName) + self._eventsChain = ROOT.TChain(self._eventsTreeName) self.RunChain = ROOT.TChain(runTreeName) if isinstance(self.fileName,list): for f in self.fileName: - self.__addFile(f) + self._addFile(f) else: - self.__addFile(self.fileName) + self._addFile(self.fileName) # Make base RDataFrame - BaseDataFrame = ROOT.RDataFrame(self.__eventsChain) + BaseDataFrame = ROOT.RDataFrame(self._eventsChain) self.BaseNode = Node('base',BaseDataFrame) self.BaseNode.children = [] # protect against memory issue when running over multiple sets in one script self.AllNodes = [self.BaseNode] @@ -133,7 +134,7 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs", createAll # Get LHAID from LHEPdfWeights branch self.lhaid = "-1" if not self.isData: - pdfbranch = self.__eventsChain.GetBranch("LHEPdfWeight") + pdfbranch = self._eventsChain.GetBranch("LHEPdfWeight") if pdfbranch != None: branch_title = pdfbranch.GetTitle() if branch_title != '': @@ -151,10 +152,7 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs", createAll self.ActiveNode = self.BaseNode # Auto create collections - self.__collectionDict = BuildCollectionDict(BaseDataFrame) - self.__builtCollections = [] - if createAllCollections: - self.__createAllCollections(silent=True) + self._collectionOrg = CollectionOrganizer(BaseDataFrame) skipHeaders = [] if 'CMSSW_BASE' not in os.environ.keys(): @@ -164,7 +162,7 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs", createAll if f.split('/')[-1] in skipHeaders: continue CompileCpp('#include "%s"\n'%f) - def __addFile(self,f): + def _addFile(self,f): '''Add file to TChains being tracked. Args: @@ -173,13 +171,13 @@ def __addFile(self,f): if f.endswith(".root"): if 'root://' not in f and f.startswith('/store/'): f='root://cms-xrd-global.cern.ch/'+f - self.__eventsChain.Add(f) + self._eventsChain.Add(f) self.RunChain.Add(f) elif f.endswith(".txt"): txt_file = open(f,"r") for l in txt_file.readlines(): thisfile = l.strip() - self.__addFile(thisfile) + self._addFile(thisfile) else: raise Exception("File name extension not supported. Please provide a single or list of .root files or a .txt file with a line-separated list of .root files to chain together.") @@ -190,7 +188,7 @@ def Close(self): None ''' self.BaseNode.Close() - self.__eventsChain.Reset() + self._eventsChain.Reset() def __str__(self): '''Call with `print()` to print a nicely formatted description @@ -253,7 +251,7 @@ def GetCollectionNames(self): Returns: list(str): Collection names. ''' - return self.__collectionDict.keys() + return self._collectionOrg.collectionDict.keys() def SetActiveNode(self,node): '''Sets the active node. @@ -405,12 +403,12 @@ def Cut(self,name,cuts,node=None,nodetype=None): if isinstance(cuts,CutGroup): for c in cuts.keys(): cut = cuts[c] - newNode = self.__collectionDefCheck(cut, newNode) + newNode = self._collectionOrg.CollectionDefCheck(cut, newNode) newNode = newNode.Cut(c,cut,nodetype=nodetype,silent=self.silent) newNode.name = cuts.name+'__'+c self.TrackNode(newNode) elif isinstance(cuts,str): - newNode = self.__collectionDefCheck(cuts, newNode) + newNode = self._collectionOrg.CollectionDefCheck(cuts, newNode) newNode = newNode.Cut(name,cuts,nodetype=nodetype,silent=self.silent) self.TrackNode(newNode) else: @@ -441,29 +439,22 @@ def Define(self,name,variables,node=None,nodetype=None): if isinstance(variables,VarGroup): for v in variables.keys(): var = variables[v] - newNode = self.__collectionDefCheck(var, newNode) + newNode = self._collectionOrg.CollectionDefCheck(var, newNode) newNode = newNode.Define(v,var,nodetype=nodetype,silent=self.silent) + self._collectionOrg.AddBranch(v) newNode.name = variables.name+'__'+v self.TrackNode(newNode) # newNode.name = variables.name elif isinstance(variables,str): - newNode = self.__collectionDefCheck(variables, newNode) + newNode = self._collectionOrg.CollectionDefCheck(variables, newNode) newNode = newNode.Define(name,variables,nodetype=nodetype,silent=self.silent) + self._collectionOrg.AddBranch(name, str(newNode.DataFrame.GetColumnType(name))) self.TrackNode(newNode) else: raise TypeError("Second argument to Define method must be a string of a single var or of type VarGroup (which provides an OrderedDict).") return self.SetActiveNode(newNode) - def __collectionDefCheck(self, action_str, node): - newNode = node - for c in self.__collectionDict.keys(): - if (c+'s' not in self.__builtCollections) and re.search(r"\b" + re.escape(c+'s') + r"\b", action_str): - print ('MAKING %ss for %s'%(c,action_str)) - newNode = self.__createCollection(c,self.__collectionDict[c],silent=True,node=newNode) - self.__builtCollections.append(c+'s') - return self.SetActiveNode(newNode) - # Applies a bunch of action groups (cut or var) in one-shot in the order they are given def Apply(self,actionGroupList,node=None,trackEach=True): '''Applies a single CutGroup/VarGroup or an ordered list of Groups to the provided node or the #ActiveNode by default. @@ -550,12 +541,6 @@ def SubCollection(self,name,basecoll,condition,skip=[]): else: if condition != '': self.Define(replacementName,'%s[%s]'%(b,name+'_idx'),nodetype='SubCollDefine') else: self.Define(replacementName,b,nodetype='SubCollDefine') - - branches_to_track = [] - for v in collBranches: - v_with_type = GetKeyValForBranch(self.DataFrame, v)[1] - if v_with_type != '': branches_to_track.append(v_with_type) - self.__trackNewCollection(name,branches_to_track) return self.ActiveNode @@ -630,30 +615,6 @@ def MergeCollections(self,name,collectionNames): self.Define('n'+name,'+'.join(['n'+n for n in collectionNames]),nodetype='MergeDefine') - self.__trackNewCollection(name,[GetKeyValForBranch(self.DataFrame,collectionNames[0]+'_'+v)[1] for v in vars_to_make]) - - def __trackNewCollection(self,name,branches): - self.__collectionDict[name] = [] - for b in branches: - self.__collectionDict[name].append(b) - - def __createCollection(self,collection,attributes,silent=True,node=None): - init_silent = self.silent - self.silent = silent - if collection+'s' not in self.__builtCollections: - self.__builtCollections.append(collection+'s') - CompileCpp(StructDef(collection,attributes)) - newNode = self.Define(collection+'s', StructObj(collection,attributes),node) - else: - raise RuntimeError('Collections `%s` already built.'%(collection+'s')) - self.silent = init_silent - return self.SetActiveNode(newNode) - - def __createAllCollections(self,silent=True): - collDict = BuildCollectionDict(self.BaseNode.DataFrame) - for c in collDict.keys(): - self.__createCollection(c,collDict[c],silent) - def CommonVars(self,collections): '''Find the common variables between collections. @@ -751,7 +712,7 @@ def AddCorrections(self,correctionList,node=None): return self.SetActiveNode(newNode) - def __checkCorrections(self,node,correctionNames,dropList): + def _checkCorrections(self,node,correctionNames,dropList): '''Starting at the provided node, will scale up the tree, grabbing all corrections added along the way. This ensures corrections from other forks of the analyzer tree are not @@ -829,7 +790,7 @@ def MakeWeightCols(self,name='',node=None,correctionNames=None,dropList=[],corre if name != '': namemod = '_'+name else: namemod = '' - correctionsToApply = self.__checkCorrections(node,correctionNames,dropList) + correctionsToApply = self._checkCorrections(node,correctionNames,dropList) # Build nominal weight first (only "weight", no "uncert") weights = {'nominal':''} @@ -1049,14 +1010,14 @@ def CalibrateVars(self,varCalibDict,evalArgs,newCollectionName,variationsFlag=Tr ``` This will apply the JES and JER calibrations and their four variations (up,down pair for each) to FatJet_pt and FatJet_mass branches - and create a new collection called "CorrectedFatJets" which will be ordered by the new pt values. Note that if you want to correct a different + and create a new collection called "CalibratedFatJets" which will be ordered by the new pt values. Note that if you want to correct a different collection (ex. AK4 based Jet collection), you need a separate payload and separate call to CalibrateVars because only one collection can be generated at a time. Also note that in this example, `jes` and `jer` are initialized with the AK8PFPuppi jets in mind. So if you'd like to apply the JES or JER calibrations to AK4 jets, you would also need to define objects like `jesAK4` and `jerAK4`. The calibrations will always be calculated as a seperate column which stores a vector named `__vec` and ordered {nominal, up, down} where "up" and "down" are the absolute weights - (ie. not relative to "nominal"). If you'd just like the weights and do not want them applied to any variable, you can provide + (ie. "up" and "down" are not relative to "nominal"). If you'd just like the weights and do not want them applied to any variable, you can provide an empty dictionary (`{}`) for the varCalibDict argument. This method will set the new active node to the one with the new collection defined. @@ -1078,7 +1039,7 @@ def CalibrateVars(self,varCalibDict,evalArgs,newCollectionName,variationsFlag=Tr ''' if node == None: node = self.ActiveNode # Type checking and create calibration branches - newNode = self.__checkCalibrations(node,varCalibDict,evalArgs) + newNode = self._checkCalibrations(node,varCalibDict,evalArgs) # Create the product of weights new_columns = OrderedDict() @@ -1111,7 +1072,7 @@ def CalibrateVars(self,varCalibDict,evalArgs,newCollectionName,variationsFlag=Tr return self.SetActiveNode(newNode) - def __checkCalibrations(self,node,varCalibDict,evalArgs): + def _checkCalibrations(self,node,varCalibDict,evalArgs): newNode = node # Type checking if not isinstance(node,Node): raise TypeError('CalibrateVar() does not support argument of type %s for node. Please provide a Node.'%(type(node))) @@ -1825,7 +1786,6 @@ def __init__(self,name,script,constructor=[],mainFunc='eval',columnList=None,isC self._constructor = constructor self._objectName = self.name self._call = None - # self.__funcNames = self.__funcInfo.keys() if not isClone: if not self._mainFunc.endswith(mainFunc): diff --git a/TIMBER/CollectionOrganizer.py b/TIMBER/CollectionOrganizer.py new file mode 100644 index 0000000..45d6ec2 --- /dev/null +++ b/TIMBER/CollectionOrganizer.py @@ -0,0 +1,222 @@ +from TIMBER.Tools.Common import CompileCpp +import re + +class CollectionOrganizer: + '''Tracks the available collections, collection attributes, and solo branches + in the dataframe while processing. The initial set of branches will be read + on startup and any new branch will be added accordingly. Collection names + are deduced from the branch name by being the string before the first underscore + (if there is an underscore). + ''' + def __init__(self, rdf): + '''Constructor + + @param rdf (RDataFrame): RDataFrame from which to organize. + ''' + self.baseBranches = [str(b) for b in rdf.GetColumnNames()] + self.generateFromRDF(rdf) + self.builtCollections = [] + + def generateFromRDF(self, rdf): + '''Generate the collection from the RDataFrame. + + @param rdf (RDataFrame): RDataFrame from which to organize. + ''' + self.collectionDict = {} + self.otherBranches = {} + + for b in self.baseBranches: + self.AddBranch(b,rdf.GetColumnType(b)) + + def parsetype(self, t): + '''Deduce the type that TIMBER needs to see for the + collection structs. If t is an RVec, deduce the internal type + of the RVec. + + @param t (str): Type name from RDataFrame.GetColumnType() + + Returns: + str: TIMBER-friendly version of the type. + ''' + if not t.startswith('ROOT::VecOps::RVec<'): + collType = False + else: + collType = str(t).replace('ROOT::VecOps::RVec<','') + if collType.endswith('>'): + collType = collType[:-1] + collType += '&' + if 'Bool_t' in collType: + collType = collType.replace('Bool_t&','std::_Bit_reference') + + if collType == '&': + collType = '' + + return collType + + def AddCollection(self, c): + '''Add a collection to tracking. + + @param c (str): Collection name only. + ''' + if c not in self.collectionDict.keys(): + self.collectionDict[c] = {'alias': False} + + def GetCollectionAttributes(self, c): + '''Get all attributes of a collection. Example, for the 'Electron' + collection, will return a list of `['pt', 'eta', ...]`. + + @param c (str): Collection name. + + Returns: + list(str): List of attributes for the collection. + ''' + return [c for c in self.collectionDict[c] if c != 'alias'] + + def AddBranch(self, b, btype=''): + '''Add a branch to track. Will deduce if it is in a collection + in which case, the attribute will be added to the tracked collection. + + @param b (str): Branch name + @param btype (str, optional): Type of branch. Defaults to '' but should only be left + this way in rare cases. + ''' + collname = b.split('_')[0] + varname = '_'.join(b.split('_')[1:]) + typeStr = self.parsetype(btype) + + if typeStr == False or varname == '' or 'n'+collname not in self.baseBranches: + self.otherBranches[b] = { + 'type': typeStr, + 'alias': False + } + elif varname != '': + self.AddCollection(collname) + self.collectionDict[collname][varname] = { + 'type': typeStr, + 'alias': False + } + + def Alias(self, alias, name): + '''Add an alias for a solo branch, collection, or collection attribute. + + @param alias (str): Alias name. + @param name (str): Full branch name or a collection name. If an alias for a + collection attribute is desired, provide the full branch name (ie. _). + + Raises: + ValueError: Entries do not exist so an alias cannot be added. + ''' + # Name is either in otherBranches, is a collection name, or is a full name _ + if name in self.otherBranches.keys(): + self.otherBranches[name]['alias'] = alias + elif name in self.collectionDict.keys(): + self.collectionDict[name]['alias'] = alias + else: + collname = name.split('_')[0] + varname = '_'.join(name.split('_')[1:]) + if collname in self.collectionDict.keys(): + if varname in self.collectionDict[collname].keys(): + self.collectionDict[collname][varname]['alias'] = alias + else: + raise ValueError('Cannot add alias `%s` because attribute `%s` does not exist in collection `%s`'%(alias,varname,collname)) + else: + raise ValueError('Cannot add alias `%s` because collection `%s` does not exist'%(alias,collname)) + + def BuildCppCollection(self,collection,node,silent=True): + '''Build the collection as a struct in C++ so that it's accessible + to the RDataFrame loop. + + @param collection (str): Collection name. + @param node (Node): Node on which to act. + @param silent (bool, optional): Whether output should be silenced. Defaults to True. + + Raises: + RuntimeError: Collection already built. + + Returns: + Node: Manipulated node with the collection struct now defined. + ''' + newNode = node + attributes = [] + for aname in self.GetCollectionAttributes(collection): + attributes.append('%s %s'%(self.collectionDict[collection][aname]['type'], aname)) + + if collection+'s' not in self.builtCollections: + self.builtCollections.append(collection+'s') + CompileCpp(StructDef(collection,attributes)) + newNode = newNode.Define(collection+'s', StructObj(collection,attributes),silent=silent) + else: + raise RuntimeError('Collections `%s` already built.'%(collection+'s')) + + return newNode + + def CollectionDefCheck(self, action_str, node): + '''Checks if a collection C++ struct is needed in the action string. + If in the string but not defined, this function builds it. Does not + apply the action. + + @param action_str (str): Action being performed on the Node/RDataFrame. + @param node (Node): Node being acted on. + + Returns: + Node: Manipulated node with the C++ struct built (the action string is not applied though). + ''' + newNode = node + for c in self.collectionDict.keys(): + if re.search(r"\b" + re.escape(c+'s') + r"\b", action_str) and (c+'s' not in self.builtCollections): + print ('MAKING %ss for %s'%(c,action_str)) + newNode = self.BuildCppCollection(c,newNode,silent=True) + return newNode + +def StructDef(collectionName, varList): + '''Defines the struct in C++/Cling memory. + + @param collectionName (str): Name of the collection to define. + @param varList (str): List of attributes of the collection to include. + + Returns: + str: C++ string defining the struct. + ''' + out_str = ''' +struct {0}Struct {{ + {1} + {0}Struct({2}) : + {3} {{ + }}; +}}; + ''' + definitions = [] + ctor_args = [] + ctor_assign = [] + for i,v in enumerate(varList): + definitions.append('%s; \n'%v) + ctor_args.append('%s'%v) + ctor_assign.append('%s(%s)'%(v.split(' ')[-1], v.split(' ')[-1])) + + out_str = out_str.format(collectionName, '\t'.join(definitions), ','.join(ctor_args),','.join(ctor_assign)) + return out_str + +def StructObj(collectionName, varList): + '''Initializes an instance of the C++ struct for the collection in C++/Cling memory. + + @param collectionName (str): Name of the collection to define. + @param varList (str): List of attributes of the collection to include. + + Returns: + str: C++ string defining the struct instance. + ''' + out_str = ''' +std::vector<{0}Struct> {0}s; +{0}s.reserve(n{0}); +for (size_t i = 0; i < n{0}; i++) {{ + {0}s.emplace_back({1}); +}} +return {0}s; +''' + attr_assignment_str = '' + print (varList) + for i,v in enumerate(varList): + varname = v.split(' ')[-1] + attr_assignment_str += '{0}_{1}[i],'.format(collectionName, varname) + out_str = out_str.format(collectionName,attr_assignment_str[:-1]) + return out_str diff --git a/TIMBER/Framework/include/TopPt_weight.h b/TIMBER/Framework/include/TopPt_weight.h index ab3da25..5383a74 100644 --- a/TIMBER/Framework/include/TopPt_weight.h +++ b/TIMBER/Framework/include/TopPt_weight.h @@ -15,11 +15,9 @@ using namespace ROOT::VecOps; * * \f[ \sqrt{e^{\alpha - \beta \cdot p_{T}^{\textrm{Gen} t}} \cdot e^{\alpha - \beta \cdot p_{T}^{\textrm{Gen} \bar{t}} }} \f]. * - * where \f$\alpha = 0.0615\f$ and \f$\beta = 0.0005\f$. See the alpha() and beta() functions - * to calculate the weights with these parameters varied. - * - * WARNING: You MUST run corr() before alpha() and beta() since these functions - * recycle information derived from corr(). + * where \f$\alpha = 0.0615\f$ and \f$\beta = 0.0005\f$. See the eval() function + * to calculate the weight plus variations of the \f$\beta\f$ parameter. The \f$\alpha\f$ parameter + * is not varied since it would only represent a flat normalization change. * */ class TopPt_weight { @@ -28,48 +26,19 @@ class TopPt_weight { ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1); public: - TopPt_weight(){}; - ~TopPt_weight(){}; - /** - * @brief Calculate the top \f$p_T\f$ reweighting value for \f$t\bar{t}\f$ simulation - * based on doing gen particle matching. The weight is calculated as - * \f[ \sqrt{e^{\alpha - \beta \cdot p_{T}^{\textrm{Gen } t}} \cdot e^{\alpha - \beta \cdot p_{T}^{\textrm{Gen } \bar{t}} }} \f]. - * where \f$\alpha = 0.0615\f$ and \f$\beta = 0.0005\f$. See the alpha() and beta() functions - * to calculate the weights with these parameters varied. - * - * @param GenPart_pdgId NanoAOD branch - * @param GenPart_statusFlags NanoAOD branch - * @param GenPart_vects Vector of ROOT::Math::PtEtaPhiMVectors (create through hardware::TLvector) - * @param jet0 - * @param jet1 - * @return RVec Will only be length 1. Stored as vector to satisfy TIMBER Correction() requirements - */ - RVec corr(RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, - ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1); /** - * @brief Calculate variations of top \f$p_T\f$ weight by varying the \f$\alpha\f$ parameter. - * The amount of variation can be changed via the scale arguement which is a - * percent change on \f$\alpha\f$. The output is the weight calculated with the variation - * divided by the nominal value. When using MakeWeightCols(), the nominal will be multiplied - * by this variation to recover the total weight. + * @brief Construct a new TopPt_weight object. No arguments. * - * @param GenPart_pdgId NanoAOD branch - * @param GenPart_statusFlags NanoAOD branch - * @param GenPart_vects Vector of ROOT::Math::PtEtaPhiMVectors (create through hardware::TLvector) - * @param jet0 - * @param jet1 - * @param scale Percent variation on \f$\alpha\f$ parameter. - * @return RVec {up, down} variations of the top \f$p_T\f$ reweighting value divided by the nominal weight. */ - RVec alpha( - RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, - ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1, float scale = 0.5); + TopPt_weight(); + ~TopPt_weight(){}; /** - * @brief Calculate variations of the top \f$p_T\f$ weight by varying the \f$\beta\f$ parameter. + * @brief Calculate the top \f$p_T\f$ reweighting value for \f$t\bar{t}\f$ simulation + * based on doing gen particle matching. Additionally, calculate variations of the top + * \f$p_T\f$ weight by varying the \f$\beta\f$ parameter. * The amount of variation can be changed via the scale arguement which is a - * percent change on \f$\beta\f$. The output is the weight calculated with the variation - * divided by the nominal value. When using MakeWeightCols(), the nominal will be multiplied - * by this variation to recover the total weight. + * percent change on \f$\beta\f$. There is no corresponding function for \f$\alpha\f$ + * because the effect is only a flat normalization change. * * @param GenPart_pdgId NanoAOD branch * @param GenPart_statusFlags NanoAOD branch @@ -77,9 +46,9 @@ class TopPt_weight { * @param jet0 * @param jet1 * @param scale Percent variation on \f$\beta\f$ parameter. - * @return RVec {up, down} variations of the top \f$p_T\f$ reweighting value divided by the nominal weight. + * @return RVec {nom, up, down} variations of the top \f$p_T\f$ reweighting value (absolute). */ - RVec beta( + RVec eval( RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1, float scale = 0.5); }; diff --git a/TIMBER/Framework/src/TopPt_weight.cc b/TIMBER/Framework/src/TopPt_weight.cc index aba0be4..fda442d 100644 --- a/TIMBER/Framework/src/TopPt_weight.cc +++ b/TIMBER/Framework/src/TopPt_weight.cc @@ -1,5 +1,7 @@ #include "../include/TopPt_weight.h" +TopPt_weight::TopPt_weight(){}; + std::vector TopPt_weight::matchingGenPt( RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1){ @@ -24,55 +26,7 @@ std::vector TopPt_weight::matchingGenPt( return {genTPt,genTbarPt}; } -RVec TopPt_weight::corr( - RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, - ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1) { - - std::vector matched = matchingGenPt(GenPart_pdgId, GenPart_statusFlags, - GenPart_vects, jet0, jet1); - float genTPt = matched[0]; - float genTbarPt = matched[1]; - - float wTPt = 1.0; - if (genTPt > 0){ - wTPt = exp(0.0615 - 0.0005*genTPt); - } - - float wTbarPt = 1.0; - if (genTbarPt > 0){ - wTbarPt = exp(0.0615 - 0.0005*genTbarPt); - } - - return {sqrt(wTPt*wTbarPt)}; -} - -RVec TopPt_weight::alpha( - RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, - ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1, float scale){ - - std::vector matched = matchingGenPt(GenPart_pdgId, GenPart_statusFlags, - GenPart_vects, jet0, jet1); - float genTPt = matched[0]; - float genTbarPt = matched[1]; - - float wTPt_up = 1.0; - float wTPt_down = 1.0; - if (genTPt > 0){ - wTPt_up = exp((1+scale)*0.0615 - 0.0005*genTPt) / exp(0.0615 - 0.0005*genTPt); - wTPt_down = exp((1-scale)*0.0615 - 0.0005*genTPt) / exp(0.0615 - 0.0005*genTPt); - } - - float wTbarPt_up = 1.0; - float wTbarPt_down = 1.0; - if (genTbarPt > 0){ - wTbarPt_up = exp((1+scale)*0.0615 - 0.0005*genTbarPt) / exp(0.0615 - 0.0005*genTbarPt); - wTbarPt_down = exp((1-scale)*0.0615 - 0.0005*genTbarPt) / exp(0.0615 - 0.0005*genTbarPt); - } - - return {sqrt(wTPt_up*wTbarPt_up),sqrt(wTPt_down*wTbarPt_down)}; -} - -RVec TopPt_weight::beta( +RVec TopPt_weight::eval( RVec GenPart_pdgId, RVec GenPart_statusFlags, RVec GenPart_vects, ROOT::Math::PtEtaPhiMVector jet0, ROOT::Math::PtEtaPhiMVector jet1, float scale){ @@ -81,19 +35,23 @@ RVec TopPt_weight::beta( float genTPt = matched[0]; float genTbarPt = matched[1]; + float wTPt = 1.0; float wTPt_up = 1.0; float wTPt_down = 1.0; if (genTPt > 0){ - wTPt_up = exp(0.0615 - (1+scale)*0.0005*genTPt) / exp(0.0615 - 0.0005*genTPt); - wTPt_down = exp(0.0615 - (1-scale)*0.0005*genTPt) / exp(0.0615 - 0.0005*genTPt); + wTPt = exp(0.0615 - 0.0005*genTPt); + wTPt_up = exp((1+scale)*0.0615 - 0.0005*genTPt); + wTPt_down = exp((1-scale)*0.0615 - 0.0005*genTPt); } + float wTbarPt = 1.0; float wTbarPt_up = 1.0; float wTbarPt_down = 1.0; - if (genTbarPt > 0){ - wTbarPt_up = exp(0.0615 - (1+scale)*0.0005*genTbarPt) / exp(0.0615 - 0.0005*genTbarPt); - wTbarPt_down = exp(0.0615 - (1-scale)*0.0005*genTbarPt) / exp(0.0615 - 0.0005*genTbarPt); + if (genTbarPt > 0){ + wTbarPt = exp(0.0615 - 0.0005*genTbarPt); + wTbarPt_up = exp((1+scale)*0.0615 - 0.0005*genTbarPt); + wTbarPt_down = exp((1-scale)*0.0615 - 0.0005*genTbarPt); } - return {sqrt(wTPt_up*wTbarPt_up),sqrt(wTPt_down*wTbarPt_down)}; + return {sqrt(wTPt*wTbarPt),sqrt(wTPt_up*wTbarPt_up),sqrt(wTPt_down*wTbarPt_down)}; } \ No newline at end of file diff --git a/test/test_Analyzer.py b/test/test_Analyzer.py index 31c7139..e16a78b 100644 --- a/test/test_Analyzer.py +++ b/test/test_Analyzer.py @@ -97,6 +97,9 @@ def test_GetFlagString(self): assert self.a.GetFlagString(['HLT_IsoMu24','HLT_IsoMu24_eta2p1','NotReal']) == '((HLT_IsoMu24==1) && (HLT_IsoMu24_eta2p1==1))' pass + def test_CollectionStruct(self): + self.a.Cut('structCut','Jets[0].pt > 0') + def test_Groups(): a = analyzer('examples/GluGluToHToTauTau.root') # CompileCpp('TIMBER/Framework/include/common.h') # Compile (via gInterpreter) commonly used c++ code @@ -113,3 +116,7 @@ def test_Groups(): a.Apply([test_vg, test_cg]) rep = a.DataFrame.Report() rep.Print() + +def test_CollectionGroup(): + a = analyzer('examples/GluGluToHToTauTau.root') + assert ('Jet' in a.GetCollectionNames())