infid.cxx

Go to the documentation of this file.

#include "DataUtil/infid.h"

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <iostream>
#include <cassert>

#include "RVersion.h"
#include "TInterpreter.h"
#include "TMath.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TString.h"
#endif

namespace FidVol {

  // component values used

  static std::string gName = "Default";

  // Near choices
  static bool   gNearFollowBeam = true;
  static double gNearR          = 1.0; 
  // z cuts were:                   1.0,     5.0
  // z choices between scint&steel  N017PVp, N084PVp
  static double gNearZData[2] = { 1.01080, 4.99059 };
  static double gNearZMC[2]   = { 1.01080, 4.99059 };

  // direction along "beam"
  static double gBeamAngleRad   = 3.34321 * TMath::DegToRad();
  static double gNearDyDz       = TMath::Tan(-gBeamAngleRad);
  static double gNearX0Beam     = 1.4828;
  static double gNearY0Beam     = 0.2384;

  // center when along "z"
  static double gNearX0Z     = 1.4885;
  static double gNearY0Z     = 0.1397;

  // Far choices
  static bool   gFarOctagon     = false;
  static bool   gFarCoilCut     = true;
  static double gFarRinner      = 0.5;
  static double gFarRouter      = TMath::Sqrt(14.0);
  // z cuts were:                  0.5,     14.3,     16.2,     28.0
  // z choices between scint&steel F008PUt, F240PUt,  F255PUt,  F452PVt
  static double gFarZData[4] = { 0.49080, 14.29300, 16.27110, 27.98270};
  static double gFarZMC[4]   = { 0.47692, 14.27860, 16.26470, 27.97240};

  //Variables to correct the position of the vtx
  //In R1.24 the trk vtx was in the downstream scintillator
  //whereas the neutrino interaction most often occurred in the
  //upstream steel plane
  //These numbers are to be subtracted from the evt/trk/shw vertex
  static double gEvtVtxZOffset=0.0;
  static double gTrkVtxZOffset=0.0392;
  static double gShwVtxZOffset=0.0;

  // constants
  //const double r_sqrt2 = 1.0 / TMath::Sqrt(2.0); // ACLiC doesn't like this
  const double r_sqrt2 = 7.07106781186547462e-01;

  // sub-functions
  bool   infid_near_z(SimFlag::SimFlag_t simflg, double z);
  bool   infid_far_z(SimFlag::SimFlag_t simflg, double z);
  bool   infid_near_circle_z(double x, double y);
  bool   infid_near_circle_beam(double x, double y, double z);
  bool   infid_far_coil(double x, double y);
  bool   infid_far_circle(double x, double y);
  bool   infid_far_octagon(double x, double y);

  bool load_setter_script(std::string name);
  bool know_setter(std::string funcname);

  bool legal_indx(int indx, int size) { return (indx>=0 && indx<size); }

}


//
// Actual definitions
//

bool infid(Detector::Detector_t det, SimFlag::SimFlag_t simflg, 
                  double x, double y, double z)
{
  if      ( Detector::kNear == det ) {
    // test z
    if ( ! FidVol::infid_near_z(simflg,z) ) return false;
    // test transverse
    if ( FidVol::gNearFollowBeam ) return FidVol::infid_near_circle_beam(x,y,z);
    else                           return FidVol::infid_near_circle_z(x,y);
  }
  else if ( Detector::kFar == det ) {
    // test z
    if ( ! FidVol::infid_far_z(simflg,z) ) return false;
    // test coil cut (if requested)
    if ( FidVol::gFarCoilCut && ! FidVol::infid_far_coil(x,y) ) return false;
    // test transverse
    if ( FidVol::gFarOctagon ) return FidVol::infid_far_octagon(x,y);
    else                       return FidVol::infid_far_circle(x,y);
  }
  // unknown detector
  return true;
}

//

