FSUDAQ_Qt6/analyzers/MCP.h
2024-06-04 11:59:18 -05:00

261 lines
7.3 KiB
C++

#ifndef MCP_h
#define MCP_h
/*********************************************
* This is online analyzer for MCP, ANL
*
* Created by Khushi @ 2024-03-27
*
* ******************************************/
#include "Analyser.h"
#include <cmath>
#include "math.h"
#include <algorithm>
//#define M_PI 3.14159265
class MCP : public Analyzer{
public:
MCP(Digitizer ** digi, unsigned int nDigi, QMainWindow * parent = nullptr): Analyzer(digi, nDigi, parent){
SetUpdateTimeInSec(2.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); //ns
//========== 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;
//Histogram2D * hXX; // X1 versus X2 : e[1] versus e[0]
//Histogram2D * hYY; // Y1 versus Y2 : e[3] versus e[2]
Histogram1D * hX; // X position:((e[0]+e[1])/(e[0]+e[1]+e[2]+e[3]))
Histogram1D * hY; // Y position:((e[2]+e[3])/((e[0]+e[1]+e[2]+e[3])))
Histogram1D * hXr; // X position angle rotated
Histogram1D * hYr; // Y position angle rotated
Histogram2D * hXY; // 2D position plot: ((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;
Histogram1D * he0; // e0: signal 0
Histogram1D * he1; // e1: signal 1
Histogram1D * he2; // e2: signal 2
Histogram1D * he3; // e3: signal 3
Histogram1D * ht; // time window
int tick2ns;
float dE, E;
unsigned long long dE_t, E_t;
float e0, e1, e2, e3;
unsigned long long t0, t1, t2, t3;
};
inline void MCP::SetUpCanvas(){
setGeometry(0, 0, 2000, 800);
//============ histograms
hPID = new Histogram2D("MCP", "E", "dE", 100, 0, 5000, 100, 0, 5000, this);
layout->addWidget(hPID, 0, 0);
hXYr = new Histogram2D("2D rot posi", "Xr position", "Yr position", 200, -0.5, 0.5, 200, -0.5, 0.5, this);
layout->addWidget(hXYr, 0, 1);
hXY = new Histogram2D("2D position plot", "X position", "Y position", 500, 0, 1, 500, 0, 1, this);
layout->addWidget(hXY, 0, 2);
ht = new Histogram1D("Time Window", "t", 50, 0, 500, this);
layout->addWidget(ht, 0, 3);
hX = new Histogram1D("X position", "X", 250, 0, 1, this);
layout->addWidget(hX, 2, 0);
hY = new Histogram1D("Y position", "Y", 250, 0, 1, this);
layout->addWidget(hY, 2, 1);
hXr = new Histogram1D("Angle rot X posi", "Xr", 250, -0.5, 0.5, this);
layout->addWidget(hXr, 2, 2);
hYr = new Histogram1D("Angle rot Y posi", "Yr", 250, -0.5, 0.5, this);
layout->addWidget(hYr, 2, 3);
// UNCOMMENT FOLLOWING 8 LINES TO SEE INDIVIDUAL SIGNALS
he0 = new Histogram1D("Signal 0", "e0", 200, 0, 8000, this);
layout->addWidget(he0, 1, 0);
he1 = new Histogram1D("Signal 1", "e1", 200, 0, 8000, this);
layout->addWidget(he1, 1, 1);
he2 = new Histogram1D("Signal 2", "e2", 200, 0, 8000, this);
layout->addWidget(he2, 1, 2);
he3 = new Histogram1D("Signal 3", "e3", 200, 0, 8000, this);
layout->addWidget(he3, 1, 3);
//
}
inline void MCP::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;
//if( event.size() < 2 ) return;
cout<< "event size " << event.size() <<endl;
e0 = 0;
e1 = 0;
e2 = 0;
e3 = 0;
t0 = 0;
t1 = 0;
t2 = 0;
t3 = 0;
for( int k = 0; k < (int) event.size(); k++ ){
//event[k].Print();
if( event[k].ch == 2 ) {dE = event[k].energy; dE_t = event[k].timestamp;}
if( event[k].ch == 2 ) {E = event[k].energy; E_t = event[k].timestamp;}
if( event[k].ch == 2 ) {e0 = event[k].energy; t0 = event[k].timestamp;}
if( event[k].ch == 3 ) {e1= event[k].energy; t1 = event[k].timestamp;}
if( event[k].ch == 4 ) {e2 = event[k].energy; t2 = event[k].timestamp;}
if( event[k].ch == 5 ) {e3= event[k].energy; t3 = event[k].timestamp;}
}
if (e0>10 && e1>10 && e2>10 && e3>10) {
float_t rotation_angle = 31.;
double_t Xr = (((e1+e2)/(e0+e1+e2+e3))-0.51)*cos(-rotation_angle*M_PI/180)-(((e2+e3)/(e0+e1+e2+e3))-0.51)*sin(-rotation_angle*M_PI/180);
double_t Yr = (((e1+e2)/(e0+e1+e2+e3))-0.51)*sin(-rotation_angle*M_PI/180)+(((e2+e3)/(e0+e1+e2+e3))-0.51)*cos(-rotation_angle*M_PI/180);
hPID->Fill(E + RandomGauss(0, 100), dE+ RandomGauss(0, 100)); // x, y
hXY->Fill(((e1+e2)/(e0+e1+e2+e3)),((e2+e3)/(e0+e1+e2+e3)));
hXYr->Fill(Xr,Yr);
hX->Fill(((e1+e2)/(e0+e1+e2+e3)));
hY->Fill(((e2+e3)/(e0+e1+e2+e3)));
hXr->Fill(Xr);
hYr->Fill(Yr);
}
he0->Fill(e0);
he1->Fill(e1);
he2->Fill(e2);
he3->Fill(e3);
ht->Fill(max(max(t0,t1),max(t2,t3))-min(min(t0,t1),min(t2,t3)));
// cout << "t0: " << t0 << endl;
// cout << "t1: " << t1 << endl;
// cout << "t2: " << t2 << endl;
// cout << "t3: " << t3 << endl;
// cout << "time window " << max(max(t0,t1),max(t2,t3))-min(min(t0,t1),min(t2,t3)) <<endl;
// cout <<"---------------------" << endl;
//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]);
}
}
}
hPID->UpdatePlot();
//hXX->UpdatePlot();//
//hYY->UpdatePlot();
hXY->UpdatePlot();
hXYr->UpdatePlot();
hX->UpdatePlot();
hY->UpdatePlot();
hXr->UpdatePlot();
hYr->UpdatePlot();
he0->UpdatePlot();
he1->UpdatePlot();
he2->UpdatePlot();
he3->UpdatePlot();
ht->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