// ---------------------------------------------------------------------
// -----                  KrattaEne header file                    -----
// -----            Created 10/11/13  by J. Lukasik                -----
// ---------------------------------------------------------------------

/** KrattaEne
 * @author J. Lukasik <Jerzy.Lukasik@ifj.edu.pl>
 * @since 10/11/13
 *   Event energy calibration for the KRATTA: parameters 
 *   container and procedures
 *
 **/


#ifndef __KRATTA_ENE_CLASS__
#define __KRATTA_ENE_CLASS__


#include "TObject.h"
#include "Rtypes.h"
#include "TClass.h"
#include "TCutG.h"
#include "TRandom3.h"
#include "TFile.h"
#include "TString.h"

#include <fstream>
#include <sstream>
#include <iostream>
#include <string>

#include "re_tab.hxx"
#include "KRATTAgeo.hxx"
#include "KRATTAtgt.hxx"

#include "mktreedstchain_sct_fitABC.h"

#include "TRootKRATEvent.h"

class RE_TAB;
class Target;
class KrattaGeo;
class ASYEvent;
class ASYFadcPeak;
class TRootKRATEvent;
class TRootKRATParticle;



class KrattaEne : public TObject{
private:
  int     run_current;
  TString mat_current;
public:

  KrattaEne();
  virtual ~KrattaEne();

  //friend std::ostream& operator<<(std::ostream&, const KrattaEne& p);
  
  //TRootKRATParticle fDat;
  //TRootKRATParticle* GetTRootKRATParticle(void ) { return &fDat;};
  
  void LoadParameters(Int_t run);
  void LoadParametersOld(Int_t run);
  void LoadFineCorrections(Int_t run);
  void LoadPD01PthrChannels(Int_t run);
  void LoadGammaLines(Int_t run);
  void LoadStarCuts(Int_t run);
  void LoadExtraProtonCuts(Int_t run);
  void LoadTime0Intervals(Int_t run);
  
  bool IsInTime(Int_t mod, Int_t ipd, Float_t time, Float_t nsigma_minus=0, Float_t nsigma_plus=0);
  Float_t DiffTime(Int_t mod, Int_t ipd, Float_t time);//returns the difference
  //between time-<time> in the units of sigma
  Float_t GetTime0Mean(Int_t mod, Int_t PDn){
  return (PDn>=0&&PDn<=2&&mod>=0&&mod<35) ? time0[mod][PDn][0] : -1;}
  Float_t GetTime0Sigma(Int_t mod, Int_t PDn){
  return (PDn>=0&&PDn<=2&&mod>=0&&mod<35) ? time0[mod][PDn][1] : -1;}
  
  //void Calibrate(Int_t run, ASYFadcPeak *peak, Int_t iz, Int_t ia);
  ///void Calibrate(Int_t run, ASYFadcPeak *peak, Int_t iz, Int_t ia, TRootKRATEvent& kevt, TRootKRATParticle& dat );
  ///void Calibrate(Int_t run, ASYFadcPeak *peak, Int_t iz, Int_t ia, TRootKRATEvent *kevt, TRootKRATParticle* dat );
  void Calibrate(  ASYFadcPeak *peak, TRootKRATEvent& kevt, TRootKRATParticle& dat);
  void Calibrate(  ASYFadcPeak *peak, TRootKRATEvent *kevt, TRootKRATParticle* dat);
  
  //void EfficencyCorrection( TRootKRATParticle* dat );
  float EfficencyCorrection( int iz, int ia, float E_kin_in_MeV );
  
  static const Int_t NNPAR=16; // number of calibration parameters (new)
  static const Int_t NPAR=20; // number of calibration parameters
  static const Int_t NMOD=35; // number of modules
  //static const Float_t cpt_fct;//    = 1.01; // skip particles 1% above pthr
  //static const Float_t uball_thck;// = 7000; // effective CsI uball thick. of ring 5
  Float_t cpt_fct;
  Float_t uball_thck;
  
  float par[NMOD][NPAR];
  float npar[NMOD][NNPAR];
  float nparaz[NMOD][NNPAR][5];
  float gamma1[NMOD][3];  // parabolic gamma cuts for CSI1
  float gamma2[NMOD][3];  // parabolic gamma cuts for CSI2
  TCutG *cutgstar1[NMOD]; // punch-through and background cut for CSI1
  TCutG *cutgstar2[NMOD][5]; // punch-through and background cut for CSI2
  TCutG *cutgpthr2[NMOD]; // punch-through cut for CSI2
  TCutG *cutgpextra[NMOD]; // extra protons in CSI1
  float time0[NMOD][3][4];// time0 parameters [MOD][PD0-2][mean,sigma,tmin,tmax]
  float corr_fct[NMOD][2][5];// energy corrections for p d t 3he 4he
  float Pthr_CsI12[NMOD][2][5];// punch through channels for p d t 3he 4he
  float Pthr_CsI12_sig[NMOD][2][5];// punch through sigma for p d t 3he 4he
  float Pthr_PD01[NMOD][2][5];// punch through channels for p d t 3he 4he
  float Pthr_PD01_sig[NMOD][2][5];// punch through sigma for p d t 3he 4he

