GeoStripNode.cxx

Go to the documentation of this file.

// $Id: GeoStripNode.cxx,v 1.17 2007/11/11 07:28:43 rhatcher Exp $
//
// GeoStripNode
//
// GeoStripNode is a strip node
//
// Author:  S. Kasahara 06/04

#include <TGeoVolume.h>
#include <TGeoMatrix.h>
#include <TGeoBBox.h>
#include <TMath.h>

#include "MessageService/MsgService.h"
#include "GeoGeometry/GeoGeometry.h"
#include "GeoGeometry/GeoStripNode.h"
#include "GeoGeometry/GeoScintPlnNode.h"
#include "GeoGeometry/GeoScintMdlNode.h"

ClassImp(GeoStripNode)
CVSID("$Id: GeoStripNode.cxx,v 1.17 2007/11/11 07:28:43 rhatcher Exp $");

//_____________________________________________________________________________
GeoStripNode::GeoStripNode(GeoGeometry* geom, TGeoVolume* stpvol, 
                           TGeoMatrix* stpmatrix, GeoScintMdlNode* mdlNode, 
                           std::string nodename,const PlexStripEndId& stripid)
  : GeoNode(geom,stpvol,stpmatrix,
            mdlNode->GetScintMdlVolume()->GetAirNode()->GetVolume(),
            mdlNode->GetGlobalPath()+"/" + 
            mdlNode->GetScintMdlVolume()->GetAirNode()->GetName() + "/" +
            nodename,nodename),fStripEndId(stripid),fScintMdlNode(mdlNode) {
  // Normal constructor
  // Node is constructed and added to volume parentvol with position specified
  // by stpmatrix.
  // Protected

  UpdateGlobalManager();

}

//_____________________________________________________________________________
GeoStripNode::~GeoStripNode() {
  // Destructor should delete any owned objects

  UpdateGlobalManager();

  if ( CountRef() ) {
    MSG("Geo",Msg::kWarning)
      << "GeoStripNode destructor " << GetSEId()
      << " still had " << CountRef()
      << " outstanding references " << endl;
  }
}

//_____________________________________________________________________________
Float_t GeoStripNode::GetHalfLength() const {
  // Return the strip's half length (along local x)

  UpdateGlobalManager();
  Float_t halflength = GetStripVolume() -> GetHalfLength();
  
  MSG("Geo",Msg::kDebug) << "GetHalfLength " << fStripEndId.AsString() 
                         << " halflength " << halflength << endl;
  
  return halflength;

}

