// Includes from TRD
#include "CbmTrdRawToDigiR.h"
#include "CbmTrdDigi.h"
#include "TMath.h"
#include "CbmTrdAddress.h"
//#include "CbmSpadicRawMessage22.h"
#include<string>


#include "assert.h"
#include <iostream>
#include <chrono>
#include <vector>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TFile.h>

using namespace std::chrono; 

//_________________________________________________________________________________
CbmTrdRawToDigiR::CbmTrdRawToDigiR()
  :TObject()
  ,fSampleMask()
  ,fElookupSmall()
  ,fElookupAsym()
  ,fElookupA()
  ,fElookupBig()
  ,fElookup() 
{
}

//_________________________________________________________________________________
CbmTrdRawToDigiR::CbmTrdRawToDigiR(std::string readfile)
  :TObject()
  ,fSampleMask()
  ,fElookupSmall()
  ,fElookupAsym()
  ,fElookupA()
  ,fElookupBig()
  ,fElookup() 
{
  SetReadFile(readfile);  
}


//_________________________________________________________________________________
CbmTrdRawToDigiR::CbmTrdRawToDigiR(Double_t cal,Double_t tau,Int_t mode)
  :TObject()
  ,fSampleMask()
  ,fElookupSmall()
  ,fElookupAsym()
  ,fElookupA()
  ,fElookupBig()
  ,fElookup() 
{
  SetCalibration(cal);
  SetTau(tau);
  SetRecoMode(mode);
  SetPars(mode,cal,tau,fSampleMask);
}



//_________________________________________________________________________________
CbmTrdRawToDigiR::CbmTrdRawToDigiR(Double_t cal,Double_t tau,std::vector<Int_t> mask)
  :TObject()
  ,fSampleMask()
  ,fElookupSmall()
  ,fElookupAsym()
  ,fElookupA()
  ,fElookupBig()
  ,fElookup() 
{
  SetCalibration(cal);
  SetTau(tau);
  SetRecoMask(mask);
  SetPars(fRecoMode,cal,tau,mask);
}


void CbmTrdRawToDigiR::SetPars(Int_t mode,Double_t cal,Double_t tau,std::vector<Int_t> mask){

  //default
  if(mode==0){
    for(UInt_t i=0;i<=mask.size();i++) fSampleMask.push_back(mask[i]);
    fCalibration = cal;
    fTau = tau;
    Double_t sum=0;
    for(UInt_t i=0;i<mask.size();i++) sum+=fCalibration*CalcResponse(mask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC));
    fEReco     = sum;
    FillLookUps();
  }
  if(mode==1){
    for(Int_t i=fMaxBin;i<=fMaxBin;i++) fSampleMask.push_back(i);
    fCalibration = cal;
    fTau = tau;
    fLookUp = 0;
    Double_t sum=0;
    for(UInt_t i=0;i<fSampleMask.size();i++) sum+=fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC));
    fEReco     = sum;
    auto start = high_resolution_clock::now(); 
    FillLookUps();
    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    if(fDebug)   std::cout<<std::endl<<std::endl<<"filling took: "<<duration.count()<<"  microseconds "<<std::endl<<fEReco<<std::endl<<std::endl<<std::endl;
  }
  if(mode==2){
    Double_t sum=0;
    for(UInt_t i=0;i<fSampleMask.size();i++) sum+=fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC));
    fEReco     = sum;
    auto start = high_resolution_clock::now(); 
    FillLookUps();
    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    if(fDebug)   std::cout<<std::endl<<std::endl<<"filling took: "<<duration.count()<<"  microseconds "<<std::endl<<fEReco<<std::endl<<std::endl<<std::endl;
  }
  
}


void CbmTrdRawToDigiR::Init(){
  //default
  Double_t sum=0;
  for(UInt_t i=0;i<fSampleMask.size();i++) sum+=fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC));
  fEReco     = sum;
  auto start = high_resolution_clock::now(); 
  if(fReadFile=="")   FillLookUps(fWriteFile);  
  else                ReadMaps(fReadFile);
  if(fLookUp==4)      FillLookUps("");  
  auto stop = high_resolution_clock::now();
  auto duration = duration_cast<microseconds>(stop - start);
  if(fDebug)  std::cout<<std::endl<<std::endl<<"filling took: "<<duration.count()<<"  microseconds "<<std::endl<<fEReco<<std::endl<<std::endl<<std::endl;
  // if(fShapingOrder == 1)     {SetMaxBin(2+fPresamples);SetMinBin(1+fPresamples);}
  // if(fShapingOrder == 2)     SetMaxBin(4+fPresamples);
  
}


