#ifndef __CINT__
#include <TROOT.h>
#include <TPad.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>
#include <TClonesArray.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TF1.h>
#include <TLegend.h>
#include <stdio.h>
#include <time.h>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <string>
#include <sstream>
#include <cstdlib>
#include <ctime>
#include <vector>
#include <algorithm>
#include "Vector.h"
#include "Particle.h"
#include "mcp.h"
#include "Detector.h"
#include "McpFunctions.h"
#include "Read_data.h"
#include "eventFunction.h"
#else
class TParticlePDG;
class TParticle;
class TLegend;
class Particle;
class Vector;
#endif

using namespace std;

int ep()
{
    // Basic decleration of the MCP:
    Vector Quartz_Dimension(0.053,0.053,0.020);
    Vector Mcp_dimension(0.053,0.053,0.030);
    double space_between_the_sub_quartz = 1.0/1000.0;
    double heigth = Quartz_Dimension.GetY();
    double ElectronicThickness = 0.05;
    double Reflective_index = 1.5;
    int number_of_Mcp = 20;
    int number_of_sub_quartz = 16; // This number should be minimum 1!!!
    double distance_between_quartz = 6.0/1000;
    char nameC[] = "T0C";
    char nameA[] = "T0A";
    vector<Vector> pore_positions; // This vector contains all pores position for both T0A and T0C
    int Output = 0;  // 1 for if the output (8*8) or 0 for intersection i all points on MCP
    int hist_ON_OFF = 0;
    
    //------------------------------------------------------------------------
    // Decleration of the MCP for T0C :
    double z_axis_sign_T0C = 1.0;
    Vector pipe_mcp_position_T0C = Vector(0.0,0.0,-0.7);
    double pipe_radius_T0C = 0.06;
    
    // The position of each MCP at T0C:
    vector<Vector> position_array_quartz_T0C;
    Array_Mcp_Quartz_Position(position_array_quartz_T0C, distance_between_quartz, heigth, pipe_mcp_position_T0C, pipe_radius_T0C);
    pore_position(position_array_quartz_T0C, heigth, pore_positions);
    //------------------------------------------------------------------------
    // Decleration of the MCP for T0A :
    double z_axis_sign_T0A = -1.0;
    Vector pipe_mcp_position_T0A = Vector(0.0,0.0,3.70);
    double pipe_radius_T0A = 0.06;
    
    // The position of each Quartz MCP at T0A:
    vector<Vector> position_array_quartz_T0A;
    Array_Mcp_Quartz_Position(position_array_quartz_T0A, distance_between_quartz, heigth, pipe_mcp_position_T0A, pipe_radius_T0A);
    pore_position(position_array_quartz_T0A, heigth, pore_positions);
    //------------------------------------------------------------------------
    //-----------------------------------------------------------------------------------------------------
    //All the pointers
    TH2F* Intersection_plot_T0A = new TH2F("MCP Intersection T0A", "MCP Intersection T0A", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    TH2F* Intersection_plot_T0C = new TH2F("MCP Intersection T0C", "MCP Intersection T0C", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    TH2F* por = new TH2F("por", "por", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    MCP* MCP_detector_T0C;
    MCP* MCP_detector_T0A;

    TLegend* leg = new TLegend(.1 ,.7 ,.3 ,.9 ," ");
    TH2F* v2_LYZ = new TH2F("v2_LYZ", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    TH2F* v2_Two = new TH2F("v2_Two", "v_{2} distribution", 100, 0, 150, 100, 0, 0.5);
    TH2F* v2_Four = new TH2F("v2_Four", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    TH2F* v2_Event = new TH2F("v2_Event", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    
    //-----------------------------------------------------------------------------------------------------
    //
    //              Here the THijing file is opened all the events are looped to be analysed.
    
    // v2_{EP}
    map< double,vector<double> > Parameter_E;
    // LYZ
    map< int,vector<double> > LYZ_Parameter_Q;
    map< int,vector<double> > LYZ_Parameter_M;
    // v2_{2}
    map< double,vector<double> > Two_Parameter_Q;
    map< double,vector<double> > Two_Parameter_M;
     // v2_{4}
    map< double,vector<double> > Four_Parameter_Q;
    map< double,vector<double> > Four_Parameter_M;
    map< double,vector<double> > Four_Parameter_Two_Q;
    
    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
    TFile *f = TFile::Open("Hijing_000244918_AA_05020_30k.root" ,"READ");
    TTree *t = static_cast<TTree*>(f->Get("T"));
    Header*       h = new Header;
    TClonesArray* p = new TClonesArray("TParticle");
    p->SetOwner(true);
    t->SetBranchAddress("header",    &(h->run));
    t->SetBranchAddress("particles", &p);
   
    for (int i = 0; i < t->GetEntries(); i++) { // Loop over each event
        vector<Vector> particles_intersection_mcp;
        t->GetEntry(i);
      
            Printf("This is event # %d with b=%8.8ffm, and %d particles, c=%5.2f perc. ",h->event, h->b, p->GetEntries(), h->c);
           for (int j = 0;   j < p->GetEntries(); j++) { // Loop over each particle GetEntriesFast
                TParticle* q = static_cast<TParticle*>(p->At(j));
                TParticlePDG* pdg = pdgDb->GetParticle(q->GetPdgCode());
                if (!pdg) {
                    Warning("Foo", "PDG (%s) %d unknown", q->GetName(), q->GetPdgCode());
                    continue;}
    //-----------------------------------------------------------------------------------------------------
    //                                      The particle detection starts from here:
                
                if ( (q->Px() !=NAN) && (q->Py() !=NAN) && (q->Pz() !=NAN) && (q->GetMass() !=0) && (pdg->Charge() !=0) ) {
                    Particle in = Make_Particle(q);
                    Vector IP = in.GetVertex();   //IP
                    
///Multithreading from here->     There is a single histogram in this section of the code, but it is inactive right now. So the whole section here do not contains histogram or ROOT functions, but only detector class functions.
                    for (unsigned int i = 0; i < position_array_quartz_T0A.size(); i++){ // Search intersection for each detectors
                        // For T0C
                        MCP_detector_T0C = new MCP(nameC, Mcp_dimension, position_array_quartz_T0C[i], Quartz_Dimension, number_of_sub_quartz,space_between_the_sub_quartz,Reflective_index, z_axis_sign_T0C);
                        MCP_detector_T0C->MCP_Intersection(in, IP, MCP_detector_T0C, Intersection_plot_T0C, Output, pore_positions, particles_intersection_mcp,0);

                        // For T0A
                        MCP_detector_T0A = new MCP(nameA, Mcp_dimension, position_array_quartz_T0A[i], Quartz_Dimension, number_of_sub_quartz,space_between_the_sub_quartz,Reflective_index, z_axis_sign_T0A);
                        MCP_detector_T0A->MCP_Intersection(in, IP, MCP_detector_T0A, Intersection_plot_T0A, Output, pore_positions, particles_intersection_mcp,0);
                        delete MCP_detector_T0C;
                        delete MCP_detector_T0A;}
///Multithreading to here->
                }
            }
            // All functions which analyse every single event should start from here:
            
            Single_Event_EventPlaneV(particles_intersection_mcp, 2, h->c, "c",  Parameter_E);
            
           Two_avg_Cumulant(particles_intersection_mcp,2,h->c, "c",Two_Parameter_Q,Two_Parameter_M);
            
           Four_avg_Cumulant(particles_intersection_mcp, 2, h->c, "c", Four_Parameter_Q, Four_Parameter_M,Four_Parameter_Two_Q);
            
            Q_theta(particles_intersection_mcp,2, h->c, "c",  LYZ_Parameter_Q, LYZ_Parameter_M);
            ///--------------------------------------------------------
        
        p->Delete();
    }


    // All functions which analyses the whole event should start from here:
    TCanvas* c_v2_b = new TCanvas("c_v2_b", "c_v2_b", 900, 900);
    LYZ(LYZ_Parameter_Q, LYZ_Parameter_M, v2_LYZ, "" );
    EventPlaneV2_vs_Parameter( Parameter_E, v2_Event, "same" );
    V_Two_Cumulant(Two_Parameter_Q, Two_Parameter_M, v2_Two, "same" );
    V_Four_Cumulant(Four_Parameter_Q, Four_Parameter_M, Four_Parameter_Two_Q, v2_Four, "same" );
    
    // Draw the Legend
    leg->SetFillColor (0) ;
    leg->AddEntry(v2_Two,"v2_{2}");
    leg->AddEntry(v2_Four,"v2_{4}");
    leg->AddEntry(v2_Event,"v2_{EP}");
    leg->AddEntry(v2_LYZ,"v2_{LYZ}");
    leg->DrawClone("Same");
    c_v2_b->SaveAs("Plots/v2_c.png");
    
   
     return 0;
}