Skip to content
Snippets Groups Projects
Commit 47c908bc authored by Stephen A. Wood's avatar Stephen A. Wood Committed by Stephen Wood
Browse files

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.
parent 7039ef27
No related branches found
No related tags found
No related merge requests found
//*-- 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)
#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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment