/////////////////////////////////////////
//                                     //
//    PandaFieldMapSim1 source         //
//                                     //
//    @2006, George Serbanut           //
//                                     //
/////////////////////////////////////////

//g++ headers
#include <iostream>
#include <fstream>
#include <vector>

//ROOT specific headers
#include "TSystem.h"
#include "TString.h"
#include "TMath.h"

//class specific headers
#include "PandaFieldMapSim1.h"
#include "CbmField.h"

//standard constructor
PFMS1::PFMS1() : CbmField("PFMS1","PANDA magnetic Field Map for Simulation v.1")
{
    fType = 2;
}

PFMS1::~PFMS1()
{
    delete (CbmField*)this;
}

void PFMS1::Init()
{
    TString tmpFN;
    TString dir = getenv("VMCWORKDIR");
    tmpFN = dir + "/bfield/bfield.map";
    ReadAsciiFile(tmpFN);
}

void PFMS1::ReadAsciiFile(const char* fileName)
{
    ifstream f;
    BField* bf;
    string buffer;
    
    f.open(fileName);
    if(f.is_open()) {
	getline(f,buffer,'\n');
	getline(f,buffer,'\n');
	while(f) {
	    bf = new BField;
	    f>>bf->r>>bf->z>>bf->theta>>bf->Br>>bf->Bz>>bf->Btheta;
	    vbf.push_back(bf);
	};
	vbf.pop_back();
	f.close();
	cout<<"B-field map loaded with "<<vbf.size()<<" values"<<endl;
    } else {
	cout<<"bfield.map not found!!!"<<endl;
	exit(1);
    };
    f.clear();
    delete bf;
    
    return;
}

void PFMS1::GetFieldValue(const Double_t Point[3], Double_t *Bfield)
{
//  Double_t unit = 10.0; //conversion T->kG
//cm:
  Double_t rmin = 0.0;
  Double_t rmax = 250.0;
  Double_t zmin = -160.0;
  Double_t zmax = 310.0;

  Double_t theta, r, z;
  double theta_a, r_a, z_a;
  int pos;

  theta = TMath::ATan2(Point[1],Point[0]);

  r = TMath::Sqrt((Point[0]*Point[0])+(Point[1]*Point[1]))/10.0;
  z = Point[2]/10.0;
  
  theta_a = (double)theta;
  r_a = (double)r;
  z_a = (double)z;
  
  Bfield[0] = 0.0;
  Bfield[1] = 0.0;
  Bfield[2] = 0.0;
  
  double bField[3];
  bField[0] = 0.0;
  bField[1] = 0.0;
  bField[2] = 0.0;

  if((r>=rmin)&&(r<=rmax))
  if((z>=zmin)&&(z<=zmax)) {
    pos = static_cast<int>(floor(r_a/2.0))*2351+static_cast<int>(floor((z+1600.0)/2.0));
    bField[0] = (vbf[pos]->Br)*cos(theta_a)*10.0;
    bField[1] = (vbf[pos]->Br)*sin(theta_a)*10.0;
    bField[2] = (vbf[pos]->Bz)*10.0;
    Bfield[0] = (Double_t)bField[0];
    Bfield[1] = (Double_t)bField[1];
    Bfield[2] = (Double_t)bField[2];
  };
  
  return;
}

Double_t PFMS1::GetBx(Double_t x, Double_t y, Double_t z)
{
    Double_t Point[3], B[3];
    
    Point[0] = x;
    Point[1] = y;
    Point[2] = z;
    
    GetFieldValue(Point, B);

    return B[0];
}

Double_t PFMS1::GetBy(Double_t x, Double_t y, Double_t z)
{
    Double_t Point[3], B[3];
    
    Point[0] = x;
    Point[1] = y;
    Point[2] = z;
    
    GetFieldValue(Point, B);

    return B[1];
}

Double_t PFMS1::GetBz(Double_t x, Double_t y, Double_t z)
{
    Double_t Point[3], B[3];
    
    Point[0] = x;
    Point[1] = y;
    Point[2] = z;
    
    GetFieldValue(Point, B);

    return B[2];
}

ClassImp(PFMS1)