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 1893 additions and 673 deletions
containers/docker/README.md deleted 100644 → 0
View file @ 36a71670
Hall A/C Software
=================
Images for running containers for all aspects of the SANE experimental
analysis.
The starting point is a pre-built image for the ROOT libraries. (ubuntu + ROOT)
Main software libraries:
- `evio`: Built from https://github.com/whit2333/hallac_evio
- `analyzer`: Hall A analyzer (podd)from https://github.com/whit2333/analyzer
- `hcana`: Hall C analyzer from https://github.com/whit2333/hcana
These are all built using the super build project `cool_halls` (https://github.com/whit2333/cool_halls)
containers/docker/deploy.env deleted 100644 → 0
View file @ 36a71670
# You have to define the values in {}
#DOCKER_REPO={account-nr}.dkr.ecr.{region}.amazonaws.com
## optional aws-cli options
#AWS_CLI_PROFILE={aws-cli-profile}
#AWS_CLI_REGION={aws-cli-region}
containers/docker/usage.sh deleted 100644 → 0
View file @ 36a71670
# INSTALL
# - copy the files deploy.env, config.env, version.sh and Makefile to your repo
# - replace the vars in deploy.env
# - define the version script
# Build the container
make build
# Build and publish the container
make release
# Publish a container to AWS-ECR.
# This includes the login to the repo
make publish
# Run the container
make run
# Build an run the container
make up
# Stop the running container
make stop
# Build the container with differnt config and deploy file
make cnf=another_config.env dpl=another_deploy.env build
\ No newline at end of file
containers/singularity/Singularity deleted 100644 → 0
View file @ 36a71670
Bootstrap: docker
From: eicweb.phy.anl.gov:4567/jlab/hallc/analyzer_software/hcana:latest
%help
Hall A/C container.
Tools:
- evio : EVIO DAQ data format
- analyzer : Hall A analyzer (podd)
- hcana : Hall C analyzer (hcana)
- root : root version used for the analyzer
- rootls, rootbrowse, root_config
%labels
Maintainer "Whitney Armstrong, Sylvester Joosten"
Version v1.0
%setup -c /bin/bash
export SINGULARITY_SHELL=/bin/bash
%environment -c /bin/bash
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include:/usr/local/include/podd:/usr/local/include/hcana
%post -c /bin/bash
echo "Hello from post"
echo "Install additional software here"
source /usr/local/bin/thisroot.sh
## libformat and nlohmann json used heavily in new generation replay scripts
## libformat
#cd /tmp && git clone https://github.com/fmtlib/fmt.git && cd fmt && \
# git checkout 5.3.0 && mkdir /tmp/build && cd /tmp/build && \
# cmake -DBUILD_SHARED_LIBS=TRUE ../fmt &&
# make -j4 install && cd /tmp && rm -r /tmp/build && rm -r /tmp/fmt
### json
# =======================
# global
# =======================
%runscript
echo "Launching a shell in the Hall A/C singularity container
exec bash
# =======================
# root
# =======================
%apprun root
root "$@"
%appenv root
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
# =======================
# analyzer
# =======================
%apprun analyzer
analyzer "$@"
%appenv analyzer
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
# =======================
# hcana
# =======================
%apphelp hcana
Run the Hall-C analyzer with same root-style arguments.
%apprun hcana
source /usr/local/bin/thisroot.sh
hcana "$@"
%appenv hcana
export DB_DIR=DBASE
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOTSYS=/usr/local
export ROOT_INCLUDE_PATH=/usr/local/include
export ROOT_INCLUDE_PATH=/usr/local/include:/usr/local/include/podd:/usr/local/include/hcana
# =======================
# root-config
# =======================
%apprun root_config
root-config "$@"
%appenv root_config
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
# =======================
# rootbrowse
# =======================
%apprun rootbrowse
rootbrowse "$@"
%appenv rootbrowse
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
# =======================
# rootls
# =======================
%apprun rootls
rootls "$@"
%appenv rootls
export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
export PATH=/usr/local/bin:${PATH}
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
podd @ 355b3ce6
Compare ecff1455...355b3ce6
Subproject commit ecff1455f2868eee65295d4aa6e772eefb2dd941
Subproject commit 355b3ce629f8f47f107e2a0fcfa7ee55b3463b6c
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
......@@ -28,16 +28,16 @@ using namespace std;
//_____________________________________________________________________________
THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus* app)
: 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()
: 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 )
......@@ -90,6 +90,11 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
fPeriodCheck = 0.0;
fCycle = 0.0;
fglHelicityScaler = 0;
fHelicityHistory = 0;
fRecommendedFreq = -1.0;
fStatus = kOK;
return fStatus;
}
......@@ -106,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 )
......@@ -120,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},
......@@ -157,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;
}
......@@ -171,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;
......@@ -197,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;
}
......@@ -266,7 +269,7 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
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;
}
......@@ -274,20 +277,6 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
fMPS = fIsMPS ? 1 : 0;
fQrt = fIsQrt ? 1 : 0; // Last of quartet
if (fglHelicityScaler) {
Int_t nhelev = fglHelicityScaler->GetNevents();
Int_t ncycles = fglHelicityScaler->GetNcycles();
if (nhelev > 0) {
for (Int_t i = 0; i < nhelev; i++) {
fScaleQuartet = (fHelicityHistory[i] & 2) != 0;
Int_t ispos = fHelicityHistory[i] & 1;
if (fScaleQuartet) {
fScalerSeed = ((fScalerSeed << 1) | ispos) & 0x3FFFFFFF;
}
}
}
}
if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
return 0;
......@@ -318,11 +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 / fTIPeriod, 1.0);
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;
......@@ -337,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
......@@ -349,14 +345,14 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
Int_t quartets_missed = (newNCycle - fFirstCycle) / 4 - (fNCycle - 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;
......@@ -389,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;
......@@ -422,12 +418,12 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
if ((fNCycle - fFirstCycle) % 4 != 0) { // Test if first in a quartet
fNQRTProblems++;
if (fNQRTProblems > 10) {
cout << "QRT Problem resetting" << endl;
_logger->warn("QRT Problem resetting");
fHaveQRT = kFALSE;
fFoundQuartet = kFALSE;
fNQRTProblems = 0;
} else {
cout << "Ignored " << fNQRTProblems << " problems" << endl;
_logger->warn("Ignored {} problems", fNQRTProblems);
}
} else {
fNQRTProblems = 0;
......@@ -442,8 +438,8 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
fNLastQuartet = fNQuartet - 1; // Make sure LoadHelicity uses
fNBits = 0;
fNQRTProblems = 0;
cout << "Phase found from QRT signal" << endl;
cout << "fFirstcycle = " << fFirstCycle << endl;
_logger->info("Phase found from QRT signal");
_logger->info("fFirstCycle = {}", fFirstCycle);
}
} else {
if (fHaveQRT) { // Using qrt signal.
......@@ -452,12 +448,12 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
if ((fNCycle - fFirstCycle) % 4 == 0) { // Shouldn't happen
fNQRTProblems++;
if (fNQRTProblems > 10) {
cout << "Shouldn't happen, cycle=" << fNCycle << "/" << fFirstCycle << endl;
_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 {
cout << "Ignored " << fNQRTProblems << " problems" << endl;
_logger->warn("Ignored {} problems", fNQRTProblems);
}
}
} else { // Presumable pre qrt signal data
......@@ -465,25 +461,22 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
if (!fFoundQuartet) {
// fFirstCycle = fNCycle - 3;
cout << "Quartet potentially found, starting at cycle " << fFirstCycle << endl;
_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
cout << "Lost quartet sync at cycle " << fNCycle << endl;
cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " "
<< fQuartet[3] << endl;
_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;
cout << "Searching for first of a quartet at cycle "
<< " " << fFirstCycle << endl;
cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " "
<< fQuartet[3] << endl;
_logger->info("Searching for first of a quartet at cycle {}: {} {} {} {}",
fFirstCycle, fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
fFirstCycle++;
}
} else {
......@@ -506,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;
......@@ -528,12 +521,11 @@ 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;
......@@ -547,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.
......@@ -558,9 +566,10 @@ 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;
......@@ -568,8 +577,8 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
// if(!fFixFirstCycle) {
if (fNQuartet - fNLastQuartet > 1) { // If we missed a quartet
if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over
cout << "fNBits = 0, missedcycles=" << missedcycles << " " << fNLastQuartet << " "
<< fNQuartet << endl;
_logger->warn("fNBits = 0, missedcycles={} {} {}", missedcycles, fNLastQuartet,
fNQuartet);
fNBits = 0;
return;
}
......@@ -586,7 +595,7 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
if (fNQuartet != fNLastQuartet) {
// Sanity check fNQuartet == fNLastQuartet+1
if (fNBits == 0) {
cout << "Start calibrating at cycle " << cyclecount << endl;
_logger->info("Start calibrating at cycle {}", cyclecount);
fRingSeed_reported = 0;
}
// Do phase stuff right here
......@@ -602,20 +611,39 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
if (fReportedHelicity == kUnknown) {
fNBits = 0;
fRingSeed_reported = 0;
} else if (fNBits == fMAXBIT) {
cout << "Seed Found " << bitset<32>(fRingSeed_reported) << " at cycle " << cyclecount
<< " with first cycle " << fFirstCycle << endl;
} else if (fNBits == fMAXBIT + 0) {
_logger->info("Seed Found {:32b} at cycle {} with first cycle {}", fRingSeed_reported,
cyclecount, fFirstCycle);
if (fglHelicityScaler) {
cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
_logger->info("Scaler Seed {:32b}", fScalerSeed);
}
Int_t backseed = GetSeed30(fRingSeed_reported);
cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec << endl;
_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);
}
}
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 {
......@@ -633,9 +661,9 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
// fRingSeed_reported << " " << fRingSeed_actual << dec
//<< endl;
if (fReportedHelicity != fPredictedHelicity) {
cout << "Helicity prediction failed " << fReportedHelicity << " " << fPredictedHelicity
<< " " << fActualHelicity << endl;
cout << bitset<32>(fRingSeed_reported) << " " << bitset<32>(fRingSeed_actual) << endl;
_logger->warn("Helicity prediction failed {} {} {}", fReportedHelicity, fPredictedHelicity,
fActualHelicity);
_logger->warn("{} {}", fRingSeed_reported, fRingSeed_actual);
fNBits = 0; // Need to reaquire seed
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
......@@ -670,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,12 +6,13 @@
// 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;
......@@ -19,113 +20,110 @@ class THcHelicityScaler;
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 void SetHelicityScaler(THcHelicityScaler *f);
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;
Bool_t fFixFirstCycle;
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 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 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;
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;
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;
Int_t* fHelicityHistory;
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
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
......@@ -23,8 +23,8 @@ using namespace std;
//____________________________________________________________________
THcHelicityReader::THcHelicityReader()
: hcana::ConfigLogging<podd2::EmptyBase>(),
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
}
//____________________________________________________________________
......@@ -45,9 +45,10 @@ Int_t THcHelicityReader::ReadDatabase(const char* /*dbfilename*/, const char* /*
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
cout << "THcHelicityReader: Helicity information from ROC 2 (SHMS)" << endl;
_logger->info("{}: Helicity information from ROC 2 (SHMS)", here);
SetROCinfo(kHel, 2, 14, 9);
SetROCinfo(kHelm, 2, 14, 8);
SetROCinfo(kMPS, 2, 14, 10);
......@@ -104,9 +105,9 @@ Int_t THcHelicityReader::ReadData(const THaEvData& evdata) {
}
// Check if ROC info is correct
if (!evdata.GetModule(fROCinfo[kTime].roc, fROCinfo[kTime].slot)) {
cout << "THcHelicityReader: ROC 2 not found" << endl;
cout << "Changing to ROC 1 (HMS)" << endl;
_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);
......@@ -119,6 +120,17 @@ Int_t THcHelicityReader::ReadData(const THaEvData& evdata) {
// 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);
// 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++;
......@@ -177,3 +189,4 @@ Int_t THcHelicityReader::SetROCinfo(EROC which, Int_t roc, Int_t slot, Int_t ind
//____________________________________________________________________
ClassImp(THcHelicityReader)
src/THcHelicityScaler.cxx
View file @ e974befb
......@@ -3,109 +3,322 @@
\brief Event handler for Hall C helicity scalers
Scalers not yet implemented. For now just picks the helicity control
bits out of the scaler words
~~~
gHaEvtHandlers->Add (new THcHelicityScaler("H","HC helicity scalers"));
~~~
To enable debugging you may try this in the setup script
~~~
THcHelcityScaler *hhelscaler = new THcHelicityScaler("H","HC helicity scalers");
hscaler->SetDebugFile("HHelScaler.txt");
gHaEvtHandlers->Add (hhelscaler);
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);
~~~
\author
\authors: S. A. Wood (saw@jlab.org)
C. Yero (cyero@jlab.org)
*/
#include "THaEvtTypeHandler.h"
#include "THcHelicityScaler.h"
#include "THaCodaData.h"
#include "THaEvData.h"
#include "THcParmList.h"
#include "THcGlobals.h"
#include "THaGlobals.h"
#include "THcHelicity.h"
#include "TNamed.h"
#include "TMath.h"
#include "TString.h"
#include "TROOT.h"
#include <cstring>
#include <bitset>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <map>
#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"
#include "Helper.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 ICURRENT = 3;
static const UInt_t ICHARGE = 4;
static const UInt_t ITIME = 5;
static const UInt_t ICUT = 6;
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)
: THaEvtTypeHandler(name,description),
fBCM_Gain(0), fBCM_Offset(0), fBCM_delta_charge(0),
fBankID(9801),
evcount(0), evcountR(0.0), ifound(0),
fUseFirstEvent(kTRUE),
fOnlySyncEvents(kFALSE), fOnlyBanks(kFALSE), fDelayedType(-1),
fClockChan(-1), fLastClock(0), fClockOverflows(0)
{
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()
{
delete [] fBCM_Gain;
delete [] fBCM_Offset;
delete [] fBCM_delta_charge;
THcHelicityScaler::~THcHelicityScaler() {
if (!TROOT::Initialized()) {
delete fScalerTree;
}
Podd::DeleteContainer(scalers);
Podd::DeleteContainer(scalerloc);
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
for (auto it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
}
Int_t THcHelicityScaler::End( THaRunBase* )
{
Int_t THcHelicityScaler::End(THaRunBase* runbase) {
// Process any delayed events in order received
cout << "THcHelicityScaler::End Analyzing " << fDelayedEvents.size() << " delayed scaler events" << endl;
for(std::vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it) {
for (std::vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end();
++it) {
UInt_t* rdata = *it;
AnalyzeBuffer(rdata,kFALSE);
AnalyzeBuffer(rdata);
}
evNumberR = -1;
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
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) {
Int_t THcHelicityScaler::ReadDatabase(const TDatime& date )
{
char prefix[2];
prefix[0]='g';
prefix[1]='\0';
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.
*
......@@ -116,67 +329,130 @@ void THcHelicityScaler::SetDelayedType(int evtype) {
*/
fDelayedType = evtype;
}
Int_t THcHelicityScaler::Analyze(THaEvData *evdata)
{
if ( !IsMyEvent(evdata->GetEvType()) ) return -1;
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;
*fDebugFile << endl << "---------------------------------- " << endl << endl;
*fDebugFile << "\nEnter THcHelicityScaler for fName = " << fName << endl;
EvDump(evdata);
}
UInt_t *rdata = (UInt_t*) evdata->GetRawDataBuffer();
if(evdata->GetEvType() == fDelayedType) { // Save this event for processing later
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];
UInt_t* datacopy = new UInt_t[evlen];
fDelayedEvents.push_back(datacopy);
memcpy(datacopy,rdata,evlen*sizeof(UInt_t));
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();
}
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,fOnlySyncEvents))) {
//
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, Bool_t onlysync)
{
fNevents = 0;
Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata) {
fNTrigsInBuf = 0;
// Parse the data, load local data arrays.
UInt_t *p = (UInt_t*) rdata;
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;
UInt_t* plast = p + *p; // Index to last word in the bank
Int_t roc = -1;
Int_t evlen = *p + 1;
ifound=0;
while(p<plast) {
Int_t ifound = 0;
while (p < plast) {
Int_t banklen = *p;
p++; // point to header
p++; // point to header
if (fDebugFile) {
*fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p-1) << endl;
*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(fRocSet.find(roc)==fRocSet.end()) { // Not a ROC with helicity scaler
p+=*(p-1)-1; // Skip to next ROC
}
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
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
......@@ -185,122 +461,991 @@ Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
// 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
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;
// cout << "BankID=" << tag << endl;
if (tag != fBankID) {
p = pnext; // Fall through to end of the above else if
// cout << " Skipping to next bank" << endl;
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 == 5) {
Int_t nevents = (banklen-2)/32;
//cout << "# of helicity events in bank:" << " " << nevents << endl;
if (nevents > 100) {
cout << "Error! Beam off for too long" << endl;
}
fNevents = nevents;
fNcycles += nevents;
// 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 nentries = 32*iev+2;
fHelicityHistory[iev] = (p[nentries-1]>>30) & 0x3;
// cout << "H: " << evNumber << endl;
}
}
// 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;
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
p = p + *(p - 1); // Skip to next bank
}
}
if (fDebugFile) {
*fDebugFile << "Finished with decoding. "<<endl;
*fDebugFile << " Found flag = "<<ifound<<endl;
*fDebugFile << "Finished with decoding. " << endl;
*fDebugFile << " Found flag = " << ifound << endl;
}
// HMS has headers which are different from SOS, but both are
// event type 0 and come here. If you found no headers, return.
if (!ifound) return 0;
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::AnalyzeHelicityScaler(UInt_t *p)
//{
//}
//_____________________________________________________________________________
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;
THaAnalysisObject::EStatus THcHelicityScaler::Init(const TDatime& date)
{
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;
fStatus = kOK;
fNormIdx = -1;
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
cout << "Howdy ! We are initializing THcHelicityScaler !! name = "
<< fName << endl;
_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;
if(eventtypes.size()==0) {
eventtypes.push_back(0); // Default Event Type
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]);
}
}
fRocSet.insert(5); // List ROCs that have helicity scalers
fRocSet.insert(8); // Should make configurable
// 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();
fNcycles = 0;
fNevents = 0;
return kOK;
}
size_t THcHelicityScaler::FindNoCase(const string& sdata, const string& skey)
{
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 = "";
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) {
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);
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
View file @ e974befb
......@@ -9,80 +9,171 @@
/////////////////////////////////////////////////////////////////////
#include "THaEvtTypeHandler.h"
#include "THcScalerEvtHandler.h"
#include "Decoder.h"
#include <string>
#include <vector>
#include <algorithm>
#include <set>
#include "TTree.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 THaEvtTypeHandler {
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, Bool_t onlysync);
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 SetHelicityDetector(THcHelicity *f) {fglHelicityDetector = f;}
virtual Int_t GetNevents() { return fNevents;}
virtual Int_t GetNcycles() { return fNcycles;}
virtual Int_t GetEvNum() { return evNumber;}
virtual Int_t* GetHelicityHistoryP() {return fHelicityHistory;}
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:
static size_t FindNoCase(const std::string& sdata, const std::string& skey);
Int_t fNumBCMs;
Double_t *fBCM_Gain;
Double_t *fBCM_Offset;
Double_t *fBCM_delta_charge;
Int_t fROC;
UInt_t fBankID;
// Helicity Scaler variables
Int_t fNevents; /* # of helicity scaler reads in last event */
Int_t fNcycles;
Int_t fHelicityHistory[200];
//
Bool_t fUseFirstEvent;
Bool_t fOnlySyncEvents;
Bool_t fOnlyBanks;
Int_t fDelayedType;
Int_t fClockChan;
UInt_t fLastClock;
Int_t fClockOverflows;
std::vector<UInt_t*> fDelayedEvents;
std::set<UInt_t> fRocSet;
THcHelicityScaler(const THcHelicityScaler& fh);
THcHelicityScaler& operator=(const THcHelicityScaler& fh);
UInt_t evcount;
Double_t evcountR;
UInt_t evNumber;
Int_t ifound;
THcHelicity *fglHelicityDetector;
ClassDef(THcHelicityScaler,0) // Scaler Event handler
//------------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
......