diff --git a/podd b/podd
index 492bc495f20d61033bc283a93260a81fd121d683..e58d6959ef72f188811ba485f53945dde13141bf 160000
--- a/podd
+++ b/podd
@@ -1 +1 @@
-Subproject commit 492bc495f20d61033bc283a93260a81fd121d683
+Subproject commit e58d6959ef72f188811ba485f53945dde13141bf
diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx
index 1dc6285783b36c9c30caec922f0e39592917c4b0..3b53d464d25481e83701a50dec76ca3a64cd4a4d 100644
--- a/src/THcAerogel.cxx
+++ b/src/THcAerogel.cxx
@@ -716,10 +716,14 @@ Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
Bool_t pulseTimeCut = adctdcdiffTime > fAdcPosTimeWindowMin[npmt] && adctdcdiffTime < fAdcPosTimeWindowMax[npmt];
// By default, the last hit within the timing cut will be considered "good"
- if (!errorFlag && pulseTimeCut) {
+ if (!errorFlag)
+ {
+ fGoodPosAdcMult.at(npmt) += 1;
+ }
+
+ if (!errorFlag && pulseTimeCut) {
fGoodPosAdcPed.at(npmt) = pulsePed;
- fGoodPosAdcMult.at(npmt) = frPosAdcPulseInt->GetEntries();
- // cout << " out = " << npmt << " " << frPosAdcPulseInt->GetEntries() << " " <<fGoodPosAdcMult.at(npmt);
+ // cout << " out = " << npmt << " " << frPosAdcPulseInt->GetEntries() << " " <<fGoodPosAdcMult.at(npmt);
fGoodPosAdcPulseInt.at(npmt) = pulseInt;
fGoodPosAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodPosAdcPulseAmp.at(npmt) = pulseAmp;
@@ -748,15 +752,18 @@ Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
Bool_t errorFlag = ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
//// Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcNegTimeWindowMin[npmt] && adctdcdiffTime < fAdcNegTimeWindowMax[npmt];
-
+ if (!errorFlag)
+ {
+ fGoodNegAdcMult.at(npmt) += 1;
+ }
+
// By default, the last hit within the timing cut will be considered "good"
if (!errorFlag && pulseTimeCut) {
fGoodNegAdcPed.at(npmt) = pulsePed;
- fGoodNegAdcMult.at(npmt) = frNegAdcPulseInt->GetEntries();
- fGoodNegAdcPulseInt.at(npmt) = pulseInt;
fGoodNegAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodNegAdcPulseAmp.at(npmt) = pulseAmp;
- fGoodNegAdcPulseTime.at(npmt) = pulseTime;
+ fGoodNegAdcPulseInt.at(npmt) = pulseInt;
+ fGoodNegAdcPulseTime.at(npmt) = pulseTime;
fGoodNegAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fNegNpe.at(npmt) = fNegGain[npmt]*fGoodNegAdcPulseInt.at(npmt);
diff --git a/src/THcDC.h b/src/THcDC.h
index c4fe61a863ae8d79db9008d9abc99236672e7bd5..80918fd04db076d448d23d98310e115080a13e55 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -189,7 +189,7 @@ protected:
// double tan_angle, sin_angle, cos_angle;
// Intermediate structure for building
- static const char MAXTRACKS = 10;
+ static const UInt_t MAXTRACKS = 10;
std::vector<THcDriftChamberPlane*> fPlanes; // List of plane objects
std::vector<THcDriftChamber*> fChambers; // List of chamber objects
diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a5e14f5bf60b03476c5c0e89c9076d67713b773d
--- /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 0000000000000000000000000000000000000000..ff31ce188e858f20247b909369dbc2cf37dec568
--- /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
+
diff --git a/src/THcHelicityReader.cxx b/src/THcHelicityReader.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..66814306b5042725dfc65620b53b9c5e6b92657a
--- /dev/null
+++ b/src/THcHelicityReader.cxx
@@ -0,0 +1,193 @@
+//*-- Author : Ole Hansen August 2006
+// Extracted from Bob Michaels' THaHelicity CVS 1.19
+////////////////////////////////////////////////////////////////////////
+//
+// THcHelicityReader
+//
+////////////////////////////////////////////////////////////////////////
+
+#include "THcHelicityReader.h"
+#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)
+{
+ // Default constructor
+
+}
+//____________________________________________________________________
+THcHelicityReader::~THcHelicityReader()
+{
+ // Destructor
+
+ // Histograms will be deleted by ROOT
+ // for( Int_t i = 0; i < NHISTR; ++i ) {
+ // delete fHistoR[i];
+ // }
+}
+
+//____________________________________________________________________
+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*/ )
+{
+
+ // 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);
+
+ fADCThreshold = 8000;
+
+ DBRequest list[] = {
+ {"helicity_adcthreshold",&fADCThreshold, kInt, 0, 1},
+ {0}
+ };
+
+ gHcParms->LoadParmValues(list, "");
+
+ return THaAnalysisObject::kOK;
+}
+
+//_____________________________________________________________________________
+void THcHelicityReader::Begin()
+{
+ // static const char* const here = "THcHelicityReader::Begin";
+ // cout<<here<<endl;
+
+ fTITime_last = 0;
+ fTITime = 0;
+ fTITime_rollovers = 0;
+
+ return;
+}
+//____________________________________________________________________
+void THcHelicityReader::End()
+{
+ // static const char* const here = "THcHelicityReader::End";
+ // cout<<here<<endl;
+
+ return;
+}
+//____________________________________________________________________
+Int_t THcHelicityReader::ReadData( const THaEvData& evdata )
+{
+ // Obtain the present data from the event for QWEAK helicity mode.
+
+ static const char* here = "THcHelicityReader::ReadData";
+
+ // std::cout<<" kHel, kTime, kRing="<< kHel<<" "<<kTime<<" "<<kRing<<endl;
+ // for (int jk=0; jk<3; jk++)
+ // {
+ // std::cout<<" which="<<jk
+ // <<" roc="<<fROCinfo[jk].roc
+ // <<" header="<<fROCinfo[jk].header
+ // <<" index="<<fROCinfo[jk].index
+ // <<endl;
+ // }
+
+ // std::cout<<" fHaveROCs="<<fHaveROCs<<endl;
+ if( !fHaveROCs ) {
+ ::Error( here, "ROC data (detector map) not properly set up." );
+ return -1;
+ }
+
+ // 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) {
+ fTITime_rollovers++;
+ }
+ fTITime = titime + fTITime_rollovers*4294967296;
+ fTITime_last = titime;
+
+ const_cast<THaEvData&>(evdata).SetEvTime(fTITime);
+
+ // 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;
+ fIsHelp = helpraw > fADCThreshold;
+ fIsHelm = helmraw > fADCThreshold;
+
+ return 0;
+}
+
+//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 )
+{
+
+ // Define source and offset of data. Normally called by ReadDatabase
+ // of the detector that is a THcHelicityReader.
+ //
+ // "which" is one of { kHel, kKelm, kQrt, kTime }.
+ // You must define at least the kHel and kTime ROCs.
+ // Returns <0 if parameter error, 0 if success
+
+ if( which<kHel || which>=kCount )
+ return -1;
+ if( roc <= 0 || roc > 255 )
+ return -2;
+
+ fROCinfo[which].roc = roc;
+ fROCinfo[which].slot = slot;
+ fROCinfo[which].index = index;
+
+ 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/THcHelicityReader.h b/src/THcHelicityReader.h
new file mode 100644
index 0000000000000000000000000000000000000000..eef01503732231f90dfbca4132c3128e4dcf591f
--- /dev/null
+++ b/src/THcHelicityReader.h
@@ -0,0 +1,70 @@
+#ifndef Podd_THcHelicityReader_h_
+#define Podd_THcHelicityReader_h_
+
+//////////////////////////////////////////////////////////////////////////
+//
+// THcHelicityReader
+//
+// Routines for decoding QWEAK helicity hardware
+//
+//////////////////////////////////////////////////////////////////////////
+
+#include "Rtypes.h"
+
+class THaEvData;
+class TDatime;
+class TH1F;
+
+class THcHelicityReader {
+
+public:
+ THcHelicityReader();
+ virtual ~THcHelicityReader();
+
+ struct ROCinfo {
+ Int_t roc; // ROC to read out
+ Int_t slot; // Headers to search for (0 = ignore)
+ Int_t index; // Index into buffer
+ };
+
+protected:
+
+ // Used by ReadDatabase
+ enum EROC { kHel = 0, kHelm, kMPS, kQrt, kTime, kCount };
+ Int_t SetROCinfo( EROC which, Int_t roc, Int_t slot, Int_t index );
+
+ virtual void Clear( Option_t* opt="" );
+ virtual Int_t ReadData( const THaEvData& evdata );
+ Int_t ReadDatabase( const char* dbfilename, const char* prefix,
+ const TDatime& date, int debug_flag = 0 );
+ void Begin();
+ void End();
+
+ ULong64_t fTITime;
+ UInt_t fTITime_last;
+ UInt_t fTITime_rollovers;
+
+ // Reported Helicity status for the event
+ Bool_t fIsMPS;
+ Bool_t fIsQrt;
+ Bool_t fIsHelp;
+ Bool_t fIsHelm;
+
+ Int_t fADCThreshold; // Threshold for On/Off of helicity signals
+
+ ROCinfo fROCinfo[kCount];
+
+ Int_t fQWEAKDebug; // Debug level
+ Bool_t fHaveROCs; // Required ROCs are defined
+ Bool_t fNegGate; // Invert polarity of gate, TO DO implement this functionality
+ static const Int_t NHISTR = 12;
+ // TH1F* fHistoR[12]; // Histograms
+
+private:
+
+ ClassDef(THcHelicityReader,0) // Helper class for reading QWEAK helicity data
+
+};
+
+#endif
+
diff --git a/src/THcRun.cxx b/src/THcRun.cxx
index 449a1f8afc868b7895a87d4daf8205a841ef7202..06f5593ac596c61e6e33317f77f9b27b7e213d62 100644
--- a/src/THcRun.cxx
+++ b/src/THcRun.cxx
@@ -40,7 +40,7 @@ THcRun::THcRun( const vector<TString>& pathList, const char* filename,
}
//_____________________________________________________________________________
-THcRun& THcRun::operator=(const THaRun& rhs)
+THcRun& THcRun::operator=(const THaRunBase& rhs)
{
// Assignment operator. Not really sure what I (saw) am doing here.
diff --git a/src/THcRun.h b/src/THcRun.h
index fa9eadacca3a0e0c8cce1c5a039acb37f1cfe7b0..0fcd7771947c87c5eb7f777123fe7915b3d07220 100644
--- a/src/THcRun.h
+++ b/src/THcRun.h
@@ -17,7 +17,7 @@ class THcRun : public THaRun {
THcRun( const THcRun& run );
THcRun( const std::vector<TString>& pathList, const char* filename,
const char* description="" );
- virtual THcRun& operator=( const THaRun& rhs );
+ THcRun& operator=( const THaRunBase& rhs );
virtual ~THcRun();
virtual void Print( Option_t* opt="" ) const;
THcParmList* GetHCParms() const { return fHcParms; }
diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
index 6a30172f1c1cbac292d59a38e18164ef8bdf9aeb..2efd325caaeca442f5dd1e0019de1cada25971a7 100644
--- a/src/THcShowerArray.cxx
+++ b/src/THcShowerArray.cxx
@@ -356,6 +356,7 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
fNumGoodAdcHits = vector<Int_t> (fNelem, 0.0);
fGoodAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0);
fGoodAdcPed = vector<Double_t> (fNelem, 0.0);
+ fGoodAdcMult = vector<Double_t> (fNelem, 0.0);
fGoodAdcPulseInt = vector<Double_t> (fNelem, 0.0);
fGoodAdcPulseAmp = vector<Double_t> (fNelem, 0.0);
fGoodAdcPulseTime = vector<Double_t> (fNelem, 0.0);
@@ -442,6 +443,7 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
{"goodAdcPulseIntRaw", "Good Raw ADC Pulse Integrals", "fGoodAdcPulseIntRaw"}, //this is defined as pulseIntRaw, NOT ADC Amplitude in FillADC_DynamicPedestal() method
{"goodAdcPed", "Good ADC Pedestals", "fGoodAdcPed"},
+ {"goodAdcMult", "Good ADC Multiplicity", "fGoodAdcMult"},
{"goodAdcPulseInt", "Good ADC Pulse Integrals", "fGoodAdcPulseInt"}, //this is defined as pulseInt, which is the pedestal subtracted pulse integrals, and is defined if threshold is passed
{"goodAdcPulseAmp", "Good ADC Pulse Amplitudes", "fGoodAdcPulseAmp"},
{"goodAdcPulseTime", "Good ADC Pulse Times", "fGoodAdcPulseTime"}, //this is defined as pulseInt, which is the pedestal subtracted pulse integrals, and is defined if threshold is passed
@@ -496,6 +498,7 @@ void THcShowerArray::Clear( Option_t* )
for (UInt_t ielem = 0; ielem < fGoodAdcPed.size(); ielem++) {
fGoodAdcPulseIntRaw.at(ielem) = 0.0;
fGoodAdcPed.at(ielem) = 0.0;
+ fGoodAdcMult.at(ielem) = 0.0;
fGoodAdcPulseInt.at(ielem) = 0.0;
fGoodAdcPulseAmp.at(ielem) = 0.0;
fGoodAdcPulseTime.at(ielem) = kBig;
@@ -871,6 +874,10 @@ void THcShowerArray::FillADC_DynamicPedestal()
Bool_t errorflag = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData();
Bool_t pulseTimeCut = (adctdcdiffTime > fAdcTimeWindowMin[npad]) && (adctdcdiffTime < fAdcTimeWindowMax[npad]);
+ if (!errorflag)
+ {
+ fGoodAdcMult.at(npad) += 1;
+ }
if (!errorflag && pulseTimeCut) {
fTotNumAdcHits++;
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
index a44592b2555390eb06ff7cd2f91a16da2ec4c9f4..dc8991af37f3456bac0e5667594409c29979effa 100644
--- a/src/THcShowerArray.h
+++ b/src/THcShowerArray.h
@@ -113,6 +113,7 @@ protected:
vector<Double_t> fGoodAdcPulseIntRaw; // [fNelem] Good Raw ADC pulse Integrals of blocks
vector<Double_t> fGoodAdcPed; // [fNelem] Event by event (FADC) good pulse pedestals
+ vector<Double_t> fGoodAdcMult;
vector<Double_t> fGoodAdcPulseInt; // [fNelem] good pedestal subtracted pulse integrals
vector<Double_t> fGoodAdcPulseAmp;
vector<Double_t> fGoodAdcPulseTime;
diff --git a/src/THcTrigDet.cxx b/src/THcTrigDet.cxx
index b4c435cb409dcd884090e3e6cb249fc9432fb71e..bb0972580b98603035196c5c35605efded757d11 100644
--- a/src/THcTrigDet.cxx
+++ b/src/THcTrigDet.cxx
@@ -255,7 +255,7 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
UInt_t good_hit=999;
for (UInt_t thit=0; thit<rawAdcHit.GetNPulses(); ++thit) {
Int_t TestTime=rawAdcHit.GetPulseTimeRaw(thit);
- if (TestTime>fAdcTimeWindowMin[cnt]&&TestTime<fAdcTimeWindowMax[cnt]) {
+ if (TestTime>fAdcTimeWindowMin[cnt]&&TestTime<fAdcTimeWindowMax[cnt]&&good_hit==999) {
good_hit=thit;
}
}
@@ -276,7 +276,7 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
UInt_t good_hit=999;
for (UInt_t thit=0; thit<rawTdcHit.GetNHits(); ++thit) {
Int_t TestTime= rawTdcHit.GetTimeRaw(thit);
- if (TestTime>fTdcTimeWindowMin[cnt]&&TestTime<fTdcTimeWindowMax[cnt]) {
+ if (TestTime>fTdcTimeWindowMin[cnt]&&TestTime<fTdcTimeWindowMax[cnt]&&good_hit==999) {
good_hit=thit;
}
}