AttenCurve.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   Store the results of the attenuation fits
/// \author  bckhouse@caltech.edu
/// \date
/////////////////////////////////////////////////////////////////////////
// C/C++ includes
#include <cassert>
#include <cstdio>
#include <iostream>
#include <set>

#include "TString.h"

// Framework includes
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvASoft includes
#include "Calibrator/func/AttenCurve.h"

#include "GeometryObjects/GeometryBase.h"

namespace calib
{
  //......................................................................
  AttenCurve* AttenCurve::Uninitialized(int det, geo::OfflineChan chan)
  {
    AttenCurve* res = new AttenCurve;
    res->det = det;
    res->chan = chan;
    res->initialized = false;
    res->center_offset = 0;
    return res;
  }

  template<class T> T FOURTH(T x){return x*x*x*x;}
  //......................................................................
  float AttenCurve::MeanPEPerCmAt(double w) const
  {
    return MeanPEPerCmAt(w,"");
  }

  //......................................................................
  float AttenCurve::MeanPEPerCmAt(double w, std::string version) const
  {
    assert(initialized);

    //Now check just once in Calibrator instead
    //if(!IsCalibrated()) return -1;

    // TODO: take account of center_offset here for muon catcher planes. Need
    // to fix the todo below at the same time.

    //The attenuation formulae has a minor correction relative to tag v6 and earlier. 
    //Backwards compatibility is ensured in calibrator by passing the tag version to 
    //this function. Users can also call this function directly without specifying a
    //tag in which case they get the latest, corrected, version of the attenuation
    //formulae.

    // Distance to travel the long way round the fibre
    double wlong;
    if(version=="v1"||version=="v2"||version=="v3"||version=="v3.1"||version=="v4"||version=="v5"||version=="v5.1"||version=="v6"){
      wlong = 1.5*cell_length+w;
    }
    else{
      wlong = cell_length+w;
    }
    
    const double expPart = coeff_exp*(exp(w/atten_length)+
                                      exp(-wlong/atten_length))+background;

    double edgePart = 1;
    if(w < -edge_low) edgePart -= coeff_low *FOURTH(w+edge_low);
    if(w > edge_high) edgePart -= coeff_high*FOURTH(w-edge_high);

    double interpPart = 1;
    if(interp_pts.size() >= 2){
      typedef std::vector<Pt>::const_iterator it_t;

      // Find the point after w
      it_t itNext = std::lower_bound(interp_pts.begin(), interp_pts.end(), w);
      if(itNext == interp_pts.end()) --itNext;
      if(itNext == interp_pts.begin()) ++itNext;
      // And the point before
      it_t itPrev = itNext; --itPrev;

      // Interpolate
      const double w0 = itPrev->w; const double y0 = itPrev->factor;
      const double w1 = itNext->w; const double y1 = itNext->factor;

      interpPart = ((w1-w)*y0+(w-w0)*y1)/(w1-w0);
    }

    return expPart*edgePart*interpPart;
  }

  //......................................................................
  bool AttenCurve::IsCalibrated() const
  {
    // This isn't actually a chisq, but more an average fractional deviation of
    // the fit from the data. 0.2 is the eyeballed position of the tails of the
    // distribution in FD data.
    return chisq >= 0 && chisq < 0.2;

    // Checking the LOWESS points is another way to achieve something
    // similar. But this implementation tends to give to much weight to
    // problems that only occur in single locations.

    // for(unsigned int n = 0; n < interp_pts.size(); ++n){
    //   // Don't let the interpolated corrections pull more than 15% off the
    //   // original fit. 15% is arbitrary, based on a desire to calibrate a
    //   // reasonable number of channels.
    //   if(fabs(interp_pts[n].factor-1) > .15){
    //     // Don't count interpolations where just the final point is off as bad
    //     if(n > 0 && n < interp_pts.size()-1)
    //       return false;
    //   }
    // }
    // return true;
  }

  //......................................................................
  void AttenCurve::AddInterpPoint(float w, float factor)
  {
    Pt pt = {w, factor};
    interp_pts.push_back(pt);
  }

  //......................................................................
  void AttenCurve::WriteToCSVs(FILE* fConsts, FILE* fPoints, bool mc) const
  {
    const int vldTime = 1; // 1970


    assert(initialized);
    fprintf(fConsts, "%g,%g,%g,%g,%g,%g,%g,%g,%d,%d\n",
            coeff_exp, atten_length, background,
            edge_low, edge_high, coeff_low, coeff_high, chisq,
            chan.ToDBValidityChan(), vldTime);

    const int N = interp_pts.size();
    for(int n = 0; n < N; ++n){
      fprintf(fPoints, "%g,%g,%d,%d\n",
              interp_pts[n].w, interp_pts[n].factor,
              1000*chan.ToDBValidityChan()+n, vldTime);
    }
  }

  //......................................................................
  std::ostream& operator<<(std::ostream& os, const AttenCurve& res)
  {
    os << "Attenuation fit results for "
       << novadaq::cnv::DetInfo::GetName(res.det)
       << ", " << res.chan << ": ";

    if(!res.initialized){
      os << " Uninitialized." << std::endl;
      return os;
    }

    os << res.atten_length << "cm " << res.coeff_exp << " amplitude on "
       << res.background << " background. Rolloffs start at -"
       << res.edge_low << "cm and +" << res.edge_high << "cm with size "
       << res.coeff_low << " and " << res.coeff_high
       << " fit quality " << res.chisq << std::endl;
    os << "Interpolation points are: [ ";
    const int N = res.interp_pts.size();
    for(int n = 0; n < N; ++n)
      os << "(" << res.interp_pts[n].w << ", "
         << res.interp_pts[n].factor << ") ";
    os << "]" << std::endl;
    return os;
  }

} // end namespace calib
//////////////////////////////////////////////////////////////////////////