# coding: utf-8 # # MAGIC Source Simulator. The macro takes a given spectral shape, # and using MAGIC performance in the range 40GeV-16TeV from # Aleksic, J., et al., 2016, Astroparticle Physics, 72, 76 # (please cite this paper if you want to use this macro in a publication) # estimates what kind of signal can be observed by MAGIC. # It computes significances of each spectral point according to # Li & Ma 1983, ApJ, 272, 317, Eq. 17 definition. # # Terminal Output: # For each estimated energy range one gets the number of excess events, # signal-to-background ratio, significance, and information if a given # bin satisfies the conditions for the detection. # # Plot Output: # Spectral points following the assumed spectrum are shown, their # error bars reflect the performance of the instrument for such a source. # Significance for each bin is given (bins without detection # have gray numbers. In the top of the plot a simple number using # information from all the shown bins is given to evaluate # if the source can be detected (roughly if sum of sigma_i / sqrt(N) >~5) # # requirements: # - python3, matplotlib, numpy, scipy # # caveats: # - the macro is operating on estimated energy doing simple # comparisons with differential in estimated energy rates seen for # Crab Nebula, for soft sources the differences in energy migration # will result in different performance than calculated here # - the threatment of extended sources is very approximate, only # increase in background is taken into account (without energy # dependence of PSF) for extended sources for extension >~0.4 deg the # dependence on the offset from the centre of the camera will further # worsen the performance w.r.t. the one calculated here # - significances are given for each differential energy bin # separatelly, but to detect a source one normally applies a cut that # keeps a broad range of energies inside resulting in better integral # sensitivity than differential one. Also, optimization of cuts for # a broad energy range usually results in somewhat better sensitivity # than what one can get by simply integrating the used here signal in # differential energy bins. As a very crude approximation for detection # capability we calculate here also a sum of significances of all the # points in SED divided by the sqrt of the number of those points # # usage: # - modify the "Assumed(x)" spectrum and "basic settings" below # and run in terminal: python3 mss.py # - The default output is a .pdf figure file. # - For display of the graph, uncomment fig.show() line. # WARNING: On MacOS fig.show() might hang unless using pythonw: # pythonw mss.py # # for questions ask Julian Sitarek (jsitarek [at] uni.lodz.pl) # for python version add: ana.babic@fer.hr # version history: # 1.1 (2017/10/16): fixed a bug which made the errorbars too small import matplotlib.pyplot as plt import numpy as np from scipy.integrate import quad # Assumed spectrum # [cm^-2 s^-1 TeV^-1] assumed source spectrum (dN/dS dt dE) as a function of energy in GeV # Some common examples below def Assumed(x): return 0.1*np.exp(-x/200)*3.39e-11*pow(x/1000., -2.51-0.21*np.log10(x/1000.)) #10% Crab with 200 GeV cut-off # return 1.e-11 * pow(x/1000, -1.9) # simple power-law # return 1.e-11 * pow(x/1000, -1.9) * np.exp(-x/300) # power-law with an exponential cut-off # return 1.e-11 * pow(x/1000, -1.9) * np.exp(-pow(x/300, 0.5)) # power-law with a non-exponential cut-off # basic settings timeh=50 # [h], time of observations extension=0. # [deg] extension radius of the source srcname ="SourceTest1" # name of the source ismidzd = False # false: performance for 0-30deg zenith angles, true: for 30-45deg zenith # you can check visibility of your source e.g. here: http://www.magic.iac.es/scheduler/ # advanced settings numoff = 3 # number of background estimation regions minev = 10.0 # minimum number of events minSBR = 0.05 # minimum ratio of excess to background PSF = 0.1 # [deg] PSF (for worsening the performance for extended sources) eplotmin = 40 # [GeV] x plot range eplotmax = 20.e3 # [GeV] x plot range yplotmin = 1.e-14 # [TeV cm^-2 s^-1] y plot range yplotmax = 1.e-9 # [TeV cm^-2 s^-1] y plot range minerror = 2 # showing only points with value > minerror * error drawsigma = True # whether to draw also sigmas on the plot #Crab Nebula [Aleksic et al. 2016] def Crab(x): return 3.39e-11*pow(x/1000., -2.51-0.21*np.log10(x/1000.)) # Calculate SED value from energy and flux, for display def CalcSED(x,flux): return 1.e-6*x*x*flux # global variables (DO NOT MODIFY) npoints = 13 def Prepare(): npoints_array = np.arange(0, npoints+1) enbins = 100.*pow(10., (npoints_array-2)*0.2) # [GeV] if ismidzd: crabrate = np.array([0,0.404836, 3.17608, 2.67108, 2.86307, 1.76124,\ 1.43988, 0.944385, 0.673335, 0.316263, 0.200331,\ 0.0991222, 0.0289831]) # [min^-1] bgdrate = np.array([1.67777, 2.91732, 2.89228, 0.542563, 0.30467, 0.0876449,\ 0.0375621, 0.0197085, 0.0111295, 0.00927459, 0.00417356,\ 0.00521696, 0.000231865]) # [min^-1] else: # 0-30 deg values crabrate = np.array([0.818446, 3.01248, 4.29046, 3.3699, 1.36207, 1.21791,\ 0.880268, 0.579754, 0.299179, 0.166192, 0.0931911,\ 0.059986, 0.017854]) # [min^-1] bgdrate = np.array([3.66424, 4.05919, 2.41479, 0.543629, 0.0660764, 0.0270313,\ 0.0132653, 0.00592351, 0.00266975, 0.00200231, 0.00141831,\ 0.00458864, 0.0016686]) # [min^-1] return npoints_array, enbins, crabrate, bgdrate #// Li & Ma eq 17. significance, function abridged from MARS def SignificanceLiMa(s, b, alpha): if b==0 and s==0: return 0 b = FLT_MIN if b==0 else b #// Guarantee that a value is output even in the case b or s == 0 s = FLT_MIN if s==0 else s #// (which is not less allowed, possible, or meaningful than #// doing it with b,s = 0.001) sumsb = s+b if sumsb<0 or alpha<=0: return -1; l = s*np.log(s/sumsb*(alpha+1)/alpha) m = b*np.log(b/sumsb*(alpha+1) ) return -1 if l+m<0 else np.sqrt((l+m)*2) def Checks(): if extension>1: print("Extension comparable to the size of the MAGIC camera cannot be simulated") return False if (numoff<=0) or (numoff>7): print("Number of OFF estimation regions must be in the range 1-7") return False if (extension >0.5) and (numoff>1): print("For large source extensions 1 OFF estimation region (numoff) should be used") return False return True ## MAIN code from mss.cpp if not Checks(): print("exiting") raise SystemExit # exit npoints_array, enbins, crabrate, bgdrate = Prepare() en = []; sed = []; dsed = []; sigmas = []; detected = []; for i, e1, e2 in zip(range(0,len(enbins)),enbins, enbins[1:]): intcrab, error = quad(Crab, e1, e2) intass, error = quad(Assumed, e1, e2) noff = bgdrate[i]*timeh*60 # number of off events noff *= (PSF*PSF+extension*extension)/(PSF*PSF) #larger integration cut due to extension dnoff = np.sqrt(noff/numoff) # error on the number of off events (computed from numoff regions) nexc = crabrate[i]*timeh*60*intass/intcrab dexc = np.sqrt(nexc + noff + dnoff*dnoff) # error on the number of excess events sigma = 0. # for tiny excesses (1.e-7 events) the function below can have # numerical problems, and either way sigma should be 0 then if nexc>0.01: sigma = SignificanceLiMa(nexc+noff, noff*numoff, 1./numoff) detect = False if sigma>=5.0 and nexc/noff>minSBR and nexc>minev: detect = True print('{0:.1f}-{1:.1f} GeV: exc. = {2:.1f} ev., SBR={3:.2f}%, sigma = {4:.1f}'\ .format(enbins[i], enbins[i+1], nexc, 100.*nexc/noff, sigma)," DETECTION" if detect else "") if nexc>minerror*dexc: tmpen = np.sqrt(enbins[i]*enbins[i+1]) en.append(tmpen) tmpsed = CalcSED(tmpen, Assumed(tmpen)) sed.append(tmpsed) dsed.append(tmpsed*dexc/nexc) sigmas.append(sigma) detected.append(detect) if len(sigmas)>0: sig_comb = sum(sigmas)/np.sqrt(len(sigmas)) print("Combined significance (using the {0:d} data points shown in the SED) = {1:.1f}"\ .format(len(sigmas), sig_comb)) if sig_comb<4: print("The source probably will not be detected by MAGIC in {0:.2f} hours".format(timeh)) elif sig_comb<6: print("The source probably might be detected by MAGIC in {0:.2f} hours".format(timeh)) else: print("The source probably will be detected by MAGIC in {0:.2f} hours".format(timeh)) # preparing reference SED graph fig, ax = plt.subplots(figsize=(7, 5)) #fig, ax = plt.subplots(figsize=(7, 5), facecolor='white') ax.set_xscale('log') ax.set_yscale('log') ax.set(xlim=(eplotmin, eplotmax), ylim=(yplotmin, yplotmax), xlabel='E [GeV]', ylabel='E$^2$ dN/dE [TeV cm$^{-2}$ s$^{-1}$]', title="$\Sigma\sigma_i / \sqrt{"+str(len(sigmas))+"} ="+str(round(sig_comb,1))+"\sigma $") x = np.logspace(np.log10(eplotmin), np.log10(eplotmax), 50) labeltext = "expected MAGIC SED ($T_{obs}$ = "+str(timeh)+" h)" plt.plot(x, 1.e-6*x*x*Assumed(x), 'lightgreen', label = srcname+" (Assumed)" ) if len(en)>0: ax.errorbar(en, sed, yerr = dsed, label = labeltext, color='0', fmt='o') plt.plot(x, 1.e-6*x*x*Crab(x), '0.3', label="Crab (Aleksic et al 2016)" ) ax.legend(loc="upper right") ax.grid(True,which="both",axis = 'x', ls="--",color='0.95') ax.grid(True,which="major",axis = 'y', ls="--",color='0.95') if drawsigma: for i in range(0,len(sigmas)): col = '0' if detected[i] else '0.75' ax.text(en[i], 2*sed[i], "{0:.1f}$\sigma$".format(sigmas[i]),\ rotation=90, size=10, ha='left', va='bottom', color=col ) fig.savefig('SED_'+srcname+'.pdf'); print ("Saving figure SED_"+srcname+".pdf") # fig.show() plt.close(fig)