#include "honlinemonhist.h"
#include "honlinehistarray.h"
#include "honlinetrendhist.h"

#include "hcategory.h"
#include "rpcdef.h"
#include "hrpcraw.h"


#include "fwdetdef.h"
#include "hstsraw.h"
#include "hstscal.h"
#include "hstsgeompar.h"
#include "hfrpcraw.h"
#include "hfrpccal.h"
#include "hfrpchit.h"
#include "hforwardcand.h"
#include "hparticletool.h"


#include "helpers.h"

#include "TList.h"
#include "TString.h"

#include <map>
#include <iostream>
using namespace std;

constexpr Float_t z_target = -124.; // mean target position

constexpr Float_t pi_value = 3.14159265359; // Pi constant value
constexpr Float_t c_value = 299.792458; // Speed of light (In mm/ns)
constexpr Float_t proton_mass = 938.272081358; // Proton mass [MeV/c^2]

//constexpr Float_t inverse_gamma2 = 0.29429; //value corresponding to 4.5 GeV beam energy
constexpr Float_t inverse_gamma2 = 0.5429;    //value corresponding to 1.58 GeV beam energy

// Cuts for elastic scattering
constexpr Float_t cut_HF_phi_diff_min = 175.0;
//constexpr Float_t cut_HF_phi_diff_max = 184.0;    //value from 4.2 GeV
constexpr Float_t cut_HF_phi_diff_max = 185.0;

// for 4.5  GeV +/- 1.5 deg
// for 1.58 GeV +/- 2.0 deg
constexpr Float_t cut_HH_phi_diff_min = 178.0;
constexpr Float_t cut_HH_phi_diff_max = 182.0;

constexpr Float_t cut_HF_tan_theta_min = inverse_gamma2 - 0.1;
constexpr Float_t cut_HF_tan_theta_max = inverse_gamma2 + 0.1;

// for 4.5  GeV +/- 0.015
// for 1.58 GeV +/- 0.025
constexpr Float_t cut_HH_tan_theta_min = inverse_gamma2 - 0.025;
constexpr Float_t cut_HH_tan_theta_max = inverse_gamma2 + 0.025;

map < TString , HOnlineMonHistAddon* > fphysicsMap;

//----------------------------------------------------------------------------------------------
// Initialize histograms
//----------------------------------------------------------------------------------------------

struct fphysicsP
{
    // 2D histograms tables (straw_module:straw_layer)
    HOnlineMonHistAddon* h_HH_and_HF_theta_P1_vs_theta_P2 = 0;
    HOnlineMonHistAddon* h_HH_and_HF_theta_P1_vs_theta_P2_el = 0;
    HOnlineMonHistAddon* h_HH_and_HF_mom_vs_theta = 0;
    HOnlineMonHistAddon* h_HH_and_HF_mom_vs_theta_el = 0;
    HOnlineMonHistAddon* h_track_start_Y_vs_X = 0;
    HOnlineMonHistAddon* h_track_start_R_vs_Z = 0;

    // 1D histograms tables (straw_module)
    HOnlineMonHistAddon* h_HF_phi_diff = 0;
    HOnlineMonHistAddon* h_HF_phi_diff_cut = 0;
    HOnlineMonHistAddon* h_HF_phi_diff_closeup = 0;
    HOnlineMonHistAddon* h_HF_phi_diff_cut_closeup = 0;
    HOnlineMonHistAddon* h_HF_tan_theta_prod = 0;
    HOnlineMonHistAddon* h_HF_tan_theta_prod_cut = 0;
    HOnlineMonHistAddon* h_HH_phi_diff = 0;
    HOnlineMonHistAddon* h_HH_phi_diff_cut = 0;
    HOnlineMonHistAddon* h_HH_phi_diff_closeup = 0;
    HOnlineMonHistAddon* h_HH_phi_diff_cut_closeup = 0;
    HOnlineMonHistAddon* h_HH_tan_theta_prod = 0;
    HOnlineMonHistAddon* h_HH_tan_theta_prod_cut = 0;
    HOnlineMonHistAddon* h_track_start_Z = 0;

    HOnlineMonHistAddon* h_track_target_proj = 0;
    HOnlineMonHistAddon* h_track_xy = 0;
    HOnlineMonHistAddon* h_track_txty = 0;

    // Trend histograms
    HOnlineMonHistAddon* h_HF_elastics_trend = 0;
    HOnlineMonHistAddon* h_HH_elastics_trend = 0;

};

fphysicsP* pFPhysics;

//----------------------------------------------------------------------------------------------
// Create histograms
//----------------------------------------------------------------------------------------------

// Array of indices for short straws

Bool_t createHistFPhysics(TList& histpool)
{
    mapHolder::setMap(fphysicsMap); // make fphysicsMap currentMap

    //----------------------------------------------------------------------------------------------
    // Defining histogram parameters
    //----------------------------------------------------------------------------------------------

    const Char_t* hists[] =
    {
        //----------------------------------------------------------------------------------------------
        // 2D histograms tables (straw_module:straw_layer)
        "FORMAT#mon TYPE#2F NAME#h_HH_and_HF_theta_P1_vs_theta_P2 TITLE#HH_HF_ThetaVsTheta ACTIVE#1 RESET#1 REFRESH#100000 BIN#90:0.0:90.0:90:0.0:90.0 SIZE#0:0 AXIS#theta_2_[deg]:theta_1_[deg]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#2F NAME#h_HH_and_HF_theta_P1_vs_theta_P2_el TITLE#HH_HF_ThetaVsTheta_EL ACTIVE#1 RESET#1 REFRESH#100000 BIN#90:0.0:90.0:90:0.0:90.0 SIZE#0:0 AXIS#theta_2_[deg]:theta_1_[deg]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#array TYPE#2F NAME#h_HH_and_HF_mom_vs_theta TITLE#HH_HF_mom_vs_theta ACTIVE#1 RESET#1 REFRESH#100000 BIN#90:0.0:90.0:175:0.0:7000.0 SIZE#1:6 AXIS#:theta_[deg]:Mom_[MeV]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#array TYPE#2F NAME#h_HH_and_HF_mom_vs_theta_el TITLE#HH_HF_mom_vs_theta_EL ACTIVE#1 RESET#1 REFRESH#100000 BIN#90:0.0:90.0:175:0.0:7000.0 SIZE#1:6 AXIS#theta_[deg]:Mom_[MeV]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",

        "FORMAT#mon TYPE#2F NAME#h_track_start_Y_vs_X TITLE#Track_start_Y_vs_X ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:-100.0:100.0:100:-100.0:100.0 SIZE#0:0 AXIS#X_[mm]:Y_[mm]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#0:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#2F NAME#h_track_start_R_vs_Z TITLE#Track_start_R_vs_Z ACTIVE#1 RESET#1 REFRESH#100000 BIN#300:-3000.0:3000.0:100:-100.0:100.0 SIZE#0:0 AXIS#Z_[mm]:X_[mm]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#0:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",

        //----------------------------------------------------------------------------------------------
        // 1D histograms tables (straw_module)
        "FORMAT#mon TYPE#1F NAME#h_HF_phi_diff TITLE#HF_Phi_Diff ACTIVE#1 RESET#1 REFRESH#100000 BIN#720:0.0:360.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HF_phi_diff_cut TITLE#HF_Phi_Diff_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#720:0.0:360.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HF_phi_diff_closeup TITLE#HF_Phi_Diff ACTIVE#1 RESET#1 REFRESH#100000 BIN#120:150.0:210.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HF_phi_diff_cut_closeup TITLE#HF_Phi_Diff_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#120:150.0:210.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HF_tan_theta_prod TITLE#HF_Tan_Theta_Prod ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:0.0:1.0:0:0.0:0.0 SIZE#0:0 AXIS#tan(theta_1)*tan(theta_2):counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HF_tan_theta_prod_cut TITLE#HF_Tan_Theta_Prod_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:0.0:1.0:0:0.0:0.0 SIZE#0:0 AXIS#tan(theta_1)*tan(theta_2):counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",

        "FORMAT#mon TYPE#1F NAME#h_HH_phi_diff TITLE#HH_Phi_Diff ACTIVE#1 RESET#1 REFRESH#100000 BIN#720:0.0:360.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HH_phi_diff_cut TITLE#HH_Phi_Diff_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#720:0.0:360.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HH_phi_diff_closeup TITLE#HH_Phi_Diff ACTIVE#1 RESET#1 REFRESH#100000 BIN#120:150.0:210.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HH_phi_diff_cut_closeup TITLE#HH_Phi_Diff_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#120:150.0:210.0:0:0.0:0.0 SIZE#0:0 AXIS#phi_[deg]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HH_tan_theta_prod TITLE#HH_Tan_Theta_Prod ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:0.0:1.0:0:0.0:0.0 SIZE#0:0 AXIS#tan(theta_1)*tan(theta_2):counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#1F NAME#h_HH_tan_theta_prod_cut TITLE#HH_Tan_Theta_Prod_CUT ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:0.0:1.0:0:0.0:0.0 SIZE#0:0 AXIS#tan(theta_1)*tan(theta_2):counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",

        "FORMAT#mon TYPE#1F NAME#h_track_start_Z TITLE#Track_start_Z ACTIVE#1 RESET#1 REFRESH#100000 BIN#100:-1000.0:1000.0:0:0.0:0.0 SIZE#0:0 AXIS#Z_[mm]:counts:no DIR#no OPT#no STATS#0 LOG#0:0:0 GRID#1:1 LINE#1:0 FILL#38:0 MARKER#0:0:0 RANGE#-99:-99",

        "FORMAT#mon TYPE#2F NAME#h_track_target_proj TITLE#Track_on_target_projection ACTIVE#1 RESET#1 REFRESH#50000 BIN#100:-50.0:50.0:100:-50.0:50.0 SIZE#0:0 AXIS#X_[mm]:Y_[mm]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#0:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#2F NAME#h_track_xy TITLE#Track_Base_XY ACTIVE#1 RESET#1 REFRESH#50000 BIN#100:-400.0:400.0:100:-400.0:400.0 SIZE#0:0 AXIS#X_[mm]:Y_[mm]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#0:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",
        "FORMAT#mon TYPE#2F NAME#h_track_txty TITLE#Track_Dir_TxTy ACTIVE#1 RESET#1 REFRESH#50000 BIN#100:-0.2:0.2:100:-0.2:0.2 SIZE#0:0 AXIS#X_[mm]:Y_[mm]:no DIR#no OPT#colz STATS#0 LOG#0:0:0 GRID#1:1 LINE#0:0 FILL#0:0 MARKER#0:0:0 RANGE#-99:-99",

        // Trend histograms
        "FORMAT#trend TYPE#1F NAME#h_HF_elastics_trend TITLE#HF_Elastics_Trend ACTIVE#1 RESET#0 REFRESH#100000 BIN#50:0.0:50.0:0:0.:0. SIZE#0:0 AXIS#trend:el_counts_per_100k_ev:no DIR#no OPT#lp STATS#0 LOG#0:0:0 GRID#0:1 LINE#1:0 FILL#0:0 MARKER#1:20:0.5 RANGE#-99:-99",
        "FORMAT#trend TYPE#1F NAME#h_HH_elastics_trend TITLE#HH_Elastics_Trend ACTIVE#1 RESET#0 REFRESH#100000 BIN#50:0.0:50.0:0:0.:0. SIZE#0:0 AXIS#trend:el_counts_per_100k_ev:no DIR#no OPT#lp STATS#0 LOG#0:0:0 GRID#0:1 LINE#1:0 FILL#0:0 MARKER#1:20:0.5 RANGE#-99:-99",
    };

    //----------------------------------------------------------------------------------------------
    // Create histograms and add them to the pool
    //----------------------------------------------------------------------------------------------

    mapHolder::createHists(sizeof(hists)/sizeof(Char_t*),hists,histpool);

    //
    pFPhysics = new fphysicsP();

    // 2D histograms tables (straw_module:straw_layer)
    pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2 = get("h_HH_and_HF_theta_P1_vs_theta_P2");
    pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2_el = get("h_HH_and_HF_theta_P1_vs_theta_P2_el");
    pFPhysics->h_HH_and_HF_mom_vs_theta = get("h_HH_and_HF_mom_vs_theta");
    pFPhysics->h_HH_and_HF_mom_vs_theta_el = get("h_HH_and_HF_mom_vs_theta_el");
    pFPhysics->h_track_start_Y_vs_X = get("h_track_start_Y_vs_X");
    pFPhysics->h_track_start_R_vs_Z = get("h_track_start_R_vs_Z");

    // 1D histograms tables (straw_module)
    pFPhysics->h_HF_phi_diff = get("h_HF_phi_diff");
    pFPhysics->h_HF_phi_diff_cut = get("h_HF_phi_diff_cut");
    pFPhysics->h_HF_phi_diff_closeup = get("h_HF_phi_diff_closeup");
    pFPhysics->h_HF_phi_diff_cut_closeup = get("h_HF_phi_diff_cut_closeup");
    pFPhysics->h_HF_tan_theta_prod = get("h_HF_tan_theta_prod");
    pFPhysics->h_HF_tan_theta_prod_cut = get("h_HF_tan_theta_prod_cut");
    pFPhysics->h_HH_phi_diff = get("h_HH_phi_diff");
    pFPhysics->h_HH_phi_diff_cut = get("h_HH_phi_diff_cut");
    pFPhysics->h_HH_phi_diff_closeup = get("h_HH_phi_diff_closeup");
    pFPhysics->h_HH_phi_diff_cut_closeup = get("h_HH_phi_diff_cut_closeup");
    pFPhysics->h_HH_tan_theta_prod = get("h_HH_tan_theta_prod");
    pFPhysics->h_HH_tan_theta_prod_cut = get("h_HH_tan_theta_prod_cut");
    pFPhysics->h_track_start_Z = get("h_track_start_Z");

    // tracking
    pFPhysics->h_track_target_proj = get("h_track_target_proj");
    pFPhysics->h_track_xy = get("h_track_xy");
    pFPhysics->h_track_txty = get("h_track_txty");

    // Trend histograms
    pFPhysics->h_HF_elastics_trend = get("h_HF_elastics_trend");
    pFPhysics->h_HH_elastics_trend = get("h_HH_elastics_trend");


    //----------------------------------------------------------------------------------------------
    // Setting colors for histogram arrays
    //----------------------------------------------------------------------------------------------

    // Colors
    HOnlineMonHistAddon* addon = 0;
    if( (addon = pFPhysics->h_HF_elastics_trend))
    {
        addon->getP()->SetLineColor(2);
        addon->getP()->SetLineWidth(2);
    }

    if( (addon = pFPhysics->h_HH_elastics_trend))
    {
        addon->getP()->SetLineColor(2);
        addon->getP()->SetLineWidth(2);
    }

    //----------------------------------------------------------------------------------------------
    // Accessing geometry container
    //----------------------------------------------------------------------------------------------

    HRuntimeDb *rtdb = HRuntimeDb::instance();//myHades -> getRuntimeDb();
    rtdb->initContainers(1);
    rtdb->print();

    return kTRUE;
}

