// ---------------------------------------------------------------------
// -----                  Swiatowid source file                    -----
// -----            Created 10/11/13  by P. Pawlowski              -----
// ---------------------------------------------------------------------


#include "Swiatowid.hxx"

namespace Kratta{

  
  enum pattern_mask{
    EMPTY   = 0b00000000,	// No signal
    PD0     = 0b00000001,	// PD0 fired
    PD1     = 0b00000010,	// PD1 fired
    PD2     = 0b00000100,	// PD2 fired
    SISI    = 0b00001000,	// SI0 vs SI1
    PUNCH   = 0b00010000,	// CsI2 punched through
    GAMMA   = 0b00100000,	// Gamma particle
    NEUTRON = 0b01000000,	// Neutron
    BCKGRD  = 0b10000000,	// background
 };

  bool pd0time(const ASYFadcPeak * peak){
    return peak->p04 >0;//  (peak->p04 > 180 && peak->p04 < 200);
  }

  bool pd1time(const ASYFadcPeak * peak){
    return peak->p14 > 0;// (peak->p14 > 190 && peak->p14 < 230);
  }

  bool pd2time(const ASYFadcPeak * peak){
    return peak->p24 > 0;// (peak->p24 > 190 && peak->p24 < 230);
  }

  bool punchthrough(const ASYFadcPeak * peak){
    ///TODO: This function should indicate the punch-through particles
    return false; //peak->p20>0;
  }

  bool gamma(const ASYFadcPeak * peak){
    ///TODO: This function should indicate the gamma particles
    return false;
  }

  bool neutron(const ASYFadcPeak * peak){
    ///TODO: This function should indicate the neutrons
    return false;
  }

  bool sisi(const ASYFadcPeak * peak){
    //return peak->p12+peak->p13<1.2*peak->p10/fct_vert[peak->mod];
   return peak->p12+peak->p13-peak->p10/fct_vert[peak->mod]<thr_sct[peak->mod];
  }

  bool background(const ASYFadcPeak * peak){
    ///TODO: This function should indicate the background particles
    return false;
  }

 /// PD0_PD0
 float pd0pd0X(const ASYFadcPeak *peak){
//    return ((peak->p09!=peak->p05)*peak->p05*peak->p09/(peak->p05-peak->p09)*log(peak->p05/peak->p09)+(peak->p09==peak->p05)*peak->p09)*1000;
    return  1000*(
	      peak->p09 == peak->p05 		?  peak->p09:
	      peak->p09 == 0 || peak->p05 == 0 	? 0:
	      peak->p05*peak->p09/(peak->p05-peak->p09)*log(peak->p05/peak->p09)
	    );
  }

float pd0pd0Y(const ASYFadcPeak *peak){
    return peak->p00;
  }
bool pd0pd0C(const ASYFadcPeak *peak){
    return peak->am1 == 0 && peak->p00 > 120;

  }

/// PD0_PD1
float pd0pd1X(const ASYFadcPeak *peak){
    return peak->p10 + peak->p12 + peak->p13;
  }
float pd0pd1Y(const ASYFadcPeak *peak){
    return peak->p00;
  }
bool pd0pd1C(const ASYFadcPeak *peak){
    return pd0time(peak) && pd1time(peak) && peak->am2==0 && peak->p00>30 && peak->am1>30
    && sisi(peak);
  }

/// PD1_PD1
float pd1pd1X(const ASYFadcPeak *peak){
    return peak->p12 + peak->p13 - peak->p10/fct_vert[peak->mod];
  }
float pd1pd1Y(const ASYFadcPeak *peak){
    return peak->p10 + peak->p10/fct_vert[peak->mod];
  }
bool  pd1pd1C(const ASYFadcPeak *peak){
    return pd1time(peak) && peak->am2==0;
  }

/// PD01_PD1
float pd01pd1X(const ASYFadcPeak *peak){
    return peak->p12 + peak->p13 - peak->p10/fct_vert[peak->mod];
  }
float pd01pd1Y(const ASYFadcPeak *peak){
    return peak->p00 +peak->p10 + peak->p10/fct_vert[peak->mod];
  }
bool  pd01pd1C(const ASYFadcPeak *peak){
    return pd1time(peak) && peak->am2==0;
  }

/// PD1_PD2
float pd1pd2X(const ASYFadcPeak *peak){
    return peak->p22 + peak->p23;
  }
float pd1pd2Y(const ASYFadcPeak *peak){
    return peak->p12 + peak->p13 - peak->p10/fct_vert[peak->mod];
  }
bool  pd1pd2C(const ASYFadcPeak *peak){
    return pd1time(peak) && pd2time(peak) &&  peak->p12 + peak->p13 - peak->p10/fct_vert[peak->mod]>60 && peak->p22 + peak->p23>30;
  }

/// PD2_PD2
float pd2pd2X(const ASYFadcPeak *peak){
    return peak->p20*20;
  }
float pd2pd2Y(const ASYFadcPeak *peak){
    return peak->p22 + peak->p23;
  }
bool  pd2pd2C(const ASYFadcPeak *peak){
    return pd2time(peak) &&  punchthrough(peak);
  }

/// FS_PD1
float fspd1X(const ASYFadcPeak *peak){
    return peak->p13;
  }
float fspd1Y(const ASYFadcPeak *peak){
    return peak->p12;
  }
bool  fspd1C(const ASYFadcPeak *peak){
    return peak->am2 == 0;
  }

/// FS_PD2
float fspd2X(const ASYFadcPeak *peak){
    return peak->p23;
  }
float fspd2Y(const ASYFadcPeak *peak){
    return peak->p22;
  }
bool  fspd2C(const ASYFadcPeak *peak){
    return true;
  }

  void RealACorr(TRootKRATParticle * part){
    static const int PARTID  = SISI | PD0 | PD1 | PD2;
    switch(part->Code & PARTID){
      case PD0 | PD1 | SISI :
	switch(part->Z){
	  case 3: part->RealA +=  0.06; break;
	  case 4: part->RealA += -0.07; break;
	  case 5: part->RealA += -0.20; break;
	  case 6: part->RealA += -0.30; break;
	  default : break;
	}; break;
      default: return;
    }
    part->A = int(part->RealA +0.5);
  }

};


///Default C-tor
Swiatowid::Swiatowid():
  Run(0)
  //slow_control(SLOW_CONTROL.c_str())
{
  /// TODO: Set in function - as a parameter of them:
  ///TString  sCalibrationDirLocalDefPath  = "calfiles/KRATTA/IdentGrids/KrattaCalib";     
  TString  sCalibrationDirLocalDefPath  = "calfiles/KRATTA/Particle_identyfication/Identyfication_grids";    

  TString sAsyeosrootMainDir = gSystem->ExpandPathName("$VMCWORKDIR");

  TString  sFullPathToKrattaCalibrDir = sAsyeosrootMainDir + "/" + sCalibrationDirLocalDefPath;
  ///SLOW_CONTROL       = KrattaCalibrationDirectory + "/slow_control.bin";  ///SLOW_CONTROL - this variable is not used. 
  ENERGY_CALIBRATION = sFullPathToKrattaCalibrDir + "/kratta_calib_%02d_r.par";

  MASK[PD0_PD0]  = sFullPathToKrattaCalibrDir + "/grid_PD0_PD0_%02d.root";         //PD0_PD0
  MASK[PD0_PD1]  = sFullPathToKrattaCalibrDir + "/grid_PD0_PD1_%02d.root";         //PD0_PD1
  MASK[PD01_PD1] = sFullPathToKrattaCalibrDir + "/grid_PD01_PD1_%02d.root";        //PD01_PD1
  MASK[PD1_PD1]  = sFullPathToKrattaCalibrDir + "/grid_PD1_PD1_%02d.root";         //PD1_PD1
  MASK[PD1_PD2]  = sFullPathToKrattaCalibrDir + "/grid_PD1_PD2_%02d.root";         //PD1_PD2
  //MASK[PD2_PD2]  = sFullPathToKrattaCalibrDir + "/grid_PD2_PD2_%02d.root";         //PD2_PD2
  //MASK[FS_PD1]   = sFullPathToKrattaCalibrDir + "/grid_FS_PD1_%02d.root";          //FS_PD1
  //MASK[FS_PD2]   = sFullPathToKrattaCalibrDir + "/grid_FS_PD2_%02d.root";          //FS_PD2,

  //Calibrator.Init(ENERGY_CALIBRATION);                                  
  
  for(int i=0; i<NSPECT; i++)
    Recognizer[i] = new KrattaRecognizer(MASK[i],"c1");

  Recognizer[0]->SetTolerance(0.1);

  for(int i=0; i<NMODULE; i++){
    Spectrum[i][PD0_PD0]  = new KrattaSpectrum("PD0_PD0",i, "p00:((p09!=p05)*p05*p09/(p05-p09)*log(p05/p09)+(p09==p05)*p09)*1000","am1==0 && p00>120", &pd0pd0X, &pd0pd0Y,&pd0pd0C);
    Spectrum[i][PD0_PD1]  = new KrattaSpectrum("PD0_PD1",i, "p00:p10+p12+p13", Form("am2==0 && p00>30 && am1>30 && p12+p13<1.2*p10/%.2f",fct_vert[i]), &pd0pd1X, &pd0pd1Y, &pd0pd1C);
    Spectrum[i][PD01_PD1] = new KrattaSpectrum("PD01_PD1",i, Form("p00+p10+p10/%.2f:p12+p13-p10/%.2f", fct_vert[i], fct_vert[i]), "am2==0", &pd01pd1X, &pd01pd1Y, &pd01pd1C);
    Spectrum[i][PD1_PD1]  = new KrattaSpectrum("PD1_PD1",i, Form("p10+p10/%.2f:p12+p13-p10/%.2f", fct_vert[i], fct_vert[i]), "am2==0", &pd1pd1X, &pd1pd1Y, &pd1pd1C);
    Spectrum[i][PD1_PD2]  = new KrattaSpectrum("PD1_PD2",i, Form("p12+p13-p10/%.2f:p22+p23",fct_vert[i]), Form("p12+p13-p10/%.2f>60&&p22+p23>30", fct_vert[i]), &pd1pd2X, &pd1pd2Y, &pd1pd2C);
    //Spectrum[i][PD2_PD2] = new KrattaSpectrum("PD2_PD2",i, "p22+p23:p20*20","p20>0", &pd2pd2X, &pd2pd2Y, &pd2pd2C);
    //Spectrum[i][FS_PD1] = new KrattaSpectrum("FS_PD1",i, "p12:p13", "am2==0", &fspd1X, &fspd1Y, &fspd1C); // Not used (2013-11-21)
    //Spectrum[i][FS_PD2] = new KrattaSpectrum("FS_PD2",i, "p22:p23", "", &fspd2X, &fspd2Y, &fspd2C); // Not used (2013-11-21)
  }

  pattern.join_mask(PD0,     &pd0time);
  pattern.join_mask(PD1,     &pd1time);
  pattern.join_mask(PD2,     &pd2time);
  pattern.join_mask(SISI,    &sisi);
  pattern.join_mask(PUNCH,   &punchthrough);
  pattern.join_mask(GAMMA,   &gamma);
  pattern.join_mask(NEUTRON, &neutron);
  pattern.join_mask(BCKGRD,  &background);
}

  ///Default D-tor
  Swiatowid::~Swiatowid(){}

  /// Analyzing a peak
  void Swiatowid::Analyze(ASYFadcPeak * peak, TRootKRATParticle * part){
    
    //part->Clear(); // It is done in Task TKratCalib
    part->Module = peak->mod;
    
    SlowControl(peak);
    SetCode(peak, part);
    Identify(peak, part);        /// <= MEMORY LEAK (>>)
    //Calibrate(peak, part);
  }