void CbmTrdRawToDigiR::FillLookUps(std::string write){
  
  if(fLookUp == 1){
    for(Int_t shift=0.;shift<CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC);shift++){      
      Double_t sum=0;
      for(Int_t i=fMinBin;i<=fMaxBin;i++) sum+=fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
      fEReco     = sum;
      Float_t temp = fCalibration*CalcResponse(fMaxBin*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
	
      for(Int_t max=0;max<fDynamicRange;max++){
	Float_t energy = max*1.0/temp;
	Int_t a = energy * fCalibration*CalcResponse(fMinBin*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
	Int_t b = energy * fCalibration*CalcResponse((fHighBin)*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);

	Int_t mina=a-fExtrapolate*a;
	Int_t maxa=a+fExtrapolate*a;
	Int_t minb=b-fExtrapolate*b;
	Int_t maxb=b+fExtrapolate*b;
	if(!(fElookupAsym[max][a][b] > 0)) fElookupAsym[max][a][b]=shift;
	fElookupSmall[shift][max]=energy;
	if(fDebug)  fQA->CreateHist("Asym",512,-12.0,500.0,512,-12.0,500.0);
	if(fDebug && max == 200 && fQA->GetCont2D("Asym",a,b) == 0. )  fQA->Fill("Asym",a,b,shift);
	if(fDebug)  fQA->CreateHist("Look",63,-0.5,62.5,512,-12.0,500.0);
	if(fDebug)  fQA->Fill("Look",shift,max,energy);
	Int_t extrapolate=0;
	while(true){
	  extrapolate++;
	  if((b - extrapolate) <= minb || (a - extrapolate) <= mina) break;
	  if(fElookupAsym[max][a-extrapolate][b-extrapolate] > 0) continue;
	  fElookupAsym[max][a-extrapolate][b-extrapolate]=fElookupAsym[max][a][b];
	  if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate,b-extrapolate,fElookupAsym[max][a][b]);
	  Int_t count=1;
	  if(!(fElookupAsym[max][a-extrapolate-count][b-extrapolate] > 0)) {
	    fElookupAsym[max][a-extrapolate-count][b-extrapolate]=fElookupAsym[max][a-extrapolate][b-extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate-count,b-extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	  count=1;
	  if(!(fElookupAsym[max][a-extrapolate+count][b-extrapolate] > 0)) {
	    fElookupAsym[max][a-extrapolate+count][b-extrapolate]=fElookupAsym[max][a-extrapolate][b-extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate+count,b-extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	}
	extrapolate=0;
	while(true){
	  extrapolate++;
	  if((b + extrapolate) >= maxb || (a + extrapolate) >= maxa) break;
	  if(fElookupAsym[max][a+extrapolate][b+extrapolate] > 0) continue;
	  fElookupAsym[max][a+extrapolate][b+extrapolate]=fElookupAsym[max][a][b];
	  if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate,b+extrapolate,fElookupAsym[max][a][b]);
	  Int_t count=1;
	  if(!(fElookupAsym[max][a+extrapolate-count][b+extrapolate] > 0)) {
	    fElookupAsym[max][a+extrapolate-count][b+extrapolate]=fElookupAsym[max][a+extrapolate][b+extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate-count,b+extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	  count=1;
	  if(!(fElookupAsym[max][a+extrapolate+count][b+extrapolate] > 0)) {
	    fElookupAsym[max][a+extrapolate+count][b+extrapolate]=fElookupAsym[max][a+extrapolate][b+extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate+count,b+extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	}
	
      }
    }
    if(write != "") WriteMaps(write);
  }  
  if(fLookUp == 2){
    for(Int_t shift=0.;shift<CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC);shift++){      
      Double_t sum=0;
      for(Int_t i=fMinBin;i<=fMaxBin;i++) sum+=fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
      fEReco     = sum;
      Float_t temp = fCalibration*CalcResponse(fMaxBin*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
      Int_t extrapolate=0;	
      for(Int_t max=0;max<fDynamicRange;max++){
	Float_t energy = max*1.0/temp;
	Int_t a = energy * fCalibration*CalcResponse(fMinBin*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
	Int_t b = energy * fCalibration*CalcResponse((fHighBin)*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
		
	Int_t mina=a-fExtrapolate*a;
	Int_t maxa=a+fExtrapolate*a;
	Int_t minb=b-fExtrapolate*b;
	Int_t maxb=b+fExtrapolate*b;
	if(!(fElookupAsym[max][a][b] > 0)) fElookupAsym[max][a][b]=shift;
	sum =0.;
	for(UInt_t i=0;i<fSampleMask.size();i++)	sum+=energy * fCalibration*CalcResponse(fSampleMask[i]*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)+shift);
	
	fElookupSmall[shift][sum]=energy;
	if(fDebug)  fQA->CreateHist("Asym",512,-12.0,500.0,512,-12.0,500.0);
	if(fDebug && max == 200 && fQA->GetCont2D("Asym",a,b) == 0. )  fQA->Fill("Asym",a,b,shift);
	if(fDebug)  fQA->CreateHist("Look",63,-0.5,62.5,512,-12.0,2000.0);
	if(fDebug)  fQA->Fill("Look",shift,sum,energy);	
	extrapolate=0;
	while(true){
	  extrapolate++;
	  if((b - extrapolate) <= minb || (a - extrapolate) <= mina) break;
	  if(fElookupAsym[max][a-extrapolate][b-extrapolate] > 0) continue;
	  fElookupAsym[max][a-extrapolate][b-extrapolate]=fElookupAsym[max][a][b];
	  if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate,b-extrapolate,fElookupAsym[max][a][b]);
	  Int_t count=1;
	  if(!(fElookupAsym[max][a-extrapolate-count][b-extrapolate] > 0)) {
	    fElookupAsym[max][a-extrapolate-count][b-extrapolate]=fElookupAsym[max][a-extrapolate][b-extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate-count,b-extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	  count=1;
	  if(!(fElookupAsym[max][a-extrapolate+count][b-extrapolate] > 0)) {
	    fElookupAsym[max][a-extrapolate+count][b-extrapolate]=fElookupAsym[max][a-extrapolate][b-extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a-extrapolate+count,b-extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	}
	extrapolate=0;
	while(true){
	  extrapolate++;
	  if((b + extrapolate) >= maxb || (a + extrapolate) >= maxa) break;
	  if(fElookupAsym[max][a+extrapolate][b+extrapolate] > 0) continue;
	  fElookupAsym[max][a+extrapolate][b+extrapolate]=fElookupAsym[max][a][b];
	  if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate,b+extrapolate,fElookupAsym[max][a][b]);
	  Int_t count=1;
	  if(!(fElookupAsym[max][a+extrapolate-count][b+extrapolate] > 0)) {
	    fElookupAsym[max][a+extrapolate-count][b+extrapolate]=fElookupAsym[max][a+extrapolate][b+extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate-count,b+extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	  count=1;
	  if(!(fElookupAsym[max][a+extrapolate+count][b+extrapolate] > 0)) {
	    fElookupAsym[max][a+extrapolate+count][b+extrapolate]=fElookupAsym[max][a+extrapolate][b+extrapolate];
	    if(fDebug && max == 200)  fQA->Fill("Asym",a+extrapolate+count,b+extrapolate,fElookupAsym[max][a][b]);
	    count++;
	  }
	}
	
      }

    }
    if(write != "") WriteMaps(write);
  }
  if(fLookUp == 4){
    for(Int_t max=0;max<fDynamicRange;max++){
      Float_t energy = max / (fCalibration*CalcResponse(fMaxBin*CbmTrdDigi::Clk(CbmTrdDigi::kSPADIC)));	
      fElookup[max]=energy;
    }
  }
  else{
    Int_t range=  fDynamicRange * fSampleMask.size();
    for(Int_t n=0;n<=range;n++){
      Float_t value = n/fEReco;
      fElookup[n]=value;
    }
  }
    
}


void CbmTrdRawToDigiR::ReadMaps(std::string file){
  if(fLookUp == 3){
    TFile f(file.data(), "OPEN");
    TProfile2D *h;
    h = (TProfile2D*)  f.Get("MAX ADC");
    for(Int_t x=1;x<=h->GetNbinsX();x++){
      for(Int_t y=1;y<=h->GetNbinsY();y++){
	fElookupSmall[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)]=h->GetBinContent(x,y);
	if(fDebug)  fQA->CreateHist("Read Small",63,-0.5,62.5,512,-12.,500.);
	if(fDebug)  fQA->Fill("Read Small",h->GetXaxis()->GetBinCenter(x),h->GetYaxis()->GetBinCenter(y),h->GetBinContent(x,y));
      }
    }
    h = (TProfile2D*)  f.Get("ASYM MAP");
    for(Int_t x=1;x<=h->GetNbinsX();x++){
      for(Int_t y=1;y<=h->GetNbinsY();y++){
	fElookupA[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)]=h->GetBinContent(x,y);
	if(fDebug)  fQA->CreateHist("Read Asym",512,0.,512.,512,-12.,500.);
	if(fDebug)  fQA->Fill("Read Asym",h->GetXaxis()->GetBinCenter(x),h->GetYaxis()->GetBinCenter(y),h->GetBinContent(x,y));
      }
    }
    f.Close();
  }
  if(fLookUp == 4){
    TFile f(file.data(), "OPEN");
    TProfile2D *h;
    h = (TProfile2D*)  f.Get("MAX ADC");
    for(Int_t x=1;x<=h->GetNbinsX();x++){
      for(Int_t y=1;y<=h->GetNbinsY();y++){
  	fElookupSmall[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)]=h->GetBinContent(x,y);
  	if(fDebug)  fQA->CreateHist("Read Small",63,-0.5,62.5,512,-12.,500.);
  	if(fDebug)  fQA->Fill("Read Small",h->GetXaxis()->GetBinCenter(x),h->GetYaxis()->GetBinCenter(y),h->GetBinContent(x,y));
      }
    }
    h = (TProfile2D*)  f.Get("ASYM MAP");
    for(Int_t x=1;x<=h->GetNbinsX();x++){
      for(Int_t y=1;y<=h->GetNbinsY();y++){
  	fElookupA[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)]=h->GetBinContent(x,y);
  	if(fDebug)  fQA->CreateHist("Read Asym",512,0.,512.,512,-12.,500.);
  	if(fDebug)  fQA->Fill("Read Asym",h->GetXaxis()->GetBinCenter(x),h->GetYaxis()->GetBinCenter(y),h->GetBinContent(x,y));
      }
    }
    f.Close();    
  }
  
}

CbmTrdDigi* CbmTrdRawToDigiR::MakeDigi(std::vector<Int_t> samples, Int_t channel,Int_t module,Int_t layer,ULong64_t time,Bool_t FN){
  Float_t digicharge=0;
  Int_t samplesum=0;
  for(size_t i=0;i<fSampleMask.size();i++){
    samplesum+=samples[fSampleMask[i]];    
  }   
  if(fReadFile == ""){
    if(fLookUp == 1)   {
      digicharge=fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samples[fMaxBin]];
    }
    if(fLookUp == 2)   {
      digicharge=fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samplesum];
    }
    else{
      digicharge=fElookup[samplesum];
    }
  }
  else{
    if(fLookUp == 3)   {
      digicharge=fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
    }
    if(fLookUp == 4 && !FN)   {
      digicharge=fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
    }
    if(fLookUp == 4 && FN)   {
      digicharge=fElookup[samples[fMaxBin]];
    }
    
  }
  
  CbmTrdDigi* digi= new CbmTrdDigi(channel, digicharge, time, 0, 0); 
  digi->SetAddress(CbmTrdAddress::GetAddress(layer, module, 0, 0, 0));
  
  return digi;
  
}

Float_t CbmTrdRawToDigiR::GetTimeShift(std::vector<Int_t> samples){

  Float_t shift=0;
  if(fReadFile == ""){
    if(fLookUp == 1)   {
      shift=fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]];
    }
    else{
      return 0.;
    }
  }
  else{
    if(fLookUp == 3)   {
      shift=fElookupA[samples[fMinBin]][samples[fMaxBin]];
    }
    if(fLookUp == 4)   {
      shift=fElookupA[samples[fMinBin]][samples[fMaxBin]];
    }
  }
    
  return shift;  
}

Double_t CbmTrdRawToDigiR::GetCharge(std::vector<Int_t> samples, Int_t shift){

  Float_t charge=0;
  if(fReadFile == ""){
    if(fLookUp == 1)   {      
      if(shift > -1) charge=fElookupSmall[shift][samples[fMaxBin]];
      else           charge=fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samples[fMaxBin]];
    }
    if(fLookUp == 2)   {      
      Int_t samplesum=0;
      for(size_t i=0;i<fSampleMask.size();i++){
	samplesum+=samples[fSampleMask[i]];    
      }   
      if(shift > -1) charge=fElookupSmall[shift][samples[fMaxBin]];
      else           charge=fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samples[fMaxBin]];
    }
    else{
      return 0.;
    }
  }
  else{
    if(fLookUp == 3)   {      
      if(shift > -1) charge=fElookupSmall[shift][samples[fMaxBin]];
      else           charge=fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
    }
    if(fLookUp == 4)   {      
      if(shift > -1) charge=fElookupSmall[shift][samples[fMaxBin]];
      else           charge=fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
    }
  }
  
  return charge;  
}


//_______________________________________________________________________________
Double_t CbmTrdRawToDigiR::CalcResponse(Double_t t)
{

  if(fShapingOrder == 1)     return (t/fTau)*TMath::Exp(-(t/fTau));
  if(fShapingOrder == 2)     return (t/fTau)*(t/fTau)*TMath::Exp(-(t/fTau));

  return 0.;
}