//----------------------------------------------------------------------------------------------
// Main part of the program
//----------------------------------------------------------------------------------------------

Bool_t fillHistFPhysics(Int_t evtCt)
{
    mapHolder::setMap(fphysicsMap); // make fphysicsMap currentMap

    //----------------------------------------------------------------------------------------------
    // Categories to access
    //----------------------------------------------------------------------------------------------

    HCategory* fForwardCand = gHades->getCurrentEvent()->getCategory(catForwardCand);
    HCategory* fParticleCand = gHades->getCurrentEvent()->getCategory(catParticleCand);

    // Accessing trigger bit
    Int_t TBit = gHades->getCurrentEvent()->getHeader()->getTBit();

    Int_t particle_cand_cnt = 0;
    Int_t forward_cand_cnt = 0;

    Int_t /*t_PT1 = 0,*/ t_PT2 = 0/*, t_PT3 = 0*/;
//     Int_t t_ToF_Rec;
    Int_t t_Sector_P1, t_Sector_P2;//, t_Sector_Pair;
//     Float_t t_z_start_P1;//, t_track_length_P1;
    Float_t t_z_start_P2;//, t_track_length_P2;
    Float_t t_theta_P1, t_phi_P1, t_mom_P1;//, t_tof_P1, t_tof_theor_P1, t_mom_P1, t_mom_theor_P1, t_beta_P1, t_beta_theor_P1, t_E_P1, t_E_theor_P1;
    Float_t t_theta_P2, t_phi_P2, t_mom_P2;//, t_tof_P2, t_tof_theor_P2, t_mom_P2, t_mom_theor_P2, t_beta_P2, t_beta_theor_P2, t_E_P2, t_E_theor_P2;
    Float_t t_phi_diff, t_tan_theta_product;

    // Temp values to make sure there's only one elastic reaction per event - if the values are not zero (an elastic reaction has already been found in this event) that means the event should be skipped without filling el. histograms
    Int_t el_Sector_P1 = 0, el_Sector_P2 = 0;
    Float_t el_theta_P1 = 0.0, /*el_phi_P1 = 0.0,*/ el_mom_P1 = 0.0;
    Float_t el_theta_P2 = 0.0, /*el_phi_P2 = 0.0,*/ el_mom_P2 = 0.0;
//     Float_t el_phi_diff = 0.0, el_tan_theta_product = 0.0;

//     Float_t bkg_phi_diff = 0.0;

    // Check if event contains HH or HF pair
    int has_HH_cnt = 0;
    int has_HF_cnt = 0;

    // Check if event contains an elastic pair
    bool has_HH_elastic = false;
    bool has_HF_elastic = false;

    // Check if event a background reference pair
    bool has_HH_bkg_ref = false;
    bool has_HF_bkg_ref = false;

    // Elastics counters
    static Int_t num_HF_el_peak = 0;
    static Int_t num_HH_el_peak = 0;
    static Int_t num_HF_bkg_ref = 0;
    static Int_t num_HH_bkg_ref = 0;

    if (fForwardCand)
    {
        Int_t vcnt = fForwardCand->getEntries();

        for (Int_t j = 0; j < vcnt; ++j)
        {
            HForwardCand* fwdetvec = HCategoryManager::getObject(fwdetvec, fForwardCand, j);

            Float_t z = fwdetvec->getBaseZ();
            if (pFPhysics->h_track_target_proj)
                pFPhysics->h_track_target_proj->getP()->Fill(
                                                             fwdetvec->getBaseX() - (z - z_target) * fwdetvec->getDirTx(),
                                                             fwdetvec->getBaseY() - (z - z_target) * fwdetvec->getDirTy());

            if (pFPhysics->h_track_xy)
                pFPhysics->h_track_xy->getP()->Fill(fwdetvec->getBaseX(), fwdetvec->getBaseY());
            if (pFPhysics->h_track_txty)
                pFPhysics->h_track_txty->getP()->Fill(fwdetvec->getDirTx(), fwdetvec->getDirTy());
        }
    }

//     if((TBit&2048)==2048)
//     {
//         t_PT1 = 1;
//     }
    if((TBit&4096)==4096)
    {
        t_PT2 = 1;
    }
//     if((TBit&8192)==8192)
//     {
//         t_PT3 = 1;
//     }

    //printf("FORWARD Containers: %#x  %#x\n", fParticleCand, fForwardCand);
    if(fForwardCand && fParticleCand && t_PT2 == 1)
    {
        particle_cand_cnt = fParticleCand->getEntries();
        forward_cand_cnt = fForwardCand->getEntries();

        if(particle_cand_cnt >= 2 && forward_cand_cnt >= 2)
        {
            HGeomVector av_PCA_fcand, av_PCA_pcand;
            av_PCA_fcand.clear();
            av_PCA_pcand.clear();


            // Vector of "HGeomVector" paris defining tracks
            std::vector<std::pair<HGeomVector, HGeomVector>> tracks_pcand;
            std::vector<std::pair<HGeomVector, HGeomVector>> tracks_fcand;

            //----------------------------------------------------------------------------------------------
            // FwDet Particles
            for (Int_t i = 0; i < forward_cand_cnt; i++)
            {
                HForwardCand* fwdetcand = 0;
                fwdetcand = HCategoryManager::getObject(fwdetcand, fForwardCand, i);

                HGeomVector base;
                HGeomVector dir;

                // Track extraction
                HParticleTool::calcSegVector(fwdetcand->getZ(), fwdetcand->getR(), fwdetcand->Phi(), fwdetcand->Theta(), base, dir);
                tracks_fcand.push_back({base, dir}); // Adding track to vector
            }

            //----------------------------------------------------------------------------------------------
            // HADES Particles
            for (Int_t i = 0; i < particle_cand_cnt; i++)
            {
                HParticleCand* fparticlecand = 0;
                fparticlecand = HCategoryManager::getObject(fparticlecand, fParticleCand, i);

                HGeomVector base;
                HGeomVector dir;

                // Track extraction
                HParticleTool::calcSegVector(fparticlecand->getZ(), fparticlecand->getR(), fparticlecand->Phi(), fparticlecand->Theta(), base, dir);
                tracks_pcand.push_back({base, dir}); // Adding track to vector
            }

            //----------------------------------------------------------------------------------------------
            // Calculating closest approach point - average from pairing tracks with each other in ForwardCand

            // Number of tracks
            auto size_fcand = tracks_fcand.size();
            Double_t pairs_num_fcand = 0.0;

            // Loop over "i"
            for (UInt_t i = 0; i < size_fcand; i++)
            {
                // Loop over "j" - starting from "i+1" to avoid pairing same tracks twice and pairing a track with itself
                for (UInt_t j = i + 1; j < size_fcand; j++)
                {
                    //auto [base1, dir1] = tracks[i]; c++17 only
                    //auto [base2, dir2] = tracks[j];
                    auto base1 = tracks_fcand[i].first;
                    auto dir1 = tracks_fcand[i].second;
                    auto base2 = tracks_fcand[j].first;
                    auto dir2 = tracks_fcand[j].second;

                    // Calculating the point of closest approach between the two tracks
                    HGeomVector PCA_fcand ;
                    PCA_fcand.setXYZ(-2000,-2000,-2000);
                    if(dir1==dir2 && dir1.getZ()==1){ // track same direction along z axis
                        //cout<<"parallel fw"<<endl;
                    } else {
                        PCA_fcand = HParticleTool::calculatePointOfClosestApproach(base1, dir1, base2, dir2);
                    }

                    // Calculating average PCA
                    av_PCA_fcand += PCA_fcand;
                    pairs_num_fcand = pairs_num_fcand + 1.0;
                }
            }

            // Number of tracks
            auto size_pcand = tracks_pcand.size();
            Double_t pairs_num_pcand = 0.0;

            // Loop over "i"
            for (UInt_t i = 0; i < size_pcand; i++)
            {
                // Loop over "j" - starting from "i+1" to avoid pairing same tracks twice and pairing a track with itself
                for (UInt_t j = i + 1; j < size_pcand; j++)
                {
                    //auto [base1, dir1] = tracks[i]; c++17 only
                    //auto [base2, dir2] = tracks[j];
                    auto base1 = tracks_pcand[i].first;
                    auto dir1 = tracks_pcand[i].second;
                    auto base2 = tracks_pcand[j].first;
                    auto dir2 = tracks_pcand[j].second;

                    // Calculating the point of closest approach between the two tracks
                    HGeomVector PCA_pcand ;
                    PCA_pcand.setXYZ(-2000,-2000,-2000);
                    if(dir1==dir2 ){ // trach same direction
                        //cout<<"parallel"<<endl;
                    } else {
                        PCA_pcand = HParticleTool::calculatePointOfClosestApproach(base1, dir1, base2, dir2);
                    }

                    // Calculating average PCA
                    av_PCA_pcand += PCA_pcand;
                    pairs_num_pcand = pairs_num_pcand + 1.0;
                }
            }

            // Calculating average PCAs
            av_PCA_fcand /= pairs_num_fcand;
            av_PCA_pcand /= pairs_num_pcand;
            if(pFPhysics->h_track_start_Y_vs_X) pFPhysics->h_track_start_Y_vs_X->getP()->Fill(av_PCA_fcand.getX(),av_PCA_fcand.getY());
            if(pFPhysics->h_track_start_R_vs_Z) pFPhysics->h_track_start_R_vs_Z->getP()->Fill(av_PCA_pcand.getZ(),TMath::Sqrt(av_PCA_fcand.getX()*av_PCA_fcand.getX() + av_PCA_fcand.getY()*av_PCA_fcand.getY() ));
            if(pFPhysics->h_track_start_Z) pFPhysics->h_track_start_Z->getP()->Fill(av_PCA_pcand.getZ());
        }
        else
        {
            //----------------------------------------------------------------------------------------------
            // Loops over all tracks
            //----------------------------------------------------------------------------------------------

            HGeomVector av_PCA;
            av_PCA.clear();

            // Vector of "HGeomVector" paris defining tracks
            std::vector<std::pair<HGeomVector, HGeomVector>> tracks;

            //----------------------------------------------------------------------------------------------
            // FwDet Particles
            for (Int_t i = 0; i < forward_cand_cnt; i++)
            {
                HForwardCand* fwdetcand = 0;
                fwdetcand = HCategoryManager::getObject(fwdetcand, fForwardCand, i);

                HGeomVector base;
                HGeomVector dir;

                // Track extraction
                HParticleTool::calcSegVector(fwdetcand->getZ(), fwdetcand->getR(), fwdetcand->Phi(), fwdetcand->Theta(), base, dir);
                tracks.push_back({base, dir}); // Adding track to vector
            }

            //----------------------------------------------------------------------------------------------
            // HADES Particles
            for (Int_t i = 0; i < particle_cand_cnt; i++)
            {
                HParticleCand* fparticlecand = 0;
                fparticlecand = HCategoryManager::getObject(fparticlecand, fParticleCand, i);

                HGeomVector base;
                HGeomVector dir;

                // Track extraction
                HParticleTool::calcSegVector(fparticlecand->getZ(), fparticlecand->getR(), fparticlecand->Phi(), fparticlecand->Theta(), base, dir);
                tracks.push_back({base, dir}); // Adding track to vector
            }

            //----------------------------------------------------------------------------------------------
            // Calculating closest approach point - average from pairing tracks with each other

            // Number of tracks
            auto s = (Int_t)tracks.size();
            Double_t pairs_num = 0.0;

            // Loop over "i"
            for (Int_t i = 0; i < s; i++)
            {
                // Loop over "j" - starting from "i+1" to avoid pairing same tracks twice and pairing a track with itself
                for (Int_t j = i + 1; j < s; j++)
                {
                    //auto [base1, dir1] = tracks[i]; c++17 only
                    //auto [base2, dir2] = tracks[j];
                    auto base1 = tracks[i].first;
                    auto dir1 = tracks[i].second;
                    auto base2 = tracks[j].first;
                    auto dir2 = tracks[j].second;

                    // Calculating the point of closest approach between the two tracks
                    HGeomVector PCA;
                    PCA.setXYZ(-2000,-2000,-2000);
                    if(dir1==dir2 ){ // trach same direction
                        //cout<<"parallel"<<endl;
                    } else {
                        PCA = HParticleTool::calculatePointOfClosestApproach(base1, dir1, base2, dir2);
                    }

                    // Calculating average PCA
                    av_PCA += PCA;
                    pairs_num = pairs_num + 1.0;
                }
            }

            // Calculating average PCA
            av_PCA /= pairs_num;
            if(pFPhysics->h_track_start_Y_vs_X) pFPhysics->h_track_start_Y_vs_X->getP()->Fill(av_PCA.getX(),av_PCA.getY());
            if(pFPhysics->h_track_start_R_vs_Z) pFPhysics->h_track_start_R_vs_Z->getP()->Fill(av_PCA.getZ(),av_PCA.getX());
            if(pFPhysics->h_track_start_Z) pFPhysics->h_track_start_Z->getP()->Fill(av_PCA.getZ());
        }


        Float_t peak_width = cut_HH_phi_diff_max - cut_HH_phi_diff_min;
        Float_t integration_offset = 0.5*peak_width;

        Float_t bkg_low_min = (cut_HH_phi_diff_min - 0.5*peak_width - integration_offset);
        Float_t bkg_low_max = (cut_HH_phi_diff_min - integration_offset);
        Float_t bkg_high_min = (cut_HH_phi_diff_max + integration_offset);
        Float_t bkg_high_max = (cut_HH_phi_diff_max + 0.5*peak_width + integration_offset);

        // Loop over ParticleCand
        for (Int_t j = 0; j < particle_cand_cnt; j++)
        {
            //----------------------------------------------------------------------------------------------
            // HADES - HADES case
            //----------------------------------------------------------------------------------------------
            
            // Loop over ParticleCand - starting from "j+1" to avoid pairing same tracks twice and pairing a track with itself
            for (Int_t k = j + 1; k < particle_cand_cnt; k++)
            {
                //----------------------------------------------------------------------------------------------
                // Particle 1 - HADES particle
                HParticleCand* fparticlecand = 0;
                fparticlecand = HCategoryManager::getObject(fparticlecand, fParticleCand, j);

                if ( !( (fparticlecand->getTofdEdx() > 1.0 || fparticlecand->getMdcdEdx() > 1.0) && fparticlecand->isFlagBit(kIsUsed))) continue;

                t_Sector_P1 = (Int_t)fparticlecand->getSector();
                t_phi_P1 = fparticlecand->getPhi(); // [deg]
                t_theta_P1 = fparticlecand->getTheta(); // [deg]
                t_mom_P1 = fparticlecand->getMomentum(); // [MeV/c]

                //----------------------------------------------------------------------------------------------
                // Particle 2 - HADES particle
                fparticlecand = HCategoryManager::getObject(fparticlecand, fParticleCand, k);

                if ( !( (fparticlecand->getTofdEdx() > 1.0 || fparticlecand->getMdcdEdx() > 1.0) && fparticlecand->isFlagBit(kIsUsed) )) continue;

                t_Sector_P2 = (Int_t)fparticlecand->getSector();
                t_phi_P2 = fparticlecand->getPhi(); // [deg]
                t_theta_P2 = fparticlecand->getTheta(); // [deg]
                t_mom_P2 = fparticlecand->getMomentum(); // [MeV/c]

                //----------------------------------------------------------------------------------------------
                // Two Partices Variables
                t_phi_diff = TMath::Abs(t_phi_P1 - t_phi_P2);  // [deg]
                t_tan_theta_product = TMath::Tan(TMath::DegToRad() * t_theta_P1)*TMath::Tan(TMath::DegToRad() * t_theta_P2); // [ ]

                //----------------------------------------------------------------------------------------------
                // Histograms filling
                ++has_HH_cnt;

                if(pFPhysics->h_HH_phi_diff) pFPhysics->h_HH_phi_diff->getP()->Fill(t_phi_diff);
                if(pFPhysics->h_HH_phi_diff_closeup) pFPhysics->h_HH_phi_diff_closeup->getP()->Fill(t_phi_diff);
                if(pFPhysics->h_HH_tan_theta_prod) pFPhysics->h_HH_tan_theta_prod->getP()->Fill(t_tan_theta_product);
                if(pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2) pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2->getP()->Fill(t_theta_P1,t_theta_P2);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta) pFPhysics->h_HH_and_HF_mom_vs_theta->getP(0,t_Sector_P1)->Fill(t_mom_P1,t_theta_P1);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta) pFPhysics->h_HH_and_HF_mom_vs_theta->getP(0,t_Sector_P2)->Fill(t_mom_P2,t_theta_P2);

                if ( t_phi_diff >= cut_HH_phi_diff_min && t_phi_diff <= cut_HH_phi_diff_max )
                {
                    if(pFPhysics->h_HH_tan_theta_prod_cut) pFPhysics->h_HH_tan_theta_prod_cut->getP()->Fill(t_tan_theta_product);
                }
                if(t_tan_theta_product >= cut_HH_tan_theta_min && t_tan_theta_product <= cut_HH_tan_theta_max)
                {
                    if(pFPhysics->h_HH_phi_diff_cut) pFPhysics->h_HH_phi_diff_cut->getP()->Fill(t_phi_diff);
                    if(pFPhysics->h_HH_phi_diff_cut_closeup) pFPhysics->h_HH_phi_diff_cut_closeup->getP()->Fill(t_phi_diff);
                }

                // Background reference selection
                if ( ( (t_phi_diff >= bkg_low_min && t_phi_diff <= bkg_low_max) || (t_phi_diff >= bkg_high_min && t_phi_diff <= bkg_high_max) ) && \
                    (t_tan_theta_product >= cut_HH_tan_theta_min && t_tan_theta_product <= cut_HH_tan_theta_max) )
                {
                    // If it's the first background reference pair found in this event
                    if (has_HH_bkg_ref == false)
                    {
                        has_HH_bkg_ref = true;
//                         bkg_phi_diff = t_phi_diff;
                    }
                }

                // Elastic scattering selection
                if ( (t_phi_diff >= cut_HH_phi_diff_min && t_phi_diff <= cut_HH_phi_diff_max) && \
                    (t_tan_theta_product >= cut_HH_tan_theta_min && t_tan_theta_product <= cut_HH_tan_theta_max) )
                {
                    // If it's the first elastic candidate found in this event
                    if (has_HH_elastic == false)
                    {
                        has_HH_elastic = true;

                        //----------------------------------------------------------------------------------------------
                        // Candidate particle 2 - HADES particle
                        el_Sector_P1 = t_Sector_P1;
//                         el_phi_P1 = t_phi_P1;
                        el_theta_P1 = t_theta_P1;
                        el_mom_P1 = t_mom_P1;

                        //----------------------------------------------------------------------------------------------
                        // Candidate particle 2 - HADES particle
                        el_Sector_P2 = t_Sector_P2;
//                         el_phi_P2 = t_phi_P2;
                        el_theta_P2 = t_theta_P2;
                        el_mom_P2 = t_mom_P2;

//                         el_phi_diff = t_phi_diff;
//                         el_tan_theta_product = t_tan_theta_product;
                    }
                }
            }

            //----------------------------------------------------------------------------------------------
            // HADES - FwDet case
            //----------------------------------------------------------------------------------------------
                                                
            // Loop over ForwardCand
            for (Int_t k = 0; k < forward_cand_cnt; k++)
            {
                //----------------------------------------------------------------------------------------------
                // Particle 1 - HADES particle
                HParticleCand* fparticlecand = 0;
                fparticlecand = HCategoryManager::getObject(fparticlecand, fParticleCand, j);

                if (!(fparticlecand->getTofdEdx() > 0.0 && fparticlecand->getSystem() == 1)) continue;

                t_Sector_P1 = (Int_t)fparticlecand->getSector();
//                 t_z_start_P1 = fparticlecand->getZ(); // [mm]
                //t_track_length_P1 = fparticlecand->getDistanceToMetaHit(); // [mm]
                t_phi_P1 = fparticlecand->getPhi(); // [deg]
                t_theta_P1 = fparticlecand->getTheta(); // [deg]
                t_mom_P1 = fparticlecand->getMomentum(); // [MeV/c]

                //----------------------------------------------------------------------------------------------
                // Particle 2 - FwDet Particle
                HForwardCand* fwdetcand = HCategoryManager::getObject(fwdetcand, fForwardCand, k);

                //if(fwdetcand->getTofRec() <= 0) t_ToF_Rec = 0;
                //else if (fwdetcand->getTofRec() > 0) t_ToF_Rec = 1;

                t_z_start_P2 = fparticlecand->getZ(); // [mm]

                // Correction for track start position (Z_START)
                fwdetcand->setStartXYZ(0.0,0.0,t_z_start_P2); // We take the track start position from the reconstructed HADES track paired with FwDet track
                fwdetcand->calcPoints(); // Recalculating ToF and Distance to account for correction
                fwdetcand->calc4vectorProperties(proton_mass); // Calculating 4 vector properties assuming particles are protons

                t_phi_P2 = fwdetcand->getPhi(); // [deg]
                t_theta_P2 = fwdetcand->getTheta(); // [deg]


                //----------------------------------------------------------------------------------------------
                // Two Partices Variables
                t_phi_diff = TMath::Abs(t_phi_P1 - t_phi_P2);  // [deg]
                t_tan_theta_product = TMath::Tan(TMath::DegToRad() * t_theta_P1)*TMath::Tan(TMath::DegToRad() * t_theta_P2); // [ ]

                //----------------------------------------------------------------------------------------------
                // Histograms filling
                ++has_HF_cnt;

                if(pFPhysics->h_HF_phi_diff) pFPhysics->h_HF_phi_diff->getP()->Fill(t_phi_diff);
                if(pFPhysics->h_HF_phi_diff_closeup) pFPhysics->h_HF_phi_diff_closeup->getP()->Fill(t_phi_diff);
                if(pFPhysics->h_HF_tan_theta_prod) pFPhysics->h_HF_tan_theta_prod->getP()->Fill(t_tan_theta_product);
                if(pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2) pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2->getP()->Fill(t_theta_P1,t_theta_P2);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta) pFPhysics->h_HH_and_HF_mom_vs_theta->getP(0,t_Sector_P1)->Fill(t_theta_P1,t_mom_P1);

                if(t_phi_diff >= cut_HF_phi_diff_min && t_phi_diff <= cut_HF_phi_diff_max)
                {
                    if(pFPhysics->h_HF_tan_theta_prod_cut) pFPhysics->h_HF_tan_theta_prod_cut->getP()->Fill(t_tan_theta_product);
                }
                if(t_tan_theta_product >= cut_HF_tan_theta_min && t_tan_theta_product <= cut_HF_tan_theta_max)
                {
                    if(pFPhysics->h_HF_phi_diff_cut) pFPhysics->h_HF_phi_diff_cut->getP()->Fill(t_phi_diff);
                    if(pFPhysics->h_HF_phi_diff_cut_closeup) pFPhysics->h_HF_phi_diff_cut_closeup->getP()->Fill(t_phi_diff);
                }

                // Background reference selection
                if ( ( (t_phi_diff >= bkg_low_min && t_phi_diff <= bkg_low_max) || (t_phi_diff >= bkg_high_min && t_phi_diff <= bkg_high_max) ) && \
                    (t_tan_theta_product >= cut_HF_tan_theta_min && t_tan_theta_product <= cut_HF_tan_theta_max) )
                {
                    // If it's the first background reference pair found in this event
                    if (has_HF_bkg_ref == false)
                    {
                        has_HF_bkg_ref = true;
//                         bkg_phi_diff = t_phi_diff;
                    }
                }

                // Elastic scattering selection
                if ( (t_phi_diff >= cut_HF_phi_diff_min && t_phi_diff <= cut_HF_phi_diff_max) && \
                    (t_tan_theta_product >= cut_HF_tan_theta_min && t_tan_theta_product <= cut_HF_tan_theta_max) )
                {
                    // If it's the first elastic candidate found in this event
                    if (has_HF_elastic == false)
                    {
                        has_HF_elastic = true;
                        //----------------------------------------------------------------------------------------------
                        // Candidate particle 2 - HADES particle
                        el_Sector_P1 = t_Sector_P1;
//                         el_phi_P1 = t_phi_P1;
                        el_theta_P1 = t_theta_P1;
                        el_mom_P1 = t_mom_P1;

                        //----------------------------------------------------------------------------------------------
                        // Candidate particle 2 - HADES particle
                        el_Sector_P2 = t_Sector_P2;
//                         el_phi_P2 = t_phi_P2;
                        el_theta_P2 = t_theta_P2;
                        el_mom_P2 = t_mom_P2;

//                         el_phi_diff = t_phi_diff;
//                         el_tan_theta_product = t_tan_theta_product;
                    }
                }
            }
        }

        //----------------------------------------------------------------------------------------------
        // HADES - HADES case - evaluating candidate event
        if (has_HH_cnt)
        {
            // Filling HH background reference histograms
            if(has_HH_bkg_ref == true && has_HF_bkg_ref == false)
            {
                num_HH_bkg_ref++;
            }
            // Filling HH elastic histograms only if there's JUST ONE elastic candidate in this event
            if(has_HH_elastic == true && has_HF_elastic == false)
            {
                num_HH_el_peak++;
                if(pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2_el) pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2_el->getP()->Fill(el_theta_P1,el_theta_P2);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta_el) pFPhysics->h_HH_and_HF_mom_vs_theta_el->getP(0,el_Sector_P1)->Fill(el_theta_P1,el_mom_P1);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta_el) pFPhysics->h_HH_and_HF_mom_vs_theta_el->getP(0,el_Sector_P2)->Fill(el_theta_P2,el_mom_P2);
            }
        }

        //----------------------------------------------------------------------------------------------
        // HADES - FwDet case - evaluating candidate event
        if (has_HF_cnt)
        {
            // Filling HH background reference histograms
            if(has_HF_bkg_ref == true && has_HH_bkg_ref == false)
            {
                num_HF_bkg_ref++;
            }

            // Filling HF elastic histograms only if there's JUST ONE elastic candidate in this event
            if(has_HF_elastic == true && has_HH_elastic == false)
            {
                num_HF_el_peak++;
                if(pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2_el) pFPhysics->h_HH_and_HF_theta_P1_vs_theta_P2_el->getP()->Fill(el_theta_P1,el_theta_P2);
                if(pFPhysics->h_HH_and_HF_mom_vs_theta_el) pFPhysics->h_HH_and_HF_mom_vs_theta_el->getP(0,el_Sector_P1)->Fill(el_theta_P1,el_mom_P1);
            }
        }

    }

    //----------------------------------------------------------------------------------------------
    // Elastic events counts trend histograms
    //----------------------------------------------------------------------------------------------

    // Elastics number found by integrating over the peak area and subtracting from it the integrals over two areas on both sides of the peak, each half of the peak area wide
    // HADES - FwDet case
    if(pFPhysics->h_HF_elastics_trend && evtCt%(pFPhysics->h_HF_elastics_trend->getRefreshRate()) == 0 && evtCt > 0)
    {
        pFPhysics->h_HF_elastics_trend ->fill(num_HF_el_peak - num_HF_bkg_ref);
        num_HF_el_peak = 0;
        num_HF_bkg_ref = 0;
    }

    // HADES - HADES case
    if(pFPhysics->h_HH_elastics_trend && evtCt%(pFPhysics->h_HH_elastics_trend->getRefreshRate()) == 0 && evtCt > 0)
    {
        pFPhysics->h_HH_elastics_trend ->fill(num_HH_el_peak - num_HH_bkg_ref);
        num_HH_el_peak = 0;
        num_HH_bkg_ref = 0;
    }

    //----------------------------------------------------------------------------------------------
    // Reset if needed
    //----------------------------------------------------------------------------------------------
    mapHolder::resetHists(evtCt);

    return kTRUE;
}