CgemGeoReadoutPlane.h

Go to the documentation of this file.

#ifndef CGEMGEOREADOUTPLANE_H
#define CGEMGEOREADOUTPLANE_H 
 
#include <string>
#include <vector>
#include <CLHEP/Units/PhysicalConstants.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Vector/ThreeVector.h>
#include "G4ThreeVector.hh"
 
using namespace std;
using namespace CLHEP;
 
class CgemGeoReadoutPlane{
 
    public:
        CgemGeoReadoutPlane(int iLayer, int iSheet,
                double rx, double rv, 
                double phi_min, double dX_strip, double dV,
                double w, double zmin, double L,
                int NXStrip, int NVStrip,
                double Xpitch, double XstripWidth, double Vpitch, double VstripWidth, double stereoAngle, 
                double midROfGap, double outROfGap);
        ~CgemGeoReadoutPlane();
 
        void print();
        
        int    getLayerId()     {return m_iLayer;};
        int    getSheetId()     {return m_iSheet;};
        double getRX()          {return m_RX;};
        double getRV()          {return m_RV;};
        double getPhimin()      {return m_Phimin;};
        double getXmin()        {return m_Xmin;};
        double getWidth()       {return m_W;};
        double getZmin()        {return m_Zmin;};
        double getLength()      {return m_L;};
        double getXPitch()      {return m_XPitch;};
        double getVPitch()      {return m_VPitch;};
        double getStereoAngle() {return m_StereoAngle;};
        double getMidRAtGap()   {return m_MidRGap;};
        double getVmax()        {return m_Vmax;};
        int    getNXstrips() const { return m_NXstrips; };
        int    getNVstrips() const { return m_NVstrips; };
        
        bool OnThePlane(double phi, double z) const;// check (phi,z) in the sensitive area or not
        
        double getX(double phi) const;// X with m_Xmin as origin
        int  getXStripID(double phi) const;// -1: beyond start, -2: beyond end
        int  getXStripID(double x, double y) const;// -1: beyond start, -2: beyond end, -3: x=y=0
        int  getXStripID(G4ThreeVector pos) const; // -1: beyond start, -2: beyond end, -3: x=y=0
        int  getClosestXStripID(double phi, double& dist) ;       // dist to the border of the closest strips, <0 inside of the strip, >0 outside of the strip
        int  getClosestVStripID(G4ThreeVector pos, double& dist) const;// dist to the border of the closest strips, <0 inside of the strip, >0 outside of the strip
        double  getDist2ClosestXStripCenter(double phi, int& id) ;       // distance to the center of the closest strip
        double  getDist2ClosestVStripCenter(G4ThreeVector pos, int& id) ;// distance to the center of the closest strip
        void getStripID(G4ThreeVector pos, int& X_ID, int& V_ID) const; //  get the closest XID/VID
        void getFiredStripID(G4ThreeVector pos1, G4ThreeVector pos2, vector<int>& vecXID, vector<int>& vecVID) const;
        
        double getVFromLocalXZ(double x,double zeta, bool checkRange=true) const;// input: local X, Zeta @ X plane with (phi_min_strip,-L/2) as origin, output: local V @ V plane
        double getVFromPhiZ(double phi, double z, bool checkRange=true) const;
        double getVFromPhiZ_nearPhiMin(double phi, double z, bool checkRange=true) const;
        double getVInNextSheetFromV(double v, double phiminNext) const { 
            double dPhi = phiminNext-m_Phimin;// assume the next m_Phimin_strip is the same
            while(dPhi<=0) dPhi+=CLHEP::twopi;
            while(dPhi>CLHEP::twopi) dPhi-=CLHEP::twopi;
            //cout<<"getVInNextSheetFromV: dphi = "<<dPhi<<endl;
            double dX = dPhi*m_Rmid;
            return v-dX*cos(m_StereoAngle); 
        }
        int  getVIDFromV(double V) const; // -1: beyond start, -2: beyond end
        int  getVIDInNextSheetFromVID(int vID, double phimin_next) const;
        
