{"cells":[{"metadata":{"_uuid":"cba85599fe3a03aa87499bc640daa83794f9e07e"},"cell_type":"markdown","source":"![abc](http://www.pbs.org/wgbh/nova/next/wp-content/uploads/2015/03/cms-inner-tracker-barrel1.jpg)\nRetrieved from http://www.pbs.org/wgbh/nova/next/physics/lhc-accidental-rainbow-universe/"},{"metadata":{"_uuid":"abb279e3d7894930dd291f80fe7cfbcc3169a189"},"cell_type":"markdown","source":"# Objective\nThe host hands us a large dataset containing records of a bunch of particle detectors and wants us to perdict what kind of particle is hitting the detector. According to Pauli Exclusion Principle, no two fermion could occupy a same quantum state meaning every particle is unique just like us. Since the information shows the particle here is fermion, the objective is to label every jiterring paricle released in each collision event. If we see the problem in a semi-classical way, each particle would has it's unique tarjectory after bouncing off from the collision. Therefore if we could figure out a model to perdict the trajectory of the particle, we are indirectly labeling the particle."},{"metadata":{"_uuid":"8db0e1aa04117216e14a762b8eeecd98ba6c500c"},"cell_type":"markdown","source":"# Dependency\nRemember to add the trackml package to the notebook\n\nHow: \n\nclick the \">\" beside the \"Commit&Run\" botton \n    \ngo to the setting, Add a custom package\n\ntype \"LAL/trackml-library\" to the GitHub user/repo and click the arrow\n\nwait until done and restart the kernal by clicking the refrashing button at the bottom"},{"metadata":{"_uuid":"a9292d4f97f769a48542fa9bee73de027098bdd5"},"cell_type":"markdown","source":"# Note\n*  particle_id 0 in the \"truth\" file represents noise particles which has a relatively high momentum.\n*  column value in the \"cells\" file represents the signal value from the detector but the host do not left too many information to us (need to figure it out).\n* The apparatus is most like a pipe placed horizontally. \n* Detectors are in a multi-layers cylindrical arragement and are coaxial with the apparatus.\n* Each detector is a square tile. (Silicon Semiconductor?)\n* x, y axis is in the cross section (transverse) plane of the apparatus while z axis represent the \"width\" of the pipe\n* a strong magnetic field in z-axis bend the particle into a helix trajectory\n* mostly the particle is circulating in the x-y plane "},{"metadata":{"_uuid":"f0343c943f14d0d3937a58894ca3666d19f7d767"},"cell_type":"markdown","source":"# Helpful Resources:\n* Announcement from the host https://www.kaggle.com/c/trackml-particle-identification/discussion/55708\n* Guideline from the host https://kaggle2.blob.core.windows.net/forum-message-attachments/321278/9331/trackml-participant-document-particle-v1.0.pdf\n* List of Elementary Particle https://en.wikipedia.org/wiki/Elementary_particle\n* Useful Background Information by Heng CherKeng https://www.kaggle.com/c/trackml-particle-identification/discussion/55726#335835\n* DL Approachs\n    * Incorporating Deep Learning https://www.kaggle.com/c/trackml-particle-identification/discussion/57503#335377\n    * cone slicing, straightening helix and fitting https://www.kaggle.com/c/trackml-particle-identification/discussion/55726#335835\n* Clustering Approchs\n    * Flattening or unrolling tracks in polar (r,s,c) coordinates https://www.kaggle.com/c/trackml-particle-identification/discussion/58078#337285\n    * metric learning + clustering https://www.kaggle.com/c/trackml-particle-identification/discussion/57931#336566\n* Visualization & Dimension Reduction\n    * Transformation Visualization (fixed) osa111 https://www.kaggle.com/osa111/transformation-visualization-fixed\n    * analyzing results of LB 0.4922 https://www.kaggle.com/c/trackml-particle-identification/discussion/57947#336192\n    * Tilted Tracks https://www.kaggle.com/c/trackml-particle-identification/discussion/57931#336566"},{"metadata":{"_cell_guid":"d359b6c8-1803-4a89-888c-5273c501ead3","_uuid":"39eee1c9f51be754aa8282953e891b33f197b593"},"cell_type":"markdown","source":"In the training data we have the following information on each **event**:\n- **Hits**: $x, y, z$ coordinates of each hit on the particle detector\n- **Particles**: Each particle's initial position ($v_x, v_y, v_z$), momentum ($p_x, p_y, p_z$), charge ($q$) and number of hits\n- **Truth**: Mapping between hits and generating particles; the particle's trajectory, momentum and the hit weight\n- **Cells**: Precise location of where each particle hit the detector and how much energy it deposited"},{"metadata":{"_cell_guid":"089ad293-ca5d-4352-9e87-35300a757f3c","_uuid":"3d4c96fa4e9f4a55f70ac8fcf0e9073c18276757"},"cell_type":"markdown","source":"# Data Exploration:\n\n#### Import `trackml-library`\nThe easiest and best way to load the data is with the [trackml-library] that was built for this purpose.\n\nUnder your kernel's *Settings* tab -> *Add a custom package* -> *GitHub user/repo* (LAL/trackml-library)\n\nRestart your  kernel'\ns session and you will be good to go.\n\n[trackml-library]: https://github.com/LAL/trackml-library"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","scrolled":true,"trusted":true},"cell_type":"code","source":"import pdb\nimport os\nimport copy\n\nimport numpy as np\nimport pandas as pd\n\nfrom trackml.dataset import load_event\nfrom trackml.randomize import shuffle_hits\nfrom trackml.score import score_event\n\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nimport seaborn as sns\n%matplotlib inline\n\n#import plotly.plotly as py\nimport plotly.graph_objs as go\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\ninit_notebook_mode(connected=True)","execution_count":2,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dcd0381765fdefdfdb035bbd7be9a39ed6d244a3"},"cell_type":"code","source":"device = pd.read_csv('../input/detectors.csv')","execution_count":3,"outputs":[]},{"metadata":{"_cell_guid":"03e8af96-85ba-4d79-8601-453a32b8942e","_kg_hide-output":false,"_uuid":"59c07c9273966e79b26594a4804f209064a94048","trusted":true,"collapsed":true},"cell_type":"code","source":"event_prefix = 'event000001000'\nhits, cells, particles, truth = load_event(os.path.join('../input/train_1', event_prefix))\n\nmem_bytes = (hits.memory_usage(index=True).sum() \n             + cells.memory_usage(index=True).sum() \n             + particles.memory_usage(index=True).sum() \n             + truth.memory_usage(index=True).sum())\nprint('{} memory usage {:.2f} MB'.format(event_prefix, mem_bytes / 2**20))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ad9c44f5-4b12-4aac-a557-d16fbbaa28be","_uuid":"967a443483052831a49bc0a019b7766c972a287e"},"cell_type":"markdown","source":"## Tracjectory Data\n\n### the truth data\nHere we have: \n* $tx, ty, tz$ global coordinates (in millimeters) of where the particles hit the detector surface\n* $tp_x, tp_y, tp_z$ the momentum component  (in $GeV/c$) of the particles in each direction [c means light speed, eV is electron volt]\n* $weight$ the weight for evaluation the predicted particle_id. \n* $particle id$ our objective"},{"metadata":{"_uuid":"d231ee4740104951068e25d12478a80bba3e9a54","trusted":false},"cell_type":"code","source":"truth.sort_values(by=\"hit_id\",inplace=True)\ntruth.head()","execution_count":4,"outputs":[]},{"metadata":{"_uuid":"7e2322a7746df2d3651c8cc452f0ca6e27841280","collapsed":true},"cell_type":"markdown","source":"### Discovey1: the host don't care about the noise particles\nWe have not penalty from the noise."},{"metadata":{"_uuid":"18ffa99687df24ad551999c0ebd7a769cd01e4ef","trusted":false},"cell_type":"code","source":"truth['hasweight'] = ~np.equal( truth.weight.values, 0 )\ntruth[['particle_id','weight','hasweight']].head()","execution_count":5,"outputs":[]},{"metadata":{"_uuid":"d293ad5cdcd6551a441eb3c369021b54a2f482f5","trusted":false},"cell_type":"code","source":"# store the magnitude of the momentum\ntruth['tp'] = np.sqrt(truth['tpx']**2+truth['tpy']**2+truth['tpz']**2)\ntruth.head()","execution_count":6,"outputs":[]},{"metadata":{"_uuid":"71720a3a5315e1ad319c3075c9ad9ee1493577f3"},"cell_type":"markdown","source":"## the particles data\nHere we have: \n* $vx, vy, vz$ global coordinates of the initial position of the particles\n* $vp_x, vp_y, vp_z$ the initial momentum component of the particles in each direction.\n* $q$ the relative electric charge w.r.t e\n* nhits the number of hit to the detectors\n* $particle id$ our objective"},{"metadata":{"_uuid":"8ba204ce3cf5ed2560ab85f73d83b4c2aae0fdf0","trusted":false},"cell_type":"code","source":"particles.head()","execution_count":7,"outputs":[]},{"metadata":{"_kg_hide-output":true,"_uuid":"ae5e8c9019ea472432805d5c6752448d76a712b8","trusted":false},"cell_type":"code","source":"# add hit_id to help the merge in future\nparticles['hit_id'] = -1\nparticles['p'] = np.sqrt(particles['px']**2+particles['py']**2+particles['pz']**2)\nparticles.head()","execution_count":8,"outputs":[]},{"metadata":{"_uuid":"70cd4fd2de3bd23f2f58bb0b8cf53b9324a9a5aa"},"cell_type":"markdown","source":"### Discovery 2: there are several particles missing in the truth data\n* Some of the particle in \"particles\" could not be found in \"truth\". \n* the only one particle in \"truth\" cannot be found in \"particles\" is particle_id 0, noise particles."},{"metadata":{"_uuid":"ef0fa941c987bca8e5e6213351674f13f3173ce9","trusted":false,"collapsed":true},"cell_type":"code","source":"init_truth = particles.rename({\n    'vx' : 'tx',\n    'vy' : 'ty',\n    'vz' : 'tz',\n    'px' : 'tpx',\n    'py' : 'tpy',\n    'pz' : 'tpz',\n    'p'  : 'tp'\n}, axis=1)\ninit_truth.drop('nhits', axis=1, inplace=True)","execution_count":9,"outputs":[]},{"metadata":{"_uuid":"9661a078701f6e9429d5eef21838ee5c392af4c1","trusted":false},"cell_type":"code","source":"# the number of particles in the particles and truth\ninituni_par = init_truth.particle_id.unique()\nuni_par = truth.particle_id.unique()\nlen(inituni_par), len(uni_par)","execution_count":10,"outputs":[]},{"metadata":{"_uuid":"b31d587e4d2f5e29f29d982eb2abca5bbbb55fa7","trusted":false},"cell_type":"code","source":"# how many particles do they share\ninter = np.intersect1d(uni_par, inituni_par)\nlen(inter)","execution_count":11,"outputs":[]},{"metadata":{"_uuid":"06f29584cb664845d9cebcf5b8187e026bf19a18","trusted":false},"cell_type":"code","source":"# what is the one can't find in the particles\nnp.setdiff1d(uni_par, inter)","execution_count":12,"outputs":[]},{"metadata":{"_uuid":"a3f0d1576c9b6fe536426425f666671fc171a8e4"},"cell_type":"markdown","source":"### Let's merge truth and particles"},{"metadata":{"_uuid":"3ec05e8c3793c847749766ee1a30aaa0446fb831","trusted":false,"collapsed":true},"cell_type":"code","source":"weight_map = truth.groupby('particle_id').first()['weight']\ncharge_map = init_truth.groupby('particle_id').first()['q']","execution_count":13,"outputs":[]},{"metadata":{"_uuid":"ea50b8cbe7be7cb50271002254affa3fd6ccc0c2","trusted":false,"collapsed":true},"cell_type":"code","source":"truth['q'] = truth.particle_id.map(charge_map)\ntruth.fillna(0, inplace=True)","execution_count":14,"outputs":[]},{"metadata":{"_uuid":"1bbd6a77a1b83b3c4c2921a483df34ed67ca8a35","trusted":false,"collapsed":true},"cell_type":"code","source":"init_truth['weight'] = init_truth.particle_id.map(weight_map)\ninit_truth.fillna(0, inplace=True)\ninit_truth['hasweight'] = ~np.equal(init_truth.weight, 0)","execution_count":15,"outputs":[]},{"metadata":{"_uuid":"8e0bf0812e1554639169b0f7de15e5d31132ab80","trusted":false},"cell_type":"code","source":"truth.set_index('hit_id',inplace=True)\ninit_truth.set_index('hit_id',inplace=True)\nfulltruth = init_truth.append(truth, sort=True)\nfulltruth.head()","execution_count":16,"outputs":[]},{"metadata":{"_uuid":"2a9dcde5b9d9855d87b991f5efc3484a8758de15"},"cell_type":"markdown","source":"### Create a feature $R_r$ or relative rotation radius from the Equation:\n\\begin{equation}\n    m\\frac{v^2}{r}=qvB\n\\end{equation}\n\\begin{equation}\nr = \\frac{mv}{qB} = \\frac{p}{qB}\n\\end{equation}\n\n\\begin{equation}\nR_r = reB = \\frac{p}{q_r}\n\\end{equation}"},{"metadata":{"_uuid":"08d760ff04d4475af2c1f04ec280c4f29350e57c","trusted":false,"collapsed":true},"cell_type":"code","source":"fulltruth['R'] =  np.sqrt( fulltruth['tpx']**2 + fulltruth['tpy']**2 ) / fulltruth['q']","execution_count":17,"outputs":[]},{"metadata":{"_uuid":"1397799f4bc82137a872af3a8609491c57b2e9bf"},"cell_type":"markdown","source":"### Discovery 3: the particle the host do not care about is the one has extremely high momentum\nnote: the host do care about the particle has large momentum (>3GeV/c)"},{"metadata":{"_kg_hide-input":false,"_kg_hide-output":true,"_uuid":"572fe0781b421db7b0c49316527daf24e0d427f0","trusted":false,"collapsed":true},"cell_type":"code","source":"tgroups = truth.groupby(\"hasweight\")\n\nunrel = tgroups.get_group(False)\nunrel_ = unrel.groupby('particle_id').first()\n\nrel = tgroups.get_group(True)\nrel_ = rel.groupby('particle_id').first()","execution_count":18,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"aaaa9d50c5d729e555bdcc5396ff2114831dc022","trusted":false},"cell_type":"code","source":"fig, axs = plt.subplots(1,2, figsize=(18,6))\naxs[0].set_title('zero weight')\nsns.distplot(unrel.tp, hist=True, kde=False, ax = axs[0] )\naxs[1].set_title('has weight')\nsns.distplot(rel.tp, hist=True, kde=False, ax= axs[1])","execution_count":19,"outputs":[]},{"metadata":{"_uuid":"7e059c6ba2dd8b29fb881b4a60d65b7ff81c192b"},"cell_type":"markdown","source":"### Discovey 4: the data is pretty noisy \n15% of the data is the noise\n\n\nthe paricle_id is the one with generally high momentum"},{"metadata":{"_uuid":"65c934f22c4b18587a0a6958431580d70415dcb4","trusted":false},"cell_type":"code","source":"par0 = truth[truth.particle_id==0]\nlen(truth), len(par0), len(par0)/len(truth)    ","execution_count":20,"outputs":[]},{"metadata":{"_uuid":"48941b3420f9771a8a318e54991f1e62621ce04c","trusted":false},"cell_type":"code","source":"par0.tp.describe()","execution_count":21,"outputs":[]},{"metadata":{"_uuid":"7f4b9e1d8a91f2119eab24c7d7b49880c97cbda8"},"cell_type":"markdown","source":"## Let's plot the initial location of the particle\nAt here the marker color indicate the magnitude of the momentum "},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"d0c8cdbd0c971a0eaa6a45d083f09dd80f43270b"},"cell_type":"code","source":"stdlayout3d = dict(\n    width=800,\n    height=700,\n        \n    autosize=False,\n    title= 'unknown',\n    scene=dict(\n        xaxis=dict(\n            title = \"unknown x\",\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        yaxis=dict(\n            title = \"unknown y\",\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        zaxis=dict(\n            title = \"unknown z\",\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        camera=dict(\n            up=dict(x=0, y=0, z=1),\n            eye=dict(x=-1.7428, y=1.0707, z=0.7100,)\n        ),\n        aspectratio = dict(x=1, y=1, z=0.7),\n        aspectmode = 'manual'\n    ),\n)\n\ndef layout_costom3d(xtitle, ytitle, ztitle, title, xrange=None, yrange=None, zrange=None):\n    layout = copy.deepcopy(stdlayout3d)\n    layout['scene']['xaxis']['title'] = xtitle \n    layout['scene']['yaxis']['title'] = ytitle\n    layout['scene']['zaxis']['title'] = ztitle\n    if xrange is not None: layout['scene']['xaxis']['range'] = xrange\n    if yrange is not None: layout['scene']['yaxis']['range'] = yrange\n    if zrange is not None: layout['scene']['zaxis']['range'] = zrange\n    layout['title'] = title\n    return layout","execution_count":186,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"6d7d3a1e83cde1aa32c20c150c413605af152107"},"cell_type":"code","source":"stdlayout2d = dict(\n    height = 800,\n    width = 800,\n    title = 'unknown',\n    yaxis = dict(),\n        #zeroline = False,\n\n    xaxis = dict(),\n        #zeroline = False,\n)    \ndef layout_costom2d(xtitle, ytitle, title, xrange=None, yrange=None, zrange=None):\n    layout = copy.deepcopy(stdlayout2d)\n    layout['xaxis']['title'] = xtitle \n    layout['yaxis']['title'] = ytitle    \n    if xrange is not None: layout['xaxis']['range'] = xrange\n    if yrange is not None: layout['yaxis']['range'] = yrange\n    layout['title'] = title\n    return layout   \n","execution_count":187,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"9ee18a23f929a3db257ed852b8afa2a27f2e045e"},"cell_type":"code","source":"xyzlayout = layout_costom3d('z axis(mm)', 'x axis(mm)', 'y axis(mm)', 'sample trajectories')\nxylayout = layout_costom2d('x axis(mm)', 'y axis(mm)', 'sample trajectories', [-1000, 1000], [-1000, 1000])","execution_count":188,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"18e1ba172fd1dedb7aff81a383275a9922b2d542"},"cell_type":"code","source":"def set_marker(pp, pp_name, isnoise, ms, cmin, cmax):\n    marker = dict(\n        size=ms,\n        symbol= \"square\" if isnoise else \"circle\",\n    )\n    #pdb.set_trace()\n    if pp is not None:\n        marker['color'] = pp\n        marker['colorscale']='Rainbow'\n        marker['colorbar']=dict(\n                title = pp_name,\n                x = 1.20\n            )\n        marker['cmin'] = cmin\n        marker['cmax'] = cmax\n    return marker\n\ndef set_line( width=1):\n    line = dict(\n        width=1\n    )","execution_count":89,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"9766e245975ea91fad051c6d592ae2d240d4b2f1"},"cell_type":"code","source":"def plotly_3d(x, y, z, pp=None, pp_name=None, isnoise=False, visible=\"legendonly\", mode=None, ms=4, cmin=0, cmax=1):     \n    marker = set_marker(pp, pp_name, isnoise, ms, cmin, cmax)        \n    trace = go.Scatter3d(\n        mode=mode,\n        visible=visible,\n        x=x, y=y, z=z,\n        marker= marker,\n        line = dict(width=1)\n    )\n    return [trace]\n\ndef plot_df(df, xyzcols, n=10, pproperty=None, pids=None, visible=\"legendonly\", mode=None, ms=2, cmin=0, cmax=1):\n    particlegroup = df.groupby('particle_id')\n    if pids is None:\n        pids = np.random.choice(df.particle_id.unique(), n)\n    \n    particles = [particlegroup.get_group(pid) for pid in pids]\n    \n    data = []\n    xc, yc, zc = xyzcols\n    for particle in particles:\n        trace=plotly_3d(\n            x=particle[xc],\n            y=particle[yc],\n            z=particle[zc],\n            pp=particle[pproperty] if pproperty is not None else None,\n            pp_name=pproperty,\n            isnoise=particle.weight.values[0] == 0,\n            visible=visible,\n            mode=mode,\n            ms=ms, cmin=cmin, cmax=cmax\n        )\n        data+=trace\n    return data","execution_count":136,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"0bba20e85516e7215a56c0476bdfeafbc1a6070b"},"cell_type":"code","source":"def plotly_2d(x, y, pp=None, pp_name=None, isnoise=False, visible=\"legendonly\", mode=None, ms=4, cmin=0, cmax=1):\n        \n    marker = set_marker(pp, pp_name, isnoise, ms, cmin, cmax)       \n    trace = go.Scatter(\n        x = x,\n        y = y,\n        mode = mode,\n        marker = marker,\n        line = dict(width=1)\n    )\n    return [trace]\n        \ndef plot_df2d(df, xycols, n=10, pproperty=None, pids=None, visible=\"legendonly\", mode=None, ms=4, cmin=0, cmax=1):\n    particlegroup = df.groupby('particle_id')\n    if pids is None:\n        pids = np.random.choice(df.particle_id.unique(), n)\n    \n    particles = [particlegroup.get_group(pid) for pid in pids]\n        \n    data = []\n    xc, yc = xycols\n    for particle in particles:\n        trace=plotly_2d(\n            x=particle[xc],\n            y=particle[yc],\n            pp=particle[pproperty] if pproperty is not None else None,\n            pp_name=pproperty,\n            isnoise=particle.weight.values[0] == 0,\n            visible=visible,\n            mode=mode,\n            ms=ms, cmin=cmin, cmax=cmax\n        )\n        data+=trace\n    return data ","execution_count":149,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"0a04d26fd608fa017110b3dbac0b26ef2a7f7d72","trusted":false},"cell_type":"code","source":"data = plotly_3d(particles.vz, particles.vx, particles.vy, particles.p, 'momentum', visible=True, ms=1, mode='markers')\nlayout = xyzlayout.copy()\nlayout['title'] = 'initial position'\niplot(dict(data=data, layout=layout), filename='local')","execution_count":104,"outputs":[]},{"metadata":{"_uuid":"754b9cb6d12be22c1ce2c4b1a5e3180592095a42"},"cell_type":"markdown","source":"## Let's plot the trajectory"},{"metadata":{"_uuid":"282cb1a87fa54ead69f50dfd76d944649e0c371d","collapsed":true},"cell_type":"markdown","source":"### plot 10 randomly picked particles "},{"metadata":{"_kg_hide-input":true,"_uuid":"ee9017580f13f30c9e894cc8f42f24ae7acfe42d","scrolled":true,"trusted":false},"cell_type":"code","source":"data = plot_df(truth, ['tz','tx','ty'], n=10, visible=True)\niplot(dict(data=data, layout=xyzlayout))","execution_count":105,"outputs":[]},{"metadata":{"_uuid":"529801fe7a43c8188956570c4322ed47e0cee673","collapsed":true},"cell_type":"markdown","source":"### plot the noise"},{"metadata":{"_kg_hide-input":true,"_uuid":"9b3d894f5b467630356cf99b26edff32d55bf0db","scrolled":true,"trusted":false},"cell_type":"code","source":"data = plot_df(truth, ['tz','tx','ty'], pids=[0], visible=True, mode='markers', ms=1)\nlayout = xyzlayout.copy()\nlayout['title'] = 'Noise'\niplot(dict(data=data, layout=layout))","execution_count":109,"outputs":[]},{"metadata":{"_uuid":"e877a3377df5dcf646ffdcebd690d0f58a30ffc0"},"cell_type":"markdown","source":"### plot the trajectory projection on xy plane"},{"metadata":{"_kg_hide-input":true,"_uuid":"11775820a618e6fddf6d660fa2d481476a27522e","trusted":false},"cell_type":"code","source":"data = plot_df2d(truth, ['tx', 'ty'], n=100, mode=None)\niplot(dict(data=data, layout=xylayout))","execution_count":127,"outputs":[]},{"metadata":{"_uuid":"36ebee2c8bfd0f8bdd9568fd9ad504b86af45672"},"cell_type":"markdown","source":"## the hits data\nHere we have: \n* $x, y, z$ global coordinates of the initial position of the particles\n* volume_id, layer_id and the module_id claim the arrgement information of each device"},{"metadata":{"_uuid":"e4d7268cfcd66de7e0494be8f2f90445b43b4b4d","trusted":false},"cell_type":"code","source":"hits.head()","execution_count":25,"outputs":[]},{"metadata":{"_uuid":"715495c155c85485d5d16b66f165e0df26e1a59e"},"cell_type":"markdown","source":"## the cells data\nHere we have: \n* $ch0, ch1$ global coordinates of the initial position of the particles\n* volume_id, layer_id and the module_id claim the arrgement information of each device\n"},{"metadata":{"_uuid":"71185906436937c56fd698e3ae01d80b9ddf4af1","trusted":false},"cell_type":"code","source":"cells.head()","execution_count":26,"outputs":[]},{"metadata":{"_uuid":"432b5b67437c860150dd9695fcb8f60258d1394c"},"cell_type":"markdown","source":"For one hit_id, there are more than two recorded data. We take the max to find the most sensitive cell."},{"metadata":{"_kg_hide-output":false,"_uuid":"4790163d4b4254035b75ab41100978e61ebf32ca","trusted":false},"cell_type":"code","source":"cells_ = cells.set_index('hit_id')\ncells_.drop(['ch0','ch1'], axis=1, inplace=True)\ncells_ = cells_.groupby('hit_id').agg('sum')\ncells_.head()","execution_count":27,"outputs":[]},{"metadata":{"_uuid":"3df0532f221f02e189bb294866f695bc3073e1c1"},"cell_type":"markdown","source":"### We merge hits, cells, particles, and truth into one dataframe"},{"metadata":{"_uuid":"6376ae0c32ba79d2a37ec54360dfa3983800fa52","trusted":false},"cell_type":"code","source":"hits_ = hits.set_index('hit_id')\ninfo = cells_.join(hits_)\ncheat = info.join(fulltruth)\ncheat.dropna(axis=0)\ncheat.R.replace(np.inf, 10000, inplace=True)\ncheat.head()","execution_count":28,"outputs":[]},{"metadata":{"_uuid":"dc7237dcc99d4d13f703923fa6992e7aa508ec8a"},"cell_type":"markdown","source":"## The distribution of signal value for different particle "},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"3fae461c07f1e106e515078547f30a4e3a6b7d08","trusted":false},"cell_type":"code","source":"pids = cheat.particle_id.unique()\nsamples = np.random.choice(pids, 2)\nfor sample in samples:\n    if sample != 0:   \n        ax = sns.distplot( cheat[cheat.particle_id==sample].value , kde=False, bins=np.linspace(0,1,10) )\nax.set_xlim([0,1])","execution_count":29,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"cbb21f1de58431651847e5bd27b3600f70baa56c","trusted":false},"cell_type":"code","source":"fig, axs =plt.subplots(1,3,figsize=(18,4))\n\nax = sns.distplot(cheat.value, ax=axs[0], kde=False, bins=np.linspace(0,1,20))\nax.set_xlim([0,1])\nax.set_title('overall')\nax = sns.distplot(cheat[cheat.particle_id == 0].value, ax=axs[1], kde=False, bins=np.linspace(0,1,20))\nax.set_xlim([0,1])\nax.set_title('noise')\nax = sns.distplot(cheat[~(cheat.particle_id == 0)].value, ax=axs[2], kde=False, bins=np.linspace(0,1,20))\nax.set_xlim([0,1])\nax.set_title('not noise');","execution_count":30,"outputs":[]},{"metadata":{"_uuid":"067bddaae566de43f7dcdde9d68fb4177cbc3eda"},"cell_type":"markdown","source":"## Let's plot the trajectory and vary the marker's color with respect to the particle property"},{"metadata":{"_uuid":"76d086418f7971d90d635ebffc644a0bfab5d1fe"},"cell_type":"markdown","source":"### w.r.t the signal value (i.e. value column)"},{"metadata":{"_kg_hide-input":true,"_uuid":"f55e9849923f552d99ebbb65f1847f11e33435d6","trusted":false},"cell_type":"code","source":"data = plot_df(cheat, ['tx', 'ty', 'tz'], n=30 ,pproperty='value', visible=True)\niplot(dict(data=data, layout=xyzlayout))","execution_count":129,"outputs":[]},{"metadata":{"_uuid":"e1781ebdbacc58e96542eb1098eb51ea125082c6"},"cell_type":"markdown","source":"#### the trajectory with high signal value (>0.5) "},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"34dc6a64270adf5ee58bad52db65ae21ed1aa27c","trusted":false},"cell_type":"code","source":"highsig = cheat.value > 0.5\nhighsigpart = cheat.particle_id[highsig].unique()\nlowsigpart = np.setdiff1d(cheat.particle_id,highsigpart)\ncheat.particle_id.nunique(), len(highsigpart), len(lowsigpart)","execution_count":34,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"20a634ae16553ffb2824b19a9fa2ddac0628e6b9","trusted":false},"cell_type":"code","source":"data = plot_df(cheat, ['tx', 'ty', 'tz'], pids=np.random.choice(highsigpart, 10), pproperty='value', visible=True)\niplot(dict(data=data, layout=xyzlayout))","execution_count":139,"outputs":[]},{"metadata":{"_uuid":"77457246ca4cbe16d298c624ca82c7405c201ce6"},"cell_type":"markdown","source":"#### trajectory with generally low signal value (>0.5) "},{"metadata":{"_kg_hide-input":true,"_uuid":"30868e7fe1822e14dc919bb059262ef7792994ad","trusted":false},"cell_type":"code","source":"data = plot_df(cheat, ['tx', 'ty', 'tz'], pids=np.random.choice(lowsigpart, 10), pproperty='value', visible=True)\niplot(dict(data=data, layout=xyzlayout))","execution_count":140,"outputs":[]},{"metadata":{"_uuid":"379ae767a0ed27e7acd7d81e605a9ddde1556cd6"},"cell_type":"markdown","source":"### w.r.t the magnitude of momentum (i.e. tp column)"},{"metadata":{"_kg_hide-input":true,"_uuid":"e7e2f67e2b9e9764e591d20e448242ee846b967c","trusted":false},"cell_type":"code","source":"data = plot_df(cheat, ['tx', 'ty', 'tz'], pproperty='tp', visible=True, cmin=0, cmax=2, ms=2)\niplot(dict(data=data, layout=xyzlayout))","execution_count":145,"outputs":[]},{"metadata":{"_uuid":"d8b1ef2f9118b99c40fc7d34f2e9da6907da53b0"},"cell_type":"markdown","source":"### w.r.t particle relative charge (i.e. q)\nWe could clearly see charge affect the bending direction in the second plot."},{"metadata":{"_kg_hide-input":true,"_uuid":"93be8fdc0db29d96f1e4acf882ecfb6739738b6d","trusted":false},"cell_type":"code","source":"data = plot_df(cheat, ['tx', 'ty', 'tz'], pproperty='q', visible=True, cmin=-1, cmax=1, ms=4)\niplot(dict(data=data, layout=xyzlayout))","execution_count":146,"outputs":[]},{"metadata":{"_uuid":"61b719fc441a7fd7609b6714eb87cabeadbaa1bf"},"cell_type":"markdown","source":"#### Projections on xy plane"},{"metadata":{"_kg_hide-input":true,"_uuid":"720f6633f94ddbb6d77e82a78c3d7f3bdd9ba121","trusted":false},"cell_type":"code","source":"data = plot_df2d(cheat, ['tx', 'ty'], pproperty='q', visible=True, cmin=-1, cmax=1, ms=6, mode='markers')\niplot(dict(data=data, layout=xylayout))","execution_count":156,"outputs":[]},{"metadata":{"_uuid":"252debf2f6293229eeef23fbd09ac74aef69c1d7"},"cell_type":"markdown","source":"### w.r.t relative rotation radius on xy plane"},{"metadata":{"_uuid":"e84344d177cbc43bdc006303b08d1ddfa3bf3524","_kg_hide-input":true,"trusted":false},"cell_type":"code","source":"data = plot_df2d(cheat, ['tx', 'ty'], pproperty='R', visible=True, cmin=-1, cmax=1, ms=6, mode='markers')\niplot(dict(data=data, layout=xylayout))","execution_count":158,"outputs":[]},{"metadata":{"_uuid":"c881481ba5bfefc23aef422b7181698a8149f48e","collapsed":true},"cell_type":"markdown","source":"## Convert the Cartesian Coordinatea to Cylindrical Coordinates\nIt provide a more straight forward visualization (trajactory: helix curve -> straight line)"},{"metadata":{"_kg_hide-input":false,"trusted":false,"collapsed":true,"_uuid":"24fb418bb836c8067541ae5f62e316f1afe2bf48"},"cell_type":"code","source":"def xyz2c(df, cols, newcols):\n    x, y, z = df[cols[0]], df[cols[1]], df[cols[2]]\n    #pdb.set_trace()\n    cr, cpsi, cz = newcols\n    r = np.sqrt( x**2 + y**2 )\n    psi = np.arctan(y/x)\n    df[cr] = r\n    df[cpsi] = psi\n    df[cz] = z","execution_count":35,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"ab6324f0e575261601e8755e0bbf20501cc368c6"},"cell_type":"code","source":"def xyz2s(df, cols, newcols):\n    x, y, z = df[cols[0]], df[cols[1]], df[cols[2]]\n    cr, cpsi, ctheta = newcols\n    r = np.sqrt( x**2 + y**2 + z**2 )\n    psi = np.arctan(y/x)\n    theta = np.arccos(z/r)\n    df[cr] = r\n    df[cpsi] = psi\n    df[ctheta] = theta","execution_count":36,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"df62cb66b222443ccc3cb0ed219ca5e6eb2375f6"},"cell_type":"code","source":"xyz2c(cheat, ['tx','ty','tz'], ['trc','tpsic','tzc'])\nxyz2s(cheat, ['tx','ty','tz'], ['trs','tpsis','tthetas'])","execution_count":163,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"37839cd547a270d7b10d85a7b71700bef3d4009d"},"cell_type":"code","source":"## Plotting the trajectory in new coordinates","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"3c851c58ab17e23eeed841ddddd05a99590d4823"},"cell_type":"code","source":"### in Spherical Coordinates","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"177a41a20b18b394ec79ec9de6305c46ccc3e0db"},"cell_type":"code","source":"spherelayout = layout_costom3d('theta (radians)','radius (mm)', 'psi (radians)', 'sample trajectories in spherical coor',xrange=[0, np.pi], zrange=[-np.pi, np.pi])","execution_count":181,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":false,"_uuid":"abda9ac2b753a092fed103960f90026f283f0648"},"cell_type":"code","source":"data = plot_df(cheat, ['tthetas', 'trs', 'tpsis'], visible=True)\niplot(dict(data=data, layout=spherelayout))","execution_count":169,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"64bf882a76f70c8de96fd4f3be9ff877b8724291"},"cell_type":"code","source":"### in Cylindricall Coordinates","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"f9e828c15c5a0bea97427014eca2d90ea609fe3f"},"cell_type":"code","source":"cylindricallayout = layout_costom3d('z (mm)','radius (mm)', 'psi (radians)', 'sample trajectories in cylindrical coor', zrange=[-np.pi, np.pi])","execution_count":189,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":false,"_uuid":"15e8d3cebf2a5c851d5f160117e09d65ffcf792a"},"cell_type":"code","source":"data = plot_df(cheat, ['tz', 'trc', 'tpsic'], visible=True)\niplot(dict(data=data, layout=cylindricallayout))","execution_count":191,"outputs":[]},{"metadata":{"_uuid":"84043ec9eb7c98aef121437c02ababc832a79275"},"cell_type":"markdown","source":"### in the r-z Plane of the Cylinderical Coordinate + detectors' location\nAccording to several discussions about transformation, visualizes the trajectory in this plane is seemingly the most promisible way. (links at the top of the notebook)\n\nWe could also plot trajectories with the location of detectors to create a better picture."},{"metadata":{"_kg_hide-input":false,"trusted":false,"collapsed":true,"_uuid":"daa573825ea18838affcb2871adef7f790c88322"},"cell_type":"code","source":"xyz2c(device, ['cx','cy','cz'], ['rc','psic','zc'])\nxyz2s(device, ['cx','cy','cz'], ['rs','psis','thetas'])","execution_count":49,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":false,"collapsed":true,"_uuid":"6845b02c63448b04ce7302e00e4718258f26a247"},"cell_type":"code","source":"zrlayout = layout_costom2d('z (mm)','radius (mm)', 'sample trajectories on rz plane')","execution_count":194,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":false,"_uuid":"d661aa5c0bcb955079197268154213c3ef9c4228"},"cell_type":"code","source":"data = plotly_2d(x=device.zc, y=device.rc, isnoise=True, ms=5, mode='markers')\ndata += plot_df2d(cheat, ['tzc','trc'])\niplot(dict(data=data, layout=zrlayout))","execution_count":197,"outputs":[]}],"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.5.2"}},"nbformat":4,"nbformat_minor":1}