FSUDAQ_Qt6/analyzers/PID.h

#ifndef PID_h
#define PID_h

/*********************************************
 * This is online analyzer for PID, ANL
 * 
 * Created by Khushi @ 2024-03-27
 * 
 * ******************************************/
#include "Analyser.h"


class PID : public Analyzer{



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


    SetUpdateTimeInSec(1.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);
 
    //========== use the influx from the Analyzer
    influx = new InfluxDB("https://fsunuc.physics.fsu.edu/influx/");
    dataBaseName = "testing"; 

    SetUpCanvas(); // see below

  };

  void SetUpCanvas();

public slots:
  void UpdateHistograms();


private:

  MultiBuilder *evtbder;

  Histogram2D * hPID;
  
  Histogram1D * hdE;    // raw dE (ch=1): ch1
  Histogram1D * hE; 	// raw E (ch=4) : ch4
  Histogram1D * hdT; 	// raw dT (ch=7): ch7
  
  Histogram1D * hTotE; 	// total energy (dE+E): ch1+ch4
  Histogram1D * hTWin; 	// coincidence time window TWin: (t4-t1)*1e9
  
  Histogram2D * hdEE; // dE versus E : ch1 versus ch4
  Histogram2D * hdEtotE; // dE versus totE : ch1 versus (ch1+ch4)
 
  Histogram2D * hdEdT; // dE versus TOF: ch1 versus (t7-t1)*1e9
 
  
  int tick2ns;
  
  float ch1, ch4, ch7;
  unsigned long long t1, t4, t7;
  

};



inline void PID::SetUpCanvas(){

  setGeometry(0, 0, 2000, 1000);  

  //============ histograms
  
  hdEE = new Histogram2D("dE vs E", "E[ch]", "dE[ch]", 100, 0, 7000, 100, 0, 4000, this);
  layout->addWidget(hdEE, 0, 0, 1, 2); 
  
  hdE = new Histogram1D("raw dE (ch=1)", "dE [ch]", 300, 0, 2500, this);
  layout->addWidget(hdE, 0, 2);
  
  hdEdT = new Histogram2D("dE vs TOF", "TOF [ns]", "dE", 100, 0, 500, 100, 0, 4000, this);
  layout->addWidget(hdEdT, 0, 3);
  
  hE = new Histogram1D("raw E (ch=4)", "E [ch]", 300, 0, 5000, this);
  layout->addWidget(hE, 0, 4);
  
  hdEtotE = new Histogram2D("dE vs TotE", "TotE[ch]", "dE[ch]", 100, 0, 8000, 100, 0, 4000, this);
  layout->addWidget(hdEtotE, 1, 0, 1, 2);
  
  hdT = new Histogram1D("raw dT (ch=7)", "dT [ch]", 300, 0, 1000, this);
  layout->addWidget(hdT, 1, 2);
  
  hTotE = new Histogram1D("total energy (dE+E)", "TotE [ch]", 300, 0, 7000, this);
  layout->addWidget(hTotE, 1, 3);
  
  hTWin = new Histogram1D("coincidence time window", "TWin [ns]", 300, 0, 500, this);
  layout->addWidget(hTWin, 1, 4);
  
}

inline void PID::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};

  //============ 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;
   
    for( int k = 0; k < (int) event.size(); k++ ){
      //event[k].Print();
      if( event[k].ch == 0 ) {ch1 = event[k].energy; t1 = event[k].timestamp;}
      if( event[k].ch == 1 ) {ch4 = event[k].energy; t4 = event[k].timestamp;}
      if( event[k].ch == 2 ) {ch7 = event[k].energy; t7 = event[k].timestamp;}
     
    }

    // printf("(E, dE) = (%f, %f)\n", E, dE);
    hPID->Fill(ch4 + RandomGauss(0, 100), ch1 + RandomGauss(0, 100)); // x, y
    hdE->Fill(ch1);    
	hE->Fill(ch4);
	hdT->Fill(ch7);
	hTotE->Fill(ch1+ch4);
	hdEE->Fill(ch4,ch1);
	hdEtotE->Fill(ch1+ch4+ RandomGauss(0, 100),ch1+ RandomGauss(0, 100));
	hdEdT->Fill((t7-t1)*1e9,ch1);
	hTWin->Fill((t4-t1)*1e9);
    
    //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(ch4, ch1), Qt::OddEvenFill) ){
        if( t1 < tMin[p] ) tMin[p] = t1;
        if( t1 > tMax[p] ) tMax[p] = t1;
        count[p] ++;
        //printf(".... %d \n", count[p]);
      }
    }
    
  }

  hPID->UpdatePlot();
  hdE->UpdatePlot();    
  hE->UpdatePlot();
  hdT->UpdatePlot();
  hTotE->UpdatePlot();
  hdEE->UpdatePlot();
  hdEtotE->UpdatePlot();
  hdEdT->UpdatePlot();
  hTWin->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