///////////////////////////////////////////////////////////////////////////////
//  Dec 17 2000 P.Murat: simple ROOT-parsable description of CalTower
// 
//  Description of the inline functions:
//  --------------------------------------
//  
// - EmCOGPhi()  : phi center of gravity by EM compartment
// - EmPmt (i)   : ADC counts in EM  PMT channel 'i'
// - Energy()    : sum of EM and HAD energies
// - HadPmt(i)   : ADC counts in HAD PMT channel 'i'
// - HadCOGPhi() : phi center of gravity by HAD compartment
//
// next 2 functions currently (Apr 15) are defined only for the central
// calorimeter
//
// - HadTdc0() : CHA time for the central
// - HadTdc1() : WHA time for the overlap between the central and the wall
//
// - Theta  () : theta of the center of the tower
// - Type   () : returns type of the tower (1-9)
//
// Oct 13 2001: v2, pack all the floats into shorts, add ValidCode
// 
// fCode: val_code&0xff + n0(TDC hits) << 8  + n1(TDC hits) << 12
///////////////////////////////////////////////////////////////////////////////
#include <Stntuple/geom/TStnCalGeometry.hh>


ClassImp(TStnCalGeometry)

//_____________________________________________________________________________
TStnCalGeometry::TStnCalGeometry() {
  Init();
}

//_____________________________________________________________________________
TStnCalGeometry::~TStnCalGeometry() {
}


//_____________________________________________________________________________
Int_t TStnCalGeometry::IEta(Float_t DetEta) {
  // return tower IEta for a given detector eta

  int   ieta;

  for (ieta=0; ieta<TOWER_NETA; ieta++) {
    if ((DetEta > fEta[ieta]) && (DetEta <= fEta[ieta+1])) break;
  }
  if (ieta == TOWER_NETA) ieta = -1;
  return ieta;
}

//_____________________________________________________________________________
Int_t TStnCalGeometry::IEta(Float_t x, Float_t y, Float_t z) {

  float R     = sqrt(x*x + y*y);
  float cotTheta;
  if (R > 0.0 ){
	cotTheta = z/R;
  } else {
	return 0;
  }
  
  for(int i =0; i<TOWER_NETA; i++){ 
    //printf("%i %f %f %f \n",i,cotTheta,fCotMin[i],fCotMax[i]);
    if(cotTheta>fCotMin[i] && cotTheta<fCotMax[i]) return i;
  }

}

//_____________________________________________________________________________
Double_t TStnCalGeometry::Phi(Int_t IEta, Int_t IPhi) {
  return  2*PI * float(IPhi+0.5)/float(TEPSEG[TETTYP[IEta]]);
}

//_____________________________________________________________________________
Int_t TStnCalGeometry::IPhi(Int_t IEta, Float_t phi) {
  
  if (phi>2*PI) { phi -= 2.*PI*int ( phi/(2.*PI)); }
  if (phi<0)    { phi += 2.*PI*int((-phi/(2.*PI))+1); }
  return  int(phi/(2*PI) * TEPSEG[TETTYP[IEta]]);

}



//_____________________________________________________________________________
void TStnCalGeometry::Init() {

  float theta;

  for (int i=0; i<TOWER_NETA; i++) {
    if (i <= 25) {
      theta      = TMath::Pi()*(1-(TOWER_THETA[i]+TOWER_THETA[i+1])/2.0/180.);
      fTheta[i] = theta;
      fEta  [i] = -TMath::Log(fabs(TMath::Tan(theta/2.)));
      fCotMin[i] = 1.0/tan( (180.-TOWER_THETA[i]  )*DEGRAD );
      fCotMax[i] = 1.0/tan( (180.-TOWER_THETA[i+1])*DEGRAD );
    }
    else if (i > 25) {
      int i1 = 51-i;
      theta      = TMath::Pi()*((TOWER_THETA[i1]+TOWER_THETA[i1+1])/2.0/180.);
      fTheta[i] = theta;
      fEta  [i] = -TMath::Log(fabs(TMath::Tan(theta/2.)));
      fCotMax[i] = 1.0/tan( (TOWER_THETA[i1]  )*DEGRAD );
      fCotMin[i] = 1.0/tan( (TOWER_THETA[i1+1])*DEGRAD );
    }
  }

  fEta[TOWER_NETA] = -fEta[0];


  for (int is=0; is<2; is++) {
    for (int it=0; it<84; it++) {
      float mpx = MPA_XPositions[is][it];
      float mpy = MPA_YPositions[is][it];
      //--- theta and eta ---
      double rtow = sqrt(pow(mpx,2) + pow(mpy,2));
      theta  = atan(rtow/580.);           // Z=580cm: MP front plane
      fEtaMP  [is][it] = -1.0*log(tan(0.5*theta));
      fThetaMP[is][it] = theta*180./TMath::Pi();
      if (is==0) {
        theta  = TMath::Pi() - theta;
	fEtaMP  [is][it] = -1.0*log(tan(0.5*theta));
	fThetaMP[is][it] = theta*180./TMath::Pi();
      }
      
      //--- phi ---
      double dphi = atan(mpy/mpx);
      if (mpx<=0)
	dphi += TMath::Pi();
      else if (mpx>0&&mpy<=0)
	dphi += 2*TMath::Pi();
      fPhiMP[is][it] = dphi;
    }
  }

}

//_____________________________________________________________________________
Int_t    TStnCalGeometry::MPITow (int ieta, int iphi)  const {
  // return MP tower number (0-83) for a given ieta,iphi

  int layer = -1;
  if (ieta <= 3) layer = ieta;
  if (ieta >=48) layer = 51-ieta;
  
  int phi = iphi;
  int tower = -1;    
  if (layer == 3) {
    int wedge = phi / 5;
    int phiInWedge = phi % 5;
    tower = 9 + 14 * wedge + phiInWedge;
  }
  else if (layer == 2) {
    int wedge = phi / 4;
    int phiInWedge = phi % 4;
    tower = 5 + 14 * wedge + phiInWedge;
  }
  else if (layer == 1) {
    int wedge = phi / 3;
    int phiInWedge = phi % 3;
    tower = 2 + 14 * wedge + phiInWedge;
  }
  else if (layer == 0) {
    int wedge = phi / 2;
    int phiInWedge = phi % 2;
    tower = 14 * wedge + phiInWedge;
  }
  return tower;
}

//_____________________________________________________________________________
void TStnCalGeometry::Clear(Option_t* opt) {
}

//_____________________________________________________________________________
void TStnCalGeometry::Print(Option_t* opt) const {
}