#ifndef _VECTOR_H
#define _VECTOR_H
#include <iosfwd>
#include <map>


// ===================================================================
class Vector
{
private:
    double fX;
    double fY;
    double fZ;
public:
    Vector(double x=0, double y=0, double z=0);
    ~Vector(void);
    double Length(void) const;
    void   SetX(double x) { fX= x; }
    void   SetY(double y) { fY= y; }
    void   SetZ(double z) { fZ= z; }
    double GetX(void) const { return fX; }
    double GetY(void) const { return fY; }
    double GetZ(void) const { return fZ; }
    
    Vector Cross(const Vector& v) const;
    double Dot(const Vector& v) const;
    
    Vector operator+(const Vector& v);
    Vector operator-(const Vector& v);
    Vector operator=(const Vector& v);
    
};



///-------------------------------------------------------------------
// Global functions
Vector operator*(const double  a, const Vector& v);
Vector operator*(const Vector& v, const Vector& u);
std::ostream& operator<<(std::ostream& o, const Vector& v);


#include <cmath>
#include <iostream>

///-------------------------------------------------------------------
Vector::Vector(double x, double y, double z)
{
    // Vector contructor. Sets the three coordinates
    fX = x;
    fY = y;
    fZ = z;
}



///-------------------------------------------------------------------
Vector::~Vector(){
}


///-------------------------------------------------------------------
double
Vector::Length(void) const
{
    // Return Length of Vector
    return sqrt(fX*fX+fY*fY+fZ*fZ);
}

///-------------------------------------------------------------------
Vector
Vector::Cross(const Vector& v) const
{
    // Cross product between this Vector ansd its argument Vector.
    return Vector(fY*v.GetZ()-fZ*v.GetY(),
                  fZ*v.GetX()-fX*v.GetZ(),
                  fX*v.GetY()-fY*v.GetX());
}

///-------------------------------------------------------------------
double
Vector::Dot(const Vector& v) const
{
    // Dot (scalar)  product between this Vector ansd its argument
    // Vector.
    return (fX*v.GetX() + fY*v.GetY() + fZ*v.GetZ());
}

///-------------------------------------------------------------------
Vector
Vector::operator+(const Vector& v)
{
    // Addition of this Vector, and argument Vector.
    return Vector(fX+v.GetX(), fY+v.GetY(), fZ+v.GetZ());
}

///-------------------------------------------------------------------
Vector
Vector::operator-(const Vector& v)
{
    // Subtraction of this Vector, and argument Vector.
    return Vector(fX-v.GetX(), fY-v.GetY(), fZ-v.GetZ());
}

///-------------------------------------------------------------------
Vector
Vector::operator=(const Vector& v)
{
    // Assignment operator. Return this Vector object.
    fX = v.GetX();
    fY = v.GetY();
    fZ = v.GetZ();
    return (*this);
}

///-------------------------------------------------------------------
Vector
operator*(double a, const Vector& v)
{
    // Global scalar muliplication operator.
    return Vector(a*v.GetX(), a*v.GetY(), a*v.GetZ());
}


///-------------------------------------------------------------------
Vector
operator*(const Vector& v, const Vector& u)
{
    // Gloabal Cross Product operator.
    return v.Cross(u);
}

///-------------------------------------------------------------------
std::ostream&
operator<<(std::ostream& o, const Vector& v)
{
    // Global put-to operator
    o << "(" << v.GetX()
    << "," << v.GetY()
    << "," << v.GetZ()
    << ")";
    
    return o;
}

// ===================================================================
#ifndef Tmap_h
#define Tmap_h
#include <TParameter.h>
#include <map>
#include <vector>
#include "TNamed.h"
#include "TCollection.h"


class Tmap : public TNamed {
    
private:
    std::map< int,std::vector<double> *> *fQ;
    std::vector<double>* fQContainer;
    
    
    public :
    Tmap(void);
    virtual ~Tmap(){};
    
    virtual void TAdd(int key,double Qnumber);
    
    virtual bool TExist(int key);
    virtual std::map< int,std::vector<double>* >* GetTmap(void) const{ return fQ; }

    Long64_t Merge(TCollection *li);

    ClassDef(Tmap,2)
};

ClassImp(Tmap)

Tmap::Tmap(void){
    fQ = new map<int, std::vector<double>*>();
    fQContainer = new vector<double>();
}


void Tmap::TAdd(int key,double Qnumber){
    
    fQContainer->push_back(Qnumber);
    fQ->insert(std::pair<int, std::vector<double>*>(key, fQContainer));
    
}


bool Tmap::TExist(int key){
    
    bool exist;
    if ( fQ->find(key) == fQ->end()){
        
        return false;
    }
    
    else{
        return true;}
}

// Merge Tmap objects in collection 'li' into itself
Long64_t Tmap::Merge(TCollection *li)
{

    Long64_t rc = 0;

    // If empty list, nothing to do
    if (!li || li->GetSize() < 1) return rc;

    // Loop over the objects in the collection
    TIter nxtm(li);
    Tmap *m = 0;
    // (Assume objects in the collection are Tmap; a check could be added ...)
    while ( (m = (Tmap *) nxtm())) {
        // link to the container
        std::map< int,std::vector<double>* >* tm = m->GetTmap();
        std::map< int,std::vector<double>* >::iterator tm_iter;
        unsigned int i = 0;
        for(tm_iter = tm->begin(); tm_iter != tm->end(); tm_iter++) {
            int key = tm_iter->first;
            std::vector<double> cont = *(tm_iter->second);
            TAdd(key, cont[i]);
            i++;
        }
        rc += i;
        std::cout<<"Merged "<<m->GetName()<<" , size: "<< m->GetTmap()->size() << "\n" ;
    }
    return rc;    
}

#endif


// ===================================================================
#include <string.h>

class Particle
{
private:
    // Data members
    string fName;     // particle name (C++ string)
    double fMass;     // Mass of the particle
    int    fCharge;   // Charge of particle
    Vector fMomentum; // Momentum vector
    Vector fVertex;
public:
    // Method members
    Particle(string name, double mass, int charge);
    ~Particle(void) {}
    void   SetMomentum(Vector& v)  { fMomentum = v; }
    void   SetVertex(Vector& v)  { fVertex = v; }
    
    string GetName(void)     const { return fName; }
    double GetMass(void)     const { return fMass; }
    int    GetCharge(void)   const { return fCharge; }
    
    double GetEnergy(void)   const;
    double GetGamma(void)    const;
    
    Vector GetVelocity(void) const;
    Vector GetVertex(void) const;
    Vector GetMomentum(void) const { return fMomentum; }
};
#include <cmath>
#include <string>

///-------------------------------------------------------------------
Particle::Particle(string name, double mass, int charge)
{
    // Particle constructor. Sets name, mass, and charge. Momentum and
    // Position are initially set to 0 vectors
    
    // No massless particles!!!!
    if (mass <= 0) return;
    
    fName = name;
    fMass = mass;
    fCharge = charge;
    
    fMomentum = Vector(0,0,0);
}

///-------------------------------------------------------------------
double
Particle::GetEnergy(void) const
{
    // Return the energy of th Particle.
    double p = fMomentum.Length();
    return sqrt(p * p + fMass * fMass);
}

///-------------------------------------------------------------------
double
Particle::GetGamma(void) const
{
    // Return the gamma factor of the Particle
    
    double p  = fMomentum.Length();
    double EE = p * p + fMass * fMass;
    return (1 / sqrt(1 - (p * p / EE) ));
}