//_____________________________________________________________________________
Float_t GeoStripNode::GetHalfThickness() const {
  // Return the strip's half thickness (along local z)

  UpdateGlobalManager();

  Float_t halfthickness = GetStripVolume() -> GetHalfThickness();
  MSG("Geo",Msg::kDebug) << "GetHalfThickness " << fStripEndId.AsString() 
                         << " halfthickness " << halfthickness << endl;
  
  return halfthickness;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::GetHalfWidth() const {
  // Return the strip's half width (along local y)

  UpdateGlobalManager();

  Float_t halfwidth = GetStripVolume() -> GetHalfWidth();
  MSG("Geo",Msg::kDebug) << "GetHalfWidth " << fStripEndId.AsString() 
                         << " halfwidth " << halfwidth << endl;
  
  return halfwidth;

}

//_____________________________________________________________________________
Float_t GeoStripNode::PartialLength(StripEnd::StripEnd_t stripend) const{
  // Return the strip's partial length at the specified end
  // If the strip end is not specified, return the total length

  UpdateGlobalManager();

  Float_t lenpart = GetStripVolume() -> GetLenPart(stripend);
  MSG("Geo",Msg::kDebug) << "PartialLength " << fStripEndId.AsString() 
                         << " end " << StripEnd::AsString(stripend) 
                         << " lenpart " << lenpart << endl;
  
  return lenpart;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::WlsBypass() const {
  // Return the length of the wls fiber along coil hole bypass

  UpdateGlobalManager();

  Float_t wlsbypass = GetStripVolume() -> GetWlsLenBypass();

  MSG("Geo",Msg::kDebug) << "WlsBypass " << fStripEndId.AsString() 
                         << " wlsbypass " << wlsbypass << endl;

  return wlsbypass;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::WlsPigtail(StripEnd::StripEnd_t stripend) const {
  // Return the length of the wls fiber "pigtail" (the portion extending
  // beyond the strip end), if either StripEnd::kEast or StripEnd::kWest
  // is specified.  Otherwise, return 0 with error message. pigtail length
  // includes the extra pigtail length associated with the scint module.

  UpdateGlobalManager();

  Float_t wlspigtail = 0;
  
  switch (stripend) {

  case StripEnd::kEast:
  case StripEnd::kWest: 
    {
      GeoScintMdlNode* mdlNode = GetScintMdlNode();
      wlspigtail = GetStripVolume()->GetWlsLen(stripend) 
        + mdlNode->GetExtraWlsFiber(stripend);
    }
    break;
  
  default:
    MSG("Geo",Msg::kWarning) << "WlsPigtail called with stripend "
                             << StripEnd::AsString(stripend) 
                             << ". Must specify kEast or kWest." << endl;
  } // end of switch

  MSG("Geo",Msg::kDebug) << "WlsPigtail " << fStripEndId.AsString() << " end " 
                         << StripEnd::AsString(stripend) 
                         << " wlspigtail " << wlspigtail << endl;

  return wlspigtail;
  
}

//_____________________________________________________________________________
Bool_t GeoStripNode::IsMirrored(StripEnd::StripEnd_t stripend) const {
  // Return true if specified stripend is mirrored, else false. 

  UpdateGlobalManager();

  Bool_t isMirrored = GetStripVolume() -> IsMirrored(stripend);

  MSG("Geo",Msg::kDebug) << "IsMirrored " << fStripEndId.AsString() << " end " 
                         << StripEnd::AsString(stripend) << " isMirrored "
                         << isMirrored << endl;

  return isMirrored;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::GetLPosRelMdl() const {
  // Return the lpos relative to module (x-pos from translation matrix)


  UpdateGlobalManager();

  const Double_t* trans = GetMatrix() -> GetTranslation();
  
  MSG("Geo",Msg::kDebug) << "GetLPosRelMdl " << fStripEndId.AsString() 
                         << " lposrelmdl " << trans[0] << endl;

  return trans[0]; 
  
}

//_____________________________________________________________________________
GeoScintPlnNode* GeoStripNode::GetScintPlnNode(void) const{
  // Return scint plane node from which this strip descends

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GetScintPlnNode " << fStripEndId.AsString() 
                         << endl;

  return (GetScintMdlNode()->GetScintPlnNode());
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::GetTPosRelMdl() const {
  // Return the tpos relative to module (y-pos from translation matrix)

  MSG("Geo",Msg::kDebug) << fStripEndId.AsString() << endl;

  UpdateGlobalManager();

  const Double_t* trans = GetMatrix() -> GetTranslation();

  return trans[1];
  
}


//_____________________________________________________________________________
Float_t GeoStripNode::GetZRotRelMdlRad() const {
  // Return the zrot relative to module in radians

  MSG("Geo",Msg::kDebug) << fStripEndId.AsString() << endl;

  UpdateGlobalManager();

  const Double_t* rotmatrix = GetMatrix() -> GetRotationMatrix();
  Float_t zrotrelmdlrad = TMath::ATan2(rotmatrix[3],rotmatrix[0]);

  return zrotrelmdlrad;
   
}

//_____________________________________________________________________________
Float_t GeoStripNode::GetTPos(Float_t orthCoord) const {
  // Get the strip's tpos offset from scint plane origin
  // orthCoord is the t-pos in the orthogonal view in ideal scint plane
  // coordinates. For example, if this strip is a U-view strip, the orthcoord
  // is the tpos in the "ideal" V-view plane coordinates.
  // This allows one to include rotation corrections.
  // if orthCoord > 100000 (cm) (default), use the center of the strip.

  MSG("Geo",Msg::kDebug) << fStripEndId.AsString() << " orthCoord "
                         << orthCoord << endl;

  UpdateGlobalManager();

  // First convert strip center to scint plane coordinates
  Double_t striplxyz[3] = {0}; // center of strip
  // Convert strip center in strip coordinates to module coordinates
  Double_t modulelxyz[3] = {0};
  GetMatrix() -> LocalToMaster(striplxyz,modulelxyz);
  // Convert strip center in module coordinates to scint plane coordinates
  Double_t scplanelxyz[3] = {0};      
  GeoScintMdlNode* mdlNode = GetScintMdlNode();
  mdlNode -> GetMatrix() -> LocalToMaster(modulelxyz,scplanelxyz);

  // tpos of center of strip
  Float_t tpos = scplanelxyz[1];
  if ( orthCoord > 100000. ) return tpos;
  
  // Adjust tpos at orthCoord position for rotation of strip about z-axis
  Float_t zrotrelmdlrad = GetZRotRelMdlRad();
  Float_t zrotrelplnrad = mdlNode->GetZRotRelPlnRad();
  GeoScintPlnNode* scplnNode = GetScintPlnNode();
  Float_t zrotrelsteelrad = scplnNode->GetZRotRelSteelRad();

  // Pick the ideal steel plane rotation and subtract it off
  Float_t zrotrelsteelideal = 0;
  Float_t lpos = orthCoord;
  if ( fStripEndId.GetDetector() == Detector::kFar ||
       fStripEndId.GetDetector() == Detector::kNear  ) {
    if ( fStripEndId.GetPlaneView() == PlaneView::kU ) {
      zrotrelsteelideal = -0.785398163; // -45 degrees
      lpos = -orthCoord; // lpos is in opposite direction of V-view tpos
    }
    else if ( fStripEndId.GetPlaneView() == PlaneView::kV ) {
      zrotrelsteelideal = +0.785398163; // +45 degrees
    }
  }
  else if ( fStripEndId.GetDetector() == Detector::kCalDet ) {
    if ( fStripEndId.GetPlaneView() == PlaneView::kU ) {
      zrotrelsteelideal = -1.570796327;
      lpos = -orthCoord; // lpos is in opposite direction of V-view tpos
    }
    else if ( fStripEndId.GetPlaneView() == PlaneView::kV ) {
      zrotrelsteelideal = 0;
    }
  }
  
  Float_t zrottot = (zrotrelsteelrad - zrotrelsteelideal) + zrotrelplnrad
                  +  zrotrelmdlrad;
  
  Float_t deltatpos = lpos * TMath::Tan(zrottot);
  tpos += deltatpos;
  
  return tpos;
  
}



//_____________________________________________________________________________
Float_t GeoStripNode::ClearFiber(StripEnd::StripEnd_t end) const {
  // Return clear fiber length of scint module to which this strip belongs

  UpdateGlobalManager();

  Float_t clearfiber = GetScintMdlNode() -> GetClearFiber(end);
  
  MSG("Geo",Msg::kDebug) << "ClearFiber " << fStripEndId.AsString() << " end "
                         << StripEnd::AsString(end) 
                         << " clearfiber " << clearfiber << endl;

  return clearfiber;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::TotalAttenuation(StripEnd::StripEnd_t /* end */,
                                       const Float_t /* alongLength */) const {
  // Dummy, but also dummy in UgliStripNode.

  UpdateGlobalManager();

  MSG("Geo",Msg::kWarning) 
    << "GeoStripNode::TotalAttenuation not implemented yet." << endl;
  return 0;
  
}

//_____________________________________________________________________________
Float_t GeoStripNode::DistanceAlong(const PlexStripEndId& orthogonalStrip) 
 const{
  // Get the distance along a strip represented by the specified orthogonal 
  // strip. Note this is relative to the strip's local coordinate system 
  // origin.
  //
  // This strip node is expressed as a line in parametric form:
  //     xyz(i) = base(i) + s*dir(i) for i= 0,1,2 (x,y,z)
  // where base is LocalToGlobal(TVector3(0,0,0)) 
  // and dir is (LocalToGlobal(TVector3(1,0,0)-base)), i.e. the 
  // unit vector along the strip length.
  // The returned value is "s", the distance along the strip length 
  // (measured from the strip local origin) defining the point on 
  // this strip with the smallest distance to any point on the orthogonal 
  // strip.
  //
  // The algorithm used here to determine the intersection or closest approach
  // of two skew lines is the same as that in R. Hatcher's UgliStripNode, 
  // which was adapted from:
  //   Graphic Gems, ed. Andrew S. Glassner, pg 304
  //   ISBN 0-12-286165-5   T385.G697 (1990)
  //
  
  UpdateGlobalManager();

  const Double_t farAway = 1.0e25;
  
  GeoStripNode* otherNode = GetGeoGeometry()->GetStripNode(orthogonalStrip);
  
  if ( !otherNode ) {
    MSG("Geo",Msg::kError)
      << "GeoStripNode::DistanceAlong(), other strip "
      << orthogonalStrip << " does not exist." << endl;
    return farAway;
  }
  
  // Define base & dir for this strip
  TVector3 baseThis = this -> GlobalPos(0.,false);
  TVector3 dirThis = this -> GlobalPos(1.,false) - baseThis; 
  
  // Define base & dir for the other strip
  TVector3 baseOther = otherNode -> GlobalPos(0.,false);
  TVector3 dirOther = otherNode -> GlobalPos(1.,false) - baseOther;
  
  TVector3 dirCross = dirThis.Cross(dirOther);
  Double_t dirCrossMag2 = dirCross.Mag2();
  if ( dirCrossMag2 == 0.0 ) {
    // exactly parallel
    return farAway;
  }

  TVector3 dp = baseOther - baseThis;
  TVector3 dpxv = dp.Cross(dirOther);
  Double_t s = dpxv.Dot(dirCross)/dirCrossMag2;
  
  return s;
  
}

//_____________________________________________________________________________
TVector3 GeoStripNode::Intersection(const PlexStripEndId& orthogonalStrip) 
const {
  // Return the 3-d "intersection" of two strips in global coordinates.  
  // The intersection is determined as the mid-point of the vector joining 
  // the strips at the distance of closest approach.  
  //
  // The strip nodes are expressed as a line in parametric form:
  //     xyz(i) = base(i) + s*dir(i) for i= 0,1,2 (x,y,z) // for this strip
  //     xyz(i) = base(i) + t*dir(i) for i= 0,1,2 (x,y,z) // for ortho strip
  // where base is LocalToGlobal(TVector3(0,0,0)) 
  // and dir is (LocalToGlobal(TVector3(1,0,0)-base)), i.e. the 
  // unit vector along the strip length.
  // The position along both strips at the distance of closest approach
  // is solved for to determine:
  //     xyz_s0(i) = base(i) + s0*dir(i) for i= 0,1,2 (x,y,z) //for this strip
  //     xyz_t0(i) = base(i) + t0*dir(i) for i= 0,1,2 (x,y,z) //for ortho strip
  // The returned mid-point is then:
  //     xyz_mid = 0.5*(xyz_s0+xyz_t0)
  //
  // The algorithm used here to determine the intersection or closest approach
  // of two skew lines is the same as that in R. Hatcher's UgliStripNode, 
  // which was adapted from:
  //   Graphic Gems, ed. Andrew S. Glassner, pg 304
  //   ISBN 0-12-286165-5   T385.G697 (1990)
  //

  UpdateGlobalManager();

  const Double_t farAway = 1.0e25;
  const TVector3 farAway3(farAway,farAway,farAway);

  GeoStripNode* otherNode = GetGeoGeometry()->GetStripNode(orthogonalStrip);
  
  if ( !otherNode ) {
    MSG("Geo",Msg::kError)
      << "GeoStripNode::DistanceAlong() other strip "
      << orthogonalStrip << " does not exist." << endl;
    return farAway3;
  }

  // Define base & dir for this strip
  TVector3 baseThis = this -> GlobalPos(0.,false);
  TVector3 dirThis = this -> GlobalPos(1.,false) - baseThis; 
  
  // Define base & dir for the other strip
  TVector3 baseOther = otherNode -> GlobalPos(0.,false);
  TVector3 dirOther = otherNode -> GlobalPos(1.,false) - baseOther;
  
  TVector3 dirCross = dirThis.Cross(dirOther);
  Double_t dirCrossMag2 = dirCross.Mag2();
  if ( dirCrossMag2 == 0.0 ) {
    // exactly parallel
    return 0.5*(baseThis+baseOther) + farAway*dirThis;
  }

  TVector3 dp = baseOther - baseThis;
  TVector3 dpxvThis = dp.Cross(dirOther);
  Double_t s0 = dpxvThis.Dot(dirCross)/dirCrossMag2;

  TVector3 xyz_s0 = baseThis + s0*dirThis;
  
  dirCross *= -1;
  dp       *= -1;
  TVector3 dpxvOther = dp.Cross(dirThis);
  Double_t t0 = dpxvOther.Dot(dirCross)/dirCrossMag2;
  
  TVector3 xyz_t0 = baseOther + t0*dirOther;
  
  return 0.5*(xyz_s0 + xyz_t0);

}

//_____________________________________________________________________________
TVector3 GeoStripNode::GlobalPos(const Float_t alongLength,
                                 const Bool_t onWLS) const {
  // Return position based on distance along centerline *or* WLS fiber
  // Position on WLS fiber is not implemented yet, as is also true for
  // UgliStripNode.

  MSG("Geo",Msg::kDebug) << fStripEndId.AsString() 
                         << " alongLength " << alongLength << endl;

  UpdateGlobalManager();

  if ( onWLS ) {
    MSG("Geo",Msg::kWarning) 
      << "GeoStripNode::GlobalPos not implemented yet with onWLS = kTRUE."
      << " Ignored." << endl;
  }
  
  return LocalToGlobal(TVector3(alongLength,0.,0.));
  
}

//_____________________________________________________________________________
void GeoStripNode::SetZRotRelMdlRad(Float_t /*radians*/) {
  // Dummy Set method

  UpdateGlobalManager();

  MSG("Geo",Msg::kWarning) 
    << "GeoStripNode::SetZRotRelMdlRad not implemented yet." << endl;
  return;
  
}

//_____________________________________________________________________________
void GeoStripNode::SetLTPosRelMdl(Float_t /*lpos*/, Float_t /*tpos*/) {
  // Dummy Set method

  UpdateGlobalManager();

  MSG("Geo",Msg::kWarning) 
    << "GeoStripNode::SetLTPosRelMdl not implemented yet." << endl;
  return;
  
}

---