#include <iostream>
#include <cstdlib>
#include <cmath>
#include <TRandom.h>

#include "IFieldManager.hxx"
#include "IOAMagneticField.hxx"
#include <ICOMETLog.hxx>
#include <HEPUnits.hxx>
#include <ICOMETEvent.hxx>
#include <ICOMETContext.hxx>
#include <IEventFolder.hxx>
#include <IOARuntimeParameters.hxx>
#include <ISIMG4Header.hxx>

static bool fMagnetCurrentInitialized = false; 

//********************************************************
COMET::IOAMagneticField::IOAMagneticField() {
//********************************************************
  Initialize();
}

//********************************************************
COMET::IOAMagneticField::~IOAMagneticField() {
//********************************************************
}

//********************************************************
void COMET::IOAMagneticField::Initialize() {
//********************************************************

  // As default set B Field error to zero
  fBErr = 0;
  //TEMPorary! These should come from a data base
  //NOTE: Don't use any of these in the GetInnerValue functions!
  fYMinC = 1.77 * unit::m;
  fYMaxC = 2.00 * unit::m;
  fZMinC = 3.55 * unit::m;
  fZMaxC = 3.78 * unit::m;
  fRMinC = 0.27 * unit::m;
  fRMaxC = 0.50 * unit::m;
  fXMinM = 1.760 * unit::m;
  fXMaxM = 2.815 * unit::m;
  fYMinM = 2.020 * unit::m;
  fYMaxM = 3.075 * unit::m;
  fZMaxM = 3.800 * unit::m;
  //fZMagnet = 327.45 * unit::cm;
  fZMagnet = 0 * unit::cm;
  fXCenl = 1.74 * unit::m;
  fYCenl = 1.60 * unit::m;
  fZCenl = 3.40 * unit::m;
  fDZCees = 0.956 * unit::m;
  fHZGaps = 0.038 * unit::m;
  fHSGap3 = 0.015 * unit::m;
  fBFMin = 0.00001 * unit::kilogauss;
  fBTMax = 4.0 * unit::kilogauss; 
  fBLMax = 7.4 * unit::kilogauss; 
  //More temporary stuff from EIKE
  fXMaxTPC = 1.12 * unit::m;
  fYMaxTPC = 1.18 * unit::m;
  fZMinTPC = -0.975 * unit::m;
  fZMaxTPC = 2.90 * unit::m;

  //fBTMax is used for the NOMAD field map. NOMAD used 0.4T as field strength 
  //so by default 0.4T is used for NOMAD. This can however be changed by 
  //SetFieldStrength. In COMET 0.2T will be used. So for the other field map
  //options 0.2T will be used. This can also be changed by SetFieldStrength.
  
  fFieldStrengthX = COMET::IOARuntimeParameters::Get().GetParameterD("CalibMagnet.FieldStrengthX") * unit::tesla;

  fFieldStrengthY = 0;
  fFieldStrengthZ = 0;


  //Initiaze values for the field map used so far by the SIMG4
  fYokeOuterWidth = 5870 * unit::mm;
  fYokeOuterHeight = 6086 * unit::mm;
  fYokeOuterLength = 7489 * unit::mm;
  fYokeInnerWidth = 3824 * unit::mm;
  fYokeInnerHeight = 4040 * unit::mm;
  fYokeInnerLength = 7489 * unit::mm;
  fCoilOuterWidth = 3824 * unit::mm;
  fCoilOuterHeight = 4020 * unit::mm;
  fCoilOuterLength = 7489 * unit::mm;
  fCoilInnerWidth = 3824 * unit::mm;
  fCoilInnerHeight = 3540 * unit::mm;
  fCoilInnerLength = 7009 * unit::mm;


  fNormalizationFromCurrent = (bool)COMET::IOARuntimeParameters::Get().GetParameterI("CalibMagnet.NormalizationFromCurrent");

  fReferenceCurrent  = COMET::IOARuntimeParameters::Get().GetParameterD("CalibMagnet.ReferenceCurrent");

  fCurrentToFieldC0  = COMET::IOARuntimeParameters::Get().GetParameterD("CalibMagnet.CurrentToFieldC0");
  fCurrentToFieldC1  = COMET::IOARuntimeParameters::Get().GetParameterD("CalibMagnet.CurrentToFieldC1");
  fCurrentToFieldC2  = COMET::IOARuntimeParameters::Get().GetParameterD("CalibMagnet.CurrentToFieldC2");

  fDefaultFieldType = COMET::IOARuntimeParameters::Get().GetParameterI("CalibMagnet.DefaultFieldType");

  fBFieldType = fDefaultFieldType; // default field type
  

  // Default the scale factor to 1.  This will get updated when the 
  // field map is first called.
  fScaleFactor = 1;  
  fLastNormalizationUpdate = -1;

  // For MC, use hard coded 7.12727e-05 value of magnetic field at origin (0,0,0).
  if (IsMC() ) {
    fScaleFactor = fFieldStrengthX/7.12727e-05;
  }
 
}

//********************************************************
void COMET::IOAMagneticField::SetFieldType(int type) {
//********************************************************
  fBFieldType = type;
}

//********************************************************
void COMET::IOAMagneticField::SetFieldError(double err) {
//********************************************************
  fBErr = err;
}

//********************************************************
void COMET::IOAMagneticField::SetFieldStrength(double bMax) {
//********************************************************
  fBLMax = fBLMax * (bMax / fBTMax);
  fBTMax = bMax;
  fFieldStrengthX = bMax;
}

//********************************************************
double COMET::IOAMagneticField::GetFieldStrength() {
//********************************************************
  if (fBFieldType == 1) return fBTMax;
  return fFieldStrengthX;
}

//********************************************************
void COMET::IOAMagneticField::SetFieldValue(double x, double y, double z) {
//********************************************************
  fFieldStrengthX = x;
  fFieldStrengthY = y;
  fFieldStrengthZ = z;

}

//********************************************************
TVector3 COMET::IOAMagneticField::GetFieldValue(const TLorentzVector& point) {
//********************************************************
  return GetFieldValue(point.Vect());
}

//********************************************************
TVector3 COMET::IOAMagneticField::GetFieldValue(const TVector3& point){ 
  //********************************************************
  //Using the current set field type
  return GetFieldValue(point, fBFieldType);
}

//********************************************************
TVector3 COMET::IOAMagneticField::GetFieldValue(const TLorentzVector& point, int fieldType) {
//********************************************************
  return GetFieldValue(point.Vect(), fieldType);
}

//********************************************************
TVector3 COMET::IOAMagneticField::GetFieldValue(const TVector3& point, int fieldType) {
//********************************************************
  // 0: No field at all
  if (fieldType == 0) return TVector3(0,0,0);
  // 1: The nominal field, i.e.the measured point-to-point field map for data (type 9) and the MC field for MC (type 5)
  if (fieldType == 1) return GetNominalFieldValue(point);
  // 2: Nominal field with some error
  if (fieldType == 2) {
    TVector3 nominal = GetNominalFieldValue(point);
    return nominal + TVector3(fBErr*gRandom->Uniform(-1.,1.),
			      fBErr*gRandom->Uniform(-1.,1.),
			      fBErr*gRandom->Uniform(-1.,1.));
  }
  // 3: "Perfect" B Field
  if (fieldType == 3) return TVector3(fFieldStrengthX, 0, 0);
  // 4: NOMAD's field map
  if (fieldType == 4) return GetNOMADFieldValue(point);
  // 5: The B Field currently used by SimG4 (will change)
  if (fieldType == 5) return GetMCFieldValue(point);
  // 6: Latest simulated field map (July 2008)
  if (fieldType == 6) return GetSimFieldValue(point);
  // 7: Measured field map (March 2010)
  if (fieldType == 7) return GetMeasFieldValue(point);
  // 8: Simple homogeneous B Field
  if (fieldType == 8) return TVector3(fFieldStrengthX, fFieldStrengthY, fFieldStrengthZ);
  // 9: Point to point map (official map)
  if (fieldType == 9) return GetP2PMeasFieldValue(point);
  // 111: NOT A REAL MAP
  // Identity map: B-Field component equals point coordinates
  // !!! This map is just for the purpose of checking
  if (fieldType == 111) return point;



  COMETError("ERROR in TVector3 IOAMagneticField::GetFieldValue: Wrong field type = "<<fieldType<<"! Exit!");;
  exit(1);
}

//To be called by SimG4
//********************************************************
void COMET::IOAMagneticField::GetFieldValue(const double pointIn[3], double *BField) {
//********************************************************
  TVector3 point(pointIn[0], pointIn[1], pointIn[2]);
  TVector3 fieldStrength = GetFieldValue(point);
  BField[0] = fieldStrength.X();
  BField[1] = fieldStrength.Y();
  BField[2] = fieldStrength.Z();
}


//********************************************************
void COMET::IOAMagneticField::UpdateScaleFactor() {
//********************************************************

  COMET::ICOMETEvent *event = COMET::IEventFolder::GetCurrentEvent();
  if(!event) return;
    
  ICOMETContext context = event->GetContext();

  // If this is still the last time as last check, then just return current value.
  if(context.GetTimeStamp() == fLastNormalizationUpdate) return;

  unsigned int partition = context.GetPartition();
  // Check for invalid or MC context;
  if( partition == 0xdeadbeef  || partition & COMET::ICOMETContext::kMCData) return;

  //  0.1909 T (for I = 2700 A for runs before March 21st 2010)
  //  0.1839 T (for I = 2600 A for runs since March 21st 2010)
  //  0.2048 T (for I = 2900 A for runs since November 2010)
  double I  = 2900; 


  // time corresonding to  10/03/21 00:30:00 JMT
  //  int t0=1269100800;  
  
  // This is only used for old data or when current = 0.0.
  double default_current = 0.0; 
  
  // The magnet current table was only filled after t=1264575356.
  // If we are dealing with data after that time, then we make database query.
  // If we are dealing with data before that time, set current to default.
  if(context.GetTimeStamp() >= 1264575356){
    
    // initialize the IMagnetCurrent the first time it is called.
    if (!fMagnetCurrentInitialized){
      fMagnetCurrent = new IMagnetCurrent( context.GetTimeStamp()-10000, context.GetTimeStamp()+10000);
      fMagnetCurrentInitialized = true;
      std::cout << "IMagnetCurrent initialized in CalibMagnet" << std::endl;
    }
    
    
    // Set the current unix time
    fMagnetCurrent->SetTime( context.GetTimeStamp() );
    
    if(COMET::ICOMETLog::GetLogLevel("CalibMagnet") >= COMET::ICOMETLog::VerboseLevel){
      std::cout << " Magnet Current = " << fMagnetCurrent->GetI() 
                << " A, at time = " << fMagnetCurrent->GetTimeStamp() << std::endl;
    }
    
    // Get the magnet current at that time
    I = fMagnetCurrent->GetI();
    
  }else{
    COMETLog("COMET::IOAMagneticField::UpdateCurrentScaleFactor() : time is before magnet current was being logged.\n"
             << "Using default of current = " << default_current << "A");
    I = default_current; // Use default for early times.
  }
  
  // Compute scale factor 
  I = I*100/100.089; 	//correction between set value for current and measured value fo current
  double I0 = fReferenceCurrent; // current at which field measurements were done
  double c0 = fCurrentToFieldC0;
  double c1 = fCurrentToFieldC1;
  double c2 = fCurrentToFieldC2;

  fScaleFactor=(c0 + c1*I + c2*I*I)/(c0 + c1*I0 + c2*I0*I0);
  fLastNormalizationUpdate = context.GetTimeStamp();

}

//********************************************************
TVector3 COMET::IOAMagneticField::GetNominalFieldValue(const TVector3& point) {
//********************************************************

  // This is the current field used by the reconstruction
  if ( IsMC() ) 
    return GetMCFieldValue(point);
  else
    return GetP2PMeasFieldValue(point);

}


//********************************************************
bool COMET::IOAMagneticField::IsMC() {
//********************************************************

  COMET::ICOMETEvent *event = COMET::IEventFolder::GetCurrentEvent();
  if(event){
    ICOMETContext context = event->GetContext();

    unsigned int partition = context.GetPartition();
    // Check for invalid or MC context;
    if( partition & COMET::ICOMETContext::kMCData)
      return true;
  }

  return false;

}

//*******************************************************************************
//******************************* NOMAD *****************************************


//********************************************************
TVector3 COMET::IOAMagneticField::GetNOMADFieldValue(const TVector3& point) {
//********************************************************


  double scale_factor = 1;

  if (fNormalizationFromCurrent){
    // scale factor based on the coil current 
    // scale factor that relates NOMAD and Meas fields  (alpha = 0.178318)
    UpdateScaleFactor();
    scale_factor = fScaleFactor*GetMeasFieldXValue(TVector3(0,0,0))/GetNOMADFieldXValue(TVector3(0,0,0));
  }
  else
    scale_factor = fFieldStrengthX/GetNOMADFieldXValue(TVector3(0,0,0));

  double Bx = scale_factor*GetNOMADFieldXValue(point);
  double By = scale_factor*GetNOMADFieldYValue(point);
  double Bz = scale_factor*GetNOMADFieldZValue(point);


  return TVector3(Bx, By, Bz);
  
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADFieldXValue(const TVector3& point) {
//********************************************************
  // Set the B_x = 0 in case the point is not inside any magnetic area
  double Bx = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet);

  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner x value
    Bx = GetNOMADInnerXValue(x,y,z);
  }
  // Otherwice, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
        double x1 = std::min(x,fXCenl);
        double y1 = std::min(y,fYCenl);
        double z1 = std::min(z,fZCenl);	
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
        double y2 = std::min(y,(fYMinC - fRMinC));
       double z2 = std::min(z,(fZMinC - fRMinC));
      double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_x value with rounded corners and soft transfer
      Bx = GetNOMADInnerXValue(x1,y1,z1) * soft;
      // If x within [fXCenl, fXMinM] make a soft transfer for B_x along x
      if (x > fXCenl) {
	double BxM = GetNOMADOuterXValue(fXMinM,y,z);
	double xRatio = (x - fXCenl) / (fXMinM - fXCenl);
	double softM = Soft01(xRatio);
	Bx = Bx + (BxM - Bx) * softM;
      }
    }
    // Otherwice, it is in the outer region
    else {
      Bx = GetNOMADOuterXValue(x,y,z);
    }
  }
  // Put Bx = 0 for very small field components
  if (fabs(Bx) <= fBFMin) Bx = 0;
  // Now ready to return B_x

  return Bx;
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADFieldYValue(const TVector3& point){ 
//********************************************************
  // Set the B_y = 0 in case the point is not inside any magnetic area
  double By = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet); 

  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner y value
    By = GetNOMADInnerYValue(x,y,z);
  }
  // Otherwise, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
        double x1 = std::min(x,fXCenl);
        double y1 = std::min(y,fYCenl);
        double z1 = std::min(z,fZCenl);	
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
        double y2 = std::min(y,(fYMinC - fRMinC));
        double z2 = std::min(z,(fZMinC - fRMinC));
      double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_y value with rounded corners and soft transfer
      By = GetNOMADInnerYValue(x1,y1,z1) * soft;
    }
    // Otherwise, it is in the outer region
    else {
      By = GetNOMADOuterYValue(x,y,z);
    }
  }
  // Put By = 0 for very small field components
  if (fabs(By) <= fBFMin) By = 0;
		
	//apply correct symmetry
	if (x>=0 && y<0) By = (-1)*By;
	else if (x<0 && y>=0) By = (-1)*By;
	
	
  return By;
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADFieldZValue(const TVector3& point) {
//********************************************************
  // Set the B_z = 0 in case the point is not inside any magnetic area
  double Bz = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet); 
  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner z value
    Bz = GetNOMADInnerZValue(x,y,z);
  }
  // Otherwice, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
        double x1 = std::min(x,fXCenl);
        double y1 = std::min(y,fYCenl);
        double z1 = std::min(z,fZCenl);
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
        double y2 = std::min(y,(fYMinC - fRMinC));
        double z2 = std::min(z,(fZMinC - fRMinC));
        double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_z value with rounded corners and soft transfer
      Bz = GetNOMADInnerZValue(x1,y1,z1) * soft;
      // If x within [fXCenl, fXMinM] make a soft transfer for B_z along x
      if (x > fXCenl) {
	double BzM = GetNOMADOuterZValue(fXMinM,y,z);
	double xRatio = (x - fXCenl) / (fXMinM - fXCenl);
	double softM = Soft01(xRatio);
	Bz = Bz + (BzM - Bz) * softM; 
      }
    }
    // Otherwice, it is in the outer region
    else {
      Bz = GetNOMADOuterZValue(x,y,z);
    }
  }

  // Put Bz = 0 for very small field components
  if (fabs(Bz) <= fBFMin) Bz = 0;

	//apply correct symmetry
	if (x>=0 && z<0) Bz = (-1)*Bz;
	else if (x<0 && z>=0) Bz = (-1)*Bz;
	
  return Bz;
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADInnerXValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bx within the measured area: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double Bx = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the NOMAD function is paremeterized assuming meters and kilogauss
  double x = x_in / unit::m;
  double y = y_in / unit::m;
  double z = z_in / unit::m;

  // Biases for x and y
  double x1 = x - 1.0;
  double y1 = y - 1.0;

  // **** Main component Bx(x,y)
  // If y is within some y limit 
  if (y <= 1.0) {
    // The main component for Bx follows
    Bx = 4.006 + 0.014*y + (-0.012 - 0.078*y*y) * x1*x1;
  } 
  // If y is not within that y limit
  else {
    // The main component for Bx follows
    Bx = 4.02 + 0.026*y1 + 0.44*y1*y1 - (0.090 + 0.372*y1 + 0.753*y1*y1)*x1*x1;
  }
  // After getting the main component for Bx we will now 
  // apply the correction factors; factorX2, factorX3 and factorX4
  // ***** factorX2(x,y) at |x|>1.3
  if (x > 1.3) {
    double dZ = 0.956;
    double z1 = z/dZ;
    if (z1 > 2.0) z1 = 2.0 + (z1-2.0)*1.06;
    int z2 = std::min(int(z1),3);
    double deltaZ = fabs(z1 - z2);
    double factorX2A = 1.0/(1.0 + (78.23*(x-1.76))*(78.23*(x-1.76))) * (x-1.30)/0.46;
    double factorX2B = 0.108;
    double factorX2C = std::min((0.1 + 7.780*(x-1.76)*(x-1.76)),0.5);
    if (z2 == 3) {
      factorX2A = factorX2A*2;
      factorX2B = factorX2B / (factorX2B + 0.5);      
    }
    double factorX2 = 1.0 - factorX2A * (1.0 - (deltaZ/factorX2B))/(1+(deltaZ/factorX2C)*(deltaZ/factorX2C));
    Bx = Bx * factorX2;
  }
  // ***** factorX3(x,z) at |z|>2.0
  if (z > 2.0) {
    double z2 = z - 2.0;
    double z3 = z2 - 1.05;
    double factorX3A = std::max(1.0 + 0.010*z2, 0.984 + 0.0314*z2);
    double z3test = std::max(0.0, z3);
    double factorX3B = 0.5 * z3test*z3test;
    double factorX3C = 0;
    if (x1 <= 0) {
      factorX3C = 1.0 + 0.00131*z2 - 0.0313*z2*z2*z2;
    } else {
      double factorX3D = 0;
      if (z2 < 0.9) factorX3C = 1.0 - 0.0522*z2 + 0.0578*z2*z2;
      if (z2 > 0.8) factorX3D = 1.056 - 1.824*z3*z3 + 10.47*z3*z3*z3*z3;
      factorX3C = std::max(factorX3C,factorX3D);
    }
    double x1test = std::max(fabs(x1/0.66),0.001);
    double xto8 = pow(x1test, 8);
    double factorX3 = (factorX3A + factorX3B*fabs(x1) + factorX3C*xto8)/(1.0 + xto8);
    Bx = Bx * factorX3;
  }
  // ***** factorX4(x,y,z) at |z|>3.25 |y|<0.65
  if ((z > 3.25)&&(y < 0.65)) {
    double z3 = z - 3.25;
    double factorX4AA = pow(33.6*z3,3)/(1.0+pow(83.0*z3,3));
    double yTest = std::min(2.58*(0.65-y), 1.0);
    double factorX4A = (factorX4AA-0.052*y)*yTest;
    double factorX4B = std::max(-0.0130+2.21*factorX4A,0.5*factorX4A);
    double factorX4F = 1.06+2.0*z3;
    double factorX4FTest = std::max(fabs(factorX4F*x),0.001);
    double factorX4 = 1.0 - (factorX4A - factorX4B*x)/(1.0 + pow(factorX4FTest,8));
    Bx = Bx * factorX4;
  }

  return Bx * unit::kilogauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADInnerYValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for By within the measured area: |x|<1.66m, |y|<1.52m, |z|<3.3m
  double By = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the NOMAD function is paremeterized assuming meters and kilogauss
  double x = x_in / unit::m;
  double y = y_in / unit::m;
  double z = z_in / unit::m;

  // Biases for x and z
  double x1 = x - 1.0;
  double z2 = z - 2.5;

  // **** Main component By(x,y,z)
  double ByB1 = - 0.5;
  if (z2 > 0) ByB1 = ByB1 - 2.48*z2*z2;
  double ByB2 = 0.662 - 1.606*ByB1;
  double ByB3 = -0.447 + 0.657*ByB1;
  double ByB = ByB1*y + ByB2*y*y + ByB3*y*y*y;
  if (y > 1.0) ByB = std::max(ByB,0.0);
  double ByC1 = std::max(-0.28,-1.375 + 0.395*z);
  double zLarge = std::max(z, 2.5);
  double ByC2 = 7.34 - 3.5*zLarge + 0.7*zLarge*zLarge;
  double ByC3 = std::max(0.81,0.158 + 0.235*z);
  double yR = y;
  if ((z > 3.2)&&(y < 0.6)) {
    yR = 0.6 - (1.0-11.5*(z-3.2))*(0.6-y);
  } 
  double ByC = y*(ByC1 + ByC2*(yR-ByC3)*(yR-ByC3));
  double ByD = -ByB - ByC;
  if (x1 < 0) ByC = -ByC;
  By = ByB*x1 + ByC*x1*x1 + ByD*x1*x1*x1;
  if (x1 > 0) By = By*0.41*y;

  return By * unit::kilogauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADInnerZValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bz within the measured area: |x|<1.66m, |y|<1.52m, |z|<3.3m
  double Bz = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the NOMAD function is paremeterized assuming meters and kilogauss
  double x = x_in / unit::m;
  double y = y_in / unit::m;
  double z = z_in / unit::m;

  // **** Bz for 0<|x|<1.1
  if (x - 1.1 < 0) {
    // for |z|>2.4
    if (z - 2.4 > 0) {
      double z2 = z - 2.4;
      double BzB2 = std::max(0.40 - 0.37*y,0.17 + 0.093*y);
      double BzB3 = 1.10 - 0.09*y;
      double BzB = 0;
      if ((x > BzB2)&&(x < BzB3)) {
	double BzB1 = (-0.72 + 0.10*y)*(1.0 + (0.43-0.38*y)*x*x); 
	BzB = BzB1*(x-BzB2)*(BzB3-x);
      }
      double BzD3 = 1.10 - 0.067*y;
      double BzD = 0;
      if (x < BzD3) {
          double BzD1 = std::max(4.85 - 1.87*y,3.49);
          double BzD2 = std::max(0.91 - 0.59*y,0.61);
	BzD = x*(BzD3-x)*BzD1*(x-BzD2);
	if (x > BzD2) BzD = BzD*(BzD3-x)/(BzD3-BzD2);
      }
      Bz += BzB*z2 + BzD*z2*z2*z2;
    }
  } 
  // **** Bz for |x|>1.1 
  else { 
    // **** Bz for |x|>1.4 (excluding the last gap)
    if (x - 1.4 > 0) {
      double dZ = 0.956;
      double z1 = int(z/dZ)*dZ;
      if (z1 < 2.5) {
	double deltaZ = z - z1;
	double BzA = pow((3.41*(x-1.4)),3);
	double BzB = 3.625 + 5.440*BzA;
	Bz += BzA*deltaZ/(1.0 + pow((BzB*deltaZ),4)) * (1.0 - 2*fabs(deltaZ)/dZ);
      }
    }
    // **** Bz for |z|>2.2
    if (z - 2.2 > 0) {
      // "background"
      double z2 = z - 2.2;
      double BzB = 0.5*(x-1.1)*(1.75-x);
      double BzD = 2.35*(x-1.1)*(1.75-x)*(x-1.39);
      Bz += BzB*z2 + BzD*z2*z2*z2;
      // "peak"
      if (x - 1.3 > 0) {
	double xs = x - 1.3;
	double zsp = 0.022*xs + 2.91*xs*xs*xs;
	double zst = 4.3 + 14.7*xs;
	Bz += zsp/(1.0 + pow((zst*(z-2.94)),4));
      }
    }
  }
  return Bz * unit::kilogauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADOuterXValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  double Bx = 0;

  if ((x > fXMaxM)||(y > fYMaxM)||(z > fZMaxM)) return Bx;
  if ((x < fXMinM)&&(y < fYMinM)) return Bx;

  double z0 = int(z/fDZCees)*fDZCees;
  double z0Limit = 2.5 * unit::m;
  if (z0 > z0Limit) z0 += fHSGap3;
  if (fabs(z-z0) <= fHZGaps) return Bx;
  
  if (y < fYMinM) {
    // Variable field area (vertical bars of C's)
    double x1 = fXMaxM - x;
    Bx = fBTMax*x1 / (fXMaxM - fXMinM);
    double rY = (fYMinM - y)/0.2;
    double rYLimit = 1.0 * unit::m;
    if (rY < rYLimit) Bx = Bx*3*rY*rY / (rYLimit + 2*rY*rY*rY);
  } else {
    // Parallel field lines area (horiz. bars of C's or angles)
    Bx = -fBLMax;
    if ((y > fYMinM)&&(x > fXMinM)) {
      // Angular area, B must be rotated
      double x1 = x - fXMinM;
      double y1 = y - fYMinM;
      double rr = sqrt(x1*x1 + y1*y1);
      if (rr < fYMaxM-fYMinM) {
	Bx = -fBLMax * y1/rr;

      } else {
	Bx = 0;
      }
    }
  }  
  return Bx;
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADOuterYValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  double By = 0;

  if ((x > fXMaxM)||(y > fYMaxM)||(z > fZMaxM)) return By;
  if ((x < fXMinM)&&(y < fYMinM)) return By;

  double z0 = int(z/fDZCees)*fDZCees;
  double z0Limit = 2.5 * unit::m;
  if (z0 > z0Limit) z0 += fHSGap3;
  if (fabs(z-z0) <= fHZGaps) return By;
  
  if (y < fYMinM) {
    // Variable field area (vertical bars of C's)
    By = fBLMax* y / fYMinM;
  } else {
    // Parallel field lines area (horiz. bars of C's or angles)
    if ((y > fYMinM)&&(x > fXMinM)) {
      // Angular area, B must be rotated
      double x1 = x - fXMinM;
      double y1 = y - fYMinM;
      double rr = sqrt(x1*x1 + y1*y1);
      if (rr < fYMaxM-fYMinM) {
	By = -fBLMax * x1/rr;
      } 
    }
  }  
  return By;
}

//********************************************************
double COMET::IOAMagneticField::GetNOMADOuterZValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}

//********************************************************
double COMET::IOAMagneticField::Soft01(double x) {
//********************************************************
  // Provides "soft jump" from 0 to 1 for x within [0,1]
  double soft01 = 0;
  if (x <= 0.0) {
    soft01 = 0.0;
  } else if (x >= 1.0) {
    soft01 = 1.0;
  } else {
    if (x <= 0.5) {
      double x1 = 3.4641*x;
      soft01 = 0.385*x1*x1*x1/(1.0 + x1*x1);
    } else {
      double x1 = 3.4641*(1.0-x);
      soft01 = 1.0 - 0.385*x1*x1*x1/(1.0 + x1*x1);
    }
  }
  return soft01;
}


//***********************************************************************
//****************** MC *************************************************

//********************************************************
TVector3 COMET::IOAMagneticField::GetMCFieldValue(const TVector3& pos) {
//********************************************************
  
  TVector3 BField(0,0,0);

   if (fabs(pos.X()) > fYokeOuterWidth/2.) return BField;
   if (fabs(pos.Y()) > fYokeOuterHeight/2.) return BField;
   if (fabs(pos.Z()) > fYokeOuterLength/2.) return BField;

   // This is the field strength we use for MC.
   // As backup we use the value from parameters file...
   double mcFieldStrengthX = fFieldStrengthX;
   //std::cout << " Field Strenght " << mcFieldStrengthX << std::endl;

   // ... but we really want to try to get the value from the MC header.
   ICOMETEvent* event = COMET::IEventFolder::GetCurrentEvent();
   if (event) {
     COMET::IHandle<COMET::ISIMG4Header> mcHeader = 
       event->Get<COMET::ISIMG4Header>("truth/mcHeader");
     if(mcHeader){
       mcFieldStrengthX = mcHeader->GetOffAxisField();
      }
   }
  
  if (fabs(pos.X()) < fCoilInnerWidth/2.
      && fabs(pos.Y()) < fCoilInnerHeight/2. 
      && fabs(pos.Z()) < fCoilInnerLength/2.) {
    // Inside the inner volume of the magnet.  
    // This will eventually come from a measured field map.
    BField.SetX(mcFieldStrengthX);
  } else if (fabs(pos.X()) < fCoilOuterWidth/2.
	     && fabs(pos.Y()) < fCoilOuterHeight/2. 
	     && fabs(pos.Z()) < fCoilOuterLength/2.) {
    // Inside the coil volume of the magnet.
    double Y = (fabs(pos.Y()) - fCoilInnerHeight/2.)
      / (fCoilOuterHeight/2. - fCoilInnerHeight/2.);
    double Z = (fabs(pos.z()) - fCoilInnerLength/2.)
      / (fCoilOuterLength/2. - fCoilInnerLength/2.);
    double depth = std::max(Y,Z);
    if (depth<0.0) depth = 0.0;
    if (depth>1.0) depth = 1.0;
    BField.SetX(mcFieldStrengthX * (1.0 - depth));
  }

  //THIS IS IMPLEMENTED IN MC BUT NEVER(!) CALLED DUE TO A #ifdef REST_OF_FIELD
//   } else if (fabs(pos.X()) < fInnerVolumeWidth/2.) {
//     // In the top or bottom part of the flux return
//     double innerArea = fInnerVolumeHeight*fInnerVolumeLength;
//     double fluxArea = (fMagnetHeight-fInnerVolumeHeight)*fMagnetLength;
//     double field = fFieldStrengthX*innerArea/fluxArea;
//     BField.SetX(-field);
//   } else if (fabs(pos.Y()) < fInnerVolumeHeight/2.) {
//     // In the side part of the flux return.
//     double innerArea = fInnerVolumeHeight*fInnerVolumeLength;
//     double fluxArea = (fMagnetHeight-fInnerVolumeHeight)*fMagnetLength;
//     double field = fFieldStrengthX*innerArea/fluxArea;
//     if (pos.X() < 0) field = -field;
//     BField.SetY(field*pos.Y()*fInnerVolumeHeight);
//   } else {
//     // In the corners of the flux return.
//     double innerArea = fInnerVolumeHeight*fInnerVolumeLength;
//     double fluxArea = (fMagnetHeight-fInnerVolumeHeight)*fMagnetLength;
//     double field = fFieldStrengthX*innerArea/fluxArea/1.414;
//     BField.SetX(-field);
//     BField.SetY(field);
//     if (pos.X() < 0.0) BField.SetY( - BField.Y() );
//   }

  return BField;
}

//***********************************************************************
//****************** SIMULATION JULY 2008********************************

//********************************************************
TVector3 COMET::IOAMagneticField::GetSimFieldValue(const TVector3& point) {
//********************************************************


  double scale_factor = 1;

  if (fNormalizationFromCurrent){
    // scale factor based on the coil current 
    // scale factor that relates Sim and Meas fields  (alpha = -0.350727)
    UpdateScaleFactor();
    scale_factor = fScaleFactor*GetMeasFieldXValue(TVector3(0,0,0))/GetSimFieldXValue(TVector3(0,0,0));
  }
  else
    scale_factor = fFieldStrengthX/GetSimFieldXValue(TVector3(0,0,0));


  double Bx = scale_factor*GetSimFieldXValue(point);
  double By = scale_factor*GetSimFieldYValue(point);
  double Bz = scale_factor*GetSimFieldZValue(point);

  return TVector3(Bx, By, Bz);
  
}

//********************************************************
double COMET::IOAMagneticField::GetSimFieldXValue(const TVector3& point) {
//********************************************************
  // Set the B_x = 0 in case the point is not inside any magnetic area
  double Bx = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet);

  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner x value
    Bx = GetSimInnerXValue(x,y,z);
  }
  // Otherwise, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
      double x1 = std::min(x,fXCenl);
      double y1 = std::min(y,fYCenl);
      double z1 = std::min(z,fZCenl);	
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
      double y2 = std::min(y,(fYMinC - fRMinC));
      double z2 = std::min(z,(fZMinC - fRMinC));
      double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_x value with rounded corners and soft transfer
      Bx = GetSimInnerXValue(x1,y1,z1) * soft;
      // If x within [fXCenl, fXMinM] make a soft transfer for B_x along x
      if (x > fXCenl) {
	double BxM = GetSimOuterXValue(fXMinM,y,z);
	double xRatio = (x - fXCenl) / (fXMinM - fXCenl);
	double softM = Soft01(xRatio);
	Bx = Bx + (BxM - Bx) * softM;
      }
    }
    // Otherwise, it is in the outer region
    else {
      Bx = GetSimOuterXValue(x,y,z);
    }
  }
  // Put Bx = 0 for very small field components
  if (fabs(Bx) <= fBFMin) Bx = 0;
  // Now ready to return B_x

  return Bx;
}

//********************************************************
double COMET::IOAMagneticField::GetSimFieldYValue(const TVector3& point) {
//********************************************************
  // Set the B_y = 0 in case the point is not inside any magnetic area
  double By = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet); 

  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner y value
    By = GetSimInnerYValue(x,y,z);
  }
  // Otherwise, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
        double x1 = std::min(x,fXCenl);
        double y1 = std::min(y,fYCenl);
        double z1 = std::min(z,fZCenl);	
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
        double y2 = std::min(y,(fYMinC - fRMinC));
        double z2 = std::min(z,(fZMinC - fRMinC));
      double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_y value with rounded corners and soft transfer
      By = GetSimInnerYValue(x1,y1,z1) * soft;
    }
    // Otherwise, it is in the outer region
    else {
      By = GetSimOuterYValue(x,y,z);
    }
  }
  // Put By = 0 for very small field components
  if (fabs(By) <= fBFMin) By = 0;
  // Now ready to return B_y
  return By;
}

//********************************************************
double COMET::IOAMagneticField::GetSimFieldZValue(const TVector3& point) {
//********************************************************
  // Set the B_z = 0 in case the point is not inside any magnetic area
  double Bz = 0;
  // Get into the "magnet frame"
  double x = fabs(point.X());
  double y = fabs(point.Y());
  double z = fabs(point.Z() - fZMagnet); 
  // Check if the point is in the central area
  if ((x <= fXCenl)&&(y <= fYCenl)&&(z <= fZCenl)) {
    // Then return the inner z value
    Bz = GetSimInnerZValue(x,y,z);
  }
  // Otherwise, check if the point is within the z-limit  
  else if (z < fZMaxM) {
    // Check if the point is between the central and magnetic area,
    // and in the "coil frame" in y and z
    if ((x < fXMinM)&&(y <= fYMaxC)&&(z <= fZMaxC)) {
      // Get a reference point to the central area
        double x1 = std::min(x,fXCenl);
        double y1 = std::min(y,fYCenl);
        double z1 = std::min(z,fZCenl);
      // Take care of rounded corners and soft transfer to B=0 from 
      // the internal surface of the coil core to the external one
        double y2 = std::min(y,(fYMinC - fRMinC));
        double z2 = std::min(z,(fZMinC - fRMinC));
        double r = sqrt((y-y2)*(y-y2) + (z-z2)*(z-z2));
      double rRatio = (fRMaxC - r) / (fRMaxC - fRMinC);
      double soft = Soft01(rRatio);
      // Get the B_z value with rounded corners and soft transfer
      Bz = GetSimInnerZValue(x1,y1,z1) * soft;
      // If x within [fXCenl, fXMinM] make a soft transfer for B_z along x
      if (x > fXCenl) {
	double BzM = GetSimOuterZValue(fXMinM,y,z);
	double xRatio = (x - fXCenl) / (fXMinM - fXCenl);
	double softM = Soft01(xRatio);
	Bz = Bz + (BzM - Bz) * softM; 
      }
    }
    // Otherwise, it is in the outer region
    else {
      Bz = GetSimOuterZValue(x,y,z);
    }
  }

  // Put Bz = 0 for very small field components
  if (fabs(Bz) <= fBFMin) Bz = 0;
  // Now ready to return B_z
  return Bz;
}

//********************************************************
double COMET::IOAMagneticField::GetSimInnerXValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bx within the measured area: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double Bx = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::cm;
  double y = y_in / unit::cm;
  double z = z_in / unit::cm;
  // Introduce auxiliary variable for absolute value of coordinates
  double xa = fabs(x);
  double ya = fabs(y);
  double za = fabs(z);

  // The parametrised area is divided into 4 regions, for each a 3rd order polynomial 
  // description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;

  if (fabs(z) <= 285) {
    if (fabs(y) <= 125) {
      // region 1: |x|<=1.25m, |y|<=1.25m, |z|<=2.85m (contains all TPC chambers)
      // statiscial error: 1.3 Gauss
      if (fabs(x) <= 125) {
	a0 = -2.0306e+03;
	a1 = 1.0269e-01; 
	a2 = -8.0973e-02;
	a3 = 4.2998e-02;
	a4 = -2.8874e-03; 
	a5 = 1.9941e-03;
	a6 = -1.9060e-03;
	a7 = 2.6655e-03; 
	a8 = -4.0723e-05;
	a9 = -2.4769e-06;
	a10 =7.9008e-06; 
	a11 = -3.0372e-06;
	a12 = 1.0050e-05;
	a13 = -4.4785e-05; 
	a14 = -2.8096e-06;
	a15 = 5.8426e-06;
	a16 = 6.3161e-06; 
	a17 = 7.8780e-07;
	a18 = 7.4948e-07;
	a19 = -1.5385e-06; 
      }
      // region 2: 1.25m<|x|<=fXCenl, |y|<=1.25m, |z|<=2.85m
      // statiscial error: 14.8 Gauss
      else {
	a0 = -8.7223e+03;
	a1 = 1.4493e+02; 
	a2 = 3.2292e-01;
	a3 = -7.3465e-01;
	a4 = -1.0446e+00; 
	a5 = -7.2493e-03;
	a6 = 1.6079e-02;
	a7 = 6.6346e-03; 
	a8 = 8.4244e-04;
	a9 = -4.3596e-03;
	a10 = 2.4986e-03; 
	a11 = 3.3442e-05;
	a12 = -7.6154e-05;
	a13 = -5.3370e-05; 
	a14 = -3.5444e-06;
	a15 = 2.7115e-05;
	a16 = -8.1217e-06; 
	a17 = -1.5485e-06;
	a18 = -1.2106e-06;
	a19 = 3.1017e-06; 
      }
    }
    // region 3: |x|<=fXCenl, 1.25m<|y|<=fYCenl, |z|<=2.85m
    // statiscial error: 8.1 Gauss
    else {
     a0 = -1.5516e+03;
     a1 = -5.1688e+00; 
     a2 = -8.6354e+00;
     a3 = 3.6899e-01;
     a4 = -1.1904e-02; 
     a5 = 9.2251e-02;
     a6 = -7.8398e-04;
     a7 = 4.6893e-02; 
     a8 = -3.5014e-03;
     a9 = -9.7560e-04;
     a10 = -4.3896e-05; 
     a11 = 1.6237e-04;
     a12 = 6.5866e-06;
     a13 = -4.7779e-04; 
     a14 = -7.4334e-06;
     a15 = 5.1944e-06;
     a16 = -3.9592e-05; 
     a17 = 9.0668e-06;
     a18 = 5.9543e-06;
     a19 = -6.4976e-07; 
    }
  } 
  // region 4 |x|<=fXCenl, |y|<=fXCenl, 2.85m<|z|<=fZCenl
  // statiscial error: 23.5 Gauss
  else {
    a0 = 0;
    a1 = 3.6599e-01; 
    a2 = -5.0177e-02;
    a3 = -1.8848e+01;
    a4 = -8.5641e-02; 
    a5 = 1.2829e-02;
    a6 = 3.2781e-02;
    a7 = 3.8802e-03; 
    a8 = -3.3750e-03;
    a9 = 5.6112e-02;
    a10 = -1.9709e-04; 
    a11 = -1.2911e-05;
    a12 = 4.5583e-04;
    a13 = -5.8187e-05; 
    a14 = -3.1428e-05;
    a15 = -1.5012e-04;
    a16 = 7.4747e-06; 
    a17 = 1.4498e-07;
    a18 = 1.0256e-05;
    a19 = -5.2350e-05; 
  }
  Bx = a0 + a1*xa + a2*ya + a3*za + a4*xa*xa + a5*xa*ya + a6*xa*za + a7*ya*ya + a8*ya*za + a9*za*za + a10*xa*xa*xa + a11*xa*xa*ya + a12*xa*xa*za + a13*xa*ya*ya + a14*xa*ya*za + a15*xa*za*za + a16*ya*ya*ya + a17*ya*ya*za + a18*ya*za*za + a19*za*za*za;
  return Bx * unit::gauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetSimInnerYValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bx within the measured area: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double By = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::cm;
  double y = y_in / unit::cm;
  double z = z_in / unit::cm;
  // Introduce auxiliary variable for absolute value of coordinates
  double xa = fabs(x);
  double ya = fabs(y);
  double za = fabs(z);

  // The parametrised area is divided into 4 regions, for each a 3rd order polynomial 
  // description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;

  if (fabs(z) <= 285) {
    if (fabs(y) <= 125) {
      // region 1: |x|<=1.25m, |y|<=1.25m, |z|<=2.85m (contains all TPC chambers)
      // statiscial error: 0.8 Gauss
      if (fabs(x) <= 125) {
	a0 = -3.1808e-01;
	a1 = -2.4453e-02; 
	a2 = 1.3890e-02;
	a3 = 2.9351e-02;
	a4 = -1.0961e-04; 
	a5 = 5.9106e-03;
	a6 = -7.2906e-05;
	a7 = -1.4767e-03; 
	a8 = -2.9495e-04;
	a9 = -2.0939e-04;
	a10 = 3.8531e-06; 
	a11 = -4.3932e-05;
	a12 = -1.4008e-06;
	a13 = 1.1395e-05; 
	a14 = 1.7887e-06;
	a15 = 8.5893e-07;
	a16 = 1.4573e-05; 
	a17 = 2.3757e-08;
	a18 = 1.1761e-06;
	a19 = 3.6597e-07; 
      }
      // region 2: 1.25m<|x|<=fXCenl, |y|<=1.25m, |z|<=2.85m
      // statiscial error: 0.4 Gauss
      else {
	a0 = -2.2174e+02;
	a1 = 4.9122e+00; 
	a2 = -9.5611e-01;
	a3 = -7.5003e-02;
	a4 = -3.5887e-02; 
	a5 = 1.5746e-02;
	a6 = 9.8987e-04;
	a7 = 5.0517e-03; 
	a8 = 2.2612e-04;
	a9 = 1.6978e-06;
	a10 = 8.6123e-05; 
	a11 = -5.6536e-05;
	a12 = -1.9365e-06;
	a13 = -4.6246e-05; 
	a14 = -3.3609e-06;
	a15 = -1.5145e-06;
	a16 = 1.7547e-05; 
	a17 = 5.2858e-08;
	a18 = 1.3548e-06;
	a19 = 4.7128e-07; 
      }
    }
    // region 3: |x|<=fXCenl, 1.25m<|y|<=fYCenl, |z|<=2.85m
    // statiscial error: 5.8 Gauss
    else {
     a0 = -1.8378e+03;
     a1 = 1.8415e+00; 
     a2 = 3.9807e+01;
     a3 = -5.5344e-02;
     a4 = 1.2887e-02; 
     a5 = -3.1304e-02;
     a6 = 3.1348e-04;
     a7 = -2.8777e-01; 
     a8 = 7.5134e-04;
     a9 = -1.7429e-04;
     a10 = 1.4881e-05; 
     a11 = -1.6263e-04;
     a12 = -3.4890e-06;
     a13 = 1.9524e-04; 
     a14 = -4.0578e-07;
     a15 = 1.4098e-06;
     a16 = 6.9818e-04; 
     a17 = -1.3875e-06;
     a18 = -8.6400e-07;
     a19 = 8.1114e-07; 
    }
  } 
  // region 4 |x|<=fXCenl, |y|<=fXCenl, 2.85m<|z|<=fZCenl
  // statiscial error: 1.8 Gauss
  else {
    a0 = 0;
    a1 = 5.8288e-01; 
    a2 = 1.4075e+00;
    a3 = -2.3401e-01;
    a4 = 5.9064e-03; 
    a5 = 2.3902e-03;
    a6 = -6.2787e-03;
    a7 = -2.9369e-03; 
    a8 = -8.8831e-03;
    a9 = 1.6506e-03;
    a10 = -6.3093e-06; 
    a11 = -5.9813e-05;
    a12 = -1.2869e-05;
    a13 = 9.3497e-06; 
    a14 = 2.2935e-05;
    a15 = 1.2548e-05;
    a16 = 1.6539e-05; 
    a17 = 4.3073e-06;
    a18 = 1.3650e-05;
    a19 = -2.8003e-06; 
  }
  By = a0 + a1*xa + a2*ya + a3*za + a4*xa*xa + a5*xa*ya + a6*xa*za + a7*ya*ya + a8*ya*za + a9*za*za + a10*xa*xa*xa + a11*xa*xa*ya + a12*xa*xa*za + a13*xa*ya*ya + a14*xa*ya*za + a15*xa*za*za + a16*ya*ya*ya + a17*ya*ya*za + a18*ya*za*za + a19*za*za*za;
  //if ((x<0 && y>0)||(x>0 && y<0)) By=-By;
  if (x*y<0) By=-By;
  return By * unit::gauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetSimInnerZValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bx within the measured area: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double Bz = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::cm;
  double y = y_in / unit::cm;
  double z = z_in / unit::cm;
  // Introduce auxiliary variable for absolute value of coordinates
  double xa = fabs(x);
  double ya = fabs(y);
  double za = fabs(z);
  
  // The parametrised area is divided into 4 regions, for each a 3rd order polynomial 
  // description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;

  if (fabs(z) <= 285) {
    if (fabs(y) <= 125) {
      // region 1: |x|<=1.25m, |y|<=1.25m, |z|<=2.85m (contains all TPC chambers)
      // statiscial error: 1.2 Gauss
      if (fabs(x) <= 125) {
	a0 = -3.6026e-01;
	a1 = 7.2131e-02; 
	a2 = 2.9800e-02;
	a3 = -7.6175e-02;
	a4 = -1.5651e-03; 
	a5 = -7.8254e-05;
	a6 = 3.0404e-04;
	a7 = -1.7688e-04; 
	a8 = -4.5884e-04;
	a9 = 8.5683e-04;
	a10 = 6.9746e-06; 
	a11 = -1.7408e-06;
	a12 = 6.2436e-06;
	a13 = 1.2929e-06; 
	a14 = 2.0103e-06;
	a15 = -6.9822e-06;
	a16 = 4.1907e-08; 
	a17 = 1.3674e-06;
	a18 = 1.1647e-06;
	a19 = -2.1283e-06; 
      }
      // region 2: 1.25m<|x|<=fXCenl, |y|<=1.25m, |z|<=2.85m
      // statiscial error: 22.5 Gauss
      else {
	a0 = -9.7738e+02;
	a1 = 1.5409e+01; 
	a2 = -5.7830e-02;
	a3 = 6.2427e+00;
	a4 = -6.9867e-02; 
	a5 = 8.7237e-04;
	a6 = -8.6529e-02;
	a7 = 1.1089e-04; 
	a8 = -1.7406e-04;
	a9 = -2.6958e-03;
	a10 = 6.8960e-05; 
	a11 = -1.8545e-06;
	a12 = 3.0188e-04;
	a13 = -2.0529e-06; 
	a14 = -2.2860e-06;
	a15 = 1.7106e-05;
	a16 = 9.0275e-08; 
	a17 = 1.9105e-06;
	a18 = 1.5764e-06;
	a19 = -6.5951e-07; 
      }
    }
    // region 3: |x|<=fXCenl, 1.25m<|y|<=fYCenl, |z|<=2.85m
    // statiscial error: 8.9 Gauss
    else {
     a0 = 3.1393e+01;
     a1 = 4.3520e-01; 
     a2 = -5.6767e-01;
     a3 = -1.9558e-01;
     a4 = -9.7754e-04; 
     a5 = -2.4753e-03;
     a6 = -3.3582e-03;
     a7 = 1.5249e-03; 
     a8 = 5.1630e-03;
     a9 = -4.9225e-04;
     a10 = 4.1915e-06; 
     a11 = -6.6473e-06;
     a12 = 9.7761e-06;
     a13 = 2.3833e-06; 
     a14 = 2.1119e-05;
     a15 = -4.0218e-06;
     a16 = 6.3975e-06; 
     a17 = -2.9604e-05;
     a18 = 7.2754e-06;
     a19 = -1.1295e-06; 
    }
  } 
  // region 4 |x|<=fXCenl, |y|<=fXCenl, 2.85m<|z|<=fZCenl
  // statiscial error: 25.4 Gauss
  else {
    a0 = 0;
    a1 = -2.0704e+01; 
    a2 = 9.1026e-01;
    a3 = 7.3106e+00;
    a4 = -1.2398e-01; 
    a5 = -3.5659e-03;
    a6 = 2.0083e-01;
    a7 = -1.1760e-03; 
    a8 = -5.6030e-03;
    a9 = -5.3729e-02;
    a10 = 5.6690e-05; 
    a11 = -1.8282e-05;
    a12 = 3.8882e-04;
    a13 = 4.4144e-06; 
    a14 = 2.1364e-05;
    a15 = -4.4458e-04;
    a16 = -3.0552e-06; 
    a17 = 6.9573e-06;
    a18 = 7.3851e-06;
    a19 = 9.8208e-05; 
  }
  Bz = a0 + a1*xa + a2*ya + a3*za + a4*xa*xa + a5*xa*ya + a6*xa*za + a7*ya*ya + a8*ya*za + a9*za*za + a10*xa*xa*xa + a11*xa*xa*ya + a12*xa*xa*za + a13*xa*ya*ya + a14*xa*ya*za + a15*xa*za*za + a16*ya*ya*ya + a17*ya*ya*za + a18*ya*za*za + a19*za*za*za;
  //if ((x<0 && z>0)||(x>0 && z<0)) Bz=-Bz;
  if (x*z<0) Bz=-Bz;
  return Bz * unit::gauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetSimOuterXValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}

//********************************************************
double COMET::IOAMagneticField::GetSimOuterYValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}

//********************************************************
double COMET::IOAMagneticField::GetSimOuterZValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}


//***********************************************************************
//****************** MEASURED MAP MARCH 2010 ****************************

// Temporary implementation (bad programming style)
// Also with bad transitions between the different regions
// Moreover it is using polynomial fits instead of more elaborated methods
// Not measured regions are covered with the NOMAD field so far


//********************************************************
TVector3 COMET::IOAMagneticField::GetMeasFieldValue(const TVector3& point) {
//********************************************************

  double scale_factor = 1;;

  if (fNormalizationFromCurrent){
    // scale factor based on the coil current 
    UpdateScaleFactor();
    scale_factor = fScaleFactor;
  }
  else
    scale_factor = fFieldStrengthX/GetMeasFieldXValue(TVector3(0,0,0));


  double Bx = scale_factor*GetMeasFieldXValue(point);
  double By = scale_factor*GetMeasFieldYValue(point);
  double Bz = scale_factor*GetMeasFieldZValue(point);

  return TVector3(Bx, By, Bz);

}

//********************************************************
double COMET::IOAMagneticField::GetMeasFieldXValue(const TVector3& point) {
//********************************************************
  // Set the B_x = 0 in case the point is not inside any magnetic area
  double Bx = 0;
  // Get into the "magnet frame"
  double x = point.X();
  double y = point.Y();
  double z = point.Z() - fZMagnet;

  // Check if the point is in the central area
  if ((fabs(x) <= fXMaxTPC) && (fabs(y) <= fYMaxTPC) && (z <= fZMaxTPC) && (z >= fZMinTPC)) {
    // Then return the inner x value
    Bx = GetMeasInnerXValue(x,y,z);
  }
  // Otherwise get normalized NOMAD value
  else {
    Bx = GetNOMADFieldXValue(point) / GetNOMADInnerXValue(0,0,0) * GetMeasInnerXValue(0,0,0);
  }

  return Bx;
}

//********************************************************
double COMET::IOAMagneticField::GetMeasFieldYValue(const TVector3& point) {
//********************************************************
  // Set the B_y = 0 in case the point is not inside any magnetic area
  double By = 0;
  // Get into the "magnet frame"
  double x = point.X();
  double y = point.Y();
  double z = point.Z() - fZMagnet;

  // Check if the point is in the central area
  if ((fabs(x) <= fXMaxTPC) && (fabs(y) <= fYMaxTPC) && (z <= fZMaxTPC) && (z >= fZMinTPC)) {
    // Then return the inner x value
    By = GetMeasInnerYValue(x,y,z);
  }
  // Otherwise get normalized NOMAD value
  else {
    By = GetNOMADFieldYValue(point) / GetNOMADInnerXValue(0,0,0) * GetMeasInnerXValue(0,0,0);
    if (x*y > 0) By=-By;
  }

  return By;
}

//********************************************************
double COMET::IOAMagneticField::GetMeasFieldZValue(const TVector3& point) {
//********************************************************
  // Set the B_z = 0 in case the point is not inside any magnetic area
  double Bz = 0;
  // Get into the "magnet frame"
  double x = point.X();
  double y = point.Y();
  double z = point.Z() - fZMagnet;

  // Check if the point is in the central area
  if ((fabs(x) <= fXMaxTPC) && (fabs(y) <= fYMaxTPC) && (z <= fZMaxTPC) && (z >= fZMinTPC)) {
    // Then return the inner x value
    Bz = GetMeasInnerZValue(x,y,z);
  }
  // Otherwise get normalized NOMAD value
  else {
    Bz = GetNOMADFieldZValue(point) / GetNOMADInnerXValue(0,0,0) * GetMeasInnerXValue(0,0,0);
    if (x*z > 0) Bz=-Bz;
  }

  return Bz;
}

