#ifndef _VECTOR_H #define _VECTOR_H #include #include // =================================================================== 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 #include ///------------------------------------------------------------------- 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 #include #include #include "TObject.h" class Tmap : public TObject { private: std::map< int,std::vector *> *fQ; std::vector* fQContainer; public : Tmap(void); virtual ~Tmap(){}; virtual void TAdd(int key,double Qnumber); virtual bool TExist(int key); virtual std::map< int,std::vector* >* GetTmap(void) const{ return fQ; } ClassDef(Tmap,2) }; ClassImp(Tmap) Tmap::Tmap(void){ fQ = new map*>(); fQContainer = new vector(); } void Tmap::TAdd(int key,double Qnumber){ fQContainer->push_back(Qnumber); fQ->insert(std::pair*>(key, fQContainer)); } bool Tmap::TExist(int key){ bool exist; if ( fQ->find(key) == fQ->end()){ return false; } else{ return true;} } #endif // =================================================================== #include 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 #include ///------------------------------------------------------------------- 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 &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 temp(position_array_quartz); position_array_quartz.clear(); if (number_of_sub_quartz == 4) { for (unsigned int i=0; i &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 &position_array_quartz, const double &heigth, 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 #include #include #include #include "TROOT.h" #include "TRandom3.h" #include "TH2F.h" #include "TCanvas.h" #include #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 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 &Angels_array); void MCP_quartz_div(vector &Single_position_quartz_array); void pore_intersection(int pore_ON_OFF, vector &pore_positions, Vector &intersect_point, TH2F *Intersection_plot, 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 &pathPhoton_array, vector &CherenkovPhoton_array); void MCP_Intersection(Particle& pname, Vector &particle_pos, MCP* name , TH2F *Intersection_plot, int pore_ON_OFF, vector &pore_positions, 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 &Angels_array){ for (unsigned int i=0; i0) { end_sign =-1.0; } else { end_sign =+1.0; } MCP_end_pos_sign = end_sign; } MCP::~MCP(){ } //------------------------------------------------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------------------------------------------------ void MCP::MCP_quartz_div(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 &pore_positions, Vector &intersect_point, TH2F *Intersection_plot, 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()<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 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 <--| // | ^ | 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 &pathPhoton_array, 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< 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 &pore_positions, 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 pathPhoton_array; 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< #include #include 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 #include #include #include #include #include #include #include #include #include typedef std::map< int,vector > mapvector; ///-------------------------------------------------------- /// Lee-Yang zeros method #include #include void Q_theta(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>::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 <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 <SetYTitle(" v_{2} "); // histogram->Draw(option);*/ } #endif