// simulation of induced charge from a track in ATLAS SCT barrel strips
// y = 0 for HV plane, y = d (285 microns)for strip plane
// x = 0 to 80 microns with strip at the center
// 2010.8.3 Taka Kondo (KEK)
//
#include "SCT1/MySCTtestmap.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGauss.h"
#include <TRandom.h>

////////////////////////////////////////////////////////

MySCTtestmap::MySCTtestmap(const std::string& name, ISvcLocator* pSvcLocator):
  AthAlgorithm(name, pSvcLocator)
{
  //---------------setting basic parameters---------------------------
   declareProperty("X_in", m_x = 50.0) ;
   declareProperty("Y_in", m_y = 50.0) ;
}

////////////////////////////////////////////////////////////////
MySCTtestmap::~MySCTtestmap() {}

//////////////////////////////////////////////////////////////////////////////
StatusCode MySCTtestmap::initialize() {

//////////////////////////////////////////////////////////////
//
  init_inducedChargesArray();
// 
/** get a handle on the NTuple and histogramming service */
  m_sc = service("THistSvc", m_thistSvc);

m_htm2= new TH1F("m_htm2","strip=-2, induced",50, 0., 25.);
m_sc = m_thistSvc->regHist("/AANT/m_htm2",m_htm2);

m_htm1= new TH1F("m_htm1","strip=-1, induced",50, 0., 25.);
m_sc = m_thistSvc->regHist("/AANT/m_htm1",m_htm1);

m_ht00= new TH1F("m_ht00","strip=0, induced",50, 0., 25.);
m_sc = m_thistSvc->regHist("/AANT/m_ht00",m_ht00);

m_htp1= new TH1F("m_htp1","strip=1, induced",50, 0., 25.);
m_sc = m_thistSvc->regHist("/AANT/m_htp1",m_htp1);

m_htp2= new TH1F("m_htp2","strip=2, induced",50, 0., 25.);
m_sc = m_thistSvc->regHist("/AANT/m_htp2",m_htp2);
   

return m_sc;
}
 
///////////////////////////////////////////////////////////////////

StatusCode MySCTtestmap::execute() {


double charge[50];
for (int istrip=-2 ; istrip<3; istrip++) {
     inducedCharges (istrip,  m_x,  m_y, charge) ;
     for (int it = 0; it<50; it++) {
       double time = it*0.5 + 0.25 ;
       double q =  charge[it];
       if(istrip== 0) m_ht00->Fill(time, q); 
       if(istrip==-1) m_htm1->Fill(time, q);
       if(istrip==-2) m_htm2->Fill(time, q);
       if(istrip== 1) m_htp1->Fill(time, q);
       if(istrip== 2) m_htp2->Fill(time, q);
     }
}
  return StatusCode::SUCCESS;
}

/////////////////////////////////////////////////////////////////////////

StatusCode MySCTtestmap::finalize() {

///////////////////////////////////////////////////////////////////////////
return StatusCode::SUCCESS;
}


//--------------------------------------------------------------------
void MySCTtestmap::inducedCharges (int& istrip, double& x, double& y, double* temp) {
   double  fx, fy ;
   int ix, iy, ix1, iy1;
   if (x > 80. || x <  0. ) return ;
   if (y < 0. || y > 285. ) return ;
   fx = x / m_mappingStep ;
   ix = int(fx);
   ix1=ix+1;
   fx = fx - ix;
   fy = (y - 0.5*m_mappingStep) / m_mappingStep ;
   iy = int(fy);
   iy1=iy+1;
   fy = fy - iy;
   if (y <= 0.5*m_mappingStep) {fy = 0.; iy=0; iy1=0;}
   if (y >= (m_iymax-0.5)*m_mappingStep){fy=1.; iy=m_iymax-1; iy1=m_iymax-1;} 
   std::cout<<"ix,iy,fx,fy"<<ix<<" "<<iy<<" "<<fx<<" "<<fy<<std::endl;
   double p00[50],p01[50],p10[50],p11[50];
   inducedChargesArray(istrip, ix,  iy,  p00);
   inducedChargesArray(istrip, ix,  iy1, p01);
   inducedChargesArray(istrip, ix1, iy,  p10);
   inducedChargesArray(istrip, ix1, iy1, p11);
   for (int i=0; i<50; i++) {
     temp[i] = p00[i]*(1.-fx)*(1.-fy) + p01[i]*(1.-fx)*fy 
               + p10[i]*fx*(1.-fy) + p11[i]*fx*fy;
   }
   return;
}

void MySCTtestmap::inducedChargesArray (int& istrip, int& ix, int& iy, double* charge) {

for (int it=0; it<50; it++) {
  if(istrip==-2) charge[it] = m_stripCharge_m2[ix][iy][it];
  if(istrip==-1) charge[it] = m_stripCharge_m1[ix][iy][it];
  if(istrip== 0) charge[it] = m_stripCharge_00[ix][iy][it];
  if(istrip== 1) charge[it] = m_stripCharge_p1[ix][iy][it];
  if(istrip== 2) charge[it] = m_stripCharge_p2[ix][iy][it];
}
return ;
}
//--------------------------------------------------------------------
//  get Charge arrays
//-------------------------------------------------------------------
double GetMappingStep();
double GetCharge150_m2(int& i, int& j, int& k);
double GetCharge150_m1(int& i, int& j, int& k);
double GetCharge150_00(int& i, int& j, int& k);
double GetCharge150_p1(int& i, int& j, int& k);
double GetCharge150_p2(int& i, int& j, int& k);

void MySCTtestmap::init_inducedChargesArray ( ) {
     m_mappingStep = GetMappingStep();
     m_ixmax = int( 80./ m_mappingStep) + 1;
     m_iymax = int( 285. / m_mappingStep) ;
     std::cout<<"m_mappingStep, m_ixmax, m_iymax= "
     << m_mappingStep<< " "<< m_ixmax << " "<< m_iymax << std::endl;
        for (int k=0; k<m_ixmax; k++) {
           for (int l=0; l<m_iymax; l++) {
              for (int m=0; m<50; m++) {
              m_stripCharge_m2[k][l][m] = GetCharge150_m2(k,l,m);
              m_stripCharge_m1[k][l][m] = GetCharge150_m1(k,l,m);
              m_stripCharge_00[k][l][m] = GetCharge150_00(k,l,m);
              m_stripCharge_p1[k][l][m] = GetCharge150_p1(k,l,m);
              m_stripCharge_p2[k][l][m] = GetCharge150_p2(k,l,m);
              }
           }
        }
 
return;
}