diff --git a/madanalysis/IOinterface/job_writer.py b/madanalysis/IOinterface/job_writer.py index 1f125067..d25e44ba 100644 --- a/madanalysis/IOinterface/job_writer.py +++ b/madanalysis/IOinterface/job_writer.py @@ -708,16 +708,20 @@ def WriteMakefiles(self, option="", **kwargs): options.has_process = True # @JACK enable usage of fastjet # If there are any root files, fastjet clusterer should not run with MA5_FASTJET_MODE Flag - options.ma5_fastjet_mode = ( - all([all([('root' not in x) for x in dataset.filenames]) for dataset in self.main.datasets]) - and self.main.archi_info.has_fastjet - and kwargs.get('ma5_fastjet_mode',True) - and self.main.fastsim.package == 'fastjet' - ) + options.ma5_fastjet_mode = self.main.archi_info.has_fastjet and self.main.archi_info.has_fjcontrib options.has_fastjet_lib = self.main.archi_info.has_fastjet options.has_fastjet_ma5lib = self.main.archi_info.has_fastjet # @JACK: fastjet_inc is required to be able to use the FJ files when FJ mode is on options.has_fastjet_inc = self.main.archi_info.has_fastjet + # @Jack: add substructure library + options.has_substructure = self.main.archi_info.has_fastjet and self.main.archi_info.has_fjcontrib + # @Jack: add HTT library + options.has_heptoptagger = ( + self.main.archi_info.has_fastjet and \ + self.main.archi_info.has_fjcontrib and \ + self.main.archi_info.has_heptoptagger + ) + options.has_root_inc = self.main.archi_info.has_root options.has_root_lib = self.main.archi_info.has_root #options.has_userpackage = True diff --git a/madanalysis/IOinterface/library_writer.py b/madanalysis/IOinterface/library_writer.py index c8e7bea1..93dc91ce 100644 --- a/madanalysis/IOinterface/library_writer.py +++ b/madanalysis/IOinterface/library_writer.py @@ -23,15 +23,10 @@ from __future__ import absolute_import -from madanalysis.selection.instance_name import InstanceName from madanalysis.IOinterface.folder_writer import FolderWriter -from madanalysis.IOinterface.job_writer import JobWriter -from string_tools import StringTools from shell_command import ShellCommand import logging -import shutil import os -import subprocess from six.moves import input @@ -126,6 +121,8 @@ def WriteMakefileForInterfaces(self,package): filename = self.path+"/SampleAnalyzer/Test/Makefile_zlib" elif package=='test_fastjet': filename = self.path+"/SampleAnalyzer/Test/Makefile_fastjet" + elif package=='test_htt': + filename = self.path+"/SampleAnalyzer/Test/Makefile_htt" elif package=='test_delphes': filename = self.path+"/SampleAnalyzer/Test/Makefile_delphes" elif package=='test_delphesMA5tune': @@ -149,6 +146,8 @@ def WriteMakefileForInterfaces(self,package): title='*zlib-interface* test' elif package=='test_fastjet': title='*fastjet-interface* test' + elif package=='test_htt': + title='*htt-interface* test' elif package=='test_delphes': title='*delphes-interface* test' elif package=='test_delphesMA5tune': @@ -173,6 +172,49 @@ def WriteMakefileForInterfaces(self,package): options.has_fastjet_inc=True # options.has_fastjet_lib=True toRemove.extend(['compilation_fastjet.log','linking_fastjet.log','cleanup_fastjet.log','mrproper_fastjet.log','../Bin/TestFastjet.log']) + elif package == "substructure": + options.has_commons=True + options.has_fastjet_inc=True + options.has_fastjet_lib=True + # @JACK: To be able to use fastjet in Ma5 data structure + options.ma5_fastjet_mode=True + options.has_fjcontrib = True + options.has_nsubjettiness = True + toRemove.extend( + ['compilation_substructure.log', + 'linking_substructure.log', + 'cleanup_substructure.log', + 'mrproper_substructure.log'] + ) + elif package == "HEPTopTagger": + options.has_commons=True + options.has_fastjet_inc=True + options.has_fastjet_lib=True + # @JACK: To be able to use fastjet in Ma5 data structure + options.ma5_fastjet_mode=True + options.has_nsubjettiness = True + options.has_substructure = True + toRemove.extend( + ['compilation_heptoptagger.log', + 'linking_heptoptagger.log', + 'cleanup_heptoptagger.log', + 'mrproper_heptoptagger.log'] + ) + elif package == "test_htt": + options.has_commons=True + options.has_fastjet_inc=True + # options.has_fastjet_lib=True + # @JACK: To be able to use fastjet in Ma5 data structure + options.ma5_fastjet_mode=True + options.has_fastjet_ma5lib = True + options.has_substructure = True + options.has_heptoptagger = True + toRemove.extend( + ['compilation_test_heptoptagger.log', + 'linking_test_heptoptagger.log', + 'cleanup_test_heptoptagger.log', + 'mrproper_test_heptoptagger.log'] + ) elif package=='configuration': toRemove.extend(['compilation.log','linking.log','cleanup.log','mrproper.log']) elif package=='commons': @@ -264,6 +306,7 @@ def WriteMakefileForInterfaces(self,package): options.has_fastjet_inc = self.main.archi_info.has_fastjet options.has_fastjet_lib = self.main.archi_info.has_fastjet options.ma5_fastjet_mode = self.main.archi_info.has_fastjet + options.has_substructure = self.main.archi_info.has_fjcontrib and self.main.archi_info.has_fastjet toRemove.extend(['compilation.log','linking.log','cleanup.log','mrproper.log']) elif package=='test_process': @@ -299,6 +342,9 @@ def WriteMakefileForInterfaces(self,package): elif package=='test_fastjet': cppfiles = ['Fastjet/*.cpp'] hfiles = ['Fastjet/*.h'] + elif package=='test_htt': + cppfiles = ['HEPTopTagger/*.cpp'] + hfiles = ['HEPTopTagger/*.h'] elif package=='test_delphes': cppfiles = ['Delphes/*.cpp'] hfiles = ['Delphes/*.h'] @@ -333,6 +379,10 @@ def WriteMakefileForInterfaces(self,package): isLibrary=False ProductName='TestFastjet' ProductPath='../Bin/' + elif package=='test_htt': + isLibrary=False + ProductName='TestHEPTopTagger' + ProductPath='../Bin/' elif package=='test_root': isLibrary=False ProductName='TestRoot' @@ -366,7 +416,7 @@ def Compile(self,ncores,package,folder): # log file name if package in ['process','commons','test','configuration']: logfile = folder+'/compilation.log' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet', "test_htt", 'test_root','test_delphes','test_delphesMA5tune']: logfile = folder+'/compilation_'+package[5:]+'.log' else: logfile = folder+'/compilation_'+package+'.log' @@ -374,7 +424,7 @@ def Compile(self,ncores,package,folder): # makefile if package in ['process','commons','test','configuration']: makefile = 'Makefile' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet', "test_htt", 'test_root','test_delphes','test_delphesMA5tune']: makefile = 'Makefile_'+package[5:] else: makefile = 'Makefile_'+package @@ -401,7 +451,7 @@ def Link(self,package,folder): # log file name if package in ['process','commons','test','configuration']: logfile = folder+'/linking.log' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet',"test_htt",'test_root','test_delphes','test_delphesMA5tune']: logfile = folder+'/linking_'+package[5:]+'.log' else: logfile = folder+'/linking_'+package+'.log' @@ -409,7 +459,7 @@ def Link(self,package,folder): # makefile if package in ['process','commons','test','configuration']: makefile = 'Makefile' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet',"test_htt",'test_root','test_delphes','test_delphesMA5tune']: makefile = 'Makefile_'+package[5:] else: makefile = 'Makefile_'+package @@ -433,7 +483,7 @@ def Clean(self,package,folder): # log file name if package in ['process','commons','configuration','test']: logfile = folder+'/cleanup.log' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet',"test_htt",'test_root','test_delphes','test_delphesMA5tune']: logfile = folder+'/cleanup_'+package[5:]+'.log' else: logfile = folder+'/cleanup_'+package+'.log' @@ -441,7 +491,7 @@ def Clean(self,package,folder): # makefile if package in ['process','commons','test','configuration']: makefile = 'Makefile' - elif package in ['test_process','test_commons','test_zlib','test_fastjet','test_root','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_fastjet',"test_htt",'test_root','test_delphes','test_delphesMA5tune']: makefile = 'Makefile_'+package[5:] else: makefile = 'Makefile_'+package @@ -465,7 +515,7 @@ def MrProper(self,package,folder): # log file name if package in ['process','commons','configuration']: logfile = folder+'/mrproper.log' - elif package in ['test_process','test_commons','test_zlib','test_root','test_fastjet','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_root','test_fastjet',"test_htt",'test_delphes','test_delphesMA5tune']: logfile = folder+'/mrproper_'+package[5:]+'.log' else: logfile = folder+'/mrproper_'+package+'.log' @@ -475,7 +525,7 @@ def MrProper(self,package,folder): # makefile if package in ['process','commons','test','configuration']: makefile = 'Makefile' - elif package in ['test_process','test_commons','test_zlib','test_root','test_fastjet','test_delphes','test_delphesMA5tune']: + elif package in ['test_process','test_commons','test_zlib','test_root','test_fastjet',"test_htt",'test_delphes','test_delphesMA5tune']: makefile = 'Makefile_'+package[5:] else: makefile = 'Makefile_'+package diff --git a/madanalysis/build/makefile_writer.py b/madanalysis/build/makefile_writer.py index 7179e673..d0dc9f66 100644 --- a/madanalysis/build/makefile_writer.py +++ b/madanalysis/build/makefile_writer.py @@ -168,6 +168,11 @@ def __init__(self): self.has_fastjet_tag = False # @JACK: cant use FASTJET_USE tag, it clashes with ROOT self.ma5_fastjet_mode = False + self.has_fjcontrib = False + self.has_substructure = False + self.has_heptoptagger = False + self.has_nsubjettiness = False + self.has_delphes_tag = False self.has_delphesMA5tune_tag = False self.has_root_tag = False @@ -363,7 +368,6 @@ def Makefile( if options.has_zlib_ma5lib: libs.extend(['-lzlib_for_ma5']) if options.has_zlib_lib: -# libs.extend(['-L'+archi_info.zlib_lib_path,'-lz']) libs.extend(['-lz']) file.write('LIBFLAGS += '+' '.join(libs)+'\n') @@ -386,14 +390,18 @@ def Makefile( if options.has_fastjet_ma5lib: libs.extend(['-lfastjet_for_ma5']) if options.has_fastjet_lib: -# libs.extend(['$(shell fastjet-config --libs --plugins)']) # --rpath=no)']) libs.extend(['$(shell $(MA5_BASE)/tools/SampleAnalyzer/ExternalSymLink/Bin/fastjet-config --libs --plugins)']) + file.write('LIBFLAGS += '+' '.join(libs)+'\n') + libs=[] + if options.has_fjcontrib: # Add fjcontrib libraries libs.extend(["-lRecursiveTools"]) # SoftDrop - libs.extend(["-lNsubjettiness"]) # Nsubjettiness libs.extend(["-lVariableR"]) # VariableR libs.extend(["-lEnergyCorrelator"]) # -lEnergyCorrelator - file.write('LIBFLAGS += '+' '.join(libs)+'\n') + if options.has_nsubjettiness: + libs.extend(["-lNsubjettiness"]) # Nsubjettiness + if options.has_nsubjettiness or options.has_fjcontrib: + file.write('LIBFLAGS += '+' '.join(libs)+'\n') # - delphes if options.has_delphes_ma5lib or options.has_delphes_lib: @@ -401,7 +409,6 @@ def Makefile( if options.has_delphes_ma5lib: libs.extend(['-ldelphes_for_ma5']) if options.has_delphes_lib: -# libs.extend(['-L'+archi_info.delphes_lib_paths[0],'-lDelphes']) libs.extend(['-lDelphes']) file.write('LIBFLAGS += '+' '.join(libs)+'\n') @@ -411,10 +418,17 @@ def Makefile( if options.has_delphesMA5tune_ma5lib: libs.extend(['-ldelphesMA5tune_for_ma5']) if options.has_delphesMA5tune_lib: -# libs.extend(['-L'+archi_info.delphesMA5tune_lib_paths[0],'-lDelphesMA5tune']) libs.extend(['-lDelphesMA5tune']) file.write('LIBFLAGS += '+' '.join(libs)+'\n') + # Substructure module + if options.has_substructure: + file.write('LIBFLAGS += -lsubstructure_for_ma5\n') + + # HEPTopTagger module + if options.has_heptoptagger: + file.write('LIBFLAGS += -lHEPTopTagger_for_ma5\n') + # - Commons if options.has_commons: libs=[] @@ -423,8 +437,6 @@ def Makefile( # - Root libs=[] - #libs.extend(['-L'+archi_info.root_lib_path, \ - # '-lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic -lEG']) # becareful: to not forget -lEG if options.has_root: libs.extend(['$(shell $(MA5_BASE)/tools/SampleAnalyzer/ExternalSymLink/Bin/root-config --libs)','-lEG']) @@ -448,6 +460,10 @@ def Makefile( libs.append('$(MA5_BASE)/tools/SampleAnalyzer/Lib/libdelphesMA5tune_for_ma5.so') if options.has_fastjet_ma5lib: libs.append('$(MA5_BASE)/tools/SampleAnalyzer/Lib/libfastjet_for_ma5.so') + if options.has_substructure: + libs.append('$(MA5_BASE)/tools/SampleAnalyzer/Lib/libsubstructure_for_ma5.so') + if options.has_heptoptagger: + libs.append('$(MA5_BASE)/tools/SampleAnalyzer/Lib/libHEPTopTagger_for_ma5.so') if len(libs)!=0: file.write('# Requirements to check before building\n') for ind in range(0,len(libs)): diff --git a/madanalysis/core/main.py b/madanalysis/core/main.py index ca92e0a5..bcbbfa0f 100644 --- a/madanalysis/core/main.py +++ b/madanalysis/core/main.py @@ -490,9 +490,11 @@ def BuildLibrary(self,forced=False): if not os.path.isfile(self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestSampleAnalyzer'): FirstUse=True - precompiler = LibraryWriter('lib',self) - if not precompiler.Run('TestSampleAnalyzer',[self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/Process/dummy_list.txt'],self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/',silent=True): - UpdateNeed=True + precompiler = LibraryWriter('lib', self) + if not precompiler.Run('TestSampleAnalyzer', + [self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/Process/dummy_list.txt'], + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/', silent=True): + UpdateNeed = True if not precompiler.CheckRun('TestSampleAnalyzer',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/',silent=True): UpdateNeed=True @@ -526,36 +528,80 @@ def BuildLibrary(self,forced=False): # |- [4] = folder # |- [5] = False=Library, True=Executable libraries = [] - libraries.append(['configuration','SampleAnalyzer configuration', 'configuration', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/PortabilityCheckup',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Configuration',True]) - libraries.append(['commons','SampleAnalyzer commons', 'commons', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libcommons_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Commons',False]) - libraries.append(['test_commons','SampleAnalyzer commons', 'test_commons', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestCommons',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['configuration', 'SampleAnalyzer configuration', 'configuration', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/PortabilityCheckup', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Configuration', True]) + libraries.append(['commons', 'SampleAnalyzer commons', 'commons', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libcommons_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Commons', False]) + libraries.append(['test_commons', 'SampleAnalyzer commons', 'test_commons', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestCommons', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # Zlib if self.archi_info.has_zlib: - libraries.append(['zlib', 'interface to zlib', 'zlib', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libzlib_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Interfaces',False]) - libraries.append(['test_zlib','interface to zlib', 'test_zlib', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestZlib',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['zlib', 'interface to zlib', 'zlib', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libzlib_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_zlib', 'interface to zlib', 'test_zlib', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestZlib', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # Fastjet if self.archi_info.has_fastjet: os.environ["FASTJET_FLAG"] = "-DMA5_FASTJET_MODE" - libraries.append(['FastJet', 'interface to FastJet', 'fastjet', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libfastjet_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Interfaces',False]) - libraries.append(['test_fastjet','interface to Fastjet', 'test_fastjet', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestFastjet',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['FastJet', 'interface to FastJet', 'fastjet', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libfastjet_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_fastjet', 'interface to Fastjet', 'test_fastjet', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestFastjet', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) + if self.archi_info.has_fjcontrib: + libraries.append(['substructure', 'interface to Jet Substructure module', 'substructure', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libsubstructure_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + if self.archi_info.has_heptoptagger: + libraries.append(['HEPTopTagger', 'interface to HEPTopTagger module', 'HEPTopTagger', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libHEPTopTagger_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_htt', 'interface to HEPTopTagger module', 'test_htt', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestHEPTopTagger', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) + else: + os.environ["FASTJET_FLAG"] = "" + # Delphes if self.archi_info.has_delphes: - libraries.append(['Delphes', 'interface to Delphes', 'delphes', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libdelphes_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Interfaces',False]) - libraries.append(['test_delphes','interface to Delphes', 'test_delphes', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestDelphes',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['Delphes', 'interface to Delphes', 'delphes', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libdelphes_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_delphes', 'interface to Delphes', 'test_delphes', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestDelphes', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # DelphesMA5tune if self.archi_info.has_delphesMA5tune: - libraries.append(['Delphes-MA5tune', 'interface to Delphes-MA5tune', 'delphesMA5tune', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libdelphesMA5tune_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Interfaces',False]) - libraries.append(['test_delphesMA5tune','interface to DelphesMA5tune', 'test_delphesMA5tune', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestDelphesMA5tune',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['Delphes-MA5tune', 'interface to Delphes-MA5tune', 'delphesMA5tune', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libdelphesMA5tune_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_delphesMA5tune', 'interface to DelphesMA5tune', 'test_delphesMA5tune', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestDelphesMA5tune', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # Root if self.archi_info.has_root: - libraries.append(['Root', 'interface to Root', 'root', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libroot_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Interfaces',False]) - libraries.append(['test_root','interface to Root', 'test_root', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestRoot',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['Root', 'interface to Root', 'root', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libroot_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Interfaces', False]) + libraries.append(['test_root', 'interface to Root', 'test_root', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestRoot', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # Process - libraries.append(['process', 'SampleAnalyzer core', 'process', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Lib/libprocess_for_ma5.so',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Process',False]) - libraries.append(['test_process','SampleAnalyzer core', 'test_process', self.archi_info.ma5dir+'/tools/SampleAnalyzer/Bin/TestSampleAnalyzer',self.archi_info.ma5dir+'/tools/SampleAnalyzer/Test/',True]) + libraries.append(['process', 'SampleAnalyzer core', 'process', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Lib/libprocess_for_ma5.so', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Process', False]) + libraries.append(['test_process', 'SampleAnalyzer core', 'test_process', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Bin/TestSampleAnalyzer', + self.archi_info.ma5dir + '/tools/SampleAnalyzer/Test/', True]) # Writing the Makefiles diff --git a/madanalysis/install/install_fastjet.py b/madanalysis/install/install_fastjet.py index 0f9a84fb..2a79a477 100644 --- a/madanalysis/install/install_fastjet.py +++ b/madanalysis/install/install_fastjet.py @@ -39,14 +39,7 @@ def __init__(self,main): self.downloaddir = self.main.session_info.downloaddir self.untardir = os.path.normpath(self.tmpdir + '/MA5_fastjet/') self.ncores = 1 -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.0.6.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.1.3.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.2.0.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.2.1.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.3.0.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.3.2.tar.gz"} - self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.3.3.tar.gz"} -# self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.4.0.tar.gz"} + self.files = {"fastjet.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fastjet-3.4.0.tar.gz"} def Detect(self): if not os.path.isdir(self.toolsdir): @@ -108,7 +101,9 @@ def Unpack(self): def Configure(self): # Input - theCommands=['./configure','--prefix='+self.installdir] + initial_env = os.environ.get("CXXFLAGS", "") + os.environ["CXXFLAGS"] = initial_env + " -std=c++11 -fPIC " + theCommands = ['./configure', '--prefix=' + self.installdir, "--enable-allplugins"] logname=os.path.normpath(self.installdir+'/configuration.log') # Execute logging.getLogger('MA5').debug('shell command: '+' '.join(theCommands)) @@ -116,6 +111,11 @@ def Configure(self): logname,\ self.tmpdir,\ silent=False) + # reset CXXFLAGS + if initial_env == "": + os.environ.pop("CXXFLAGS") + else: + os.environ["CXXFLAGS"] = initial_env # return result if not ok: logging.getLogger('MA5').error('impossible to configure the project. For more details, see the log file:') diff --git a/madanalysis/install/install_fastjetcontrib.py b/madanalysis/install/install_fastjetcontrib.py index 02fb7000..11af544c 100644 --- a/madanalysis/install/install_fastjetcontrib.py +++ b/madanalysis/install/install_fastjetcontrib.py @@ -40,13 +40,9 @@ def __init__(self,main): self.downloaddir = self.main.session_info.downloaddir self.untardir = os.path.normpath(self.tmpdir + '/MA5_fastjetcontrib/') self.ncores = 1 -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.012.tar.gz"} -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.017.tar.gz"} -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.024.tar.gz"} -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.027.tar.gz"} -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.039.tar.gz"} -# self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.042.tar.gz"} - self.files = {"fastjetcontrib.tar.gz" : "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.045.tar.gz"} + self.files = { + "fastjetcontrib.tar.gz": "http://madanalysis.irmp.ucl.ac.be/raw-attachment/wiki/WikiStart/fjcontrib-1.048.tar.gz" + } def GetNcores(self): self.ncores = InstallService.get_ncores(self.main.archi_info.ncores,\ @@ -85,7 +81,9 @@ def Unpack(self): def Configure(self): # Input - theCommands=['./configure','--fastjet-config='+self.bindir] + # TODO: figure out how to give `-std=c++11 -fPIC` together to CXXFLAGS + # using " or ' doesn't work on linux systems + theCommands = ['./configure', '--fastjet-config=' + self.bindir, 'CXXFLAGS=-fPIC'] logname=os.path.normpath(self.installdir+'/configuration_contrib.log') # Execute logging.getLogger('MA5').debug('shell command: '+' '.join(theCommands)) diff --git a/madanalysis/install/install_heptoptagger.py b/madanalysis/install/install_heptoptagger.py new file mode 100644 index 00000000..6ad855dd --- /dev/null +++ b/madanalysis/install/install_heptoptagger.py @@ -0,0 +1,161 @@ +################################################################################ +# +# Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +# The MadAnalysis development team, email: +# +# This file is part of MadAnalysis 5. +# Official website: +# +# MadAnalysis 5 is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# MadAnalysis 5 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with MadAnalysis 5. If not, see +# +################################################################################ + +from __future__ import absolute_import +from madanalysis.install.install_service import InstallService +from shell_command import ShellCommand +import os, json +import logging + + +class InstallHEPTopTagger: + def __init__(self, main): + self.main = main + self.installdir = os.path.normpath(self.main.archi_info.ma5dir + "/tools/HEPTopTagger/") + self.toolsdir = os.path.normpath(self.main.archi_info.ma5dir + "/tools") + self.tmpdir = self.main.session_info.tmpdir + self.downloaddir = self.main.session_info.downloaddir + self.untardir = os.path.normpath(self.tmpdir + "/MA5_heptoptagger/") + self.github_repo = "https://api.github.com/repos/MadAnalysis/HEPTopTagger/releases/latest" + self.meta = {} + + def Detect(self): + if not os.path.isdir(self.toolsdir): + logging.getLogger("MA5").debug("The folder '" + self.toolsdir + "' is not found") + return False + if not os.path.isdir(self.installdir): + logging.getLogger("MA5").debug("The folder " + self.installdir + "' is not found") + return False + return True + + def Remove(self, question=True): + from madanalysis.IOinterface.folder_writer import FolderWriter + return FolderWriter.RemoveDirectory(self.installdir, question) + + def CreatePackageFolder(self): + if not InstallService.create_tools_folder(self.toolsdir): + return False + if not InstallService.create_package_folder(self.toolsdir, "HEPTopTagger"): + return False + return True + + def CreateTmpFolder(self): + ok = InstallService.prepare_tmp(self.untardir, self.downloaddir) + if ok: + self.tmpdir = self.untardir + return ok + + def Download(self): + # Checking the link to GitHub + theCommands=['curl', '-s', self.github_repo] + logging.getLogger('MA5').debug('shell command: '+' '.join(theCommands)) + ok, out = ShellCommand.ExecuteWithLog( + theCommands, os.path.normpath(self.installdir + "/github_meta.log"), ".", silent=False + ) + if not ok: + return False + + self.meta = json.loads(out) + + # Downloading the package from GitHub + logging.getLogger('MA5').debug(f" -> HEPTopTagger {self.meta['tag_name']}") + logging.getLogger('MA5').debug(f" -> Published at {self.meta['published_at']}") + theCommands = [ + 'curl', '-s', "-L", "--create-dirs", "-o", + os.path.join(self.downloaddir, self.meta['tag_name'] + ".tar.gz"), + self.meta["tarball_url"] + ] + ok, out = ShellCommand.ExecuteWithLog( + theCommands, os.path.normpath(self.installdir + "/curl.log"), ".", silent=False + ) + logging.getLogger('MA5').debug(f"is download ok? {ok} :: {out}") + + return ok + + def Unpack(self): + # Logname + logname = os.path.normpath(self.installdir + "/unpack.log") + # Unpacking the tarball + logging.getLogger("MA5").debug(f"Unpack : {self.meta['tag_name']+'.tar.gz'}") + ok, packagedir = InstallService.untar( + logname, self.downloaddir, self.installdir, self.meta['tag_name']+'.tar.gz' + ) + logging.getLogger("MA5").debug(f"Unpack : {packagedir} is ok? {ok}") + if not ok: + return False + + # Checking where all files are + from glob import glob + content = glob(self.installdir + "/*") + if len(content) == 0: + return False + + found = False + main_folder = None + for file in content: + if os.path.isdir(file) and "MadAnalysis-HEPTopTagger" in file: + content = glob(os.path.join(file, "*")) + found = True + main_folder = file + break + + if not found: + return False + + # Copying all files where relevant + import shutil + for htt_file in content: + logging.getLogger("MA5").debug(f"copy {htt_file} to {self.installdir}") + shutil.copyfile(htt_file, os.path.join(self.installdir, os.path.basename(htt_file))) + + if main_folder is not None: + try: + shutil.rmtree(main_folder, ignore_errors=True) + except Exception as err: + logging.getLogger('MA5').debug(err) + + with open(os.path.join(self.installdir, "metadata.json"), "w") as meta: + json.dump(self.meta, meta, indent=4) + + # Ok: returning the good folder + self.tmpdir = packagedir + return True + + def Check(self): + # Check HTT files + for htt_file in ["HEPTopTagger.hh", "HEPTopTagger.cc"]: + if not os.path.isfile(os.path.join(self.installdir, htt_file)): + logging.getLogger("MA5").error(f"{htt_file} is missing.") + self.display_log() + return False + + return True + + def display_log(self): + logging.getLogger("MA5").error("More details can be found into the log files:") + logging.getLogger("MA5").error(" - " + os.path.normpath(self.installdir + "/github_meta.log")) + logging.getLogger("MA5").error(" - " + os.path.normpath(self.installdir + "/curl.log")) + logging.getLogger("MA5").error(" - " + os.path.normpath(self.installdir + "/unpack.log")) + + def NeedToRestart(self): + return True diff --git a/madanalysis/install/install_manager.py b/madanalysis/install/install_manager.py index 08110b37..270532d0 100644 --- a/madanalysis/install/install_manager.py +++ b/madanalysis/install/install_manager.py @@ -24,7 +24,6 @@ from __future__ import absolute_import import logging -import sys from string_tools import StringTools from chronometer import Chronometer @@ -49,6 +48,13 @@ def Execute(self, rawpackage): elif package=='fastjet-contrib': from madanalysis.install.install_fastjetcontrib import InstallFastjetContrib installer=InstallFastjetContrib(self.main) + elif package=='heptoptagger': + if self.main.archi_info.has_fastjet and self.main.archi_info.has_fjcontrib: + from madanalysis.install.install_heptoptagger import InstallHEPTopTagger + installer=InstallHEPTopTagger(self.main) + else: + self.logger.error("HEPTopTagger requires FastJet and the FJContrib libraries to be installed. Installation skipped") + return True elif package in ['delphes', 'delphesma5tune']: if self.main.archi_info.has_root: from madanalysis.install.install_delphes import InstallDelphes diff --git a/madanalysis/interpreter/cmd_install.py b/madanalysis/interpreter/cmd_install.py index e8e77439..ad8a73af 100644 --- a/madanalysis/interpreter/cmd_install.py +++ b/madanalysis/interpreter/cmd_install.py @@ -124,6 +124,8 @@ def UpdatePaths(): if installer.Execute('fastjet')==False: return False return installer.Execute('fastjet-contrib') + elif args[0]=='HEPTopTagger': + return installer.Execute('HEPTopTagger') elif args[0]=='gnuplot': return installer.Execute('gnuplot') elif args[0]=='matplotlib': @@ -223,7 +225,7 @@ def complete(self,text,args,begidx,endidx): else: output = ["samples","zlib","fastjet", "delphes", "delphesMA5tune",\ "gnuplot", "matplotlib", "root" , "numpy", "PAD", "PADForMA5tune",\ - "PADForSFS", "pyhf"] + "PADForSFS", "pyhf", "HEPTopTagger"] return self.finalize_complete(text,output) diff --git a/madanalysis/system/architecture_info.py b/madanalysis/system/architecture_info.py index d86dbc2e..32f07cbc 100644 --- a/madanalysis/system/architecture_info.py +++ b/madanalysis/system/architecture_info.py @@ -41,6 +41,8 @@ def __init__(self): # Is there optional package? self.has_root = False self.has_fastjet = False + self.has_fjcontrib = False + self.has_heptoptagger = False self.has_zlib = False self.has_delphes = False self.has_delphesMA5tune = False diff --git a/madanalysis/system/checkup.py b/madanalysis/system/checkup.py index 4f8ee9c8..462536f4 100644 --- a/madanalysis/system/checkup.py +++ b/madanalysis/system/checkup.py @@ -321,6 +321,10 @@ def CheckOptionalProcessingPackages(self): return False if not self.checker.Execute('fastjet'): return False + if not self.checker.Execute('fastjet-contrib'): + return False + if not self.checker.Execute('HEPTopTagger'): + return False if not self.checker.Execute('root'): return False diff --git a/madanalysis/system/detect_fastjetcontrib.py b/madanalysis/system/detect_fastjetcontrib.py new file mode 100644 index 00000000..4447c3db --- /dev/null +++ b/madanalysis/system/detect_fastjetcontrib.py @@ -0,0 +1,150 @@ +################################################################################ +# +# Copyright (C) 2012-2020 Jack Araz, Eric Conte & Benjamin Fuks +# The MadAnalysis development team, email: +# +# This file is part of MadAnalysis 5. +# Official website: +# +# MadAnalysis 5 is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# MadAnalysis 5 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with MadAnalysis 5. If not, see +# +################################################################################ + + +from __future__ import absolute_import +import logging, os +from shell_command import ShellCommand +from madanalysis.enumeration.detect_status_type import DetectStatusType + + +class DetectFastjetContrib: + def __init__(self, archi_info, user_info, session_info, debug): + self.archi_info = archi_info + self.user_info = user_info + self.session_info = session_info + self.debug = debug + self.name = "FastJet Contrib" + self.mandatory = False + self.force = False + self.lib_paths = [] + self.bin_file = "" + self.bin_path = "" + self.version = "" + self.logger = logging.getLogger("MA5") + self.required = ["RecursiveTools", "Nsubjettiness", "VariableR", "EnergyCorrelator"] + + def IsItVetoed(self): + if self.user_info.fastjet_veto: + self.logger.debug("user setting: veto on FastJet") + return True + else: + self.logger.debug("no user veto") + return False + + def detect_libs(self, fjconfig): + """ + Detect fastjet contrib libraries + """ + msg = "" + + theCommands = [fjconfig, "--libs"] + ok, out, err = ShellCommand.ExecuteWithCapture(theCommands, "./") + if not ok: + return False, "The fastjet-config executable does not work properly." + + libpath = [x[2:] for x in out.split() if x.startswith("-L")][0] + + not_found = [] + for libs in self.required: + if not any( + [f"lib{libs}{ext}" in os.listdir(libpath) for ext in [".a", ".dylib", ".so"]] + ): + not_found.append(libs) + if len(not_found) != 0: + return ( + DetectStatusType.UNFOUND, + "The following FJContrib libraries cannot be found: " + ", ".join(not_found), + ) + + # Ok + return DetectStatusType.FOUND, msg + + def ManualDetection(self): + msg = "" + + # User setting + if self.user_info.fastjet_bin_path == None: + return DetectStatusType.UNFOUND, msg + # Ok + return self.detect_libs(self.user_info.fastjet_bin_path + "/fastjet-config") + + def ToolsDetection(self): + msg = "" + + filename = os.path.normpath(self.archi_info.ma5dir + "/tools/fastjet/bin/fastjet-config") + self.logger.debug("Look for FastJet in the folder here:" + filename + " ...") + if os.path.isfile(filename): + self.logger.debug("-> found") + self.bin_file = filename + self.bin_path = os.path.dirname(self.bin_file) + else: + self.logger.debug("-> not found") + return DetectStatusType.UNFOUND, msg + + return self.detect_libs(self.archi_info.ma5dir + "/tools/fastjet/bin/fastjet-config") + + def AutoDetection(self): + msg = "" + + # Trying to call fastjet-config with which + result = ShellCommand.Which("fastjet-config", mute=True) + if len(result) == 0: + msg = "The FastJet package is not found." + return DetectStatusType.UNFOUND, msg + islink = os.path.islink(result[0]) + if not islink: + self.bin_file = os.path.normpath(result[0]) + else: + self.bin_file = os.path.normpath(os.path.realpath(result[0])) + self.bin_path = os.path.dirname(self.bin_file) + + # Debug mode + if self.debug: + self.logger.debug( + " which: " + str(result[0]) + " [is it a link? " + str(islink) + "]" + ) + if islink: + self.logger.debug(" -> " + os.path.realpath(result[0])) + + # Which all + if self.debug: + result = ShellCommand.Which("fastjet-config", all=True, mute=True) + if len(result) == 0: + msg = "The FastJet package is not found." + return DetectStatusType.UNFOUND, msg + self.logger.debug(" which-all: ") + for file in result: + self.logger.debug(" - " + str(file)) + + return self.detect_libs("fastjet-config") + + def ExtractInfo(self): + return True + + def SaveInfo(self): + # archi_info + self.archi_info.has_fjcontrib = True + + # Ok + return True diff --git a/madanalysis/system/detect_gpp.py b/madanalysis/system/detect_gpp.py index 7d662a80..a3c344a9 100644 --- a/madanalysis/system/detect_gpp.py +++ b/madanalysis/system/detect_gpp.py @@ -68,6 +68,23 @@ def AutoDetection(self): if self.debug: self.logger.debug(" which: " + str(result[0])) + # Check GCC version + try: + command = ["g++ -dumpversion > .gcc_version"] + result = ShellCommand.Execute(command, self.archi_info.ma5dir, shell=True) + with open(os.path.join(self.archi_info.ma5dir, ".gcc_version"), "r") as f: + result = f.read() + os.remove(os.path.join(self.archi_info.ma5dir, ".gcc_version")) + gcc_version = [int(x) for x in result.split(".")] + if (gcc_version[0] < 8 and not self.archi_info.isMac) or (gcc_version[0] < 9 and self.archi_info.isMac): + msg = "MadAnalysis 5 requires " + self.archi_info.isMac*"clang version 9 " + \ + (not self.archi_info.isMac)*"GCC version 8 " + "or above." + \ + f" Current version is " + ".".join([str(x) for x in gcc_version]) + self.logger.error(msg) + return DetectStatusType.UNFOUND, msg + except Exception as err: + self.logger.debug("Problem with compiler version detection:: " + str(err)) + # Check C++ version cpp = open(os.path.join(self.archi_info.ma5dir,"cxxtest.cc"), 'w') cpp.write("int main() { return 0; }") diff --git a/madanalysis/system/detect_heptoptagger.py b/madanalysis/system/detect_heptoptagger.py new file mode 100644 index 00000000..a95e205a --- /dev/null +++ b/madanalysis/system/detect_heptoptagger.py @@ -0,0 +1,102 @@ +################################################################################ +# +# Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +# The MadAnalysis development team, email: +# +# This file is part of MadAnalysis 5. +# Official website: +# +# MadAnalysis 5 is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# MadAnalysis 5 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with MadAnalysis 5. If not, see +# +################################################################################ + + +from __future__ import absolute_import +import logging, os +from madanalysis.enumeration.detect_status_type import DetectStatusType + + +class DetectHEPTopTagger: + def __init__(self, archi_info, user_info, session_info, debug): + self.archi_info = archi_info + self.user_info = user_info + self.session_info = session_info + self.debug = debug + self.name = "HEPTopTagger" + self.mandatory = False + self.force = False + self.lib_paths = [] + self.bin_file = "" + self.bin_path = "" + self.version = "" + self.logger = logging.getLogger("MA5") + + def IsItVetoed(self): + if self.user_info.fastjet_veto: + self.logger.debug("user setting: veto on FastJet") + return True + else: + self.logger.debug("no user veto") + return False + + def detect_files(self): + """ + Detect HEPTopTagger files + """ + msg = "" + + if not self.archi_info.has_fastjet: + logging.getLogger("MA5").debug(f" -> FastJet not found") + return DetectStatusType.UNFOUND, "FastJet not found." + if not self.archi_info.has_fjcontrib: + logging.getLogger("MA5").debug(f" -> FastJet contrib not found.") + return DetectStatusType.UNFOUND, "FastJet contrib not found." + + if not os.path.isdir(os.path.join(self.archi_info.ma5dir, "tools", "HEPTopTagger")): + logging.getLogger("MA5").debug( + f" -> The {os.path.join(self.archi_info.ma5dir, 'tools', 'HEPTopTagger')} folder does not exist." + ) + return ( + DetectStatusType.UNFOUND, + f"The {os.path.join(self.archi_info.ma5dir, 'tools', 'HEPTopTagger')} folder does not exist.", + ) + + # Check HTT files + for htt_file in ["HEPTopTagger.hh", "HEPTopTagger.cc"]: + if not os.path.isfile( + os.path.join(self.archi_info.ma5dir, "tools", "HEPTopTagger", htt_file) + ): + logging.getLogger("MA5").debug(f" -> {htt_file} is missing.") + return DetectStatusType.UNFOUND, f"{htt_file} is missing." + # Ok + return DetectStatusType.FOUND, msg + + def ManualDetection(self): + return self.detect_files() + + def ToolsDetection(self): + return self.detect_files() + + def AutoDetection(self): + return self.detect_files() + + def ExtractInfo(self): + return True + + def SaveInfo(self): + # archi_info + self.archi_info.has_heptoptagger = True + + # Ok + return True diff --git a/madanalysis/system/detect_manager.py b/madanalysis/system/detect_manager.py index 24a7ba99..93e80d2a 100644 --- a/madanalysis/system/detect_manager.py +++ b/madanalysis/system/detect_manager.py @@ -69,6 +69,9 @@ def Execute(self, rawpackage): elif package=='fastjet-contrib': from madanalysis.system.detect_fastjetcontrib import DetectFastjetContrib checker=DetectFastjetContrib(self.archi_info, self.user_info, self.session_info, self.debug) + elif package=='heptoptagger': + from madanalysis.system.detect_heptoptagger import DetectHEPTopTagger + checker=DetectHEPTopTagger(self.archi_info, self.user_info, self.session_info, self.debug) elif package=='delphes': from madanalysis.system.detect_delphes import DetectDelphes checker=DetectDelphes(self.archi_info, self.user_info, self.session_info, self.debug) diff --git a/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.cpp b/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.cpp new file mode 100644 index 00000000..8cb7c7b8 --- /dev/null +++ b/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.cpp @@ -0,0 +1,258 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2020 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// Fastjet headers +#ifdef MA5_FASTJET_MODE +#include "fastjet/PseudoJet.hh" +#endif + +// SampleAnalyzer headers +#include "SampleAnalyzer/Commons/DataFormat/RecEventFormat.h" + +namespace MA5 { + /// Clearing all information + void RecEventFormat::Reset() + { + PrimaryJetID_ = "Ma5Jet"; + photons_.clear(); + electrons_.clear(); + muons_.clear(); + taus_.clear(); + jetcollection_.clear(); + fatjets_.clear(); + emptyjet_.clear(); +#ifdef MA5_FASTJET_MODE + // pseudojet : only available when fastjet is in use (code efficiency) + input_hadrons_.clear(); +#endif + towers_ok_=false; + towers_.clear(); + tracks_ok_=false; + tracks_.clear(); + vertices_ok_=false; + vertices_.clear(); + EFlowTracks_ok_=false; + EFlowTracks_.clear(); + EFlowPhotons_ok_=false; + EFlowPhotons_.clear(); + EFlowNeutralHadrons_ok_=false; + EFlowNeutralHadrons_.clear(); + genjets_.clear(); + MET_.Reset(); + MHT_.Reset(); + TET_ = 0.; + THT_ = 0.; + Meff_ = 0.; + MCHadronicTaus_.clear(); + MCMuonicTaus_.clear(); + MCElectronicTaus_.clear(); + MCBquarks_.clear(); + MCCquarks_.clear(); + } + + //======================// + //===== Jet Tools ======// + //======================// + + // Remove an element from jet collection + void RecEventFormat::Remove_Collection(std::string id) + { + if (hasJetID(id)) jetcollection_.erase(id); + else ERROR << "Remove_Collection:: '" << id << "' does not exist." << endmsg; + } + + /// Change Jet ID + void RecEventFormat::ChangeJetID(std::string previous_id, std::string new_id) + { + if (!hasJetID(new_id) && hasJetID(previous_id)) + { + auto it = jetcollection_.find(previous_id); + std::swap(jetcollection_[new_id],it->second); + jetcollection_.erase(it); + } + else + { + if (hasJetID(new_id)) + ERROR << "ChangeJetID:: '" << new_id << "' already exists." << endmsg; + if (!hasJetID(previous_id)) + ERROR << "ChangeJetID:: '" << previous_id << "' does not exist." << endmsg; + } + } + + // Get the list of jet collection IDs + const std::vector RecEventFormat::GetJetIDs() const + { + std::vector keys; + keys.reserve(jetcollection_.size()); + for (auto &key: jetcollection_) + keys.emplace_back(key.first); + return keys; + } + + + // Add a new hadron to be clustered. (for code efficiency) + void RecEventFormat::AddHadron(MCParticleFormat& v, MAuint32& idx) + { +#ifdef MA5_FASTJET_MODE + fastjet::PseudoJet input; + input.reset(v.px(), v.py(), v.pz(), v.e()); + input.set_user_index(idx); + input_hadrons_.push_back(input); +#endif + } + +#ifdef MA5_FASTJET_MODE + // Get hadrons to cluster (code efficiency) + std::vector& RecEventFormat::cluster_inputs() {return input_hadrons_;} +#endif + + //======================// + //=== Get New Object ===// + //======================// + + /// Giving a new photon entry + RecPhotonFormat* RecEventFormat::GetNewPhoton() + { + photons_.push_back(RecPhotonFormat()); + return &photons_.back(); + } + + /// Giving a new electron entry + RecLeptonFormat* RecEventFormat::GetNewElectron() + { + electrons_.push_back(RecLeptonFormat()); + (&electrons_.back())->setElectronId(); + return &electrons_.back(); + } + + /// Giving a new muon entry + RecLeptonFormat* RecEventFormat::GetNewMuon() + { + muons_.push_back(RecLeptonFormat()); + (&muons_.back())->setMuonId(); + return &muons_.back(); + } + + /// Giving a new tower entry + RecTowerFormat* RecEventFormat::GetNewTower() + { + towers_.push_back(RecTowerFormat()); + return &towers_.back(); + } + + /// Giving a new EFlowTrack entry + RecTrackFormat* RecEventFormat::GetNewEFlowTrack() + { + EFlowTracks_.push_back(RecTrackFormat()); + return &EFlowTracks_.back(); + } + + /// Giving a new EFlowPhoton entry + RecParticleFormat* RecEventFormat::GetNewEFlowPhoton() + { + EFlowPhotons_.push_back(RecParticleFormat()); + return &EFlowPhotons_.back(); + } + + /// Giving a new EFlowNeutralHadron entry + RecParticleFormat* RecEventFormat::GetNewEFlowNeutralHadron() + { + EFlowNeutralHadrons_.push_back(RecParticleFormat()); + return &EFlowNeutralHadrons_.back(); + } + + /// Giving a new tau entry + RecTauFormat* RecEventFormat::GetNewTau() + { + taus_.push_back(RecTauFormat()); + return &taus_.back(); + } + + /// Giving a new primary jet entry + RecJetFormat* RecEventFormat::GetNewJet() + { + std::pair< std::map >::iterator, bool> new_jet; + new_jet = jetcollection_.insert(std::make_pair(PrimaryJetID_,std::vector() )); + new_jet.first->second.push_back(RecJetFormat()); + return &(new_jet.first->second.back()); + } + + /// Giving a new primary jet entry + void RecEventFormat::CreateEmptyJetAccesor() + { + std::pair< std::map >::iterator, bool> new_jet; + new_jet = jetcollection_.insert(std::make_pair(PrimaryJetID_,std::vector() )); + } + + // Get a new jet entry with specific ID + RecJetFormat* RecEventFormat::GetNewJet(std::string id) + { + std::pair< std::map >::iterator,bool> new_jet; + new_jet = jetcollection_.insert(std::make_pair(id,std::vector() )); + new_jet.first->second.push_back(RecJetFormat()); + return &(new_jet.first->second.back()); + } + + // Create an empty jet accessor with specific id + void RecEventFormat::CreateEmptyJetAccesor(std::string id) + { + std::pair< std::map >::iterator,bool> new_jet; + new_jet = jetcollection_.insert(std::make_pair(id,std::vector() )); + } + + /// Giving a new fat jet entry + RecJetFormat* RecEventFormat::GetNewFatJet() + { + fatjets_.push_back(RecJetFormat()); + return &fatjets_.back(); + } + + /// Giving a new gen jet entry + RecJetFormat* RecEventFormat::GetNewGenJet() + { + genjets_.push_back(RecJetFormat()); + return &genjets_.back(); + } + + /// Giving a new track entry + RecTrackFormat* RecEventFormat::GetNewTrack() + { + tracks_.push_back(RecTrackFormat()); + return &tracks_.back(); + } + + /// Giving a new vertex entry + RecVertexFormat* RecEventFormat::GetNewVertex() + { + vertices_.push_back(RecVertexFormat()); + return &vertices_.back(); + } + + /// Giving a pointer to the Missing Transverse Energy + RecParticleFormat* RecEventFormat::GetNewMet() { return &MET_; } + + /// Giving a pointer to the Missing Transverse Energy + RecParticleFormat* RecEventFormat::GetNewMht() { return &MHT_; } + + +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.h b/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.h index 9b463862..2b295631 100644 --- a/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.h +++ b/tools/SampleAnalyzer/Commons/DataFormat/RecEventFormat.h @@ -25,11 +25,6 @@ #ifndef RecEventFormat_h #define RecEventFormat_h -// Fastjet headers -#ifdef MA5_FASTJET_MODE -#include "fastjet/PseudoJet.hh" -#endif - // STL headers #include #include @@ -46,6 +41,13 @@ #include "SampleAnalyzer/Commons/DataFormat/RecVertexFormat.h" #include "SampleAnalyzer/Commons/Service/LogService.h" +#ifdef MA5_FASTJET_MODE + namespace fastjet { + class PseudoJet; + } +#endif + + namespace MA5 { @@ -112,6 +114,9 @@ namespace MA5 std::vector input_hadrons_; #endif + /// Create an empty jet + std::vector emptyjet_; + /// Collection of generated jets std::vector genjets_; @@ -205,8 +210,7 @@ namespace MA5 return jetcollection_.at(PrimaryJetID_); else { - static const std::vector empty; - return empty; + return emptyjet_; } } @@ -220,8 +224,7 @@ namespace MA5 return jetcollection_.at(id); else { - static const std::vector empty; - return empty; + return emptyjet_; } } @@ -294,8 +297,7 @@ namespace MA5 return jetcollection_.at(PrimaryJetID_); else { - static std::vector empty; - return empty; + return emptyjet_; } } @@ -309,8 +311,7 @@ namespace MA5 return jetcollection_.at(id); else { - static std::vector empty; - return empty; + return emptyjet_; } } @@ -365,91 +366,29 @@ namespace MA5 {return MCCquarks_;} /// Clearing all information - void Reset() - { - PrimaryJetID_ = "Ma5Jet"; - photons_.clear(); - electrons_.clear(); - muons_.clear(); - taus_.clear(); - jetcollection_.clear(); - fatjets_.clear(); -#ifdef MA5_FASTJET_MODE - // pseudojet : only available when fastjet is in use (code efficiency) - input_hadrons_.clear(); -#endif - towers_ok_=false; - towers_.clear(); - tracks_ok_=false; - tracks_.clear(); - vertices_ok_=false; - vertices_.clear(); - EFlowTracks_ok_=false; - EFlowTracks_.clear(); - EFlowPhotons_ok_=false; - EFlowPhotons_.clear(); - EFlowNeutralHadrons_ok_=false; - EFlowNeutralHadrons_.clear(); - genjets_.clear(); - MET_.Reset(); - MHT_.Reset(); - TET_ = 0.; - THT_ = 0.; - Meff_ = 0.; - MCHadronicTaus_.clear(); - MCMuonicTaus_.clear(); - MCElectronicTaus_.clear(); - MCBquarks_.clear(); - MCCquarks_.clear(); - } + void Reset(); // Initialize PrimaryJetID void SetPrimaryJetID(std::string v) {PrimaryJetID_ = v;} // Remove an element from jet collection - void Remove_Collection(std::string id) - { - if (hasJetID(id)) jetcollection_.erase(id); - else ERROR << "Remove_Collection:: '" << id << "' does not exist." << endmsg; - } + void Remove_Collection(std::string id); /// Change Jet ID - void ChangeJetID(std::string previous_id, std::string new_id) - { - if (!hasJetID(new_id) && hasJetID(previous_id)) - { - auto it = jetcollection_.find(previous_id); - std::swap(jetcollection_[new_id],it->second); - jetcollection_.erase(it); - } - else - { - if (hasJetID(new_id)) - ERROR << "ChangeJetID:: '" << new_id << "' already exists." << endmsg; - if (!hasJetID(previous_id)) - ERROR << "ChangeJetID:: '" << previous_id << "' does not exist." << endmsg; - } - } + void ChangeJetID(std::string previous_id, std::string new_id); // Get the list of jet collection IDs - const std::vector GetJetIDs() const - { - std::vector keys; - keys.reserve(jetcollection_.size()); - for (auto &key: jetcollection_) - keys.emplace_back(key.first); - return keys; - } + const std::vector GetJetIDs() const; // Check if collection ID exists MAbool hasJetID(std::string id) { return (jetcollection_.find(id) != jetcollection_.end()); } -#ifdef MA5_FASTJET_MODE // Add a new hadron to be clustered. (for code efficiency) - void AddHadron(fastjet::PseudoJet v) {input_hadrons_.push_back(v);} + void AddHadron(MCParticleFormat& v, MAuint32& idx); +#ifdef MA5_FASTJET_MODE // Get hadrons to cluster (code efficiency) - std::vector& cluster_inputs() {return input_hadrons_;} + std::vector& cluster_inputs(); #endif /// Displaying data member values. @@ -474,130 +413,58 @@ namespace MA5 } /// Giving a new photon entry - RecPhotonFormat* GetNewPhoton() - { - photons_.push_back(RecPhotonFormat()); - return &photons_.back(); - } + RecPhotonFormat* GetNewPhoton(); /// Giving a new electron entry - RecLeptonFormat* GetNewElectron() - { - electrons_.push_back(RecLeptonFormat()); - (&electrons_.back())->setElectronId(); - return &electrons_.back(); - } + RecLeptonFormat* GetNewElectron(); /// Giving a new muon entry - RecLeptonFormat* GetNewMuon() - { - muons_.push_back(RecLeptonFormat()); - (&muons_.back())->setMuonId(); - return &muons_.back(); - } + RecLeptonFormat* GetNewMuon(); /// Giving a new tower entry - RecTowerFormat* GetNewTower() - { - towers_.push_back(RecTowerFormat()); - return &towers_.back(); - } + RecTowerFormat* GetNewTower(); /// Giving a new EFlowTrack entry - RecTrackFormat* GetNewEFlowTrack() - { - EFlowTracks_.push_back(RecTrackFormat()); - return &EFlowTracks_.back(); - } + RecTrackFormat* GetNewEFlowTrack(); /// Giving a new EFlowPhoton entry - RecParticleFormat* GetNewEFlowPhoton() - { - EFlowPhotons_.push_back(RecParticleFormat()); - return &EFlowPhotons_.back(); - } + RecParticleFormat* GetNewEFlowPhoton(); /// Giving a new EFlowNeutralHadron entry - RecParticleFormat* GetNewEFlowNeutralHadron() - { - EFlowNeutralHadrons_.push_back(RecParticleFormat()); - return &EFlowNeutralHadrons_.back(); - } + RecParticleFormat* GetNewEFlowNeutralHadron(); /// Giving a new tau entry - RecTauFormat* GetNewTau() - { - taus_.push_back(RecTauFormat()); - return &taus_.back(); - } + RecTauFormat* GetNewTau(); /// Giving a new primary jet entry - RecJetFormat* GetNewJet() - { - std::pair< std::map >::iterator, bool> new_jet; - new_jet = jetcollection_.insert(std::make_pair(PrimaryJetID_,std::vector() )); - new_jet.first->second.push_back(RecJetFormat()); - return &(new_jet.first->second.back()); - } + RecJetFormat* GetNewJet(); /// Giving a new primary jet entry - void CreateEmptyJetAccesor() - { - std::pair< std::map >::iterator, bool> new_jet; - new_jet = jetcollection_.insert(std::make_pair(PrimaryJetID_,std::vector() )); - } + void CreateEmptyJetAccesor(); // Get a new jet entry with specific ID - RecJetFormat* GetNewJet(std::string id) - { - std::pair< std::map >::iterator,bool> new_jet; - new_jet = jetcollection_.insert(std::make_pair(id,std::vector() )); - new_jet.first->second.push_back(RecJetFormat()); - return &(new_jet.first->second.back()); - } + RecJetFormat* GetNewJet(std::string id); // Create an empty jet accessor with specific id - void CreateEmptyJetAccesor(std::string id) - { - std::pair< std::map >::iterator,bool> new_jet; - new_jet = jetcollection_.insert(std::make_pair(id,std::vector() )); - } + void CreateEmptyJetAccesor(std::string id); /// Giving a new fat jet entry - RecJetFormat* GetNewFatJet() - { - fatjets_.push_back(RecJetFormat()); - return &fatjets_.back(); - } + RecJetFormat* GetNewFatJet(); /// Giving a new gen jet entry - RecJetFormat* GetNewGenJet() - { - genjets_.push_back(RecJetFormat()); - return &genjets_.back(); - } + RecJetFormat* GetNewGenJet(); /// Giving a new track entry - RecTrackFormat* GetNewTrack() - { - tracks_.push_back(RecTrackFormat()); - return &tracks_.back(); - } + RecTrackFormat* GetNewTrack(); /// Giving a new vertex entry - RecVertexFormat* GetNewVertex() - { - vertices_.push_back(RecVertexFormat()); - return &vertices_.back(); - } + RecVertexFormat* GetNewVertex(); /// Giving a pointer to the Missing Transverse Energy - RecParticleFormat* GetNewMet() - { return &MET_; } + RecParticleFormat* GetNewMet(); /// Giving a pointer to the Missing Transverse Energy - RecParticleFormat* GetNewMht() - { return &MHT_; } + RecParticleFormat* GetNewMht(); }; diff --git a/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.cpp b/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.cpp new file mode 100644 index 00000000..62d97519 --- /dev/null +++ b/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.cpp @@ -0,0 +1,80 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" + +#ifdef MA5_FASTJET_MODE + +// FastJet headers +#include "fastjet/PseudoJet.hh" + +namespace MA5 { + // return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) + // that would be obtained when running the algorithm with the given dcut. + std::vector RecJetFormat::exclusive_subjets(MAfloat32 dcut) const + { + std::vector output_jets; + try { + if (!RecJetFormat::has_exclusive_subjets()) + throw EXCEPTION_ERROR("This jet structure does not contain exclusive_subjets", + "Exclusive subjets only exist for jets clustered through an exclusive algorithm.", 0); + } catch (const std::exception& err) { + MANAGE_EXCEPTION(err); + return output_jets; + } + std::vector current_jets = fastjet::sorted_by_pt(pseudojet_.exclusive_subjets(dcut)); + output_jets.reserve(current_jets.size()); + for (auto &jet: current_jets) + { + RecJetFormat * NewJet = new RecJetFormat(jet); + output_jets.push_back(NewJet); + } + return output_jets; + } + + // return the list of subjets obtained by unclustering the supplied jet down to nsub subjets. + std::vector RecJetFormat::exclusive_subjets(MAint32 nsub) const + { + std::vector output_jets; + try { + if (!RecJetFormat::has_exclusive_subjets()) + throw EXCEPTION_ERROR("This jet structure does not contain exclusive_subjets", + "Exclusive subjets only exist for jets clustered through an exclusive algorithm.", 0); + } catch (const std::exception& err) { + MANAGE_EXCEPTION(err); + return output_jets; + } + std::vector current_jets = fastjet::sorted_by_pt(pseudojet_.exclusive_subjets(nsub)); + output_jets.reserve(current_jets.size()); + for (auto &jet: current_jets) + { + RecJetFormat * NewJet = new RecJetFormat(jet); + output_jets.push_back(NewJet); + } + return output_jets; + } + + //returns true if the PseudoJet has support for exclusive subjets + MAbool RecJetFormat::has_exclusive_subjets() const { return pseudojet_.has_exclusive_subjets(); } +} +#endif diff --git a/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.h b/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.h index c6e16e31..e36745cf 100644 --- a/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.h +++ b/tools/SampleAnalyzer/Commons/DataFormat/RecJetFormat.h @@ -48,7 +48,7 @@ namespace MA5 class ClusterBase; class Pruner; class Nsubjettiness; - class SoftFrop; + class SoftDrop; class Filter; class EnergyCorrelator; } @@ -79,7 +79,7 @@ namespace MA5 friend class Substructure::ClusterBase; friend class Substructure::Pruner; friend class Substructure::Nsubjettiness; - friend class Substructure::SoftFrop; + friend class Substructure::SoftDrop; friend class Substructure::Filter; friend class Substructure::EnergyCorrelator; @@ -213,6 +213,16 @@ namespace MA5 // Accessor for pseudojets const fastjet::PseudoJet& pseudojet() const {return pseudojet_;} + // return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) + // that would be obtained when running the algorithm with the given dcut. + std::vector exclusive_subjets(MAfloat32 dcut) const; + + // return the list of subjets obtained by unclustering the supplied jet down to nsub subjets. + std::vector exclusive_subjets(MAint32 nsub) const; + + //returns true if the PseudoJet has support for exclusive subjets + MAbool has_exclusive_subjets() const; + // Add one pseudojet void setPseudoJet (const fastjet::PseudoJet& v) {pseudojet_=v;} void setPseudoJet (MALorentzVector& v) diff --git a/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.cpp b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.cpp new file mode 100644 index 00000000..8e9a0dbf --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.cpp @@ -0,0 +1,177 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#include "HEPTopTagger/HEPTopTagger.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h" + +namespace MA5 { + namespace Substructure { + + HTT::~HTT() { delete _tagger; } + + //============================// + // Initialization // + //============================// + + void HTT::Initialize(InputParameters& param) + { + _tagger = new fastjet::HEPTopTagger::HEPTopTagger(); + + // Optimal R + _tagger->do_optimalR(param.do_optimalR); + _tagger->set_optimalR_min(param.optimalR_min); + _tagger->set_optimalR_step(param.optimalR_step); + _tagger->set_optimalR_threshold(param.optimalR_threshold); + + // Candidate selection + fastjet::HEPTopTagger::Mode mode; + if (param.mode == HTT::EARLY_MASSRATIO_SORT_MASS) + mode = fastjet::HEPTopTagger::EARLY_MASSRATIO_SORT_MASS; + else if (param.mode == HTT::LATE_MASSRATIO_SORT_MASS) + mode = fastjet::HEPTopTagger::LATE_MASSRATIO_SORT_MASS; + else if (param.mode == HTT::EARLY_MASSRATIO_SORT_MODDJADE) + mode = fastjet::HEPTopTagger::EARLY_MASSRATIO_SORT_MODDJADE; + else if (param.mode == HTT::LATE_MASSRATIO_SORT_MODDJADE) + mode = fastjet::HEPTopTagger::LATE_MASSRATIO_SORT_MODDJADE; + else + mode = fastjet::HEPTopTagger::TWO_STEP_FILTER; + _tagger->set_mode(mode); + _tagger->set_mt(param.top_mass); + _tagger->set_mw(param.W_mass); + _tagger->set_top_mass_range(param.Mtop_min, param.Mtop_max); + _tagger->set_fw(param.fw); + _tagger->set_mass_ratio_range(param.mass_ratio_range_min, param.mass_ratio_range_max); + _tagger->set_mass_ratio_cut(param.m23cut, param.m13cutmin, param.m13cutmax); + + // Filtering + _tagger->set_filtering_n(param.filt_N); + _tagger->set_filtering_R(param.filtering_R); + _tagger->set_filtering_minpt_subjet(param.filtering_minpt); + + fastjet::JetAlgorithm algo_ = fastjet::antikt_algorithm; + if (param.filtering_algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; + else if (param.filtering_algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; + _tagger->set_filtering_jetalgorithm(algo_); + + // Reclustering + algo_ = fastjet::antikt_algorithm; + if (param.reclustering_algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; + else if (param.reclustering_algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; + _tagger->set_reclustering_jetalgorithm(algo_); + + // MassDrop + _tagger->set_mass_drop_threshold(param.mass_drop); + _tagger->set_mass_drop_threshold(param.mass_drop); + + // Pruning + _tagger->set_pruning_rcut_factor(param.prun_rcut); + _tagger->set_pruning_zcut(param.prun_zcut); + } + + //====================// + // Execute // + //====================// + + // Method to run top tagger + void HTT::Execute(const RecJetFormat *jet) { _tagger->run(jet->pseudojet()); } + + //======================// + // Accessors // + //======================// + + // accessor to the top jet + const RecJetFormat * HTT::top() const + { + RecJetFormat * NewJet = new RecJetFormat(const_cast(_tagger->t())); + return NewJet; + } + + // accessor to the bottom jet inside the top + const RecJetFormat * HTT::b() const + { + RecJetFormat * NewJet = new RecJetFormat(const_cast(_tagger->b())); + return NewJet; + } + + // accessor to the W jet inside the top + const RecJetFormat * HTT::W() const + { + RecJetFormat * NewJet = new RecJetFormat(const_cast(_tagger->W())); + return NewJet; + } + + // accessor to the leading subjet inside the W jet + const RecJetFormat * HTT::W1() const + { + RecJetFormat * NewJet = new RecJetFormat(const_cast(_tagger->W1())); + return NewJet; + } + + // accessor to second leading subjet inside the W jet + const RecJetFormat * HTT::W2() const + { + RecJetFormat * NewJet = new RecJetFormat(const_cast(_tagger->W2())); + return NewJet; + } + + // accessor to all PT-ordered subjets + std::vector HTT::subjets() const + { + std::vector output; + output.reserve(3); + RecJetFormat * j1 = new RecJetFormat(const_cast(_tagger->j1())); + RecJetFormat * j2 = new RecJetFormat(const_cast(_tagger->j2())); + RecJetFormat * j3 = new RecJetFormat(const_cast(_tagger->j3())); + output.push_back(j1); + output.push_back(j2); + output.push_back(j3); + return output; + } + + // print tagger information + void HTT::get_info() const { _tagger->get_info(); } + + // print tagger settings + void HTT::get_settings() const { _tagger->get_setting(); } + + // accessor to the pruned mass + MAfloat32 HTT::pruned_mass() const { return _tagger->pruned_mass(); } + + // accessor to the unfiltered mass + MAfloat32 HTT::unfiltered_mass() const { return _tagger->unfiltered_mass(); } + + // accessor to the difference between the reconstructed top mass and the true top mass + MAfloat32 HTT::delta_top() const { return _tagger->delta_top(); } + + // Is a given jet top-tagged + MAbool HTT::is_tagged() const { return _tagger->is_tagged(); } + + // Is the top mass window requirement satisfied? + MAbool HTT::is_maybe_top() const { return _tagger->is_maybe_top(); } + + // Are 2D mass plane requirements satisfied? + MAbool HTT::is_masscut_passed() const { return _tagger->is_masscut_passed(); } + } +} diff --git a/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h new file mode 100644 index 00000000..f56b3cbb --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h @@ -0,0 +1,168 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef MADANALYSIS5_HTT_H +#define MADANALYSIS5_HTT_H + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/Commons.h" + +namespace fastjet { + namespace HEPTopTagger { + class HEPTopTagger; + } +} + +namespace MA5 { + namespace Substructure { + class HTT { + + protected: + fastjet::HEPTopTagger::HEPTopTagger* _tagger; + + public: + + enum Mode { + EARLY_MASSRATIO_SORT_MASS, // applies 2D mass plane requirements then select the candidate which minimizes |m_cand-mt| + LATE_MASSRATIO_SORT_MASS, // selects the candidate which minimizes |m_cand-mt| + EARLY_MASSRATIO_SORT_MODDJADE, // applies the 2D mass plane requirements then select the candidate with highest jade distance + LATE_MASSRATIO_SORT_MODDJADE, // selects the candidate with highest modified jade distance + TWO_STEP_FILTER // only analyzes the candidate built with the highest pT(t) after unclustering + }; + + struct InputParameters { + Mode mode = EARLY_MASSRATIO_SORT_MASS; // execution mode + + MAbool do_optimalR = true; // initialize optimal R or set to fixed R + // optimal R parameters + MAfloat32 optimalR_min = 0.5; // min jet size + MAfloat32 optimalR_step = 0.1; // step size + MAfloat32 optimalR_threshold = 0.2; // step size + + // massdrop - unclustering + MAfloat32 mass_drop = 0.8; + MAfloat32 max_subjet = 30.; // set_max_subjet_mass + + // filtering + MAuint32 filt_N = 5; // set_nfilt + MAfloat32 filtering_R = 0.3; // max subjet distance for filtering + MAfloat32 filtering_minpt = 0.; // min subjet pt for filtering + // jet algorithm for filtering + Algorithm filtering_algorithm = Algorithm::cambridge; + + // Reclustering + // reclustering jet algorithm + Algorithm reclustering_algorithm = Algorithm::cambridge; + + //top mass range + MAfloat32 top_mass = 172.3; + MAfloat32 W_mass = 80.4; + MAfloat32 Mtop_min = 150.; + MAfloat32 Mtop_max = 200.; //set_top_range(min,max) + + // set top mass ratio range + MAfloat32 fw = 0.15; + MAfloat32 mass_ratio_range_min = (1.-fw)*W_mass/top_mass; + MAfloat32 mass_ratio_range_max = (1.+fw)*W_mass/top_mass; + + //mass ratio cuts + MAfloat32 m23cut = 0.35; + MAfloat32 m13cutmin = 0.2; + MAfloat32 m13cutmax = 1.3; + + // pruning + MAfloat32 prun_zcut = 0.1; // set_prun_zcut + MAfloat32 prun_rcut = .5; // set_prun_rcut + }; + + // Constructor without arguments + HTT() {} + + // Destructor + ~HTT(); + + //============================// + // Initialization // + //============================// + + HTT(HTT::InputParameters& param) { Initialize(param); } + + void Initialize(HTT::InputParameters& param); + + //====================// + // Execute // + //====================// + + // Method run top tagger. + void Execute(const RecJetFormat *jet); + + //======================// + // Accessors // + //======================// + + // accessor to top jet + const RecJetFormat * top() const; + + //accessor to bottom jet + const RecJetFormat * b() const; + + //accessor to W jet + const RecJetFormat * W() const; + + //accessor to leading subjet from W + const RecJetFormat * W1() const; + + //accessor to second leading subjet from W + const RecJetFormat * W2() const; + + // accessor to PT ordered subjets + std::vector subjets() const; + + // print tagger information + void get_info() const; + + // print tagger settings + void get_settings() const; + + // accessor to pruned mass + MAfloat32 pruned_mass() const; + + // accessor to unfiltered mass + MAfloat32 unfiltered_mass() const; + + // accessor to delta top + MAfloat32 delta_top() const; + + // Is given jet tagged + MAbool is_tagged() const; + + // top mass window requirement passed? + MAbool is_maybe_top() const; + + // 2D mass plane requirements passed? + MAbool is_masscut_passed() const; + }; + } +} + +#endif //MADANALYSIS5_HTT_H diff --git a/tools/SampleAnalyzer/Interfaces/HEPTopTagger/README.md b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/README.md new file mode 100644 index 00000000..ee3e0bcb --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/HEPTopTagger/README.md @@ -0,0 +1,107 @@ +# HEPTopTagger for MadAnalysis 5 +This module is based on the `Substructure` module of MadAnalysis 5, and is thus +constructed following the same principle. The corresponding class consists of an +`Initialize` and an `Execute` function. However, due to the complexity of the +HEPTopTgger algorithm, the `Execute` function does not return anything. Instead +it is equipped with various accessors. +For details see `tools/SampleAnalyzer/Interfaces/substructure/README.md`. + +## Installation +Simply call the `install HEPTopTagger` command through the MadAnalysis 5 command +line interface. Please make sure that the `FastJet` and `FastJet contrib` +packages have already been installed, and that they are available from +MadAnalysis 5. + +## Examples +```c++ +#ifndef analysis_test_h +#define analysis_test_h + +#include "SampleAnalyzer/Process/Analyzer/AnalyzerBase.h" + +namespace MA5 +{ + + class test : public AnalyzerBase + { + INIT_ANALYSIS(test,"test") + + private: + Substructure::HTT tagger; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize filtering + Substructure::InputParameters parameters; + parameters.mode = Substructure::HTT::EARLY_MASSRATIO_SORT_MASS; // execution mode + + parameters.do_optimalR = true; // initialize optimal R or set to fixed R + // optimal R parameters + parameters.optimalR_min = 0.5; // min jet size + parameters.optimalR_step = 0.1; // step size + parameters.optimalR_threshold = 0.2; // step size + + // massdrop - unclustering + parameters.mass_drop = 0.8; + parameters.max_subjet = 30.; // set_max_subjet_mass + + // filtering + parameters.filt_N = 5; // set_nfilt + parameters.filtering_R = 0.3; // max subjet distance for filtering + parameters.filtering_minpt = 0.; // min subjet pt for filtering + // jet algorithm for filtering + parameters.filtering_algorithm = Substructure::cambridge; + + // Reclustering + // reclustering jet algorithm + parameters.reclustering_algorithm = Substructure::Algorithm::cambridge; + + //top mass range + parameters.top_mass = 172.3; + parameters.W_mass = 80.4; + parameters.Mtop_min = 150.; + parameters.Mtop_max = 200.; //set_top_range(min,max) + + // set top mass ratio range + parameters.fw = 0.15; + parameters.mass_ratio_range_min = (1-parameters.fw)*parameters.W_mass/parameters.top_mass; + parameters.mass_ratio_range_max = (1+parameters.fw)*parameters.W_mass/parameters.top_mass; + + //mass ratio cuts + parameters.m23cut = 0.35; + parameters.m13cutmin = 0.2; + parameters.m13cutmax = 1.3; + + // pruning + parameters.prun_zcut = 0.1; // set_prun_zcut + parameters.prun_rcut = .5; // set_prun_rcut + tagger.Initialize(parameters); + return true; + } + virtual void Finalize(const SampleFormat& summary, const std::vector& files){} + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + + if (AK08.size() > 0) + { + taggedJet = tagger.Execute(AK08[0]); + if (taggedJet->is_tagged()) + INFO << "Top pT = " << taggedJet.top()->pt() << endmsg; + } + return true; + } + }; +} + +#endif +``` diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/ClusterBase.h b/tools/SampleAnalyzer/Interfaces/fastjet/ClusterBase.h deleted file mode 100644 index 6ddd09a1..00000000 --- a/tools/SampleAnalyzer/Interfaces/fastjet/ClusterBase.h +++ /dev/null @@ -1,268 +0,0 @@ -//////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks -// The MadAnalysis development team, email: -// -// This file is part of MadAnalysis 5. -// Official website: -// -// MadAnalysis 5 is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// MadAnalysis 5 is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with MadAnalysis 5. If not, see -// -//////////////////////////////////////////////////////////////////////////////// - -#ifndef MADANALYSIS5_CLUSTERBASE_H -#define MADANALYSIS5_CLUSTERBASE_H - -// STL headers -#include -#include - -// FastJet headers -#include "fastjet/ClusterSequence.hh" -#include "fastjet/PseudoJet.hh" - -// SampleAnalyser headers -#include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" -#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" -#include "SampleAnalyzer/Commons/DataFormat/EventFormat.h" -#include "SampleAnalyzer/Commons/Service/PDGService.h" -#include "SampleAnalyzer/Commons/Service/LogService.h" - -namespace fastjet -{ - class JetDefinition; -} - -using namespace std; - -namespace MA5{ - namespace Substructure { - - // Accessor for jet clustering algorithms - enum Algorithm {antikt, cambridge, kt}; - - class Recluster; - class SoftDrop; - class Pruner; - class Filter; - - class ClusterBase { - - friend class Recluster; - friend class SoftDrop; - friend class Pruner; - friend class Filter; - - //--------------------------------------------------------------------------------- - // data members - //--------------------------------------------------------------------------------- - protected: - - // External parameters - MAfloat32 ptmin_; // minimum transverse momentum - MAbool isExclusive_; // if false return a vector of all jets (in the sense of the inclusive algorithm) - // with pt >= ptmin. Time taken should be of the order of the number of jets - // returned. if True return a vector of all jets (in the sense of the exclusive - // algorithm) that would be obtained when running the algorithm with the given ptmin. - - /// Jet definition - fastjet::JetDefinition* JetDefinition_; - fastjet::JetDefinition::Plugin* JetDefPlugin_; - MAbool isPlugin_; - MAbool isClustered_; - - // Shared Cluster sequence - std::shared_ptr clust_seq; - - public: - - /// Constructor without argument - ClusterBase() {} - - /// Destructor - virtual ~ClusterBase() {} - - // Set the Jet definition using algorithm and radius input - void SetJetDef(Algorithm algorithm, MAfloat32 radius) - { - isPlugin_ = false; - JetDefinition_ = new fastjet::JetDefinition(__get_clustering_algorithm(algorithm), radius); - } - - //=======================// - // Execution // - //=======================// - - // Wrapper for event based execution - void Execute(const EventFormat& event, std::string JetID) - { __execute(const_cast(event), JetID); } - - // Execute with a single jet. This method reclusters the given jet using its constituents - std::vector Execute(const RecJetFormat *jet) - { - std::vector reclustered_jets = __cluster(jet->pseudojet().constituents()); - return __transform_jets(reclustered_jets); - } - - // Execute with a single jet. This method reclusters the given jet using its constituents by filtering - // reclustered events with respect to the initial jet - template - std::vector Execute(const RecJetFormat *jet, Func func) - { - std::vector output_jets; - std::vector reclustered_jets = __cluster(jet->pseudojet().constituents()); - - for (auto &recjet: reclustered_jets) - { - RecJetFormat *NewJet = __transform_jet(recjet); - if (func(jet, const_cast(NewJet))) output_jets.push_back(NewJet); - } - - return output_jets; - } - - // Execute with a list of jets. This method reclusters the given collection - // of jets by combining their constituents - std::vector Execute(std::vector &jets) - { - std::vector constituents; - for (auto &jet: jets) - { - std::vector current_constituents = jet->pseudojet().constituents(); - constituents.reserve(constituents.size() + current_constituents.size()); - constituents.insert( - constituents.end(), current_constituents.begin(), current_constituents.end() - ); - } - - std::vector reclustered_jets = __cluster(constituents); - - return __transform_jets(reclustered_jets); - } - - // Handler for clustering step - void cluster(const RecJetFormat *jet) - { - if (isPlugin_) - { - fastjet::JetDefinition jetDefinition(JetDefPlugin_); - clust_seq.reset(new fastjet::ClusterSequence(jet->pseudojet().constituents(), jetDefinition)); - } else { - clust_seq.reset(new fastjet::ClusterSequence( - jet->pseudojet().constituents(), - const_cast(*JetDefinition_) - )); - } - isClustered_ = true; - } - - // return a vector of all jets when the event is clustered (in the exclusive sense) to exactly njets. - // If there are fewer than njets particles in the ClusterSequence the function just returns however many - // particles there were. - std::vector exclusive_jets_up_to(MAint32 njets) - { - if (!isClustered_) throw EXCEPTION_ERROR("No clustered jet available", "", 1); - std::vector clustered_jets = fastjet::sorted_by_pt( - clust_seq->exclusive_jets_up_to(njets) - ); - return __transform_jets(clustered_jets); - } - - private: - - // Generic clustering method - std::vector __cluster(std::vector particles) - { - if (isPlugin_) - { - fastjet::JetDefinition jetDefinition(JetDefPlugin_); - clust_seq.reset(new fastjet::ClusterSequence(particles, jetDefinition)); - } else { - clust_seq.reset(new fastjet::ClusterSequence( - particles, const_cast(*JetDefinition_) - )); - } - isClustered_ = true; - std::vector jets; - if (isExclusive_) jets = clust_seq->exclusive_jets(ptmin_); - else jets = clust_seq->inclusive_jets(ptmin_); - - return fastjet::sorted_by_pt(jets); - } - - // Method to transform pseudojet into recjetformat - static RecJetFormat * __transform_jet(fastjet::PseudoJet jet) - { - RecJetFormat * NewJet = new RecJetFormat(jet); - return NewJet; - } - - static std::vector __transform_jets(std::vector jets) - { - std::vector output_jets; - output_jets.reserve(jets.size()); - for (auto &jet: jets) - output_jets.push_back(__transform_jet(jet)); - return output_jets; - } - - // Method to get jet algorithm - static fastjet::JetAlgorithm __get_clustering_algorithm(Substructure::Algorithm algorithm) - { - fastjet::JetAlgorithm algo_; - if (algorithm == Substructure::antikt) algo_ = fastjet::antikt_algorithm; - else if (algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; - else if (algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; - else throw EXCEPTION_ERROR("Unknown algorithm","",1); - return algo_; - } - - // Execute with the Reconstructed event. This method creates a new Jet in RecEventFormat which - // can be accessed via JetID. The algorithm will only be executed if a unique JetID is given - MAbool __execute(EventFormat& myEvent, std::string JetID) - { - try { - if (myEvent.rec()->hasJetID(JetID)) - throw EXCEPTION_ERROR("Substructure::ClusterBase - Jet ID `" + JetID + \ - "` already exits. Skipping execution.","",1); - } catch (const std::exception& err) { - MANAGE_EXCEPTION(err); - return false; - } - - std::vector jets = __cluster(myEvent.rec()->cluster_inputs()); - - std::vector output_jets; - output_jets.reserve(jets.size()); - - for (auto &jet: jets) { - output_jets.emplace_back(jet); - std::vector constituents = clust_seq->constituents(jet); - output_jets.back().Constituents_.reserve(constituents.size()); - output_jets.back().ntracks_ = 0; - for (auto &constit: constituents) { - output_jets.back().Constituents_.emplace_back(constit.user_index()); - if (PDG->IsCharged(myEvent.mc()->particles()[constit.user_index()].pdgid())) - output_jets.back().ntracks_++; - } - } - myEvent.rec()->jetcollection_.insert(std::make_pair(JetID, output_jets)); - isClustered_ = false; - return true; - } - }; - } -} - -#endif //MADANALYSIS5_CLUSTERBASE_H diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Nsubjettiness.h b/tools/SampleAnalyzer/Interfaces/fastjet/Nsubjettiness.h deleted file mode 100644 index 3cc0ee56..00000000 --- a/tools/SampleAnalyzer/Interfaces/fastjet/Nsubjettiness.h +++ /dev/null @@ -1,193 +0,0 @@ -//////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks -// The MadAnalysis development team, email: -// -// This file is part of MadAnalysis 5. -// Official website: -// -// MadAnalysis 5 is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// MadAnalysis 5 is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with MadAnalysis 5. If not, see -// -//////////////////////////////////////////////////////////////////////////////// - -#ifndef MADANALYSIS5_NSUBJETTINESS_H -#define MADANALYSIS5_NSUBJETTINESS_H - -// STL headers -#include -#include - -// FastJet headers -#include "fastjet/contrib/Nsubjettiness.hh" - -// SampleAnalyser headers -#include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" -#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" - -using namespace MA5; -using namespace std; - -namespace MA5 { - namespace Substructure { - class Nsubjettiness { - //--------------------------------------------------------------------------------- - // data members - //--------------------------------------------------------------------------------- - protected: - MAint32 order_; - fastjet::contrib::AxesDefinition* axesdef_; - fastjet::contrib::MeasureDefinition* measuredef_; - - public: - - enum AxesDef { - KT_Axes, - CA_Axes, - AntiKT_Axes, // (R0) - WTA_KT_Axes, - WTA_CA_Axes, - // GenKT_Axes, // (p, R0 = infinity) - // WTA_GenKT_Axes, // (p, R0 = infinity) - // GenET_GenKT_Axes, // (delta, p, R0 = infinity) - Manual_Axes, - OnePass_KT_Axes, - OnePass_CA_Axes, - OnePass_AntiKT_Axes, // (R0) - OnePass_WTA_KT_Axes, - OnePass_WTA_CA_Axes, - // OnePass_GenKT_Axes, // (p, R0 = infinity) - // OnePass_WTA_GenKT_Axes, // (p, R0 = infinity) - // OnePass_GenET_GenKT_Axes, // (delta, p, R0 = infinity) - // OnePass_Manual_Axes, - // MultiPass_Axes, // (NPass) (currently only defined for KT_Axes) - // MultiPass_Manual_Axes, // (NPass) - // Comb_GenKT_Axes, // (nExtra, p, R0 = infinity) - // Comb_WTA_GenKT_Axes, // (nExtra, p, R0 = infinity) - // Comb_GenET_GenKT_Axes, // (nExtra, delta, p, R0 = infinity) - }; - - enum MeasureDef { - NormalizedMeasure, // (beta,R0) - UnnormalizedMeasure, // (beta) - NormalizedCutoffMeasure, // (beta,R0,Rcutoff) - UnnormalizedCutoffMeasure, // (beta,Rcutoff) - }; - - /// Constructor without argument - Nsubjettiness() {} - - /// Destructor - virtual ~Nsubjettiness() {} - - //============================// - // Initialization // - //============================// - // Initialize the parameters of the algorithm. Initialization includes multiple if conditions - // Hence it would be optimum execution to initialize the algorithm during the initialisation - // of the analysis - - // Constructor with arguments - Nsubjettiness( - MAint32 order, - AxesDef axesdef, - MeasureDef measuredef, - MAfloat32 beta, - MAfloat32 R0, - MAfloat32 Rcutoff=std::numeric_limits::max() - ) - { Initialize(order, axesdef, measuredef, beta, R0, Rcutoff); } - - void Initialize( - MAint32 order, - AxesDef axesdef, - MeasureDef measuredef, - MAfloat32 beta, - MAfloat32 R0, - MAfloat32 Rcutoff=std::numeric_limits::max() - ) - { - if (axesdef == Substructure::Nsubjettiness::KT_Axes) - axesdef_ = new fastjet::contrib::KT_Axes(); - else if (axesdef == Substructure::Nsubjettiness::CA_Axes) - axesdef_ = new fastjet::contrib::CA_Axes(); - else if (axesdef == Substructure::Nsubjettiness::AntiKT_Axes) - axesdef_ = new fastjet::contrib::AntiKT_Axes(R0); - else if (axesdef == Substructure::Nsubjettiness::WTA_KT_Axes) - axesdef_ = new fastjet::contrib::WTA_KT_Axes(); - else if (axesdef == Substructure::Nsubjettiness::WTA_CA_Axes) - axesdef_ = new fastjet::contrib::WTA_CA_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::GenKT_Axes) -// axesdef_ = new fastjet::contrib::GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::WTA_GenKT_Axes) -// axesdef_ = new fastjet::contrib::WTA_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::GenET_GenKT_Axes) -// axesdef_ = new fastjet::contrib::GenET_GenKT_Axes(); - else if (axesdef == Substructure::Nsubjettiness::Manual_Axes) - axesdef_ = new fastjet::contrib::Manual_Axes(); - else if (axesdef == Substructure::Nsubjettiness::OnePass_KT_Axes) - axesdef_ = new fastjet::contrib::OnePass_KT_Axes(); - else if (axesdef == Substructure::Nsubjettiness::OnePass_CA_Axes) - axesdef_ = new fastjet::contrib::OnePass_CA_Axes(); - else if (axesdef == Substructure::Nsubjettiness::OnePass_AntiKT_Axes) - axesdef_ = new fastjet::contrib::OnePass_AntiKT_Axes(R0); - else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_KT_Axes) - axesdef_ = new fastjet::contrib::OnePass_WTA_KT_Axes(); - else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_CA_Axes) - axesdef_ = new fastjet::contrib::OnePass_WTA_CA_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::OnePass_GenKT_Axes) -// axesdef_ = new fastjet::contrib::OnePass_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_GenKT_Axes) -// axesdef_ = new fastjet::contrib::OnePass_WTA_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::OnePass_GenET_GenKT_Axes) -// axesdef_ = new fastjet::contrib::OnePass_GenET_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::OnePass_Manual_Axes) -// axesdef_ = new fastjet::contrib::OnePass_Manual_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::MultiPass_Axes) -// axesdef_ = new fastjet::contrib::MultiPass_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::MultiPass_Manual_Axes) -// axesdef_ = new fastjet::contrib::MultiPass_Manual_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::Comb_GenKT_Axes) -// axesdef_ = new fastjet::contrib::Comb_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::Comb_WTA_GenKT_Axes) -// axesdef_ = new fastjet::contrib::Comb_WTA_GenKT_Axes(); -// else if (axesdef == Substructure::Nsubjettiness::Comb_GenET_GenKT_Axes) -// axesdef_ = new fastjet::contrib::Comb_GenET_GenKT_Axes(); - - if (measuredef == Substructure::Nsubjettiness::NormalizedCutoffMeasure) - measuredef_ = new fastjet::contrib::NormalizedCutoffMeasure(beta, R0, Rcutoff); - else if (measuredef == Substructure::Nsubjettiness::NormalizedMeasure) - measuredef_ = new fastjet::contrib::NormalizedMeasure(beta, R0); - else if (measuredef == Substructure::Nsubjettiness::UnnormalizedMeasure) - measuredef_ = new fastjet::contrib::UnnormalizedMeasure(beta); - else if (measuredef == Substructure::Nsubjettiness::UnnormalizedCutoffMeasure) - measuredef_ = new fastjet::contrib::UnnormalizedCutoffMeasure(beta, Rcutoff); - - order_ = order; - } - - //=======================// - // Execution // - //=======================// - - // Method to calculate nsub for a given jet with respect to initialization parameters - MAdouble64 Execute(const RecJetFormat *jet) - { - fastjet::contrib::Nsubjettiness nsubjettiness(order_, *axesdef_, *measuredef_); - return nsubjettiness(jet->pseudojet()); - } - }; - } -} - -#endif //MADANALYSIS5_NSUBJETTINESS_H diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Cluster.cpp b/tools/SampleAnalyzer/Interfaces/substructure/Cluster.cpp new file mode 100644 index 00000000..cb0cfebc --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Cluster.cpp @@ -0,0 +1,36 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#include "SampleAnalyzer/Interfaces/substructure/Cluster.h" + +namespace MA5 { + namespace Substructure { + void Cluster::Initialize( + Algorithm algorithm, MAfloat32 radius, MAfloat32 ptmin, MAbool isExclusive + ) + { + SetJetDef(algorithm, radius); + ptmin_ = ptmin < 0. ? 0. : ptmin; isExclusive_ = isExclusive; + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Cluster.h b/tools/SampleAnalyzer/Interfaces/substructure/Cluster.h similarity index 84% rename from tools/SampleAnalyzer/Interfaces/fastjet/Cluster.h rename to tools/SampleAnalyzer/Interfaces/substructure/Cluster.h index 91e76375..6b6887ce 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/Cluster.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/Cluster.h @@ -25,9 +25,7 @@ #define MADANALYSIS5_CLUSTER_H // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" - -using namespace std; +#include "SampleAnalyzer/Interfaces/substructure/ClusterBase.h" namespace MA5 { namespace Substructure{ @@ -42,7 +40,7 @@ namespace MA5 { Cluster() {} /// Destructor - virtual ~Cluster() {} + ~Cluster() {} //============================// // Initialization // @@ -55,11 +53,9 @@ namespace MA5 { Cluster(Algorithm algorithm, MAfloat32 radius, MAfloat32 ptmin=0., MAbool isExclusive = false) { Initialize(algorithm, radius, ptmin, isExclusive); } - void Initialize(Algorithm algorithm, MAfloat32 radius, MAfloat32 ptmin=0., MAbool isExclusive = false) - { - SetJetDef(algorithm, radius); - ptmin_ = ptmin < 0. ? 0. : ptmin; isExclusive_ = isExclusive; - } + void Initialize( + Algorithm algorithm, MAfloat32 radius, MAfloat32 ptmin=0., MAbool isExclusive = false + ); }; } diff --git a/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.cpp b/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.cpp new file mode 100644 index 00000000..cd928858 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.cpp @@ -0,0 +1,209 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/ClusterSequence.hh" +#include "fastjet/PseudoJet.hh" + +#include "SampleAnalyzer/Interfaces/substructure/ClusterBase.h" + +namespace MA5 { + namespace Substructure { + + //=======================// + // Initialize // + //=======================// + + // Set the Jet definition using algorithm and radius input + void ClusterBase::SetJetDef(Algorithm algorithm, MAfloat32 radius) + { + isPlugin_ = false; + JetDefinition_ = new fastjet::JetDefinition(__get_clustering_algorithm(algorithm), radius); + } + + //=======================// + // Execution // + //=======================// + + // Wrapper for event based execution + void ClusterBase::Execute(const EventFormat& event, std::string JetID) + { __execute(const_cast(event), JetID); } + + // Execute with a single jet. This method reclusters the given jet using its constituents + std::vector ClusterBase::Execute(const RecJetFormat *jet) + { + std::vector reclustered_jets = __cluster(jet->pseudojet().constituents()); + return __transform_jets(reclustered_jets); + } + + // Execute with a single jet. This method reclusters the given jet using its constituents by filtering + // reclustered events with respect to the initial jet + template + std::vector ClusterBase::Execute(const RecJetFormat *jet, Func func) + { + std::vector output_jets; + std::vector reclustered_jets = __cluster(jet->pseudojet().constituents()); + + for (auto &recjet: reclustered_jets) + { + RecJetFormat *NewJet = __transform_jet(recjet); + if (func(jet, const_cast(NewJet))) output_jets.push_back(NewJet); + } + + return output_jets; + } + + // Execute with a list of jets. This method reclusters the given collection + // of jets by combining their constituents + std::vector ClusterBase::Execute(std::vector &jets) + { + std::vector constituents; + for (auto &jet: jets) + { + std::vector current_constituents = jet->pseudojet().constituents(); + constituents.reserve(constituents.size() + current_constituents.size()); + constituents.insert( + constituents.end(), current_constituents.begin(), current_constituents.end() + ); + } + + std::vector reclustered_jets = __cluster(constituents); + + return __transform_jets(reclustered_jets); + } + + // Handler for clustering step + void ClusterBase::cluster(const RecJetFormat *jet) + { + if (isPlugin_) + { + fastjet::JetDefinition jetDefinition(JetDefPlugin_); + clust_seq.reset(new fastjet::ClusterSequence(jet->pseudojet().constituents(), jetDefinition)); + } else { + clust_seq.reset(new fastjet::ClusterSequence( + jet->pseudojet().constituents(), + const_cast(*JetDefinition_) + )); + } + isClustered_ = true; + } + + // return a vector of all jets when the event is clustered (in the exclusive sense) to exactly njets. + // If there are fewer than njets particles in the ClusterSequence the function just returns however many + // particles there were. + std::vector ClusterBase::exclusive_jets_up_to(MAint32 njets) + { + if (!isClustered_) throw EXCEPTION_ERROR("No clustered jet available", "", 1); + std::vector clustered_jets = fastjet::sorted_by_pt( + clust_seq->exclusive_jets_up_to(njets) + ); + return __transform_jets(clustered_jets); + } + + //=======================// + // Private Functions // + //=======================// + + // Generic clustering method + std::vector ClusterBase::__cluster(std::vector particles) + { + if (isPlugin_) + { + fastjet::JetDefinition jetDefinition(JetDefPlugin_); + clust_seq.reset(new fastjet::ClusterSequence(particles, jetDefinition)); + } else { + clust_seq.reset(new fastjet::ClusterSequence( + particles, const_cast(*JetDefinition_) + )); + } + isClustered_ = true; + std::vector jets; + if (isExclusive_) jets = clust_seq->exclusive_jets(ptmin_); + else jets = clust_seq->inclusive_jets(ptmin_); + + return fastjet::sorted_by_pt(jets); + } + + // Method to transform pseudojet into recjetformat + RecJetFormat * ClusterBase::__transform_jet(fastjet::PseudoJet jet) const + { + RecJetFormat * NewJet = new RecJetFormat(jet); + return NewJet; + } + + // Transform pseudojets into RecJetFormat + std::vector ClusterBase::__transform_jets(std::vector jets) const + { + std::vector output_jets; + output_jets.reserve(jets.size()); + for (auto &jet: jets) + output_jets.push_back(__transform_jet(jet)); + return output_jets; + } + + // Method to get jet algorithm + fastjet::JetAlgorithm ClusterBase::__get_clustering_algorithm(Substructure::Algorithm algorithm) const + { + fastjet::JetAlgorithm algo_; + if (algorithm == Substructure::antikt) algo_ = fastjet::antikt_algorithm; + else if (algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; + else if (algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; + else throw EXCEPTION_ERROR("Unknown algorithm","",1); + return algo_; + } + + // Execute with the Reconstructed event. This method creates a new Jet in RecEventFormat which + // can be accessed via JetID. The algorithm will only be executed if a unique JetID is given + MAbool ClusterBase::__execute(EventFormat& myEvent, std::string JetID) + { + try { + if (myEvent.rec()->hasJetID(JetID)) + throw EXCEPTION_ERROR("Substructure::ClusterBase - Jet ID `" + JetID + \ + "` already exits. Skipping execution.","",1); + } catch (const std::exception& err) { + MANAGE_EXCEPTION(err); + return false; + } + + std::vector jets = __cluster(myEvent.rec()->cluster_inputs()); + + std::vector output_jets; + output_jets.reserve(jets.size()); + + for (auto &jet: jets) { + output_jets.emplace_back(jet); + std::vector constituents = clust_seq->constituents(jet); + output_jets.back().Constituents_.reserve(constituents.size()); + output_jets.back().ntracks_ = 0; + for (auto &constit: constituents) { + output_jets.back().Constituents_.emplace_back(constit.user_index()); + if (PDG->IsCharged(myEvent.mc()->particles()[constit.user_index()].pdgid())) + output_jets.back().ntracks_++; + } + } + myEvent.rec()->jetcollection_.insert(std::make_pair(JetID, output_jets)); + isClustered_ = false; + return true; + } + } +} diff --git a/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.h b/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.h new file mode 100644 index 00000000..a5f1a155 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/ClusterBase.h @@ -0,0 +1,140 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef MADANALYSIS5_CLUSTERBASE_H +#define MADANALYSIS5_CLUSTERBASE_H + +// STL headers +#include +#include + +// FastJet headers +#include "fastjet/ClusterSequence.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Commons/DataFormat/EventFormat.h" +#include "SampleAnalyzer/Commons/Service/PDGService.h" +#include "SampleAnalyzer/Commons/Service/LogService.h" +#include "SampleAnalyzer/Interfaces/substructure/Commons.h" + +namespace fastjet +{ + class JetDefinition; + class PseudoJet; +} + +namespace MA5{ + namespace Substructure { + + class Recluster; + + class ClusterBase { + + friend class Recluster; + + //--------------------------------------------------------------------------------- + // data members + //--------------------------------------------------------------------------------- + protected: + + // External parameters + MAfloat32 ptmin_; // minimum transverse momentum + MAbool isExclusive_; // if false return a vector of all jets (in the sense of the inclusive algorithm) + // with pt >= ptmin. Time taken should be of the order of the number of jets + // returned. if True return a vector of all jets (in the sense of the exclusive + // algorithm) that would be obtained when running the algorithm with the given ptmin. + + /// Jet definition + fastjet::JetDefinition* JetDefinition_; + fastjet::JetDefinition::Plugin* JetDefPlugin_; + MAbool isPlugin_; + MAbool isClustered_; + + // Shared Cluster sequence + std::shared_ptr clust_seq; + + public: + + /// Constructor without argument + ClusterBase() {} + + /// Destructor + virtual ~ClusterBase() + { + // clean heap allocation + delete JetDefinition_; + delete JetDefPlugin_; + } + + // Set the Jet definition using algorithm and radius input + void SetJetDef(Algorithm algorithm, MAfloat32 radius); + + //=======================// + // Execution // + //=======================// + + // Wrapper for event based execution + virtual void Execute(const EventFormat& event, std::string JetID); + + // Execute with a single jet. This method reclusters the given jet using its constituents + std::vector Execute(const RecJetFormat *jet); + + // Execute with a single jet. This method reclusters the given jet using its constituents by filtering + // reclustered events with respect to the initial jet + template + std::vector Execute(const RecJetFormat *jet, Func func); + + // Execute with a list of jets. This method reclusters the given collection + // of jets by combining their constituents + virtual std::vector Execute(std::vector &jets); + + // Handler for clustering step + void cluster(const RecJetFormat *jet); + + // return a vector of all jets when the event is clustered (in the exclusive sense) to exactly njets. + // If there are fewer than njets particles in the ClusterSequence the function just returns however many + // particles there were. + std::vector exclusive_jets_up_to(MAint32 njets); + + private: + + // Generic clustering method + std::vector __cluster(std::vector particles); + + // Method to transform pseudojet into recjetformat + RecJetFormat * __transform_jet(fastjet::PseudoJet jet) const; + + // Transform pseudojets into RecJetFormat + std::vector __transform_jets(std::vector jets) const; + + // Method to get jet algorithm + fastjet::JetAlgorithm __get_clustering_algorithm(Substructure::Algorithm algorithm) const; + + // Execute with the Reconstructed event. This method creates a new Jet in RecEventFormat which + // can be accessed via JetID. The algorithm will only be executed if a unique JetID is given + MAbool __execute(EventFormat& myEvent, std::string JetID); + }; + } +} + +#endif //MADANALYSIS5_CLUSTERBASE_H diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Commons.h b/tools/SampleAnalyzer/Interfaces/substructure/Commons.h new file mode 100644 index 00000000..baa83ce9 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Commons.h @@ -0,0 +1,40 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef MADANALYSIS5_COMMONS_H +#define MADANALYSIS5_COMMONS_H + +// SampleAnalyser headers +#include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" +#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" + +namespace MA5{ + namespace Substructure { + + // Accessor for jet clustering algorithms + enum Algorithm {antikt, cambridge, kt}; + + } +} + +#endif //MADANALYSIS5_COMMONS_H diff --git a/tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.cpp b/tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.cpp new file mode 100644 index 00000000..872fe09d --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.cpp @@ -0,0 +1,66 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/contrib/EnergyCorrelator.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.h" + +namespace MA5 { + namespace Substructure { + + EnergyCorrelator::~EnergyCorrelator() + { + delete _EC; + } + + void EnergyCorrelator::Initialize( + MAuint32 N, + MAfloat32 beta, + EnergyCorrelator::Measure measure, + EnergyCorrelator::Strategy strategy + ) + { + fastjet::contrib::EnergyCorrelator::Measure measure_; + fastjet::contrib::EnergyCorrelator::Strategy strategy_; + + if (measure == EnergyCorrelator::Measure::pt_R) + measure_ = fastjet::contrib::EnergyCorrelator::Measure::pt_R; + else if (measure == EnergyCorrelator::Measure::E_theta) + measure_ = fastjet::contrib::EnergyCorrelator::Measure::E_theta; + else + measure_ = fastjet::contrib::EnergyCorrelator::Measure::E_inv; + + if (strategy == EnergyCorrelator::Strategy::storage_array) + strategy_ = fastjet::contrib::EnergyCorrelator::Strategy::storage_array; + else strategy_ = fastjet::contrib::EnergyCorrelator::Strategy::slow; + + _EC = new fastjet::contrib::EnergyCorrelator(N, beta, measure_, strategy_); + } + + // Method to execute with a single jet + MAdouble64 EnergyCorrelator::Execute(const RecJetFormat* jet) const + { return (*_EC)(jet->pseudojet()); } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/EnergyCorrelator.h b/tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.h similarity index 77% rename from tools/SampleAnalyzer/Interfaces/fastjet/EnergyCorrelator.h rename to tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.h index 2f308f1a..9afc2d5d 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/EnergyCorrelator.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.h @@ -28,15 +28,18 @@ #include #include -// FastJet headers -#include "fastjet/contrib/EnergyCorrelator.hh" - // SampleAnalyser headers #include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" #include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" using namespace std; +namespace fastjet { + namespace contrib { + class EnergyCorrelator; + } +} + namespace MA5 { namespace Substructure { class EnergyCorrelator { @@ -70,7 +73,7 @@ namespace MA5 { EnergyCorrelator() {} /// Destructor - virtual ~EnergyCorrelator() {} + ~EnergyCorrelator(); /// constructs an N-point correlator with angular exponent beta, /// using the specified choice of energy and angular measure as well @@ -87,27 +90,10 @@ namespace MA5 { MAfloat32 beta, EnergyCorrelator::Measure measure = EnergyCorrelator::Measure::pt_R, EnergyCorrelator::Strategy strategy = EnergyCorrelator::Strategy::storage_array - ) - { - fastjet::contrib::EnergyCorrelator::Measure measure_; - fastjet::contrib::EnergyCorrelator::Strategy strategy_; - - if (measure == EnergyCorrelator::Measure::pt_R) - measure_ = fastjet::contrib::EnergyCorrelator::Measure::pt_R; - else if (measure == EnergyCorrelator::Measure::E_theta) - measure_ = fastjet::contrib::EnergyCorrelator::Measure::E_theta; - else - measure_ = fastjet::contrib::EnergyCorrelator::Measure::E_inv; - - if (strategy == EnergyCorrelator::Strategy::storage_array) - strategy_ = fastjet::contrib::EnergyCorrelator::Strategy::storage_array; - else strategy_ = fastjet::contrib::EnergyCorrelator::Strategy::slow; - - _EC = new fastjet::contrib::EnergyCorrelator(N, beta, measure_, strategy_); - } + ); // Method to execute with a single jet - MAdouble64 Execute(const RecJetFormat* jet) { return (*_EC)(jet->pseudojet()); } + MAdouble64 Execute(const RecJetFormat* jet) const; }; } } diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Filter.cpp b/tools/SampleAnalyzer/Interfaces/substructure/Filter.cpp new file mode 100644 index 00000000..1f09546e --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Filter.cpp @@ -0,0 +1,100 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/tools/Filter.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/Filter.h" + +namespace MA5 { + namespace Substructure { + + Filter::~Filter() + { + delete JetDefinition_; + delete JetFilter_; + } + + //============================// + // Initialization // + //============================// + + void Filter::Initialize(Algorithm algorithm, MAfloat32 radius, Selector selector, MAfloat32 rho) + { + fastjet::JetAlgorithm algo_ = fastjet::antikt_algorithm; + if (algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; + else if (algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; + + JetDefinition_ = new fastjet::JetDefinition(algo_, radius); + Rfilt_=-1.; rho_=rho; + init_filter(selector, true); + } + + void Filter::Initialize(MAfloat32 Rfilt, Selector selector, MAfloat32 rho) + { + Rfilt_=Rfilt; rho_=rho; + init_filter(selector,false); + } + + void Filter::init_filter(Selector selector, MAbool isJetDefined) + { + if (isJetDefined) + JetFilter_ = new fastjet::Filter(*JetDefinition_, selector.selector_, rho_); + else + JetFilter_ = new fastjet::Filter(Rfilt_, selector.selector_, rho_); + } + + //=======================// + // Execution // + //=======================// + + // Method to filter a given jet. + const RecJetFormat* Filter::Execute(const RecJetFormat *jet) const + { + fastjet::PseudoJet filtered_jet = (*JetFilter_)(jet->pseudojet()); + RecJetFormat * NewJet = new RecJetFormat(filtered_jet); + return NewJet; + } + + // Method to filter all the jets in a vector + std::vector Filter::Execute(std::vector &jets) const + { + std::vector output_jets; + output_jets.reserve(jets.size()); + for (auto &jet: jets) + output_jets.push_back(Execute(jet)); + + std::sort( + output_jets.begin(), + output_jets.end(), + [](const RecJetFormat *j1, const RecJetFormat *j2) + { + return (j1->pt() > j2->pt()); + } + ); + + return output_jets; + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Filter.h b/tools/SampleAnalyzer/Interfaces/substructure/Filter.h similarity index 66% rename from tools/SampleAnalyzer/Interfaces/fastjet/Filter.h rename to tools/SampleAnalyzer/Interfaces/substructure/Filter.h index 09de5782..a233fcfd 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/Filter.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/Filter.h @@ -24,14 +24,13 @@ #ifndef MADANALYSIS5_FILTER_H #define MADANALYSIS5_FILTER_H -// FastJet headers -#include "fastjet/tools/Filter.hh" - // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" -#include "SampleAnalyzer/Interfaces/fastjet/Selector.h" +#include "SampleAnalyzer/Interfaces/substructure/Commons.h" +#include "SampleAnalyzer/Interfaces/substructure/Selector.h" -using namespace std; +namespace fastjet { + class Filter; +} namespace MA5 { namespace Substructure { @@ -54,11 +53,9 @@ namespace MA5 { // data members //--------------------------------------------------------------------------------- protected: - fastjet::Selector selector_; // the Selector applied to compute the kept subjets /// Jet definition fastjet::JetDefinition* JetDefinition_; // the jet definition applied to obtain the subjets - MAbool isJetDefined_; MAfloat32 rho_; // if non-zero, backgruond-subtract each subjet befor selection MAfloat32 Rfilt_; // the filtering radius @@ -73,7 +70,7 @@ namespace MA5 { Filter() {} /// Destructor - virtual ~Filter() {} + ~Filter(); //============================// // Initialization // @@ -86,60 +83,22 @@ namespace MA5 { Filter(MAfloat32 Rfilt, Selector selector, MAfloat32 rho=0.) { Initialize(Rfilt, selector, rho); } - void Initialize(Algorithm algorithm, MAfloat32 radius, Selector selector, MAfloat32 rho=0.) - { - JetDefinition_ = new fastjet::JetDefinition( - ClusterBase().__get_clustering_algorithm(algorithm), radius - ); - selector_ = selector.__get(); Rfilt_=-1.; rho_=rho; isJetDefined_=true; - init_filter(); - } - - void Initialize(MAfloat32 Rfilt, Selector selector, MAfloat32 rho=0.) - { - selector_ = selector.__get(); Rfilt_=Rfilt; rho_=rho; isJetDefined_=false; - init_filter(); - } + void Initialize(Algorithm algorithm, MAfloat32 radius, Selector selector, MAfloat32 rho=0.); + void Initialize(MAfloat32 Rfilt, Selector selector, MAfloat32 rho=0.); //=======================// // Execution // //=======================// // Method to filter a given jet. - const RecJetFormat* Execute(const RecJetFormat *jet) - { - fastjet::PseudoJet filtered_jet = (*JetFilter_)(jet->pseudojet()); - return ClusterBase().__transform_jet(filtered_jet); - } + const RecJetFormat* Execute(const RecJetFormat *jet) const; // Method to filter all the jets in a vector - std::vector Execute(std::vector &jets) - { - std::vector output_jets; - for (auto &jet: jets) - output_jets.push_back(Execute(jet)); - - std::sort( - output_jets.begin(), - output_jets.end(), - [](const RecJetFormat *j1, const RecJetFormat *j2) - { - return (j1->pt() > j2->pt()); - } - ); - - return output_jets; - } + std::vector Execute(std::vector &jets) const; private: - void init_filter() - { - if (isJetDefined_) - JetFilter_ = new fastjet::Filter(*JetDefinition_, selector_, rho_); - else - JetFilter_ = new fastjet::Filter(Rfilt_, selector_, rho_); - } + void init_filter(Selector selector, MAbool isJetDefined); }; diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.cpp b/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.cpp new file mode 100644 index 00000000..9472d904 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.cpp @@ -0,0 +1,125 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/contrib/Nsubjettiness.hh" + +#include "SampleAnalyzer/Interfaces/substructure/Nsubjettiness.h" + +namespace MA5 { + namespace Substructure { + + Nsubjettiness::~Nsubjettiness() + { + // clean heap allocation + delete axesdef_; + delete measuredef_; + } + + //============================// + // Initialization // + //============================// + // Initialize the parameters of the algorithm. Initialization includes multiple if conditions + // Hence it would be optimum execution to initialize the algorithm during the initialisation + // of the analysis + + void Nsubjettiness::Initialize( + MAint32 order, + AxesDef axesdef, + MeasureDef measuredef, + MAfloat32 beta, + MAfloat32 R0, + MAfloat32 Rcutoff + ) + { + if (axesdef == Substructure::Nsubjettiness::KT_Axes) + axesdef_ = new fastjet::contrib::KT_Axes(); + else if (axesdef == Substructure::Nsubjettiness::CA_Axes) + axesdef_ = new fastjet::contrib::CA_Axes(); + else if (axesdef == Substructure::Nsubjettiness::AntiKT_Axes) + axesdef_ = new fastjet::contrib::AntiKT_Axes(R0); + else if (axesdef == Substructure::Nsubjettiness::WTA_KT_Axes) + axesdef_ = new fastjet::contrib::WTA_KT_Axes(); + else if (axesdef == Substructure::Nsubjettiness::WTA_CA_Axes) + axesdef_ = new fastjet::contrib::WTA_CA_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::GenKT_Axes) +// axesdef_ = new fastjet::contrib::GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::WTA_GenKT_Axes) +// axesdef_ = new fastjet::contrib::WTA_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::GenET_GenKT_Axes) +// axesdef_ = new fastjet::contrib::GenET_GenKT_Axes(); + else if (axesdef == Substructure::Nsubjettiness::Manual_Axes) + axesdef_ = new fastjet::contrib::Manual_Axes(); + else if (axesdef == Substructure::Nsubjettiness::OnePass_KT_Axes) + axesdef_ = new fastjet::contrib::OnePass_KT_Axes(); + else if (axesdef == Substructure::Nsubjettiness::OnePass_CA_Axes) + axesdef_ = new fastjet::contrib::OnePass_CA_Axes(); + else if (axesdef == Substructure::Nsubjettiness::OnePass_AntiKT_Axes) + axesdef_ = new fastjet::contrib::OnePass_AntiKT_Axes(R0); + else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_KT_Axes) + axesdef_ = new fastjet::contrib::OnePass_WTA_KT_Axes(); + else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_CA_Axes) + axesdef_ = new fastjet::contrib::OnePass_WTA_CA_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::OnePass_GenKT_Axes) +// axesdef_ = new fastjet::contrib::OnePass_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::OnePass_WTA_GenKT_Axes) +// axesdef_ = new fastjet::contrib::OnePass_WTA_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::OnePass_GenET_GenKT_Axes) +// axesdef_ = new fastjet::contrib::OnePass_GenET_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::OnePass_Manual_Axes) +// axesdef_ = new fastjet::contrib::OnePass_Manual_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::MultiPass_Axes) +// axesdef_ = new fastjet::contrib::MultiPass_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::MultiPass_Manual_Axes) +// axesdef_ = new fastjet::contrib::MultiPass_Manual_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::Comb_GenKT_Axes) +// axesdef_ = new fastjet::contrib::Comb_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::Comb_WTA_GenKT_Axes) +// axesdef_ = new fastjet::contrib::Comb_WTA_GenKT_Axes(); +// else if (axesdef == Substructure::Nsubjettiness::Comb_GenET_GenKT_Axes) +// axesdef_ = new fastjet::contrib::Comb_GenET_GenKT_Axes(); + + if (measuredef == Substructure::Nsubjettiness::NormalizedCutoffMeasure) + measuredef_ = new fastjet::contrib::NormalizedCutoffMeasure(beta, R0, Rcutoff); + else if (measuredef == Substructure::Nsubjettiness::NormalizedMeasure) + measuredef_ = new fastjet::contrib::NormalizedMeasure(beta, R0); + else if (measuredef == Substructure::Nsubjettiness::UnnormalizedMeasure) + measuredef_ = new fastjet::contrib::UnnormalizedMeasure(beta); + else if (measuredef == Substructure::Nsubjettiness::UnnormalizedCutoffMeasure) + measuredef_ = new fastjet::contrib::UnnormalizedCutoffMeasure(beta, Rcutoff); + + order_ = order; + } + + //=======================// + // Execution // + //=======================// + + // Method to calculate nsub for a given jet with respect to initialization parameters + MAdouble64 Nsubjettiness::Execute(const RecJetFormat *jet) const + { + fastjet::contrib::Nsubjettiness nsubjettiness(order_, *axesdef_, *measuredef_); + return nsubjettiness(jet->pseudojet()); + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.h b/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.h new file mode 100644 index 00000000..908aae81 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Nsubjettiness.h @@ -0,0 +1,132 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef MADANALYSIS5_NSUBJETTINESS_H +#define MADANALYSIS5_NSUBJETTINESS_H + +// STL headers +#include + +// SampleAnalyser headers +#include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" +#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" + +using namespace std; + +namespace fastjet { + namespace contrib { + class AxesDefinition; + class MeasureDefinition; + } +} + +namespace MA5 { + namespace Substructure { + class Nsubjettiness { + //--------------------------------------------------------------------------------- + // data members + //--------------------------------------------------------------------------------- + protected: + MAint32 order_; + fastjet::contrib::AxesDefinition* axesdef_; + fastjet::contrib::MeasureDefinition* measuredef_; + + public: + + enum AxesDef { + KT_Axes, + CA_Axes, + AntiKT_Axes, // (R0) + WTA_KT_Axes, + WTA_CA_Axes, + // GenKT_Axes, // (p, R0 = infinity) + // WTA_GenKT_Axes, // (p, R0 = infinity) + // GenET_GenKT_Axes, // (delta, p, R0 = infinity) + Manual_Axes, + OnePass_KT_Axes, + OnePass_CA_Axes, + OnePass_AntiKT_Axes, // (R0) + OnePass_WTA_KT_Axes, + OnePass_WTA_CA_Axes, + // OnePass_GenKT_Axes, // (p, R0 = infinity) + // OnePass_WTA_GenKT_Axes, // (p, R0 = infinity) + // OnePass_GenET_GenKT_Axes, // (delta, p, R0 = infinity) + // OnePass_Manual_Axes, + // MultiPass_Axes, // (NPass) (currently only defined for KT_Axes) + // MultiPass_Manual_Axes, // (NPass) + // Comb_GenKT_Axes, // (nExtra, p, R0 = infinity) + // Comb_WTA_GenKT_Axes, // (nExtra, p, R0 = infinity) + // Comb_GenET_GenKT_Axes, // (nExtra, delta, p, R0 = infinity) + }; + + enum MeasureDef { + NormalizedMeasure, // (beta,R0) + UnnormalizedMeasure, // (beta) + NormalizedCutoffMeasure, // (beta,R0,Rcutoff) + UnnormalizedCutoffMeasure, // (beta,Rcutoff) + }; + + /// Constructor without argument + Nsubjettiness() {} + + /// Destructor + ~Nsubjettiness(); + + //============================// + // Initialization // + //============================// + // Initialize the parameters of the algorithm. Initialization includes multiple if conditions + // Hence it would be optimum execution to initialize the algorithm during the initialisation + // of the analysis + + // Constructor with arguments + Nsubjettiness( + MAint32 order, + AxesDef axesdef, + MeasureDef measuredef, + MAfloat32 beta, + MAfloat32 R0, + MAfloat32 Rcutoff=std::numeric_limits::max() + ) + { Initialize(order, axesdef, measuredef, beta, R0, Rcutoff); } + + void Initialize( + MAint32 order, + AxesDef axesdef, + MeasureDef measuredef, + MAfloat32 beta, + MAfloat32 R0, + MAfloat32 Rcutoff=std::numeric_limits::max() + ); + + //=======================// + // Execution // + //=======================// + + // Method to calculate nsub for a given jet with respect to initialization parameters + MAdouble64 Execute(const RecJetFormat *jet) const; + }; + } +} + +#endif //MADANALYSIS5_NSUBJETTINESS_H diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Pruner.cpp b/tools/SampleAnalyzer/Interfaces/substructure/Pruner.cpp new file mode 100644 index 00000000..85e38a34 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Pruner.cpp @@ -0,0 +1,85 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/tools/Pruner.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/Pruner.h" + +namespace MA5 { + namespace Substructure { + + Pruner::~Pruner() + { + delete JetDefinition_; + } + + // Initialize with all the arguments. Note that if R <= 0 max allowable radius will be used + void Pruner::Initialize( + Substructure::Algorithm algorithm, MAfloat32 R, MAfloat32 zcut, MAfloat32 Rcut_factor + ) + { + fastjet::JetAlgorithm algo_ = fastjet::antikt_algorithm; + if (algorithm == Substructure::cambridge) algo_ = fastjet::cambridge_algorithm; + else if (algorithm == Substructure::kt) algo_ = fastjet::kt_algorithm; + + if (R<=0.) + JetDefinition_ = new fastjet::JetDefinition(algo_, fastjet::JetDefinition::max_allowable_R); + else + JetDefinition_ = new fastjet::JetDefinition(algo_, R); + + zcut_ = zcut; Rcut_factor_=Rcut_factor; + } + + //=======================// + // Execution // + //=======================// + + // Method to prune a given jet with respect to initialization parameters + const RecJetFormat * Pruner::Execute(const RecJetFormat *jet) const + { + fastjet::PseudoJet pruned_jet = __prune(jet->pseudojet()); + RecJetFormat * NewJet = new RecJetFormat(pruned_jet); + return NewJet; + } + + // Method to prune each given jet individually with respect to initialization parameters + std::vector Pruner::Execute(std::vector &jets) const + { + std::vector output_jets; + output_jets.reserve(jets.size()); + for (auto &jet: jets) + output_jets.push_back(Execute(jet)); + return output_jets; + } + + fastjet::PseudoJet Pruner::__prune(fastjet::PseudoJet jet) const + { + fastjet::Pruner pruner( + *const_cast(JetDefinition_), zcut_, Rcut_factor_ + ); + return pruner(jet); + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Pruner.h b/tools/SampleAnalyzer/Interfaces/substructure/Pruner.h similarity index 71% rename from tools/SampleAnalyzer/Interfaces/fastjet/Pruner.h rename to tools/SampleAnalyzer/Interfaces/substructure/Pruner.h index 953b0bf4..2dbf8db6 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/Pruner.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/Pruner.h @@ -24,13 +24,12 @@ #ifndef MADANALYSIS5_PRUNER_H #define MADANALYSIS5_PRUNER_H -// FastJet headers -#include "fastjet/tools/Pruner.hh" - // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" +#include "SampleAnalyzer/Interfaces/substructure/Commons.h" -using namespace std; +namespace fastjet { + class JetDefinition; +} namespace MA5 { namespace Substructure { @@ -52,7 +51,7 @@ namespace MA5 { Pruner() {} /// Destructor - virtual ~Pruner() {} + ~Pruner(); //============================// // Initialization // @@ -73,47 +72,23 @@ namespace MA5 { { Initialize(algorithm, -1., zcut, Rcut_factor); } // Initialize with all the arguments. Note that if R <= 0 max allowable radius will be used - void Initialize(Substructure::Algorithm algorithm, MAfloat32 R, MAfloat32 zcut, MAfloat32 Rcut_factor) - { - fastjet::JetAlgorithm algo_ = ClusterBase().__get_clustering_algorithm(algorithm); - - if (R<=0.) - JetDefinition_ = new fastjet::JetDefinition(algo_, fastjet::JetDefinition::max_allowable_R); - else - JetDefinition_ = new fastjet::JetDefinition(algo_, R); - - zcut_ = zcut; Rcut_factor_=Rcut_factor; - } + void Initialize( + Substructure::Algorithm algorithm, MAfloat32 R, MAfloat32 zcut, MAfloat32 Rcut_factor + ); //=======================// // Execution // //=======================// // Method to prune a given jet with respect to initialization parameters - const RecJetFormat * Execute(const RecJetFormat *jet) - { - fastjet::PseudoJet pruned_jet = __prune(jet->pseudojet()); - return ClusterBase().__transform_jet(pruned_jet); - } + const RecJetFormat * Execute(const RecJetFormat *jet) const; // Method to prune each given jet individually with respect to initialization parameters - std::vector Execute(std::vector &jets) - { - std::vector output_jets; - for (auto &jet: jets) - output_jets.push_back(Execute(jet)); - return output_jets; - } + std::vector Execute(std::vector &jets) const; private: - fastjet::PseudoJet __prune(fastjet::PseudoJet jet) - { - fastjet::Pruner pruner( - *const_cast(JetDefinition_), zcut_, Rcut_factor_ - ); - return pruner(jet); - } + fastjet::PseudoJet __prune(fastjet::PseudoJet jet) const; }; } diff --git a/tools/SampleAnalyzer/Interfaces/substructure/README.md b/tools/SampleAnalyzer/Interfaces/substructure/README.md new file mode 100644 index 00000000..cc624b57 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/README.md @@ -0,0 +1,456 @@ +# Usage of the substructure module shipped with MadAnalysis 5 + +The Substructure module of MadAnalysis 5 has been designed to be used within the +[SFS framework](https://arxiv.org/abs/2006.09387), in which jet clustering +information is preserved and used throughout the analysis execution. This allows +users to use the `Substructure` interface to the `FastJet` libraries without +having to deal with how jets are treated internally (in particular with respect +to memory management). The`Substructure` interface uses the clustering history +to link `RecJetFormat` objects to `FastJet`'s `PseudoJet` objects, and converts +the two types of object seamlesly in the background. + +Each class that has been constructed under the `Substructure` namespace has been +equiped with an `Initialize` and an `Execute` function. The `Initialize` +function creates a memory allocation for the corresponding `FastJet` (or +`fjcontib` backend). + +**Note:** Initialising any given of the `Substructure` classes during the +execution of the analysis can slow down the code and might create memory +allocation problems. The best practice requires to **always** initialise the +used substructure classes in the analysis initialisation routine. Examples are +provided below. + +Depending on the nature of the jet substructure tool considered, the `Execute` +function can take as an argument a `const RecJetFormat *` pointer or a +`std::vector` vector of pointers. The former possibility +is designed to work with a single jet, whereas the latter is designed to work +with a collection of jets. The precise action of the `Execute` function depends +on the substructure tools to be used. The `Execute` function, as suggested by +its name, needs to be executed during the analysis execution. + +## Outline +* [Examples](#examples) + * [SoftDrop](#softdrop) + * [Filtering](#filtering) + * [Energy Correlator](#energy-correlator) + * [Nsubjettiness](#nsubjettiness) + * [Clustering](#clustering) + * [VariableR](#variabler) + * [Pruning](#pruning) + * [Reclustering](#reclustering) + +# Examples +The examples below are designed to illutrate only analysis header files, but can +easily be generalised to be used separately in a header and a source file. We +only made this choice as it is easier to examplify the functionalities through a +single file. To create an environment compatible with the substructure module, +we begin with writing a file named `substructure_example.ma5` including the +commands shown below: +``` +# Define standard AK04 jet +set main.fastsim.package = fastjet +set main.fastsim.algorithm = antikt +set main.fastsim.radius = 0.4 +set main.fastsim.JetID = AK04 +set main.fastsim.bjet_id.matching_dr = 0.3 +set main.fastsim.bjet_id.exclusive = true +set main.fastsim.exclusive_id = false +set main.fastsim.ptmin = 20.0 + +# Define AK08 jet +define jet_algorithm AK08 antikt radius=0.8 ptmin=200 + +# Define a VariableR algorithm name "varR" +define jet_algorithm varR VariableR rho=2000 +set varR.minR = 0 +set varR.maxR = 2 +set varR.ptmin = 20 +``` +This will generate a MadAnalysis 5 workspace which includes a definition of a +primary jet class named `AK04`, and clustered according to the `antikt` +algorithm with the parameters `R=0.4` and `ptmin=20` GeV. Furthermore, the +commands lead to the initialisation of two additional jet classes named `AK08` +and `varR`. In the former case, jets are clustered with the `antikt` algorithm +with parameters `R=0.8` and `ptmin=200` GeV. The latter definition leads to a +jet definition through the `VariableR` algorithm. All available `jet_algorithm` +definitions are shown with their default values in the table below. + +| Algorithm | Parameters & Default values | +|:---------------------:|:-----------------------------------------------------------------------------------| +| `antikt`, `cambridge` | `radius=0.4`, `ptmin=5.` | +| `genkt` | `radius=0.4`, `ptmin=5.`, `exclusive=False`, `p=-1` | +| `kt` | `radius=0.4`, `ptmin=5.`, `exclusive=False` | +| `gridjet` | `ymax=3.`, `ptmin=5.` | +| `cdfjetclu` | `radius=0.4`, `ptmin=5.`, `overlap=0.5`, `seed=1.`, `iratch=0.` | +| `cdfmidpoint` | `radius=0.4`, `ptmin=5.`, `overlap=0.5`, `seed=1.`, `iratch=0.`, `areafraction=1.` | +| `siscone` | `radius=0.4`, `ptmin=5.`, `overlap=0.5`, `input_ptmin=5.`, `npassmax=1.` | +| `VariableR` | `rho=2000.`, `minR=0.`, `maxR=2.`, `ptmin=20.` `exclusive=False` `clustertype=CALIKE` `strategy=Best` | + +The above MadAnalysis 5 script can be executed as: +```bash +cd madanalysis5 +./bin/ma5 -Re my_workspace test substructure_example.ma5 +``` +This will create a MadAnalysus 5 workspace in the Expert mode of the code. The +workspace is named `my_workspace` and the analysis file is named `test`. The +`test.h` and `test.cpp` files can be found under + `my_workspace/Build/SampleAnalyzer/User/Analyzer`. +The structure of `test.h` reads: + +`my_workspace/Build/SampleAnalyzer/User/Analyzer/test.h`: +```c++ +#ifndef analysis_test_h +#define analysis_test_h + +#include "SampleAnalyzer/Process/Analyzer/AnalyzerBase.h" + +namespace MA5 +{ +class test : public AnalyzerBase +{ + INIT_ANALYSIS(test,"test") + + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters); + virtual void Finalize(const SampleFormat& summary, const std::vector& files); + virtual bool Execute(SampleFormat& sample, const EventFormat& event); + + private: +}; +} + +#endif +``` +For the purposes of this example `test.cpp` is not be used and can hence be +removed. We will only modify `test.h` from now on. + +[Back to top](#outline) + +## SoftDrop: + +First we need to create a private `softdrop` object to be able to access it +thorough the different functions of the `test` class +```c++ +private: + Substructure::SoftDrop softDrop; +``` +Then during the `Initialize` function, it can be called as +```c++ +virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) +{ + // Initialize SoftDrop + MAfloat32 z_cut = 0.10; + MAfloat32 beta = 2.0; + softDrop.Initialize(beta, z_cut); + return true; +} +``` +Note that it is possible to define a `SoftDrop` object through a pointer, + `Substructure::SoftDrop* softDrop;` +In this case, the class should be initialised as: + `softDrop = new Substructure::SoftDrop(beta, z_cut);`. + +The `SoftDrop` object can then be further used during the execution of the code. +```c++ +virtual bool Execute(SampleFormat& sample, const EventFormat& event) +{ + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + + if (AK08.size() > 0) + { + INFO << "SoftDrop on leading AK08 jet:" << endmsg; + const RecJet softdrop_jet = softDrop.Execute(AK08[0]); + INFO << "M_SD = " << softdrop_jet->m() << " M(j_1) = " << AK08[0]->m() << endmsg; + } + return true; +} +``` +Note that if `SoftDrop` is defined as a pointer it needs to be executed as + `const RecJet softdrop_jet = softDrop->Execute(AK08[0]);`. +Here the `filter` function cleans the "AK08" jet collection defined +[above](#examples) and removes from it any jet with pT<200 GeV and |eta|>2.5. In +the following we apply the `softDrop` method on the leading `AK08` jet, which +results in a `const RecJetFormat *` object, the `RecJet` definition used here +having been defined in the `SampleAnalyzer` backend. + +**Note:** All substructure classes are desinged in the same way. Therefore, from +now on we only show simple examples of usage without any detailed explanations. + +[Back to top](#outline) + +## Filtering +```c++ + private: + Substructure::Filter JetFilter; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize filtering + MAfloat32 Rfilt = 0.2; + INFO << "Initializing Filter" << endmsg; + JetFilter.Initialize(Rfilt, Substructure::SelectorPtFractionMin(0.03)); + return true; + } + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + + if (AK08.size() > 0) + { + INFO << "Filter the leading AK08 jet:" << endmsg; + const RecJet filteredJet = JetFilter.Execute(AK08[0]); + INFO << "pT(j_filt) = " << filteredJet->pt() << " pT(j_1) = " << AK08[0]->pt() << endmsg; + } + + // RecJets shorthand for `std::vector` + RecJets AK08_filtered = JetFilter.Execute(AK08); + for (MAuint32 i=0; i 0) + { + INFO << "EnergyCorrelator:" << endmsg; + MAdouble64 ec = EC.Execute(AK08[0]); + INFO << "EC = " << ec << endmsg; + } + return true; + } +``` +[Back to top](#outline) +## Nsubjettiness +```c++ + private: + Substructure::Nsubjettiness nsub; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize Nsubjettiness + MAfloat32 beta = 0.1, R0 = 0.2; + nsub.Initialize( + 1., + Substructure::Nsubjettiness::KT_Axes, + Substructure::Nsubjettiness::NormalizedMeasure, + beta, + R0, + ); + return true; + } + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + if (AK08.size() > 0) + { + MAdouble64 tau1 = nsub.Execute(AK08[0]); + INFO << "tau1 = " << tau1 << endmsg; + } + return true; + } +``` +[Back to top](#outline) +## Clustering +```c++ + private: + Substructure::Cluster cluster; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize Cluster + cluster.Initialize(Substructure::cambridge, 0.8, 200., true); + return true; + } + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + if (AK08.size() > 0) + { + RecJets clusteredFromLeading = cluster.Execute(AK08[0]); + RecJets CAJets = cluster.Execute(AK08); + RecJets CAJetsFiltered = cluster.Execute( + AK08[0], [](const RecJet mainjet, const RecJet subjet) + { return (subjet->pt() > mainjet->pt()*0.05);} + ); + } + return true; + } +``` +[Back to top](#outline) +## VariableR +```c++ + private: + Substructure::VariableR varR; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize Energy Correlator + varR.Initialize( + 2000., // mass scale for effective radius (i.e. R ~ rho/pT) + 0., //minimum jet radius + 2., // maximum jet radius + Substructure::VariableR::CALIKE, + Substructure::VariableR::Best, + 200.,// Minimum pT + false // isexclusive + ); + return true; + } + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + if (AK08.size() > 0) + { + RecJets clusteredFromLeading = varR.Execute(AK08[0]); + RecJets varRJets = varR.Execute(AK08); + RecJets varRJetsFiltered = varR.Execute( + AK08[0], [](const RecJet mainjet, const RecJet subjet) + { return (subjet->pt() > mainjet->pt()*0.05);} + ); + } + return true; + } +``` +[Back to top](#outline) +## Pruning +```c++ + private: + Substructure::Pruner pruner; + public: + virtual bool Initialize(const MA5::Configuration& cfg, const std::map& parameters) + { + // Initializing PhysicsService for MC + PHYSICS->mcConfig().Reset(); + AddDefaultHadronic(); + AddDefaultInvisible(); + + // Initializing PhysicsService for RECO + PHYSICS->recConfig().Reset(); + + // Initialize pruner + pruner.Initialize(Substructure::cambridge, 0.2, 0.2, 0.3); + return true; + } + virtual bool Execute(SampleFormat& sample, const EventFormat& event) + { + // Get antikt R=0.8 jets with pT > 200 |eta| < 2.5 + RecJets AK08 = filter(event.rec()->jets("AK08"), 200., 2.5); + + if (AK08.size() > 0) + { + INFO << "Prune the leading AK08 jet:" << endmsg; + const RecJet prunedjet = pruner.Execute(AK08[0]); + INFO << "pT(j_filt) = " << prunedjet->pt() << " pT(j_1) = " << AK08[0]->pt() << endmsg; + } + + // RecJets shorthand for `std::vector` + RecJets AK08_pruned = pruner.Execute(AK08); + for (MAuint32 i=0; i 0) + { + const RecJet reclusteredJet = recluster.Execute(AK08[0]); + RecJets reclusteredJets = recluster.Execute(AK08); + } + return true; + } +``` +[Back to top](#outline) diff --git a/tools/SampleAnalyzer/Interfaces/substructure/Recluster.cpp b/tools/SampleAnalyzer/Interfaces/substructure/Recluster.cpp new file mode 100644 index 00000000..0b8bfbd9 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/Recluster.cpp @@ -0,0 +1,86 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/contrib/Recluster.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/Recluster.h" + +namespace MA5 { + namespace Substructure { + + //=======================// + // Execution // + //=======================// + + // Method to recluster a given jet. Returns only the hardest reclustered jet. + const RecJetFormat* Recluster::Execute(const RecJetFormat *jet) + { + fastjet::contrib::Recluster recluster(*JetDefinition_); + fastjet::PseudoJet reclustered_jet = recluster(jet->pseudojet()); + return __transform_jet(reclustered_jet); + } + + // Method to recluster each jet in a given vector + std::vector Recluster::Execute(std::vector &jets) + { + std::vector output_jets; + for (auto &jet: jets) + output_jets.push_back(Execute(jet)); + + std::sort( + output_jets.begin(), + output_jets.end(), + [](const RecJetFormat *j1, const RecJetFormat *j2) + { + return (j1->pt() > j2->pt()); + } + ); + + return output_jets; + } + + //=============================// + // NOT IMPLEMENTED // + //=============================// + + void Recluster::Execute(const EventFormat& event, std::string JetID) + { + try { throw EXCEPTION_ERROR("This method has not been implemented", "", 1); } + catch (const std::exception& err) { + MANAGE_EXCEPTION(err); + } + } + + template + std::vector Recluster::Execute(const RecJetFormat *jet, Func func) + { + try { throw EXCEPTION_ERROR("This method has not been implemented", "", 1); } + catch (const std::exception& err) { + MANAGE_EXCEPTION(err); + return std::vector(); + } + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Recluster.h b/tools/SampleAnalyzer/Interfaces/substructure/Recluster.h similarity index 64% rename from tools/SampleAnalyzer/Interfaces/fastjet/Recluster.h rename to tools/SampleAnalyzer/Interfaces/substructure/Recluster.h index 07042fc4..2c7e19ab 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/Recluster.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/Recluster.h @@ -24,11 +24,8 @@ #ifndef MADANALYSIS5_RECLUSTER_H #define MADANALYSIS5_RECLUSTER_H -// FastJet headers -#include "fastjet/contrib/Recluster.hh" - // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" +#include "SampleAnalyzer/Interfaces/substructure/ClusterBase.h" using namespace std; @@ -45,7 +42,7 @@ namespace MA5 { Recluster() {} /// Destructor - virtual ~Recluster() {} + ~Recluster() {} //============================// // Initialization // @@ -64,53 +61,19 @@ namespace MA5 { //=======================// // Method to recluster a given jet. Returns only the hardest reclustered jet. - const RecJetFormat* Execute(const RecJetFormat *jet) - { - fastjet::contrib::Recluster recluster(*JetDefinition_); - fastjet::PseudoJet reclustered_jet = recluster(jet->pseudojet()); - return __transform_jet(reclustered_jet); - } + const RecJetFormat* Execute(const RecJetFormat *jet); // Method to recluster each jet in a given vector - std::vector Execute(std::vector &jets) - { - std::vector output_jets; - for (auto &jet: jets) - output_jets.push_back(Execute(jet)); - - std::sort( - output_jets.begin(), - output_jets.end(), - [](const RecJetFormat *j1, const RecJetFormat *j2) - { - return (j1->pt() > j2->pt()); - } - ); - - return output_jets; - } + std::vector Execute(std::vector &jets) override; //=============================// // NOT IMPLEMENTED // //=============================// - void Execute(const EventFormat& event, std::string JetID) - { - try { throw EXCEPTION_ERROR("This method has not been implemented", "", 1); } - catch (const std::exception& err) { - MANAGE_EXCEPTION(err); - } - } + void Execute(const EventFormat& event, std::string JetID) override; template - std::vector Execute(const RecJetFormat *jet, Func func) - { - try { throw EXCEPTION_ERROR("This method has not been implemented", "", 1); } - catch (const std::exception& err) { - MANAGE_EXCEPTION(err); - return std::vector(); - } - } + std::vector Execute(const RecJetFormat *jet, Func func); }; } } diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/Selector.h b/tools/SampleAnalyzer/Interfaces/substructure/Selector.h similarity index 100% rename from tools/SampleAnalyzer/Interfaces/fastjet/Selector.h rename to tools/SampleAnalyzer/Interfaces/substructure/Selector.h diff --git a/tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.cpp b/tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.cpp new file mode 100644 index 00000000..9984af68 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.cpp @@ -0,0 +1,75 @@ +////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +////////////////////////////////////////////////////// + +// FastJet headers +#include "fastjet/contrib/SoftDrop.hh" + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/SoftDrop.h" + +namespace MA5 { + namespace Substructure { + + SoftDrop::~SoftDrop() { delete softDrop_; } + + //============================// + // Initialization // + //============================// + + void SoftDrop::Initialize(MAfloat32 beta, MAfloat32 symmetry_cut, MAfloat32 R0) + { softDrop_ = new fastjet::contrib::SoftDrop(beta, symmetry_cut, R0); } + + //=======================// + // Execution // + //=======================// + + // Execute with a single jet + const RecJetFormat * SoftDrop::Execute(const RecJetFormat *jet) const + { + fastjet::PseudoJet sd_jet = (*softDrop_)(jet->pseudojet()); + RecJetFormat * NewJet = new RecJetFormat(sd_jet); + return NewJet; + } + + // Execute with a list of jets + std::vector SoftDrop::Execute(std::vector &jets) const + { + std::vector output_jets; + output_jets.reserve(jets.size()); + for (auto &jet: jets) + output_jets.push_back(Execute(jet)); + + // Sort with respect to jet pT + std::sort( + output_jets.begin(), + output_jets.end(), + [](const RecJetFormat *j1, const RecJetFormat *j2) + { + return (j1->pt() > j2->pt()); + } + ); + + return output_jets; + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/SoftDrop.h b/tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.h similarity index 77% rename from tools/SampleAnalyzer/Interfaces/fastjet/SoftDrop.h rename to tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.h index 96241ffb..9e6e2880 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/SoftDrop.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/SoftDrop.h @@ -24,13 +24,15 @@ #ifndef MADANALYSIS5_SOFTDROP_H #define MADANALYSIS5_SOFTDROP_H -// FastJet headers -#include "fastjet/contrib/SoftDrop.hh" - // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" +#include "SampleAnalyzer/Commons/Base/PortableDatatypes.h" +#include "SampleAnalyzer/Commons/DataFormat/RecJetFormat.h" -using namespace std; +namespace fastjet { + namespace contrib { + class SoftDrop; + } +} namespace MA5 { namespace Substructure { @@ -71,7 +73,7 @@ namespace MA5 { SoftDrop() {} // Destructor - virtual ~SoftDrop() {} + ~SoftDrop(); //============================// // Initialization // @@ -85,40 +87,17 @@ namespace MA5 { ) { Initialize(beta, symmetry_cut, R0); } - void Initialize(MAfloat32 beta, MAfloat32 symmetry_cut, MAfloat32 R0=1.) - { softDrop_ = new fastjet::contrib::SoftDrop(beta, symmetry_cut, R0); } + void Initialize(MAfloat32 beta, MAfloat32 symmetry_cut, MAfloat32 R0=1.); //=======================// // Execution // //=======================// // Execute with a single jet - const RecJetFormat * Execute(const RecJetFormat *jet) - { - fastjet::PseudoJet sd_jet = (*softDrop_)(jet->pseudojet()); - return ClusterBase().__transform_jet(sd_jet); - } + const RecJetFormat * Execute(const RecJetFormat *jet) const; // Execute with a list of jets - std::vector Execute(std::vector &jets) - { - std::vector output_jets; - for (auto &jet: jets) - output_jets.push_back(Execute(jet)); - - // Sort with respect to jet pT - std::sort( - output_jets.begin(), - output_jets.end(), - [](const RecJetFormat *j1, const RecJetFormat *j2) - { - return (j1->pt() > j2->pt()); - } - ); - - return output_jets; - } - + std::vector Execute(std::vector &jets) const; }; } } diff --git a/tools/SampleAnalyzer/Interfaces/substructure/VariableR.cpp b/tools/SampleAnalyzer/Interfaces/substructure/VariableR.cpp new file mode 100644 index 00000000..a6aeb1e7 --- /dev/null +++ b/tools/SampleAnalyzer/Interfaces/substructure/VariableR.cpp @@ -0,0 +1,80 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +// SampleAnalyser headers +#include "SampleAnalyzer/Interfaces/substructure/VariableR.h" + +// FastJet headers +#include "fastjet/contrib/VariableRPlugin.hh" + +namespace MA5 { + namespace Substructure { + // Initialization method + void VariableR::Initialize( + MAfloat32 rho, + MAfloat32 minR, + MAfloat32 maxR, + Substructure::VariableR::ClusterType clusterType, + Substructure::VariableR::Strategy strategy, + MAfloat32 ptmin, + MAbool isExclusive + ) + { + fastjet::contrib::VariableRPlugin::ClusterType clusterType_ = fastjet::contrib::VariableRPlugin::AKTLIKE; + // whether to use CA-like, kT-like, + // or anti-kT-like distance measure + // (this value is the same as the p exponent in + // generalized-kt, with anti-kt = -1.0, CA = 0.0, and + // kT = 1.0) + if (clusterType == Substructure::VariableR::CALIKE) + clusterType_ = fastjet::contrib::VariableRPlugin::CALIKE; + else if (clusterType == Substructure::VariableR::KTLIKE) + clusterType_ = fastjet::contrib::VariableRPlugin::KTLIKE; + + fastjet::contrib::VariableRPlugin::Strategy strategy_; + // decodes which algorithm to apply for the clustering + + if (strategy == Substructure::VariableR::Best) + strategy_ = fastjet::contrib::VariableRPlugin::Best; + else if (strategy == Substructure::VariableR::N2Tiled) + strategy_ = fastjet::contrib::VariableRPlugin::N2Tiled; + else if (strategy == Substructure::VariableR::N2Plain) + strategy_ = fastjet::contrib::VariableRPlugin::N2Plain; + else if (strategy == Substructure::VariableR::NNH) + strategy_ = fastjet::contrib::VariableRPlugin::NNH; + else if (strategy == Substructure::VariableR::Native) + strategy_ = fastjet::contrib::VariableRPlugin::Native; + + ptmin_ = ptmin; isExclusive_ = isExclusive; + + JetDefPlugin_ = new fastjet::contrib::VariableRPlugin( + rho, minR, maxR, clusterType_, false, strategy_ + ); + isPlugin_ = true; + + /// Note that pre-clustering is deprecated and will likely be + /// removed in a future releasse of this contrib. + /// (precluster = false at the moment) + } + } +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Interfaces/fastjet/VariableR.h b/tools/SampleAnalyzer/Interfaces/substructure/VariableR.h similarity index 55% rename from tools/SampleAnalyzer/Interfaces/fastjet/VariableR.h rename to tools/SampleAnalyzer/Interfaces/substructure/VariableR.h index 20ff8045..c9824b98 100644 --- a/tools/SampleAnalyzer/Interfaces/fastjet/VariableR.h +++ b/tools/SampleAnalyzer/Interfaces/substructure/VariableR.h @@ -24,32 +24,12 @@ #ifndef MADANALYSIS5_VARIABLER_H #define MADANALYSIS5_VARIABLER_H -// FastJet headers -#include "fastjet/contrib/VariableRPlugin.hh" - // SampleAnalyser headers -#include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" - -using namespace std; +#include "SampleAnalyzer/Interfaces/substructure/ClusterBase.h" namespace MA5 { namespace Substructure { class VariableR : public ClusterBase { - - //--------------------------------------------------------------------------------- - // data members - //--------------------------------------------------------------------------------- - protected: - fastjet::contrib::VariableRPlugin::ClusterType clusterType_; // whether to use CA-like, kT-like, - // or anti-kT-like distance measure - // (this value is the same as the p exponent in - // generalized-kt, with anti-kt = -1.0, CA = 0.0, and - // kT = 1.0) - - fastjet::contrib::VariableRPlugin::Strategy strategy_; - // decodes which algorithm to apply for the clustering - - public: enum ClusterType {CALIKE, KTLIKE, AKTLIKE}; @@ -66,7 +46,7 @@ namespace MA5 { VariableR() {} /// Destructor - virtual ~VariableR() {} + ~VariableR() {} //============================// // Initialization // @@ -96,37 +76,7 @@ namespace MA5 { Substructure::VariableR::Strategy strategy = Substructure::VariableR::Best, MAfloat32 ptmin = 0., MAbool isExclusive = false - ) - { - if (clusterType == Substructure::VariableR::CALIKE) - clusterType_ = fastjet::contrib::VariableRPlugin::CALIKE; - else if (clusterType == Substructure::VariableR::KTLIKE) - clusterType_ = fastjet::contrib::VariableRPlugin::KTLIKE; - else if (clusterType == Substructure::VariableR::AKTLIKE) - clusterType_ = fastjet::contrib::VariableRPlugin::AKTLIKE; - - if (strategy == Substructure::VariableR::Best) - strategy_ = fastjet::contrib::VariableRPlugin::Best; - else if (strategy == Substructure::VariableR::N2Tiled) - strategy_ = fastjet::contrib::VariableRPlugin::N2Tiled; - else if (strategy == Substructure::VariableR::N2Plain) - strategy_ = fastjet::contrib::VariableRPlugin::N2Plain; - else if (strategy == Substructure::VariableR::NNH) - strategy_ = fastjet::contrib::VariableRPlugin::NNH; - else if (strategy == Substructure::VariableR::Native) - strategy_ = fastjet::contrib::VariableRPlugin::Native; - - ptmin_ = ptmin; isExclusive_ = isExclusive; - - JetDefPlugin_ = new fastjet::contrib::VariableRPlugin( - rho, minR, maxR, clusterType_, false, strategy_ - ); - isPlugin_ = true; - - /// Note that pre-clustering is deprecated and will likely be - /// removed in a future releasse of this contrib. - /// (precluster = false at the moment) - } + ); }; } } diff --git a/tools/SampleAnalyzer/Process/Analyzer/AnalyzerBase.h b/tools/SampleAnalyzer/Process/Analyzer/AnalyzerBase.h index b44bf29c..38dbf929 100644 --- a/tools/SampleAnalyzer/Process/Analyzer/AnalyzerBase.h +++ b/tools/SampleAnalyzer/Process/Analyzer/AnalyzerBase.h @@ -42,15 +42,16 @@ // FJcontrib tools #ifdef MA5_FASTJET_MODE - #include "SampleAnalyzer/Interfaces/fastjet/SoftDrop.h" - #include "SampleAnalyzer/Interfaces/fastjet/Cluster.h" - #include "SampleAnalyzer/Interfaces/fastjet/Recluster.h" - #include "SampleAnalyzer/Interfaces/fastjet/Nsubjettiness.h" - #include "SampleAnalyzer/Interfaces/fastjet/VariableR.h" - #include "SampleAnalyzer/Interfaces/fastjet/Pruner.h" - #include "SampleAnalyzer/Interfaces/fastjet/Selector.h" - #include "SampleAnalyzer/Interfaces/fastjet/Filter.h" - #include "SampleAnalyzer/Interfaces/fastjet/EnergyCorrelator.h" + #include "SampleAnalyzer/Interfaces/substructure/SoftDrop.h" + #include "SampleAnalyzer/Interfaces/substructure/Cluster.h" + #include "SampleAnalyzer/Interfaces/substructure/Recluster.h" + #include "SampleAnalyzer/Interfaces/substructure/Nsubjettiness.h" + #include "SampleAnalyzer/Interfaces/substructure/VariableR.h" + #include "SampleAnalyzer/Interfaces/substructure/Pruner.h" + #include "SampleAnalyzer/Interfaces/substructure/Selector.h" + #include "SampleAnalyzer/Interfaces/substructure/Filter.h" + #include "SampleAnalyzer/Interfaces/substructure/EnergyCorrelator.h" + #include "SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h" #endif // STL headers @@ -60,7 +61,6 @@ #include #include - // initializing MACRO #define INIT_ANALYSIS(CLASS,NAME) public: CLASS() {setName(NAME);} virtual ~CLASS() {} private: @@ -76,540 +76,538 @@ #define RecTrack MA5::RecTrackFormat * #define RecTracks std::vector -namespace MA5 -{ - -class AnalyzerBase -{ - - // ------------------------------------------------------------- - // data members - // ------------------------------------------------------------- - public : - - /// name of the analysis - std::string name_; - - /// Weighted events mode - MAbool weighted_events_; - - /// A RS manager is associated with each analysis - RegionSelectionManager manager_; - std::string outputdir_; +namespace MA5 { + class AnalyzerBase + { - /// Writer SAF - SAFWriter out_; + // ------------------------------------------------------------- + // data members + // ------------------------------------------------------------- + public : + + /// name of the analysis + std::string name_; + + /// Weighted events mode + MAbool weighted_events_; + + /// A RS manager is associated with each analysis + RegionSelectionManager manager_; + std::string outputdir_; + + /// Writer SAF + SAFWriter out_; - // options - std::map options_; - - // parameters - std::map parameters_; - - // ------------------------------------------------------------- - // method members - // ------------------------------------------------------------- - public : - - /// Constructor without argument - AnalyzerBase() - { name_="unknown"; outputdir_=""; } - - /// Destructor - virtual ~AnalyzerBase() - { } - - /// Initialize (common part to all analyses) - MAbool PreInitialize(const std::string& outputName, - const Configuration* cfg) - { - weighted_events_ = !cfg->IsNoEventWeight(); - out_.Initialize(cfg,outputName.c_str()); - options_ = cfg->Options(); - return true; - } - - /// Initialize (specific to the analysis) - virtual MAbool Initialize(const Configuration& cfg, - const std::map& parameters)=0; - - MAbool Initialize(const Configuration& cfg) - { - return Initialize(cfg, parameters_); - } - - /// PreFinalize - void PreFinalize(const SampleFormat& summary, - const std::vector& samples) - { - out_.WriteHeader(summary); - out_.WriteFiles(samples); - } - - /// Finalize - virtual void Finalize(const SampleFormat& summary, - const std::vector& samples)=0; - - /// PostFinalize - void PostFinalize(const SampleFormat& summary, - const std::vector& samples) - { - // Closing output file - out_.WriteFoot(summary); - out_.Finalize(); - } - - /// Execute - MAbool PreExecute(const SampleFormat& mySample, - const EventFormat& myEvent) - { - PHYSICS->Id->SetFinalState(myEvent.mc()); - PHYSICS->Id->SetInitialState(myEvent.mc()); - return true; - } - - virtual MAbool Execute(SampleFormat& mySample, - const EventFormat& myEvent)=0; - - /// Accessor to analysis name - const std::string name() const {return name_;} - - /// Accessor to the manager - RegionSelectionManager *Manager() { return &manager_; } - - /// Mutator to analysis name - void setName(const std::string& Name) {name_=Name;} - - /// Accessor to the output directory name - const std::string Output() const {return outputdir_;} - - /// Mutator to the output directory name - void SetOutputDir(const std::string &name) {outputdir_=name;} - - - SAFWriter& out() - { return out_; } - - // Set command line options - void SetOptions(std::map options) {options_=options;} - - // Set parameters, initialized in main.cpp - void SetParameters(std::map params) {parameters_=params;} - - // Accessor to parameters - const std::map GetParameters() {return parameters_;} - - /// Get an option for this analysis instance as a string. - std::string getOption(std::string optname) const - { - if ( options_.find(optname) != options_.end() ) - return options_.find(optname)->second; - return ""; - } - - /// Get an option for this analysis instance converted to a specific type - /// The return type is given by the specified @a def value, or by an explicit template - /// type-argument, e.g. getOption("FOO", 3). - template - T getOption(std::string optname, T def) const { - if (options_.find(optname) == options_.end()) return def; - std::stringstream ss; - ss << options_.find(optname)->second; - T ret; - ss >> ret; - return ret; - } - - /// overload for literal character strings (which don't play well with stringstream) - /// Note this isn't a template specialisation, because we can't return a non-static - /// char*, and T-as-return-type is built into the template function definition. - std::string getOption(std::string optname, const char* def) { - return getOption(optname, def); - } - - - void AddDefaultHadronic() - { - // definition of the multiparticle "hadronic" - PHYSICS->mcConfig().AddHadronicId(-20543); - PHYSICS->mcConfig().AddHadronicId(-20533); - PHYSICS->mcConfig().AddHadronicId(-20523); - PHYSICS->mcConfig().AddHadronicId(-20513); - PHYSICS->mcConfig().AddHadronicId(-20433); - PHYSICS->mcConfig().AddHadronicId(-20423); - PHYSICS->mcConfig().AddHadronicId(-20413); - PHYSICS->mcConfig().AddHadronicId(-20323); - PHYSICS->mcConfig().AddHadronicId(-20313); - PHYSICS->mcConfig().AddHadronicId(-20213); - PHYSICS->mcConfig().AddHadronicId(-10543); - PHYSICS->mcConfig().AddHadronicId(-10541); - PHYSICS->mcConfig().AddHadronicId(-10533); - PHYSICS->mcConfig().AddHadronicId(-10531); - PHYSICS->mcConfig().AddHadronicId(-10523); - PHYSICS->mcConfig().AddHadronicId(-10521); - PHYSICS->mcConfig().AddHadronicId(-10513); - PHYSICS->mcConfig().AddHadronicId(-10511); - PHYSICS->mcConfig().AddHadronicId(-10433); - PHYSICS->mcConfig().AddHadronicId(-10431); - PHYSICS->mcConfig().AddHadronicId(-10423); - PHYSICS->mcConfig().AddHadronicId(-10421); - PHYSICS->mcConfig().AddHadronicId(-10413); - PHYSICS->mcConfig().AddHadronicId(-10411); - PHYSICS->mcConfig().AddHadronicId(-10323); - PHYSICS->mcConfig().AddHadronicId(-10321); - PHYSICS->mcConfig().AddHadronicId(-10313); - PHYSICS->mcConfig().AddHadronicId(-10311); - PHYSICS->mcConfig().AddHadronicId(-10213); - PHYSICS->mcConfig().AddHadronicId(-10211); - PHYSICS->mcConfig().AddHadronicId(-5554); - PHYSICS->mcConfig().AddHadronicId(-5544); - PHYSICS->mcConfig().AddHadronicId(-5542); - PHYSICS->mcConfig().AddHadronicId(-5534); - PHYSICS->mcConfig().AddHadronicId(-5532); - PHYSICS->mcConfig().AddHadronicId(-5524); - PHYSICS->mcConfig().AddHadronicId(-5522); - PHYSICS->mcConfig().AddHadronicId(-5514); - PHYSICS->mcConfig().AddHadronicId(-5512); - PHYSICS->mcConfig().AddHadronicId(-5503); - PHYSICS->mcConfig().AddHadronicId(-5444); - PHYSICS->mcConfig().AddHadronicId(-5442); - PHYSICS->mcConfig().AddHadronicId(-5434); - PHYSICS->mcConfig().AddHadronicId(-5432); - PHYSICS->mcConfig().AddHadronicId(-5424); - PHYSICS->mcConfig().AddHadronicId(-5422); - PHYSICS->mcConfig().AddHadronicId(-5414); - PHYSICS->mcConfig().AddHadronicId(-5412); - PHYSICS->mcConfig().AddHadronicId(-5403); - PHYSICS->mcConfig().AddHadronicId(-5401); - PHYSICS->mcConfig().AddHadronicId(-5342); - PHYSICS->mcConfig().AddHadronicId(-5334); - PHYSICS->mcConfig().AddHadronicId(-5332); - PHYSICS->mcConfig().AddHadronicId(-5324); - PHYSICS->mcConfig().AddHadronicId(-5322); - PHYSICS->mcConfig().AddHadronicId(-5314); - PHYSICS->mcConfig().AddHadronicId(-5312); - PHYSICS->mcConfig().AddHadronicId(-5303); - PHYSICS->mcConfig().AddHadronicId(-5301); - PHYSICS->mcConfig().AddHadronicId(-5242); - PHYSICS->mcConfig().AddHadronicId(-5232); - PHYSICS->mcConfig().AddHadronicId(-5224); - PHYSICS->mcConfig().AddHadronicId(-5222); - PHYSICS->mcConfig().AddHadronicId(-5214); - PHYSICS->mcConfig().AddHadronicId(-5212); - PHYSICS->mcConfig().AddHadronicId(-5203); - PHYSICS->mcConfig().AddHadronicId(-5201); - PHYSICS->mcConfig().AddHadronicId(-5142); - PHYSICS->mcConfig().AddHadronicId(-5132); - PHYSICS->mcConfig().AddHadronicId(-5122); - PHYSICS->mcConfig().AddHadronicId(-5114); - PHYSICS->mcConfig().AddHadronicId(-5112); - PHYSICS->mcConfig().AddHadronicId(-5103); - PHYSICS->mcConfig().AddHadronicId(-5101); - PHYSICS->mcConfig().AddHadronicId(-4444); - PHYSICS->mcConfig().AddHadronicId(-4434); - PHYSICS->mcConfig().AddHadronicId(-4432); - PHYSICS->mcConfig().AddHadronicId(-4424); - PHYSICS->mcConfig().AddHadronicId(-4422); - PHYSICS->mcConfig().AddHadronicId(-4414); - PHYSICS->mcConfig().AddHadronicId(-4412); - PHYSICS->mcConfig().AddHadronicId(-4403); - PHYSICS->mcConfig().AddHadronicId(-4334); - PHYSICS->mcConfig().AddHadronicId(-4332); - PHYSICS->mcConfig().AddHadronicId(-4324); - PHYSICS->mcConfig().AddHadronicId(-4322); - PHYSICS->mcConfig().AddHadronicId(-4314); - PHYSICS->mcConfig().AddHadronicId(-4312); - PHYSICS->mcConfig().AddHadronicId(-4303); - PHYSICS->mcConfig().AddHadronicId(-4301); - PHYSICS->mcConfig().AddHadronicId(-4232); - PHYSICS->mcConfig().AddHadronicId(-4224); - PHYSICS->mcConfig().AddHadronicId(-4222); - PHYSICS->mcConfig().AddHadronicId(-4214); - PHYSICS->mcConfig().AddHadronicId(-4212); - PHYSICS->mcConfig().AddHadronicId(-4203); - PHYSICS->mcConfig().AddHadronicId(-4201); - PHYSICS->mcConfig().AddHadronicId(-4132); - PHYSICS->mcConfig().AddHadronicId(-4122); - PHYSICS->mcConfig().AddHadronicId(-4114); - PHYSICS->mcConfig().AddHadronicId(-4112); - PHYSICS->mcConfig().AddHadronicId(-4103); - PHYSICS->mcConfig().AddHadronicId(-4101); - PHYSICS->mcConfig().AddHadronicId(-3334); - PHYSICS->mcConfig().AddHadronicId(-3324); - PHYSICS->mcConfig().AddHadronicId(-3322); - PHYSICS->mcConfig().AddHadronicId(-3314); - PHYSICS->mcConfig().AddHadronicId(-3312); - PHYSICS->mcConfig().AddHadronicId(-3303); - PHYSICS->mcConfig().AddHadronicId(-3224); - PHYSICS->mcConfig().AddHadronicId(-3222); - PHYSICS->mcConfig().AddHadronicId(-3214); - PHYSICS->mcConfig().AddHadronicId(-3212); - PHYSICS->mcConfig().AddHadronicId(-3203); - PHYSICS->mcConfig().AddHadronicId(-3201); - PHYSICS->mcConfig().AddHadronicId(-3122); - PHYSICS->mcConfig().AddHadronicId(-3114); - PHYSICS->mcConfig().AddHadronicId(-3112); - PHYSICS->mcConfig().AddHadronicId(-3103); - PHYSICS->mcConfig().AddHadronicId(-3101); - PHYSICS->mcConfig().AddHadronicId(-2224); - PHYSICS->mcConfig().AddHadronicId(-2214); - PHYSICS->mcConfig().AddHadronicId(-2212); - PHYSICS->mcConfig().AddHadronicId(-2203); - PHYSICS->mcConfig().AddHadronicId(-2114); - PHYSICS->mcConfig().AddHadronicId(-2112); - PHYSICS->mcConfig().AddHadronicId(-2103); - PHYSICS->mcConfig().AddHadronicId(-2101); - PHYSICS->mcConfig().AddHadronicId(-1114); - PHYSICS->mcConfig().AddHadronicId(-1103); - PHYSICS->mcConfig().AddHadronicId(-545); - PHYSICS->mcConfig().AddHadronicId(-543); - PHYSICS->mcConfig().AddHadronicId(-541); - PHYSICS->mcConfig().AddHadronicId(-535); - PHYSICS->mcConfig().AddHadronicId(-533); - PHYSICS->mcConfig().AddHadronicId(-531); - PHYSICS->mcConfig().AddHadronicId(-525); - PHYSICS->mcConfig().AddHadronicId(-523); - PHYSICS->mcConfig().AddHadronicId(-521); - PHYSICS->mcConfig().AddHadronicId(-515); - PHYSICS->mcConfig().AddHadronicId(-513); - PHYSICS->mcConfig().AddHadronicId(-511); - PHYSICS->mcConfig().AddHadronicId(-435); - PHYSICS->mcConfig().AddHadronicId(-433); - PHYSICS->mcConfig().AddHadronicId(-431); - PHYSICS->mcConfig().AddHadronicId(-425); - PHYSICS->mcConfig().AddHadronicId(-423); - PHYSICS->mcConfig().AddHadronicId(-421); - PHYSICS->mcConfig().AddHadronicId(-415); - PHYSICS->mcConfig().AddHadronicId(-413); - PHYSICS->mcConfig().AddHadronicId(-411); - PHYSICS->mcConfig().AddHadronicId(-325); - PHYSICS->mcConfig().AddHadronicId(-323); - PHYSICS->mcConfig().AddHadronicId(-321); - PHYSICS->mcConfig().AddHadronicId(-315); - PHYSICS->mcConfig().AddHadronicId(-313); - PHYSICS->mcConfig().AddHadronicId(-311); - PHYSICS->mcConfig().AddHadronicId(-215); - PHYSICS->mcConfig().AddHadronicId(-213); - PHYSICS->mcConfig().AddHadronicId(-211); - PHYSICS->mcConfig().AddHadronicId(111); - PHYSICS->mcConfig().AddHadronicId(113); - PHYSICS->mcConfig().AddHadronicId(115); - PHYSICS->mcConfig().AddHadronicId(130); - PHYSICS->mcConfig().AddHadronicId(211); - PHYSICS->mcConfig().AddHadronicId(213); - PHYSICS->mcConfig().AddHadronicId(215); - PHYSICS->mcConfig().AddHadronicId(221); - PHYSICS->mcConfig().AddHadronicId(223); - PHYSICS->mcConfig().AddHadronicId(225); - PHYSICS->mcConfig().AddHadronicId(310); - PHYSICS->mcConfig().AddHadronicId(311); - PHYSICS->mcConfig().AddHadronicId(313); - PHYSICS->mcConfig().AddHadronicId(315); - PHYSICS->mcConfig().AddHadronicId(321); - PHYSICS->mcConfig().AddHadronicId(323); - PHYSICS->mcConfig().AddHadronicId(325); - PHYSICS->mcConfig().AddHadronicId(331); - PHYSICS->mcConfig().AddHadronicId(333); - PHYSICS->mcConfig().AddHadronicId(335); - PHYSICS->mcConfig().AddHadronicId(411); - PHYSICS->mcConfig().AddHadronicId(413); - PHYSICS->mcConfig().AddHadronicId(415); - PHYSICS->mcConfig().AddHadronicId(421); - PHYSICS->mcConfig().AddHadronicId(423); - PHYSICS->mcConfig().AddHadronicId(425); - PHYSICS->mcConfig().AddHadronicId(431); - PHYSICS->mcConfig().AddHadronicId(433); - PHYSICS->mcConfig().AddHadronicId(435); - PHYSICS->mcConfig().AddHadronicId(441); - PHYSICS->mcConfig().AddHadronicId(443); - PHYSICS->mcConfig().AddHadronicId(445); - PHYSICS->mcConfig().AddHadronicId(511); - PHYSICS->mcConfig().AddHadronicId(513); - PHYSICS->mcConfig().AddHadronicId(515); - PHYSICS->mcConfig().AddHadronicId(521); - PHYSICS->mcConfig().AddHadronicId(523); - PHYSICS->mcConfig().AddHadronicId(525); - PHYSICS->mcConfig().AddHadronicId(531); - PHYSICS->mcConfig().AddHadronicId(533); - PHYSICS->mcConfig().AddHadronicId(535); - PHYSICS->mcConfig().AddHadronicId(541); - PHYSICS->mcConfig().AddHadronicId(543); - PHYSICS->mcConfig().AddHadronicId(545); - PHYSICS->mcConfig().AddHadronicId(551); - PHYSICS->mcConfig().AddHadronicId(553); - PHYSICS->mcConfig().AddHadronicId(555); - PHYSICS->mcConfig().AddHadronicId(1103); - PHYSICS->mcConfig().AddHadronicId(1114); - PHYSICS->mcConfig().AddHadronicId(2101); - PHYSICS->mcConfig().AddHadronicId(2103); - PHYSICS->mcConfig().AddHadronicId(2112); - PHYSICS->mcConfig().AddHadronicId(2114); - PHYSICS->mcConfig().AddHadronicId(2203); - PHYSICS->mcConfig().AddHadronicId(2212); - PHYSICS->mcConfig().AddHadronicId(2214); - PHYSICS->mcConfig().AddHadronicId(2224); - PHYSICS->mcConfig().AddHadronicId(3101); - PHYSICS->mcConfig().AddHadronicId(3103); - PHYSICS->mcConfig().AddHadronicId(3112); - PHYSICS->mcConfig().AddHadronicId(3114); - PHYSICS->mcConfig().AddHadronicId(3122); - PHYSICS->mcConfig().AddHadronicId(3201); - PHYSICS->mcConfig().AddHadronicId(3203); - PHYSICS->mcConfig().AddHadronicId(3212); - PHYSICS->mcConfig().AddHadronicId(3214); - PHYSICS->mcConfig().AddHadronicId(3222); - PHYSICS->mcConfig().AddHadronicId(3224); - PHYSICS->mcConfig().AddHadronicId(3303); - PHYSICS->mcConfig().AddHadronicId(3312); - PHYSICS->mcConfig().AddHadronicId(3314); - PHYSICS->mcConfig().AddHadronicId(3322); - PHYSICS->mcConfig().AddHadronicId(3324); - PHYSICS->mcConfig().AddHadronicId(3334); - PHYSICS->mcConfig().AddHadronicId(4101); - PHYSICS->mcConfig().AddHadronicId(4103); - PHYSICS->mcConfig().AddHadronicId(4112); - PHYSICS->mcConfig().AddHadronicId(4114); - PHYSICS->mcConfig().AddHadronicId(4122); - PHYSICS->mcConfig().AddHadronicId(4132); - PHYSICS->mcConfig().AddHadronicId(4201); - PHYSICS->mcConfig().AddHadronicId(4203); - PHYSICS->mcConfig().AddHadronicId(4212); - PHYSICS->mcConfig().AddHadronicId(4214); - PHYSICS->mcConfig().AddHadronicId(4222); - PHYSICS->mcConfig().AddHadronicId(4224); - PHYSICS->mcConfig().AddHadronicId(4232); - PHYSICS->mcConfig().AddHadronicId(4301); - PHYSICS->mcConfig().AddHadronicId(4303); - PHYSICS->mcConfig().AddHadronicId(4312); - PHYSICS->mcConfig().AddHadronicId(4314); - PHYSICS->mcConfig().AddHadronicId(4322); - PHYSICS->mcConfig().AddHadronicId(4324); - PHYSICS->mcConfig().AddHadronicId(4332); - PHYSICS->mcConfig().AddHadronicId(4334); - PHYSICS->mcConfig().AddHadronicId(4403); - PHYSICS->mcConfig().AddHadronicId(4412); - PHYSICS->mcConfig().AddHadronicId(4414); - PHYSICS->mcConfig().AddHadronicId(4422); - PHYSICS->mcConfig().AddHadronicId(4424); - PHYSICS->mcConfig().AddHadronicId(4432); - PHYSICS->mcConfig().AddHadronicId(4434); - PHYSICS->mcConfig().AddHadronicId(4444); - PHYSICS->mcConfig().AddHadronicId(5101); - PHYSICS->mcConfig().AddHadronicId(5103); - PHYSICS->mcConfig().AddHadronicId(5112); - PHYSICS->mcConfig().AddHadronicId(5114); - PHYSICS->mcConfig().AddHadronicId(5122); - PHYSICS->mcConfig().AddHadronicId(5132); - PHYSICS->mcConfig().AddHadronicId(5142); - PHYSICS->mcConfig().AddHadronicId(5201); - PHYSICS->mcConfig().AddHadronicId(5203); - PHYSICS->mcConfig().AddHadronicId(5212); - PHYSICS->mcConfig().AddHadronicId(5214); - PHYSICS->mcConfig().AddHadronicId(5222); - PHYSICS->mcConfig().AddHadronicId(5224); - PHYSICS->mcConfig().AddHadronicId(5232); - PHYSICS->mcConfig().AddHadronicId(5242); - PHYSICS->mcConfig().AddHadronicId(5301); - PHYSICS->mcConfig().AddHadronicId(5303); - PHYSICS->mcConfig().AddHadronicId(5312); - PHYSICS->mcConfig().AddHadronicId(5314); - PHYSICS->mcConfig().AddHadronicId(5322); - PHYSICS->mcConfig().AddHadronicId(5324); - PHYSICS->mcConfig().AddHadronicId(5332); - PHYSICS->mcConfig().AddHadronicId(5334); - PHYSICS->mcConfig().AddHadronicId(5342); - PHYSICS->mcConfig().AddHadronicId(5401); - PHYSICS->mcConfig().AddHadronicId(5403); - PHYSICS->mcConfig().AddHadronicId(5412); - PHYSICS->mcConfig().AddHadronicId(5414); - PHYSICS->mcConfig().AddHadronicId(5422); - PHYSICS->mcConfig().AddHadronicId(5424); - PHYSICS->mcConfig().AddHadronicId(5432); - PHYSICS->mcConfig().AddHadronicId(5434); - PHYSICS->mcConfig().AddHadronicId(5442); - PHYSICS->mcConfig().AddHadronicId(5444); - PHYSICS->mcConfig().AddHadronicId(5503); - PHYSICS->mcConfig().AddHadronicId(5512); - PHYSICS->mcConfig().AddHadronicId(5514); - PHYSICS->mcConfig().AddHadronicId(5522); - PHYSICS->mcConfig().AddHadronicId(5524); - PHYSICS->mcConfig().AddHadronicId(5532); - PHYSICS->mcConfig().AddHadronicId(5534); - PHYSICS->mcConfig().AddHadronicId(5542); - PHYSICS->mcConfig().AddHadronicId(5544); - PHYSICS->mcConfig().AddHadronicId(5554); - PHYSICS->mcConfig().AddHadronicId(10111); - PHYSICS->mcConfig().AddHadronicId(10113); - PHYSICS->mcConfig().AddHadronicId(10211); - PHYSICS->mcConfig().AddHadronicId(10213); - PHYSICS->mcConfig().AddHadronicId(10221); - PHYSICS->mcConfig().AddHadronicId(10223); - PHYSICS->mcConfig().AddHadronicId(10311); - PHYSICS->mcConfig().AddHadronicId(10313); - PHYSICS->mcConfig().AddHadronicId(10321); - PHYSICS->mcConfig().AddHadronicId(10323); - PHYSICS->mcConfig().AddHadronicId(10331); - PHYSICS->mcConfig().AddHadronicId(10333); - PHYSICS->mcConfig().AddHadronicId(10411); - PHYSICS->mcConfig().AddHadronicId(10413); - PHYSICS->mcConfig().AddHadronicId(10421); - PHYSICS->mcConfig().AddHadronicId(10423); - PHYSICS->mcConfig().AddHadronicId(10431); - PHYSICS->mcConfig().AddHadronicId(10433); - PHYSICS->mcConfig().AddHadronicId(10441); - PHYSICS->mcConfig().AddHadronicId(10443); - PHYSICS->mcConfig().AddHadronicId(10511); - PHYSICS->mcConfig().AddHadronicId(10513); - PHYSICS->mcConfig().AddHadronicId(10521); - PHYSICS->mcConfig().AddHadronicId(10523); - PHYSICS->mcConfig().AddHadronicId(10531); - PHYSICS->mcConfig().AddHadronicId(10533); - PHYSICS->mcConfig().AddHadronicId(10541); - PHYSICS->mcConfig().AddHadronicId(10543); - PHYSICS->mcConfig().AddHadronicId(10551); - PHYSICS->mcConfig().AddHadronicId(10553); - PHYSICS->mcConfig().AddHadronicId(20113); - PHYSICS->mcConfig().AddHadronicId(20213); - PHYSICS->mcConfig().AddHadronicId(20223); - PHYSICS->mcConfig().AddHadronicId(20313); - PHYSICS->mcConfig().AddHadronicId(20323); - PHYSICS->mcConfig().AddHadronicId(20333); - PHYSICS->mcConfig().AddHadronicId(20413); - PHYSICS->mcConfig().AddHadronicId(20423); - PHYSICS->mcConfig().AddHadronicId(20433); - PHYSICS->mcConfig().AddHadronicId(20443); - PHYSICS->mcConfig().AddHadronicId(20513); - PHYSICS->mcConfig().AddHadronicId(20523); - PHYSICS->mcConfig().AddHadronicId(20533); - PHYSICS->mcConfig().AddHadronicId(20543); - PHYSICS->mcConfig().AddHadronicId(20553); - PHYSICS->mcConfig().AddHadronicId(100443); - PHYSICS->mcConfig().AddHadronicId(100553); - PHYSICS->mcConfig().AddHadronicId(9900440); - PHYSICS->mcConfig().AddHadronicId(9900441); - PHYSICS->mcConfig().AddHadronicId(9900443); - PHYSICS->mcConfig().AddHadronicId(9900551); - PHYSICS->mcConfig().AddHadronicId(9900553); - PHYSICS->mcConfig().AddHadronicId(9910441); - PHYSICS->mcConfig().AddHadronicId(9910551); - } - - void AddDefaultInvisible() - { - // definition of the multiparticle "invisible" - PHYSICS->mcConfig().AddInvisibleId(-16); - PHYSICS->mcConfig().AddInvisibleId(-14); - PHYSICS->mcConfig().AddInvisibleId(-12); - PHYSICS->mcConfig().AddInvisibleId(12); - PHYSICS->mcConfig().AddInvisibleId(14); - PHYSICS->mcConfig().AddInvisibleId(16); - PHYSICS->mcConfig().AddInvisibleId(1000022); - PHYSICS->mcConfig().AddInvisibleId(1000039); - } - - protected : - -}; + // options + std::map options_; + + // parameters + std::map parameters_; + + // ------------------------------------------------------------- + // method members + // ------------------------------------------------------------- + public : + + /// Constructor without argument + AnalyzerBase() + { name_="unknown"; outputdir_=""; } + + /// Destructor + virtual ~AnalyzerBase() + { } + + /// Initialize (common part to all analyses) + MAbool PreInitialize(const std::string& outputName, + const Configuration* cfg) + { + weighted_events_ = !cfg->IsNoEventWeight(); + out_.Initialize(cfg,outputName.c_str()); + options_ = cfg->Options(); + return true; + } + + /// Initialize (specific to the analysis) + virtual MAbool Initialize(const Configuration& cfg, + const std::map& parameters)=0; + + MAbool Initialize(const Configuration& cfg) + { + return Initialize(cfg, parameters_); + } + + /// PreFinalize + void PreFinalize(const SampleFormat& summary, + const std::vector& samples) + { + out_.WriteHeader(summary); + out_.WriteFiles(samples); + } + + /// Finalize + virtual void Finalize(const SampleFormat& summary, + const std::vector& samples)=0; + + /// PostFinalize + void PostFinalize(const SampleFormat& summary, + const std::vector& samples) + { + // Closing output file + out_.WriteFoot(summary); + out_.Finalize(); + } + + /// Execute + MAbool PreExecute(const SampleFormat& mySample, + const EventFormat& myEvent) + { + PHYSICS->Id->SetFinalState(myEvent.mc()); + PHYSICS->Id->SetInitialState(myEvent.mc()); + return true; + } + + virtual MAbool Execute(SampleFormat& mySample, + const EventFormat& myEvent)=0; + + /// Accessor to analysis name + const std::string name() const {return name_;} + + /// Accessor to the manager + RegionSelectionManager *Manager() { return &manager_; } + + /// Mutator to analysis name + void setName(const std::string& Name) {name_=Name;} + + /// Accessor to the output directory name + const std::string Output() const {return outputdir_;} + + /// Mutator to the output directory name + void SetOutputDir(const std::string &name) {outputdir_=name;} + + + SAFWriter& out() + { return out_; } + + // Set command line options + void SetOptions(std::map options) {options_=options;} + + // Set parameters, initialized in main.cpp + void SetParameters(std::map params) {parameters_=params;} + + // Accessor to parameters + const std::map GetParameters() {return parameters_;} + + /// Get an option for this analysis instance as a string. + std::string getOption(std::string optname) const + { + if ( options_.find(optname) != options_.end() ) + return options_.find(optname)->second; + return ""; + } + + /// Get an option for this analysis instance converted to a specific type + /// The return type is given by the specified @a def value, or by an explicit template + /// type-argument, e.g. getOption("FOO", 3). + template + T getOption(std::string optname, T def) const { + if (options_.find(optname) == options_.end()) return def; + std::stringstream ss; + ss << options_.find(optname)->second; + T ret; + ss >> ret; + return ret; + } + + /// overload for literal character strings (which don't play well with stringstream) + /// Note this isn't a template specialisation, because we can't return a non-static + /// char*, and T-as-return-type is built into the template function definition. + std::string getOption(std::string optname, const char* def) { + return getOption(optname, def); + } + + + void AddDefaultHadronic() + { + // definition of the multiparticle "hadronic" + PHYSICS->mcConfig().AddHadronicId(-20543); + PHYSICS->mcConfig().AddHadronicId(-20533); + PHYSICS->mcConfig().AddHadronicId(-20523); + PHYSICS->mcConfig().AddHadronicId(-20513); + PHYSICS->mcConfig().AddHadronicId(-20433); + PHYSICS->mcConfig().AddHadronicId(-20423); + PHYSICS->mcConfig().AddHadronicId(-20413); + PHYSICS->mcConfig().AddHadronicId(-20323); + PHYSICS->mcConfig().AddHadronicId(-20313); + PHYSICS->mcConfig().AddHadronicId(-20213); + PHYSICS->mcConfig().AddHadronicId(-10543); + PHYSICS->mcConfig().AddHadronicId(-10541); + PHYSICS->mcConfig().AddHadronicId(-10533); + PHYSICS->mcConfig().AddHadronicId(-10531); + PHYSICS->mcConfig().AddHadronicId(-10523); + PHYSICS->mcConfig().AddHadronicId(-10521); + PHYSICS->mcConfig().AddHadronicId(-10513); + PHYSICS->mcConfig().AddHadronicId(-10511); + PHYSICS->mcConfig().AddHadronicId(-10433); + PHYSICS->mcConfig().AddHadronicId(-10431); + PHYSICS->mcConfig().AddHadronicId(-10423); + PHYSICS->mcConfig().AddHadronicId(-10421); + PHYSICS->mcConfig().AddHadronicId(-10413); + PHYSICS->mcConfig().AddHadronicId(-10411); + PHYSICS->mcConfig().AddHadronicId(-10323); + PHYSICS->mcConfig().AddHadronicId(-10321); + PHYSICS->mcConfig().AddHadronicId(-10313); + PHYSICS->mcConfig().AddHadronicId(-10311); + PHYSICS->mcConfig().AddHadronicId(-10213); + PHYSICS->mcConfig().AddHadronicId(-10211); + PHYSICS->mcConfig().AddHadronicId(-5554); + PHYSICS->mcConfig().AddHadronicId(-5544); + PHYSICS->mcConfig().AddHadronicId(-5542); + PHYSICS->mcConfig().AddHadronicId(-5534); + PHYSICS->mcConfig().AddHadronicId(-5532); + PHYSICS->mcConfig().AddHadronicId(-5524); + PHYSICS->mcConfig().AddHadronicId(-5522); + PHYSICS->mcConfig().AddHadronicId(-5514); + PHYSICS->mcConfig().AddHadronicId(-5512); + PHYSICS->mcConfig().AddHadronicId(-5503); + PHYSICS->mcConfig().AddHadronicId(-5444); + PHYSICS->mcConfig().AddHadronicId(-5442); + PHYSICS->mcConfig().AddHadronicId(-5434); + PHYSICS->mcConfig().AddHadronicId(-5432); + PHYSICS->mcConfig().AddHadronicId(-5424); + PHYSICS->mcConfig().AddHadronicId(-5422); + PHYSICS->mcConfig().AddHadronicId(-5414); + PHYSICS->mcConfig().AddHadronicId(-5412); + PHYSICS->mcConfig().AddHadronicId(-5403); + PHYSICS->mcConfig().AddHadronicId(-5401); + PHYSICS->mcConfig().AddHadronicId(-5342); + PHYSICS->mcConfig().AddHadronicId(-5334); + PHYSICS->mcConfig().AddHadronicId(-5332); + PHYSICS->mcConfig().AddHadronicId(-5324); + PHYSICS->mcConfig().AddHadronicId(-5322); + PHYSICS->mcConfig().AddHadronicId(-5314); + PHYSICS->mcConfig().AddHadronicId(-5312); + PHYSICS->mcConfig().AddHadronicId(-5303); + PHYSICS->mcConfig().AddHadronicId(-5301); + PHYSICS->mcConfig().AddHadronicId(-5242); + PHYSICS->mcConfig().AddHadronicId(-5232); + PHYSICS->mcConfig().AddHadronicId(-5224); + PHYSICS->mcConfig().AddHadronicId(-5222); + PHYSICS->mcConfig().AddHadronicId(-5214); + PHYSICS->mcConfig().AddHadronicId(-5212); + PHYSICS->mcConfig().AddHadronicId(-5203); + PHYSICS->mcConfig().AddHadronicId(-5201); + PHYSICS->mcConfig().AddHadronicId(-5142); + PHYSICS->mcConfig().AddHadronicId(-5132); + PHYSICS->mcConfig().AddHadronicId(-5122); + PHYSICS->mcConfig().AddHadronicId(-5114); + PHYSICS->mcConfig().AddHadronicId(-5112); + PHYSICS->mcConfig().AddHadronicId(-5103); + PHYSICS->mcConfig().AddHadronicId(-5101); + PHYSICS->mcConfig().AddHadronicId(-4444); + PHYSICS->mcConfig().AddHadronicId(-4434); + PHYSICS->mcConfig().AddHadronicId(-4432); + PHYSICS->mcConfig().AddHadronicId(-4424); + PHYSICS->mcConfig().AddHadronicId(-4422); + PHYSICS->mcConfig().AddHadronicId(-4414); + PHYSICS->mcConfig().AddHadronicId(-4412); + PHYSICS->mcConfig().AddHadronicId(-4403); + PHYSICS->mcConfig().AddHadronicId(-4334); + PHYSICS->mcConfig().AddHadronicId(-4332); + PHYSICS->mcConfig().AddHadronicId(-4324); + PHYSICS->mcConfig().AddHadronicId(-4322); + PHYSICS->mcConfig().AddHadronicId(-4314); + PHYSICS->mcConfig().AddHadronicId(-4312); + PHYSICS->mcConfig().AddHadronicId(-4303); + PHYSICS->mcConfig().AddHadronicId(-4301); + PHYSICS->mcConfig().AddHadronicId(-4232); + PHYSICS->mcConfig().AddHadronicId(-4224); + PHYSICS->mcConfig().AddHadronicId(-4222); + PHYSICS->mcConfig().AddHadronicId(-4214); + PHYSICS->mcConfig().AddHadronicId(-4212); + PHYSICS->mcConfig().AddHadronicId(-4203); + PHYSICS->mcConfig().AddHadronicId(-4201); + PHYSICS->mcConfig().AddHadronicId(-4132); + PHYSICS->mcConfig().AddHadronicId(-4122); + PHYSICS->mcConfig().AddHadronicId(-4114); + PHYSICS->mcConfig().AddHadronicId(-4112); + PHYSICS->mcConfig().AddHadronicId(-4103); + PHYSICS->mcConfig().AddHadronicId(-4101); + PHYSICS->mcConfig().AddHadronicId(-3334); + PHYSICS->mcConfig().AddHadronicId(-3324); + PHYSICS->mcConfig().AddHadronicId(-3322); + PHYSICS->mcConfig().AddHadronicId(-3314); + PHYSICS->mcConfig().AddHadronicId(-3312); + PHYSICS->mcConfig().AddHadronicId(-3303); + PHYSICS->mcConfig().AddHadronicId(-3224); + PHYSICS->mcConfig().AddHadronicId(-3222); + PHYSICS->mcConfig().AddHadronicId(-3214); + PHYSICS->mcConfig().AddHadronicId(-3212); + PHYSICS->mcConfig().AddHadronicId(-3203); + PHYSICS->mcConfig().AddHadronicId(-3201); + PHYSICS->mcConfig().AddHadronicId(-3122); + PHYSICS->mcConfig().AddHadronicId(-3114); + PHYSICS->mcConfig().AddHadronicId(-3112); + PHYSICS->mcConfig().AddHadronicId(-3103); + PHYSICS->mcConfig().AddHadronicId(-3101); + PHYSICS->mcConfig().AddHadronicId(-2224); + PHYSICS->mcConfig().AddHadronicId(-2214); + PHYSICS->mcConfig().AddHadronicId(-2212); + PHYSICS->mcConfig().AddHadronicId(-2203); + PHYSICS->mcConfig().AddHadronicId(-2114); + PHYSICS->mcConfig().AddHadronicId(-2112); + PHYSICS->mcConfig().AddHadronicId(-2103); + PHYSICS->mcConfig().AddHadronicId(-2101); + PHYSICS->mcConfig().AddHadronicId(-1114); + PHYSICS->mcConfig().AddHadronicId(-1103); + PHYSICS->mcConfig().AddHadronicId(-545); + PHYSICS->mcConfig().AddHadronicId(-543); + PHYSICS->mcConfig().AddHadronicId(-541); + PHYSICS->mcConfig().AddHadronicId(-535); + PHYSICS->mcConfig().AddHadronicId(-533); + PHYSICS->mcConfig().AddHadronicId(-531); + PHYSICS->mcConfig().AddHadronicId(-525); + PHYSICS->mcConfig().AddHadronicId(-523); + PHYSICS->mcConfig().AddHadronicId(-521); + PHYSICS->mcConfig().AddHadronicId(-515); + PHYSICS->mcConfig().AddHadronicId(-513); + PHYSICS->mcConfig().AddHadronicId(-511); + PHYSICS->mcConfig().AddHadronicId(-435); + PHYSICS->mcConfig().AddHadronicId(-433); + PHYSICS->mcConfig().AddHadronicId(-431); + PHYSICS->mcConfig().AddHadronicId(-425); + PHYSICS->mcConfig().AddHadronicId(-423); + PHYSICS->mcConfig().AddHadronicId(-421); + PHYSICS->mcConfig().AddHadronicId(-415); + PHYSICS->mcConfig().AddHadronicId(-413); + PHYSICS->mcConfig().AddHadronicId(-411); + PHYSICS->mcConfig().AddHadronicId(-325); + PHYSICS->mcConfig().AddHadronicId(-323); + PHYSICS->mcConfig().AddHadronicId(-321); + PHYSICS->mcConfig().AddHadronicId(-315); + PHYSICS->mcConfig().AddHadronicId(-313); + PHYSICS->mcConfig().AddHadronicId(-311); + PHYSICS->mcConfig().AddHadronicId(-215); + PHYSICS->mcConfig().AddHadronicId(-213); + PHYSICS->mcConfig().AddHadronicId(-211); + PHYSICS->mcConfig().AddHadronicId(111); + PHYSICS->mcConfig().AddHadronicId(113); + PHYSICS->mcConfig().AddHadronicId(115); + PHYSICS->mcConfig().AddHadronicId(130); + PHYSICS->mcConfig().AddHadronicId(211); + PHYSICS->mcConfig().AddHadronicId(213); + PHYSICS->mcConfig().AddHadronicId(215); + PHYSICS->mcConfig().AddHadronicId(221); + PHYSICS->mcConfig().AddHadronicId(223); + PHYSICS->mcConfig().AddHadronicId(225); + PHYSICS->mcConfig().AddHadronicId(310); + PHYSICS->mcConfig().AddHadronicId(311); + PHYSICS->mcConfig().AddHadronicId(313); + PHYSICS->mcConfig().AddHadronicId(315); + PHYSICS->mcConfig().AddHadronicId(321); + PHYSICS->mcConfig().AddHadronicId(323); + PHYSICS->mcConfig().AddHadronicId(325); + PHYSICS->mcConfig().AddHadronicId(331); + PHYSICS->mcConfig().AddHadronicId(333); + PHYSICS->mcConfig().AddHadronicId(335); + PHYSICS->mcConfig().AddHadronicId(411); + PHYSICS->mcConfig().AddHadronicId(413); + PHYSICS->mcConfig().AddHadronicId(415); + PHYSICS->mcConfig().AddHadronicId(421); + PHYSICS->mcConfig().AddHadronicId(423); + PHYSICS->mcConfig().AddHadronicId(425); + PHYSICS->mcConfig().AddHadronicId(431); + PHYSICS->mcConfig().AddHadronicId(433); + PHYSICS->mcConfig().AddHadronicId(435); + PHYSICS->mcConfig().AddHadronicId(441); + PHYSICS->mcConfig().AddHadronicId(443); + PHYSICS->mcConfig().AddHadronicId(445); + PHYSICS->mcConfig().AddHadronicId(511); + PHYSICS->mcConfig().AddHadronicId(513); + PHYSICS->mcConfig().AddHadronicId(515); + PHYSICS->mcConfig().AddHadronicId(521); + PHYSICS->mcConfig().AddHadronicId(523); + PHYSICS->mcConfig().AddHadronicId(525); + PHYSICS->mcConfig().AddHadronicId(531); + PHYSICS->mcConfig().AddHadronicId(533); + PHYSICS->mcConfig().AddHadronicId(535); + PHYSICS->mcConfig().AddHadronicId(541); + PHYSICS->mcConfig().AddHadronicId(543); + PHYSICS->mcConfig().AddHadronicId(545); + PHYSICS->mcConfig().AddHadronicId(551); + PHYSICS->mcConfig().AddHadronicId(553); + PHYSICS->mcConfig().AddHadronicId(555); + PHYSICS->mcConfig().AddHadronicId(1103); + PHYSICS->mcConfig().AddHadronicId(1114); + PHYSICS->mcConfig().AddHadronicId(2101); + PHYSICS->mcConfig().AddHadronicId(2103); + PHYSICS->mcConfig().AddHadronicId(2112); + PHYSICS->mcConfig().AddHadronicId(2114); + PHYSICS->mcConfig().AddHadronicId(2203); + PHYSICS->mcConfig().AddHadronicId(2212); + PHYSICS->mcConfig().AddHadronicId(2214); + PHYSICS->mcConfig().AddHadronicId(2224); + PHYSICS->mcConfig().AddHadronicId(3101); + PHYSICS->mcConfig().AddHadronicId(3103); + PHYSICS->mcConfig().AddHadronicId(3112); + PHYSICS->mcConfig().AddHadronicId(3114); + PHYSICS->mcConfig().AddHadronicId(3122); + PHYSICS->mcConfig().AddHadronicId(3201); + PHYSICS->mcConfig().AddHadronicId(3203); + PHYSICS->mcConfig().AddHadronicId(3212); + PHYSICS->mcConfig().AddHadronicId(3214); + PHYSICS->mcConfig().AddHadronicId(3222); + PHYSICS->mcConfig().AddHadronicId(3224); + PHYSICS->mcConfig().AddHadronicId(3303); + PHYSICS->mcConfig().AddHadronicId(3312); + PHYSICS->mcConfig().AddHadronicId(3314); + PHYSICS->mcConfig().AddHadronicId(3322); + PHYSICS->mcConfig().AddHadronicId(3324); + PHYSICS->mcConfig().AddHadronicId(3334); + PHYSICS->mcConfig().AddHadronicId(4101); + PHYSICS->mcConfig().AddHadronicId(4103); + PHYSICS->mcConfig().AddHadronicId(4112); + PHYSICS->mcConfig().AddHadronicId(4114); + PHYSICS->mcConfig().AddHadronicId(4122); + PHYSICS->mcConfig().AddHadronicId(4132); + PHYSICS->mcConfig().AddHadronicId(4201); + PHYSICS->mcConfig().AddHadronicId(4203); + PHYSICS->mcConfig().AddHadronicId(4212); + PHYSICS->mcConfig().AddHadronicId(4214); + PHYSICS->mcConfig().AddHadronicId(4222); + PHYSICS->mcConfig().AddHadronicId(4224); + PHYSICS->mcConfig().AddHadronicId(4232); + PHYSICS->mcConfig().AddHadronicId(4301); + PHYSICS->mcConfig().AddHadronicId(4303); + PHYSICS->mcConfig().AddHadronicId(4312); + PHYSICS->mcConfig().AddHadronicId(4314); + PHYSICS->mcConfig().AddHadronicId(4322); + PHYSICS->mcConfig().AddHadronicId(4324); + PHYSICS->mcConfig().AddHadronicId(4332); + PHYSICS->mcConfig().AddHadronicId(4334); + PHYSICS->mcConfig().AddHadronicId(4403); + PHYSICS->mcConfig().AddHadronicId(4412); + PHYSICS->mcConfig().AddHadronicId(4414); + PHYSICS->mcConfig().AddHadronicId(4422); + PHYSICS->mcConfig().AddHadronicId(4424); + PHYSICS->mcConfig().AddHadronicId(4432); + PHYSICS->mcConfig().AddHadronicId(4434); + PHYSICS->mcConfig().AddHadronicId(4444); + PHYSICS->mcConfig().AddHadronicId(5101); + PHYSICS->mcConfig().AddHadronicId(5103); + PHYSICS->mcConfig().AddHadronicId(5112); + PHYSICS->mcConfig().AddHadronicId(5114); + PHYSICS->mcConfig().AddHadronicId(5122); + PHYSICS->mcConfig().AddHadronicId(5132); + PHYSICS->mcConfig().AddHadronicId(5142); + PHYSICS->mcConfig().AddHadronicId(5201); + PHYSICS->mcConfig().AddHadronicId(5203); + PHYSICS->mcConfig().AddHadronicId(5212); + PHYSICS->mcConfig().AddHadronicId(5214); + PHYSICS->mcConfig().AddHadronicId(5222); + PHYSICS->mcConfig().AddHadronicId(5224); + PHYSICS->mcConfig().AddHadronicId(5232); + PHYSICS->mcConfig().AddHadronicId(5242); + PHYSICS->mcConfig().AddHadronicId(5301); + PHYSICS->mcConfig().AddHadronicId(5303); + PHYSICS->mcConfig().AddHadronicId(5312); + PHYSICS->mcConfig().AddHadronicId(5314); + PHYSICS->mcConfig().AddHadronicId(5322); + PHYSICS->mcConfig().AddHadronicId(5324); + PHYSICS->mcConfig().AddHadronicId(5332); + PHYSICS->mcConfig().AddHadronicId(5334); + PHYSICS->mcConfig().AddHadronicId(5342); + PHYSICS->mcConfig().AddHadronicId(5401); + PHYSICS->mcConfig().AddHadronicId(5403); + PHYSICS->mcConfig().AddHadronicId(5412); + PHYSICS->mcConfig().AddHadronicId(5414); + PHYSICS->mcConfig().AddHadronicId(5422); + PHYSICS->mcConfig().AddHadronicId(5424); + PHYSICS->mcConfig().AddHadronicId(5432); + PHYSICS->mcConfig().AddHadronicId(5434); + PHYSICS->mcConfig().AddHadronicId(5442); + PHYSICS->mcConfig().AddHadronicId(5444); + PHYSICS->mcConfig().AddHadronicId(5503); + PHYSICS->mcConfig().AddHadronicId(5512); + PHYSICS->mcConfig().AddHadronicId(5514); + PHYSICS->mcConfig().AddHadronicId(5522); + PHYSICS->mcConfig().AddHadronicId(5524); + PHYSICS->mcConfig().AddHadronicId(5532); + PHYSICS->mcConfig().AddHadronicId(5534); + PHYSICS->mcConfig().AddHadronicId(5542); + PHYSICS->mcConfig().AddHadronicId(5544); + PHYSICS->mcConfig().AddHadronicId(5554); + PHYSICS->mcConfig().AddHadronicId(10111); + PHYSICS->mcConfig().AddHadronicId(10113); + PHYSICS->mcConfig().AddHadronicId(10211); + PHYSICS->mcConfig().AddHadronicId(10213); + PHYSICS->mcConfig().AddHadronicId(10221); + PHYSICS->mcConfig().AddHadronicId(10223); + PHYSICS->mcConfig().AddHadronicId(10311); + PHYSICS->mcConfig().AddHadronicId(10313); + PHYSICS->mcConfig().AddHadronicId(10321); + PHYSICS->mcConfig().AddHadronicId(10323); + PHYSICS->mcConfig().AddHadronicId(10331); + PHYSICS->mcConfig().AddHadronicId(10333); + PHYSICS->mcConfig().AddHadronicId(10411); + PHYSICS->mcConfig().AddHadronicId(10413); + PHYSICS->mcConfig().AddHadronicId(10421); + PHYSICS->mcConfig().AddHadronicId(10423); + PHYSICS->mcConfig().AddHadronicId(10431); + PHYSICS->mcConfig().AddHadronicId(10433); + PHYSICS->mcConfig().AddHadronicId(10441); + PHYSICS->mcConfig().AddHadronicId(10443); + PHYSICS->mcConfig().AddHadronicId(10511); + PHYSICS->mcConfig().AddHadronicId(10513); + PHYSICS->mcConfig().AddHadronicId(10521); + PHYSICS->mcConfig().AddHadronicId(10523); + PHYSICS->mcConfig().AddHadronicId(10531); + PHYSICS->mcConfig().AddHadronicId(10533); + PHYSICS->mcConfig().AddHadronicId(10541); + PHYSICS->mcConfig().AddHadronicId(10543); + PHYSICS->mcConfig().AddHadronicId(10551); + PHYSICS->mcConfig().AddHadronicId(10553); + PHYSICS->mcConfig().AddHadronicId(20113); + PHYSICS->mcConfig().AddHadronicId(20213); + PHYSICS->mcConfig().AddHadronicId(20223); + PHYSICS->mcConfig().AddHadronicId(20313); + PHYSICS->mcConfig().AddHadronicId(20323); + PHYSICS->mcConfig().AddHadronicId(20333); + PHYSICS->mcConfig().AddHadronicId(20413); + PHYSICS->mcConfig().AddHadronicId(20423); + PHYSICS->mcConfig().AddHadronicId(20433); + PHYSICS->mcConfig().AddHadronicId(20443); + PHYSICS->mcConfig().AddHadronicId(20513); + PHYSICS->mcConfig().AddHadronicId(20523); + PHYSICS->mcConfig().AddHadronicId(20533); + PHYSICS->mcConfig().AddHadronicId(20543); + PHYSICS->mcConfig().AddHadronicId(20553); + PHYSICS->mcConfig().AddHadronicId(100443); + PHYSICS->mcConfig().AddHadronicId(100553); + PHYSICS->mcConfig().AddHadronicId(9900440); + PHYSICS->mcConfig().AddHadronicId(9900441); + PHYSICS->mcConfig().AddHadronicId(9900443); + PHYSICS->mcConfig().AddHadronicId(9900551); + PHYSICS->mcConfig().AddHadronicId(9900553); + PHYSICS->mcConfig().AddHadronicId(9910441); + PHYSICS->mcConfig().AddHadronicId(9910551); + } + + void AddDefaultInvisible() + { + // definition of the multiparticle "invisible" + PHYSICS->mcConfig().AddInvisibleId(-16); + PHYSICS->mcConfig().AddInvisibleId(-14); + PHYSICS->mcConfig().AddInvisibleId(-12); + PHYSICS->mcConfig().AddInvisibleId(12); + PHYSICS->mcConfig().AddInvisibleId(14); + PHYSICS->mcConfig().AddInvisibleId(16); + PHYSICS->mcConfig().AddInvisibleId(1000022); + PHYSICS->mcConfig().AddInvisibleId(1000039); + } + + protected : + + }; } diff --git a/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.cpp b/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.cpp index 204e772f..d54fdcee 100644 --- a/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.cpp +++ b/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.cpp @@ -31,8 +31,7 @@ #include "SampleAnalyzer/Process/JetClustering/NullSmearer.h" #ifdef MA5_FASTJET_MODE - // FastJet headers - #include + #include "SampleAnalyzer/Interfaces/substructure/VariableR.h" #endif using namespace MA5; @@ -40,7 +39,7 @@ using namespace MA5; /// Set isolation cones for tracks, e, mu, photon based on tower objects template void SetConeRadius( - std::vector cone_radius, std::vector& objects, MCParticleFormat part, MAbool addself=false + std::vector cone_radius, std::vector& objects, MCParticleFormat part, MAbool addself=false ) { for (MAuint32 iR=0; iR& options) { - // algo defined ? - if (algo_==0) return false; - - // configure tagger - myBtagger_ = new bTagger(); - myCtagger_ = new cTagger(); - myTautagger_ = new TauTagger(); - mySmearer_ = new NullSmearer(); - mySmearer_->Initialize(true); - - // Loop over options - for (std::map::const_iterator - it=options.begin();it!=options.end();it++) - { - std::string key = ClusterAlgoBase::Lower(it->first); - MAbool result=false; - - // exclusive_id - if (key=="exclusive_id") + // algo defined ? + if (algo_==0) return false; + + // configure tagger + myBtagger_ = new bTagger(); + myCtagger_ = new cTagger(); + myTautagger_ = new TauTagger(); + mySmearer_ = new NullSmearer(); + mySmearer_->Initialize(true); + + // Loop over options + for (std::map::const_iterator + it=options.begin();it!=options.end();it++) { - MAint32 tmp=0; - std::stringstream str; - str << it->second; - str >> tmp; - try - { - if (tmp==1) ExclusiveId_=true; - else if (tmp==0) ExclusiveId_=false; - else throw EXCEPTION_WARNING("'exclusive_id' must be equal to 0 or 1. Using default value 'exclusive_id' = "+CONVERT->ToString(ExclusiveId_),"",0); - } - catch(const std::exception& e) - { - MANAGE_EXCEPTION(e); - } - result=true; - } + std::string key = ClusterAlgoBase::Lower(it->first); + MAbool result=false; - // b-tagging - else if (key.find("bjet_id.")==0) - { - result=myBtagger_->SetParameter(key.substr(8),it->second,"bjet_id."); - } + // exclusive_id + if (key=="exclusive_id") + { + MAint32 tmp=0; + std::stringstream str; + str << it->second; + str >> tmp; + try + { + if (tmp==1) ExclusiveId_=true; + else if (tmp==0) ExclusiveId_=false; + else throw EXCEPTION_WARNING("'exclusive_id' must be equal to 0 or 1. Using default value 'exclusive_id' = "+CONVERT->ToString(ExclusiveId_),"",0); + } + catch(const std::exception& e) + { + MANAGE_EXCEPTION(e); + } + result=true; + } - // c-tagging - // else if (key.find("cjet_id.")==0) - // { - // result=myCtagger_->SetParameter(key.substr(8),it->second,"cjet_id."); - // } + // b-tagging + else if (key.find("bjet_id.")==0) + { + result=myBtagger_->SetParameter(key.substr(8),it->second,"bjet_id."); + } - // tau-tagging - else if (key.find("tau_id.")==0) - { - result=myTautagger_->SetParameter(key.substr(7),it->second,"tau_id."); - } + // c-tagging + // else if (key.find("cjet_id.")==0) + // { + // result=myCtagger_->SetParameter(key.substr(8),it->second,"cjet_id."); + // } - // clustering algo - else if (key.find("cluster.")==0) - { - result=algo_->SetParameter(key.substr(8),it->second); - } + // tau-tagging + else if (key.find("tau_id.")==0) + { + result=myTautagger_->SetParameter(key.substr(7),it->second,"tau_id."); + } - // Primary Jet ID - else if (key == "jetid") - { - JetID_ = it->second; - result = true; - } + // clustering algo + else if (key.find("cluster.")==0) + { + result=algo_->SetParameter(key.substr(8),it->second); + } - // Isolation cone radius for tracker - else if (key.find("isolation")==0) - { - std::stringstream str(it->second); - for (MAfloat64 tmp; str >> tmp;) + // Primary Jet ID + else if (key == "jetid") { - if (tmp>0. && key.substr(10) == "track.radius") isocone_track_radius_.push_back(tmp); - if (tmp>0. && key.substr(10) == "electron.radius") isocone_electron_radius_.push_back(tmp); - if (tmp>0. && key.substr(10) == "muon.radius") isocone_muon_radius_.push_back(tmp); - if (tmp>0. && key.substr(10) == "photon.radius") isocone_photon_radius_.push_back(tmp); - if (str.peek() == ',' || str.peek() == ' ') str.ignore(); + JetID_ = it->second; + result = true; } - result = true; - } - // Other - try - { - if (!result) throw EXCEPTION_WARNING("Parameter = "+key+" unknown. It will be skipped.","",0); - } - catch(const std::exception& e) - { - MANAGE_EXCEPTION(e); - } + // Isolation cone radius for tracker + else if (key.find("isolation")==0) + { + std::stringstream str(it->second); + for (MAfloat64 tmp; str >> tmp;) + { + if (tmp>0. && key.substr(10) == "track.radius") isocone_track_radius_.push_back(tmp); + if (tmp>0. && key.substr(10) == "electron.radius") isocone_electron_radius_.push_back(tmp); + if (tmp>0. && key.substr(10) == "muon.radius") isocone_muon_radius_.push_back(tmp); + if (tmp>0. && key.substr(10) == "photon.radius") isocone_photon_radius_.push_back(tmp); + if (str.peek() == ',' || str.peek() == ' ') str.ignore(); + } + result = true; + } - } + // Other + try + { + if (!result) throw EXCEPTION_WARNING("Parameter = "+key+" unknown. It will be skipped.","",0); + } + catch(const std::exception& e) + { + MANAGE_EXCEPTION(e); + } - // configure algo - algo_->Initialize(); + } + // configure algo + algo_->Initialize(); - return true; + + return true; } @@ -177,11 +176,11 @@ MAbool JetClusterer::Initialize(const std::map& options // ----------------------------------------------------------------------------- void JetClusterer::Finalize() { - if (algo_!=0) delete algo_; - if (myBtagger_!=0) delete myBtagger_; - if (myCtagger_!=0) delete myCtagger_; - if (myTautagger_!=0) delete myTautagger_; - if (mySmearer_!=0) delete mySmearer_; + if (algo_!=0) delete algo_; + if (myBtagger_!=0) delete myBtagger_; + if (myCtagger_!=0) delete myCtagger_; + if (myTautagger_!=0) delete myTautagger_; + if (mySmearer_!=0) delete mySmearer_; } @@ -190,11 +189,11 @@ void JetClusterer::Finalize() // ----------------------------------------------------------------------------- void JetClusterer::GetFinalState(const MCParticleFormat* part, std::set& finalstates) { - for (MAuint32 i=0; idaughters().size(); i++) - { - if (PHYSICS->Id->IsFinalState(part->daughters()[i])) finalstates.insert(part->daughters()[i]); - else return GetFinalState(part->daughters()[i],finalstates); - } + for (MAuint32 i=0; idaughters().size(); i++) + { + if (PHYSICS->Id->IsFinalState(part->daughters()[i])) finalstates.insert(part->daughters()[i]); + else return GetFinalState(part->daughters()[i],finalstates); + } } @@ -203,11 +202,11 @@ void JetClusterer::GetFinalState(const MCParticleFormat* part, std::setdaughters().size(); i++) - { - if (part->daughters()[i]->pdgid()==part->pdgid()) return false; - } - return true; + for (MAuint32 i=0; idaughters().size(); i++) + { + if (part->daughters()[i]->pdgid()==part->pdgid()) return false; + } + return true; } // ----------------------------------------------------------------------------- @@ -215,335 +214,327 @@ MAbool JetClusterer::IsLast(const MCParticleFormat* part, EventFormat& myEvent) // ----------------------------------------------------------------------------- MAbool JetClusterer::Execute(SampleFormat& mySample, EventFormat& myEvent) { - // Safety - if (mySample.mc()==0 || myEvent.mc()==0) return false; - if (mySample.rec()==0) mySample.InitializeRec(); - if (myEvent.rec() ==0) myEvent.InitializeRec(); + // Safety + if (mySample.mc()==0 || myEvent.mc()==0) return false; + if (mySample.rec()==0) mySample.InitializeRec(); + if (myEvent.rec() ==0) myEvent.InitializeRec(); - // Reseting the reconstructed event - myEvent.rec()->Reset(); + // Reseting the reconstructed event + myEvent.rec()->Reset(); - // Veto - std::vector vetos(myEvent.mc()->particles().size(),false); - std::set vetos2; + // Veto + std::vector vetos(myEvent.mc()->particles().size(),false); + std::set vetos2; - // Filling the dataformat with electron/muon - for (MAuint32 i=0;iparticles().size();i++) - { - const MCParticleFormat& part = myEvent.mc()->particles()[i]; - MAuint32 absid = std::abs(part.pdgid()); + // Filling the dataformat with electron/muon + for (MAuint32 i=0;iparticles().size();i++) + { + const MCParticleFormat& part = myEvent.mc()->particles()[i]; + MAuint32 absid = std::abs(part.pdgid()); - // Rejecting particle with a null pt (initial state ?) - if (part.pt()<1e-10) continue; + // Rejecting particle with a null pt (initial state ?) + if (part.pt()<1e-10) continue; - // Run particle propagator - if (mySmearer_->isPropagatorOn() && part.mothers().size()>0) - mySmearer_->ParticlePropagator(const_cast(&part)); + // Run particle propagator + if (mySmearer_->isPropagatorOn() && part.mothers().size()>0) + mySmearer_->ParticlePropagator(const_cast(&part)); -#ifdef MA5_FASTJET_MODE - // @JACK delphes based analyses already has tracks - // Set up tracks as charged FS particles OR charged interstate particles with nonzero ctau - if (PDG->IsCharged(part.pdgid()) && part.mothers().size()>0) - { - // Minimum tracking requirement is around 0.5 mm see ref. 1007.1988 - if (part.ctau() > 0. || PHYSICS->Id->IsFinalState(part)) + // @JACK delphes based analyses already has tracks + // Set up tracks as charged FS particles OR charged interstate particles with nonzero ctau + if (PDG->IsCharged(part.pdgid()) && part.mothers().size()>0 && algo_!=0) { - // Reminder: -1 is reserved for the tracks - MCParticleFormat smeared_track = mySmearer_->Execute(&part, -1); - if (smeared_track.pt() > 1e-5) + // Minimum tracking requirement is around 0.5 mm see ref. 1007.1988 + if (part.ctau() > 0. || PHYSICS->Id->IsFinalState(part)) { - RecTrackFormat * track = myEvent.rec()->GetNewTrack(); - MALorentzVector trk_mom; - trk_mom.SetPtEtaPhiM(smeared_track.pt(), - smeared_track.eta(), - smeared_track.phi(),0.0); - track->setMomentum(trk_mom); - track->setD0(smeared_track.d0()); - track->setDZ(smeared_track.dz()); - track->setD0Approx(smeared_track.d0_approx()); - track->setDZApprox(smeared_track.dz_approx()); - MAdouble64 ctau = PHYSICS->Id->IsFinalState(part) ? 0.0 : part.mothers()[0]->ctau(); - MALorentzVector new_vertex(part.mothers()[0]->decay_vertex().X(), - part.mothers()[0]->decay_vertex().Y(), - part.mothers()[0]->decay_vertex().Z(), ctau); - track->setProductionVertex(new_vertex); - track->setClosestApproach(smeared_track.closest_approach()); - track->setMc(&(part)); - track->SetCharge(PDG->GetCharge(part.pdgid()) / 3.); + // Reminder: -1 is reserved for the tracks + MCParticleFormat smeared_track = mySmearer_->Execute(&part, -1); + if (smeared_track.pt() > 1e-5) + { + RecTrackFormat * track = myEvent.rec()->GetNewTrack(); + MALorentzVector trk_mom; + trk_mom.SetPtEtaPhiM(smeared_track.pt(), + smeared_track.eta(), + smeared_track.phi(),0.0); + track->setMomentum(trk_mom); + track->setD0(smeared_track.d0()); + track->setDZ(smeared_track.dz()); + track->setD0Approx(smeared_track.d0_approx()); + track->setDZApprox(smeared_track.dz_approx()); + MAdouble64 ctau = PHYSICS->Id->IsFinalState(part) ? 0.0 : part.mothers()[0]->ctau(); + MALorentzVector new_vertex(part.mothers()[0]->decay_vertex().X(), + part.mothers()[0]->decay_vertex().Y(), + part.mothers()[0]->decay_vertex().Z(), ctau); + track->setProductionVertex(new_vertex); + track->setClosestApproach(smeared_track.closest_approach()); + track->setMc(&(part)); + track->SetCharge(PDG->GetCharge(part.pdgid()) / 3.); + } } } - } -#endif - // Treating intermediate particles - if (PHYSICS->Id->IsInterState(part)) - { - // rejecting not interesting particles - if (absid!=5 && absid!=4 && absid!=15) continue; - - // keeping the last particle with the same id in the decay chain - if (!IsLast(&part, myEvent)) continue; - - // looking for b quarks - if (absid==5) - { - MAbool found=false; - for (MAuint32 j=0;jMCBquarks_.size();j++) - { - if (myEvent.rec()->MCBquarks_[j]==&(part)) - {found=true; break;} - } - if (!found) myEvent.rec()->MCBquarks_.push_back(&(part)); - } - - // looking for c quarks - else if (absid==4) - { - MAbool found=false; - for (MAuint32 j=0;jMCCquarks_.size();j++) + // Treating intermediate particles + if (PHYSICS->Id->IsInterState(part)) { - if (myEvent.rec()->MCCquarks_[j]==&(part)) - {found=true; break;} - } - if (!found) myEvent.rec()->MCCquarks_.push_back(&(part)); - } - - // looking for taus - else if (absid==15) - { - // rejecting particle if coming from hadronization - if (LOOP->ComingFromHadronDecay(&part,mySample,myEvent.mc()->particles().size())) continue; - - // Looking taus daughters id - MAbool leptonic = true; - MAbool muonic = false; - MAbool electronic = false; - for (MAuint32 j=0;jpdgid()); - if (pdgid==13) muonic=true; - else if (pdgid==11) electronic=true; - else if (pdgid!=22 /*photons*/ && - !(pdgid>=11 && pdgid<=16) /*neutrinos*/) - leptonic=false; - } - if (!leptonic) {muonic=false; electronic=false;} + // rejecting not interesting particles + if (absid!=5 && absid!=4 && absid!=15) continue; - // Saving taus decaying into muons (only one copy) - if (muonic) - { - MAbool found=false; - for (MAuint32 j=0;jMCMuonicTaus_.size();j++) - { - if (myEvent.rec()->MCMuonicTaus_[j]==&(part)) - {found=true; break;} - } - if (!found) myEvent.rec()->MCMuonicTaus_.push_back(&(part)); - } + // keeping the last particle with the same id in the decay chain + if (!IsLast(&part, myEvent)) continue; - // Saving taus decaying into electrons (only one copy) - else if (electronic) - { - MAbool found=false; - for (MAuint32 j=0;jMCElectronicTaus_.size();j++) - { - if (myEvent.rec()->MCElectronicTaus_[j]==&(part)) - {found=true; break;} - } - if (!found) myEvent.rec()->MCElectronicTaus_.push_back(&(part)); - } + // looking for b quarks + if (absid==5) + { + MAbool found=false; + for (MAuint32 j=0;jMCBquarks_.size();j++) + { + if (myEvent.rec()->MCBquarks_[j]==&(part)) + {found=true; break;} + } + if (!found) myEvent.rec()->MCBquarks_.push_back(&(part)); + } - // Saving taus decaying into hadrons (only copy) - else - { - MAbool found=false; - for (MAuint32 j=0;jMCHadronicTaus_.size();j++) - { - if (myEvent.rec()->MCHadronicTaus_[j]==&(part)) - {found=true; break;} - } - if (!found) - { - // Saving the hadrons in MC container - myEvent.rec()->MCHadronicTaus_.push_back(&(part)); - - // Applying efficiency - if (!myTautagger_->IsIdentified()) continue; - - // Smearing the hadronic taus - MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); - // If smeared pt is zero, no need to count the particle but it still needs - // to be vetoed for jet clustering. - if (smeared.pt() > 1e-10) + // looking for c quarks + else if (absid==4) { - // Creating reco hadronic taus - RecTauFormat* myTau = myEvent.rec()->GetNewTau(); - if (part.pdgid()>0) myTau->setCharge(-1); - else myTau->setCharge(+1); - myTau->setMomentum(smeared.momentum()); - myTau->setD0(smeared.d0()); - myTau->setDZ(smeared.dz()); - myTau->setD0Approx(smeared.d0_approx()); - myTau->setDZApprox(smeared.dz_approx()); - myTau->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), - part.mothers()[0]->decay_vertex().Y(), - part.mothers()[0]->decay_vertex().Z(),0.0)); - myTau->setClosestApproach(smeared.closest_approach()); - myTau->setMc(&part); - myTau->setDecayMode(PHYSICS->GetTauDecayMode(myTau->mc())); - if (myTau->DecayMode()<=0) myTau->setNtracks(0); // ERROR case - else if (myTau->DecayMode()==7 || - myTau->DecayMode()==9) myTau->setNtracks(3); // 3-Prong - else myTau->setNtracks(1); // 1-Prong + MAbool found=false; + for (MAuint32 j=0;jMCCquarks_.size();j++) + { + if (myEvent.rec()->MCCquarks_[j]==&(part)) + {found=true; break;} + } + if (!found) myEvent.rec()->MCCquarks_.push_back(&(part)); } - // Searching final state - GetFinalState(&part,vetos2); - } - } - } - } + // looking for taus + else if (absid==15) + { + // rejecting particle if coming from hadronization + if (LOOP->ComingFromHadronDecay(&part,mySample,myEvent.mc()->particles().size())) continue; + + // Looking taus daughters id + MAbool leptonic = true; + MAbool muonic = false; + MAbool electronic = false; + for (MAuint32 j=0;jpdgid()); + if (pdgid==13) muonic=true; + else if (pdgid==11) electronic=true; + else if (pdgid!=22 /*photons*/ && + !(pdgid>=11 && pdgid<=16) /*neutrinos*/) + leptonic=false; + } + if (!leptonic) {muonic=false; electronic=false;} - // Keeping only final states - else if (PHYSICS->Id->IsFinalState(part)) - { - // keeping only electron, muon and photon - if (absid==22 || absid==11 || absid==13) - { - // rejecting particle if coming from hadronization - if ( !(ExclusiveId_ && LOOP->ComingFromHadronDecay(&part,mySample))) - { - - // Muons - if (absid==13) - { - vetos[i]=true; - - // Smearing its momentum - MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); - if (smeared.pt() <= 1e-10) continue; - - RecLeptonFormat * muon = myEvent.rec()->GetNewMuon(); - muon->setMomentum(smeared.momentum()); - muon->setD0(smeared.d0()); - muon->setDZ(smeared.dz()); - muon->setD0Approx(smeared.d0_approx()); - muon->setDZApprox(smeared.dz_approx()); - muon->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), - part.mothers()[0]->decay_vertex().Y(), - part.mothers()[0]->decay_vertex().Z(),0.0)); - muon->setClosestApproach(smeared.closest_approach()); - muon->setMc(&(part)); - if (part.pdgid()==13) muon->SetCharge(-1); - else muon->SetCharge(+1); - } - - // Electrons - else if (absid==11) - { - vetos[i]=true; - - // Smearing the electron momentum - MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); - if (smeared.pt() <= 1e-10) continue; - - RecLeptonFormat * elec = myEvent.rec()->GetNewElectron(); - elec->setMomentum(smeared.momentum()); - elec->setD0(smeared.d0()); - elec->setDZ(smeared.dz()); - elec->setD0Approx(smeared.d0_approx()); - elec->setDZApprox(smeared.dz_approx()); - elec->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), - part.mothers()[0]->decay_vertex().Y(), - part.mothers()[0]->decay_vertex().Z(),0.0)); - elec->setClosestApproach(smeared.closest_approach()); - elec->setMc(&(part)); - if (part.pdgid()==11) elec->SetCharge(-1); - else elec->SetCharge(+1); - } - - // Photons - else if (absid==22) - { - if (!LOOP->IrrelevantPhoton(&part,mySample)) + // Saving taus decaying into muons (only one copy) + if (muonic) + { + MAbool found=false; + for (MAuint32 j=0;jMCMuonicTaus_.size();j++) + { + if (myEvent.rec()->MCMuonicTaus_[j]==&(part)) + {found=true; break;} + } + if (!found) myEvent.rec()->MCMuonicTaus_.push_back(&(part)); + } + + // Saving taus decaying into electrons (only one copy) + else if (electronic) { - vetos[i]=true; - - // Smearing the photon momentum - MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); - if (smeared.pt() <= 1e-10) continue; - - RecPhotonFormat * photon = myEvent.rec()->GetNewPhoton(); - photon->setMomentum(smeared.momentum()); - photon->setD0(smeared.d0()); - photon->setDZ(smeared.dz()); - photon->setD0Approx(smeared.d0_approx()); - photon->setDZApprox(smeared.dz_approx()); - photon->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), - part.mothers()[0]->decay_vertex().Y(), - part.mothers()[0]->decay_vertex().Z(),0.0)); - photon->setClosestApproach(smeared.closest_approach()); - photon->setMc(&(part)); + MAbool found=false; + for (MAuint32 j=0;jMCElectronicTaus_.size();j++) + { + if (myEvent.rec()->MCElectronicTaus_[j]==&(part)) + {found=true; break;} + } + if (!found) myEvent.rec()->MCElectronicTaus_.push_back(&(part)); } - } + + // Saving taus decaying into hadrons (only copy) + else + { + MAbool found=false; + for (MAuint32 j=0;jMCHadronicTaus_.size();j++) + { + if (myEvent.rec()->MCHadronicTaus_[j]==&(part)) + {found=true; break;} + } + if (!found) + { + // Saving the hadrons in MC container + myEvent.rec()->MCHadronicTaus_.push_back(&(part)); + + // Applying efficiency + if (!myTautagger_->IsIdentified()) continue; + + // Smearing the hadronic taus + MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); + // If smeared pt is zero, no need to count the particle but it still needs + // to be vetoed for jet clustering. + if (smeared.pt() > 1e-10) + { + // Creating reco hadronic taus + RecTauFormat* myTau = myEvent.rec()->GetNewTau(); + if (part.pdgid()>0) myTau->setCharge(-1); + else myTau->setCharge(+1); + myTau->setMomentum(smeared.momentum()); + myTau->setD0(smeared.d0()); + myTau->setDZ(smeared.dz()); + myTau->setD0Approx(smeared.d0_approx()); + myTau->setDZApprox(smeared.dz_approx()); + myTau->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), + part.mothers()[0]->decay_vertex().Y(), + part.mothers()[0]->decay_vertex().Z(),0.0)); + myTau->setClosestApproach(smeared.closest_approach()); + myTau->setMc(&part); + myTau->setDecayMode(PHYSICS->GetTauDecayMode(myTau->mc())); + if (myTau->DecayMode()<=0) myTau->setNtracks(0); // ERROR case + else if (myTau->DecayMode()==7 || + myTau->DecayMode()==9) myTau->setNtracks(3); // 3-Prong + else myTau->setNtracks(1); // 1-Prong + } + + // Searching final state + GetFinalState(&part,vetos2); + } + } + } } - } -#ifdef MA5_FASTJET_MODE - // Putting the good inputs into the containter - // Good inputs = - final state - // - visible - // - if exclusiveID=1: particles not vetoed - // - if exclusiveID=0: all particles except muons - if (PHYSICS->Id->IsInvisible(part)) continue; - - // ExclusiveId mode - if (ExclusiveId_) - { - if (vetos[i]) continue; - if (vetos2.find(&part)!=vetos2.end()) continue; - } - // NonExclusive Id mode - else if (std::abs(part.pdgid())==13) continue; + // Keeping only final states + else if (PHYSICS->Id->IsFinalState(part)) + { + // keeping only electron, muon and photon + if (absid==22 || absid==11 || absid==13) + { + // rejecting particle if coming from hadronization + if ( !(ExclusiveId_ && LOOP->ComingFromHadronDecay(&part,mySample))) + { - // Smearer module returns a smeared MCParticleFormat object - // Default: NullSmearer, that does nothing - // Reminder: 0 is reserved for the jet constituents - MCParticleFormat smeared = mySmearer_->Execute(&part, 0); - if (smeared.pt() <= 1e-10) continue; + // Muons + if (absid==13) + { + vetos[i]=true; + + // Smearing its momentum + MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); + if (smeared.pt() <= 1e-10) continue; + + RecLeptonFormat * muon = myEvent.rec()->GetNewMuon(); + muon->setMomentum(smeared.momentum()); + muon->setD0(smeared.d0()); + muon->setDZ(smeared.dz()); + muon->setD0Approx(smeared.d0_approx()); + muon->setDZApprox(smeared.dz_approx()); + muon->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), + part.mothers()[0]->decay_vertex().Y(), + part.mothers()[0]->decay_vertex().Z(),0.0)); + muon->setClosestApproach(smeared.closest_approach()); + muon->setMc(&(part)); + if (part.pdgid()==13) muon->SetCharge(-1); + else muon->SetCharge(+1); + } + + // Electrons + else if (absid==11) + { + vetos[i]=true; + + // Smearing the electron momentum + MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); + if (smeared.pt() <= 1e-10) continue; + + RecLeptonFormat * elec = myEvent.rec()->GetNewElectron(); + elec->setMomentum(smeared.momentum()); + elec->setD0(smeared.d0()); + elec->setDZ(smeared.dz()); + elec->setD0Approx(smeared.d0_approx()); + elec->setDZApprox(smeared.dz_approx()); + elec->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), + part.mothers()[0]->decay_vertex().Y(), + part.mothers()[0]->decay_vertex().Z(),0.0)); + elec->setClosestApproach(smeared.closest_approach()); + elec->setMc(&(part)); + if (part.pdgid()==11) elec->SetCharge(-1); + else elec->SetCharge(+1); + } + + // Photons + else if (absid==22) + { + if (!LOOP->IrrelevantPhoton(&part,mySample)) + { + vetos[i]=true; + + // Smearing the photon momentum + MCParticleFormat smeared = mySmearer_->Execute(&part, static_cast(absid)); + if (smeared.pt() <= 1e-10) continue; + + RecPhotonFormat * photon = myEvent.rec()->GetNewPhoton(); + photon->setMomentum(smeared.momentum()); + photon->setD0(smeared.d0()); + photon->setDZ(smeared.dz()); + photon->setD0Approx(smeared.d0_approx()); + photon->setDZApprox(smeared.dz_approx()); + photon->setProductionVertex(MALorentzVector(part.mothers()[0]->decay_vertex().X(), + part.mothers()[0]->decay_vertex().Y(), + part.mothers()[0]->decay_vertex().Z(),0.0)); + photon->setClosestApproach(smeared.closest_approach()); + photon->setMc(&(part)); + } + } + } + } + // Putting the good inputs into the containter + // Good inputs = - final state + // - visible + // - if exclusiveID=1: particles not vetoed + // - if exclusiveID=0: all particles except muons + if (PHYSICS->Id->IsInvisible(part) || algo_==0) continue; + + // ExclusiveId mode + if (ExclusiveId_) + { + if (vetos[i]) continue; + if (vetos2.find(&part)!=vetos2.end()) continue; + } - // Filling good particle for clustering - fastjet::PseudoJet input; - input.reset(smeared.px(), smeared.py(), smeared.pz(), smeared.e()); - input.set_user_index(i); - myEvent.rec()->AddHadron(input); -#endif + // NonExclusive Id mode + else if (std::abs(part.pdgid())==13) continue; + + // Smearer module returns a smeared MCParticleFormat object + // Default: NullSmearer, that does nothing + // Reminder: 0 is reserved for the jet constituents + MCParticleFormat smeared = mySmearer_->Execute(&part, 0); + if (smeared.pt() <= 1e-10) continue; + // Filling good particle for clustering + myEvent.rec()->AddHadron(smeared, i); + } } - } - // Sorting the objecfts after smearing - if (mySmearer_->isElectronSmearerOn()) - std::sort(myEvent.rec()->electrons_.begin(), myEvent.rec()->electrons_.end(), - [](RecLeptonFormat const & lep1, RecLeptonFormat const & lep2){ return lep1.pt() > lep2.pt(); }); - if (mySmearer_->isMuonSmearerOn()) - std::sort(myEvent.rec()->muons_.begin(), myEvent.rec()->muons_.end(), - [](RecLeptonFormat const & lep1, RecLeptonFormat const & lep2){ return lep1.pt() > lep2.pt(); }); - if (mySmearer_->isTauSmearerOn()) - std::sort(myEvent.rec()->taus_.begin(), myEvent.rec()->taus_.end(), - [](RecTauFormat const & ta1, RecTauFormat const & ta2){ return ta1.pt() > ta2.pt(); }); - if (mySmearer_->isPhotonSmearerOn()) - std::sort(myEvent.rec()->photons_.begin(), myEvent.rec()->photons_.end(), - [](RecPhotonFormat const & ph1, RecPhotonFormat const & ph2){ return ph1.pt() > ph2.pt(); }); - - // Set Primary Jet ID - myEvent.rec()->SetPrimaryJetID(JetID_); - // Launching the clustering - // -> Filling the collection: myEvent->rec()->jets() - algo_->Execute(mySample,myEvent,mySmearer_); + // Sorting the objecfts after smearing + if (mySmearer_->isElectronSmearerOn()) + std::sort(myEvent.rec()->electrons_.begin(), myEvent.rec()->electrons_.end(), + [](RecLeptonFormat const & lep1, RecLeptonFormat const & lep2){ return lep1.pt() > lep2.pt(); }); + if (mySmearer_->isMuonSmearerOn()) + std::sort(myEvent.rec()->muons_.begin(), myEvent.rec()->muons_.end(), + [](RecLeptonFormat const & lep1, RecLeptonFormat const & lep2){ return lep1.pt() > lep2.pt(); }); + if (mySmearer_->isTauSmearerOn()) + std::sort(myEvent.rec()->taus_.begin(), myEvent.rec()->taus_.end(), + [](RecTauFormat const & ta1, RecTauFormat const & ta2){ return ta1.pt() > ta2.pt(); }); + if (mySmearer_->isPhotonSmearerOn()) + std::sort(myEvent.rec()->photons_.begin(), myEvent.rec()->photons_.end(), + [](RecPhotonFormat const & ph1, RecPhotonFormat const & ph2){ return ph1.pt() > ph2.pt(); }); + + // Set Primary Jet ID + myEvent.rec()->SetPrimaryJetID(JetID_); + // Launching the clustering + // -> Filling the collection: myEvent->rec()->jets() + algo_->Execute(mySample,myEvent,mySmearer_); #ifdef MA5_FASTJET_MODE - // Cluster additional jets separately. In order to save time Execute function + // Cluster additional jets separately. In order to save time Execute function // saves hadron inputs into memory and that configuration is used for the rest // of the jets. for (auto &collection_item: cluster_collection_) @@ -553,48 +544,48 @@ MAbool JetClusterer::Execute(SampleFormat& mySample, EventFormat& myEvent) #endif - // shortcut for TET & THT - MAfloat64 & TET = myEvent.rec()->TET(); - // MAfloat64 & THT = myEvent.rec()->THT(); - RecParticleFormat* MET = &(myEvent.rec()->MET()); - RecParticleFormat* MHT = &(myEvent.rec()->MHT()); + // shortcut for TET & THT + MAfloat64 & TET = myEvent.rec()->TET(); + // MAfloat64 & THT = myEvent.rec()->THT(); + RecParticleFormat* MET = &(myEvent.rec()->MET()); + RecParticleFormat* MHT = &(myEvent.rec()->MHT()); - // End - if (ExclusiveId_) - { - for (MAuint32 i=0;ielectrons().size();i++) - { - (*MET) -= myEvent.rec()->electrons()[i].momentum(); - TET += myEvent.rec()->electrons()[i].pt(); - } - for (MAuint32 i=0;iphotons().size();i++) + // End + if (ExclusiveId_) { - (*MET) -= myEvent.rec()->photons()[i].momentum(); - TET += myEvent.rec()->photons()[i].pt(); + for (MAuint32 i=0;ielectrons().size();i++) + { + (*MET) -= myEvent.rec()->electrons()[i].momentum(); + TET += myEvent.rec()->electrons()[i].pt(); + } + for (MAuint32 i=0;iphotons().size();i++) + { + (*MET) -= myEvent.rec()->photons()[i].momentum(); + TET += myEvent.rec()->photons()[i].pt(); + } + for (MAuint32 i=0;itaus().size();i++) + { + (*MET) -= myEvent.rec()->taus()[i].momentum(); + TET += myEvent.rec()->taus()[i].pt(); + } } - for (MAuint32 i=0;itaus().size();i++) + + for (MAuint32 i=0;imuons().size();i++) { - (*MET) -= myEvent.rec()->taus()[i].momentum(); - TET += myEvent.rec()->taus()[i].pt(); + (*MET) -= myEvent.rec()->muons()[i].momentum(); + TET += myEvent.rec()->muons()[i].pt(); } - } - - for (MAuint32 i=0;imuons().size();i++) - { - (*MET) -= myEvent.rec()->muons()[i].momentum(); - TET += myEvent.rec()->muons()[i].pt(); - } - MET->momentum().SetPz(0.); - MET->momentum().SetE(MET->momentum().Pt()); - MHT->momentum().SetPz(0.); - MHT->momentum().SetE(MHT->momentum().Pt()); + MET->momentum().SetPz(0.); + MET->momentum().SetE(MET->momentum().Pt()); + MHT->momentum().SetPz(0.); + MHT->momentum().SetE(MHT->momentum().Pt()); - myBtagger_->Execute(mySample,myEvent); - myTautagger_->Execute(mySample,myEvent); + myBtagger_->Execute(mySample,myEvent); + myTautagger_->Execute(mySample,myEvent); #ifdef MA5_FASTJET_MODE - // Setup isolation cones + // Setup isolation cones if (isocone_track_radius_.size() > 0 || isocone_electron_radius_.size() > 0 || \ isocone_muon_radius_.size() > 0 || isocone_photon_radius_.size() > 0) { @@ -617,5 +608,153 @@ MAbool JetClusterer::Execute(SampleFormat& mySample, EventFormat& myEvent) } #endif - return true; + return true; } + +// Load additional Jets +MAbool JetClusterer::LoadJetConfiguration(std::map options) +{ +#ifdef MA5_FASTJET_MODE + std::string new_jetid; + std::string algorithm; + if (options.find("algorithm") == options.end()) + { + ERROR << "Jet configuration needs to have `algorithm` option. Jet configuration will be ignored." << endmsg; + return true; + } + else algorithm = options["algorithm"]; + if (options.find("JetID") == options.end()) + { + ERROR << "Jet configuration needs to have `JetID` option. Jet configuration will be ignored." << endmsg; + return true; + } + if (substructure_collection_.find(options["JetID"]) != substructure_collection_.end() || \ + cluster_collection_.find(options["JetID"]) != cluster_collection_.end() ) + { + ERROR << "Jet ID " + options["JetID"] + \ + " already exists. Jet configuration will be ignored." << endmsg; + return true; + } + + if (algorithm != "VariableR") + { + std::map clustering_params; + + // decide if its good to keep this jet + ClusterAlgoBase* new_algo; + // Loop over options + for (const auto &it: options) + { + std::string key = ClusterAlgoBase::Lower(it.first); + if (key=="jetid") + { + // Check if JetID is used before + new_jetid = it.second; + continue; + } + + // Find the clustering algorithm + if (key=="algorithm") + { + if (it.second == "antikt") new_algo = new ClusterAlgoStandard("antikt"); + else if (it.second == "cambridge") new_algo = new ClusterAlgoStandard("cambridge"); + else if (it.second == "genkt") new_algo = new ClusterAlgoStandard("genkt"); + else if (it.second == "kt") new_algo = new ClusterAlgoStandard("kt"); + else if (it.second == "siscone") new_algo = new ClusterAlgoSISCone(); + else if (it.second == "cdfmidpoint") new_algo = new ClusterAlgoCDFMidpoint(); + else if (it.second == "cdfjetclu") new_algo = new ClusterAlgoCDFJetClu(); + else if (it.second == "gridjet") new_algo = new ClusterAlgoGridJet(); + else { + ERROR << "Unknown algorithm : " << it.second << ". It will be ignored." << endmsg; + return true; + } + continue; + } + // clustering algo -> keep the previous syntax + else if (key.find("cluster.")==0) + { + clustering_params.insert(std::pair(key.substr(8),it.second)); + continue; + } + + // Other + try + { + throw EXCEPTION_WARNING("Parameter = "+key+" unknown. It will be skipped.","",0); + } + catch(const std::exception& e) + { + MANAGE_EXCEPTION(e); + return false; + } + } + + cluster_collection_.insert(std::pair(new_jetid,new_algo)); + for (const auto &it: clustering_params) + cluster_collection_[new_jetid]->SetParameter(it.first, it.second); + std::string algoname = cluster_collection_[new_jetid]->GetName(); + std::string params = cluster_collection_[new_jetid]->GetParameters(); + INFO << " - Adding Jet ID : " << new_jetid << endmsg; + INFO << " with algo : " << algoname << ", " << params << endmsg; + cluster_collection_[new_jetid]->Initialize(); + } + else if (algorithm == "VariableR") + { + for (std::string key: {"rho", "minR", "maxR", "PTmin", "clustertype", "strategy", "exclusive"}) + { + if (options.find("cluster."+key) == options.end()) + { + ERROR << "Option 'cluster." + key + "' is missing. VariableR clustering will be ignored." << endmsg; + return true; + } + } + MAfloat32 rho = std::stof(options["cluster.rho"]); + MAfloat32 minR = std::stof(options["cluster.minR"]); + MAfloat32 maxR = std::stof(options["cluster.maxR"]); + MAfloat32 ptmin = std::stof(options["cluster.PTmin"]); + MAbool isExclusive = (options["cluster.exclusive"] == "1"); + + Substructure::VariableR::ClusterType ctype = Substructure::VariableR::AKTLIKE; + if (options["cluster.clustertype"] == "CALIKE") ctype = Substructure::VariableR::CALIKE; + else if (options["cluster.clustertype"] == "KTLIKE") ctype = Substructure::VariableR::KTLIKE; + else if (options["cluster.clustertype"] == "AKTLIKE") ctype = Substructure::VariableR::AKTLIKE; + + Substructure::VariableR::Strategy strategy = Substructure::VariableR::Best; + if (options["cluster.strategy"] == "Best") strategy = Substructure::VariableR::Best; + else if (options["cluster.strategy"] == "N2Tiled") strategy = Substructure::VariableR::N2Tiled; + else if (options["cluster.strategy"] == "N2Plain") strategy = Substructure::VariableR::N2Plain; + else if (options["cluster.strategy"] == "NNH") strategy = Substructure::VariableR::NNH; + else if (options["cluster.strategy"] == "Native") strategy = Substructure::VariableR::Native; + + Substructure::VariableR* variableR; + variableR = new Substructure::VariableR(rho, minR, maxR, ctype, strategy, ptmin, isExclusive); + + substructure_collection_.insert( + std::pair(options["JetID"], variableR) + ); + + std::string exclusive = isExclusive ? "True" : "False"; + INFO << " - Adding Jet ID : " << options["JetID"] << endmsg; + INFO << " with algo : VariableR" << ", " + << "rho = " << options["cluster.rho"] << ", " + << "minR = " << options["cluster.minR"] << ", " + << "maxR = " << options["cluster.maxR"] << ", " + << "ptmin = " << options["cluster.PTmin"] << ", \n" + << " " + << "isExclusive = " << exclusive << ", " + << "clustertype = " << options["cluster.clustertype"] << ", " + << "strategy = " << options["cluster.strategy"] + << endmsg; + } + else + { + ERROR << "Unknown algorithm: " << algorithm << endmsg; + return false; + } + + return true; +#else + ERROR << "FastJet has not been enabled. Can not add jets to the analysis." << endmsg; + return true; +#endif +} \ No newline at end of file diff --git a/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.h b/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.h index 5c60da8d..87911de3 100644 --- a/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.h +++ b/tools/SampleAnalyzer/Process/JetClustering/JetClusterer.h @@ -42,8 +42,7 @@ #include "SampleAnalyzer/Interfaces/fastjet/ClusterAlgoCDFMidpoint.h" #include "SampleAnalyzer/Interfaces/fastjet/ClusterAlgoCDFJetClu.h" #include "SampleAnalyzer/Interfaces/fastjet/ClusterAlgoGridJet.h" - #include "SampleAnalyzer/Interfaces/fastjet/VariableR.h" - #include "SampleAnalyzer/Interfaces/fastjet/ClusterBase.h" + #include "SampleAnalyzer/Interfaces/substructure/ClusterBase.h" #endif // STL headers @@ -113,6 +112,7 @@ namespace MA5 algo_ = algo; #ifdef MA5_FASTJET_MODE cluster_collection_.clear(); + substructure_collection_.clear(); #endif myBtagger_ = 0; myCtagger_ = 0; @@ -153,152 +153,7 @@ namespace MA5 } // Load additional Jets - MAbool LoadJetConfiguration(std::map options) - { -#ifdef MA5_FASTJET_MODE - - std::string new_jetid; - std::string algorithm; - if (options.find("algorithm") == options.end()) - { - ERROR << "Jet configuration needs to have an `algorithm` option. Jet configuration ignored." << endmsg; - return true; - } - else algorithm = options["algorithm"]; - if (options.find("JetID") == options.end()) - { - ERROR << "Jet configuration needs to have a `JetID` option. Jet configuration ignored." << endmsg; - return true; - } - if (substructure_collection_.find(options["JetID"]) != substructure_collection_.end() || \ - cluster_collection_.find(options["JetID"]) != cluster_collection_.end() ) - { - ERROR << "Jet ID " + options["JetID"] + " already defined. Jet configuration ignored." << endmsg; - return true; - } - - if (algorithm != "VariableR") - { - std::map clustering_params; - - // decide if its good to keep this jet - ClusterAlgoBase* new_algo; - // Loop over options - for (const auto &it: options) - { - std::string key = ClusterAlgoBase::Lower(it.first); - if (key=="jetid") - { - // Check if JetID is used before - new_jetid = it.second; - continue; - } - - // Find the clustering algorithm - if (key=="algorithm") - { - if (it.second == "antikt") new_algo = new ClusterAlgoStandard("antikt"); - else if (it.second == "cambridge") new_algo = new ClusterAlgoStandard("cambridge"); - else if (it.second == "genkt") new_algo = new ClusterAlgoStandard("genkt"); - else if (it.second == "kt") new_algo = new ClusterAlgoStandard("kt"); - else if (it.second == "siscone") new_algo = new ClusterAlgoSISCone(); - else if (it.second == "cdfmidpoint") new_algo = new ClusterAlgoCDFMidpoint(); - else if (it.second == "cdfjetclu") new_algo = new ClusterAlgoCDFJetClu(); - else if (it.second == "gridjet") new_algo = new ClusterAlgoGridJet(); - else { - ERROR << "Unknown algorithm " << it.second << " ignored." << endmsg; - return true; - } - continue; - } - // clustering algo -> keep the previous syntax - else if (key.find("cluster.")==0) - { - clustering_params.insert(std::pair(key.substr(8),it.second)); - continue; - } - - // Other - try - { - throw EXCEPTION_WARNING("Parameter = "+key+" unknown and thus skipped.","",0); - } - catch(const std::exception& e) - { - MANAGE_EXCEPTION(e); - return false; - } - } - - cluster_collection_.insert(std::pair(new_jetid,new_algo)); - for (const auto &it: clustering_params) - cluster_collection_[new_jetid]->SetParameter(it.first, it.second); - std::string algoname = cluster_collection_[new_jetid]->GetName(); - std::string params = cluster_collection_[new_jetid]->GetParameters(); - INFO << " - Adding Jet ID : " << new_jetid << endmsg; - INFO << " with algo : " << algoname << ", " << params << endmsg; - cluster_collection_[new_jetid]->Initialize(); - } - else if (algorithm == "VariableR") - { - for (std::string key: {"rho", "minR", "maxR", "PTmin", "clustertype", "strategy", "exclusive"}) - { - if (options.find("cluster."+key) == options.end()) - { - ERROR << "Option 'cluster." + key + "' is missing. VariableR clustering ignored." << endmsg; - return true; - } - } - MAfloat32 rho = std::stof(options["cluster.rho"]); - MAfloat32 minR = std::stof(options["cluster.minR"]); - MAfloat32 maxR = std::stof(options["cluster.maxR"]); - MAfloat32 ptmin = std::stof(options["cluster.PTmin"]); - MAbool isExclusive = (options["cluster.exclusive"] == "1"); - - Substructure::VariableR::ClusterType ctype = Substructure::VariableR::AKTLIKE; - if (options["cluster.clustertype"] == "CALIKE") ctype = Substructure::VariableR::CALIKE; - else if (options["cluster.clustertype"] == "KTLIKE") ctype = Substructure::VariableR::KTLIKE; - else if (options["cluster.clustertype"] == "AKTLIKE") ctype = Substructure::VariableR::AKTLIKE; - - Substructure::VariableR::Strategy strategy = Substructure::VariableR::Best; - if (options["cluster.strategy"] == "Best") strategy = Substructure::VariableR::Best; - else if (options["cluster.strategy"] == "N2Tiled") strategy = Substructure::VariableR::N2Tiled; - else if (options["cluster.strategy"] == "N2Plain") strategy = Substructure::VariableR::N2Plain; - else if (options["cluster.strategy"] == "NNH") strategy = Substructure::VariableR::NNH; - else if (options["cluster.strategy"] == "Native") strategy = Substructure::VariableR::Native; - - Substructure::VariableR* variableR; - variableR = new Substructure::VariableR(rho, minR, maxR, ctype, strategy, ptmin, isExclusive); - - substructure_collection_.insert( - std::pair(options["JetID"], variableR) - ); - - std::string exclusive = isExclusive ? "True" : "False"; - INFO << " - Adding Jet ID : " << options["JetID"] << endmsg; - INFO << " with algo : VariableR" << ", " - << "rho = " << options["cluster.rho"] << ", " - << "minR = " << options["cluster.minR"] << ", " - << "maxR = " << options["cluster.maxR"] << ", " - << "ptmin = " << options["cluster.PTmin"] << ", \n" - << " " - << "isExclusive = " << exclusive << ", " - << "clustertype = " << options["cluster.clustertype"] << ", " - << "strategy = " << options["cluster.strategy"] - << endmsg; - } - else - { - ERROR << "Unknown algorithm: " << algorithm << endmsg; - return false; - } - - return true; -#else - ERROR << "FastJet has not been enabled. Can not add jets to the analysis." << endmsg; - return true; -#endif - } + MAbool LoadJetConfiguration(std::map options); /// Accessor to the jet clusterer name std::string GetName() diff --git a/tools/SampleAnalyzer/Test/HEPTopTagger/Test.cpp b/tools/SampleAnalyzer/Test/HEPTopTagger/Test.cpp new file mode 100644 index 00000000..b82dddc4 --- /dev/null +++ b/tools/SampleAnalyzer/Test/HEPTopTagger/Test.cpp @@ -0,0 +1,44 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2012-2022 Jack Araz, Eric Conte & Benjamin Fuks +// The MadAnalysis development team, email: +// +// This file is part of MadAnalysis 5. +// Official website: +// +// MadAnalysis 5 is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// MadAnalysis 5 is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with MadAnalysis 5. If not, see +// +//////////////////////////////////////////////////////////////////////////////// + +#include "SampleAnalyzer/Interfaces/HEPTopTagger/HTT.h" + +using namespace MA5; + +// ----------------------------------------------------------------------- +// main program +// ----------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + std::cout << "BEGIN-SAMPLEANALYZER-TEST" << std::endl; + MA5::Substructure::HTT tagger; + MA5::Substructure::HTT::InputParameters param; + param.do_optimalR = false; + param.reclustering_algorithm = Substructure::kt; + INFO << "initializing HTT " << endmsg; + tagger.Initialize(param); + INFO << "HTT initialized" << endmsg; + tagger.get_settings(); + std::cout << "END-SAMPLEANALYZER-TEST" << std::endl; + return 0; +}