#include "CbmTrdParModGas.h"

#include <FairLogger.h>         // for Logger, LOG

#include <TAxis.h>              // for TAxis
#include <TGenericClassInfo.h>  // for TGenericClassInfo
#include <TH2.h>                // for TH2F
#include <TMath.h>              // for Exp
#include <TMathBase.h>          // for Abs
#include <TObjArray.h>          // for TObjArray
#include <TObjString.h>         // for TObjString

#include <stdio.h>              // for printf

#define VERBOSE 0  

// binding energy in keV for 'K' 'L' and 'M' shells of Ar and Xe 
Float_t CbmTrdParModGas::fgkBindingEnergy[2][NSHELLS] = {
  {34.5,  5.1, 1.1}, // Xe
  { 3.2, 0.25, 0.1}  // Ar
};
// binding energy in keV for 'K' 'L' and 'M' shells of Ar and Xe 
Float_t CbmTrdParModGas::fgkBR[2][NSHELLS-1] = {
  {0.11, 0.02}, // Xe
  {0.11, 0.02}  // Ar
};
// average energy to produce one electron-ion pair
Float_t CbmTrdParModGas::fgkWi[3] = {
  23.1, // Xe
  25.8, // Ar 
  33    // CO2
};
// Bucharest detector gas gain parametrization based on 55Fe measurements with ArCO2(80/20) 
Float_t CbmTrdParModGas::fgkGGainUaPar[2]  = {-10.1676,  8.3745};
Float_t CbmTrdParModGas::fgkE0 = 866.1047;// energy offset in ADC ch @ 0 keV

//___________________________________________________________________
CbmTrdParModGas::CbmTrdParModGas(const char* title)
  : CbmTrdParMod("CbmTrdParModGas", title)
  ,fConfig(0)
  ,fUa(0)
  ,fUd(0)
  ,fDw(0.3)
  ,fGasThick(0.6)
  ,fPercentCO2(0.2)
  ,fDriftMap(nullptr)
  ,fFileNamePID()
{
  TString s(title);
  TString name; 
//  Int_t val;
  TObjArray *so = s.Tokenize("/");
  for(Int_t ie(0); ie<so->GetEntries(); ie+=2){
    name = ((TObjString*)(*so)[ie])->String();
    if(name.EqualTo("Module")) fModuleId = ((TObjString*)(*so)[ie+1])->String().Atoi();
    else if(name.EqualTo("Ua")) fUa = ((TObjString*)(*so)[ie+1])->String().Atoi();
    else if(name.EqualTo("Ud")) fUd = ((TObjString*)(*so)[ie+1])->String().Atoi();
    else if(name.EqualTo("Gas")){ 
      TString gas = ((TObjString*)(*so)[ie+1])->String();
      if(gas.EqualTo("Ar")) SetNobleGasType(1);
      else if(gas.EqualTo("Xe")) SetNobleGasType(0);
      else {
        LOG(warn) << GetName() <<":: gas type \""<<gas<<"\" not defined. Default to Xe." ;
        SetNobleGasType(0);
      } 
    }
  }
  so->Delete(); delete so;

  if(VERBOSE) printf("Module[%2d] U[%4d %3d]\n", fModuleId, fUa, fUd);
}

//___________________________________________________________________
CbmTrdParModGas::~CbmTrdParModGas()
{
//  if(fDriftMap) delete fDriftMap;
}  

//_______________________________________________________________________________________________
Float_t CbmTrdParModGas::EkevFC(Float_t ekev) const
{
  // Convert energy deposit to no of primary ionizations and apply gas gain.
  // Currently gas gain is evalauted from 55Fe spectrum analysis on ArCO2(80/20)  
  Int_t     gasId=GetNobleGasType()-1;
  Float_t   wi=(1.-fPercentCO2)*fgkWi[gasId]+fPercentCO2*fgkWi[2];
  
  //gas gain
  // G = G[ev->ADC] * wi[ArCO2]/C[mV->ADC]/A[fC->mV]/e
  // G[ev->ADC] : measured gain based on 55Fe spectrum (expo)
  // wi[ArCO2]  : average energy to produce a ele-ion pair in mixture (27.24 ev)
  // C[mV->ADC] : FASP out [2V] to ADC range [4096 ch] (2) 
  // A[fC->mV]  : FASP gain from CADENCE (6)
  // e : 1.6e-4 [fC] electric charge
  Double_t  gain = 170.25*TMath::Exp(fgkGGainUaPar[0]+fgkGGainUaPar[1]*fUa*1.e-3)/12.;
  // for Xe correct Ar gain measurements TODO
  if(gasId==0) gain *= 0.6;
  
  Float_t efC = gain*ekev*0.16 / wi;  
  if(VERBOSE) printf("        ua[V]=%d gain[%5.2e] wi[eV]=%5.2f :: E[keV]=%6.3f E[fC]=%6.2f\n", fUa, gain, wi, ekev, efC);
  return efC;
  
//   Double_t eadc = fgkE0 + gain * ekev;
//   //printf("E = %fkeV %fADC\n", ekev, eadc);  
//   // apply FASP gain -> should be done without intermediate ADC conversion TODO 
//   Double_t sFASP = eadc/2.; // FASP signal [mV]; FASP amplification 1
//   // FASP gaincharacteristic -> should be defined elsewhere
//   // data based on CADENCE simulations
//   Double_t s0FASP = 10, gFASP = 6;
//   return (sFASP-s0FASP)/gFASP;
}

//_______________________________________________________________________________________________
Int_t CbmTrdParModGas::GetShellId(const Char_t shell) const
{
/** Return index of atomic shell. 
 * shell name can be 'K', 'L' and 'M'
 */  
  switch(shell){
    case 'K': return 0;
    case 'L': return 1;
    case 'M': return 2;
    default:
      LOG(warn) << GetName() <<"::GetShellId: Atomic shell : "<<shell<<" not defined for gas "<<(GetNobleGasType()==2?"Ar":"Xe");
      return -1;
  }
}