        double getPhiMin_strip() const {return m_Phimin_strip+m_Phimin;} 
        double getCentralXFromXID(int X_ID) const 
        {   //X_ID=X_ID-m_NXstrips*m_iSheet; 
            return (X_ID+0.5)*m_XPitch+m_Xmin_strip; 
        }// central X with m_Xmin as origin
        double getPhiFromXID(int X_ID) const; /* from -pi to pi */
        double getCentralVFromVID(int V_ID) const 
        { 
            //V_ID=V_ID-m_NVstrips*m_iSheet; 
            return (V_ID+0.5)*m_VPitch+m_Vmin_strip; 
        }
        double getVStripLength(int V_ID) const
        {
            double vStripLegth(0.);
            if(V_ID>=0&&V_ID<m_NVstrips) 
            {
                double V=getCentralVFromVID(V_ID);
                double Lz=V/fabs(sin(m_StereoAngle));
                if(Lz<=m_L) 
                {
                    vStripLegth=Lz/cos(m_StereoAngle);
                }
                else if(V/cos(m_StereoAngle)<=m_NXstrips*m_XPitch)
                {
                    vStripLegth=m_L/cos(m_StereoAngle);
                }
                else vStripLegth=(m_Vtot-V)/fabs(sin(m_StereoAngle)*cos(m_StereoAngle));
            }
            return vStripLegth;
        }
 
        double getZFromXV(double X, double V, int checkXRange=1, int checkVRange=1) const;// X with m_Xmin as origin
        double getZFromPhiV(double phi, double V, int checkXRange=1) const;
 
        double incidentAngleX(Hep3Vector momentum, double phi) const;
        double incidentAngleV(Hep3Vector momentum, double phi) const;
        double dAngleAcute(double phi1, double phi2) const;
 
        /*      bool intersection(int iX, int iV) */
        /*      { */
        /*    double V=getCentralVFromVID(iV); */
        /*    double Xmin=(V/sin(m_StereoAngle)-m_L); */
        /*      } */
 
    private:
        int    m_iLayer; // layer id: 0, 1 ,2
        int    m_iSheet; // sheet id: 0, 1
        double m_RX;     // radius (mm) for X strips
        double m_RV;     // radius (mm) for V strips
        double m_Rmid;   // middle radius (mm) for anode
        double m_MidRGap;// middle radius of the drift region (mm)
        double m_OutRGap;// outer  radius of the drift region (mm)
        double m_Zmin;   // starting position of the readout plane in Z (mm)
        double m_Zmax;   // ending   position of the readout plane in Z (mm)
        double m_L;      // length in z (mm)
 
        double m_W;           // width in X (mm) excluding gap area (to be connected by bands)
        double m_Phimin;      // starting position of the readout plane in phi (radian) (-pi to pi)
        double m_Xmin;        // starting position of the readout plane in X (mm) = m_Phimin*m_RX 
        double m_Phimin_strip;// position of the first strip in phi (radian) center -0.5*pitch (-pi,pi) with m_Phimin as origin
        double m_Xmin_strip;  // position of the first strip in X (mm) = m_Phimin_strip*m_Rmid  (center-0.5*pitch) with m_Xmin as origin
        double m_Phimax_strip;// position of the last  strip in phi (radian) center +0.5*pitch (-pi,pi) with m_Phimin as origin
        double m_Xmax_strip;  // position of the last  strip in X (mm) = m_Xmin_strip+m_XPitch*m_NXstrips (center+0.5*pitch) with m_Xmin as origin
 
        //double m_Vmin;
        double m_Vtot;         // total V in the sensitive rectangle area (m_XPitch*m_NXstrips, L)
        double m_Vmin_strip;   // position of the first V strip in V (mm) center-0.5*pitch with respect to the corner of the sensitive area
        double m_Vmax;         // position of the last  V strip in V (mm) center+0.5*pitch with respect to the corner of the sensitive area
 
        //double m_Wdead;  // dead area (mm)
 
        double m_XPitch;      // pitch in X (mm)
        double m_VPitch;      // pitch in V (mm)
        double m_XstripWidth; // width of X strips (mm)
        double m_VstripWidth; // width of V strips (mm)
        double m_StereoAngle; // Stereo angle (rad)
        double m_k; // +1 or -1, sign of the stereo angle
 
        int    m_NXstrips;
        int    m_NVstrips;
};
 