//********************************************************
double COMET::IOAMagneticField::GetMeasInnerXValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bx within the measured area of NOMAD: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double Bx = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::mm;
  double y = y_in / unit::mm;
  double z = z_in / unit::mm;

  // The parametrised area is divided into 4 regions (TPC1, TPC2, TPC3 and the rest)
  // For each region a 3rd order polynomial description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;


  if (-1120 <= x && x <= 1120) { 
    if (0 <= y && y <= 1180) {//FOR THE MOMENT ASSUME SYMMETRY IN Y UNTIL FIT METHOD HAS BEEN IMPROVED (so that arm 3 can also be included)
      //region of TPC1: |x|<=1.12m, |y|<=1.18m, -0.975m<z<=0.325m
      if (-975 < z && z <= 325) {
	a0 = 7.121860E+02;
	a1 = -6.508397E-03;
	a2 = 3.553099E-03;
	a3 = -5.768018E-04;
	a4 = 1.836052E-05;
	a5 = -1.439496E-06;
	a6 = 9.355738E-07;
	a7 = -1.029321E-05;
	a8 = -4.275516E-07;
	a9 = 1.674928E-07;
	a10 = -1.042519E-08;
	a11 = -8.960122E-09;
	a12 = -6.927128E-10;
	a13 = 2.132294E-08;
	a14 = 2.239468E-10;
	a15 = 1.194714E-09;
	a16 = -9.815048E-11;
	a17 = 2.283092E-11;
	a18 = -3.102793E-10;
	a19 = 2.917436E-10;
	Bx = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
      // region of TPC2: |x|<=1.12m, |y|<=1.18m, 0.325m<z<=1.625m
      else if (325 < z && z <= 1625) {
	a0 = 7.122429E+02;
	a1 = -4.568123E-03;
	a2 = 3.435124E-03;
	a3 = -1.173135E-03;
	a4 = 1.728083E-05;
	a5 = -9.632399E-07;
	a6 = -4.271543E-06;
	a7 = -1.064613E-05;
	a8 = 4.601862E-07;
	a9 = 1.337533E-06;
	a10 =-9.865123E-09;
	a11 = -9.475482E-09;
	a12 = 5.577070E-10;
	a13 = 2.134384E-08;
	a14 = -2.861562E-10;
	a15 = 2.366843E-09;
	a16 = -1.584219E-11;
	a17 = 2.170653E-10;
	a18 = -1.499117E-10;
	a19 = -6.790860E-10;
	Bx = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
      // region of TPC3 |x|<=1.12m, |y|<=1.18m, 1.625m<z<= 2.90m
      else if (1625 < z && z <= 2900) {
	a0 = 6.737949E+02;
	a1 = 2.353285E-02;
	a2 = -7.567251E-03;
	a3 = 5.009049E-02;
	a4 = 7.011087E-05;
	a5 = 2.649338E-05;
	a6 = -6.446029E-05;
	a7 = -1.009779E-05;
	a8 = 3.432304E-06;
	a9 = -1.692050E-05;
	a10 = -2.262848E-08;
	a11 = -6.288378E-09;
	a12 = -1.685814E-08;
	a13 = 2.042897E-08;
	a14 = -1.625052E-08;
	a15 = 2.373631E-08;
	a16 = -1.597465E-10;
	a17 = 1.622080E-10;
	a18 = 1.482099E-09;
	a19 = 5.262672E-10;
	Bx = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
    }//end if (y == ...)
    else if (y < 0 &&  y >= -1180) {
      Bx = GetMeasInnerXValue(x,fabs(y),z) / unit::gauss;//RECURSIVE!!!BUT IT WORKS 
    }
  }// end if (x == ...)
  
  return Bx * unit::gauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetMeasInnerYValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for By within the measured area of NOMAD: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double By = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::mm;
  double y = y_in / unit::mm;
  double z = z_in / unit::mm;

  // The parametrised area is divided into 4 regions (TPC1, TPC2, TPC3 and the rest)
  // For each region a 3rd order polynomial description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;

  if (-1120 <= x && x <= 1120) { 
    if (0 <= y && y <= 1180) {//FOR THE MOMENT ASSUME SYMMETRY IN Y UNTIL FIT METHOD HAS BEEN IMPROVED (so that arm 3 can also be included)
      //region of TPC1: |x|<=1.12m, |y|<=1.18m, -0.975m<z<=0.325m
      if (-975 < z && z <= 325) {
	a0 = 3.317998E-02; 
	a1 = 1.782106E-03;
	a2 = -3.794425E-04;
	a3 = -5.113176E-05;
	a4 = -7.117799E-06;
	a5 = 9.601580E-06; 
	a6 = -1.999435E-07;
	a7 = 3.135181E-07; 
	a8 = 4.068269E-07; 
	a9 = -1.721770E-07;
	a10 = 5.623354E-09; 
	a11 = -9.587882E-09;
	a12 = 2.688968E-10; 
	a13 = -2.532657E-09;
	a14 = -3.505973E-10;
	a15 = -1.890408E-11;
	a16 = 2.241078E-09; 
	a17 = 5.678665E-10;
	a18 = 5.662437E-10; 
	a19 = -1.862113E-12;
	By = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
      // region of TPC2: |x|<=1.12m, |y|<=1.18m, 0.325m<z<=1.625m
      else if (325 < z && z <= 1625) {
	a0 = 1.007469E-01; 
	a1 = 1.961923E-04; 
	a2 = -6.501191E-04;
	a3 = 7.086530E-04; 
	a4 = -4.821129E-06;
	a5 = 9.948425E-06; 
	a6 = 1.077218E-06; 
	a7 = 9.603443E-07; 
	a8 = -4.608350E-07;
	a9 = -1.029192E-06;
	a10 = 4.698984E-09; 
	a11 = -8.982391E-09;
	a12 = -1.083649E-09;
	a13 = -3.688314E-09;
	a14 = 1.450533E-10; 
	a15 = 5.585552E-11; 
	a16 = 2.514698E-09; 
	a17 = -4.684959E-10;
	a18 = 2.235947E-10; 
	a19 = 3.112712E-10; 
	By = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
      // region of TPC3 |x|<=1.12m, |y|<=1.18m, 1.625m<z<= 2.90m
      else if (1625 < z && z <= 2900) {
	a0 = 1.289923E+01;  
	a1 = 2.182300E-03; 
	a2 = 4.077481E-03; 
	a3 = -2.220091E-02;
	a4 = -2.304038E-07;
	a5 = 2.286387E-05; 
	a6 = -2.283726E-06;
	a7 = -3.585013E-06;
	a8 = -5.296133E-06;
	a9 = 1.216952E-05; 
	a10 =4.434555E-11; 
	a11 =-1.645308E-08;
	a12 =7.991628E-10; 
	a13 =-2.877473E-09;
	a14 = -1.809401E-09;
	a15 = 2.807819E-10; 
	a16 = 1.218047E-09; 
	a17 = 3.363846E-09; 
	a18 = 2.954316E-10; 
	a19 = -2.105540E-09;
	By = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
    }//end if (y == ...)
    else if (y >= -1180) {
      By = -GetMeasInnerYValue(x,fabs(y),z) / unit::gauss;//RECURSIVE!!!BUT IT WORKS
    }
  }// end if (x == ...)
  
  return By * unit::gauss;  // get back to internal unit system
}


//********************************************************
double COMET::IOAMagneticField::GetMeasInnerZValue(double x_in, double y_in, double z_in) {
//********************************************************
  // Gives parametrisation for Bz within the measured area of NOMAD: |x|<1.66m, |y|<1.50m, |z|<3.3m
  double Bz = 0;
  // Within this function I am forced to leave the internal unit system 
  // since the simulation is paremeterized assuming centimeters and gauss
  double x = x_in / unit::mm;
  double y = y_in / unit::mm;
  double z = z_in / unit::mm;

  // The parametrised area is divided into 4 regions (TPC1, TPC2, TPC3 and the rest)
  // For each region a 3rd order polynomial description will be given with the parameters a0...an
  double a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19;

  if (-1120 <= x && x <= 1120) { 
    if (0 <= y && y <= 1180) {//FOR THE MOMENT ASSUME SYMMETRY IN Y UNTIL FIT METHOD HAS BEEN IMPROVED (so that arm 3 can also be included)
      //region of TPC1: |x|<=1.12m, |y|<=1.18m, -0.975m<z<=0.325m
      if (-975 < z && z <= 325) {
	a0 = -6.966551E-01;
	a1 = -6.502174E-03;
	a2 = -8.932913E-04;
	a3 = 6.264287E-04; 
	a4 = 1.193333E-05; 
	a5 = 7.229467E-07;  
	a6 = 8.026402E-07; 
	a7 = 9.576926E-07; 
	a8 = -1.026951E-06;
	a9 = 3.469261E-07; 
	a10 = -6.036586E-09;
	a11 = -5.485606E-10;
	a12 = -2.578502E-09;
	a13 = 1.335251E-10; 
	a14 = 1.241602E-09; 
	a15 = -1.521325E-10;
	a16 = -4.122226E-10;
	a17 = 2.787803E-10; 
	a18 = -5.449095E-10;
	a19 = 4.237943E-11;
	Bz = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
      // region of TPC2: |x|<=1.12m, |y|<=1.18m, 0.325m<z<=1.625m
      else if (325 < z && z <= 1625) {
	a0 = -1.416051E+00;
	a1 = -8.860277E-03;
	a2 = -7.235186E-08;
	a3 = 3.005751E-03; 
	a4 = 1.951060E-05;
	a5 = -5.179616E-07;
	a6 = -1.343574E-06;
	a7 = 4.249685E-07; 
	a8 = -7.519501E-07;
	a9 = -2.527755E-06;
	a10 = -9.118536E-09;
	a11 = 1.225402E-10; 
	a12 = -2.672049E-09;
	a13 = -4.247607E-11;
	a14 = 1.655568E-10; 
	a15 = 1.545208E-09;
	a16 = -3.454769E-10;
	a17 = 1.793908E-11; 
	a18 = 3.132915E-10; 
	a19 = 7.701537E-10; 
	Bz = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);      
      }
      // region of TPC3 |x|<=1.12m, |y|<=1.18m, 1.625m<z<= 2.90m
      else if (1625 < z && z <= 2900) {
	a0 = 2.299212E+01;
	a1 = -3.735329E-02;
	a2 = -7.819633E-03;
	a3 = -2.084685E-02;
	a4 = 1.703276E-05; 
	a5 = -8.131301E-06;
	a6 = 3.099678E-05; 
	a7 = -5.151532E-06;
	a8 = 1.489178E-05; 
	a9 = 2.293184E-06; 
	a10 = 9.400963E-10; 
	a11 = 3.856641E-09; 
	a12 = -1.055542E-08;
	a13 = 4.942807E-10; 
	a14 = 1.498134E-09; 
	a15 = -5.136526E-09;
	a16 = 9.942945E-11; 
	a17 = 2.601885E-09;
	a18 = -5.441056E-09;
	a19 = 1.157495E-09;
	Bz = (a0 + a1*x + a2*y + a3*z + a4*x*x + a5*x*y + a6*x*z + a7*y*y + a8*y*z + a9*z*z + a10*x*x*x + a11*x*x*y + a12*x*x*z + a13*x*y*y + a14*x*y*z + a15*x*z*z + a16*y*y*y + a17*y*y*z + a18*y*z*z + a19*z*z*z);
      }
    }//end if (y == ...)
    else if (y >= -1180) {
      Bz = GetMeasInnerZValue(x,fabs(y),z) / unit::gauss;//RECURSIVE!!!BUT IT WORKS
    }
  }//end if (x == ...)
   
  return Bz * unit::gauss;  // get back to internal unit system
}

//********************************************************
double COMET::IOAMagneticField::GetMeasOuterXValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}

//********************************************************
double COMET::IOAMagneticField::GetMeasOuterYValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}

//********************************************************
double COMET::IOAMagneticField::GetMeasOuterZValue(double x, double y, double z) {
//********************************************************
  // Temporary definition of the magnetic field outside measured and parametrized volume
  return 0;
}


//********************************************************
TVector3 COMET::IOAMagneticField::GetP2PMeasFieldValue(const TVector3& point) {
//********************************************************
  //std::cout << " GetP2PValue " << std::endl;
  double X = point.X();
  double Y = point.Y();
  double Z = point.Z();

  if (fNormalizationFromCurrent){
    // scale factor based on the coil current
    UpdateScaleFactor();
  }
 

  //std::cout << " Scale factor " << fscale_factor << std::endl;
  //std::cout << " field stength " << fFieldStrengthX << std::endl;
  //std::cout << " P2P orgin X val " << GetP2PMeasFieldValue(TVector3(0,0,0)).X() << std::endl;

  TVector3 p000, p001, p010, p011, p100, p101, p110, p111;
  int arm_nr = 0;
  
  // define different regions and how to get the value for each
  // some border conditions are redundant in order to see directly what are the exact dimesions of each region
  //0. volume outside the mapped region
  if (fabs(X)>1083 || (Y<-1059.7 || Y>1175.3) || (Z<-887.9 || Z>2962.1)) {
    return COMET::IOAMagneticField::GetNOMADFieldValue(point);
  }
  //1. region which is covered by arm2
  else if (fabs(X)<=1083 && (Y>=-804.7 && Y<=1175.3) && (Z>=-887.9 && Z<=2962.1)) {
    arm_nr = 2;
  }
  //2. region covered by arm3 exlusively
  else if (fabs(X)<=855 && (Y>=-1059.7 && Y<=-861.7) && (Z>=-686.5 && Z<=2962.1)) {
    arm_nr = 3;
  }
  //3. extension of arm3 in z-region
  else if (fabs(X)<=855 && (Y>=-1059.7 && Y<=-861.7) && (Z>=-887.9 && Z<-686.5)) {
    TVector3 ptemp (X, Y, -686.4);
    return COMET::IOAMagneticField::GetP2PMeasFieldValue(ptemp);
  }
  //4. regions below arm2 and out of reach of arm3 in x-direction
  else if (fabs(X)>855 && (Y>=-1059.7 && Y<=-861.7) && (Z>=-887.9 && Z<=2962.1)) {
    TVector3 ptemp (855, Y, Z);
    return COMET::IOAMagneticField::GetP2PMeasFieldValue(ptemp);
  }
  //5. region right below arm2, make a transition of the values from arm2 to the values derived from arm3
  else if (fabs(X)<=1083 && (Y>-861.7 && Y<-804.7) && (Z>=-887.9 && Z<=2962.1)) {
    TVector3 ptemp_arm2 (X, -804.7, Z);
    TVector3 ptemp_arm3 (X, -861.7, Z);
    float weight_arm2 = (Y + 861.7)/(861.7-804.7);//between 0 and 1
    return weight_arm2*(COMET::IOAMagneticField::GetP2PMeasFieldValue(ptemp_arm2)) + (1 - weight_arm2)*(COMET::IOAMagneticField::GetP2PMeasFieldValue(ptemp_arm3));
  }
  //ERROR: everything should be covered by now  -> making sure you are not exceeding the measurement limits
  else {
    std::cerr << "ERROR: Point (" << X << ", " << Y << ", " << Z << ") not associated with a field value.Take NOMAD value and continue." << std::endl;
    return COMET::IOAMagneticField::GetNOMADFieldValue(point);
  }

  //make linear interpolation in case of arm_nr=2 or arm_nr=3, i.e. we are within a well defined mapped region
  /*
  std::map<TVector3, TVector3, classcomp> *ptrMapB = 0;
  std::set<float> *ptrNodesX = 0;
  std::set<float> *ptrNodesY = 0;
  std::set<float> *ptrNodesZ = 0;
  std::set<float>::iterator it;
  

  
  TVector3 magfieldarm2[39][41][78];
  COMET::IFieldManager::GetBMapArr(2,magfieldarm2);
  ptrMapB = COMET::IFieldManager::Get().GetBMap(arm_nr);
  ptrNodesX = COMET::IFieldManager::Get().GetNodesX(arm_nr);
  ptrNodesY = COMET::IFieldManager::Get().GetNodesY(arm_nr);
  ptrNodesZ = COMET::IFieldManager::Get().GetNodesZ(arm_nr);
  
  //Make sure you are not exceeding the measurement limits
  if (X < *ptrNodesX->begin() || Y < *ptrNodesY->begin() || Z < *ptrNodesZ->begin()
      || X > *(--ptrNodesX->end()) || Y > *(--ptrNodesY->end()) || Z > *(--ptrNodesZ->end())) {
    std::cerr << "ERROR: Point (" << X << ", " << Y << ", " << Z << ") exceeded limit of measurement region. Take NOMAD value and continue." << std::endl;
    return COMET::IOAMagneticField::GetNOMADFieldValue(point);
  } 
  */
  
  int xidx,yidx,zidx;

  int xmax,ymax,zmax;
  float dx,dy,dz;
 
  float xstart,ystart,zstart ;
  
  TVector3 *magfieldarm;
  ymax = 41;
  zmax = 78;

  if(arm_nr == 2){
    //magfieldarm[39][41][78];
    //fieldman->GetBMapArr(arm_nr,magfieldarm);
    magfieldarm = COMET::IFieldManager::Get().GetBMapArr(arm_nr);
    dx = 57.0;
    dy = 49.5;
    dz = 50.0;
    xstart = -1083.0;
    ystart = -804.70001220703125;
    zstart = -887.9000244140625;
    xmax = 39;

    }
  
  else if(arm_nr == 3){
    magfieldarm = COMET::IFieldManager::Get().GetBMapArr(arm_nr); 
    //COMET::IFieldManager::GetBMapArr(2,magfieldarm); //BADBADBD
    //COMET::IFieldManager::GetBMap2Arr(arm_nr,magfieldarm);
    dx = 171.0;
    dy = 49.5;
    dz = 50.0;

    xstart = -855.0;
    ystart = -1059.699951171875;
    zstart = -686.5;
    xmax =11;
  }
  
  
  xidx = (int) floor( (X-xstart)/dx );
  yidx = (int) floor( (Y-ystart)/dy ); 
  zidx = (int) floor( (Z-zstart)/dz ); 

  /*
  xlow = xstart+xidx*dx;
  ylow = ystart+yidx*dy;
  zlow = zstart+zidx*dz; 
  
  xhig = xstart+(xidx+1)*dx;
  yhig = ystart+(yidx+1)*dy;
  zhig = zstart+(zidx+1)*dz; 
  */

  //magfieldarm[zidx*41+yidx*39+xidx];
  //std::cout << magfieldarm[zidx*41+yidx*39+xidx] << std::endl;
  float Bx_final=0;
  float By_final=0;
  float Bz_final=0;
  
  for(int xval = xidx; (xval<xidx+2) && (xval<xmax); xval++){
    for(int yval = yidx; (yval<yidx+2) && (yval<ymax); yval++){
      for(int zval = zidx; (zval<zidx+2) &&(zval<zmax); zval++){
	
	
	if(  magfieldarm[zval*1599+yval*39+xval].X() == 0){
	  std::cout << " x " << xval << " y " << yval << " z " << zval << std::endl;
	}
	
	Bx_final += magfieldarm[zval*1599+yval*39+xval].X()* (dx-fabs(X-(xstart+xval*dx)))* (dy-fabs(Y-(ystart+yval*dy)))*(dz-fabs(Z-(zstart+zval*dz))) / (dx*dy*dz);
	By_final += magfieldarm[zval*1599+yval*39+xval].Y()* (dx-fabs(X-(xstart+xval*dx)))* (dy-fabs(Y-(ystart+yval*dy)))*(dz-fabs(Z-(zstart+zval*dz))) / (dx*dy*dz);
	Bz_final += magfieldarm[zval*1599+yval*39+xval].Z()* (dx-fabs(X-(xstart+xval*dx)))* (dy-fabs(Y-(ystart+yval*dy)))*(dz-fabs(Z-(zstart+zval*dz))) / (dx*dy*dz);
      }
    }
  }

 
  
  TVector3 B_Vec;
  B_Vec.SetX(fScaleFactor*Bx_final * unit::gauss);
  B_Vec.SetY(fScaleFactor*By_final * unit::gauss);
  B_Vec.SetZ(fScaleFactor*Bz_final * unit::gauss);

  //std::cout << " fscale_factor*Bx_final * unit::gauss " << fscale_factor*Bx_final * unit::gauss << std::endl;
 
  //std::cout << Bx_final * unit::gauss << std::endl;
 return B_Vec;
}