-
Sylvester Joosten authoredmerged in changes from upstream, including major changes to the THcHodoscope class. I had to manually propagate some of these changes because of the differences in implementation for this classs.
Sylvester Joosten authoredmerged in changes from upstream, including major changes to the THcHodoscope class. I had to manually propagate some of these changes because of the differences in implementation for this classs.
THcHelicity.cxx 22.82 KiB
//*-- Author : Julie Roche, November 2010
// this is a modified version of THaGOHelicity.C
////////////////////////////////////////////////////////////////////////
//
// THcHelicity
//
// Helicity of the beam from QWEAK electronics in delayed mode
// +1 = plus, -1 = minus, 0 = unknown
//
// Also supports in-time mode with delay = 0
//
////////////////////////////////////////////////////////////////////////
#include "THcHelicity.h"
#include "TH1F.h"
#include "THaApparatus.h"
#include "THaEvData.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "TMath.h"
#include <iostream>
using namespace std;
//_____________________________________________________________________________
THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus* app)
: hcana::ConfigLogging<THaHelicityDet>(name, description, app), fnQrt(-1), fHelDelay(8),
fMAXBIT(30) {
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i] = 0;
// memset(fHbits, 0, sizeof(fHbits));
}
//_____________________________________________________________________________
THcHelicity::THcHelicity()
: hcana::ConfigLogging<THaHelicityDet>(), fnQrt(-1), fHelDelay(8), fMAXBIT(30) {
// Default constructor for ROOT I/O
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i] = 0;
}
//_____________________________________________________________________________
THcHelicity::~THcHelicity() {
DefineVariables(kDelete);
// for( Int_t i = 0; i < NHIST; ++i ) {
// delete fHisto[i];
// }
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
// Call `Setup` before everything else.
Setup(GetName(), GetTitle());
fFirstEvProcessed = kFALSE;
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
fLastMPSTime = 0;
fFoundMPS = kFALSE;
// Call initializer for base class.
// This also calls `ReadDatabase` and `DefineVariables`.
EStatus status = THaDetector::Init(date);
if (status) {
fStatus = status;
return fStatus;
}
fEvNumCheck = 0;
fDisabled = kFALSE;
fStatus = kOK;
return fStatus;
}
//_____________________________________________________________________________
void THcHelicity::Setup(const char* name, const char* description) {
// Prefix for parameters in `param` file.
string kwPrefix = string(GetApparatus()->GetName()) + "_" + name;
std::transform(kwPrefix.begin(), kwPrefix.end(), kwPrefix.begin(), ::tolower);
fKwPrefix = kwPrefix;
}
//_____________________________________________________________________________
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 )
// return st;
// Read readout parameters (ROC addresses etc.)
Int_t st = THcHelicityReader::ReadDatabase(GetDBFileName(), GetPrefix(), date, fQWEAKDebug);
if (st != kOK)
return st;
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;
DBRequest list[] = {// {"_hsign", &fSign, kInt, 0, 1},
{"helicity_delay", &fHelDelay, kInt, 0, 1},
{"helicity_freq", &fFreq, kDouble, 0, 1},
// {"helicity_seed", &fRingSeed_reported_initial, kInt, 0, 1},
// {"helicity_cycle", &fFirstCycle, kInt, 0, 1},
{0}};
gHcParms->LoadParmValues(list, "");
fMAXBIT = 30;
fTIPeriod = 250000000.0 / fFreq;
// maximum of event in the pattern, for now we are working with quartets
// Int_t localpattern[4]={1,-1,-1,1};
// careful, the first value here should always +1
// for(int i=0;i<fQWEAKNPattern;i++)
// {
// fPatternSequence.push_back(localpattern[i]);
// }
HWPIN = kTRUE;
fQuartet[0] = fQuartet[1] = fQuartet[2] = fQuartet[3] = 0;
if (fFirstCycle >= 0 && fRingSeed_reported_initial != 0) {
// Set the seed for predicted reported and predicted actual
} else {
// Initialize mode to find quartets and then seed
}
_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;
}
//_____________________________________________________________________________
void THcHelicity::MakePrefix() { THaDetector::MakePrefix(); }
//_____________________________________________________________________________
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;
fIsSetup = (mode == kDefine);
// Define standard variables from base class
THaHelicityDet::DefineVariables(mode);
const RVarDef var[] = {{"nqrt", "position of cycle in quartet", "fnQrt"},
{"hel", "actual helicity for event", "fActualHelicity"},
{"helrep", "reported helicity for event", "fReportedHelicity"},
{"helpred", "predicted reported helicity for event", "fPredictedHelicity"},
{"mps", "In MPS blanking period", "fMPS"},
{0}};
// cout << "Calling THcHelicity DefineVarsFromList" << endl;
_logger->info("Calling THcHelicity DefineVarsFromList");
return DefineVarsFromList(var, mode);
}
//_____________________________________________________________________________
void THcHelicity::PrintEvent(Int_t evtnum) {
cout << " ++++++ THcHelicity::Print ++++++\n";
cout << " +++++++++++++++++++++++++++++++++++++\n";
return;
}
//_____________________________________________________________________________
Int_t THcHelicity::Begin(THaRunBase*) {
THcHelicityReader::Begin();
// fHisto[0] = new TH1F("hel.seed","hel.seed",32,-1.5,30.5);
// fHisto[1] = new TH1F("hel.error.code","hel.error.code",35,-1.5,33.5);
return 0;
}
//_____________________________________________________________________________
// void THcHelicity::FillHisto()
//{
// fHisto[0]->Fill(fRing_NSeed);
// fHisto[1]->Fill(fErrorCode);
// return;
//}
//_____________________________________________________________________________
void THcHelicity::SetErrorCode(Int_t error) {
// used as a control for the helciity computation
// 2^0: if the reported number of events in a pattern is larger than fQWEAKNPattern
// 2^1: if the offset between the ring reported value and TIR value is not fOffsetTIRvsRing
// 2^2: if the reported time in the ring is 0
// 2^3: if the predicted reported helicity doesn't match the reported helicity in the ring // 2^4: if the helicity cannot be computed using the SetHelicity routine
// 2^5: if seed is being gathered
if (fErrorCode == 0)
fErrorCode = (1 << error);
// only one reported error at the time
return;
}
//_____________________________________________________________________________
void THcHelicity::Clear(Option_t* opt) {
// Clear event-by-event data
THaHelicityDet::Clear(opt);
THcHelicityReader::Clear(opt);
fEvtype = 0;
fQrt = 0;
fErrorCode = 0;
return;
}
//_____________________________________________________________________________
Int_t THcHelicity::Decode(const THaEvData& evdata) {
// Decode Helicity data.
// Return 1 if helicity was assigned, 0 if not, <0 if error.
Int_t err = ReadData(evdata); // from THcHelicityReader class
if (err) {
Error(Here("THcHelicity::Decode"), "Error decoding helicity data.");
return err;
}
fReportedHelicity = (fIsHelp ? (fIsHelm ? kUnknown : kPlus) : (fIsHelm ? kMinus : kUnknown));
fMPS = fIsMPS ? 1 : 0;
if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
return 0;
}
if (fDisabled) {
fActualHelicity = kUnknown;
return 0;
}
fEvNumCheck++;
Int_t evnum = evdata.GetEvNum();
if (fEvNumCheck != evnum) {
_logger->info("THcHelicity: Missed {} events at event {}.", evnum - fEvNumCheck, evnum);
_logger->info(" Disabling helicity decoding for rest of run.");
_logger->info(
" Make sure \"RawDecode_master in cuts file accepts all physics events.");
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
fLastMPSTime = (fTITime + fLastMPSTime + missed * fTIPeriod) / 2;
fIsNewCycle = kTRUE;
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
} else {
fLastMPSTime = (fLastMPSTime + fTITime - missed * fTIPeriod) / 2;
}
// If there is a skip, pass it off to next non MPS event
// Need to also check here for missed MPS's
// cout << "Found MPS" << endl;
// check for Nint((time-last)/period) > 1
} else {
fFoundMPS = kTRUE;
fLastMPSTime = fTITime;
}
} else if (fFoundMPS) { //
if (fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods
missed = TMath::Nint(floor((fTITime - fLastMPSTime) / fTIPeriod));
if (missed > 1) {
// cout << "Missed " << missed << " MPSes" << endl;
Int_t newNCycle = fNCycle + missed - 1; // How many cycles really missed
Int_t quartets_missed = (newNCycle - fFirstCycle) / 4 - (fNCycle - fFirstCycle) / 4;
for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
int quartetphase = (newNCycle - fFirstCycle) % 4;
// cout << " " << 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) {
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];
} else {
fQuartet[0] = fQuartet[1] = -fReportedHelicity;
fQuartet[2] = -fQuartet[1];
}
} else {
fActualHelicity = kUnknown;
fQuartet[0] = fReportedHelicity;
fQuartet[1] = 0;
}
}
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];
fQuartet[0] = fReportedHelicity;
fNCycle++;
if ((fNCycle - fFirstCycle) % 4 == 3) { // Test if last in a quartet
if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
if (!fFoundQuartet) {
// fFirstCycle = fNCycle - 3;
_logger->info("Quartet potentially found, starting at cycle {} - event {}",
fFirstCycle, evdata.GetEvNum());
// cout << "Quartet potentially found, starting at cycle " << fFirstCycle << " - event
// "
// << evdata.GetEvNum() << endl;
fFoundQuartet = kTRUE;
}
} else {
if (fNCycle - fFirstCycle > 4) { // Not at start of run. Reset
_logger->warn("Lost quartet sync at cycle {} - event {}", fNCycle, evdata.GetEvNum());
_logger->warn("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
// cout << "Lost quartet sync at cycle " << fNCycle << " - event " <<
// evdata.GetEvNum()
// << endl;
// cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " <<
// fQuartet[3]
// << endl;
fFirstCycle += 4 * ((fNCycle - fFirstCycle) / 4); // Update, but don't change phase
}
fFoundQuartet = kFALSE;
fNBits = 0;
_logger->info("Searching for first of a quartet at cycle {} - event {}", fFirstCycle,
evdata.GetEvNum());
// cout << "Searching for first of a quartet at cycle "
// << " " << fFirstCycle << " - event " << evdata.GetEvNum() << endl;
// cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " << fQuartet[3]
// << endl;
_logger->info("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
fFirstCycle++;
}
}
// Load the actual helicity. Calibrate if not calibrated.
fActualHelicity = kUnknown;
LoadHelicity(fReportedHelicity, fNCycle, missed);
fLastReportedHelicity = fReportedHelicity;
fIsNewCycle = kFALSE;
// cout << fTITime/250000000.0 << " " << fNCycle << " " << fReportedHelicity << endl;
// cout << fNCycle << ": " << fReportedHelicity << " "
// << 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;
fPredictedHelicity = kUnknown;
fFirstEvTime = fTITime;
fLastEvTime = fTITime; fLastMPSTime = fTITime; // Not necessarily during the MPS
fNCycle = 0;
fFirstEvProcessed = kTRUE;
fFoundMPS = kFALSE;
fFoundQuartet = kFALSE;
fIsNewCycle = kFALSE;
fNBits = 0;
}
// Some sanity checks
if (fActualHelicity < -5) {
_logger->info("Actual Helicity never got defined");
}
if (fNBits < fMAXBIT) {
if (fActualHelicity == -1 || fActualHelicity == 1) {
cout << "Helicity of " << fActualHelicity << " reported prematurely at cycle " << fNCycle
<< endl;
}
}
fLastActualHelicity = fActualHelicity;
return 0;
}
//_____________________________________________________________________________
Int_t THcHelicity::End(THaRunBase*) {
// End of run processing. Write histograms.
THcHelicityReader::End();
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i]->Write();
return 0;
}
//_____________________________________________________________________________
void THcHelicity::SetDebug(Int_t level) {
// Set debug level of this detector as well as the THcHelicityReader
// helper class.
THaHelicityDet::SetDebug(level);
fQWEAKDebug = level;
}
//_____________________________________________________________________________
void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles) {
// static const char* const here = "THcHelicity::LoadHelicity";
int quartetphase = (cyclecount - fFirstCycle) % 4;
fnQrt = quartetphase;
if (missedcycles > 1) { // If we missed windows
if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over
fNBits = 0;
return;
}
}
if (!fFoundQuartet) { // Wait until we have found quad phase before starting
return; // to calibrate
}
if (quartetphase == 0) { // Start of a quad
if (fNBits < fMAXBIT) {
if (fNBits == 0) {
_logger->info("Start calibrating at cycle {}", cyclecount);
// cout << "Start calibrating at cycle " << cyclecount << endl;
fRingSeed_reported = 0;
}
if (fReportedHelicity == kPlus) {
fRingSeed_reported = ((fRingSeed_reported << 1) | 1) & 0x3FFFFFFF;
} else {
fRingSeed_reported = (fRingSeed_reported << 1) & 0x3FFFFFFF;
} fNBits++;
if (fReportedHelicity == kUnknown) {
fNBits = 0;
fRingSeed_reported = 0;
} else if (fNBits == fMAXBIT) {
_logger->info("Seed Found {} at cycle {} with first cycle {}", fRingSeed_reported,
cyclecount, fFirstCycle);
// cout << "Seed Found " << hex << fRingSeed_reported << dec << " at cycle " << cyclecount
// << " with first cycle " << fFirstCycle << endl;
Int_t backseed = GetSeed30(fRingSeed_reported);
_logger->info("Seed at cycle {} should be {}", fFirstCycle, backseed);
// cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec <<
// endl;
}
fActualHelicity = kUnknown;
} else if (fNBits >= fMAXBIT) {
fRingSeed_reported = RanBit30(fRingSeed_reported);
if (fNBits == fMAXBIT) {
fRingSeed_actual = fRingSeed_reported;
for (Int_t i = 0; i < fHelDelay / 4; i++) {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fNBits++;
} else {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fActualHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
// if(fTITime/250000000.0 > 380.0) cout << fTITime/250000000.0 << " " << fNCycle << " "
// << hex <<
// fRingSeed_reported << " " << fRingSeed_actual << dec
//<< endl;
if (fReportedHelicity != fPredictedHelicity) {
_logger->warn("Helicity prediction failed {} {} {}", fReportedHelicity, fPredictedHelicity,
fActualHelicity);
// cout << "Helicity prediction failed " << fReportedHelicity << " "
// << fPredictedHelicity << " " << fActualHelicity << endl;
// cout << hex << fRingSeed_reported << " " << fRingSeed_actual << dec << endl;
fNBits = 0; // Need to reaquire seed
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
}
}
fQuartetStartHelicity = fActualHelicity;
fQuartetStartPredictedHelicity = fPredictedHelicity;
} else { // Not the beginning of a quad
if (fNBits >= fMAXBIT) {
fActualHelicity =
(quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity : -fQuartetStartHelicity;
fPredictedHelicity = (quartetphase == 0 || quartetphase == 3)
? fQuartetStartPredictedHelicity
: -fQuartetStartPredictedHelicity;
}
}
return;
}
//_____________________________________________________________________________
Int_t THcHelicity::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;
if (ranseed <= 0) {
if (fQWEAKDebug > 1)
std::cerr << "ranseed must be greater than zero!" << "\n";
newbit = 0;
}
ranseed = ((ranseed << 1) | newbit) & 0x3FFFFFFF;
// here ranseed is changed
if (fQWEAKDebug > 1) {
cout << "THcHelicity::RanBit30, newbit=" << newbit << "\n";
}
return ranseed;
}
//_____________________________________________________________________________
Int_t THcHelicity::GetSeed30(Int_t currentseed)
/* Back track the seed by 30 samples */
{
#if 1
Int_t seed = currentseed;
for (Int_t i = 0; i < 30; i++) {
UInt_t bit1 = (seed & 0x00000001) != 0;
UInt_t bit8 = (seed & 0x00000080) != 0;
UInt_t bit29 = (seed & 0x10000000) != 0;
UInt_t bit30 = (seed & 0x20000000) != 0;
UInt_t newbit30 = (bit30 ^ bit29 ^ bit8 ^ bit1) & 0x1;
seed = (seed >> 1) | (newbit30 << 29);
}
#else
Int_t bits = currentseed;
Int_t seed = 0;
for (Int_t i = 0; i < 30; i++) {
Int_t val;
// XOR at virtual position 0 and 29
if (i == 0) {
val = ((bits & (1 << (i))) != 0) ^ ((bits & (1 << (i + 29))) != 0);
} else {
val = ((bits & (1 << (i))) != 0) ^ ((seed & (1 << (i - 1))) != 0);
}
if (i <= 1) {
val = ((bits & (1 << (1 - i))) != 0) ^ val;
} else {
val = ((seed & (1 << (i - 2))) != 0) ^ val;
}
if (i <= 22) {
val = ((bits & (1 << (i - 22))) != 0) ^ val;
} else {
val = ((seed & (1 << (i - 23))) != 0) ^ val;
}
seed |= (val << i);
}
#endif
return seed;
}
ClassImp(THcHelicity)