#define Calibration_cxx

#include <TH2.h>
#include <TF1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <fstream>
#include <utility>
#include <algorithm>
#include <TProfile.h>
#include <TVector3.h>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TGraphErrors.h"
#include "Calibration.h"

int padID = 0;

SX3 sx3_contr;
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
// TVector3 anodeIntersection;
std::map<int, std::pair<double, double>> slopeInterceptMap;

bool HitNonZero;
bool sx3ecut;
bool qqqEcut;

TH2F *hSX3FvsB;
TH2F *hSX3FvsB_g;
TH2F *hSX3;
TH1F *hZProj;
TH2F *hsx3IndexVE;
TH2F *hqqqIndexVE;
TH2F *hqqqIndexVE_cal;
TH2F *hsx3Coin;
TH2F *hqqqCoin;
TH2F *hqqqPolar;

TCutG *cut;
TCutG *cut1;

// Gain arrays

const int MAX_SX3 = 24;
const int MAX_UP = 4;
const int MAX_DOWN = 4;
const int MAX_BK = 4;
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double backGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool backGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
double frontGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
double uvdslope[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
double qqqGain[MAX_QQQ][MAX_BK][MAX_UP] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_BK][MAX_UP] = {{{false}}};

void Calibration::Begin(TTree * /*tree*/)
{
    TString option = GetOption();

    hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
    hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
    hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
    hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
    hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
    hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
    hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000);
    hqqqIndexVE_cal = new TH2F("hqqqIndexVE_cal", "QQQ index vs Energy (calibrated); QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000);
    hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
    hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4 * 2 * 16, 0, 4 * 2 * 16, 4 * 2 * 16, 0, 4 * 2 * 16);

    hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16 * 4, -TMath::Pi(), TMath::Pi(), 16, 10, 50);

    sx3_contr.ConstructGeo();
    pw_contr.ConstructGeo();
    // ----------------------- Load Back Gains
    {
        std::string filename = "sx3_GainMatchback.txt";
        std::ifstream infile(filename);
        if (!infile.is_open())
        {
            std::cerr << "Error opening " << filename << "!" << std::endl;
        }
        else
        {
            int id, bk, u, d;
            double gain;
            while (infile >> id >> bk >> u >> d >> gain)
            {
                backGain[id][bk][u][d] = gain;
                backGainValid[id][bk][u][d] = (gain > 0);
            }
            infile.close();
            std::cout << "Loaded back gains from " << filename << std::endl;
        }
    }

    // ----------------------- Load Front Gains
    {
        std::string filename = "sx3_GainMatchfront.txt";
        std::ifstream infile(filename);
        if (!infile.is_open())
        {
            std::cerr << "Error opening " << filename << "!" << std::endl;
        }
        else
        {
            int id, bk, u, d;
            double gain;
            while (infile >> id >> bk >> u >> d >> gain)
            {
                frontGain[id][bk][u][d] = gain;
                frontGainValid[id][bk][u][d] = (gain > 0);
            }
            infile.close();
            std::cout << "Loaded front gains from " << filename << std::endl;
        }
    }

    // ----------------------- Load QQQ Gains
    {
        std::string filename = "qqq_GainMatch.txt";
        std::ifstream infile(filename);
        if (!infile.is_open())
        {
            std::cerr << "Error opening " << filename << "!" << std::endl;
        }
        else
        {
            int det, ring, wedge;
            double gain;
            while (infile >> det >> ring >> wedge >> gain)
            {
                qqqGain[det][ring][wedge] = gain;
                qqqGainValid[det][ring][wedge] = (gain > 0);
            }
            infile.close();
            std::cout << "Loaded QQQ gains from " << filename << std::endl;
        }
    }

    SX3 sx3_contr;
}
Bool_t Calibration::Process(Long64_t entry)
{

    // if ( entry > 100 ) return kTRUE;

    hitPos.Clear();
    HitNonZero = false;

    // if( entry > 1) return kTRUE;
    // printf("################### ev : %llu \n", entry);

    b_sx3Multi->GetEntry(entry);
    b_sx3ID->GetEntry(entry);
    b_sx3Ch->GetEntry(entry);
    b_sx3E->GetEntry(entry);
    b_sx3T->GetEntry(entry);
    b_qqqMulti->GetEntry(entry);
    b_qqqID->GetEntry(entry);
    b_qqqCh->GetEntry(entry);
    b_qqqE->GetEntry(entry);
    b_qqqT->GetEntry(entry);
    b_pcMulti->GetEntry(entry);
    b_pcID->GetEntry(entry);

    sx3.CalIndex();
    qqq.CalIndex();
    pc.CalIndex();

    // sx3.Print();

    // ########################################################### Raw data
    //  //======================= SX3
    sx3ecut = false;
    std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
    for (int i = 0; i < sx3.multi; i++)
    {
        ID.push_back(std::pair<int, int>(sx3.id[i], i));
        hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);

        if (sx3.e[i] > 100)
        {
            sx3ecut = true;
        }

        for (int j = i + 1; j < sx3.multi; j++)
        {
            hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
        }
    }

    // --- safe SX3 handling (replace your existing block that builds sx3ID) ---
    if (ID.size() > 0)
    {
        std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
                  { return a.first < b.first; });

        std::vector<std::pair<int, int>> sx3ID;
        sx3ID.push_back(ID[0]);
        bool found = false;
        for (size_t i = 1; i < ID.size(); i++)
        {
            if (ID[i].first == sx3ID.back().first)
            {
                sx3ID.push_back(ID[i]);
                if (sx3ID.size() >= 3)
                {
                    found = true;
                }
            }
            else
            {
                if (!found)
                {
                    sx3ID.clear();
                    sx3ID.push_back(ID[i]);
                    found = false;
                }
            }
        }

        if (found)
        {
            // initialize to sentinel values
            int sx3ChUp = -1;
            int sx3ChDn = -1;
            int sx3ChBk = -1;
            float sx3EUp = 0.0f;
            float sx3EDn = 0.0f;

            // collect channels/energies
            for (size_t i = 0; i < sx3ID.size(); i++)
            {
                int index = sx3ID[i].second;
                int ch = sx3.ch[index];
                float e = sx3.e[index];

                if (ch < 8) // front channels
                {
                    // you used even/odd to denote down/up — keep that convention
                    if ((ch % 2) == 0) // down
                    {
                        sx3ChDn = ch;
                        sx3EDn = e;
                    }
                    else // up
                    {
                        sx3ChUp = ch;
                        sx3EUp = e;
                    }
                }
                else // back channels (assuming back channels are 8..11 or so)
                {
                    sx3ChBk = ch; // store as raw channel number; adapt if you index bk differently
                }
            }

            // Basic sanity checks before using indices:
            bool haveFrontPair = (sx3ChUp >= 0 && sx3ChDn >= 0);
            bool haveBack = (sx3ChBk >= 0);

            // convert raw channel numbers to array indices if needed:
            // example: if back channels are stored with ch>=8 and bk index is ch-8:
            int bk_index = (haveBack ? (sx3ChBk - 8) : -1);
            int up_index = (haveFrontPair ? sx3ChUp : -1);
            int dn_index = (haveFrontPair ? sx3ChDn : -1);

            // bounds checks against your array sizes
            auto sx3Id = sx3ID[0].first;
            auto inRange = [&](int id, int bk, int up, int dn)
            {
                if (id < 0 || id >= MAX_SX3)
                    return false;
                if (bk < 0 || bk >= MAX_BK)
                    return false;
                if (up < 0 || up >= MAX_UP)
                    return false;
                if (dn < 0 || dn >= MAX_DOWN)
                    return false;
                return true;
            };

            // apply calibration safely, fallback to raw ADC if no valid gain
            double calibEUp = sx3EUp;
            double calibEDn = sx3EDn;
            double calibEBack = 0.0;

            if (haveFrontPair && haveBack && inRange(sx3Id, bk_index, up_index, dn_index))
            {
                // If you stored front gains indexed by [id][bk][up][down]
                if (frontGainValid[sx3Id][bk_index][up_index][dn_index])
                {
                    calibEUp = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EUp;
                    // calibEDn = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EDn;
                }
                if (backGainValid[sx3Id][bk_index][up_index][dn_index])
                {
                    calibEBack = backGain[sx3Id][bk_index][up_index][dn_index] * sx3.e[sx3ID[0].second /* if you have back index */];
                    // Note: if back is from different index in sx3ID vector, you should find the exact index for the back hit
                }
            }
            else
            {
                // partial information: try best-effort per-channel checks
                if (haveFrontPair && up_index >= 0 && dn_index >= 0 && sx3Id >= 0 && sx3Id < MAX_SX3)
                {
                    if (frontGainValid[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)])
                    {
                        // attempt with default bk=0 if that makes sense in your geometry
                        calibEUp = frontGain[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)] * sx3EUp;
                        // calibEDn = frontGain[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)] * sx3EDn;
                    }
                }
                // keep calibEBack==0 if unavailable
            }

            // Only call CalSX3Pos if we have reasonable energies (avoid calling with zeros/uninitialized)
            if (haveFrontPair && (calibEUp > 0.0) && haveBack)
            {
                // find exact back energy value from sx3 entries if you tracked it above
                float backEnergyRaw = 0.0f;
                // locate the back index in sx3ID if needed
                for (size_t k = 0; k < sx3ID.size(); ++k)
                {
                    int idx = sx3ID[k].second;
                    if (sx3.ch[idx] >= 8)
                    {
                        backEnergyRaw = sx3.e[idx];
                        break;
                    }
                }
                // use calibrated back if available else raw
                double backEnergyToUse = (calibEBack > 0.0 ? calibEBack : backEnergyRaw);
                hSX3->Fill(sx3ChDn + 4, sx3ChBk);
                hSX3->Fill(sx3ChUp, sx3ChBk);

                // Fill the histogram for the front vs back
                hSX3FvsB->Fill(sx3EUp + sx3EDn, calibEBack);

                sx3_contr.CalSX3Pos(sx3Id, sx3ChUp, sx3ChDn, sx3ChBk, static_cast<float>(calibEUp), static_cast<float>(calibEDn));
                hitPos = sx3_contr.GetHitPos();
                HitNonZero = true;
            }
        } // found
    }

    // //======================= QQQ
    for (int i = 0; i < qqq.multi; i++)
    {

        int det = qqq.id[i]; // detector ID (0–3)
        int ch = qqq.ch[i];  // raw channel (0–31)

        // Separate ring vs wedge channel
        int ring = -1;
        int wedge = -1;
        if (ch < 16)
        { // wedge
            wedge = ch;
        }
        else
        { // ring
            ring = ch - 16;
        }

        double Ecal = qqq.e[i]; // default = raw
        if (ring >= 0 && wedge >= 0 && qqqGainValid[det][ring][wedge])
        {
            Ecal *= qqqGain[det][ring][wedge];
        }
        // for( int j = 0; j < pc.multi; j++){
        // if(pc.index[j]==4){
        hqqqIndexVE_cal->Fill(qqq.index[i], Ecal);
        hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]);

        // }
        // printf("QQQ ID : %d, ch : %d, e : %d \n", qqq.id[i], qqq.ch[i], qqq.e[i]);
        if (qqq.e[i] > 100)
        {
            qqqEcut = true;
        }
        // }
        for (int j = 0; j < qqq.multi; j++)
        {
            if (j == i)
                continue;
            hqqqCoin->Fill(qqq.index[i], qqq.index[j]);
        }

        // }

        for (int j = i + 1; j < qqq.multi; j++)
        {
            // if( qqq.used[i] == true ) continue;

            // if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0  ){ // must be same detector and wedge and ring
            if (qqq.id[i] == qqq.id[j])
            { // must be same detector

                int chWedge = -1;
                int chRing = -1;
                if (qqq.ch[i] < qqq.ch[j])
                {
                    chRing = qqq.ch[j] - 16;
                    chWedge = qqq.ch[i];
                }
                else
                {
                    chRing = qqq.ch[i];
                    chWedge = qqq.ch[j] - 16;
                }
                // printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);

                double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
                double rho = 50. + 40. / 16. * (chRing + 0.5);
                // if(qqq.e[i]>50){
                hqqqPolar->Fill(theta, rho);
                // }
                // qqq.used[i] = true;
                // qqq.used[j] = true;

                if (!HitNonZero)
                {
                    double x = rho * TMath::Cos(theta);
                    double y = rho * TMath::Sin(theta);
                    hitPos.SetXYZ(x, y, 23 + 75 + 30);
                    HitNonZero = true;
                }
            }
        }
    }

    return kTRUE;
}
void Calibration::Terminate()
{
}