#!/usr/bin/env python """ Calcrate is a utility to access instrument resonse functions and fluxes and to compute binned sensitivity estimates. The type output is determined by -m option: -m energy : print table vs E of Aeff, Flux, efficiencies, and rates -m mrf : print a table of the total signal, background and model rejection factor (mrf) for each declination, and cut. -m opt : optimize the cone and energy cut (best MRF if -l MRF or best MDP if -l MDP) and print table with results -m gopt : show plots of rates & mrf vs cut-values options can be combined; e.g. -m energy:mrf:opt to do all three (*) For the variables -d,-a and -e, you can indicate a range of values; e.g: -d 80:80:20 will do all declinations from -80 to 80 with steps of 20 degrees. examples: calcrate.py -F 1e-4*(E**-2) calcrate.py -F 1e-4*(E**-2) -m mrf -d -80:80:10 calcrate.py -F "1e-4*(E**-2)" -m gopt -d -40:-35:999 -a 1:1.5:999 -e 0:5:999 -w 0.6 -u 0 #RXJ LOI calcrate.py -m mrf -d -40:-35:999 -F "3.41e-10*pow(x[0]/1e3,-1.8)*exp(-sqrt(x[0]/3.7e3))" -g True -a 1:1.5:999 -e 0:5:999 -w 0.6 -u 0 calcrate.py -m mrf -d -40:-35:999 -F "16.8e-11*pow(x[0]/1e3,-1.72)*exp(-sqrt(x[0]/2.1e3))" -a 1:1.5:999 -e 0:5:999 -w 0.6 -u 0 -t numu #OPT calcrate.py -m opt -d -40:-35:999 -F "3.41e-10*pow(x[0]/1e3,-1.8)*exp(-sqrt(x[0]/3.7e3))" -g True -a 1:1.5:0.1 -e 0:5:1 -w 0.6 -u 0 #RXJ AGATA PAPER calcrate.py -m mrf -d -40:-35:999 -F "2.3e-10*pow(x[0]/1e3,-2.06)*exp(-x[0]/12.9e3)" -g True -a 1:1.5:999 -e 0:5:999 -w 0.6 -u 0 calcrate.py -m mrf -d -40:-35:999 -F "8.9e-11*pow(x[0]/1e3,-2.06)*exp(-x[0]/8.04e3)" -a 1:1.5:999 -e 0:5:999 -w 0.6 -u 0 -t numu #VELAX python -i calcrate.py -m mrf -d -45.6:80:999 -F "1.46e-10*pow(x[0]/1e3,-1.32)*exp(-x[0]/1.4e4)" -g True -a 1:9:999 -e 0:99:999 -w 0.8 -u 0 python -i calcrate.py -m mrf -d -45.6:80:999 -F "7.2e-11 *pow(x[0]/1e3,-1.36)*exp(-x[0]/7e3)" -a 1:100:999 -e 0:1:999 -w 0.8 -u 0 python -i calcrate.py -m mrf -d -45.6:80:999 -F "7.2e-11 *pow(x[0]/1e3,-1.36)*exp(-x[0]/7e3)" -a 1:100:999 -e 0:1:999 -w 0.8 -u 0 -t numu #HAWC1907 python -i calcrate.py -m mrf -d 6.32:80:999 -F "0.95e-12*pow(x[0]/1e4,-2.46-0.11*log10(x[0]/1e4))" -g True -a 1:9:999 -e 0:99:999 -w 0.67 -u 1 options: -f : RDF/ingredients ROOT files default=$AADIR/data/aanet_detresponse.cut2.v1.0.0.root -m : output-type (see above) default=energy -F : neutrino flux expression [1/(GeV m2 s sr)] default=(1e-4*(E**-2)) -g : flux expression refers to gamma-rays default=False -T : livetime (s) default=31557600 -YR : years of observation default=1:10:999 -t : Neutrino flavour interactions default=numuCC,anumuCC -c : observational channels default=track -p : point source mode (flux in /GeVm2s) default=True -d : declination(range) for point-mode (degrees) default=0:99:999 -a : search cone radius(range) for point-mode (degrees) default=1.0:99:999 -e : minimal log of reconstructed energy (range) default=0:99:999 -w : source width (degrees) default=0.0 -u : source shape (if width>0); 0=disk, 1=gaus default=0 -z : confidence level for limits default=0.90 -s : significance level for discovery default=3 -b : 1-b statistical power for discovery default=0.50 -y : output style (str,html,latex) default=str -h : help message and exit default=False -i : Python interative mode (prompt when done) default=False -FC : use Feldman Cousins intervals default=False -o : path+filename for plots default=interact -X : asimov: nexps,nbins,minloge,cone_trk,cone_casc default=10000,5,0,0.293,3.56 -l : MRF or MDP default:MRF """ #-f : RDF/ingredients ROOT files default=$AADIR/data/aanet_detresponse.cut2.v1.0.0.root,$AADIR/data/aanet_detresponse.cut2.v1.0.0.25.root import aa, ROOT, sys, itertools from ana.search.limits import * from ana import search from util import aautil from ROOT import Table, track, shower, nue, anue, anumu, numu from math import * from util import pyroot_util from array import array import numpy as np options = aa.Options( __doc__.split("options:")[1], sys.argv[1:] ) chans = [ ROOT.channel_from_name( c ) for c in options.c.split(",") ] flav_ints = [ ROOT.flavor_interaction_from_name( f ) for f in options.t.split(",") ] print (" channels " + ",".join( ROOT.channel_name( c ) for c in chans ) ) print (" interaction flavors " + ",".join( ROOT.flavor_interaction_name( f ) for f in flav_ints ) ) flavs=ROOT.get_flavor_from_flavor_interaction(flav_ints) print (" flavors " + ",".join( ROOT.flavor_name( f ) for f in flavs ) ) files=options.f.split(",") print(files) detress = [] index=0 for f in files: print (f) detress.append( ROOT.FullDetResponse( "detresponse"+str(index), f , flav_ints , chans, options.w, options.u )) index=index+1 expr = options.F.replace("E","x") if options.F.endswith(".txt") : # read a text file with flux. First column is Energy, second flux in units GeV-1 m-2 s-1 f = open( options.F ) loge_s, flux_s = [], [] for line in f : if line.startswith("#") : continue energy, flux = map( float, line.split() ) loge_s.append(log10(energy)) flux_s.append(flux) graph = ROOT.TGraph( len(loge_s), array('d',loge_s), array('d', flux_s ) ) ROOT.gInterpreter.Declare(""" struct GraphFlux: public Flux { TGraph h; GraphFlux( TGraph& g ) : h(g) {} virtual double dNdE( int pdg_type, double loge ) { return h.Eval( loge ); } virtual void Print() { cout << "graphFlux:" << h.GetName() << endl; } }; """) F = ROOT.GraphFlux( graph ) else : # it is an expression expr = options.F.replace("E","x") if not options.g : F = ROOT.NuExprFlux( 1, expr ) else : print ("gamma-ray flux : ", expr ) F = ROOT.NuFromGammaFlux( expr ) point = options.p livetime = options.T significance=(1-ROOT.TMath.Erf(options.s/sqrt(2))) # 2 side convention #year = options.YR outfile = options.o limit =options.l if outfile !="interact": ROOT.gROOT.SetBatch() try : Table.__str__ = { "html" : Table.html, "str" : Table.str, "latex": Table.latex }[ options.y ] except : print ("invalid option", options.y ,"for table output style.") raise def torange( opt ) : "string like 1:10:3 to range(1,10,3)" if type(opt) in (float,int) : return [opt] return aautil.frange(*map(float, opt.split(":"))) decl_range = list( map ( lambda a:a*pi/180 , torange( options.d ) ) ) cone_range = list( map ( lambda a:a*pi/180 , torange( options.a ) ) ) year_range = torange( options.YR ) #print(year_range,len(year_range),type(year_range)) minloge_range = torange( options.e ) ranges = list( itertools.product( year_range,decl_range, cone_range, minloge_range ) ) funcs = {} def print_energy_table(): "print extensive table of Aeff,Flux,Rate vs Energy" for detres in detress: for chan in chans : for flav in flavs : print (flav, flavs) for year,decl,cone,minloge in ranges : table = detres.get(flav,chan).compute_rates_table( F, livetime, point, decl , cone, minloge ) print(table) funcs['energy'] = print_energy_table """ def sensi_plot_decl(nsig_vec,nbkg_vec,decl_list,year,plot): decls=array('d',decl_list) declinations=array('d') for ddd in decls: declinations.append(ddd*180./pi) lim_vec, mrf_vec = array( 'd' ), array( 'd' ) n3s_vec, mdp_vec = array( 'd' ), array( 'd' ) if outfile != "interact" : out_txt = open(options.o+".dat", "w") for i in range(len(declinations)): nbkg=nbkg_vec[i]*year nsig=nsig_vec[i]*year decl=declinations[i] #print("debug ",year,nsig_vec[i],nsig,nbkg_vec[i],nbkg,decl) lim=mean_limit( nbkg , options.z,options.FC ) lim_vec.append(lim) mrf = lim/nsig mrf_vec.append(mrf) n3s=get_mu_lds(nbkg,significance,options.b) n3s_vec.append(n3s) mdp = n3s/nsig mdp_vec.append(mdp) print(declinations[i],",",nsig,",",nbkg,",",lim,",",mrf,",",n3s,",",mdp) if outfile != "interact" : out_txt.write(str(declinations[i])+","+str(nsig)+","+str(nbkg)+","+str(lim)+","+str(mrf)+","+str(n3s)+","+str(mdp)+"\n") if(plot): gr_n90=ROOT.TGraph(len(declinations),declinations,lim_vec) gr_n90.GetXaxis().SetTitle("declination [degree]") gr_n90.GetXaxis().CenterTitle() gr_n90.GetYaxis().SetTitle("n90%CL") gr_n90.GetYaxis().CenterTitle() gr_sensi=ROOT.TGraph(len(declinations),declinations,mrf_vec) gr_sensi.GetXaxis().SetTitle("declination [degree]") gr_sensi.GetXaxis().CenterTitle() gr_sensi.GetYaxis().SetTitle("MRF") gr_sensi.GetYaxis().CenterTitle() gr_n3sigma=ROOT.TGraph(len(declinations),declinations,n3s_vec) gr_n3sigma.GetXaxis().SetTitle("declination [degree]") gr_n3sigma.GetXaxis().CenterTitle() gr_n3sigma.GetYaxis().SetTitle("n3#sigma") gr_n3sigma.GetYaxis().CenterTitle() gr_discovery=ROOT.TGraph(len(declinations),declinations,mdp_vec) gr_discovery.GetXaxis().SetTitle("declination [degree]") gr_discovery.GetXaxis().CenterTitle() gr_discovery.GetYaxis().SetTitle("MDP") gr_discovery.GetYaxis().CenterTitle() c1 = ROOT.TCanvas('c1','c1',800,800) ROOT.SetOwnership(c1, False ) c1.Divide(2,2) c1.cd(1) gr_n90.DrawClone("APL") c1.cd(2) gr_sensi.DrawClone("APL") c1.cd(3) gr_n3sigma.DrawClone("APL") c1.cd(4) gr_discovery.DrawClone("APL") if outfile != "interact" : print("save the plots in ",outfile+"_sensi_decl.png") c1.SaveAs(outfile+"_sensi_decl.png") if outfile != "interact": out_txt.close() return declinations,lim_vec, mrf_vec """ """ def sensi_plot_year(nsig,nbkg,year_list,plot): year_vec=array('d',year_list) lim_vec, mrf_vec = array( 'd' ), array( 'd' ) n3s_vec, mdp_vec = array( 'd' ), array( 'd' ) if outfile != "interact" : out_txt = open(options.o+".dat", "w") for y in year_vec: nsigy=nsig*y nbkgy=nbkg*y lim=mean_limit( nbkgy , options.z,options.FC ) lim_vec.append(lim) mrf = lim/nsigy mrf_vec.append(mrf) n3s=get_mu_lds(nbkgy,significance,options.b) n3s_vec.append(n3s) mdp = n3s/nsigy mdp_vec.append(mdp) print(y,",",nsigy,",",nbkgy,",",lim,",",mrf,",",n3s,",",mdp) if outfile != "interact" : out_txt.write(str(y)+","+str(nsigy)+","+str(nbkgy)+","+str(lim)+","+str(mrf)+","+str(n3s)+","+str(mdp)+"\n") if(plot): gr_n90=ROOT.TGraph(len(year_vec),year_vec,lim_vec) gr_n90.GetXaxis().SetTitle("nyears") gr_n90.GetXaxis().CenterTitle() gr_n90.GetYaxis().SetTitle("n90%CL") gr_n90.GetYaxis().CenterTitle() gr_sensi=ROOT.TGraph(len(year_vec),year_vec,mrf_vec) gr_sensi.GetXaxis().SetTitle("nyears") gr_sensi.GetXaxis().CenterTitle() gr_sensi.GetYaxis().SetTitle("MRF") gr_sensi.GetYaxis().CenterTitle() gr_n3sigma=ROOT.TGraph(len(year_vec),year_vec,n3s_vec) gr_n3sigma.GetXaxis().SetTitle("nyears") gr_n3sigma.GetXaxis().CenterTitle() gr_n3sigma.GetYaxis().SetTitle("n3#sigma") gr_n3sigma.GetYaxis().CenterTitle() gr_discovery=ROOT.TGraph(len(year_vec),year_vec,mdp_vec) gr_discovery.GetXaxis().SetTitle("nyears") gr_discovery.GetXaxis().CenterTitle() gr_discovery.GetYaxis().SetTitle("MDP") gr_discovery.GetYaxis().CenterTitle() c1 = ROOT.TCanvas('c1','c1',800,800) ROOT.SetOwnership(c1, False ) c1.Divide(2,2) c1.cd(1) gr_n90.DrawClone("APL") c1.cd(2) gr_sensi.DrawClone("APL") c1.cd(3) gr_n3sigma.DrawClone("APL") c1.cd(4) gr_discovery.DrawClone("APL") if outfile != "interact" : print("save the plots in ",outfile+"_sensi_year.png") c1.SaveAs(outfile+"_sensi_year.png") if outfile != "interact": out_txt.close() return year_vec,lim_vec, mrf_vec """ def get_sig_bg( chan, F, livetime, point, decl , cone, minloge, year,FC=False ) : sig, bg = 0.0, 0.0 for detres in detress: for flav in flavs : sig_,bg_,_ = detres.get(flav,chan).compute_rates( F, livetime, point, decl , cone, minloge ) sig+=sig_*year bg+=bg_*year lim = mean_limit( bg , options.z,FC ) mrf = lim/sig mu_lds=get_mu_lds(bg,significance,options.b) mdp=mu_lds/sig return sig,bg,lim,mrf,mu_lds,mdp """ def plot_year(): "sensitivity and discovery plots vs year" if len(cone_range)>1 : print("Only 1 alpha value supported in plot mode" ) exit() if len(minloge_range)>1 : print("Only 1 loge_min value supported in plot mode" ) exit() if len(decl_range)>1 : print("Only 1 decl value supported in plot(year) mode" ) exit() if len(year_range)<2 : print("years array not provided, no plot available" ) exit() for chan in chans : table = ROOT.Table("decl","cone","minloge","total-sig","total-bg","meanlimit","MRF","n limit discovery","MDP") for decl,cone,minloge in ranges : # todo: decl does not make sense for diffuse sig, bg, lim, mrf,mu_lds,mdp = get_sig_bg( chan, F, livetime, point, decl, cone, minloge, options.FC ) table.add(decl).add(cone * 180/pi ).add(minloge).add(sig).add(bg).add(lim).add( lim/sig ).add(mu_lds).add( mu_lds/sig ) sensi_plot_year(sig,bg,year_range,True) flavnames = [ ROOT.flavor_name( f ) for f in flavs ] table.title = " results for " + expr + " " + "+".join( flavnames ) + " " + ROOT.channel_name(chan) print( table ) funcs['plot_year'] = plot_year """ """ def plot_decl(): "sensitivity and discovery plots vs decl" if len(cone_range)>1 : print("Only 1 alpha value supported in plot mode" ) exit() if len(minloge_range)>1 : print("Only 1 loge_min value supported in plot mode" ) exit() if len(year_range)>1 : print("Only 1 year value supported in plot(decl) mode" ) exit() if len(decl_range)<2 : print("decl array not provided, no plot available" ) exit() if point==False : print("decl does not make sense for diffuse" ) exit() nsig_vec = [] nbkg_vec = [] for chan in chans : table = ROOT.Table("decl","cone","minloge","total-sig","total-bg","meanlimit","MRF","n limit discovery","MDP") for year,decl,cone,minloge in ranges : # todo: decl does not make sense for diffuse sig, bg, lim, mrf,mu_lds,mdp = get_sig_bg( chan, F, livetime, point, decl, cone, minloge, options.FC ) table.add(decl).add(cone * 180/pi ).add(minloge).add(sig).add(bg).add(lim).add( lim/sig ).add(mu_lds).add( mu_lds/sig ) nsig_vec.append(sig) nbkg_vec.append(bg) sensi_plot_decl(nsig_vec,nbkg_vec,decl_range,year,True) flavnames = [ ROOT.flavor_name( f ) for f in flavs ] table.title = " results for " + expr + " " + "+".join( flavnames ) + " " + ROOT.channel_name(chan) print( table ) funcs['plot_decl'] = plot_decl """ def print_mrf_table(): "show total signal and background for each channel (added over flavours)" if outfile != "interact" : out_txt = open(options.o+".dat", "w") for chan in chans : table = ROOT.Table("year","decl","cone","minloge","total-sig","total-bg","meanlimit","MRF","n limit discovery","MDP") for year,decl,cone,minloge in ranges : # todo: decl does not make sense for diffuse sig, bg, lim, mrf,mu_lds,mdp = get_sig_bg( chan, F, livetime, point, decl, cone, minloge, year,options.FC ) table.add(year).add(decl).add(cone * 180/pi ).add(minloge).add(sig).add(bg).add(lim).add( lim/sig ).add(mu_lds).add( mu_lds/sig ) if outfile != "interact" : out_txt.write(str(year)+" "+str(decl)+" "+str(sig)+" "+str(bg)+" "+str(lim)+" "+str(mrf)+" "+str(mu_lds)+" "+str(mdp)+" "+str(cone* 180/pi)+" "+ str(minloge)+"\n") flavnames = [ ROOT.flavor_name( f ) for f in flavs ] #table.title = " results for " + F.expression() + " " + "+".join( flavnames ) + " " + ROOT.channel_name(chan) table.title = " results for " + expr + " " + "+".join( flavnames ) + " " + ROOT.channel_name(chan) print( table ) if outfile != "interact" : out_txt.close() #plot_year(table) #plot_decl(table) funcs['mrf'] = print_mrf_table def print_optimized_mrf_table(): "optimize the cone and logE cut for best MRF and print results for each declination" for chan in chans : table = ROOT.Table("year","decl","cone(opt)","minlogE(opt)","total-sig","total_bg","mean limit","MRF","n limit discovery","MDP") if outfile != "interact" : out_txt = open(options.o+"_optimize_"+limit+".dat", "w") def score(cone,minloge) : sig, bg, lim, mrf,mu_lds,mdp = get_sig_bg( chan, F, livetime, point, decl, cone, minloge,year,options.FC ) if limit=="MRF": return mrf else: return mdp for year in year_range: for decl in decl_range: cone, minloge = pyroot_util.minimize( score, [ 0.1 *pi/180, 0.1 ], [[0.1 *pi/180,2.*pi/180],[0.,5.]],tollerance = 1e-8 ) sig,bg,lim,mrf,mu_lds,mdp = get_sig_bg( chan, F, livetime, point, decl, cone, minloge,year,options.FC ) table.add( year ).add( decl ).add( cone * 180/pi ).add( minloge ).add( sig ).add( bg ).add(lim).add( mrf ).add(mu_lds).add( mdp ) if outfile != "interact" : out_txt.write(str(year)+" "+str(decl)+" "+str(sig)+" "+str(bg)+" "+str(lim)+" "+str(mrf)+" "+str(mu_lds)+" "+str(mdp)+" "+str(cone* 180/pi)+" "+ str(minloge)+"\n") print (table) if outfile != "interact" : out_txt.close() funcs['opt'] = print_optimized_mrf_table def show_optimization_plots(): "show plots vs cone and energy cut" chan = chans[0] decl = decl_range[0] year = year_range[0] binning = (30, 0.2, 2.0, 100, 0.0, 2.0) hsig = ROOT.TH2D("hsig", "sig", *binning ) hbg = ROOT.TH2D("hbg", "bg", *binning ) hlim = ROOT.TH2D("hlim", "lim", *binning ) hmrf = ROOT.TH2D("hmrf", "mrf", *binning ) hmu_lds = ROOT.TH2D("hmu_lds", "lim", *binning ) hmdp = ROOT.TH2D("hmdp", "mrf", *binning ) for bx in range( 1, hsig.GetNbinsX() + 1 ) : for by in range( 1, hsig.GetNbinsY() + 1 ) : cone, minloge = hsig.GetXaxis().GetBinCenter( bx ), hsig.GetYaxis().GetBinCenter( by ) sig,bg, lim, mrf,mu_lds, mdp = get_sig_bg( chan, F, livetime, point, decl, cone / 180 * pi , minloge, year, options.FC ) hsig.SetBinContent( bx,by, sig ) hbg.SetBinContent( bx,by, bg ) hlim.SetBinContent( bx,by, lim ) hmrf.SetBinContent( bx,by, mrf ) hmu_lds.SetBinContent( bx,by,mu_lds ) hmdp.SetBinContent( bx,by, mdp ) c = ROOT.TCanvas('c','c',800,800) ROOT.SetOwnership(c, False ) c.Divide(3,2) c.cd(1) hsig.DrawCopy("colz") c.cd(2) hbg.DrawCopy("colz") c.cd(3) hlim.DrawCopy("colz") c.cd(4) hmrf.DrawCopy("colz") c.cd(5) hmu_lds.DrawCopy("colz") c.cd(6) hmdp.DrawCopy("colz") print("debug outfile: ", outfile) if outfile != "interact" : print("save the plots in ",outfile+"_opti.png") c.SaveAs(outfile+"_opti.png") funcs['gopt'] = show_optimization_plots # functions for asimov estimation def get_sig_bg_range( chan, F, livetime, point, decl, year, cone_min, cone_max, minloge, maxloge ) : "get signal and background for a certain cone bin and energy bin" sig, bg = 0.0, 0.0 for detres in detress: for flav in flavs : sig_,bg_,_ = detres.get(flav,chan).compute_rates_angle_energy_range( F, livetime, point, decl, cone_min, cone_max, minloge, maxloge ) sig+=sig_*year bg+=bg_*year return sig,bg def cone_range_equal_size( angle_low, angle_high, angle_bins ): "return radius edges of equal cone size bins" result = {} cone_area_bin = 2*pi*(1-cos(angle_high))/angle_bins # https://en.wikipedia.org/wiki/Solid_angle for i in range(1, angle_bins + 1): cone_area_low = cone_area_bin*(i-1) cone_area_high = cone_area_bin*(i) result[i] = [ acos(1 - cone_area_low/(2*pi)), acos(1 - cone_area_high/(2*pi)) ] return result def generate_energy_bins( energy_low, energy_high, energy_bins ): "return edges of equal energy bins" result = {} bin_size = (energy_high-energy_low)/energy_bins for i in range(1, energy_bins + 1): result[i] = [ energy_low + (i-1)*bin_size, energy_low + (i)*bin_size ] return result def asimov(): "asimov sensitivy estimate" "call using options.X = nexps,nbins,minloge,cone_trk,cone_casc" "nexps used to cross check results" parse = options.X.split(",") nexps = int(parse[0]) nbins = int(parse[1]) minloge = float(parse[2]) cone = { ROOT.channel_from_name("track"):float(parse[3])*pi/180, ROOT.channel_from_name("shower"):float(parse[4])*pi/180 } maxloge = 8 for year in year_range: for decl in decl_range: total_sig = {} total_bg = {} s = ROOT.vector('double')() b = ROOT.vector('double')() for chan in chans : energy_bins = generate_energy_bins(minloge, maxloge, nbins) # store min and max of bin angle_bins = cone_range_equal_size( 0.0, cone[chan], nbins ) total_sig[chan], total_bg[chan],_,_,_,_ = get_sig_bg( chan, F, livetime, point, decl, cone[chan], minloge, year, options.FC ) total_sig_check = 0 total_bg_check = 0 for energy_bin in energy_bins: for angle_bin in angle_bins: print(energy_bin, angle_bin) sig_, bg_ = get_sig_bg_range( chan, F, livetime, point, decl, year, angle_bins[angle_bin][0], angle_bins[angle_bin][1], energy_bins[energy_bin][0], energy_bins[energy_bin][1] ) s.push_back(sig_) b.push_back(bg_) total_sig_check+=sig_ total_bg_check+=bg_ print("channel", ROOT.channel_name(chan) ) print("total sig in all bins", total_sig[chan], "cross check with get_sig_bg", total_sig_check ) print("total bg in all bins", total_bg[chan], "cross check with get_sig_bg", total_bg_check ) p = ROOT.AsimovHypothesis(s,b) # limit setting mu_limit = p.find_mu( ROOT.Z_from_pvalue( 1 - options.z ), False ) # p value = 1 - cl print(20*"-", "limit setting", 20*"-") print("asimov estimation: mu_limit", mu_limit) if nexps > 0: p_value = p.p_value( mu_limit, mu_limit, 0.0, nexps) # cross check result print("corresponding p from pseudo-exp:", p_value) print(5*"-", "tune mu limit to obtain right p: {}".format( 1- options.z ), 5*"-") while p_value > 1 - options.z: mu_limit+=0.01 p_value = p.p_value( mu_limit, mu_limit, 0.0, nexps) print("mu_limit", mu_limit, "p_value", p_value) # discovery mu_discovery = p.find_mu( options.s, True ) print(20*"-", "discovery", 20*"-") print("asimov estimation: mu_discovery", mu_discovery) if nexps > 0: p_value = p.p_value( 0.0, 0.0, mu_discovery, nexps) print("corresponding p from pseudo-exp:", p_value) print(5*"-", "tune mu to obtain right p: {}".format( ROOT.pvalue_from_Z( options.s ) ), 5*"-") while p_value > ROOT.pvalue_from_Z( options.s ): mu_discovery+=0.01 p_value = p.p_value( 0.0, 0.0, mu_discovery, nexps) print("mu_discovery", mu_discovery, "p_value", p_value) funcs['asimov'] = asimov modes = options.m.split(":") for mode in modes : funcs[mode]()