///-------------------------------------------------------------------
Vector
Particle::GetVertex(void) const
{
    return fVertex;
}


///-------------------------------------------------------------------
Vector
Particle::GetVelocity(void) const
{
    return (1.0/(GetMass()*GetGamma()))*GetMomentum();
}

// ===================================================================
struct Header {
    unsigned int    run;
    unsigned int    event;
    unsigned int    ntgt;
    unsigned int    nproj;
    unsigned int    nbin;
    unsigned int    type;
    double  ipx;
    double  ipy;
    double  ipz;
    double  b;
    double  c;
    double  phir;
    unsigned int    nsnp;
    unsigned int    nsnt;
    unsigned int    nspp;
    unsigned int    nspt;
};

// ===================================================================
void quart_div(int number_of_sub_quartz,
               vector<Vector> &position_array_quartz,
               double &quartz_hight,
               Vector pipe_mcp_position){
    
    //This function divides quartz in 4 equal pieces (square):
    
    
    //Now we consider the case where the there is sub quartz and define there positions:
    double space_between_the_sub_quartz = 1.0/1000.0;
    vector<Vector> temp(position_array_quartz);
    position_array_quartz.clear();
    if (number_of_sub_quartz == 4) {
        for (unsigned int i=0; i<temp.size(); i++) {
            // Defining the start position fo the first sub quartz:
            
            //              ------------
            //             |     |      |
            //             |  a  |  b   |
            //             |------------|
            //             |     |      |               Front view of the MCP
            //             |  c  |  d   |
            //              ------------
            
            double start_q_x = (temp[i]).GetX() - ((quartz_hight - space_between_the_sub_quartz)/  number_of_sub_quartz);
            double start_q_y =(temp[i]).GetY() + ((quartz_hight - space_between_the_sub_quartz)/ number_of_sub_quartz);
            
            position_array_quartz.push_back(Vector(start_q_x, start_q_y, pipe_mcp_position.GetZ())); //a
            
            double space = (quartz_hight - space_between_the_sub_quartz)/ number_of_sub_quartz;
            
            position_array_quartz.push_back(Vector(start_q_x + space + space + space_between_the_sub_quartz, start_q_y, pipe_mcp_position.GetZ()));  //b
            
            
            position_array_quartz.push_back(Vector(start_q_x, start_q_y - space - space_between_the_sub_quartz -space, pipe_mcp_position.GetZ())); //c
            
            
            position_array_quartz.push_back(Vector(start_q_x + space + space_between_the_sub_quartz + space,     start_q_y - space - space_between_the_sub_quartz -space, pipe_mcp_position.GetZ())); // d
            
            
        }
        
    }
    
    else if (number_of_sub_quartz == 16){
        
        for (unsigned int i=0; i<temp.size(); i++) {
            
            Vector upper_left_corner = Vector(    ((temp[i]).GetX() - (quartz_hight/2)),    ((temp[i]).GetY() + (quartz_hight/2)),   pipe_mcp_position.GetZ());
            
            double xa = (quartz_hight - 3*space_between_the_sub_quartz)/  (sqrt(number_of_sub_quartz)*2) ;
            
            double ya = (quartz_hight - 3*space_between_the_sub_quartz)/  (sqrt(number_of_sub_quartz)*2);
            
            
            Vector start_point( upper_left_corner.GetX() + xa,        upper_left_corner.GetY()-ya,       upper_left_corner.GetZ());
            
            
            double spacy = 2*xa + space_between_the_sub_quartz;
            
            for (int y=0; y<4; y++) {
                for (int x=0; x<4; x++){
                    Vector position_4_4 = Vector(  start_point.GetX() + (x*(spacy)),       start_point.GetY() - (y*(spacy)),      start_point.GetZ());
                    
                    position_array_quartz.push_back(position_4_4);
                }
            }
            
        }
        
        
    }
}



///-------------------------------------------------------------------
void Array_Mcp_Quartz_Position(vector<Vector> &position_array_quartz,
                               double &distance_between_quartz,
                               double quartz_hight,
                               Vector &pipe_mcp_position,
                               double pipe_diameter){
    
    //Here I define the first 8 detectors:
    
    //                      8
    //                      9
    //              1   5       14  18
    //                      10
    //                      11
    //
    //
    
    //       ---------------Christian kig her ......(De første 8 detector bliver defineret her)-----
    
    const double distance_first_mcp = ((pipe_diameter/2) + (0.5*distance_between_quartz + quartz_hight*0.5));
    const double distance_second_mcp = distance_first_mcp + (distance_between_quartz+quartz_hight);
    
    for (int i = 0; i<4; i++) {
        int pos_x = int(cos(double((M_PI/2.0)*i)));
        int pos_y = int(sin(double((M_PI/2.0)*i)));
        
        position_array_quartz.push_back(Vector(distance_first_mcp*pos_x, distance_first_mcp*pos_y,  pipe_mcp_position.GetZ()));
        position_array_quartz.push_back(Vector(distance_second_mcp*pos_x, distance_second_mcp*pos_y,  pipe_mcp_position.GetZ()));
    }
    
    
    //The MCP are divided in 4 parts with three detector each:
    
    //          part2       8    3
    //                     [2]   1   2     <-- part1
    //              1   5       [0]  [1]
    //                      10
    //          part3       11    part4
    //
    //
    
    for (int i = 0; i<1; i++) {
        //part1
        Vector Vector1 = Vector( position_array_quartz[2].GetX()+(distance_between_quartz+quartz_hight),    position_array_quartz[0].GetY()+(distance_between_quartz+quartz_hight),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector1);
        
        Vector Vector2 = Vector( Vector1.GetX()+(distance_between_quartz+quartz_hight),  Vector1.GetY(),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector2);
        
        Vector Vector3 = Vector( Vector1.GetX(),  Vector1.GetY()+(distance_between_quartz+quartz_hight),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector3);
        
        //part2
        Vector Vector4 = Vector( -(position_array_quartz[2].GetX()+(distance_between_quartz+quartz_hight)),    position_array_quartz[0].GetY()+(distance_between_quartz+quartz_hight),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector4);
        
        Vector Vector5 = Vector(- (Vector1.GetX()+(distance_between_quartz+quartz_hight)),  Vector1.GetY(),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector5);
        
        Vector Vector6 = Vector(- (Vector1.GetX()),  Vector1.GetY()+(distance_between_quartz+quartz_hight),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector6);
        
        //part3
        Vector Vector7 = Vector( -(position_array_quartz[2].GetX()+(distance_between_quartz+quartz_hight)),   -(position_array_quartz[0].GetY()+(distance_between_quartz+quartz_hight)),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector7);
        
        Vector Vector8 = Vector(-( Vector1.GetX()+(distance_between_quartz+quartz_hight)),  -Vector1.GetY(),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector8);
        
        Vector Vector9 = Vector(-( Vector1.GetX()), -( Vector1.GetY()+(distance_between_quartz+quartz_hight)),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector9);
        
        //part4
        Vector Vector10 = Vector( +(position_array_quartz[2].GetX()+(distance_between_quartz+quartz_hight)),   -(position_array_quartz[0].GetY()+(distance_between_quartz+quartz_hight)),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector10);
        
        Vector Vector11 = Vector(+( Vector1.GetX()+(distance_between_quartz+quartz_hight)),  -Vector1.GetY(),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector11);
        
        Vector Vector12 = Vector(+ Vector1.GetX(), -( Vector1.GetY()+(distance_between_quartz+quartz_hight)),   pipe_mcp_position.GetZ());
        
        position_array_quartz.push_back(Vector12);
        
    }
    
    
    
    
}

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

void pore_position(vector<Vector> &position_array_quartz, const double &heigth, vector<Vector> &pore_positions){
    
    //This function initialize all the pore positions for all detetctors.
    
    
    // upper_left_corner  ---->  ------------
    //                          | . . . . .  |
    //                          | . . . . .  |   <-- 8 by 8 pores
    //                          | . . . . .  |
    //                          | . . . . .  |
    //                          | . . . . .  |
    //                           ------------
    
    
    for (unsigned int pos = 0 ; pos<position_array_quartz.size(); pos++) {
        
        Vector position = position_array_quartz[pos];
        
        Vector upper_left_corner = Vector(    (position.GetX()-heigth/2.0),       (position.GetY()+heigth/2.0),      position.GetZ());
        
        
        Vector start_point(upper_left_corner.GetX()+0.0064,        upper_left_corner.GetY()-0.0064,       upper_left_corner.GetZ());
        
        
        for (int pore_y=0; pore_y<8; pore_y++) {
            for (int pore_x=0; pore_x<8; pore_x++){
                Vector pore = Vector(  start_point.GetX() + (pore_x*(0.0065)),       start_point.GetY() - (pore_y*(0.0065)),      start_point.GetZ());
                
                pore_positions.push_back(pore);
            }
        }
        
    }
}
// ===================================================================
#include <cmath>
#include <iostream>
#include <iomanip>
#include <vector>
#include "TROOT.h"
#include "TRandom3.h"
#include "TH2F.h"
#include "TCanvas.h"
#include <TMath.h>
#include "TProofServ.h"

class MCP
{
private:
    Vector MCP_Dimension;
    Vector MCP_QuartDimension;
    Vector MCP_Position;
    int MCP_NumberSubQuartz;
    double MCP_Space_betw_quartz;
    double MCP_Quartz_Reflective_index;
    double MCP_z_axis_sign;
    double MCP_end_pos_sign;
    vector<double> MCP_Incident_Angle;
    
public:
    Vector GetMcpPos(void) const { return MCP_Position; }
    
    MCP(char* name,Vector& dimension, const Vector &position, Vector& QuartzDimension, int NumberSubQuartz, double space_between_the_sub_quartz, double Reflective_index, double z_axis_sign);
    ~MCP(void);
    
    Vector Recurssiv_wall_Intersection(Vector &path_photon, Vector &Cherenkov_photon, MCP &name, Vector &Last_wall_intersection);
    void Append_angles(double a);
    void Get_angles(vector<double> &Angels_array);
    
    void MCP_quartz_div(vector<Vector> &Single_position_quartz_array);
    
    
    void pore_intersection(int pore_ON_OFF,   vector<Vector> &pore_positions,  Vector &intersect_point,   TH2F *Intersection_plot, vector<Vector> &event_plan_array, int &hist_ON_OFF);
    
    Vector normal_vector_initilizer(Vector &point, Vector &back_point, double thickness_of_mcp);
    
    Vector Quartz_wall_Intersection(MCP* name, Vector &Last_wall_intersection, Vector &path_photon, Vector &direction_vector, int reflection_on_off, Vector &Cherenkov_photon);
    
    Vector MCP_end_Intersection_ref(Vector &Photon_momentum, Vector &particle_pos);
    Vector MCP_end_Intersection_for_particle( Particle& pname, Vector &particle_pos, MCP* name);
    
    Vector MCP_end_Intersection_for_Photon(Vector &Photon_momentum, Vector &particle_pos, MCP* name);
    
    void Cherenkov_Photon_Generation( Particle& pname, Vector &first_intersection, Vector &second_intersection, vector<Vector> &pathPhoton_array, vector<Vector> &CherenkovPhoton_array);
    
    void MCP_Intersection(Particle& pname, Vector &particle_pos, MCP* name , TH2F *Intersection_plot, int pore_ON_OFF, vector<Vector> &pore_positions, vector<Vector> &event_plan_array, int hist_ON_OFF);
    
    
    
    
    
};


//____________________________________________________________
// Functions for critical angle:
void MCP::Append_angles(double a){
    MCP_Incident_Angle.push_back(a);
}

void MCP::Get_angles(vector<double> &Angels_array){
    
    for (unsigned int i=0; i<MCP_Incident_Angle.size(); i++){
        Angels_array.push_back(MCP_Incident_Angle[i]);
    }
    
}
//____________________________________________________________

MCP::MCP(char* name,Vector& dimension,const Vector &position, Vector& QuartzDimension, int NumberSubQuartz, double space_between_the_sub_quartz, double Reflective_index, double z_axis_sign)//: Detector(name, dimension, 0)
{
    MCP_Dimension = dimension;
    MCP_Position = position;
    MCP_QuartDimension = QuartzDimension;
    MCP_Space_betw_quartz = space_between_the_sub_quartz;
    MCP_NumberSubQuartz = NumberSubQuartz;
    MCP_Quartz_Reflective_index = Reflective_index;
    MCP_z_axis_sign = z_axis_sign;
    
    int end_sign =1.0;
    if (MCP_z_axis_sign>0) {
        end_sign =-1.0;
    }
    else {
        end_sign =+1.0;
    }
    MCP_end_pos_sign = end_sign;
}

