CherenkovCalc.cxx

Go to the documentation of this file.

#include <string.h>
#include <signal.h>
#include <vector>

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TString.h"
#include "TFile.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TMinuit.h"
#include "TRandom3.h"

#include "consts.h"
#include "CherenkovCalc.h"
#include "Plot.h"

// Remove for a quick and inaccurate fit
#define FULL

const bool float_ndfx = false,
           float_ndfy = false,

           float_fdfx = true,
           float_fdfy = true,

           float_ec   = false,

           float_attb = false,
           float_att1 = false,
           float_att2 = false,

           // So far most fits have preferred a value of 1, so leave it
           // Based on Kurarray data, expect ~0.99 efficiency.
           float_rngeff = false;

static TRandom3 R(0); // 0 means random seed

using std::vector;

// We copy the pointers to the samples into these static variables so that
// MINUIT's objective function can get at them.
static vector<CherenkovCalc*> sfdSamples, sndSamples;

static const unsigned int npar = 10;

//#define USEPROTONWEIGHT

#ifdef USEPROTONWEIGHT
  // Weight true proton hits in the MC by this amount, to crudely adjust for
  // cross section uncertainties.
  const double PROTONWEIGHT = 1;
#endif

// The minimum accepted reconstructed path through the cell.
const double MINPATHLENGTH = 1; // cm

// For speed, don't accept MC events after this many, per sample.
static const unsigned int maxMChitsPERw = 1e7;

// In order to prevent one sample from overpowering the rest (primarily this
// is a concern for FD cosmics), limit the number of data events per sample.
// A discussion of statistics philosophy should appear here.
static const unsigned int maxDATAhitsPERw = 1e6;

// How many hits an MC bin has to have using the initial parameters for that
// bin to be used
static const int stat_req = 1000;

// Warn if a bin's contents goes below this level as the fit proceeds
static const int stat_warn_level = 100;

// If improvement is less than this after SCAN or SIMPLEX, or whatnot, don't
// trigger another round of trying to improve.
static const double epsilonll = 1e-4;

// Project 'in' along XZ to 'out'.  The user is responsible for making the
// number of X and Z bins match between the two histograms.  All overflow bins
// are used.  This is a replacement for ROOT's insane ProjectionX().
static void saneProjectionXZ(TH2D * out, TH3D * in)
{
  double sumsum = 0;
  for(int bz = 0; bz <= in->GetNbinsZ()+1; bz++){
    for(int b = 0; b <= in->GetNbinsX()+1; b++){
      double sum = 0;
      for(int by = 0; by <= in->GetNbinsY()+1; by++)
        sum += in->GetBinContent(b, by, bz);
      out->SetBinContent(b, bz, sum);
      sumsum += sum;
    }
  }
  out->SetEntries(sumsum);
}

CherenkovCalc::CherenkovCalc(const char * const dataFileName,
                             const char * const mcFileName,
                             const char * const dirName,
                             const char * const plot_name_,
                             const vector<double> & binX,
                             const vector<double> & binY,
                             const vector<double> & binZ,
                             const bool isFD_)
: isFD(isFD_), plot_name(plot_name_)
{
  if(isFD) cellhalflength = 1549.4/2;
  else     cellhalflength = 399.28/2;

  nattensteps = (cellhalflength*4 + maxpigtail)/atten_step + 1;

  _dataFile = new TFile(dataFileName);
  _mcFile   = new TFile(mcFileName);

  if(_dataFile == NULL || _dataFile->IsZombie()){
    fprintf(stderr, "Failed to open %s\n", dataFileName);
    _exit(1);
  }
  if(_mcFile == NULL || _mcFile->IsZombie()){
    fprintf(stderr, "Failed to open %s\n", mcFileName);
    _exit(1);
  }

  _dirName = dirName;

  if(binX.size() < 2 || binY.size() < 2 || binZ.size() < 2){
    fprintf(stderr, "Must have at least two bin boundaries per dimension\n");
    _exit(1);
  }

  ndistbins = nx = binX.size()-1;
  npebins   = ny = binY.size()-1;
  nwbins    = nz = binZ.size()-1;

  distbins = binX;
  pebins = binY;
  wbins = binZ;

  // For smearing
  extpebins = binY;
  extpebins.insert(extpebins.begin(), pebins[0] - (pebins[1] - pebins[0]));
  extpebins.push_back(pebins[ny] + pebins[ny] - pebins[ny-1]);

  const char * const title = ";Distance To Stop (cm);PE/Path (cm^{-1})";

  _mcX   = new TH3D(Form("mcX_%s_%s", mcFileName, dirName),
                    title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
  _mcY   = new TH3D(Form("mcY_%s_%s", mcFileName, dirName),
                    title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
  _dataX = new TH3D("dataX", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
  _dataY = new TH3D("dataY", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);

  _maskX = new TH3D("maskX", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
  _maskY = new TH3D("maskY", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);

  _mcXdist   = new TH2D("mcXdist",   "", nx, &binX[0], nz, &binZ[0]);
  _mcYdist   = new TH2D("mcYdist",   "", nx, &binX[0], nz, &binZ[0]);
  _dataXdist = new TH2D("dataXdist", "", nx, &binX[0], nz, &binZ[0]);
  _dataYdist = new TH2D("dataYdist", "", nx, &binX[0], nz, &binZ[0]);
}

void CherenkovCalc::SetInitialParams(const params & mcparams,
                                     const enum DET det)
{
  _Init.Fx = det == ND? mcparams.NDFx: mcparams.FDFx;
  _Init.Fy = det == ND? mcparams.NDFy: mcparams.FDFy;
  _Init.attB = mcparams.attB;
  _Init.inv_att1 = 1/mcparams.att1;
  _Init.inv_att2 = 1/mcparams.att2;
  _Init.Ys = mcparams.Ys;
  _Init.Ec = mcparams.Ec;
  _Init.rngeff = mcparams.rngeff;
}

// Copied from Geo::getPigtail
double CherenkovCalc::getpigtail(const int cell, const bool isX)
{
  int cellid = cell % kCellsPerModule;
  // In vertical planes, high cell numbers have longer fibres.
  // In horizontal planes it's the opposite, so flip it round.
  if(!isX) cellid = kCellsPerModule-cellid-1;
  // This should really never happen, but just to be safe...
  if(cellid < 0 || cellid >= kCellsPerModule) return 100;
  // Email from Tom Chase 2011-04-29
  // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
  // different happens.
  return kPigtails[cellid];
}

// Pure function that returns attenuation
static double _exact_atten(const double attB, const double inv_att1,
                           const double inv_att2, const double dist,
                           const double rngeff, const bool longway)
{
  return (longway?rngeff:1)*
         (attB * exp(-dist*inv_att1) + (1-attB) * exp(-dist*inv_att2));
}

// Calculate the exact attenaution for a given distance 'dist' and the
// current parameters.  Intended only for use to fill the cache.
double CherenkovCalc::exact_atten(const double dist, const bool longway)
{
  return _exact_atten(_Shift.attB, _Shift.inv_att1, _Shift.inv_att2,
                      dist, _Shift.rngeff, longway);
}

// Because calculating attenuation exactly is expensive, make a table
// of attenuation to the nearest 'atten_step' (5cm as I write this).
void CherenkovCalc::init_atten_cache()
{
  atten_cache_short.resize(nattensteps);
  atten_cache_long .resize(nattensteps);
  for(int d = 0; d < nattensteps; d++){
    atten_cache_short[d] = exact_atten(d*atten_step, false);
    atten_cache_long [d] = exact_atten(d*atten_step, true );
  }
}

// Make a TTreeReader given a file and the name of the desired tree.
// Be tolerant of name changes over time and try a few.
//
// Implemented with goto's for fun!
static TTreeReader * makereader(const char * const dir, TFile * file)
{
  TTreeReader * reader = new TTreeReader(Form("cerenkov/%s", dir), file);
  if(!reader->IsInvalid()) goto ok;

  reader->SetTree(Form("cerenkovselection/%s", dir), file);
  if(!reader->IsInvalid()) goto ok;

  reader->SetTree(Form("cerenkovresponse/%s", dir), file);
  if(!reader->IsInvalid()) goto ok;

  if(strcmp(dir, "cerenkovMuonTree")) goto fail;

  reader->SetTree("cerenkovresponse/cerenkovTree", file);
  if(!reader->IsInvalid()) goto ok;

  goto fail;

  ok:

  if(reader->GetEntries(false) == 0){
    fprintf(stderr, "Empty input from %s\n", file->GetName());
    _exit(1);
  }
  return reader;

  fail:

  fprintf(stderr, "Couldn't get tree from %s, %s\n", file->GetName(), dir);
  _exit(1);
}

static inline double clamp(const double v, const double mi, const double ma)
{
  return v < mi? mi: v > ma? ma: v;
}

struct binw{
  double prev, curr, next;
};

#define SMEAR

static binw smearbin(
#ifndef SMEAR
  __attribute__((unused)) const vector<double> & bins,
  __attribute__((unused)) const double val,
  __attribute__((unused)) const int bin
#else
  const vector<double> & bins, const double val, const int bin
#endif
)
{
  binw ans;

  // This is the regular way of binning!
#ifndef SMEAR
  ans.prev = ans.next = 0;
  ans.curr = 1;
  return ans;
#else
  const double binlowedge = bins[bin];
  const double binhigedge = bins[bin+1];
  const double binfrac = clamp((   val     - binlowedge)/
                               (binhigedge - binlowedge), -0.5, 1.5);
  if(binfrac > 0.5){
    ans.prev = 0;
    ans.next = binfrac - 0.5;
    ans.curr = 1 - ans.next;
  }
  else{
    ans.prev = 0.5 - binfrac;
    ans.curr = 1 - ans.prev;
    ans.next = 0;
  }
#endif

#if 0
  printf("%f %d:\t%f %f %f %f \n", val, bin, binfrac, ans.prev, ans.curr, ans.next);
#endif

  return ans;
}

static void Add(TH3D * h, const int bx, const int by, const int bz,
                const double w)
{
  if(w == 0) return;
  const int gbin = h->GetBin(bx, by, bz);
  h->SetBinContent(gbin, h->GetBinContent(gbin) + w);
}

// Get the zero-indexed bin number that 'pe' is in. Optimize by checking
// if it is the same as the previous answer. This works because we
// sorted by the initial PE.
static unsigned int getpebin(const double pe,
                             const std::vector<double> & bins)
{
  static unsigned int prev = 0;
  if(prev < bins.size()-1 && bins[prev] < pe && bins[prev+1] > pe)
    return prev;

  // Don't try to optimize by restricting the search range.  The pe value may
  // have been modified so that it is out of range.
  return prev = std::lower_bound(bins.begin(), bins.end(), pe)
                - bins.begin() - 1;
}

bool CherenkovCalc::fill_no_overflow(TH3D * h, const double x,
                                     const double y, const double z,
                                     const double weight)
{
  const double xlo = distbins[0];
  const double xhi = distbins[nx]; // sic, because nx is 1 less than the size
  const double zlo = wbins[0];
  const double zhi = wbins[nz];

  // Do not check y here; that's handled below
  if(!(x < xhi && x >= xlo && z < zhi && z >= zlo)) return false;

  const int xbin = h->GetXaxis()->FindBin(x);
  const int ybin = getpebin(y, extpebins);
  const int zbin = h->GetZaxis()->FindBin(z);

  // extpebins is indexed like ROOT histograms with overflow: 0 to ny+1
  const binw bw = smearbin(extpebins, y, ybin);

  if(ybin-1 >= 1)                      Add(h,xbin,ybin-1,zbin,bw.prev*weight);
  if(ybin   >= 1 && ybin   <= (int)ny) Add(h,xbin,ybin  ,zbin,bw.curr*weight);
  if(               ybin+1 <= (int)ny) Add(h,xbin,ybin+1,zbin,bw.next*weight);

  return true;
}

bool CherenkovCalc::ok_w(const float w, __attribute__((unused)) const bool isX)
{
  // From looking at the curves at
  // https://nusoft.fnal.gov/nova/calibration/Attenuation_Calibration/SecondAna/Period3/Epoch3c
  // Want to avoid the roll-off at the ends, even if we think we get
  // some of it right.

//#define RESTRICTFDW

  if(isFD){
#ifdef RESTRICTFDW
     if(!isX) return w > fdwmin && w < fdwmax;
     const double cen = 500;
     return abs(w - cen) < 100;
#else
     return w > fdwmin && w < fdwmax;
#endif
  }
  else     return w > ndwmin && w < ndwmax;
}

static bool goodfdplane(const int plane)
{
//#define BYDIBLOCK

#ifdef BYDIBLOCK
  const int diblock = 14;
  return plane >= 64*(diblock-1) &&
         plane <  64*diblock;
#else
  return plane >= 64;
#endif
}

static bool goodfdcell(__attribute__((unused)) const bool isX,
                       __attribute__((unused)) const int cell)
{
  if(isX) return true;

#define ALL

#if defined ALL
  return true;
#elif defined BOTTOMTHIRD
  return cell                   < 32 * 4;
#elif defined MIDDLETHIRD
  return cell >= 32 * 4 && cell < 32 * 8;
#elif defined TOPTHIRD
  return cell >= 32 * 8                 ;
#endif
}

bool CherenkovCalc::accepthit(const bool isFD, const bool isX, const int plane,
                              const int cell,
                              const float pathlength, const float disttostop,
                              const float w, const float PE)
{
  // remove first two blocks of FD, which have lower pseudocumene scintillator
  if(isFD && !goodfdplane(plane)) return false;

  if(isFD && !goodfdcell(isX, cell)) return false;

  if(pathlength <= MINPATHLENGTH) return false;
  if(disttostop > distbins[distbins.size()-1]) return false;
  if(!ok_w(w, isX)) return false;

  // Data and MC disagree on where the threshold is, so remove the
  // disagreeing part. Not sure if this is also true for the ND, but ND
  // hits aren't mostly pasted up against the threshold, so it's not
  // such an obvious problem.
  if(isFD && PE < 30) return false;

  return true;
}

void CherenkovCalc::FillDataHistos()
{
  TTreeReader * reader = makereader(_dirName, _dataFile);

  TTreeReaderValue<float> rvPathLength(*reader, "pathLength");
  TTreeReaderValue<float> rvPE        (*reader, "pe");
  TTreeReaderValue<float> rvDistToStop(*reader, "distFromEnd");
  TTreeReaderValue<char>  rvIsX       (*reader, "isX");
  TTreeReaderValue<float> rvw         (*reader, "w");
  TTreeReaderValue<int>   rvPlane     (*reader, "plane");
  TTreeReaderValue<int>   rvCell      (*reader, "cell");

  _dataX->Reset();
  _dataY->Reset();
  _dataXdist->Reset();
  _dataYdist->Reset();

  printf("%9u %s %s\n", (unsigned)reader->GetEntries(false),
         _dataFile->GetName(), _dirName);
  fflush(stdout);

  unsigned int got[nwbins] = { 0 };
  bool capped[nwbins] = { false };
  unsigned int ncapped = 0;
  while (reader->Next()) {
    // remove first 2 blocks of FD, which have lower pseudocumene scintillator
    if(!accepthit(isFD, *rvIsX, *rvPlane, *rvCell, *rvPathLength, *rvDistToStop, *rvw, *rvPE))
      continue;

    // Zero-indexed bin number
    const unsigned int wBin =
      std::lower_bound(wbins.begin(), wbins.end(), *rvw) - wbins.begin() - 1;

    if(wBin >= wbins.size()) continue; // shouldn't be possible

    if(got[wBin] >= maxDATAhitsPERw){
      if(!capped[wBin]){
        capped[wBin] = true;
        ncapped++;
      }
      if(nwbins == ncapped) break;
      else continue;
    }

    if(fill_no_overflow((*rvIsX)?_dataX:_dataY,
                        *rvDistToStop, (*rvPE)/(*rvPathLength), *rvw))
      got[wBin]++;
  }

  for(unsigned int i = 0; i < nwbins; i++)
    printf("%9u data hits accepted in w bin %d%s\n", got[i],
      i, capped[i]?" (capped)": Form(" (want %.1fx more)",
      (float)maxDATAhitsPERw/got[i]));
  printf("\n");

  delete reader;

  saneProjectionXZ(_dataXdist, _dataX);
  saneProjectionXZ(_dataYdist, _dataY);

  _dataFile->Close();
}

static TH3D * invmask(TH3D * mask)
{
  TH3D * notmask = (TH3D*)mask->Clone("notX");
  for(int i = 1; i <= mask->GetNbinsX(); i++)
    for(int j = 1; j <= mask->GetNbinsY(); j++)
      for(int k = 1; k <= mask->GetNbinsZ(); k++)
        notmask->SetBinContent(i, j, k, !mask->GetBinContent(i, j, k));
  return notmask;
}

TH3D * CherenkovCalc::invmaskX()
{
  return invmask(_maskX);
}

TH3D * CherenkovCalc::invmaskY()
{
  return invmask(_maskY);
}

// Return true if any change was made
static bool fill_mask(TH3D * mask, TH3D * h)
{
  printf("New mask for %s\n", h->GetName());
  bool changed = false;
  for(int k = 1; k <= h->GetNbinsZ(); k++){
    for(int i = 1; i <= h->GetNbinsX(); i++){
      for(int j = 1; j <= h->GetNbinsY(); j++){
        const bool newval = h->GetBinContent(i, j, k) > stat_req;
        const bool oldval = mask->GetBinContent(i, j, k);
        printf("%c", newval?'.':'X');
        if(oldval != newval){
          changed = true;
          mask->SetBinContent(i, j, k, newval);
        }
      }
      puts("");
    }
    puts("");
  }
  return changed;
}


bool CherenkovCalc::MakeMasks()
{
  const bool changedx = fill_mask(_maskX, _mcX);
  const bool changedy = fill_mask(_maskY, _mcY);
  return changedx || changedy;
}

static bool compare_dat(const dat & a, const dat & b)
{
  return a.PE_PathLength < b.PE_PathLength;
}

double CherenkovCalc::peWeight(const bool isX, const dat & d)
{
  const double attenafter = atten_cache_short[d.distshort]
                          + atten_cache_long [d.distlong ];

  const double NphoShift = (isX?_Shift.Fx:_Shift.Fy)*
    (_Shift.Ys*d.ScintPho + _Shift.Ec*d.CherenkovPho);

  return NphoShift * attenafter * d.InvTotPhoAtten;
}

// Speeds up the whole program by 40% compared to naive use of TH2::Fill().
void CherenkovCalc::FastMCFill(TH3D * hist,
                               const std::vector<dat> & dats,
                               const bool isX)
{
  double newhist[nx][ny][nz];
  memset(newhist, 0, sizeof(double)*nx*ny*nz);

  for(const dat & d : dats){
    const double weight = peWeight(isX, d);

    if(weight == 0) continue;

    // Apply a very rough threshold so that low energy hits drop out entirely
    // Of course, this can only kill hits, not revive them.
    const double threshold = 17.;
    const double scaled_pe = d.PE * weight;
    if(scaled_pe < threshold) continue;

    const double scaled_pe_pathLength = d.PE_PathLength * weight;

    // Range checked when we cached the hits: don't check again
    const int xbin = d.DistBin, zbin = d.wBin;

    // Returns zero-indexed answer for extpebins, and subtract one to get
    // zero-indexed answer for pebins.
    const int ybin = getpebin(scaled_pe_pathLength, extpebins) - 1;

    if(!(ybin >= -1 && ybin < (int)ny+1 /* for overflow smear bins */)) continue;

    // Simple smearing over this bin and one adjacent bin so that events don't
    // move abruptly and MINUIT can work better.  Only needed for y (pe/path).
    const binw bw = smearbin(extpebins, scaled_pe_pathLength, ybin+1 /* to index into extpebins */);

    #ifdef USEPROTONWEIGHT
    if(d.proton && PROTONWEIGHT != 1.){
      if(ybin-1>=0)                newhist[xbin][ybin-1][zbin] += bw.prev*PROTONWEIGHT;
      if(ybin  >=0&&ybin  <(int)ny)newhist[xbin][ybin  ][zbin] += bw.curr*PROTONWEIGHT;
      if(           ybin+1<(int)ny)newhist[xbin][ybin+1][zbin] += bw.next*PROTONWEIGHT;
    }
    else
    #endif
    {
      if(ybin-1>=0)                newhist[xbin][ybin-1][zbin] += bw.prev;
      if(ybin  >=0&&ybin  <(int)ny)newhist[xbin][ybin  ][zbin] += bw.curr;
      if(           ybin+1<(int)ny)newhist[xbin][ybin+1][zbin] += bw.next;
    }
  }

  double tot = 0;
  for(unsigned int i = 0; i < nx; i++){
    for(unsigned int j = 0; j < ny; j++){
      for(unsigned int k = 0; k < nz; k++){
        tot += newhist[i][j][k];
        hist->SetBinContent(i+1, j+1, k+1, newhist[i][j][k]);
      }
    }
  }
  hist->SetEntries(tot);
}

void CherenkovCalc::fill_hit_cache_and_init_things()
{
  TTreeReader * reader = makereader(_dirName, _mcFile);

  TTreeReaderValue<float> rvPathLength (*reader, "pathLength");
  TTreeReaderValue<float> rvPE         (*reader, "pe");
  TTreeReaderValue<float> rvDistToStop (*reader, "distFromEnd");
  TTreeReaderValue<char>  rvIsX        (*reader, "isX");
  TTreeReaderValue<float> rvCerenkovPho(*reader, "cerenkovPho");
  TTreeReaderValue<float> rvEBirks     (*reader, "eBirks");
  TTreeReaderValue<float> rvw          (*reader, "w");
  TTreeReaderValue<char>  rvIsProton   (*reader, "isProton");
  TTreeReaderValue<int>   rvPlane      (*reader, "plane");
  TTreeReaderValue<int>   rvCell       (*reader, "cell");

  _mcX->Reset();
  _mcY->Reset();

  printf("%9u %s %s\n", (unsigned)reader->GetEntries(false),
      _mcFile->GetName(), _dirName);

  unsigned int got[nwbins] = { 0 };
  unsigned int gotall = 0, trueprotons = 0, truemuons = 0;
  bool capped[nwbins] = { false };
  unsigned int ncapped = 0;

  printf("Filling hit cache:\n");
  fflush(stdout);
  while (reader->Next()) {
    // remove first two blocks of FD, which have lower pseudocumene scintillator
    if(!accepthit(isFD, *rvIsX, *rvPlane, *rvCell, *rvPathLength, *rvDistToStop, *rvw, *rvPE)) continue;

    if(*rvEBirks <= 0) continue;

    dat d;
    d.PE            =  *rvPE;
    d.PE_PathLength = (*rvPE)/(*rvPathLength);
    d.CherenkovPho  = (*rvCerenkovPho)*_Init.Ec;

    // zero-indexed bins
    d.DistBin = std::lower_bound(distbins.begin(), distbins.end(), *rvDistToStop)
                - distbins.begin() - 1;
    if(d.DistBin < 0 || d.DistBin >= (int)nx) continue;
    d.wBin    = std::lower_bound(wbins.begin(), wbins.end(), *rvw)
                - wbins.begin() - 1;
    if(d.wBin < 0 || d.wBin >= (int)nz) continue;
    d.ScintPho      = (*rvEBirks)*1000*_Init.Ys;
    d.proton = *rvIsProton;

    const double distshort =   cellhalflength - (*rvw) + getpigtail(*rvCell, *rvIsX);
    const double distlong  = 3*cellhalflength + (*rvw) + getpigtail(*rvCell, *rvIsX);

    const double attenbefore = _exact_atten(_Init.attB, _Init.inv_att1, _Init.inv_att2,
                                 distshort, _Init.rngeff, false) +
                    _exact_atten(_Init.attB, _Init.inv_att1, _Init.inv_att2,
                                 distlong , _Init.rngeff, true );

    d.InvTotPhoAtten = 1/((d.ScintPho + d.CherenkovPho)*attenbefore);
    if(d.InvTotPhoAtten <= 0) continue;

    d.distshort = clamp(0, distshort/atten_step + 0.5, nattensteps-1);
    d.distlong  = clamp(0, distlong /atten_step + 0.5, nattensteps-1);

    if(got[d.wBin] >= maxMChitsPERw){
      if(!capped[d.wBin]){
        capped[d.wBin] = true;
        ncapped++;
      }
      if(nwbins == ncapped) break;
      else continue;
    }

    if(*rvIsX) datsx.push_back(d);
    else       datsy.push_back(d);

    if(fill_no_overflow((*rvIsX)?_mcX:_mcY, *rvDistToStop,
                        (*rvPE)/(*rvPathLength), *rvw,
                        #ifdef USEPROTONWEIGHT
                          (*rvIsProton)?PROTONWEIGHT:1)){
                        #else
                          1)){
                        #endif
      got[d.wBin]++;
      gotall++;
      if(*rvIsProton) trueprotons++;
    }

    if(gotall%1000000 == 0)
      printf("Accepted %d hits and counting...\n", gotall);
  }

  printf("Accepted %.0f%% true protons\n", 100.*trueprotons/gotall);
  printf("Accepted %.0f%% true muons\n"  , 100.*truemuons  /gotall);

  delete reader;

  printf("Sorting hits\n");
  // Optimize by sorting by PE/pathlength and then starting each bin search
  // by guessing that the PE/pathlength bin is the same as the previous.
  std::sort(datsx.begin(), datsx.end(), compare_dat);
  std::sort(datsy.begin(), datsy.end(), compare_dat);
  printf("Hit cache ready\n");

  for(unsigned int i = 0; i < nwbins; i++)
    printf("%9u  MC  hits accepted in w bin %d%s\n", got[i],
      i, capped[i]?" (capped)": Form(" (want %.1fx more)",
      (float)maxMChitsPERw/got[i]));
  printf("\n");

  saneProjectionXZ(_mcXdist, _mcX);
  saneProjectionXZ(_mcYdist, _mcY);

  MakeMasks();
}

void CherenkovCalc::FillMCHistos(const detparams & mcparams,
                                 const bool changed_x, const bool changed_y)
{
  _Shift.Fx = mcparams.Fx/_Init.Fx;
  _Shift.Fy = mcparams.Fy/_Init.Fy;
  _Shift.Ys = mcparams.Ys/_Init.Ys;
  _Shift.Ec = mcparams.Ec/_Init.Ec;

  _Shift.attB = mcparams.attB;
  _Shift.inv_att1 = mcparams.inv_att1;
  _Shift.inv_att2 = mcparams.inv_att2;
  _Shift.rngeff = mcparams.rngeff;
  init_atten_cache();

  if(datsx.empty()) fill_hit_cache_and_init_things();

  if(changed_x) FastMCFill(_mcX, datsx, true);
  if(changed_y) FastMCFill(_mcY, datsy, false);

  // Do not need to fill projections again because they are unchanged.
}

static bool complain_about_problems = true;

static double
llike_hist(const TH3D * const hObs,     const TH3D * const hMC,
           const TH2D * const hObsDist, const TH2D * const hMCDist,
           const TH3D * const mask)
{
  double llike = 0;
  int problems = 0;
  for(int iZ = 1; iZ <= hMC->GetNbinsZ(); ++iZ) {
    for(int iX = 1; iX <= hMC->GetNbinsX(); ++iX) {
      const double data_at_dist = hObsDist->GetBinContent(iX, iZ);
      const double   mc_at_dist =  hMCDist->GetBinContent(iX, iZ);
      const double mc_scale = data_at_dist/mc_at_dist;

      for(int iY = 1; iY <= hMC->GetNbinsY(); ++iY) {
        if(mask->GetBinContent(iX, iY, iZ) == 0) continue;

        const double data = hObs->GetBinContent(iX, iY, iZ);
        const double rawmc = hMC->GetBinContent(iX, iY, iZ);

        if(rawmc < stat_warn_level){
          problems++;
          if(rawmc <= 0) continue;
        }

        // normalize MC to match data
        const double mc = rawmc * mc_scale;

        llike += mc - data + (data <= 0?0:data * log (data/mc));
      }
    }
  }
  if(complain_about_problems && problems)
    printf("%d low-stat MC bin%s in %s!\n",
           problems, problems == 1?"":"s", hMC->GetName());

  return llike;
}

double CherenkovCalc::CalcChi2(const detparams & mcparams,
                               const bool changed_x, const bool changed_y)
{
  FillMCHistos(mcparams, changed_x, changed_y);
  return llike_hist(_dataX, _mcX, _dataXdist, _mcXdist, _maskX) +
         llike_hist(_dataY, _mcY, _dataYdist, _mcYdist, _maskY);
}

CherenkovFitter::CherenkovFitter(const params & mcparams)
  : _Init(mcparams)
{
}

void CherenkovFitter::AddFDSample(CherenkovCalc* calc)
{
  calc->SetInitialParams(_Init, FD);
  calc->FillDataHistos();
  _fdSamples.push_back(calc);
}

void CherenkovFitter::AddNDSample(CherenkovCalc* calc)
{
  calc->SetInitialParams(_Init, ND);
  calc->FillDataHistos();
  _ndSamples.push_back(calc);
}

// Mechanism for the user to request plotting in the middle of the run
static CherenkovFitter * thefitter = NULL;
static bool plot = false;

static void intermediate_plot(const int signal)
{
  if(signal != SIGUSR1) return;
  plot = true;
}

// Keep track of the best values independently of MINUIT, because if
// MINUIT doesn't think a fit has converged, it may not leave the
// parameters at the best point seen so far (yes, this certainly
// happens).
static double lowest_llike = 1e300;
static double bestpars[npar] = { 0 };

// The objective function for MINUIT.
static void objfcn(__attribute__((unused)) int & np,
                   __attribute__((unused)) double * gin,
                                           double & llike,
                                           double * par,
                                           int flag)
{
  // taking the absolute value is a nice way to prevent non-physical nonsense
  // without using SET LIMITS.  The fit will only hit negative values when it is
  // going crazy, and this can give it a way to get back home.
  detparams FD, ND;
  static detparams old_FD, old_ND;
  FD.Fx = fabs(par[0]);
  FD.Fy = fabs(par[1]);
  ND.Fx = fabs(par[2]);
  ND.Fy = fabs(par[3]);

  FD.Ys       = ND.Ys       = fabs(par[4]);
  FD.Ec       = ND.Ec       = fabs(par[5]);
  FD.attB     = ND.attB     = fabs(par[6]);
  FD.inv_att1 = ND.inv_att1 = 1/fabs(par[7]);
  FD.inv_att2 = ND.inv_att2 = 1/fabs(par[8] + par[7]);
  FD.rngeff   = ND.rngeff   = fabs(par[9]);

  double ndllike = 0, fdllike = 0;
  static double old_ndllike = 0, old_fdllike = 0;

  const bool fdchanged = old_FD != FD;
  const bool ndchanged = old_ND != ND;

  if(fdchanged){
    const bool fdchanged_x = !old_FD.equalsx(FD);
    const bool fdchanged_y = !old_FD.equalsy(FD);

    for(CherenkovCalc* calc : sfdSamples)
      fdllike += calc->CalcChi2(FD, fdchanged_x, fdchanged_y);
    old_fdllike = fdllike;
  }
  else{
    fdllike = old_fdllike;
  }

  if(ndchanged){
    const bool ndchanged_x = !old_ND.equalsx(ND);
    const bool ndchanged_y = !old_ND.equalsy(ND);

    for(CherenkovCalc* calc : sndSamples)
      ndllike += calc->CalcChi2(ND, ndchanged_x, ndchanged_y);
    old_ndllike = ndllike;
  }
  else{
    ndllike = old_ndllike;
  }

  old_FD = FD;
  old_ND = ND;

  llike = ndllike + fdllike;

  {
    if(llike < lowest_llike){
      lowest_llike = llike;
      printf("lowest llike so far: %.9f%s\n", llike,
        lowest_llike - llike < epsilonll?" (barely)":"");
      lowest_llike = llike;
      memcpy(bestpars, par, sizeof(double)*npar);
    }
    printf(".");
    if(flag == 3) printf("\n");
    fflush(stdout);
  }
  if(plot){
    static int setn = 0;
    thefitter->Plot(Form("intermediate%d", setn));
    thefitter->Plot("intermediatelatest"); // inefficient, but ok...
    setn++;
    plot = false;
  }
}

// Get the zero-indexed value 'i'.
static double getpar(TMinuit & mn, const int i)
{
  double val, err;
  mn.GetParameter(i, val, err);
  return val;
}

// Get the zero-indexed error 'i'.
double geterr(TMinuit & mn, const int i)
{
  double val, err;
  mn.GetParameter(i, val, err);
  return err;
}

// Check if the value of the objective function at the current
// parameters is less than, equal to, or greater than the last time this
// function was called. If it greater, then reset the parameters to
// where they were before. This is useful because sometimes MIGRAD or
// SIMPLEX(?) get stuck and leave things worse than before.
void mincheck(TMinuit & mn, const bool reset = false)
{
  static double oldllike = 0;

  if(oldllike == 0 || reset){
    printf("\nllike = %8.1f\n", lowest_llike);
    lowest_llike = 1e300;
    memset(bestpars, 0, sizeof(double)*npar);
  }
  else if(oldllike < lowest_llike){
    const int prevprint = mn.fISW[4];
    mn.Command("SET PRINT -1");
    if(prevprint >= 0)
      printf("\nllike = %8.1f, %8.5f MORE than before (%8.1f). Resetting.\n",
             lowest_llike, lowest_llike - oldllike, oldllike);
    for(unsigned int i = 1; i <= npar; i++)
      mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
    mn.Command(Form("SET PRINT %d", prevprint));
    return;
  }
  else if(lowest_llike == 0){
    printf("\nllike = 0 which is nonsense. Resetting.\n");
    for(unsigned int i = 1; i <= npar; i++)
      mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
    return;
  }
  else if(oldllike == lowest_llike){
    printf("\nllike = %8.1f, which is the same as last time\n", lowest_llike);
  }
  else if(oldllike - lowest_llike > epsilonll){
    printf("\nllike = %8.1f, which is %f less than last time\n",
           lowest_llike, oldllike - lowest_llike);
  }
  else{
    printf("\nllike = %8.1f, which is barely (%g) less than last time\n",
           lowest_llike, oldllike - lowest_llike);
  }

  if(lowest_llike < mn.fAmin){
    const int prevprint = mn.fISW[4];
    mn.Command("SET PRINT -1");
    if(prevprint >= 0)
      printf("\nNote: best = %f < fAmin = %f, using best pars instead\n",
             lowest_llike, mn.fAmin);
    for(unsigned int i = 1; i <= npar; i++)
      mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
    mn.Command(Form("SET PRINT %d", prevprint));
  }

  oldllike = mn.fAmin;
  for(unsigned int i = 0; i < npar; i++) bestpars[i] = getpar(mn, i);
}

void CherenkovFitter::Plot(const char * const tag)
{
  for(CherenkovCalc * sample : _fdSamples){
    TwoDplots(*sample, Form("%s_%sfit", sample->plot_name, tag));
    Flatplots(*sample, Form("%s_3fitflat", sample->plot_name), tag);
  }
  for(CherenkovCalc * sample : _ndSamples){
    TwoDplots(*sample, Form("%s_%sfit", sample->plot_name, tag));
    Flatplots(*sample, Form("%s_3fitflat", sample->plot_name), tag);
  }
}

__attribute__((unused)) static void Console(TMinuit & mn)
{
  std::string in;
  printf("MINIUT> ");
  while(getline(std::cin, in)){
    if(in == "quit" || in == "exit" || in == ".q") break;
    mn.Command(in.c_str());
    printf("MINIUT> ");
  }
}

static void setarand(TMinuit & mn, const int p, const double min,
                     const double max)
{
  const double v = min + (max-min)*R.Rndm();
  mn.Command(Form("SET PAR %d %f", p, v));
}

#ifndef FULL
__attribute__((unused))
#endif
static void randomparameters(TMinuit & mn)
{
  if(float_fdfx)   setarand(mn, 1, 0.3, 0.6);
  if(float_fdfy)   setarand(mn, 2, 0.3, 0.6);
  if(float_ndfx)   setarand(mn, 3, 0.3, 0.6);
  if(float_ndfy)   setarand(mn, 4, 0.3, 0.6);

  if(float_ec)     setarand(mn, 6, 0.7, 0.9);
  if(float_attb)   setarand(mn, 7, 0.1, 0.9);
  if(float_att1)   setarand(mn, 8, 110, 400);
  if(float_att2)   setarand(mn, 9, 110, 900);
  if(float_rngeff) setarand(mn, 10, 0.98, 1.0);
}

void CherenkovFitter::Fit(const params & mcparams)
{
  sfdSamples = _fdSamples; // So the MINUIT FCN can get at them
  sndSamples = _ndSamples;

  signal(SIGUSR1, intermediate_plot);
  thefitter = this; // so the signal handler can get at it.  Oh boy.

  TMinuit mn(npar);
  mn.SetFCN(objfcn);

  mn.Command("SET ERR 0.5"); // Since we're using log likelihoods

  // Make SCAN plots smaller than default
  mn.Command("SET WIDTH 80");
  mn.Command("SET LINES 40");

  int ierr;
  mn.mnparm( 1 - 1, "FDFx", mcparams.FDFx, 0.01, 0.2, 2, ierr);
  mn.mnparm( 2 - 1, "FDFy", mcparams.FDFy, 0.01, 0.2, 2, ierr);
  mn.mnparm( 3 - 1, "NDFx", mcparams.NDFx, 0.01, 0, 2, ierr);
  mn.mnparm( 4 - 1, "NDFy", mcparams.NDFy, 0.01, 0, 2, ierr);
  mn.mnparm( 5 - 1, "Ys",   mcparams.Ys  ,   10, 0, 0, ierr);
  mn.mnparm( 6 - 1, "Ec",   mcparams.Ec  ,  0.1, 0, 0, ierr);
  mn.mnparm( 7 - 1, "attB", mcparams.attB, 0.01, 0, 1, ierr);
  mn.mnparm( 8 - 1, "att1", mcparams.att1, 10, 100, 2000, ierr);

  // Define this parameter to be the difference between the two attenuation
  // parameters.  This is a handy way to keep them from becoming degenerate
  // and spoiling the fit
  mn.mnparm( 9 - 1, "Da12", mcparams.att2 - mcparams.att1, 10, 100, 2000, ierr);
  mn.mnparm(npar-1, "rngeff", mcparams.rngeff, 0.01, 0, 1.5, ierr);
  if(ierr > 0){
    fprintf(stderr, "Failed to set up parameters\n");
    _exit(1);
  }

  mn.Command("SET STRATEGY 2");

  // Reduce MINUIT's expectations of accuracy because the objective function is
  // not going to be very smooth since it involves events migrating over bin
  // boundaries.  Experimentally, near a minimum, parameter changes of epsilon
  // result in FCN value changes of ~1E-6 of the FCN value.
  mn.Command("SET EPS 1E-6");

  mn.fGraphicsMode = false;  // So SCAN prints a plot

  // Only three of {x, y, Ys} should be free.  Which ones is arbitrary. It's
  // easiest to compare to previous results if we fix Ys.
  mn.Command("SET PAR 5 3151.04");
  mn.Command("FIX 5");

  if(!float_fdfx) mn.Command("FIX 1");
  if(!float_fdfy) mn.Command("FIX 2");
  if(!float_ndfx) mn.Command("FIX 3");
  if(!float_ndfy) mn.Command("FIX 4");

  if(!float_ec)   mn.Command("FIX 6");
  if(!float_attb) mn.Command("FIX 7");
  if(!float_att1) mn.Command("FIX 8");
  if(!float_att2) mn.Command("FIX 9");
  if(!float_rngeff) mn.Command("FIX 10");

#ifdef FULL
  mn.Command("CALL");
  for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
  for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();
  mn.Command("CALL");

  Plot("1ralsopre");

  // In case input parameters are crazy, do quick random search
  printf("WHAT WE ARE DOING: Trying some random parameters\n");

  // Start in random places instead of using any prior information
  #define NOHELP

  #ifndef NOHELP
    mincheck(mn, true);
  #endif

  const int prevprint = mn.fISW[4];
  mn.Command("SET PRINT -1");
  complain_about_problems = false;
  for(int i = 0; i < 1000; i++){
    randomparameters(mn);
    mn.Command("CALL");
    mincheck(mn,
      #ifdef NOHELP
        i == 0
      #else
        false
      #endif
    );
  }
  complain_about_problems = true;
  mn.Command(Form("SET PRINT %d", prevprint));

  for(int cornot = 0; cornot < 2; cornot++){
    if(float_ec){
      if(cornot == 0) mn.Command("FIX 6");
      if(cornot == 1) mn.Command("RELEASE 6");
    }
    if(cornot == 1 && !float_ec) break;

    printf("WHAT WE ARE DOING: Set masks at current value\n");
    mn.Command("CALL");
    for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
    for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();
    mn.Command("CALL");

    #if 0
    Console(mn);
    #endif

    mincheck(mn, true);

    printf("WHAT WE ARE DOING: Do a little SIMPLEX to get in a nice spot\n");
    mn.Command("SIMPLEX 100");
    mincheck(mn);

    printf("WHAT WE ARE DOING: See if MIGRAD helps\n");
    mn.Command("MIGRAD");
    mincheck(mn);

    for(int i = 0; i < 2; i++){
      printf("WHAT WE ARE DOING: Reset masks at current value and minimize again\n");
      mn.Command("CALL");
      bool changed = false;
      for (CherenkovCalc* calc : sfdSamples) changed |= calc->MakeMasks();
      for (CherenkovCalc* calc : sndSamples) changed |= calc->MakeMasks();
      mn.Command("CALL");

      if(!changed){
        printf("No change in masks, so no need to minimize again\n");
        break;
      }

      mincheck(mn, true);

      mn.Command("SIMPLEX");
      mincheck(mn);
      mn.Command("SEEK 300 8");
      mincheck(mn);
      mn.Command("MIGRAD");
      mincheck(mn);
      mn.Command("IMPROVE");
      mincheck(mn);
    }

    while(true){
      {
        const double oldmin = mn.fAmin;

        printf("WHAT WE ARE DOING: Scan all to make sure things are sensible\n");
        const double r = 0.02, lo = 1-r, hi = 1+r;
        if(float_fdfx)  mn.Command(Form("SCAN  1 81 %f %f",getpar(mn,0)*lo,getpar(mn,0)*hi));
        if(float_fdfy)  mn.Command(Form("SCAN  2 81 %f %f",getpar(mn,1)*lo,getpar(mn,1)*hi));
        if(float_ndfx)  mn.Command(Form("SCAN  3 81 %f %f",getpar(mn,2)*lo,getpar(mn,2)*hi));
        if(float_ndfy)  mn.Command(Form("SCAN  4 81 %f %f",getpar(mn,3)*lo,getpar(mn,3)*hi));

        if(float_ec)    mn.Command(Form("SCAN  6 81 %f %f",getpar(mn,5)*lo,getpar(mn,5)*hi));
        if(float_attb)  mn.Command(Form("SCAN  7 81 %f %f",getpar(mn,6)*lo,getpar(mn,6)*hi));
        if(float_att1)  mn.Command(Form("SCAN  8 81 %f %f",getpar(mn,7)*lo,getpar(mn,7)*hi));
        if(float_att2)  mn.Command(Form("SCAN  9 81 %f %f",getpar(mn,8)*lo,getpar(mn,8)*hi));
        if(float_rngeff)mn.Command(Form("SCAN 10 81 %f %f",getpar(mn,9)*lo,getpar(mn,9)*hi));

        const double newmin = mn.fAmin;
        const bool significantreduction = oldmin - newmin > epsilonll;
        if(newmin < oldmin)
          printf("SCANs reduced llike by %g%s\n", oldmin - newmin,
                 significantreduction?"":" -- not significant");
        if(!significantreduction) break;
      }

      {
        printf("WHAT WE ARE DOING: found new min, so SIMPLEX, SEEK, MIGRAD, IMPROVE\n");
        const double oldmin = mn.fAmin;
        mn.Command("SIMPLEX");
        mincheck(mn);
        mn.Command("SEEK 300 8");
        mincheck(mn);
        mn.Command("MIGRAD");
        mincheck(mn);
        mn.Command("IMPROVE");
        mincheck(mn);
        const double newmin = mn.fAmin;
        const bool significantreduction = oldmin - newmin > epsilonll;
        if(newmin < oldmin)
          printf("SIMPLEX, etc. reduced llike by %g%s\n", oldmin - newmin,
                 significantreduction?"":" -- not significant");
        if(!significantreduction) break;
      }
    }
  }

  printf("WHAT WE ARE DOING: Call once more at min to make plots right\n");
  mn.Command("CALL");
#else
  mn.Command("SET PAR 1  6.41795e-01");
  mn.Command("SET PAR 2  5.53340e-01");
  mn.Command("SET PAR 3  5.09120e-01");
  mn.Command("SET PAR 4  4.94346e-01");
  mn.Command("SET PAR 5  3.15104e+03");
  mn.Command("SET PAR 6  8.10680e-01");
  mn.Command("SET PAR 7  4.46867e-01");
  mn.Command("SET PAR 8  2.03410e+02");
  mn.Command("SET PAR 9  5.48778e+02");
  mn.Command("SET PAR 10 1.00000e+00 ");

  const double besta1 = getpar(mn, 8-1), bestda = getpar(mn, 9-1);

  printf("WHAT WE ARE DOING: Reset masks at current value\n");
  mn.Command("CALL");
  for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
  for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();

  // Helps?
  mincheck(mn, true);
  mn.Command("SIMPLEX 30");
  mincheck(mn);
#if 0
  mn.Command("MIGRAD");
  mincheck(mn);
#endif

  mn.Command("FIX 8 9");

  mn.Command("FIX 5"); // Fix Ys as in the regular fit
  mn.Command("FIX 10"); // not messing with this at the moment
  mn.Command("SET ERR 5"); // just to cut off SIMPLEX earlier without
                           // changing the number of calls

  const int nside = 28;
  const double a1step = 0.003, dastep = 0.005/4;
  printf("contcenter2: %f %f %d %f %f\n",
         besta1, bestda, nside, a1step, dastep);
  for(int i = nside; i >= -nside; i--){
    for(int j = -nside; j <= nside; j++){
      const double a1rat = 1 + i*a1step;
      const double darat = 1 + j*dastep;
      const double a1 = besta1*a1rat;
      const double da = bestda*darat;
      mn.Command("SET PAR 1  6.41795e-01");
      mn.Command("SET PAR 2  5.53340e-01");
      mn.Command("SET PAR 3  5.09120e-01");
      mn.Command("SET PAR 4  4.94346e-01");
      mn.Command("SET PAR 5  3.15104e+03");
      mn.Command("SET PAR 6  8.10680e-01");
      mn.Command("SET PAR 7  4.46867e-01");

      mn.Command(Form("SET PAR 8 %f", a1));
      mn.Command(Form("SET PAR 9 %f", da));

      mn.Command("FIX 1 2 3 4 5 6 7 8 9");

      // Have to do this every time, because sometimes MINUIT
      // decides 0 isn't good enough, and switches to 1, and it
      // is a permanent switch!
      mn.Command("SET STRATEGY 0");
      mn.Command("CALL");
#if 0
      printf("cont: %10f %10f %10f\n", a1, da, mn.fAmin);
#endif
      mincheck(mn, true);
#if 1
      mn.Command("RELEASE 7");
      mn.Command("SIMPLEX 10");
      mn.Command("MIGRAD 100");
      mincheck(mn);
#endif
#if 1
      printf("cont2: %10f %10f %10f\n", a1, da, mn.fAmin);
#endif
#if 0
      mn.Command("RELEASE 1 2 3 4");
      mn.Command("SIMPLEX 10");
      //mn.Command("MIGRAD 50");
      mincheck(mn);

      mn.Command("RELEASE 6");
      mn.Command("SIMPLEX 10");
      mn.Command("MIGRAD 50");
      mincheck(mn);

      printf("cont3: %10f %10f %10f\n", a1, da, lowest_llike);
#endif
    }
  }

  Console(mn);

  #if 0
  mincheck(mn, true);
  mn.Command("SIMPLEX");

  printf("WHAT WE ARE DOING: Scan all to make sure things are sensible\n");
  mn.Command(Form("SCAN  1 21 %f %f", getpar(mn, 0)*0.98, getpar(mn, 0)*1.02));
  mn.Command(Form("SCAN  2 21 %f %f", getpar(mn, 1)*0.98, getpar(mn, 1)*1.02));
  mn.Command(Form("SCAN  3 21 %f %f", getpar(mn, 2)*0.98, getpar(mn, 2)*1.02));
  mn.Command(Form("SCAN  4 21 %f %f", getpar(mn, 3)*0.98, getpar(mn, 3)*1.02));
  mn.Command(Form("SCAN  6 21 %f %f", getpar(mn, 5)*0.98, getpar(mn, 5)*1.02));
  mn.Command(Form("SCAN  7 21 %f %f", getpar(mn, 6)*0.98, getpar(mn, 6)*1.02));
  mn.Command(Form("SCAN  8 21 %f %f", getpar(mn, 7)*0.98, getpar(mn, 7)*1.02));
  mn.Command(Form("SCAN  9 21 %f %f", getpar(mn, 8)*0.98, getpar(mn, 8)*1.02));
  mn.Command(Form("SCAN 10 21 %f %f", getpar(mn, 9)*0.98, getpar(mn, 9)*1.02));

  printf("WHAT WE ARE DOING: Call once more at min to make plots right\n");
  mn.Command("CALL");
  #endif
#endif

  printf("\n\nFD Fx: %.4f   Fy: %.4f\n"
         "ND Fx: %.4f   Fy: %.4f\n"
         "Ys: %.2f (by definition)\n"
         "Ec: %.4f\n"
         "attB: %6.4f  att1: %6.1f att2: %6.1f\n",
         getpar(mn, 0), getpar(mn, 1), getpar(mn, 2), getpar(mn, 3),
         getpar(mn, 4), getpar(mn, 5), getpar(mn, 6), getpar(mn, 7),
         getpar(mn, 7) + getpar(mn, 8));

  Plot("2post");
}