#include "PamVMCNDDig.h"
#include <TRandom.h>
#include <TMath.h>


#define DEBUG 0

using TMath::Log;
using TMath::Exp;
using TMath::Abs;

ClassImp(PamVMCNDDig)

void PamVMCNDDig::Digitize(){

  if(DEBUG)  cout<<"Starting ND Digitization...."<<endl;

  Int_t NdN = 0;

  const Double_t Showercut = 10.e-6; // 10 us
  const Double_t Window = 200.e-6; //200 us
  const Double_t Tau = 17.e-6; // 17 us. Time when amplitude of signal reduce
                               // in e times
  const Double_t SigmaTau = 4.e-6; //Sigma of fall time (obtained from 
                                   //Lebedev's oscillograms)
  const Double_t Tresh = 130; //According Lebedev th = 130 mV
  const Double_t keVtomV = Tresh/200.; //We suppose that 130 mV = 200 keV... everything linear
  

  Double_t deadtime[36];
  Double_t Time,U, Uback;
  Int_t pos;

  for (Int_t i = 0; i<36; i++) deadtime[i] = Showercut;

  TClonesArray * hc = (TClonesArray*)fhitscolmap.GetValue("NDTI");

  
  PamVMCDetectorHit * hit = 0;
  PamVMCDetectorHit * hit1 = 0;
  PamVMCDetectorHit *hit0 =0;
   

  if(hc){
    
    Int_t HitsNum = hc->GetEntriesFast();

    //sorting by time and put results into vector fpHits
    if (HitsNum > 1){
      for(Int_t i=0; i<hc->GetEntriesFast(); i++){
	Time = 0.;
	for(Int_t j=0; j<hc->GetEntriesFast(); j++){
	  hit = (PamVMCDetectorHit*)hc->At(j);
	  if((hit->GetTOF()>Time) && (!IsInsideVector(hit))){
	    hit1 = hit;
	    Time = hit->GetTOF();
	  }
	}
	
	fpHits.push_back(hit1);
      }
    }
    
    //now we separate hits by tubes

    for(Int_t i = fpHits.size()-1; i>=0; i--){
      hit = (PamVMCDetectorHit*)fpHits.at(i);
      pos = hit->GetPOS()-1;
      tubeHits[pos].push_back(hit);
    }
    fpHits.clear();// we don't need it

    //From this moment all hits sorted by time and present in tubes

    for(Int_t i = 0; i<NTUBES; i++){
      if( tubeHits[i].size() > 0 ){
	for(UInt_t j = 0; j<tubeHits[i].size(); j++){
	    Uback = 0.;
	    hit =(PamVMCDetectorHit*)tubeHits[i].at(j);
	    U = keVtomV*(hit->GetEREL())*1.e6; //initial amplitude of this signal;
	    Time = hit->GetTOF();
	    for(UInt_t k = 1; k<=j; k++){
	      hit0 = (PamVMCDetectorHit*)tubeHits[i].at(k-1);
	      //hit1 = (PamVMCDetectorHit*)tubeHits[i].at(j);
	      if (tauHits[i].at(k-1)> 0)
		Uback += keVtomV*(hit0->GetEREL())*1.e6*Exp((hit0->GetTOF()-/*hit1->GetTOF()*/hit->GetTOF())/Abs(tauHits[i].at(k-1)));	  
	    }
	      
	    tauHits[i].push_back(frandom->Gaus(Tau,SigmaTau));
	    if (DEBUG){
	      cout<<"++++++HIT+++++++"<<endl;
	      cout<<"TubeNo="<<hit->GetPOS()<<endl
		  <<"Particle="<<((TParticlePDG*)(TDatabasePDG::Instance()->GetParticle(hit->GetPDG())))->GetName()<<endl
		  <<"    Erel="<<hit->GetEREL()*1.0e6<<" (keV)"<<endl
		  <<"    Time="<<hit->GetTOF()*1.e6<<" (us)"<<endl
		  <<"    Patn="<<hit->GetPATH()*10.<<" (mm)"<<endl
		  <<"     Tau="<<tauHits[i].at(j)*1.e6<<" (us)"<<endl
		  <<"       U="<<U<<" (mV)"<<endl
		  <<"    Uback="<<Uback<<" (mV)"<<endl;
	      if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut))
		cout<<" Particle detected as NEUTRON"<<endl;
	    }
	    if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut)) NdN++;

	  } 
      }
    }

   if(DEBUG)
     cout<<"==== Total nuber of detected neutrons for for event ==="<<NdN<<" ==="<<endl;

      }  
    
	  
    
  //Cleaning of everything...
  for (Int_t i = 0; i<NTUBES; i++){
    tauHits[i].clear();
    tubeHits[i].clear();
  }
  
  fData.clear();
  
  if (DEBUG) cout<<"REGISTERED: "<< NdN <<" neutrons"<<endl;
  for(Int_t i=0;i<3;i++){
    fData.push_back(0x0000);
    fData.push_back(0x010F);
  }
  fData.at(0) = ((0xFF00 & (256*NdN)));
}