From 47c908bc1da9a0c984b7543aee6c44254ab5442e Mon Sep 17 00:00:00 2001 From: "Stephen A. Wood" <saw@jlab.org> Date: Wed, 7 Nov 2018 13:01:42 -0500 Subject: [PATCH] Determine the beam helicity for each event. Undoes delayed helicity reporting by learning seed of pseudo random sequence used by injector helicity electronics. This detector class should be added to the Trigger apparatus. Delivers the global variables "hel" and "helrep". "helrep" is the reported delayed helicity and "hel" is the actual helicity. + and - helicity are reported as +1 and -1. If the event is during a MPS settling period, then the helicity variables are zero. For the first few seconds of a run, the "hel" variable will be zero as the pseudo random sequence seed has not been determined. (If the trigger class is "T" and this detector is called "helicity", then these variables will be T.helicity.hel and T.helicity.helrep. By default, it is assumed that there is delayed reporting of 8 cycles and that Quartets are used. If delayed helicity reporting is not in use, then set the parameter "helicity_delay" to zero. --- src/THcHelicity.cxx | 536 ++++++++++++++++++++++++++++++++++++++++++++ src/THcHelicity.h | 109 +++++++++ 2 files changed, 645 insertions(+) create mode 100644 src/THcHelicity.cxx create mode 100644 src/THcHelicity.h diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx new file mode 100644 index 0000000..a5e14f5 --- /dev/null +++ b/src/THcHelicity.cxx @@ -0,0 +1,536 @@ +//*-- 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 "THaApparatus.h" +#include "THaEvData.h" +#include "THcGlobals.h" +#include "THcParmList.h" +#include "TH1F.h" +#include "TMath.h" +#include <iostream> + +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() + : 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; + + // Call initializer for base class. + // This also calls `ReadDatabase` and `DefineVariables`. + EStatus status = THaDetector::Init(date); + if (status) { + fStatus = status; + return fStatus; + } + + 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 ) +{ + + 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 + } + + return kOK; +} + +//_____________________________________________________________________________ + +void THcHelicity::MakePrefix() +{ + THaDetector::MakePrefix(); +} +//_____________________________________________________________________________ +Int_t THcHelicity::DefineVariables( EMode mode ) +{ + // Initialize global variables + + 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; + 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(fFirstEvProcessed) { // Normal processing + Int_t missed = 0; + // Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0; + if(fIsMPS) { + fFoundMPS = kTRUE; + Int_t missed = TMath::Nint(floor((fTITime-fLastMPSTime)/fTIPeriod)); + // cout << fTITime/250000000.0 << " " << fNCycle << " MPS " << fReportedHelicity <<endl; + if(missed <= 1) { + fLastMPSTime = fTITime; + fIsNewCycle = kTRUE; + fActualHelicity = kUnknown; + fPredictedHelicity = kUnknown; + } // 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 if (fFoundMPS) { // + if(fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods + Int_t 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); + } + // 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 (((fNCycle - fFirstCycle)%2)==1) { + // fQuartet[2] = fQuartet[3] = -fQuartet[0]; + // fQuartet[1] = fQuartet[0]; + fQuartet[0] = fReportedHelicity; + fQuartet[1] = fQuartet[2] = -fQuartet[0]; + } else { + // fQuartet[1] = fQuartet[2] = -fQuartet[0]; + // fQuartet[3] = fQuartet[0]; + fQuartet[0] = fQuartet[1] = -fReportedHelicity; + fQuartet[2] = -fQuartet[1]; + } + } else { + fQuartet[0] = fReportedHelicity; + fQuartet[1] = 0; + } + } + fLastMPSTime += missed*fTIPeriod; + fIsNewCycle = kTRUE; + fLastReportedHelicity = fReportedHelicity; + } + 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; + 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 + 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; + cout << "Searching for first of a quartet at cycle " << " " << fFirstCycle + << " - event " << evdata.GetEvNum() << endl; + cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " + << fQuartet[3] << endl; + 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 { + cout << "Initializing" << endl; + 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; + } + + 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) { + 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) { + cout << "Seed Found " << hex << fRingSeed_reported << dec << " at cycle " << cyclecount << " with first cycle " << fFirstCycle << endl; + Int_t backseed = GetSeed30(fRingSeed_reported); + 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) { + cout << "Helicity prediction failed " << fReportedHelicity << " " + << fPredictedHelicity << " " << fActualHelicity << 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) + diff --git a/src/THcHelicity.h b/src/THcHelicity.h new file mode 100644 index 0000000..ff31ce1 --- /dev/null +++ b/src/THcHelicity.h @@ -0,0 +1,109 @@ +#ifndef Podd_THcHelicity_h_ +#define Podd_THcHelicity_h_ + +//////////////////////////////////////////////////////////////////////// +// +// THcHelicity +// +// Helicity of the beam - from QWEAK electronics in delayed mode +// +//////////////////////////////////////////////////////////////////////// + +#include "THaHelicityDet.h" +#include "THcHelicityReader.h" + +class TH1F; + +class THcHelicity : public THaHelicityDet, public THcHelicityReader { + +public: + + THcHelicity( const char* name, const char* description, + THaApparatus* a = NULL ); + THcHelicity(); + virtual ~THcHelicity(); + + virtual EStatus Init(const TDatime& date); + virtual void MakePrefix(); + + virtual Int_t Begin( THaRunBase* r=0 ); + virtual void Clear( Option_t* opt = "" ); + virtual Int_t Decode( const THaEvData& evdata ); + virtual Int_t End( THaRunBase* r=0 ); + virtual void SetDebug( Int_t level ); + virtual Bool_t HelicityValid() const { return fValidHel; } + + void PrintEvent(Int_t evtnum); + +protected: + 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; + Double_t fFreq; + + ULong64_t fTIPeriod; // Reversal period in TI time units + + Bool_t fFirstEvProcessed; + Int_t fLastReportedHelicity; + ULong64_t fFirstEvTime; + ULong64_t fLastEvTime; + ULong64_t fLastMPSTime; + Int_t fReportedHelicity; + Int_t fMPS; + Int_t fPredictedHelicity; + Int_t fActualHelicity; + Int_t fQuartetStartHelicity; + Int_t fQuartetStartPredictedHelicity; + Bool_t fFoundMPS; + Bool_t fFoundQuartet; // True if quartet phase probably found. + Bool_t fIsNewCycle; + Int_t fNCycle; // Count of # of helicity cycles + Int_t fQuartet[4]; // For finding and checking quartet pattern + Int_t fNBits; + Int_t fnQrt; // Position in quartet + + Int_t fRingSeed_reported; + Int_t fRingSeed_actual; + + + // Offset between the ring reported value and the reported value + Int_t fHelDelay; + // delay of helicity (# windows) + 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 fQrt; + Int_t fTSettle; + Bool_t fValidHel; + + Int_t fHelicityLastTIR; + Int_t fPatternLastTIR; + void SetErrorCode(Int_t error); + Double_t fErrorCode; + + Int_t fEvtype; // Current CODA event type + + static const Int_t NHIST = 2; + TH1F* fHisto[NHIST]; + + virtual Int_t DefineVariables( EMode mode = kDefine ); + virtual Int_t ReadDatabase( const TDatime& date ); + + ClassDef(THcHelicity,0) // Beam helicity from QWEAK electronics in delayed mode + +}; + +#endif + -- GitLab