CgemGeoReadoutPlane.cxx

Go to the documentation of this file.

#include "CgemGeomSvc/CgemGeoReadoutPlane.h"
 
using namespace std;
using namespace CLHEP;
 
CgemGeoReadoutPlane::CgemGeoReadoutPlane(int iLayer, int iSheet,
        double rx, double rv, 
        double phiMin, 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){
 
    m_iLayer = iLayer;
    m_iSheet = iSheet;
 
    m_RX = rx;
    m_RV = rv;
    m_Rmid = 0.5*(rx+rv);
    
    m_W = w;
    m_Zmin = zmin;
    m_Zmax = zmin+L;
    m_L = L;
    
    m_XPitch = Xpitch;            
    m_XstripWidth = XstripWidth;  
    m_VPitch = Vpitch;            
    m_VstripWidth = VstripWidth;  
    //m_XPitch = Xpitch*m_RX/m_Rmid;            // scaling from plane to cylinder
    //m_XstripWidth = XstripWidth*m_RX/m_Rmid;  // scaling from plane to cylinder
    //m_VPitch = Vpitch*m_RV/m_Rmid;            // scaling from plane to cylinder
    //m_VstripWidth = VstripWidth*m_RV/m_Rmid;  // scaling from plane to cylinder
    m_StereoAngle = stereoAngle / 180.*CLHEP::pi; //degree -> radian
    m_k = stereoAngle/fabs(stereoAngle); // sign of the stereo angle
    m_MidRGap = midROfGap;
    m_OutRGap = outROfGap;
    m_NXstrips = NXStrip;
    m_NVstrips = NVStrip;
 
    m_Phimin = phiMin;
    while(m_Phimin<-M_PI) m_Phimin+=2*M_PI;
    while(m_Phimin>M_PI)  m_Phimin-=2*M_PI;
    //m_Xmin = m_Phimin*rx;
    m_Xmin = m_Phimin*m_Rmid;
 
    //double dPhi_sensitive = NXStrip*m_XPitch/m_Rmid;
    //m_Phimin_strip = dX_strip/m_Rmid+m_Phimin;
    m_Phimin_strip = dX_strip/m_Rmid;
    //while(m_Phimin_strip<0)      m_Phimin_strip+=2*M_PI;
    //while(m_Phimin_strip>2*M_PI) m_Phimin_strip-=2*M_PI;
    m_Xmin_strip = m_Phimin_strip*m_Rmid;
    m_Xmax_strip = m_Xmin_strip+NXStrip*m_XPitch;
    m_Phimax_strip = m_Xmax_strip/m_Rmid;
    //if(NXStrip*m_XPitch>2*M_PI*rx/) cout<<"CgemGeoReadoutPlane : "<<endl;
 
    //m_Vmax = m_L*fabs(sin(m_StereoAngle))+m_W*m_RV/m_RX*cos(m_StereoAngle);
    //m_Vmax = m_L*fabs(sin(m_StereoAngle))+m_W*cos(m_StereoAngle);
    //m_NXstrips = floor(m_W / m_XPitch);
    //m_NVstrips = floor(m_Vmax / m_VPitch);
    //m_Vmin = 0.;
    
    m_Vtot = NXStrip*m_XPitch*cos(m_StereoAngle)+m_L*fabs(sin(m_StereoAngle));
    dV=(m_Vtot-NVStrip*m_VPitch)/2.0;// FIXME?
    m_Vmin_strip = dV;
    m_Vmax = m_Vmin_strip+NVStrip*m_VPitch;
 
 
    cout<<"    dX = "<<dX_strip<<endl;
    print();
}
 
//Long Peixun's update: Add destructor
CgemGeoReadoutPlane::~CgemGeoReadoutPlane()
{}
 
