SIREN_HNLAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "715b07eb-e2e4-4977-9572-810409037c00",
   "metadata": {},
   "outputs": [],
   "source": [
    "import awkward as ak\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "lumi=3000\n",
    "\n",
    "# Remove empty sublists within each list\n",
    "def clean_array(array):\n",
    "    return np.array(np.squeeze(ak.Array([[sublist for sublist in inner_list if len(sublist) > 0] for inner_list in array])))\n",
    "                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b0990db4-ff5d-437b-a7e1-4122febefa7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = ak.from_parquet(\"hnl_test2.parquet\")\n",
    "events_to_inject = len(data)\n",
    "in_fid = np.all(np.array(data.in_fiducial),axis=-1)\n",
    "muon_flag = np.abs(data.secondary_types) == 13\n",
    "hnl_flag = np.abs(data.primary_type) == 5914\n",
    "mu_vertex = clean_array(data.vertex[hnl_flag])\n",
    "muon_momenta = clean_array(data.secondary_momenta[muon_flag])\n",
    "muon_momentum = np.linalg.norm(muon_momenta[:,:,1:],axis=-1)\n",
    "mu_dir = muon_momenta[:,:,1:] / np.expand_dims(muon_momentum,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bc4c2b62-e211-47e9-ba8e-7b6ce0cddebe",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "dot() got an unexpected keyword argument 'axis'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[16], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmu_dir\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m:\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43mmu_dir\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m:\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mTypeError\u001b[0m: dot() got an unexpected keyword argument 'axis'"
     ]
    }
   ],
   "source": [
    "np.dot(mu_dir[:,0,:],mu_dir[:,1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "909542e1-3eab-4949-ac7d-6b0deea94e82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "139"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(data.weights[in_fid])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e2de7e62-37b3-4a98-be68-921285de08b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>[{event_weight: 1e-48, event_gen_time: 0.151, event_weight_time: 0.00248, ...},\n",
       " {event_weight: 2.36e-26, event_gen_time: 0.000716, ...},\n",
       " {event_weight: 1.33e-65, event_gen_time: 0.00056, event_weight_time: ..., ...},\n",
       " {event_weight: 2.36e-21, event_gen_time: 0.00055, event_weight_time: ..., ...},\n",
       " {event_weight: 3.92e-21, event_gen_time: 0.000528, ...},\n",
       " {event_weight: 2.01e-37, event_gen_time: 0.000534, ...},\n",
       " {event_weight: 1.21e-27, event_gen_time: 0.000517, ...},\n",
       " {event_weight: 1.17e-92, event_gen_time: 0.000523, ...},\n",
       " {event_weight: 1.34e-17, event_gen_time: 0.000512, ...},\n",
       " {event_weight: 5.73e-19, event_gen_time: 0.000557, ...},\n",
       " ...,\n",
       " {event_weight: 1.95e-23, event_gen_time: 0.000517, ...},\n",
       " {event_weight: 2.06e-84, event_gen_time: 0.000506, ...},\n",
       " {event_weight: 5.98e-28, event_gen_time: 0.000515, ...},\n",
       " {event_weight: 5.84e-34, event_gen_time: 0.000517, ...},\n",
       " {event_weight: 8.39e-23, event_gen_time: 0.00051, event_weight_time: ..., ...},\n",
       " {event_weight: 3.08e-17, event_gen_time: 0.000504, ...},\n",
       " {event_weight: 1.02e-95, event_gen_time: 0.000512, ...},\n",
       " {event_weight: nan, event_gen_time: 0.000514, event_weight_time: 0.00246, ...},\n",
       " {event_weight: 1.79e-28, event_gen_time: 0.000519, ...}]\n",
       "--------------------------------------------------------------------------------\n",
       "type: 1000 * {\n",
       "    event_weight: float64,\n",
       "    event_gen_time: float64,\n",
       "    event_weight_time: float64,\n",
       "    event_global_time: float64,\n",
       "    num_interactions: int64,\n",
       "    vertex: var * var * float64,\n",
       "    primary_initial_position: var * var * float64,\n",
       "    in_fiducial: var * bool,\n",
       "    primary_type: var * int64,\n",
       "    target_type: var * int64,\n",
       "    num_secondaries: var * int64,\n",
       "    secondary_types: var * var * int64,\n",
       "    primary_momentum: var * var * float64,\n",
       "    secondary_momenta: var * var * var * float64,\n",
       "    parent_idx: var * int64,\n",
       "    num_daughters: var * int64,\n",
       "    int_probs: var * float64,\n",
       "    PDG: float64,\n",
       "    bjorken_x: float64,\n",
       "    bjorken_y: float64,\n",
       "    energy: float64,\n",
       "    hPDG: float64,\n",
       "    thx: float64,\n",
       "    thy: float64,\n",
       "    wgt: float64,\n",
       "    weights: float64\n",
       "}</pre>"
      ],
      "text/plain": [
       "<Array [{event_weight: 1e-48, ...}, ..., {...}] type='1000 * {event_weight:...'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.se"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e4650203-acbf-481d-a79f-ad57a12fa72a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "KeyboardInterrupt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "data = ak.from_parquet(\"Data/SIREN/Output/UNDINE_LHCb_North/LHC13_SIBYLL_light_14_NC.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "4320051c-30bf-40d5-a676-fe72d03a3c1f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>[0.0697,\n",
       " 0.03,\n",
       " 0.117,\n",
       " 0.164,\n",
       " 0.0889,\n",
       " 0.0697,\n",
       " 0.0197,\n",
       " 0.0258,\n",
       " 0.0582,\n",
       " 0.295,\n",
       " ...,\n",
       " 0.421,\n",
       " 0.066,\n",
       " 0.226,\n",
       " 0.0741,\n",
       " 0.0142,\n",
       " 0.697,\n",
       " 0.616,\n",
       " 0.0928,\n",
       " 0.0212]\n",
       "----------------------\n",
       "type: 138398 * float64</pre>"
      ],
      "text/plain": [
       "<Array [0.0697, 0.03, 0.117, ..., 0.0928, 0.0212] type='138398 * float64'>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.weights[in_fid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "aad2c51b-69ac-4512-9cd8-4af36f10ab4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35929.91449331461"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(data.weights[in_fid])#int_probs[data.in_fiducial]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f4bc6e0c-4919-43e1-aa74-e5c3c5e7146d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import siren\n",
    "import os\n",
    "# Now include DIS interaction\n",
    "cross_section_model = \"HNLDISSplines\"\n",
    "\n",
    "xsfiledir = siren.utilities.get_cross_section_model_path(cross_section_model)\n",
    "\n",
    "# Cross Section Model\n",
    "target_type = siren.dataclasses.Particle.ParticleType.Nucleon\n",
    "\n",
    "nutype = \"nu\"\n",
    "primary = 14\n",
    "primary_type = (siren.dataclasses.Particle.ParticleType)(primary)\n",
    "\n",
    "m4 = \"0000300\"\n",
    "Ue4 = 0\n",
    "Umu4 = 0.01\n",
    "Utau4 = 0\n",
    "\n",
    "HNL_DIS_xs = siren.interactions.HNLDISFromSpline(\n",
    "    os.path.join(xsfiledir, \"M_0000000MeV/dsdxdy-%s-N-nc-GRV98lo_patched_central.fits\"%nutype),\n",
    "    os.path.join(xsfiledir, \"M_%sMeV/sigma-%s-N-nc-GRV98lo_patched_central.fits\"%(m4,nutype)),\n",
    "    float(m4)*1e-3,\n",
    "    [Ue4,Umu4,Utau4],\n",
    "    siren.utilities.Constants.isoscalarMass,\n",
    "    1,\n",
    "    [primary_type],\n",
    "    [target_type],\n",
    ")\n",
    "xsfiledir = \"/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/pweigel/cross_sections/20241017\"\n",
    "\n",
    "# Cross Section Model\n",
    "target_type = siren.dataclasses.Particle.ParticleType.Nucleon\n",
    "nutype = \"neutrino\"\n",
    "nuflavor = \"muon\"\n",
    "xs_mode = \"NC\"\n",
    "minQ2 = 1\n",
    "\n",
    "DIS_xs = siren.interactions.DISFromSpline(\n",
    "        os.path.join(xsfiledir, \"wcg24b_dsdxdy_%s_%s_%s_isoscalar.fits\"%(xs_mode,nuflavor,nutype)),\n",
    "        os.path.join(xsfiledir, \"wcg24b_sigma_%s_%s_%s_isoscalar.fits\"%(xs_mode,nuflavor,nutype)),\n",
    "        1,siren.utilities.Constants.isoscalarMass,minQ2,\n",
    "        [primary_type],\n",
    "        [target_type], \"cm\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b624bc2c-f619-41b6-b52c-ab14ce2b9339",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9182853408097895e-36"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "DIS_xs.TotalCrossSection(primary_type,1000) # cm^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "325a80e2-6229-4c79-b66c-30deedeae0ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9821680062720344e-40"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HNL_DIS_xs.TotalCrossSection(primary_type,1000)# * 1e18 * siren.utilities.Constants.hbarc**2 * 1e4 # cm^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6671cce4-55cc-4038-97d9-43b1db44029b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lienv",
   "language": "python",
   "name": "lienv"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}