BfldCanvas.cxx

Go to the documentation of this file.

#include "BField/Viz/BfldCanvas.h"

#include <iostream>
#include <map>
using namespace std;

#include "TMath.h"
#include "TStyle.h"
#include "TPolyLine.h"
#include "TEllipse.h"
#include "TStopwatch.h"

#include "Conventions/Munits.h"


int    BfldCanvas::fgVerbose  = 0;
double BfldCanvas::fgKGMin    =  0.0;    // 10.0, 15.0 ..
double BfldCanvas::fgKGMax    = 22.0;    // 20.0, 18.5, 24.0 ...
int    BfldCanvas::fgArrowGridX = 40;
BfldInterpMethod::InterpMethod_t BfldCanvas::fgDefaultInterpMethod 
                              = BfldInterpMethod::kDefault;


BfldCanvas::BfldCanvas()
  : fCanvas(0)
  , fScatter(0)
  , fPhase(0)
  , fIsCoil(0)
  , fVldc()
  , fZ(-999)
  , fScaleMax(1.0)
  , fArrowScale(1.0)
  , fInterpMethod(fgDefaultInterpMethod)
{
  // default ctor for ROOTCINT
}

BfldCanvas::~BfldCanvas()
{
  if (fCanvas)  delete fCanvas;
  if (fScatter) delete fScatter;
  if (fPhase)   delete fPhase;
}

BfldCanvas::BfldCanvas(VldContext& vldc, std::string tag, int iscoil)
  : fCanvas(0)
  , fScatter(0)
  , fPhase(0)
  , fIsCoil(iscoil)
  , fVldc(vldc)
  , fZ(-999)
  , fScaleMax(1.0)
  , fArrowScale(1.0)
  , fInterpMethod(fgDefaultInterpMethod)
{
  // some initial defaults
  double xmin = -5.0*Munits::m, xmax = +5.0*Munits::m;
  double ymin = -5.0*Munits::m, ymax = +5.0*Munits::m;
  double scalea = 1.0;   // scale for arrows
  int  xpixels = 800, ypixels = 800;
  int  nx=500, ny=500;  // grid for color map
  int  nxa=40, nya=40;  // grid for arrows

  DetermineLimits(xmin,xmax,ymin,ymax,xpixels,ypixels);
  DetermineColorGrid(nx,ny);
  DetermineArrowGrid(nxa,nya,scalea);
  
  std::string cname = tag;   cname += "_canvas";
  std::string sname = tag;
  std::string pname = tag;   pname += "_phase";
  fCanvas  = new TCanvas(cname.c_str(),cname.c_str(),xpixels,ypixels);
  fScatter = new TH2F(sname.c_str(),sname.c_str(),nx,xmin,xmax,ny,ymin,ymax);
  fPhase   = new TH2F(pname.c_str(),pname.c_str(),nxa,xmin,xmax,nya,ymin,ymax);
  AdjustLabels();
  SetScaleMax();
  
}

void BfldCanvas::SetScaleMax(double scale) 
{ 
  fScaleMax = scale;
  if ( fScatter ) {
    fScatter->SetMinimum(fgKGMin*Munits::kilogauss);
    fScatter->SetMaximum(fgKGMax*Munits::kilogauss*fScaleMax);
  }
}

void BfldCanvas::SetTitle(std::string title)
{
  std::string name = title.substr(0,title.find(" "));

  //static int icnt = 0; ++icnt;
  //std::string nameuniq = Form("_%d",icnt);
  static std::map<std::string,int> namecnt;
  int icnt = namecnt[title]++;
  std::string nameuniq = "";
  if ( icnt > 0 ) nameuniq = Form("_%d",icnt);

  std::string cname = name;
  cname += "_canvas";
  fCanvas->SetTitle(cname.c_str());
  cname += nameuniq;
  fCanvas->SetName(cname.c_str());

  fScatter->SetTitle(""); // (title.c_str());
  title += "_phase";
  fPhase->SetTitle(title.c_str());
  
  name += nameuniq;
  fScatter->SetName(name.c_str());
  name += "_phase";
  fPhase->SetName(name.c_str());

}

void BfldCanvas::Fill(BField& bfld)
{
  Fill(&bfld,0);
}
void BfldCanvas::Fill(BfldHandler& bhandler)
{
  Fill(0,&bhandler);
}

void BfldCanvas::Draw()
{
  fCanvas->Draw();
  fCanvas->cd();
  AdjustMargins();

  fScatter->SetContour(gStyle->GetNumberOfColors());
  fScatter->Draw("COLZ");

  if ( fArrowScale > 0 ) DrawArrows();

  // coil region limits
  //draw_border(0.298); 
  //draw_border(0.448); 

  fCanvas->Update();

}

void BfldCanvas::DetermineLimits(double& xmin, double& xmax,
                                 double& ymin, double& ymax,
                                 int& xpixels, int& ypixels)
{
  xmin = -5.0*Munits::m; xmax = +5.0*Munits::m;
  ymin = -5.0*Munits::m; ymax = +5.0*Munits::m;

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

  if ( fIsCoil > 0 ) {

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

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

  xpixels = ypixels = 800; // canvas size 
  xpixels = TMath::Nint(ypixels * (xmax-xmin)/(ymax-ymin));
}

void BfldCanvas::DetermineColorGrid(int& nx, int& ny)
{ 
  nx  = ny  = 500;  // grid for color map 
}

void BfldCanvas::DetermineArrowGrid(int& nxa, int& nya, double& scalea)
{
  nxa = nya = 40;
  scalea = 1.0; 

  if ( fIsCoil ) {
    // change arrow size, grid
    nxa = nya = 25;
    scalea = 2.5;
    
    nxa = nya = 40;
    scalea = 2.0;
  }

  double xmin, xmax, ymin, ymax;
  int xpixels, ypixels;
  DetermineLimits(xmin,xmax,ymin,ymax,xpixels,ypixels);

  scalea = 1.5*5.0*(Munits::cm)*(xmax-xmin)/10.0*(Munits::m);

  scalea *= fArrowScale;     // scale for arrows size

  // change number of phase histogram bins in X
  nxa = TMath::Nint((float)(nxa) * (float)fgArrowGridX/40.0);

  // keep aspect ratio
  nya = TMath::Nint((float)(nxa) * (ymax-ymin)/(xmax-xmin));

} 

void BfldCanvas::AdjustLabels()
{
  fScatter->SetStats(false);
  fScatter->SetDirectory(0);
  fScatter->SetXTitle("x [ m ]");
  fScatter->SetYTitle("y [ m ]");
  fScatter->SetZTitle("B [ T ]");
  double tsize = 0.045, lsize = 0.035;
  fScatter->GetXaxis()->SetTitleSize(tsize);
  fScatter->GetYaxis()->SetTitleSize(tsize);
  fScatter->GetZaxis()->SetTitleSize(tsize);
  fScatter->GetXaxis()->SetLabelSize(lsize);
  fScatter->GetYaxis()->SetLabelSize(lsize);
  fScatter->GetZaxis()->SetLabelSize(lsize);
}

void BfldCanvas::AdjustMargins()
{
  fCanvas->SetTopMargin(0.1625); // 0.16875);   // 0.175);   // 0.075
  fCanvas->SetBottomMargin(0.10);   // 0.070   0.145
  fCanvas->SetLeftMargin(0.10);     // 0.06 
  fCanvas->SetRightMargin(0.1875); // 0.19375); //0.20);  // 0.110   0.170  0.025 (for z)
}

void BfldCanvas::Fill(BField* bfld, BfldHandler* bhandler)
{

  TVector3 xyz, bxyz;
  
  xyz.SetZ(fZ);
  
  if (fgVerbose) cout << " start loop over grid " << endl;
  
  // fill the TH2F's, first the color grid, then the phase
  for (int pass=0; pass<2; ++pass) {
    TH2F* h2 = ((pass==0) ? fScatter : fPhase );
    double x, y;
    int ix, iy;

    // don't rely on functions for limit/grid info 
    // take it from scatter plot

    TAxis* xaxis = h2->GetXaxis();
    int    nx    = xaxis->GetNbins();
    double xmin  = xaxis->GetXmin();
    double xmax  = xaxis->GetXmax();

    TAxis* yaxis = h2->GetYaxis();
    int    ny    = yaxis->GetNbins();
    double ymin  = yaxis->GetXmin();
    double ymax  = yaxis->GetXmax();

    if (fgVerbose) {
      cout << " x " << xmin << " " << xmax << " nx " << nx << endl;
      cout << " y " << ymin << " " << ymax << " ny " << ny << endl;
    }

    //double dx = (xmax-xmin)/(double)nx;
    //double dy = (ymax-ymin)/(double)ny;
  
    for (ix=1; ix<=nx; ++ix) {
      x = xaxis->GetBinCenter(ix);
      xyz.SetX(x);
      for (iy=1; iy<=ny; ++iy) {
        y = yaxis->GetBinCenter(iy);
        xyz.SetY(y);
        
        //int ihashbin = h2->GetBin(ix,iy);

        if ( bfld ) {
          bxyz = bfld->GetBField(xyz);
        } else if ( bhandler ) {
          bxyz = bhandler->GetBField(xyz,false);
        }
        
        double bmag  = bxyz.Mag();
        double value = bmag;
        if ( pass == 1 ) {
          double phi   = ( (bmag > 0.0) ? bxyz.Phi() : 42 );
          value = phi;
        }
        h2->Fill(x,y,value);
      
#ifdef RWHDEBUG
        double r2   = x*x + y*y;
        double r2lo = 0.09 *  0.09;
        double 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 (fgVerbose) cout << " pos " << xyz << " B " << bxyz << endl;
        }
#endif
        
      } // done y
    } // done x

  } // done pass
  
}

void BfldCanvas::DrawArrows()
{

  int i,j;
  double scalea;
  DetermineArrowGrid(i,j,scalea); // only use "scalea"

  TAxis* xaxis = fPhase->GetXaxis();
  int    nx    = xaxis->GetNbins();
  //double xmin  = xaxis->GetXmin();
  //double xmax  = xaxis->GetXmax();

  TAxis* yaxis = fPhase->GetYaxis();
  int    ny    = yaxis->GetNbins();
  //double ymin  = yaxis->GetXmin();
  //double ymax  = yaxis->GetXmax();

  double x, y;
  int ix, iy;

  //double dx = (xmax-xmin)/(double)nx;
  //double dy = (ymax-ymin)/(double)ny;
  
  for (ix=1; ix<=nx; ++ix) {
    x = xaxis->GetBinCenter(ix);
    for (iy=1; iy<=ny; ++iy) {
      y = yaxis->GetBinCenter(iy);

      int ihashbin = fPhase->GetBin(ix,iy);
      double phi = fPhase->GetBinContent(ihashbin);
      if ( phi > TMath::TwoPi() ) continue;  // don't draw nonsense

      double bx = TMath::Cos(phi);
      double by = TMath::Sin(phi);
      DrawCheveron(x,y,bx,by,scalea);
   
    } // done y
  } // done x

}

void BfldCanvas::DrawCheveron(double x0, double y0, double dx, double dy, 
                              double scale)
{
  static double xp[] = { -0.16, -0.50, 0.60, -0.50, -0.16};
  static double yp[] = {  0.00,  0.25, 0.00, -0.25,  0.00};
  double x[5], y[5];
  double len = TMath::Sqrt(dx*dx+dy*dy);
  if (len<=0) return;
  double c = dx/len;
  double 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(1);
  cheveron.DrawPolyLine(5,x,y);
  //static TEllipse  circle;
  //circle.SetLineWidth(1);
  //circle.SetFillStyle(0);  // transparent fill
  //circle.DrawEllipse(x0,y0,scale,0.,0.,360.,0.,"");
}

std::string BfldCanvas::GetZoomInterpString()
{
  std::string result = ( GetIsCoil() ? "_coil" : "_full");
  int im = GetInterpMethod();
  if ( im != (int)BfldInterpMethod::kDefault )
    result += Form("_im%d",im);
  return result;
}

---