#include "Digitizer.h"

using namespace std;

void Digitizer::DigitizeS4(){
  Int_t DEBUG=0;
  // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
  TString ciao,modo="ns";
  Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
  Float_t E0,E1=1e-6,Ert,X,Y,Z,x,y,z,V[3],Xs[2],Ys[2],Zs[2],Yp[6],q,w,p=0.1,l,l0=500;
  Xs[0]=-24.1;
  Xs[1]=24.1;
  Ys[0]=-24.1;
  Ys[1]=24.1;
  Zs[0]=-0.5;
  Zs[1]=0.5;
  Yp[0]=-20.;
  Yp[2]=-1.;
  Yp[4]=17.;
  for(i=0;i<3;i++)
    Yp[2*i+1]=Yp[2*i]+3;
  srand(time(NULL));
  // --- activate branches:
  //   fhBookTree->SetBranchStatus("Nthtof",1); //modified by E.Vannuccini 03/08
  //   fhBookTree->SetBranchStatus("Ipltof",1);
  //   fhBookTree->SetBranchStatus("Ipaddle",1);

  //   fhBookTree->SetBranchStatus("Xintof",1);
  //   fhBookTree->SetBranchStatus("Yintof",1);
  //   fhBookTree->SetBranchStatus("Xouttof",1);
  //   fhBookTree->SetBranchStatus("Youttof",1);

  //   fhBookTree->SetBranchStatus("Ereltof",1);
  //   fhBookTree->SetBranchStatus("Timetof",1);
  NdFT=0;
  Ert=0;
  for(i=0;i<Nthtof;i++){
    if(Ipltof[i]!=6) continue;
    Ert+=Ereltof[i];
    if(modo=="ns") continue;
    NdF=Int_t(Ereltof[i]/E1);
    NdFT=0;
    X=Xintof[i];
    Y=Yintof[i];
    Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
    //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
    for(j=0;j<NdF;j++){
      q=(Float_t)random()/(Float_t)0x7fffffff;
      w=(Float_t)random()/(Float_t)0x7fffffff;
      // cout<<"qw "<<q<<" "<<w<<endl;
      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;
          //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
          //cin>>ciao;
        }
        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;
              //cout<<NdFT<<endl;
              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;
        //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
	//cin>>ciao;
      }
    }
  }
  Ert=Ert/0.002;
  q=(Float_t)(random())/(Float_t)0x7fffffff;
  w=0.7;
  //E0=(Float_t)(4064./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;
  fDataS4[0]=S4v;//0xf028;
  fDataS4[1]=0xd800;
  fDataS4[2]=0x0300;
  //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
  //cin>>ciao;
}