AttenFit_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    AttenFit_module.cc
// \brief   Fit attenuation behaviour from AttenHists files
// \author  Christopher Backhouse, Dan Pershey - bckhouse@caltech.edu, pershey@caltech.edu
////////////////////////////////////////////////////////////////////////

// C/C++ includes
#include <cmath>
#include <iostream>

// ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TString.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TProfile.h"
#include "TFitResult.h"
#include "TFitResultPtr.h"
#include "TStyle.h"
#include "TSystem.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/SubRun.h"

// NOvASoft includes
#include "GeometryObjects/CellGeo.h"
#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "Calibration/func/ThresholdCorrMap.h"
#include "CalibrationDataProducts/AttenHists.h"
#include "CalibrationDataProducts/AttenProfiles.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"
#include "RunHistory/service/RunHistoryService.h"
#include "Utilities/func/MathUtil.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "Calibrator/func/AttenCurve.h"

namespace calib
{
  // Parameters for fits, should be moved to .fcl file
  int kStatsLimit = 1;
  //const double kTrueShortAtten = 289.5;

  class AttenFit: public art::EDAnalyzer
  {
  public:
    explicit AttenFit(const fhicl::ParameterSet& pset);
    ~AttenFit();

    void reconfigure(const fhicl::ParameterSet& pset);
    virtual void beginRun(const art::Run& run);
    virtual void analyze(const art::Event& evt);
    virtual void endSubRun(const art::SubRun& sr);
    virtual void endRun(const art::Run& run);
    virtual void endJob();

  protected:
    int MinBlock() const;
    int MaxBlock() const;
    int GetView(int plane, art::ServiceHandle<geo::Geometry> geom);
    void ConsolidateViews();

    caldp::AttenHists GetChannelHists(bool isData, int plane, int cell);
    caldp::AttenProfiles GetChannelProfiles(bool isData, int plane, int cell);

    TH2* GetBestHistogram(caldp::AttenHists v,caldp::AttenProfiles u);
    TH1* MedianProfile(TH2* f);
    TH1* MeanProfile(TH2* f);
    TH1* TruncatedMeanProfile(TH2* h);
    void ApplyThresholdCorrection(int view, int cell, TH1* prof);
    TF1* expFit(TH1* prof, TString opt, geo::OfflineChan chan, calib::AttenCurve*& res);
    TF1* rolloffFit(TH1* prof, TString opt, geo::OfflineChan chan, calib::AttenCurve*& res);
    TGraph* lowessFit(TH1* prof, int plane, calib::AttenCurve*& res);
    void fit_channel_prof(TH1* atten, int plane, int cell, calib::AttenCurve*& res);
    void DrawDetectorEdges(double range);
    void DrawDetectorEdges(double rangedown, double rangeup);

    double FitQuality(TH1* prof, TF1* fit) const;

    /// Ensure every cell has a CSV entry even if it couldn't be fit
    void WriteDummyCSVRow(FILE* f, int plane, int cell) const;

    double fInitEdgeX, fInitEdgeY;
    double fInitEdgeMuXMin, fInitEdgeMuXMax;
    std::string fAttenHistsMapLabel;
    std::string fPlotsDirectory;
    std::string fThresholdCorrMapFile;
    bool fProfileMode;
    bool fConsolidateViewsMode;
    bool fIncludeRolloffFit;
    bool fIncludeLowessFit;
    int fPlaneMin, fPlaneMax;
    bool fLookInSubRuns;

    bool fIsMC;

    ThresholdCorrMap* fThreshCorrMap;

    /// This is where \ref GetChannelHist stores its histograms. Map is from block
    caldp::AttenHistsMap fChannelMap;
    std::map<int, caldp::AttenProfilesMap> fChannelMapProf;
  };

  //......................................................................
  AttenFit::AttenFit(const fhicl::ParameterSet& pset)
    : EDAnalyzer(pset)
  {
    reconfigure(pset);

    fThreshCorrMap = 0;

    gStyle->SetOptStat(0);
  }

  //......................................................................
  AttenFit::~AttenFit()
  {
    delete fThreshCorrMap;
  }

  //......................................................................
  void AttenFit::reconfigure(const fhicl::ParameterSet& pset)
  {
    fAttenHistsMapLabel = pset.get<std::string>("AttenHistsMapLabel");
    fThresholdCorrMapFile = pset.get<std::string>("ThresholdCorrMapFile");
    fPlotsDirectory = pset.get<std::string>("PlotsDirectory");
    fProfileMode = pset.get<bool>("ProfileMode");
    fConsolidateViewsMode = pset.get<bool>("ConsolidateViewsMode");
    fIncludeRolloffFit = pset.get<bool>("IncludeRolloffFit");
    fIncludeLowessFit = pset.get<bool>("IncludeLowessFit");
    fPlaneMin = pset.get<int>("PlaneMin");
    fPlaneMax = pset.get<int>("PlaneMax");
    fLookInSubRuns = pset.get<bool>("LookInSubRuns");
  }

