#ifndef MCPandPSD_h
#define MCPandPSD_h

/*********************************************
 * This is online analyzer for RASIOR, ANL
 * 
 * Created by Ryan @ 2023-10-16
 * 
 * ******************************************/
#include "Analyser.h"
#include <cmath>
#include "math.h"
#include <algorithm>


class MCPandPSD : public Analyzer{

public:
  MCPandPSD(Digitizer ** digi, unsigned int nDigi, QMainWindow * parent = nullptr): Analyzer(digi, nDigi, parent){


    SetUpdateTimeInSec(4.0);

    RedefineEventBuilder({0}); // only builder for the 0-th digitizer.
    tick2ns = digi[0]->GetTick2ns();

    SetBackwardBuild(false, 100); // using normal building (acceding in time) or backward building, int the case of backward building, default events to be build is 100. 
    evtbder = GetEventBuilder();
    evtbder->SetTimeWindow(500);
 
    SetDatabase("https://localhost:8086", "testing", "zKhzKk4Yhf1l9QU-yE2GsIZ1RazqUgoW3NlF8LJqq_xDMwatOJwg1sKrjgq36uLEsQf8Fmn4sJALP7Kkilk14A==");

    SetUpCanvas(); // see below

  };

  void SetUpCanvas();

public slots:
  void UpdateHistograms();


private:

  MultiBuilder *evtbder;

  //Histogram2D * hPID;
  
  Histogram2D * hXYE; // 2D energy plot: e[2]+e[3] versus e[0]+e[1]
  
  Histogram1D * hX; // X position:((e[0]-e[1])/(e[0]+e[1]))
  Histogram1D * hY; // Y position:((e[2]-e[3])/(e[2]+e[3]))
  
  Histogram1D * hXmcp; // X position
  Histogram1D * hYmcp; // Y position
  
  Histogram2D * hXY; // 2D position plot: ((e[2]-e[3])/(e[2]+e[3])) versus ((e[0]-e[1])/(e[0]+e[1]))
  Histogram2D * hXYMCP; // 2D position plot for MCP: ((e[2]+e[3])/((e[0]+e[1]+e[2]+e[3]))) versus ((e[0]+e[1])/(e[0]+e[1]+e[2]+e[3]))
  Histogram2D * hXYr; // 2D position plot rotated for MCP:
  
  
  Histogram2D * hXEdE1; //X energy versus dE signal 1
  Histogram2D * hYEdE1; //Y energy versus dE signal 1
  
  Histogram2D * hXEdE2; //X energy versus dE signal 2
  Histogram2D * hYEdE2; //Y energy versus dE signal 2
  /*
  Histogram1D * he0; // e0: signal 0 from PSD
  Histogram1D * he1; // e1: signal 1 from PSD
  Histogram1D * he2; // e2: signal 2 from PSD
  Histogram1D * he3; // e3: signal 3 from PSD
  
  Histogram1D * hmcp0; // s0: signal 0 from MCP
  Histogram1D * hmcp1; // s1: signal 1 from MCP
  Histogram1D * hmcp2; // s2: signal 2 from MCP
  Histogram1D * hmcp3; // s3: signal 3 from MCP
  */
  
  
  int tick2ns;
  
  //float dE, E;
  //unsigned long long dE_t, E_t;
  
  float e0, e1, e2, e3, dE1, dE2;
  unsigned long long t0, t1, t2, t3, dE1_t, dE2_t;
  
  float s0, s1, s2, s3;
  unsigned long long s_t0, s_t1, s_t2, s_t3;
  
};