void print_infid()
{
  std::cout << "Current infid settings: \"" << FidVol::gName 
            << "\"" << std::endl;
  std::cout << "Near:";
  std::cout << " cylinder radius " << FidVol::gNearR;
  if ( FidVol::gNearFollowBeam ) {
    std::cout << ", x0=" << FidVol::gNearX0Beam 
              << ", y0=" << FidVol::gNearY0Beam
              << ", follow beam" << std::endl;
    std::cout << "      dy/dz " << FidVol::gNearDyDz
              << " (angle " << FidVol::gBeamAngleRad << " radians)" 
              << std::endl;
  } else {
    std::cout << ", x0=" << FidVol::gNearX0Z 
              << ", y0=" << FidVol::gNearY0Z 
              << ", along z axis" << std::endl;
  }
  std::cout << "      z limits data: " 
            << Form("%8.5f",FidVol::gNearZData[0]) << " " 
            << Form("%8.5f",FidVol::gNearZData[1]) << std::endl
            << "                 MC: " 
            << Form("%8.5f",FidVol::gNearZMC[0]) << " " 
            << Form("%8.5f",FidVol::gNearZMC[1]) << std::endl;

  std::cout << "Far:";
  std::cout << "  " << (FidVol::gFarOctagon?"octagon inscribed":"circle") 
            << " radius " << FidVol::gFarRouter << ", ";
  if ( FidVol::gFarCoilCut ) 
    std::cout << "coil cut " << FidVol::gFarRinner;
  else
    std::cout << "no coil cut";     
  std::cout << std::endl;
  std::cout << "      z limits data: " 
            << Form("%8.5f",FidVol::gFarZData[0]) << " " 
            << Form("%8.5f",FidVol::gFarZData[1]) << " "
            << Form("%8.5f",FidVol::gFarZData[2]) << " "
            << Form("%8.5f",FidVol::gFarZData[3]) << std::endl
            << "                 MC: " 
            << Form("%8.5f",FidVol::gFarZMC[0]) << " " 
            << Form("%8.5f",FidVol::gFarZMC[1]) << " "
            << Form("%8.5f",FidVol::gFarZMC[2]) << " " 
            << Form("%8.5f",FidVol::gFarZMC[3]) << std::endl;

  std::cout << "Evt/Trk/Shw vertex offsets:"
            << " " << Form("%8.5f",FidVol::gEvtVtxZOffset)
            << "/" << Form("%8.5f",FidVol::gTrkVtxZOffset)
            << "/" << Form("%8.5f",FidVol::gShwVtxZOffset)
            << std::endl;

}

bool FidVol::load_setter_script(std::string script)
{
  const char* paths = 
    ".:~:$SRT_PRIVATE_CONTEXT/DataUtil:$SRT_PUBLIC_CONTEXT/DataUtil";

  if ( script.find(".C") == std::string::npos ) script += ".C";
  const char* filename = gSystem->Which(paths,script.c_str());
  if ( filename ) {
    //std::cout << "found " << filename << std::endl;
    int errcode = gROOT->LoadMacro(filename);
    // gSystem->Which() returns allocated space, give it back when done
    delete [] filename;
    if ( errcode == 0 ) {
      //std::cout << "load macro" << filename << std::endl;
      return true;
    } else {
      std::cout << "failed loading " << script 
                << ", status = " << filename << std::endl;
    }
  }
  return false;
}

bool FidVol::know_setter(std::string funcname)
{
  // determine if we know how to call this function
  TInterpreter* interp = gROOT->GetInterpreter();
  void* func = interp->GetInterfaceMethod(0,funcname.c_str(),"");
  return ( (func) ? true : false );
}