inline bool CgemGeoReadoutPlane::OnThePlane(double phi, double z) const
{
    double dPhi=phi-m_Phimin;// [0, 2pi) with m_Phimin as origin
    while(dPhi<0)       dPhi+=CLHEP::twopi;
    while(dPhi>=2*M_PI)  dPhi-=CLHEP::twopi;
    if(dPhi>=m_Phimin_strip && dPhi<=m_Xmax_strip/m_Rmid && z>=m_Zmin && z<=m_Zmax) return true;
    else return false;
}
 
inline double CgemGeoReadoutPlane::getVFromLocalXZ(double x, double zeta, bool checkRange) const 
{
    double V = -9999.0;
    if( (!checkRange) || (x>=0. && x<=m_NXstrips*m_XPitch && zeta>=0 && zeta<=m_L) )
    {
        //double X = x/m_RX*m_RV;// X at V plane
        //if(m_StereoAngle>0) return zeta*sin(m_StereoAngle) + X*cos(m_StereoAngle);
        //else return (zeta-m_L)*sin(m_StereoAngle)+X*cos(m_StereoAngle);
        double z0=m_Zmax; if(m_k<0) z0=m_Zmin;
        //return X*cos(m_StereoAngle)+((m_k-1.0)*m_L/2.0-zeta)*sin(m_StereoAngle);
        V = x*cos(m_StereoAngle)+(z0-zeta-m_Zmin)*sin(m_StereoAngle);
    } 
    return V;
}
 
inline double CgemGeoReadoutPlane::getVFromPhiZ(double phi, double z, bool checkRange) const
{
    double dPhi=phi-m_Phimin;
    while(dPhi<0) dPhi+=CLHEP::twopi;
    while(dPhi>=CLHEP::twopi) dPhi-=CLHEP::twopi;
    double X=dPhi*m_Rmid;// X with m_Xmin as origin
    double x=X-m_Xmin_strip;// x with m_Xmin_strip as origin
    double zeta=z-m_Zmin;// local zeta [0,L]
    return getVFromLocalXZ(x,zeta, checkRange);
}
 
inline double CgemGeoReadoutPlane::getVFromPhiZ_nearPhiMin(double phi, double z, bool checkRange) const
{
    double dPhi=phi-m_Phimin;
    while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
    while(dPhi>=CLHEP::pi) dPhi-=CLHEP::twopi;
    double X=dPhi*m_Rmid;// X with m_Xmin as origin
    double x=X-m_Xmin_strip;// x with m_Xmin_strip as origin
    double zeta=z-m_Zmin;// local zeta [0,L]
    return getVFromLocalXZ(x,zeta, checkRange);
}
 
inline int CgemGeoReadoutPlane::getVIDInNextSheetFromVID(int vID, double phimin_next) const {
    int vId_next = -1;
    //vID=vID-m_NVstrips*m_iSheet;
    if( (vID>=0) && (vID<m_NVstrips) ){
        double v = getCentralVFromVID(vID);
        double v2 = getVInNextSheetFromV(v, phimin_next);
        //if(v2>=0.) return floor(v2/m_VPitch);
        vId_next = getVIDFromV(v2);// assume the next sheet with the same configuration
        //if(m_iLayer>0) { vId_next=vId_next+(1-2*m_iSheet)*m_NVstrips;// with assumption Nsheet=2 for layerID>0 FIXME}
    }
    return vId_next;
}
 
inline double CgemGeoReadoutPlane::getZFromPhiV(double phi, double V, int checkXRange) const {
    phi = phi-m_Phimin;
    while(phi<0)             phi+=CLHEP::twopi;
    while(phi>=CLHEP::twopi) phi-=CLHEP::twopi;
    return getZFromXV(phi*m_Rmid, V, checkXRange);
}
 
inline int  CgemGeoReadoutPlane::getVIDFromV(double V) const{
        double Vid = floor((V-m_Vmin_strip)/m_VPitch);
        if(Vid<0)           Vid=-1;
        if(Vid>=m_NVstrips) Vid=-2;
        //Vid=Vid+m_NVstrips*m_iSheet;
        return Vid;
}
#endif  /* CGEMGEOREADOUTPLANE_H */