pha2.py

rsps = []
phas = []
baks = []

for interval in intervals:
  
  this_rsp = my_instrument_response_set.weight_by_counts(...)
  this_pha = [this is what you are doing]
  this_bak = [this is what you are doing]
  
  rsps.append(this_rsp)
  phas.append(this_pha)
  baks.append(this_bak)
 
# Now write the RSP file
ebounds = EBOUNDS(...) # You can get what you need from the first matrix
ogip_rsps = [SPECRESP_MATRIX(this_rsp.monte_carlo_energies, this_rsp.ebounds, this_rsp.matrix) for this_rsp in rsps]
rsp2_file = FITSFile(fits_extensions=ogip_rsps)

# and so on