diff --git a/.gitignore b/.gitignore index 0e45a3f..2eef895 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,9 @@ DataReaderScript data +splitpole.C +splitpole.h + *.d *.pcm diff --git a/ClassData.h b/ClassData.h index 6a481d7..b3327fb 100644 --- a/ClassData.h +++ b/ClassData.h @@ -658,7 +658,7 @@ inline int Data::DecodePHADualChannelBlock(unsigned int ChannelMask, bool fastDe Energy[channel][DataIndex[channel]] = energy; Timestamp[channel][DataIndex[channel]] = timeStamp; - if(extra2Option == 0 || extra2Option == 2 ) fineTime[channel][DataIndex[channel]] = (extra2 & 0x07FF ); + if(extra2Option == 2 ) fineTime[channel][DataIndex[channel]] = (extra2 & 0x03FF ); PileUp[channel][DataIndex[channel]] = pileUp; NumEventsDecoded[channel] ++; TotNumEvents[channel] ++; @@ -850,6 +850,7 @@ inline int Data::DecodePSDDualChannelBlock(unsigned int ChannelMask, bool fastDe Energy[channel][DataIndex[channel]] = Qshort; Energy2[channel][DataIndex[channel]] = Qlong; Timestamp[channel][DataIndex[channel]] = timeStamp; + if( extraOption == 2 ) fineTime[channel][DataIndex[channel]] = extra & 0x3FF; if( !fastDecode ) { if( hasDualTrace ){ diff --git a/DigiSettingsPanel.cpp b/DigiSettingsPanel.cpp index a43e887..bc8ed46 100644 --- a/DigiSettingsPanel.cpp +++ b/DigiSettingsPanel.cpp @@ -2066,7 +2066,9 @@ void DigiSettingsPanel::SetUpPSDChannel(){ SetUpSpinBox(sbFixedBaseline[ID][ch], "", tabLayout, ch + 1, 7, DPP::PSD::FixedBaseline, ch); connect(cbBaseLineAvg[ID][ch], &RComboBox::currentIndexChanged, this, [=](){ - sbFixedBaseline[ID][ch]->setEnabled( cbBaseLineAvg[ID][ch]->currentData().toInt() == 0); + for( int jj = 0; jj < 16 ; jj++ ){ + sbFixedBaseline[ID][jj]->setEnabled( cbBaseLineAvg[ID][jj]->currentData().toInt() == 0); + } }); } @@ -2891,12 +2893,16 @@ void DigiSettingsPanel::SyncAllChannelsTab_PSD(){ SyncComboBox(cbVetoMode); SyncComboBox(cbVetoStep); + for( int jj = 0; jj < 16 ; jj++ ){ + sbFixedBaseline[ID][jj]->setEnabled( cbBaseLineAvg[ID][jj]->currentData().toInt() == 0); + } + } void DigiSettingsPanel::UpdatePSDSetting(){ enableSignalSlot = false; - //printf("------ %s \n", __func__); + // printf("------ %s \n", __func__); for(int ch = 0; ch < digi[ID]->GetNChannels(); ch ++){ @@ -2917,7 +2923,6 @@ void DigiSettingsPanel::UpdatePSDSetting(){ UpdateSpinBox(sbNumEventAgg[ID][ch], DPP::NumberEventsPerAggregate_G, ch); UpdateComboBox(cbDynamicRange[ID][ch], DPP::InputDynamicRange, ch); - uint32_t dpp = digi[ID]->GetSettingFromMemory(DPP::DPPAlgorithmControl, ch); UpdateComboBoxBit(cbPolarity[ID][ch], dpp, DPP::Bit_DPPAlgorithmControl_PSD::Polarity); @@ -2966,6 +2971,7 @@ void DigiSettingsPanel::UpdatePSDSetting(){ UpdateSpinBox(sbVetoWidth[ID][ch], DPP::VetoWidth, ch); UpdateComboBoxBit(cbVetoStep[ID][ch], vetoBit, DPP::Bit_VetoWidth::VetoStep); + } enableSignalSlot = true; diff --git a/FSUDAQ.cpp b/FSUDAQ.cpp index c00e4f0..919e6ee 100644 --- a/FSUDAQ.cpp +++ b/FSUDAQ.cpp @@ -316,10 +316,12 @@ void MainWindow::OpenDataPath(){ leDataPath->setText(fileDialog.selectedFiles().at(0)); rawDataPath = leDataPath->text(); bnStartACQ->setEnabled(true); + bnStartACQ->setStyleSheet("background-color: green;"); }else{ leDataPath->clear(); rawDataPath = ""; bnStartACQ->setEnabled(false); + bnStartACQ->setStyleSheet(""); } SaveProgramSettings(); @@ -551,7 +553,9 @@ void MainWindow::OpenDigitizers(){ LogMsg(QString("Done. Opened %1 digitizer(s).").arg(nDigi)); WaitForDigitizersOpen(false); + bnStartACQ->setStyleSheet("background-color: green;"); bnStopACQ->setEnabled(false); + bnStopACQ->setStyleSheet(""); if( rawDataPath == "" ) { chkSaveData->setChecked(false); @@ -622,6 +626,7 @@ void MainWindow::WaitForDigitizersOpen(bool onOff){ bnOpenScaler->setEnabled(!onOff); bnStartACQ->setEnabled(!onOff); bnStopACQ->setEnabled(!onOff); + bnStopACQ->setStyleSheet(""); chkSaveData->setEnabled(!onOff); bnCanvas->setEnabled(!onOff); bnAnalyzer->setEnabled(!onOff); @@ -849,7 +854,9 @@ void MainWindow::StartACQ(){ if( canvas != nullptr ) histThread->start(); bnStartACQ->setEnabled(false); + bnStartACQ->setStyleSheet(""); bnStopACQ->setEnabled(true); + bnStopACQ->setStyleSheet("background-color: red;"); bnOpenScope->setEnabled(false); if( onlineAnalyzer ) onlineAnalyzer->StartThread(); @@ -918,7 +925,9 @@ void MainWindow::StopACQ(){ lbScalarACQStatus->setText("ACQ Off"); bnStartACQ->setEnabled(true); + bnStartACQ->setStyleSheet("background-color: green;"); bnStopACQ->setEnabled(false); + bnStopACQ->setStyleSheet(""); bnOpenScope->setEnabled(true); {//^=== elog and database @@ -1054,7 +1063,9 @@ void MainWindow::OpenScope(){ connect(scope, &Scope::SendLogMsg, this, &MainWindow::LogMsg); connect(scope, &Scope::CloseWindow, this, [=](){ bnStartACQ->setEnabled(true); + bnStartACQ->setStyleSheet("background-color: green;"); bnStopACQ->setEnabled(false); + bnStopACQ->setStyleSheet(""); }); connect(scope, &Scope::TellACQOnOff, this, [=](bool onOff){ if( scope ) { @@ -1089,7 +1100,9 @@ void MainWindow::OpenScope(){ } bnStartACQ->setEnabled(false); + bnStartACQ->setStyleSheet(""); bnStopACQ->setEnabled(false); + bnStopACQ->setStyleSheet(""); chkSaveData->setChecked(false); } diff --git a/FSUDAQ_Qt6.pro b/FSUDAQ_Qt6.pro index b160239..3cb01d1 100644 --- a/FSUDAQ_Qt6.pro +++ b/FSUDAQ_Qt6.pro @@ -38,6 +38,7 @@ HEADERS += ClassData.h \ MultiBuilder.h \ Analyser.h \ qcustomplot.h \ + Isotope.h \ SplitPoleAnalyzer.h SOURCES += ClassDigitizer.cpp \ DigiSettingsPanel.cpp \ diff --git a/Histogram1D.h b/Histogram1D.h index c1b3010..da2cd6e 100644 --- a/Histogram1D.h +++ b/Histogram1D.h @@ -3,6 +3,8 @@ #include "qcustomplot.h" +#define MaxNHist 10 + //^============================================== //^============================================== class Histogram1D : public QCustomPlot{ @@ -11,6 +13,7 @@ public: Histogram1D(QString title, QString xLabel, int xbin, double xmin, double xmax, QWidget * parent = nullptr) : QCustomPlot(parent){ for( int i = 0; i < 3; i ++) txt[i] = nullptr; + nData = 1; Rebin(xbin, xmin, xmax); xAxis->setLabel(xLabel); @@ -35,7 +38,7 @@ public: connect(xAxis, SIGNAL(rangeChanged(QCPRange)), xAxis2, SLOT(setRange(QCPRange))); connect(yAxis, SIGNAL(rangeChanged(QCPRange)), yAxis2, SLOT(setRange(QCPRange))); - graph(0)->setData(xList, yList); + graph(0)->setData(xList, yList[0]); //setInteractions(QCP::iRangeDrag | QCP::iRangeZoom | QCP::iSelectPlottables); //setInteractions( QCP::iRangeDrag | QCP::iRangeZoom ); @@ -62,7 +65,7 @@ public: connect(this, &QCustomPlot::mouseMove, this, [=](QMouseEvent *event){ double x = this->xAxis->pixelToCoord(event->pos().x()); double bin = (x - xMin)/dX; - double z = yList[2*qFloor(bin) + 1]; + double z = yList[0][2*qFloor(bin) + 1]; QString coordinates = QString("Bin: %1, Value: %2").arg(qFloor(bin)).arg(z); QToolTip::showText(event->globalPosition().toPoint(), coordinates, this); @@ -176,8 +179,22 @@ public: double GetXMin() const {return xMin;} double GetXMax() const {return xMax;} + void SetColor(QColor color, unsigned short ID = 0) { + graph(ID)->setPen(QPen(color)); + QColor haha = color; + haha.setAlpha(20); + graph(ID)->setBrush(QBrush(haha)); + } + void AddDataList(QString title, QColor color){ + nData ++; + addGraph(); + graph(nData - 1)->setName(title); + SetColor(color, nData-1); + yList[nData-1] = yList[0]; + } + void UpdatePlot(){ - graph(0)->setData(xList, yList); + for( int ID = 0 ; ID < nData; ID ++) graph(ID)->setData(xList, yList[ID]); xAxis->setRangeLower(xMin); xAxis->setRangeUpper(xMax); yAxis->setRangeLower(0); @@ -186,7 +203,9 @@ public: } void Clear(){ - for( int i = 0; i <= yList.count(); i++) yList[i] = 0; + for( int ID = 0 ; ID < nData; ID ++) { + for( int i = 0; i <= yList[ID].count(); i++) yList[ID][i] = 0; + } yMax = 0; txt[0]->setText("Under Flow : 0"); txt[1]->setText("Total Entry : 0"); @@ -202,13 +221,15 @@ public: dX = (xMax - xMin)/(xBin); xList.clear(); - yList.clear(); + for( int i = 0 ; i < nData ; i ++) yList[i].clear(); for( int i = 0; i <= xBin; i ++ ){ xList.append(xMin + i*dX-(dX)*0.000001); xList.append(xMin + i*dX); - yList.append(0); - yList.append(0); + for( int ID = 0 ; ID < nData; ID ++ ){ + yList[ID].append(0); + yList[ID].append(0); + } } yMax = 0; @@ -222,34 +243,36 @@ public: if( txt[2] ) txt[2]->setText("Over Flow : 0"); } - void Fill(double value){ - totalEntry ++; - txt[1]->setText("Total Entry : "+ QString::number(totalEntry)); + void Fill(double value, unsigned int ID = 0){ + if( ID == 0 ){ + totalEntry ++; + txt[1]->setText("Total Entry : "+ QString::number(totalEntry)); - if( value < xMin ) { - underFlow ++; - txt[0]->setText("Under Flow : "+ QString::number(underFlow)); - return; - } - if( value > xMax ) { - overFlow ++; - txt[2]->setText("Over Flow : "+ QString::number(overFlow)); - return; + if( value < xMin ) { + underFlow ++; + txt[0]->setText("Under Flow : "+ QString::number(underFlow)); + return; + } + if( value > xMax ) { + overFlow ++; + txt[2]->setText("Over Flow : "+ QString::number(overFlow)); + return; + } } double bin = (value - xMin)/dX; int index1 = 2*qFloor(bin) + 1; int index2 = index1 + 1; - if( 0 <= index1 && index1 <= 2*xBin) yList[index1] += 1; - if( 0 <= index1 && index2 <= 2*xBin) yList[index2] += 1; + if( 0 <= index1 && index1 <= 2*xBin) yList[ID][index1] += 1; + if( 0 <= index1 && index2 <= 2*xBin) yList[ID][index2] += 1; - if( yList[index1] > yMax ) yMax = yList[index1]; + if( yList[ID][index1] > yMax ) yMax = yList[ID][index1]; } - void Print(){ + void Print(unsigned int ID = 0){ for( int i = 0; i < xList.count(); i++){ - printf("%f %f\n", xList[i], yList[i]); + printf("%f %f\n", xList[i], yList[ID][i]); } } @@ -266,12 +289,15 @@ private: int underFlow; int overFlow; - QVector xList, yList; + unsigned short nData; + QVector xList; + QVector yList[MaxNHist]; QCPItemText * txt[3]; bool usingMenu; + }; #endif \ No newline at end of file diff --git a/Isotope.h b/Isotope.h new file mode 100644 index 0000000..d2aa94b --- /dev/null +++ b/Isotope.h @@ -0,0 +1,533 @@ +/*********************************************************************** + * + * This is Isotope.h, To extract the isotope mass from massXX.txt + * + *------------------------------------------------------- + * created by Ryan (Tsz Leung) Tang, Nov-18, 2018 + * email: goluckyryan@gmail.com + * ********************************************************************/ + + +#ifndef ISOTOPE_H +#define ISOTOPE_H + +#include +#include +#include +#include +#include +#include +#include //atoi +#include +using namespace std; + +const double mp = 938.2720813; // MeV/c^2 +const double mn = 939.56542052; // MeV/c^2 +const double amu = 931.0; + +const string massData="mass20.txt"; + +// about the mass**.txt +// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu +// mass excess uncertaintly +// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons +// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons + +class Isotope { +public: + int A, Z; + double Mass, MassError, BEA; + string Name, Symbol; + string dataSource; + + Isotope(){findHeliosPath();}; + Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; + Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; + + void SetIso(int a, int z); + void SetIsoByName(string name); + + double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal + double CalSp2(int a, int z); // this is for (a,z) nucleus removal + + double CalBeta(double T){ + //double Etot = Mass + T; + double gamma = 1. + T/Mass; + double beta = sqrt(1. - 1. / gamma / gamma ) ; + return beta; + } + + void Print(); + void ListShell(); + +private: + void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; + void FindMassByName(string name); // give Z, mass, massError, BEA; + + int TwoJ(int nShell); + string Orbital(int nShell); + int magic(int i){ + switch (i){ + case 0: return 2; break; + case 1: return 8; break; + case 2: return 20; break; + case 3: return 28; break; + case 4: return 40; break; + case 5: return 50; break; + case 6: return 82; break; + case 7: return 128; break; + } + return 0; + } + + int magicShellID(int i){ + switch (i){ + case 0: return 0; break; + case 1: return 2; break; + case 2: return 5; break; + case 3: return 6; break; + case 4: return 9; break; + case 5: return 10; break; + case 6: return 15; break; + case 7: return 21; break; + } + return 0; + } + + int fileStartLine; + int fileEndLine; + int lineMass050_099; + int lineMass100_149; + int lineMass150_199; + int lineMass200; + + + void setFileLines(){ + fileStartLine = 37; + fileEndLine = 3594; + + lineMass050_099 = 466; + lineMass100_149 = 1160; + lineMass150_199 = 1994; + lineMass200 = 2774; + } + + char * heliosPath; + bool isFindOnce; + + void findHeliosPath(){ + heliosPath = getenv("HELIOSSYS"); + if( heliosPath ){ + dataSource = heliosPath; + dataSource += "/analysis" + massData; + }else{ + dataSource = massData; + } + } + +}; + +inline void Isotope::SetIso(int a, int z){ + this->A = a; + this->Z = z; + FindMassByAZ(a,z); +} + +inline void Isotope::SetIsoByName(string name){ + FindMassByName(name); +} + +inline void Isotope::FindMassByAZ(int A, int Z){ + string line; + int lineNum=0; + int list_A, list_Z; + + ifstream myfile; + int flag=0; + + setFileLines(); + + int numLineStart = fileStartLine; + int numLineEnd = fileEndLine; + + if ( A >= 50 && A < 100) numLineStart = lineMass050_099; + if ( A >=100 && A < 150) numLineStart = lineMass100_149; + if ( A >=150 && A < 200) numLineStart = lineMass150_199; + if ( A >=200 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_Z = atoi((line.substr(10,5)).c_str()); + list_A = atoi((line.substr(15,5)).c_str()); + + if ( A == list_A && Z == list_Z) { + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + break; + } + + } + } + + if( this->Name == "1H" ) this->Name = "p"; + if( this->Name == "2H" ) this->Name = "d"; + if( this->Name == "3H" ) this->Name = "t"; + if( this->Name == "4He" ) this->Name = "a"; + + myfile.close(); + }else { + printf("Unable to open %s\n", dataSource.c_str()); + } +} + +inline void Isotope::FindMassByName(string name){ + + // done seperate the Mass number and the name + if( name == "n" ) { + this->Name = "1n"; + this->BEA = 0; + this->Mass = mn; + this->MassError = 0; + this->Name = "n"; + this->A = 1; + this->Z = 0; + return; + } + if( name == "p" ) name = "1H"; + if( name == "d" ) name = "2H"; + if( name == "t" ) name = "3H"; + if( name == "a" ) name = "4He"; + + string temp = name; + int lastDigit = 0; + + for(int i=0; temp[i]; i++){ + if(temp[i] == '0') lastDigit = i; + if(temp[i] == '1') lastDigit = i; + if(temp[i] == '2') lastDigit = i; + if(temp[i] == '3') lastDigit = i; + if(temp[i] == '4') lastDigit = i; + if(temp[i] == '5') lastDigit = i; + if(temp[i] == '6') lastDigit = i; + if(temp[i] == '7') lastDigit = i; + if(temp[i] == '8') lastDigit = i; + if(temp[i] == '9') lastDigit = i; + } + + this->Symbol = temp.erase(0, lastDigit +1); + //check is Symbol is 2 charaters, if not, add " " at the end + if( this->Symbol.length() == 1 ){ + this->Symbol = this->Symbol + " "; + } + + + temp = name; + int len = temp.length(); + temp = temp.erase(lastDigit+1, len); + + this->A = atoi(temp.c_str()); + //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); + + // find the nucleus in the data + string line; + int lineNum=0; + int list_A; + string list_symbol; + + ifstream myfile; + int flag=0; + + setFileLines(); + + int numLineStart = fileStartLine; + int numLineEnd = fileEndLine; + + if ( A >= 50 && A < 100) numLineStart = lineMass050_099; + if ( A >=100 && A < 150) numLineStart = lineMass100_149; + if ( A >=150 && A < 200) numLineStart = lineMass150_199; + if ( A >=200 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_symbol = line.substr(20,2); + list_A = atoi((line.substr(15,5)).c_str()); + + //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str()); + + if ( this->A == list_A && this->Symbol == list_symbol) { + this->Z = atoi((line.substr(10,5)).c_str()); + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << this->A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > this->A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + break; + } + + } + } + myfile.close(); + }else { + printf("Unable to open %s\n", dataSource.c_str()); + } +} + +inline double Isotope::CalSp(int Np, int Nn){ + Isotope nucleusD(A - Np - Nn, Z - Np); + + if( nucleusD.Mass != -404){ + return nucleusD.Mass + Nn*mn + Np*mp - this->Mass; + }else{ + return -404; + } +} + +inline double Isotope::CalSp2(int a, int z){ + Isotope nucleusD(A - a , Z - z); + Isotope nucleusS(a,z); + + if( nucleusD.Mass != -404 && nucleusS.Mass != -404){ + return nucleusD.Mass + nucleusS.Mass - this->Mass; + }else{ + return -404; + } +} + +inline int Isotope::TwoJ(int nShell){ + + switch(nShell){ + case 0: return 1; break; // 0s1/2 + case 1: return 3; break; // 0p3/2 + case 2: return 1; break; // 0p1/2 -- 8 + case 3: return 5; break; // 0d5/2 + case 4: return 1; break; // 1s1/2 + case 5: return 3; break; // 0d3/2 -- 20 + case 6: return 7; break; // 0f7/2 -- 28 + case 7: return 3; break; // 1p3/2 + case 8: return 1; break; // 1p1/2 + case 9: return 5; break; // 0f5/2 -- 40 + case 10: return 9; break; // 0g9/2 -- 50 + case 11: return 7; break; // 0g7/2 + case 12: return 5; break; // 1d5/2 + case 13: return 11; break; // 0h11/2 + case 14: return 3; break; // 1d3/2 + case 15: return 1; break; // 2s1/2 -- 82 + case 16: return 9; break; // 0h9/2 + case 17: return 7; break; // 1f7/2 + case 18: return 13; break; // 0i13/2 + case 19: return 3; break; // 2p3/2 + case 20: return 5; break; // 1f5/2 + case 21: return 1; break; // 1p1/2 -- 126 + case 22: return 9; break; // 1g9/2 + case 23: return 11; break; // 0i11/2 + case 24: return 15; break; // 0j15/2 + case 25: return 5; break; // 2d5/2 + case 26: return 1; break; // 3s1/2 + case 27: return 3; break; // 2d3/2 + case 28: return 7; break; // 1g7/2 + } + + return 0; +} + +inline string Isotope::Orbital(int nShell){ + + switch(nShell){ + case 0: return "0s1 "; break; // + case 1: return "0p3 "; break; // + case 2: return "0p1 "; break; //-- 8 + case 3: return "0d5 "; break; // + case 4: return "1s1 "; break; // + case 5: return "0d3 "; break; //-- 20 + case 6: return "0f7 "; break; //-- 28 + case 7: return "1p3 "; break; // + case 8: return "1p1 "; break; // + case 9: return "0f5 "; break; //-- 40 + case 10: return "0g9 "; break; //-- 50 + case 11: return "0g7 "; break; // + case 12: return "1d5 "; break; // + case 13: return "0h11"; break; // + case 14: return "1d3 "; break; // + case 15: return "2s1 "; break; //-- 82 + case 16: return "0h9 "; break; // + case 17: return "1f7 "; break; // + case 18: return "0i13"; break; // + case 19: return "2p3 "; break; // + case 20: return "1f5 "; break; // + case 21: return "1p1 "; break; //-- 126 + case 22: return "1g9 "; break; // + case 23: return "0i11"; break; // + case 24: return "0j15"; break; // + case 25: return "2d5 "; break; // + case 26: return "3s1 "; break; // + case 27: return "2d3 "; break; // + case 28: return "1g7 "; break; // + } + + return "nan"; +} + +inline void Isotope::ListShell(){ + + if( Mass < 0 ) return; + + int n = A-Z; + int p = Z; + + int k = min(n,p); + int nMagic = 0; + for( int i = 0; i < 7; i++){ + if( magic(i) < k && k <= magic(i+1) ){ + nMagic = i; + break; + } + } + + int coreShell = magicShellID(nMagic-1); + int coreZ1 = magic(nMagic-1); + int coreZ2 = magic(nMagic); + + Isotope core1( 2*coreZ1, coreZ1); + Isotope core2( 2*coreZ2, coreZ2); + + printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str()); + printf(" || "); + int t = max(n,p); + int nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if( nShell > coreShell) { + printf("%4s", Orbital(nShell).c_str()); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + } + t = t - occ; + nShell++; + }while( t > 0 && nShell < 29); + for( int i = 1; i <= 6; i++){ + if (nShell < 28) { + printf("%4s,", Orbital(nShell).c_str()); + }else if( nShell == 28 ) { + printf("%4s", Orbital(nShell).c_str()); + } + nShell ++; + } + if (nShell < 29) printf("%4s", Orbital(nShell).c_str()); + printf("\n"); + + + printf(" Z = %3d || ", p); + nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if( nShell > coreShell ){ + if( p > occ ) { + printf("%-4d", occ); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + }else{ + printf("%-4d", p); + } + } + p = p - occ; + nShell++; + }while( p > 0 && nShell < 29); + printf("\n"); + + printf(" N = %3d || ", n); + nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if ( nShell > coreShell ){ + if( n > occ ) { + printf("%-4d", occ); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + }else{ + printf("%-4d", n); + } + } + n = n - occ; + nShell++; + }while( n > 0 && nShell < 29); + printf("\n"); + + printf("------------------ \n"); +} + +inline void Isotope::Print(){ + + if (Mass > 0){ + + findHeliosPath(); + + printf(" using mass data : %s \n", dataSource.c_str()); + printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.); + printf(" total BE : %12.5f MeV\n",BEA*A/1000.); + printf(" mass in amu : %12.5f u\n",Mass/amu); + printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu); + printf("-------------- Seperation energy \n"); + printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1)); + printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1)); + printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n", CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2)); + printf(" S4p: %8.4f| S4n: %8.4f| \n", CalSp(4, 0), CalSp(0, 4)); + + }else{ + printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A); + } + + +} + +#endif diff --git a/MultiBuilder.cpp b/MultiBuilder.cpp index f15e62c..f423c23 100644 --- a/MultiBuilder.cpp +++ b/MultiBuilder.cpp @@ -207,7 +207,7 @@ void MultiBuilder::BuildEvents(bool isFinal, bool skipTrace, bool verbose){ eventBuilt = 0; //======= Start building event - EventMember em; + Hit em; do{ eventIndex ++; @@ -239,6 +239,8 @@ void MultiBuilder::BuildEvents(bool isFinal, bool skipTrace, bool verbose){ em.ch = ch; em.energy = data[k]->Energy[ch][nextIndex[k][ch]]; em.timestamp = time; + em.fineTime = data[k]->fineTime[ch][nextIndex[k][ch]]; + if( !skipTrace ) em.trace = data[k]->Waveform1[ch][nextIndex[k][ch]]; if( typeList[k] == V1730_DPP_PSD_CODE ) em.energy2 = data[k]->Energy2[ch][nextIndex[k][ch]]; @@ -259,7 +261,7 @@ void MultiBuilder::BuildEvents(bool isFinal, bool skipTrace, bool verbose){ } - std::sort(events[eventIndex].begin(), events[eventIndex].end(), [](const EventMember& a, const EventMember& b) { + std::sort(events[eventIndex].begin(), events[eventIndex].end(), [](const Hit& a, const Hit& b) { return a.timestamp < b.timestamp; }); @@ -308,7 +310,7 @@ void MultiBuilder::BuildEventsBackWard(int maxNumEvent, bool verbose){ //========== build event eventBuilt = 0; - EventMember em; + Hit em; do{ eventIndex ++; if( eventIndex >= MaxNEvent ) eventIndex = 0; @@ -354,7 +356,7 @@ void MultiBuilder::BuildEventsBackWard(int maxNumEvent, bool verbose){ if( timeWindow == 0 ) break; } - std::sort(events[eventIndex].begin(), events[eventIndex].end(), [](const EventMember& a, const EventMember& b) { + std::sort(events[eventIndex].begin(), events[eventIndex].end(), [](const Hit& a, const Hit& b) { return a.timestamp < b.timestamp; }); diff --git a/MultiBuilder.h b/MultiBuilder.h index ae2fa44..b7a12cb 100644 --- a/MultiBuilder.h +++ b/MultiBuilder.h @@ -6,7 +6,7 @@ #define MaxNEvent 5000 // circular -class EventMember{ +class Hit{ public: int sn; unsigned short bd; @@ -18,7 +18,7 @@ public: std::vector trace; - EventMember(){ + Hit(){ Clear(); } @@ -34,7 +34,7 @@ public: } void Print(){ - printf("(%2d, %2d)[%3d] %6d %10llu, %5ld\n", bd, ch, sn, energy, timestamp, trace.size()); + printf("(%2d, %2d)[%3d] %6d %10llu, %6d, %5ld\n", bd, ch, sn, energy, timestamp, fineTime, trace.size()); } }; @@ -67,7 +67,7 @@ public: long eventIndex; long eventBuilt; // reset once call BuildEvents() long totalEventBuilt; - std::vector events[MaxNEvent]; + std::vector events[MaxNEvent]; private: std::vector typeList; diff --git a/RegisterAddress.h b/RegisterAddress.h index b7f1d2c..8de2632 100644 --- a/RegisterAddress.h +++ b/RegisterAddress.h @@ -70,7 +70,13 @@ class Reg{ uint32_t ActualAddress(int ch = -1){ if( address == 0x8180 ) return (ch < 0 ? address : (address + 4*(ch/2))); - if( address < 0x8000 ) return (ch < 0 ? (address + 0x7000) : (address + (ch << 8)) ); + if( address < 0x8000 ){ + if( group ) { + if( ch < 0 ) return address + 0x7000; + return address + ((ch % 2 == 0 ? ch : ch - 1) << 8) ; + } + return (ch < 0 ? (address + 0x7000) : (address + (ch << 8)) ); + } if( address >= 0x8000 ) return address; return 0; } diff --git a/Scope.cpp b/Scope.cpp index 5f163e0..fa7a574 100644 --- a/Scope.cpp +++ b/Scope.cpp @@ -19,6 +19,8 @@ Scope::Scope(Digitizer ** digi, unsigned int nDigi, ReadDataThread ** readDataTh enableSignalSlot = false; + isACQStarted = false; + plot = new RChart(); for( int i = 0; i < MaxNumberOfTrace; i++) { dataTrace[i] = new QLineSeries(); @@ -105,6 +107,10 @@ Scope::Scope(Digitizer ** digi, unsigned int nDigi, ReadDataThread ** readDataTh connect(cbScopeDigi, &RComboBox::currentIndexChanged, this, [=](int index){ if( !enableSignalSlot ) return; + + bool saveACQStartStatus = isACQStarted; + if( isACQStarted) StopScope(); + ID = index; tick2ns = digi[ID]->GetTick2ns(); //---setup cbScopeCh @@ -118,11 +124,20 @@ Scope::Scope(Digitizer ** digi, unsigned int nDigi, ReadDataThread ** readDataTh ReadSettingsFromBoard(); + if( saveACQStartStatus )StartScope(); + }); connect(cbScopeCh, &RComboBox::currentIndexChanged, this, [=](){ if( !enableSignalSlot ) return; + + bool saveACQStartStatus = isACQStarted; + if( isACQStarted) StopScope(); + ReadSettingsFromBoard(); + + if( saveACQStartStatus )StartScope(); + }); @@ -202,7 +217,9 @@ Scope::Scope(Digitizer ** digi, unsigned int nDigi, ReadDataThread ** readDataTh layout->setColumnStretch(5, 1); bnScopeStart->setEnabled(true); + bnScopeStart->setStyleSheet("background-color: green;"); bnScopeStop->setEnabled(false); + bnScopeStop->setStyleSheet(""); UpdatePanelFromMomeory(); @@ -259,12 +276,16 @@ void Scope::StartScope(){ updateTraceThread->start(); bnScopeStart->setEnabled(false); + bnScopeStart->setStyleSheet(""); bnScopeStop->setEnabled(true); + bnScopeStop->setStyleSheet("background-color: red;"); EnableControl(false); TellACQOnOff(true); + isACQStarted = true; + } void Scope::StopScope(){ @@ -292,12 +313,16 @@ void Scope::StopScope(){ emit UpdateOtherPanels(); bnScopeStart->setEnabled(true); + bnScopeStart->setStyleSheet("background-color: green;"); bnScopeStop->setEnabled(false); + bnScopeStop->setStyleSheet(""); EnableControl(true); TellACQOnOff(false); + isACQStarted = false; + // printf("----- end of %s\n", __func__); } @@ -351,6 +376,13 @@ void Scope::UpdateScope(){ } digiMTX[ID].unlock(); + if( data->TriggerRate[ch] == 0 ){ + dataTrace[0]->clear(); + dataTrace[1]->clear(); + dataTrace[2]->clear(); + dataTrace[3]->clear(); + } + plot->axes(Qt::Horizontal).first()->setRange(0, tick2ns * traceLength * factor); emit UpdateScaler(); @@ -457,6 +489,7 @@ void Scope::SetUpSpinBox(RSpinBox * &sb, QString str, int row, int col, const Re digiMTX[ID].lock(); digi[ID]->WriteRegister(para, value, ch); + digiMTX[ID].unlock(); if( digi[ID]->GetErrorCode() == CAEN_DGTZ_Success ){ SendLogMsg(msg + " | OK."); @@ -610,7 +643,8 @@ void Scope::SetUpPSDPanel(){ SetUpSpinBox(sbGateOffset, "Gate Offset [ns] ", rowID, 6, DPP::PSD::GateOffset); rowID ++; //============================================================= - SetUpComboBoxSimple(cbAnaProbe1, "Ana. Probe ", rowID, 0); + SetUpSpinBox(sbThreshold, "Threshold [LSB] ", rowID, 0, DPP::PSD::TriggerThreshold); + SetUpComboBoxSimple(cbAnaProbe1, "Ana. Probe ", rowID, 2); for( int i = 0; i < (int) DPP::Bit_BoardConfig::ListAnaProbe_PSD.size(); i++){ cbAnaProbe1->addItem(QString::fromStdString(DPP::Bit_BoardConfig::ListAnaProbe_PSD[i].first), DPP::Bit_BoardConfig::ListAnaProbe_PSD[i].second); } diff --git a/Scope.h b/Scope.h index b482ee2..2fc696e 100644 --- a/Scope.h +++ b/Scope.h @@ -71,6 +71,7 @@ private: Digitizer ** digi; unsigned short nDigi; unsigned short ID; // the id of digi, index of cbScopeDigi + bool isACQStarted; int tick2ns; bool traceOn[MaxNDigitizer]; diff --git a/SplitPoleAnalyzer.h b/SplitPoleAnalyzer.h index 14b6a2d..db8b5fe 100644 --- a/SplitPoleAnalyzer.h +++ b/SplitPoleAnalyzer.h @@ -15,10 +15,142 @@ * >make * * ******************************************/ - - #include "Analyser.h" +#include "Isotope.h" +namespace ChMap{ + + const short ScinR = 0; + const short ScinL = 1; + const short dFR = 9; + const short dFL = 8; + const short dBR = 10; + const short dBL = 11; + const short Cathode = 7; + const short AnodeF = 13; + const short AnodeB = 15; + +}; + +const double c = 299.792458; // mm/ns +const double pi = M_PI; +const double deg2rad = pi/180.; + +class SplitPoleHit{ + +public: + SplitPoleHit(){ + + target.SetIso(12, 6); + beam.SetIso(2,1); + recoil.SetIso(1,1); + + Bfield = 0.76; // Tesla + angleDegree = 20; // degree + beamKE = 16; // MeV + + heavyRecoil.SetIso(target.A + beam.A - recoil.A, target.Z + beam.Z - recoil.Z); + + double Q = target.Mass + beam.Mass - recoil.Mass - heavyRecoil.Mass; + + double haha1 = sqrt(beam.Mass + beamKE + recoil.Mass)/(recoil.Mass + heavyRecoil.Mass) / cos(angleDegree * deg2rad); + double haha2 = ( beamKE * ( heavyRecoil.Mass + beam.Mass) + heavyRecoil.Mass * Q) / (recoil.Mass + heavyRecoil.Mass); + + double recoilKE = pow(haha1 + sqrt(haha1*haha1 + haha2), 2); + + printf("Q value : %f \n", Q); + printf("proton enegry : %f \n", recoilKE); + + double recoilP = sqrt( recoilKE* ( recoilKE + 2*recoil.Mass)); + double rho = recoilP/(target.Z * Bfield * c); // in m + double haha = sqrt( recoil.Mass * beam.Mass * beamKE / recoilKE ); + double k = haha * sin(angleDegree * deg2rad) / ( recoil.Mass + heavyRecoil.Mass - haha * cos(angleDegree * deg2rad)); + + const double SPS_DISPERSION = 1.96; // x-position/rho + const double SPS_MAGNIFICATION = 0.39; // in x-position + + zOffset = -1000.0 * rho * k * SPS_DISPERSION * SPS_MAGNIFICATION; + + printf("rho: %f m; z-offset: %f mm\n", rho, zOffset); + + Clear(); + } + + unsigned int eSR; unsigned long long tSR; + unsigned int eSL; unsigned long long tSL; + unsigned int eFR; unsigned long long tFR; + unsigned int eFL; unsigned long long tFL; + unsigned int eBR; unsigned long long tBR; + unsigned int eBL; unsigned long long tBL; + unsigned int eCath; unsigned long long tCath; + unsigned int eAF; unsigned long long tAF; + unsigned int eAB; unsigned long long tAB; + + float eSAvg; + float x1, x2, theta; + float xAvg; + + void Clear(){ + eSR = 0; tSR = 0; + eSL = 0; tSL = 0; + eFR = 0; tFR = 0; + eFL = 0; tFL = 0; + eBR = 0; tBR = 0; + eBL = 0; tBL = 0; + eCath = 0; tCath = 0; + eAF = 0; tAF = 0; + eAB = 0; tAB = 0; + + eSAvg = -1; + x1 = NAN; + x2 = NAN; + theta = NAN; + xAvg = NAN; + } + + void CalData(){ + + if( eSR > 0 && eSL > 0 ) eSAvg = (eSR + eSL)/2; + if( eSR > 0 && eSL == 0 ) eSAvg = eSR; + if( eSR == 0 && eSL > 0 ) eSAvg = eSL; + + if( tFR > 0 && tFL > 0 ) x1 = (tFL - tFR)/2./2.1; + if( tBR > 0 && tBL > 0 ) x2 = (tBL - tBR)/2./1.98; + + if( !std::isnan(x1) && !std::isnan(x2)) { + + if( x2 > x1 ) { + theta = atan((x2-x1)/36.0); + }else if(x2 < x1){ + theta = pi + atan((x2-x1)/36.0); + }else{ + theta = pi * 0.5; + } + + double w1 = 0.5 - zOffset/4.28625; + xAvg = w1 * x1 + (1-w1)* x2; + + } + + } + +private: + + Isotope target; + Isotope beam; + Isotope recoil; + Isotope heavyRecoil; + + double Bfield; + double angleDegree; + double beamKE; + + double zOffset; + +}; + +//^=========================================== +//^=========================================== class SplitPole : public Analyzer{ Q_OBJECT public: @@ -37,6 +169,9 @@ public: dataBaseName = "testing"; SetUpCanvas(); + + hit.Clear(); + } /// ~SplitPole(); // comment out = defalt destructor @@ -51,26 +186,39 @@ private: MultiBuilder *evtbder; // declaie histograms - Histogram2D * h2; + Histogram2D * hPID; + Histogram1D * h1; + Histogram1D * h1g; + Histogram1D * hMulti; int tick2ns; + SplitPoleHit hit; + }; inline void SplitPole::SetUpCanvas(){ - setGeometry(0, 0, 1600, 800); + setGeometry(0, 0, 1600, 1000); // the "this" make the histogram a child of the SplitPole class. When SplitPole destory, all childs destory as well. - h2 = new Histogram2D("Split Pole PID", "x", "y", 100, 0, 10000, 100, 0, 10000, this); + hPID = new Histogram2D("Split Pole PID", "Scin-L", "Anode-Font", 100, 0, 2000, 100, 0, 2000, this); //layout is inheriatge from Analyzer - layout->addWidget(h2, 0, 0); + layout->addWidget(hPID, 0, 0, 2, 1); - h1 = new Histogram1D("Spectrum", "x", 400, 0, 10000, this); + h1 = new Histogram1D("Spectrum", "x", 100, 0, 2000, this); + h1->SetColor(Qt::darkGreen); + h1->AddDataList("Test", Qt::red); // add another histogram in h1, Max Data List is 10 layout->addWidget(h1, 0, 1); + hMulti = new Histogram1D("Multiplicity", "", 10, 0, 10, this); + layout->addWidget(hMulti, 1, 1); + + h1g = new Histogram1D("Spectrum (gated)", "x", 100, 0, 2000, this); + layout->addWidget(h1g, 2, 0, 1, 2); + } @@ -83,7 +231,7 @@ inline void SplitPole::UpdateHistograms(){ if( eventBuilt == 0 ) return; //============ Get the cut list, if any - QList cutList = h2->GetCutList(); + QList cutList = hPID->GetCutList(); const int nCut = cutList.count(); unsigned long long tMin[nCut] = {0xFFFFFFFFFFFFFFFF}, tMax[nCut] = {0}; unsigned int count[nCut]={0}; @@ -94,38 +242,54 @@ inline void SplitPole::UpdateHistograms(){ if(eventStart < 0 ) eventStart += MaxNEvent; for( long i = eventStart ; i <= eventIndex; i ++ ){ - unsigned short e1 = 0, e2 = 0; - unsigned long long t1 = 0, t2 = 0; - std::vector event = evtbder->events[i]; + std::vector event = evtbder->events[i]; //printf("-------------- %ld\n", i); + + hMulti->Fill((int) event.size()); + //if( event.size() < 9 ) return; + if( event.size() == 0 ) return; + + hit.Clear(); + for( int k = 0; k < (int) event.size(); k++ ){ //event[k].Print(); - if( event[k].ch == 9 ) {e1 = event[k].energy; t1 = event[k].timestamp;} - if( event[k].ch == 10 ) {e2 = event[k].energy; t2 = event[k].timestamp;} + if( event[k].ch == ChMap::ScinR ) {hit.eSR = event[k].energy; hit.tSR = event[k].timestamp;} + if( event[k].ch == ChMap::ScinL ) {hit.eSL = event[k].energy; hit.tSL = event[k].timestamp;} + if( event[k].ch == ChMap::dFR ) {hit.eFR = event[k].energy; hit.tFR = event[k].timestamp;} + if( event[k].ch == ChMap::dFL ) {hit.eFL = event[k].energy; hit.tFL = event[k].timestamp;} + if( event[k].ch == ChMap::dBR ) {hit.eBL = event[k].energy; hit.tBL = event[k].timestamp;} + if( event[k].ch == ChMap::dBL ) {hit.eBL = event[k].energy; hit.tBL = event[k].timestamp;} + if( event[k].ch == ChMap::Cathode ) {hit.eCath = event[k].energy; hit.tCath = event[k].timestamp;} + if( event[k].ch == ChMap::AnodeF ) {hit.eAF = event[k].energy; hit.tAF = event[k].timestamp;} + if( event[k].ch == ChMap::AnodeB ) {hit.eAB = event[k].energy; hit.tAB = event[k].timestamp;} } - if( e1 == 0 ) continue; - if( e2 == 0 ) continue; + hit.CalData(); - h2->Fill(e1, e2); - h1->Fill(e1); + hPID->Fill(hit.eSL, hit.eSR); // x, y - //check events inside any Graphical cut and extract the rate, using t1 only + h1->Fill(hit.eSL); + h1->Fill(hit.eSR, 1); + + //check events inside any Graphical cut and extract the rate, using tSR only for(int p = 0; p < cutList.count(); p++ ){ if( cutList[p].isEmpty() ) continue; - if( cutList[p].containsPoint(QPointF(e1, e2), Qt::OddEvenFill) ){ - if( t1 < tMin[p] ) tMin[p] = t1; - if( t1 > tMax[p] ) tMax[p] = t1; + if( cutList[p].containsPoint(QPointF(hit.eSL, hit.eSR), Qt::OddEvenFill) ){ + if( hit.tSR < tMin[p] ) tMin[p] = hit.tSR; + if( hit.tSR > tMax[p] ) tMax[p] = hit.tSR; count[p] ++; //printf(".... %d \n", count[p]); + if( p == 0 ) h1g->Fill(hit.eSR); } } } - h2->UpdatePlot(); + hPID->UpdatePlot(); h1->UpdatePlot(); + hMulti->UpdatePlot(); + h1g->UpdatePlot(); - QList cutNameList = h2->GetCutNameList(); + QList 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