  Target *tgt;
  KrattaGeo *kgeo;

//-->   tgt|air|win|ddl0|sil0|ddl1|sil1|myl0|csi1|myl1|csi2|myl2|sil2|ddl2
  float tgt_thickness;
  float tgt_thck[NMOD];
  float air_thck; // um (40 cm of air)
  float win_thck; // module window 100 um of Cu    
  float dl0_thck; // 15 um of Si dead layer for PD0 
  float dl2_thck; // 15 um of Si dead layer for PD2 

  float si0_thck[NMOD]; // um of Si active layer for PD0     
  float dl1_thck[NMOD]; // um of Si dead   layer for PD1     
  float si1_thck[NMOD]; // um of Si active layer for PD1     
  float si2_thck[NMOD]; // um of Si active layer for PD2     

  float myl0_thck[NMOD];  // um of Mylar betw PD1 CSI1      
  float myl1_thck[NMOD];  // um of Mylar betw CSI1 CSI2     
  float myl2_thck [NMOD]; // um of Mylar betw CSI2 PD2      

  float csi1_thck[NMOD]; // um CSI1 thickness 	
  float csi2_thck[NMOD]; // um CSI2 thickness 	

  float fct_si0[NMOD]; //channels/MeV PD0
  float fct_si1[NMOD]; //channels/MeV PD1
  float fct_si2[NMOD]; //channels/MeV PD2

  float a11[NMOD]; //Parlog parameters for CSI1
  float a21[NMOD]; //Parlog parameters for CSI1
  float a31[NMOD]; //Parlog parameters for CSI1
  float a41[NMOD]; //Parlog parameters for CSI1

  float a12[NMOD]; //Parlog parameters for CSI2
  float a22[NMOD]; //Parlog parameters for CSI2
  float a32[NMOD]; //Parlog parameters for CSI2
  float a42[NMOD]; //Parlog parameters for CSI2

  float na11[NMOD][5]; //new parameters for CSI1
  float na21[NMOD][5]; //new parameters for CSI1
  float nof1[NMOD][5]; //new parameters for CSI1
  float na12[NMOD][5]; //new parameters for CSI2
  float na22[NMOD][5]; //new parameters for CSI2
  float nof2[NMOD][5]; //new parameters for CSI2
  RE_TAB *re_tgt; // target material
  RE_TAB *re_air; // air
  RE_TAB *re_cu;  // Cu
  RE_TAB *re_si;  // Si
  RE_TAB *re_myl; // Mylar
  RE_TAB *re_csi; // CsI
  
  
  //static const char *dat_dir;
  //static const char *TABS;
  //static const float fct_vert[NMOD];
  //static const float thr_sct[NMOD]; // separation between ionization/scintillation
                                    // may be module and Z dependent
  // SK change: 
  char dat_dir[200];
  char TABS[200];
  //static const float fct_vert[NMOD];
  //static const float thr_sct[NMOD];
  float fct_vert[NMOD];
  float thr_sct[NMOD];

  float H2E(float fH, int iz, int ia, float a1,float a2,float a3,float a4, float fct_si);
  float H2E_GAMMA(float fH, float a1, float a2,float a3,float a4, float fct_si);
  float E2H(float e_csi, int iz, int ia,float a1,float a2,float a3,float a4, float fct_si);
  float H2E_PTHR(float fH, int iz, int ia, float fa1, float fa2,
	         float fct_si, float csi2_thcki, RE_TAB *re_csi1,float off);
  float H2E_PTHRI(float fH, int iz, int ia, float fa1, float fa2,
	         float fct_si, float csi2_thcki, RE_TAB *re_csi1,float off);

  float SI01_PTHR(float ch_sil0, float ch_sil1,int iz, int ia, 
                  float fct_si0i, float fct_si1i, 
	          float si0_thcki,float dl1_thcki,float si1_thcki,
	          RE_TAB *re_si1);

  float SI0_PTHR(float ch_sil0, int iz, int ia, 
                  float fct_si0i,  
	          float si0_thcki,
	          RE_TAB *re_si1);
  float NH2E( float fH, int iz, int ia, float a1,float a2,float off, float fct_si);
  float NH2EI(float fH, int iz, int ia, float a1,float a2,float off, float fct_si);
  float NE2H( float fE, int iz, int ia, float a1,float a2,float off, float fct_si);
  float NE2HI(float fE, int iz, int ia, float a1,float a2,float off, float fct_si);
  
  ClassDef(KrattaEne,1)  //energy calibration

};

#endif