#ifndef TOAMagneticField_h
#define TOAMagneticField_h

#include "TVector3.h"
#include "TLorentzVector.h"
#include "IMagnetCurrent.hxx"
#include <map>
#include <set>


namespace COMET {
  class IOAMagneticField;
};

///    \brief comet magnetic field (current nominal field is from SIMFIELD)
///    \author Christian Hansen, hansen@uvic.ca
///
///    IOAMagneticField
///   
///    This assumes the magnet is centered at the origin.
///    IOAMagneticField provides maps of the magnetic field in the 
///    COMET detector. The type of map to use is controlled by either 
///    a parameter in the data base or as an input to a set function. 
///    The current nominal field map (type 1) is the field map measured
///    in September 2009. It is the put as type 9 for later reference. 
///    An other unhomogenic field map (chosen as type 4) was 
///    given to us by NOMAD, the experiment that used the UA1 magnet 
///    before T2K. It is important to note that NOMAD used this magnet with
///    0.4 tesla as the field strength and not 0.2 tesla as COMET.
///    This can however be changed by the SetFieldStrength function.
///    We do not know if this NOMAD field map orgins from measurements 
///    or simulations. The f77 code has been translated to c++ and can 
///    now be used to get the B field in any point in the COMET experiment.
///    There are several other possible B field type options:
///    Without any B-Field (type 0).
///    The nominal field with an added randomly error (type 2)
///    where the strength of the error can be set by the SetFieldError function.
///    A "perfect" B field, i.e. only an x component (type 3)
///    The field map currently used by the SimG4 (type 5).
///    The results of a simulation with OPERA/TOSC (type 6)
///    An obsolete map using polynomial fits to the measured data (type 7)
///    A homogenous map where you can set all three directions of the B vector (type 8)
///    A map giving back the position in space. This intended to be used for
///    sanity checks exclusively, hence the rather uncommon number (type 111)


class COMET::IOAMagneticField {
 public:

  /// Constructor, default: using the B-Field type requested in DB
  IOAMagneticField();
  /// Destructor
  ~IOAMagneticField();

  /*!
   * Overwriting the B-Field type requested in DB
   *  - 0; A zero B-field
   *  - 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).  
   *  - 2; Nominal field with an error (error strength can be set by SetFieldError)
   *  - 3; "Perfect" B Field. Only Bx = field strength
   *  - 4; NOMAD field map. Note that the NOMAD field
   *       map has by default 0.4T as the field strength. This can be changed by SetFieldStrength.
   *  - 5; The B Field currently used by SimG4 (will change)
   *  - 6; The SIMFIELD simulated field map 
   *  - 7; Obsolete map using polyomial fits to the mapped data.
   *  - 8; Simple homogenous field where one can set each of the B-field coordinates.
   *  - 9; A point to point field map of the measurements done in Sep 2009. This map uses a linear interpoaltion between the meaurement points.
   *  - 111; This is not a real map, it is just there for the purpose of testing the geometrical implementation of the map.
  */
  void SetFieldType(int fieldType);
  /// Returns current field type
  int GetFieldType() {return fBFieldType;};

  /// Set the maximal field strength. Default is 0.4T for NOMAD and 0.2T for the other map options.
  void SetFieldStrength(double bMax);
  /// Returns current field strength (depending on what the current field type is)
  double GetFieldStrength();

  /// Set the field error strength (only used for the B-Field type 2)
  void SetFieldError(double err);
  /// Returns current field error
  double GetFieldError() {return fBErr;};

  // method to set a homogeneous field value
  void SetFieldValue(double x, double y, double z);

  /// Method to return the B field at a position (using the current B field type)
  TVector3 GetFieldValue(const TLorentzVector& point);
  /// Method to return the B field at a position (using the current B field type)
  TVector3 GetFieldValue(const TVector3& point);
  /// Method to return the B field at a position for a specific B field type
  TVector3 GetFieldValue(const TLorentzVector& point, int fieldType);
  /// Method to return the B field at a position for a specific B field type
  TVector3 GetFieldValue(const TVector3& point, int fieldType);
  /// Method to be called by SimG4
  void GetFieldValue(const double point[3], double *BField);

  /// Methods to return the B field at a position from a special field map 
  /// These do not care about the current field type
  TVector3 GetNominalFieldValue(const TVector3& point);
  TVector3 GetNominalWithErrorFieldValue(const TVector3& point) {return GetFieldValue(point, 2);};
  TVector3 GetPerfectFieldValue(const TVector3& point) {return GetFieldValue(point, 3);};
  TVector3 GetNOMADFieldValue(const TVector3& point);
  TVector3 GetSIMG4FieldValue(const TVector3& point) {return GetFieldValue(point, 5);};
  TVector3 GetSimFieldValue(const TVector3& point);
  TVector3 GetMeasFieldValue(const TVector3& point); 
  TVector3 GetHomogeneousFieldValue(const TVector3& point) {return GetFieldValue(point, 8);}  
  //Point to point (P2P) MEASURED MAP August 2010
  TVector3 GetP2PMeasFieldValue(const TVector3& point);

