2024-06-04 12:59:18 -04:00
# 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 ) ;
2024-09-17 10:44:53 -04:00
influx = new InfluxDB ( ) ;
2024-09-14 18:17:00 -04:00
SetDatabase ( " https://localhost:8086 " , " testing " , " zKhzKk4Yhf1l9QU-yE2GsIZ1RazqUgoW3NlF8LJqq_xDMwatOJwg1sKrjgq36uLEsQf8Fmn4sJALP7Kkilk14A== " ) ;
2024-06-04 12:59:18 -04:00
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