{"cells":[{"metadata":{"_cell_guid":"1f10c4b1-1d72-4714-8a65-6f393ae2a98e","_uuid":"ab05eab5e414bdc025ee4bad6f6625ea5986d012"},"cell_type":"markdown","source":"In this notebook, I'm looking at a few things in the small training set to try to understand what is in the data.\n\nIf you want to understand more about how detectors work and how particles interact with them, I would strongly suggest at least going over the relevant review articles from the Particle Data Group (PDG), which are all available for free online. There are also several textbooks on detectors available if you're really interested in knowing more details. Unfortunately, there aren't so many easily accessible resources for learning in detail how things like track reconstruction (the topic of this competition) are implemented in practice."},{"metadata":{"_cell_guid":"fb798923-26a9-46d2-9d02-a74197110110","_uuid":"a59eb58b2d38181cf7c55420e994a4556d378c1d","collapsed":true,"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a7dbda84-7ffd-4def-bdf4-a92821ff67ae","_uuid":"e438abe4edc178e364d3a0525c38c67d935cee4a"},"cell_type":"markdown","source":"First, let's read in one of the events."},{"metadata":{"_cell_guid":"126aa4b5-1750-4014-8a28-2b80c60865e2","_uuid":"4814881776b9b7d6592915fb597e1db7c339f38b","collapsed":true,"trusted":true},"cell_type":"code","source":"def get_event_data(path='../input/train_1/event000001000'):\n    cells = pd.read_csv(path+'-cells.csv')\n    hits = pd.read_csv(path+'-hits.csv',index_col=0)\n    particles = pd.read_csv(path+'-particles.csv',index_col=0)\n    truth = pd.read_csv(path+'-truth.csv',index_col=0)\n    return cells, hits, particles, truth","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f3d59bbc-cb2e-4c16-9289-213a960feb20","_uuid":"108d5864bbec20875dd1a9e4b7913fad8ada9b63","collapsed":true,"trusted":true},"cell_type":"code","source":"cells, hits, particles, truth = get_event_data()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b257d2e5-6687-4e4d-b6a3-fc0ff6384e76","_uuid":"e287752ecca5d6d347f1877a96d0f4982babfaf8"},"cell_type":"markdown","source":"Now, let's print the head of each of the dataframes."},{"metadata":{"_cell_guid":"0cae66cf-4342-464f-9089-d3c5764c05cf","_uuid":"7b1c6bdd5ece83d89da96464e3e9d6773011eee8","collapsed":true,"trusted":true},"cell_type":"code","source":"truth.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"21f1d500-d4f9-49a2-a5f6-8bd4dc0bb82b","_uuid":"95ce31149966eecc23e4168a7485835dbd9992ea"},"cell_type":"markdown","source":"We see that the particle IDs are in the form of a very long integer. My guess is that the number encodes some information about the particle but it may also just be some sort of randomized or hashed number to keep participants from learning any information from that. \n\nA couple of the first hits have no corresponding particle ID. The momentum suggests that these would have PeV-scale energies, which is obviously too high to be created at the LHC. These hits have weight 0, so they should be excluded from the data."},{"metadata":{"_cell_guid":"cd9c38e4-8f8b-4f9c-ab3c-b22a9bd7f4ea","_uuid":"6bd945d93336f0806a7fc53552bd288c07bbdd17","collapsed":true,"trusted":true},"cell_type":"code","source":"particles.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"fa40fd37-acb6-4311-b2f2-4e6152648f2c","_uuid":"63d57ccd80a256dc6676d019533c3cdf1ca99238"},"cell_type":"markdown","source":"This looks more reasonable."},{"metadata":{"_cell_guid":"d5e93bce-9285-4a5e-ad89-7306786cd24e","_uuid":"95f9f87ea6b23856b17db6185fb4f637e855a6a1","collapsed":true,"trusted":true},"cell_type":"code","source":"hits.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"5395c701-00bb-4a34-831c-5204f45b2d6e","_uuid":"fc4e63bfe688626d99359b5be250e0615c67b4e5"},"cell_type":"markdown","source":"This also looks fine. It gives us some mappings from the hit IDs to positions and detector IDs."},{"metadata":{"_cell_guid":"98773908-6b27-49ba-aa4b-f68c4d9601cc","_uuid":"6f8c7d87201966f08460d0ee5fe92196e9c76c3c","collapsed":true,"trusted":true},"cell_type":"code","source":"cells.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"4fd9087a-cad4-4854-b409-bba2834d613b","_uuid":"2927b4b507813092f89aa7dec1141ba5932b185c"},"cell_type":"markdown","source":"And this gives us channel and energy loss information. A single hit can cover multiple channels. We might want to try to study things like the resolution for different types of channels. It's also not clear exactly how to use the signal. Do we add all the signals for a given hit?"},{"metadata":{"_cell_guid":"651aa754-570a-4756-88b2-319abd1d6867","_uuid":"44742b5ca8873c9ecbe5c103fb550901ef2323b0"},"cell_type":"markdown","source":"## Understanding the Geometry\n\nFirst, we want to understand the geometry. We do have a geometry file available, but we can also plot out hte hits. First, let's look at the 2D projections for everything. The volume ID provides some information about what piece of the detctor we hit."},{"metadata":{"_cell_guid":"851c54ed-5755-4102-a172-f498f1cc79b3","_uuid":"66db9737578efd95a487086b597a152a9fdc58f8","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nplt.scatter(hits.x,hits.y,marker='.',alpha=0.1,c=hits.volume_id)\nplt.xlabel('x')\nplt.ylabel('y')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b986bb3b-bb64-478e-af0e-3c98d94f3f09","_uuid":"511a5e3dfc3b89fb97ba77bfaa84aeb733eef3fc"},"cell_type":"markdown","source":"Here, we see from the colors that there are 3 different volumes included here. We also see that there are a series of concentric rings with some random scattered hits in between the rings. There are gaps between the different volumes. The concentric-ring type geometry is very common in collider physics, as it matches the cylindrical symmetry of collisions well."},{"metadata":{"_cell_guid":"98b679c5-50a6-422c-af0c-31f9a2be6013","_uuid":"277748c3f567e6d9ba8fb0c78bd8162db949dc59","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nplt.scatter(hits.x,hits.z,marker='.',alpha=0.1,c=hits.volume_id)\nplt.xlabel('x')\nplt.ylabel('z')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"be3595d9-8038-43ca-bce0-da6bc5e4c816","_uuid":"e7ebda968a8c902acb1595fab270b3acc73a5ffe"},"cell_type":"markdown","source":"The $xz$ projection shows a bit more about what's going on. The random scattering between the rings in the $xy$ projection are from vertically-aligned layers at both ends of the tracker. (Basically, endcap layers.)\n\nWe can also clearly see the distinction between the inner tracker and the outer tracker. I haven't checked this but it's likely that much of the outer tracker is made of strip detectors while the inner tracker may be made only of pixels.\n\nIn this projection, the beam is going in the vertical dimension. The region around $z=0$ consists of concentric circles around the beam pipe.\n\nNow we can start zooming in. First, remove the endcap layers."},{"metadata":{"_cell_guid":"7ec3700d-f5f0-4334-a9c6-39ff53f2b72e","_uuid":"75ea3b9c6e455ce9dfee864910ac8ca21687c50b","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nplt.scatter(hits[np.abs(hits.z)<1000].x,hits[np.abs(hits.z)<1000].y,marker='.',alpha=0.1)\nplt.xlabel('x')\nplt.ylabel('y')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"7cfd9430-032f-477b-9225-11cbcd29fc77","_uuid":"ce7f5ae8d906790c7a5f2dc2b9671d9ce60050ef"},"cell_type":"markdown","source":"As we might expect, we now see the outer tracker as just a series of concentric rings around the beam pipe."},{"metadata":{"_cell_guid":"f2ad643f-0b68-4169-a121-40d666049ed6","_uuid":"a7ce3907f4c81efcfcd137aa7a802c65f6aa8e1b","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nhits_inner = hits[(np.hypot(hits.x,hits.y)<200) & (np.abs(hits.z)<1500)]\nplt.scatter(hits_inner.x,hits_inner.y,marker='.',alpha=0.1,c=hits_inner.volume_id)\nplt.xlabel('x')\nplt.ylabel('y')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"8acb5128-8410-4207-a6f2-7b88ec3d49a4","_uuid":"4e61f0a08aa8c26a66a754dfdf89094c52189746"},"cell_type":"markdown","source":"Looking closer at the inner tracker, we see that the innermost layer is only around 2.5 cm from the center. Having detectors so close to the beam is critical to many measurements, as that is what allows for good vertex resolution.\n\nSome interesting types of events, such as tau lepton and b quark production events can be identified because the tau or b decays a measurable distance from the original collision vertex. This can only be seen with a very good vertex position resolution.\n\nOne unfortunate effect of this is that the harsh environment slowly kills the detectors as they succumb to radiation damage. Modules in trackers typically need to be periodically replaced or there will end up with dead regions of the detector."},{"metadata":{"_cell_guid":"15983c80-c23e-4d05-9ed4-76974fd01b83","_uuid":"a188f44e4e03ef9989603b7976266e0b3e5f8827","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nplt.scatter(hits_inner.x,hits_inner.z,marker='.',alpha=0.1,c=hits_inner.volume_id)\nplt.xlabel('x')\nplt.ylabel('z')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"8db7ecfd-c811-4f0e-a76d-563d92911089","_uuid":"d01d5a6be6d0ad25fa720fc407c25938e0452946"},"cell_type":"markdown","source":"The inner tracker looks a lot like the outer tracker (but smaller obviously). There are a series of concentric rings right near the origin and then a series of endcap layers for very forward/backward-going tracks."},{"metadata":{"_cell_guid":"9b38d444-7838-4916-b11d-623f71f93e4f","_uuid":"e4aaf217243a407f37eecf170f132e973e3a8962","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,10))\nhits_inner2 = hits[(np.hypot(hits.x,hits.y)<200) & (np.abs(hits.z)<500)]\nplt.scatter(hits_inner2.x,hits_inner2.y,marker='.',alpha=1,c=hits_inner2.layer_id)\nplt.xlabel('x')\nplt.ylabel('y')\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"4f1794a1-0bb2-46bf-98fc-dd02c294f54d","_uuid":"727fb49da1ed9c730b5beaf31048dce9dbdde6a1"},"cell_type":"markdown","source":"Now looking only at the concentric rings, I've colored things by the layer ID. It is fairly simple to select hits only in certain layers using the volume and layer ID fields.\n\nWhen zoomed in this far, we also see something interesting. While previously, it might have looked like there was some aliasing in the images, here we clearly see that the concentric circles are not really circles. Instead, they are build from a series of planes. The actual silicon detectors are small flat rectangles, so this shows how the rings are created."},{"metadata":{"_cell_guid":"8f592c89-a106-49ca-a592-e5137794bb05","_uuid":"934c0d2620c70b9f20f6ee7011dec1ed0374e44d"},"cell_type":"markdown","source":"## Energy Loss Distribution\n\nNow, let's select all the hits from a particular layer. "},{"metadata":{"_cell_guid":"45104195-17b8-4f45-906a-57b5e1cd80ef","_uuid":"ddcf75185f5c46a525db32278a696f291b5027f6","collapsed":true,"trusted":true},"cell_type":"code","source":"hits_layer3 = hits[(np.abs(hits.z)<500)&(np.hypot(hits.x,hits.y)<190)&(np.hypot(hits.x,hits.y)>150)]\nhits_layer3.describe()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"96c4abe5-86cd-4637-9696-efb4e321f74c","_uuid":"2ce80881f2c449c98ddc2b652d602572cc349751"},"cell_type":"markdown","source":"This is the outermost concentric layer of the inner tracker, which has volume_id=8 and layer_id=8. \n\nI can merge the hits with the cells data to bring in signal information. Silicon detectors measure induced currents from particle-hole pairs created by ionization in the depletion region of the detector (typically a pn junction). \n\nOf course, there is a lot more that goes into taking a raw signal and getting something useful. The total energy lost is miniscule, and what I really might care about is how much energy is lost per amount of stuff the particle traveled through. To do this, I should add in the particle momentum information and use the angle of incidence and the detector thickness to normalize the results.\n\nThe dataset also does not include truth information about the particle type beyond the charge, so particle identification (PID) might not be feasible here."},{"metadata":{"_cell_guid":"cebc7270-0c18-4c54-aa20-fbf57f3fcd6e","_uuid":"7197c0b344b2cae22e6f1cd2b93ead592af729c5","collapsed":true,"trusted":true},"cell_type":"code","source":"hits_with_energy = cells.merge(hits_layer3,how='inner',left_on='hit_id',right_index=True)\nfig = plt.figure(figsize=(12,6))\nax = fig.add_subplot(121)\nax.hist(hits_with_energy.value,bins=100)\nax.set_xlabel('Cell Signal Size')\nax = fig.add_subplot(122)\nax.hist(hits_with_energy.groupby('hit_id').value.sum(),bins=100,range=(0,1.2))\nax.set_xlabel('Summed Hit Signal')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d32b1e1d-d058-4dfb-9715-b19e99c00780","_uuid":"25adffbcdb0cfed82b0087d702bced018ba66644"},"cell_type":"markdown","source":"This is quite interesting. Typically, as relativistic particles travel through thin detectors, the energy loss due to ionization follows a Landau distribution (think a Gaussian with a very long positive tail).\n\nHere, if we look at the individual cell signals, we see a nice peak with a long positive tail, but we also see a shoulder at low energy loss.\n\nIf we add all the cells together, we lose the shoulder, but also see the signal distribution mostly cut off around 1. Maybe this is related to how cells are combined into hits."},{"metadata":{"_cell_guid":"69dd0a98-1ed8-4476-87f0-91c1affba080","_uuid":"9cf3c85ec06c13d12c5007c543953397e050124a","collapsed":true,"trusted":true},"cell_type":"code","source":"hits_with_energy = cells.merge(hits[(hits.volume_id==8)&(hits.layer_id==6)],how='inner',left_on='hit_id',right_index=True)\nfig = plt.figure(figsize=(12,6))\nax = fig.add_subplot(121)\nax.hist(hits_with_energy.value,bins=100,range=(0,0.25))\nax.set_xlabel('Cell Signal Size')\n\nax = fig.add_subplot(122)\nax.hist(hits_with_energy.groupby('hit_id').value.sum(),bins=100,range=(0,1.4))\nax.set_xlabel('Summed Hit Signal')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1aabf602-6134-46d4-97bb-0087b825dd0b","_uuid":"9218892a478fe23acdf77d291a2e680e978be679"},"cell_type":"markdown","source":"We see something similar if we go inward one layer, but now the summed hit signal distribution extends quite a bit farther. This may just be an effect from (1) maybe using different detector or readout hardware and (2) having a higher particle density."},{"metadata":{"_cell_guid":"60383638-9100-4bf5-91ee-1b2223e67456","_uuid":"d6ea9b9d03d9551c8283a90ecbb1603062232a51","collapsed":true,"trusted":true},"cell_type":"code","source":"hits_with_energy = cells.merge(hits[(hits.volume_id==8)&(hits.layer_id==4)],how='inner',left_on='hit_id',right_index=True)\nfig = plt.figure(figsize=(12,6))\nax = fig.add_subplot(121)\nax.hist(hits_with_energy.value,bins=100,range=(0,0.25))\nax = fig.add_subplot(122)\nax.hist(hits_with_energy.groupby('hit_id').value.sum(),bins=100,range=(0,2))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b5c7d6ae-74fc-43a0-859d-e575fc51324e","_uuid":"d71257b22777439fbe2b7aca269d810734057f44"},"cell_type":"markdown","source":"And we see the same thing as we continue to move closer to the beam."},{"metadata":{"_cell_guid":"d3f84059-2ea0-4ac3-b123-18985770c8c8","_uuid":"ab9a302cb8977bb13b4b6e71393f42f45376e20d"},"cell_type":"markdown","source":"## Particle Truth Information\n\nThe particles file includes truth information from the particles when they are created. We don't have much beyond basic kinematics and a charge. First, let's add in some common kinematic variables used in high energy physics. The are:\n\n  - $p_{t}$: The transverse momentum (i.e. momentum perpendicular to the beam)\n  - $E$: Energy. Technically, this is just the total momentum since we don't have the mass, but in the relativistic limit, energy and the magnitude of the momentum are identical.\n  - $\\eta$: The pseudorapidity. This is an alternative measure of the direction with respect to the beam. High rapidity means very close to the beam direction\n  - $\\phi$: Azimuthal angle. Things are very close to cylindrically symmetric here, so we expect that this should not matter when large amounts of data are aggregated. If the data show significant dependence on this, then we'll need to understand why."},{"metadata":{"_cell_guid":"57ce4949-90f0-491c-9de1-f7b3891d8e2f","_uuid":"7bd20a2086874cfc1434e61470287d22b4af0b94","collapsed":true,"trusted":true},"cell_type":"code","source":"def add_kinematics(df):\n    df['pt'] = np.hypot(df.px,df.py)\n    df['E'] = np.hypot(df.pt,df.pz)\n    df['eta'] = -0.5 * np.log( (df.E + df.pz)/(df.E-df.pz))\n    df['phi'] = np.arctan2(df.vy,df.vx)\nadd_kinematics(particles)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2a81618b-f86b-4b9c-b759-e08e43ad1bd0","_uuid":"45fc75886542d54ba3a44eb34a32909c31833704","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.hist(particles.nhits,bins=20,range=(0,20))\nplt.xlabel('Number of hits')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"203ce59b-9fbb-4fa4-8973-646a6bcf94a1","_uuid":"47628333ff6816ea31dd491b373502d7d9a69fee"},"cell_type":"markdown","source":"Particles mostly either have very few or around 12 hits. The number 12 is probably just related to how the layers are set up in the detector. Low momentum or very forward particles will tend to miss everything."},{"metadata":{"_cell_guid":"408dc1de-650a-4be7-aac3-cd793aca69c0","_uuid":"bec3e00cc48dfebd455f87501ce6a52b9f398d82","collapsed":true,"trusted":true},"cell_type":"code","source":"particles.q.value_counts()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a726f8b4-218c-420b-8a9a-90b5774591ff","_uuid":"d28fd8005fc24d0c58e21f1ad7d5418e51ad4ae9"},"cell_type":"markdown","source":"There are substantially more positively-charged particles than negatively charged ones. I guess we expect somewhat more positively charged particles since the LHC runs proton-proton collisions rather than proton-antiproton like at the Tevatron or the S$\\rm p\\bar{p}$S collider. I'm not sure if we expect this much of a difference though."},{"metadata":{"_cell_guid":"9229f3e1-b81a-48d8-88f5-1a7bd76b0b3a","_uuid":"90e9b837b2f1b672113fabc2cbf4c12cbaf108f7","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.hist(particles.pt,bins=100,range=(0,10))\nplt.xlabel('Transverse Momentum [GeV]')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"9a1c872c-77b5-433e-ba58-d70956ab0000","_uuid":"887a7821601f88f4cdd1813fc27b4fcafd6ae9d1"},"cell_type":"markdown","source":"The transverse momentum distribution looks pretty exponential, with few particles having $p_{t}>2$ GeV. One thing we can try to look at here is the momentum balance. We can sort by transverse momentum first and then check the balance by looking at the \"missing $p_{t}$\" for different subsets of particles.\n\nNote that this won't really work here since I haven't tried to work out particle ancestry. If one high-$p_{t}$ particle is the parent of another (due to something like a decay process), I could be double counting the momentum.\n\nBut if you want to figure out things like when decays and inelastic interactions happen, or at least select only things originating right near the beam, maybe you will see something. This is complicated by the fact that there are going  to be many interactions per collision and that most collisions just lead to a bunch of junk being sprayed everywhere. You might miss a lot of momentum due to particles staying in the beam pipe. Ideally, you can find the particles associated with each interaction vertex and then look at the momentum balance for each vertex that looks like it might be intersting."},{"metadata":{"_cell_guid":"d9f9acae-bce7-41da-b863-23a87d9c3e8c","_uuid":"25ebbc7d48727ccf4e04a593cb7b78847f56b64a","collapsed":true,"trusted":true},"cell_type":"code","source":"particles_pt_sort = particles.sort_values(by='pt')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1b2d2433-82ac-4b3e-a53a-bd8943cfded3","_uuid":"d8ebff4eeb7d06f0a2a4acd00b4d3bf1c17502f9","collapsed":true,"trusted":true},"cell_type":"code","source":"px_summed = particles_pt_sort.py.cumsum()\npy_summed =particles_pt_sort.px.cumsum()\npt_summed = np.hypot(px_summed,py_summed)\nprint(pt_summed.values[-20:])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c92e6184-57dc-49bf-add0-103eb4bb50e3","_uuid":"9032bccc3532ac43c859f2577c8ad2fe78313913"},"cell_type":"markdown","source":"This first set shows what happens if I look at the total momentum balance excluding some of the highest $p_{t}$ particles. We can see that the momentum balance is quite far off but actually is reduced significantly by those particles. But again, I really need to do this in a more careful way to really say anything about this."},{"metadata":{"_cell_guid":"6a4de40f-d087-4ad9-9c4c-da98a6fb19d7","_uuid":"a361e799db2674dea2bdc9e3928a8a63015ea937","collapsed":true,"trusted":true},"cell_type":"code","source":"particles_pt_sort = particles.sort_values(by='pt',ascending=False)\npx_summed = particles_pt_sort.py.cumsum()\npy_summed =particles_pt_sort.px.cumsum()\npt_summed = np.hypot(px_summed,py_summed)\nprint(pt_summed.values[:20])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ff9ec495-0ea5-4f90-88c3-dfd0d1e0cbad","_uuid":"0745893dd9bd80b66b16ba1882845ae39df9fc40"},"cell_type":"markdown","source":"And this shows the momentum balance for only the few highest $p_{t}$ particles. Again, it looks pretty high but I am summing up everything from a number of vertices and may be double-counting momentum from decay and inelastic interaction events."},{"metadata":{"_cell_guid":"29573d80-5bc5-40c7-9a8c-ebf056b2c573","_uuid":"c67ad3e9e126e3818b3059dfb32f555825e5b809","collapsed":true,"trusted":true},"cell_type":"code","source":"particles_pt_sort.head(20)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"9e1fd0bb-dc08-453b-ac73-54c406f707cf","_uuid":"4c8760a099413dcb25b670505e1ebfe4e4e5b4b9"},"cell_type":"markdown","source":"We can see that there are some very high transverse momentum particles, but many of them also have a large longitudinal momentum. The kinds of things that are really interesting to analyses for things like electroweak physics are particles with high transverse momentum and not much longitudinal momentum. A low pseudorapidity and high energy means that a lot of momentum had to have been transferred somewhere."},{"metadata":{"_cell_guid":"aca6d7f2-18ee-4502-afac-3033c8637f24","_uuid":"fc2b9cbe6fbcc570f61adc5e0a95277d080f844e"},"cell_type":"markdown","source":"## Particles Near/In the Beam\n\nNow let's take particles with fairly high $p_t$ (1 GeV or more) and zoom in on the beam."},{"metadata":{"_cell_guid":"78a666a4-de5d-475f-8e12-04b74efb5549","_uuid":"35d130057bdbd9fa409320d504401fb60c8718fa","collapsed":true,"trusted":true},"cell_type":"code","source":"particles_high_pt = particles[particles.pt>1]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"72daa878-9c3e-4062-963f-921f4948d5e2","_uuid":"dc47bab2ba9240e1407bec4037f5a4d8a39fa2fa","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.scatter(particles_high_pt.vx,particles_high_pt.vy,alpha=0.3)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"dfdcc7e6-5d9e-476c-bacd-b79cf6687898","_uuid":"4ec1c326711eb83f8b8b4180239a853ed0094cdf","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.scatter(particles_high_pt.vx,particles_high_pt.vy,alpha=0.3)\nplt.xlim([-50,50])\nplt.ylim([-50,50])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3705695d-2eb3-405b-bfaa-694751357b52","_uuid":"63875085690e0c6d6e1d42e6e8653742efb69c78","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.scatter(particles_high_pt.vz,particles_high_pt.vx,alpha=0.3)\nplt.ylim([-5,5])\nplt.xlim([-25,25])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2052933d-21b2-47a8-80d5-fd12d633750a","_uuid":"a629f4a4ab3e640ea93dd27b5b7aadec692e6e73"},"cell_type":"markdown","source":"I didn't put axis labels on everything but the distance scales are in millimeters. The collisions here are all happening within around 1 cm of the origin in the $z$ (beam) direction and a fraction of a cm in the transverse direction.\n\n## Invariant Mass\n\nOne way to look for decays of interesting particles is to calculate the invariant mass of sets of particles that you're interested in. Often, a particle will appear as a peak in a invariant mass distribution. For example, if you select a bunch of $Z\\rightarrow\\mu^+\\mu^-$ decays, you should see a nice peak centered at 91.2 GeV (well, it can get a bit more complicated than that due to various effects). The invariant mass is nice because it corrects for things like boosted frames, but it can be hard to interpret if you don't reconstruct all the particles in the final state. Here, we also don't know how particles are selected in the simulation, so we could accidentally look at too many particles. \n\nFor now, let's look at everything with $p_{t}>1$ GeV."},{"metadata":{"_cell_guid":"520826fd-890a-490b-9705-46b1f6b75e1c","_uuid":"0887ef30f4f22635adf55a64571259c205e5b642","collapsed":true,"trusted":true},"cell_type":"code","source":"invariant_mass = []\nfor ids,gp in particles[particles.pt>1].groupby(['vx','vy','vz']):\n    if gp.vx.count()==1:\n        continue\n    px_tot = gp.px.sum()\n    py_tot = gp.py.sum()\n    pz_tot = gp.pz.sum()\n    E_tot = gp.E.sum()\n    invariant_mass.append( np.sqrt(E_tot*E_tot-px_tot*px_tot-py_tot*py_tot-pz_tot*pz_tot) )\n    \n","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"002a4ae7-bbc3-4a47-b724-e495a9cbf8b9","_uuid":"59a686a73f38445ecd4d4d5580e98b8e16a00100","collapsed":true,"trusted":true},"cell_type":"code","source":"plt.hist(invariant_mass,bins=100,range=(1,2000))\nplt.xlabel('Invariant Mass [GeV]')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1911f980-96cc-464e-8e51-eab44123eb9a","_uuid":"1937a0acde837ca414fc178d2b971f9c22fe4f81"},"cell_type":"markdown","source":"We don't see much, but that's not surprising. Most of the particles here are probably from uninteresting sprays of hadrons or are secondary particles. It's actually not going to be particularly common to see truly interesting events, and probably we'd need more than just the tracker to tell how interesting something is.\n\n## Back to particles near the beam\n"},{"metadata":{"_cell_guid":"f93d91c7-fdcb-4d52-b86f-bffab3aadb3d","_uuid":"85d143dcb8b4ee66605399f868955aa19d30c602","collapsed":true,"trusted":true},"cell_type":"code","source":"particles_near_center = particles[(np.hypot(particles['vx'],particles['vy'])<0.05)\n                                  &(np.abs(particles['vz'])<25)]\nplt.scatter(particles_near_center.vz,particles_near_center.vx,alpha=0.3)\nplt.ylim([-0.1,0.1])\nplt.xlim([-25,25])\nplt.show() ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"5bc4f431-12ed-4098-beb7-92ff40fc1a36","_uuid":"fe3ca5ce22118ded225649101f7b8782126b6660","collapsed":true,"trusted":true},"cell_type":"code","source":"\nplt.hist(particles_near_center.drop_duplicates(subset=['vx','vy','vz']).vz,bins=25)\nplt.title('Particle Position Near Beam')\nplt.xlabel('Z [mm]')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"df7251da-29ef-4589-92cb-68bdb7ec7e57","_uuid":"11fb4358350a00be5f643c6a3951c1a6f72402d3"},"cell_type":"markdown","source":"We see that if we select only things right near the beam, there is a distribution peaked near 0 but with a width of order 5 mm. It doesn't really look Gaussian, but we probably want to add in more data to get a nicer profile and to define a cleaner selection.\n\n## $\\eta-\\phi$ Distributions\n\nFinally, let's look at the direction distributions of particles originating near the beam"},{"metadata":{"_cell_guid":"34b3913e-0cfa-4b06-b868-ca7d572859df","_uuid":"58fc7b43115aa34bdc06b752b89698a8caaa618f","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\nplt.hist2d(particles_near_center['eta'],particles_near_center['phi'],bins=(50,50))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"7db122a2-5bb9-42ab-a53f-907f7ef909cc","_uuid":"ecfc61edcece4f2c21626e9b11da21ed2aacd42c"},"cell_type":"markdown","source":"We see a couple potentially worrying things here. (Again, it would be good to get a cleaner selection). First, there is a big peak in the distribution. With just one event, maybe this is due to some correlated particles (like a jet or shower) and low statistics.\n\nWe also see clear horizontal bands, indicating that some azimuthal angles seem to be favored."},{"metadata":{"_cell_guid":"b05a1743-4095-494e-ac86-3e141c73632b","_uuid":"41a87118c5e9e13f1d0e5c6304502e5b6f7d4485","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\n\nplt.hist2d(particles_near_center['eta'],particles_near_center['phi'],bins=(50,50),range=([-4,4],[1,3]))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"30730378-62d4-4189-9a13-fd30f9ea37b6","_uuid":"5df56425e16e6829214eae5f3ba72adefad60808","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\n\nparticles_tmp = particles_near_center[particles_near_center.E>1]\nplt.hist2d(particles_tmp['eta'],particles_tmp['phi'],bins=(50,50),range=([-4,4],[1,3]))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f4f3ac38-4de7-4207-bc11-cc3682408d6f","_uuid":"42cf79734029fb93bfdab4ae1963643661ce1588"},"cell_type":"markdown","source":"We see that these effects persits even as we zoom in or make some minimum cut on the transverse momentum. So, let's bring in more data."},{"metadata":{"_cell_guid":"0bdaafc2-eb8c-4b67-b869-3793f3c8fd9b","_uuid":"62315b42bf74ac234a7f3cb51d20f8b4c7925be1","collapsed":true,"trusted":true},"cell_type":"code","source":"def get_particle_data(path='../input/train_1/event0000010{:02}-particles.csv'):\n    df = pd.DataFrame()\n    for i in range(100):\n        data = pd.read_csv(path.format(i),index_col=0)\n        df = pd.concat([df,data])\n    return df\nparticle_all  = get_particle_data()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b6b7217c-7948-45e0-bf1d-8be774d1707d","_uuid":"8a6de7322c598ea4d256cc4910583bf19f92ad14","collapsed":true,"trusted":true},"cell_type":"code","source":"add_kinematics(particle_all)    ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e98c56e9-1a80-4d15-8093-0e7f81c8c688","_uuid":"fecac54a414404e561bd94c2f4cbbdfeb0b18e07","collapsed":true,"trusted":true},"cell_type":"code","source":"def near_center(df):\n    return df[(np.hypot(df['vx'],df['vy'])<0.05)\n               &(np.abs(df['vz'])<25)]\ncenter = near_center(particle_all)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d9c342c5-9cac-4615-b158-f3b1876a4697","_uuid":"f778151c8ab39d53ee0346f97dabe6db0f66e807","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\n\nplt.hist2d(center['eta'],center['phi'],bins=(50,50))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c3b26d88-f5ca-41f5-a4b3-7bb00d390180","_uuid":"b25028b4744841ce3375c8ea27681c8f1de49996"},"cell_type":"markdown","source":"The big peak has gone away clearly see that certain azimuthal angles still appear to be favored. So, it might be important to figure out if this is some weird effect from how I chose my particles (maybe double counting some things?) or if there is some reason for this."},{"metadata":{"_cell_guid":"0362531e-d8bb-490c-b54b-094edc28e0d3","_uuid":"3a4512bfbd8e9833044c55f07929a690ba35b001","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\n\nplt.hist2d(center[center.E>1]['eta'],center[center.E>1]['phi'],bins=(50,50))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c1cda1eb-18d4-41d7-a748-7f60ae7c228d","_uuid":"d9540e73dfb9bb5bb51f5a8f86d2bf90749d617c"},"cell_type":"markdown","source":"Things still persist as we select on high energy. High energies are dominated by forward-going tracks."},{"metadata":{"_cell_guid":"c22a5c2a-8f4f-4579-8a52-0ede3308a318","_uuid":"dfd7fcbc2334a769443fa6a91ef78106160fe200","collapsed":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8,6))\n\nplt.hist2d(center[center.E>5]['eta'],center[center.E>5]['phi'],bins=(50,50))\nplt.xlabel(r'$\\eta$')\nplt.ylabel(r'$\\phi$')\nplt.colorbar()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"bf45ffaa-b6f2-4379-8108-ea552ead3efd","_uuid":"b63bf4f30ff96be74ced3f219bb17928c8b49357"},"cell_type":"markdown","source":"And that's it for now. There's a lot more explore in the data before even attempting any model building. It would definitely be good to start building some diagnostic tools based on the truth information to identify various kinds of events. I also haven't even looked at trying to build tracks, but a proper model will have to somehow account for the fact that tracks mostly move in helical paths due to a magnetic field that is probably in the geometry and also the fact that multiple scattering can lead to significant changes in particle directions (especially for electrons/positrons). So, it won't be enough to just look for straight lines. Additionally, there are various processes that can create new particles far from the beam, and a model will need to try to identify these as well."}],"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}