/// Getting particle code
void Swiatowid::SetCode(ASYFadcPeak * peak, TRootKRATParticle * part){
  part->Code = pattern.set(peak);
}

/// Identifying a particle
void Swiatowid::Identify(ASYFadcPeak * peak, TRootKRATParticle * part)const{
  
    if(part->Code & GAMMA){
	     part->Z = 0;
	     part->A = 0;
	     return;
    }
    if(part->Code & NEUTRON){
       part->Z = 0;
       part->A = 1;
       return;
    }

    static const int PARTID  = PD0 | PD1 | PD2;
    // switch(part->Code & PARTID){
    // 	case  0  |  0  |  0  :	break;
    // 	case PD0 |  0  |  0  : 	Recognize(PD0_PD0,peak,part); break;
    // 	case  0  | PD1 |  0  :	Recognize(PD1_PD1,peak,part); break;
    // 	case PD0 | PD1 |  0  : 	if(part->Code & SISI)
    // 				                    Recognize(PD0_PD1,peak,part);
    // 				                  else
    // 				                    Recognize(PD01_PD1,peak,part);
    // 				                  //Recognize(PD1_PD1,peak,part);
    // 				                  break;
    // 	case  0  |  0  | PD2 :
    // 	case PD0 |  0  | PD2 : 	break;
    // 	case  0  | PD1 | PD2 :
    // 	case PD0 | PD1 | PD2 :  Recognize(PD1_PD2,peak,part); break;
    // };
    switch(part->Code & PARTID){
      case  0  |  0  |  0  :	return;
      case PD0 |  0  |  0  : 	part->Spect = PD0_PD0; break;
      case  0  | PD1 |  0  :	part->Spect = PD1_PD1; break;
      case PD0 | PD1 |  0  : 	if(part->Code & SISI)
                                part->Spect = PD0_PD1;
                              else
                                part->Spect = PD01_PD1;
                              //Recognize(PD1_PD1,peak,part);
                              break;
      case  0  |  0  | PD2 :
      case PD0 |  0  | PD2 : 	return;
      case  0  | PD1 | PD2 :
      case PD0 | PD1 | PD2 :  part->Spect = PD1_PD2; break;
    };
    Recognize(peak, part);     /// <= MEMORY LEAK (>>>)
    RealACorr(part);
  } 

