#if  !defined(__CINT__) || defined(__MAKECINT__)
#include <iostream>
#include <iomanip>
using namespace std;

#include "TVector3.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TArrow.h"
#include "TFile.h"
#include "TMath.h"
#include "TPolyLine.h"
#include "TStopwatch.h"

#include "Validity/VldContext.h"
#include "BField/BField.h"
#include "BField/BfldDbiPlaneMap.h"
#include "BField/BfldLoanPool.h"
#include "Conventions/Munits.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliSteelPlnHandle.h"
#include "DatabaseInterface/DbiResultPtr.h"
#endif

// forward declarations
std::ostream& operator<<(std::ostream& os, const TVector3& tv3);
void do_one_z(TCanvas* c, TH2F* s, Double_t z, BField& bfld, bool verbose);
void draw_cheveron(Double_t x0, Double_t y0, 
                   Double_t dx, Double_t dy, Double_t scale);
Double_t rescale(Double_t value);
void adjust_margins();
void adjust_labels(TH2F* hist);
void dump_planemaps(VldContext& vldc);

// some defaults
const Float_t kgauss_min_default     =  0.0;  // 10.0
//const Float_t kgauss_max_default     = 24.0;  // 18.5 or 24.0
const Float_t kgauss_max_default     = 20.0;  // 18.5
const Float_t kgauss_max_gap_default =  2.0; 

// some globals
UInt_t  nx=500, ny=500;  // grid for color map
//UInt_t  nx=100, ny=100;  // grid for color map
//UInt_t  nx=50, ny=50;  // grid for color map
UInt_t  nxa=60, nya=60;  // grid for arrows
Double_t scalea = 1.0;   // scale for arrows
Float_t xmin = -5.0*Munits::m, xmax = +5.0*Munits::m;
Float_t ymin = -5.0*Munits::m, ymax = +5.0*Munits::m;

// main function
void bfld_slices(Detector::Detector_t detector = Detector::kFar, 
                 Int_t ipln0=0, Int_t ipln1=12,
                 bool coil_region = true,
                 bool doprint = false,
                 Float_t kilogauss_min = kgauss_min_default, 
                 Float_t kilogauss_max = kgauss_max_default,
                 Float_t kilogauss_max_gap = kgauss_max_gap_default,
                 bool do_gap = true,
                 bool write_root_file = false)
{

  bool verbose = false;
  
  BfldDbiPlaneMap::SetDefensiveUnpkg(true);

  if (verbose) cout << "adjust the limits" << endl;

  xmin = -5.0*Munits::m; xmax = +5.0*Munits::m;
  ymin = -5.0*Munits::m; ymax = +5.0*Munits::m;

  if (Detector::kNear == detector) {
      xmin = -2.8*Munits::m, xmax = +4.0*Munits::m;
      ymin = -2.4*Munits::m, ymax = +2.4*Munits::m;
  }

  const char* region_name="";
  if (coil_region) {
      region_name = "_coil";

      xmin = -0.35*Munits::m, xmax = +0.35*Munits::m;
      ymin = -0.35*Munits::m, ymax = +0.35*Munits::m;

      xmin = -0.55*Munits::m, xmax = +0.55*Munits::m;
      ymin = -0.55*Munits::m, ymax = +0.55*Munits::m;

      // change arrow size, grid
      nxa = nya = 25;
      scalea = 2.5;

      nxa = nya = 40;
      scalea = 2.0;

      //xmin =  0.15*Munits::m, xmax =  0.23*Munits::m;
      //ymin = -0.12*Munits::m, ymax = -0.04*Munits::m;

      //xmin =  0.66*Munits::m, xmax =  0.84*Munits::m;
      //ymin = -0.49*Munits::m, ymax = -0.31*Munits::m;
      //kilogauss_min = 0.9/Munits::kilogauss;
      //kilogauss_max = 1.6/Munits::kilogauss;
  }

  Float_t xpixels = 800, ypixels = xpixels * (ymax-ymin)/(xmax-xmin);
  nya = (float)(nxa) * (ymax-ymin)/(xmax-xmin);
  
  if (false) { nxa = nya = 0; }

  VldTimeStamp now;
  
  gStyle->SetPalette(1,(Int_t*)0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetFrameFillColor(0);
  gStyle->SetFrameFillStyle(0);
  gStyle->SetPadColor(0);
  gStyle->SetTextSize(0.015);

  SimFlag::SimFlag_t sflg = SimFlag::kData;
  //SimFlag::SimFlag_t sflg = SimFlag::kMC;
  //if (Detector::kNear == detector ) sflg = SimFlag::kReroot;

  VldContext vldc(detector,sflg,now);
  const char* detname = Detector::AsString(detector);

  dump_planemaps(vldc);

  PlexPlaneId plnid[3];
  plnid[0] = PlexPlaneId(detector,ipln0,kTRUE);
  plnid[1] = PlexPlaneId(detector,ipln1,kTRUE);
  plnid[2] = PlexPlaneId(detector,ipln0,kTRUE);

  TString fname   = Form("DetMap_%s_pln%d_pln%d.root",detname,ipln0,ipln1);

  if (verbose) cout << " process " << vldc << endl;

  UgliGeomHandle ugh(vldc);

  if (verbose) cout << " build bfield " << vldc << endl;

  // do BfldLoanPool config *before* building BField object
  BfldLoanPool* bfldpool = BfldLoanPool::Instance();
  bfldpool->Set("UseDCSCoilDir=0"); // just so we don't screw sign up
  bfldpool->Set("DoLocalTransform=2"); // tranform both position and field
  bfldpool->Set("ApplyBdotScale=1");
  if ( do_gap ) {
    cout << "Default BfldLoanPool Configuration" << endl;
    bfldpool->Print();
    bfldpool->Set("IgnoreUseEverywhere=1");
    bfldpool->Set("RequireInZTest=3");
    bfldpool->Set("DoInterPlaneField=1");
  }
  bfldpool->Update(); // don't forget to signal it to update
  cout << "New BfldLoanPool Configuration" << endl;
  bfldpool->Print();

  //bool dointerp = true;

  BField bfld(vldc);
  bfld.SetInterpMethod(BfldInterpMethod::kDefault);
  //bfld.SetInterpMethod(BfldInterpMethod::kClosest);
  //bfld.SetInterpMethod(BfldInterpMethod::kBilinear);
  //bfld.SetInterpMethod(BfldInterpMethod::kPlanar);
  //bfld.SetInterpMethod(BfldInterpMethod::kPlanarVec);

  //cout << "bfld.GetRequireInZTest()  " << bfld.GetRequireInZTest() << endl;
  //if ( bfld.GetRequireInZTest() < 2 ) {
  //  cout << " HEY! That's not supposed to be <2 " << endl;
  //  return;
  //}

  if (verbose) cout << " create TH2F " << endl;

  TString pngname[3], psname[3];
  char what[1024];
  TCanvas* c[3];
  TH2F* s[3];
  Double_t z[3];

  int nslice = ( ( do_gap && coil_region ) ? 3 : 2 );
  for (int i=0; i<nslice; ++i) {
    if ( i != 2 ) {
      sprintf(what,"pln%d",plnid[i].GetPlane());
      UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(plnid[i]);
      z[i] = steelpln.GetZ0();
    } else {
      sprintf(what,"gap%d_%d",plnid[0].GetPlane(),plnid[1].GetPlane());
      //z[i] = 0.5*(z[0]+z[1]);
      z[i] = z[0] + 0.5*(z[1]-z[0])/(ipln1-ipln0);
    }
    TString cname   = Form("DetMap_%s_%s%sCanvas",detname,what,region_name);
    TString ctitle  = Form("DetMap %s %s %s Canvas",detname,what,region_name);
    TString hname   = Form("DetMap_%s_%s%s",detname,what,region_name);
    TString htitle  = Form("DetMap %s %s %s",detname,what,region_name);
    pngname[i] = Form("DetMap_%s_%s%s.png",detname,what,region_name);
    psname[i] = Form("DetMap_%s_%s%s.ps",detname,what,region_name);

    if (verbose) cout << " create canvas " << cname << endl;
    cout << "hname: " << hname << endl;
    c[i] = new TCanvas(cname,ctitle,xpixels,ypixels);  
    s[i] = new TH2F(hname,htitle,nx,xmin,xmax,ny,ymin,ymax);
    adjust_labels(s[i]);

    if ( i != 2 ) {
      s[i]->SetMaximum(kilogauss_max*Munits::kilogauss);
    } else {
      s[i]->SetMaximum(kilogauss_max_gap*Munits::kilogauss);
    }
    //s[i]->SetMinimum(0.0);
    s[i]->SetMinimum(kilogauss_min*Munits::kilogauss);

    do_one_z(c[i],s[i],z[i],bfld,verbose);
  }

  if ( coil_region || ! coil_region ) {
    // hackery here ....
    TCanvas* cdiffer = new TCanvas("differCanvas","differCanvas",xpixels,ypixels);
    TH2F*    sdiffer = new TH2F("differ","differ",nx,xmin,xmax,ny,ymin,ymax);
    adjust_labels(sdiffer);

    TCanvas* cfrac = new TCanvas("fracCanvas","fracCanvas",xpixels,ypixels);
    TH2F*    sfrac = new TH2F("frac","frac",nx,xmin,xmax,ny,ymin,ymax);
    adjust_labels(sfrac);

    TH2F* s1 = s[1];
    TH2F* s0 = s[0];
    //*sdiffer = *s1 - *s0;
    //sdiffer->SetMinimum(-kilogauss_max*Munits::kilogauss);
    Float_t dx = (xmax-xmin)/(float)nx;
    Float_t dy = (ymax-ymin)/(float)ny;
    Float_t x, y;
    for (UInt_t ix=0; ix<nx; ++ix) {
      x = xmin + ((float)ix+0.5)*dx;
      for (UInt_t iy=0; iy<ny; ++iy) {
        y = ymin + ((float)iy+0.5)*dy;
        Int_t indxbin = sdiffer->FindBin(x,y);
        Float_t val1 = s1->GetBinContent(indxbin);
        Float_t val0 = s0->GetBinContent(indxbin);
        Float_t absdiff = TMath::Abs(val1 - val0);
        sdiffer->Fill(x,y,absdiff);
        Float_t ave = 0.5*TMath::Abs(val1 + val0);
        Float_t ferr = 0.0;
        if (ave != 0.0) ferr = (val1-val0)/ave;
        sfrac->Fill(x,y,ferr);
      }
    }
    cdiffer->Draw();
    cdiffer->cd();
    adjust_margins();
    sdiffer->SetMaximum(0.2*Munits::tesla);
    sdiffer->SetContour(gStyle->GetNumberOfColors());
    sdiffer->Draw("COLZ");
    cdiffer->Update();

    cfrac->Draw();
    cfrac->cd();
    adjust_margins();
    sfrac->SetMaximum( 0.2);
    sfrac->SetMinimum(-0.2);
    sfrac->SetContour(gStyle->GetNumberOfColors());
    sfrac->Draw("COLZ");
    cfrac->Update();
  }

  if (write_root_file) {
      TFile hfile(fname,"RECREATE");
      for (int j=0; j<3; ++j) s[j]->Write();
      hfile.Write();
      hfile.Close();
  }
  if (doprint) {
    for (int k=0; k<3; ++k) {
      c[k]->Print(pngname[k]);
      c[k]->Print(psname[k]);
    }
  }

  //dump_planemaps(vldc);

}

void do_one_z(TCanvas* c, TH2F* s, Double_t z, BField& bfld, bool verbose)
{
  
  TStopwatch sw;
  sw.Start();

  if (verbose) cout << " lookup z within steel plane 0 " << endl;
  
  Float_t x,y;
  TVector3 xyz, bxyz;
  
  xyz.SetZ(z);
  
  if (verbose) cout << " start loop over grid " << endl;
  
  Float_t dx = (xmax-xmin)/(float)nx;
  Float_t dy = (ymax-ymin)/(float)ny;
  
  for (UInt_t ix=0; ix<nx; ++ix) {
    x = xmin + ((float)ix+0.5)*dx;
    //      cout << " x = " << x << " x0 " << x0 << " dx " << dx 
    //           << " ix " << ix << endl;
    
    xyz.SetX(x);
    for (UInt_t iy=0; iy<ny; ++iy) {
      y = ymin + ((float)iy+0.5)*dy;
      xyz.SetY(y);
      
      bxyz = bfld.GetBField(xyz);
      Double_t value = bxyz.Mag();
      //value = rescale(value);
      s->Fill(x,y,value);
      
      Float_t r2   = x*x + y*y;
      Float_t r2lo = 0.09 *  0.09;
      Float_t r2hi = 0.12 *  0.12;
      //if (TMath::Abs(x)<10.*Munits::cm && TMath::Abs(y)<10.*Munits::cm ) {
      if ( false && r2lo < r2 && r2 < r2hi  ) {
        if (verbose) cout << " pos " << xyz << " B " << bxyz << endl;
      }
      
    } // done y
  } // done x
  
  if (verbose) cout << " draw the canvas " << endl;
  
  //c->Draw();
  //c->Divide(3,3);
  //c->cd(1);

  c->Draw();
  c->cd();
  adjust_margins();
  s->SetContour(gStyle->GetNumberOfColors());
  //s->Draw("BOX");
  if (true) {
    s->Draw("COLZ");
  } else {
    // test of contour plots for printing grey scale
    s->SetMinimum(0.25*20*Munits::kilogauss);
    s->Draw("cont1");
    c->Update();
    return; // skip arrows
  }

  if (verbose) cout << " start drawing arrows " << endl;

  TArrow arrow;
  bool keepscale = false;  // if true arrow reflect B magnitude
  arrow.SetFillStyle(0);  // transparent fill
  arrow.SetLineWidth(0.5);

  //oldarrow nx = ny = 50;
  //nx = ny = 60;
  dx = (xmax-xmin)/(float)nxa;
  dy = (ymax-ymin)/(float)nya;
  float dxy = TMath::Max(dx,dy);
  dx = dy = dxy;

  for (UInt_t ix=0; ix<nxa; ++ix) {
    x = xmin + ((float)ix+0.5)*dx;
    if (x>xmax) continue;

    //      cout << " x = " << x << " x0 " << x0 << " dx " << dx 
    //           << " ix " << ix << endl;
    
    xyz.SetX(x);
    for (UInt_t iy=0; iy<nya; ++iy) {
      y = ymin + ((float)iy+0.5)*dy;
      if (y>ymax) continue;

      xyz.SetY(y);
      
      bxyz = bfld.GetBField(xyz);
      
      if (false && TMath::Abs(x)>3.5 && TMath::Abs(y)>3.5) {
        if (verbose) cout << xyz << " B = " << bxyz 
                          << " mag " << bxyz.Mag() << endl;
      }
      
      //oldarrow Float_t bs = 0.08;
      //Float_t bs = 0.15;
      //Float_t head = 0.006; 
      Float_t bs = 0.09;
      Float_t head = 0.005; 
      if (!keepscale) {
        bs = 0.1;
        head = 0.006;
        if (bxyz.Mag2()>0) bxyz.SetMag(1.0);
      }
      
      Float_t bx = bxyz.X();
      Float_t by = bxyz.Y();
      //arrow.DrawArrow(x,y,x+bx*bs,y+by*bs,head,"|>");
      Float_t arrowScale = 6.5*(Munits::cm)*(xmax-xmin)/10.0*(Munits::m);
      arrowScale *= scalea;
      draw_cheveron(x,y,bx,by,arrowScale);

    } // done y
  } // done x
  
  if (verbose) cout << " arrows done, update canvas for " << endl;
  c->Update();  

  if (verbose && false) {
      xyz.SetY(0.0);
      xyz.SetX(-1*Munits::m);
      bxyz = bfld.GetBField(xyz);
      cout << "xyz=" << xyz << " B=" << bxyz << endl;
      xyz.SetY(0.0);
      xyz.SetX(-1*Munits::m);
      bxyz = bfld.GetBField(xyz);
      cout << "xyz=" << xyz << " B=" << bxyz << endl;
      xyz.SetY(0.0);
      xyz.SetX(+1*Munits::m);
      bxyz = bfld.GetBField(xyz);
      cout << "xyz=" << xyz << " B=" << bxyz << endl;
      xyz.SetY(-1*Munits::m);
      xyz.SetX(-1*Munits::m);
      bxyz = bfld.GetBField(xyz);
      cout << "xyz=" << xyz << " B=" << bxyz << endl;
  }

  sw.Stop();
  sw.Print();

}

std::ostream& operator<<(std::ostream& os, const TVector3& tv3)
{
  int w=9, p=4; // 1mm precision
  os << "[" << setw(w)   << setprecision(p) << tv3.x()
     << "," << setw(w)   << setprecision(p) << tv3.y()
     << "," << setw(w+1) << setprecision(p) << tv3.z()
     << "]";
 
  return os;
}


void draw_cheveron(Double_t x0, Double_t y0, Double_t dx, Double_t dy, Double_t scale)
{
    static Double_t xp[] = { -0.16, -0.50, 0.60, -0.50, -0.16};
    static Double_t yp[] = {  0.00,  0.25, 0.00, -0.25,  0.00};
    Double_t x[5], y[5];
    Double_t len = TMath::Sqrt(dx*dx+dy*dy);
    if (len<=0) return;
    Double_t c = dx/len;
    Double_t s = dy/len;
    // rotate, scale, translate
    for (int i=0; i<5; ++i) {
        x[i] = (xp[i]*c - yp[i]*s)*scale + x0;
        y[i] = (xp[i]*s + yp[i]*c)*scale + y0;
    }
    static TPolyLine cheveron;
    cheveron.SetLineWidth(0.25);
    cheveron.DrawPolyLine(5,x,y);
    //static TEllipse  circle;
    //circle.SetLineWidth(0.25);
    //circle.SetFillStyle(0);  // transparent fill
    //circle.DrawEllipse(x0,y0,scale,0.,0.,360.,0.,"");
}

Double_t rescale(Double_t value)
{
  // make scale nonlinear
  // map range [0.0,0.8] to [0.0,0.2]
  // map range [0.8,1.5] to [0.2,1.5]
  // map range [1.5,2.0] to [1.5,2.0]

  Double_t x[2] = { 0.8, 1.5 };
  Double_t y[2] = { 0.2, 1.5 };

  if ( value > x[1] ) return value;
  if ( value > x[0] )
    return y[0] + (value-x[0]) * (y[1]-y[0]) / (x[1]-x[0]); 
  else
    return value * y[0] / x[0];
}

void adjust_labels(TH2F* hist)
{
    hist->SetStats(false);
    hist->SetDirectory(0);
    hist->SetXTitle("x (m)");
    hist->SetYTitle("y (m)");
    hist->GetXaxis()->SetTitleSize(0.025);
    hist->GetYaxis()->SetTitleSize(0.025);
    hist->GetXaxis()->SetLabelSize(0.025);
    hist->GetYaxis()->SetLabelSize(0.025);
    hist->GetZaxis()->SetLabelSize(0.025);
}

void adjust_margins()
{
  gPad->SetTopMargin(0.075);
  gPad->SetBottomMargin(0.070);
  gPad->SetRightMargin(0.110);
  gPad->SetLeftMargin(0.06);
}

void dump_planemaps(VldContext& vldc)
{
  DbiResultPtr<BfldDbiPlaneMap> tblrows(vldc);
  unsigned int nrows = tblrows.GetNumRows();
  for (unsigned int irow=0; irow < nrows; ++irow) {
    const BfldDbiPlaneMap* onerow = tblrows.GetRow(irow);
    bool printit = false;
    int ipln = onerow->GetPlaneId().GetPlane();
    if (vldc.GetDetector() == Detector::kFar) {
      int lastsm0 = PlexPlaneId::LastPlaneFarSM0();
      int lastsm1 = PlexPlaneId::LastPlaneFarSM1();
      if (ipln < 12 || 
           (ipln > lastsm0-12 && ipln < lastsm0+12 ) || 
          ipln > lastsm1-12 ) printit = true;
    } else if ( vldc.GetDetector() == Detector::kNear ) {
      int lastspec = PlexPlaneId::LastPlaneNearSpect();
      if (ipln < 12 || ipln > lastspec-12 ) printit = true;
    }
    if (printit) onerow->Print();
  }
}