# # 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: # - ROOT (https://root.cern.ch), the macro was tested on 5.34/34 version # # 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 "basic settings" below and run in terminal: # root mss.cpp+ # # for questions ask Julian Sitarek (jsitarek [at] uni.lodz.pl) # # version history: # 1.2 (2020/07/24): a number of additional features: # - EBL attenuation with Dominguez+11 model # - SUMT performance (only for low zenith) # - pulsar mode # 1.1 (2017/10/16): fixed a bug which made the errorbars too smallget_ipython().run_line_magic('matplotlib', 'inline') import matplotlib.pyplot as plt import numpy as np from scipy.integrate import quad from scipy import interpolate # Source settings # [cm^-2 s^-1 TeV^-1] assumed (intrinsic) source spectrum (dN/dS dt dE) as a function of energy in GeV # you can create your own, some common examples below def Assumed(x): return 2.7e-9*pow(x/100., -3.18) # simple power-law #return 1.0e-11*pow(x/1000., -1.9)*np.exp(-x/100) # power-law with an exponential cut-off #return 1.0e-11*pow(x/1000., -1.9)*np.exp(-pow(x/300, 0.5)) # power-law with a non-exponential cut-off #return 0.007*3.39e-11*pow(x/1000., -2.51-0.21*np.log10(x/1000.)) # Crab Nebula [Aleksic et al. 2016] #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# More examples of assumed spectrum: # # const Double_t redshift=0.94; # const TString assumedstr="2.7e-9 * pow(x/100, -3.18)"; // simple power-law # const TString srcname="PKS1441"; // name of the source # # const TString assumedstr="2.28e-9 * pow(x/32.1, -5.62)"; // simple power-law # const TString srcname="Geminga"; // name of the source # # const TString assumedstr="4.e-9 * pow(x/100., -2.6)"; // simple power-law # const TString srcname="PG 1553"; // name of the source # # const TString assumedstr="1.4e-12 * pow(x/1310., -2.74)"; // simple power-law # const TString srcname="J1835-069"; // name of the source # # const TString assumedstr="7.9e-11 * pow(x/175, -3.26)"; // simple power-law # const Double_t redshift=0.360; // redshift of the source (for the EBL absorption), if -1 then no absorption # const TString srcname="PKS1510low"; // name of the source # basic settings timeh=20 # [h], time of observations extension=0. # [deg] extension radius of the source # source settings redshift = -1 # redshift of the source (for the EBL absorption), if -1 then no absorption srcname ="Test" # 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/ isSUMT = False # 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) offsetdegrad = 1. # degradation factor if observations are taken at higher offset then 0.4 deg from camera center 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# pulsar mode settings pulsarmode = False # if true the background is reduced to on phase (see below) and SBR cut is ignored pulsarOnRange = 0.092 # range of ON phases used for pulsar mode pulsarOffRange = 0.25 # range of OFF phases used for pulsar mode# global variables (DO NOT MODIFY) npoints = 13 pathebl = "dominguez_ebl_tau.txt" # path with EBL model of Dominguez+11 #Crab Nebula [Aleksic et al. 2016] def Crab(x): return 3.39e-11*pow(x/1000., -2.51-0.21*np.log10(x/1000.)) #const TString crabstr="3.39e-11 * pow(x/1000, -2.51-0.21*log10(x/1000))"; //Crab Nebula [Aleksic et al. 2016] # Calculate SED value from energy and flux, for display def CalcSED(x,flux): return 1.e-6*x*x*flux def LoadEBL(): f1 = open(pathebl, 'r') line = f1.readline() firstlineEBL = list(map(float, line.split())) nz = firstlineEBL[0] zz = firstlineEBL[1:] #print(nz, zz) f1.close() eblbody = np.loadtxt(pathebl,delimiter=' ', skiprows=1) energies = eblbody[:,0]*1.e3 #in GeV #en*=1.e3; // GeV taus = eblbody[:,1:] if len(taus>0): return zz, energies, taus def FluxObs(xx, z, fluxint): en = xx EBLz, EBLen, EBLtaus = LoadEBL() ftau = interpolate.interp2d(EBLz, EBLen, EBLtaus, kind='cubic') tau = ftau(z,en) atten = np.exp(-tau) return fluxint*atten def Prepare(): npoints_array = np.arange(0, npoints+1) enbins = 100.*pow(10., (npoints_array-2)*0.2) # [GeV] if isSUMT: if ismidzd: crabrate = np.array([]) bgdrate = np.array([]) # to be implemented else: crabrate = np.array([1.39684, 3.12657, 3.09145, 2.40268, 1.32915, 0.86180, 0.51666, 0.31533, 0.16207, 0.09279, 0.04624, 0.02345, 0.00874]) # [min^-1] # 1st point ~50GeV, last point ~12TeV bgdrate = np.array([ 3.33321, 3.24046, 1.32361, 0.406588, 0.091944, 0.032226, 0.007277, 0.003123, 0.001487, 0.001464, 0.001231, 0.001152, 0.000957]) # [min^-1] else: if ismidzd: # Aleksic et al 2016 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], 1st point below threshold !!! 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], 1st point below threshold !!! else: # 0-30 deg values, Aleksic et al 2016 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 if (isSUMT and ismidzd): print("Sorry not implemented yet :-(") return False if (offsetdegrad>1.00001): print("No cheating! the performance degradation ({0}) should not be larger then 1".format(offsetdegrad)) return False if pulsarmode: if (pulsarOnRange<=0 or pulsarOnRange>=1): print("Pulsar mode ON phase range is {0}, and it should be in range (0,1)".format(pulsarOnRange)) return False if (pulsarOffRange<=0 or pulsarOffRange>=1): print("Pulsar mode OFF phase range is {0}, and it should be in range (0,1)".format(pulsarOffRange)) return False if (redshift > 0): print("Do you really want to observe a pulsar at redshift of {0} ??".format(redshift)) return True ## MAIN code mss.cpp if not Checks(): print("exiting") #return 0 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) if redshift>0: intass, error = quad(lambda x: FluxObs(x,redshift,Assumed(x)), e1, e2) else: 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 noffon=0 if pulsarmode: noffon = noff*pulsarOnRange # number of bgd events in ON phase noff *= pulsarOffRange # number of bgd events in OFF phase dnoff = np.sqrt(noff)*pulsarOnRange/pulsarOffRange #ignoring numoff for pulsars and scaling for the phase difference dexc = np.sqrt(nexc + noffon + dnoff*dnoff) # error on the number of excess events # for tiny excesses (1.e-7 events) the function below can have numerical problems, and either way sigma should be 0 then sigma = 0 # added sept2021 if nexc>0.01: if (pulsarmode): sigma = SignificanceLiMa(nexc+noffon, noff, pulsarOnRange/pulsarOffRange) noff = noffon # needed later for SBR else: sigma = SignificanceLiMa(nexc+noff, noff*numoff, 1./numoff) detect = False if pulsarmode: if ((sigma>=5.0) and (nexc>minev)): detect = True else: if sigma>=5.0 and nexc/noff>minSBR and nexc>minev: detect = True print('{0:.1f}-{1:.1f} GeV: exc. = {2:.1f}+-{3:.1f} ev., SBR={4:.2f}%, sigma = {5:.1f}' .format(enbins[i], enbins[i+1], nexc, dexc, 100.*nexc/noff, sigma)," DETECTION" if detect else "") print(nexc, minerror, dexc) if nexc>minerror*dexc: tmpen = np.sqrt(enbins[i]*enbins[i+1]) en.append(tmpen) if redshift>0: tmpsed = CalcSED(tmpen, FluxObs(tmpen,redshift,Assumed(tmpen))) else: tmpsed = CalcSED(tmpen, Assumed(tmpen)) sed.append(float(tmpsed)) dsed.append(float(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 ",timeh, " hours") elif sig_comb<6: print("The source probably might be detected by MAGIC in ",timeh, " hours") else: print("The source probably will be detected by MAGIC in ",timeh, " hours") # preparing reference SED graph fig, ax = plt.subplots(figsize=(8, 6)) 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))+"$") 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*Crab(x), '0.3', label="Crab (Aleksic et al 2016)" ) if redshift>0: plt.plot(x, 1.e-6*x*x*[float(FluxObs(ee,redshift,Assumed(ee))) for ee in x], 'g', label = srcname+" (Assumed, z={0:.2f})".format(redshift) ) else: plt.plot(x, 1.e-6*x*x*Assumed(x), 'g', label = srcname+" (Assumed)" ) if len(en)>0: ax.errorbar(en, sed, yerr = dsed, label = labeltext, color='0', fmt='o') 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'); plt.show() plt.close(fig)