//
//  McpProof.cpp
//  Speciale
//
//  Created by Zaki G on 20/06/16.
//  Copyright © 2016 Zaki G. All rights reserved.
//




//
//  EventPlane.cpp
//  Speciale
//
//  Created by Zaki G on 08/04/16.
//  Copyright © 2016 Zaki G. All rights reserved.
//


#ifndef __CINT__
#include <TStyle.h>
#include "TSelector.h"
#include "TNamed.h"
#include <Riostream.h>
#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 <TString.h>

#include <TProfile.h>
#include <TStopwatch.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TF1.h>
#include <TLegend.h>
#include <sys/time.h>
#include <stdio.h>
#include <time.h>

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <utility>
#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"
#include "ProofEventProc.h"
#else
class TParticlePDG;
class TParticle;
class TLegend;
class Particle;
class Vector;
#endif

using namespace std;



///_____________________________________________________________________________
void ProofEventProc::Begin(TTree *tree){
    
    t = tree;

    h = new Header;
    p = new TClonesArray("TParticle");
    p->SetOwner(true);
    t->SetBranchAddress("header",    &(h->run));
    t->SetBranchAddress("particles", &p);
    pdgDb = TDatabasePDG::Instance();
    
   
}

///_____________________________________________________________________________
void ProofEventProc::SlaveBegin(void){
    
    

    
    
    //TString option = GetOption();
    
    //All the pointers
     Intersection_plot_T0A = new TH2F("MCP Intersection T0A", "MCP Intersection T0A", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    fOutput->Add(Intersection_plot_T0A);
    
    Intersection_plot_T0C = new TH2F("MCP Intersection T0C", "MCP Intersection T0C", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    fOutput->Add(Intersection_plot_T0C);
    
    //TH2F* por = new TH2F("por", "por", 5000, -0.20, 0.20, 5000, -0.20, 0.20);
    
    // TH2F* NParticle_vs_paramter = new TH2F("NParticle_vs_paramter", "N Particle vs N Participants", 500, 0, 150,500, 0, 160000);
    
    //TH1F* CP_hist = new TH1F("CP_hist", "CP_hist", 500, 0, 500000);
    
    //TH2F* dN_dPhi_dEta = new TH2F("dN_dPhi_dEta", "dN_dPhi_dEta", 100, 0, 100, 100, 0, 1000000);
    
    //TH2F* v2_Two = new TH2F("v2_Two", "v2_Two", 50, 0, 20, 50, 0, 0.5);
    //TH2F* v2_Four = new TH2F("v2_Four", "v2_Four", 50, 0, 20, 50, 0, 0.5);
    // Histogram for :
    
    //TH2F* dN_dPhi_dEta = new TH2F("dN_dPhi_dEta", "dN_dPhi_dEta", 100, 0, 7, 100, -5, 5);
    
    //TH2F* hist = new TH2F("v2_b", "v2 vs b", 500, 0, 20, 500, 0, 0.5);
    
    //TH1F* v2_distribution = new TH1F("v2", "v_2 v2_distribution", 400, 0.0, 0.8);
    
    //TH2F* Phi_Eta = new TH2F("Phi_Eta", "Phi_Eta", 100, 0, 7, 100, -5, 5);
    //TH2F* v2_Event = new TH2F("v2_Event", "v2_Event", 50, 0, 20, 50, 0, 0.5);
    
    
     leg = new TLegend(.1 ,.7 ,.3 ,.9 ," ");
    fOutput->Add(leg);
     v2_LYZ = new TH2F("v2_LYZ", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    fOutput->Add(v2_LYZ);
    v2_Two = new TH2F("v2_Two", "v_{2} distribution", 100, 0, 150, 100, 0, 0.5);
    fOutput->Add(v2_Two);
     v2_Four = new TH2F("v2_Four", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    fOutput->Add(v2_Four);
    v2_Event = new TH2F("v2_Event", "v_{2} Distribution", 100, 0, 150, 100, 0, 0.5);
    fOutput->Add(v2_Event);
    
}

///_____________________________________________________________________________
Bool_t ProofEventProc::Process(Long64_t entry )
{
    MCP* MCP_detector_T0C;
    MCP* MCP_detector_T0A;
    
    // 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);
    //------------------------------------------------------------------------
    //
    //              Here the THijing file is opened all the events are looped to be analysed.
    // vector<Vector> particles_intersection_mcp;
    //map< double,vector<int> > AvgCHParameter_M;
    
    //map< double,vector<int> > GlauberParameter_M;
    
    // 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;
    
    //-----------------------------------------------------------------------------------------------------
    // Creating particles
    vector<MCP*> MCP_array_T0C;
    vector<MCP*> MCP_array_T0A;
    
    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_array_T0C.push_back(MCP_detector_T0C);
        
        // 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_array_T0A.push_back(MCP_detector_T0A);
        
    }
    //-----------------------------------------------------------------------------------------------------
  
        vector<Vector> particles_intersection_mcp;
    t->GetEntry(entry);
    if (true) {
    //if (h->event==1) {
        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->GetEntriesFast(); j++) { /// Loop over each particle GetEntriesFast
            TParticle* q = static_cast<TParticle*>(p->At(j));
            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) /*&& ((q->Px() <0.0)&& (q->Px() > -0.02) && (q->Py() > -0.15) &&(q->Py() < -0.12))*/ ) {
                Particle in = Make_Particle(q);
                Vector IP = in.GetVertex();   //IP
                for (unsigned int i = 0; i < MCP_array_T0A.size(); i++){ // Search intersection for each detectors
                    
                    MCP_array_T0C[i]->MCP_Intersection(in, IP, MCP_array_T0C[i], Intersection_plot_T0C, Output, pore_positions, particles_intersection_mcp,1);
                    
                    MCP_array_T0A[i]->MCP_Intersection(in, IP, MCP_array_T0A[i], Intersection_plot_T0A, Output, pore_positions, particles_intersection_mcp,1);
                }
            }
        }
        /// 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);
         //NumberParticles_vs_parameter(particles_intersection_mcp,h->nsnp, NParticle_vs_paramter);*/
        particles_intersection_mcp.clear();
        ///--------------------------------------------------------
    }

    p->Delete();
    return true;
}

///_____________________________________________________________________________
void ProofEventProc::Terminate()
{
    /// All functions which analyses the whole event should start from here:
    // For T0A
    TCanvas* c_T0A= new TCanvas("c_T0A","T0A",700,700);
    Intersection_plot_T0A->SetMarkerStyle(20);
    Intersection_plot_T0A->SetMarkerSize(0.01);
    Intersection_plot_T0A->GetXaxis()->SetLabelSize(0.02);
    Intersection_plot_T0A->GetYaxis()->SetLabelSize(0.02);
    Intersection_plot_T0A->SetXTitle("x-axis");
    Intersection_plot_T0A->SetYTitle("y-axis");
    Intersection_plot_T0A->GetXaxis()->SetTitleOffset(1);
    Intersection_plot_T0A->GetYaxis()->SetTitleOffset(1);
    Intersection_plot_T0A->Draw();
    c_T0A->SaveAs("Plots/Intersection_plot_T0A.png");
    
    
    
}