diff --git a/README.md b/README.md
index e91ff7a..fa1bd53 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,7 @@
 [![PyPI version](https://badge.fury.io/py/lvec.svg)](https://badge.fury.io/py/lvec)
 [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
 
-> ⚠️ This project is a work in progress
+>  This project is a work in progress
 
 A Python package for seamless handling of Lorentz vectors, 2D vectors, and 3D vectors in HEP analysis, bridging the gap between Scikit-HEP and ROOT ecosystems.
 
@@ -19,6 +19,11 @@ LVec aims to simplify HEP analysis by providing a unified interface for working
 pip install lvec
 ```
 
+For JIT acceleration support (recommended):
+```bash
+pip install lvec numba
+```
+
 For development installation:
 ```bash
 git clone https://github.com/MohamedElashri/lvec
@@ -64,6 +69,30 @@ stats = v.cache_stats
 print(f"Cache hit ratio: {v.cache_hit_ratio}")
 ```
 
+### JIT Acceleration
+```python
+from lvec import LVec, is_jit_available, enable_jit
+import numpy as np
+
+# Check if JIT is available (requires numba)
+print(f"JIT acceleration available: {is_jit_available()}")
+
+# Create a large array of particles
+px = np.random.normal(0, 10, 1_000_000)
+py = np.random.normal(0, 10, 1_000_000)
+pz = np.random.normal(0, 10, 1_000_000)
+E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)  # pion mass ~ 0.14 GeV
+
+# Process with JIT enabled (default if numba is installed)
+vectors = LVec(px, py, pz, E)
+_ = vectors.pt  # JIT-accelerated calculation
+_ = vectors.mass  # JIT-accelerated calculation
+
+# Disable JIT for debugging
+enable_jit(False)
+_ = vectors.pt  # Now uses non-JIT implementation
+```
+
 ### Reference Frames
 ```python
 from lvec import LVec, Frame
@@ -202,6 +231,7 @@ vectors_ak = LVec(ak.Array(px), ak.Array(py), ak.Array(pz), ak.Array(E))
 - Python >= 3.10
 - NumPy >= 1.20.0
 - Awkward >= 2.0.0
+- Numba (optional, for JIT acceleration)
 
 For development:
 - pytest >= 7.0.0
diff --git a/docs/README.md b/docs/README.md
index be0136e..9669fd1 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -9,6 +9,7 @@
 7. [Advanced Usage](#advanced-usage)
 8. [Performance Considerations](#performance-considerations)
 9. [Reference Frames](#reference-frames)
+10. [JIT Acceleration](#jit-acceleration)
 
 ## Overview
 
@@ -18,6 +19,7 @@ LVec is a Python package designed for High Energy Physics (HEP) analysis, provid
 - Python 3.10+ (due to numpy requirement)
 - NumPy (required)
 - Awkward Array (optional, for Awkward array support)
+- Numba (optional, for JIT acceleration)
 
 ## Installation
 
@@ -831,6 +833,75 @@ with uproot.recreate("output.root") as f:
    selected = vectors[mask]
    ```
 
+## JIT Acceleration
+
+LVec provides JIT (Just-In-Time) compilation support to accelerate computationally intensive operations, particularly when working with large arrays of particles.
+
+### How JIT Acceleration Works
+
+The JIT acceleration feature in LVec uses `Numba` to compile critical parts of the code directly to optimized machine code. This provides significant performance improvements for operations like calculating masses and transverse momentum.
+
+JIT acceleration is fully compatible with the caching system, providing optimal performance through both compilation and caching of results:
+1. First calculation: JIT-accelerated computation → cached result
+2. Subsequent calculations: Retrieved from cache (fastest)
+
+### Using JIT Acceleration
+
+```python
+from lvec import LVec, is_jit_available, enable_jit
+import numpy as np
+
+# Check if JIT is available (requires numba)
+if is_jit_available():
+    print("JIT acceleration enabled")
+else:
+    print("JIT acceleration not available, install numba for better performance")
+    
+# Create vectors (JIT is automatically used if available)
+px = np.random.normal(0, 10, 1_000_000)  # Large array is better JIT benefit
+py = np.random.normal(0, 10, 1_000_000)
+pz = np.random.normal(0, 10, 1_000_000)
+E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)  # pion mass ~ 0.14 GeV
+
+vectors = LVec(px, py, pz, E)
+
+# Perform calculations (automatically JIT-accelerated if available)
+mass = vectors.mass  # JIT accelerated
+pt = vectors.pt      # JIT accelerated
+
+# Disable JIT if needed (for debugging or testing)
+enable_jit(False)
+```
+
+### JIT Configuration
+
+- **Global enabling/disabling**: Use `enable_jit(True/False)` to control JIT compilation globally
+- **Availability check**: Use `is_jit_available()` to check if JIT is available (requires numba package)
+- **Automatic fallback**: If numba is not installed or JIT is disabled, the code automatically falls back to standard implementations
+
+### Supported Operations
+
+The following operations benefit from JIT acceleration:
+
+| Operation | Description | Typical Speedup |
+|-----------|-------------|----------------|
+| `pt`      | Transverse momentum | 2-3x |
+| `p`       | Total momentum | 2-3x |
+| `mass`    | Invariant mass | 1.5-3x |
+| `eta`     | Pseudorapidity | 1.5-2x |
+| `phi`     | Azimuthal angle | 1-1.5x |
+| Constructors | Creation from pt,eta,phi,m | 1.5-2x |
+
+Note: Speedup factors depend on array size, hardware, and specific operation. JIT works best with large NumPy arrays.
+
+### Best Practices
+
+1. **Use larger arrays**: JIT compilation has overhead, so benefits increase with larger dataset sizes
+2. **Batch operations**: Apply operations to full arrays, not in Python loops
+3. **NumPy vs Awkward**: JIT acceleration works with NumPy arrays, while Awkward arrays use standard implementations
+4. **First-run overhead**: First execution includes compilation time; subsequent runs are much faster
+5. **JIT+Cache synergy**: Use both for optimal performance - JIT for fast computation, caching for reuse
+
 ## Reference Frames
 
 The `Frame` class provides a powerful abstraction for handling reference frame transformations in relativistic calculations.
@@ -892,4 +963,3 @@ In the center-of-mass frame, the total momentum should be zero:
 total_p_cm = sum(p.to_frame(cm_frame) for p in particles)
 print(f"Total momentum in CM: ({total_p_cm.px}, {total_p_cm.py}, {total_p_cm.pz})")
 # Should be very close to (0, 0, 0)
-```
diff --git a/examples/12_jit_acceleration.py b/examples/12_jit_acceleration.py
new file mode 100644
index 0000000..d41d031
--- /dev/null
+++ b/examples/12_jit_acceleration.py
@@ -0,0 +1,585 @@
+#!/usr/bin/env python3
+"""
+JIT Acceleration Demo
+====================
+
+This example demonstrates how to use the JIT (Just-In-Time) compilation 
+feature in the lvec package to accelerate performance-critical calculations.
+
+It shows:
+1. How to check if JIT compilation is available
+2. How to enable/disable JIT compilation globally
+3. Performance comparison between JIT-enabled and JIT-disabled code
+4. How JIT acceleration works with the caching system
+5. Best practices for JIT acceleration
+6. Pre-compilation to eliminate first-run overhead
+7. Adaptive batch processing thresholds
+
+Note: This example requires the numba package to demonstrate JIT acceleration.
+If numba is not installed, the example will show how lvec gracefully falls back
+to non-JIT implementations.
+"""
+
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+from lvec import (LVec, Vector3D, is_jit_available, enable_jit, 
+                 precompile_jit_functions, get_jit_batch_threshold)
+
+def check_jit_status():
+    """Check if JIT compilation is available and enabled."""
+    print("\n=== JIT Acceleration Status ===")
+    
+    if is_jit_available():
+        print(" JIT compilation is available")
+        print("   Using numba for accelerated computations")
+        print(f"   Current batch threshold: {get_jit_batch_threshold():,}")
+    else:
+        print(" JIT compilation is not available")
+        print("   To enable JIT, install numba: pip install numba")
+    
+    # Try enabling JIT (will only work if numba is installed)
+    original_status = is_jit_available()
+    enable_jit(True)
+    new_status = is_jit_available()
+    
+    if original_status != new_status:
+        print(f"JIT status changed: {original_status} -> {new_status}")
+    else:
+        print(f"JIT status unchanged: {new_status}")
+
+def demonstrate_precompilation():
+    """Demonstrate the benefits of pre-compilation."""
+    print("\n=== JIT Pre-compilation Demo ===")
+    
+    if not is_jit_available():
+        print("JIT is not available, skipping pre-compilation demo")
+        return
+    
+    # Create test data
+    n_particles = 1_000_000
+    rng = np.random.RandomState(42)
+    px = rng.normal(0, 10, n_particles)
+    py = rng.normal(0, 10, n_particles)
+    pz = rng.normal(0, 10, n_particles)
+    E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)
+    
+    # First, disable JIT to reset any compiled functions
+    enable_jit(False)
+    enable_jit(True)
+    
+    print("\n1. Without pre-compilation (first run includes compilation overhead)")
+    # Create vectors and time the first access to properties
+    vectors = LVec(px, py, pz, E)
+    
+    start_time = time.time()
+    _ = vectors.pt
+    pt_time = time.time() - start_time
+    print(f"   First pt access: {pt_time:.6f} seconds")
+    
+    start_time = time.time()
+    _ = vectors.mass
+    mass_time = time.time() - start_time
+    print(f"   First mass access: {mass_time:.6f} seconds")
+    
+    # Now reset and use pre-compilation
+    enable_jit(False)
+    enable_jit(True)
+    
+    print("\n2. With pre-compilation (compilation done ahead of time)")
+    # Pre-compile JIT functions
+    start_time = time.time()
+    precompile_jit_functions()
+    precompile_time = time.time() - start_time
+    print(f"   Pre-compilation time: {precompile_time:.6f} seconds")
+    
+    # Create vectors and time the first access to properties
+    vectors = LVec(px, py, pz, E)
+    
+    start_time = time.time()
+    _ = vectors.pt
+    pt_time_precompiled = time.time() - start_time
+    print(f"   First pt access: {pt_time_precompiled:.6f} seconds")
+    
+    start_time = time.time()
+    _ = vectors.mass
+    mass_time_precompiled = time.time() - start_time
+    print(f"   First mass access: {mass_time_precompiled:.6f} seconds")
+    
+    # Calculate speedup
+    if pt_time > 0 and pt_time_precompiled > 0:
+        speedup = pt_time / pt_time_precompiled
+        print(f"\nSpeedup for first pt access: {speedup:.1f}x")
+    
+    if mass_time > 0 and mass_time_precompiled > 0:
+        speedup = mass_time / mass_time_precompiled
+        print(f"Speedup for first mass access: {speedup:.1f}x")
+    
+    return {
+        "without_precompile": {"pt": pt_time, "mass": mass_time},
+        "with_precompile": {"pt": pt_time_precompiled, "mass": mass_time_precompiled},
+        "precompile_time": precompile_time
+    }
+
+def benchmark_jit_performance():
+    """Benchmark the performance difference between JIT-enabled and disabled code."""
+    print("\n=== JIT Performance Benchmark ===")
+    
+    # Create a large dataset for benchmarking
+    n_particles = 20_000_000  # Increased from 10M to 20M for better demonstration
+    print(f"Creating test dataset with {n_particles:,} particles...")
+    
+    # Generate random particle data
+    rng = np.random.RandomState(42)  # Fixed seed for reproducibility
+    px = rng.normal(0, 10, n_particles)
+    py = rng.normal(0, 10, n_particles)
+    pz = rng.normal(0, 10, n_particles)
+    E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)  # pion mass ~ 0.14 GeV
+    
+    # Define operations to benchmark
+    def benchmark_operations(vectors):
+        _ = vectors.pt      # Transverse momentum 
+        _ = vectors.eta     # Pseudorapidity
+        _ = vectors.phi     # Azimuthal angle
+        _ = vectors.mass    # Invariant mass
+        _ = vectors.p       # Total momentum
+        
+        # Additional operations that use multiple properties
+        _ = vectors.pt * np.cosh(vectors.eta)  # pz calculation
+        _ = vectors.pt**2 + vectors.pz**2      # p^2 calculation
+
+    # Run benchmark with JIT enabled vs disabled
+    results = {}
+    
+    for jit_enabled in [False, True]:  # Test non-JIT first, then JIT
+        # Set JIT status
+        enable_jit(jit_enabled)
+        status = "enabled" if jit_enabled else "disabled"
+        print(f"\nRunning benchmark with JIT {status}...")
+        
+        if jit_enabled:
+            # Pre-compile JIT functions to avoid compilation overhead
+            precompile_jit_functions()
+            print(f"  Pre-compiled JIT functions")
+            print(f"  Current batch threshold: {get_jit_batch_threshold():,}")
+        
+        # Create vectors (JIT setting affects vector creation too)
+        start_time = time.time()
+        vectors = LVec(px, py, pz, E)
+        create_time = time.time() - start_time
+        print(f"  Vector creation: {create_time:.4f} seconds")
+        
+        # First run (cold start)
+        start_time = time.time()
+        benchmark_operations(vectors)
+        cold_time = time.time() - start_time
+        print(f"  First run (cold): {cold_time:.4f} seconds")
+        
+        # Clear cache to ensure fair comparison
+        vectors.clear_cache()
+        
+        # Second run (cold again due to cache clear)
+        start_time = time.time()
+        benchmark_operations(vectors)
+        second_time = time.time() - start_time
+        print(f"  Second run (cold): {second_time:.4f} seconds")
+        
+        # Third run (should be warm with cached values)
+        start_time = time.time()
+        benchmark_operations(vectors)
+        warm_time = time.time() - start_time
+        print(f"  Third run (warm): {warm_time:.4f} seconds")
+        
+        results[jit_enabled] = {
+            "create": create_time,
+            "cold": cold_time,
+            "second": second_time,
+            "warm": warm_time
+        }
+    
+    # Calculate and display speedup
+    if False in results and True in results:
+        print("\n=== JIT Speedup Summary ===")
+        for phase in ["cold", "second"]:
+            speedup = results[False][phase] / results[True][phase]
+            print(f"  JIT speedup for {phase} run: {speedup:.2f}x")
+    
+    return results
+
+def benchmark_jit_with_caching():
+    """Show how JIT acceleration works together with the caching system."""
+    print("\n=== JIT with Caching System ===")
+    
+    # Enable JIT for this test
+    enable_jit(True)
+    precompile_jit_functions()
+    
+    # Create test data
+    n_particles = 1_000_000  # Increased from 100K to 1M
+    rng = np.random.RandomState(42)
+    px = rng.normal(0, 10, n_particles)
+    py = rng.normal(0, 10, n_particles)
+    pz = rng.normal(0, 10, n_particles)
+    E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)
+    
+    vectors = LVec(px, py, pz, E)
+    
+    # Reset cache counters
+    vectors._cache.reset_counters()
+    
+    print("\n1. First access (JIT-accelerated calculation + cache store)")
+    start_time = time.time()
+    _ = vectors.mass  # This should use JIT
+    first_time = time.time() - start_time
+    print(f"   Time: {first_time:.6f} seconds")
+    
+    print("\n2. Second access (cache retrieval, no calculation)")
+    start_time = time.time()
+    _ = vectors.mass  # This should use cache
+    cached_time = time.time() - start_time
+    print(f"   Time: {cached_time:.6f} seconds")
+    
+    # Print cache stats
+    stats = vectors._cache.get_stats()
+    print(f"\nCache hit ratio for mass: {stats['properties'].get('mass', {}).get('hit_ratio', 0):.2%}")
+    
+    # Show the speedup from caching vs JIT vs both
+    if first_time > 0:
+        cache_speedup = first_time / cached_time
+        print(f"\nSpeedup from caching: {cache_speedup:.1f}x")
+    
+    print("\n3. Clearing cache and disabling JIT")
+    vectors._cache.clear_cache()
+    enable_jit(False)
+    
+    print("\n4. Access with JIT disabled (non-JIT calculation + cache store)")
+    start_time = time.time()
+    _ = vectors.mass  # This should use non-JIT implementation
+    nojit_time = time.time() - start_time
+    print(f"   Time: {nojit_time:.6f} seconds")
+    
+    # Calculate JIT speedup if available
+    if nojit_time > 0 and first_time > 0:
+        jit_speedup = nojit_time / first_time
+        print(f"\nSpeedup from JIT: {jit_speedup:.1f}x")
+    
+    return first_time, cached_time, nojit_time
+
+def benchmark_vector_operations():
+    """Benchmark JIT performance for different vector operations."""
+    print("\n=== JIT Performance Across Different Operations ===")
+    
+    # Create test data
+    n_particles = 5_000_000  # Increased from 500K to 5M
+    rng = np.random.RandomState(42)
+    px = rng.normal(0, 10, n_particles)
+    py = rng.normal(0, 10, n_particles)
+    pz = rng.normal(0, 10, n_particles)
+    E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)
+    
+    # Define operations to benchmark
+    operations = {
+        "pt": lambda v: v.pt,
+        "eta": lambda v: v.eta,
+        "phi": lambda v: v.phi,
+        "mass": lambda v: v.mass,
+        "p": lambda v: v.p,
+    }
+    
+    # Run benchmarks
+    results = {}
+    
+    for jit_enabled in [False, True]:  # Test non-JIT first, then JIT
+        enable_jit(jit_enabled)
+        status = "enabled" if jit_enabled else "disabled"
+        print(f"\nRunning with JIT {status}...")
+        
+        if jit_enabled:
+            # Pre-compile JIT functions to avoid compilation overhead
+            precompile_jit_functions()
+        
+        # Create vectors
+        vectors = LVec(px, py, pz, E)
+        
+        # Benchmark each operation
+        op_results = {}
+        for name, op in operations.items():
+            # Clear cache to ensure cold start
+            vectors.clear_cache()
+            
+            # Time the operation
+            start_time = time.time()
+            _ = op(vectors)
+            op_time = time.time() - start_time
+            
+            print(f"  {name}: {op_time:.6f} seconds")
+            op_results[name] = op_time
+        
+        results[jit_enabled] = op_results
+    
+    # Calculate speedups
+    if False in results and True in results:
+        print("\n=== Operation-specific JIT Speedups ===")
+        for op in operations:
+            speedup = results[False][op] / results[True][op]
+            print(f"  {op}: {speedup:.2f}x faster with JIT")
+    
+    return results
+
+def benchmark_array_sizes():
+    """Benchmark JIT performance across different array sizes to show adaptive thresholds."""
+    print("\n=== JIT Performance Across Array Sizes ===")
+    
+    # Test different array sizes
+    sizes = [1_000, 10_000, 100_000, 1_000_000, 10_000_000]
+    operations = ["pt", "mass"]
+    
+    results = {"sizes": sizes, "operations": {}}
+    
+    # Enable JIT and pre-compile
+    enable_jit(True)
+    precompile_jit_functions()
+    
+    for op in operations:
+        print(f"\nBenchmarking {op} calculation across array sizes:")
+        op_results = {"jit": [], "nojit": [], "speedup": []}
+        
+        for size in sizes:
+            print(f"  Array size: {size:,}")
+            
+            # Generate data
+            rng = np.random.RandomState(42)
+            px = rng.normal(0, 10, size)
+            py = rng.normal(0, 10, size)
+            pz = rng.normal(0, 10, size)
+            E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)
+            
+            # Test with JIT enabled
+            enable_jit(True)
+            vectors_jit = LVec(px, py, pz, E)
+            vectors_jit.clear_cache()
+            
+            start_time = time.time()
+            if op == "pt":
+                _ = vectors_jit.pt
+            elif op == "mass":
+                _ = vectors_jit.mass
+            jit_time = time.time() - start_time
+            
+            # Test with JIT disabled
+            enable_jit(False)
+            vectors_nojit = LVec(px, py, pz, E)
+            vectors_nojit.clear_cache()
+            
+            start_time = time.time()
+            if op == "pt":
+                _ = vectors_nojit.pt
+            elif op == "mass":
+                _ = vectors_nojit.mass
+            nojit_time = time.time() - start_time
+            
+            # Calculate speedup
+            speedup = nojit_time / jit_time if jit_time > 0 else 0
+            
+            print(f"    JIT: {jit_time:.6f}s, Non-JIT: {nojit_time:.6f}s, Speedup: {speedup:.2f}x")
+            
+            op_results["jit"].append(jit_time)
+            op_results["nojit"].append(nojit_time)
+            op_results["speedup"].append(speedup)
+        
+        results["operations"][op] = op_results
+    
+    return results
+
+def plot_benchmark_results(jit_results, operation_results, precompile_results=None, array_size_results=None):
+    """Plot benchmark results for visualization."""
+    plt.figure(figsize=(15, 10))
+    
+    # Plot 1: JIT vs non-JIT for different phases
+    plt.subplot(2, 2, 1)
+    phases = ['cold', 'second', 'warm']
+    jit_times = [jit_results[True][p] for p in phases]
+    nojit_times = [jit_results[False][p] for p in phases]
+    
+    x = np.arange(len(phases))
+    width = 0.35
+    
+    plt.bar(x - width/2, nojit_times, width, label='JIT Disabled')
+    plt.bar(x + width/2, jit_times, width, label='JIT Enabled')
+    
+    plt.xlabel('Benchmark Phase')
+    plt.ylabel('Time (seconds)')
+    plt.title('JIT vs non-JIT Performance')
+    plt.xticks(x, phases)
+    plt.legend()
+    
+    # Add speedup annotations
+    for i, (nojit, jit) in enumerate(zip(nojit_times, jit_times)):
+        speedup = nojit / jit if jit > 0 else 0
+        plt.text(i, max(nojit, jit) + 0.05, f'{speedup:.1f}x', 
+                ha='center', va='bottom', fontweight='bold')
+    
+    # Plot 2: Operation-specific performance
+    plt.subplot(2, 2, 2)
+    ops = list(operation_results[True].keys())
+    jit_op_times = [operation_results[True][op] for op in ops]
+    nojit_op_times = [operation_results[False][op] for op in ops]
+    
+    x = np.arange(len(ops))
+    
+    plt.bar(x - width/2, nojit_op_times, width, label='JIT Disabled')
+    plt.bar(x + width/2, jit_op_times, width, label='JIT Enabled')
+    
+    plt.xlabel('Operation')
+    plt.ylabel('Time (seconds)')
+    plt.title('Operation-specific JIT Performance')
+    plt.xticks(x, ops)
+    plt.legend()
+    
+    # Add speedup annotations
+    for i, (nojit, jit) in enumerate(zip(nojit_op_times, jit_op_times)):
+        speedup = nojit / jit if jit > 0 else 0
+        plt.text(i, max(nojit, jit) + 0.01, f'{speedup:.1f}x', 
+                ha='center', va='bottom', fontweight='bold')
+    
+    # Plot 3: Precompilation benefits or speedup factors
+    plt.subplot(2, 2, 3)
+    
+    if precompile_results:
+        # Plot precompilation benefits
+        operations = ["pt", "mass"]
+        without_precompile = [precompile_results["without_precompile"][op] for op in operations]
+        with_precompile = [precompile_results["with_precompile"][op] for op in operations]
+        
+        x = np.arange(len(operations))
+        
+        plt.bar(x - width/2, without_precompile, width, label='Without Precompilation')
+        plt.bar(x + width/2, with_precompile, width, label='With Precompilation')
+        
+        plt.xlabel('Operation')
+        plt.ylabel('Time (seconds)')
+        plt.title('Precompilation Benefits')
+        plt.xticks(x, operations)
+        plt.legend()
+        
+        # Add speedup annotations
+        for i, (wo, w) in enumerate(zip(without_precompile, with_precompile)):
+            speedup = wo / w if w > 0 else 0
+            plt.text(i, max(wo, w) + 0.001, f'{speedup:.1f}x', 
+                    ha='center', va='bottom', fontweight='bold')
+    else:
+        # Calculate speedups
+        phase_speedups = [nojit_times[i]/jit_times[i] if jit_times[i] > 0 else 0 
+                         for i in range(len(phases))]
+        op_speedups = [nojit_op_times[i]/jit_op_times[i] if jit_op_times[i] > 0 else 0 
+                      for i in range(len(ops))]
+        
+        categories = phases + ops
+        speedups = phase_speedups + op_speedups
+        colors = ['blue']*len(phases) + ['green']*len(ops)
+        
+        plt.bar(categories, speedups, color=colors)
+        plt.axhline(y=1.0, color='r', linestyle='-', alpha=0.3)
+        
+        plt.xlabel('Category')
+        plt.ylabel('Speedup Factor (higher is better)')
+        plt.title('JIT Speedup Factors')
+        plt.xticks(rotation=45)
+        
+        # Add value labels
+        for i, v in enumerate(speedups):
+            plt.text(i, v + 0.1, f'{v:.1f}x', ha='center', va='bottom')
+    
+    # Plot 4: Array size performance or cache vs JIT
+    plt.subplot(2, 2, 4)
+    
+    if array_size_results:
+        # Plot array size performance
+        sizes = array_size_results["sizes"]
+        op = "mass"  # Choose one operation to display
+        
+        if op in array_size_results["operations"]:
+            speedups = array_size_results["operations"][op]["speedup"]
+            
+            plt.semilogx(sizes, speedups, 'o-', linewidth=2, markersize=8)
+            plt.axhline(y=1.0, color='r', linestyle='--', alpha=0.3)
+            
+            plt.xlabel('Array Size')
+            plt.ylabel('Speedup Factor (higher is better)')
+            plt.title(f'JIT Speedup vs Array Size for {op}')
+            plt.grid(True, alpha=0.3)
+            
+            # Add value labels
+            for i, (size, speedup) in enumerate(zip(sizes, speedups)):
+                plt.text(size, speedup + 0.1, f'{speedup:.1f}x', 
+                        ha='center', va='bottom')
+    else:
+        # This assumes benchmark_jit_with_caching has been run
+        # and returned the times
+        try:
+            first_time, cached_time, nojit_time = benchmark_jit_with_caching()
+            
+            times = [nojit_time, first_time, cached_time]
+            labels = ['non-JIT', 'JIT', 'Cached']
+            
+            plt.bar(labels, times)
+            plt.yscale('log')  # Log scale to show the dramatic difference with caching
+            
+            plt.xlabel('Access Method')
+            plt.ylabel('Time (seconds, log scale)')
+            plt.title('JIT vs Caching Performance')
+            
+            # Add speedup annotations
+            jit_speedup = nojit_time / first_time if first_time > 0 else 0
+            cache_speedup = first_time / cached_time if cached_time > 0 else 0
+            total_speedup = nojit_time / cached_time if cached_time > 0 else 0
+            
+            plt.text(0, nojit_time, f'{nojit_time:.6f}s', ha='center', va='bottom')
+            plt.text(1, first_time, f'{first_time:.6f}s\n({jit_speedup:.1f}x faster)', ha='center', va='bottom')
+            plt.text(2, cached_time, f'{cached_time:.6f}s\n({total_speedup:.1f}x faster)', ha='center', va='bottom')
+            
+        except Exception as e:
+            plt.text(0.5, 0.5, f"Error generating cache vs JIT plot: {e}", 
+                    ha='center', va='center', transform=plt.gca().transAxes)
+    
+    plt.tight_layout()
+    plt.savefig('jit_benchmark_results.pdf')
+    print("\nBenchmark plot saved as 'jit_benchmark_results.pdf'")
+    plt.show()
+
+def main():
+    """Run the JIT acceleration demonstration."""
+    print("=" * 80)
+    print("LVec JIT Acceleration Demonstration")
+    print("=" * 80)
+    
+    # Check if JIT is available
+    check_jit_status()
+    
+    # Demonstrate precompilation benefits
+    precompile_results = demonstrate_precompilation()
+    
+    # Run the benchmarks
+    jit_results = benchmark_jit_performance()
+    operation_results = benchmark_vector_operations()
+    
+    # Benchmark different array sizes to show adaptive thresholds
+    array_size_results = benchmark_array_sizes()
+    
+    # Plot benchmark results for visualization
+    plot_benchmark_results(jit_results, operation_results, 
+                          precompile_results, array_size_results)
+    
+    print("\n" + "=" * 80)
+    print("JIT Acceleration Best Practices:")
+    print("=" * 80)
+    print("1. JIT works best with large NumPy arrays (>100K elements)")
+    print("2. Use precompile_jit_functions() to eliminate first-run compilation overhead")
+    print("3. For maximum performance, combine JIT with the caching system")
+    print("4. JIT only works with NumPy arrays, not with Awkward arrays")
+    print("5. Enable JIT globally with: from lvec import enable_jit; enable_jit(True)")
+    print("6. The adaptive threshold automatically optimizes batch processing for your hardware")
+    print("=" * 80)
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/README.md b/examples/README.md
index 10d5467..253dea3 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -6,11 +6,15 @@ This directory contains example scripts demonstrating how to use the LVec packag
 Make sure you have LVec installed along with its dependencies:
 ```bash
 pip install lvec
-pip install uproot numpy awkward
+```
+
+For JIT acceleration examples:
+```bash
+pip install numba
 ```
 
 ## Data Sample
-The examples use a simulated Z→μμ decay sample. To generate the sample:
+The examples use a simulated Z decay sample. To generate the sample:
 
 ```bash
 python create_test_data.py
@@ -18,8 +22,8 @@ python create_test_data.py
 
 This will create `samples/physics_data.root` containing:
 - Mother particle (Z boson): px, py, pz, E
-- Daughter 1 (μ⁺): px, py, pz, E
-- Daughter 2 (μ⁻): px, py, pz, E
+- Daughter 1 (+): px, py, pz, E
+- Daughter 2 (): px, py, pz, E
 
 ## Available Examples
 
@@ -108,7 +112,7 @@ Shows how to:
 
 ### 8. LHCb Analysis (`08_lhcb_data.py`)
 Demonstrates how to:
-- Work with real LHCb open data (B→hhh decay)
+- Work with real LHCb open data (Bhhh decay)
 - Calculate two-body invariant masses
 - Create publication-quality plots with LHCb style
 - Handle multiple particle combinations
@@ -131,7 +135,7 @@ pip install uproot matplotlib
 Uses LHCb open data from:
 - Dataset: B2HHH_MagnetDown.root (just a sample)
 - Source: [CERN Open Data](https://opendata.cern.ch/record/4900)
-- Description: B→hhh decay data collected by LHCb in 2011 at √s = 7 TeV
+- Description: Bhhh decay data collected by LHCb in 2011 at s = 7 TeV
 
 #### Running the Example
 ```bash
@@ -176,19 +180,94 @@ total_cm = p1_cm + p2_cm
 print(f"Total momentum in CM: px={total_cm.px:.6f}, py={total_cm.py:.6f}, pz={total_cm.pz:.6f}")
 ```
 
-Sample output:
+### 11. Cache Optimization Demo (`11_cache_optimization_demo.py`)
+Demonstrates advanced cache optimization features:
+- Setting Time-To-Live (TTL) for cached properties
+- Configuring LRU (Least Recently Used) eviction policies
+- Monitoring cache performance
+
+```python
+# Configure cache with size limit and TTL
+vector = LVec(px, py, pz, E, max_cache_size=100, default_ttl=60)
+
+# Check cache size and clear expired entries
+print(f"Cache size: {vector.cache_size}")
+expired_count = vector.clear_expired()
 ```
-=== Particles in center-of-mass frame ===
-Particle 1: px=0.000000, py=0.000000, pz=17.329527, E=22.919697, mass=15.000000
-Particle 2: px=0.000000, py=0.000000, pz=-17.329527, E=21.801663, mass=13.228757
 
-Total momentum in CM: px=0.000000, py=0.000000, pz=0.000000
-Total energy in CM: E=44.721360
-Invariant mass in CM: √s=44.721360
+### 12. JIT Acceleration (`12_jit_acceleration.py`)
+Demonstrates how to use JIT (Just-In-Time) compilation to speed up calculations:
+```python
+from lvec import LVec, is_jit_available, enable_jit
+import numpy as np
+
+# Check if JIT is available
+print(f"JIT available: {is_jit_available()}")
+
+# Create particle vectors with JIT acceleration
+px = np.random.normal(0, 10, 10_000_000)  # Large dataset for better JIT benefit
+py = np.random.normal(0, 10, 10_000_000)
+pz = np.random.normal(0, 10, 10_000_000)
+E = np.sqrt(px**2 + py**2 + pz**2 + 0.14**2)  # pion mass ~ 0.14 GeV
+
+vectors = LVec(px, py, pz, E)  # JIT-enabled by default
+
+# Access properties (JIT-accelerated)
+pt = vectors.pt
+mass = vectors.mass
+
+# Disable JIT for debugging or comparison
+enable_jit(False)
 ```
 
+#### Running the Example
+```bash
+python 12_jit_acceleration.py
+```
 
+This will:
+1. Check if JIT acceleration is available (requires numba)
+2. Run a series of benchmarks comparing JIT and non-JIT performance
+3. Generate performance visualization plots
+4. Demonstrate how JIT works with the caching system
+
+#### Sample Output
+```
+=== JIT Performance Benchmark ===
+Creating test dataset with 10,000,000 particles...
+
+Running benchmark with JIT enabled...
+  Vector creation: 0.0049 seconds
+  First run (cold): 0.7907 seconds  # Includes JIT compilation time
+  Second run (cold): 0.4795 seconds
+  Third run (warm): 0.0886 seconds
+
+Running benchmark with JIT disabled...
+  First run (cold): 0.5160 seconds
+  Second run (cold): 0.5197 seconds
+  Third run (warm): 0.0911 seconds
+
+=== JIT Performance Across Different Operations ===
+Running operation benchmarks with JIT enabled...
+  pt   : 0.000514 seconds
+  eta  : 0.003842 seconds
+  phi  : 0.013952 seconds
+  mass : 0.003011 seconds
+  p    : 0.000760 seconds
+
+Running operation benchmarks with JIT disabled...
+  pt   : 0.001242 seconds (2.4x slower)
+  eta  : 0.006445 seconds (1.7x slower)
+  phi  : 0.013638 seconds (1.0x slower)
+  mass : 0.003630 seconds (1.2x slower)
+  p    : 0.001861 seconds (2.5x slower)
+```
 
+#### Dependencies
+Required for this example:
+```bash
+pip install numba matplotlib
+```
 ## Running the Examples
 
 Run each example individually:
@@ -200,6 +279,8 @@ python 04_boost_frame.py
 python 08_lhcb_data.py
 python 09_cache_performance.py
 python 10_reference_frames.py
+python 11_cache_optimization_demo.py
+python 12_jit_acceleration.py
 ```
 
 ## Expected Output
@@ -219,7 +300,7 @@ Original mass: 91.20 GeV
 Reconstructed mass: 0.22 GeV
 Mass resolution: 0.017 GeV
 
-Average ΔR between daughters: 0.000
+Average R between daughters: 0.000
 ```
 
 ### Physics Selections
@@ -229,8 +310,8 @@ Total events: 1000
 Selected events: 625
 
 Selected Z properties:
-Mass mean: 0.22 ± 0.02 GeV
-pT mean: 75.25 ± 15.55 GeV
+Mass mean: 0.22  0.02 GeV
+pT mean: 75.25  15.55 GeV
 ```
 
 ### Reference Frames
@@ -255,7 +336,7 @@ Azimuthal angle (phi): [1.32581766 1.19028995 1.10714872]
 
 Vector Operations:
 v1 + v2: ([6. 8.], [10. 12.])
-v1 · v2 (dot product): [26. 44.]
+v1  v2 (dot product): [26. 44.]
 
 === 2D Vector Operations with Awkward Backend ===
 
@@ -265,7 +346,7 @@ Azimuthal angle (phi): [1.33, 1.19, 1.11]
 
 Vector Operations:
 v1 + v2: ([6, 8], [10, 12])
-v1 · v2 (dot product): [26, 44]
+v1  v2 (dot product): [26, 44]
 ```
 
 ### 3D Vectors
@@ -281,8 +362,8 @@ Polar angle (theta): [0.5323032  0.59247462 0.64052231]
 
 Vector Operations:
 v1 + v2: ([ 8. 10.], [12. 14.], [16. 18.])
-v1 · v2 (dot product): [ 89. 128.]
-v1 × v2 (cross product): ([-12. -12.], [24. 24.], [-12. -12.])
+v1  v2 (dot product): [ 89. 128.]
+v1  v2 (cross product): ([-12. -12.], [24. 24.], [-12. -12.])
 
 === 3D Vector Operations with Awkward Backend ===
 
@@ -294,8 +375,8 @@ Polar angle (theta): [0.532, 0.592, 0.641]
 
 Vector Operations:
 v1 + v2: ([8, 10], [12, 14], [16, 18])
-v1 · v2 (dot product): [89, 128]
-v1 × v2 (cross product): ([-12, -12], [24, 24], [-12, -12])
+v1  v2 (dot product): [89, 128]
+v1  v2 (cross product): ([-12, -12], [24, 24], [-12, -12])
 ```
 
 ### LHCb Analysis
diff --git a/lvec/__init__.py b/lvec/__init__.py
index a294921..fdefabc 100644
--- a/lvec/__init__.py
+++ b/lvec/__init__.py
@@ -4,5 +4,43 @@
 from .exceptions import ShapeError
 from .frame import Frame
 
-__all__ = ['LVec', 'Vector2D', 'Vector3D', 'Frame', 'ShapeError']
+# Import JIT functions and settings
+try:
+    from .jit import (enable_jit, is_jit_available, precompile_jit_functions,
+                     BATCH_THRESHOLD, MIN_BATCH_THRESHOLD, MAX_BATCH_THRESHOLD)
+except ImportError:
+    # Create dummy functions if JIT is not available
+    def enable_jit(enable=True):
+        """Dummy function when Numba is not available."""
+        import warnings
+        warnings.warn("Numba not installed, JIT compilation not available")
+        return False
+        
+    def is_jit_available():
+        """Check if JIT compilation is available."""
+        return False
+        
+    def precompile_jit_functions():
+        """Dummy function when Numba is not available."""
+        return False
+        
+    # Dummy threshold values
+    BATCH_THRESHOLD = 100_000
+    MIN_BATCH_THRESHOLD = 10_000
+    MAX_BATCH_THRESHOLD = 1_000_000
+
+def get_jit_batch_threshold():
+    """
+    Get the current batch processing threshold for JIT operations.
+    
+    Returns
+    -------
+    int
+        Current threshold value. Arrays larger than this will use parallel batch processing.
+    """
+    return BATCH_THRESHOLD
+
+__all__ = ['LVec', 'Vector2D', 'Vector3D', 'Frame', 'ShapeError', 
+           'enable_jit', 'is_jit_available', 'precompile_jit_functions',
+           'get_jit_batch_threshold']
 __version__ = '0.1.4'
diff --git a/lvec/jit.py b/lvec/jit.py
new file mode 100644
index 0000000..4f950be
--- /dev/null
+++ b/lvec/jit.py
@@ -0,0 +1,583 @@
+"""
+JIT compilation support for performance-critical operations in LVec.
+
+This module provides Numba-accelerated versions of core computation functions
+while maintaining compatibility with LVec's caching system. Functions gracefully 
+fall back to non-JIT versions when Numba is not available.
+"""
+
+import numpy as np
+import time
+import warnings
+from functools import lru_cache
+
+try:
+    import numba as nb
+    HAS_NUMBA = True
+except ImportError:
+    HAS_NUMBA = False
+
+# Configuration
+JIT_ENABLED = HAS_NUMBA  # Global flag to enable/disable JIT
+OPTIMIZE_NUMBA = True    # Apply additional optimizations specific to Numba
+PARALLEL = True          # Enable parallel execution where possible
+FASTMATH = True          # Enable fast math optimizations (may reduce precision slightly)
+
+# Adaptive batch processing
+# Initial threshold - will be adjusted based on runtime performance
+BATCH_THRESHOLD = 100_000
+MIN_BATCH_THRESHOLD = 10_000
+MAX_BATCH_THRESHOLD = 1_000_000
+BATCH_THRESHOLD_ADJUSTMENT_FACTOR = 1.5  # How much to adjust threshold up/down
+
+# Performance tracking for adaptive threshold
+_performance_history = {}
+
+# Function to check if JIT is available and enabled
+def is_jit_available():
+    """Check if JIT compilation is available and enabled."""
+    return HAS_NUMBA and JIT_ENABLED
+
+def enable_jit(enable=True):
+    """Enable or disable JIT compilation globally."""
+    global JIT_ENABLED
+    if enable and not HAS_NUMBA:
+        warnings.warn("Numba not installed, JIT compilation not available")
+        return False
+    JIT_ENABLED = enable
+    return JIT_ENABLED
+
+# ------------------- JIT implementations of core computation functions -------------------
+
+if HAS_NUMBA:
+    # Numba-optimized functions for NumPy arrays
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_pt_np(px, py):
+        """JIT-optimized pt calculation for NumPy arrays."""
+        return np.sqrt(px*px + py*py)
+        
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_p_np(px, py, pz):
+        """JIT-optimized momentum calculation for NumPy arrays."""
+        return np.sqrt(px*px + py*py + pz*pz)
+        
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_mass_np(E, px, py, pz):
+        """JIT-optimized mass calculation for NumPy arrays."""
+        p_squared = px*px + py*py + pz*pz
+        m_squared = E*E - p_squared
+        
+        # Handle numerical precision issues that can lead to small negative mass-squared values
+        # For physical particles, m_squared should be positive
+        abs_m_squared = np.abs(m_squared)
+        return np.sqrt(abs_m_squared)
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_eta_np(p, pz):
+        """JIT-optimized eta calculation for NumPy arrays."""
+        epsilon = 1e-10  # Small constant to avoid division by zero
+        return 0.5 * np.log((p + pz) / (p - pz + epsilon))
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_phi_np(px, py):
+        """JIT-optimized phi calculation for NumPy arrays."""
+        return np.arctan2(py, px)
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_compute_p4_from_ptepm_np(pt, eta, phi, m):
+        """JIT-optimized conversion from pt,eta,phi,m to p4 for NumPy arrays."""
+        px = pt * np.cos(phi)
+        py = pt * np.sin(phi)
+        pz = pt * np.sinh(eta)
+        p = np.sqrt(px*px + py*py + pz*pz)
+        E = np.sqrt(p*p + m*m)
+        return px, py, pz, E
+    
+    # Batch processing functions for large arrays
+    @nb.njit(cache=True, fastmath=FASTMATH, parallel=PARALLEL)
+    def _jit_batch_compute_pt_np(px, py):
+        """Batch-optimized pt calculation for large NumPy arrays."""
+        result = np.empty_like(px)
+        n = len(px)
+        for i in nb.prange(n):
+            result[i] = np.sqrt(px[i]*px[i] + py[i]*py[i])
+        return result
+    
+    @nb.njit(cache=True, fastmath=FASTMATH, parallel=PARALLEL)
+    def _jit_batch_compute_p_np(px, py, pz):
+        """Batch-optimized momentum calculation for large NumPy arrays."""
+        result = np.empty_like(px)
+        n = len(px)
+        for i in nb.prange(n):
+            result[i] = np.sqrt(px[i]*px[i] + py[i]*py[i] + pz[i]*pz[i])
+        return result
+    
+    @nb.njit(cache=True, fastmath=FASTMATH, parallel=PARALLEL)
+    def _jit_batch_compute_mass_np(E, px, py, pz):
+        """Batch-optimized mass calculation for large NumPy arrays."""
+        result = np.empty_like(E)
+        n = len(E)
+        for i in nb.prange(n):
+            p_squared = px[i]*px[i] + py[i]*py[i] + pz[i]*pz[i]
+            m_squared = E[i]*E[i] - p_squared
+            # Handle numerical precision issues
+            abs_m_squared = abs(m_squared)
+            result[i] = np.sqrt(abs_m_squared)
+        return result
+    
+    # Vectorized versions that use NumPy's optimized operations directly
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_vectorized_compute_pt_np(px, py):
+        """Vectorized pt calculation using NumPy operations."""
+        return np.sqrt(px**2 + py**2)
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_vectorized_compute_p_np(px, py, pz):
+        """Vectorized momentum calculation using NumPy operations."""
+        return np.sqrt(px**2 + py**2 + pz**2)
+    
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_vectorized_compute_mass_np(E, px, py, pz):
+        """Vectorized mass calculation using NumPy operations."""
+        p_squared = px**2 + py**2 + pz**2
+        m_squared = E**2 - p_squared
+        return np.sqrt(np.abs(m_squared))
+    
+    # Optimized phi calculation that might perform better with JIT
+    @nb.njit(cache=True, fastmath=FASTMATH)
+    def _jit_optimized_phi_np(px, py):
+        """Optimized phi calculation that avoids arctan2 for common cases."""
+        result = np.empty_like(px)
+        n = len(px)
+        
+        for i in range(n):
+            x, y = px[i], py[i]
+            
+            # Handle special cases
+            if x == 0.0:
+                if y > 0.0:
+                    result[i] = np.pi/2
+                elif y < 0.0:
+                    result[i] = -np.pi/2
+                else:
+                    result[i] = 0.0
+                continue
+                
+            # Use arctan2 for general case
+            result[i] = np.arctan2(y, x)
+            
+        return result
+
+# Function to determine the optimal implementation based on array size
+def _select_optimal_implementation(array_size, operation):
+    """
+    Select the optimal implementation based on array size and operation.
+    
+    Parameters
+    ----------
+    array_size : int
+        Size of the input array
+    operation : str
+        Operation type ('pt', 'p', 'mass', 'eta', 'phi')
+        
+    Returns
+    -------
+    str
+        Implementation type ('vectorized', 'batch', 'standard')
+    """
+    # For very small arrays, use standard implementation
+    if array_size < MIN_BATCH_THRESHOLD:
+        return 'standard'
+    
+    # For medium-sized arrays, use vectorized implementation
+    if array_size < BATCH_THRESHOLD:
+        return 'vectorized'
+    
+    # For large arrays, use batch implementation
+    return 'batch'
+
+# Adaptive threshold adjustment
+def _update_batch_threshold(operation, implementation, array_size, execution_time):
+    """
+    Update the batch threshold based on performance data.
+    
+    Parameters
+    ----------
+    operation : str
+        Operation type ('pt', 'p', 'mass', 'eta', 'phi')
+    implementation : str
+        Implementation used ('vectorized', 'batch', 'standard')
+    array_size : int
+        Size of the input array
+    execution_time : float
+        Execution time in seconds
+    """
+    global BATCH_THRESHOLD
+    
+    # Only track performance for arrays near the threshold
+    if not (MIN_BATCH_THRESHOLD <= array_size <= MAX_BATCH_THRESHOLD):
+        return
+    
+    # Initialize performance history for this operation if needed
+    if operation not in _performance_history:
+        _performance_history[operation] = {'vectorized': [], 'batch': []}
+    
+    # Add performance data
+    if implementation in ('vectorized', 'batch'):
+        _performance_history[operation][implementation].append((array_size, execution_time))
+    
+    # Only adjust threshold if we have enough data points
+    if (len(_performance_history[operation]['vectorized']) >= 3 and 
+        len(_performance_history[operation]['batch']) >= 3):
+        
+        # Calculate average performance for each implementation
+        vec_sizes = [s for s, _ in _performance_history[operation]['vectorized']]
+        vec_times = [t for _, t in _performance_history[operation]['vectorized']]
+        batch_sizes = [s for s, _ in _performance_history[operation]['batch']]
+        batch_times = [t for _, t in _performance_history[operation]['batch']]
+        
+        # Find average performance per element
+        if vec_sizes and vec_times:
+            vec_perf = sum(t/s for s, t in zip(vec_sizes, vec_times)) / len(vec_sizes)
+        else:
+            vec_perf = float('inf')
+            
+        if batch_sizes and batch_times:
+            batch_perf = sum(t/s for s, t in zip(batch_sizes, batch_times)) / len(batch_sizes)
+        else:
+            batch_perf = float('inf')
+        
+        # Adjust threshold based on which is faster
+        if vec_perf < batch_perf:
+            # Vectorized is faster, increase threshold
+            new_threshold = min(MAX_BATCH_THRESHOLD, 
+                               int(BATCH_THRESHOLD * BATCH_THRESHOLD_ADJUSTMENT_FACTOR))
+        else:
+            # Batch is faster, decrease threshold
+            new_threshold = max(MIN_BATCH_THRESHOLD, 
+                               int(BATCH_THRESHOLD / BATCH_THRESHOLD_ADJUSTMENT_FACTOR))
+        
+        # Update threshold and reset history
+        if new_threshold != BATCH_THRESHOLD:
+            BATCH_THRESHOLD = new_threshold
+            _performance_history[operation] = {'vectorized': [], 'batch': []}
+
+# ------------------- JIT-aware wrapper functions -------------------
+
+# Helper for the non-JIT implementations
+def _get_original_compute_pt():
+    """Get the original compute_pt function without JIT."""
+    from lvec.utils import _compute_pt_original
+    return _compute_pt_original
+
+def _get_original_compute_p():
+    """Get the original compute_p function without JIT."""
+    from lvec.utils import _compute_p_original
+    return _compute_p_original
+
+def _get_original_compute_mass():
+    """Get the original compute_mass function without JIT."""
+    from lvec.utils import _compute_mass_original
+    return _compute_mass_original
+
+def _get_original_compute_eta():
+    """Get the original compute_eta function without JIT."""
+    from lvec.utils import _compute_eta_original
+    return _compute_eta_original
+
+def _get_original_compute_phi():
+    """Get the original compute_phi function without JIT."""
+    from lvec.utils import _compute_phi_original
+    return _compute_phi_original
+
+def _get_original_compute_p4_from_ptepm():
+    """Get the original compute_p4_from_ptepm function without JIT."""
+    from lvec.utils import _compute_p4_from_ptepm_original
+    return _compute_p4_from_ptepm_original
+
+def jit_compute_pt(px, py, lib):
+    """
+    JIT-optimized pt calculation with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    px, py : scalar or array-like
+        X and Y components of momentum
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    scalar or array-like
+        Transverse momentum with the same type as input
+    """
+    # Only use JIT for NumPy arrays, fall back for scalars and awkward arrays
+    if lib == 'np' and is_jit_available() and not isinstance(px, (float, int)):
+        array_size = getattr(px, 'size', len(px) if hasattr(px, '__len__') else 0)
+        
+        # Select implementation based on array size
+        start_time = time.time()
+        implementation = _select_optimal_implementation(array_size, 'pt')
+        
+        if implementation == 'batch':
+            result = _jit_batch_compute_pt_np(px, py)
+        elif implementation == 'vectorized':
+            result = _jit_vectorized_compute_pt_np(px, py)
+        else:
+            result = _jit_compute_pt_np(px, py)
+            
+        # Update adaptive threshold
+        execution_time = time.time() - start_time
+        _update_batch_threshold('pt', implementation, array_size, execution_time)
+            
+        # Convert to scalar if both inputs were scalars (shouldn't happen but just in case)
+        if isinstance(px, (float, int)) and isinstance(py, (float, int)):
+            return float(result)
+        return result
+    else:
+        # Let the original function handle awkward arrays and scalars
+        return _get_original_compute_pt()(px, py, lib)
+
+def jit_compute_p(px, py, pz, lib):
+    """
+    JIT-optimized momentum calculation with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    px, py, pz : scalar or array-like
+        X, Y, and Z components of momentum
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    scalar or array-like
+        Total momentum with the same type as input
+    """
+    if lib == 'np' and is_jit_available() and not isinstance(px, (float, int)):
+        array_size = getattr(px, 'size', len(px) if hasattr(px, '__len__') else 0)
+        
+        # Select implementation based on array size
+        start_time = time.time()
+        implementation = _select_optimal_implementation(array_size, 'p')
+        
+        if implementation == 'batch':
+            result = _jit_batch_compute_p_np(px, py, pz)
+        elif implementation == 'vectorized':
+            result = _jit_vectorized_compute_p_np(px, py, pz)
+        else:
+            result = _jit_compute_p_np(px, py, pz)
+            
+        # Update adaptive threshold
+        execution_time = time.time() - start_time
+        _update_batch_threshold('p', implementation, array_size, execution_time)
+            
+        # Convert to scalar if all inputs were scalars
+        if isinstance(px, (float, int)) and isinstance(py, (float, int)) and isinstance(pz, (float, int)):
+            return float(result)
+        return result
+    else:
+        # Let the original function handle awkward arrays and scalars
+        return _get_original_compute_p()(px, py, pz, lib)
+
+def jit_compute_mass(E, px, py, pz, lib):
+    """
+    JIT-optimized mass calculation with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    E : scalar or array-like
+        Energy
+    px, py, pz : scalar or array-like
+        Momentum components
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    scalar or array-like
+        Mass with the same type as input
+    """
+    if lib == 'np' and is_jit_available() and not isinstance(E, (float, int)):
+        array_size = getattr(E, 'size', len(E) if hasattr(E, '__len__') else 0)
+        
+        # Select implementation based on array size
+        start_time = time.time()
+        implementation = _select_optimal_implementation(array_size, 'mass')
+        
+        if implementation == 'batch':
+            result = _jit_batch_compute_mass_np(E, px, py, pz)
+        elif implementation == 'vectorized':
+            result = _jit_vectorized_compute_mass_np(E, px, py, pz)
+        else:
+            result = _jit_compute_mass_np(E, px, py, pz)
+            
+        # Update adaptive threshold
+        execution_time = time.time() - start_time
+        _update_batch_threshold('mass', implementation, array_size, execution_time)
+            
+        # Convert to scalar if all inputs were scalars
+        if isinstance(E, (float, int)) and isinstance(px, (float, int)) and isinstance(py, (float, int)) and isinstance(pz, (float, int)):
+            return float(result)
+        return result
+    else:
+        # For awkward arrays, we need a different approach since we calculate
+        # mass differently in the original function
+        from lvec.utils import compute_p
+        p = compute_p(px, py, pz, lib)
+        return _get_original_compute_mass()(E, p, lib)
+
+def jit_compute_eta(p, pz, lib):
+    """
+    JIT-optimized eta calculation with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    p : scalar or array-like
+        Total momentum
+    pz : scalar or array-like
+        Z component of momentum
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    scalar or array-like
+        Pseudorapidity with the same type as input
+    """
+    if lib == 'np' and is_jit_available() and not isinstance(p, (float, int)):
+        result = _jit_compute_eta_np(p, pz)
+        # Convert to scalar if all inputs were scalars
+        if isinstance(p, (float, int)) and isinstance(pz, (float, int)):
+            return float(result)
+        return result
+    else:
+        # Let the original function handle awkward arrays and scalars
+        return _get_original_compute_eta()(p, pz, lib)
+
+def jit_compute_phi(px, py, lib):
+    """
+    JIT-optimized phi calculation with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    px, py : scalar or array-like
+        X and Y components of momentum
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    scalar or array-like
+        Azimuthal angle with the same type as input
+    """
+    if lib == 'np' and is_jit_available() and not isinstance(px, (float, int)):
+        array_size = getattr(px, 'size', len(px) if hasattr(px, '__len__') else 0)
+        
+        # For phi, we'll try the optimized implementation for large arrays
+        if array_size > BATCH_THRESHOLD:
+            result = _jit_optimized_phi_np(px, py)
+        else:
+            result = _jit_compute_phi_np(px, py)
+            
+        # Convert to scalar if all inputs were scalars
+        if isinstance(px, (float, int)) and isinstance(py, (float, int)):
+            return float(result)
+        return result
+    else:
+        # Let the original function handle awkward arrays and scalars
+        return _get_original_compute_phi()(px, py, lib)
+
+def jit_compute_p4_from_ptepm(pt, eta, phi, m, lib):
+    """
+    JIT-optimized conversion from pt,eta,phi,m to p4 with backend awareness.
+    
+    Falls back to non-JIT implementation for non-NumPy backends or when JIT is disabled.
+    
+    Parameters
+    ----------
+    pt : scalar or array-like
+        Transverse momentum
+    eta : scalar or array-like
+        Pseudorapidity
+    phi : scalar or array-like
+        Azimuthal angle
+    m : scalar or array-like
+        Mass
+    lib : str
+        Backend library ('np' or 'ak')
+        
+    Returns
+    -------
+    tuple
+        (px, py, pz, E) with the same type as input
+    """
+    if lib == 'np' and is_jit_available() and not isinstance(pt, (float, int)):
+        return _jit_compute_p4_from_ptepm_np(pt, eta, phi, m)
+    else:
+        # Let the original function handle awkward arrays and scalars
+        return _get_original_compute_p4_from_ptepm()(pt, eta, phi, m, lib)
+
+# Pre-compilation function to eliminate first-run overhead
+def precompile_jit_functions():
+    """
+    Pre-compile all JIT functions to eliminate first-run compilation overhead.
+    
+    This function should be called once at module import time or when JIT is enabled.
+    """
+    if not is_jit_available():
+        return False
+    
+    try:
+        # Create small test arrays
+        size = 10
+        px = np.ones(size)
+        py = np.ones(size)
+        pz = np.ones(size)
+        E = np.ones(size) * 2
+        p = np.ones(size) * np.sqrt(3)
+        pt = np.ones(size) * np.sqrt(2)
+        eta = np.zeros(size)
+        phi = np.zeros(size)
+        m = np.ones(size)
+        
+        # Pre-compile standard functions
+        _ = _jit_compute_pt_np(px, py)
+        _ = _jit_compute_p_np(px, py, pz)
+        _ = _jit_compute_mass_np(E, px, py, pz)
+        _ = _jit_compute_eta_np(p, pz)
+        _ = _jit_compute_phi_np(px, py)
+        _ = _jit_compute_p4_from_ptepm_np(pt, eta, phi, m)
+        
+        # Pre-compile vectorized functions
+        _ = _jit_vectorized_compute_pt_np(px, py)
+        _ = _jit_vectorized_compute_p_np(px, py, pz)
+        _ = _jit_vectorized_compute_mass_np(E, px, py, pz)
+        
+        # Pre-compile batch functions
+        _ = _jit_batch_compute_pt_np(px, py)
+        _ = _jit_batch_compute_p_np(px, py, pz)
+        _ = _jit_batch_compute_mass_np(E, px, py, pz)
+        
+        # Pre-compile optimized phi
+        _ = _jit_optimized_phi_np(px, py)
+        
+        return True
+    except Exception as e:
+        warnings.warn(f"Failed to pre-compile JIT functions: {e}")
+        return False
+
+# Attempt to pre-compile if JIT is available
+if HAS_NUMBA and JIT_ENABLED:
+    precompile_jit_functions()
diff --git a/lvec/lvec.py b/lvec/lvec.py
index e6bb947..05f3781 100644
--- a/lvec/lvec.py
+++ b/lvec/lvec.py
@@ -10,6 +10,14 @@
 from .frame import Frame
 import numpy as np
 
+# Import JIT functions if available
+try:
+    from .jit import (jit_compute_pt, jit_compute_p, jit_compute_mass,
+                     jit_compute_eta, jit_compute_phi, jit_compute_p4_from_ptepm,
+                     is_jit_available)
+    HAS_JIT = True
+except ImportError:
+    HAS_JIT = False
 
 class LVec:
     """
@@ -82,7 +90,13 @@ def from_ptepm(cls, pt, eta, phi, m, max_cache_size=None, default_ttl=None):
         """Create from pt, eta, phi, mass."""
         # First convert to arrays and get the lib type
         pt, eta, phi, m, lib = ensure_array(pt, eta, phi, m)
-        px, py, pz, E = compute_p4_from_ptepm(pt, eta, phi, m, lib)
+        
+        # Use JIT if available for NumPy arrays
+        if lib == 'np' and HAS_JIT and is_jit_available() and not isinstance(pt, (float, int)):
+            px, py, pz, E = jit_compute_p4_from_ptepm(pt, eta, phi, m, lib)
+        else:
+            px, py, pz, E = compute_p4_from_ptepm(pt, eta, phi, m, lib)
+            
         return cls(px, py, pz, E, max_cache_size=max_cache_size, default_ttl=default_ttl)
     
     @classmethod
@@ -215,109 +229,63 @@ def _p_squared(self):
     @property
     def pt(self):
         """Transverse momentum."""
-        return self._get_cached('pt', 
-                             lambda: backend_sqrt(self._pt_squared(), self._lib),
-                             ['px', 'py'])
-    
+        def _compute_pt():
+            # Use JIT-optimized function when available for NumPy arrays
+            if HAS_JIT and self._lib == 'np' and is_jit_available():
+                return jit_compute_pt(self._px, self._py, self._lib)
+            else:
+                return compute_pt(self._px, self._py, self._lib)
+        return self._get_cached('pt', _compute_pt)
+        
     @property
     def p(self):
-        """Total momentum magnitude."""
-        return self._get_cached('p', 
-                             lambda: backend_sqrt(self._p_squared(), self._lib),
-                             ['px', 'py', 'pz'])
+        """Magnitude of the 3-momentum."""
+        def _compute_p():
+            # Use JIT-optimized function when available for NumPy arrays
+            if HAS_JIT and self._lib == 'np' and is_jit_available():
+                return jit_compute_p(self._px, self._py, self._pz, self._lib)
+            else:
+                return compute_p(self._px, self._py, self._pz, self._lib)
+        return self._get_cached('p', _compute_p)
     
     @property
     def mass(self):
-        """Invariant mass."""
-        def calc_mass():
-            # E² - p² calculation optimized with intermediate results
-            E_squared = self._E**2
-            p_squared = self._p_squared()
-            m_squared = E_squared - p_squared
-            
-            # Handle numerical precision issues (for physical particles, m² should be ≥ 0)
-            if isinstance(m_squared, (float, int)):
-                if m_squared < 0:
-                    if abs(m_squared) < 1e-10:  # Small tolerance for numerical precision
-                        m_squared = 0
-                    # Large negative values are a problem
-                    elif m_squared < -1e-6:
-                        import warnings
-                        warnings.warn(f"Mass calculation resulted in negative m²={m_squared}")
-                        m_squared = abs(m_squared)
+        """Invariant mass of the Lorentz vector."""
+        def _compute_mass():
+            # Use direct JIT-optimized mass calculation when available
+            # This avoids unnecessary intermediate calculations
+            if HAS_JIT and self._lib == 'np' and is_jit_available():
+                return jit_compute_mass(self._E, self._px, self._py, self._pz, self._lib)
             else:
-                # Array case
-                if self._lib == 'np':
-                    if (m_squared < 0).any():
-                        # Check for numerically significant negative values
-                        problematic = (m_squared < -1e-6)
-                        if problematic.any():
-                            import warnings
-                            warnings.warn("Mass calculation resulted in negative m² values")
-                        # Take absolute value to ensure physical result
-                        m_squared = np.where(m_squared < 0, abs(m_squared), m_squared)
-                elif self._lib == 'ak':
-                    import awkward as ak
-                    if ak.any(m_squared < 0):
-                        problematic = (m_squared < -1e-6)
-                        if ak.any(problematic):
-                            import warnings
-                            warnings.warn("Mass calculation resulted in negative m² values")
-                        # Take absolute value with the appropriate backend
-                        m_squared = ak.where(m_squared < 0, abs(m_squared), m_squared)
+                # Fall back to the standard computation
+                p = self.p  # This uses the cached p value
+                return compute_mass(self._E, p, self._lib)
             
-            return backend_sqrt(m_squared, self._lib)
-        
-        return self._get_cached('mass', calc_mass, ['px', 'py', 'pz', 'E'])
+        return self._get_cached('mass', _compute_mass)
     
     @property
     def phi(self):
-        """Azimuthal angle."""
-        return self._get_cached('phi',
-                             lambda: backend_atan2(self._py, self._px, self._lib),
-                             ['px', 'py'])
+        """Azimuthal angle in radians."""
+        def _compute_phi():
+            # Use JIT-optimized function when available for NumPy arrays
+            if HAS_JIT and self._lib == 'np' and is_jit_available():
+                return jit_compute_phi(self._px, self._py, self._lib)
+            else:
+                return compute_phi(self._px, self._py, self._lib)
+        return self._get_cached('phi', _compute_phi)
     
     @property
     def eta(self):
         """Pseudorapidity."""
-        def calc_eta():
-            p = self.p  # This will use our cached property
-            pz = self._pz
-            # Use a small epsilon to avoid division by zero
-            epsilon = 1e-10
-            
-            if self._lib == 'np':
-                # For NumPy, handle scalar and array cases
-                if isinstance(p, (float, int)) and isinstance(pz, (float, int)):
-                    if abs(p - abs(pz)) < epsilon:
-                        # When p ≈ |pz|, the particle is along the beam axis
-                        return float('inf') if pz >= 0 else float('-inf')
-                    # Numerically stable formula for eta
-                    return 0.5 * np.log((p + pz) / (p - pz))
-                else:
-                    # Array case with mask for special values
-                    near_beam = (abs(p - abs(pz)) < epsilon)
-                    # First compute the standard formula
-                    result = 0.5 * np.log((p + pz) / (p - pz + epsilon))
-                    # Then correct special cases (if any)
-                    if near_beam.any():
-                        result = np.where(near_beam & (pz >= 0), np.inf, result)
-                        result = np.where(near_beam & (pz < 0), -np.inf, result)
-                    return result
+        def _compute_eta():
+            # Use JIT-optimized function when available for NumPy arrays
+            if HAS_JIT and self._lib == 'np' and is_jit_available():
+                p = self.p  # Use cached value
+                return jit_compute_eta(p, self._pz, self._lib)
             else:
-                # For Awkward arrays
-                import awkward as ak
-                # Define a small epsilon to avoid division by zero
-                near_beam = (abs(p - abs(pz)) < epsilon)
-                # Compute the standard formula
-                result = 0.5 * backend_log((p + pz) / (p - pz + epsilon), self._lib)
-                # Correct special cases
-                if ak.any(near_beam):
-                    result = ak.where(near_beam & (pz >= 0), np.inf, result)
-                    result = ak.where(near_beam & (pz < 0), -np.inf, result)
-                return result
-                
-        return self._get_cached('eta', calc_eta, ['px', 'py', 'pz'])
+                p = self.p  # Use cached value
+                return compute_eta(p, self._pz, self._lib)
+        return self._get_cached('eta', _compute_eta)
     
     def __add__(self, other):
         """Add two Lorentz vectors."""
diff --git a/lvec/utils.py b/lvec/utils.py
index f775793..3e28fb5 100644
--- a/lvec/utils.py
+++ b/lvec/utils.py
@@ -6,6 +6,15 @@
                       backend_where)
 from .exceptions import ShapeError, BackendError
 
+# Import JIT-optimized functions if available
+try:
+    from .jit import (jit_compute_pt, jit_compute_p, jit_compute_mass,
+                     jit_compute_eta, jit_compute_phi, jit_compute_p4_from_ptepm,
+                     is_jit_available)
+    HAS_JIT = True
+except ImportError:
+    HAS_JIT = False
+
 def ensure_array(*args):
     """
     Convert inputs to consistent array type.
@@ -132,6 +141,121 @@ def check_shapes(*arrays):
                 shapes=shape_info
             )
 
+def _compute_pt_original(px, py, lib):
+    """
+    Original compute_pt implementation without JIT.
+    """
+    # Use the appropriate backend for the square root
+    pt = backend_sqrt(px*px + py*py, lib)
+    
+    # Convert to scalar only if both inputs were scalars
+    if isinstance(px, (float, int)) and isinstance(py, (float, int)):
+        return float(pt)
+    return pt
+
+def _compute_p_original(px, py, pz, lib):
+    """
+    Original compute_p implementation without JIT.
+    """
+    # Use the appropriate backend for the square root
+    p = backend_sqrt(px*px + py*py + pz*pz, lib)
+    
+    # Convert to scalar only if all inputs were scalars
+    if isinstance(px, (float, int)) and isinstance(py, (float, int)) and isinstance(pz, (float, int)):
+        return float(p)
+    return p
+
+def _compute_mass_original(E, p, lib):
+    """
+    Original compute_mass implementation without JIT.
+    """
+    # Calculate m² = E² - p²
+    m_squared = E*E - p*p
+    
+    # For proper physics, m² should be positive or zero
+    # Small negative values can occur due to numerical precision issues
+    negative_m_squared = False
+    
+    if isinstance(m_squared, (float, int)):
+        negative_m_squared = m_squared < 0
+        if negative_m_squared:
+            import warnings
+            warnings.warn(f"Negative mass-squared detected ({m_squared}), "
+                         "taking absolute value")
+    else:
+        if lib == 'np':
+            negative_m_squared = (m_squared < 0).any()
+            if negative_m_squared:
+                import warnings
+                warnings.warn(f"Negative mass-squared detected in array, "
+                             "taking absolute value")
+        elif lib == 'ak':
+            import awkward as ak
+            negative_m_squared = ak.any(m_squared < 0)
+            if negative_m_squared:
+                import warnings
+                warnings.warn(f"Negative mass-squared detected in array, "
+                             "taking absolute value")
+                
+    # Take the absolute value of m² if negative
+    if negative_m_squared:
+        if isinstance(m_squared, (float, int)):
+            m_squared = abs(m_squared)
+        elif lib == 'np':
+            m_squared = np.abs(m_squared)
+        elif lib == 'ak':
+            m_squared = ak.abs(m_squared)
+            
+    # Compute the square root
+    m = backend_sqrt(m_squared, lib)
+    
+    # Convert to scalar only if both inputs were scalars
+    if isinstance(E, (float, int)) and isinstance(p, (float, int)):
+        return float(m)
+    return m
+
+def _compute_eta_original(p, pz, lib):
+    """
+    Original compute_eta implementation without JIT.
+    """
+    # Small constant to avoid division by zero
+    epsilon = 1e-10
+    
+    # Compute eta using the numerically stable formula
+    eta = 0.5 * backend_log((p + pz) / (p - pz + epsilon), lib)
+    
+    # Convert to scalar only if both inputs were scalars
+    if isinstance(p, (float, int)) and isinstance(pz, (float, int)):
+        return float(eta)
+    return eta
+
+def _compute_phi_original(px, py, lib):
+    """
+    Original compute_phi implementation without JIT.
+    """
+    # Use the appropriate backend for arc-tangent
+    phi = backend_atan2(py, px, lib)
+    
+    # Convert to scalar only if both inputs were scalars
+    if isinstance(px, (float, int)) and isinstance(py, (float, int)):
+        return float(phi)
+    return phi
+
+def _compute_p4_from_ptepm_original(pt, eta, phi, m, lib):
+    """
+    Original compute_p4_from_ptepm implementation without JIT.
+    """
+    # Calculate Cartesian components
+    px = pt * backend_cos(phi, lib)
+    py = pt * backend_sin(phi, lib)
+    pz = pt * backend_sinh(eta, lib)
+    
+    # Calculate total momentum and energy
+    p = backend_sqrt(px*px + py*py + pz*pz, lib)
+    E = backend_sqrt(p*p + m*m, lib)
+    
+    return px, py, pz, E
+
 def compute_pt(px, py, lib):
     """
     Compute transverse momentum.
@@ -148,13 +272,12 @@ def compute_pt(px, py, lib):
     scalar or array-like
         Transverse momentum with the same type as input
     """
-    # Use the appropriate backend for the square root
-    pt = backend_sqrt(px*px + py*py, lib)
+    # Use JIT version if available and applicable
+    if HAS_JIT and lib == 'np':
+        return jit_compute_pt(px, py, lib)
     
-    # Convert to scalar only if both inputs were scalars
-    if isinstance(px, (float, int)) and isinstance(py, (float, int)):
-        return float(pt)
-    return pt
+    # Fall back to original implementation
+    return _compute_pt_original(px, py, lib)
 
 def compute_p(px, py, pz, lib):
     """
@@ -172,14 +295,12 @@ def compute_p(px, py, pz, lib):
     scalar or array-like
         Total momentum with the same type as input
     """
-    # Use the appropriate backend for the square root
-    p = backend_sqrt(px*px + py*py + pz*pz, lib)
+    # Use JIT version if available and applicable
+    if HAS_JIT and lib == 'np':
+        return jit_compute_p(px, py, pz, lib)
     
-    # Convert to scalar only if all inputs were scalars
-    if isinstance(px, (float, int)) and isinstance(py, (float, int)) and isinstance(pz, (float, int)):
-        return float(p)
-    return p
-
+    # Fall back to original implementation
+    return _compute_p_original(px, py, pz, lib)
 
 def compute_mass(E, p, lib):
     """
@@ -203,50 +324,11 @@ def compute_mass(E, p, lib):
     scalar or array-like
         Mass with the same type as input
     """
-    import warnings
-    m2 = E*E - p*p
-    
-    # Check for negative m² values
-    if isinstance(m2, (float, int)):
-        if m2 < 0:
-            warnings.warn(f"Negative m² value encountered: {m2}. Taking absolute value, but this may indicate unphysical particles or numerical issues.")
-            m2 = abs(m2)
-    else:
-        # For array inputs, check if any values are negative
-        has_negative = False
-        
-        if lib == 'ak' and HAS_AWKWARD:
-            # For Awkward arrays
-            try:
-                import awkward as ak
-                # Use ak.any for Awkward arrays
-                has_negative = ak.any(m2 < 0)
-            except Exception:
-                # Fallback if ak.any doesn't work
-                has_negative = any(m2 < 0)
-        else:
-            # For NumPy arrays
-            if hasattr(m2, 'any'):
-                has_negative = (m2 < 0).any()
-            else:
-                # For other iterable types
-                try:
-                    has_negative = any(m2 < 0)
-                except TypeError:
-                    # Not iterable
-                    has_negative = False
-            
-        if has_negative:
-            warnings.warn("Negative m² values encountered. Taking absolute values, but this may indicate unphysical particles or numerical issues.")
-            m2 = backend_where(m2 < 0, abs(m2), m2, lib)
+    # If we have JIT and full p4 components available in the caller, 
+    # the caller should use jit_compute_mass directly
     
-    m = backend_sqrt(m2, lib)
-    
-    # Convert to scalar only if all inputs were scalars
-    if isinstance(E, (float, int)) and isinstance(p, (float, int)):
-        return float(m)
-    return m
-
+    # Fall back to original implementation
+    return _compute_mass_original(E, p, lib)
 
 def compute_eta(p, pz, lib):
     """
@@ -270,14 +352,12 @@ def compute_eta(p, pz, lib):
     scalar or array-like
         Pseudorapidity with the same type as input
     """
-    epsilon = 1e-10
-    eta = 0.5 * backend_log((p + pz) / (p - pz + epsilon), lib)
+    # Use JIT version if available and applicable
+    if HAS_JIT and lib == 'np':
+        return jit_compute_eta(p, pz, lib)
     
-    # Convert to scalar only if all inputs were scalars
-    if isinstance(p, (float, int)) and isinstance(pz, (float, int)):
-        return float(eta)
-    return eta
-
+    # Fall back to original implementation
+    return _compute_eta_original(p, pz, lib)
 
 def compute_phi(px, py, lib):
     """
@@ -295,12 +375,12 @@ def compute_phi(px, py, lib):
     scalar or array-like
         Azimuthal angle with the same type as input
     """
-    phi = backend_atan2(py, px, lib)
+    # Use JIT version if available and applicable
+    if HAS_JIT and lib == 'np':
+        return jit_compute_phi(px, py, lib)
     
-    # Convert to scalar only if all inputs were scalars
-    if isinstance(px, (float, int)) and isinstance(py, (float, int)):
-        return float(phi)
-    return phi
+    # Fall back to original implementation
+    return _compute_phi_original(px, py, lib)
 
 def compute_p4_from_ptepm(pt, eta, phi, m, lib):
     """
@@ -324,12 +404,9 @@ def compute_p4_from_ptepm(pt, eta, phi, m, lib):
     tuple
         (px, py, pz, E) with the same type as input
     """
-    px = pt * backend_cos(phi, lib)
-    py = pt * backend_sin(phi, lib)
-    pz = pt * backend_sinh(eta, lib)
-    E = backend_sqrt(pt*pt * backend_cosh(eta, lib)**2 + m*m, lib)
+    # Use JIT version if available and applicable
+    if HAS_JIT and lib == 'np':
+        return jit_compute_p4_from_ptepm(pt, eta, phi, m, lib)
     
-    # Convert to scalar only if all inputs were scalars
-    if all(isinstance(x, (float, int)) for x in [pt, eta, phi, m]):
-        return float(px), float(py), float(pz), float(E)
-    return px, py, pz, E
\ No newline at end of file
+    # Fall back to original implementation
+    return _compute_p4_from_ptepm_original(pt, eta, phi, m, lib)
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 28a6475..42fdd47 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -30,6 +30,14 @@ dev = [
     "pytest>=7.0.0",
     "uproot>=5.0.0",
 ]
+performance = [
+    "numba>=0.61.2",
+]
+all = [
+    "numba>=0.61.2",
+    "pytest>=7.0.0",
+    "uproot>=5.0.0",
+]
 
 [project.urls]
 Homepage = "https://github.com/MohamedElashri/lvec"