void choose_infid_set(std::string setname, bool doassert)
{
  //
  // Call function that will configure "infid" parameters
  // Function name could be any of:
  //    infid_set_<setname>()
  //    <setname>()
  // And can be defined in any of:
  //    infid_sets.C
  //    infid_set_<setname>.C
  //    <setname>.C
  // Found in any of the directories (searched in this order):
  //     . ~ $SRT_PRIVATE_CONTEXT/DataUtil $SRT_PUBLIC_CONTEXT/DataUtil
  // If doassert is true and function can't be found, then assert()
  //
  // std::cout << "choose_infid_set(\"" << setname << "\")" << std::endl;
  // one-time load of generic macro that that has infid setter functions
  static bool loadGeneric = FidVol::load_setter_script("infid_sets.C");
  if ( ! loadGeneric ) {
    std::cout << "failed to find \"infid_sets.C\"" << std::endl;
  }

  const char* patterns[] = { "infid_set_%s", "%s" };
  int npatterns = sizeof(patterns)/sizeof(const char*);
  bool done = false;

  for (int i = 0; i < npatterns; ++i) {
    std::string funcname = Form(patterns[i],setname.c_str());
    if ( ! FidVol::know_setter(funcname) ) {
      // not known, try loading script if one can be found
      FidVol::load_setter_script(funcname);
    }
    if ( FidVol::know_setter(funcname) ) {
      funcname += "();";
      int psuccess = 0;
      gROOT->ProcessLine(funcname.c_str(),&psuccess);
      if ( psuccess == 0 ) {
        std::cout << "Successfully ran function \"" << funcname 
                  << "\" for set \"" << setname << "\"" << std::endl;
        done = true;
        break; // exit loop
      } else {
        std::cout << "Failed running macro \"" << funcname << "\" for set \"" 
                  << setname << "\"" << std::endl;
        if (doassert) assert(0);
      }
    } else {
      std::cout << "Found no setter function " << funcname 
                << "() for \"" << setname << "\"" << std::endl;
    }
  } // end loop over patterns
  if ( ! done ) {
    if ( doassert ) assert(0);
    else {
      std::cout << "choose_infid_set(\"" << setname << "\") FAILED" 
                << ", but user chose not to assert()" << std::endl;
    }
  }
  print_infid();
}

//
//

bool FidVol::infid_near_z(SimFlag::SimFlag_t simflg, double z)
{
  double* zcuts = gNearZData;
  if ( SimFlag::kMC == simflg ) zcuts = gNearZMC;
  if ( z < zcuts[0] ) return false;
  if ( z > zcuts[1] ) return false;
  return true;
}

bool FidVol::infid_far_z(SimFlag::SimFlag_t simflg, double z)
{
  double* zcuts = gFarZData;
  if ( SimFlag::kMC == simflg ) zcuts = gFarZMC;
  if ( z < zcuts[0] ) return false;
  if ( z > zcuts[3] ) return false;
  if ( z > zcuts[1] && z < zcuts[2] ) return false;
  return true;
}

double FidVol::infid_near_radius(double x, double y, double z)
{
  double xx = x, yy = y;
  if ( gNearFollowBeam ) {
    xx = xx - gNearX0Beam;
    yy = yy - gNearY0Beam - z*gNearDyDz;
  } else {
    xx = xx - gNearX0Z;
    yy = yy - gNearY0Z;
  }
  double r2 = xx*xx + yy*yy;
  return TMath::Sqrt(r2);
}

bool FidVol::infid_near_circle_z(double x, double y)
{
  double xx = x - gNearX0Z;
  double yy = y - gNearY0Z;
  double r2 = xx*xx + yy*yy;
  if ( r2 > gNearR*gNearR ) return false;
  return true;
}

bool FidVol::infid_near_circle_beam(double x, double y, double z)
{
  double xx = x - gNearX0Beam;
  double yy = y - gNearY0Beam - z*gNearDyDz;
  double r2 = xx*xx + yy*yy;
  if ( r2 > gNearR*gNearR ) return false;
  return true;
}

bool FidVol::infid_far_coil(double x, double y)
{
  double r2 = x*x + y*y;
  if ( r2 < gFarRinner*gFarRinner ) return false;
  return true;
}

bool FidVol::infid_far_circle(double x, double y)
{
  double r2 = x*x + y*y;
  if ( r2 > gFarRouter*gFarRouter ) return false;
  return true;
}