/// Recognizing a particle
// void Swiatowid::Recognize(spectrum spect, ASYFadcPeak * peak, TRootKRATParticle * part, KrattaGrid * grid)const{
//         float x = Spectrum[peak->mod][spect]->GetX(peak);
//         float y = Spectrum[peak->mod][spect]->GetY(peak);
// 	part->Xident = x;
// 	part->Yident = y;
// 	part->Spectrum = spect;
// 	if(!grid) Recognizer[spect]->Recognize(part,x,y);
// 	else KrattaRecognizer::SetUp(part, x, y, grid);
// }

/// Recognizing a particle
void Swiatowid::Recognize(ASYFadcPeak * peak, TRootKRATParticle * part, KrattaGrid * grid)const{
  
	part->Xident = Spectrum[part->Module][part->Spect]->GetX(peak);
	part->Yident = Spectrum[part->Module][part->Spect]->GetY(peak);
	//std::cout << "[Swiatowid::Recognize:] KrattaGrid * grid = " <<  * grid << std::endl;
  //printf("[Swiatowid::Recognize:] KrattaGrid *grid=%d \n", grid);
  
  if(!grid)
  {
    Recognizer[part->Spect]->Recognize(part);         /// <= MEMORY LEAK (>>>>)
  }
	else
  {
    KrattaRecognizer::SetUp(part, grid);
  }
 }

/// Calibrating energy
//void Swiatowid::Calibrate(ASYFadcPeak * peak, TRootKRATParticle * part)const{
//  ///TODO:NEW INTERFACE:Calibrator.Calibrate(run, part, peak);
//  Calibrator.Calibrate(part, peak);
//}

/// Slow control
void Swiatowid::SlowControl(ASYFadcPeak * peak)const{}