{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"},"cell_type":"markdown","source":"# Motivation: Measurements in the Pixel detector\n\nThe innermost componets of tracking detectors (``volume_id = 7,8,9``)  for high energy physics are often Silicon Pixel detector, as it is is modelled in this dataset.  A pixelated readout structure allows to detect the *single* or *group of* pixels that are traversed by the particle and thus receive a signal. Evidently, this is strongly dependent on the incident angle of the particle into the detection sensor, as illustrated below:\n\n<img src=\"https://asalzbur.web.cern.ch/asalzbur/work/tml/PixelModule.png\" width=600>\n\nThe markers here are:\n- **black** dots : pixel center positions\n- **magenta** dot: reconstructed cluster position\n- **red** dot: true intersection of particle with sensor mid-surface\n\n\nEach single module has thus a channel system which is two-dimensional to the local coordinates. These are the channel indizes ``ch0`` and ``ch1``, which indicate which pixel(s) on the module have been crossed by a particle. When a particle crosses a pixel, it induces charge by *ionisation*, this does not alter the charge of the traversing particle, but reduces its energy (and thus momentum) slightly. The value of the charge can be read out in the pixel detector and is stored for each channel identifiec by ``(ch0,ch1)`` as the approprated ``value``. The following shows such a pixel module channel schema with the came cluster as above, the coloring of the traversed pixels indicates their ``value`` (charge).\n\n<img src=\"https://asalzbur.web.cern.ch/asalzbur/work/tml/PixelChannels.png\" width=600>\n\n## Position information\n\nIn a first stage, pixels that are adjunct are grouped together into *clusters*, an operation that has already been done in the presented dataset. The particle intersection with the module can then be rather precisely determined by taking the pixel position - and, as it it the case in the dataset - using the charge information of the individual pixels that contribute to the cluster. \n"},{"metadata":{"collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"60d62720b94bbf2c32713ad0295fe0e37db8ee4f","_cell_guid":"f8e31dda-67f7-4182-b576-f927b778374c"},"cell_type":"markdown","source":"Import the `trackml` library for convenient data loading."},{"metadata":{"collapsed":true,"_uuid":"6d305847511c75d68ff5b780eb3cc3d099f617ca","_cell_guid":"4a521e26-8939-4990-a183-1d279e37d121","trusted":false},"cell_type":"code","source":"import trackml\nfrom trackml.dataset import load_event, load_dataset\nfrom trackml.randomize import shuffle_hits\nfrom trackml.score import score_event","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"6b3a43c7a44658ed0c18b043fb17e9d69bd97289","_cell_guid":"e94f8a63-a64e-4909-a387-2b057948c136"},"cell_type":"markdown","source":"Load a specific event ..."},{"metadata":{"_uuid":"317df99e0eb4962ad97fe2eeac496361153e74df","_cell_guid":"a411352f-c3c9-4bcb-837e-64754f261735","trusted":false,"collapsed":true},"cell_type":"code","source":"# One event of 8850\nevent_id = 'event000001000'\n# \"All methods either take or return pandas.DataFrame objects\"\nhits, cells, particles, truth = load_event('../input/train_1/'+event_id)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"c6492d01f181a5e5e9471728ccd69aa699c0efb2","_cell_guid":"7f6ee0d3-4aa2-4e32-b8b2-60b52461a2af"},"cell_type":"markdown","source":"Let us now pick a single *Pixel cluster*, I've picked a rather large one ... "},{"metadata":{"collapsed":true,"_uuid":"89f7b9776a0a6ac3a5092bbec6941b0d1a3529c6","_cell_guid":"4fbff9d4-0e31-4170-9c62-b9354d8419a2","trusted":false},"cell_type":"code","source":"h_id = 19144 \npixel_cluster = cells[ cells['hit_id']==h_id ]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fe0b3c77399d5d398918959163ddc1baf808067d","scrolled":true,"_cell_guid":"700cc44d-bb80-4261-9b60-1ad1d69d3f06","trusted":false,"collapsed":true},"cell_type":"code","source":"len(pixel_cluster)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"380123ce576759cecf3e9f272c80b24338f0ac75","_cell_guid":"a6130c72-a186-4222-a18f-20a644a464eb"},"cell_type":"markdown","source":"This cluster has  1 to many individual pixels contributing, let's display the cluster "},{"metadata":{"collapsed":true,"_uuid":"f46db315307aaa04c458009b004e5a02006c4bb8","_cell_guid":"ece67463-a064-48fa-b05e-f3c5326ffcd1","trusted":false},"cell_type":"code","source":"# a function that calculates the cluster size and makes a pixel matrix\ndef pixel_matrix(pixel_cluster, show=False):\n    # cluster size\n    min0 = min(pixel_cluster['ch0'])\n    max0 = max(pixel_cluster['ch0'])\n    min1 = min(pixel_cluster['ch1'])\n    max1 = max(pixel_cluster['ch1'])\n    # the matrix\n    matrix = np.zeros(((max1-min1+3),(max0-min0+3)))\n    for pixel in pixel_cluster.values :\n        i0 = int(pixel[1]-min0+1)\n        i1 = int(pixel[2]-min1+1)\n        value = pixel[3]\n        matrix[i1][i0] = value \n    # return the matris\n    if show :\n        fig = plt.figure()\n        ax = fig.add_subplot(1,1,1)\n        ax.set_aspect('equal')\n        plt.imshow(matrix, interpolation='nearest', cmap=plt.cm.YlOrRd)\n        plt.colorbar()\n        plt.show()\n    return matrix, max0-min0+1, max1-min1+1","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c20d4316b6c8b8afc09aedfc438b55efa642e214","_cell_guid":"e4f2440f-caf4-445c-ba1e-2f69425f7b43","trusted":false,"collapsed":true},"cell_type":"code","source":"cluster,width,length = pixel_matrix(pixel_cluster,True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2d533dcc73863cfcad1763b76f54fe1051e7c4a7","_cell_guid":"53725a9e-6458-4527-a0cb-dbd918da3086"},"cell_type":"markdown","source":"The cluster size in `u` and `v` direction is:"},{"metadata":{"_uuid":"d5a99f96c3fe0d6e7e73bbc78844c326d8a50a14","_cell_guid":"c37ece43-f44d-42d7-b78f-f0e51fa4b9bc","trusted":false,"collapsed":true},"cell_type":"code","source":"print(width,length)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9d13165f92b6edb41531b39f443ba8e4d6240371","_cell_guid":"3f67d513-fd3c-4f1c-be00-ab1c615105a3"},"cell_type":"markdown","source":"You can almost *see* the particle's trajectory through the silicon, the coloring corresonds to the charge of the individual pixel. As you can see, the edge pixels which are not fully traversed by the particle have only little charge, which is what you expect, because the particle did traverse less of Silicon and thus induces less charge.\n\nLet us now see how well the position is estamated for this measurement:"},{"metadata":{"collapsed":true,"_uuid":"2527803ca75bed8dce9b43036de5cec599269ecb","_cell_guid":"a8d92365-3117-434f-8c07-2fb5f62d40ca","trusted":false},"cell_type":"code","source":"pixel_hit = hits[ hits['hit_id']==h_id ]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6cb6961e1324c1ea64af3edb2cc638e9afd90449","_cell_guid":"c7640e3f-70af-4e0f-87a3-7fae179e9908","trusted":false,"collapsed":true},"cell_type":"code","source":"print(pixel_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"76848bd1186c6e157f26988746d7fe51c758c050","_cell_guid":"efde8470-633a-4053-b25d-9f98f4364a4d"},"cell_type":"markdown","source":"In the truth file we can find the **truth** position of the intersection:"},{"metadata":{"collapsed":true,"_uuid":"1bebacb3855ccb71e722a8b6d4474a7cfa7dc024","_cell_guid":"98493e76-fbd0-4aa0-b15f-5fe7eb58857f","trusted":false},"cell_type":"code","source":"truth_pixel_hit = truth[ truth['hit_id']==h_id]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7917205905967ccde8d2c9e8e310d2b38b282589","_cell_guid":"b921082c-d410-481d-b429-97eab14475e1","trusted":false,"collapsed":true},"cell_type":"code","source":"print(truth_pixel_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"22c947a964296dc94e8bd1c60b9c30e2376068d9","_cell_guid":"d033321f-31b5-4be7-8dcf-e289de3bf68b"},"cell_type":"markdown","source":"The reconstruction did really well!! Compare ``(x,y,z)`` with ``(tx,ty,tz)`` , the values are **really** close, that's a good detector.\nHow did we come from the local cluster position on the surface, which was made of the local positions of the individual cells - to the global positions ? \n\nThis is with the help of the geometry description. \n\nEach module is unequile positioned in space through a **center** position and a **rotation** that transforms the local coordinate system ``(u,v,w)`` to the global coordinate system ``(x,y,z)``. For the dataset, ``ch0`` is measured in ``u`` and ``ch1`` is measured in the ``v`` direction, while the ``w`` direction is the ``thickness`` of the module as shown above.\n\n<img src=\"https://asalzbur.web.cern.ch/asalzbur/work/tml/localToGlobal.png\" width=600>\n\nIt's time to load the detector descrption. \nWe load the full detector and then picke the module our pixel cluster in question as measured at:"},{"metadata":{"collapsed":true,"_uuid":"b9e46fb3534227217acb8e0feb7f5032465e1909","_cell_guid":"775714aa-b594-4850-9c2a-67dd86eb61a9","trusted":false},"cell_type":"code","source":"detector = pd.read_csv('./../input/detectors.csv')\n# method to retrieve the according module associated to a hit\ndef retrieve_module(detector,hit) :\n    volume = detector[ detector['volume_id']==hit.volume_id.data[0] ]\n    layer  = volume[ volume['layer_id']==hit.layer_id.data[0] ]\n    module = layer[ layer['module_id']== hit.module_id.data[0] ]\n    return module\n# get the one for our example\nmodule = retrieve_module(detector,pixel_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"63e963cf88a6c342724d003b5e721b76818cb2aa","_cell_guid":"adc2cdb1-cb19-47ff-9ce9-6c88b37c28f0"},"cell_type":"markdown","source":"We now have the full detector information of this detection module:\n- its *translation* given by the module center position `(cx,cy,cz)`\n- its *rotation* with reference to the global coordinate system `(rot_xy, ..., rot_zw)`\n- the  module thicknes  `thickness = 2 * module_t`\n- the module dimension (rectangular): `2 * module_minh` in local `u` and `2 * mondule_hv` in locl `v`\n- the measurement segmentation, i.e. the pixel size `pitch_u` in `u` and `pitch_v` in `v`\n\n## Direction information\n\nIn the pixel cluster shape there's obviously some directional information decoded about the particle, to access this, we can compare the *expected* cluster size and the *measured* ones.\nTo calculate the expected cluster shape of a track hypothesis on a sensor, you need to express the *momentum direction* in the local coordinate system of the module. This is done by applying the inverse rotation to the *global direction*:\n\n`direction_local(_hit) = rotation.inverse() * direction_global_(hit)`\n\nWe will take the truth direction here to demonstrate - at the ``hit`` position, in global coordinates."},{"metadata":{"collapsed":true,"_uuid":"92ff6cf7a8da03ead5472b7f363dabe47cfe767f","_cell_guid":"8a8e0dcc-21e6-47b1-a42b-fcbfa6a90f4f","trusted":false},"cell_type":"code","source":"# method to build and nomralize direction vector from the \ndef direction_vector(ipx, ipy, ipz) :\n    # the absolute momentum for normalization\n    p = np.sqrt(ipx*ipx+ipy*ipy+ipz*ipz)\n    # build the direction vector - to be used with the matrix \n    direction = [[ipx/p], [ipy/p], [ipz/p]]\n    return direction","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1481cdb19404bf9dc5182c1ba7c09d473cd5854a","_cell_guid":"6d97e6b3-d510-469a-83f1-017eb452f232","trusted":false,"collapsed":true},"cell_type":"code","source":"# get the truth direction at the module, it's more accurate than the starting position\ndirection_global_hit = direction_vector(truth_pixel_hit.tpx.data[0],truth_pixel_hit.tpy.data[0],truth_pixel_hit.tpz.data[0])\nprint(direction_global_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f8830a4f98661c092c58a6511317e078d521820f","_cell_guid":"c2644862-bae4-4ecc-b448-e7a51cd013f1"},"cell_type":"markdown","source":"Let's do a comparison how the global momentum direction at particle creation:"},{"metadata":{"_uuid":"9af008add42d2e4d33631b7d7c9e17c8f3e388b3","_cell_guid":"c8972ac9-908a-49dd-b728-26a0ce132333","trusted":false,"collapsed":true},"cell_type":"code","source":"# get the truth particle information\nparticle = particles[ particles['particle_id'] == truth_pixel_hit.particle_id.data[0] ]\nprint(particle)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b8c0d4a855d6a43123f0f025fae7a459dfaabf72","_cell_guid":"72947206-e4b1-4de7-8b68-e8ebcd472149","trusted":false,"collapsed":true},"cell_type":"code","source":"# build the direction vector with the start momentum\ndirection_global_start = direction_vector(particle.px.data[0],particle.py.data[0],particle.pz.data[0])\nprint(direction_global_start)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f25d5fc8a366531427baa64fa1bd0c95085eb3a2","_cell_guid":"9fbbca67-95d4-4f82-9865-7afd8114d401"},"cell_type":"markdown","source":"In global coordinates, the *polar* angle ``theta_global(_hit)`` is (from truth information):"},{"metadata":{"_uuid":"99c5c7edb8168f87636409b64841d7bf7e7818ca","_cell_guid":"55155e0d-d358-4ebd-b860-cd7adc137d12","trusted":false,"collapsed":true},"cell_type":"code","source":"# extract phi and theta from a direciton vector\ndef phi_theta(dx,dy,dz) :\n    dr  = np.sqrt(dx*dx+dy*dy)\n    phi = np.arctan2(dy,dx)\n    theta = np.arctan2(dr,dz)\n    return phi, theta\n# get thet and phi\nphi_hit, theta_hit = phi_theta(direction_global_hit[0][0],\n                               direction_global_hit[1][0],\n                               direction_global_hit[2][0])\nprint(theta_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6231c26fa465d9c5a407741fd1cbe27f9fc9e5ee","_cell_guid":"2cfb1547-b449-4057-93e6-ff67a893383a"},"cell_type":"markdown","source":"From this, the expected cluster length is given from trigenometry, we using the module thicknes ( = `2*module_t`) to calucate the path in the Silicon wafer.\nWhen staying in the global system, we can only assume that the module is parallel to the ``z``-axis, so we can only get an estimate. There's little we can say about the `u` direction, if we do not know how the moudle is oriented."},{"metadata":{"_uuid":"c223ff43e7314731cdc14ef7435fb51780708eba","_cell_guid":"9e54fdb5-5c96-40b0-a4a5-c8a59cbb958e","trusted":false,"collapsed":true},"cell_type":"code","source":"# get the length of the cluster in v direction\ncluster_length_v_hit = np.abs(2.*module.module_t.data[0]/np.tan(theta_hit))\ncluster_size_v_hit   = cluster_length_v_hit/module.pitch_v.data[0]\nprint(cluster_length_v_hit,cluster_size_v_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b84ea4f0768d8cfdf6eb2337d3d306f418b3d927","_cell_guid":"271909f1-ae26-4f92-9005-ccd07979000a"},"cell_type":"markdown","source":"Not bad, we estimated the cluster size rather accurately only from the global information of the particle and the knowledge that the detector element is from the barrel detector (this, you can see from the `volume_id`).\n\n\n### Refining with geometry information\n\nSo far, we have not used the `module` information in order to get the actual incident angle into the module (which is the direct cause of the cluster size), hence this information should help to increase the accuracy of our prediction."},{"metadata":{"collapsed":true,"_uuid":"1744a32027591bb6d73abcd7b419c286bef97758","_cell_guid":"20b0e8b7-7783-497a-b5d1-ea3da1517acb","trusted":false},"cell_type":"code","source":"# function to extract the rotation matrix (and its inverse) from module dataframe\ndef extract_rotation_matrix(module) :\n    rot_matrix = np.matrix( [[ module.rot_xu.data[0], module.rot_xv.data[0], module.rot_xw.data[0]],\n                            [  module.rot_yu.data[0], module.rot_yv.data[0], module.rot_yw.data[0]],\n                            [  module.rot_zu.data[0], module.rot_zv.data[0], module.rot_zw.data[0]]])\n    return rot_matrix, np.linalg.inv(rot_matrix)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d4d8e789f1e6c0adb4436fa38118729ec32c8646","_cell_guid":"88e4c99e-2f16-4f4e-a4c7-f302e1fa5fc4"},"cell_type":"markdown","source":"Let's extract the rotation matrix and it's inverse from the module then, and transform the global direction into a local direction:"},{"metadata":{"_uuid":"7ddf490ff869742f205e551a3b2a8ff9eb61510e","_cell_guid":"a497c7f9-94f8-4f91-bf1d-6d9574198d14","trusted":false,"collapsed":true},"cell_type":"code","source":"module_matrix, module_matrix_inv = extract_rotation_matrix(module)\nprint (module_matrix)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1bdd0897ec45b5eaf70df0b3a0e9645a2e40a724","_cell_guid":"83004ddf-43c6-49f4-a0be-1ad676c196cd"},"cell_type":"markdown","source":"And let's have a look at the inverse martrix as well:"},{"metadata":{"_uuid":"22b1bd91226aa1d95df471f71e0442968b29e965","_cell_guid":"a62401a3-ad05-4342-b27d-853f46db2785","trusted":false,"collapsed":true},"cell_type":"code","source":"direction_local_hit =  module_matrix_inv*direction_global_hit\nprint(direction_local_hit)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd88720f40f9f2248c7cd651b42051b5e75daa9d","_cell_guid":"23afd10c-9224-41dc-86a4-c10144a90e32"},"cell_type":"markdown","source":"As you can see, the local momentum direction has the biggest component along the local `v` axis, this is not surprising giving that the cluster size in `v` is about \nNow let's see, what *predicted cluster sizes we get*, first we need to calculate the `phi` and `theta` in the local coordinate frame:"},{"metadata":{"_uuid":"5f3f25807e8f620e894c2bfed31a02ce125c3ce7","_cell_guid":"7b9315e8-d2f1-48d1-8d80-42225983ace4","trusted":false,"collapsed":true},"cell_type":"code","source":"# theta is defined as the arctan of the radial vs the longitudinal components\n# phi is defined as the acran of the two transvese components\nphi_local,theta_local = phi_theta(direction_local_hit[0][0],\n                                  direction_local_hit[1][0],\n                                  direction_local_hit[2][0])\nprint(phi_local,theta_local)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f5c32d5fc3507aa6d87f173bae7cf20e9526284f","_cell_guid":"4b9144a6-3e6f-4824-a9ce-6815e439943a"},"cell_type":"markdown","source":"From the module thickness, we can get the full path length in the silicon:"},{"metadata":{"_uuid":"d64c17f215f33fa30f9382db421e5150a96e2fdd","_cell_guid":"f96f4eaa-7c13-4513-802d-e87c87d28901","trusted":false,"collapsed":true},"cell_type":"code","source":"path_in_silicon = 2*module.module_t.data[0]/np.cos(theta_local)\nprint(path_in_silicon)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1eebf2fff6d006d8ac8f73f88c5e94a53cb44a5b","_cell_guid":"2735d942-1a10-4876-bece-90f8e2698ffc"},"cell_type":"markdown","source":"Which finally allows us to calculate the path length in `u` and `v` and thus the cluster sizes:\n"},{"metadata":{"_uuid":"872abe3d08440b1b13588d93245e68930d634d0f","_cell_guid":"0a3fcdfe-a697-4fb4-9372-66a6f1085ed6","trusted":false,"collapsed":true},"cell_type":"code","source":"# calculate the component in u and v\npath_component_u = path_in_silicon*np.sin(theta_local)*np.cos(phi_local)\npath_component_v = path_in_silicon*np.sin(theta_local)*np.sin(phi_local)\ncluster_size_in_u = path_component_u/module.pitch_u.data[0]\ncluster_size_in_v = path_component_v/module.pitch_v.data[0]\n# print the cluster size \nprint(cluster_size_in_u, cluster_size_in_v)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"2e245437a7c60b67f03dc5a68151116a5074ebcb","_cell_guid":"bed7fe67-e60c-438d-8423-e2e36f24291a"},"cell_type":"markdown","source":"\nWe have `reconstructed` the cluster size quite accurately in `v` and somewhat okish in `u` - why is this ?\nThe cluster is rather long in `v` and thus we do not really converned if the particle entered really at the outermost extend in the `v` direction, while for the `u` direction this makes some difference, indeed. \n"},{"metadata":{"collapsed":true,"_uuid":"f1bc1772dc2ead8ef7257d26332f346627d4dc52","_cell_guid":"68114e15-efd2-4a8e-9bf3-971b291a311c","trusted":false},"cell_type":"code","source":"","execution_count":null,"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}