ThresholdCorrMap.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file   ThresholdCorrMap.cxx
// \brief  Parameterization for use in threshold and shadowing corrections
// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "Calibration/func/ThresholdCorrMap.h"

#include "TFile.h"
#include "TTree.h"

#include <cassert>
#include <cmath>

#include "cetlib/search_path.h"
#include "cetlib_except/exception.h"

namespace calib
{
  //......................................................................
  ThresholdCorrMap::ThresholdCorrMap(novadaq::cnv::DetId det,
                                     std::string const& fileName,
                                     bool FiberBrightnessMode)
  {
    std::cout << "ThresholdCorrMap fileName: " << fileName << std::endl;

    /*    std::string fname;
    std::string filePath("Calibration/");
    filePath += fileName;
    cet::search_path sp("FW_SEARCH_PATH");
    sp.find_file(filePath.c_str(), fname);

    std::cout << "sp.find_file(" << filePath.c_str() << ", " << fname << std::endl;
    */

    switch(det){
    case novadaq::cnv::kNEARDET:
      if(fileName.find("nd") == std::string::npos)//was filePath
          throw cet::exception("ThresholdCorrMap")
            << "threshold correction map file " << fileName //was filePath
          << " does not appear to be valid for current detector"
          << det;
        break;
    case novadaq::cnv::kFARDET:
      if(fileName.find("fd") == std::string::npos)//was filePath
          throw cet::exception("ThresholdCorrMap")
            << "threshold correction map file " << fileName //was filePath
          << " does not appear to be valid for current detector"
          << det;
      break;
    case novadaq::cnv::kTESTBEAM:
      if(fileName.find("tb") == std::string::npos)//was filePath
        throw cet::exception("ThresholdCorrMap")
          << "threshold correction map file " << fileName //was filePath
          << " does not appear to be valid for current detector"
          << det;
      break;

    default:
      assert(0 && "Unhandled detector, tb test");
      }
    assert(!fileName.empty());//was fname

    TFile corrFile(fileName.c_str());//was fname
    assert(!corrFile.IsZombie());

    TTree* tr = (TTree*)corrFile.Get("tree");
    int view;
    int fiberbrightness=0;
    tr->SetBranchAddress("view", &view);
    if(FiberBrightnessMode) tr->SetBranchAddress("fiberbrightness", &fiberbrightness);
    Coeffs cs;
    double cutw, cutwmax;
    tr->SetBranchAddress("Na",&fNa);
    tr->SetBranchAddress("Nb",&fNb);
    tr->SetBranchAddress("Nflag",&fNflag);
    tr->SetBranchAddress("wpow", &cs.wpow);
    tr->SetBranchAddress("cellpow", &cs.cellpow);
    tr->SetBranchAddress("modpow", &cs.modpow);
    tr->SetBranchAddress("coeff", &cs.coeff);
    tr->SetBranchAddress("dcellMax", &fCellMax);
    tr->SetBranchAddress("dcutW", &cutw);
    tr->SetBranchAddress("dcutWMax", &cutwmax);

    for(int n = 0; n < tr->GetEntries(); ++n){
      tr->GetEntry(n);
      fCoeffs[view][fiberbrightness].push_back(cs);
      fCutW[view] = cutw;
      fCutWMax[view] = cutwmax;
    }

    corrFile.Close();
  }

  //......................................................................
  ThresholdCorrMap::~ThresholdCorrMap()
  {
  }

  namespace{
    /// Raise a double to a (small) integer
    inline double ipow(double a, int b)
    {
      double ret = 1;
      for(int i = 0; i < b; ++i) ret *= a;
      return ret;
    }
  }

  //......................................................................
  double ThresholdCorrMap::ThresholdCorrAt(int view, int cell, double w) const
  {  
    return ThresholdCorrAt(view,cell,w,0);
  }

  //......................................................................
  double ThresholdCorrMap::ThresholdCorrAt(int view, int cell, double w, int fiberbrightness) const
  {
    
    double kMaxW = fCutWMax[view];
    double kCutW = fCutW[view];
    if(w < -kMaxW) w = -kMaxW;
    if(w > +kCutW) w = +kCutW;

    const int mod = cell%32; // Position within module
    double ret = 0;
   
    if(fCoeffs[view][fiberbrightness].size()){ 
      for(int a = 0; a < fNa; ++a){
        for(int b = 0; b < fNb; ++b){
          ret += fCoeffs[view][fiberbrightness][fNb*a+b].coeff*pow(w/kMaxW, a)*pow(cell/fCellMax, b);
        }
      }
      ret+=fCoeffs[view][fiberbrightness][fNa*fNb+mod].coeff;
    } 
    std::cout << "ThresholdCorrMapAt(" << view << "," << cell << "," << w << "," << fiberbrightness << ") returns " << ret << std::endl;

    return ret;

  }
}