#define rcpicker_cxx
// This is a utility class for interactive viewing of rocking curve fit
// results produced by the rcfitter TSelector.
//
// To use this macro, try the following session from within the directory
// containing your fit results root file, entitled "pristine_007_results.root"
// in this example.  This macro is located in the parent dir.
//
//  Root > TFile *_file0 = TFile::Open("pristine_007_results.root");
//  Root > .L rcfitter.C++
//  Root > .L rcpicker.C++
//  Root > TString rootfile("pristine_007_rocking_curves.root");
//  Root > TString logfile("pristine/pristine");
//  Root > rcpicker rcp("pristine",7,rootfile,logfile);
//  Root > _file0->cd();
//  Root > rcp.PickOn(hmu);
//

//#define USING_FLICAM_CAMERA 1
//#define USING_ANDOR_CAMERA_2012 1
//#define USING_ANDOR_CAMERA_2013 1
#define USING_ANDOR_CAMERA_2014 1
//#define USING_ANDOR_CAMERA_2015 1
//#define USING_ANDOR_CAMERA_2016 1
//#define USING_HAMA_CAMERA_2016 1

#define USING_SI331_MONO 1
//#define USING_SI220_MONO_REVERSE 1
//#define USING_SI220_MONO_FORWARD 1

#define SLOW_STARTUP_IN_PICKON 1

//#define SCAN_COMBINE_STEPS 4

#include "rcpicker.h"
#include "rcfitter.h"

#include <iostream>
#include <fstream>

#include <TROOT.h>
#include <TStyle.h>
#include <TObjString.h>
#include <TMath.h>
#include <TFile.h>
#include <TF1.h>

Double_t rcpicker::fAngleScaleFactor;
Double_t rcpicker::fDistanceScaleFactor;
Double_t rcpicker::fPixelAspectRatio;
Double_t rcpicker::fDispersionSlope;
TCanvas *rcpicker::fCanvas2;
TCanvas *rcpicker::fCanvas1;
TFile *rcpicker::fFile;
TTree *rcpicker::fTree;
TH1D *rcpicker::fHist;
TH1D* rcpicker::fShow;

rcpicker::rcpicker(const char *prefix, Int_t seq, const char *rootfile,
                                                  const char *elogfile)
 : fPrefix(prefix), fSeq(seq),
   fRootfile(rootfile), fElogfile(elogfile)
{
   fFile = TFile::Open(fRootfile,"READ");
   if (!fFile) {
      std::cerr << "Error - unable to open rctree file" << std::endl;
      return;
   }
   fTree = (TTree *)fFile->Get("rctree");
   if (!fTree) {
      std::cerr << "Error - unable to find rctree in input file" << std::endl;
      return;
   }
   fHist = new TH1D();
   fShow = new TH1D();
   fTree->SetBranchAddress("rchist", &fHist);
   fTree->GetEntry(0);
   fDistanceScaleFactor = 0;
   fPixelAspectRatio = 1;
   fAngleScaleFactor = 0;
}

rcpicker::~rcpicker()
{
   delete fFile;   // this auto-deletes the remaining objects below 
   // delete fTree;
   // delete fHist;
   // delete fShow;
}

Int_t rcpicker::PickOn(TH2D *rcimage)
{
   GetMicroradPerStep();
   TH2D *hbase  = (TH2D*)gDirectory->Get("hbase");
   TH2D *hamp   = (TH2D*)gDirectory->Get("hamp");
   TH2D *hmu    = (TH2D*)gDirectory->Get("hmu");
   TH2D *hsigma = (TH2D*)gDirectory->Get("hsigma");
   TH2D *hmean  = (TH2D*)gDirectory->Get("hmean");
   TH2D *hrms   = (TH2D*)gDirectory->Get("hrms");
   TH2D *hmax   = (TH2D*)gDirectory->Get("hmax");
   TH2D *hpeak  = (TH2D*)gDirectory->Get("hpeak");
   TH1D *rcurve  = (TH1D*)gDirectory->Get("rcurve");
   if (rcimage->GetXaxis()->GetXmax() == rcimage->GetNbinsX()) {
      Normalize(hbase, hamp, hmu, hsigma, hmean, hrms, hmax, hpeak);
   }
   else {
      fShow = rcurve;
   }
   Draw(rcimage);
   fCanvas2->DeleteExec("rcPointPicker");
   fCanvas2->AddExec("rcPointPicker","rcpicker::PointPicker();");
   fCanvas1 = new TCanvas("c1","single-pixel rocking curve",700,10,600,500);
   fCanvas1->SetTopMargin(0.12);
   fCanvas1->SetLeftMargin(0.15);
   fShow->GetXaxis()->SetTitle("theta (#mu^{}rad)");
   fShow->GetYaxis()->SetTitle("intensity (arb. units)");
   fShow->GetYaxis()->SetTitleOffset(1.6);
   Int_t nbins = fShow->GetNbinsX();
   fShow->GetXaxis()->SetLimits(0,nbins*fAngleScaleFactor);
   PointFit(fShow);
   fShow->Draw();
   fCanvas1->Update();
   std::cout << "rcPointPicker initialized - move the cursor over the image"
             << std::endl;
   return 0;
}

