#!/usr/bin/env python #Plotting of Likelihood spaces Jordan """ usage : aashowerfit.py -f input_file Aashowerfit is a reconstruction program for cascade events. The options can be set via the command line (see below), or via the following environment variables AA_INFILE, AA_OUTFILE, AA_DETFILE, AA_LIBDIR, AA_NEVE, and AA_VERB. Environment variables are intended for use in batch jobs and take precident over command line arugments. The reconstructed tracks are appended to the Evt.trks container in the aanet event. This means that previous reconstructions are not overwritten by default (see option -c). Aashowerfit will write the following tracks (in order): o 1 starting point of the M-estimator fit , rec_type = 101, rec_stages = {1} o 1 (M-est) pos+time fit , rec-type = 101, rec_stages = {1,2} o 11 sub-optimal likelihood fits , rec-type = 101, rec_stages = {1,2,3} o 1 optimal likelihood fit (best fit) , rec-type = 101, rec_stages = {1,2,3,4} (In fact the 12 likelihood fits are simply sorted by likelihood, but the best (i.e. last) one receives a diffent value of rec_stages for easy identificaion in analysis.) The likelihood fit tracks (rec_stage>=3) have: fitinf[0] = uncorrected energy, fitinf[1] = number of hits used in the fit options: -f -input_file : Input file in any format aanet can read. default=/sps/km3net/users/jseneca/aanet/data/km3_sim/km3_v4_nueCC_1.JTE.aareco.root -i -id : Event id number default=0 -l -ll_files: default=./rec_files/ -n : process this many events default=100000000 -e : only process events w/ MC-Enu greater than this default=0 -v : verbose (0,1 or 2) default=0 -t : report timing informtion every N events default=100 -s : Show plots default=False """ print "aashowerfit" import ROOT,os,sys,collections, array from math import * from copy import copy import matplotlib.pyplot as pl import numpy as np import pickle try : import aa except ImportError: print "Failed to import aa; Make sure AADIR is in your PYHTONPATH!" sys.exit() import rec options = aa.Options( __doc__, sys.argv[1:] ) #options.read_env() print options chosen_id = options.i #Energy increases in the file from ROOT import vector, Hit, Det, Evt, Trk, Vec, EventFile, Timer def todeg(rad): return rad*180./pi def torad(deg): return deg*pi/180. #----------------------------------------------------- # Fancy-smancy setting of members #----------------------------------------------------- Trk.set = aa.__set #----------------------------------------- # input & output file #----------------------------------------- basefilename = os.path.basename(options.f) f = ROOT.EventFile( options.f ) print("Filename: ", options.f) #=================================================================================== # end of setup #=================================================================================== for evt in f: if evt.id < chosen_id: continue print("Found event: ", evt.id) fits_3d = [] fits_5d = [] fit_eval = [] LLs = [] try: print("Trying to open file: %s" %(options.l + basefilename + "%i_LLs" %(evt.id) ) ) load_LLs = open(options.l + basefilename + "%i_LLs" %(evt.id), 'r') LLs = pickle.load(load_LLs) except IOError as error: print("ERROR IN LOADING FILES: ", error) with Timer("1 plotting"): Lspace_3d = False zzzscale = 0.005 zzscale = 1 zscale = 10 pscale = 50 neut_D = evt.mc_trks[0].dir.normalize() neut_index = [neut_D.phi()*pscale/pi + pscale, neut_D.theta()*pscale/pi] neut_E = evt.mc_trks[0].E lep_D = evt.mc_trks[1].dir.normalize() reco3d_D = fits_3d[-1].dir.normalize() reco3d_E = fits_3d[-1].E reco5d_D = fits_5d[-1].dir.normalize() reco5d_E = fits_5d[-1].E startDir = Vec(0.0, 0.0, 0.0) startDir.set_angles(evt.mc_trks[0].dir.theta() + 0.1, evt.mc_trks[0].dir.phi() + 0.1) minLL_D = Vec(0.0, 0.0, 0.0) min_eval_trks = fit_eval[0] min_evals = fit_eval[1] #use reco Energy max_LL = 0.0 min_LL = 999999.0 Zmin_LL = 999999.0 ZZmin_LL = 999999.0 true_LL = LLs[0] reco5d_LL = LLs[1] LL_5d = np.zeros((pscale, 2*pscale)) ThetaArray = np.linspace(0, pi, pscale) PhiArray = np.linspace(-pi, pi, 2*pscale) try: pickled_file = open(options.l + basefilename + "%iLL_TESTF" %(evt.id), 'r') LL_5d = pickle.load(pickled_file) except IOError: print("Could not open LL") ZLL_5d = np.zeros((pscale, pscale)) ZThetaArray = np.linspace(neut_D.theta() - torad(zscale), neut_D.theta() + torad(zscale), pscale) ZPhiArray = np.linspace(neut_D.phi() - torad(zscale), neut_D.phi() + torad(zscale), pscale) try: pickled_file = open(options.l + basefilename + "%iZLL_TESTF" %(evt.id), 'r') ZLL_5d = pickle.load(pickled_file) except IOError: print("Could not open ZLL") Zmin_LL = ZLL_5d.min() minLL_indices = np.unravel_index(ZLL_5d.argmin(), ZLL_5d.shape) minLL_D.set_angles(ZThetaArray[minLL_indices[0]], ZPhiArray[minLL_indices[1]]) #extra zoomed version ZZLL_5d = np.zeros((pscale, pscale)) ZZThetaArray = np.linspace(reco5d_D.theta() - torad(zzscale), reco5d_D.theta() + torad(zzscale), pscale) ZZPhiArray = np.linspace(reco5d_D.phi() - torad(zzscale), reco5d_D.phi() + torad(zzscale), pscale) try: pickled_file = open(options.l + basefilename + "%iZZLL_TESTF" %(evt.id), 'r') ZZLL_5d = pickle.load(pickled_file) except IOError: print("Could not open reco ZZLL") #extra super zoomed version ZZZLL_5d = np.zeros((pscale, pscale)) ZZZThetaArray = np.linspace(reco5d_D.theta() - torad(zzzscale), reco5d_D.theta() + torad(zzzscale), pscale) ZZZPhiArray = np.linspace(reco5d_D.phi() - torad(zzzscale), reco5d_D.phi() + torad(zzzscale), pscale) try: pickled_file = open(options.l + basefilename + "%iZZZLL_TESTF" %(evt.id), 'r') ZZZLL_5d = pickle.load(pickled_file) except IOError: print("Could not open reco ZZZLL") #1d version one_LL_5d = np.zeros(pscale * pscale) one_Theta = pi/2. one_PhiArray = np.linspace(-pi, pi, pscale*pscale) try: pickled_file = open(options.l + basefilename + "%ione_LL_TESTF" %(evt.id), 'r') one_LL_5d = pickle.load(pickled_file) except IOError: print("Could not open reco one_LL") #plot starting direction vs fitted print("--------------fitted direction: ", reco3d_D) print("Number of recos: %d" %len(fits_5d)) print("True LL: %f, ZminLL: %f, reco_LL: %f" %(true_LL, Zmin_LL, reco5d_LL)) print("Neut dir theta: %f, phi: %f" %(neut_D.theta(), neut_D.phi())) print("Lep dir theta: %f, phi: %f" %(lep_D.theta(), lep_D.phi())) print("start reco dir theta: %f, phi: %f" %(startDir.theta(), startDir.phi())) print("3d reco dir theta: %f, phi: %f" %(reco3d_D.theta(), reco3d_D.phi())) print("5d minLL dir theta: %f, phi: %f" %(minLL_D.theta(), minLL_D.phi())) print("5d reco dir theta: %f, phi: %f" %(reco5d_D.theta(), reco5d_D.phi())) print("Neut E: %f" %(neut_E)) print("3d reco E: %f" %(reco3d_E)) print("5d reco E: %f" %(reco5d_E)) res3dAng = (180/pi)*acos(neut_D.dot( reco3d_D ) ) res5dAng = (180/pi)*acos(neut_D.dot( reco5d_D ) ) res3dE = (reco3d_E - neut_E) / neut_E res5dE = (reco5d_E - neut_E ) / neut_E print("3dAng reco diff: %f" %res3dAng ) print("5dAng reco diff: %f" %res5dAng ) print("3dE reco diff: %f" %res3dE ) print("5dE reco diff: %f" %res5dE ) pl.plot(todeg( neut_D.phi() ), todeg(neut_D.theta()), marker='*', markersize=7, color='k') pl.plot(todeg( lep_D.phi() ), todeg(lep_D.theta()), marker='x', markersize=7, color='b') pl.plot(todeg( startDir.phi() ), todeg(startDir.theta()), marker='x', markersize=5, color='k') pl.plot(todeg( reco5d_D.phi()) , todeg(reco5d_D.theta()), marker='o', markersize=7, color='k') pl.plot(todeg( minLL_D.phi() ) , todeg(minLL_D.theta()), marker='*', markersize=5, color='r') pl.plot(todeg(reco3d_D.phi()) , todeg(reco3d_D.theta()), marker='o', markersize=4, color='g') for itrk, trk in enumerate(min_eval_trks): pl.plot( todeg(trk.dir.phi()), todeg(trk.dir.theta()), marker='.', markersize=3, color='r') pl.imshow(LL_5d, cmap='YlOrRd', origin='lower', extent=[-180, 180, 0, 180], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('Full 5d LL space, evt %i' %evt.id) pl.colorbar() pl.savefig(options.l + basefilename + "LL_TESTF_evt_%i" %evt.id) if options.s: pl.show() pl.clf() zmin_theta = todeg( neut_D.theta() ) - zscale zmax_theta = todeg( neut_D.theta() ) + zscale zmin_phi = todeg( neut_D.phi() ) - zscale zmax_phi = todeg( neut_D.phi() ) + zscale pl.plot(todeg( neut_D.phi() ), todeg(neut_D.theta()), marker='*', markersize=7, color='k') pl.plot(todeg( lep_D.phi() ), todeg(lep_D.theta()), marker='x', markersize=7, color='b') pl.plot(todeg( startDir.phi() ), todeg(startDir.theta()), marker='x', markersize=5, color='k') pl.plot(todeg( reco5d_D.phi()) , todeg(reco5d_D.theta()), marker='o', markersize=7, color='k') pl.plot(todeg( minLL_D.phi() ) , todeg(minLL_D.theta()), marker='*', markersize=5, color='r') pl.plot(todeg(reco3d_D.phi()) , todeg(reco3d_D.theta()), marker='o', markersize=4, color='g') for itrk, trk in enumerate(min_eval_trks): pl.plot( todeg(trk.dir.phi()), todeg(trk.dir.theta()), marker='.', markersize=3, color='r') pl.imshow(ZLL_5d, cmap='YlOrRd', origin='lower', extent=[zmin_phi, zmax_phi, zmin_theta, zmax_theta], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('%ix%i deg zoom, evt %i' %(zscale, zscale, evt.id) ) pl.colorbar() pl.savefig(options.l + basefilename + "ZLL_TESTF_evt_%i" %evt.id) if options.s: pl.show() pl.clf() zzmin_theta = todeg(reco5d_D.theta() ) - zzscale zzmax_theta = todeg(reco5d_D.theta() ) + zzscale zzmin_phi = todeg(reco5d_D.phi() ) - zzscale zzmax_phi = todeg(reco5d_D.phi() ) + zzscale pl.plot(todeg( neut_D.phi() ), todeg(neut_D.theta()), marker='*', markersize=7, color='k') pl.plot(todeg( lep_D.phi() ), todeg(lep_D.theta()), marker='x', markersize=7, color='b') pl.plot(todeg( startDir.phi() ), todeg(startDir.theta()), marker='x', markersize=5, color='k') pl.plot(todeg( reco5d_D.phi()) , todeg(reco5d_D.theta()), marker='o', markersize=7, color='k') pl.plot(todeg( minLL_D.phi() ) , todeg(minLL_D.theta()), marker='*', markersize=5, color='r') pl.plot(todeg(reco3d_D.phi()) , todeg(reco3d_D.theta()), marker='o', markersize=4, color='g') for itrk, trk in enumerate(min_eval_trks): pl.plot( todeg(trk.dir.phi()), todeg(trk.dir.theta()), marker='.', markersize=3, color='r') pl.imshow(ZZLL_5d, cmap='YlOrRd', origin='lower', extent=[zzmin_phi, zzmax_phi, zzmin_theta, zzmax_theta], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('%.1fx%.1f deg zoom, evt %i' %(zzscale, zzscale, evt.id) ) pl.colorbar() pl.savefig(options.l + basefilename + "ZZLL_TESTF_reco_evt_%i" %evt.id) if options.s: pl.show() pl.clf() zzzmin_theta = todeg(reco5d_D.theta() ) - zzzscale zzzmax_theta = todeg(reco5d_D.theta() ) + zzzscale zzzmin_phi = todeg(reco5d_D.phi() ) - zzzscale zzzmax_phi = todeg(reco5d_D.phi() ) + zzzscale pl.plot(todeg( neut_D.phi() ), todeg(neut_D.theta()), marker='*', markersize=7, color='k') pl.plot(todeg( lep_D.phi() ), todeg(lep_D.theta()), marker='x', markersize=7, color='b') pl.plot(todeg( startDir.phi() ), todeg(startDir.theta()), marker='x', markersize=5, color='k') pl.plot(todeg( reco5d_D.phi()) , todeg(reco5d_D.theta()), marker='o', markersize=7, color='k') pl.plot(todeg( minLL_D.phi() ) , todeg(minLL_D.theta()), marker='*', markersize=5, color='r') pl.plot(todeg(reco3d_D.phi()) , todeg(reco3d_D.theta()), marker='o', markersize=4, color='g') for itrk, trk in enumerate(min_eval_trks): pl.plot( todeg(trk.dir.phi()), todeg(trk.dir.theta()), marker='.', markersize=3, color='r') pl.imshow(ZZZLL_5d, cmap='YlOrRd', origin='lower', extent=[zzzmin_phi, zzzmax_phi, zzzmin_theta, zzzmax_theta], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('%.4fx%.4f deg zoom, evt %i' %(zzzscale, zzzscale, evt.id) ) pl.colorbar() pl.savefig(options.l + basefilename + "ZZZLL_TESTF_reco_evt_%i" %evt.id) if options.s: pl.show() pl.clf() marked_one_LL_5d = np.zeros( len(one_LL_5d) ) + one_LL_5d.mean() print("shape one ll: ", np.shape(one_LL_5d) ) for i, ill in enumerate(one_LL_5d): if i == 0 or i == len(one_LL_5d) - 1: continue grad_forward = one_LL_5d[i + 1] - ill grad_backward = one_LL_5d[i - 1] - ill if grad_forward == 0 or grad_backward == 0 or abs( log( abs( grad_forward/grad_backward ) ) ) > 3: marked_one_LL_5d[i] = one_LL_5d[i] pl.plot(todeg(one_PhiArray), one_LL_5d, 'b') pl.plot(todeg(one_PhiArray), marked_one_LL_5d, 'r') pl.xlabel(r'phi [deg]') pl.ylabel(r'LL') pl.title('1d #theta = #pi/2, evt %i' %evt.id) pl.legend() pl.savefig(options.l + basefilename + "one_LL_TESTF_reco_evt_%i" %evt.id) if options.s: pl.show() pl.clf() microscale = 0.004 micromin_theta = todeg(reco5d_D.theta() ) - microscale micromax_theta = todeg(reco5d_D.theta() ) + microscale micromin_phi = todeg(reco5d_D.phi() ) - microscale micromax_phi = todeg(reco5d_D.phi() ) + microscale Minuit_arr = np.zeros((pscale, pscale)) - 0.01 for itrk, trk in enumerate(min_eval_trks): itheta = ( todeg(trk.dir.theta()) - micromin_theta) * pscale / ( micromax_theta - micromin_theta ) iphi = ( todeg(trk.dir.phi()) - micromin_phi ) * pscale / ( micromax_phi - micromin_phi ) if abs(itheta) <= pscale and abs(iphi) <= pscale: Minuit_arr[itheta][iphi] = min_evals[itrk] - reco5d_LL pl.imshow(Minuit_arr, cmap='YlOrRd', origin='lower', extent=[micromin_phi, micromax_phi, micromin_theta, micromax_theta], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('0.004 x 0.004 deg zoom, evt %i' %evt.id) pl.colorbar() #pl.savefig(options.l + basefilename + "microLL_reco_evt_%i" %evt.id) if options.s: pl.show() pl.clf() nanoscale = 0.0001 nanomin_theta = todeg(reco5d_D.theta() ) - nanoscale nanomax_theta = todeg(reco5d_D.theta() ) + nanoscale nanomin_phi = todeg(reco5d_D.phi() ) - nanoscale nanomax_phi = todeg(reco5d_D.phi() ) + nanoscale Minuit_arr = np.zeros((pscale, pscale)) - 0.0004 for itrk, trk in enumerate(min_eval_trks): itheta = ( todeg(trk.dir.theta()) - nanomin_theta) * pscale / ( nanomax_theta - nanomin_theta ) iphi = ( todeg(trk.dir.phi()) - nanomin_phi ) * pscale / ( nanomax_phi - nanomin_phi ) if abs(itheta) <= pscale and abs(iphi) <= pscale: Minuit_arr[itheta][iphi] = min_evals[itrk] - reco5d_LL pl.imshow(Minuit_arr, cmap='YlOrRd', origin='lower', extent=[nanomin_phi, nanomax_phi, nanomin_theta, nanomax_theta], interpolation='none') #, vmin=true_LL*0.99, vmax=true_LL*1.01) pl.xlabel(r'phi [deg]') pl.ylabel(r'theta [deg]') pl.title('0.0001 x 0.0001 deg zoom, evt %i' %evt.id) pl.colorbar() #pl.savefig(options.l + basefilename + "nanoLL_reco_evt_%i" %evt.id) if options.s: pl.show() pl.clf() print "bye"