Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jlab/hallc/analyzer_software/hcana
  • whit/hcana
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2303 additions and 515 deletions
src/PRadException.cxx deleted 100644 → 0
View file @ dfb590c3
//============================================================================//
// A exception class for PRad Event Viewer //
// //
// Chao Peng //
// 02/27/2016 //
//============================================================================//
#include "PRadException.h"
using namespace std;
PRadException::PRadException(const string &typ, const string &txt, const string &aux)
: title(typ), text(txt), auxText(aux)
{
}
PRadException::PRadException(PRadExceptionType typ, const string &txt, const string &aux)
: type(typ), text(txt), auxText(aux)
{
}
PRadException::PRadException(PRadExceptionType typ, const string &txt, const string &file, const string &func, int line)
: type(typ), text(txt)
{
ostringstream oss;
oss << " evioException occured in file " << file << ", function " << func << ", line " << line;
auxText=oss.str();
}
string PRadException::FailureDesc(void) const
{
ostringstream oss;
oss << text << endl
<< auxText;
return(oss.str());
}
string PRadException::FailureType(void) const
{
if(!title.empty())
return title;
string oss;
switch(type)
{
case ET_CONNECT_ERROR:
oss = "ET CONNECT ERROR";
break;
case ET_CONFIG_ERROR:
oss = "ET CONFIG ERROR";
break;
case ET_STATION_CONFIG_ERROR:
oss = "ET STATION CONFIG ERROR";
break;
case ET_STATION_CREATE_ERROR:
oss = "ET STATION CREATE ERROR";
break;
case ET_STATION_ATTACH_ERROR:
oss = "ET ATTACH ERROR";
break;
case ET_READ_ERROR:
oss = "ET READ ERROR";
break;
case ET_PUT_ERROR:
oss = "ET PUT ERROR";
break;
case HIGH_VOLTAGE_ERROR:
oss = "HIGH VOLTAGE SYSTEM ERROR";
break;
default:
oss = "UNKNOWN ERROR";
break;
}
return(oss);
}
const char* PRadException::what() const noexcept
{
string failure = FailureType() + ": " + FailureDesc();
return failure.c_str();
}
src/PRadException.h deleted 100644 → 0
View file @ dfb590c3
#ifndef PRAD_EXCEPTION_H
#define PRAD_EXCEPTION_H
#include <stdlib.h>
#include <string.h>
#include <exception>
#include <string>
#include <sstream>
class PRadException : public std::exception
{
public:
enum PRadExceptionType
{
UNKNOWN_ERROR,
ET_CONNECT_ERROR,
ET_CONFIG_ERROR,
ET_STATION_CONFIG_ERROR,
ET_STATION_CREATE_ERROR,
ET_STATION_ATTACH_ERROR,
ET_READ_ERROR,
ET_PUT_ERROR,
HIGH_VOLTAGE_ERROR,
};
PRadException(const std::string &typ, const std::string &txt = "", const std::string &aux = "");
PRadException(PRadExceptionType typ = UNKNOWN_ERROR, const std::string &txt = "", const std::string &aux = "");
PRadException(PRadExceptionType typ, const std::string &txt, const std::string &file, const std::string &func, int line);
virtual ~PRadException(void) {}
virtual std::string FailureDesc(void) const;
virtual std::string FailureType(void) const;
const char *what() const noexcept;
public:
PRadExceptionType type; // exception type
std::string title;
std::string text; // primary text
std::string auxText; // auxiliary text
};
#endif
src/THcAerogel.cxx
View file @ e974befb
......@@ -39,7 +39,8 @@ THcAerogel::THcAerogel( const char* name, const char* description,
THaNonTrackingDetector(name,description,apparatus),
fPresentP(0),
fAdcPosTimeWindowMin(0), fAdcPosTimeWindowMax(0), fAdcNegTimeWindowMin(0),
fAdcNegTimeWindowMax(0), fRegionValue(0), fPosGain(0), fNegGain(0),
fAdcNegTimeWindowMax(0),fPedNegDefault(0),fPedPosDefault(0),
fRegionValue(0), fPosGain(0), fNegGain(0),
frPosAdcPedRaw(0), frPosAdcPulseIntRaw(0), frPosAdcPulseAmpRaw(0),
frPosAdcPulseTimeRaw(0), frPosAdcPed(0), frPosAdcPulseInt(0),
frPosAdcPulseAmp(0), frPosAdcPulseTime(0), frNegAdcPedRaw(0),
......@@ -59,7 +60,8 @@ THcAerogel::THcAerogel( const char* name, const char* description,
THcAerogel::THcAerogel( ) :
THaNonTrackingDetector(),
fAdcPosTimeWindowMin(0), fAdcPosTimeWindowMax(0), fAdcNegTimeWindowMin(0),
fAdcNegTimeWindowMax(0), fRegionValue(0), fPosGain(0), fNegGain(0),
fAdcNegTimeWindowMax(0),
fPedNegDefault(0),fPedPosDefault(0),fRegionValue(0), fPosGain(0), fNegGain(0),
frPosAdcPedRaw(0), frPosAdcPulseIntRaw(0), frPosAdcPulseAmpRaw(0),
frPosAdcPulseTimeRaw(0), frPosAdcPed(0), frPosAdcPulseInt(0),
frPosAdcPulseAmp(0), frPosAdcPulseTime(0), frNegAdcPedRaw(0),
......@@ -111,6 +113,8 @@ void THcAerogel::DeleteArrays()
delete [] fAdcPosTimeWindowMax; fAdcPosTimeWindowMax = 0;
delete [] fAdcNegTimeWindowMin; fAdcNegTimeWindowMin = 0;
delete [] fAdcNegTimeWindowMax; fAdcNegTimeWindowMax = 0;
delete [] fPedNegDefault; fPedNegDefault = 0;
delete [] fPedPosDefault; fPedPosDefault = 0;
// 6 GeV variables
delete fPosTDCHits; fPosTDCHits = NULL;
......@@ -296,6 +300,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
fAdcPosTimeWindowMax = new Double_t [fNelem];
fAdcNegTimeWindowMin = new Double_t [fNelem];
fAdcNegTimeWindowMax = new Double_t [fNelem];
fPedNegDefault = new Int_t [fNelem];
fPedPosDefault = new Int_t [fNelem];
DBRequest list[]={
{"aero_num_regions", &fNRegions, kInt},
......@@ -315,6 +321,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
{"aero_adcPosTimeWindowMax", fAdcPosTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem), 1},
{"aero_adcNegTimeWindowMin", fAdcNegTimeWindowMin, kDouble, static_cast<UInt_t>(fNelem), 1},
{"aero_adcNegTimeWindowMax", fAdcNegTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem), 1},
{"aero_PedNegDefault", fPedNegDefault, kInt, static_cast<UInt_t>(fNelem), 1},
{"aero_PedPosDefault", fPedPosDefault, kInt, static_cast<UInt_t>(fNelem), 1},
{"aero_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
{"aero_debug_adc", &fDebugAdc, kInt, 0, 1},
{"aero_six_gev_data", &fSixGevData, kInt, 0, 1},
......@@ -331,6 +339,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
fAdcNegTimeWindowMin[ip] = -1000.;
fAdcPosTimeWindowMax[ip] = 1000.;
fAdcNegTimeWindowMax[ip] = 1000.;
fPedNegDefault[ip] = 0.;
fPedPosDefault[ip] = 0.;
}
fSixGevData = 0; // Set 6 GeV data parameter to false unless set in parameter file
......@@ -646,6 +656,21 @@ Int_t THcAerogel::Decode( const THaEvData& evdata )
if (rawPosAdcHit.GetPulseAmpRaw(thit) > 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 0);
if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 1);
if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawPosAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawPosAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawPosAdcHit.GetAdcTopC();
Double_t AdcToV = rawPosAdcHit.GetAdcTomV();
if (fPedPosDefault[npmt-1] !=0) {
Double_t tPulseInt = AdcToC*(rawPosAdcHit.GetPulseIntRaw(thit) - fPedPosDefault[npmt-1]*PeakPedRatio);
((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(npmt, tPulseInt);
((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, fPedPosDefault[npmt-1]);
((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(npmt, float(fPedPosDefault[npmt-1])/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(npmt, 0.);
}
++nrPosAdcHits;
fTotNumAdcHits++;
......@@ -669,6 +694,21 @@ Int_t THcAerogel::Decode( const THaEvData& evdata )
if (rawNegAdcHit.GetPulseAmpRaw(thit) > 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 0);
if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 1);
if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawNegAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawNegAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawNegAdcHit.GetAdcTopC();
Double_t AdcToV = rawNegAdcHit.GetAdcTomV();
if (fPedNegDefault[npmt-1] !=0) {
Double_t tPulseInt = AdcToC*(rawNegAdcHit.GetPulseIntRaw(thit) - fPedNegDefault[npmt-1]*PeakPedRatio);
((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(npmt, tPulseInt);
((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, fPedNegDefault[npmt-1]);
((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(npmt, float(fPedNegDefault[npmt-1])/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(npmt, 0.);
}
++nrNegAdcHits;
fTotNumAdcHits++;
fTotNumNegAdcHits++;
......@@ -690,6 +730,8 @@ Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
{
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
Double_t OffsetTime = 0.0;
if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
//cout << " starttime = " << StartTime << endl;
// Loop over the elements in the TClonesArray
for(Int_t ielem = 0; ielem < frPosAdcPulseInt->GetEntries(); ielem++) {
......@@ -700,18 +742,13 @@ Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
//// Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcPosTimeWindowMin[npmt] && adctdcdiffTime < fAdcPosTimeWindowMax[npmt];
// By default, the last hit within the timing cut will be considered "good"
if (!errorFlag)
{
fGoodPosAdcMult.at(npmt) += 1;
}
if (!errorFlag && pulseTimeCut) {
if (pulseTimeCut) {
fGoodPosAdcPed.at(npmt) = pulsePed;
// cout << " out = " << npmt << " " << frPosAdcPulseInt->GetEntries() << " " <<fGoodPosAdcMult.at(npmt);
fGoodPosAdcPulseInt.at(npmt) = pulseInt;
......@@ -738,17 +775,13 @@ Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
//// Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcNegTimeWindowMin[npmt] && adctdcdiffTime < fAdcNegTimeWindowMax[npmt];
if (!errorFlag)
{
fGoodNegAdcMult.at(npmt) += 1;
}
// By default, the last hit within the timing cut will be considered "good"
if (!errorFlag && pulseTimeCut) {
if (pulseTimeCut) {
fGoodNegAdcPed.at(npmt) = pulsePed;
fGoodNegAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodNegAdcPulseAmp.at(npmt) = pulseAmp;
......
src/THcAerogel.h
View file @ e974befb
......@@ -73,6 +73,8 @@ class THcAerogel : public THaNonTrackingDetector, public THcHitList {
Double_t *fAdcPosTimeWindowMax;
Double_t *fAdcNegTimeWindowMin;
Double_t *fAdcNegTimeWindowMax;
Int_t* fPedNegDefault;
Int_t* fPedPosDefault;
Double_t fAdcTdcOffset;
Double_t *fRegionValue;
......
src/THcCherenkov.cxx
View file @ e974befb
......@@ -129,27 +129,20 @@ void THcCherenkov::DeleteArrays() {
fGain = NULL;
// 6 Gev variables
delete[] fPedSum;
fPedSum = NULL;
delete[] fPedSum2;
fPedSum2 = NULL;
delete[] fPedLimit;
fPedLimit = NULL;
delete[] fPedMean;
fPedMean = NULL;
delete[] fPedCount;
fPedCount = NULL;
delete[] fPed;
fPed = NULL;
delete[] fThresh;
fThresh = NULL;
delete[] fAdcTimeWindowMin;
fAdcTimeWindowMin = 0;
delete[] fAdcTimeWindowMax;
fAdcTimeWindowMax = 0;
delete[] fRegionValue;
fRegionValue = 0;
delete [] fPedSum; fPedSum = NULL;
delete [] fPedSum2; fPedSum2 = NULL;
delete [] fPedLimit; fPedLimit = NULL;
delete [] fPedMean; fPedMean = NULL;
delete [] fPedCount; fPedCount = NULL;
delete [] fPed; fPed = NULL;
delete [] fThresh; fThresh = NULL;
delete [] fPedDefault; fPedDefault = 0;
delete [] fAdcTimeWindowMin; fAdcTimeWindowMin = 0;
delete [] fAdcTimeWindowMax; fAdcTimeWindowMax = 0;
delete [] fRegionValue; fRegionValue = 0;
}
//_____________________________________________________________________________
......@@ -212,39 +205,48 @@ Int_t THcCherenkov::ReadDatabase(const TDatime& date) {
// << GetName() << " with " << fNelem << " PMTs" << endl;
// 6 GeV pedestal paramters
fPedLimit = new Int_t[fNelem];
fGain = new Double_t[fNelem];
fPedMean = new Double_t[fNelem];
fAdcTimeWindowMin = new Double_t[fNelem];
fAdcTimeWindowMax = new Double_t[fNelem];
fPedLimit = new Int_t[fNelem];
fGain = new Double_t[fNelem];
fPedMean = new Double_t[fNelem];
fAdcTimeWindowMin = new Double_t[fNelem];
fAdcTimeWindowMax= new Double_t[fNelem];
fPedDefault= new Int_t[fNelem];
// Region parameters
fRegionsValueMax = fNRegions * 8;
fRegionValue = new Double_t[fRegionsValueMax];
fAdcGoodElem = new Int_t[fNelem];
fAdcPulseAmpTest = new Double_t[fNelem];
DBRequest list[] = {{"_ped_limit", fPedLimit, kInt, (UInt_t)fNelem, optional},
{"_adc_to_npe", fGain, kDouble, (UInt_t)fNelem},
{"_red_chi2_min", &fRedChi2Min, kDouble},
{"_red_chi2_max", &fRedChi2Max, kDouble},
{"_beta_min", &fBetaMin, kDouble},
{"_beta_max", &fBetaMax, kDouble},
{"_enorm_min", &fENormMin, kDouble},
{"_enorm_max", &fENormMax, kDouble},
{"_dp_min", &fDpMin, kDouble},
{"_dp_max", &fDpMax, kDouble},
{"_mirror_zpos", &fMirrorZPos, kDouble},
{"_npe_thresh", &fNpeThresh, kDouble},
{"_debug_adc", &fDebugAdc, kInt, 0, 1},
{"_adcTimeWindowMin", fAdcTimeWindowMin, kDouble, (UInt_t)fNelem, 1},
{"_adcTimeWindowMax", fAdcTimeWindowMax, kDouble, (UInt_t)fNelem, 1},
{"_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
{"_region", &fRegionValue[0], kDouble, (UInt_t)fRegionsValueMax},
{"_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
{0}};
for (Int_t i = 0; i < fNelem; i++) {
fAdcTimeWindowMin[i] = -1000.;
fAdcTimeWindowMax[i] = 1000.;
DBRequest list[]={
{"_ped_limit", fPedLimit, kInt, (UInt_t) fNelem, optional},
{"_adc_to_npe", fGain, kDouble, (UInt_t) fNelem},
{"_red_chi2_min", &fRedChi2Min, kDouble},
{"_red_chi2_max", &fRedChi2Max, kDouble},
{"_beta_min", &fBetaMin, kDouble},
{"_beta_max", &fBetaMax, kDouble},
{"_enorm_min", &fENormMin, kDouble},
{"_enorm_max", &fENormMax, kDouble},
{"_dp_min", &fDpMin, kDouble},
{"_dp_max", &fDpMax, kDouble},
{"_mirror_zpos", &fMirrorZPos, kDouble},
{"_npe_thresh", &fNpeThresh, kDouble},
{"_debug_adc", &fDebugAdc, kInt, 0, 1},
{"_adcTimeWindowMin", fAdcTimeWindowMin, kDouble,(UInt_t) fNelem,1},
{"_adcTimeWindowMax", fAdcTimeWindowMax, kDouble, (UInt_t) fNelem,1},
{"_PedDefault", fPedDefault, kInt, (UInt_t) fNelem,1},
{"_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
{"_region", &fRegionValue[0], kDouble, (UInt_t) fRegionsValueMax},
{"_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
{0}
};
for (Int_t i=0;i<fNelem;i++) {
fAdcTimeWindowMin[i]=-1000.;
fAdcTimeWindowMax[i]=1000.;
fPedDefault[i]=0;
}
fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file
fAdcTdcOffset = 0.0;
......@@ -298,32 +300,36 @@ Int_t THcCherenkov::DefineVariables(EMode mode) {
} // end debug statement
RVarDef vars[] = {
{"adcCounter", "ADC counter numbers", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"},
{"numGoodAdcHits", "Number of Good ADC Hits Per PMT",
"fNumGoodAdcHits"}, // Cherenkov occupancy
{"totNumGoodAdcHits", "Total Number of Good ADC Hits",
"fTotNumGoodAdcHits"}, // Cherenkov multiplicity
{"numTracksMatched", "Number of Tracks Matched Per Region", "fNumTracksMatched"},
{"numTracksFired", "Number of Tracks that Fired Per Region", "fNumTracksFired"},
{"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"},
{"totNumTracksFired", "Total Number of Tracks that Fired", "fTotNumTracksFired"},
{"xAtCer", "Track X at Cherenkov mirror", "fXAtCer"},
{"yAtCer", "Track Y at Cherenkov mirror", "fYAtCer"},
{"npe", "Number of PEs", "fNpe"},
{"npeSum", "Total Number of PEs", "fNpeSum"},
{"goodAdcPed", "Good ADC pedestals", "fGoodAdcPed"},
{"goodAdcMult", "Good ADC Multiplicity", "fGoodAdcMult"},
{"goodAdcHitUsed", "Good ADC Hit Used", "fGoodAdcHitUsed"},
{"goodAdcPulseInt", "Good ADC pulse integrals", "fGoodAdcPulseInt"},
{"goodAdcPulseIntRaw", "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"},
{"goodAdcPulseAmp", "Good ADC pulse amplitudes", "fGoodAdcPulseAmp"},
{"goodAdcPulseTime", "Good ADC pulse times", "fGoodAdcPulseTime"},
{"goodAdcTdcDiffTime", "Good Hodo Start - ADC pulse times", "fGoodAdcTdcDiffTime"},
{0}};
{"adcCounter", "ADC counter numbers", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"},
{"numGoodAdcHits", "Number of Good ADC Hits Per PMT", "fNumGoodAdcHits"}, // Cherenkov occupancy
{"totNumGoodAdcHits", "Total Number of Good ADC Hits", "fTotNumGoodAdcHits"}, // Cherenkov multiplicity
{"numTracksMatched", "Number of Tracks Matched Per Region", "fNumTracksMatched"},
{"numTracksFired", "Number of Tracks that Fired Per Region", "fNumTracksFired"},
{"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"},
{"totNumTracksFired", "Total Number of Tracks that Fired", "fTotNumTracksFired"},
{"xAtCer", "Track X at Cherenkov mirror", "fXAtCer"},
{"yAtCer", "Track Y at Cherenkov mirror", "fYAtCer"},
{"npe", "Number of PEs", "fNpe"},
{"npeSum", "Total Number of PEs", "fNpeSum"},
{"goodAdcPed", "Good ADC pedestals", "fGoodAdcPed"},
{"goodAdcMult", "Good ADC Multiplicity", "fGoodAdcMult"},
{"goodAdcHitUsed", "Good ADC Hit Used", "fGoodAdcHitUsed"},
{"goodAdcPulseInt", "Good ADC pulse integrals", "fGoodAdcPulseInt"},
{"goodAdcPulseIntRaw", "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"},
{"goodAdcPulseAmp", "Good ADC pulse amplitudes", "fGoodAdcPulseAmp"},
{"goodAdcPulseTime", "Good ADC pulse times", "fGoodAdcPulseTime"},
{"goodAdcTdcDiffTime", "Good Hodo Start - ADC pulse times", "fGoodAdcTdcDiffTime"},
{"RefTime", "Raw ADC RefTime (chan) ", "fRefTime"}, // Raw reference time
{ 0 }
};
return DefineVarsFromList(vars, mode);
}
......@@ -360,6 +366,7 @@ void THcCherenkov::Clear(Option_t* opt) {
fYAtCer = 0.0;
fNpeSum = 0.0;
fRefTime=kBig;
frAdcPedRaw->Clear();
frAdcPulseIntRaw->Clear();
......@@ -406,7 +413,12 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
}
fNhits = DecodeToHitList(evdata, !present);
if (gHaCuts->Result("Pedestal_event")) {
//THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
// cout << "Cerenkov Event num = " << evdata.GetEvNum() << " spec = " << app->GetName() << endl;
if(gHaCuts->Result("Pedestal_event")) {
AccumulatePedestals(fRawHitList);
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return (0);
......@@ -421,14 +433,16 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
UInt_t nrAdcHits = 0;
_waveforms.clear();
while (ihit < fNhits) {
THcCherenkovHit* hit = (THcCherenkovHit*)fRawHitList->At(ihit);
Int_t npmt = hit->fCounter;
THcRawAdcHit& rawAdcHit = hit->GetRawAdcHitPos();
_waveforms.push_back({rawAdcHit.GetSampleBuffer()});
while(ihit < fNhits) {
THcCherenkovHit* hit = (THcCherenkovHit*) fRawHitList->At(ihit);
Int_t npmt = hit->fCounter;
THcRawAdcHit& rawAdcHit = hit->GetRawAdcHitPos();
if (rawAdcHit.GetNPulses() >0 && rawAdcHit.HasRefTime()) {
fRefTime=rawAdcHit.GetRefTime() ;
}
//if (rawAdcHit.GetNPulses()>0) cout << "Cer npmt = " << " ped = " << rawAdcHit.GetPed() << endl;
for (UInt_t thit = 0; thit < rawAdcHit.GetNPulses(); thit++) {
......@@ -455,6 +469,23 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
if (rawAdcHit.GetPulseAmpRaw(thit) <= 0)
((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 1);
if (rawAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawAdcHit.GetAdcTopC();
Double_t AdcToV = rawAdcHit.GetAdcTomV();
if (fPedDefault[npmt-1] !=0) {
Double_t tPulseInt = AdcToC*(rawAdcHit.GetPulseIntRaw(thit) - fPedDefault[npmt-1]*PeakPedRatio);
((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(npmt, tPulseInt);
((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(npmt, fPedDefault[npmt-1]);
((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(npmt, float(fPedDefault[npmt-1])/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(npmt, 0.);
}
++nrAdcHits;
fTotNumAdcHits++;
fNumAdcHits.at(npmt - 1) = npmt;
......@@ -471,7 +502,9 @@ Int_t THcCherenkov::ApplyCorrections(void) { return (0); }
Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
for(Int_t ipmt = 0; ipmt < fNelem; ipmt++) {
Double_t OffsetTime = 0.0;
if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
for(Int_t ipmt = 0; ipmt < fNelem; ipmt++) {
fAdcPulseAmpTest[ipmt] = -1000.;
fAdcGoodElem[ipmt]=-1;
}
......@@ -480,17 +513,18 @@ Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
Int_t npmt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
Double_t pulseTime = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
Bool_t errorFlag = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin[npmt] && adctdcdiffTime < fAdcTimeWindowMax[npmt];
if (!errorFlag)
{
fGoodAdcMult.at(npmt) += 1;
}
if (!errorFlag && pulseTimeCut && pulseAmp > fAdcPulseAmpTest[npmt]) {
fAdcGoodElem[npmt]=ielem;
fAdcPulseAmpTest[npmt] = pulseAmp;
}
fGoodAdcMult.at(npmt) += 1;
if (!errorFlag) {
if (pulseTimeCut && pulseAmp > fAdcPulseAmpTest[npmt]) {
fAdcGoodElem[npmt]=ielem;
fAdcPulseAmpTest[npmt] = pulseAmp;
}
} else {
if (pulseTimeCut) fAdcGoodElem[npmt]=ielem;
}
}
// Loop over the npmt
for(Int_t npmt = 0; npmt < fNelem; npmt++) {
......@@ -501,17 +535,18 @@ Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
Double_t pulseIntRaw = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
// By default, the last hit within the timing cut will be considered "good"
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
fGoodAdcPed.at(npmt) = pulsePed;
fGoodAdcHitUsed.at(npmt) = ielem + 1;
fGoodAdcPulseInt.at(npmt) = pulseInt;
fGoodAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodAdcPulseAmp.at(npmt) = pulseAmp;
fGoodAdcPulseTime.at(npmt) = pulseTime;
fGoodAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fNpe.at(npmt) = fGain[npmt] * fGoodAdcPulseInt.at(npmt);
fGoodAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fNpe.at(npmt) = fGain[npmt]*pulseInt;
fNpeSum += fNpe.at(npmt);
fTotNumGoodAdcHits++;
......
src/THcCherenkov.h
View file @ e974befb
......@@ -69,6 +69,7 @@ public:
Int_t fTotNumGoodAdcHits;
Int_t fTotNumTracksMatched;
Int_t fTotNumTracksFired;
Double_t fRefTime;
Double_t fNpeSum;
Double_t* fGain;
......@@ -102,6 +103,7 @@ public:
Double_t fNpeThresh;
Double_t* fAdcTimeWindowMin;
Double_t* fAdcTimeWindowMax;
Int_t* fPedDefault;
Double_t fAdcTdcOffset;
Double_t* fRegionValue;
......
src/THcCoinTime.cxx
View file @ e974befb
......@@ -65,6 +65,14 @@ void THcCoinTime::Clear( Option_t* opt )
fROC2_ePosCoinTime=kBig;
fROC1_RAW_CoinTime=kBig;
fROC2_RAW_CoinTime=kBig;
fTRIG1_ePosCoinTime=kBig;
fTRIG4_ePosCoinTime=kBig;
fTRIG1_ePiCoinTime=kBig;
fTRIG4_ePiCoinTime=kBig;
fTRIG1_eKCoinTime=kBig;
fTRIG4_eKCoinTime=kBig;
fTRIG1_epCoinTime=kBig;
fTRIG4_epCoinTime=kBig;
}
//_____________________________________________________________________________
......@@ -263,6 +271,9 @@ Int_t THcCoinTime::Process( const THaEvData& evdata )
Double_t hms_ypfp = theHMSTrack->GetPhi();
Double_t HMS_FPtime = theHMSTrack->GetFPTime();
if (SHMS_FPtime==-2000 || HMS_FPtime==-2000) return 1;
if (SHMS_FPtime==-1000 || HMS_FPtime==-1000) return 1;
//Get raw TDC Times for HMS/SHMS (3/4 trigger)
pTRIG1_TdcTime_ROC1 = fCoinDet->Get_CT_Trigtime(0); //SHMS
pTRIG4_TdcTime_ROC1 = fCoinDet->Get_CT_Trigtime(1); //HMS
......
src/THcDC.h
View file @ e974befb
......@@ -227,7 +227,7 @@ protected:
public:
THcDriftChamberPlane* GetPlane(unsigned int i_plane) {
if(i_plane < fNPlanes) {
if(static_cast<int>(i_plane) < fNPlanes) {
return fPlanes[i_plane];
}
return nullptr;
......
src/THcDriftChamberPlane.cxx
View file @ e974befb
......@@ -298,6 +298,7 @@ Int_t THcDriftChamberPlane::DefineVariables( EMode mode )
{"dist","Drift distancess",
"fHits.THcDCHit.GetDist()"},
{"nhit", "Number of hits", "GetNHits()"},
{"RefTime", "TDC reference time", "fTdcRefTime"},
{ 0 }
};
......@@ -351,7 +352,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fHits->Clear();
fRawHits->Clear();
fTdcRefTime = kBig;
Int_t nrawhits = rawhits->GetLast()+1;
fNRawhits=0;
Int_t ihit = nexthit;
......@@ -365,6 +366,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
Int_t wireNum = hit->fCounter;
THcDCWire* wire = GetWire(wireNum);
Bool_t First_Hit_In_Window = kTRUE;
if (hit->GetRawTdcHit().HasRefTime()) fTdcRefTime = hit->GetRawTdcHit().GetRefTime();
for(UInt_t mhit=0; mhit<hit->GetRawTdcHit().GetNHits(); mhit++) {
fNRawhits++;
/* Sort into early, late and ontime */
......@@ -391,6 +393,7 @@ Int_t THcDriftChamberPlane::SubtractStartTime()
{
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
if (StartTime == -1000) StartTime = 0.0;
for(Int_t ihit=0;ihit<GetNHits();ihit++) {
THcDCHit *thishit = (THcDCHit*) fHits->At(ihit);
Double_t temptime= thishit->GetTime()-StartTime;
......
src/THcDriftChamberPlane.h
View file @ e974befb
......@@ -86,6 +86,8 @@ protected:
TClonesArray* fHits;
TClonesArray* fRawHits;
TClonesArray* fWires;
Double_t fTdcRefTime;
Int_t fVersion;
Int_t fWireOrder;
......
src/THcHallCSpectrometer.cxx
View file @ e974befb
......@@ -142,6 +142,10 @@ Int_t THcHallCSpectrometer::DefineVariables( EMode mode )
fIsSetup = ( mode == kDefine );
RVarDef vars[] = {
{ "tr.betachisq", "Chi2 of beta", "fTracks.THaTrack.GetBetaChi2()"},
{ "tr.GoodPlane4", "Flag for track hitting hodo plane 4", "fTracks.THaTrack.GetGoodPlane4()"},
{ "tr.GoodPlane3", "Flag for track hitting hodo plane 3", "fTracks.THaTrack.GetGoodPlane3()"},
{ "tr.fptime", "Track hodo focal plane time", "fTracks.THaTrack.GetFPTime()"},
{ "tr.npmt", "Track number of hodo PMTs hit", "fTracks.THaTrack.GetNPMT()"},
{ "tr.PruneSelect", "Prune Select ID", "fPruneSelect"},
{ "present", "Trigger Type includes this spectrometer", "fPresent"},
{ 0 }
......
src/THcHelicity.cxx
View file @ e974befb
//*-- Author : Julie Roche, November 2010
// this is a modified version of THaGOHelicity.C
////////////////////////////////////////////////////////////////////////
//
// THcHelicity
//
// Helicity of the beam from QWEAK electronics in delayed mode
// Helicity of the beam in delayed mode using quartets
// +1 = plus, -1 = minus, 0 = unknown
//
// Also supports in-time mode with delay = 0
......@@ -17,24 +15,29 @@
#include "THaApparatus.h"
#include "THaEvData.h"
#include "THcGlobals.h"
#include "THcHelicityScaler.h"
#include "THcParmList.h"
#include "TMath.h"
#include <bitset>
#include <iostream>
#include "THaHelicityDet.h"
#include "hcana/Logger.h"
using namespace std;
//_____________________________________________________________________________
THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus* app)
: hcana::ConfigLogging<THaHelicityDet>(name, description, app), fnQrt(-1), fHelDelay(8),
fMAXBIT(30) {
: THaHelicityDet(name, description, app), fnQrt(-1), fHelDelay(8), fMAXBIT(30) {
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i] = 0;
// memset(fHbits, 0, sizeof(fHbits));
fglHelicityScaler = 0;
fHelicityHistory = 0;
}
//_____________________________________________________________________________
THcHelicity::THcHelicity()
: hcana::ConfigLogging<THaHelicityDet>(), fnQrt(-1), fHelDelay(8), fMAXBIT(30) {
THcHelicity::THcHelicity() : fnQrt(-1), fHelDelay(8), fMAXBIT(30) {
// Default constructor for ROOT I/O
// for( Int_t i = 0; i < NHIST; ++i )
......@@ -72,6 +75,26 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
fEvNumCheck = 0;
fDisabled = kFALSE;
fLastHelpCycle = -1;
fQuadPattern[0] = fQuadPattern[1] = fQuadPattern[2] = fQuadPattern[3] = 0;
fQuadPattern[4] = fQuadPattern[5] = fQuadPattern[6] = fQuadPattern[7] = 0;
fHelperHistory = 0;
fHelperQuartetHistory = 0;
fScalerSeed = 0;
fNBits = 0;
lastispos = 0;
fLastReportedHelicity = kUnknown;
fFixFirstCycle = kFALSE;
fHaveQRT = kFALSE;
fNQRTProblems = 0;
fPeriodCheck = 0.0;
fCycle = 0.0;
fglHelicityScaler = 0;
fHelicityHistory = 0;
fRecommendedFreq = -1.0;
fStatus = kOK;
return fStatus;
}
......@@ -88,7 +111,6 @@ void THcHelicity::Setup(const char* name, const char* description) {
Int_t THcHelicity::ReadDatabase(const TDatime& date) {
_logger->info("In THcHelicity::ReadDatabase");
// cout << "In THcHelicity::ReadDatabase" << endl;
// Read general HelicityDet database values (e.g. fSign)
// Int_t st = THaHelicityDet::ReadDatabase( date );
// if( st != kOK )
......@@ -102,9 +124,12 @@ Int_t THcHelicity::ReadDatabase(const TDatime& date) {
fSign = 1; // Default helicity sign
fRingSeed_reported_initial = 0; // Initial see that should predict reported
// helicity of first quartet.
fFirstCycle = -1; // First Cycle that starts a quad (0 to 3)
fFreq = 29.5596;
fHelDelay = 8;
fFirstCycle = -1; // First Cycle that starts a quad (0 to 3)
fNLastQuartet = -1;
fNQuartet = 0;
// fFreq = 29.5596;
fFreq = 120.0007547169;
fHelDelay = 8;
DBRequest list[] = {// {"_hsign", &fSign, kInt, 0, 1},
{"helicity_delay", &fHelDelay, kInt, 0, 1},
......@@ -139,8 +164,6 @@ Int_t THcHelicity::ReadDatabase(const TDatime& date) {
_logger->info(
"Helicity decoder initialized with frequency of {} Hz and reporting delay of {} cycles.",
fFreq, fHelDelay);
// cout << "Helicity decoder initialized with frequency of " << fFreq
// << " Hz and reporting delay of " << fHelDelay << " cycles." << endl;
return kOK;
}
......@@ -153,8 +176,6 @@ Int_t THcHelicity::DefineVariables(EMode mode) {
// Initialize global variables
_logger->info("Called THcHelicity::DefineVariables with mode == {}", mode);
// cout << "Called THcHelicity::DefineVariables with mode == "
// << mode << endl;
if (mode == kDefine && fIsSetup)
return kOK;
......@@ -168,8 +189,10 @@ Int_t THcHelicity::DefineVariables(EMode mode) {
{"helrep", "reported helicity for event", "fReportedHelicity"},
{"helpred", "predicted reported helicity for event", "fPredictedHelicity"},
{"mps", "In MPS blanking period", "fMPS"},
{"pcheck", "Period check", "fPeriodCheck"},
{"cycle", "Helicity Cycle", "fCycle"},
{"qrt", "Last cycle of quartet", "fQrt"},
{0}};
// cout << "Calling THcHelicity DefineVarsFromList" << endl;
_logger->info("Calling THcHelicity DefineVarsFromList");
return DefineVarsFromList(var, mode);
}
......@@ -177,9 +200,9 @@ Int_t THcHelicity::DefineVariables(EMode mode) {
void THcHelicity::PrintEvent(Int_t evtnum) {
cout << " ++++++ THcHelicity::Print ++++++\n";
_logger->info(" ++++++ THcHelicity::Print ++++++");
cout << " +++++++++++++++++++++++++++++++++++++\n";
_logger->info(" +++++++++++++++++++++++++++++++++++++\n");
return;
}
......@@ -231,21 +254,29 @@ void THcHelicity::Clear(Option_t* opt) {
return;
}
//_____________________________________________________________________________
void THcHelicity::SetHelicityScaler(THcHelicityScaler* f) {
fglHelicityScaler = f;
fHelicityHistory = fglHelicityScaler->GetHelicityHistoryP();
}
//_____________________________________________________________________________
Int_t THcHelicity::Decode(const THaEvData& evdata) {
// Decode Helicity data.
// Return 1 if helicity was assigned, 0 if not, <0 if error.
evnum = evdata.GetEvNum();
Int_t err = ReadData(evdata); // from THcHelicityReader class
if (err) {
Error(Here("THcHelicity::Decode"), "Error decoding helicity data.");
_logger->error("THcHelicity::Decode : Error decoding helicity data.");
return err;
}
fReportedHelicity = (fIsHelp ? (fIsHelm ? kUnknown : kPlus) : (fIsHelm ? kMinus : kUnknown));
fMPS = fIsMPS ? 1 : 0;
fQrt = fIsQrt ? 1 : 0; // Last of quartet
if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
return 0;
......@@ -276,9 +307,17 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
Int_t missed = 0;
// Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0;
if (fIsMPS) {
fPeriodCheck = fmod(fTITime / (250000000.0 / fFreq) - fPeriodCheckOffset, 1.0);
fCycle = (fTITime / fTIPeriod);
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
if (fFoundMPS) {
if (fRecommendedFreq < 0.0) {
if (TMath::Abs(fPeriodCheck - 0.5) > 0.25) {
fRecommendedFreq = fFreq * (1 - (fPeriodCheck - 0.5) / fCycle);
}
}
missed = TMath::Nint(fTITime / fTIPeriod - fLastMPSTime / fTIPeriod);
if (missed < 1) { // was <=1
fLastMPSTime = (fTITime + fLastMPSTime + missed * fTIPeriod) / 2;
......@@ -293,8 +332,9 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
// cout << "Found MPS" << endl;
// check for Nint((time-last)/period) > 1
} else {
fFoundMPS = kTRUE;
fLastMPSTime = fTITime;
fFoundMPS = kTRUE;
fLastMPSTime = fTITime;
fPeriodCheckOffset = (fPeriodCheck - .5);
}
} else if (fFoundMPS) { //
if (fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods
......@@ -303,20 +343,23 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
// cout << "Missed " << missed << " MPSes" << endl;
Int_t newNCycle = fNCycle + missed - 1; // How many cycles really missed
Int_t quartets_missed = (newNCycle - fFirstCycle) / 4 - (fNCycle - fFirstCycle) / 4;
for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
int quartetphase = (newNCycle - fFirstCycle) % 4;
int quartetphase = (newNCycle - fFirstCycle) % 4;
// cout << " " << fNCycle << " " << newNCycle << " " << fFirstCycle << " " <<
// quartets_missed << " " << quartetphase << endl; cout << "Cycles " << fNCycle << "
// " << newNCycle << " " << fFirstCycle
// quartets_missed << " " << quartetphase << endl;
// cout << "Cycles " << fNCycle << " " << newNCycle << " " << fFirstCycle
// << " skipped " << quartets_missed << " quartets" << endl;
fNCycle = newNCycle;
// Need to reset fQuartet to reflect where we are based on the current
// reported helicity. So we don't fail quartet testing.
// But only do this if we are calibrated.
if (fNBits >= fMAXBIT) {
if (fNBits >= fMAXBIT + 0) {
for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
// cout << "Advancing seed A " << fNBits << endl;
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
fActualHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity
......@@ -342,7 +385,7 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
fIsNewCycle = kTRUE;
fLastReportedHelicity = fReportedHelicity;
} else { // No missed periods. Get helicities from rings
if (fNBits >= fMAXBIT) {
if (fNBits >= fMAXBIT + 0) {
int quartetphase = (fNCycle - fFirstCycle) % 4;
fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
......@@ -356,49 +399,98 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
}
}
if (fIsNewCycle) {
// cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
// Int_t predictedScalerSeed = RanBit30(fScalerSeed);
// cout << "Predct Seed " << bitset<32>(predictedScalerSeed) << endl;
// cout << "Event Seed " << bitset<32>(fRingSeed_reported) << endl;
fQuartet[3] = fQuartet[2];
fQuartet[2] = fQuartet[1];
fQuartet[1] = fQuartet[0];
fQuartet[0] = fReportedHelicity;
fNCycle++;
if ((fNCycle - fFirstCycle) % 4 == 3) { // Test if last in a quartet
if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
if (!fFoundQuartet) {
// fFirstCycle = fNCycle - 3;
_logger->info("Quartet potentially found, starting at cycle {} - event {}",
fFirstCycle, evdata.GetEvNum());
// cout << "Quartet potentially found, starting at cycle " << fFirstCycle << " - event
// "
// << evdata.GetEvNum() << endl;
fFoundQuartet = kTRUE;
}
} else {
if (fNCycle - fFirstCycle > 4) { // Not at start of run. Reset
_logger->warn("Lost quartet sync at cycle {} - event {}", fNCycle, evdata.GetEvNum());
_logger->warn("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
// cout << "Lost quartet sync at cycle " << fNCycle << " - event " <<
// evdata.GetEvNum()
// << endl;
// cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " <<
// fQuartet[3]
// << endl;
fFirstCycle += 4 * ((fNCycle - fFirstCycle) / 4); // Update, but don't change phase
// cout << "XX: " << fNCycle << " " << fReportedHelicity << " " << fIsQrt <<endl;
if (fIsQrt) { //
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
if (fHaveQRT) { // Should already have phase set
if ((fNCycle - fFirstCycle) % 4 != 0) { // Test if first in a quartet
fNQRTProblems++;
if (fNQRTProblems > 10) {
_logger->warn("QRT Problem resetting");
fHaveQRT = kFALSE;
fFoundQuartet = kFALSE;
fNQRTProblems = 0;
} else {
_logger->warn("Ignored {} problems", fNQRTProblems);
}
} else {
fNQRTProblems = 0;
}
fFoundQuartet = kFALSE;
}
if (!fHaveQRT) {
fHaveQRT = kTRUE;
fFoundQuartet = kTRUE;
fFirstCycle = fNCycle; //
fNQuartet = (fNCycle - fFirstCycle) / 4;
fNLastQuartet = fNQuartet - 1; // Make sure LoadHelicity uses
fNBits = 0;
_logger->info("Searching for first of a quartet at cycle {} - event {}", fFirstCycle,
evdata.GetEvNum());
// cout << "Searching for first of a quartet at cycle "
// << " " << fFirstCycle << " - event " << evdata.GetEvNum() << endl;
// cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " << fQuartet[3]
// << endl;
_logger->info("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
fFirstCycle++;
fNQRTProblems = 0;
_logger->info("Phase found from QRT signal");
_logger->info("fFirstCycle = {}", fFirstCycle);
}
} else {
if (fHaveQRT) { // Using qrt signal.
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
if ((fNCycle - fFirstCycle) % 4 == 0) { // Shouldn't happen
fNQRTProblems++;
if (fNQRTProblems > 10) {
_logger->warn("Shouldn't happen, cycle= {} / {}", fNCycle, fFirstCycle);
fHaveQRT = kFALSE; // False until a new QRT seen
fNBits = 0; // Reset
fNLastQuartet = fNQuartet; // Make sure LoadHelicity does not use
} else {
_logger->warn("Ignored {} problems", fNQRTProblems);
}
}
} else { // Presumable pre qrt signal data
if ((fNCycle - fFirstCycle) % 4 == 3) { // Test if last in a quartet
if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
if (!fFoundQuartet) {
// fFirstCycle = fNCycle - 3;
_logger->warn("Quartet potentially found, starting at cycle {}", fFirstCycle);
fNQuartet = (fNCycle - fFirstCycle) / 4;
fNLastQuartet = fNQuartet - 1; // Make sure LoadHelicity uses
fFoundQuartet = kTRUE;
}
} else {
if (fNCycle - fFirstCycle > 4) { // Not at start of run. Reset
_logger->warn("Lost quartet sync at cycle {}: {} {} {} {}", fNCycle, fQuartet[0],
fQuartet[1], fQuartet[2], fQuartet[3]);
fFirstCycle +=
4 * ((fNCycle - fFirstCycle) / 4); // Update, but don't change phase
}
fFoundQuartet = kFALSE;
fNBits = 0;
_logger->info("Searching for first of a quartet at cycle {}: {} {} {} {}",
fFirstCycle, fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
fFirstCycle++;
}
} else {
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
}
}
}
// Load the actual helicity. Calibrate if not calibrated.
fActualHelicity = kUnknown;
// Here if we know we missed some earlier in the quartet, we need
// to make sure we get here and call LoadHelicity for the missing
// Cycles, reducing the missed count for each extra loadhelicity
// call that we make.
LoadHelicity(fReportedHelicity, fNCycle, missed);
fLastReportedHelicity = fReportedHelicity;
fIsNewCycle = kFALSE;
......@@ -407,11 +499,11 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
// << fPredictedHelicity << " " << fActualHelicity << endl;
}
// Ignore until a MPS Is found
} else { // No MPS found yet
fActualHelicity = kUnknown;
}
} else {
// cout << "Initializing" << endl;
_logger->info("Initializing Helicity");
fLastReportedHelicity = fReportedHelicity;
fActualHelicity = kUnknown;
......@@ -429,18 +521,16 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
// Some sanity checks
if (fActualHelicity < -5) {
_logger->info("Actual Helicity never got defined");
_logger->warn("Actual Helicity never got defined");
}
if (fNBits < fMAXBIT) {
if (fActualHelicity == -1 || fActualHelicity == 1) {
cout << "Helicity of " << fActualHelicity << " reported prematurely at cycle " << fNCycle
<< endl;
_logger->warn("Helicity of {} reported prematurely at cycle {}", fActualHelicity, fNCycle);
}
}
fLastActualHelicity = fActualHelicity;
return 0;
}
//_____________________________________________________________________________
Int_t THcHelicity::End(THaRunBase*) {
// End of run processing. Write histograms.
......@@ -449,10 +539,26 @@ Int_t THcHelicity::End(THaRunBase*) {
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i]->Write();
if (fRecommendedFreq < 0.0) {
fRecommendedFreq = fFreq * (1 - (fPeriodCheck - 0.5) / fCycle);
}
if (TMath::Abs(1 - fRecommendedFreq / fFreq) >= 0.5e-6) {
cout << "------------- HELICITY DECODING ----------------------" << endl;
cout << "Actual helicity reversal frequency differs from \"helicity_freq\" value" << endl;
cout << "If there are helicity decoding errors beyond the start of the run, " << endl;
streamsize ss = cout.precision();
cout.precision(10);
cout << "try replacing helicity_freq value of " << fFreq << " with " << fRecommendedFreq
<< endl;
cout << "If that still gives helicity errors, try " << 0.9999999 * fRecommendedFreq << endl;
cout.precision(ss);
cout << "------------------------------------------------------" << endl;
}
return 0;
}
//_____________________________________________________________________________
//_____________________________________________________________________________
void THcHelicity::SetDebug(Int_t level) {
// Set debug level of this detector as well as the THcHelicityReader
// helper class.
......@@ -460,61 +566,94 @@ void THcHelicity::SetDebug(Int_t level) {
THaHelicityDet::SetDebug(level);
fQWEAKDebug = level;
}
//_____________________________________________________________________________
//_____________________________________________________________________________
void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles) {
// static const char* const here = "THcHelicity::LoadHelicity";
int quartetphase = (cyclecount - fFirstCycle) % 4;
fnQrt = quartetphase;
if (missedcycles > 1) { // If we missed windows
if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over
// if(!fFixFirstCycle) {
if (fNQuartet - fNLastQuartet > 1) { // If we missed a quartet
if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over
_logger->warn("fNBits = 0, missedcycles={} {} {}", missedcycles, fNLastQuartet,
fNQuartet);
fNBits = 0;
return;
}
}
// }
if (!fFoundQuartet) { // Wait until we have found quad phase before starting
return; // to calibrate
}
if (quartetphase == 0) { // Start of a quad
if (fNBits < fMAXBIT) {
// LoadHelicity is called first event of each cycle.
// But only do it's thing if it is first cycle of quartet.
// Need to instead have it do it's thing on the first event of a quartet.
// But only for seed building. Current logic is OK once seed is found.
if (fNBits < fMAXBIT) {
if (fNQuartet != fNLastQuartet) {
// Sanity check fNQuartet == fNLastQuartet+1
if (fNBits == 0) {
_logger->info("Start calibrating at cycle {}", cyclecount);
// cout << "Start calibrating at cycle " << cyclecount << endl;
fRingSeed_reported = 0;
}
if (fReportedHelicity == kPlus) {
// Do phase stuff right here
if ((fReportedHelicity == kPlus && (quartetphase == 0 || quartetphase == 3)) ||
(fReportedHelicity == kMinus && (quartetphase == 1 || quartetphase == 2))) {
fRingSeed_reported = ((fRingSeed_reported << 1) | 1) & 0x3FFFFFFF;
// cout << " " << fNQuartet << " 1" << endl;
} else {
fRingSeed_reported = (fRingSeed_reported << 1) & 0x3FFFFFFF;
// cout << " " << fNQuartet << " 0" << endl;
}
fNBits++;
if (fReportedHelicity == kUnknown) {
fNBits = 0;
fRingSeed_reported = 0;
} else if (fNBits == fMAXBIT) {
_logger->info("Seed Found {} at cycle {} with first cycle {}", fRingSeed_reported,
} else if (fNBits == fMAXBIT + 0) {
_logger->info("Seed Found {:32b} at cycle {} with first cycle {}", fRingSeed_reported,
cyclecount, fFirstCycle);
// cout << "Seed Found " << hex << fRingSeed_reported << dec << " at cycle " << cyclecount
// << " with first cycle " << fFirstCycle << endl;
if (fglHelicityScaler) {
_logger->info("Scaler Seed {:32b}", fScalerSeed);
}
Int_t backseed = GetSeed30(fRingSeed_reported);
_logger->info("Seed at cycle {} should be {}", fFirstCycle, backseed);
// cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec <<
// endl;
}
fActualHelicity = kUnknown;
} else if (fNBits >= fMAXBIT) {
fRingSeed_reported = RanBit30(fRingSeed_reported);
if (fNBits == fMAXBIT) {
_logger->info("Seed at cycle {} should be {:#x}", fFirstCycle, backseed);
// Create the "actual seed"
fRingSeed_actual = fRingSeed_reported;
for (Int_t i = 0; i < fHelDelay / 4; i++) {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fNBits++;
} else {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
} else if (fglHelicityScaler && fNBits > 2) { // Try the scalers
if (fglHelicityScaler->IsSeedGood()) {
Int_t scalerseed = fglHelicityScaler->GetReportedSeed();
fRingSeed_reported = RanBit30(scalerseed);
_logger->info(" -- Getting seed from scalers -- ");
cout << " Seed Found " << bitset<32>(fRingSeed_reported) << " at cycle " << cyclecount
<< " with first cycle " << fFirstCycle << endl;
cout << "Scaler Seed " << bitset<32>(scalerseed) << endl;
// Create the "actual seed"
fRingSeed_actual = fRingSeed_reported;
for (Int_t i = 0; i < fHelDelay / 4; i++) {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
fNBits = fMAXBIT + 0;
}
}
fActualHelicity = kUnknown;
} // Need to change this to build seed even when not at start of quartet
} else {
if (quartetphase == 0) {
// If quartetphase !=, the seeds will alread have been advanced
// except that we won't have made the initial
// cout << "Advancing seed B " << fNBits << endl;
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
fActualHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
// if(fTITime/250000000.0 > 380.0) cout << fTITime/250000000.0 << " " << fNCycle << " "
......@@ -524,24 +663,18 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
if (fReportedHelicity != fPredictedHelicity) {
_logger->warn("Helicity prediction failed {} {} {}", fReportedHelicity, fPredictedHelicity,
fActualHelicity);
// cout << "Helicity prediction failed " << fReportedHelicity << " "
// << fPredictedHelicity << " " << fActualHelicity << endl;
// cout << hex << fRingSeed_reported << " " << fRingSeed_actual << dec << endl;
_logger->warn("{} {}", fRingSeed_reported, fRingSeed_actual);
fNBits = 0; // Need to reaquire seed
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
}
fQuartetStartHelicity = fActualHelicity;
fQuartetStartPredictedHelicity = fPredictedHelicity;
}
fQuartetStartHelicity = fActualHelicity;
fQuartetStartPredictedHelicity = fPredictedHelicity;
} else { // Not the beginning of a quad
if (fNBits >= fMAXBIT) {
fActualHelicity =
(quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity : -fQuartetStartHelicity;
fPredictedHelicity = (quartetphase == 0 || quartetphase == 3)
? fQuartetStartPredictedHelicity
: -fQuartetStartPredictedHelicity;
}
fActualHelicity =
(quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity : -fQuartetStartHelicity;
fPredictedHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartPredictedHelicity
: -fQuartetStartPredictedHelicity;
}
return;
}
......@@ -565,7 +698,7 @@ Int_t THcHelicity::RanBit30(Int_t ranseed) {
ranseed = ((ranseed << 1) | newbit) & 0x3FFFFFFF;
// here ranseed is changed
if (fQWEAKDebug > 1) {
cout << "THcHelicity::RanBit30, newbit=" << newbit << "\n";
_logger->info("THcHelicity::RanBit30, newbit={}", newbit);
}
return ranseed;
......
src/THcHelicity.h
View file @ e974befb
......@@ -6,107 +6,124 @@
// THcHelicity
//
// Helicity of the beam - from QWEAK electronics in delayed mode
//
//
////////////////////////////////////////////////////////////////////////
#include "hcana/Logger.h"
#include "THaHelicityDet.h"
#include "THcHelicityReader.h"
#include "hcana/Logger.h"
class TH1F;
class THcHelicityScaler;
class THcHelicity : public hcana::ConfigLogging<THaHelicityDet>, public THcHelicityReader {
class THcHelicity : public THaHelicityDet, public THcHelicityReader {
public:
THcHelicity( const char* name, const char* description,
THaApparatus* a = NULL );
THcHelicity(const char* name, const char* description, THaApparatus* a = NULL);
THcHelicity();
virtual ~THcHelicity();
virtual EStatus Init(const TDatime& date);
virtual void MakePrefix();
virtual void MakePrefix();
virtual Int_t Begin( THaRunBase* r=0 );
virtual void Clear( Option_t* opt = "" );
virtual Int_t Decode( const THaEvData& evdata );
virtual Int_t End( THaRunBase* r=0 );
virtual void SetDebug( Int_t level );
virtual Bool_t HelicityValid() const { return fValidHel; }
virtual Int_t Begin(THaRunBase* r = 0);
virtual void Clear(Option_t* opt = "");
virtual Int_t Decode(const THaEvData& evdata);
virtual Int_t End(THaRunBase* r = 0);
virtual void SetDebug(Int_t level);
virtual void SetHelicityScaler(THcHelicityScaler* f);
void PrintEvent(Int_t evtnum);
protected:
void Setup(const char* name, const char* description);
void Setup(const char* name, const char* description);
std::string fKwPrefix;
void FillHisto();
void LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles);
Int_t RanBit30(Int_t ranseed);
Int_t GetSeed30(Int_t currentseed);
// Fixed Parameters
Int_t fRingSeed_reported_initial;
Int_t fFirstCycle;
Int_t fRingSeed_reported_initial;
Int_t fFirstCycle;
Bool_t fFixFirstCycle;
Double_t fFreq;
Double_t fRecommendedFreq;
Double_t fTIPeriod; // Reversal period in TI time units
Double_t fTIPeriod; // Reversal period in TI time units
Double_t fPeriodCheck;
Double_t fPeriodCheckOffset;
Double_t fCycle;
Bool_t fFirstEvProcessed;
Int_t fLastReportedHelicity;
Bool_t fFirstEvProcessed;
Int_t fLastReportedHelicity;
Long64_t fFirstEvTime;
Long64_t fLastEvTime;
Long64_t fLastMPSTime;
Int_t fReportedHelicity;
Int_t fMPS;
Int_t fPredictedHelicity;
Int_t fActualHelicity;
Int_t fQuartetStartHelicity;
Int_t fQuartetStartPredictedHelicity;
Bool_t fFoundMPS;
Bool_t fFoundQuartet; // True if quartet phase probably found.
Bool_t fIsNewCycle;
Int_t fNCycle; // Count of # of helicity cycles
Int_t fQuartet[4]; // For finding and checking quartet pattern
Int_t fNBits;
Int_t fnQrt; // Position in quartet
Int_t fReportedHelicity;
Int_t fMPS;
Int_t fPredictedHelicity;
Int_t fActualHelicity;
Int_t fQuartetStartHelicity;
Int_t fQuartetStartPredictedHelicity;
Bool_t fFoundMPS;
Bool_t fFoundQuartet; // True if quartet phase probably found.
Bool_t fIsNewCycle;
Int_t fNCycle; // Count of # of helicity cycles
Int_t fNQuartet; // Quartet count
Int_t fNLastQuartet;
Int_t fQuartet[4]; // For finding and checking quartet pattern
Int_t fNBits;
Int_t fnQrt; // Position in quartet
Bool_t fHaveQRT; // QRT signal exists
Int_t fNQRTProblems;
Int_t fRingSeed_reported;
Int_t fRingSeed_actual;
// Offset between the ring reported value and the reported value
Int_t fHelDelay;
Int_t fHelDelay;
// delay of helicity (# windows)
Int_t fMAXBIT;
//number of bit in the pseudo random helcity generator
Int_t fMAXBIT;
// number of bit in the pseudo random helcity generator
std::vector<Int_t> fPatternSequence; // sequence of 0 and 1 in the pattern
Int_t fQWEAKNPattern; // maximum of event in the pattern
Bool_t HWPIN;
Int_t fQWEAKNPattern; // maximum of event in the pattern
Bool_t HWPIN;
Int_t fQrt;
Int_t fTSettle;
Bool_t fValidHel;
Int_t fHelicityLastTIR;
Int_t fPatternLastTIR;
void SetErrorCode(Int_t error);
void SetErrorCode(Int_t error);
Double_t fErrorCode;
Int_t fEvtype; // Current CODA event type
Int_t fLastActualHelicity;
Int_t fEvNumCheck;
Bool_t fDisabled;
static const Int_t NHIST = 2;
TH1F* fHisto[NHIST];
virtual Int_t DefineVariables( EMode mode = kDefine );
virtual Int_t ReadDatabase( const TDatime& date );
ClassDef(THcHelicity,0) // Beam helicity from QWEAK electronics in delayed mode
Int_t fEvtype; // Current CODA event type
Int_t fLastActualHelicity;
Int_t fEvNumCheck;
Bool_t fDisabled;
static const Int_t NHIST = 2;
TH1F* fHisto[NHIST];
virtual Int_t DefineVariables(EMode mode = kDefine);
virtual Int_t ReadDatabase(const TDatime& date);
THcHelicityScaler* fglHelicityScaler = nullptr;
Int_t* fHelicityHistory = nullptr;
Int_t fLastHelpCycle;
Int_t fScaleQuartet;
Int_t fQuadPattern[8];
Int_t fHelperHistory;
Int_t fHelperQuartetHistory;
Int_t fScalerSeed;
Int_t lastispos;
Int_t evnum;
Int_t fThisScaleHel;
Int_t fLastScaleHel;
Int_t fLastLastScaleHel;
ClassDef(THcHelicity, 0) // Beam helicity from QWEAK electronics in delayed mode
};
#endif
......
src/THcHelicityReader.cxx
View file @ e974befb
......@@ -7,30 +7,28 @@
////////////////////////////////////////////////////////////////////////
#include "THcHelicityReader.h"
#include "TError.h"
#include "TH1F.h"
#include "THaAnalysisObject.h" // For LoadDB
#include "THaEvData.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "TMath.h"
#include "TError.h"
#include "VarDef.h"
#include "THaAnalysisObject.h" // For LoadDB
#include "hcana/Logger.h"
#include <iostream>
#include <vector>
#include "TH1F.h"
using namespace std;
//____________________________________________________________________
THcHelicityReader::THcHelicityReader()
: fTITime(0), fTITime_last(0), fTITime_rollovers(0),
fHaveROCs(kFALSE)
{
: hcana::ConfigLogging<podd2::EmptyBase>(), fTITime(0), fTITime_last(0), fTITime_rollovers(0),
fHaveROCs(kFALSE) {
// Default constructor
}
//____________________________________________________________________
THcHelicityReader::~THcHelicityReader()
{
THcHelicityReader::~THcHelicityReader() {
// Destructor
// Histograms will be deleted by ROOT
......@@ -40,34 +38,26 @@ THcHelicityReader::~THcHelicityReader()
}
//____________________________________________________________________
void THcHelicityReader::Clear( Option_t* )
{
fIsMPS = fIsQrt = fIsHelp = fIsHelm = kFALSE;
}
void THcHelicityReader::Clear(Option_t*) { fIsMPS = fIsQrt = fIsHelp = fIsHelm = kFALSE; }
//_____________________________________________________________________________
Int_t THcHelicityReader::ReadDatabase( const char* /*dbfilename*/,
const char* /*prefix*/,
const TDatime& /*date*/,
int /*debug_flag*/ )
{
Int_t THcHelicityReader::ReadDatabase(const char* /*dbfilename*/, const char* /*prefix*/,
const TDatime& /*date*/, int /*debug_flag*/) {
// Eventually get these from the parameter file
const static char* const here = "THcHelicityReader::ReadDatabase";
// SHMS settings see https://logbooks.jlab.org/entry/3614445
_logger->info("{}: Helicity information from ROC 2 (SHMS)", here);
SetROCinfo(kHel, 2, 14, 9);
SetROCinfo(kHelm, 2, 14, 8);
SetROCinfo(kMPS, 2, 14, 10);
SetROCinfo(kQrt, 2, 14, 7); // Starting about run 5818
SetROCinfo(kTime, 2, 21, 2);
SetROCinfo(kHel,2,14,9);
SetROCinfo(kHelm,2,14,8);
SetROCinfo(kMPS,2,14,10);
SetROCinfo(kQrt,2,14,7);
SetROCinfo(kTime,2,21,2);
fADCThreshold = 8000;
DBRequest list[] = {
{"helicity_adcthreshold",&fADCThreshold, kInt, 0, 1},
{0}
};
DBRequest list[] = {{"helicity_adcthreshold", &fADCThreshold, kInt, 0, 1}, {0}};
gHcParms->LoadParmValues(list, "");
......@@ -75,28 +65,25 @@ Int_t THcHelicityReader::ReadDatabase( const char* /*dbfilename*/,
}
//_____________________________________________________________________________
void THcHelicityReader::Begin()
{
void THcHelicityReader::Begin() {
// static const char* const here = "THcHelicityReader::Begin";
// cout<<here<<endl;
fTITime_last = 0;
fTITime = 0;
fTITime_last = 0;
fTITime = 0;
fTITime_rollovers = 0;
return;
}
//____________________________________________________________________
void THcHelicityReader::End()
{
void THcHelicityReader::End() {
// static const char* const here = "THcHelicityReader::End";
// cout<<here<<endl;
return;
}
//____________________________________________________________________
Int_t THcHelicityReader::ReadData( const THaEvData& evdata )
{
Int_t THcHelicityReader::ReadData(const THaEvData& evdata) {
// Obtain the present data from the event for QWEAK helicity mode.
static const char* here = "THcHelicityReader::ReadData";
......@@ -107,27 +94,48 @@ Int_t THcHelicityReader::ReadData( const THaEvData& evdata )
// std::cout<<" which="<<jk
// <<" roc="<<fROCinfo[jk].roc
// <<" header="<<fROCinfo[jk].header
// <<" index="<<fROCinfo[jk].index
// <<" index="<<fROCinfo[jk].index
// <<endl;
// }
// std::cout<<" fHaveROCs="<<fHaveROCs<<endl;
if( !fHaveROCs ) {
::Error( here, "ROC data (detector map) not properly set up." );
if (!fHaveROCs) {
::Error(here, "ROC data (detector map) not properly set up.");
return -1;
}
// Check if ROC info is correct
if (!evdata.GetModule(fROCinfo[kTime].roc, fROCinfo[kTime].slot)) {
_logger->info("{}: ROC 2 not found, changing to ROC 1 (HMS)", here);
SetROCinfo(kHel, 1, 18, 9);
SetROCinfo(kHelm, 1, 18, 8);
SetROCinfo(kMPS, 1, 18, 10);
SetROCinfo(kQrt, 1, 18, 7);
SetROCinfo(kTime, 1, 21, 2);
}
// Get the TI Data
// Int_t fTIType = evData.GetData(fTICrate, fTISlot, 0, 0);
// Int_t fTIEvNum = evData.GetData(fTICrate, fTISlot, 1, 0);
UInt_t titime = (UInt_t) evdata.GetData(fROCinfo[kTime].roc,
fROCinfo[kTime].slot,
fROCinfo[kTime].index, 0);
//cout << fTITime_last << " " << titime << endl;
if(titime < fTITime_last) {
UInt_t titime =
(UInt_t)evdata.GetData(fROCinfo[kTime].roc, fROCinfo[kTime].slot, fROCinfo[kTime].index, 0);
// Check again if ROC info is correct
if (titime == 0 && fTITime_last == 0) {
_logger->info("{}: ROC 2 not found, changing to ROC 1 (HMS)", here);
SetROCinfo(kHel, 1, 18, 9);
SetROCinfo(kHelm, 1, 18, 8);
SetROCinfo(kMPS, 1, 18, 10);
SetROCinfo(kQrt, 1, 18, 7);
SetROCinfo(kTime, 1, 21, 2);
titime =
(UInt_t)evdata.GetData(fROCinfo[kTime].roc, fROCinfo[kTime].slot, fROCinfo[kTime].index, 0);
}
// cout << fTITime_last << " " << titime << endl;
if (titime < fTITime_last) {
fTITime_rollovers++;
}
fTITime = titime + fTITime_rollovers*4294967296;
fTITime = titime + fTITime_rollovers * 4294967296;
fTITime_last = titime;
const_cast<THaEvData&>(evdata).SetEvTime(fTITime);
......@@ -135,36 +143,26 @@ Int_t THcHelicityReader::ReadData( const THaEvData& evdata )
// Get the helicity control signals. These are from the pedestals
// acquired by FADC channels.
Int_t helpraw = evdata.GetData(Decoder::kPulsePedestal,
fROCinfo[kHel].roc,
fROCinfo[kHel].slot,
fROCinfo[kHel].index, 0);
Int_t helmraw = evdata.GetData(Decoder::kPulsePedestal,
fROCinfo[kHelm].roc,
fROCinfo[kHelm].slot,
fROCinfo[kHelm].index, 0);
Int_t mpsraw = evdata.GetData(Decoder::kPulsePedestal,
fROCinfo[kMPS].roc,
fROCinfo[kMPS].slot,
fROCinfo[kMPS].index, 0);
Int_t qrtraw = evdata.GetData(Decoder::kPulsePedestal,
fROCinfo[kQrt].roc,
fROCinfo[kQrt].slot,
fROCinfo[kQrt].index, 0);
fIsQrt = qrtraw > fADCThreshold;
fIsMPS = mpsraw > fADCThreshold;
Int_t helpraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kHel].roc, fROCinfo[kHel].slot,
fROCinfo[kHel].index, 0);
Int_t helmraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kHelm].roc, fROCinfo[kHelm].slot,
fROCinfo[kHelm].index, 0);
Int_t mpsraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kMPS].roc, fROCinfo[kMPS].slot,
fROCinfo[kMPS].index, 0);
Int_t qrtraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kQrt].roc, fROCinfo[kQrt].slot,
fROCinfo[kQrt].index, 0);
fIsQrt = qrtraw > fADCThreshold;
fIsMPS = mpsraw > fADCThreshold;
fIsHelp = helpraw > fADCThreshold;
fIsHelm = helmraw > fADCThreshold;
return 0;
}
//TODO: this should not be needed once LoadDB can fill fROCinfo directly
// TODO: this should not be needed once LoadDB can fill fROCinfo directly
//____________________________________________________________________
Int_t THcHelicityReader::SetROCinfo( EROC which, Int_t roc,
Int_t slot, Int_t index )
{
Int_t THcHelicityReader::SetROCinfo(EROC which, Int_t roc, Int_t slot, Int_t index) {
// Define source and offset of data. Normally called by ReadDatabase
// of the detector that is a THcHelicityReader.
......@@ -173,18 +171,19 @@ Int_t THcHelicityReader::SetROCinfo( EROC which, Int_t roc,
// You must define at least the kHel and kTime ROCs.
// Returns <0 if parameter error, 0 if success
if( which<kHel || which>=kCount )
if (which < kHel || which >= kCount)
return -1;
if( roc <= 0 || roc > 255 )
if (roc <= 0 || roc > 255)
return -2;
fROCinfo[which].roc = roc;
fROCinfo[which].slot = slot;
fROCinfo[which].index = index;
fROCinfo[which].roc = roc;
fROCinfo[which].slot = slot;
fROCinfo[which].index = index;
_logger->info("SetROCInfo: {} (roc: {}, slot: {}, index: {})", which, fROCinfo[which].roc,
fROCinfo[which].slot, fROCinfo[which].index);
fHaveROCs = (fROCinfo[kHel].roc > 0 && fROCinfo[kTime].roc > 0 && fROCinfo[kMPS].roc);
//cout << "SetROCInfo: " << which << " " << fROCinfo[kHel].roc << " " << fROCinfo[kTime].roc << endl;
fHaveROCs = ( fROCinfo[kHel].roc > 0 && fROCinfo[kTime].roc > 0 );
return 0;
}
......
src/THcHelicityReader.h
View file @ e974befb
......@@ -10,60 +10,59 @@
//////////////////////////////////////////////////////////////////////////
#include "Rtypes.h"
#include "THaHelicityDet.h"
#include "hcana/Logger.h"
class THaEvData;
class TDatime;
class TH1F;
class THcHelicityReader {
class THcHelicityReader : public hcana::ConfigLogging<podd2::EmptyBase> {
public:
THcHelicityReader();
virtual ~THcHelicityReader();
struct ROCinfo {
Int_t roc; // ROC to read out
Int_t slot; // Headers to search for (0 = ignore)
Int_t index; // Index into buffer
Int_t roc; // ROC to read out
Int_t slot; // Headers to search for (0 = ignore)
Int_t index; // Index into buffer
};
protected:
protected:
// Used by ReadDatabase
enum EROC { kHel = 0, kHelm, kMPS, kQrt, kTime, kCount };
Int_t SetROCinfo( EROC which, Int_t roc, Int_t slot, Int_t index );
Int_t SetROCinfo(EROC which, Int_t roc, Int_t slot, Int_t index);
virtual void Clear( Option_t* opt="" );
virtual Int_t ReadData( const THaEvData& evdata );
Int_t ReadDatabase( const char* dbfilename, const char* prefix,
const TDatime& date, int debug_flag = 0 );
void Begin();
void End();
virtual void Clear(Option_t* opt = "");
virtual Int_t ReadData(const THaEvData& evdata);
Int_t ReadDatabase(const char* dbfilename, const char* prefix, const TDatime& date,
int debug_flag = 0);
void Begin();
void End();
ULong64_t fTITime;
UInt_t fTITime_last;
UInt_t fTITime_rollovers;
UInt_t fTITime_last;
UInt_t fTITime_rollovers;
// Reported Helicity status for the event
Bool_t fIsMPS;
Bool_t fIsQrt;
Bool_t fIsHelp;
Bool_t fIsHelm;
Int_t fADCThreshold; // Threshold for On/Off of helicity signals
Int_t fADCThreshold; // Threshold for On/Off of helicity signals
ROCinfo fROCinfo[kCount];
ROCinfo fROCinfo[kCount];
Int_t fQWEAKDebug; // Debug level
Bool_t fHaveROCs; // Required ROCs are defined
Bool_t fNegGate; // Invert polarity of gate, TO DO implement this functionality
Int_t fQWEAKDebug; // Debug level
Bool_t fHaveROCs; // Required ROCs are defined
Bool_t fNegGate; // Invert polarity of gate, TO DO implement this functionality
static const Int_t NHISTR = 12;
// TH1F* fHistoR[12]; // Histograms
private:
ClassDef(THcHelicityReader,0) // Helper class for reading QWEAK helicity data
ClassDef(THcHelicityReader, 0) // Helper class for reading QWEAK helicity data
};
#endif
......
src/THcHelicityScaler.cxx 0 → 100644
View file @ e974befb
/** \class THcHelicityScaler
\ingroup Base
\brief Event handler for Hall C helicity scalers
~~~
~~~
THcHelcityScaler *phelscaler = new THcHelicityScaler("P","HC helicity scalers");
phelscaler->SetDebugFile("PHelScaler.txt");
phelscaler->SetROC(8); // 5 for HMS defaults to 8 for SHMS
phelscaler->SetBankID(9801); // Will default to this
gHaEvtHandlers->Add (phelscaler);
~~~
\authors: S. A. Wood (saw@jlab.org)
C. Yero (cyero@jlab.org)
*/
#include <bitset>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <map>
#include <sstream>
#include "TMath.h"
#include "TNamed.h"
#include "TROOT.h"
#include "TString.h"
//#include "fmt/core.h"
#include "nlohmann/json.hpp"
#include "GenScaler.h"
#include "Scaler3801.h"
#include "THaCodaData.h"
#include "THaEvData.h"
#include "THaEvtTypeHandler.h"
#include "THaGlobals.h"
#include "THcGlobals.h"
#include "THcHelicity.h"
#include "THcHelicityScaler.h"
#include "THcParmList.h"
#include "Helper.h"
#include "THaVarList.h"
#include "VarDef.h"
using namespace std;
using namespace Decoder;
static const UInt_t ICOUNT = 1;
static const UInt_t IRATE = 2;
static const UInt_t ICURRENT = 3;
static const UInt_t ICHARGE = 4;
static const UInt_t ITIME = 5;
static const UInt_t ICUT = 6;
static const UInt_t MAXCHAN = 32;
static const UInt_t defaultDT = 4;
THcHelicityScaler::THcHelicityScaler(const char* name, const char* description)
: hcana::ConfigLogging<THaEvtTypeHandler>(name, description),
fBankID(9801),
fUseFirstEvent(kTRUE),
fDelayedType(-1),
fBCM_Gain(0),
fBCM_Offset(0),
fBCM_SatOffset(0),
fBCM_SatQuadratic(0),
fBCM_delta_charge(0),
evcount(0),
evcountR(0.0),
ifound(0),
fNormIdx(-1),
fNormSlot(-1),
dvars(0),
dvarsFirst(0),
fScalerTree(0),
fOnlyBanks(kFALSE),
fClockChan(-1),
fLastClock(0) {
fRocSet.clear();
fModuleSet.clear();
//---------------------------------------------------------------------------
fROC = -1;
fNScalerChannels = 32;
AddEvtType(1);
AddEvtType(2);
AddEvtType(4);
AddEvtType(5);
AddEvtType(6);
AddEvtType(7);
AddEvtType(129);
SetDelayedType(129);
}
THcHelicityScaler::~THcHelicityScaler() {
if (!TROOT::Initialized()) {
delete fScalerTree;
}
Podd::DeleteContainer(scalers);
Podd::DeleteContainer(scalerloc);
for (auto it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
}
Int_t THcHelicityScaler::End(THaRunBase* runbase) {
// Process any delayed events in order received
for (std::vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end();
++it) {
UInt_t* rdata = *it;
AnalyzeBuffer(rdata);
}
evNumberR = -1;
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
// Write the helicity variables to scaler tree
if (fScalerTree) {
fScalerTree->Write();
}
// Compute Scaler Asymmetries
/*
for(Int_t i=0;i<fNScalerChannels;i++) {
if(fScalerSums[i]>0.5) {
fScalAsymmetry[i] = (fHScalers[0][i]-fHScalers[1][i])/fScalerSums[i];
fScalAsymmetryError[i] = 2*TMath::Sqrt(fHScalers[0][i]*fHScalers[1][i]
*fScalerSums[i])
/(fScalerSums[i]*fScalerSums[i]);
} else {
fScalAsymmetry[i] = 0.0;
fScalAsymmetryError[i] = 0.0;
}
}
*/
//------Compute Time Asymmetries-------
// C.Y. 12/15/2020 Time Asymmetry / Error Calculations (with scaler_current cut)
if (fQuartetCount <= 1) {
fTimeAsymmetry = -1000.;
fTimeAsymmetryError = 0.0;
} else {
fTimeAsymmetry = fTimeAsymSum / fQuartetCount; // normalize asymmetry to total number of
// quartets
if (fTimeAsymSum2 >= fQuartetCount * TMath::Power(fTimeAsymmetry, 2)) {
fTimeAsymmetryError =
TMath::Sqrt((fTimeAsymSum2 - fQuartetCount * TMath::Power(fTimeAsymmetry, 2)) /
(fQuartetCount * (fQuartetCount - 1)));
} else {
fTimeAsymmetryError = 0.0;
}
}
//------Compute Scaler Asymmetries-------
// C.Y. 12/15/2020 Scaler Asymmetry / Error Calculations (with scaler_current cut)
for (Int_t i = 0; i < fNScalerChannels; i++) {
if (fQuartetCount <= 1) {
fScalAsymmetry[i] = -1000.;
fScalAsymmetryError[i] = 0.0;
} else {
fScalAsymmetry[i] =
fScalAsymSum[i] / fQuartetCount; // normalize asymmetry to total number of quartets
if (fScalAsymSum2[i] >= fQuartetCount * TMath::Power(fScalAsymmetry[i], 2)) {
fScalAsymmetryError[i] =
TMath::Sqrt((fScalAsymSum2[i] - fQuartetCount * TMath::Power(fScalAsymmetry[i], 2)) /
(fQuartetCount * (fQuartetCount - 1)));
} else {
fScalAsymmetryError[i] = 0.0;
}
}
}
// json dump of helicity charge info
const int run_number = runbase->GetNumber();
j_helicity["run_number"] = run_number;
//------Compute Charge Asymmetries-------
// C.Y. 12/14/2020 Charge Asymmetry / Error Calculations (with scaler_current cut)
// Set the helicity scaler module channels for each BCM
std::map<std::string, Int_t> bcmindex;
bcmindex["BCM1_Hel.scal"] = 0;
bcmindex["BCM2_Hel.scal"] = 2;
bcmindex["Unser_Hel.scal"] = 6;
bcmindex["BCM4A_Hel.scal"] = 10;
bcmindex["BCM4B_Hel.scal"] = 4;
bcmindex["BCM4C_Hel.scal"] = 12;
// bcmindex["1MHz_Hel.scal"] = 8;
// cout << endl << "---------------------- Beam Charge Asymmetries ---------------------- " <<
// endl; cout << " BCM Total Charge Beam ON Beam ON Asymmetry" <<
// endl; cout << " Name Charge Asymmetry Charge Asymmetry Error" <<
// endl;
for (Int_t i = 0; i < fNumBCMs; i++) {
if (fQuartetCount <= 1) {
fChargeAsymmetry[i] = -1000.;
fChargeAsymmetryError[i] = 0.0;
} else {
fChargeAsymmetry[i] =
fChargeAsymSum[i] / fQuartetCount; // normalize charge asymmetry to total number of
// quartets (as the sum is for every quartet)
if (fChargeAsymSum2[i] >= fQuartetCount * TMath::Power(fChargeAsymmetry[i], 2)) {
fChargeAsymmetryError[i] = TMath::Sqrt(
(fChargeAsymSum2[i] - fQuartetCount * TMath::Power(fChargeAsymmetry[i], 2)) /
(fQuartetCount * (fQuartetCount - 1)));
} else {
fChargeAsymmetryError[i] = 0.0;
}
}
j_helicity[fBCM_Name[i]] = {{"charge", fChargeSum[i]},
{"charge_asymmetry", fChargeAsymmetry[i]},
{"charge_asymmetry_error", fChargeAsymmetryError[i]}};
// printf("%6s %12.2f %12.8f %12.2f %12.8f %12.8f\n",fBCM_Name[i].c_str(),fCharge[i],
// fChargeAsymmetry[i],fChargeSum[i],asy,asyerr);
}
// Compute +/- helicity Times (no BCM cut)
Double_t pclock = fHScalers[0][fClockChan];
Double_t mclock = fHScalers[1][fClockChan];
fTimePlus = pclock / fClockFreq;
fTimeMinus = mclock / fClockFreq;
fTime = (pclock + mclock) / fClockFreq;
// if(pclock+mclock>0) {
// fTimeAsymmetry = (pclock-mclock)/(pclock+mclock);
//} else {
// fTimeAsymmetry = 0.0;
//}
// printf("TIME(s)%12.2f %12.8f %12.2f\n",fTime, fTimeAsymmetry, fTimeSum);
j_helicity["clock"] = {{"time", fTime},
{"time_asymmetry", fTimeAsymmetry},
{"time_asymmetry_error", fTimeAsymmetryError}};
// Compute Helicity Trigger Asymmetries (no BCM cut)
if (fNTriggersPlus + fNTriggersMinus > 0) {
fTriggerAsymmetry =
((Double_t)(fNTriggersPlus - fNTriggersMinus)) / (fNTriggersPlus + fNTriggersMinus);
} else {
fTriggerAsymmetry = 0.0;
}
j_helicity["triggers"] = {{"N_triggers", fNTriggersPlus + fNTriggersMinus},
{"triggers_asymmetry", fTriggerAsymmetry}};
// write output to
return 0;
}
Int_t THcHelicityScaler::ReadDatabase(const TDatime& date) {
char prefix[2];
prefix[0] = 'g';
prefix[1] = '\0';
fNumBCMs = 0;
DBRequest list[] = {{"NumBCMs", &fNumBCMs, kInt, 0, 1}, {0}};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
if (fNumBCMs > 0) {
fBCM_Gain.resize(fNumBCMs);
fBCM_Offset.resize(fNumBCMs);
fBCM_SatOffset.resize(fNumBCMs);
fBCM_SatQuadratic.resize(fNumBCMs);
fBCM_delta_charge.resize(fNumBCMs);
string bcm_namelist;
DBRequest list2[] = {{"BCM_Gain", &fBCM_Gain[0], kDouble, (UInt_t)fNumBCMs},
{"BCM_Offset", &fBCM_Offset[0], kDouble, (UInt_t)fNumBCMs},
{"BCM_SatQuadratic", &fBCM_SatQuadratic[0], kDouble, (UInt_t)fNumBCMs, 1},
{"BCM_SatOffset", &fBCM_SatOffset[0], kDouble, (UInt_t)fNumBCMs, 1},
{"BCM_Names", &bcm_namelist, kString},
{"BCM_Current_threshold", &fbcm_Current_Threshold, kDouble, 0, 1},
{"BCM_Current_threshold_index", &fbcm_Current_Threshold_Index, kInt, 0, 1},
{0}};
fbcm_Current_Threshold = 0.0;
fbcm_Current_Threshold_Index = 0;
for (Int_t i = 0; i < fNumBCMs; i++) {
fBCM_SatOffset[i] = 0.;
fBCM_SatQuadratic[i] = 0.;
}
gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
vector<string> bcm_names = vsplit(bcm_namelist);
for (Int_t i = 0; i < fNumBCMs; i++) {
fBCM_Name.push_back(bcm_names[i] + "_Hel.scal");
fBCM_delta_charge[i] = 0.;
}
}
fTotalTime = 0.;
fPrevTotalTime = 0.;
fDeltaTime = -1.;
return kOK;
}
void THcHelicityScaler::SetDelayedType(int evtype) {
/**
* \brief Delay analysis of this event type to end.
*
* Final scaler events generated in readout list end routines may not
* come in order in the data stream. If the event type of a end routine
* scaler event is set, then the event contents will be saved and analyzed
* at the end of the analysis so that time ordering of scaler events is preserved.
*/
fDelayedType = evtype;
}
Int_t THcHelicityScaler::Analyze(THaEvData* evdata) {
// C.Y. | THcScalerEvtHandler::Analyze() uses this flag (which is forced to be 1),
// but as to why, it is beyond me. For consistency, I have also used it here.
Int_t lfirst = 1;
if (!IsMyEvent(evdata->GetEvType()))
return -1;
if (fDebugFile) {
*fDebugFile << endl << "---------------------------------- " << endl << endl;
*fDebugFile << "\nEnter THcHelicityScaler for fName = " << fName << endl;
EvDump(evdata);
}
if (lfirst && !fScalerTree) {
lfirst = 0;
// Assign a name to the helicity scaler tree
TString sname1 = "TSHel";
TString sname2 = sname1 + fName;
TString sname3 = fName + " Scaler Data";
if (fDebugFile) {
*fDebugFile << "\nAnalyze 1st time for fName = " << fName << endl;
*fDebugFile << sname2 << " " << sname3 << endl;
}
// Create Scaler Tree
fScalerTree = new TTree(sname2.Data(), sname3.Data());
fScalerTree->SetAutoSave(200000000);
TString name, tinfo;
name = "evcount";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &evcountR, tinfo.Data(), 4000);
name = "evNumber";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &evNumberR, tinfo.Data(), 4000);
// C.Y. 12/02/2020 Added actual helicity to be stored in scaler tree
name = "actualHelicity";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &actualHelicityR, tinfo.Data(), 4000);
// C.Y. 12/09/2020 Added quartetphase to be stored in scaler tree
name = "quartetPhase";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &quartetphaseR, tinfo.Data(), 4000);
for (size_t i = 0; i < scalerloc.size(); i++) {
name = scalerloc[i]->name;
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &dvars[i], tinfo.Data(), 4000);
}
}
UInt_t* rdata = (UInt_t*)evdata->GetRawDataBuffer();
if (evdata->GetEvType() == fDelayedType) { // Save this event for processing later
Int_t evlen = evdata->GetEvLength();
UInt_t* datacopy = new UInt_t[evlen];
fDelayedEvents.push_back(datacopy);
memcpy(datacopy, rdata, evlen * sizeof(UInt_t));
return 1;
}
else { // A normal event
if (fDebugFile)
*fDebugFile << "\n\nTHcHelicityScaler :: Debugging event type " << dec << evdata->GetEvType()
<< " event num = " << evdata->GetEvNum() << endl
<< endl;
evNumber = evdata->GetEvNum();
evNumberR = evNumber;
Int_t ret;
if ((ret = AnalyzeBuffer(rdata))) {
if (fDebugFile)
*fDebugFile << "scaler tree ptr 1 " << fScalerTree << endl;
if (fDebugFile)
*fDebugFile << "ret = " << ret << endl;
}
return ret;
}
}
Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata) {
fNTrigsInBuf = 0;
// Parse the data, load local data arrays.
UInt_t* p = (UInt_t*)rdata;
UInt_t* plast = p + *p; // Index to last word in the bank
Int_t roc = -1;
Int_t evlen = *p + 1;
Int_t ifound = 0;
while (p < plast) {
Int_t banklen = *p;
p++; // point to header
if (fDebugFile) {
*fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p - 1) << endl;
}
if ((*p & 0xff00) == 0x1000) { // Bank Containing banks
if (evlen - *(p - 1) > 1) { // Don't use overall event header
roc = (*p >> 16) & 0xf;
if (fDebugFile)
*fDebugFile << "ROC: " << roc << " " << evlen << " " << *(p - 1) << hex << " " << *p
<< dec << endl;
// cout << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p << dec <<
// endl;
if (roc != fROC) { // Not a ROC with helicity scaler
p += *(p - 1) - 1; // Skip to next ROC
}
}
p++; // Now pointing to a bank in the bank
} else if (((*p & 0xff00) == 0x100) && (*p != 0xC0000100)) {
// Bank containing integers. Look for scalers
// This is either ROC bank containing integers or
// a bank within a ROC containing data from modules of a single type
// Look for banks with the helicity scaler bank ID (9801)
// Assume that very first word is a scaler header
// At any point in the bank where the word is not a matching
// header, we stop.
UInt_t tag = (*p >> 16) & 0xffff; // Bank ID (ROC #)
// UInt_t num = (*p) & 0xff;
UInt_t* pnext = p + banklen; // Next bank
p++; // First data word
// If the bank is not a helicity scaler bank
// or it is not one of the ROC containing helcity scaler data
// skip to the next bank
// cout << "BankID=" << tag << endl;
if (tag != fBankID) {
p = pnext; // Fall through to end of the above else if
// cout << " Skipping to next bank" << endl;
} else {
// This is a helicity scaler bank
if (roc == fROC) {
Int_t nevents = (banklen - 2) / fNScalerChannels;
// cout << "# of helicity events in bank:" << " " << nevents << endl;
if (nevents > 100) {
cout << "Error! Beam off for too long" << endl;
}
fNTrigsInBuf = 0;
// Save helcitiy and quad info for THcHelicity
for (Int_t iev = 0; iev < nevents; iev++) { // find number of helicity events in each bank
Int_t index = fNScalerChannels * iev + 1;
// C.Y. 11/26/2020 This methods extracts the raw helicity information
AnalyzeHelicityScaler(p + index);
}
}
}
while (p < pnext) {
Int_t nskip = 0;
if (fDebugFile) {
*fDebugFile << "Scaler Header: " << hex << *p << dec;
}
if (nskip == 0) {
if (fDebugFile) {
*fDebugFile << endl;
}
break; // Didn't find a matching header
}
p = p + nskip;
}
p = pnext;
} else {
p = p + *(p - 1); // Skip to next bank
}
}
if (fDebugFile) {
*fDebugFile << "Finished with decoding. " << endl;
*fDebugFile << " Found flag = " << ifound << endl;
}
if (!ifound)
return 0;
return 1;
}
Int_t THcHelicityScaler::AnalyzeHelicityScaler(UInt_t* p) {
const static char* const here = "THcHelicityScaler::AnalyzeHelicityScaler";
Int_t hbits = (p[0] >> 30) & 0x3; // quartet and helcity bits in scaler word
Bool_t isquartet = (hbits & 2) != 0;
Int_t ispos = hbits & 1;
Int_t actualhelicity = 0;
fHelicityHistory[fNTrigsInBuf] = hbits;
fNTrigsInBuf++;
fNTriggers++;
Int_t quartetphase = (fNTriggers - fFirstCycle) % 4;
if (fFirstCycle >= -10) {
if (quartetphase == 0) {
Int_t predicted = RanBit30(fRingSeed_reported);
fRingSeed_reported = ((fRingSeed_reported << 1) | ispos) & 0x3FFFFFFF;
// Check if ringseed_predicted agrees with reported if(fNBits>=30)
if (fNBits >= 30 && predicted != fRingSeed_reported) {
_param_logger->warn("{}: Helicity Prediction Failed, Reported {:32b}, Predicted {:32b}",
here, fRingSeed_reported, predicted);
}
fNBits++;
if (fNBits == 30) {
_param_logger->info("{}: A {:32b} found at cycle {}", here, fRingSeed_reported, fNTriggers);
}
} else if (quartetphase == 3) {
if (!isquartet) {
_param_logger->warn("{}: Quartet bit expected but not set ({})", here, fNTriggers);
fNBits = 0;
fRingSeed_reported = 0;
fRingSeed_actual = 0;
fFirstCycle = -100;
}
}
} else { // First cycle not yet identified
if (isquartet) { // Helicity and quartet signal for next set of scalers
fFirstCycle = fNTriggers - 3;
quartetphase = (fNTriggers - fFirstCycle) % 4;
// Helicity at start of quartet is same as last of quartet, so we can start filling the seed
fRingSeed_reported = ((fRingSeed_reported << 1) | ispos) & 0x3FFFFFFF;
fNBits++;
if (fNBits == 30) {
_param_logger->info("{}: B {:32b} fount at cycle {}", here, fRingSeed_reported, fNTriggers);
}
}
}
if (fNBits >= 30) {
fRingSeed_actual = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
#define DELAY9
#ifdef DELAY9
if (quartetphase == 3) {
fRingSeed_actual = RanBit30(fRingSeed_actual);
actualhelicity = (fRingSeed_actual & 1) ? +1 : -1;
} else {
actualhelicity = (fRingSeed_actual & 1) ? +1 : -1;
if (quartetphase == 0 || quartetphase == 1) {
actualhelicity = -actualhelicity;
}
}
quartetphase = (quartetphase + 1) % 4;
#else
actualhelicity = (fRingSeed_actual & 1) ? +1 : -1;
if (quartetphase == 1 || quartetphase == 2) {
actualhelicity = -actualhelicity;
}
#endif
} else {
fRingSeed_actual = 0;
}
// C.Y. 12/09/2020 : Pass quartet phase to scaler tree varable
quartetphaseR = quartetphase;
// C.Y. 11/26/2020 The block of code below is used to extract the helicity information from each
// channel of the helicity scaler onto a variable to be stored in the scaler tree. For each
// channel, each helicity state (+, -, or undefined) is stored in a single varibale. Each helicity
// state is tagged to the corresponding scaler read via 'actualHelicityR' which is stored below..
// C.Y. Assign actual helicity value to a variable to be written to tree (the value may be (0,
// +hel (+1), or -hel(-1)) Where actualhelicity=0 is NOT the MPS. It is zero when the actual
// helicity has not been determined.
actualHelicityR = actualhelicity;
// C.Y. 11/26/2020 Loop over all 32 scaler channels for a specific helicity scaler module (SIS
// 3801)
for (Int_t i = 0; i < fNScalerChannels; i++) {
// C.Y. 11/26/2020 the count expression below gets the scaler raw helicity information (+, -, or
// MPS helicity states) for the ith channel
Int_t count = p[i] & 0xFFFFFF; // Bottom 24 bits equivalent of scalers->Decode()
fScalerChan[i] =
count; // pass the helicity raw information to each helicity scaler channel array element
}
// C.Y. 11/26/2020 The block of code below is used to get a cumulative sum of +/- helicity used
// for calculation of the cumulative beam charge asymmetry and other quantities in the ::End()
// method
if (actualhelicity != 0) {
// index=0 (+ helicity), index=1 (- helicity)
Int_t hindex = (actualhelicity > 0) ? 0 : 1;
// C.Y. 11/24/2020 increment the counter for either '+' or '-' helicity states
(actualhelicity > 0) ? (fNTriggersPlus++) : (fNTriggersMinus++);
// C.Y. 11/24/2020 Loop over all 32 scaler channels for a specific helicity scaler module (SIS
// 3801)
for (Int_t i = 0; i < fNScalerChannels; i++) {
Int_t count = p[i] & 0xFFFFFF; // Bottom 24 bits equivalent of scalers->Decode()
// Increment either the '+' (hindex=0) or '-' (hindex=1) helicity counts for each [i] scaler
// channel channel of a given module
fHScalers[hindex][i] += count;
fScalerSums[i] += count;
}
}
// Set the helicity scaler clock to define the time
fDeltaTime =
fScalerChan[fClockChan] /
fClockFreq; // total clock counts / clock_frequency (1MHz) for a specific scaler read interval
fTotalTime = fPrevTotalTime + fDeltaTime; // cumulative scaler time
if (fDeltaTime == 0) {
_param_logger->error("{}: ******************* Severe Warning ****************************",
here);
_param_logger->error("{}: In THcHelicityScaler have found fDeltaTime is zero !!", here);
_param_logger->error("{}: ******************* Alert DAQ experts ***************************",
here);
if (fDebugFile)
*fDebugFile << " In THcHelicityScaler have found fDeltaTime is zero !! " << endl;
}
fPrevTotalTime =
fTotalTime; // set the current total time to the previous time for the upcoming read
// C.Y. Nov 27, 2020 : (below) code to write the helicity raw data to a variable
// and map the variable to the scaler location
Double_t scal_current = 0;
// Loop over each scaler variable from the map
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t idx = scalerloc[i]->index;
size_t ichan = scalerloc[i]->ichan;
// ANALYZE 1ST SCALER READ SEPARATE (There is no previous read before this one)
if (evcount == 0) {
// Loop over the scaler variable (ivar), helicity slot (idx), and slot channel (ichan) -->
// idx=0 (only one helicity module per spectrometer arm)
if ((ivar < scalerloc.size()) && (idx < scalers.size()) && (ichan < MAXCHAN)) {
// If fUseFirstEvent=kTRUE (do not skip 1st read)
if (fUseFirstEvent) {
// Write scaler counts
if (scalerloc[ivar]->ikind == ICOUNT) {
dvars[ivar] = fScalerChan[ichan];
dvarsFirst[ivar] = 0.0;
}
// Write scaler time
if (scalerloc[ivar]->ikind == ITIME) {
dvars[ivar] = fDeltaTime;
dvarsFirst[ivar] = 0.0;
}
// Write scaler rate
if (scalerloc[ivar]->ikind == IRATE) {
dvars[ivar] = (fScalerChan[ichan]) / fDeltaTime;
dvarsFirst[ivar] = 0.0;
}
// Write either scaler current or scaler charge
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
// Get BCM index
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
// Calculate and write scaler current
if (scalerloc[ivar]->ikind == ICURRENT) {
// set default to zero
dvars[ivar] = 0.;
// Check bcm index and calculate current in a temporary placeholder, "cur_temp"
if (bcm_ind != -1) {
Double_t cur_temp =
((fScalerChan[ichan]) / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp = cur_temp + fBCM_SatQuadratic[bcm_ind] *
TMath::Max(cur_temp - fBCM_SatOffset[i], 0.0);
// Assign the calculated scaler current to dvars
dvars[ivar] = cur_temp;
}
// Check which bcm index to place a bcm current threshold cut later on. Assign the the
// beam current read to "scal_current" for later use.
if (bcm_ind == fbcm_Current_Threshold_Index) {
scal_current = dvars[ivar];
}
}
// Calculate andd write scaler charge
if (scalerloc[ivar]->ikind == ICHARGE) {
// Check bcm index and calculate current in a temporary placeholder, "cur_temp", to
// determine the charge later on.
if (bcm_ind != -1) {
Double_t cur_temp =
((fScalerChan[ichan]) / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
// Calculate the charge for this scaler read
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
// Assign the charge to dvars
dvars[ivar] = fBCM_delta_charge[bcm_ind];
}
}
}
}
// If user has set fUseFirstEvent=kFALSE (then use dvarsFirst to store the information for
// the 1st event, and leave dvars at 0.0) (The same calculations as above, but assign
// dvars=0.0, so that it does not count when written to the scaler tree)
else { // ifnotusefirstevent
if (scalerloc[ivar]->ikind == ICOUNT) {
dvarsFirst[ivar] = fScalerChan[ichan];
dvars[ivar] = 0.0;
}
if (scalerloc[ivar]->ikind == ITIME) {
dvarsFirst[ivar] = fDeltaTime;
dvars[ivar] = 0.0;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvarsFirst[ivar] = (fScalerChan[ichan]) / fDeltaTime;
dvars[ivar] = 0.0;
}
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvarsFirst[ivar] = 0.0;
dvars[ivar] = 0.0;
if (bcm_ind != -1) {
Double_t cur_temp =
((fScalerChan[ichan]) / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp + fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[i], 0.0), 2.0);
dvarsFirst[ivar] = cur_temp;
}
if (bcm_ind == fbcm_Current_Threshold_Index)
scal_current = dvarsFirst[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
Double_t cur_temp =
((fScalerChan[ichan]) / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
dvarsFirst[ivar] = fBCM_delta_charge[bcm_ind];
}
}
}
}
} else {
_param_logger->warn("{}: ERROR: incorrect index {} {} {}", here, ivar, idx, ichan);
}
}
// Calculate Scaler Quantities for ALL OTHER SCALER READS (OTHER THAN THE FIRST READ)
else { // evcount != 0
if ((ivar < scalerloc.size()) && (idx < scalers.size()) && (ichan < MAXCHAN)) {
if (scalerloc[ivar]->ikind == ICOUNT) {
dvars[ivar] = fScalerChan[ichan];
}
if (scalerloc[ivar]->ikind == ITIME) {
dvars[ivar] = fDeltaTime;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvars[ivar] = fScalerChan[ichan] / fDeltaTime;
}
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvars[ivar] = 0.;
if (bcm_ind != -1) {
if (fDeltaTime > 0) {
Double_t cur_temp =
(fScalerChan[ichan] / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
dvars[ivar] = cur_temp;
}
}
if (bcm_ind == fbcm_Current_Threshold_Index)
scal_current = dvars[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
dvars[ivar] = 0.0;
if (fDeltaTime > 0) {
Double_t cur_temp =
(fScalerChan[ichan] / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
}
dvars[ivar] = fBCM_delta_charge[bcm_ind];
}
}
}
} else {
_param_logger->warn("{}: ERROR: incorrect index {} {} {}", here, ivar, idx, ichan);
}
}
}
// Analyze Scaler Reads ONLY FOR SCAL_CURRENTS > BCM_CURRENT THRESHOLD
//(These variables have the name structure: varibaleName_Cut)
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t ichan = scalerloc[i]->ichan;
if (scalerloc[ivar]->ikind == ICUT + ICOUNT) {
if (scal_current > fbcm_Current_Threshold) {
dvars[ivar] = fScalerChan[ichan];
}
}
if (scalerloc[ivar]->ikind == ICUT + ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scal_current > fbcm_Current_Threshold && bcm_ind != -1) {
dvars[ivar] = fBCM_delta_charge[bcm_ind];
}
}
if (scalerloc[ivar]->ikind == ICUT + ITIME) {
if (scal_current > fbcm_Current_Threshold) {
dvars[ivar] = fDeltaTime;
}
}
}
//--------------------------------------
// C.Y. 12/13/2020 : Calculate the asymmetries / errors at the end of each quartet (S.A. Wood
// approach)
if (actualhelicity != 0) {
// Reset fHaveCycle to kFALSE at the start of each quartet (e.g. - + + -, + - - +
if (quartetphase == 0) {
fHaveCycle[0] = fHaveCycle[1] = fHaveCycle[2] = fHaveCycle[3] = kFALSE;
}
// Check if BCM scaler current is above set threshold
if (scal_current > fbcm_Current_Threshold &&
(quartetphase == 0 || fHaveCycle[max(quartetphase - 1, 0)])) {
fHaveCycle[quartetphase] = kTRUE;
// Loop over each BCM to get the charge for each cycle of a quartet
for (Int_t i = 0; i < fNumBCMs; i++) {
fChargeCycle[quartetphase][i] = fBCM_delta_charge[i];
}
// Loop over each Scaler Channel to the the counts for each cycle of a quartet
for (Int_t i = 0; i < fNScalerChannels; i++) {
fScalCycle[quartetphase][i] = fScalerChan[i];
}
// Set the time for each cycle of the quartet
fTimeCycle[quartetphase] = fDeltaTime;
}
// Compute asymmetries at the end of each quartet
if (quartetphase == 3 && fHaveCycle[3]) {
// Loop over BCMs
for (Int_t i = 0; i < fNumBCMs; i++) {
// compute charge asymmetry for each quartet at the end of said quartet
Double_t asy =
actualhelicity *
(fChargeCycle[0][i] + fChargeCycle[3][i] - fChargeCycle[1][i] - fChargeCycle[2][i]) /
(fChargeCycle[0][i] + fChargeCycle[3][i] + fChargeCycle[1][i] + fChargeCycle[2][i]);
fChargeSum[i] +=
fChargeCycle[0][i] + fChargeCycle[1][i] + fChargeCycle[2][i] + fChargeCycle[3][i];
// keep track of sums for proper error calculation
fChargeAsymSum[i] += asy;
fChargeAsymSum2[i] += asy * asy;
}
//-------
// Loop over Scaler Channels
for (Int_t i = 0; i < fNScalerChannels; i++) {
// compute scaler asymmetry for each quartet at the end of said quartet
Double_t asy = actualhelicity *
(fScalCycle[0][i] + fScalCycle[3][i] - fScalCycle[1][i] - fScalCycle[2][i]) /
(fScalCycle[0][i] + fScalCycle[3][i] + fScalCycle[1][i] + fScalCycle[2][i]);
fScalSum[i] += fScalCycle[0][i] + fScalCycle[1][i] + fScalCycle[2][i] + fScalCycle[3][i];
// keep track of sums for proper error calculation
fScalAsymSum[i] += asy;
fScalAsymSum2[i] += asy * asy;
}
//-------
// compute time asymmetry for each quartet at the end of said quartet
Double_t asy = actualhelicity *
(fTimeCycle[0] + fTimeCycle[3] - fTimeCycle[1] - fTimeCycle[2]) /
(fTimeCycle[0] + fTimeCycle[3] + fTimeCycle[1] + fTimeCycle[2]);
// keep track of total time
fTimeSum += fTimeCycle[0] + fTimeCycle[1] + fTimeCycle[2] + fTimeCycle[3];
// keep track of sums for proper error calculation
fTimeAsymSum += asy;
fTimeAsymSum2 += asy * asy;
//------
// keep track of the total number of quartets
fQuartetCount++;
}
}
//--------------------------------------
// increment scaler reads
evcount = evcount + 1;
evcountR = evcount;
// clear Genscaler scalers
for (size_t j = 0; j < scalers.size(); j++)
scalers[j]->Clear("");
// C.Y. 12/02/2020 Fill Scaler Tree Here
if (fScalerTree) {
fScalerTree->Fill();
}
return (0);
}
//_____________________________________________________________________________
Int_t THcHelicityScaler::RanBit30(Int_t ranseed) {
UInt_t bit7 = (ranseed & 0x00000040) != 0;
UInt_t bit28 = (ranseed & 0x08000000) != 0;
UInt_t bit29 = (ranseed & 0x10000000) != 0;
UInt_t bit30 = (ranseed & 0x20000000) != 0;
UInt_t newbit = (bit30 ^ bit29 ^ bit28 ^ bit7) & 0x1;
ranseed = ((ranseed << 1) | newbit) & 0x3FFFFFFF;
return ranseed;
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHelicityScaler::Init(const TDatime& date) {
const static char* const here = "THcHelicityScaler::Init";
ReadDatabase(date);
const int LEN = 200;
char cbuf[LEN];
fStatus = kOK;
fNormIdx = -1;
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
_param_logger->info("Howdy! We are initializing THcHelicityScaler !! name = {}", fName.Data());
if (eventtypes.size() == 0) {
eventtypes.push_back(0); // Default Event Type
}
if (fROC < 0) {
fROC = 8; // Default to SHMS crate
}
// C.Y. 11/26/2020 : Read In and Parse the variables in the helicity scaler map file
TString dfile;
dfile = fName + "scaler.txt";
TString sname0 = "Scalevt";
TString sname;
sname = "hel" + fName + sname0; // This should be: 'helPScalevt' or 'helHScalevt'
// Open helicity scaler .dat map file
FILE* fi = OpenFile(sname.Data(), date);
if (!fi) {
_param_logger->warn("{}: Cannot find db file for {} scaler event handler", here, fName);
return kFileError;
}
size_t minus1 = -1;
size_t pos1;
string scomment = "#";
string svariable = "variable";
string smap = "map";
vector<string> dbline;
while (fgets(cbuf, LEN, fi) != NULL) {
std::string sin(cbuf);
std::string sinput(sin.substr(0, sin.find_first_of("#")));
if (fDebugFile)
*fDebugFile << "string input " << sinput << endl;
dbline = vsplit(sinput);
if (dbline.size() > 0) {
pos1 = FindNoCase(dbline[0], scomment);
if (pos1 != minus1)
continue;
pos1 = FindNoCase(dbline[0], svariable);
if (pos1 != minus1 && dbline.size() > 4) {
string sdesc = "";
for (size_t j = 5; j < dbline.size(); j++)
sdesc = sdesc + " " + dbline[j];
UInt_t islot = atoi(dbline[1].c_str());
UInt_t ichan = atoi(dbline[2].c_str());
UInt_t ikind = atoi(dbline[3].c_str());
if (fDebugFile)
*fDebugFile << "add var " << dbline[1] << " desc = " << sdesc << " islot= " << islot
<< " " << ichan << " " << ikind << endl;
TString tsname(dbline[4].c_str());
TString tsdesc(sdesc.c_str());
AddVars(tsname, tsdesc, islot, ichan, ikind);
// add extra scaler which is cut on the current
if (ikind == ICOUNT || ikind == ITIME || ikind == ICHARGE) {
tsname = tsname + "Cut";
AddVars(tsname, tsdesc, islot, ichan, ICUT + ikind);
}
}
pos1 = FindNoCase(dbline[0], smap);
if (fDebugFile)
*fDebugFile << "map ? " << dbline[0] << " " << smap << " " << pos1 << " "
<< dbline.size() << endl;
if (pos1 != minus1 && dbline.size() > 6) {
Int_t imodel, icrate, islot, inorm;
UInt_t header, mask;
char cdum[20];
sscanf(sinput.c_str(), "%s %d %d %d %x %x %d \n", cdum, &imodel, &icrate, &islot, &header,
&mask, &inorm);
if ((fNormSlot >= 0) && (fNormSlot != inorm))
_param_logger->warn("{}: contradictory norm slot {} {}", here, fNormSlot, inorm);
fNormSlot =
inorm; // slot number used for normalization. This variable is not used but is checked.
Int_t clkchan = -1;
Double_t clkfreq = 1;
if (dbline.size() > 8) {
clkchan = atoi(dbline[7].c_str());
clkfreq = 1.0 * atoi(dbline[8].c_str());
fClockChan = clkchan;
fClockFreq = clkfreq;
}
if (fDebugFile) {
*fDebugFile << "map line " << dec << imodel << " " << icrate << " " << islot << endl;
*fDebugFile << " header 0x" << hex << header << " 0x" << mask << dec << " " << inorm
<< " " << clkchan << " " << clkfreq << endl;
}
switch (imodel) {
case 3801: // Hall C Helicity Scalers (SIS 3801)
scalers.push_back(new Scaler3801(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
}
if (scalers.size() > 0) {
UInt_t idx = scalers.size() - 1;
// Headers must be unique over whole event, not
// just within a ROC
scalers[idx]->SetHeader(header, mask);
// The normalization slot has the clock in it, so we automatically recognize it.
// fNormIdx is the index in scaler[] and
// fNormSlot is the slot#, checked for consistency
if (clkchan >= 0) {
scalers[idx]->SetClock(defaultDT, clkchan, clkfreq);
_param_logger->info("{}: Setting scaler clock ... channel = {} ... freq = {}", here,
clkchan, clkfreq);
if (fDebugFile)
*fDebugFile << "Setting scaler clock ... channel = " << clkchan
<< " ... freq = " << clkfreq << endl;
fNormIdx = idx;
if (islot != fNormSlot)
_param_logger->warn("{}: contradictory norm slot {}", here, islot);
}
}
}
}
} // end while loop
// can't compare UInt_t to Int_t (compiler warning), so do this
nscalers = 0;
for (size_t i = 0; i < scalers.size(); i++)
nscalers++;
// need to do LoadNormScaler after scalers created and if fNormIdx found
if (fDebugFile)
*fDebugFile << "fNormIdx = " << fNormIdx << endl;
if ((fNormIdx >= 0) && fNormIdx < nscalers) {
for (Int_t i = 0; i < nscalers; i++) {
if (i == fNormIdx)
continue;
scalers[i]->LoadNormScaler(scalers[fNormIdx]);
}
}
// Called after AddVars() has been called
DefVars();
// Verify that the slots are not defined twice
for (UInt_t i1 = 0; i1 < scalers.size() - 1; i1++) {
for (UInt_t i2 = i1 + 1; i2 < scalers.size(); i2++) {
if (scalers[i1]->GetSlot() == scalers[i2]->GetSlot())
_param_logger->warn("{}: same slot defined twice", here);
}
}
// Identify indices of scalers[] vector to variables.
for (UInt_t i = 0; i < scalers.size(); i++) {
for (UInt_t j = 0; j < scalerloc.size(); j++) {
if (scalerloc[j]->islot == static_cast<UInt_t>(scalers[i]->GetSlot()))
scalerloc[j]->index = i;
}
}
if (fDebugFile)
*fDebugFile << "THcHelicityScaler:: Name of scaler bank " << fName << endl;
for (size_t i = 0; i < scalers.size(); i++) {
if (fDebugFile) {
*fDebugFile << "Scaler # " << i << endl;
scalers[i]->SetDebugFile(fDebugFile);
scalers[i]->DebugPrint(fDebugFile);
}
}
//------Initialize Helicity Variables / Arrays-----
fNTriggers = 0;
fNTrigsInBuf = 0;
fFirstCycle = -100;
fRingSeed_reported = 0;
fRingSeed_actual = 0;
fNBits = 0;
fNTriggersPlus = fNTriggersMinus = 0;
fHScalers[0].resize(fNScalerChannels); //+ helicity scaler counts
fHScalers[1].resize(fNScalerChannels); //- helicity scaler counts
fScalerSums.resize(fNScalerChannels); // total helicity scaler counts (hel=0, excluded)
fScalerChan.resize(fNScalerChannels); // C.Y. 11/26/2020 : added array to store helicity
// information per channel
for (Int_t i = 0; i < fNScalerChannels; i++) {
fHScalers[0][i] = 0.0;
fHScalers[1][i] = 0.0;
fScalerSums[i] = 0.0;
fScalerChan[i] = 0.0;
}
fTimePlus = fTimeMinus = 0.0;
fTriggerAsymmetry = 0.0;
//------C.Y.: 12/13/2020 Initialize Variables for quartet-by-quartet asymmetries (include BCM
// current threshold cut)--------- (and error calculation)
for (Int_t i = 0; i < 4; i++) {
fChargeCycle[i].resize(fNumBCMs);
fScalCycle[i].resize(fNScalerChannels);
fHaveCycle[i] = kFALSE;
for (Int_t j = 0; j < fNumBCMs; j++) {
fChargeCycle[i][j] = 0.0;
}
for (Int_t j = 0; j < fNScalerChannels; j++) {
fScalCycle[i][j] = 0.0;
}
fTimeCycle[i] = 0.0;
}
// Initialize quartet counter
fQuartetCount = 0.0;
// Initialize variables for time asymmetry calculation
fTimeSum = 0.0;
fTimeAsymmetry = 0.0;
fTimeAsymmetryError = 0.0;
fTimeAsymSum = 0.0;
fTimeAsymSum2 = 0.0;
// Initialize variables for charge asymmetry calculation
fChargeSum.resize(fNumBCMs);
fChargeAsymmetry.resize(fNumBCMs);
fChargeAsymmetryError.resize(fNumBCMs);
fChargeAsymSum.resize(fNumBCMs);
fChargeAsymSum2.resize(fNumBCMs);
for (Int_t i = 0; i < fNumBCMs; i++) {
fChargeSum[i] = 0.0;
fChargeAsymmetry[i] = 0.0;
fChargeAsymSum[i] = 0.0;
fChargeAsymSum2[i] = 0.0;
}
// Initialize variables for scaler asymmetry calculation
fScalSum.resize(fNScalerChannels);
fScalAsymmetry.resize(fNScalerChannels);
fScalAsymmetryError.resize(fNScalerChannels);
fScalAsymSum.resize(fNScalerChannels);
fScalAsymSum2.resize(fNScalerChannels);
for (Int_t i = 0; i < fNScalerChannels; i++) {
fScalSum[i] = 0.0;
fScalAsymmetry[i] = 0.0;
fScalAsymSum[i] = 0.0;
fScalAsymSum2[i] = 0.0;
}
//------------------------------------------------------------------------------------------
// Call MakeParms() to define variables to be used in report file
MakeParms();
return kOK;
}
void THcHelicityScaler::MakeParms() {
// Put Various helicity scaler results in gHcParms so they can be included in results.
// no bcm cut required
gHcParms->Define(Form("g%s_hscaler_plus[%d]", fName.Data(), fNScalerChannels),
"Plus Helcity Scalers", fHScalers[0].front());
gHcParms->Define(Form("g%s_hscaler_minus[%d]", fName.Data(), fNScalerChannels),
"Minus Helcity Scalers", fHScalers[1].front());
// gHcParms->Define(Form("g%s_hscaler_sum[%d]",fName.Data(),fNScalerChannels),
// "Helcity Scalers Sum",*fScalerSums);
gHcParms->Define(Form("g%s_hscaler_triggers", fName.Data()), "Total Helicity Scaler Triggers",
fNTriggers);
gHcParms->Define(Form("g%s_hscaler_triggers_plus", fName.Data()),
"Positive Helicity Scaler Triggers", fNTriggersPlus);
gHcParms->Define(Form("g%s_hscaler_triggers_minus", fName.Data()),
"Negative Helicity Scaler Triggers", fNTriggersMinus);
gHcParms->Define(Form("g%s_hscaler_trigger_asy", fName.Data()), "Helicity Trigger Asymmetry",
fTriggerAsymmetry);
// gHcParms->Define(Form("g%s_hscaler_time_plus",fName.Data()),
// "Positive Helicity Time",fTimePlus);
// gHcParms->Define(Form("g%s_hscaler_time_minus",fName.Data()),
// "Negative Helicity Time",fTimeMinus);
// bcm cut
gHcParms->Define(Form("g%s_hscaler_sum[%d]", fName.Data(), fNScalerChannels),
"Helcity Scalers Sum", fScalSum[0]);
gHcParms->Define(Form("g%s_hscaler_asy[%d]", fName.Data(), fNScalerChannels),
"Helicity Scaler Asymmetry[%d]", fScalAsymmetry[0]);
gHcParms->Define(Form("g%s_hscaler_asyerr[%d]", fName.Data(), fNScalerChannels),
"Helicity Scaler Asymmetry Error[%d]", fScalAsymmetryError[0]);
gHcParms->Define(Form("g%s_hscaler_charge[%d]", fName.Data(), fNumBCMs), "Helicity Gated Charge",
fChargeSum[0]);
gHcParms->Define(Form("g%s_hscaler_charge_asy[%d]", fName.Data(), fNumBCMs),
"Helicity Gated Charge Asymmetry", fChargeAsymmetry[0]);
gHcParms->Define(Form("g%s_hscaler_charge_asyerr[%d]", fName.Data(), fNumBCMs),
"Helicity Gated Charge Asymmetry Error", fChargeAsymmetryError[0]);
gHcParms->Define(Form("g%s_hscaler_time", fName.Data()), "Helicity Gated Time (sec)", fTimeSum);
gHcParms->Define(Form("g%s_hscaler_time_asy", fName.Data()), "Helicity Gated Time Asymmetry",
fTimeAsymmetry);
gHcParms->Define(Form("g%s_hscaler_time_asyerr", fName.Data()),
"Helicity Gated Time Asymmetry Error", fTimeAsymmetryError);
}
//--------------------C.Y. Sep 20, 2020 : Added Utility Functions--------------------
void THcHelicityScaler::AddVars(TString name, TString desc, UInt_t islot, UInt_t ichan,
UInt_t ikind) {
if (fDebugFile)
*fDebugFile << "C.Y. | Calling THcHelicityScaler::AddVars() " << endl;
// need to add fName here to make it a unique variable. (Left vs Right HRS, for example)
TString name1 = fName + name;
TString desc1 = fName + desc;
// We don't yet know the correspondence between index of scalers[] and slots.
// Will put that in later.
scalerloc.push_back(new HCScalerLoc(name1, desc1, 0, islot, ichan, ikind, scalerloc.size()));
}
void THcHelicityScaler::DefVars() {
const char* const here = "THcHelicityScaler::DefVars()";
if (fDebugFile)
*fDebugFile << "C.Y. | Calling THcHelicityScaler::DefVars() " << endl;
// called after AddVars has finished being called.
Nvars = scalerloc.size();
if (Nvars == 0)
return;
dvars.resize(Nvars); // dvars is a member of this class
dvarsFirst.resize(Nvars); // dvarsFirst is a member of this class
dvars = {0.};
dvarsFirst = {0.};
if (gHaVars) {
if (fDebugFile)
*fDebugFile << "THcScalerEVtHandler:: Have gHaVars " << gHaVars << endl;
} else {
_param_logger->warn("{}: No gHaVars ?! Well, that's a problem !!", here);
return;
}
if (fDebugFile)
*fDebugFile << "THcHelicityScaler:: scalerloc size " << scalerloc.size() << endl;
const Int_t* count = 0;
for (size_t i = 0; i < scalerloc.size(); i++) {
gHaVars->DefineByType(scalerloc[i]->name.Data(), scalerloc[i]->description.Data(), &dvars[i],
kDouble, count);
}
}
size_t THcHelicityScaler::FindNoCase(const string& sdata, const string& skey) {
// Find iterator of word "sdata" where "skey" starts. Case insensitive.
string sdatalc, skeylc;
sdatalc = "";
skeylc = "";
for (string::const_iterator p = sdata.begin(); p != sdata.end(); ++p) {
sdatalc += tolower(*p);
}
for (string::const_iterator p = skey.begin(); p != skey.end(); ++p) {
skeylc += tolower(*p);
}
if (sdatalc.find(skeylc, 0) == string::npos)
return -1;
return sdatalc.find(skeylc, 0);
};
void THcHelicityScaler::WriteJson(const std::string& path) const {
std::ofstream ofile{
fmt::format("{}/hel_scalers_{}.json", path, j_helicity["run_number"].get<int>())};
ofile << std::setw(4) << j_helicity << "\n";
}
//---------------------------------------------------------------------------------
ClassImp(THcHelicityScaler)
src/THcHelicityScaler.h 0 → 100644
View file @ e974befb
#ifndef THcHelicityScaler_
#define THcHelicityScaler_
/////////////////////////////////////////////////////////////////////
//
// THcHelicityScaler
// class to handle Helicity scaler events
//
/////////////////////////////////////////////////////////////////////
#include "THaEvtTypeHandler.h"
#include "THcScalerEvtHandler.h"
#include "Decoder.h"
#include "THaRunBase.h"
#include "TString.h"
#include "TTree.h"
#include "hcana/Logger.h"
#include <algorithm>
#include <cstring>
#include <nlohmann/json.hpp>
#include <set>
#include <string>
#include <vector>
class THcHelicity;
class HCScalerLoc;
class THcHelicityScaler : public hcana::ConfigLogging<THaEvtTypeHandler> {
public:
THcHelicityScaler(const char*, const char*);
virtual ~THcHelicityScaler();
Int_t Analyze(THaEvData* evdata);
Int_t AnalyzeBuffer(UInt_t* rdata);
Int_t AnalyzeHelicityScaler(UInt_t* p);
virtual EStatus Init(const TDatime& run_time);
virtual Int_t ReadDatabase(const TDatime& date);
virtual Int_t End(THaRunBase* r = 0);
virtual void SetUseFirstEvent(Bool_t b = kFALSE) { fUseFirstEvent = b; }
virtual void SetDelayedType(int evtype);
virtual void SetROC(Int_t roc) { fROC = roc; }
virtual void SetBankID(Int_t bankid) { fBankID = bankid; }
virtual void SetNScalerChannels(Int_t n) { fNScalerChannels = n; }
virtual Int_t GetNevents() { return fNTrigsInBuf; }
virtual Int_t GetNcycles() { return fNTriggers; }
virtual Int_t GetEvNum() { return evNumber; }
virtual Int_t* GetHelicityHistoryP() { return fHelicityHistory; }
virtual Int_t GetReportedSeed() { return fRingSeed_reported; }
virtual Int_t GetReportedActual() { return fRingSeed_actual; }
virtual Bool_t IsSeedGood() { return fNBits >= 30; }
Double_t GetPlusCharge(const std::string& name);
Double_t GetMinusCharge(const std::string& name);
// utility function to write out json helicity file
void WriteJson(const std::string& path) const;
private:
//------------C.Y. Sep 20, 2020 :: Added Utility Function Prototypes----------------
void AddVars(TString name, TString desc, UInt_t iscal, UInt_t ichan, UInt_t ikind);
void DefVars();
static size_t FindNoCase(const std::string& sdata, const std::string& skey);
std::vector<Decoder::GenScaler*> scalers;
std::vector<HCScalerLoc*> scalerloc;
Int_t RanBit30(Int_t ranseed);
void MakeParms();
UInt_t fBankID;
// Helicity Scaler variables
Int_t fNTrigsInBuf; /* # of helicity scaler reads in last event */
Int_t fNTriggers;
Int_t fFirstCycle;
Int_t fHelicityHistory[200];
//
Bool_t fUseFirstEvent;
Int_t fDelayedType;
Int_t fRingSeed_reported;
Int_t fRingSeed_actual;
Int_t fNBits;
Int_t fNTriggersPlus;
Int_t fNTriggersMinus;
std::vector<Double_t> fHScalers[2];
Int_t fGateCount[2];
std::vector<Double_t> fScalerSums, fAsymmetry, fAsymmetryError;
std::vector<Double_t> fCharge;
Double_t fTime;
Double_t fTimePlus;
Double_t fTimeMinus;
Double_t fTriggerAsymmetry;
//---- C.Y.: 12/14/2020 Variables for quartet-by-quartet asymmetry/error calculations ----
Bool_t fHaveCycle[4];
Int_t fQuartetCount; // keep track of number of quartets
// quartet-by-quartet time asymmetry variables
Double_t fTimeCycle[4];
Double_t fTimeSum;
Double_t fTimeAsymmetry;
Double_t fTimeAsymmetryError;
Double_t fTimeAsymSum;
Double_t fTimeAsymSum2;
// quartet-by-quartet scaler counts asymmetry variables
std::vector<Double_t> fScalCycle[4];
std::vector<Double_t> fScalSum; // reminder: need to initialize
std::vector<Double_t> fScalAsymmetry;
std::vector<Double_t> fScalAsymmetryError;
std::vector<Double_t> fScalAsymSum;
std::vector<Double_t> fScalAsymSum2;
// quartet-by-quartet charge asymmetry variables
std::vector<Double_t> fChargeCycle[4];
std::vector<Double_t> fChargeSum;
std::vector<Double_t> fChargeAsymmetry;
std::vector<Double_t> fChargeAsymmetryError;
std::vector<Double_t> fChargeAsymSum;
std::vector<Double_t> fChargeAsymSum2;
//----------------------
//----C.Y. Nov 26, 2020----
std::vector<Double_t> fScalerChan;
std::vector<UInt_t*> fDelayedEvents;
Int_t fROC;
Int_t fNScalerChannels; // Number of scaler channels/event
Int_t fNumBCMs;
std::vector<Double_t> fBCM_Gain, fBCM_Offset;
std::vector<std::string> fBCM_Name;
//---C.Y. Sep 2020 : Added additional BCM-related variables--
std::vector<Double_t> fBCM_SatOffset;
std::vector<Double_t> fBCM_SatQuadratic;
std::vector<Double_t> fBCM_delta_charge;
Double_t fTotalTime;
Double_t fDeltaTime;
Double_t fPrevTotalTime;
Double_t fbcm_Current_Threshold;
Double_t fClockFreq;
Int_t fbcm_Current_Threshold_Index;
//----C.Y. Sep 20, 2020 : Added additional variables-----
// (required by utility functions and scaler tree output)
UInt_t evcount;
Double_t evcountR;
UInt_t evNumber;
Double_t evNumberR;
Double_t actualHelicityR;
Double_t quartetphaseR;
Int_t Nvars, ifound, fNormIdx, fNormSlot, nscalers;
std::vector<Double_t> dvars;
std::vector<Double_t> dvarsFirst;
TTree* fScalerTree;
Bool_t fOnlySyncEvents;
Bool_t fOnlyBanks;
Int_t fClockChan;
UInt_t fLastClock;
std::set<UInt_t> fRocSet;
std::set<UInt_t> fModuleSet;
//--------------------------------------------------------
// json file with helicity info
nlohmann::json j_helicity;
THcHelicityScaler(const THcHelicityScaler& fh);
THcHelicityScaler& operator=(const THcHelicityScaler& fh);
ClassDef(THcHelicityScaler, 0) // Scaler Event handler
};
#endif
src/THcHitList.cxx
View file @ e974befb
......@@ -263,10 +263,14 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
fRefIndexMaps[i].hashit = kFALSE;
Bool_t goodreftime=kFALSE;
Int_t reftime = 0;
Int_t prevtime = 0;
Int_t difftime = 0;
for(Int_t ihit=0; ihit<nrefhits; ihit++) {
reftime = evdata.GetData(Decoder::kPulseTime,fRefIndexMaps[i].crate,
fRefIndexMaps[i].slot, fRefIndexMaps[i].channel,ihit);
reftime += 64*timeshift;
if (ihit != 0) difftime=reftime-prevtime;
prevtime = reftime;
if(reftime >= fADC_RefTimeCut) {
goodreftime = kTRUE;
break;
......@@ -274,6 +278,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
}
if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
fRefIndexMaps[i].reftime = reftime;
fRefIndexMaps[i].refdifftime = difftime;
fRefIndexMaps[i].hashit = kTRUE;
}
} else { // Assume this is a TDC
......@@ -285,9 +290,13 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
// then fTDC_RefTimeCut
Bool_t goodreftime=kFALSE;
Int_t reftime = 0;
Int_t prevtime = 0;
Int_t difftime = 0;
for(Int_t ihit=0; ihit<nrefhits; ihit++) {
reftime = evdata.GetData(fRefIndexMaps[i].crate,fRefIndexMaps[i].slot,
fRefIndexMaps[i].channel,ihit);
if( ihit != 0) difftime=reftime-prevtime;
prevtime=reftime;
if(reftime >= fTDC_RefTimeCut) {
goodreftime = kTRUE;
break;
......@@ -295,6 +304,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
}
if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
fRefIndexMaps[i].reftime = reftime;
fRefIndexMaps[i].refdifftime = difftime;
fRefIndexMaps[i].hashit = kTRUE;
}
}
......@@ -360,8 +370,12 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
Int_t nrefhits = evdata.GetNumHits(d->crate,d->slot,d->refchan);
Bool_t goodreftime=kFALSE;
Int_t reftime=0;
Int_t prevtime=0;
Int_t difftime=0;
for(Int_t ihit=0; ihit<nrefhits; ihit++) {
reftime = evdata.GetData(d->crate, d->slot, d->refchan, ihit);
if (ihit != 0 ) difftime=reftime-prevtime;
prevtime = reftime;
if(reftime >= fTDC_RefTimeCut) {
goodreftime = kTRUE;
break;
......@@ -371,6 +385,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
// hits make the RefTimeCut
if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
rawhit->SetReference(signal, reftime);
rawhit->SetReferenceDiff(signal, difftime);
} else if (!suppresswarnings) {
cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
" missing for (" << d->crate << ", " << d->slot <<
......@@ -381,6 +396,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
if(d->refindex >=0 && d->refindex < fNRefIndex) {
if(fRefIndexMaps[d->refindex].hashit) {
rawhit->SetReference(signal, fRefIndexMaps[d->refindex].reftime);
rawhit->SetReferenceDiff(signal, fRefIndexMaps[d->refindex].refdifftime);
} else {
if(!suppresswarnings) {
cout << "HitList(event=" << evdata.GetEvNum() << "): refindex " << d->refindex <<
......@@ -399,12 +415,17 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
if (fPSE125) {
if(!fHaveFADCInfo) {
fNSA = fPSE125->GetNSA(d->crate);
fNSB = fPSE125->GetNSB(d->crate);
fNPED = fPSE125->GetNPED(d->crate);
//fNSA = fPSE125->GetNSA(d->crate);
//fNSB = fPSE125->GetNSB(d->crate);
//fNPED = fPSE125->GetNPED(d->crate);
// Hard coding for now so we don't need event typ 125 at the start of run.
fNSA = 26;
fNSB = 3;
fNPED = 4;
fHaveFADCInfo = kTRUE;
}
// Set F250 parameters.
}
//std::cout << fNSA << ", " << fNSB << ", " << fNPED << "\n";
// Set F250 parameters.
rawhit->SetF250Params(fNSA, fNSB, fNPED);
}
......@@ -446,6 +467,8 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
d->crate, d->slot, d->refchan);
Bool_t goodreftime=kFALSE;
Int_t reftime = 0;
Int_t prevtime = 0;
Int_t difftime = 0;
timeshift=0;
if(fTISlot>0) { // Get the trigger time for this module
if(fTrigTimeShiftMap.find(d->slot)
......@@ -460,6 +483,8 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
for(Int_t ihit=0; ihit<nrefhits; ihit++) {
reftime = evdata.GetData(Decoder::kPulseTime, d->crate, d->slot, d->refchan, ihit);
reftime += 64*timeshift;
if (ihit != 0) difftime=reftime-prevtime;
prevtime=reftime;
if(reftime >= fADC_RefTimeCut) {
goodreftime=kTRUE;
break;
......@@ -469,6 +494,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
// hits make the RefTimeCut
if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
rawhit->SetReference(signal, reftime);
rawhit->SetReferenceDiff(signal, difftime);
} else if (!suppresswarnings) {
#ifndef SUPPRESSMISSINGADCREFTIMEMESSAGES
cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
......@@ -481,6 +507,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
if(d->refindex >=0 && d->refindex < fNRefIndex) {
if(fRefIndexMaps[d->refindex].hashit) {
rawhit->SetReference(signal, fRefIndexMaps[d->refindex].reftime);
rawhit->SetReferenceDiff(signal, fRefIndexMaps[d->refindex].refdifftime);
} else {
if(!suppresswarnings) {
#ifndef SUPPRESSMISSINGADCREFTIMEMESSAGES
......
src/THcHitList.h
View file @ e974befb
......@@ -73,6 +73,7 @@ protected:
Int_t slot;
Int_t channel;
Int_t reftime;
Int_t refdifftime;
};
std::vector<RefIndexMap> fRefIndexMaps;
// Should this be a sparse list instead in case user
......
src/THcHodoHit.h
View file @ e974befb
......@@ -34,6 +34,9 @@ public:
Double_t GetNegADCpeak() const { return fNegADC_Peak; }
Double_t GetPosADCtime() const { return fPosADC_Time; }
Double_t GetNegADCtime() const { return fNegADC_Time; }
Double_t GetPosADCCorrtime() const { return fPosADC_CorrTime; }
Double_t GetNegADCCorrtime() const { return fNegADC_CorrTime; }
Double_t GetCalcPosition() const { return fCalcPosition; }
Int_t GetPosTDC() const { return fPosTDC; }
Int_t GetNegTDC() const { return fNegTDC; }
Double_t GetPosCorrectedTime() const { return fPosCorrectedTime;}
......@@ -46,8 +49,9 @@ public:
Int_t GetPaddleNumber() const { return fPaddleNumber; }
Double_t GetPaddleCenter() const { return fPaddleCenter; }
void SetCorrectedTimes(Double_t pos, Double_t neg, Double_t) {
void SetCorrectedTimes(Double_t pos, Double_t neg) {
fPosCorrectedTime = pos; fNegCorrectedTime = neg;
fHasCorrectedTimes = kFALSE;
}
void SetCorrectedTimes(Double_t pos, Double_t neg,
Double_t postof, Double_t negtof,
......@@ -75,6 +79,15 @@ public:
void SetNegADCtime( Double_t ptime) {
fNegADC_Time =ptime;
}
void SetPosADCCorrtime( Double_t ptime) {
fPosADC_CorrTime =ptime;
}
void SetNegADCCorrtime( Double_t ptime) {
fNegADC_CorrTime =ptime;
}
void SetCalcPosition( Double_t calcpos) {
fCalcPosition =calcpos;
}
protected:
static const Double_t kBig; //!
......@@ -86,6 +99,9 @@ protected:
Double_t fNegADC_Peak; // ADC peak amplitude
Double_t fPosADC_Time; // ADC time
Double_t fNegADC_Time; // ADC time
Double_t fPosADC_CorrTime; // ADC time
Double_t fNegADC_CorrTime; // ADC time
Double_t fCalcPosition; // Position along paddle calculated by time diff
Int_t fPaddleNumber;
Double_t fPosCorrectedTime; // Pulse height corrected time
......