<repo_name>mjadiaz/BL-SSM<file_sep>/HepAid/submit
#!/bin/bash

#SBATCH --job-name=Scan
#SBATCH --ntasks-per-node=16
#SBATCH --time=4:00:00

source /home/mjad1g20/.bashrc
#module load conda/py3-latest
source activate HepAid

#conda env list
python --version
which python
#python BLSSM.py
#python BLSSM_HEscan_DoublePeak.py
python BLSSM_mg5.py<file_sep>/The-model.md
## [ Double Higgs peak in the minimal SUSY B − L model](https://arxiv.org/pdf/1409.7837.pdf).

We want calculate in the BLSSM model, the processes $h^{\prime} \rightarrow Z Z \rightarrow 4 l$ and $h^{\prime} \rightarrow \gamma \gamma$ with the following benchmark in the low energy scale:

![[LowEnergyParameterPoint.png | 450]]


## Useful resources

-  [ Di-photon decay of a light Higgs state in the BLSSM](https://arxiv.org/pdf/2012.04952.pdf)
-  [ Exploring new models in all detail with SARAH](https://arxiv.org/pdf/1503.04200.pdf)
-  [Mass sprecturm for the minimal SUSY *B-L* model ](https://arxiv.org/pdf/1112.4600.pdf)
-  [Beyond-MSSM Higgs sectors](https://arxiv.org/pdf/1409.7182.pdf)
-  [The Higgs sector of the minimal SUSY B − L model](https://arxiv.org/pdf/1504.05328.pdf)<file_sep>/Codes/submit_delphes
#!/bin/bash

#SBATCH --job-name=Delphes
#SBATCH --ntasks-per-node=4


module load conda/py3-latest
source activate root_env

python run_delphes.py
<file_sep>/Codes/read_lhe.py
    #%%
import pylhe
import math
import numpy as np
import pandas as pd
import os
import subprocess
import gzip
import shutil


from compute_processes import process, madgraph_path, param_card_name, param_card_path, output_name
#%%

#process = [['p p > z,( z > l+ l- a $a,( a > l+ l-))', 'ppza4l_test_2']]
def invariant_mass(p1, p2, p3, p4):
    return math.sqrt(
        sum(
            (1 if mu == "e" else -1) * (getattr(p1, mu) + getattr(p2, mu) + getattr(p3, mu) + getattr(p4, mu)) ** 2
            for mu in ["e", "px", "py", "pz"]
        )
    )


def PT(p1):
    return math.sqrt(
        sum((getattr(p1, mu)) ** 2 for mu in ["px", "py"])
    )
# %%


if __name__ == '__main__':

    cut_pt = 5
    cut = True

    # Create directory for store the data
    if not os.path.exists('double_peak_data'):
        os.makedirs('double_peak_data')

    
        

    for p in range(3):

        lhe_path = madgraph_path+output_name+process[p][1]+'/Events/run_01'
        lhe_file  = lhe_path+'/unweighted_events.lhe'

        # Unzip the LHE.gz files 
        with gzip.open(lhe_path+'/unweighted_events.lhe.gz', 'rb') as f_in:
            with open(lhe_file, 'wb') as f_out:
                shutil.copyfileobj(f_in, f_out)

        # Get the invariant mass data for each process
        invariant_mass_vec = np.array([])       #  Vector to store the invariant mass
        for e in pylhe.readLHE(lhe_file):
            # block raising an exception
            m = invariant_mass(e.particles[-1], e.particles[-2],  e.particles[-3], e.particles[-4])
            
            PT_event_list = [PT(e.particles[-1]),PT(e.particles[-2]),PT(e.particles[-3]),PT(e.particles[-4])]
            above_cut_test =  all(x > cut_pt for x in PT_event_list)
            if cut == True:
                if above_cut_test == True:
                    print('all above cut in: ', PT_event_list)
                    invariant_mass_vec=np.append(invariant_mass_vec,m)
                # print('at least one bellow cut in: ', PT_event_list)
                # print(invariant_mass_vec)
            else:
                invariant_mass_vec=np.append(invariant_mass_vec,m)

        df = pd.DataFrame(invariant_mass_vec,columns=['imass'])
        print(df.head(), len(df))
        df.to_csv('double_peak_data/'+process[p][1]+'.csv')
        print(lhe_path, ' Done')


<file_sep>/README.md


### Model Files
The Sarah file was taken from [B-L-SSM – Trac – Hepforge (sarah.hepforge.org)](https://sarah.hepforge.org/trac/wiki/B-L-SSM).  This model was created by L.Basso, F.Staub. The UFO and SPHENO files where created in Mathematica by SARAH by an `generate_model.m` through: 
```mathematica
Needs["SARAH`","/mainfs/home/mjad1g20/HEP/SARAH/SARAH-4.14.5/SARAH.m"]

Start["B-L-SSM"]
ModelOutput[EWSB]

MakeSPheno[]
MakeUFO[]
MakeCHep[]

Quit[];
```


- There is an option to [create the CalcHep](https://gitlab.in2p3.fr/goodsell/sarah/-/wikis/CalcHep/CompHep) file to run it with spheno, that might be helpful. 
	- `RunSPhenoViaCalcHep`: writes C code to run SPhenofrom the graphical interface of CalcHepto calculate the spectrum on the fly.

- There are more options that should be looked.

#### Useful links

- [Running Jupyter Notebook on a remote server](https://docs.anaconda.com/anaconda/user-guide/tasks/remote-jupyter-notebook/)


<file_sep>/Codes/pySpheno.py
import os 
import re
import subprocess

import time




class pySpheno:
    def __init__(self, spheno_dir):
        self._dir=spheno_dir
        self.models = pySpheno._models_in_dir(self)
    
    def _models_in_dir(self):
        models = []
        for f in os.listdir(self._dir+'/bin'):
            if not (re.search(r'SPheno\w+',f)==None):
                m = re.search('(?<=SPheno).*','SPhenoBLSSM')[0]
                models.append(m)
        return models

    def run(self, model, in_file, out_file):
        #subprocess.call([self._dir+'/bin'+'/SPheno'+model, in_file])
       
        run = subprocess.run([self._dir+'/bin'+'/SPheno'+model, in_file, out_file], capture_output=True, text=True)        
        if 'Finished' in run.stdout:
            print(run.stdout)                
        else:
            print('Error')
            

       
        


spheno_dir = '/home/mjad1g20/HEP/SPHENO/SPheno-3.3.8'
in_file = '/mainfs/home/mjad1g20/HEP/SPHENO/SPheno-3.3.8/LesHouches.in.BLSSM_low'
out_file = '/mainfs/home/mjad1g20/HEP/WorkArea/BLSSM-Double-Peak/Codes/spheno_cards/param_card.dat'


spheno = pySpheno(spheno_dir)
spheno.run(model='BLSSM', in_file=in_file, out_file=out_file)


<file_sep>/Codes/submit
#!/bin/bash

#SBATCH --job-name=MadRun
#SBATCH --ntasks-per-node=16

module load conda/py3-latest

#conda env list
source activate root_env

python --version
python compute_processes.py
<file_sep>/Codes/run_delphes.py
import os
import subprocess
import gzip
import shutil

#from compute_processes  import process, output_name, madgraph_path

process=    [['p p > z > l+ l- a $a , a > l+ l-', 'ppza4l'],
             ['p p > h1 > z z, (z > l+ l-, z > l+ l-)', 'pph14l'],
             ['p p > h2 > z z, (z > l+ l-, z > l+ l-)', 'pph24l'],
             ['p p > h1 > a a','pph1aa'],
             ['p p > h2 > a a','pph2aa']]

madgraph_path = '/home/mjad1g20/HEP/MG5/MG5_aMC_v3_1_1'  
output_name = '/double_peak/'

delphes_path = '/home/mjad1g20/HEP/DELPHES/Delphes-3.5.0'
delphes_card = '/mainfs/home/mjad1g20/HEP/DELPHES/Delphes-3.5.0/cards/delphes_card_CMS.tcl'
DelphesHepMC2 = '/mainfs/home/mjad1g20/HEP/DELPHES/Delphes-3.5.0/DelphesHepMC2'


#print(DelphesHepMC2)

if not os.path.exists(delphes_path+output_name):
    print('Create directory for Delphes output: {}'.format(output_name))
    os.makedirs(delphes_path+output_name)

def gunzip(gzfile, gunzipfile):
    with gzip.open(gzfile, 'rb') as f_in:
        with open(gunzipfile, 'wb') as f_out:
            shutil.copyfileobj(f_in, f_out)

for p in process:
    _, p_name = p

    hepmc_path = madgraph_path+output_name+p_name+'/Events/run_01/tag_1_pythia8_events.hepmc'
    delphes_output = delphes_path+output_name+p_name+'.root'

    print(hepmc_path)
    
    if os.path.exists(delphes_output):
        print('Removing existing .root output for the process')
        os.remove(delphes_output)
    
    if not os.path.exists(hepmc_path):
        print('Unzip hepmc.gz for DelphesHepMC2')
        gunzip(hepmc_path+'.gz',hepmc_path)
    
    print('Running DelphesHepMC2')
    #run_delphes = subprocess.run([DelphesHepMC2, delphes_card, delphes_output, hepmc_path], capture_output=True, text=True)
    subprocess.run([DelphesHepMC2, delphes_card, delphes_output, hepmc_path])


    # if os.path.exists(delphes_path+output_name+p_name+'_output.txt'):
    #     os.remove(delphes_path+output_name+p_name+'_output.txt')
    # with open(delphes_path+output_name+p_name+'_output.txt','w') as f:
    #     f.write(run_delphes.stdout)
    #     f.write(run_delphes.stderr)
    
    print('Removing hepmc file ... t0o heavy')
    os.remove(hepmc_path)



    

    <file_sep>/HepAid/BLSSM_mg5.py
import HepRead
import HepTools
import os

model = 'BLSSM'
work_dir = '/scratch/mjad1g20/HEP/WorkArea/BLSSM_Work'
spheno_dir = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8'
reference_lhs = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8/BLSSM/Input_Files/LesHouches.in.BLSSM'
madgraph_dir = '/scratch/mjad1g20/HEP/MG5_aMC_v3_1_1'

# Initiating Madgraph
mg5 = HepTools.Madgraph(madgraph_dir, work_dir)

# Selected param_card.dat
param_card_path = os.path.join(work_dir,'SPhenoBLSSM_output/SPheno.spc.BLSSM_HEscan_770')

# UFO model name in madraph_dir
ufo_model = 'BLSSM-2'

# Creating the Madgraph script
script = HepRead.MG5Script(work_dir,ufo_model)

processes = [['p p > h1, h1 > a a', 'pph1aa'],
             ['p p > h2, h2 > a a', 'pph2aa'],
             ['p p > a a / h1 h2', 'ppaa_bg']]

for p in processes:
    process, syntax = p
    script.process(process)
    script.output(syntax)
    script.launch(syntax)
    script.param_card(param_card_path)
    script.delphes_card('/mainfs/scratch/mjad1g20/HEP/MG5_aMC_v3_1_1/Delphes/cards/delphes_card_CMS.tcl')
    script.write()
    mg5.run()
<file_sep>/HepAid/HepTools.py
import os 
import re
import subprocess

from scipy.stats import loguniform
import matplotlib.pyplot as plt
import numpy as np


import HepRead



class Spheno:
    '''
    ## Spheno class
    -   Identifies what models are available in the work direction. \n 
        and store all the models in .model_list  \n
        
    -   Call .run(in, out) to run Spheno with input and output files. \n
        The output file will be stored in a folder called SPhenoModel_output \n
        inside work_dir.  
    Todo:
    - Figure out how to redirect the additional files that spheno creates in a run.
    '''
    def __init__(self, spheno_dir, work_dir, model=None, input_lhs=None):
        self._dir=spheno_dir
        self.model_list = Spheno._models_in_dir(self)
        self._model = Spheno._model_init(self, model)
        self.work_dir = work_dir
  

    def _model_init(self, model):
        if model == None:
            if len(self.model_list) == 1:
                return self.model_list[0]
            else:
                return print('Define model. \nModels available: {}'.format(self.model_list))
        else:
            print(f'{model} model activated.')
            return model


            
    def _models_in_dir(self):
        models_in_dir = []
        for f in os.listdir(self._dir+'/bin'):
            if not (re.search(r'SPheno\w+',f)==None):
                m = re.search('(?<=SPheno).*',f)[0]
                models_in_dir.append(m)
        return models_in_dir

    def run(self, in_file_name, out_file_name, mode='local'):
        
        out_dir = os.path.join(self.work_dir, 'SPheno'+self._model+'_output')
        in_file = os.path.join(self.work_dir, 'SPheno'+self._model+'_input',in_file_name)

        if not(os.path.exists(out_dir)):
            os.makedirs(out_dir)
        file_dir=os.path.join(out_dir,out_file_name)

        print(f'Save {out_file_name} in :{file_dir}')
        if mode == 'local':
            run = subprocess.run([self._dir+'/bin'+'/SPheno'+self._model, in_file, file_dir], capture_output=True,  text=True)        
            if 'Finished' in run.stdout:
                print(run.stdout) 
                return  file_dir            
            else:
                print('Parameters point error!')
                print(run.stdout)
                return None
        elif mode == 'cluster':
            print('Implement cluster mode')




class Scanner:
    '''
    Class containing useful tools to do scans with SPheno and Madgraph.
    Maybe it will a bunch of functions.
    To do: 
    - Implement general scanner function (for loops)
    - How to store these functions Class or what...
    '''
    
    
    def rlog_array(min, max, n_points, show=False):
        '''
        Creates a random array distributed uniformly in log scale. \n
        - from min to max with n_points number of points. \n
        - turn show=True to see the histogram of the distribution.
        '''
        ar = loguniform.rvs(min, max, size=n_points)

        if show == True:
            plt.hist(np.log10(ar),density=True, histtype='step',color='orange')
            plt.hist(np.log10(ar),density=True, histtype='stepfilled', alpha=0.2, color='orange')
            plt.xlabel('Variable [log]')
            plt.show()
        
        return ar
    def rlog_float(min, max):
        '''
        Creates a random float distributed uniformly in log scale.
        '''
        ar = loguniform.rvs(min, max, size=1)
        return ar[0]


    def runiform_array(min,max, n_points, show=False):
        '''
        Create a uniform random array from min to mas with n_points.
        - show= True to see the distribution histogram.
        '''
        ar=np.random.uniform(min, max ,size = n_points)

        if show == True:
            plt.hist(ar,density=True, histtype='step',color='orange')
            plt.hist(ar,density=True, histtype='stepfilled', alpha=0.2, color='orange')
            plt.xlabel('Variable')
            plt.show()
        
        return ar

    def runiform_float(min,max):
        '''
        Create a uniform random float within min and max.
        '''
        ar=np.random.uniform(min, max ,size = 1)        
        return ar[0]


class Madgraph:
    '''
    Basicaly just run madgraph with in the script mode with the scripts created \n
    by the MG5Script class. \n
    Todo:
    - Figure out how to print some output from madgraph.
    '''
    def __init__(self, madgraph_dir, work_dir):
        self._dir = madgraph_dir
        self.work_dir = work_dir       
        
    def run(self, input_file = 'MG5Script.txt', mode='local'):
        '''
        Run madgraph with an script named MG5Script.txt (created by the MG5Script class) in within work_dir. \n
        Change input_file to change to another script within work_dir.
        '''
        if mode == 'local':
            subprocess.run([os.path.join(self._dir,'bin/mg5_aMC'), os.path.join(self.work_dir,input_file)])    
        elif mode == 'cluster':
            print('Implement cluster mode.')
            
class HiggsBounds:
    def __init__(self, higgs_bounds_dir, work_dir, model=None, neutral_higgs=None, charged_higgs=None):
        self._dir = higgs_bounds_dir
        self.work_dir = work_dir
        self.model = model
        self.neutral_higgs = neutral_higgs
        self.charged_higgs = charged_higgs
           
        
    def run(self):
        '''
        Runs HiggsBounds with the last point calculated by SPheno. Returns the path for HiggsBounds_results.dat if SPheno and HiggsBounds succeed. If there is an Error in one them it returns None and prints the error output.
        '''
        run = subprocess.run([os.path.join(self._dir, 'HiggsBounds'), 'LandH', 'effC', str(self.neutral_higgs), str(self.charged_higgs), self.work_dir+'/'], capture_output=True,  text=True)        
        if 'finished' in run.stdout:
            print(run.stdout) 
            return  os.path.join(self.work_dir, 'HiggsBounds_results.dat')            
        else:
            print('HiggsBound not finished!')
            print(run.stdout)
            return None

<file_sep>/Codes/events_invariant_mass_plot.py
#%%
import seaborn as sns # for data visualization
import pandas as pd # for data analysis
import matplotlib.pyplot as plt # for data visualization
import numpy as np

from compute_processes import process, process_1, madgraph_path, param_card_name, param_card_path, output_name

# %%
sns.set(style='darkgrid',)
sns.__version__

# %%
# df_mass1 = pd.read_csv('double_peak_data/'+process[0][1]+'.csv')
# print(process[0][1])
df_mass1 = pd.read_csv('double_peak_data/'+process_1[0][1]+'.csv')
print(process_1[0][1])

df_mass2 = pd.read_csv('double_peak_data/'+process[1][1]+'.csv')
print(process[1][1])

df_mass3 = pd.read_csv('double_peak_data/'+process[2][1]+'.csv')
print(process[2][1])

#%%
ppza4l = df_mass1.imass.to_numpy()
pph14l = df_mass2.imass.to_numpy()
pph24l = df_mass3.imass.to_numpy()

pph14l_vec = [r'$pp\rightarrow h1 \rightarrow ZZ \rightarrow 4l$']*len(pph14l)
pph24l_vec = [r'$pp \rightarrow h2 \rightarrow ZZ \rightarrow 4l$']*len(pph24l)
ppza4l_vec  = [r'$pp \rightarrow Z \rightarrow 2l\gamma * \rightarrow 4l$']*len(ppza4l)
df1 = pd.DataFrame(dict(m = pph14l, channel = pph14l_vec), columns=['m', 'channel'])
df2 = pd.DataFrame(dict(m = pph24l, channel = pph24l_vec), columns=['m', 'channel'])
df3 = pd.DataFrame(dict(m = ppza4l, channel = ppza4l_vec), columns=['m', 'channel'])

df = pd.concat([df1,df2,df3])
df
# %%
sns.histplot(
    df, x="m", hue="channel", element="step",
)

plt.xlim(0,200)
#plt.yscale('log')
#plt.xlabel(r'$M_{\mu^{+}\mu^{-}} $')
plt.xlabel(r'$M_{4ll} $')
plt.ylabel('Events')
plt.show()
# %%


# %%
<file_sep>/HepAid/HepRead.py
import re 
import os
from shutil import copy
import awkward as ak
import numpy as np
from tabulate import tabulate



#########################################
# Classes for reading LesHouches files. #
# Focusing on Spheno.                   #
#########################################

class BlockLine:
    '''## Line:
        Each line can be of a different category:
        - block_header
        - on_off
        - value
        - matrix_value '''
    def __init__(self, entries, line_category):
        self.entries = entries
        self.line_category = line_category
        self.line_format = BlockLine.fline(line_category)
    def fline(cat):
        ''' Text format of each line'''
        if cat == 'block_header':
            line_format = '{:6s} {:20s}  {:13s}'
        elif cat == 'on_off':
            line_format = '{:6s} {:18s}  {:13s}'
        elif cat == 'value':
            line_format = '{:6s} {:18s}  {:13s}'
        elif cat == 'matrix_value':
            line_format = '{:3s}{:3s} {:18s}  {:13s}'
        return line_format


class Block:
    '''
    ## Block
    It holds each line of a block.\n
    Call .show() to print a block. \n
    Call .set(parameter_number, value) to change the parameter value in the instance.
    '''
    def __init__(self, block_name, block_comment=None, category=None):
        self.block_name = block_name
        self.block_comment = block_comment
        self.block_body = []
        self.category = category

    def show(self):
        '''
        Print block information in the LesHouches format.
        '''
        print('{} {}   {:10s}'.format('Block',self.block_name,self.block_comment))
        for b in self.block_body:
            print(b.line_format.format(*b.entries))
    
        
    
    def set(self,option, param_value):
        '''
        Call set(option, param_value) method to modify the option N with a parameter_value \n.
        -option = Can be an int or a list [i, j]. \n
        -param_value = Can be an int (on/off) or a float. \n
        '''
        for line in self.block_body:
            if (line.line_category == 'matrix_value'):
                if (option == [int(line.entries[0]),int(line.entries[1])]):
                    line.entries[2] = '{:E}'.format(param_value)
                    print('{} setted to : {}'.format(line.entries[-1], line.entries[1]))
                    break
            if (line.line_category == 'value') & (option == int(line.entries[0])):              
                line.entries[1] = '{:E}'.format(param_value) 
                print('{} setted to : {}'.format(line.entries[-1], line.entries[1])) 
                break                  
            elif (line.line_category == 'on_off') & (option == int(line.entries[0])):
                if isinstance(param_value, int):
                    line.entries[1] = '{}'.format(param_value)
                    print('{} setted to : {}'.format(line.entries[-1], line.entries[1]))
                else:
                    print('param_value={} is not integer'.format(param_value))
                break



class LesHouches:
    '''
    ## LesHuches
    Read a LesHouches file and stores each block in block classes. \n
    - To get all the names of the blocks call .block_list. \n
    - work_dir is the directory where all the outputs will be saved. \n
    - The new LesHouches files will be saved in a folder called SPhenoMODEL_input since is the input for spheno.
    '''
    def __init__(self, file_dir, work_dir, model):
        self.file_dir = file_dir
        print(f'Reading LesHouches from : {file_dir}')
        self._blocks = LesHouches.read_leshouches(file_dir)
        self.block_list = [name.block_name for name in self._blocks]
        self.work_dir = work_dir
        self.model = model
    

    def block(self, name):
        block = LesHouches.find_block(name.upper(), self._blocks)
        return block  


    def find_block(name, block_list):
            try:
                if isinstance(name, str):
                    for b in block_list:
                        if b.block_name == name:
                            b_found = b
                            break
                        else:
                            None
                return b_found
            except:
                print('block not found')

    def read_leshouches(file_dir):
        block_list = []
        paterns =   dict(   block_header= r'(?P<block>BLOCK)\s+(?P<block_name>\w+)\s+(?P<comment>#.*)',
                            on_off= r'(?P<index>\d+)\s+(?P<on_off>-?\d+\.?)\s+(?P<comment>#.*)',
                            value= r'(?P<index>\d+)\s+(?P<value>-?\d+\.\d+E.\d+)\s+(?P<comment>#.*)',
                            matrix_value= r'(?P<i>\d+)\s+(?P<j>\d+)\s+(?P<value>-?\d+\.\d+E.\d+)\s+(?P<comment>#.*)')



        with open(file_dir, 'r') as f:
            for line in f:
                # Match in 3 groups a pattern like: '1000001     2.00706278E+02   # Sd_1'
                m_block = re.match(paterns['block_header'], line.upper().strip()) 
                if not(m_block == None):

                    if m_block.group('block_name') in ['MODSEL','SPHENOINPUT','DECAYOPTIONS']:
                        block_list.append(Block(    block_name=m_block.group('block_name'), 
                                                    block_comment=m_block.group('comment'),
                                                    category= 'spheno_data' ))
                        in_block, block_from = m_block.group('block_name'), 'spheno_data'

                    else:
                        block_list.append(Block(        block_name=m_block.group('block_name'), 
                                                        block_comment=m_block.group('comment'),
                                                        category= 'parameters_data'))
                        in_block, block_from = m_block.group('block_name'), 'parameters_data'

                m_body =  re.match(paterns['on_off'], line.strip())
                if not(m_body == None):            
                    LesHouches.find_block(in_block,block_list).block_body.append(BlockLine(list(m_body.groups()),'on_off'))
                m_body =  re.match(paterns['value'], line.strip())
                if not(m_body == None):            
                    LesHouches.find_block(in_block,block_list).block_body.append(BlockLine(list(m_body.groups()),'value'))
                m_body =  re.match(paterns['matrix_value'], line.strip())
                if not(m_body == None):            
                    LesHouches.find_block(in_block,block_list).block_body.append(BlockLine(list(m_body.groups()), 'matrix_value'))
        return block_list

  
    def new_file(self, new_file_name):
        '''
        Writes a new LesHouches file with the blocks defined in the instance. \n
        Possibly with new values for the parameters and options.
        '''
        new_file_dir = os.path.join(self.work_dir, 'SPheno'+self.model+'_input')

        if not(os.path.exists(new_file_dir)):
            os.makedirs(new_file_dir)
        file_dir=os.path.join(new_file_dir,new_file_name)
        print(f'Writing new LesHouches in :{file_dir}')
        
        with open(file_dir,'w+') as f:
            for block in self._blocks:
                head = '{} {}  {:10s}'.format('Block',block.block_name,block.block_comment)+'\n'
                f.write(head)
                for b in block.block_body:
                    f.write(b.line_format.format(*b.entries)+'\n')

################################### 
# Classes for reading SLHA files. #  
# Focusing on madgraph.           #
###################################

class Particle():
    def __init__(self, pid, ufo_name=None, total_width=None, mass=None, comment=None, decays=None):
        self.pid = pid
        self._ufo_name = ufo_name  # To link it with madgraph but maybe it is not necessary.
        self.total_width = total_width
        self.mass = mass
        self.comment = comment 

        # How to store the decays
        self.decays = decays
    
    def find_particle(id, particles_list):
        try:
            if isinstance(id, int):        
                for p in particles_list:
                    if p.pid == id:
                        p_found = p
                        break
                    else:
                        None
            elif isinstance(id, str):
                for p in particles_list:
                    if p.comment == id:
                        p_found = p
                        break
                    else:
                        None
            return p_found
        except:
            print('No particle found')

    def show(self):
        #particle = Particle.find_particle(pid, particles_list)
        n_data = lambda n: [self.decays[n][i] for i in self.decays.fields]
        n_decays = ak.num(self.decays,axis=0)
        data = [n_data(n) for n in range(n_decays)]
        headers=[field.upper() for field in self.decays.fields]
        print(tabulate([[self.mass,self.total_width]], headers=['MASS GeV', "TOTAL WIDTH GeV"],tablefmt="fancy_grid"))
        print(tabulate(data, headers=headers,tablefmt="fancy_grid"))


class Slha:
    '''
    ## Read SLHA
    Read the given model file in SLHA format. It stores PID of all the particles of the models in .model_particles \n
    in a list of tuples (PID, NAME) to have a sense of what particles are in the model. \n
    The internal list _particles contain each particle saved as a Particle class so that we can use \n
    all the internal methods and properties of this class, for example: 
    - .particle('25').show() ---> It display all the information about the particle.
    - .particle('25').mass    
    \n
    This class focuses on extracting information not writing the file because we always can \n
    change the parameters inside madgraph by the command: set parameter = value.

    '''
    def __init__(self, model_file):
        self.model_file = model_file
        self._particles, self.model_particles = Slha.read_file(model_file)
    
    def read_file(model_file):
        '''
        Read all the particles of the model saving it in the particles object list. \n
        It creates a Particle class for each particle of the model reading it \n
        from the MASS block of the SLHA file.
        '''
        with open(model_file, 'r') as f:
            slha_file = f.read() # Read the whole file
            slha_blocks = re.split('BLOCK|Block|block|DECAY|Decay|decay', slha_file) #  Separate the text file into parameter blocks and decay blocks
            
            particles = []      # Initiate the Partcle objects list.
            

            for block in slha_blocks:
                block = re.split('\n',block)    # Split block in lines
                # Create the particles of the model from the MASS block and populate the
                # particles list with Particles objects.
                if "MASS" in block[0].split()[0].upper():
                    for line in block:
                        line = re.split('\s+',line.strip()) 
                        #print(l)
                        if line[0].isnumeric():                            #  To skip to the first particle information, since the first element is the PID (numeric)
                            particles.append(Particle(int(line[0])))       #  Append to particles list a Particle object with the PID given by l[0]
                            particles[-1].mass = float(line[1])            #  Assign the mass for the given Particle object
                            particles[-1].comment = line[3]                #  Assign the ufo_name wich is the comment (last element of the line).
                                                                            #  Note: Not always is the ufo_name.
                            #print(particles[-1].pid, particles[-1].mass, particles[-1].comment) # Just to check
                    break                

            # Another foor loop because we want to create the particles first            
            for block in slha_blocks:
                block = re.split('\n',block)    # Split block in lines again            

                # If the first element on the first line is numeric correspond to decay 
                # block with that ID value else is a parameter block.

                if block[0].split()[0].isnumeric(): 
                    tmp_pid = int(block[0].split()[0]) # Temporal particle id for a decay block
                    decays = ak.ArrayBuilder() #  Initiate the decays as an awkward array builder for the current block

                    for line in block: 


                        # Match in 3 groups a pattern like: '1000001     2.00706278E+02   # Sd_1'
                        match = re.match(r'(?P<pid>\d+)\s+(?P<width>\d+\.\d+E.\d+)\s+#(?P<comment>.*)', line.strip()) 
                        if not(match == None):
                            #print('PID: {}, Width: {:E}, Comment: {}'.format(int(match.group('pid')),float(match.group('width')) ,match.group("comment").strip()))
                            Particle.find_particle(tmp_pid,particles).total_width = float(match.group('width'))
                            #print(Particle.find_particle(tmp_pid,particles).pid, Particle.find_particle(tmp_pid,particles).total_width)
                        
                        # Match in 4 groups a pattern like: '2.36724485E-01    2            2   -1000024   # BR(Sd_1 -> Fu_1 Cha_1 )'
                        match = re.match(r'(?P<br>\d+\.\d+E.\d+)\s+(?P<nda>\d+)\s+(?P<fps>.*?)\s+#(?P<comment>.*)', line.strip())

                        if not(match == None):
                            tmp_decay = ak.Array([{ "br": float(match.group('br')),
                                                    "nda": int(match.group('nda')), 
                                                    "fps": [int(p) for p in match.group('fps').split()], 
                                                    "comment": match.group('comment').strip()}])
                            decays.append(tmp_decay, at=0)

                        #print(tmp_decay)                    
                    Particle.find_particle(tmp_pid,particles).decays = decays.snapshot()

        def particle_ids(particles_list):
            ''' Returns a simple list holding all the particle ID of the model'''
            pid_list=[(p.pid, p.comment) for p in particles_list]
            #comment_list = [p.comment for p in particles_list]
            #pid_list.sort()
            return pid_list

        particles_list = particle_ids(particles) 

        return particles, particles_list



    def particle(self, pid):
        particle = Particle.find_particle(pid, self._particles)
        return particle

#######################################
# Class for writing a Madgraph script #
#######################################

class MG5Script:
    '''
    Create a object containing mostly* all the necesary commands to compute a process in madgraph with the \n
    text-file-input/script mode. A default or template input file can be preview with the show() method within this class.
    '''
    def __init__(self,work_dir, ufo_model):
        self.work_dir = work_dir
        self.ufo_model = ufo_model 
        self._default_input_file() 

    def import_model(self):
        out = 'import model {} --modelname'.format(self.ufo_model)
        self._import_model = out
      
        
    def define_multiparticles(self,syntax = ['p = g d1 d2  u1 u2  d1bar d2bar u1bar u2bar', 'l+ = e1bar e2bar', 'l- = e1 e2']):
        out = []
        if not(syntax==None):
            [out.append('define {}'.format(i)) for i in syntax]         
            self._define_multiparticles = out            
        else:
            self._define_multiparticles = None

    def process(self,syntax = 'p p > h1 > a a'):
        '''
        Example: InputFile.process('p p > h1 > a a')
        '''
        out = 'generate {}'.format(syntax)
        self._process = out

        
    
    def add_process(self,syntax = None):
        '''
        Example: InputFile.add_process('p p > h2 > a a')
        '''
        out = []
        if not(syntax==None):
            [out.append('add process {}'.format(i)) for i in syntax]         
            self._add_process= out       
        else:
            self._add_process= None
    
    def output(self, name='pph1aa'):
        output_dir = os.path.join(self.work_dir,name)
        out = 'output {}'.format(output_dir)
        self._output = out

        

    def launch(self, name='pph1aa'):
        launch_dir = os.path.join(self.work_dir,name)
        out = 'launch {}'.format(launch_dir)
        self._launch = out

        
    def shower(self,shower='Pythia8'):
        '''
        Call .shower('OFF') to deactivate shower effects.
        '''
        out ='shower={}'.format(shower)
        self._shower = out
 
         
    def detector(self,detector='Delphes'):
        '''
        Call .detector('OFF') to deactivate detector effects.
        '''
        out ='detector={}'.format(detector)
        self._detector = out

               
    def param_card(self, path=None):
        if path==None:
            self._param_card = None
        else:
            self._param_card = path
 
            
    def delphes_card(self, path=None):
        if path==None:
           self._delphes_card = None
        else:
            self._delphes_card = path

            
    def set_parameters(self,set_param=None):
        if set_param==None:
            self._set_parameters = None
        else:
            out = ['set {}'.format(i) for i in set_param]
            self._set_parameters = out

            
    def _default_input_file(self):
        self.import_model()
        self.define_multiparticles()
        self.process()
        self.add_process()
        self.output()
        self.launch()
        self.shower()
        self.detector()
        self.param_card()
        self.delphes_card()
        self.set_parameters()

    def show(self):
        '''
        Print the current MG5InputFile
        '''
        write = [self._import_model, self._define_multiparticles, self._process, self._add_process, self._output, self._launch, self._shower,
        self._detector, '0', self._param_card, self._delphes_card, self._set_parameters, '0'] 
        for w in write:
            if not(w==None):
                if isinstance(w,str):
                    print(w)
                elif isinstance(w,list):
                    [print(i) for i in w]
    def write(self):
        '''
        Write a new madgraph script as MG5Script.txt used internally by madgraph.
        '''
        write = [self._import_model, self._define_multiparticles, self._process, self._add_process, self._output, self._launch, self._shower,
        self._detector, '0', self._param_card, self._delphes_card, self._set_parameters, '0']
        f = open(os.path.join(self.work_dir,'MG5Script.txt'),'w+')
        for w in write:
            if not(w==None):
                if isinstance(w,str):
                    f.write(w+'\n')
                elif isinstance(w,list):
                    [f.write(i+'\n') for i in w]
        f.close()
        return 
        

##########################################
# Class for reading a HiggsBounds output #
##########################################


class HiggsBoundsResults:
    def __init__(self, work_dir, model=None):
        self.work_dir = work_dir
        self.model = model
    
    def read(self, direct_path=None):
        '''
        Read HiggsBounds_results.dat and outputs a tuple with: \n
        (HBresult, chan, obsratio, ncomb)
        '''
        if direct_path:
            file_path = direct_path
        else:
            file_path = os.path.join(self.work_dir, 'HiggsBounds_results.dat')

        with open(file_path , 'r') as f:
            for line in f:
                line = line.strip()
                if ('HBresult' in line) & ('chan' in line):
                    line = line.split()[1:]
                    index_hbresult = line.index('HBresult')
                    index_chan = line.index('chan')
                    index_obsratio = line.index('obsratio')
                    index_ncomb = line.index('ncomb')
                    #print(line)
                    _ = f.readline()
                    for subline in f:
                        subline = subline.split()
                        #print(subline[index_hbresult],subline[index_chan],subline[index_obsratio],subline[index_ncomb])
                        break
        return (subline[index_hbresult],subline[index_chan],subline[index_obsratio],subline[index_ncomb])

    def save(self, output_name,in_spheno_output = True):
        if in_spheno_output:
            copy(os.path.join(self.work_dir, 'HiggsBounds_results.dat'), os.path.join(self.work_dir,'SPheno'+self.model+'_output' ,'HiggsBounds_results_'+str(output_name)+'.dat'))
        pass                

<file_sep>/Codes/compute_processes.py
#%%
import os
import subprocess
import random
#%%
############################
# Paths and processes info #
############################

# MadGraph process syntax and launch/output name
#process=    [['p p > h1 > z z, z  > l+ l-,z > l+ l-', 'pph12z4l']]
#process=    [['p p > h2 > a a', 'pph2aa']]
process=    [['p p > h1 > a a','pph1aa'],
             ['p p > h2 > a a','pph2aa'],
             ['p p > h1 > z z, (z > l+ l-, z > l+ l-)', 'pph14l'],
             ['p p > h2 > z z, (z > l+ l-, z > l+ l-)', 'pph24l'],
             ['p p > z > l+ l- a $a , a > l+ l-', 'ppza4l']]
# wget https://launchpad.net/mg5amcnlo/3.0/3.1.x/+download/MG5_aMC_v3.1.1.tar.gz

madgraph_path = '/scratch/mjad1g20/HEP/MG5_aMC_v3_1_1'  


param_card_name_2 = '/SPheno.spc.BLSSM'
param_card_path_2 = '/home/mjad1g20/HEP/SPHENO/SPheno-3.3.8'
param_card_name = '/param_card.dat'
param_card_path = '/home/mjad1g20/HEP/WorkArea/BLSSM-Double-Peak/Codes/Cards'

delphes_card = '/mainfs/home/mjad1g20/HEP/WorkArea/BLSSM-Double-Peak/Codes/Cards/delphes_card_CMS.dat'

output_name = '/double_peak/'

#################################
# Compute processes in Madgraph #
#################################
if __name__ == '__main__':

    #for p in range(len(process)):
    for p in process:
        mg5_syntax, p_name = p

        

        f = open(madgraph_path+'/compute_process.txt',"w+")
        f.write('import model BLSSM-2 --modelname\n')
        f.write('define p = g d1 d2  u1 u2  d1bar d2bar u1bar u2bar\n')
        f.write('define l+ = e1bar e2bar\n')
        f.write('define l- = e1 e2\n')
        f.write('generate '+mg5_syntax+'\n')

        while os.path.exists(madgraph_path+output_name+p_name):
            p_name = p_name+'_'+str(random.randint(1,100))

        f.write('output '+madgraph_path+output_name+p_name+'\n')
        f.write('launch '+madgraph_path+output_name+p_name+'\n')
        f.write('shower=Pythia8\n')
        f.write('detector=Delphes\n')
        f.write('0\n')
        f.write(param_card_path_2+param_card_name_2+'\n')
        f.write(delphes_card+'\n')
        f.write('set mh1 = 125\n')
        f.write('set mh2 = 136\n')
        f.write('set ebeam1 = 4000.0\n')
        f.write('set ebeam2 = 4000.0\n')
        f.write('0\n')
        f.close()

        subprocess.call([madgraph_path+'/bin/mg5_aMC', madgraph_path+'/compute_process.txt'])
        break
    
    os.remove(madgraph_path+'/compute_process.txt')
    
<file_sep>/HepAid/BLSSM_HEscan.py
import numpy as np

import HepRead
import HepTools
import os

from HepTools import Scanner as sc

model = 'BLSSM'
work_dir = '/scratch/mjad1g20/HEP/WorkArea/BLSSM_Work'
spheno_dir = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8'
reference_lhs = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8/BLSSM/Input_Files/LesHouches.in.BLSSM'
madgraph_dir = '/scratch/mjad1g20/HEP/MG5_aMC_v3_1_1'


SPhenoBLSSM = HepTools.Spheno(spheno_dir, work_dir ,model)
lhs  = HepRead.LesHouches(reference_lhs, work_dir, model)

#print('Blocks in High Energy LesHouches:')
#print('---------------------------------')
#print(lhs.block_list)

#lhs.block('MODSEL').show()
#lhs.block('SMINPUTS').show()

# Fixed parameters
lhs.block('MINPAR').set(8,2500)
lhs.block('MINPAR').set(7,1.15)
lhs.block('MINPAR').set(3,5)
lhs.block('MINPAR').show()
[lhs.block('extpar').set(i,0) for i in [11,12,13,14]]
lhs.block('extpar').set(3,2500)
lhs.block('extpar').set(10,1.15)
lhs.block('extpar').set(1,0.55)
lhs.block('extpar').set(2,-0.12)
lhs.block('extpar').show()

option = [13, 16,78,520,67,60] 
value  = [3,0,1,0,1,1]
for i,j in zip(option, value):
    lhs.block('sphenoinput').set(i,j)
#lhs.block('sphenoinput').show()

option = [1,2,3,4,5,6,7,8,9,1001,1002,1003,1004,1005,1006,1007,1008,1009,1010,1011,1012,1013] 
value  = np.zeros(len(option),dtype = int)
for i,j in zip(option, value):
    lhs.block('decayoptions').set(i,int(j))
                                  
#lhs.block('decayoptions').show()


# Scanned parameters 
m0 = lambda low, high: sc.runiform_float(low,high)
m12 = lambda low, high: sc.runiform_float(low,high)
a0 = lambda low, high: sc.runiform_float(low,high)
tanbeta = lambda low, high: sc.runiform_float(low,high)
N_POINTS = 1300

points = 0
index = []
for i in range(N_POINTS):
    
    lhs.block('MINPAR').set(1,m0(100,1000))
    lhs.block('MINPAR').set(2,m12(1000,4500))
    lhs.block('MINPAR').set(5,a0(1000,4000))
    lhs.block('MINPAR').set(3,tanbeta(1,60))
    new_lhs_name = 'LesHouches.in.BLSSM_HEscan_'+str(i)
    out_spheno_name = 'SPheno.spc.BLSSM_HEscan_'+str(i)
    # Creating the new LesHouches file, and runing Spheno.
    lhs.new_file(new_lhs_name)
    param_card_path = SPhenoBLSSM.run(new_lhs_name,out_spheno_name)
    # Reading spheno output
    if not(param_card_path==None):
        points = points + 1
        index.append(i)
        
        
        
with open('number_of_points.txt','w+') as f:
    f.write(f'{points} \n')
    for i in index:
        f.write(f'{i} \n')
    
        





<file_sep>/HepAid/BLSSM_HEscan_DoublePeak.py
import numpy as np

import HepRead
import HepTools
import os

from HepTools import Scanner as sc

model = 'BLSSM'
work_dir = '/scratch/mjad1g20/HEP/WorkArea/BLSSM_DoublePeak'
spheno_dir = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8'
reference_lhs = '/scratch/mjad1g20/HEP/SPHENO/SPheno-3.3.8/BLSSM/Input_Files/LesHouches.in.BLSSM'
madgraph_dir = '/scratch/mjad1g20/HEP/MG5_aMC_v3_1_1'


SPhenoBLSSM = HepTools.Spheno(spheno_dir, work_dir ,model)
lhs  = HepRead.LesHouches(reference_lhs, work_dir, model)


# Fixed parameters
option, value = [8,7,3], [1700,1.15,5]
for i,j in zip(option, value):
    lhs.block('MINPAR').set(i,j)

[lhs.block('EXTPAR').set(i,0) for i in [11,12,13,14]]
option, value = [3,10,1,2], [1700,1.15,0.55,-0.12]
for i,j in zip(option, value):
    lhs.block('EXTPAR').set(i,j)

# Spheno options to run with Madgraph
option = [13, 16,78,520,67,60] 
value  = [3,0,1,0,1,1]
for i,j in zip(option, value):
    lhs.block('sphenoinput').set(i,j)

# Turning off loop decay options because they are already defined
option = [1,2,3,4,5,6,7,8,9,1001,1002,1003,1004,1005,1006,1007,1008,1009,1010,1011,1012,1013] 
value  = np.zeros(len(option),dtype = int)
for i,j in zip(option, value):
    lhs.block('decayoptions').set(i,int(j))
                                  



# Scanned parameters 
m0 = lambda low, high: sc.runiform_float(low,high)
m12 = lambda low, high: sc.runiform_float(low,high)
a0 = lambda low, high: sc.runiform_float(low,high)
N_POINTS = 1000

points = 0
index = []
for i in range(N_POINTS):
    
    lhs.block('MINPAR').set(1,m0(100,1000))
    lhs.block('MINPAR').set(2,m12(2000,4500))
    lhs.block('MINPAR').set(5,- a0(1000,4000))
    new_lhs_name = 'LesHouches.in.BLSSM_HEscan_'+str(i)
    out_spheno_name = 'SPheno.spc.BLSSM_HEscan_'+str(i)
    # Creating the new LesHouches file, and runing Spheno.
    lhs.new_file(new_lhs_name)
    param_card_path = SPhenoBLSSM.run(new_lhs_name,out_spheno_name)
    # Reading spheno output
    if not(param_card_path==None):
        points = points + 1
        index.append(i)
        
        
        
with open('number_of_points.txt','w+') as f:
    f.write(f'{points} \n')
    for i in index:
        f.write(f'{i} \n')
    
        