//_______________________________________________________________________________________________
Float_t CbmTrdParModGas::GetBindingEnergy(const Char_t shell, Bool_t main) const
{
  Int_t gasId=GetNobleGasType()-1;
  Int_t shellId = GetShellId(shell);
  if(shellId<0) return 0;
  
  if(!main) return fgkBindingEnergy[gasId][shellId];
  else{
    if(shellId<NSHELLS-1) return fgkBindingEnergy[gasId][shellId+1];
    else{
      LOG(warn) << GetName() <<"::GetBindingEnergy: Request atomic shell : "<<shellId+1<<" not defined for gas "<<(gasId?"Ar":"Xe");
      return 0;
    }
  }
  return 0;
}

//_______________________________________________________________________________________________
Float_t CbmTrdParModGas::GetNonIonizingBR(const Char_t shell) const
{
  Int_t gasId=GetNobleGasType()-1;
  Int_t shellId = GetShellId(shell);
  if(shellId<0) return 0;

  return fgkBR[gasId][shellId];  
}

//_______________________________________________________________________________________________
Char_t CbmTrdParModGas::GetPEshell(Float_t Ex) const
{
  const Char_t shellName[NSHELLS] = {'K', 'L', 'M'};
  Int_t gasId=GetNobleGasType()-1;
  for(Int_t ishell(0); ishell<NSHELLS; ishell++){
    if(Ex<fgkBindingEnergy[gasId][ishell]) continue;
    return shellName[ishell];
  }  
  LOG(debug) << GetName() <<"::GetPEshell: Ex[keV] "<<Ex<<" less than highes atomic shell binding energy : "<<fgkBindingEnergy[gasId][NSHELLS-1]<<" for gas "<<(gasId?"Ar":"Xe");
  return 0;
}

//_______________________________________________________________________________________________
Double_t CbmTrdParModGas::GetDriftTime(Double_t y0, Double_t z0) const
{
  const TAxis *ay(fDriftMap->GetXaxis()), *az(fDriftMap->GetYaxis());
  Int_t by(ay->FindBin(y0)), 
        bz(az->FindBin(z0));
  Double_t tmin(fDriftMap->GetBinContent(by, bz));      
  if(VERBOSE) printf("GetDriftTime :: Start @ dt=%3d [ns]\n", Int_t(tmin));
  return tmin;
}

//_______________________________________________________________________________________________
void CbmTrdParModGas::Print(Option_t */*opt*/) const
{
  printf("%s @ %4d ", GetName(), fModuleId);
  printf("Type[%s] ", GetDetName());
  printf("%s[%4.1f%%] Ua[V]=%d Ud[V]=%d ", GetNobleGasName(), 1e2*GetNobleGas(), fUa, fUd);
  printf("Pid Type[%d] DB[%s]\n", GetPidType(), fFileNamePID.Data());
}

//_______________________________________________________________________________________________
Double_t CbmTrdParModGas::ScanDriftTime(Double_t y0, Double_t z0, Double_t dzdy, Double_t dy) const
{
  Double_t y1=y0+dy, z1=z0+dzdy*dy,
           dw(fDw), dwh(0.5*dw);
//  Double_t dhh(fGasThick);
  
  if(VERBOSE) printf("ScanDriftTime :: Try : [%7.4f %7.4f] => [%7.4f %7.4f] dzdy[%5.2f] dy[%5.2f]\n", y0, z0, y1, z1, dzdy, dy);
  while(y1<-dwh-1e-3){
    y0+=dw;y1+=dw;
  }
  while(y1>dwh+1.e-3){
    y0-=dw; y1-=dw;
  }
  TH2F *h = fDriftMap;  
  y1=y0; z1=z0;
  const TAxis *ay(fDriftMap->GetXaxis()), *az(fDriftMap->GetYaxis());
  Int_t by(ay->FindBin(y1)), 
        bz(az->FindBin(z1)),
        nby(ay->GetNbins()),
        nbz(az->GetNbins());
  Float_t dyStep = ay->GetBinWidth(1), 
          tmin(500), tmax(0), tc(0);
  while(by>0 && by<=nby && bz<=nbz){
    bz=az->FindBin(z1);
    tc=h->GetBinContent(by, bz);
    //if(VERBOSE) printf("ScanDriftTime :: Do  : y0(%7.4f), z0(%7.4f), by(%3d), bz(%3d) td[ns]=%5.1f  t0[ns]=%5.1f \n", y1, z1, by, bz, tc, tmin);
    
    if(tc>0){
      if(tc<tmin) tmin=tc;
      if(tc>tmax) tmax=tc;
    }
    z1+=TMath::Abs(dzdy)*dyStep;
    if(dzdy>0){ y1+=dyStep; by++;}
    else{       y1-=dyStep; by--;}
  }
  if(VERBOSE) printf("ScanDriftTime :: Start @ dt=%3d [ns]\n", Int_t(tmin));

  return tmin;
}

//___________________________________________________________________
void CbmTrdParModGas::SetDriftMap(TH2F *hm, TDirectory *d)
{
/**  
 * Load drift map in the module and get ownership
 */

  if(VERBOSE) printf("CbmTrdParModGas::SetDriftMap : Module[%2d] U[%4d %3d]\n", fModuleId, fUa, fUd);
  
  fDriftMap = (TH2F*)hm->Clone(Form("trdDM%02d", fModuleId));
  fDriftMap->SetTitle(GetTitle());
  fDriftMap->SetDirectory(d);
}

ClassImp(CbmTrdParModGas)