#ifndef Geom_C
#define Geom_C

// Provides rootish Int_t and Double_t type types
#if !defined (__CINT__) || defined (__MAKECINT__)
#include "Rtypes.h"
#endif

#include "TMath.h"
#include "TString.h"

// std::pair is provided in utility.
#include <utility>
#include <iostream>
#include <cassert>

namespace Geom
{
  // This array maps an integer index 0-35 corresponding to DAQ channels in the V1742
  // and indexes in digiChans to other integers which correspond to physical things in
  // the prototype.  Most notably if chMap[i] = k for k in [0,27], then that i is mapped
  // to a drift chamber sense wire.  The inverse mapping (below) is much more useful.
  const Int_t chMap[36] = { 0, 1, 2, 3, 4, 5, 6, 7,32, 8, 9,10,11,12,13,14,15,33,
                           16,17,18,19,20,21,22,23,34,24,25,26,27,28,29,30,31,35};

  // This array maps an integer index 0-28 corresponding to actual sense
  // wires to an index 0-35 corresponding to channels in the DAQ system 
  // (and in the digiChan array).  So sense wire #10 is in "DAQ channel" chInvMap[10].
  const Int_t chInvMap[28] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 
                              16, 18, 19, 20, 21, 22, 23, 24, 25, 27, 28, 29, 30};

  //vector<Int_t> chMap_v = chMap; // Only works in C++0x or 11, apparently.

  // This is a map of the actual sense wire coordinates in cm.
  const Double_t chPos[28][2] = {{-2.1,4.9},{-0.7,4.9},{0.7,4.9},{2.1,4.9},
                                 {-1.4,3.5},{0.0,3.5},{1.4,3.5},
                                 {-2.1,2.1},{-0.7,2.1},{0.7,2.1},{2.1,2.1},
                                 {-1.4,0.7},{0.0,0.7},{1.4,0.7},
                                 {-2.1,-0.7},{-0.7,-0.7},{0.7,-0.7},{2.1,-0.7},
                                 {-1.4,-2.1},{0.0,-2.1},{1.4,-2.1},
                                 {-2.1,-3.5},{-0.7,-3.5},{0.7,-3.5},{2.1,-3.5},
                                 {-1.4,-4.9},{0.0,-4.9},{1.4,-4.9}};

  // This map is used to fill the cell occupancy TH2.  It is basically
  // a set of coordinates by sublayer.
  extern const Int_t chCoords[28][2] = {{0,0},{1,0},{2,0},{3,0},
                                        {0,1},{1,1},{2,1},
                                        {0,2},{1,2},{2,2},{3,2},
                                        {0,3},{1,3},{2,3},
                                        {0,4},{1,4},{2,4},{3,4},
                                        {0,5},{1,5},{2,5},
                                        {0,6},{1,6},{2,6},{3,6},
                                        {0,7},{1,7},{2,7}};

  // Obtained by fitting the derivative of the histogram of all the t0
  // measurements in runs:
  // 101,102,104,105,110,112,113,115,117,119,121,155,170,177,179,411
  const Double_t shift = +0.05e+02;                                  
  extern const Double_t t0_Offset = 1.59229e+02 + shift;
  
  const Double_t FieldPos[][2] = { // Field wires first
                                  {-0.7,-0.0},{-0.7,-1.4},{-0.7,-2.1},{-0.7,-2.8},
                                  {-0.7,-4.2},{-0.7,-4.9},{-0.7,-5.6},{-0.7,0.7},
                                  {-0.7,1.4},{-0.7,2.8},{-0.7,3.5},{-0.7,4.2},
                                  {-0.7,5.6},{-0.7,6.3}, 
                                  {-1.4,-0.0},{-1.4,-0.7},{-1.4,-1.4},{-1.4,-2.8},
                                  {-1.4,-3.5},{-1.4,-4.2},{-1.4,-5.6},{-1.4,-6.3},
                                  {-1.4,1.4},{-1.4,2.1},{-1.4,2.8},{-1.4,4.2},
                                  {-1.4,4.9},{-1.4,5.6},
                                  {-2.1,-0.0},{-2.1,-1.4},{-2.1,-2.1},{-2.1,-2.8},
                                  {-2.1,-4.2},{-2.1,-4.9},{-2.1,-5.6},{-2.1,0.7},
                                  {-2.1,1.4},{-2.1,2.8},{-2.1,3.5},{-2.1,4.2},
                                  {-2.1,5.6},{-2.1,6.3},
                                  {-2.8,-0.0},{-2.8,-0.7},{-2.8,-1.4},{-2.8,-2.8},
                                  {-2.8,-3.5},{-2.8,-4.2},{-2.8,-5.6},{-2.8,-6.3},
                                  {-2.8,1.4},{-2.8,2.1},{-2.8,2.8},{-2.8,4.2},
                                  {-2.8,4.9},{-2.8,5.6},
                                  {-3.5,-0.0},{-3.5,-1.4},{-3.5,-2.1},{-3.5,-2.8},
                                  {-3.5,-4.2},{-3.5,-4.9},{-3.5,-5.6}, {-3.5,0.7},
                                  {-3.5,1.4},{-3.5,2.8},{-3.5,3.5},{-3.5,4.2},{-3.5,5.6},
                                  {0.0,-0.0},{0.0,-0.7},{0.0,-1.4},{0.0,-2.8},
                                  {0.0,-3.5},{0.0,-4.2},{0.0,-5.6},{0.0,-6.3},
                                  {0.0,1.4},{0.0,2.1},{0.0,2.8},{0.0,4.2},{0.0,4.9},{0.0,5.6},
                                  {0.7,-0.0},{0.7,-1.4},{0.7,-2.1},{0.7,-2.8},
                                  {0.7,-4.2},{0.7,-4.9},{0.7,-5.6},
                                  {0.7,0.7},{0.7,1.4},{0.7,2.8},{0.7,3.5},{0.7,4.2},
                                  {0.7,5.6},{0.7,6.3},
                                  {1.4,-0.0},{1.4,-0.7},{1.4,-1.4},{1.4,-2.8},
                                  {1.4,-3.5},{1.4,-4.2},{1.4,-5.6},{1.4,-6.3},
                                  {1.4,1.4},{1.4,2.1},{1.4,2.8},{1.4,4.2},{1.4,4.9},{1.4,5.6},
                                  {2.1,-0.0},{2.1,-1.4},{2.1,-2.1},{2.1,-2.8},
                                  {2.1,-4.2},{2.1,-4.9},{2.1,-5.6},
                                  {2.1,0.7},{2.1,1.4},{2.1,2.8},{2.1,3.5},
                                  {2.1,4.2},{2.1,5.6},{2.1,6.3},
                                  {2.8,-0.0},{2.8,-0.7},{2.8,-1.4},{2.8,-2.8},
                                  {2.8,-3.5},{2.8,-4.2},{2.8,-5.6},{2.8,-6.3},
                                  {2.8,1.4},{2.8,2.1},{2.8,2.8},{2.8,4.2},{2.8,4.9},{2.8,5.6},
                                  {3.5,-0.0},{3.5,-1.4},{3.5,-2.1},{3.5,-2.8},
                                  {3.5,-4.2},{3.5,-4.9},{3.5,-5.6},
                                  {3.5,0.7},{3.5,1.4},{3.5,2.8},{3.5,3.5},{3.5,4.2},{3.5,5.6},
                                  // Then guard wires
                                  {-0.7,-6.3},{-1.4,6.3},{-2.1,-6.3},{-2.8,-2.1},
                                  {-2.8,-4.9},{-2.8,0.7},{-2.8,3.5},{-2.8,6.3},
                                  {-3.5,-0.7},{-3.5,-3.5},{-3.5,2.1},{-3.5,4.9},
                                  {0.0,6.3},{0.7,-6.3},{1.4,6.3},{2.1,-6.3},
                                  {2.8,-2.1},{2.8,-4.9},{2.8,0.7},{2.8,3.5},
                                  {2.8,6.3},{3.5,-0.7},{3.5,-3.5},{3.5,2.1},{3.5,4.9}
                                };

  std::pair<Double_t,Double_t> field_coords(Int_t N)
  {
    if(N < 0 or N > 176)
    {
      std::cout << "Error, field and guard wire index must be within 0-176, returning NaNs." << std::endl;
      return std::make_pair(TMath::QuietNaN(), TMath::QuietNaN());
    }
    return std::make_pair(FieldPos[N][0],FieldPos[N][1]);
  }

  std::pair<Double_t,Double_t> wire_coords(Int_t N)
  { // This returns the wire coordinates in cm
    // for channel N.
    
    // The 28 channels are arranged as such:
    //  00  01  02  03  \ Two sublayers 
    //    04  05  06    / form a superlayer
    //  07  08  09  10  - Upper sublayer
    //    11  12  13    - lower sublayer
    //  14  15  16  17
    //    18  19  20
    //  21  22  23  24
    //    25  26  27

    // x = 0 is in the middle of the
    // chamber (intersecting the centres of
    // cells 5, 12, 19, 26).
    // y = 0 is in the middle of the chamber,
    // at the boundary between cells 11, 12, 13 and
    // 14, 15, 16, 17.

    // Each cell is 1.4cm wide, and the
    // lower sublayer is offset from the upper
    // by 0.7cm.

    // Each superlayer is 2.8cm tall.
    
    const Int_t SL = N/7;    // Superlayer number
    const Int_t ISL = N % 7; // Index within superlayer
    const Int_t UL = ISL > 3;// 0 or 1 if in upper or lower sublayer
    const Int_t IL = ISL % 4;// Index within sublayer
    return std::make_pair(-2.1+UL*0.7+IL*1.4, 4.9-SL*2.8-UL*1.4);
  }

  void check_Geom()
  { // This assertion loop checks that the channel map
    // and inverse map actually match.  It also checks
    // that the generated wire coordinates match the hand-
    // typed ones.
    for(Int_t i=0;i<28;i++)
    {
      assert(chMap[chInvMap[i]] == i);
      assert(TMath::Abs(chPos[i][0] - wire_coords(i).first) < 0.0001 &&
             TMath::Abs(chPos[i][1] - wire_coords(i).second) < 0.0001);
    }
  }

  bool Crosses_Cell(int i, double x0, double theta)
  {
    const double L = chPos[i][0] - 0.7;
    const double R = chPos[i][0] + 0.7;
    const double B = chPos[i][1] - 0.7;
    const double T = chPos[i][1] + 0.7;
    const double m = TMath::Tan(theta); // This m is the slope of the track/
    auto y = [&](double x_arg){return m*(x_arg-x0);};
    auto x = [&](double y_arg){return x0+y_arg/m;};
    const double yL = y(L);
    const double yR = y(R);
    const double xB = x(B);
    const double xT = x(T);
    unsigned int crossings = 0;
    if(B <= yL and yL <= T) crossings++;
    if(B <= yR and yR <= T) crossings++;
    if(L <= xB and xB <= R) crossings++;
    if(L <= xT and xT <= R) crossings++;
    // if(crossings != 0 and crossings != 2) 
    //   {
    // 	std::cout << "In checking whether track at x0: "
    // 		  << TString::Format("%g, theta: %g crosses cell %d, %d",
    // 				     x0,theta,i,crossings)
    // 		  << " crossings were found." << std::endl;
    //   }
    assert(crossings < 4);
    return crossings > 0;
  }
  
  std::vector<Int_t> Cells_Crossed(double x0, double theta)
  {
    std::vector<Int_t> crossed;
    for(int i = 0;i<28;i++)
      {
        if(Crosses_Cell(i,x0,theta)) crossed.push_back(i);
      }
    return crossed; // used to be std::move(crossed);, but apparently this is better.
  }

  // double DistanceFromWire_old(int chan, const double * pars)
  // { // This function returns a signed quantity.  
  //   // The magnitude is the actual
  //   // distance from the sense wire #chan to the point of closest approach (POCA)
  //   // along the track.  
  //   // The sign is the sign of the difference between the sense wire x coordinate and 
  //   // the POCA x coordinate.  In order words: + for tracks of the right, - for
  //   // tracks on the left.  

  //   const double x0 = pars[0];
  //   const double theta = pars[1];
    
  //   const double sintheta = TMath::Sin(theta);
  //   // Here m is the slope of a line orthogonal to the track,
  //   // used for the shortest distance between the track and the sense wire.
  //   const double m = 1.0/TMath::Tan(theta); 
    
  //   // Sense wire coordinates.
  //   const double xi = Geom::chPos[chan][0];
  //   const double yi = Geom::chPos[chan][1];
    
  //   // Coordinates of POCA on proposed track.
  //   const double yp = sintheta*sintheta*(yi-m*(x0-xi));
  //   const double xp = m*yp+x0;
    
  //   // Distance of POCA on proposed track from sense wire.
  //   const double D = TMath::Sqrt((xp-xi)*(xp-xi) + (yp-yi)*(yp-yi));
  //   return TMath::Sign(D,xp-xi);
  // }

  inline double DistanceFromWire(int chan, const double * pars)
  {
    const double x0 = pars[0];
    const double theta = pars[1];
    const double xi = Geom::chPos[chan][0];
    const double yi = Geom::chPos[chan][1];
    
    return yi*TMath::Cos(theta) + (x0-xi)*TMath::Sin(theta);
  }

  enum Side {LEFT, RIGHT, BOTH};
  const TString SideStrings[] = {"LEFT","RIGHT","BOTH"};
  
} // Ends Geom namespace.
#endif