 private: 
  int fBFieldType; 
  double fBErr; //Only used for field type 3
  double fFieldStrengthX; //Default is 0.4T for NOMAD and 0.2T for the other field types 

  /// only use for homogeneous field definition (type 8)
  double fFieldStrengthY;
  double fFieldStrengthZ;
  
  /// Flag to decide whether to normalize the measured maps using magnet current.
  bool fNormalizationFromCurrent;
  /// The last time that normalization was updated.
  int fLastNormalizationUpdate;
  /// The scale factor to use for normalizing measured maps.
  double fScaleFactor;

  // Method for updating the scale factor based on the magnet current.
  // This method checks that the time is different from the last call;
  // if not, then uses cached result.  Also, only changes scale factor for
  // data; MC is left unchanged.
  void UpdateScaleFactor();

  
  /// Class for holding the retrieving the current magnet current.
  IMagnetCurrent* fMagnetCurrent;

  
  // current at which field measurements were done
  double fReferenceCurrent;

	//parameters of fitting function describing the B dependce on I
  double fCurrentToFieldC0;
  double fCurrentToFieldC1;
  double fCurrentToFieldC2;

  int fDefaultFieldType;


  void Initialize();
  
  bool IsMC();

  /// Measured field with polynomial fits (as of March 2010)
  double GetMeasFieldXValue(const TVector3& point);
  double GetMeasFieldYValue(const TVector3& point);
  double GetMeasFieldZValue(const TVector3& point);
  double GetMeasInnerXValue(double x_in, double y_in, double z_in);
  double GetMeasInnerYValue(double x_in, double y_in, double z_in);
  double GetMeasInnerZValue(double x_in, double y_in, double z_in);
  double GetMeasOuterXValue(double x, double y, double z);
  double GetMeasOuterYValue(double x, double y, double z);
  double GetMeasOuterZValue(double x, double y, double z);

  /// SIMULATED FIELD with polynomial fits (2008 by Eike Frank)
  double GetSimFieldXValue(const TVector3& point);
  double GetSimFieldYValue(const TVector3& point);
  double GetSimFieldZValue(const TVector3& point);
  double GetSimInnerXValue(double x_in, double y_in, double z_in);
  double GetSimInnerYValue(double x_in, double y_in, double z_in);
  double GetSimInnerZValue(double x_in, double y_in, double z_in);
  double GetSimOuterXValue(double x, double y, double z);
  double GetSimOuterYValue(double x, double y, double z);
  double GetSimOuterZValue(double x, double y, double z);  

  /// NOMAD
  /// The field map x, y and z components
  double GetNOMADFieldXValue(const TVector3& point);
  double GetNOMADFieldYValue(const TVector3& point);
  double GetNOMADFieldZValue(const TVector3& point);
  /// The field x, y and z components in the inner, measured region
  double GetNOMADInnerXValue(double x, double y, double z);
  double GetNOMADInnerYValue(double x, double y, double z);
  double GetNOMADInnerZValue(double x, double y, double z);
  /// The field x, y and z components in the oute, unmeasured region
  double GetNOMADOuterXValue(double x, double y, double z);
  double GetNOMADOuterYValue(double x, double y, double z);
  double GetNOMADOuterZValue(double x, double y, double z);
  /// Soft transfer between 0 to 1 for x within [0,1]
  double Soft01(double x);
  /// Member Variablestpcfitstatus_fgd1_ccqe1binPNC.root
  double fYMinC;
  double fYMaxC;
  double fZMinC;
  double fZMaxC;
  double fRMinC;
  double fRMaxC;
  double fXMinM;
  double fXMaxM;
  double fYMinM;
  double fYMaxM;
  double fZMaxM;
  double fZMagnet;
  double fXCenl;
  double fYCenl;
  double fZCenl;
  double fDZCees;
  double fHZGaps;
  double fHSGap3;
  double fBFMin;
  double fBTMax;
  double fBLMax;
  //additional from EIKE
  double fXMaxTPC;
  double fYMaxTPC;
  double fZMinTPC;
  double fZMaxTPC;

  //FROM SIMG4
  TVector3 GetMCFieldValue(const TVector3& pos);
  double fYokeOuterWidth; 
  double fYokeOuterHeight; 
  double fYokeOuterLength;
  double fYokeInnerWidth; 
  double fYokeInnerHeight; 
  double fYokeInnerLength;
  double fCoilOuterWidth; 
  double fCoilOuterHeight; 
  double fCoilOuterLength;
  double fCoilInnerWidth; 
  double fCoilInnerHeight; 
  double fCoilInnerLength;

};
  
#endif // #ifdef TOAMagneticField_hxx