{"cells":[{"metadata":{"_uuid":"4bca6f1dcac7f319922fa52e5fdeab2dad9b3788"},"cell_type":"markdown","source":"# Introduction to TrackML Challenge"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","collapsed":true,"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport matplotlib.pylab as plt\nfrom mpl_toolkits import mplot3d\nfrom mpl_toolkits.mplot3d import Axes3D\n%matplotlib inline\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport trackml\nfrom trackml.dataset import load_event","execution_count":26,"outputs":[]},{"metadata":{"_uuid":"62d9dac2d403fa8ebfa6969dc1006ba7840ba13b"},"cell_type":"markdown","source":"### Look for the information out of the event by reading all the information (hits/cells/particles/truths) from the event."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"cFirstEvent=1010\ncEventDataDir='../input/train_1'\ndef getPath(pDataDir,pEventID) : \n    return '%s/event%09d' % (pDataDir, pEventID)\n\n\nhits, cells, particles, truth = load_event(getPath(cEventDataDir,cFirstEvent))\nparticles.head()","execution_count":27,"outputs":[]},{"metadata":{"_uuid":"af9ee46417f5799289d5369c7a2cdc9324465744"},"cell_type":"markdown","source":"### Now simply look at the hit information in the (x,y,z) coordinate system. We can also do something like only look at hits in a given layer/volume/module. This block of codes returns a plot showing the hits in (y,x) and (r,z) coordinates for a given volume id."},{"metadata":{"trusted":true,"_uuid":"4e3797e86d732013d2ead630262f3f1b37ec9bce"},"cell_type":"code","source":"from pandas.plotting import scatter_matrix\nimport matplotlib.pyplot as plt\nimport math as math\n# Load in convex hull method\nfrom scipy.stats.stats import pearsonr\nfrom scipy.spatial import ConvexHull\n#circle\nfrom scipy import optimize\n\nnMinHits=5\n#draw track + hits \ndef getTrackParameters(pIndex) : \n    dataFrame = pd.DataFrame(particles)\n    \n# start with those that have 5 hits \ndef getTracks(sampleSize) : \n    dataFrame = pd.DataFrame(particles)\n    dataFrame = dataFrame[dataFrame['nhits']>=nMinHits]\n    # get unique list of particle IDs \n    particle_IDs = np.random.choice(dataFrame.particle_id.unique(),sampleSize)\n    print(particle_IDs)\n    dataFrame = pd.DataFrame(truth)\n    df_truth = dataFrame[dataFrame['particle_id'].isin(particle_IDs)]\n    return df_truth\n\ndef getHitsFromTracks(df_truth, sampleSize) : \n    dataFrame = pd.DataFrame(hits)\n    df_hits = dataFrame[dataFrame['hit_id'].isin(df_truth.hit_id)]\n    return df_hits\n\ndef getOtherHits(df_truth, sampleSize) : \n    dataFrame = pd.DataFrame(hits)\n    df_hits = dataFrame[dataFrame['hit_id'].isin(df_truth.hit_id)== False]\n    return  df_hits.sample(n=sampleSize)\n\n#return truths for a given particle \ndef getTruth(pTruths, particleID) :\n    dataFrame = pd.DataFrame(pTruths)\n    df_t = dataFrame[dataFrame['particle_id'] == particleID]\n    return df_t\n\n\n#return hits in a given volume \ndef getHitsForVolume(pHits, pVolumeID) : \n    dataFrame = pd.DataFrame(pHits)\n    df_v = dataFrame[dataFrame['volume_id'] == pVolumeID]\n    #df_v = df_v[df_v['layer_id'] < 6]\n    return df_v\n\n#return hits in a given volume \ndef getHitsForVolume_perLayer(pHits, pVolumeID, pLayerID) : \n    dataFrame = pd.DataFrame(pHits)\n    df_v = dataFrame[dataFrame['volume_id'] == pVolumeID]\n    df_v = df_v[df_v['layer_id'] == pLayerID]\n    return df_v\n\n# make things look familiar...\n#plots hits in (x,y) [cartesian] and (z,r) coordinate system [cylindrical]\ndef showHitsForVolume(pHits, pVolumeID) : \n    df_v = getHitsForVolume(pHits,pVolumeID)   \n    #now estimate r-coordinate (in x,y plane)\n    r = (df_v.x**2 + df_v.y**2)**0.5\n    phi = np.arctan(df_v.y/df_v.x)\n    plt.figure(1)\n    plt.subplot(121)\n    plt.plot(df_v.x,df_v.y, 'bs')\n    plt.xlabel('x [cm]')\n    plt.ylabel('y [cm]')\n\n    plt.subplot(122)\n    plt.plot(df_v.z,r, 'bs')\n    plt.xlabel('z [cm]')\n    plt.ylabel('r [cm]')\n    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=1.55, hspace=0.25, wspace=0.35)\n    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=1.55, hspace=0.25, wspace=0.35)\n    \n    return plt\n\ndef showHitsForVolume_perLayer(pHits, pVolumeID, pLayerID) : \n    df_v = getHitsForVolume_perLayer(pHits,pVolumeID,pLayerID)   \n    #now estimate r-coordinate (in x,y plane)\n    r = (df_v.x**2 + df_v.y**2)**0.5\n    phi = np.arctan(df_v.y/df_v.x)\n    plt.figure(1)\n    plt.subplot(121)\n    plt.plot(df_v.x,df_v.y, 'bs')\n    plt.xlabel('x [cm]')\n    plt.ylabel('y [cm]')\n\n    plt.subplot(122)\n    plt.plot(df_v.z,r, 'bs')\n    plt.xlabel('z [cm]')\n    plt.ylabel('r [cm]')\n    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=1.55, hspace=0.25, wspace=0.35)\n    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=1.55, hspace=0.25, wspace=0.35)\n    \n    return plt\n\ndef showHitsForParticle(pTruth,particleID) : \n    df_t = getTruth(pTruth,particleID)\n    r = (df_t.tx**2 + df_t.ty**2)**0.5\n    plt.figure(1)\n    plt.subplot(121)\n    plt.plot(df_t.tx,df_t.ty, 'bs')\n    plt.xlabel('x [cm]')\n    plt.ylabel('y [cm]')\n    \n    plt.subplot(122)\n    plt.plot(df_t.tz,r, 'bs')\n    plt.xlabel('z [cm]')\n    plt.ylabel('r [cm]')\n    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=1.55, hspace=0.25, wspace=0.35)    \n    return plt\n\ndef draw(x,y) : \n    plt.figure(1)\n    plt.plot(x,y, 'bs')\n    plt.xlabel('x [cm]')\n    plt.ylabel('y [cm]')\n    \n    return plt \n\nnTrueTracks=1\nnFakeHits=5\ndh = pd.DataFrame(hits)\ndh = dh[np.fabs(dh['z']) < 1]\nd_t = getTracks(nTrueTracks)\nd_ht = getHitsFromTracks(d_t,nTrueTracks)\nd_hf = getOtherHits(d_t,nFakeHits)\nr_ht = np.sqrt(d_ht.x**2 + d_ht.y**2)\nd_ht['r'] = r_ht\nr_hf = np.sqrt(d_hf.x**2 + d_hf.y**2)\nd_hf['r'] = r_hf\nd = pd.concat([d_ht, d_hf])\nplt.plot(dh.x,dh.y,'or')","execution_count":54,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"8f9a2cb862102b21337337f790c461b89e163553"},"cell_type":"code","source":"wdir = os.getcwd()\n\nhits_cols = \"hit_id,x,y,z,volume_id,layer_id,module_id,event_name\"\nparticle_cols = \"particle_id,vx,vy,vz,px,py,pz,q,nhits,event_name\"\ntruth_cols = \"hit_id,particle_id,tx,ty,tz,tpx,tpy,tpz,weight,event_name\"\ncells_cols = \"hit_id,ch0,ch1,value,event_name\"\n\nhits_df = pd.DataFrame(columns = hits_cols.split(\",\"))\nparticle_df = pd.DataFrame(columns=particle_cols.split(\",\"))\ntruth_df =  pd.DataFrame(columns = truth_cols.split(\",\"))\ncells_df = pd.DataFrame(columns= cells_cols.split(','))","execution_count":29,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3c62ecf4510be2d1c3c43584f9a0e68d884737cb"},"cell_type":"code","source":"def calc_R(xc, yc):\n    \"\"\" calculate the distance of each 2D points from the center (xc, yc) \"\"\"\n    return np.sqrt((x-xc)**2 + (y-yc)**2)\n\ndef f_2(c):\n    \"\"\" calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) \"\"\"\n    Ri = calc_R(*c)\n    return Ri - Ri.mean()\n\nx = d['x']\ny = d['y']\nx_m = np.mean(x)\ny_m = np.mean(y)\n\ncenter_estimate = x_m,y_m\ncenter_2, ier = optimize.leastsq(f_2, center_estimate)\n\nxc_2, yc_2 = center_2\nRi_2       = calc_R(*center_2)\nR_2        = Ri_2.mean()\nresidu_2   = sum((Ri_2 - R_2)**2)\nprint(center_2)\nxC = np.linspace((np.min(x)-0.1*R_2), (np.max(x)+0.1*R_2), 100)\nyC = np.linspace((np.min(y)-0.1*R_2), (np.max(y)+0.1*R_2), 100)\nX, Y = np.meshgrid(xC,yC)\nF = (X-xc_2)**2 + (Y-yc_2)**2 - R_2**2\nplt.plot(x, y, 'ok')\nplt.show()","execution_count":30,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"39bdb0e52705e09984243021a22ba3fa86d0bad4"},"cell_type":"code","source":"hits, cells, particles, truth = load_event(getPath(cEventDataDir,cFirstEvent))\nhits.head()","execution_count":31,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"96192b4bb15714da05c885fd9cd422557a734eab"},"cell_type":"code","source":"hits.describe()","execution_count":32,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ebed105b4a77e485da4e01d0e7cce4e03919e262"},"cell_type":"code","source":"cells.head()","execution_count":33,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3944df676773faa7f3a880e062f2d7adabce5b4a"},"cell_type":"code","source":"cells.describe()","execution_count":34,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"843d4c1c3a68e6c8b6d7fad6fb2c356e1d94c624"},"cell_type":"code","source":"particles[(particles['q'] != -1) & (particles['q'] != 1)]","execution_count":35,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a52e9e6441cd44f100eef802fc12e580ebf5e400"},"cell_type":"code","source":"truth.head()","execution_count":36,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ee7264c95e09f4ba78c51797a0eec0abb20b3083"},"cell_type":"code","source":"truth.describe()","execution_count":37,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c9ef5c46dd46b99428f76bee704723887fd0eda9"},"cell_type":"code","source":"track = truth[truth['particle_id'] == 4503737066323968]\nfig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\nfor hit_id, hit in track.iterrows():\n    ax.scatter(hit.tx, hit.ty, hit.tz)","execution_count":38,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"7a30190d4cebe8d8487716f4ee219dee5658d3e3"},"cell_type":"code","source":"def calc_curvature(data_fr):\n    x = data_fr.tx\n    y = data_fr.ty\n    z = data_fr.tz\n    ddx  = np.diff(np.diff(x))\n    ddy  = np.diff(np.diff(y))\n    ddz  = np.diff(np.diff(z))\n#     take the mean curvature (not the sum) to avoid bias \n#     since some particles generate more hits and others less\n    return np.sqrt(ddx**2 + ddy**2 + ddz**2).mean() ","execution_count":39,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a8179fa00c297bbe0fb47076356e43d2e08fe7b0"},"cell_type":"code","source":"df  = pd.merge(hits_df,truth_df,how = 'left', on = ['hit_id','event_name'])\ndf = df[df['particle_id']!= 0] # drop particle 0 \ngrouped = df.groupby(['event_name','particle_id'])\ncurvatures = grouped.apply(calc_curvature)","execution_count":42,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2c565c40204d5a84458d5cbb67a3a8b656b50403"},"cell_type":"code","source":"import seaborn as sns\ng = sns.jointplot(hits.x, hits.y,  s=1, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n    v = hits[hits.volume_id == volume]\n    plt.scatter(v.x, v.y, s=3, label='volume {}'.format(volume))\n\nplt.xlabel('X (mm)')\nplt.ylabel('Y (mm)')\nplt.legend()\nplt.show()","execution_count":57,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"da643d6c4752c5088a34645ac2a4eeafb6e6e9eb"},"cell_type":"code","source":"g = sns.jointplot(hits.z, hits.y, s=1, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n    v = hits[hits.volume_id == volume]\n    plt.scatter(v.z, v.y, s=3, label='volume {}'.format(volume))\n\nplt.xlabel('Z (mm)')\nplt.ylabel('Y (mm)')\nplt.legend()\nplt.show()","execution_count":58,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"952a23981f5a0d8d0101c53c2df1b631bfa34edb"},"cell_type":"code","source":"hits_sample = hits.sample(8000)\nsns.pairplot(hits_sample, hue='volume_id', size=8)\nplt.show()","execution_count":59,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}