{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "from math import *\n", "import os, aa, rec\n", "import ROOT\n", "\n", "from ipywidgets import interact, interactive, fixed, interact_manual, Output, VBox, HBox\n", "import ipywidgets as widgets\n", "import IPython\n", "from IPython.display import clear_output, display\n", "from array import array\n", "from ctypes import string_at\n", "from copy import copy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# If you don't already have the pdf-data files, get them from Lyon with\n", "# cp username@cca7.in2p3.fr:/pbs/throng/km3net/src/Jpp/trunk/data/J*p.dat datadir\n", "\n", "datadir = os.environ['JPP_DIR']+'/data/'\n", "track_files = [ datadir + 'J%dp.dat'% n for n in range( 1,7)]\n", "shower_files = [datadir + 'J14p.dat', datadir + 'J13p.dat']\n", "aa_shower_file = [\"/sps/km3net/users/heijboer/he/output/genhen_JUNE_2___.hist.root\"]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "TTS = 0.0 #1.3\n", "jpp_shower_pdf = ROOT.JppShower( *shower_files, True, TTS,\n", " -100, 900, False, 20)\n", "#file_direct, file_scattered, TTS\n", "#jpp_shower_pdf = ROOT.JppShowerNpe( ROOT.JppShowerPdf( *shower_files ) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_args( logE, D, cd, angle_0, angle_1, dt, elong_steps = 0, icoord = 0):\n", " return { '0': (10**logE, D, cd, [angle_0, cos(angle_0)][icoord], angle_1, dt),\n", " 'S': (10**logE, D, cd, [angle_0, cos(angle_0)][icoord], angle_1, dt),\n", " 'E': (10**logE, D, cd, [angle_0, cos(angle_0)][icoord], angle_1, dt, elong_steps) }\n", "\n", "def get_pdf_info_dic( pdf, icoord, args ):\n", " return {\"jpp\": (\"jpp\", getattr(*( pdf, ['dnpe_dt', 'dnpe_dt_rad' ][icoord] )), args['0'], \"Jpp 4d\"),\n", " \"jpp_shower_max\": (\"jpp_shower_max\", getattr(*( pdf, ['dnpe_dt_S', 'dnpe_dt_S_rad'][icoord] )), args['S'], \"Jpp 4d_{shower max.}\"),\n", " \"jpp_elongated\": (\"jpp_elongated\", getattr(*( pdf, ['dnpe_dt_E', 'dnpe_dt_E_rad'][icoord] )), args['E'], \"Jpp 4d_{E elong.}\"),\n", " }\n", "\n", "def get_qvals( quantity, x0, x1, icoord ):\n", " #[index, mint, maxt, legend, xlabel]\n", " return { 'E' : [0, x0, 9*x1, [0.1, 0.6, 0.48, 0.9],'E [GeV]' ],\n", " 'D' : [1, x0, 10**(3*x1), [0.5, 0.6, 0.9 , 0.9],'D [m]' ],\n", " 'cd' : [2, x0, x1, [0.1, 0.6, 0.48, 0.9],'cos(#theta)' ],\n", " 'angle_0' : [3, [1,0][icoord]+x0, [pi,1][icoord]*x1, [0.1, 0.6, 0.48, 0.9],['#theta_{PMT}', 'cos_{PMT}'][icoord]],\n", " 'angle_1' : [4, 1 + x0, 2*pi*x1, [0.1, 0.6, 0.48, 0.9],['#phi_{PMT}' , 'irad' ][icoord]],\n", " 'dt' : [5, 100*x0, 1000*x1, [0.1, 0.6, 0.48, 0.9],'dt']\n", " }[quantity]\n", " \n", " \n", "color_scheme = {\n", " \"jpp\": [ 1, 34],\n", " \"jpp_shower_max\": [ 209, 211],\n", " \"jpp_elongated\": [97, 93],\n", " \"aanet\": [65, 65],\n", "}\n", "\n", "color_scheme[\"j\"] = color_scheme[\"o\"] = color_scheme[\"jpp\"] #j for jpp, o for original\n", "color_scheme[\"s\"] = color_scheme[\"c\"] = color_scheme[\"jpp_shower_max\"] #s for shower max, c for corrected\n", "color_scheme[\"e\"] = color_scheme[\"jpp_elongated\"] #e for elongated\n", "color_scheme[\"a\"] = color_scheme[\"aanet\"]\n", "\n", "\n", "\n", " \n", " \n", " \n", "def get_color( pdf_string, table = color_scheme): \n", " return table[pdf_string]\n", "\n", "def is_kinked(i0, i1, i2, threshold = 1):\n", " grad_0 = i1 - i0\n", " grad_1 = i2 - i1\n", " if(grad_0 == 0 or grad_1 == 0):\n", " return False\n", " return abs( log( abs( grad_0/grad_1 ) ) ) > threshold\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fill_histograms( histlist, b,\n", " pdf_info, quantity, qvals, logE, D, cd, dt,\n", " smax_coord= False,\n", " logx = False ):\n", " \n", " F = pdf_info[1] \n", "\n", " iq = histlist[0].GetBinCenter(b)\n", "\n", " new_D = copy(D)\n", " new_cd = copy(cd)\n", " new_dt = copy(dt)\n", " if smax_coord: #treat D_s, cd_s, dt_s and iq, as if coordinates to shower max, and perform transformations from that frame.\n", " D_s = copy(D)\n", " cd_s = copy(cd)\n", " dt_s = copy(dt)\n", " if quantity == \"D\" : \n", " iq = ROOT.get_vertex_D ( iq, cd_s, 10**logE )\n", " new_cd = ROOT.get_vertex_cd( iq, cd_s, 10**logE )\n", " new_dt = ROOT.get_vertex_dt( iq, cd_s, dt_s, 10**logE )\n", " elif quantity == \"cd\":\n", " new_D = ROOT.get_vertex_D ( D_s, iq, 10**logE )\n", " iq = ROOT.get_vertex_cd( D_s, iq, 10**logE )\n", " new_dt = ROOT.get_vertex_dt( D_s, iq, dt_s, 10**logE )\n", " elif quantity == \"dt\" : \n", " new_D = ROOT.get_vertex_D ( D_s, cd_s, 10**logE )\n", " new_cd = ROOT.get_vertex_cd( D_s, cd_s, 10**logE )\n", " iq = ROOT.get_vertex_dt( D_s, cd_s, iq, 10**logE )\n", " else:\n", " new_D = ROOT.get_vertex_D ( D_s, cd_s, 10**logE )\n", " new_cd = ROOT.get_vertex_cd( D_s, cd_s, 10**logE )\n", " new_dt = ROOT.get_vertex_dt( D_s, cd_s, dt_s, 10**logE )\n", "\n", " if quantity == \"E\" and logx: iq = 10**iq\n", "\n", "\n", " this_args = list(pdf_info[2])\n", "\n", " this_args[1] = new_D\n", " this_args[2] = new_cd\n", " this_args[5] = new_dt\n", " this_args[qvals[0]] = iq\n", "\n", " for icomp in range(3):\n", " histlist[icomp].SetBinContent(b, F(icomp, *this_args))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_pdf( pdf0 = \"jpp\",\n", " coord = \"Jpp\", \n", " quantity = \"cd\",\n", " smax_coord = False,\n", " logE = 6,\n", " D = 100,\n", " cd = 0.8,\n", " angle_0 = pi/2,\n", " angle_1 = pi/2,\n", " dt = 0,\n", " elong_steps = 10,\n", " nbins = 200 ,\n", " logy = False, logx = False,\n", " x0 = -1.0, x1 = 1.0,\n", " logy0 = 0.0, logy1 = 0.0,\n", " norm = False):\n", " \"\"\"\n", " Function for interactive plotting of the 4d Jpp Shower PDF using either Jpp (theta, phi) or Aanet (cos_pmt, irad) coordinates.\n", " Parameters:\n", " quantity: the quantity to be plotted on the x-axis (E, D, cd, angle_0, angle_1)\n", " logE: logarithm of energy [GeV]\n", " cd: cosine of reception angle\n", " coord: coordinates to be used (Jpp or Aanet)\n", " angle_0: theta of PMT in Jpp coordinates, PMT angle with vertex-DOM vector in Aanet coordinates\n", " angle_1: phi of PMT in Jpp coordinates, radians around PMT angle in Aanet coordinate\n", " \n", " nbins: number of bins to be computed\n", " logy: set y-axis to log scale\n", " x0, x1: x-axis range\n", " logy0, logy1: log of y-axis range\n", " norm: normalize integral to 1\n", " \"\"\"\n", " \n", " icoord = {\"Jpp\": 0, \"Aanet\": 1}[coord]\n", "\n", " args = get_args( logE, D, cd, angle_0, angle_1, dt, elong_steps, icoord )\n", " \n", " pdf_info = get_pdf_info_dic( jpp_shower_pdf, icoord, args )[ pdf0 ] \n", " \n", " qvals = get_qvals( quantity, x0, x1, icoord )\n", "\n", " histlist = [ ROOT.TH1D(\"h\"+str(i),\\\n", " \"JPP Pdf;\"+qvals[-1]+\";\",\\\n", " nbins,\\\n", " qvals[1],\\\n", " qvals[2]) for i in range( 3 ) ]\n", " \n", " for b in range(1, histlist[0].GetNbinsX()+1):\n", " fill_histograms( histlist, b, pdf_info, quantity, qvals, logE, D, cd, dt, smax_coord, logx)\n", " \n", " if norm:\n", " pdf_int = histlist[0].Integral()\n", " for i in range(3):\n", " histlist[i].Scale( 1.0 / float(pdf_int) )\n", "\n", " \n", " if ROOT.gPad:\n", " ROOT.gPad.SetLogy(logy)\n", " if quantity != \"E\": ROOT.gPad.SetLogx(logx)\n", " if logy : histlist[0].SetMinimum(1e-8)\n", "\n", " histlist[0].SetStats(False)\n", " histlist[0].GetXaxis().SetTitle(qvals[-1])\n", " histlist[0].GetXaxis().SetRangeUser(qvals[1], qvals[2])\n", " if logy1 > 0 :\n", " histlist[0].GetYaxis().SetRangeUser(10**logy0, 10**logy1)\n", " histlist[0].Draw(\"HISTL\")\n", " \n", " legendbox1 = [ qvals[3][0], qvals[3][1] + 0.15, qvals[3][2], qvals[3][3] ]\n", " legendbox2 = [ qvals[3][0], qvals[3][1], qvals[3][2], qvals[3][1] + 0.15 ]\n", " \n", " legend = ROOT.TLegend(*legendbox1)\n", " \n", " parleg = ROOT.TLegend(*legendbox2)\n", " parleg.AddEntry(\"\", \"10^{%i} GeV, %i m, cos#theta:%.2f, %s:%.2f, %s:%.2f, %s:%.2f\" %(logE, D, cd, ['#theta_{PMT}','cos_{PMT}'][icoord], [angle_0, cos(angle_0)][icoord], ['#phi_{PMT}','irad'][icoord], angle_1, 'dt', dt), \"\")\n", " parleg.SetMargin(0.05)\n", " parleg.SetBorderSize(0)\n", " parleg.SetFillColorAlpha(0, 0)\n", " \n", " histlist[0].SetLineColor(1)\n", " histlist[0].SetLineWidth(3)\n", " \n", " \n", " for i in range(3): \n", " col = [get_color('jpp')[0], get_color('jpp')[1], get_color('jpp')[1]][i]\n", " sty = [1, 2, 1][i] # dashed = scattered, solid = direct\n", " ltext = [\"Total\", \"1-scattered\", \"Direct\"][i]\n", " \n", " histlist[i].SetLineWidth(2)\n", " histlist[i].SetLineStyle(sty)\n", " histlist[i].SetLineColor(col)\n", " histlist[i].Draw(\"HISTLsame\")\n", "\n", " legend.AddEntry(histlist[i], ltext,\"l\")\n", " \n", " legend.Draw();\n", " parleg.Draw();\n", " \n", " with out:\n", " \n", " # create png in memory \n", " image = ROOT.TImage.Create()\n", " image.FromPad( ROOT.gPad )\n", " x = array('l', [0])\n", " y = array('i', [0])\n", " image.GetImageBuffer(x,y, ROOT.TImage.kPng)\n", " s = string_at( x[0],y[0]) \n", " img = IPython.display.Image( s, format='png' )\n", " \n", " display(img)\n", " del x,y,s,img \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "#compare Jpp and Jpp_e\n", "iplot = interactive(plot_pdf,\n", " pdf0=['jpp', 'jpp_shower_max', 'jpp_elongated'],\n", " coord=['Jpp','Aanet'],\n", " quantity=['E','D','cd', 'angle_0', 'angle_1', 'dt'],\n", " logE = (2.0,8.0), D = (2.0,500), cd = (-1,1,0.01),\n", " angle_0 = ( 0.0, pi, 0.02 ),\n", " angle_1 = ( 0.0, pi, 0.02 ),\n", " dt = ( -100.0, 1000.0, 0.1 ),\n", " elong_steps = (1, 100),\n", " nbins=(10,1000), logy = False,\n", " x0 = (-1, 1, 0.001), x1 = (-1, 1, 0.001),\n", " logy0 = (0.0, 5, 0.01), logy1 = (0.0, 5, 0.01) )\n", "\n", "# interactive plot has his own output widget, note we use it ('out') in plot_pdf\n", "leftbox = VBox( iplot.children[: int(len(iplot.children) / 2)])\n", "rightbox = VBox( iplot.children[int(len(iplot.children) / 2) : -1] )\n", "display( HBox([leftbox, rightbox]) )\n", "\n", "out = iplot.children[-1]\n", "out.layout = {\n", " 'width': '600px',\n", " 'height': '500px',\n", " 'border': '3px solid blue'\n", " }\n", "iplot\n", "display(out)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }