Quick Start Guide
=================

This guide will help you get started with ixpepy for IXPE data analysis.

Basic Workflow
--------------

1. Download IXPE Data
~~~~~~~~~~~~~~~~~~~~~

Download observations by source name or observation ID:

.. code-block:: python

    from ixpepy.io.download_ixpe_data import download_source, download_obsid
    
    # Download all observations for a source
    obsids = download_source("Crab")
    
    # Download observation data
    download_obsid("01001099", destination="/path/to/data")

2. Process the data with ixpeobssim or HEASOFT
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This is usually done externally by using `ixpeobssim <https://ixpeobssim.readthedocs.io/en/latest/>`_ or tools provided
by HEASARC. With `ixpeobssim` the standard procedure consists in running
`xpselect` to select the source and background regions, and then `xpbin`
to bin the data with the `PHA1`, `PHA1Q` and `PHA1U` algorithms and
produce the Stokes spectra. 

If ixpeobssim is installed and available in the same environment as ixpepy,
the processing can happen in the same python script:

.. code-block:: python

    from ixpeobssim.bin import xpselect, xpbin
    from ixpeobssim.irf import DEFAULT_IRF_NAME

    selections = xpselect.xpselect(filelist=datafiles, ra=ra, dec=dec, 
                                   rad=rad)

    bin_pha1 = xpbin.xpbin(filelist=selections, algorithm='PHA1', 
        irfname=DEFAULT_IRF_NAME)
    bin_pha1u = xpbin.xpbin(filelist=selections, algorithm='PHA1U', 
        irfname=DEFAULT_IRF_NAME)
    bin_pha1q = xpbin.xpbin(filelist=selections, algorithm='PHA1Q', 
        irfname=DEFAULT_IRF_NAME)
    file_list  = bin_pha1 + bin_pha1u + bin_pha1q


3. Load the Stokes Spectra Data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Load the IXPE spectral data files for analysis:

.. code-block:: python

    from ixpepy.ixpeOGIPSpectrumLike import ixpeOGIPSpectrumLike
    
    # Initialize spectrum manager
    ixpe = ixpeOGIPSpectrumLike(energy_range='2-8')
    
    # Load files (9 files: 3 DUs × 3 Stokes parameters)
    source_name = 'source'
    file_list = []
    for du in [1, 2, 3]:  
        for stokes in ['', 'q', 'u']:
            file_list = os.path.join('.', '%s_det%s_pha1%s.fits' % (source_name, du, stokes))
            obs_list.append(file_list)
    ixpe.load_file_list(file_list)

4. Perform Spectral Analysis
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Use 3ML to fit a spectral model:

.. code-block:: python

    from threeML import *
    from astromodels import *
    
    # Define a spectral model
    spectrum = Powerlaw()
    polarization = StokesPolarization(Q=Constant(), U=Constant())
    source = PointSource(source_name, 0, 0, spectral_shape=spectrum,  polarization=polarization)
    model = Model(source)
    
    # Create joint likelihood and fit
    datalist = ixpe.get_datalist()
    jl = JointLikelihood(model, datalist)

    # Fit the model
    results = jl.fit()

5. Visualize Results
~~~~~~~~~~~~~~~~~~~~

Calculate polarization properties:

.. code-block:: python

    from ixpepy.utils.stokes import get_polarization_angle, get_polarization_degree

    u = results.get_variates('Crab.spectrum.main.polarization.U.Constant.k')
    q = results.get_variates('Crab.spectrum.main.polarization.Q.Constant.k')

    # Calculate polarization angle
    pol_angle = get_polarization_angle(q, u)
    print(f"Polarization angle: {pol_angle:.1f} degrees")
    
    # Calculate polarization degree
    pol_degree = get_polarization_degree(q, u)
    print(f"Polarization degree: {pol_degree:.2%}")

Plot the spectral fit and polarization:

.. code-block:: python

    from ixpepy.io.plotting import plot_I_spectrum, plot_QU_spectra, plot_QU
    
    # Plot I, Q, U spectra with model
    plot_I_spectrum(jk)
    plot_QU_spectra(jl)

    # Plot Q-U plane
    plot_QU(q, u)

For detailed examples and tutorials, see the :doc:`examples` section.

Advanced Features
-----------------

Rebinning Spectra
~~~~~~~~~~~~~~~~~

Rebin spectra to have minimum counts per bin:

.. code-block:: python

    ixpe.rebin_on_source(min_number_of_counts=25)

Effective Area Corrections
~~~~~~~~~~~~~~~~~~~~~~~~~~~

Fit effective area corrections for detector calibration:

.. code-block:: python

    # After creating JointLikelihood
    jl = JointLikelihood(model, ixpe.get_datalist())
    ixpe.fit_area_correction(model, min_value=0.8, max_value=1.2)
    jl.fit()