bool FidVol::infid_far_octagon(double x, double y)
{
  if ( TMath::Abs(x) > gFarRouter ) return false;
  if ( TMath::Abs(y) > gFarRouter ) return false;

  double u = (  x + y ) * r_sqrt2;
  double v = ( -x + y ) * r_sqrt2;

  if ( TMath::Abs(u) > gFarRouter ) return false;
  if ( TMath::Abs(v) > gFarRouter ) return false;

  return true;
}


std::string FidVol::getName() { return gName; }
void FidVol::setName(std::string name) { gName = name; }

bool   FidVol::getNearFollowBeam() { return gNearFollowBeam; }
void   FidVol::setNearFollowBeam(bool follow) { gNearFollowBeam = follow; }
double FidVol::getNearR() { return gNearR; }
void   FidVol::setNearR(double r) { gNearR = r; }
double FidVol::getNearZData(int indx)
{ if (!legal_indx(indx,2)) assert(0); return gNearZData[indx]; }
void   FidVol::setNearZData(int indx, double z)
{ if (!legal_indx(indx,2)) assert(0); gNearZData[indx] = z; }
double FidVol::getNearZMC(int indx)
{ if (!legal_indx(indx,2)) assert(0); return gNearZMC[indx]; }
void   FidVol::setNearZMC(int indx, double z)
{ if (!legal_indx(indx,2)) assert(0); gNearZMC[indx] = z; }
double FidVol::getBeamAngleRad() { return gBeamAngleRad; }
void   FidVol::setBeamAngleRad(double angle)
{ gBeamAngleRad = angle; gNearDyDz = TMath::Tan(-angle); }
double FidVol::getNearDyDz() { return gNearDyDz; }
void   FidVol::setNearDyDz(double dydz)
{ gNearDyDz = dydz; gBeamAngleRad = -TMath::ATan(dydz); }
double FidVol::getNearX0Beam() { return gNearX0Beam; }
void   FidVol::setNearX0Beam(double x0) { gNearX0Beam = x0; }
double FidVol::getNearY0Beam() { return gNearY0Beam; }
void   FidVol::setNearY0Beam(double y0) { gNearY0Beam = y0; }
double FidVol::getNearX0Z() { return gNearX0Z; }
void   FidVol::setNearX0Z(double x0) { gNearX0Z = x0; }
double FidVol::getNearY0Z() { return gNearY0Z; }
void   FidVol::setNearY0Z(double y0) { gNearY0Z = y0; }
  
bool   FidVol::getFarOctagon() { return gFarOctagon; }
void   FidVol::setFarOctagon(bool octagon) { gFarOctagon = octagon; }
bool   FidVol::getFarCoilCut() { return gFarCoilCut; }
void   FidVol::setFarCoilCut(bool coilcut) { gFarCoilCut = coilcut; }
double FidVol::getFarRinner() { return gFarRinner; }
void   FidVol::setFarRinner(double r) { gFarRinner = r; }
double FidVol::getFarRouter() { return gFarRouter; }
void   FidVol::setFarRouter(double r) { gFarRouter = r; }
double FidVol::getFarZData(int indx)
{ if (!legal_indx(indx,4)) assert(0); return gFarZData[indx]; }
void   FidVol::setFarZData(int indx, double z)
{ if (!legal_indx(indx,4)) assert(0); gFarZData[indx] = z; }
double FidVol::getFarZMC(int indx)
{ if (!legal_indx(indx,4)) assert(0); return gFarZMC[indx]; }
void   FidVol::setFarZMC(int indx, double z)
{ if (!legal_indx(indx,4)) assert(0); gFarZMC[indx] = z; }


double FidVol::getEvtVtxZOffset() { return gEvtVtxZOffset; }
void   FidVol::setEvtVtxZOffset(double zoff) { gEvtVtxZOffset = zoff; }
double FidVol::getTrkVtxZOffset() { return gTrkVtxZOffset; }
void   FidVol::setTrkVtxZOffset(double zoff) { gTrkVtxZOffset = zoff; }
double FidVol::getShwVtxZOffset() { return gShwVtxZOffset; }
void   FidVol::setShwVtxZOffset(double zoff) { gShwVtxZOffset = zoff; }

---