Skip to content

Commit

Permalink
Added vertex information to track projections in StFcsDb (star-bnl#565)
Browse files Browse the repository at this point in the history 
Changed functions in StFcsDb that are related to projecting a g2t track
onto the FCS detector plane and now can also take vertex information
from g2t vertex.
  • Loading branch information
@dkapukchyan
dkapukchyan authored Aug 11, 2023
1 parent d14cd9d commit 3a748c7
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 55 deletions.
95 changes: 59 additions & 36 deletions StRoot/StFcsDbMaker/StFcsDb.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -706,70 +706,93 @@ double StFcsDb::getProjectedDistance(StFcsPoint* ecal, StFcsCluster* hcal, doub
return sqrt(dX*dX + dY*dY);
};

StThreeVectorD StFcsDb::projectTrackToEcal(const g2t_track_st* g2ttrk) const
StThreeVectorD StFcsDb::getNormal(int det) const
{
double phi = atan2(g2ttrk->p[1],g2ttrk->p[0]);
double theta = 2.0*atan(exp(-1.0*g2ttrk->eta));
if( g2ttrk->p[0]<0 ){ return projectToDet(0,theta,phi); } // North side for negative px
else{ return projectToDet(1,theta,phi); } // South side for positive px, since px==0 does not hit detector choose south side for that case
double detangle = getDetectorAngle(det)*M_PI/180.0;
if( det%2==0 ){ detangle *= -1.0; } //North side use negative angle
return StThreeVectorD( sin(detangle), 0 ,cos(detangle) );
double planenormal[3] = {sin(detangle),0,cos(detangle)};
}

StThreeVectorD StFcsDb::projectTrackToHcal(const g2t_track_st* g2ttrk) const
StThreeVectorD StFcsDb::projectTrackToEcal(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert) const
{
double phi = atan2(g2ttrk->p[1],g2ttrk->p[0]);
double theta = 2.0*atan(exp(-1.0*g2ttrk->eta));
if( g2ttrk->p[0]<0 ){ return projectToDet(2,theta,phi); } // North side for negative px
else{ return projectToDet(3,theta,phi); } // South side for positive px, since px==0 does not hit detector choose south side for that case
int det=0; // North side for negative px
// South side for positive px, since px==0 does not hit detector choose south side for that case
if( g2ttrk->p[2]>=0 && g2ttrk->p[0]>=0 ){ det=1; }
if( g2ttrk->p[2]<0 && g2ttrk->p[0]<0 ){ det=1; }
return projectTrack(det,g2ttrk,g2tvert,0);
}

StThreeVectorD StFcsDb::projectTrackToEcalSMax(const g2t_track_st* g2ttrk) const
StThreeVectorD StFcsDb::projectTrackToHcal(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert) const
{
double phi = atan2(g2ttrk->p[1],g2ttrk->p[0]);
double theta = 2.0*atan(exp(-1.0*g2ttrk->eta));
if( g2ttrk->p[0]<0 ){ return projectToShowerMax(0,theta,phi); } // North side for negative px
else{ return projectToShowerMax(1,theta,phi); } // South side for positive px, since px==0 does not hit detector choose south side for that case
int det = 2; // North side for negative px
// South side for positive px, since px==0 does not hit detector choose south side for that case
if( g2ttrk->p[2]>=0 && g2ttrk->p[0]>=0 ){ det=3; }
if( g2ttrk->p[2]<0 && g2ttrk->p[0]<0 ){ det=3; }
return projectTrack(det,g2ttrk,g2tvert,0);
}

StThreeVectorD StFcsDb::projectTrackToHcalSMax(const g2t_track_st* g2ttrk) const
StThreeVectorD StFcsDb::projectTrackToEcalSMax(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert) const
{
double phi = atan2(g2ttrk->p[1],g2ttrk->p[0]);
double theta = 2.0*atan(exp(-1.0*g2ttrk->eta));
if( g2ttrk->p[0]<0 ){ return projectToShowerMax(2,theta,phi); } // North side for negative px
else{ return projectToShowerMax(3,theta,phi); } // South side for positive px, since px==0 does not hit detector choose south side for that case
int det=0; // North side for negative px
// South side for positive px, since px==0 does not hit detector choose south side for that case
if( g2ttrk->p[2]>=0 && g2ttrk->p[0]>=0 ){ det=1; }
if( g2ttrk->p[2]<0 && g2ttrk->p[0]<0 ){ det=1; }
return projectTrack(det,g2ttrk,g2tvert,-1);
}

StThreeVectorD StFcsDb::projectToDet(int det, double azimuth, double polar, double xvertex, double yvertex, double zvertex) const
StThreeVectorD StFcsDb::projectTrackToHcalSMax(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert) const
{
return projectLine(det,azimuth,polar,0,xvertex,yvertex,zvertex);
int det = 2; // North side for negative px
// South side for positive px, since px==0 does not hit detector choose south side for that case
if( g2ttrk->p[2]>=0 && g2ttrk->p[0]>=0 ){ det=3; }
if( g2ttrk->p[2]<0 && g2ttrk->p[0]<0 ){ det=3; }
return projectTrack(det,g2ttrk,g2tvert,-1);
}

StThreeVectorD StFcsDb::projectToShowerMax(int det, double azimuth, double polar, double xvertex, double yvertex, double zvertex) const
StThreeVectorD StFcsDb::projectTrack(int det, const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert, double showermaxz) const
{
return projectLine(det,azimuth,polar,-1,xvertex,yvertex,zvertex);
double linedir[3] = {g2ttrk->p[0],g2ttrk->p[1],g2ttrk->p[2]};
if( g2tvert!=0 ){
int vertind = g2ttrk->start_vertex_p - 1; //To correct for offset by one between fortran array and c array. 0 start index means it was generated at the starting vertex
std::cout << "+++++ DEBUG: vertind = " << vertind << " +++++" << std::endl;
double linestart[3] = {g2tvert[vertind].ge_x[0],g2tvert[vertind].ge_x[1],g2tvert[vertind].ge_x[2]};
if( vertind >= 0 ){//Since start index==0 means no start then vertind<0 will default to using origin
return projectLine(det, linedir, linestart, showermaxz);
}
}
double zero[3] = {0,0,0};
return projectLine(det,linedir,zero,showermaxz);
}

StThreeVectorD StFcsDb::getNormal(int det) const
StThreeVectorD StFcsDb::projectLine(int det, StThreeVectorD &linedirection, StThreeVectorD &lineorigin, double showermaxz) const
{
double detangle = getDetectorAngle(det)*M_PI/180.0;
if( det%2==0 ){ detangle *= -1.0; } //North side use negative angle
return StThreeVectorD( sin(detangle), 0 ,cos(detangle) );
double planenormal[3] = {sin(detangle),0,cos(detangle)};
double linedir[3] = {linedirection.x(),linedirection.y(),linedirection.z()};
double linestart[3] = {lineorigin.x(),lineorigin.y(),lineorigin.z()};
return projectLine(det,linedir,linestart,showermaxz);
}

StThreeVectorD StFcsDb::projectLine(int det, double azimuth, double polar, double showermaxz, double xvertex, double yvertex, double zvertex) const
StThreeVectorD StFcsDb::projectLine(int det, double* linedirection, double* lineorigin, double showermaxz) const
{
if( showermaxz<0 ){ showermaxz = getShowerMaxZ(det); } //when negative use default showermax
if( det%2==0 ){ //North side is negative x-axis
if( linedirection[2]>=0 && linedirection[0]>=0 ){ LOG_WARN << "Incorrect Det" << endm; }
if( linedirection[2]<0 && linedirection[0]<=0 ){ LOG_WARN << "Incorrect Det" << endm; }
}
else{ //South side is positive x-axis
if( linedirection[2]>=0 && linedirection[0]<0 ){ LOG_WARN << "Incorrect Det" << endm; }
if( linedirection[2]<0 && linedirection[0]>=0 ){ LOG_WARN << "Incorrect Det" << endm; } //(If x and z direction are both negative then also projects to south side which is positive x-axis
}
double detangle = getDetectorAngle(det)*M_PI/180.0;
if( det%2==0 ){ detangle *= -1.0; } //North side use negative angle
StThreeVectorD xyzoff = getDetectorOffset(det,showermaxz);
StThreeVectorD planenormal = getNormal(det);
double linedir[3] = {cos(polar)*sin(azimuth),sin(polar)*sin(azimuth),cos(azimuth)};//This is the direction of a line formed for a given azimuth and polar angle in xyx coordinates where line starts at origin
//Solution of intersection of line and plane where line starts at origin and plane has some normal and has a point at the detector offset, "t" is the free parameter in the parametric equation of the line.
//Solution of intersection of line and plane where line has direction {xdir,ydir,zdir}*t and starts at {xorigin,yorigin,zorigin} and a plane that has some normal with a point on the plane as the detector offset; "t" is the free parameter in the parametric equation of the line.
double tintersection =
(planenormal.x()*(xyzoff.x()-xvertex)+planenormal.y()*(xyzoff.y()-yvertex)+planenormal.z()*(xyzoff.z()-zvertex)) /
(planenormal.x()*linedir[0]+planenormal.y()*linedir[1]+planenormal.z()*linedir[2]);
return StThreeVectorD(linedir[0]*tintersection,linedir[1]*tintersection,linedir[2]*tintersection);
(planenormal.x()*(xyzoff.x()-lineorigin[0])+planenormal.y()*(xyzoff.y()-lineorigin[1])+planenormal.z()*(xyzoff.z()-lineorigin[2])) /
(planenormal.x()*linedirection[0]+planenormal.y()*linedirection[1]+planenormal.z()*linedirection[2]);

return StThreeVectorD( linedirection[0]*tintersection+lineorigin[0], linedirection[1]*tintersection+lineorigin[1], linedirection[2]*tintersection+lineorigin[2] );
}

//! get coordinates of center of the cell STAR frame from StFcsHit
Expand Down
36 changes: 17 additions & 19 deletions StRoot/StFcsDbMaker/StFcsDb.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,7 @@
#include "tables/St_fcsPresValley_Table.h"
#include "tables/St_vertexSeed_Table.h"
#include "tables/St_g2t_track_Table.h"
#include "tables/St_g2t_vertex_Table.h"
class StFcsHit;
class StFcsCluster;
class StFcsPoint;
Expand Down Expand Up @@ -282,28 +283,25 @@ class StFcsDb : public TDataSet {
const g2t_track_st* getPrimaryG2tTrack(StFcsHit* h, g2t_track_st* g2ttrk, float& fraction, int& ntrk, unsigned int order=0);
const g2t_track_st* getPrimaryG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, float& fraction, int& ntrk, unsigned int order=0);

StThreeVectorD projectTrackToEcal(const g2t_track_st* g2ttrk) const;
StThreeVectorD projectTrackToHcal(const g2t_track_st* g2ttrk) const;
StThreeVectorD projectTrackToEcalSMax(const g2t_track_st* g2ttrk) const; //!< SMax = Shower Max Z
StThreeVectorD projectTrackToHcalSMax(const g2t_track_st* g2ttrk) const; //!< SMax = Shower Max Z
StThreeVectorD projectToDet(int det,double azimuth, double polar, double xvertex=0, double yvertex=0, double zvertex=0) const;
StThreeVectorD projectToShowerMax(int det,double azimuth, double polar, double xvertex=0, double yvertex=0, double zvertex=0) const;
StThreeVectorD projectTrackToEcal(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert=0) const;//!< project a g2t track to Ecal with a given track and vertex. If no vertex given assume a vertex of (0,0,0)
StThreeVectorD projectTrackToHcal(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert=0) const;//!< project a g2t track to Hcal with a given track and vertex. If no vertex given assume a vertex of (0,0,0)
StThreeVectorD projectTrackToEcalSMax(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert=0) const; //!< SMax = Shower Max Z
StThreeVectorD projectTrackToHcalSMax(const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert=0) const; //!< SMax = Shower Max Z
StThreeVectorD projectTrack(int det, const g2t_track_st* g2ttrk, const g2t_vertex_st* g2tvert, double showermaxz=-1) const;//!< Generic g2t track projection function but #det and #showermaxz needs to be specified; if #det or #showermaxz not known use corresponding #projectTrackToEcal(), #projectTrackToHcal(), #projectTrackToEcalSMax, #projectTrackToHcalSMax instead
StThreeVectorD projectLine(int det, StThreeVectorD &linedirection, StThreeVectorD &lineorigin, double showermaxz=-1) const;//!< Like #projectLine(det, double*, double*, double) except use StThreeVectorD for line direction and origin
/**@brief XYZ of a projected line to the FCS detector plane
Get the STAR XYZ that corresponds to the intersection of the FCS plane and a line given by the direction of some azimuthal angle and polar angle. You can also specify the origin of the line in the global STAR XYZ coordinates.
Note: if you choose the wrong detector for a given polar angle you will get the wrong projection so please use #projectTrackToEcal(), #projectTrackToHcal(), #projectTrackToEcalSMax(), or #projectTrackToHcalSMax() to take into account which polar angle goes to which detector.
@param det detector id to project to (needed for correct angle).
@param azimuth azimuthal angle (angle from positive z-axis) of the line in radians
@param polar polar angle (angle in x-y plane) of the line in radians
@param zhowermax depth in z of the detector to project the line to. Negative values use #getShowerMaxZ()
@param xvertex STAR global x coordinate of the origin of the line
@param yvertex STAR global y coordinate of the origin of the line
@param zvertex STAR global z coordinate of the origin of the line
*/
StThreeVectorD projectLine(int det, double azimuth, double polar, double zshowermax=-1, double xvertex=0, double yvertex=0, double zvertex=0) const;

Note: if you choose the wrong detector for a given polar angle you will get the wrong projection so please use #projectTrackToEcal(), #projectTrackToHcal(), #projectTrackToEcalSMax(), or #projectTrackToHcalSMax() to take into account which polar angle goes to which detector.
@param det detector id to project line to (needed for correct angle).
@param linedirection size 3 array indicating the direction of the line to project onto #det, first element x-direction, second element y-direction, third element z-direction
@param lineorigin size 3 array indication the origin of the line to project onto #det, first element is x-coordinate, second element is y-coordinate, third element is z-coordinate
@param zhowermax depth in z of the detector to project the line to. Any negative value will use #getShowerMaxZ()
*/
StThreeVectorD projectLine(int det, double* linedirection, double* lineorigin, double showermaxz=-1) const;

private:
int mDbAccess=1; //! enable(1) or disabe(0) DB access
int mRun=0; //! run#
Expand Down

0 comments on commit 3a748c7

Please sign in to comment.