inline void MCPandPSD::SetUpCanvas(){

  setGeometry(0, 0, 1500, 2000);  

  //============ histograms
  //hPID = new Histogram2D("RAISOR2", "E", "dE", 100, 0, 11000, 100, 0, 11000, this);
  //layout->addWidget(hPID, 0, 0); 
  
  hXY = new Histogram2D("2D position plot PSD_E", "X position", "Y position", 200, -1, 1, 200, -1, 1, this);
  layout->addWidget(hXY, 0, 0);
  
  hX = new Histogram1D("X position", "X", 300, -1, 1, this);
  layout->addWidget(hX, 0, 1);
  
  hY = new Histogram1D("Y position", "Y", 300, -1, 1, this);
  layout->addWidget(hY, 0, 2);
  /*
  he0 = new Histogram1D("PSD_E 0", "e0", 200, 0, 8000, this);
  layout->addWidget(he0, 0, 1);
  
  he1 = new Histogram1D("PSD_E 1", "e1", 200, 0, 8000, this);
  layout->addWidget(he1, 0, 2);
  
  he2 = new Histogram1D("PSD_E 2", "e2", 200, 0, 8000, this);
  layout->addWidget(he2, 0, 3);
 
  he3 = new Histogram1D("PSD_E 3", "e3", 200, 0, 8000, this);
  layout->addWidget(he3, 0, 4);
  */
  hXYMCP = new Histogram2D("2D position MCP", "X position", "Y position", 500, 0, 1, 500, 0,  1, this);
  layout->addWidget(hXYMCP, 1, 1);
  
  hXYr = new Histogram2D("2D rot pos MCP", "Xr position", "Yr position", 200, -0.5, 0.5, 200, -0.5, 0.5, this);
  layout->addWidget(hXYr, 1, 0);
  /*
  hmcp0 = new Histogram1D("MCP 0", "s0", 200, 0, 8000, this);
  layout->addWidget(hmcp0, 1, 1);
  
  hmcp1 = new Histogram1D("MCP 1", "s1", 200, 0, 8000, this);
  layout->addWidget(hmcp1, 1, 2);
  
  hmcp2 = new Histogram1D("MCP 2", "s2", 200, 0, 8000, this);
  layout->addWidget(hmcp2, 1, 3);
 
  hmcp3 = new Histogram1D("MCP 3", "s3", 200, 0, 8000, this);
  layout->addWidget(hmcp3, 1, 4);
  */
  hXmcp = new Histogram1D("X pos rot MCP", "X", 250, -0.5, 0.5, this);
  layout->addWidget(hXmcp, 1, 2);
  
  hYmcp = new Histogram1D("Y pos rot MCP", "Y", 250, -0.5, 0.5, this);
  layout->addWidget(hYmcp, 1, 3);
  
  hXEdE1 = new Histogram2D("X energy versus dE signal 1", "Ex", "dE signal 1", 100, 0, 8000, 100, 0, 8000, this);
  layout->addWidget(hXEdE1, 2, 0);
  
  hYEdE1 = new Histogram2D("Y energy versus dE signal 1", "Ey", "dE signal 1", 100, 0, 8000, 100, 0, 8000, this);
  layout->addWidget(hYEdE1, 2, 1);
  
  hXEdE2 = new Histogram2D("X energy versus dE signal 2", "Ex", "dE signal 2", 100, 0, 8000, 100, 0, 8000, this);
  layout->addWidget(hXEdE2, 2, 2);
  
  hYEdE2 = new Histogram2D("Y energy versus dE signal 2", "Ey", "dE signal 2", 100, 0, 8000, 100, 0, 8000, this);
  layout->addWidget(hYEdE2, 2, 3);
  
  hXYE = new Histogram2D("2D energy plot", "Ex", "Ey", 100, 0, 8000, 100, 0, 8000, this);
  layout->addWidget(hXYE, 0, 3);
  


}

