#include "PamVMCS4Dig.h"
#include <TRandom.h>

ClassImp(PamVMCS4Dig)


void PamVMCS4Dig::Digitize(){

#ifdef DIG_DEBUG
  cout<<"Starting S4 Digitization...."<<endl;
#endif
  //This is simple approximation, reflection + atennuation..
  //I hope in future we will use for generation and propagation of 
  //optical photons abilitiees of montecarlo...




  Int_t DEBUG=0;
  // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
  Int_t i,j,t,NdF,pmt,S4;
  Float_t E0,X,Y,Z,x,y,z,V[3],q,w,l;

  const Float_t Xs[2] = {-24.1, 24.1};
  const Float_t Ys[2] = {-24.1, 24.1};
  const Float_t Zs[2] = {-0.5, 0.5};
  const Float_t Yp[6] = {-20., -17., -1., 2., 17., 20.}; //PMT positions

  const Float_t E1 = 1e-6; //Energy which we spend to produce 1 photon
  const Float_t p = 0.1;
  const Int_t l0 = 500;

  const Int_t S4p = 32; //pedestal

  UShort_t S4v = 0; // ADC count

  Int_t NdFT = 0;
  
  Float_t Ert = 0.;

  TClonesArray * hc = (TClonesArray*)fhitscolmap.GetValue("S4");
  PamVMCDetectorHit * hit = 0;

  if(hc){
    for(i=0;i<hc->GetEntriesFast();i++){

      hit = (PamVMCDetectorHit*)hc->At(i);

      Ert += hit->GetEREL();
      NdF=Int_t(hit->GetEREL()/E1);
      NdFT=0;
      X=hit->GetXIN();
      Y=hit->GetYIN();
      Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
      for(j=0;j<NdF;j++){
	q=(Float_t)random()/(Float_t)0x7fffffff;
	w=(Float_t)random()/(Float_t)0x7fffffff;
	V[0]=p*cos(6.28318*q);
	V[1]=p*sin(6.28318*q);
	V[2]=p*(2.*w-1.);
	pmt=0;
	x=X;
	y=Y;
	z=Z;
	while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
	  l=0;
	  while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
	    x+=V[0];
	    y+=V[1];
	    z+=V[2];
	    l+=p;
	  }
	  if((x<Xs[0]+p || x>Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p)){
	    for(t=0;t<3;t++){
	      if(y>=Yp[2*t] && y<Yp[2*t+1]){
		if(pmt==0)NdFT++;
		pmt=1;
		break;
	      }
	    }
	    if(pmt==1)break;
	    V[0]=-V[0];
	  }
	  q=(Float_t)random()/(Float_t)0x7fffffff;
	  w=1-exp(-l/l0);
	  if(q<w)break;
	  q=(Float_t)random()/(Float_t)0x7fffffff;
	  w=0.5;
	  if(q<w)break;
	  if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];
	  if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];
	  x+=V[0];
	  y+=V[1];
	  z+=V[2];
	  l=0;
	}
      }
    }
    Ert=Ert/0.002;
    q=(Float_t)(random())/(Float_t)0x7fffffff;
    w=0.7;
    E0=4064./7.;
    if(Ert<1) S4=0;
    else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
    i=S4/4;
    if(S4%4==0)
      S4v=S4+S4p;
    else if(S4%4==1){
      if(q<w) S4v=S4-1+S4p;
      else S4v=S4+1+S4p;
    } else if(S4%4==2) S4v=S4+S4p;
    else if(S4%4==3){
      if(q<w) S4v=S4+1+S4p;
      else S4v=S4-1+S4p;
    }
    if (DEBUG)
      cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
    
    fData.clear();
    
    fData.push_back(S4v);
    fData.push_back(0xd800);
    fData.push_back(0x0300); 
  }
}