void rcpicker::PointPicker()
{
   TObject *select = fCanvas2->GetSelected();
   if (!select) {
      return;
   }
   if (!select->InheritsFrom("TH2")) {
      return;
   }
   if (fCanvas2->GetEvent() == 1) {
      fCanvas2->DeleteExec("rcPointPicker");
      return;
   }
   int px = fCanvas2->GetEventX();
   int py = fCanvas2->GetEventY();
   Double_t x = fCanvas2->AbsPixeltoX(px);
   Double_t y = fCanvas2->AbsPixeltoY(py);
   Int_t ix = ((TH2D *)select)->GetXaxis()->FindBin(x);
   Int_t iy = ((TH2D *)select)->GetYaxis()->FindBin(y);
   Int_t width = ((TH2D *)select)->GetNbinsX();
   fTree->GetEntry(width*iy+ix);
   fCanvas1->cd();
   delete fShow;
   fShow = (TH1D *)fHist->Clone("hshow");
   Int_t nbins = fShow->GetNbinsX();
   fShow->GetXaxis()->SetLimits(0,nbins*fAngleScaleFactor);
   PointFit(fShow);
   fShow->GetXaxis()->SetTitle("theta (#mu^{}rad)");
   fShow->GetYaxis()->SetTitle("intensity (arb. units)");
   fShow->GetYaxis()->SetTitleOffset(1.6);
   fShow->Draw();
   fCanvas1->Update();
   fCanvas2->cd();
}

void rcpicker::PointFit(TH1D *rchist)
{
   rchist->ResetStats();
   if (rchist->GetEntries() > 0) {
      rcfitter f;
#if SCAN_COMBINE_STEPS
      rchist->Rebin(SCAN_COMBINE_STEPS);
#endif
      TFitResultPtr result = f.FitRC(rchist);
      rchist->GetFunction("rcshape")->ResetBit(1<<9);
   }
}

Double_t rcpicker::GetMicronsPerPixel()
{
#if USING_FLICAM_CAMERA
   return fDistanceScaleFactor = 6000./1024;
#elif USING_ANDOR_CAMERA_2012
   return fDistanceScaleFactor = 2.132; // +/-0.02, November 2012 zoom
#elif USING_ANDOR_CAMERA_2013
   return fDistanceScaleFactor = 3.250; // +/-0.02, May, 2013 zoom
#elif USING_ANDOR_CAMERA_2014
   return fDistanceScaleFactor = 3.50; // +/-0.02, May, 2014 zoom
#elif USING_ANDOR_CAMERA_2015
   return fDistanceScaleFactor = 4.08; // +/-0.02, February, 2015 zoom
#elif USING_ANDOR_CAMERA_2016
   return fDistanceScaleFactor = 4.30; // +/-0.02, February, 2016 zoom
#elif USING_HAMA_CAMERA_2016
   return fDistanceScaleFactor = 9.00; // +/-0.02, August, 2016 zoom
#else
   return 999;
#endif
}

Double_t rcpicker::GetPixelAspectRatio()
{
#if USING_FLICAM_CAMERA
   return 1;
#elif USING_ANDOR_CAMERA_2012
   return 1;
#elif USING_ANDOR_CAMERA_2013
   return 1;
#elif USING_ANDOR_CAMERA_2014
   return 1;
#elif USING_ANDOR_CAMERA_2015
   return 1;
#elif USING_ANDOR_CAMERA_2016
   return 1;
#elif USING_HAMA_CAMERA_2016
   return 1.06;    // +/- 1%, include diffraction half-angle projection factor
#else
   return 999;
#endif
}
Double_t rcpicker::GetMicroradPerStep()
{
   if (fAngleScaleFactor) {
      return fAngleScaleFactor;
   }
   std::ifstream scanlog(fElogfile);
   if (scanlog.is_open()) {
      TString line;
      TString pat;
      pat.Form("#S %d",fSeq);
      while (line.ReadLine(scanlog,true)) {
         if (line.BeginsWith(pat)) {
            TObjArray *arg = line.Tokenize(" ");
            TObjString *begin = (TObjString *)(*arg)[4];
            TObjString *end = (TObjString *)(*arg)[5];
            TObjString *steps = (TObjString *)(*arg)[6];
            Double_t begin_deg = begin->GetString().Atof();
            Double_t end_deg = end->GetString().Atof();
            Double_t step_count = steps->GetString().Atoi();
            fAngleScaleFactor = (end_deg-begin_deg)/step_count
                               *(TMath::Pi()/180)*1e6;
            delete arg;
         }
      }
      if (fAngleScaleFactor == 0) {
         std::cout << "looked for but did not find line with pattern /"
                   << pat << "/" << std::endl;
      }
      return fAngleScaleFactor;
   }
 
   // look for a plain text log file that contains this information
   std::ifstream scanlog2(fElogfile + ".info");
   if (scanlog2.is_open()) {
      TString line;
      TString pat;
      pat.Form("scan %d",fSeq);
      while (line.ReadLine(scanlog2,true)) {
         if (line.BeginsWith(pat)) {
            TObjArray *arg = line.Tokenize(" ");
            TObjString *begin = (TObjString *)(*arg)[3];
            TObjString *end = (TObjString *)(*arg)[6];
            TObjString *steps = (TObjString *)(*arg)[9];
            TObjString *delta = (TObjString *)(*arg)[12];
            Double_t begin_deg = begin->GetString().Atof();
            Double_t end_deg = end->GetString().Atof();
            Double_t step_count = steps->GetString().Atoi();
            Double_t delta_deg = delta->GetString().Atof();
            fAngleScaleFactor = (end_deg-begin_deg)/step_count
                               *(TMath::Pi()/180)*1e6;
            fAngleScaleFactor = delta_deg * (TMath::Pi()/180)*1e6;
            delete arg;
         }
      }
      if (fAngleScaleFactor == 0) {
         std::cout << "looked for but did not find line with pattern /"
                   << pat << "/" << std::endl;
      }
      return fAngleScaleFactor;
   }

   std::cerr << "steptomicrorad error - cannot open file " << fElogfile
             << std::endl;
   return 0;
}

Double_t rcpicker::GetDispersionSlope()
{
#if USING_SI331_MONO
   return fDispersionSlope = -0.5;  // (urad/mm)
#elif USING_SI220_MONO_REVERSE
   return fDispersionSlope = -92.;  // (urad/mm) CLS 2016, up-bounce
#elif USING_SI220_MONO_FORWARD
   return fDispersionSlope = 20.5;  // (urad/mm) CLS 2019, down-bounce
//   return fDispersionSlope = 8.2;  // (urad/mm) CLS 2017, down-bounce
//   return fDispersionSlope = 21.0;  // (urad/mm) CLS 2016, down-bounce
//  return fDispersionSlope = 30.0;  // (urad/mm) CLS 2016, down-bounce
#else
   return 999;
#endif
}

Int_t rcpicker::RemoveDispersion(TH2D *hmu, Double_t dthetady)
{
   double dispersion_slope = GetDispersionSlope(); // in urad/mm
   if (dthetady > 0)
      dispersion_slope = dthetady;
   Int_t nbinx = hmu->GetNbinsX();
   Int_t nbiny = hmu->GetNbinsY();
   for (Int_t iy=1; iy <= nbiny; ++iy) {
      double ydisp = hmu->GetYaxis()->GetBinCenter(iy) -
                     hmu->GetYaxis()->GetBinCenter(nbiny/2);
      double dmu = dispersion_slope * ydisp;
      for (Int_t ix=1; ix <= nbinx; ++ix) {
         Double_t mu = hmu->GetBinContent(ix,iy);
         if (mu > 0) { 
            // This correction is + because images are inverted in y
            hmu->SetBinContent(ix,iy,mu + dmu);
         }
      }
   }
   return 0;
}

Int_t rcpicker::Normalize(TH2D *hbase, TH2D *hamp, TH2D *hmu, TH2D *hsigma,
                          TH2D *hmean, TH2D *hrms, TH2D *hmax, TH2D *hpeak)
{
   std::cout << "Normalizing..." << std::flush;
   if (!hbase) {
      hbase = (TH2D *)gROOT->FindObject("hbase");
   }
   if (!hbase) {
      std::cerr << "rcpicker::Normalize error - hbase histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hamp) {
      hamp = (TH2D *)gROOT->FindObject("hamp");
   }
   if (!hamp) {
      std::cerr << "rcpicker::Normalize error - hamp histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hmu) {
      hmu = (TH2D *)gROOT->FindObject("hmu");
   }
   if (!hmu) {
      std::cerr << "rcpicker::Normalize error - hmu histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hsigma) {
      hsigma = (TH2D *)gROOT->FindObject("hsigma");
   }
   if (!hsigma) {
      std::cerr << "rcpicker::Normalize error - hsigma histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hmean) {
      hmean = (TH2D *)gROOT->FindObject("hmean");
   }
   if (!hmean) {
      std::cerr << "rcpicker::Normalize error - hmean histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hrms) {
      hrms = (TH2D *)gROOT->FindObject("hrms");
   }
   if (!hrms) {
      std::cerr << "rcpicker::Normalize error - hrms histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hmax) {
      hmax = (TH2D *)gROOT->FindObject("hmax");
   }
   if (!hmax) {
      std::cerr << "rcpicker::Normalize error - hmax histogram not found" 
                << std::endl;
      return 1;
   }
   if (!hpeak) {
      hpeak = (TH2D *)gROOT->FindObject("hpeak");
   }
   if (!hpeak) {
      std::cerr << "rcpicker::Normalize error - hpeak histogram not found" 
                << std::endl;
      return 1;
   }

   delete fShow;
   fTree->GetEntry(0);
   fShow = (TH1D *)fHist->Clone("hshow");
   TString title;
   title.Form("whole crystal rocking curve of sample %s scan %d",
              fPrefix.Data(),fSeq);
   fShow->SetNameTitle("rcurve",title);
   fShow->Reset();

   Double_t angleScale;
   angleScale = GetMicroradPerStep();
   Double_t distanceScale;
   distanceScale = GetMicronsPerPixel()/1000;
   Double_t uvratio = GetPixelAspectRatio();
   Int_t nbinx = hbase->GetNbinsX();
   Int_t nbiny = hbase->GetNbinsY();
   Int_t rescale_needed = (hbase->GetXaxis()->GetXmax() == nbinx);
   if (rescale_needed) {
      hbase->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hbase->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hamp->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hamp->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hmu->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hmu->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hsigma->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hsigma->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hmean->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hmean->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hrms->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hrms->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hmax->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hmax->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
      hpeak->GetXaxis()->SetLimits(0,nbinx*distanceScale);
      hpeak->GetYaxis()->SetLimits(0,nbiny*distanceScale*uvratio);
   }

#if SLOW_STARTUP_IN_PICKON
   double dispersion_slope = GetDispersionSlope();
   for (Int_t iy=1,e=0; iy <= nbiny; ++iy) {
      double ydisp = hmu->GetYaxis()->GetBinCenter(iy) -
                     hmu->GetYaxis()->GetBinCenter(nbiny/2);
      double dmu = dispersion_slope * ydisp;
      for (Int_t ix=1; ix <= nbinx; ++ix,++e) {
         fTree->GetEntry(e);
         for (Int_t it=1; it <= fShow->GetNbinsX(); ++it) {
            double theta = fShow->GetXaxis()->GetBinCenter(it);
            double sum = fShow->GetBinContent(it);
            sum += fHist->Interpolate(theta - dmu);
            fShow->SetBinContent(it, sum);
         }
         if (rescale_needed) {
            Int_t steps = fHist->GetNbinsX()-1;
            Double_t amp = hamp->GetBinContent(ix,iy);
            Double_t base = hbase->GetBinContent(ix,iy);
            Double_t mu = hmu->GetBinContent(ix,iy);
            Double_t sigma = hsigma->GetBinContent(ix,iy);
            Double_t mean = hmean->GetBinContent(ix,iy);
            Double_t rms = hrms->GetBinContent(ix,iy);
            Double_t max = hmax->GetBinContent(ix,iy);
            Double_t peak = hpeak->GetBinContent(ix,iy);
            if (mu < 0 || mu > steps || amp*sigma < 1 || base == 0 || max < 1)
            { 
               hbase->SetBinContent(ix,iy,0);
               hamp->SetBinContent(ix,iy,0);
               hmu->SetBinContent(ix,iy,0);
               hsigma->SetBinContent(ix,iy,0);
            }
            else {
               hmu->SetBinContent(ix,iy,mu*angleScale + dmu);
               hsigma->SetBinContent(ix,iy,sigma*angleScale);
               hmean->SetBinContent(ix,iy,mean*angleScale);
               hrms->SetBinContent(ix,iy,rms*angleScale);
               hpeak->SetBinContent(ix,iy,peak*angleScale);
            }
         }
      }
   }

   // number of contours encodes whether units are microradians (=100)
   // or arbitrary units of intensity (=99)
   hbase->SetContour(99);
   hamp->SetContour(99);
   hmu->SetContour(100);
   hsigma->SetContour(100);
   hmean->SetContour(100);
   hrms->SetContour(100);
   hmax->SetContour(99);
   hrms->SetContour(100);
#endif
   std::cout << "done" << std::endl;
   return 0;
}

Int_t rcpicker::SliceX(TH2D *rcimage, Int_t ix)
{
   Draw(rcimage);
   Int_t nbinx = rcimage->GetNbinsX();
   Int_t nbiny = rcimage->GetNbinsY();
   if (ix < 0 || ix >= nbinx) {
      std::cerr << "rcpicker::SliceX error - argument ix must be between "
                << 0 << " and " << nbinx << std::endl;
      return 1;
   }
   delete fShow;
   Double_t distanceScale = GetMicronsPerPixel()/1000;
   Double_t uvratio = GetPixelAspectRatio();
   TString title;
   title.Form("slice of %s through column %d",rcimage->GetTitle(),ix);
   fShow = new TH1D("xslice",title,nbiny,0,nbiny*distanceScale*uvratio);
   for (Int_t iy=1; iy <= nbiny; ++iy) {
      fShow->SetBinContent(iy,rcimage->GetBinContent(ix,iy));
   }
   if (fCanvas1) {
      fCanvas1->cd();
   }
   else {
      fCanvas1 = new TCanvas("c1","single-pixel rocking curve",
                             610,10,600,500);
   }
   fCanvas1->SetLeftMargin(0.15);
   fCanvas1->SetTopMargin(0.12);
   fShow->Draw();
   fCanvas1->Update();
   return 0;
}

Int_t rcpicker::SliceY(TH2D *rcimage, Int_t iy)
{
   Draw(rcimage);
   Int_t nbinx = rcimage->GetNbinsX();
   Int_t nbiny = rcimage->GetNbinsY();
   if (iy < 0 || iy >= nbiny) {
      std::cerr << "rcpicker::SliceY error - argument iy must be between "
                << 0 << " and " << nbiny << std::endl;
      return 1;
   }
   delete fShow;
   Double_t distanceScale = GetMicronsPerPixel()/1000;
   TString title;
   title.Form("slice of %s through row %d",rcimage->GetTitle(),iy);
   fShow = new TH1D("yslice",title,nbinx,0,nbinx*distanceScale);
   for (Int_t ix=1; ix <= nbinx; ++ix) {
      fShow->SetBinContent(ix,rcimage->GetBinContent(ix,iy));
   }
   if (fCanvas1) {
      fCanvas1->cd();
   }
   else {
      fCanvas1 = new TCanvas("c1","single-pixel rocking curve",
                             610,10,600,500);
   }
   fCanvas1->SetLeftMargin(0.15);
   fCanvas1->SetTopMargin(0.12);
   fShow->Draw();
   fCanvas1->Update();
   return 0;
}

Int_t rcpicker::Draw(TH2D *rcimage)
{
   gStyle->SetPalette(1,0);
   gStyle->SetNumberContours(50);
   gStyle->SetOptStat(0);
   gStyle->SetOptFit(1);
   gStyle->SetStatX(0.95);
   gStyle->SetStatY(0.93);
   gStyle->SetStatW(0.18);
   gStyle->SetStatH(0.15);
   fCanvas2 = (TCanvas*)gROOT->FindObject("c2");
   if (fCanvas2 == 0) {
      double imagewidth = rcimage->GetXaxis()->GetNbins();
      double imageheight = rcimage->GetYaxis()->GetNbins();
      double plotheight = 400;
      double plotwidth = plotheight*imagewidth/imageheight;
      double leftmargin = 50;
      double rightmargin = 90;
      double topmargin = 50;
      double bottommargin = 50;
      double canvaswidth = plotwidth+leftmargin+rightmargin;
      double canvasheight = plotheight+topmargin+bottommargin;
      fCanvas2 = new TCanvas("c2","rocking curve display",50,10,
                             int(canvaswidth),int(canvasheight));
      fCanvas2->SetLeftMargin(leftmargin/canvaswidth);
      fCanvas2->SetRightMargin(rightmargin/canvaswidth);
      fCanvas2->SetTopMargin(topmargin/canvasheight);
      fCanvas2->SetBottomMargin(bottommargin/canvasheight);
   }
   fCanvas2->cd();
   rcimage->GetXaxis()->SetTitle("u (mm)");
   rcimage->GetYaxis()->SetTitle("v (mm)");
   rcimage->GetYaxis()->SetTitleOffset(1.1);
   if (rcimage->GetContour() == 100) {
      rcimage->GetZaxis()->SetTitle("#mu^{}rad");
      rcimage->GetZaxis()->SetTitleOffset(1.5);
   }
   TString title;
   title.Form("%s, %s scan %d", rcimage->GetTitle(), fPrefix.Data(), fSeq);
   rcimage->SetTitle(title);
   rcimage->Draw("colz");
   fCanvas2->Update();
   return 0;
}

TH2D *rcpicker::Crop(TH2D* rcimage, int *ixmin, int *ixmax,
                                    int *iymin, int *iymax)
{
   int nbinx = rcimage->GetNbinsX();
   int nbiny = rcimage->GetNbinsY();
   if (*ixmin == 1 && *ixmax == nbinx && *iymin == 1 && *iymax == nbiny) {
      return rcimage;
   }
   if (*ixmin == 0 || *ixmax == 0 || *iymin == 0 || *iymax == 0) {
      int xmin=1e9, xmax=0;
      int ymin=1e9, ymax=0;
      for (int ix = 1; ix <= nbinx; ++ix) {
         for (int iy = 1; iy <= nbiny; ++iy) {
            if (rcimage->GetBinContent(ix,iy) != 0) {
               xmin = (ix < xmin)? ix : xmin;
               xmax = (ix > xmax)? ix : xmax;
               ymin = (iy < ymin)? iy : ymin;
               ymax = (iy > ymax)? iy : ymax;
            }
         }
      }
      *ixmin = (*ixmin == 0)? xmin : *ixmin;
      *ixmax = (*ixmax == 0)? xmax : *ixmax;
      *iymin = (*iymin == 0)? ymin : *iymin;
      *iymax = (*iymax == 0)? ymax : *iymax;
   }
   double xmin = rcimage->GetXaxis()->GetBinLowEdge(*ixmin);
   double xmax = rcimage->GetXaxis()->GetBinUpEdge(*ixmax);
   double ymin = rcimage->GetYaxis()->GetBinLowEdge(*iymin);
   double ymax = rcimage->GetYaxis()->GetBinUpEdge(*iymax);
   nbinx = *ixmax - *ixmin + 1;
   nbiny = *iymax - *iymin + 1;
   TString name(rcimage->GetName());
   TString title(rcimage->GetTitle());
   rcimage->SetName("temp_dummy");
   TH2D *hcrop = new TH2D(name, title, nbinx, xmin, xmax, nbiny, ymin, ymax);
   for (int ix = 1; ix <= nbinx; ++ix) {
      for (int iy = 1; iy <= nbiny; ++iy) {
         double x = hcrop->GetXaxis()->GetBinCenter(ix);
         double y = hcrop->GetYaxis()->GetBinCenter(iy);
         double z = rcimage->GetBinContent(rcimage->GetXaxis()->FindBin(x),
                                           rcimage->GetYaxis()->FindBin(y));
         hcrop->SetBinContent(ix, iy, z);
      }
   }
   hcrop->GetXaxis()->SetTitle(rcimage->GetXaxis()->GetTitle());
   hcrop->GetYaxis()->SetTitle(rcimage->GetYaxis()->GetTitle());
   hcrop->GetYaxis()->SetTitleOffset(1.4);
   delete rcimage;
   return hcrop;
}