inline void MCPandPSD::UpdateHistograms(){

  if( this->isVisible() == false ) return;
  
  BuildEvents(false); // call the event builder to build events

  //============ Get events, and do analysis
  long eventBuilt = evtbder->eventBuilt;
  if( eventBuilt == 0 ) return;

  //============ Get the cut list, if any
  /*
  QList<QPolygonF> cutList = hPID->GetCutList();
  const int nCut = cutList.count();
  unsigned long long tMin[nCut] = {0xFFFFFFFFFFFFFFFF}, tMax[nCut] = {0};
  unsigned int count[nCut]={0};
  
  QList<QPolygonF> cutList1 = hXEdE1->GetCutList();
  const int nCut1 = cutList1.count();
  unsigned long long tMin1[nCut1] = {0xFFFFFFFFFFFFFFFF}, tMax1[nCut1] = {0};
  unsigned int count1[nCut1]={0};
  
  QList<QPolygonF> cutList2 = hYEdE1->GetCutList();
  const int nCut2 = cutList2.count();
  unsigned long long tMin2[nCut2] = {0xFFFFFFFFFFFFFFFF}, tMax2[nCut2] = {0};
  unsigned int count2[nCut2]={0};
  
  QList<QPolygonF> cutList3 = hXY->GetCutList();
  const int nCut3 = cutList3.count();
  unsigned long long tMin3[nCut3] = {0xFFFFFFFFFFFFFFFF}, tMax3[nCut3] = {0};
  unsigned int count3[nCut3]={0};
  */ 
  //============ Processing data and fill histograms
  long eventIndex = evtbder->eventIndex;
  long eventStart = eventIndex - eventBuilt + 1;
  if(eventStart < 0 ) eventStart += MaxNEvent;
  
  for( long i = eventStart ; i <= eventIndex; i ++ ){
    std::vector<Hit> event = evtbder->events[i];
    //printf("-------------- %ld\n", i);

    if( event.size() == 0 ) return;
    
    
    if( event.size() == 0 ) return;
    //if( event.size() < 2 ) return;
    cout<< "event size " << event.size() <<endl;
    
    s0 = 0;
	s1 = 0;
	s2 = 0;
	s3 = 0;
    
    s_t0 = 0;
	s_t1 = 0;
	s_t2 = 0;
	s_t3 = 0;
    
    for( int k = 0; k < (int) event.size(); k++ ){
      //event[k].Print();
      
      if( event[k].ch == 2 ) {s0 = event[k].energy; s_t0 = event[k].timestamp;}	//
      if( event[k].ch == 3 ) {s1= event[k].energy;  s_t1 = event[k].timestamp;} //    The 4 output signals from the 
      if( event[k].ch == 4 ) {s2 = event[k].energy; s_t2 = event[k].timestamp;} //			MCP detector
      if( event[k].ch == 5 ) {s3= event[k].energy;  s_t3 = event[k].timestamp;} //

      if( event[k].ch == 10 ) {e0 = event[k].energy; t0 = event[k].timestamp;} 		//
      if( event[k].ch == 11 ) {e1= event[k].energy;  t1 = event[k].timestamp;} 		//  The 4 output signals from the
      if( event[k].ch == 12 ) {e2 = event[k].energy; t2 = event[k].timestamp;} 		//  position sensitive E detector
      if( event[k].ch == 13 ) {e3= event[k].energy;  t3 = event[k].timestamp;} 		//
      
      if( event[k].ch == 14 ) {dE1 = event[k].energy; dE1_t = event[k].timestamp;} 	//  The 2 output signals from the
      if( event[k].ch == 15 ) {dE2= event[k].energy;  dE2_t = event[k].timestamp;} 	//	    square dE detector
    }

	if (s0>10 && s1>10 && s2>10 && s3>10) {
	float_t rotation_angle = 31.;
	double_t Xr = (((s1+s2)/(s0+s1+s2+s3))-0.51)*cos(-rotation_angle*M_PI/180)-(((s2+s3)/(s0+s1+s2+s3))-0.51)*sin(-rotation_angle*M_PI/180);
    double_t Yr = (((s1+s2)/(s0+s1+s2+s3))-0.51)*sin(-rotation_angle*M_PI/180)+(((s2+s3)/(s0+s1+s2+s3))-0.51)*cos(-rotation_angle*M_PI/180);
    // printf("(E, dE) = (%f, %f)\n", E, dE);
    //hPID->Fill(E + RandomGauss(0, 100), dE+ RandomGauss(0, 100)); // x, y
    
	hXY->Fill(((e0-e1)/(e0+e1)),((e3-e2)/(e2+e3)));
	hXYMCP->Fill(((s1+s2)/(s0+s1+s2+s3)),((s2+s3)/(s0+s1+s2+s3)));
	
	hX->Fill(((e0-e1)/(e0+e1)));
	hY->Fill(((e3-e2)/(e2+e3)));
	hXmcp->Fill(Xr);
	hYmcp->Fill(Yr);
	hXEdE1->Fill((e0+e1),dE1);
	hYEdE1->Fill(e2+e3,dE1);
	hXEdE2->Fill(e0+e1,dE2);
	hYEdE2->Fill(e2+e3,dE2);
	hXYE->Fill(e0+e1,e2+e3);
	hXYr->Fill(Xr,Yr);
	}
	/*
	he0->Fill(e0);
	he1->Fill(e1);
	he2->Fill(e2);
	he3->Fill(e3);
	
	hmcp0->Fill(s0);
	hmcp1->Fill(s1);
	hmcp2->Fill(s2);
	hmcp3->Fill(s3);
	*/ 
    //check events inside any Graphical cut and extract the rate
    /*
    for(int p = 0; p < cutList.count(); p++ ){ 
      if( cutList[p].isEmpty() ) continue;
      if( cutList[p].containsPoint(QPointF(E, dE), Qt::OddEvenFill) ){
        if( dE_t < tMin[p] ) tMin[p] = dE_t;
        if( dE_t > tMax[p] ) tMax[p] = dE_t;
        count[p] ++;
        //printf(".... %d \n", count[p]);
      }
    }
    
    for(int p = 0; p < cutList1.count(); p++ ){ 
      if( cutList1[p].isEmpty() ) continue;
      if( cutList1[p].containsPoint(QPointF((e0+e1), dE1), Qt::OddEvenFill) ){
        if( dE1_t < tMin1[p] ) tMin1[p] = dE1_t;
        if( dE1_t > tMax1[p] ) tMax1[p] = dE1_t;
        count1[p] ++;
        //printf("hXX.... %d \n", count1[p]);
      }
    }
    
    for(int p = 0; p < cutList2.count(); p++ ){ 
      if( cutList2[p].isEmpty() ) continue;
      if( cutList2[p].containsPoint(QPointF((e2+e3), dE1), Qt::OddEvenFill) ){
        if( dE1_t < tMin2[p] ) tMin2[p] = dE1_t;
        if( dE1_t > tMax2[p] ) tMax2[p] = dE1_t;
        count2[p] ++;
        //printf("hXX.... %d \n", count2[p]);
      }
    }
    
    for(int p = 0; p < cutList3.count(); p++ ){ 
      if( cutList3[p].isEmpty() ) continue;
      if( cutList3[p].containsPoint(QPointF(((e0-e1)/(e0+e1)), ((e2-e3)/(e2+e3))), Qt::OddEvenFill) ){
        if( ((t2-t3)/(t2+t3)) < tMin3[p] ) tMin3[p] = ((t2-t3)/(t2+t3));
        if( ((t2-t3)/(t2+t3)) > tMax3[p] ) tMax3[p] = ((t2-t3)/(t2+t3));
        count3[p] ++;
        //printf("hXX.... %d \n", count3[p]);
      }
    }
   */ 
  }

  //hPID->UpdatePlot();
  hXY->UpdatePlot();
  hXYr->UpdatePlot();
  hXYMCP->UpdatePlot();
  
  hX->UpdatePlot();
  hY->UpdatePlot();
  hXmcp->UpdatePlot();
  hYmcp->UpdatePlot();
  hXEdE1->UpdatePlot();
  hYEdE1->UpdatePlot();
  hXEdE2->UpdatePlot();
  hYEdE2->UpdatePlot();
  hXYE->UpdatePlot();
  /*
  he0->UpdatePlot();
  he1->UpdatePlot();
  he2->UpdatePlot();
  he3->UpdatePlot();
  
  hmcp0->UpdatePlot();
  hmcp1->UpdatePlot();
  hmcp2->UpdatePlot();
  hmcp3->UpdatePlot();
  */
  //========== output to Influx
  /*
  QList<QString> cutNameList = hPID->GetCutNameList();
  for( int p = 0; p < cutList.count(); p ++){
    if( cutList[p].isEmpty() ) continue;
    double dT = (tMax[p]-tMin[p]) * tick2ns / 1e9; // tick to sec
    double rate = count[p]*1.0/(dT);
    //printf("%llu %llu, %f %d\n", tMin[p], tMax[p], dT, count[p]);
    //printf("%10s | %d | %f Hz \n", cutNameList[p].toStdString().c_str(), count[p], rate);  
    
    influx->AddDataPoint("Cut,name=" + cutNameList[p].toStdString()+ " value=" + std::to_string(rate));
    influx->WriteData(dataBaseName);
    influx->ClearDataPointsBuffer();
  }
  */ 
}


#endif