#ifndef __CINT__
#include <iostream>
#include <iomanip>
using namespace std;

#include "TVector3.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TArrow.h"

#include "Validity/VldContext.h"
#include "BField/BField.h"
#include "Conventions/Munits.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliSteelPlnHandle.h"
#endif

// forward declarations
std::ostream& operator<<(std::ostream& os, const TVector3& tv3);
void draw_cheveron(Double_t x0, Double_t y0, 
                   Double_t dx, Double_t dy, Double_t scale);
Double_t rescale(Double_t value);

// some defaults
const Float_t kgauss_min_default =  0.0;  // 10.0
const Float_t kgauss_max_default = 20.0;  // 18.5

// main function
void bfld_imap(Int_t imap_everywhere, 
               Detector::Detector_t detector = Detector::kFar,
               bool coil_region = false,
               bool doprint = false,
               Float_t kilogauss_min = kgauss_min_default, 
               Float_t kilogauss_max = kgauss_max_default,
               bool write_root_file = false)
{

  bool verbose = false;

  VldTimeStamp now;
  
  gStyle->SetPalette(1,(Int_t*)0);

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

  VldContext vldc(detector,sflg,now);
  TString cname   = Form("%sDetMap%dCanvas",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString ctitle  = Form("%sDet map %d Canvas",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString hname   = Form("%sDetMap%d",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString htitle  = Form("%sDet map %d",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString fname   = Form("%sDetMap%d.root",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString pngname = Form("%sDetMap%d.png",
                         Detector::AsString(detector),
                         imap_everywhere);
  TString psname  = Form("%sDetMap%d.ps",
                         Detector::AsString(detector),
                         imap_everywhere);

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

  UgliGeomHandle ugh(vldc);

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

  bool dointerp = true;

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

  if (verbose) cout << " create canvas " << cname << endl;

  UInt_t  nx=500, ny=500;
  UInt_t  nxa=75, nya=75;
  Float_t xmin = -5.0*Munits::m, xmax = +5.0*Munits::m;
  Float_t 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;
  }

  if (coil_region) {
      xmin = -0.35*Munits::m, xmax = +0.35*Munits::m;
      ymin = -0.35*Munits::m, ymax = +0.35*Munits::m;

      //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 = 1200, ypixels = xpixels * (ymax-ymin)/(xmax-xmin);
  nya = (float)(nxa) * (ymax-ymin)/(xmax-xmin);
  
  if (false) { nxa = nya = 0; }

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

  TCanvas* c = new TCanvas(cname,ctitle,xpixels,ypixels);  
  TH2F* s = new TH2F(hname,htitle,nx,xmin,xmax,ny,ymin,ymax);
  s->SetStats(false);
  s->SetDirectory(0);
  s->SetXTitle("x (m)");
  s->SetYTitle("y (m)");
  
  if (verbose) cout << " lookup z within steel plane 0 " << endl;
  
  Float_t x,y;
  TVector3 xyz, bxyz;
  
  PlexPlaneId pln0steel(detector,0,kTRUE);  // steel
  if (verbose) cout << " PlexPlaneId " << pln0steel << endl;
  UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(pln0steel);
  if (verbose) cout << " steelpln " << endl;
  xyz.SetZ(steelpln.GetZ0());
  
  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 ( 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();
  s->SetMaximum(kilogauss_max*Munits::kilogauss);
  s->SetMinimum(0.0);
  s->SetMinimum(kilogauss_min*Munits::kilogauss);
  s->SetContour(gStyle->GetNumberOfColors());
  //s->Draw("BOX");
  s->Draw("COLZ");
  
  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 = 5.0*(Munits::cm)*(xmax-xmin)/10.0*(Munits::m);
      draw_cheveron(x,y,bx,by,arrowScale);

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

  if (verbose) {
      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;
  }

  if (write_root_file) {
      TFile hfile(fname,"RECREATE");
      s->Write();
      hfile.Write();
      hfile.Close();
  }
  if (doprint) {
      c->Print(pngname);
      c->Print(psname);
  }

}

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];
}