From 76bbcc81ebf0a842f1365817780e7d409ddf1415 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Thu, 11 Apr 2019 11:58:37 -0400
Subject: [PATCH] Improve THcHelicity Previously the helicity code was not
always setting the actual helicity. It would then just use the helicity
from the last event. This was almost always OK. With this change, the
helicity (T.helicity.hel) should be always defined. It is now sufficient
to use just this variable when sorting events into plus and minus helicity.
The helicity value will be zero if the event occured during the settle
period, or if the helicity decoder has not got enough data to have a seed.
---
src/THcHelicity.cxx | 72 ++++++++++++++++++++++++++++++++++++++-------
src/THcHelicity.h | 3 ++
2 files changed, 64 insertions(+), 11 deletions(-)
diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx
index 2bb7418..6cfdd36 100644
--- a/src/THcHelicity.cxx
+++ b/src/THcHelicity.cxx
@@ -73,6 +73,9 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
return fStatus;
}
+ fEvNumCheck = 0;
+ fDisabled = kFALSE;
+
fStatus = kOK;
return fStatus;
}
@@ -140,6 +143,9 @@ Int_t THcHelicity::ReadDatabase( const TDatime& date )
// Initialize mode to find quartets and then seed
}
+ cout << "Helicity decoder initialized with frequency of " << fFreq
+ << " Hz and reporting delay of " << fHelDelay << " cycles." << endl;
+
return kOK;
}
@@ -256,11 +262,32 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
fActualHelicity = fIsMPS?kUnknown:fReportedHelicity;
return 0;
}
- Long64_t lastlastmpstime = fLastMPSTime;
+
+ if(fDisabled) {
+ fActualHelicity = kUnknown;
+ return 0;
+ }
+
+ fEvNumCheck++;
+ Int_t evnum = evdata.GetEvNum();
+ if(fEvNumCheck!=evnum) {
+ cout << "THcHelicity: Missed " << evnum-fEvNumCheck << " events at event " << evnum << endl;
+ cout << " Disabling helicity decoding for rest of run." << endl;
+ cout << " Make sure \"RawDecode_master in cuts file accepts all physics events." <<endl;
+ fDisabled = kTRUE;
+ fActualHelicity = kUnknown;
+ return 0;
+ }
+
+ fActualHelicity = -10.0;
if(fFirstEvProcessed) { // Normal processing
+ // cout << evnum << " " << fNCycle << " " << fIsMPS << " " << fFoundMPS << " " << fTITime << " "
+ // << fLastMPSTime << " " << fNBits << endl;
Int_t missed = 0;
// Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0;
if(fIsMPS) {
+ fActualHelicity = kUnknown;
+ fPredictedHelicity = kUnknown;
if(fFoundMPS) {
missed = TMath::Nint(fTITime/fTIPeriod-fLastMPSTime/fTIPeriod);
if(missed < 1) { // was <=1
@@ -292,12 +319,6 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
}
int quartetphase = (newNCycle-fFirstCycle)%4;
// cout << " " << fNCycle << " " << newNCycle << " " << fFirstCycle << " " << quartets_missed << " " << quartetphase << endl;
- fQuartetStartHelicity = (fRingSeed_actual&1)?kPlus:kMinus;
- fQuartetStartPredictedHelicity = (fRingSeed_reported&1)?kPlus:kMinus;
- fActualHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartHelicity:-fQuartetStartHelicity;
- fPredictedHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartPredictedHelicity:-fQuartetStartPredictedHelicity;
// cout << "Cycles " << fNCycle << " " << newNCycle << " " << fFirstCycle
// << " skipped " << quartets_missed << " quartets" << endl;
fNCycle = newNCycle;
@@ -305,6 +326,13 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
// reported helicity. So we don't fail quartet testing.
// But only do this if we are calibrated.
if(fNBits >= fMAXBIT) {
+ fQuartetStartHelicity = (fRingSeed_actual&1)?kPlus:kMinus;
+ fQuartetStartPredictedHelicity = (fRingSeed_reported&1)?kPlus:kMinus;
+ fActualHelicity = (quartetphase==0||quartetphase==3)?
+ fQuartetStartHelicity:-fQuartetStartHelicity;
+ fPredictedHelicity = (quartetphase==0||quartetphase==3)?
+ fQuartetStartPredictedHelicity:-fQuartetStartPredictedHelicity;
+
if (((fNCycle - fFirstCycle)%2)==1) {
fQuartet[0] = fReportedHelicity;
fQuartet[1] = fQuartet[2] = -fQuartet[0];
@@ -313,6 +341,7 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
fQuartet[2] = -fQuartet[1];
}
} else {
+ fActualHelicity = kUnknown;
fQuartet[0] = fReportedHelicity;
fQuartet[1] = 0;
}
@@ -320,6 +349,18 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
fLastMPSTime += missed*fTIPeriod;
fIsNewCycle = kTRUE;
fLastReportedHelicity = fReportedHelicity;
+ } else { // No missed periods. Get helicities from rings
+ if(fNBits>=fMAXBIT) {
+ int quartetphase = (fNCycle-fFirstCycle)%4;
+ fQuartetStartHelicity = (fRingSeed_actual&1)?kPlus:kMinus;
+ fQuartetStartPredictedHelicity = (fRingSeed_reported&1)?kPlus:kMinus;
+ fActualHelicity = (quartetphase==0||quartetphase==3)?
+ fQuartetStartHelicity:-fQuartetStartHelicity;
+ fPredictedHelicity = (quartetphase==0||quartetphase==3)?
+ fQuartetStartPredictedHelicity:-fQuartetStartPredictedHelicity;
+ } else {
+ fActualHelicity = 0;
+ }
}
if(fIsNewCycle) {
fQuartet[3]=fQuartet[2]; fQuartet[2]=fQuartet[1]; fQuartet[1]=fQuartet[0];
@@ -361,6 +402,8 @@ 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;
@@ -377,10 +420,17 @@ Int_t THcHelicity::Decode( const THaEvData& evdata )
fIsNewCycle = kFALSE;
fNBits = 0;
}
- // cout << setprecision(9) << "HEL " << fTITime/250000000.0 << " " << fNCycle << "(" << (fNCycle-fFirstCycle)%4 << "): "
- // << fMPS << " " << fReportedHelicity << " "
- // << fPredictedHelicity << " " << fActualHelicity << " " << lastlastmpstime/250000000.0 << endl;
- // bitset<32>(v)
+
+ // Some sanity checks
+ if(fActualHelicity < -5) {
+ cout << "Actual Helicity never got defined" << endl;
+ }
+ if(fNBits < fMAXBIT) {
+ if(fActualHelicity == -1 || fActualHelicity == 1) {
+ cout << "Helicity of " << fActualHelicity << " reported prematurely at cycle " << fNCycle << endl;
+ }
+ }
+ fLastActualHelicity = fActualHelicity;
return 0;
}
diff --git a/src/THcHelicity.h b/src/THcHelicity.h
index adc598e..cbfa7f5 100644
--- a/src/THcHelicity.h
+++ b/src/THcHelicity.h
@@ -94,6 +94,9 @@ protected:
Double_t fErrorCode;
Int_t fEvtype; // Current CODA event type
+ Int_t fLastActualHelicity;
+ Int_t fEvNumCheck;
+ Bool_t fDisabled;
static const Int_t NHIST = 2;
TH1F* fHisto[NHIST];
--
GitLab