void CgemGeoReadoutPlane::print(){
    cout<<"construct one readout-plane with: "<<endl;
    cout<<"    iLayer "<<m_iLayer<<", iSheet "<<m_iSheet<<endl;
    cout<<"    Rx, Rv, Xmin, width = "<< m_RX<<", "<<m_RV<<", "<<m_Xmin<<", "<<m_W<<endl;
    cout<<"    Zmin, L = "<< m_Zmin <<", "<<m_L<<endl;
    cout<<"    XPitch, XstripWidth ="<<m_XPitch<<", "<<m_XstripWidth<<endl;
    cout<<"    sensitive X (Xmin as origin): "<<m_Xmin_strip<<" to "<<m_Xmax_strip<<endl;
    cout<<"    sensitive range in phi (phimin as origin): "<<m_Phimin_strip<<" to "<<m_Xmax_strip/m_Rmid<<endl;
    //cout<<"    VPitch, VstripWidth, StereoAngle, Vmax ="<<m_VPitch<<", "<<m_VstripWidth<<", "<<m_StereoAngle<<", "<<m_Vmax<<endl;
    cout<<"    VPitch, VstripWidth, StereoAngle, k ="<<m_VPitch<<", "<<m_VstripWidth<<", "<<m_StereoAngle<<", "<<m_k<<endl;
    cout<<"    Vtot, Vmin, Vmax = "<<m_Vtot<<", "<<m_Vmin_strip<<", "<<m_Vmax<<endl;
    cout<<"    N_Xstrips, N_Vstrips="<<m_NXstrips<<", "<<m_NVstrips<<endl;
    cout<<"    MidRGap="<<m_MidRGap<<", OutRGap="<<m_OutRGap<<endl;
}
 
double CgemGeoReadoutPlane::getX(double phi) const
{
    double dPhi = phi - m_Phimin;
    while(dPhi<0)              dPhi += CLHEP::twopi;
    while(dPhi>=CLHEP::twopi)  dPhi -= CLHEP::twopi;
    if(dPhi>m_Phimax_strip) {
        double dphi_start = phi-(m_Phimin+m_Phimin_strip);
        while(dphi_start> 0)        dphi_start-=2.*M_PI;
        while(dphi_start<=-2.*M_PI) dphi_start+=2.*M_PI;
        //double dphi_end   = phi-(m_Phimin+m_Phimax_strip);
        //while(dphi_end< 0)      dphi_end+=2.*M_PI;
        //while(dphi_end>=2*M_PI) dphi_end-=2.*M_PI;
        double dphi_end = dPhi-m_Phimax_strip;
        if(dphi_start+dphi_end>0) dPhi=dphi_start+m_Phimin_strip;
    }
    double X = m_Rmid*dPhi;// in X-zeta plane (X-V plane)
    return X;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int CgemGeoReadoutPlane::getXStripID(double phi) const {
    double X = getX(phi);
    //cout<<"X="<<X<<", Xmin_strip="<<m_Xmin_strip<<endl;
    X -= m_Xmin_strip;
    //cout<<"X="<<X<<endl;
    int xID = floor(X/m_XPitch);
    //cout<<"getXStripID: "<<xID<<endl;
    if(xID<0) xID=-1;
    //if(xID>=m_NXstrips) xID=m_NXstrips-1;
    if(xID>=m_NXstrips) xID=-2;
    //if(xID>=0&&xID<m_NXstrips) return xID;
    //else return -1;
    return xID;
}
 
int CgemGeoReadoutPlane::getClosestXStripID(double phi, double& dist) 
{
    //cout<<"getClosestXStripID: phi0="<<phi<<endl;
    //cout<<"getClosestXStripID: phi="<<phi<<endl;
    //cout<<"getXStripID: "<<xID<<endl;
    int xID = getXStripID(phi);
    int xID2 = xID;
    if(xID2==-1) xID2=0;
    if(xID2==-2) xID2=m_NXstrips-1;
    double X = getX(phi);
    dist = fabs(X-getCentralXFromXID(xID2))-0.5*m_XstripWidth;
    //cout<<"getClosestXStripID: m_Xmin="<<m_Xmin<<", m_RX="<<m_RX<<", X="<<X
    //  <<", getCentralXFromXID(xID)="<<getCentralXFromXID(xID)
    //  <<", m_XstripWidth="<<m_XstripWidth<<endl;
    return xID;
}
 
double CgemGeoReadoutPlane::getDist2ClosestXStripCenter(double phi, int& id) 
{
    //cout<<"getClosestXStripID: phi0="<<phi<<endl;
    //phi = phi - m_Xmin/m_RX;
    //while(phi<0) phi += CLHEP::twopi;
    //while(phi>CLHEP::twopi)  phi-=CLHEP::twopi;
    //cout<<"getClosestXStripID: phi="<<phi<<endl;
    //double X = m_RX * phi;// in X-zeta plane (X-V plane)
    //if(X == m_W) X -= 0.001 * m_XPitch;
    //int xID = floor((X-(m_Xmin_strip-m_Xmin))/m_XPitch);
    //cout<<"getXStripID: "<<xID<<endl;
    //if(xID==m_NXstrips) xID=m_NXstrips-1;
    //if(xID<0) xID=0;
    //if(xID>=m_NXstrips) xID=m_NXstrips-1;
    double dist=-9999.;
    int xID = getXStripID(phi);
    id=xID;
    int xID2 = xID;
    if(xID2==-1) xID2=0;
    if(xID2==-2) xID2=m_NXstrips-1;
    if(xID2>=0&&xID2<m_NXstrips) {
        double X = getX(phi);
        dist = X-getCentralXFromXID(xID2);
        //cout<<"getClosestXStripID: m_Xmin="<<m_Xmin<<", m_RX="<<m_RX<<", X="<<X
        //  <<", getCentralXFromXID(xID)="<<getCentralXFromXID(xID)
        //  <<", m_XstripWidth="<<m_XstripWidth<<endl;
    }
    return dist;
}
 
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int CgemGeoReadoutPlane::getXStripID(double x, double y) const {
    if(x!=0||y!=0){
        double phi=atan2(y,x);
        return getXStripID(phi);
    } else{
        return -3;
    }
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int CgemGeoReadoutPlane::getXStripID(G4ThreeVector pos) const {
    double x = pos.x();
    double y = pos.y();
    return getXStripID(x,y);
}
 
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int CgemGeoReadoutPlane::getClosestVStripID(G4ThreeVector pos, double& dist) const
{
    double z = pos.z();
    dist = -9999.;
    int V_ID=-99;
    if(z < m_Zmin){
        V_ID = -3;
    } else if(z > m_Zmin+m_L){
        V_ID = -4;
    } else {
        double phi = pos.phi(); 
        double v   = getVFromPhiZ(phi, z);
        //V_ID = floor((v-m_Vmin_strip) / m_VPitch);
        V_ID = getVIDFromV(v);
        int vid2 = V_ID;
        if(V_ID==-1) vid2=0;
        if(V_ID==-2) vid2=m_NVstrips-1;
        if(v>=0.) dist       = fabs(v-getCentralVFromVID(vid2))-0.5*m_VstripWidth;
    }
 
    return V_ID;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double CgemGeoReadoutPlane::getDist2ClosestVStripCenter(G4ThreeVector pos, int& id) {
    double z = pos.z();
    double dist=-9999.;
    int V_ID=-99;
    if(z < m_Zmin){
        V_ID = -3;
    } else if(z > m_Zmin+m_L){
        V_ID = -4;
    } else {
        double phi = pos.phi(); 
        double v   = getVFromPhiZ(phi, z);
        V_ID = getVIDFromV(v);
        int vid2 = V_ID;
        if(V_ID==-1) vid2=0;
        if(V_ID==-2) vid2=m_NVstrips-1;
        dist       = v-getCentralVFromVID(vid2);
    }
    id=V_ID;
 
    return dist;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void CgemGeoReadoutPlane::getStripID(G4ThreeVector pos, int& X_ID, int& V_ID) const {
    /*
    double z = pos.z();
    if(z < m_Zmin){
        V_ID = -2;
    } else if(z > m_Zmin+m_L){
        V_ID = -1;
    } else if(X_ID >= 0){
        double x = pos.x();
        double y = pos.y();
        double phi = atan2(y,x);
        phi = phi - m_Xmin/m_RX;
        while(phi < 0)             phi+=CLHEP::twopi;
        while(phi > CLHEP::twopi)  phi-=CLHEP::twopi;
        double X = m_RX * phi;// in X-zeta plane (X-V plane) @ m_RX
        double zeta = z - m_Zmin;
        V_ID = floor((getVFromLocalXZ(X, zeta)-m_Vmin_strip) / m_VPitch);
    } else{
        V_ID = -3;
    }
    */
    
    X_ID = getXStripID(pos);
    double dist=-9999.;
    V_ID = getClosestVStripID(pos, dist);
}
 
double CgemGeoReadoutPlane::getPhiFromXID(int X_ID) const { // from -pi to pi
    double X=getCentralXFromXID(X_ID);// central X with m_Xmin as origin
    double phi=(X+m_Xmin)/m_Rmid;
    while(phi<-M_PI) phi+=CLHEP::twopi;
    while(phi>M_PI)  phi-=CLHEP::twopi;
    return phi;
}
 
double CgemGeoReadoutPlane::getZFromXV(double X, double V, int checkXRange, int checkVRange) const {
    double localX = X-m_Xmin_strip;// with m_Xmin_strip as origin
    double z = -9999.0;
    if( (checkVRange==0||((V>=m_Vmin_strip)&&(V<m_Vmax))) && ((checkXRange==0)||((localX>=0.0)&&(localX<=m_Xmax_strip-m_Xmin_strip))) ){
        //double localXOfVstrip = V / cos(m_StereoAngle);
        //double zeta = (localXOfVstrip - localX) / tan(m_StereoAngle);
        //if( m_StereoAngle < 0 ) zeta += m_L;
        //if( (zeta>=0) && (zeta<=m_L) ) 
        //z = zeta + m_Zmin;
        double z0=m_Zmax; if(m_k<0) z0=m_Zmin;
        z = localX/tan(m_StereoAngle) - V/sin(m_StereoAngle) + z0;
    }
    return z;
}
 
void CgemGeoReadoutPlane::getFiredStripID(G4ThreeVector pos1, G4ThreeVector pos2,
        vector<int>& vecXID, vector<int>& vecVID) const {
    //cout << "get fired strip id" << endl;
    bool printDebug = false;
    double z[2];
    z[0] = pos1.z();
    z[1] = pos2.z();
    if( (z[0]<m_Zmin && z[1]<m_Zmin) || (z[0]>(m_Zmin+m_L) && z[1]>(m_Zmin+m_L)) ) return;
 
    double phi[2];
    double phimin=-CLHEP::pi;
    phi[0]=pos1.phi();
    phi[1]=pos2.phi();
    if(fabs(phi[1]-phi[0])>CLHEP::pi){
        phimin=0.0;
        for(int i=0; i<2; i++) if(phi[i]<0.0) phi[i]+=CLHEP::twopi;
    }
    if(phi[0]>phi[1]){
        double tmp=phi[0];
        phi[0]=phi[1];
        phi[1]=tmp;
        tmp=z[1];
        z[1]=z[0];
        z[0]=tmp;
    }
    if(printDebug) cout<<"phimin="<<phimin<<", phi[0,1]="<<phi[0]<<","<<phi[1]<<", z[0,1]="<<z[0]<<","<<z[1]<<endl;
    for(int i=0; i<2; i++){
        if(z[i]<m_Zmin){
            double phiTmp=(phi[1]-phi[0])/(z[1]-z[0])*(m_Zmin-z[0])+phi[0];
            phi[i]=phiTmp;
            z[i]=m_Zmin;
        }
        if(z[i]>(m_Zmin+m_L)){
            double phiTmp=(phi[1]-phi[0])/(z[1]-z[0])*(m_Zmin+m_L-z[0])+phi[0];
            phi[i]=phiTmp;
            z[i]=m_Zmin+m_L;
        }   
    }
    double phiActiveMin=m_Phimin+m_Phimin_strip+10e-10;// 2019-05-05 wangll, +10e-10 2019-10-30
    double phiActiveMax=m_Phimin+m_Phimax_strip-10e-10;// 2019-05-05 wangll, -10e-10 2019-10-30
    if(printDebug) cout<<"phimin="<<phimin<<", phiActiveMin="<<phiActiveMin<<", phiActiveMax="<<phiActiveMax<<", phi[0,1]="<<phi[0]<<","<<phi[1]<<", z[0,1]="<<z[0]<<","<<z[1]<<endl;
    while(phiActiveMax>=CLHEP::twopi) phiActiveMax-=CLHEP::twopi;// [0,2pi]
    while(phiActiveMin>=CLHEP::twopi) phiActiveMin-=CLHEP::twopi;// [0,2pi]
    while(phiActiveMax<0            ) phiActiveMax+=CLHEP::twopi;// [0,2pi]
    while(phiActiveMin<0            ) phiActiveMin+=CLHEP::twopi;// [0,2pi]
    if(phimin==-CLHEP::pi) 
    {
        /*if(phiActiveMin>CLHEP::pi||phiActiveMax>CLHEP::pi) {
          phiActiveMin-=CLHEP::twopi;
          phiActiveMax-=CLHEP::twopi;
          }*/
        while(phiActiveMin>CLHEP::pi) phiActiveMin-=CLHEP::twopi;
        while(phiActiveMax>CLHEP::pi) phiActiveMax-=CLHEP::twopi;
    }// [-pi,pi]
    if(printDebug) cout<<"phimin="<<phimin<<", phiActiveMin="<<phiActiveMin<<", phiActiveMax="<<phiActiveMax<<", phi[0,1]="<<phi[0]<<","<<phi[1]<<", z[0,1]="<<z[0]<<","<<z[1]<<endl;
 
    int iCase=-1;
    if(phiActiveMin<=phi[0])
    {
        if(phiActiveMax<=phiActiveMin) iCase=1;
        else if(phiActiveMax<=phi[0])  iCase=2;
        else if(phiActiveMax<=phi[1])  iCase=3;
        else                           iCase=4;
    }
    else if(phiActiveMin<=phi[1])
    {
        if(phiActiveMax<=phi[0])             iCase=5;
        else if(phiActiveMax<=phiActiveMin)  iCase=6;
        else if(phiActiveMax<=phi[1])        iCase=7;
        else                                 iCase=8;
    }
    else
    {
        if(phiActiveMax<=phi[0])             iCase=9;
        else if(phiActiveMax<=phi[1])        iCase=10;
        else if(phiActiveMax<=phiActiveMin)  iCase=11;
        else                                 iCase=12;
    }
    if(printDebug) cout<<"iCase = "<<iCase<<endl;
 
    double zBoundary[2];
    zBoundary[0]=z[0]+(z[1]-z[0])/(phi[1]-phi[0])*(phiActiveMin-phi[0]);
    zBoundary[1]=z[0]+(z[1]-z[0])/(phi[1]-phi[0])*(phiActiveMax-phi[0]);
    switch(iCase)
    {
        case 2:
        case 9:
        case 12:
            return;// no cluster
        case 1:
        case 4:
        case 6:
        case 11:
            break;// no change in points
        case 3:
        case 10:
            phi[1]=phiActiveMax;
            z[1]=zBoundary[1];
            break;
        case 5:
        case 8:
            phi[0]=phiActiveMin;
            z[0]=zBoundary[0];
            break;
        case 7:
            phi[0]=phiActiveMin;
            phi[1]=phiActiveMax;
            z[0]=zBoundary[0];
            z[1]=zBoundary[1];
            break;
        default:
            cout<<"Error in CgemGeoReadoutPlane::getFiredStripID: iCase = "<<iCase<<endl;
            return;
    }
    if(printDebug) cout<<"phimin="<<phimin<<", phiActiveMin="<<phiActiveMin<<", phiActiveMax="<<phiActiveMax<<", phi[0,1]="<<phi[0]<<","<<phi[1]<<", z[0,1]="<<z[0]<<","<<z[1]<<endl;
 
    G4ThreeVector pos_1,pos_2;
    pos_1.setRhoPhiZ(m_Rmid,phi[0],z[0]);
    pos_2.setRhoPhiZ(m_Rmid,phi[1],z[1]);
    int XID_min(-9),XID_max(-9);
    int VID_min(-9),VID_max(-9);
    getStripID(pos_1,XID_min,VID_min);
    getStripID(pos_2,XID_max,VID_max);
    if( iCase!=6 && VID_max<VID_min ) {
        int tmp=VID_max;
        VID_max=VID_min;
        VID_min=tmp;
    }
    if(printDebug) {
        cout<<"on the sheet: phi "<<phi[0]<<" to "<<phi[1]<<", z "<<z[0]<<" to "<<z[1]<<endl;
        cout<<"              XID "<<XID_min<<" to "<<XID_max<<", VID "<<VID_min<<" to "<<VID_max<<endl;
    }
 
    if(iCase!=6)
    {
        for(int i=0; i<=abs(XID_min-XID_max);i++) vecXID.push_back(XID_min+i);
        for(int i=0; i<=abs(VID_min-VID_max);i++) vecVID.push_back(VID_min+i);
    }
    else
    {
        pos_1.setRhoPhiZ(m_Rmid,phiActiveMax,zBoundary[1]);
        pos_2.setRhoPhiZ(m_Rmid,phiActiveMin,zBoundary[0]);
        int XID_1(-9),XID_2(-9);
        int VID_1(-9),VID_2(-9);
        getStripID(pos_1,XID_1,VID_1);
        getStripID(pos_2,XID_2,VID_2);
        for(int i=0; i<=abs(XID_min-XID_1);i++) vecXID.push_back(XID_min+i);
        for(int i=0; i<=abs(XID_max-XID_2);i++) vecXID.push_back(XID_2+i);
        int VID_0=min(VID_1,VID_min);
        for(int i=0; i<=abs(VID_min-VID_1);i++) vecVID.push_back(VID_0+i);
        VID_0=min(VID_2,VID_max);
        for(int i=0; i<=abs(VID_max-VID_2);i++) vecVID.push_back(VID_0+i);
    }
}
 
double CgemGeoReadoutPlane::incidentAngleX(Hep3Vector momentum, double phi) const
{
    double phiTrk = momentum.phi();
    /*double dPhi = phiTrk-phi;
    while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
    while(dPhi> CLHEP::pi) dPhi-=CLHEP::twopi;
    if(dPhi<-0.5*CLHEP::pi) dPhi+=CLHEP::pi;
    else if(dPhi>0.5*CLHEP::pi) dPhi-=CLHEP::pi;
    */
    double dPhi = dAngleAcute(phiTrk,phi);
    return dPhi;
}
 
double CgemGeoReadoutPlane::incidentAngleV(Hep3Vector momentum, double phi) const
{
    Hep3Vector zaxis(0,0,1);
    momentum.rotate(zaxis, -phi);
    Hep3Vector xaxis(1,0,0);
    momentum.rotate(xaxis, m_StereoAngle);
    double phiTrk = momentum.phi();
    double dPhi = dAngleAcute(phiTrk,0);
    return dPhi;
}
 
double CgemGeoReadoutPlane::dAngleAcute(double phi1, double phi2) const
{
    double dPhi = phi1-phi2;
    while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
    while(dPhi> CLHEP::pi) dPhi-=CLHEP::twopi;
    if(dPhi<-0.5*CLHEP::pi) dPhi+=CLHEP::pi;
    else if(dPhi>0.5*CLHEP::pi) dPhi-=CLHEP::pi;
    return dPhi;
}