OpParamAction.cxx

Go to the documentation of this file.

// OpParamAction class implementation : Ben Jones, MIT 2013
// See header file for detailed description.

#include "larsim/LegacyLArG4/OpParamAction.h"

#include "cetlib_except/exception.h"

#include <cmath>

namespace larg4
{

  //-----------------------------------------------
  //  OpParamAction base class methods
  //-----------------------------------------------



  OpParamAction::OpParamAction()
  {
  }

  OpParamAction::~OpParamAction()
  {
  }

  double OpParamAction::GetAttenuationFraction(G4ThreeVector /*PhotonDirection*/, G4ThreeVector /*PhotonPosition*/)
  {
    return 0;
  }



  //-----------------------------------------------
  //  SimpleWireplaneAction Methods
  //-----------------------------------------------


  SimpleWireplaneAction::SimpleWireplaneAction(TVector3 WireDirection, TVector3 PlaneNormal, double WirePitch, double WireDiameter)
  {
    fWireDirection = G4ThreeVector(WireDirection.X(),
                                   WireDirection.Y(),
                                   WireDirection.Z()).unit();
    fPlaneNormal   = G4ThreeVector(PlaneNormal.X(),
                                   PlaneNormal.Y(),
                                   PlaneNormal.Z()).unit();
    fDPRatio       = WireDiameter/WirePitch;
  }

  //-----------------------------------------------

  SimpleWireplaneAction::~SimpleWireplaneAction()
  {
  }

  //-----------------------------------------------
  // An ideal simple wireplane attenuates the light by a fraction
  //  Wire_Diameter / (Pitch * cos theta) where theta is the angle
  //  of incident light projected into the plane perpendicular to the
  //  wires.  The photon position is not used.
  //
  double SimpleWireplaneAction::GetAttenuationFraction(G4ThreeVector PhotonDirection, G4ThreeVector /*PhotonPosition*/)
  {
    G4ThreeVector ProjDirection = PhotonDirection - fWireDirection*(fWireDirection.dot(PhotonDirection));
    double CosTheta = std::abs(fPlaneNormal.dot(ProjDirection));
    if(CosTheta < fDPRatio)
      return 0;
    else
      return (1.0 - fDPRatio /  CosTheta);
   }



  //-----------------------------------------------
  //  OverlaidWireplanesAction Methods
  //-----------------------------------------------


  OverlaidWireplanesAction::OverlaidWireplanesAction(std::vector<std::vector<double> > InputVectors, int Orientation)
  {

    G4ThreeVector WireBasis1, WireBasis2;

    if(Orientation==0)
      {
        fPlaneNormal=G4ThreeVector(1,0,0);
        WireBasis1 =G4ThreeVector(0,1,0);
        WireBasis2 =G4ThreeVector(0,0,1);
      }
    else if(Orientation==1)
      {
        fPlaneNormal=G4ThreeVector(0,1,0);
        WireBasis1 =G4ThreeVector(1,0,0);
        WireBasis2 =G4ThreeVector(0,0,1);
      }
    else if(Orientation==2)
      {
        fPlaneNormal=G4ThreeVector(0,0,1);
        WireBasis1 =G4ThreeVector(0,1,0);
        WireBasis2 =G4ThreeVector(0,0,1);
      }
    else
      {
        throw cet::exception("OpParamAction") << "Unrecognized wireplane orientation. Options are 1=Xdrift, 2=Ydrift, 3=Zdrift\n";
      }
    for(size_t i=0; i!=InputVectors.size(); ++i)
      {
        if(InputVectors.at(i).size()!=3)
          {
            throw cet::exception("OpParamAction") << "Unrecognized wireplane parameter format. Expected vector(3)'s with v[0] = wire angle, v[1] = wire pitch, v[2] = wire diameter\n";
          }
        else
          {
            double theta = InputVectors[i][0]*3.142/180.;
            fWireDirections.push_back(cos(theta)*WireBasis1 + sin(theta)*WireBasis2);
            fDPRatios.push_back(InputVectors[i][2]/InputVectors[i][1]);
          }
      }
  }


  //-----------------------------------------------

  OverlaidWireplanesAction::~OverlaidWireplanesAction()
  {
  }


  //-----------------------------------------------

  double OverlaidWireplanesAction::GetAttenuationFraction(G4ThreeVector PhotonDirection, G4ThreeVector /*PhotonPosition*/)
  {

    double AttenFraction=1.;

    for(size_t i=0; i!=fWireDirections.size(); ++i)
      {
        G4ThreeVector ProjDirection = PhotonDirection - fWireDirections.at(i) * (fWireDirections.at(i).dot(PhotonDirection.unit()));

          // fWireDirections.at(i).cross(PhotonDirection.cross(fWireDirections.at(i)).unit());
        double CosTheta =  (ProjDirection.mag() > 0 ? std::abs(fPlaneNormal.dot(ProjDirection))/ProjDirection.mag() : 1.0);
        if(CosTheta < fDPRatios.at(i))
          return 0;
        else
          AttenFraction *= (1.0 - fDPRatios.at(i) / CosTheta);
      }
    return AttenFraction;
  }





}