MCP::~MCP(){
    
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

void MCP::MCP_quartz_div(vector<Vector> &Single_position_quartz_array){
    // Single_position_quartz_array contains all the quartz positions for single MCP
    
    //This function divides quartz in 4 equal pieces (square):
    
    double quartz_hight = MCP_Dimension.GetY();
    double pipe_mcp_position_z_axis = MCP_Position.GetZ();
    
    //Now we consider the case where the there is sub quartz and define there positions:
    double space_between_the_sub_quartz = MCP_Space_betw_quartz;//1.0/1000.0;
    
    if (MCP_NumberSubQuartz == 1) {
        Single_position_quartz_array.push_back(MCP_Position);
    }
    
    else if (MCP_NumberSubQuartz == 4) {
        // Defining the start position fo the first sub quartz:
        //              ------------
        //             |     |      |
        //             |  a  |  b   |
        //             |------------|
        //             |     |      |               Front view of the MCP
        //             |  c  |  d   |
        //              ------------
        
        double start_q_x = MCP_Position.GetX() - (space_between_the_sub_quartz/2.0)- ((quartz_hight - space_between_the_sub_quartz)/  MCP_NumberSubQuartz);
        double start_q_y =MCP_Position.GetY()+ (space_between_the_sub_quartz/2.0) + ((quartz_hight - space_between_the_sub_quartz)/ MCP_NumberSubQuartz);
        
        Single_position_quartz_array.push_back(Vector(start_q_x, start_q_y, pipe_mcp_position_z_axis)); //a
        
        double space = (quartz_hight - space_between_the_sub_quartz)/ MCP_NumberSubQuartz;
        
        Single_position_quartz_array.push_back(Vector(start_q_x + space + space + space_between_the_sub_quartz, start_q_y, pipe_mcp_position_z_axis));  //b
        
        
        Single_position_quartz_array.push_back(Vector(start_q_x, start_q_y - space - space_between_the_sub_quartz -space, pipe_mcp_position_z_axis)); //c
        
        
        Single_position_quartz_array.push_back(Vector(start_q_x + space + space_between_the_sub_quartz + space,     start_q_y - space - space_between_the_sub_quartz -space, pipe_mcp_position_z_axis)); // d
        
    }
    
    else if (MCP_NumberSubQuartz == 16){
        
        Vector upper_left_corner = Vector(    (MCP_Position.GetX() - (quartz_hight/2.0)),    (MCP_Position.GetY() + (quartz_hight/2.0)),   pipe_mcp_position_z_axis);
        
        double xa = (quartz_hight - 3.0*space_between_the_sub_quartz)/  (sqrt(MCP_NumberSubQuartz)*2.0) ;
        
        double ya = (quartz_hight - 3.0*space_between_the_sub_quartz)/  (sqrt(MCP_NumberSubQuartz)*2.0);
        
        
        Vector start_point( upper_left_corner.GetX() + xa,        upper_left_corner.GetY()-ya,       upper_left_corner.GetZ());
        
        
        double spacy = 2.0*xa + space_between_the_sub_quartz;
        
        for (int y=0; y<4; y++) {
            for (int x=0; x<4; x++){
                Vector position_4_4 = Vector(  start_point.GetX() + (x*(spacy)),       start_point.GetY() - (y*(spacy)),      start_point.GetZ());
                
                Single_position_quartz_array.push_back(position_4_4);
            }
        }
    }
}


//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

void MCP::pore_intersection(int pore_ON_OFF,   vector<Vector> &pore_positions,  Vector &intersect_point,   TH2F *Intersection_plot, vector<Vector> &event_plan_array, int &hist_ON_OFF){
    //This functions compare all the actual intersection om the MCP surface, and choose the point
    // which is close to one of the pores and plot the results.
    
    if (    (pore_ON_OFF == 1) &&  (intersect_point.GetX()!= 0) && (intersect_point.GetY()!= 0) && (TMath::IsNaN(intersect_point.GetX())!= 1) && (TMath::IsNaN(intersect_point.GetY())!= 1) ) {
        for (unsigned int i=0; i < pore_positions.size(); i++){
            double x_distance = abs( pore_positions[i].GetX() -  intersect_point.GetX());
            double y_distance = abs( pore_positions[i].GetY() -  intersect_point.GetY());
            
            
            if (    (x_distance <= 0.0065/2.0)  &&  (y_distance <= 0.0065/2.0)    ) {
                if (hist_ON_OFF == 1) {
                    Intersection_plot->Fill(pore_positions[i].GetX(),  pore_positions[i].GetY());
                }
                event_plan_array.push_back(pore_positions[i]);
            }
        }
    }
    
    else if ((pore_ON_OFF == 0) &&  (intersect_point.GetX()!= 0) && (intersect_point.GetY()!= 0) && (TMath::IsNaN(intersect_point.GetX())!= 1) && (TMath::IsNaN(intersect_point.GetY())!= 1)  ){
        //        cout << "intersect_point: "<<  intersect_point.GetX()<<endl;
        
        
        if (hist_ON_OFF == 1) {
            Intersection_plot->Fill(intersect_point.GetX(),  intersect_point.GetY());
        }
        
        event_plan_array.push_back(intersect_point);
        
    }
    
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

Vector MCP::normal_vector_initilizer(Vector &point, Vector &back_point, double thickness_of_mcp){
    
    Vector normal_vector = Vector(point-back_point).Cross(Vector(0,0, (back_point.GetZ()+(+1*MCP_end_pos_sign)*thickness_of_mcp) - back_point.GetZ() ));
    return normal_vector;
}

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
Vector MCP::Quartz_wall_Intersection(MCP* name, Vector &Last_wall_intersection, Vector &path_photon, Vector &direction_vector, int reflection_on_off, Vector &Cherenkov_photon) //The direction vector is the vector against the normal vectors of the walls, and the Cherenkov_photon is only called to
{
    vector<Vector> Single_position_quartz_array; // The container of all the sub quartz positions.
    MCP_quartz_div(Single_position_quartz_array);
    Vector MCP_QuartDim;
    if(MCP_NumberSubQuartz > 1){
        MCP_QuartDim = Vector(   (MCP_QuartDimension.GetX() - (sqrt(MCP_NumberSubQuartz)-1)*MCP_Space_betw_quartz)/  (sqrt(MCP_NumberSubQuartz)),       (MCP_QuartDimension.GetY() - (sqrt(MCP_NumberSubQuartz)-1)*MCP_Space_betw_quartz)/  (sqrt(MCP_NumberSubQuartz)),         MCP_QuartDimension.GetZ());
    }
    else if(MCP_NumberSubQuartz==1){
        MCP_QuartDim = Vector(   MCP_QuartDimension.GetX(),     MCP_QuartDimension.GetY(),         MCP_QuartDimension.GetZ());
    }
    
    
    Vector Sub_Quartz_position;
    for (unsigned int i=0; i<Single_position_quartz_array.size(); i++) { // Loop over every sub quartz of the single MCP
        
        //cout << "0.5"<< endl;
        //cout << "path_photon"<< path_photon<<endl;
        
        Vector w_path = Single_position_quartz_array[i] - path_photon;
        
        if ((abs(w_path.GetX()) < (MCP_QuartDim.GetX()/2.0)) && (abs(w_path.GetY()) < (MCP_QuartDim.GetX()/2.0)))
        {Sub_Quartz_position = Single_position_quartz_array[i];
        }
    }
    
    if (Sub_Quartz_position.GetX()==0 && Sub_Quartz_position.GetY()==0)
    {
        return Vector(NAN,NAN,NAN);
    }
    
    Vector return_intersection_vector;
    
    Vector Quartz_wall_Intersection; //The intersection point on one of the walls
    
    // Function for finding the intersection with the quartz walls and it can also calculated the reflected photon vector:
    //              ------------
    //             |            |
    //             |            |
    //             |            |
    //             |            |               Front view of the MCP
    //             |            |
    //              ------------
    Vector detector_posi = Sub_Quartz_position; // The position of the detector
    
    // Front view of the MCP:
    Vector Quartz_LDcorner = Vector(double(detector_posi.GetX()-(MCP_QuartDim.GetX()/2)),   detector_posi.GetY() - MCP_QuartDim.GetY()/2,  detector_posi.GetZ() ); //The position of the Quartz down left corner.
    
    Vector Quartz_LUcorner = Vector(Quartz_LDcorner.GetX(),   Quartz_LDcorner.GetY() + MCP_QuartDim.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up left corner.
    
    Vector Quartz_RUcorner = Vector(Quartz_LDcorner.GetX()+(MCP_QuartDim.GetX()),   Quartz_LDcorner.GetY() + MCP_QuartDim.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up right corner.
    
    Vector Quartz_RDcorner = Vector(Quartz_LDcorner.GetX()+(MCP_QuartDim.GetX()),   Quartz_LDcorner.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up right corner.
    
    Vector QuartzSurfacePoints_array[4] = {Quartz_LDcorner, Quartz_RDcorner, Quartz_RUcorner, Quartz_LDcorner};
    
    //Defining the four normal vectors for the quartz side walls and roofs:
    //              ------------
    //             |      |     |
    //             |      v     |
    //             |-->      <--|
    //             |      ^     |               Front view of the MCP
    //             |      |     |
    //              ------------
    
    //Walls normal vector:
    //Left :
    Vector normal_left = (MCP_end_pos_sign) *normal_vector_initilizer(Quartz_LUcorner, Quartz_LDcorner, MCP_QuartDim.GetZ());
    
    //Right:
    Vector normal_right =  (-1.0*MCP_end_pos_sign) *normal_vector_initilizer(Quartz_RUcorner, Quartz_RDcorner, MCP_QuartDim.GetZ());
    
    //Roof normal vectors:
    //Top:
    Vector normal_top = (MCP_end_pos_sign)*normal_vector_initilizer(Quartz_RUcorner, Quartz_LUcorner, MCP_QuartDim.GetZ());
    
    //Bottom:
    Vector normal_bottom = (-1.0*MCP_end_pos_sign) *normal_vector_initilizer(Quartz_RDcorner, Quartz_LDcorner, MCP_QuartDim.GetZ());
    
    // Putting the normal vectors in a array
    Vector normal_walls_roof[4]={normal_left,normal_right,normal_top,normal_bottom};
    
    // Finding the center position of the walls and roof:
    //                   * <--center position of the roof
    //              ------------
    //             |            |
    //             |            |
    //            *|            |*
    //             |            |               Front view of the MCP
    //             |            |
    //              ------------
    //                   *
    
    // Husk for tegne på hvor meget man skal den går tilbage i z aksen
    Vector Left = Quartz_LDcorner + Vector(0,   MCP_QuartDim.GetY()/2.0 , (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ())/2.0 ) ; //Left //sign
    Vector Right = Quartz_RDcorner+ Vector(0,   MCP_QuartDim.GetY()/2.0 ,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) ); // right
    Vector Top = Quartz_LUcorner  + Vector(MCP_QuartDim.GetX()/2.0 ,   0,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) ); //Top
    Vector Bottom = Quartz_LDcorner   + Vector(MCP_QuartDim.GetX()/2.0 ,   0,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) );//Bottom
    
    Vector quartz_surface_position[4]={Left,Right,Top,Bottom}; // position of the walls and roof
    
    //Looping over which wall it hits
    //cout << "3"<< endl;
    for (int i=0; i<=3; i++){
        
        double dot_normal = direction_vector.Dot(normal_walls_roof[i]);
        
        if (dot_normal < 0.0){ //the dot product should be less than zero so the angle is stump (over 90 degree)
            
            //The intersection point on one of the walls
            Vector Quartz_wall_Intersection = path_photon-((( path_photon - QuartzSurfacePoints_array[i] ).Dot(normal_walls_roof[i])/dot_normal)*direction_vector);
            
            Vector w_quartz = Quartz_wall_Intersection - quartz_surface_position[i];
            
            if ((abs(w_quartz.GetY()) < (MCP_QuartDim.GetY())/2.0) && (abs(w_quartz.GetZ()) < (MCP_QuartDim.GetZ())/2.0)){
                
                if (reflection_on_off==1){ //If the reflected momentum is needed:
                    // Refelction vector with in the quartz walls for the Cherenkov photons
                    //Vector Reflection_vector_quartz = Cherenkov_photon +(-2.0*(Cherenkov_photon.Dot(normal_walls_roof[i]))*normal_walls_roof[i]);
                    
                    Vector n  = normal_walls_roof[i];
                    
                    Vector eC =(1.0/Cherenkov_photon.Length())*Cherenkov_photon;
                    
                    Vector en =(1.0/n.Length())*n;
                    
                    Vector Reflection_vector_quartz = eC +(-2.0*(eC.Dot(en))*en);
                    char myword[] = "a";
                    
                    MCP MCP_reflection(myword, MCP_QuartDim, detector_posi, MCP_QuartDim, 1,0,0, MCP_z_axis_sign);
                    Vector RefInter = MCP_reflection.MCP_end_Intersection_ref(Reflection_vector_quartz, Quartz_wall_Intersection);
                    
                    if ((  (RefInter.GetX()!= 0) && (RefInter.GetY()!= 0)  )) {
                        Last_wall_intersection = Quartz_wall_Intersection;
                        return_intersection_vector = Reflection_vector_quartz;
                        
                        return Reflection_vector_quartz;
                        
                    }
                    
                    else{
                        
                        
                        return_intersection_vector = Recurssiv_wall_Intersection(Quartz_wall_Intersection, Reflection_vector_quartz, MCP_reflection, Last_wall_intersection);
                        
                        return return_intersection_vector;
                    }
                    
                }
                
                else{return_intersection_vector = Quartz_wall_Intersection;
                    return return_intersection_vector;
                }
            }
        }
    }
    
    return return_intersection_vector;
}

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

Vector MCP::Recurssiv_wall_Intersection(Vector &path_photon, Vector &Cherenkov_photon, MCP &name, Vector &Last_wall_intersection){
    Vector MCP_QuartDim = name.MCP_QuartDimension;
    Vector return_intersection_vector;
    // Function for finding the intersection with the quartz walls and it can also calculated the reflected photon vector:
    //              ------------
    //             |            |
    //             |            |
    //             |            |
    //             |            |               Front view of the MCP
    //             |            |
    //              ------------
    
    //Vector detector_posi = Single_position_quartz_array; // The position of the detector
    
    Vector detector_posi = name.MCP_Position; // The position of the detector
    
    // Front view of the MCP:
    Vector Quartz_LDcorner = Vector(double(detector_posi.GetX()-(MCP_QuartDim.GetX()/2.0)),   detector_posi.GetY() - MCP_QuartDim.GetY()/2.0,  detector_posi.GetZ() ); //The position of the Quartz down left corner.
    
    Vector Quartz_LUcorner = Vector(Quartz_LDcorner.GetX(),   Quartz_LDcorner.GetY() + MCP_QuartDim.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up left corner.
    
    Vector Quartz_RUcorner = Vector(Quartz_LDcorner.GetX()+(MCP_QuartDim.GetX()),   Quartz_LDcorner.GetY() + MCP_QuartDim.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up right corner.
    
    Vector Quartz_RDcorner = Vector(Quartz_LDcorner.GetX()+(MCP_QuartDim.GetX()),   Quartz_LDcorner.GetY(),  Quartz_LDcorner.GetZ() ); //The position of the Quartz up right corner.
    
    Vector QuartzSurfacePoints_array[4] = {Quartz_LDcorner, Quartz_RDcorner, Quartz_RUcorner, Quartz_LDcorner};
    
    
    //Defining the four normal vectors for the quartz side walls and roofs:
    //              ------------
    //             |      |     |
    //             |      v     |
    //             |-->      <--|
    //             |      ^     |               Front view of the MCP
    //             |      |     |
    //              ------------
    //Walls normal vector:
    //Left :
    Vector normal_left = (MCP_end_pos_sign)*normal_vector_initilizer(Quartz_LUcorner, Quartz_LDcorner, MCP_QuartDim.GetZ());
    
    //Right:
    Vector normal_right = (-1.0*MCP_end_pos_sign) * normal_vector_initilizer(Quartz_RUcorner, Quartz_RDcorner, MCP_QuartDim.GetZ());
    
    //Roof normal vectors:
    //Top:
    Vector normal_top = (MCP_end_pos_sign)*normal_vector_initilizer(Quartz_RUcorner, Quartz_LUcorner, MCP_QuartDim.GetZ());
    
    //Bottom:
    Vector normal_bottom = (-1.0*MCP_end_pos_sign) *normal_vector_initilizer(Quartz_RDcorner, Quartz_LDcorner, MCP_QuartDim.GetZ());
    
    // Putting the normal vectors in a array
    Vector normal_walls_roof[4]={normal_left,normal_right,normal_top,normal_bottom};
    
    // Finding the center position of the walls and roof:
    //                   * <--center position of the roof
    //              ------------
    //             |            |
    //             |            |
    //            *|            |*
    //             |            |               Front view of the MCP
    //             |            |
    //              ------------
    //                   *
    
    // Husk for tegne på hvor meget man skal den går tilbage i z aksen
    Vector Left = Quartz_LDcorner + Vector(0,   MCP_QuartDim.GetY()/2.0 , (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ())/2.0 ) ; //Left //sign
    Vector Right = Quartz_RDcorner+ Vector(0,   MCP_QuartDim.GetY()/2.0 ,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) ); // right
    Vector Top = Quartz_LUcorner  + Vector(MCP_QuartDim.GetX()/2.0 ,   0,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) ); //Top
    Vector Bottom = Quartz_LDcorner   + Vector(MCP_QuartDim.GetX()/2.0 ,   0,  (+1.0*MCP_end_pos_sign) *(MCP_QuartDim.GetZ()/2.0) );//Bottom
    
    Vector quartz_surface_position[4]={Left,Right,Top,Bottom}; // position of the walls and roof
    
    //Looping over which wall it hits
    for (int i=0; i<=3; i++){
        double dot_normal = Cherenkov_photon.Dot(normal_walls_roof[i]);
        if (dot_normal < 0.0){ //the dot product should be less than zero so the angle is stump (over 90 degree)
            //The intersection point on one of the walls
            Vector Quartz_wall_Intersection = path_photon-((( path_photon - QuartzSurfacePoints_array[i] ).Dot(normal_walls_roof[i])/dot_normal)*Cherenkov_photon);
            
            Vector w_quartz = Quartz_wall_Intersection - quartz_surface_position[i];
            
            if ((abs(w_quartz.GetY()) < (MCP_QuartDim.GetY())/2.0) && (abs(w_quartz.GetZ()) < (MCP_QuartDim.GetZ())/2.0)){
                
                //Vector Reflection_vector_quartz = Cherenkov_photon +(-2.0*(Cherenkov_photon.Dot(normal_walls_roof[i]))*normal_walls_roof[i]);
                Vector n  = normal_walls_roof[i];
                Vector eC =(1/Cherenkov_photon.Length())*Cherenkov_photon;
                Vector en =(1/n.Length())*n;
                
                Vector Reflection_vector_quartz = eC +(-2.0*(eC.Dot(en))*en);
                
                Vector RefInter = name.MCP_end_Intersection_ref(Reflection_vector_quartz, Quartz_wall_Intersection);
                
                if ((  (RefInter.GetX()!= 0) && (RefInter.GetY()!= 0)  )) {
                    Last_wall_intersection = Quartz_wall_Intersection;
                    return_intersection_vector =Reflection_vector_quartz;
                    
                    return return_intersection_vector;
                }
                else{
                    
                    return_intersection_vector = Recurssiv_wall_Intersection(Quartz_wall_Intersection, Reflection_vector_quartz, name, Last_wall_intersection);
                    
                    return return_intersection_vector;
                }
            }
        }
    }
    return return_intersection_vector;
}


//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
void MCP::Cherenkov_Photon_Generation( Particle& pname, Vector &first_intersection, Vector &second_intersection, vector<Vector> &pathPhoton_array, vector<Vector> &CherenkovPhoton_array)
{ //The first_intersection is the position where the particle first intersects the mcp and the second_intersection is the end intersection point of the particle.
    // Now we generates cherenkov photons along the particle path.//
    
    Vector particle_momentum = pname.GetMomentum();
    double particle_energy =  pname.GetEnergy();
    
    Vector velocity_vector = pname.GetVelocity();
    
    
    double v = velocity_vector.Length()* copysign(1.0,particle_momentum.GetZ());
    
    if ((1.0/velocity_vector.Length()) < (1.0/MCP_Quartz_Reflective_index*1.0)) {
        cout<< "return "<< "b = " << (1.0/ velocity_vector.Length()) <<  " n^-1" << (1.0/MCP_Quartz_Reflective_index*1.0) << endl;
        return;
    }
    
    double distance_travelled = (first_intersection-second_intersection).Length();
    
    
    
    double dE = (1.0/(300.0*pow(10.0,-9.0)))    -       (1.0/(500.0*pow(10.0,-9.0)));
    
    
    int Number_particles =abs( dE*2.0*M_PI*( 1.0 - (1.0/ pow((MCP_Quartz_Reflective_index*(1.0/v)),2.0))) *distance_travelled*(1.0/137.0));
    
    //cout<< "Number_particles" << Number_particles<<endl;
    
    // The end_value is a constant, which describes the position where the particle is at the end of the quart plan.
    double end_value =((second_intersection.GetZ() - first_intersection.GetZ()) / particle_momentum.GetZ());
    
    // Generating the path for the particles where the photon should emits as Cherenkov light.
    
    srand((unsigned)time(0) );
    for (int n=0; n<Number_particles; n++) // Cherenkov photons per particle
    {
        
        
        double random_number_d =  (double(rand())/(RAND_MAX))*end_value; // generate numbers between 0 and end_value.
        
        double path_photon_x =first_intersection.GetX() + random_number_d*particle_momentum.GetX();
        double path_photon_y =first_intersection.GetY() + random_number_d*particle_momentum.GetY();
        double path_photon_z =first_intersection.GetZ() + random_number_d*particle_momentum.GetZ();
        
        
        Vector path_photon(path_photon_x,path_photon_y,path_photon_z);
        
        double random_angle = ((double)rand()/RAND_MAX)* (2.0*M_PI); //Generate number between 0 and 2*Pi
        
        
        double Cherenkov_Angle = acos((1.0)/(MCP_Quartz_Reflective_index*v));
        double x = (pname.GetMomentum()).Length()*sin(Cherenkov_Angle)*cos(random_angle); //random angle is from 0-2pi
        double y = (pname.GetMomentum()).Length()*sin(Cherenkov_Angle)*sin(random_angle);
        double z = (pname.GetMomentum()).Length()*cos(Cherenkov_Angle);
        Vector Cherenkov_photon(x,y,z);
        
        pathPhoton_array.push_back(path_photon);
        CherenkovPhoton_array.push_back(Cherenkov_photon);
    }
    return;
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

Vector MCP::MCP_end_Intersection_ref(Vector &Photon_momentum, Vector &particle_pos)
{
    // Intersection with the end surface of the quartz//
    
    // end_position is the position of the end plan of the quartz.
    Vector front_quartz_position = MCP_Position;
    
    Vector end_position(front_quartz_position.GetX(),front_quartz_position.GetY(),front_quartz_position.GetZ()+(MCP_Dimension.GetZ()*MCP_end_pos_sign)); //0.020*
    
    Vector normal_end_quartz(0, 0, MCP_QuartDimension.GetX()*MCP_QuartDimension.GetY()*MCP_z_axis_sign);
    
    double dot_end_quartz = Photon_momentum.Dot(normal_end_quartz);
    
    if (dot_end_quartz > 0.0)
    {
        return Vector(0,0,0);
        
    }
    
    // inter_end_quartz is the intersection point on the end surface of the quartz
    Vector inter_end_quartz = particle_pos-(((particle_pos-end_position).Dot(normal_end_quartz)/dot_end_quartz)*Photon_momentum);
    
    Vector w_end_quartz = inter_end_quartz - end_position;
    
    
    // If the intersection of the plane falls outside window, return a
    // Zero vector.
    if ((abs(w_end_quartz.GetX()) > (MCP_QuartDimension.GetX()/2.0)) ||( abs(w_end_quartz.GetY()) > (MCP_QuartDimension.GetY()/2.0)))
    {
        return Vector(0,0,0);
    }
    
    return inter_end_quartz;
}

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
Vector MCP::MCP_end_Intersection_for_particle( Particle& pname, Vector &particle_pos, MCP* name)
{
    // Intersection with the end surface of the quartz//
    
    // end_position is the position of the end plan of the quartz.
    Vector front_quartz_position = name->MCP_Position;
    
    Vector end_position(front_quartz_position.GetX(),front_quartz_position.GetY(),front_quartz_position.GetZ()+(MCP_QuartDimension.GetZ()*(1.0*MCP_end_pos_sign)));
    
    Vector particle_momentum = pname.GetMomentum();
    
    Vector normal_end_quartz(0, 0, MCP_QuartDimension.GetX()*MCP_QuartDimension.GetY()*MCP_z_axis_sign);
    
    double dot_end_quartz = particle_momentum.Dot(normal_end_quartz);
    
    if (dot_end_quartz > 10e-18)
    {
        
        return Vector(0,0,0);
    }
    
    // inter_end_quartz is the intersection point on the end surface of the quartz
    Vector inter_end_quartz = particle_pos-(((particle_pos-end_position).Dot(normal_end_quartz)/dot_end_quartz)*particle_momentum);
    
    Vector w_end_quartz = inter_end_quartz - end_position;
    
    // If the intersection of the plane falls outside window, return a
    // Zero vector.
    if ( (abs(w_end_quartz.GetX()) > ( MCP_QuartDimension.GetX())/2.0) || ( abs(w_end_quartz.GetY()) > (MCP_QuartDimension.GetY())/2.0))
    {
        
        return Vector(0,0,0);
    }
    
    return inter_end_quartz;
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------


Vector MCP::MCP_end_Intersection_for_Photon(Vector &Photon_momentum, Vector &particle_pos, MCP* name)
{
    
    if (    (Photon_momentum.GetX()==0 && Photon_momentum.GetY()==0 )    ||  (Photon_momentum.GetX()!=Photon_momentum.GetX() && Photon_momentum.GetY()!=Photon_momentum.GetY() )               ) {
        return Vector(0,0,0);
    }
    // Intersection with the end surface of the quartz//
    
    // end_position is the position of the end plan of the quartz.
    Vector front_quartz_position = name->MCP_Position;
    
    Vector end_position(front_quartz_position.GetX(),front_quartz_position.GetY(),front_quartz_position.GetZ()+(MCP_QuartDimension.GetZ()*MCP_end_pos_sign));
    
    Vector normal_end_quartz(0, 0, MCP_QuartDimension.GetX()*MCP_QuartDimension.GetY()*MCP_z_axis_sign);
    
    double dot_end_quartz = Photon_momentum.Dot(normal_end_quartz);
    
    if (dot_end_quartz > 0.0)
    {
        return Vector(0,0,0);
    }
    
    // inter_end_quartz is the intersection point on the end surface of the quartz
    Vector inter_end_quartz = particle_pos-(((particle_pos-end_position).Dot(normal_end_quartz)/dot_end_quartz)*Photon_momentum);
    
    Vector w_end_quartz = inter_end_quartz - end_position;
    
    // If the intersection of the plane falls outside window, return a
    // Zero vector.
    if ((abs(w_end_quartz.GetX()) > (MCP_QuartDimension.GetX()/2.0)) ||( abs(w_end_quartz.GetY()) > (MCP_QuartDimension.GetY()/2.0)))
    {
        return Vector(0,0,0);
    }
    return inter_end_quartz;
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

void MCP::MCP_Intersection(Particle& pname, Vector &particle_pos, MCP* name, TH2F *Intersection_plot, int pore_ON_OFF, vector<Vector> &pore_positions, vector<Vector> &event_plan_array, int hist_ON_OFF)
{
    
    //Intersection with the front surface of the MCP which is the quartz.//
    
    Vector particle_momentum = pname.GetMomentum(); //Momentum of the particle
    
    Vector front_quartz_position = name->MCP_Position; // position of the MCP on the FIT  (Husk nedeunder MCP flad position)
    
    Vector normal_front_quartz(0,0,MCP_QuartDimension.GetY()*MCP_QuartDimension.GetX()*MCP_z_axis_sign);
    
    double dot_front_quartz = particle_momentum.Dot(normal_front_quartz);
    
    if (dot_front_quartz > 10e-18)
    {
        return;
    }
    
    // inter_front_quartz is the intersection point on the quartz
    Vector inter_front_quartz =  particle_pos -((( particle_pos - front_quartz_position).Dot(normal_front_quartz)/dot_front_quartz)*particle_momentum);
    
    Vector w = inter_front_quartz - front_quartz_position;
    
    // If the intersection of the plane falls outside window, return a
    
    if ((abs(w.GetX()) > (name->MCP_QuartDimension.GetX()/2.0)) || (abs(w.GetY()) > (name->MCP_QuartDimension.GetY()/2.0)))
    {
        return;
    }
    
    //---------------------------------------------------------------------
    vector<Vector> pathPhoton_array;
    vector<Vector> CherenkovPhoton_array;
    
    
    
    //---------------------------------------------------------------------
    
    // Intersection with the end://
    // Here it checks that it will hit the walls of the quartz first or the end panel of the quartz:
    int ref = 0;
    int dir = 0;
    //The particle intersection with the MCP directly
    Vector intersection_with_end_MCP = MCP_end_Intersection_for_particle(pname, particle_pos, name);
    //cout << "intersection_with_end_MCP"<< intersection_with_end_MCP<<endl;
    
    Vector zero_vector(0,0,0);
    Vector Last_wall_intersection;
    
    if (intersection_with_end_MCP.GetX()==0 && intersection_with_end_MCP.GetY()==0){
        Vector zero_vector(0,0,0);
        Vector particle_intersection_with_wall = Quartz_wall_Intersection(name, Last_wall_intersection, inter_front_quartz, particle_momentum, 0, zero_vector);
        
        if (particle_intersection_with_wall.GetX() != particle_intersection_with_wall.GetX() ) {
            
            return;}
        
        Cherenkov_Photon_Generation(pname, inter_front_quartz, particle_intersection_with_wall, pathPhoton_array, CherenkovPhoton_array);
        
        for (unsigned int i = 0; i < pathPhoton_array.size(); i++)
        {
            
            Vector photon_intersection_with_wall = Quartz_wall_Intersection(name, Last_wall_intersection, pathPhoton_array[i], CherenkovPhoton_array[i], 0, zero_vector);
            
            if ( ((photon_intersection_with_wall.GetX()!= 0) && (photon_intersection_with_wall.GetY()!= 0) && (photon_intersection_with_wall.GetZ()!= 0)) &&  ((TMath::IsNaN(photon_intersection_with_wall.GetX())!= 1) && (TMath::IsNaN(photon_intersection_with_wall.GetX())!= 1))  ){
                //If the photon interc. with the wall:
                Vector reflected_photon = Quartz_wall_Intersection(name, Last_wall_intersection, pathPhoton_array[i] , CherenkovPhoton_array[i], 1, CherenkovPhoton_array[i]);
                
                
                Vector intersection_with_slut_MCP = MCP_end_Intersection_for_Photon(reflected_photon, Last_wall_intersection, name);
                pore_intersection(pore_ON_OFF, pore_positions, intersection_with_slut_MCP, Intersection_plot, event_plan_array, hist_ON_OFF);
            }
            
            else if ((TMath::IsNaN(photon_intersection_with_wall.GetX())!= 1) && (TMath::IsNaN(photon_intersection_with_wall.GetY())!= 1)){
                Vector intersection_with_slut_MCP = MCP_end_Intersection_for_Photon(CherenkovPhoton_array[i], pathPhoton_array[i], name);
                pore_intersection(pore_ON_OFF, pore_positions, intersection_with_slut_MCP, Intersection_plot, event_plan_array, hist_ON_OFF);
            }
        }
    }
    
    
    else {
        //Generate Cherenkov photon:
        Cherenkov_Photon_Generation(pname, inter_front_quartz, intersection_with_end_MCP, pathPhoton_array, CherenkovPhoton_array);
        
        Vector zero_vector(0,0,0);
        Vector Last_wall_intersection;
        
        for (unsigned int i = 0; i < pathPhoton_array.size(); i++)
        {
            Vector photon_intersection_with_wall = Quartz_wall_Intersection(name, Last_wall_intersection, pathPhoton_array[i], CherenkovPhoton_array[i], 0, zero_vector);
            
            
            
            //cout << "photon_intersection_with_wall"<< photon_intersection_with_wall<<endl;
            
            if ( ((photon_intersection_with_wall.GetX()!= 0) && (photon_intersection_with_wall.GetY()!= 0) && (photon_intersection_with_wall.GetZ()!= 0)) &&  ((TMath::IsNaN(photon_intersection_with_wall.GetX())!= 1) && (TMath::IsNaN(photon_intersection_with_wall.GetY())!= 1))  )
            {  ref+=1;                 //If the photon interc. with the wall:
                // Here the reflected photon is found:
                
                Vector reflected_photon = Quartz_wall_Intersection(name, Last_wall_intersection, pathPhoton_array[i] , CherenkovPhoton_array[i], 1, CherenkovPhoton_array[i]);
                
                
                Vector intersection_with_slut_MCP = MCP_end_Intersection_for_Photon(reflected_photon, Last_wall_intersection, name);
                
                pore_intersection(pore_ON_OFF, pore_positions, intersection_with_slut_MCP, Intersection_plot, event_plan_array,hist_ON_OFF);
            }
            
            else if ((TMath::IsNaN(photon_intersection_with_wall.GetX())!= 1) && (TMath::IsNaN(photon_intersection_with_wall.GetY())!= hist_ON_OFF))
            { dir+=1;
                Vector intersection_with_slut_MCP = MCP_end_Intersection_for_Photon(CherenkovPhoton_array[i], pathPhoton_array[i], name);
                pore_intersection(pore_ON_OFF, pore_positions, intersection_with_slut_MCP, Intersection_plot, event_plan_array,hist_ON_OFF);
            }
            
        }
  
    }
    pathPhoton_array.clear();
    CherenkovPhoton_array.clear();
    
    
    
    
}
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>


Particle Make_Particle(TParticle* q){
    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
    string name = (q)->GetName();
    double x = (q)->Px();
    double y = (q)->Py();
    double z = (q)->Pz();
    double vx = (q)->Vx();
    double vy = (q)->Vy();
    double vz = (q)->Vz();
    double Mass =(q)->GetMass();
    
    TParticlePDG* pdg = pdgDb->GetParticle(q->GetPdgCode());
    
    double charge =(pdg->Charge())*(1.0/3.0);
    //cout << "charge: " << charge << endl;
    Vector momentum(x,y,z);
    Vector Vertex(vx,vy,vz);
    Particle in(name, Mass, charge);
    in.SetMomentum(momentum);
    in.SetVertex(Vertex);
    
    return in;
}

#include <TTree.h>
#include <TClonesArray.h>
#include <TSystem.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2.h>
#include <TF1.h>
#include <TVector3.h>
#include <TFormula.h>
#include <TParameter.h>
typedef std::map< int,vector<double> > mapvector;


///--------------------------------------------------------
///                      Lee-Yang zeros method
#include <TComplex.h>
#include <TNtuple.h>

void Q_theta(vector<Vector> &particles_intersection_mcp, int n_order, double parameter, string TypeParameter, Tmap *LYZ_Parameter_Q, Tmap *LYZ_Parameter_M){
    
    int M = particles_intersection_mcp.size();
    double Qx = 0.0;
    double Qy = 0.0;
    for (int i = 0; i < M; i++) {
        double phi    = TMath::ATan2(particles_intersection_mcp[i].GetY(),particles_intersection_mcp[i].GetX());
        Qx+= TMath::Cos(n_order * phi);
        Qy+= TMath::Sin(n_order * phi);
    }
    
    
    
    int MapParameter = 0.0;
    if (TypeParameter == "c") {
        double percent_0 = 0.0;
        double percent_1 = 0.0;
        int PercentParameter = 0;
        for (int i=0; i<33; i++) {
            percent_1+=3;
            PercentParameter += 3;
            if (parameter>percent_0 && parameter<= percent_1 ) {
                break;
            }
            percent_0 = percent_1;
        }
        MapParameter = PercentParameter;
    }
    else if(TypeParameter == " "){
        MapParameter = parameter;
    }
    
    
    
     LYZ_Parameter_Q->TAdd(MapParameter, Qx);
 
    
    
    
    /// Info during PROOF:
    TString msgde = TString::Format("Value = %f" ,LYZ_Parameter_Q->GetTmap()->at(3)->at(0) );
     if (gProofServ) gProofServ->SendAsynMessage(msgde);
    
}


void LYZ(Tmap *LYZ_Parameter_Q, Tmap *LYZ_Parameter_M, TH2F *histogram, Option_t * option ){
    
    //TString tlv = TString::Format("on on on  = %d" , 12);
    //if (gProofServ) gProofServ->SendAsynMessage(tlv);
    
    /* TString msgde = TString::Format("on on on  = %f" , LYZ_Parameter_Q->GetTmap()->at(3)->at(0));
     if (gProofServ) gProofServ->SendAsynMessage(msgde);*/
    
    
    /*for (map<int,vector<double>>::iterator it=LYZ_Parameter_Q->GetTmap().begin(); it!=LYZ_Parameter_Q->GetTmap().end(); ++it){
     
     
     
     int N_event = LYZ_Parameter_Q->GetTmap()[it->first].size();
     
     
     
     double M_avg=0;
     double M_sum=0;
     for (int i = 0; i < N_event; i++) {
     M_sum += LYZ_Parameter_M->GetTmap()[it->first][i];
     }
     M_avg = M_sum/N_event;
     
     string equation = "[0]*TComplex::Abs(";
     for (int i = 0; i < N_event; i++) {
     double Q = LYZ_Parameter_Q->GetTmap()[it->first][i];
     equation+="+TComplex(cos((";
     equation+=to_string(Q);
     equation+= ")*x),sin((";
     equation+=to_string(Q);
     equation+=")*x))";
     }
     equation += ")";
     const char *cstr = equation.c_str();
     
     
     
     TF1 *G = new TF1("G",cstr,0, 0.2);
     G->SetParameter(0,(1.0/N_event*1.0));
     G->Draw();
     cout << "cstr = "<< cstr <<endl;
     double r_min1 = G->GetMinimumX(0, 0.01, 1.E-10, 100, false);
     
     double v = 2.40483/(r_min1*M_avg);
     
     double r_min2 = G->GetX (0, 0, 0.01, 1.E-10, 100, false);
     //cout << "equation = "<< equation <<endl;
     //cout << "r_min1 = "<< r_min1 <<endl;
     //cout << "r_min2 = "<< r_min2 <<endl;
     histogram->Fill(it->first, v);
     //delete G;
     }
     histogram->SetMarkerStyle(21);
     histogram->SetMarkerSize(0.5);
     histogram->SetMarkerColor(3);
     histogram->SetXTitle("Centrality [%]");
     histogram->SetYTitle(" v_{2} ");
     // histogram->Draw(option);*/
    
}
#endif