  //......................................................................
  void AttenFit::beginRun(const art::Run& run)
  {
    art::ServiceHandle<geo::Geometry> geom;
    fInitEdgeX = geom->DetHalfWidth();
    fInitEdgeY = geom->DetHalfHeight();
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      fInitEdgeX-=50;
      fInitEdgeY-=50;
      fInitEdgeMuXMin=-150;
      fInitEdgeMuXMax=50;
    }

    art::ServiceHandle<nova::dbi::RunHistoryService> rh;
    fIsMC = (rh->GetDataType() == nova::dbi::kMCOnly);
    if(!fThreshCorrMap) fThreshCorrMap = new ThresholdCorrMap(geom->DetId(),fThresholdCorrMapFile);
  }

  //......................................................................
  void AttenFit::analyze(const art::Event& evt)
  {
    // If we get given any events, make sure they match what we guessed for the
    // run.
    assert(fIsMC != evt.isRealData());
  }

  //......................................................................
  void AttenFit::endSubRun(const art::SubRun& sr)
  {
    if(!fLookInSubRuns) return;

    if(fProfileMode){
      for(int block = MinBlock(); block <= MaxBlock(); ++block){
        art::Handle<caldp::AttenProfilesMap> profsmap;
        sr.getByLabel(fAttenHistsMapLabel, TString::Format("block%d", block).Data(), profsmap);
        if(profsmap.failedToGet()) continue;

        // Accumulate across runs
        if(fChannelMapProf.find(block) == fChannelMapProf.end())
          fChannelMapProf.insert(std::make_pair(block, caldp::AttenProfilesMap(profsmap->MinW(), profsmap->MaxW())));
        auto it = fChannelMapProf.find(block);
        it->second += *profsmap;

        // No one else will ever want this. Save the memory
        sr.removeCachedProduct(profsmap);
      }
    }
    else{
      art::Handle<caldp::AttenHistsMap> histsmap;
      sr.getByLabel(fAttenHistsMapLabel, histsmap);

      // Accumulate across runs
      fChannelMap += *histsmap;
    }
  }

  //......................................................................
  void AttenFit::endRun(const art::Run& run)
  {
    // Profiles could be stored in subruns or runs, depending whether summing
    // has been done.
    for(int block = MinBlock(); block <= MaxBlock(); ++block){
      art::Handle<caldp::AttenProfilesMap> profsmap;
      run.getByLabel(fAttenHistsMapLabel, TString::Format("block%d", block).Data(), profsmap);
      
      if(profsmap.failedToGet()) continue;

      // Accumulate across runs
      if(fChannelMapProf.find(block) == fChannelMapProf.end())
        fChannelMapProf.insert(std::make_pair(block, caldp::AttenProfilesMap(profsmap->MinW(), profsmap->MaxW())));
      auto it = fChannelMapProf.find(block);
      it->second += *profsmap;

      // No one else will ever want this. Save the memory
      run.removeCachedProduct(profsmap);
    }
  }

  //......................................................................
  void AttenFit::endJob()
  {
    if(fConsolidateViewsMode) ConsolidateViews();
    
    art::ServiceHandle<geo::Geometry> geom;
    
    TFile* outfile = new TFile("Atten_Fit_Results.root","RECREATE");
    
    unsigned int plane_max = fPlaneMax;
    bool isND=geom->DetId() == novadaq::cnv::kNEARDET;
    if(fPlaneMax < 0) plane_max = geom->NPlanes()-1;
    if(fConsolidateViewsMode) plane_max = 1;
    if(fConsolidateViewsMode&&isND) plane_max = 3;
    
    const TString mcStr = fIsMC ? "_mc" : "";
    FILE* fConstsCSV = fopen("calib_atten_consts"+mcStr+".csv","w");
    FILE* fPointsCSV = fopen("calib_atten_points"+mcStr+".csv","w");
    
    fprintf(fConstsCSV, "#coeff_exp,atten_length,background,edge_low,edge_high,coeff_low,coeff_high,chisq,channel,tv\n");
    fprintf(fPointsCSV, "#w,factor,channel,tv\n");
    
    for (unsigned int plane = fPlaneMin; plane <= plane_max; plane++){
      mf::LogVerbatim("AttenFit") << "Fitting plane "
                                  << plane << "/" << plane_max;
      const int block = plane/32;
      int view = GetView(plane,geom);
      if(fConsolidateViewsMode) view = plane;
      auto it = fChannelMapProf.find(block);
      std::cout<<"it == fChannelMapProf.end()"<<(it == fChannelMapProf.end())<<std::endl;
      if(it == fChannelMapProf.end()) continue;
      caldp::AttenProfilesMap& profmap = it->second;

      TDirectory* dir = 0;
      const unsigned int cell_max = geom->Plane(plane)->Ncells();
      std::cout<<"CellMax: "<<cell_max<<std::endl;
      for(unsigned int cell = 0; cell < cell_max; ++cell){
        mf::LogVerbatim("AttenFit") << "  Fitting cell "
                                    << cell << "/" << cell_max;

        TH1* prof = 0;
        if(fProfileMode){
          if(!profmap.HasProfiles(plane, cell)){
            mf::LogWarning("AttenFit") << "Nothing to fit for plane "
                                       << plane << ", cell " << cell;
            WriteDummyCSVRow(fConstsCSV, plane, cell);
            continue;
          }       
          caldp::AttenProfiles v = profmap.GetProfiles(plane, cell);

          prof = v.WPE_corr_xy.ToTH1F();
          // if(plane==2&&cell==33){

          //   std::cout<<"Integral "<<prof->Integral()<<std::endl;
          //   std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl;
          //   std::cout<<"GetEffectiveEntries "<<prof->GetEffectiveEntries()<<std::endl;

          //   std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl;

          //   std::cout<<"TotalEntries z "<<v.WPE_corr_z.TotalEntries()<<std::endl;

          //   std::cout<<"TotalEntries "<<v.WPE_corr.TotalEntries()<<std::endl;
            
          //   assert(0);
          // }

          // TODO: figure out where these bad axis limits come from
          
          if(prof->GetXaxis()->GetXmin() == 0 ||
             prof->GetXaxis()->GetXmax() == 0 ||
             v.WPE_corr_xy.TotalEntries() < kStatsLimit){
            // std::cout<<"Stats Fail "<<v.WPE_corr_xy.TotalEntries()<<" "<<kStatsLimit<<std::endl;
            
            // std::cout<<"Integral "<<prof->Integral()<<std::endl;
            // std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl;     
            // std::cout<<"GetEffectiveEntries "<<prof->GetEffectiveEntries()<<std::endl;
            
            // std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl;
            
            // std::cout<<"TotalEntries z "<<v.WPE_corr_z.TotalEntries()<<std::endl;
            
            // std::cout<<"TotalEntries "<<v.WPE_corr.TotalEntries()<<std::endl;
          
            prof = v.WPE_corr_z.ToTH1F();
            
            if(prof->GetXaxis()->GetXmin() == 0 ||
               prof->GetXaxis()->GetXmax() == 0 ||
               v.WPE_corr_z.TotalEntries() < kStatsLimit){
              std::cout<<"Stats Double Fail "<<v.WPE_corr_z.TotalEntries()<<" "<<kStatsLimit<<std::endl;
              WriteDummyCSVRow(fConstsCSV, plane, cell);
              //continue;
            }
          }
        }
        else{
          caldp::AttenHists v = fChannelMap.GetHists(plane, cell);
          caldp::AttenProfiles u = profmap.GetProfiles(plane, cell);
          TH2* h_WPE = GetBestHistogram(v,u);
          if(!h_WPE){
            WriteDummyCSVRow(fConstsCSV, plane, cell);
            continue;
          }
          prof = MedianProfile(h_WPE);
        }

        // Never evaluated for NDOS
        if(geom->DetId() == novadaq::cnv::kFARDET ||
           geom->DetId() == novadaq::cnv::kNEARDET){
          ApplyThresholdCorrection(view, cell, prof);
        }

        calib::AttenCurve* res = calib::AttenCurve::Uninitialized(geom->DetId(), geo::OfflineChan(plane, cell));
        fit_channel_prof(prof, plane, cell, res);
        outfile->cd();
        if(!dir) dir = outfile->mkdir(Form("plane_%03d", plane));
        dir->cd();
        prof->Write(Form("Atten_Fit_%03d_%03d",plane,cell));
        // Encoding to indicate it's the backup averaged constants
        if(fConsolidateViewsMode && !fIsMC)
          res->chan = geo::OfflineChan(1000+view, cell);
        res->WriteToCSVs(fConstsCSV, fPointsCSV, fIsMC);

        delete res;
        delete prof;
      }
    }
    fclose(fConstsCSV);
    fclose(fPointsCSV);
  }

  //......................................................................
  int AttenFit::MinBlock() const
  {
    return fPlaneMin/32;
  }

  //......................................................................
  int AttenFit::MaxBlock() const
  {
    if(fPlaneMax < 0) return 27;
    return fPlaneMax/32;
  }
  //......................................................................
  int AttenFit::GetView(int plane, art::ServiceHandle<geo::Geometry> geom){//get the view for the current combination of view and plane (x,y,mux, or muy)
    int view=0;
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      int firstPlane=geom->FirstPlaneInMuonCatcher();
      int isMuonCatcher=firstPlane <= plane;
      view = geom->Plane(plane)->View()+(2*isMuonCatcher);
    }
    else{
      view = geom->Plane(plane)->View();
    } 
    return view;
  }
  //......................................................................
  void AttenFit::ConsolidateViews()
  {
    art::ServiceHandle<geo::Geometry> geom;

    // Has consolidation already been done by an earlier step?
    bool isConsolidated = true;
    bool isND=geom->DetId() == novadaq::cnv::kNEARDET;
    if(fChannelMapProf.find(0) != fChannelMapProf.end()){
      const std::vector<geo::OfflineChan> chans = fChannelMapProf.find(0)->second.Channels();
      for(geo::OfflineChan chan: chans){
        if((chan.Plane() > 3 && isND) || (chan.Plane() > 1 && isND==false)){
          isConsolidated = false;
          break;
        }
      }
    }

    if(isConsolidated){
      mf::LogInfo("AttenFit") << "Input file appears to already be consolidated by view. Skipping that step.";
      return;
    }

    // Need to do it for ourselves
    mf::LogInfo("AttenFit") << "Input file is not consolidated by view. Consolidating";

    // All x-planes will go into plane 0 of this, y to 1.
    caldp::AttenProfilesMap* sum = 0;

    for(int block = MinBlock(); block <= MaxBlock(); ++block){
      auto it = fChannelMapProf.find(block);
      if(it == fChannelMapProf.end()) continue;
      caldp::AttenProfilesMap& profmap = it->second;
      const std::vector<geo::OfflineChan> chans = profmap.Channels();

      for(geo::OfflineChan chan: chans){
        const caldp::AttenProfiles& profs = profmap.GetProfiles(chan);
        if(!sum) sum = new caldp::AttenProfilesMap(profmap.MinW(),
                                                   profmap.MaxW());

        const int view = GetView(chan.Plane(),geom);
        
        sum->GetProfiles(view, chan.Cell()) += profs;
      } // end for chan

      profmap.Clear();
    } // end for block

    fChannelMapProf.insert(std::make_pair(0, *sum));
    delete sum;
  }

  //......................................................................
  double AttenFit::FitQuality(TH1* prof, TF1* fit) const
  {
    // chisq lets bad fits to low stats data through, and cuts out good(ish)
    // fits to high stats data. Try this metric, mean fractional deviation of
    // the data from the fit (in quadrature).

    double xmin, xmax;
    fit->GetRange(xmin, xmax);
    int tot = 0;
    double diff = 0;
//std::cout << "xmin = " << xmin << " , xmax = " << xmax << std::endl;
    for(int i = 0; i < prof->GetNbinsX()+2; ++i){
      const double x = prof->GetXaxis()->GetBinCenter(i);
//std::cout << "i = " << i << " , x = " << x <<  " ,xmin = " << xmin << " , xmax = " << xmax << " , Fit value = " << fit->Eval(x) << " , Bin content = " << prof->GetBinContent(i) << ", tot = " << tot << " and diff = " << diff << std::endl;
      if(x < xmin || x > xmax) {
//              std::cout << "i = " << i << " , x = " << x <<  " ,xmin = " << xmin << " , xmax = " << xmax << " , Fit value = " << fit->Eval(x) << " , Bin content = " << prof->GetBinContent(i) << ", tot = " << tot << " and diff = " << diff << std::endl;
                continue;
        }
      ++tot;
      diff += util::sqr((fit->Eval(x)-prof->GetBinContent(i))/fit->Eval(x));
    }

    const double ret = sqrt(diff/tot);
    if(std::isnan(ret) || std::isinf(ret)) {
//      std::cout << "xmin = " << xmin << " , xmax = " << xmax << " , tot = " << tot << " , diff = " << diff << " and ret = " << ret << std::endl;
        return 2;
        }
    return ret;
  }

  //......................................................................
  TH2* AttenFit::GetBestHistogram(caldp::AttenHists v, caldp::AttenProfiles u)
  {
    if(!v.WPE) return 0;
    if(u.WPE_corr_xy.TotalEntries() > kStatsLimit) return v.WPE_corr_xy;
    if(u.WPE_corr_z.TotalEntries() > kStatsLimit) return v.WPE_corr_z;
    if(u.WPE_corr.TotalEntries() > kStatsLimit) return v.WPE_corr;
    return 0;
  }

  //......................................................................
  TH1* AttenFit::MeanProfile(TH2* h)
  {
    TProfile* p = h->ProfileX();

    TAxis* xax = h->GetXaxis();
    TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
    ret->GetXaxis()->SetTitle(xax->GetTitle());
    for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
      ret->SetBinContent(ix, p->GetBinContent(ix));
      ret->SetBinError(ix, p->GetBinError(ix));
    }
    return ret;
  }

  //......................................................................
  TH1* AttenFit::MedianProfile(TH2* h)
  {
    art::ServiceHandle<geo::Geometry> geom;
    TAxis* xax = h->GetXaxis();
    TAxis* yax = h->GetYaxis();
    TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
    ret->GetXaxis()->SetTitle(xax->GetTitle());
    for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
      double tot = 0;
      // Make sure to include overflow bins
      for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
        tot += h->GetBinContent(ix, iy);
      }
      double seen = 0;
      for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
        seen += h->GetBinContent(ix, iy);
        if(seen > tot/2){
          const double z = h->GetBinContent(ix, iy);
          // Interpolate
          const double frac = 1-(seen-tot/2)/z;
          ret->SetBinContent(ix, yax->GetBinLowEdge(iy)+frac*yax->GetBinWidth(iy));
          // Intro of http://www.jstor.org/stable/2286545
          ret->SetBinError(ix, sqrt(tot)*yax->GetBinWidth(iy)/(2*z));
          break;
        }
      } // end for iy
    } // end for ix
    return ret;
  }

  //......................................................................
  TH1* AttenFit::TruncatedMeanProfile(TH2* h)
  {
    TAxis* xax = h->GetXaxis();
    TAxis* yax = h->GetYaxis();
    TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
    ret->GetXaxis()->SetTitle(xax->GetTitle());
    for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
      double tot = 0;
      // Make sure to include overflow bins
      for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
        tot += h->GetBinContent(ix, iy);
      }

      double avg = 0;
      double rms = 0;
      double seen = 0;
      double inc = 0;
      for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
        const double z = h->GetBinContent(ix, iy);
        if(!z) continue;

        seen += z;

        // Interpolation
        double weight_lo = (seen-(.25*tot))/z;
        double weight_hi = 1-(seen-(.75*tot))/z;
        if(weight_lo > 1) weight_lo = 1;
        if(weight_hi > 1) weight_hi = 1;
        double weight = weight_lo+weight_hi-1;
        if(weight < 0) weight = 0;

        inc += z*weight;
        const double c = yax->GetBinCenter(iy);
        avg += z*weight*c;
        rms += z*weight*util::sqr(c);
      } // end for iy

      // In principle this should be problematic because of floating point errors
      assert(inc == tot/2);

      if(inc){
        avg /= inc;
        rms /= inc;
        rms -= util::sqr(avg);
        assert(!std::isinf(avg) && !std::isnan(avg));
        assert(!std::isinf(rms) && !std::isnan(rms));
        assert(rms >= 0);
        ret->SetBinContent(ix, avg);
        ret->SetBinError(ix, sqrt(rms)/sqrt(tot/2));
      }
      else{
        ret->SetBinContent(ix, 0);
        ret->SetBinError(ix, 1e10);
      }
    }

    return ret;
  }

  //......................................................................
  void AttenFit::ApplyThresholdCorrection(int view, int cell, TH1* prof)
  {
//std::cout << "////////////////////////////// Apply Threshold correction /////////////////////////////" << std::endl;
    for(int i = 0; i < prof->GetNbinsX()+2; ++i){
      const double w = prof->GetXaxis()->GetBinCenter(i);
      const double y = prof->GetBinContent(i);
      const double corr = fThreshCorrMap->ThresholdCorrAt(view, cell, w);
//std::cout << "i = " << i << " , Bin center = " << w << " , Bincontent = " << y << " , threshold corr = " << corr << " , view = " << view << " , cell " << cell << " , new Bin content " << y/corr << " , old bin error = "<< prof->GetBinError(i) << " and new bin error " << prof->GetBinError(i)/corr << std::endl;
      prof->SetBinContent(i, y/corr);
      prof->SetBinError(i, prof->GetBinError(i)/corr);
    }
//std::cout << "////////////////////////////// Done Threshold correction /////////////////////////////" << std::endl;
  }

  //......................................................................
  TF1* AttenFit::expFit(TH1* prof, TString opt, geo::OfflineChan chan,
                        calib::AttenCurve*& res)
  {
    art::ServiceHandle<geo::Geometry> geom;

    int view = GetView(chan.Plane(),geom);
    double range = geom->Plane(chan.Plane())->Cell(chan.Cell())->HalfL();
    if(fConsolidateViewsMode){ 
      view = (chan.Plane());
      range = (view == geo::kX) ? geom->DetHalfWidth() : geom->DetHalfHeight();
    }
    res->cell_length = 2*range;
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      range -= 50; 
    }
    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      res->cell_length =range+50+fInitEdgeMuXMax;
    }
    res->coeff_low = 0;
    res->coeff_high = 0;
    res->edge_low = 0;
    res->edge_high = 0;

    struct ExpFitAdaptor
    {
      ExpFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
      double operator()(double* xs, double* ps)
      {
        fRes->coeff_exp = ps[0];
        fRes->atten_length = ps[1];
        fRes->background = ps[2];
        return fRes->MeanPEPerCmAt(xs[0]);
      }
      calib::AttenCurve* fRes;
    };
    ExpFitAdaptor* adapt = new ExpFitAdaptor(res);

    TString fitStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
    TString fit_name = fitStrArr[view];
    double rangeup=range;
    if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
      rangeup=fInitEdgeMuXMax;
    }

    TF1* fit = new TF1(fit_name, adapt, -range, +rangeup, 3);
    // Empirically, the fits end up close to here
    fit->SetParameter(0, 11);
    if(geom->DetId() == novadaq::cnv::kNEARDET){
      fit->SetParameter(0, 43);
      fit->SetParLimits(0, 0, 1000);
      fit->SetParLimits(1, 0, 1000);
      fit->SetParLimits(2, 0, 1000);
    }    
    if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
      fit->SetParameter(0, 70);
    }
    fit->SetParLimits(0, 0, 1000);
    fit->SetParLimits(1, 0, 1000);
    fit->SetParLimits(2, 0, 1000);

    fit->SetParameter(1, 400);//kTrueShortAtten);
    fit->SetParameter(2, 7.5);
    fit->SetParName(0, "norm");
    fit->SetParName(1, "atten");
    fit->SetParName(2, "bkgd");

    fit->SetRange(-range, +rangeup);

    res->chisq = 0; // Need to count as calibrated during fit
    TFitResultPtr fr = prof->Fit(fit, "RS"+opt);
    //std::cout<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell()<<std::endl;
    //std::cout<<fit->GetParameter(0)<<" "<<fit->GetParameter(1)<<" "<<fit->GetParameter(2)<<std::endl;
    
    if(FitQuality(prof, fit)>0.2 && view==2){
      //std::cout<<"Default "<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell()<<std::endl;
      //70 400 7.5
      fit->SetParameter(0, 43);
      fit->SetParameter(1, 452);//kTrueShortAtten);                                                                                                                                                                                          
      fit->SetParameter(2, 7.6);
      res->atten_length = 73;
      res->coeff_exp    = 452;
      res->background   = 7.6;
      res->chisq        = FitQuality(prof, fit);
      //res->chisq        = 0;

      return fit;
    }
    
    res->atten_length = fit->GetParameter(0);
    res->coeff_exp    = fit->GetParameter(1);
    res->background   = fit->GetParameter(2);
    res->chisq        = FitQuality(prof, fit);


    //    const int dof = prof->FindBin(+range)-prof->FindBin(-range)-3;
    //    res->chisq        =(fr == 0) ? fr->Chi2()/dof : -1;
    
    return fit;
  }

  //......................................................................
  TF1* AttenFit::rolloffFit(TH1* prof, TString opt, geo::OfflineChan chan,
                            calib::AttenCurve*& res)
  {
    struct RolloffFitAdaptor
    {
      RolloffFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
      double operator()(double* xs, double* ps)
      {
        fRes->coeff_low = ps[0];
        fRes->coeff_high = ps[1];
        fRes->edge_low = ps[2];
        fRes->edge_high = ps[3];
        return fRes->MeanPEPerCmAt(xs[0]);
      }
      calib::AttenCurve* fRes;
    };
    RolloffFitAdaptor* adapt = new RolloffFitAdaptor(res);

    art::ServiceHandle<geo::Geometry> geom;
    int view = GetView(chan.Plane(),geom);
    if(fConsolidateViewsMode) view = (chan.Plane());

    const double edge = (view == geo::kX) ? fInitEdgeX : fInitEdgeY;
    const double range = (view == geo::kX) ? geom->DetHalfWidth() : geom->DetHalfHeight();
    double edgeup=edge;
    double rangeup=range;
    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      edgeup=fInitEdgeMuXMax;
      rangeup=fInitEdgeMuXMax+50;
    }
    
    TF1* fit = new TF1(view == geo::kX ? "rollFitX" : "rollFitY",
                       adapt, -range, +rangeup, 4);
    fit->SetParName(0, "coeffL");
    fit->SetParName(1, "coeffR");
    fit->SetParName(2, "edgeL");
    fit->SetParName(3, "edgeR");
    fit->SetParameter(0, 0);
    fit->SetParameter(1, 0);
    fit->SetParameter(2, -edge);
    fit->SetParameter(3, edgeup);
    fit->SetRange(-range, +rangeup);

    if(fIncludeRolloffFit&&view!=2){
      res->chisq = 0; // Need to count as calibrated during fit
      TFitResultPtr fr = prof->Fit(fit, "RS"+opt);

      res->coeff_low  = fit->GetParameter(0);
      res->coeff_high = fit->GetParameter(1);
      res->edge_low   = fit->GetParameter(2);
      res->edge_high  = fit->GetParameter(3);

      res->chisq      = FitQuality(prof, fit);

      //      const int dof = prof->FindBin(+range)-prof->FindBin(-range)-7;
      //      res->chisq      = (fr == 0) ? fr->Chi2()/dof : -1;
    }

    return fit;
  }

  //......................................................................
  TGraph* AttenFit::lowessFit(TH1* prof, int plane, calib::AttenCurve*& res)
  {
    art::ServiceHandle<geo::Geometry> geom;
    int view = GetView(plane, geom);
    TGraph* ret = new TGraph;
    double x0 = (view == geo::kX) ? -geom->DetHalfHeight() : -geom->DetHalfWidth();
    double x1 = (view == geo::kX) ? +geom->DetHalfHeight() : +geom->DetHalfWidth();
    const int kNumControlPts = 20;
    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      x0 =  -geom->DetHalfHeight();
      x1 = fInitEdgeMuXMax+50;
    }
    // Range of measurements that influence each control point
    const double kRange = 1.5*(x1-x0)/kNumControlPts;

    for(int n = 0; n < kNumControlPts; ++n){
      const double x = x0+double(n)/(kNumControlPts-1)*(x1-x0);

      std::vector<double> xs, ys, ws;

      for(int i = 0; i < prof->GetNbinsX()+2; ++i){
        const double xi = prof->GetBinCenter(i);
        if(xi < x0 || xi > x1) continue;
        if(fabs(xi-x) > kRange) continue;

        xs.push_back(xi);
        ys.push_back(prof->GetBinContent(i));
        // Apparently "tricube" weights are traditional
        ws.push_back(util::cube(1-util::cube(fabs(xi-x)/kRange)));
      }

      double m, c;

      util::LinFit(xs, ys, ws, m, c);

      const double y = m*x+c;

      ret->SetPoint(n, x, y);

      res->AddInterpPoint(x, y);
    }
    ret->SetLineWidth(2);
    ret->SetMarkerColor(kRed);
    ret->SetLineColor(kRed);
    ret->SetMarkerStyle(kFullCircle);
    return ret;
  }

  //......................................................................
  void AttenFit::DrawDetectorEdges(double range)
  {
    for(int sign = -1; sign <= +1; sign += 2){
      TGraph* g = new TGraph;
      g->SetPoint(0, sign*range, 0);
      g->SetPoint(1, sign*range, 1e4);
      g->SetLineStyle(7);
      g->Draw("l same");
    }
  }

  //......................................................................
  void AttenFit::DrawDetectorEdges(double rangedown, double rangeup)
  {
      TGraph* g = new TGraph;
      g->SetPoint(0,  rangedown, 0);
      g->SetPoint(1, rangedown, 1e4);
      g->SetLineStyle(7);
      g->Draw("l same");
      TGraph* g2 = new TGraph;
      g2->SetPoint(0,  rangeup, 0);
      g2->SetPoint(1, rangeup, 1e4);
      g2->SetLineStyle(7);
      g2->Draw("l same");



  }

  //......................................................................
  void AttenFit::fit_channel_prof(TH1* atten, int plane, int cell, calib::AttenCurve*& res)
  {

    art::ServiceHandle<geo::Geometry> geom;
    int view = GetView(plane,geom);
    if(fConsolidateViewsMode) view = (plane);


    const TString dataMCStr = fIsMC ? "_mc" : "_data";
    const TString viewStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
    const TString viewStr =viewStrArr[view];

    res->initialized = true;

    double cell_length = (view == geo::kX) ? 2*geom->DetHalfHeight() : 2*geom->DetHalfWidth();
    res->cell_length = cell_length;
    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      res->cell_length =(cell_length/2)+50+fInitEdgeMuXMax;
    }

    TString suffix = dataMCStr+viewStr;
    if(cell >= 0)
      suffix += TString::Format("_%03d_%03d", plane, cell);

    TString quiet;
    if(cell >= 0) quiet = "Q";

    // TString title = "NDOS cosmic ";
    TString title;
    if (geom->DetId() == novadaq::cnv::kFARDET) title = "FD cosmic ";
    if (geom->DetId() == novadaq::cnv::kNEARDET) title = "ND cosmic ";
    if (geom->DetId() == novadaq::cnv::kNDOS) title = "NDOS cosmic ";
    if (fIsMC) title += "Monte Carlo"; else title += "data";
   
    TString viewStrArrVertHor[4]={"vertical","horizontal","vertical","horizontal"};
    TString viewStrVertHor =viewStrArrVertHor[view];
    
    title += TString::Format(" - plane %d (%s), cell %d",
                             plane,
                             viewStrVertHor.Data(),
                             cell);

    new TCanvas;

    TF1* fit = expFit(atten, quiet, geo::OfflineChan(plane, cell), res);

    atten->Draw();
    atten->SetTitle(title+";Distance from center (cm);Mean PE / cm");

    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      DrawDetectorEdges(fInitEdgeMuXMin-50,fInitEdgeMuXMax+50);
    }
    else{
      DrawDetectorEdges(cell_length/2);
    }

    const double attenErr = fit->GetParError(0); /* unused */

    atten->GetYaxis()->SetRangeUser(0, 120);

    TLegend* leg = new TLegend(.3, .15, .7, .35);
    leg->SetFillStyle(0);

    leg->AddEntry((TObject*)0, TString::Format("Attenuation length = %d#pm%.1lf cm", int(res->atten_length), attenErr), "");
    leg->AddEntry((TObject*)0, TString::Format("Background = %.1lf", res->background), "");
    leg->Draw();

    std::string plotsDirectory = fPlotsDirectory;
    if(!fPlotsDirectory.empty()){
      plotsDirectory = TString::Format("%s/%03d", fPlotsDirectory.c_str(), plane);
      gSystem->mkdir(plotsDirectory.c_str());
      gPad->Print(plotsDirectory+"/atten"+suffix+".eps");
    }

    new TCanvas;
    TH1D* profCorr = (TH1D*)atten->Clone();
    profCorr->GetYaxis()->SetTitle("Ratio to exponential fit");
    fit->SetRange(-1.1*cell_length, +1.1*cell_length);
    profCorr->Divide(fit);
    profCorr->Draw();

    TF1* rollFit = rolloffFit(atten, quiet, geo::OfflineChan(plane, cell), res);

    if(res->coeff_exp == 0)
      mf::LogWarning("AttenFit") << "Fit failed: plane = " << plane << ", cell = " << cell;

    profCorr->GetYaxis()->SetRangeUser(.4, 1.1);

    leg = new TLegend(.3, .15, .7, .35);
    leg->SetFillStyle(0);

    leg->AddEntry((TObject*)0, TString::Format("Left edge = %.1lf cm", res->edge_low), "");
    leg->AddEntry((TObject*)0, TString::Format("Right edge = %.1lf cm", res->edge_high), "");
    leg->Draw();

    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      DrawDetectorEdges(fInitEdgeMuXMin-50,fInitEdgeMuXMax+50);
    }
    else{
      DrawDetectorEdges(cell_length/2);
    }

    if(!plotsDirectory.empty()){
      gPad->Print(plotsDirectory+"/rolloff"+suffix+".eps");
    }


    new TCanvas;
    TH1* profFullCorr = (TH1*)atten->Clone();
    profFullCorr->GetYaxis()->SetTitle("Ratio to full fit");
    rollFit->SetRange(-800, +800);
    if(geom->DetId() == novadaq::cnv::kNEARDET){
    rollFit->SetRange(-250, +250);
    }
    profFullCorr->Divide(rollFit);
    profFullCorr->GetYaxis()->SetRangeUser(0.9, 1.1);
    profFullCorr->Draw();

    TGraph* one = new TGraph;
    one->SetPoint(0, -1e4, 1);
    one->SetPoint(1, +1e4, 1);
    one->SetLineStyle(7);
    one->Draw("l same");

    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      DrawDetectorEdges(fInitEdgeMuXMin-50,fInitEdgeMuXMax+50);
    }
    else{
      DrawDetectorEdges(cell_length/2);
    }

    if(fIncludeLowessFit){
      TGraph* lowess = lowessFit(profFullCorr, plane, res);
      lowess->Draw("lp same");
    }

    res->initialized = true;

    if(!plotsDirectory.empty()){
      gPad->Print(plotsDirectory+"/residual"+suffix+".eps");
    }


    new TCanvas;
    TH1* plotProf = (TH1*)atten->Clone();
    plotProf->Draw();

    struct FullFitAdaptor
    {
      FullFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
      double operator()(double* xs, double* ps)
      {
        return fRes->MeanPEPerCmAt(xs[0]);
      }
      calib::AttenCurve* fRes;
    } adapt(res);

    const double range = cell_length/2;
    double rangeup=range;
    if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
      rangeup=fInitEdgeMuXMax+50;
    }
    TF1* fullFit = new TF1("fullFit", adapt, -range, +rangeup, 0);
    fullFit->SetLineColor(kBlue);
    fullFit->SetRange(-range, +rangeup);
    fullFit->Draw("same");

    plotProf->SetLineColor(kBlack);
    plotProf->GetYaxis()->SetTitle("Mean PE / cm");
    plotProf->GetYaxis()->SetRangeUser(0, 150);

    // HACKs for making a pretty plot for blessing
    plotProf->GetXaxis()->CenterTitle();
    plotProf->GetYaxis()->CenterTitle();
    plotProf->GetYaxis()->SetTitleSize(plotProf->GetXaxis()->GetTitleSize());
    plotProf->GetYaxis()->SetTitleOffset(plotProf->GetXaxis()->GetTitleOffset());

    if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
      DrawDetectorEdges(fInitEdgeMuXMin-50,fInitEdgeMuXMax+50);
    }
    else{
      DrawDetectorEdges(cell_length/2);
    }

    if(!plotsDirectory.empty()){
      gPad->Print(plotsDirectory+"/totfit"+suffix+".eps");
      gPad->Print(plotsDirectory+"/totfit"+suffix+".png");
    }

    new TCanvas;
    TH1* resid = new TH1F("", ";W (cm);Ratio to total fit", atten->GetNbinsX(), atten->GetXaxis()->GetXmin(), atten->GetXaxis()->GetXmax());

    for(int n = 0; n < atten->GetNbinsX()+2; ++n){
      resid->SetBinContent(n, atten->GetBinContent(n)/res->MeanPEPerCmAt(atten->GetBinCenter(n)));
    }
    resid->GetYaxis()->SetRangeUser(.9, 1.1);
    resid->Draw();
    one->Draw("l same");

    if(!plotsDirectory.empty()){
      gPad->Print(plotsDirectory+"/residual_final"+suffix+".eps");
    }

    delete fit;
    delete rollFit;
    delete profCorr;
    delete profFullCorr;
    delete one;
    delete leg;
    delete plotProf;
    delete fullFit;
    delete resid;
  }

  //......................................................................
  void AttenFit::WriteDummyCSVRow(FILE* f, int plane, int cell) const
  {
    const int vldTime = 1; // 1970

    geo::OfflineChan chan(plane, cell);

    if(fConsolidateViewsMode && !fIsMC){
      assert(plane == 0 || plane == 1 || plane == 2 || plane == 3);
      chan = geo::OfflineChan(1000+plane, cell);
    }

    fprintf(f, "-1,-1,-1,-1,-1,-1,-1,1000,%d,%d\n",
            chan.ToDBValidityChan(), vldTime);
  }

  DEFINE_ART_MODULE(AttenFit)

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