////////////////////////////////////////////////////////////////////////
//
// THcHelicity
//
// Helicity of the beam in delayed mode using quartets
// +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 "THcHelicityScaler.h"
#include "THcParmList.h"
#include "TMath.h"
#include <bitset>
#include <iostream>
#include "THaHelicityDet.h"
#include "hcana/Logger.h"
using namespace std;
//_____________________________________________________________________________
THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus* app)
: 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()
: 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;
fLastHelpCycle = -1;
fQuadPattern[0] = fQuadPattern[1] = fQuadPattern[2] = fQuadPattern[3] = 0;
fQuadPattern[4] = fQuadPattern[5] = fQuadPattern[6] = fQuadPattern[7] = 0;
fHelperHistory = 0;
fHelperQuartetHistory = 0;
fScalerSeed = 0;
fNBits = 0;
lastispos = 0;
fLastReportedHelicity = kUnknown;
fFixFirstCycle = kFALSE;
fHaveQRT = kFALSE;
fNQRTProblems = 0;
fPeriodCheck = 0.0;
fCycle = 0.0;
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"},
{"pcheck", "Period check", "fPeriodCheck"},
{"cycle", "Helicity Cycle", "fCycle"},
{"qrt", "Last cycle of quartet", "fQrt"},
{0}};
_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;
}
//_____________________________________________________________________________
void THcHelicity::SetHelicityScaler(THcHelicityScaler* f) {
fglHelicityScaler = f;
fHelicityHistory = fglHelicityScaler->GetHelicityHistoryP();
}
//_____________________________________________________________________________
Int_t THcHelicity::Decode(const THaEvData& evdata) {
// Decode Helicity data.
// Return 1 if helicity was assigned, 0 if not, <0 if error.
evnum = evdata.GetEvNum();
Int_t err = ReadData(evdata); // from THcHelicityReader class
if (err) {
Error(Here("THcHelicity::Decode"), "Error decoding helicity data.");
return err;
}
fReportedHelicity = (fIsHelp ? (fIsHelm ? kUnknown : kPlus) : (fIsHelm ? kMinus : kUnknown));
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;
}
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) {
fPeriodCheck = fmod(fTITime / fTIPeriod, 1.0);
fCycle = (fTITime / fTIPeriod);
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;
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) {
for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
// cout << "Advancing seed A " << fNBits << endl;
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
fActualHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity
: -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) {
// cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
// Int_t predictedScalerSeed = RanBit30(fScalerSeed);
// cout << "Predct Seed " << bitset<32>(predictedScalerSeed) << endl;
// cout << "Event Seed " << bitset<32>(fRingSeed_reported) << endl;
fQuartet[3] = fQuartet[2];
fQuartet[2] = fQuartet[1];
fQuartet[1] = fQuartet[0];
fQuartet[0] = fReportedHelicity;
fNCycle++;
// cout << "XX: " << fNCycle << " " << fReportedHelicity << " " << fIsQrt <<endl;
if (fIsQrt) { //
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
if (fHaveQRT) { // Should already have phase set
if ((fNCycle - fFirstCycle) % 4 != 0) { // Test if first in a quartet
fNQRTProblems++;
if (fNQRTProblems > 10) {
cout << "QRT Problem resetting" << endl;
fHaveQRT = kFALSE;
fFoundQuartet = kFALSE;
fNQRTProblems = 0;
} else {
cout << "Ignored " << fNQRTProblems << " problems" << endl;
}
} else {
fNQRTProblems = 0;
}
}
if (!fHaveQRT) {
fHaveQRT = kTRUE;
fFoundQuartet = kTRUE;
fFirstCycle = fNCycle; //
fNQuartet = (fNCycle - fFirstCycle) / 4;
fNLastQuartet = fNQuartet - 1; // Make sure LoadHelicity uses
fNBits = 0;
fNQRTProblems = 0;
cout << "Phase found from QRT signal" << endl;
cout << "fFirstcycle = " << fFirstCycle << endl;
}
} else {
if (fHaveQRT) { // Using qrt signal.
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
if ((fNCycle - fFirstCycle) % 4 == 0) { // Shouldn't happen
fNQRTProblems++;
if (fNQRTProblems > 10) {
cout << "Shouldn't happen, cycle=" << fNCycle << "/" << fFirstCycle << endl;
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;
}
}
} else { // Presumable pre qrt signal data
if ((fNCycle - fFirstCycle) % 4 == 3) { // Test if last in a quartet
if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
if (!fFoundQuartet) {
// fFirstCycle = fNCycle - 3;
cout << "Quartet potentially found, starting at cycle " << fFirstCycle << endl;
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;
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;
fFirstCycle++;
}
} else {
fNLastQuartet = fNQuartet;
fNQuartet = (fNCycle - fFirstCycle) / 4;
}
}
}
// Load the actual helicity. Calibrate if not calibrated.
fActualHelicity = kUnknown;
// Here if we know we missed some earlier in the quartet, we need
// to make sure we get here and call LoadHelicity for the missing
// Cycles, reducing the missed count for each extra loadhelicity
// call that we make.
LoadHelicity(fReportedHelicity, fNCycle, missed);
fLastReportedHelicity = fReportedHelicity;
fIsNewCycle = kFALSE;
// 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(!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;
fNBits = 0;
return;
}
}
// }
if (!fFoundQuartet) { // Wait until we have found quad phase before starting
return; // to calibrate
}
// LoadHelicity is called first event of each cycle.
// But only do it's thing if it is first cycle of quartet.
// Need to instead have it do it's thing on the first event of a quartet.
// But only for seed building. Current logic is OK once seed is found.
if (fNBits < fMAXBIT) {
if (fNQuartet != fNLastQuartet) {
// Sanity check fNQuartet == fNLastQuartet+1
if (fNBits == 0) {
cout << "Start calibrating at cycle " << cyclecount << endl;
fRingSeed_reported = 0;
}
// Do phase stuff right here
if ((fReportedHelicity == kPlus && (quartetphase == 0 || quartetphase == 3)) ||
(fReportedHelicity == kMinus && (quartetphase == 1 || quartetphase == 2))) {
fRingSeed_reported = ((fRingSeed_reported << 1) | 1) & 0x3FFFFFFF;
// cout << " " << fNQuartet << " 1" << endl;
} else {
fRingSeed_reported = (fRingSeed_reported << 1) & 0x3FFFFFFF;
// cout << " " << fNQuartet << " 0" << endl;
}
fNBits++;
if (fReportedHelicity == kUnknown) {
fNBits = 0;
fRingSeed_reported = 0;
} else if (fNBits == fMAXBIT) {
cout << "Seed Found " << bitset<32>(fRingSeed_reported) << " at cycle " << cyclecount
<< " with first cycle " << fFirstCycle << endl;
if (fglHelicityScaler) {
cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
}
Int_t backseed = GetSeed30(fRingSeed_reported);
cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec << endl;
// Create the "actual seed"
fRingSeed_actual = fRingSeed_reported;
for (Int_t i = 0; i < fHelDelay / 4; i++) {
fRingSeed_actual = RanBit30(fRingSeed_actual);
}
}
fActualHelicity = kUnknown;
} // Need to change this to build seed even when not at start of quartet
} else {
if (quartetphase == 0) {
// If quartetphase !=, the seeds will alread have been advanced
// except that we won't have made the initial
// cout << "Advancing seed B " << fNBits << endl;
fRingSeed_reported = RanBit30(fRingSeed_reported);
fRingSeed_actual = RanBit30(fRingSeed_actual);
fActualHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
fPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
// if(fTITime/250000000.0 > 380.0) cout << fTITime/250000000.0 << " " << fNCycle << " "
// << hex <<
// 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;
fNBits = 0; // Need to reaquire seed
fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
}
fQuartetStartHelicity = fActualHelicity;
fQuartetStartPredictedHelicity = fPredictedHelicity;
}
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)