#ifndef _CAL_TOWER_HH_
#define _CAL_TOWER_HH_

// Modified by X. Rouby for integration in Delphes
// !!! the numbering of the towers (iEta, iPhi) is **NOT** following any convention of the corresponding experiment
// v3 -- not working

#include <cmath>
#include <iostream>

#ifndef M_PI
#define M_PI  3.141592653589793238462643383279502884197 
#endif

const double pi = acos(-1);

static bool einmal=true;
/*
 // CDF data : 22 towers. step=2.7° at the beginning, and after, it changes
const unsigned int ntower = 22;
const double TOWER_THETA[ntower+1] = {3.000,  5.700,  8.400, 11.100, 13.800, 16.500, 19.200, 21.900, 24.600, 27.300, 30.000,// step=2.7°
          33.524,  36.822, 40.261, 43.614, 47.436, 51.790, 56.735, 62.310, 68.516, 75.297, 82.526, 90.000 };
const double tower_eta_edges[ntower+1] = {
    0.0,   0.130817437415,   0.259479460739,   0.384075299436,   0.503273260393,   0.616250691646,   0.72264587494,   0.822472442947,   
    0.916007799535,   1.00361552815,   1.10000635195,   1.19999982651,   1.31695789692,   1.41531379676,   1.52309848418, 1.64260787639,  
    1.77704423284,  1.93106912741,  2.1118548488,  2.33129451223,  2.61134904017,  3.00008274261,  3.64253335732 };
const double tower_dphi[ntower] = { // list of the tower size in phi (in degrees) 
	7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 15.,
	15., 15., 15., 15., 15., 7.5, 7.5, 7.5,
	7.5, 7.5, 7.5, 7.5, 7.5, 7.5 };
*/

// CMS data
const unsigned int ntower = 40;
const double tower_eta_edges[ntower+1] = { // list of the edges of each tower, in eta, for eta>0, assuming a symmetric detector in eta<0 
    0., // lower limit of the most central tower
    0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
    1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
    4.350, 4.525, 
    4.700,  // lower limit of the most forward tower
    5.000}; // higher limit of the most forward tower

const double tower_dphi[ntower] = { // list of the tower size in phi (in degrees)
     5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 10,
    10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; 


class CalTower
{
 public:
  double Et,eta,phi;
  int iEta,iPhi;

  CalTower(): Et(0), eta(0), phi(0), iEta(-1), iPhi(-1) {}

  CalTower(const double Et0, const double eta0, const double phi0): Et(Et0), eta(eta0), phi(phi0)
  {

	if (fabs(eta) > tower_eta_edges[ntower]) {iEta = -1; iPhi = -1; } // outside the detector
        else {  
		const unsigned int nn = ntower*2; 
		double all_towers[nn];
		unsigned int i =0;
		for(unsigned int j= ntower-1; i<ntower; j--,i++) { all_towers[i] = -tower_eta_edges[j];}
		for(unsigned int j= 1; j <ntower+1; j++, i++) { all_towers[i] =  tower_eta_edges[j]; }
		for(i = 0; i< nn; i++) { if(eta < all_towers[i]) {iEta=i; break;} }

		for(unsigned int j=0; j<ntower; j++) {
			int wedge = int(360./tower_dphi[j]);
			iPhi = 0; int(phi/(2.*pi)* wedge)%wedge;
		}
	}
  }

  CalTower(const double Et0, const double eta0, const double phi0, const int iEta0, const int iPhi0): Et(Et0), eta(eta0), phi(phi0), iEta(iEta0), iPhi(iPhi0) {}

  CalTower(const CalTower& c): Et(c.Et), eta(c.eta), phi(c.phi), iEta(c.iEta), iPhi(c.iPhi) {}

  bool isEqual(const CalTower& c) {
    return Et == c.Et && eta == c.eta && phi == c.phi && iEta == c.iEta && iPhi == c.iPhi;
  }

// converts theta angle [degrees] into pseudorapidity
 double pseudorapidity(const double theta) {
         return -log(tan( (theta/2.) * (pi/180.) ));
 }

};

/***** 
 e.g.  CDF
  iEta==4  si -3.642 < eta < -3.000
  iEta==47 si  3.000 < eta <  3.642
  iEta==-1 si eta < -3.642 ou si eta > 3.642

 i   theta    pseudorapidity     iEta
--------------------------------------
 0   3.0      -3.00008274261        4
 1   5.7      -2.61134904017        5
 2   8.4      -2.33129451223        6
 3   11.1     -2.1118548488         7
 4   13.8     -1.93106912741        8
 5   16.5     -1.77704423284        9
 6   19.2     -1.64260787639       10
 7   21.9     -1.52309848418       11
 8   24.6     -1.41531379676       12
 9   27.3     -1.31695789692       13
10   30.0     -1.19999982651       14
11   33.524   -1.10000635195       15
12   36.822   -1.00361552815       16
13   40.261   -0.916007799535      17
14   43.614   -0.822472442947      18
15   47.436   -0.72264587494       19
16   51.79    -0.616250691646      20
17   56.735   -0.503273260393      21
18   62.31    -0.384075299436      22
19   68.516   -0.259479460739      23
20   75.297   -0.130817437415      24
21   82.526    1.11022302463e-16   25


 i   theta    pseudorapidity     iEta
--------------------------------------
 0   3.0       3.00008274261      47
 1   5.7       2.61134904017      46
 2   8.4       2.33129451223      45
 3   11.1      2.1118548488       44
 4   13.8      1.93106912741      43
 5   16.5      1.77704423284      42
 6   19.2      1.64260787639      41
 7   21.9      1.52309848418      40
 8   24.6      1.41531379676      39
 9   27.3      1.31695789692      38
10   30.0      1.19999982651      37
11   33.524    1.10000635195      36
12   36.822    1.00361552815      35
13   40.261    0.916007799535     34
14   43.614    0.822472442947     33
15   47.436    0.72264587494      32
16   51.79     0.616250691646     31
17   56.735    0.503273260393     30
18   62.31     0.384075299436     29
19   68.516    0.259479460739     28
20   75.297    0.130817437415     27
21   82.526    1.1102230246e-16   26
*****/

#endif