diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx index ae335b1cd918af23fceee1b5283ac78cfd32134d..80b2610d3f87c627bf45259f8540aba88ad18f24 100644 --- a/src/THcHelicity.cxx +++ b/src/THcHelicity.cxx @@ -1,10 +1,8 @@ -//*-- Author : Julie Roche, November 2010 -// this is a modified version of THaGOHelicity.C //////////////////////////////////////////////////////////////////////// // // THcHelicity // -// Helicity of the beam from QWEAK electronics in delayed mode +// Helicity of the beam in delayed mode using quartets // +1 = plus, -1 = minus, 0 = unknown // // Also supports in-time mode with delay = 0 @@ -17,8 +15,10 @@ #include "THaApparatus.h" #include "THaEvData.h" #include "THcGlobals.h" +#include "THcHelicityScaler.h" #include "THcParmList.h" #include "TMath.h" +#include <bitset> #include <iostream> using namespace std; @@ -72,6 +72,21 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) { 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; } @@ -161,6 +176,7 @@ Int_t THcHelicity::DefineVariables(EMode mode) { fIsSetup = (mode == kDefine); // Define standard variables from base class +<<<<<<< HEAD THaHelicityDet::DefineVariables(mode); const RVarDef var[] = {{"nqrt", "position of cycle in quartet", "fnQrt"}, @@ -172,6 +188,21 @@ Int_t THcHelicity::DefineVariables(EMode mode) { // cout << "Calling THcHelicity DefineVarsFromList" << endl; _logger->info("Calling THcHelicity DefineVarsFromList"); return DefineVarsFromList(var, mode); +======= + 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}}; + cout << "Calling THcHelicity DefineVarsFromList" << endl; + return DefineVarsFromList(var, mode); +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 } //_____________________________________________________________________________ @@ -231,12 +262,18 @@ void THcHelicity::Clear(Option_t* opt) { 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) { @@ -244,10 +281,33 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) { return err; } +<<<<<<< HEAD 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; +======= + 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; +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 return 0; } @@ -275,8 +335,15 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) { // << fLastMPSTime << " " << fNBits << endl; Int_t missed = 0; // Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0; +<<<<<<< HEAD if (fIsMPS) { fActualHelicity = kUnknown; +======= + if (fIsMPS) { + fPeriodCheck = fmod(fTITime / fTIPeriod, 1.0); + fCycle = (fTITime / fTIPeriod); + fActualHelicity = kUnknown; +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 fPredictedHelicity = kUnknown; if (fFoundMPS) { missed = TMath::Nint(fTITime / fTIPeriod - fLastMPSTime / fTIPeriod); @@ -296,6 +363,7 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) { fFoundMPS = kTRUE; fLastMPSTime = fTITime; } +<<<<<<< HEAD } else if (fFoundMPS) { // if (fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods missed = TMath::Nint(floor((fTITime - fLastMPSTime) / fTIPeriod)); @@ -405,6 +473,172 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) { // cout << fTITime/250000000.0 << " " << fNCycle << " " << fReportedHelicity << endl; // cout << fNCycle << ": " << fReportedHelicity << " " // << fPredictedHelicity << " " << fActualHelicity << endl; +======= + } 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; +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 } // Ignore until a MPS Is found } else { // No MPS found yet @@ -440,7 +674,6 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) { fLastActualHelicity = fActualHelicity; return 0; } - //_____________________________________________________________________________ Int_t THcHelicity::End(THaRunBase*) { // End of run processing. Write histograms. @@ -467,6 +700,7 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m int quartetphase = (cyclecount - fFirstCycle) % 4; fnQrt = quartetphase; +<<<<<<< HEAD if (missedcycles > 1) { // If we missed windows if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over fNBits = 0; @@ -530,8 +764,87 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m fNBits = 0; // Need to reaquire seed fActualHelicity = kUnknown; fPredictedHelicity = kUnknown; +======= + // 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; +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 } + fQuartetStartHelicity = fActualHelicity; + fQuartetStartPredictedHelicity = fPredictedHelicity; } +<<<<<<< HEAD fQuartetStartHelicity = fActualHelicity; fQuartetStartPredictedHelicity = fPredictedHelicity; } else { // Not the beginning of a quad @@ -542,6 +855,12 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m ? fQuartetStartPredictedHelicity : -fQuartetStartPredictedHelicity; } +======= + fActualHelicity = + (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity : -fQuartetStartHelicity; + fPredictedHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartPredictedHelicity + : -fQuartetStartPredictedHelicity; +>>>>>>> 77b916167fbc44bfd82e8c4fe2b063c386e89d99 } return; } diff --git a/src/THcHelicity.h b/src/THcHelicity.h index 70c651e529ef26e0b829080c8226018f2e47b18d..65c1bb00a6d14b3e8b1f88c24f4c4602eb268834 100644 --- a/src/THcHelicity.h +++ b/src/THcHelicity.h @@ -14,6 +14,7 @@ #include "hcana/Logger.h" class TH1F; +class THcHelicityScaler; class THcHelicity : public hcana::ConfigLogging<THaHelicityDet>, public THcHelicityReader { @@ -32,7 +33,7 @@ public: 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; } + virtual void SetHelicityScaler(THcHelicityScaler *f); void PrintEvent(Int_t evtnum); @@ -48,10 +49,14 @@ protected: // Fixed Parameters Int_t fRingSeed_reported_initial; Int_t fFirstCycle; + Bool_t fFixFirstCycle; Double_t fFreq; Double_t fTIPeriod; // Reversal period in TI time units + Double_t fPeriodCheck; + Double_t fCycle; + Bool_t fFirstEvProcessed; Int_t fLastReportedHelicity; Long64_t fFirstEvTime; @@ -67,9 +72,13 @@ protected: Bool_t fFoundQuartet; // True if quartet phase probably found. Bool_t fIsNewCycle; Int_t fNCycle; // Count of # of helicity cycles + Int_t fNQuartet; // Quartet count + Int_t fNLastQuartet; Int_t fQuartet[4]; // For finding and checking quartet pattern Int_t fNBits; Int_t fnQrt; // Position in quartet + Bool_t fHaveQRT; // QRT signal exists + Int_t fNQRTProblems; Int_t fRingSeed_reported; Int_t fRingSeed_actual; @@ -86,11 +95,7 @@ protected: Int_t fQrt; - Int_t fTSettle; - Bool_t fValidHel; - Int_t fHelicityLastTIR; - Int_t fPatternLastTIR; void SetErrorCode(Int_t error); Double_t fErrorCode; @@ -105,6 +110,20 @@ protected: virtual Int_t DefineVariables( EMode mode = kDefine ); virtual Int_t ReadDatabase( const TDatime& date ); + THcHelicityScaler *fglHelicityScaler; + Int_t* fHelicityHistory; + Int_t fLastHelpCycle; + Int_t fScaleQuartet; + Int_t fQuadPattern[8]; + Int_t fHelperHistory; + Int_t fHelperQuartetHistory; + Int_t fScalerSeed; + Int_t lastispos; + Int_t evnum; + Int_t fThisScaleHel; + Int_t fLastScaleHel; + Int_t fLastLastScaleHel; + ClassDef(THcHelicity,0) // Beam helicity from QWEAK electronics in delayed mode }; diff --git a/src/THcHelicityReader.cxx b/src/THcHelicityReader.cxx index ffc8eadb561688023972c295682970519bac9f8f..3319da87bf263a931801463e536b3d7fdbe432e0 100644 --- a/src/THcHelicityReader.cxx +++ b/src/THcHelicityReader.cxx @@ -7,30 +7,26 @@ //////////////////////////////////////////////////////////////////////// #include "THcHelicityReader.h" +#include "TError.h" +#include "TH1F.h" +#include "THaAnalysisObject.h" // For LoadDB #include "THaEvData.h" #include "THcGlobals.h" #include "THcParmList.h" #include "TMath.h" -#include "TError.h" #include "VarDef.h" -#include "THaAnalysisObject.h" // For LoadDB #include <iostream> #include <vector> -#include "TH1F.h" using namespace std; //____________________________________________________________________ THcHelicityReader::THcHelicityReader() - : fTITime(0), fTITime_last(0), fTITime_rollovers(0), - fHaveROCs(kFALSE) -{ + : fTITime(0), fTITime_last(0), fTITime_rollovers(0), fHaveROCs(kFALSE) { // Default constructor - } //____________________________________________________________________ -THcHelicityReader::~THcHelicityReader() -{ +THcHelicityReader::~THcHelicityReader() { // Destructor // Histograms will be deleted by ROOT @@ -40,34 +36,25 @@ THcHelicityReader::~THcHelicityReader() } //____________________________________________________________________ -void THcHelicityReader::Clear( Option_t* ) -{ - - fIsMPS = fIsQrt = fIsHelp = fIsHelm = kFALSE; - -} +void THcHelicityReader::Clear(Option_t*) { fIsMPS = fIsQrt = fIsHelp = fIsHelm = kFALSE; } //_____________________________________________________________________________ -Int_t THcHelicityReader::ReadDatabase( const char* /*dbfilename*/, - const char* /*prefix*/, - const TDatime& /*date*/, - int /*debug_flag*/ ) -{ +Int_t THcHelicityReader::ReadDatabase(const char* /*dbfilename*/, const char* /*prefix*/, + const TDatime& /*date*/, int /*debug_flag*/) { // Eventually get these from the parameter file - SetROCinfo(kHel,2,14,9); - SetROCinfo(kHelm,2,14,8); - SetROCinfo(kMPS,2,14,10); - SetROCinfo(kQrt,2,14,7); - SetROCinfo(kTime,2,21,2); - + // SHMS settings see https://logbooks.jlab.org/entry/3614445 + cout << "THcHelicityReader: Helicity information from ROC 2 (SHMS)" << endl; + SetROCinfo(kHel, 2, 14, 9); + SetROCinfo(kHelm, 2, 14, 8); + SetROCinfo(kMPS, 2, 14, 10); + SetROCinfo(kQrt, 2, 14, 7); // Starting about run 5818 + SetROCinfo(kTime, 2, 21, 2); + fADCThreshold = 8000; - - DBRequest list[] = { - {"helicity_adcthreshold",&fADCThreshold, kInt, 0, 1}, - {0} - }; + + DBRequest list[] = {{"helicity_adcthreshold", &fADCThreshold, kInt, 0, 1}, {0}}; gHcParms->LoadParmValues(list, ""); @@ -75,28 +62,25 @@ Int_t THcHelicityReader::ReadDatabase( const char* /*dbfilename*/, } //_____________________________________________________________________________ -void THcHelicityReader::Begin() -{ +void THcHelicityReader::Begin() { // static const char* const here = "THcHelicityReader::Begin"; // cout<<here<<endl; - fTITime_last = 0; - fTITime = 0; + fTITime_last = 0; + fTITime = 0; fTITime_rollovers = 0; return; } //____________________________________________________________________ -void THcHelicityReader::End() -{ +void THcHelicityReader::End() { // static const char* const here = "THcHelicityReader::End"; // cout<<here<<endl; return; } //____________________________________________________________________ -Int_t THcHelicityReader::ReadData( const THaEvData& evdata ) -{ +Int_t THcHelicityReader::ReadData(const THaEvData& evdata) { // Obtain the present data from the event for QWEAK helicity mode. static const char* here = "THcHelicityReader::ReadData"; @@ -107,27 +91,37 @@ Int_t THcHelicityReader::ReadData( const THaEvData& evdata ) // std::cout<<" which="<<jk // <<" roc="<<fROCinfo[jk].roc // <<" header="<<fROCinfo[jk].header - // <<" index="<<fROCinfo[jk].index + // <<" index="<<fROCinfo[jk].index // <<endl; // } - + // std::cout<<" fHaveROCs="<<fHaveROCs<<endl; - if( !fHaveROCs ) { - ::Error( here, "ROC data (detector map) not properly set up." ); + if (!fHaveROCs) { + ::Error(here, "ROC data (detector map) not properly set up."); return -1; } + // Check if ROC info is correct + if (!evdata.GetModule(fROCinfo[kTime].roc, fROCinfo[kTime].slot)) { + cout << "THcHelicityReader: ROC 2 not found" << endl; + cout << "Changing to ROC 1 (HMS)" << endl; + SetROCinfo(kHel, 1, 18, 9); + SetROCinfo(kHelm, 1, 18, 8); + SetROCinfo(kMPS, 1, 18, 10); + SetROCinfo(kQrt, 1, 18, 7); + SetROCinfo(kTime, 1, 21, 2); + } + // Get the TI Data // Int_t fTIType = evData.GetData(fTICrate, fTISlot, 0, 0); // Int_t fTIEvNum = evData.GetData(fTICrate, fTISlot, 1, 0); - UInt_t titime = (UInt_t) evdata.GetData(fROCinfo[kTime].roc, - fROCinfo[kTime].slot, - fROCinfo[kTime].index, 0); - //cout << fTITime_last << " " << titime << endl; - if(titime < fTITime_last) { + UInt_t titime = + (UInt_t)evdata.GetData(fROCinfo[kTime].roc, fROCinfo[kTime].slot, fROCinfo[kTime].index, 0); + // cout << fTITime_last << " " << titime << endl; + if (titime < fTITime_last) { fTITime_rollovers++; } - fTITime = titime + fTITime_rollovers*4294967296; + fTITime = titime + fTITime_rollovers * 4294967296; fTITime_last = titime; const_cast<THaEvData&>(evdata).SetEvTime(fTITime); @@ -135,36 +129,26 @@ Int_t THcHelicityReader::ReadData( const THaEvData& evdata ) // Get the helicity control signals. These are from the pedestals // acquired by FADC channels. - Int_t helpraw = evdata.GetData(Decoder::kPulsePedestal, - fROCinfo[kHel].roc, - fROCinfo[kHel].slot, - fROCinfo[kHel].index, 0); - Int_t helmraw = evdata.GetData(Decoder::kPulsePedestal, - fROCinfo[kHelm].roc, - fROCinfo[kHelm].slot, - fROCinfo[kHelm].index, 0); - Int_t mpsraw = evdata.GetData(Decoder::kPulsePedestal, - fROCinfo[kMPS].roc, - fROCinfo[kMPS].slot, - fROCinfo[kMPS].index, 0); - Int_t qrtraw = evdata.GetData(Decoder::kPulsePedestal, - fROCinfo[kQrt].roc, - fROCinfo[kQrt].slot, - fROCinfo[kQrt].index, 0); - - fIsQrt = qrtraw > fADCThreshold; - fIsMPS = mpsraw > fADCThreshold; + Int_t helpraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kHel].roc, fROCinfo[kHel].slot, + fROCinfo[kHel].index, 0); + Int_t helmraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kHelm].roc, fROCinfo[kHelm].slot, + fROCinfo[kHelm].index, 0); + Int_t mpsraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kMPS].roc, fROCinfo[kMPS].slot, + fROCinfo[kMPS].index, 0); + Int_t qrtraw = evdata.GetData(Decoder::kPulsePedestal, fROCinfo[kQrt].roc, fROCinfo[kQrt].slot, + fROCinfo[kQrt].index, 0); + + fIsQrt = qrtraw > fADCThreshold; + fIsMPS = mpsraw > fADCThreshold; fIsHelp = helpraw > fADCThreshold; fIsHelm = helmraw > fADCThreshold; return 0; } -//TODO: this should not be needed once LoadDB can fill fROCinfo directly +// TODO: this should not be needed once LoadDB can fill fROCinfo directly //____________________________________________________________________ -Int_t THcHelicityReader::SetROCinfo( EROC which, Int_t roc, - Int_t slot, Int_t index ) -{ +Int_t THcHelicityReader::SetROCinfo(EROC which, Int_t roc, Int_t slot, Int_t index) { // Define source and offset of data. Normally called by ReadDatabase // of the detector that is a THcHelicityReader. @@ -173,21 +157,21 @@ Int_t THcHelicityReader::SetROCinfo( EROC which, Int_t roc, // You must define at least the kHel and kTime ROCs. // Returns <0 if parameter error, 0 if success - if( which<kHel || which>=kCount ) + if (which < kHel || which >= kCount) return -1; - if( roc <= 0 || roc > 255 ) + if (roc <= 0 || roc > 255) return -2; - fROCinfo[which].roc = roc; - fROCinfo[which].slot = slot; - fROCinfo[which].index = index; + fROCinfo[which].roc = roc; + fROCinfo[which].slot = slot; + fROCinfo[which].index = index; + + _det_logger->info("SetROCInfo: {} (roc: {}, slot: {}, index: {})", which, + fROCinfo[which].roc, fROCinfo[which].slot, fROCinfo[which].index]); + fHaveROCs = (fROCinfo[kHel].roc > 0 && fROCinfo[kTime].roc > 0 && fROCinfo[kMPS].roc); - //cout << "SetROCInfo: " << which << " " << fROCinfo[kHel].roc << " " << fROCinfo[kTime].roc << endl; - fHaveROCs = ( fROCinfo[kHel].roc > 0 && fROCinfo[kTime].roc > 0 ); - return 0; } //____________________________________________________________________ ClassImp(THcHelicityReader) - diff --git a/src/THcHelicityScaler.cxx b/src/THcHelicityScaler.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4c482faa2437c5820df69d5f69a4a76659216bb0 --- /dev/null +++ b/src/THcHelicityScaler.cxx @@ -0,0 +1,306 @@ +/** \class THcHelicityScaler + \ingroup Base + +\brief Event handler for Hall C helicity scalers + +Scalers not yet implemented. For now just picks the helicity control +bits out of the scaler words + +~~~ + gHaEvtHandlers->Add (new THcHelicityScaler("H","HC helicity scalers")); +~~~ +To enable debugging you may try this in the setup script +~~~ + THcHelcityScaler *hhelscaler = new THcHelicityScaler("H","HC helicity scalers"); + hscaler->SetDebugFile("HHelScaler.txt"); + gHaEvtHandlers->Add (hhelscaler); +~~~ +\author +*/ + +#include "THaEvtTypeHandler.h" +#include "THcHelicityScaler.h" +#include "THaCodaData.h" +#include "THaEvData.h" +#include "THcParmList.h" +#include "THcGlobals.h" +#include "THaGlobals.h" +#include "THcHelicity.h" +#include "TNamed.h" +#include "TMath.h" +#include "TString.h" +#include "TROOT.h" +#include <cstring> +#include <cstdio> +#include <cstdlib> +#include <iostream> +#include <sstream> +#include <map> +#include <iterator> +#include "THaVarList.h" +#include "VarDef.h" +#include "Helper.h" + +using namespace std; +using namespace Decoder; + +static const UInt_t ICOUNT = 1; +static const UInt_t IRATE = 2; +static const UInt_t ICURRENT = 3; +static const UInt_t ICHARGE = 4; +static const UInt_t ITIME = 5; +static const UInt_t ICUT = 6; +static const UInt_t MAXCHAN = 32; +static const UInt_t defaultDT = 4; + +THcHelicityScaler::THcHelicityScaler(const char *name, const char* description) + : THaEvtTypeHandler(name,description), + fBCM_Gain(0), fBCM_Offset(0), fBCM_delta_charge(0), + fBankID(9801), + evcount(0), evcountR(0.0), ifound(0), + fUseFirstEvent(kTRUE), + fOnlySyncEvents(kFALSE), fOnlyBanks(kFALSE), fDelayedType(-1), + fClockChan(-1), fLastClock(0), fClockOverflows(0) +{ + fRocSet.clear(); +} + +THcHelicityScaler::~THcHelicityScaler() +{ + delete [] fBCM_Gain; + delete [] fBCM_Offset; + delete [] fBCM_delta_charge; + + for( vector<UInt_t*>::iterator it = fDelayedEvents.begin(); + it != fDelayedEvents.end(); ++it ) + delete [] *it; + fDelayedEvents.clear(); +} + +Int_t THcHelicityScaler::End( THaRunBase* ) +{ + // Process any delayed events in order received + + cout << "THcHelicityScaler::End Analyzing " << fDelayedEvents.size() << " delayed scaler events" << endl; + for(std::vector<UInt_t*>::iterator it = fDelayedEvents.begin(); + it != fDelayedEvents.end(); ++it) { + UInt_t* rdata = *it; + AnalyzeBuffer(rdata,kFALSE); + } + + for( vector<UInt_t*>::iterator it = fDelayedEvents.begin(); + it != fDelayedEvents.end(); ++it ) + delete [] *it; + fDelayedEvents.clear(); + + return 0; +} + + +Int_t THcHelicityScaler::ReadDatabase(const TDatime& date ) +{ + char prefix[2]; + prefix[0]='g'; + prefix[1]='\0'; + + return kOK; +} +void THcHelicityScaler::SetDelayedType(int evtype) { + /** + * \brief Delay analysis of this event type to end. + * + * Final scaler events generated in readout list end routines may not + * come in order in the data stream. If the event type of a end routine + * scaler event is set, then the event contents will be saved and analyzed + * at the end of the analysis so that time ordering of scaler events is preserved. + */ + fDelayedType = evtype; +} + +Int_t THcHelicityScaler::Analyze(THaEvData *evdata) +{ + + if ( !IsMyEvent(evdata->GetEvType()) ) return -1; + + if (fDebugFile) { + *fDebugFile << endl << "---------------------------------- "<<endl<<endl; + *fDebugFile << "\nEnter THcHelicityScaler for fName = "<<fName<<endl; + EvDump(evdata); + } + + UInt_t *rdata = (UInt_t*) evdata->GetRawDataBuffer(); + + if(evdata->GetEvType() == fDelayedType) { // Save this event for processing later + Int_t evlen = evdata->GetEvLength(); + UInt_t *datacopy = new UInt_t[evlen]; + fDelayedEvents.push_back(datacopy); + memcpy(datacopy,rdata,evlen*sizeof(UInt_t)); + return 1; + } else { // A normal event + if (fDebugFile) *fDebugFile<<"\n\nTHcHelicityScaler :: Debugging event type "<<dec<<evdata->GetEvType()<< " event num = " << evdata->GetEvNum() << endl<<endl; + evNumber=evdata->GetEvNum(); + Int_t ret; + if((ret=AnalyzeBuffer(rdata,fOnlySyncEvents))) { + // + } + return ret; + } + +} +Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync) +{ + fNevents = 0; + + // Parse the data, load local data arrays. + UInt_t *p = (UInt_t*) rdata; + + UInt_t *plast = p+*p; // Index to last word in the bank + Int_t roc = -1; + Int_t evlen = *p+1; + + ifound=0; + + while(p<plast) { + Int_t banklen = *p; + p++; // point to header + + if (fDebugFile) { + *fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p-1) << endl; + } + if((*p & 0xff00) == 0x1000) { // Bank Containing banks + if(evlen-*(p-1) > 1) { // Don't use overall event header + roc = (*p>>16) & 0xf; + if(fDebugFile) *fDebugFile << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p << dec << endl; +// cout << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p << dec << endl; + if(fRocSet.find(roc)==fRocSet.end()) { // Not a ROC with helicity scaler + p+=*(p-1)-1; // Skip to next ROC + } + } + p++; // Now pointing to a bank in the bank + } else if (((*p & 0xff00) == 0x100) && (*p != 0xC0000100)) { + // Bank containing integers. Look for scalers + // This is either ROC bank containing integers or + // a bank within a ROC containing data from modules of a single type + // Look for banks with the helicity scaler bank ID (9801) + // Assume that very first word is a scaler header + // At any point in the bank where the word is not a matching + // header, we stop. + UInt_t tag = (*p>>16) & 0xffff; // Bank ID (ROC #) + // UInt_t num = (*p) & 0xff; + UInt_t *pnext = p+banklen; // Next bank + p++; // First data word + // If the bank is not a helicity scaler bank + // or it is not one of the ROC containing helcity scaler data + // skip to the next bank + //cout << "BankID=" << tag << endl; + + if (tag != fBankID) { + p = pnext; // Fall through to end of the above else if + // cout << " Skipping to next bank" << endl; + + } else { + // This is a helicity scaler bank + if (roc == 5) { + Int_t nevents = (banklen-2)/32; + //cout << "# of helicity events in bank:" << " " << nevents << endl; + if (nevents > 100) { + cout << "Error! Beam off for too long" << endl; + } + fNevents = nevents; + fNcycles += nevents; + + // Save helcitiy and quad info for THcHelicity + for (Int_t iev = 0; iev < nevents; iev++) { // find number of helicity events in each bank + Int_t nentries = 32*iev+2; + fHelicityHistory[iev] = (p[nentries-1]>>30) & 0x3; + // cout << "H: " << evNumber << endl; + } + } + } + + while(p < pnext) { + Int_t nskip = 0; + if(fDebugFile) { + *fDebugFile << "Scaler Header: " << hex << *p << dec; + } + if(nskip == 0) { + if(fDebugFile) { + *fDebugFile << endl; + } + break; // Didn't find a matching header + } + p = p + nskip; + } + p = pnext; + } else { + p = p+*(p-1); // Skip to next bank + } + + } + + if (fDebugFile) { + *fDebugFile << "Finished with decoding. "<<endl; + *fDebugFile << " Found flag = "<<ifound<<endl; + } + + // HMS has headers which are different from SOS, but both are + // event type 0 and come here. If you found no headers, return. + + if (!ifound) return 0; + + return 1; + + +} + +//Int_t THcHelicityScaler::AnalyzeHelicityScaler(UInt_t *p) +//{ +//} + + +THaAnalysisObject::EStatus THcHelicityScaler::Init(const TDatime& date) +{ + + ReadDatabase(date); + + fStatus = kOK; + + for( vector<UInt_t*>::iterator it = fDelayedEvents.begin(); + it != fDelayedEvents.end(); ++it ) + delete [] *it; + fDelayedEvents.clear(); + + cout << "Howdy ! We are initializing THcHelicityScaler !! name = " + << fName << endl; + + if(eventtypes.size()==0) { + eventtypes.push_back(0); // Default Event Type + } + + fRocSet.insert(5); // List ROCs that have helicity scalers + fRocSet.insert(8); // Should make configurable + + fNcycles = 0; + fNevents = 0; + + return kOK; +} + +size_t THcHelicityScaler::FindNoCase(const string& sdata, const string& skey) +{ + // Find iterator of word "sdata" where "skey" starts. Case insensitive. + string sdatalc, skeylc; + sdatalc = ""; skeylc = ""; + for (string::const_iterator p = + sdata.begin(); p != sdata.end(); ++p) { + sdatalc += tolower(*p); + } + for (string::const_iterator p = + skey.begin(); p != skey.end(); ++p) { + skeylc += tolower(*p); + } + if (sdatalc.find(skeylc,0) == string::npos) return -1; + return sdatalc.find(skeylc,0); +}; + +ClassImp(THcHelicityScaler) diff --git a/src/THcHelicityScaler.h b/src/THcHelicityScaler.h new file mode 100644 index 0000000000000000000000000000000000000000..d467483469fe74d69251f5e5de30371e4ae38314 --- /dev/null +++ b/src/THcHelicityScaler.h @@ -0,0 +1,88 @@ +#ifndef THcHelicityScaler_ +#define THcHelicityScaler_ + +///////////////////////////////////////////////////////////////////// +// +// THcHelicityScaler +// class to handle Helicity scaler events +// +///////////////////////////////////////////////////////////////////// + +#include "THaEvtTypeHandler.h" +#include "Decoder.h" +#include <string> +#include <vector> +#include <algorithm> +#include <set> +#include "TTree.h" +#include "TString.h" +#include <cstring> + +class THcHelicity; + +class THcHelicityScaler : public THaEvtTypeHandler { + +public: + + THcHelicityScaler(const char*, const char*); + virtual ~THcHelicityScaler(); + + Int_t Analyze(THaEvData *evdata); + Int_t AnalyzeBuffer(UInt_t *rdata, Bool_t onlysync); + Int_t AnalyzeHelicityScaler(UInt_t *p); + + virtual EStatus Init( const TDatime& run_time); + virtual Int_t ReadDatabase(const TDatime& date ); + virtual Int_t End( THaRunBase* r=0 ); + + virtual void SetUseFirstEvent(Bool_t b = kFALSE) {fUseFirstEvent = b;} + virtual void SetDelayedType(int evtype); + virtual void SetROC(Int_t roc) {fROC=roc;} + virtual void SetBankID(Int_t bankid) {fBankID=bankid;} + virtual void SetHelicityDetector(THcHelicity *f) {fglHelicityDetector = f;} + virtual Int_t GetNevents() { return fNevents;} + virtual Int_t GetNcycles() { return fNcycles;} + virtual Int_t GetEvNum() { return evNumber;} + virtual Int_t* GetHelicityHistoryP() {return fHelicityHistory;} + +private: + + static size_t FindNoCase(const std::string& sdata, const std::string& skey); + + Int_t fNumBCMs; + Double_t *fBCM_Gain; + Double_t *fBCM_Offset; + Double_t *fBCM_delta_charge; + + Int_t fROC; + UInt_t fBankID; + // Helicity Scaler variables + Int_t fNevents; /* # of helicity scaler reads in last event */ + Int_t fNcycles; + Int_t fHelicityHistory[200]; + // + Bool_t fUseFirstEvent; + Bool_t fOnlySyncEvents; + Bool_t fOnlyBanks; + Int_t fDelayedType; + Int_t fClockChan; + UInt_t fLastClock; + Int_t fClockOverflows; + + std::vector<UInt_t*> fDelayedEvents; + std::set<UInt_t> fRocSet; + + THcHelicityScaler(const THcHelicityScaler& fh); + THcHelicityScaler& operator=(const THcHelicityScaler& fh); + + UInt_t evcount; + Double_t evcountR; + UInt_t evNumber; + Int_t ifound; + THcHelicity *fglHelicityDetector; + + ClassDef(THcHelicityScaler,0) // Scaler Event handler + +}; + +#endif