src/THcScintillatorPlane.h

@@ -63,6 +63,8 @@ class THcScintillatorPlane : public THaSubDetector {
   void SetFpTime(Double_t f) {fFptime=f;};
   void SetNGoodHits(Int_t ng) {fNGoodHits=ng;};
   void SetHitDistance(Double_t f) {fHitDistance=f;}; // Distance between track and hit paddle
+  void SetScinYPos(Double_t f) {fScinYPos=f;}; // Scint Average Y position at plane (cm)
+  void SetScinXPos(Double_t f) {fScinXPos=f;}; // Scint Average X position at plane (cm)
   void SetTrackXPosition(Double_t f) {fTrackXPosition=f;}; // Distance track X position at plane
   void SetTrackYPosition(Double_t f) {fTrackYPosition=f;}; // Distance track Y position at plane
   void SetNumberClusters(Int_t nclus) {fNumberClusters=nclus;}; // number of paddle 
@@ -181,7 +183,17 @@ class THcScintillatorPlane : public THaSubDetector {
   vector<Double_t>  fGoodDiffDistTrack;
 
   Int_t fDebugAdc;
+  Double_t fPosTdcRefTime;
+  Double_t fPosAdcRefTime;
+  Double_t fNegTdcRefTime;
+  Double_t fNegAdcRefTime;
+  Double_t fPosTdcRefDiffTime;
+  Double_t fPosAdcRefDiffTime;
+  Double_t fNegTdcRefDiffTime;
+  Double_t fNegAdcRefDiffTime;
   Double_t fHitDistance;
+  Double_t fScinXPos;
+  Double_t fScinYPos;
   Double_t fTrackXPosition;
   Double_t fTrackYPosition;
   Int_t fCosmicFlag; //
@@ -272,7 +284,10 @@ class THcScintillatorPlane : public THaSubDetector {
   Double_t *fNegPed;
   Double_t *fNegSig;
   Double_t *fNegThresh;
-
+  //
+   Int_t        fEventType;
+  Int_t        fEventNum;
+ 
   //
   Int_t fNScinGoodHits; // number of hits for which both ends of the paddle fired in time!
   Double_t fpTime; // the original code only has one fpTime per plane!

src/THcShower.cxx

@@ -36,6 +36,7 @@ THcShower::THcShower( const char* name, const char* description,
   hcana::ConfigLogging<THaNonTrackingDetector>(name,description,apparatus),
   fPosAdcTimeWindowMin(0), fNegAdcTimeWindowMin(0),
   fPosAdcTimeWindowMax(0), fNegAdcTimeWindowMax(0),
+  fPedPosDefault(0),fPedNegDefault(0),
   fShPosPedLimit(0), fShNegPedLimit(0), fPosGain(0), fNegGain(0),
   fClusterList(0), fLayerNames(0), fLayerZPos(0), BlockThick(0),
   fNBlocks(0), fXPos(0), fYPos(0), fZPos(0), fPlanes(0), fArray(0)
@@ -53,6 +54,7 @@ THcShower::THcShower( ) :
   hcana::ConfigLogging<THaNonTrackingDetector>(),
   fPosAdcTimeWindowMin(0), fNegAdcTimeWindowMin(0),
   fPosAdcTimeWindowMax(0), fNegAdcTimeWindowMax(0),
+  fPedPosDefault(0),fPedNegDefault(0),
   fShPosPedLimit(0), fShNegPedLimit(0), fPosGain(0), fNegGain(0),
   fClusterList(0), fLayerNames(0), fLayerZPos(0), BlockThick(0),
   fNBlocks(0), fXPos(0), fYPos(0), fZPos(0), fPlanes(0), fArray(0)
@@ -388,6 +390,8 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
   fNegAdcTimeWindowMin = new Double_t [fNTotBlocks];
   fPosAdcTimeWindowMax = new Double_t [fNTotBlocks];
   fNegAdcTimeWindowMax = new Double_t [fNTotBlocks];
+  fPedPosDefault = new Int_t [fNTotBlocks];
+  fPedNegDefault = new Int_t [fNTotBlocks];
 
   DBRequest list[]={
     {"cal_pos_cal_const", hcal_pos_cal_const, kDouble, fNTotBlocks},
@@ -400,6 +404,8 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
     {"cal_neg_AdcTimeWindowMin", fNegAdcTimeWindowMin, kDouble, static_cast<UInt_t>(fNTotBlocks),1},
     {"cal_pos_AdcTimeWindowMax", fPosAdcTimeWindowMax, kDouble, static_cast<UInt_t>(fNTotBlocks),1},
     {"cal_neg_AdcTimeWindowMax", fNegAdcTimeWindowMax, kDouble, static_cast<UInt_t>(fNTotBlocks),1},
+    {"cal_PedNegDefault", fPedNegDefault, kInt, static_cast<UInt_t>(fNTotBlocks),1},
+    {"cal_PedPosDefault", fPedNegDefault, kInt, static_cast<UInt_t>(fNTotBlocks),1},
     {"cal_min_peds", &fShMinPeds, kInt,0,1},
     {0}
   };
@@ -410,6 +416,8 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
     fNegAdcTimeWindowMin[ip] = -1000.;
     fPosAdcTimeWindowMax[ip] = 1000.;
     fNegAdcTimeWindowMax[ip] = 1000.;
+    fPedNegDefault[ip] = 0;
+    fPedPosDefault[ip] = 0;
    }
 
   gHcParms->LoadParmValues((DBRequest*)&list, prefix);
@@ -617,6 +625,8 @@ void THcShower::DeleteArrays()
   delete [] fNegAdcTimeWindowMin; fNegAdcTimeWindowMin = 0;
   delete [] fPosAdcTimeWindowMax; fPosAdcTimeWindowMax = 0;
   delete [] fNegAdcTimeWindowMax; fNegAdcTimeWindowMax = 0;
+  delete [] fPedNegDefault; fPedNegDefault = 0;
+  delete [] fPedPosDefault; fPedPosDefault = 0;
   delete [] fShPosPedLimit; fShPosPedLimit = 0;
   delete [] fShNegPedLimit; fShNegPedLimit = 0;
   delete [] fPosGain; fPosGain = 0;

src/THcShower.h

@@ -73,6 +73,17 @@ public:
     return ( Side == 0 ? fPosGain[nelem] : fNegGain[nelem]);
   }
 
+  Double_t GetPedDefault(Int_t NBlock, Int_t NLayer, Int_t Side) {
+    if (Side!=0&&Side!=1) {
+      cout << "*** Wrong Side in GetPedDefault:" << Side << " ***" << endl;
+      return -1;
+    }
+    Int_t nelem = 0;
+    for (Int_t i=0; i<NLayer; i++) nelem += fNBlocks[i];
+    nelem += NBlock;
+    return ( Side == 0 ? fPedPosDefault[nelem] : fPedNegDefault[nelem] );
+  }
+
   Double_t GetWindowMin(Int_t NBlock, Int_t NLayer, Int_t Side) {
     if (Side!=0&&Side!=1) {
       cout << "*** Wrong Side in GetWindowMin:" << Side << " ***" << endl;
@@ -189,7 +200,9 @@ protected:
   Double_t* fNegAdcTimeWindowMin;
   Double_t* fPosAdcTimeWindowMax;
   Double_t* fNegAdcTimeWindowMax;
-  Double_t fAdcTdcOffset;
+  Int_t* fPedPosDefault;
+  Int_t* fPedNegDefault;
+ Double_t fAdcTdcOffset;
 
   Int_t fAnalyzePedestals;   // Flag for pedestal analysis.
 

src/THcShowerArray.cxx

@@ -105,6 +105,7 @@ THcShowerArray::~THcShowerArray()
 
   delete [] fAdcTimeWindowMin; fAdcTimeWindowMin = 0;
   delete [] fAdcTimeWindowMax; fAdcTimeWindowMax = 0;
+  delete [] fPedDefault; fPedDefault = 0;
 }
 
 //_____________________________________________________________________________
@@ -280,6 +281,7 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
 
   fAdcTimeWindowMin = new Double_t [fNelem];
   fAdcTimeWindowMax = new Double_t [fNelem];
+  fPedDefault = new Int_t [fNelem];
 
   DBRequest list1[]={
     {"cal_arr_ped_limit", fPedLimit, kInt, static_cast<UInt_t>(fNelem),1},
@@ -287,12 +289,14 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
     {"cal_arr_gain_cor",  cal_arr_gain_cor,  kDouble, static_cast<UInt_t>(fNelem)},
     {"cal_arr_AdcTimeWindowMin", fAdcTimeWindowMin, kDouble, static_cast<UInt_t>(fNelem),1},
     {"cal_arr_AdcTimeWindowMax", fAdcTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem),1},
+    {"cal_arr_PedDefault", fPedDefault, kInt, static_cast<UInt_t>(fNelem),1},
     {0}
   };
 
    for(Int_t ip=0;ip<fNelem;ip++) {
     fAdcTimeWindowMin[ip] = -1000.;
     fAdcTimeWindowMax[ip] = 1000.;
+    fPedDefault[ip] = 0;
    }
 
   gHcParms->LoadParmValues((DBRequest*)&list1, prefix);
@@ -873,6 +877,8 @@ void THcShowerArray::FillADC_DynamicPedestal()
 {
   Double_t StartTime = 0.0;
   if( fglHod ) StartTime = fglHod->GetStartTime();
+   Double_t OffsetTime = 0.0;
+   if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
   for (Int_t ielem=0;ielem<frAdcPulseInt->GetEntries();ielem++) {
     
     Int_t npad           = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
@@ -881,15 +887,10 @@ void THcShowerArray::FillADC_DynamicPedestal()
     Double_t pulseInt    = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData();
     Double_t pulseAmp    = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
     Double_t pulseTime   = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
-    Double_t adctdcdiffTime = StartTime-pulseTime;
-    Bool_t errorflag     = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData();
+    Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
     Bool_t pulseTimeCut  = (adctdcdiffTime > fAdcTimeWindowMin[npad]) &&  (adctdcdiffTime < fAdcTimeWindowMax[npad]);
-
-    if (!errorflag)
-      {
 	fGoodAdcMult.at(npad) += 1;
-      }
-    if (!errorflag && pulseTimeCut) {
+    if (pulseTimeCut) {
       
       fTotNumAdcHits++;
       fGoodAdcPulseIntRaw.at(npad) = pulseIntRaw;
@@ -982,6 +983,24 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
       } else {
 	((THcSignalHit*) frAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(padnum,1);
       }
+
+      if (rawAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawAdcHit.GetAdcTomV();
+	if (fPedDefault[padnum-1] !=0) {
+	  Double_t tPulseInt = AdcToC*(rawAdcHit.GetPulseIntRaw(thit) - fPedDefault[padnum-1]*PeakPedRatio);
+	  ((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(padnum, tPulseInt);
+          ((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(padnum, fPedDefault[padnum-1]);
+          ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(padnum, float(fPedDefault[padnum-1])/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(padnum, 0.);
+	
+      }
+
+
       ++nrAdcHits;
     }
     ihit++;

src/THcShowerArray.h

@@ -150,6 +150,7 @@ protected:
   static const Int_t kADCSampIntDynPed=3;
   Double_t *fAdcTimeWindowMin ;
   Double_t *fAdcTimeWindowMax ;
+  Int_t *fPedDefault ;
   Double_t fAdcThreshold ;
   Double_t fAdcTdcOffset;
 

src/THcShowerHit.h

@@ -7,6 +7,7 @@
 #include <iterator>
 #include <iostream>
 #include <memory>
+#include <vector>
 #include "TMath.h"
 
 using namespace std;

src/THcShowerPlane.cxx

@@ -612,6 +612,21 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
       } else {
 	((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,1);
       }
+     if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawPosAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawPosAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawPosAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawPosAdcHit.GetAdcTomV();
+	Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,0);
+	if (PedDefaultTemp !=0) {
+	  Double_t tPulseInt = AdcToC*(rawPosAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
+	  ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(padnum, tPulseInt);
+          ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, PedDefaultTemp);
+          ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, 0.);	
+      }
       ++nrPosAdcHits;
       fTotNumAdcHits++;
       fTotNumPosAdcHits++;
@@ -637,6 +652,21 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
       } else {
 	((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,1);
       }
+     if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawNegAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawNegAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawNegAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawNegAdcHit.GetAdcTomV();
+	Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,1);
+	if (PedDefaultTemp !=0) {
+	  Double_t tPulseInt = AdcToC*(rawNegAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
+	  ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(padnum, tPulseInt);
+          ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, PedDefaultTemp);
+          ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, 0.);	
+      }
       ++nrNegAdcHits;
       fTotNumAdcHits++;
       fTotNumNegAdcHits++;
@@ -744,6 +774,8 @@ void THcShowerPlane::FillADC_DynamicPedestal()
 {
   Double_t StartTime = 0.0;
   if( fglHod ) StartTime = fglHod->GetStartTime();
+   Double_t OffsetTime = 0.0;
+   if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
   for (Int_t ielem=0;ielem<frNegAdcPulseInt->GetEntries();ielem++) {
    Int_t    npad         = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
    Double_t pulseInt     = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
@@ -751,18 +783,11 @@ void THcShowerPlane::FillADC_DynamicPedestal()
     Double_t pulseAmp     = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
     Double_t pulseIntRaw  = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
     Double_t pulseTime    = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
-    Double_t adctdcdiffTime = StartTime-pulseTime;
+    Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
     Double_t threshold    = ((THcSignalHit*) frNegAdcThreshold->ConstructedAt(ielem))->GetData();
-    Bool_t   errorflag    = ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
     Bool_t   pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,1)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,1) );
-
-    
- 
-    if (!errorflag)
-      {
 	fGoodNegAdcMult.at(npad) += 1;
-      }
-    if (!errorflag && pulseTimeCut) {
+    if (pulseTimeCut) {
       fGoodNegAdcPulseIntRaw.at(npad) =pulseIntRaw;
 
       if(fGoodNegAdcPulseIntRaw.at(npad) >  threshold && fGoodNegAdcPulseInt.at(npad)==0) {
@@ -793,18 +818,10 @@ void THcShowerPlane::FillADC_DynamicPedestal()
     Double_t pulseInt     = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData();
     Double_t pulseIntRaw  = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
     Double_t pulseTime    = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
-     Double_t adctdcdiffTime = StartTime-pulseTime;
-   Bool_t   errorflag    = ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
+     Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
    Bool_t   pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,0)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,0) );
-
-
-
-    if (!errorflag)
-      {
 	fGoodPosAdcMult.at(npad) += 1;
-      }
-    
-    if (!errorflag && pulseTimeCut) {
+    if (pulseTimeCut) {
       fGoodPosAdcPulseIntRaw.at(npad) = pulseIntRaw;
 
       if(fGoodPosAdcPulseIntRaw.at(npad) >  threshold && fGoodPosAdcPulseInt.at(npad)==0) {

src/THcSpacePoint.h

@@ -29,7 +29,7 @@ public:
   };
 
   void SetXY(Double_t x, Double_t y) {fX = x; fY = y;};
-  void Clear(Option_t* opt="") {fNHits=0; fNCombos=0; fHits.clear();};
+  void Clear(Option_t* /* opt */ ="") {fNHits=0; fNCombos=0; fHits.clear();};
   void AddHit(THcDCHit* hit) {
     Hit newhit;
     newhit.dchit = hit;

src/THcTrigDet.cxx

@@ -211,7 +211,8 @@ void THcTrigDet::Clear(Option_t* opt) {
   THaAnalysisObject::Clear(opt);
 
   // Reset all data.
-  for (int i=0; i<fNumAdc; ++i) {
+  fTdcRefTime = kBig;
+ for (int i=0; i<fNumAdc; ++i) {
     fAdcPedRaw[i] = 0;
     fAdcPulseIntRaw[i] = 0;
     fAdcPulseAmpRaw[i] = 0;
@@ -272,14 +273,21 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
     }
     else if (hit->fPlane == 2) {
       THcRawTdcHit rawTdcHit = hit->GetRawTdcHit();
-
+    if (rawTdcHit.GetNHits() >0 && rawTdcHit.HasRefTime() && fTdcRefTime == kBig) fTdcRefTime=rawTdcHit.GetRefTime() ;
       UInt_t good_hit=999;
+      UInt_t closest_hit=999;
+      Int_t TimeDiff=1000000;
            for (UInt_t thit=0; thit<rawTdcHit.GetNHits(); ++thit) {
 	    Int_t TestTime= rawTdcHit.GetTimeRaw(thit);
+	    if (abs(TestTime-fTdcTimeWindowMin[cnt]) < TimeDiff) {
+	      closest_hit=thit;
+	      TimeDiff=abs(TestTime-fTdcTimeWindowMin[cnt]);
+	    }
 	    if (TestTime>=fTdcTimeWindowMin[cnt]&&TestTime<=fTdcTimeWindowMax[cnt]&&good_hit==999) {
 	      good_hit=thit;
 	    }
 	   }
+	   if (good_hit == 999 and closest_hit != 999) good_hit=closest_hit;
 	 if (good_hit<rawTdcHit.GetNHits()) {
       fTdcTimeRaw[cnt] = rawTdcHit.GetTimeRaw(good_hit);
       fTdcTime[cnt] = rawTdcHit.GetTime(good_hit)*fTdcChanperNS+fTdcOffset;
@@ -396,7 +404,16 @@ Int_t THcTrigDet::DefineVariables(THaAnalysisObject::EMode mode) {
   std::vector<TString> adcPulseIntTitle(fNumAdc), adcPulseIntVar(fNumAdc);
   std::vector<TString> adcPulseAmpTitle(fNumAdc), adcPulseAmpVar(fNumAdc);
   std::vector<TString> adcMultiplicityTitle(fNumAdc), adcMultiplicityVar(fNumAdc);
-
+  
+  TString RefTimeTitle= "TdcRefTime";
+   TString RefTimeVar= "fTdcRefTime";
+   RVarDef entryRefTime {
+      RefTimeTitle.Data(),
+      RefTimeTitle.Data(),
+      RefTimeVar.Data()
+    };
+     vars.push_back(entryRefTime);
+ 
   for (int i=0; i<fNumAdc; ++i) {
     adcPedRawTitle.at(i) = fAdcNames.at(i) + "_adcPedRaw";
     adcPedRawVar.at(i) = TString::Format("fAdcPedRaw[%d]", i);

src/THcTrigDet.h

@@ -82,6 +82,7 @@ class THcTrigDet : public THaDetector, public THcHitList {
 
     Int_t fTdcMultiplicity[fMaxTdcChannels];
     Int_t fAdcMultiplicity[fMaxAdcChannels];
+    Double_t fTdcRefTime;
 
     TString fSpectName;
     std::vector<Int_t> eventtypes;

src/include/HallC_LinkDef.h

 // Preamble to HallC_LinkDef.h file
 
 #ifdef __CINT__
- 
+
 #pragma link off all globals;
 #pragma link off all classes;
 #pragma link off all functions;
- 
+
 #pragma link C++ nestedclass;
 #pragma link C++ nestedtypedef;
 
@@ -19,14 +19,12 @@
 #pragma link C++ class hallc::data::DriftChamber+;
 
 #pragma link C++ class hallc::data::PulseWaveForm+;
-#pragma link C++ class std::vector<hallc::data::PulseWaveForm>+;
-
+#pragma link C++ class std::vector < hallc::data::PulseWaveForm>+;
 
-#pragma link C++ class podd2::HitLogging<TObject>+;
+#pragma link C++ class podd2::HitLogging < TObject>+;
 
 #pragma link C++ global gHcParms;
 #pragma link C++ global gHcDetectorMap;
- 
 
 #pragma link C++ class Decoder::Scaler9001+;
 #pragma link C++ class Decoder::Scaler9250+;
@@ -55,6 +53,7 @@
 #pragma link C++ class THcHodoEff+;
 #pragma link C++ class THcHelicity+;
 #pragma link C++ class THcHelicityReader+;
+#pragma link C++ class THcHelicityScaler+;
 #pragma link C++ class THcHodoHit+;
 #pragma link C++ class THcHodoscope+;
 #pragma link C++ class THcInterface+;
@@ -93,6 +92,5 @@
 #pragma link C++ class hcana::Scandalizer+;
 #pragma link C++ class hcana::TrackingEfficiency+;
 
-#pragma link C++ class THcOnlineRun+;
 // Postamble for HallC_Linkdef.h file
 #endif

src/include/hcana/HallC_Data.h

@@ -15,7 +15,7 @@ namespace hallc {
       static const UInt_t MaxNSamples = 511;
       // From THcRawAdcHit.h
       PulseWaveForm() {}
-      PulseWaveForm(Int_t* buf, Int_t size = MaxNSamples) { std::copy_n(buf, MaxNSamples, std::begin(_buffer)); }
+      PulseWaveForm(Int_t* buf, Int_t size = MaxNSamples) { std::copy_n(buf, size, std::begin(_buffer)); }
       virtual ~PulseWaveForm() {}
 
       void ZeroBuffer() { std::fill(std::begin(_buffer), std::end(_buffer), 0); }

tests/README.md

-Hcana tests
-===========
-
-## Current tests:
-
-- ep elastic tests
-
-## Tests to add:
-
-- Start time check
-- Focal plane times test
-- "Tracking efficiency" test
-- Invariant mass check (jpsi)

tests/elastic_coin_replay.cxx

-#include <fmt/core.h>
-#include <fmt/ostream.h>
-#include <vector>
-
-#include "TString.h"
-
-R__LOAD_LIBRARY(libHallC.so)
-#include "hcana/HallC_Data.h"
-#include "DecData.h"
-//R__LOAD_LIBRARY(libScandalizer.so)
-//#include "monitor/DetectorDisplay.h"
-//#include "monitor/DisplayPlots.h"
-//#include "monitor/MonitoringDisplay.h"
-//#include "scandalizer/PostProcessors.h"
-//#include "scandalizer/ScriptHelpers.h"
-//
-//#include "THaPostProcess.h"
-//#include "monitor/ExperimentMonitor.h"
-//#include "scandalizer/PostProcessors.h"
-#include "THcAnalyzer.h"
-#include "THaCut.h"
-#include "THcGlobals.h"
-#include "THcHallCSpectrometer.h"
-#include "THcDetectorMap.h"
-#include "THcCherenkov.h"
-#include "THcDC.h"
-#include "THcHodoscope.h"
-#include "THcParmList.h"
-#include "THaGoldenTrack.h"
-#include "THcHodoEff.h"
-#include "THcScalerEvtHandler.h"
-#include "THcShower.h"
-#include "THcReactionPoint.h"
-#include "THcExtTarCor.h"
-#include "THcRasteredBeam.h"
-#include "THcRun.h"
-#include "THcCoinTime.h"
-#include "THcConfigEvtHandler.h"
-#include "THcTrigDet.h"
-#include "THcTrigApp.h"
-#include "THcSecondaryKine.h"
-#include "THcAerogel.h"
-#include "THcPrimaryKine.h"
-#include "THaReactionPoint.h"
-
-
-
-void elastic_coin_replay(Int_t RunNumber = 0, Int_t MaxEvent = -1) {
-  using namespace std;
-
-  // Get RunNumber and MaxEvent if not provided.
-  if( RunNumber<=0 ) {
-    std::exit(-1);
-  }
-
-  // Create file name patterns.
-  const char* RunFileNamePattern = "coin_all_%05d.dat";
-  vector<TString> pathList;
-  pathList.push_back(".");
-  pathList.push_back("./raw");
-  pathList.push_back("./raw/../raw.copiedtotape");
-  pathList.push_back("./cache");
-
-  //const char* RunFileNamePattern = "raw/coin_all_%05d.dat";
-  const char* ROOTFileNamePattern = "ROOTfiles/coin_replay_production_%d_%d.root";
-  
-  // Load global parameters
-  gHcParms->Define("gen_run_number", "Run Number", RunNumber);
-  gHcParms->AddString("g_ctp_database_filename", "DBASE/COIN/standard.database");
-  gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
-  gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
-  gHcParms->Load(gHcParms->GetString("g_ctp_kinematics_filename"), RunNumber);
-  // Load params for COIN trigger configuration
-  gHcParms->Load("PARAM/TRIG/tcoin.param");
-  // Load fadc debug parameters
-  gHcParms->Load("PARAM/HMS/GEN/h_fadc_debug.param");
-  gHcParms->Load("PARAM/SHMS/GEN/p_fadc_debug.param");
-
-  // const char* CurrentFileNamePattern = "low_curr_bcm/bcmcurrent_%d.param";
-  // gHcParms->Load(Form(CurrentFileNamePattern, RunNumber));
-
-  // Load the Hall C detector map
-  gHcDetectorMap = new THcDetectorMap();
-  gHcDetectorMap->Load("MAPS/COIN/DETEC/coin.map");
-
-     // Dec data
-   gHaApps->Add(new Podd::DecData("D","Decoder raw data"));
-  //=:=:=:=
-  // SHMS 
-  //=:=:=:=
-
-  // Set up the equipment to be analyzed.
-  THcHallCSpectrometer* SHMS = new THcHallCSpectrometer("P", "SHMS");
-  SHMS->SetEvtType(1);
-  SHMS->AddEvtType(4);
-  SHMS->AddEvtType(5);
-  SHMS->AddEvtType(6);
-  SHMS->AddEvtType(7);
-  gHaApps->Add(SHMS);
-  // Add Noble Gas Cherenkov to SHMS apparatus
-  THcCherenkov* pngcer = new THcCherenkov("ngcer", "Noble Gas Cherenkov");
-  SHMS->AddDetector(pngcer);
-  // Add drift chambers to SHMS apparatus
-  THcDC* pdc = new THcDC("dc", "Drift Chambers");
-  SHMS->AddDetector(pdc);
-  // Add hodoscope to SHMS apparatus
-  THcHodoscope* phod = new THcHodoscope("hod", "Hodoscope");
-  SHMS->AddDetector(phod);
-  // Add Heavy Gas Cherenkov to SHMS apparatus
-  THcCherenkov* phgcer = new THcCherenkov("hgcer", "Heavy Gas Cherenkov");
-  SHMS->AddDetector(phgcer);
-  // Add Aerogel Cherenkov to SHMS apparatus
-  THcAerogel* paero = new THcAerogel("aero", "Aerogel");
-  SHMS->AddDetector(paero);
-  // Add calorimeter to SHMS apparatus
-  THcShower* pcal = new THcShower("cal", "Calorimeter");
-  SHMS->AddDetector(pcal);
-
-  // THcBCMCurrent* hbc = new THcBCMCurrent("H.bcm", "BCM current check");
-  // gHaPhysics->Add(hbc);
-
-  // Add rastered beam apparatus
-  THaApparatus* pbeam = new THcRasteredBeam("P.rb", "Rastered Beamline");
-  gHaApps->Add(pbeam);
-  // Add physics modules
-  // Calculate reaction point
-  THcReactionPoint* prp = new THcReactionPoint("P.react", "SHMS reaction point", "P", "P.rb");
-  gHaPhysics->Add(prp);
-  // Calculate extended target corrections
-  THcExtTarCor* pext = new THcExtTarCor("P.extcor", "HMS extended target corrections", "P", "P.react");
-  gHaPhysics->Add(pext);
-  // Calculate golden track quantites
-  THaGoldenTrack* pgtr = new THaGoldenTrack("P.gtr", "SHMS Golden Track", "P");
-  gHaPhysics->Add(pgtr);
-  // Calculate the hodoscope efficiencies
-  THcHodoEff* peff = new THcHodoEff("phodeff", "SHMS hodo efficiency", "P.hod");
-  gHaPhysics->Add(peff);   
-
-  // Add event handler for scaler events
-  THcScalerEvtHandler* pscaler = new THcScalerEvtHandler("P", "Hall C scaler event type 1");
-  pscaler->AddEvtType(1);
-  pscaler->AddEvtType(4);
-  pscaler->AddEvtType(5);
-  pscaler->AddEvtType(6);
-  pscaler->AddEvtType(7);
-  pscaler->AddEvtType(129);
-  pscaler->SetDelayedType(129);
-  pscaler->SetUseFirstEvent(kTRUE);
-  gHaEvtHandlers->Add(pscaler);
-
-  //=:=:=
-  // HMS 
-  //=:=:=
-
-  // Set up the equipment to be analyzed.
-  THcHallCSpectrometer* HMS = new THcHallCSpectrometer("H", "HMS");
-  HMS->SetEvtType(2);
-  HMS->AddEvtType(4);
-  HMS->AddEvtType(5);
-  HMS->AddEvtType(6);
-  HMS->AddEvtType(7);
-  gHaApps->Add(HMS);
-  // Add drift chambers to HMS apparatus
-  THcDC* hdc = new THcDC("dc", "Drift Chambers");
-  HMS->AddDetector(hdc);
-  // Add hodoscope to HMS apparatus
-  THcHodoscope* hhod = new THcHodoscope("hod", "Hodoscope");
-  HMS->AddDetector(hhod);
-  // Add Cherenkov to HMS apparatus
-  THcCherenkov* hcer = new THcCherenkov("cer", "Heavy Gas Cherenkov");
-  HMS->AddDetector(hcer);
-  // Add Aerogel Cherenkov to HMS apparatus
-  // THcAerogel* haero = new THcAerogel("aero", "Aerogel");
-  // HMS->AddDetector(haero);
-  // Add calorimeter to HMS apparatus
-  THcShower* hcal = new THcShower("cal", "Calorimeter");
-  HMS->AddDetector(hcal);
-
-  // Add rastered beam apparatus
-  THaApparatus* hbeam = new THcRasteredBeam("H.rb", "Rastered Beamline");
-  gHaApps->Add(hbeam);  
-  // Add physics modules
-  // Calculate reaction point
-  THcReactionPoint* hrp = new THcReactionPoint("H.react", "HMS reaction point", "H", "H.rb");
-  gHaPhysics->Add(hrp);
-  // Calculate extended target corrections
-  THcExtTarCor* hext = new THcExtTarCor("H.extcor", "HMS extended target corrections", "H", "H.react");
-  gHaPhysics->Add(hext);
-  // Calculate golden track quantities
-  THaGoldenTrack* hgtr = new THaGoldenTrack("H.gtr", "HMS Golden Track", "H");
-  gHaPhysics->Add(hgtr);
-  // Calculate the hodoscope efficiencies
-  THcHodoEff* heff = new THcHodoEff("hhodeff", "HMS hodo efficiency", "H.hod");
-  gHaPhysics->Add(heff);
-
-  // Add event handler for scaler events
-  THcScalerEvtHandler *hscaler = new THcScalerEvtHandler("H", "Hall C scaler event type 4");  
-  hscaler->AddEvtType(2);
-  hscaler->AddEvtType(4);
-  hscaler->AddEvtType(5);
-  hscaler->AddEvtType(6);
-  hscaler->AddEvtType(7);
-  hscaler->AddEvtType(129);
-  hscaler->SetDelayedType(129);
-  hscaler->SetUseFirstEvent(kTRUE);
-  gHaEvtHandlers->Add(hscaler);
-
-  //=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
-  // Kinematics Modules
-  //=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
-
-  // Add Physics Module to calculate primary (scattered electrons) beam kinematics
-  THcPrimaryKine* hkin_primary = new THcPrimaryKine("H.kin.primary", "HMS Single Arm Kinematics", "H", "H.rb");
-  gHaPhysics->Add(hkin_primary);
-  // Add Physics Module to calculate secondary (scattered hadrons) beam kinematics
-  THcSecondaryKine* pkin_secondary = new THcSecondaryKine("P.kin.secondary", "SHMS Single Arm Kinematics", "P", "H.kin.primary");
-  gHaPhysics->Add(pkin_secondary);
-  
-  //=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
-  // Global Objects & Event Handlers
-  //=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
-
-  // Add trigger apparatus
-  THaApparatus* TRG = new THcTrigApp("T", "TRG");
-  gHaApps->Add(TRG);
-  // Add trigger detector to trigger apparatus
-  THcTrigDet* coin = new THcTrigDet("coin", "Coincidence Trigger Information");
-  // Suppress missing reference time warnings for these event types
-  coin->SetEvtType(1);
-  coin->AddEvtType(2);
-  TRG->AddDetector(coin); 
-
-  
-  //Add coin physics module THcCoinTime::THcCoinTime (const char *name, const char* description, const char* hadArmName, 
-  // const char* elecArmName, const char* coinname) :
-  THcCoinTime* coinTime = new THcCoinTime("CTime", "Coincidende Time Determination", "P", "H", "T.coin");
-  gHaPhysics->Add(coinTime);
-
-  // Add event handler for prestart event 125.
-  THcConfigEvtHandler* ev125 = new THcConfigEvtHandler("HC", "Config Event type 125");
-  gHaEvtHandlers->Add(ev125);
-  // Add event handler for EPICS events
-  THaEpicsEvtHandler* hcepics = new THaEpicsEvtHandler("epics", "HC EPICS event type 180");
-  gHaEvtHandlers->Add(hcepics);
- 
-  // Set up the analyzer - we use the standard one,
-  // but this could be an experiment-specific one as well.
-  // The Analyzer controls the reading of the data, executes
-  // tests/cuts, loops over Acpparatus's and PhysicsModules,
-  // and executes the output routines.
-  THcAnalyzer* analyzer = new THcAnalyzer;
-
-  // A simple event class to be output to the resulting tree.
-  // Creating your own descendant of THaEvent is one way of
-  // defining and controlling the output.
-  THaEvent* event = new THaEvent;
-
-  // Define the run(s) that we want to analyze.
-  // We just set up one, but this could be many.
-  THcRun* run = new THcRun( pathList, Form(RunFileNamePattern, RunNumber) );
-
-  // Set to read in Hall C run database parameters
-  run->SetRunParamClass("THcRunParameters");
-  
-  // Eventually need to learn to skip over, or properly analyze the pedestal events
-  run->SetEventRange(1, MaxEvent); // Physics Event number, does not include scaler or control events.
-  run->SetNscan(1);
-  run->SetDataRequired(0x7);
-  run->Print();
-
-  // Define the analysis parameters
-  TString ROOTFileName = Form(ROOTFileNamePattern, RunNumber, MaxEvent);
-  analyzer->SetCountMode(2);  // 0 = counter is # of physics triggers
-                              // 1 = counter is # of all decode reads
-                              // 2 = counter is event number
-
-  analyzer->SetEvent(event);
-  // Set EPICS event type
-  analyzer->SetEpicsEvtType(180);
-  // Define crate map
-  analyzer->SetCrateMapFileName("MAPS/db_cratemap.dat");
-  // Define output ROOT file
-  analyzer->SetOutFile(ROOTFileName.Data());
-  // Define DEF-file+
-  analyzer->SetOdefFile("DEF-files/COIN/PRODUCTION/coin_production_hElec_pProt.def");
-  // Define cuts file
-  analyzer->SetCutFile("DEF-files/COIN/PRODUCTION/CUTS/coin_production_cuts.def");  // optional
-  // File to record accounting information for cuts
-  analyzer->SetSummaryFile(Form("REPORT_OUTPUT/COIN/PRODUCTION/summary_production_%d_%d.report", RunNumber, MaxEvent));  // optional
-  // Start the actual analysis.
-  analyzer->Process(run);
-  // Create report file from template
-  analyzer->PrintReport("TEMPLATES/COIN/PRODUCTION/coin_production.template",
-  			Form("REPORT_OUTPUT/COIN/PRODUCTION/replay_coin_production_%d_%d.report", RunNumber, MaxEvent));  // optional
-
-}

tests/elastic_test.cxx

-#include "nlohmann/json.hpp"
-#include <cmath>
-#include <iostream>
-
-#include "ROOT/RDataFrame.hxx"
-#include "ROOT/RVec.hxx"
-
-#include "Math/Vector3D.h"
-#include "Math/Vector4D.h"
-#include "Math/VectorUtil.h"
-#include "TCanvas.h"
-#include "TLatex.h"
-#include "TStyle.h"
-#include "TSystem.h"
-R__LOAD_LIBRARY(libMathMore.so)
-R__LOAD_LIBRARY(libGenVector.so)
-
-R__LOAD_LIBRARY(libfmt.so)
-#include "fmt/core.h"
-
-#include "THStack.h"
-
-#ifdef __cpp_lib_filesystem
-#include <filesystem>
-namespace fs = std::filesystem;
-#else
-#include <experimental/filesystem>
-namespace fs = std::experimental::filesystem;
-#endif
-
-using RDFNode = decltype(ROOT::RDataFrame{0}.Filter(""));
-using Histo1DProxy =
-    decltype(ROOT::RDataFrame{0}.Histo1D(ROOT::RDF::TH1DModel{"", "", 128u, 0., 0.}, ""));
-
-struct RDFInfo {
-  RDFNode&          df;
-  const std::string title;
-  RDFInfo(RDFNode& df, std::string_view title) : df{df}, title{title} {}
-};
-
-constexpr const double M_P     = .938272;
-constexpr const double M_P2    = M_P * M_P;
-constexpr const double M_pion  = 0.139;
-constexpr const double M_pion2 = M_pion * M_pion;
-constexpr const double M_e     = .000511;
-
-// =================================================================================
-// Cuts
-// =================================================================================
-std::string goodTrackSHMS = "P.gtr.dp > -10 && P.gtr.dp < 22";
-std::string goodTrackHMS  = "H.gtr.dp > -8 && H.gtr.dp < 8";
-
-std::string piCutSHMS = "P.cal.etottracknorm<1.0";
-//std::string piCutSHMS = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && P.cal.etottracknorm<1.0";
-
-std::string eCutHMS = "H.cal.etottracknorm > 0.50 && H.cal.etottracknorm < 2. && "
-                      "H.cer.npeSum > 1.";
-
-std::string epiCut = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && "
-                     "H.cer.npeSum > 1.0 && H.cal.etottracknorm > 0.6 && "
-                     "H.cal.etottracknorm < 2.0 && P.cal.etottracknorm<1.0";
-
-using Pvec3D = ROOT::Math::XYZVector;
-using Pvec4D = ROOT::Math::PxPyPzMVector;
-
-// =================================================================================
-// reconstruction
-// =================================================================================
-auto p_pion = [](double px, double py, double pz) {
-  return Pvec4D{px * 0.996, py * 0.996, pz * 0.996, M_e};
-};
-auto p_electron = [](double px, double py, double pz) {
-  return Pvec4D{px * 0.994, py * 0.994, pz * 0.994, M_e};
-};
-auto t = [](const double Egamma, Pvec4D& jpsi) {
-  Pvec4D beam{0, 0, Egamma, 0};
-  return (beam - jpsi).M2();
-};
-
-bool root_file_exists(std::string rootfile) {
-  bool found_good_file = false;
-  if (!gSystem->AccessPathName(rootfile.c_str())) {
-    TFile file(rootfile.c_str());
-    if (file.IsZombie()) {
-      std::cout << rootfile << " is a zombie.\n";
-      std::cout
-          << " Did your replay finish?  Check that the it is done before running this script.\n";
-      return false;
-      // return;
-    } else {
-      std::cout << " using : " << rootfile << "\n";
-      return true;
-    }
-  }
-  return false;
-}
-
-void elastic_test(int RunNumber = 6012, int nevents = 50000, int prompt = 0, int update = 0,
-                        int default_count_goal = 10000, int redo_timing = 0) {
-
-  // ===============================================================================================
-  // Initialization
-  // ===============================================================================================
-  using nlohmann::json;
-  json j;
-  {
-    std::ifstream json_input_file("db2/run_list_coin.json");
-    try {
-      json_input_file >> j;
-    } catch (json::parse_error) {
-      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
-      std::quick_exit(-127);
-    }
-  }
-
-  auto runnum_str = std::to_string(RunNumber);
-  if (j.find(runnum_str) == j.end()) {
-    std::cout << "Run " << RunNumber << " not found in db2/run_list_coin.json\n";
-    std::cout << "Check that run number and replay exists. \n";
-    std::cout << "If problem persists please contact Sylvester (217-848-0565)\n";
-  }
-  double P0_shms_setting = j[runnum_str]["spectrometers"]["shms_momentum"].get<double>();
-  double P0_shms         = std::abs(P0_shms_setting);
-
-  bool found_good_file = false;
-
-  std::string rootfile =
-      fmt::format("full_online/coin_replay_production_{}_{}.root", RunNumber, nevents);
-  found_good_file = root_file_exists(rootfile.c_str());
-  if (!found_good_file) {
-    rootfile =
-        fmt::format("ROOTfiles_volatile/coin_replay_production_{}_{}.root", RunNumber, nevents);
-    found_good_file = root_file_exists(rootfile.c_str());
-  }
-  if (!found_good_file) {
-    rootfile = fmt::format("ROOTfiles_csv/coin_replay_production_{}_{}.root", RunNumber, nevents);
-    found_good_file = root_file_exists(rootfile.c_str());
-  }
-  if (!found_good_file) {
-    rootfile = fmt::format("ROOTfiles/coin_replay_production_{}_{}.root", RunNumber, nevents);
-    found_good_file = root_file_exists(rootfile.c_str());
-  }
-  if (!found_good_file) {
-    std::cout << " Error: suitable root file not found\n";
-    return;
-  }
-
-  // ===============================================================================================
-  // Dataframe
-  // ===============================================================================================
-
-  ROOT::EnableImplicitMT(24);
-
-  //---------------------------------------------------------------------------
-  // Detector tree
-  ROOT::RDataFrame d("T", rootfile);
-
-  //// SHMS Scaler tree
-  //ROOT::RDataFrame d_sh("TSP", rootfile);
-  //// int N_scaler_events = *(d_sh.Count());
-
-  auto d_coin = d.Filter("fEvtHdr.fEvtType == 4");
-
-  // Good track cuts
-  auto dHMSGoodTrack  = d_coin.Filter(goodTrackHMS);
-  auto dSHMSGoodTrack = d_coin.Filter(goodTrackSHMS);
-  auto dCOINGoodTrack = dHMSGoodTrack.Filter(goodTrackSHMS)
-                            .Define("p_electron", p_electron, {"H.gtr.px", "H.gtr.py", "H.gtr.pz"})
-                            .Define("p_pion", p_pion, {"P.gtr.px", "P.gtr.py", "P.gtr.pz"});
-  // PID cuts
-  auto dHMSEl  = dHMSGoodTrack.Filter(eCutHMS);
-  auto dSHMSEl = dSHMSGoodTrack.Filter(piCutSHMS);
-  auto dCOINEl = dCOINGoodTrack.Filter(eCutHMS + " && " + piCutSHMS);
-                     //.Filter(
-                     //    [=](double npe, double dp) {
-                     //      double p_track = P0_shms * (100.0 + dp) / 100.0;
-                     //      // no cerenkov cut needed when momentum is below 2.8 GeV/c
-                     //      if (p_track < 2.8) {
-                     //        return true;
-                     //      }
-                     //      return npe > 1.0;
-                     //    },
-                     //    {"P.hgcer.npeSum", "P.gtr.dp"});
-
-  // Timing cuts
-  // Find the timing peak
-  // Find the coin peak
-  double coin_peak_center = 0;
-  if (redo_timing) {
-    auto h_coin_time =
-        dCOINEl.Histo1D({"coin_time", "coin_time", 8000, 0, 1000}, "CTime.ePositronCoinTime_ROC2");
-    h_coin_time->DrawClone();
-    int coin_peak_bin = h_coin_time->GetMaximumBin();
-    coin_peak_center  = h_coin_time->GetBinCenter(coin_peak_bin);
-    std::cout << "COINCIDENCE time peak found at: " << coin_peak_center << std::endl;
-  } else {
-    //coin_peak_center = 43.0; // run 7240-7241
-    coin_peak_center = 23.0; // run 6012
-    std::cout << "COINCIDENCE time peak: using pre-calculated value at: " << coin_peak_center
-              << std::endl;
-    ;
-  }
-  // timing cut lambdas
-  // TODO: evaluate timing cut and offset for random background
-  auto timing_cut = [=](double coin_time) { return std::abs(coin_time - coin_peak_center) < 2.; };
-  // anti-timing set to 5x width of regular
-  auto anti_timing_cut = [=](double coin_time) {
-    return std::abs(coin_time - coin_peak_center - 28.) < 10.;
-  };
-
-  // timing counts
-  auto dHMSElInTime  = dHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
-  auto dHMSElRandom  = dHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});
-  auto dSHMSElInTime = dSHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
-  auto dSHMSElRandom = dSHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});
-
-  auto dCOINElInTime = dCOINEl.Filter(timing_cut, {"CTime.ePiCoinTime_ROC2"});
-  auto dCOINElRandom = dCOINEl.Filter(anti_timing_cut, {"CTime.ePiCoinTime_ROC2"});
-
-  // Output root file
-  //auto out_file =
-  //    new TFile(fmt::format("monitoring/{}/good_csv_counter.root", RunNumber).c_str(), "UPDATE");
-  //out_file->cd();
-
-  // =========================================================================================
-  // Histograms
-  // =========================================================================================
-  // 2D correlations
-  auto hTheta2DNoCuts = d_coin.Histo2D(
-      {"theta2D", "No cuts;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
-      "P.gtr.th", "H.gtr.th");
-  auto hTheta2DTracking = dCOINGoodTrack.Histo2D(
-      {"theta2D", "Cuts: tracking;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
-      "P.gtr.th", "H.gtr.th");
-  auto hTheta2DPID =
-      dCOINEl.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1,
-                       .1, 50, -.1, .1},
-                      "P.gtr.th", "H.gtr.th");
-  auto hTheta2DTiming =
-      dCOINElInTime.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50,
-                             -.1, .1, 50, -.1, .1},
-                            "P.gtr.th", "H.gtr.th");
-  // timing
-  auto hCoinTimeNoCuts =
-      d_coin.Histo1D({"coin_time.NoCuts", "No Cuts;coin_time;counts", 8000, 0, 1000},
-                     "CTime.ePositronCoinTime_ROC2");
-  auto hCoinTimeTracking = dCOINGoodTrack.Histo1D(
-      {"coin_time.Tracking", "Cuts: Tracking;coin_time;counts", 8000, 0, 1000},
-      "CTime.ePositronCoinTime_ROC2");
-  auto hCoinTimePID =
-      dCOINEl.Histo1D({"coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
-                      "CTime.ePositronCoinTime_ROC2");
-  auto hCoinTimeTiming = dCOINElInTime.Histo1D(
-      {"coin_time.Timing", "Cuts: Tracking+PID+Timing;coin_time;counts", 8000, 0, 1000},
-      "CTime.ePositronCoinTime_ROC2");
-
-  auto hRandCoinTimePID = dCOINElRandom.Histo1D(
-      {"rand_coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
-      "CTime.ePositronCoinTime_ROC2");
-
-  // P.gtr.dp
-  auto hPdpNoCuts =
-      d_coin.Histo1D({"P.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
-  auto hPdpTracking = dSHMSGoodTrack.Histo1D(
-      {"P.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
-  auto hPdpPID = dSHMSEl.Histo1D(
-      {"P.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
-  auto hPdpTiming = dSHMSElInTime.Histo1D(
-      {"P.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
-      "P.gtr.dp");
-  // P.gtr.th
-  auto hPthNoCuts = d_coin.Histo1D(
-      {"P.gtr.th.NoCuts", "No Cuts;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
-  auto hPthTracking = dSHMSGoodTrack.Histo1D(
-      {"P.gtr.th.Tracking", "Cuts: Tracking;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
-  auto hPthPID = dSHMSEl.Histo1D(
-      {"P.gtr.th.PID", "Cuts: Tracking+PID;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
-  auto hPthTiming = dSHMSElInTime.Histo1D(
-      {"P.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{SHMS};counts", 200, -0.1, 0.1},
-      "P.gtr.th");
-  // P.gtr.ph
-  auto hPphNoCuts =
-      d_coin.Histo1D({"P.gtr.ph.NoCuts", "No Cuts;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
-  auto hPphTracking = dSHMSGoodTrack.Histo1D(
-      {"P.gtr.ph.Tracking", "Cuts: Tracking;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
-  auto hPphPID = dSHMSEl.Histo1D(
-      {"P.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
-  auto hPphTiming = dSHMSElInTime.Histo1D(
-      {"P.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{SHMS};counts", 200, -0.1, 0.1},
-      "P.gtr.ph");
-  // P.gtr.y
-  auto hPyNoCuts =
-      d_coin.Histo1D({"P.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "P.gtr.y");
-  auto hPyTracking = dSHMSGoodTrack.Histo1D(
-      {"P.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "P.gtr.y");
-  auto hPyPID =
-      dSHMSEl.Histo1D({"P.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "P.gtr.y");
-  auto hPyTiming = dSHMSElInTime.Histo1D(
-      {"P.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "P.gtr.y");
-  // P.cal.etottracknorm
-  auto hPcalEPNoCuts =
-      d_coin.Histo1D({"P.cal.etottracknorm.NoCuts", "No Cuts;SHMS E/P;counts", 200, -.5, 1.5},
-                     "P.cal.etottracknorm");
-  auto hPcalEPTracking = dSHMSGoodTrack.Histo1D(
-      {"P.cal.etottracknorm.Tracking", "Cuts: Tracking;SHMS E/P;counts", 200, -.5, 1.5},
-      "P.cal.etottracknorm");
-  auto hPcalEPPID = dSHMSEl.Histo1D(
-      {"P.cal.etottracknorm.PID", "Cuts: Tracking+PID;SHMS E/P;counts", 200, -.5, 1.5},
-      "P.cal.etottracknorm");
-  auto hPcalEPAll = dCOINElInTime.Histo1D(
-      {"P.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;SHMS E/P;counts", 200, -.5, 1.5},
-      "P.cal.etottracknorm");
-  // P.ngcer.npeSum
-  auto hPcerNpheNoCuts = d_coin.Histo1D(
-      {"P.ngcer.npeSum.NoCuts", "No Cuts;SHMS NGC #phe;counts", 200, -5, 76}, "P.ngcer.npeSum");
-  auto hPcerNpheTracking = dSHMSGoodTrack.Histo1D(
-      {"P.ngcer.npeSum.Tracking", "Cuts: Tracking;SHMS NGC #phe;counts", 200, -5, 76},
-      "P.ngcer.npeSum");
-  auto hPcerNphePID = dSHMSEl.Histo1D(
-      {"P.ngcer.npeSum.PID", "Cuts: Tracking+PID;SHMS NGC #phe;counts", 200, -5, 76},
-      "P.ngcer.npeSum");
-  auto hPcerNpheAll = dCOINElInTime.Histo1D(
-      {"P.ngcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS NGC #phe;counts", 200, -5, 76},
-      "P.ngcer.npeSum");
-  // P.hgcer.npeSum
-  auto hPhgcerNpheNoCuts = d_coin.Histo1D(
-      {"P.hgcer.npeSum.NoCuts", "No Cuts;SHMS HGC #phe;counts", 200, -5, 76}, "P.hgcer.npeSum");
-  auto hPhgcerNpheTracking = dSHMSGoodTrack.Histo1D(
-      {"P.hgcer.npeSum.Tracking", "Cuts: Tracking;SHMS HGC #phe;counts", 200, -5, 76},
-      "P.hgcer.npeSum");
-  auto hPhgcerNphePID = dSHMSEl.Histo1D(
-      {"P.hgcer.npeSum.PID", "Cuts: Tracking+PID;SHMS HGC #phe;counts", 200, -5, 76},
-      "P.hgcer.npeSum");
-  auto hPhgcerNpheAll = dCOINElInTime.Histo1D(
-      {"P.hgcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS HGC #phe;counts", 200, -5, 76},
-      "P.hgcer.npeSum");
-  // H.cal.etottracknorm
-  auto hHcalEPNoCuts =
-      d_coin.Histo1D({"H.cal.etottracknorm.NoCuts", "No Cuts;HMS E/P;counts", 200, -.5, 1.5},
-                     "H.cal.etottracknorm");
-  auto hHcalEPTracking = dHMSGoodTrack.Histo1D(
-      {"H.cal.etottracknorm.Tracking", "Cuts: Tracking;HMS E/P;counts", 200, -.5, 1.5},
-      "H.cal.etottracknorm");
-  auto hHcalEPPID = dHMSEl.Histo1D(
-      {"H.cal.etottracknorm.PID", "Cuts: Tracking+PID;HMS E/P;counts", 200, -.5, 1.5},
-      "H.cal.etottracknorm");
-  auto hHcalEPAll = dCOINElInTime.Histo1D(
-      {"H.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;HMS E/P;counts", 200, -.5, 1.5},
-      "H.cal.etottracknorm");
-  // H.cer.npeSum
-  auto hHcerNpheNoCuts = d_coin.Histo1D(
-      {"H.cer.npeSum.NoCuts", "No Cuts;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
-  auto hHcerNpheTracking = dHMSGoodTrack.Histo1D(
-      {"H.cer.npeSum.Tracking", "Cuts: Tracking;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
-  auto hHcerNphePID = dHMSEl.Histo1D(
-      {"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
-      "H.cer.npeSum");
-  auto hHcerNpheAll = dCOINElInTime.Histo1D(
-      {"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
-      "H.cer.npeSum");
-  // H.gtr.dp
-  auto hHdpNoCuts =
-      d_coin.Histo1D({"H.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
-  auto hHdpTracking = dHMSGoodTrack.Histo1D(
-      {"H.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
-  auto hHdpPID = dHMSEl.Histo1D(
-      {"H.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
-  auto hHdpTiming = dHMSElInTime.Histo1D(
-      {"H.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
-      "H.gtr.dp");
-  // H.gtr.th
-  auto hHthNoCuts = d_coin.Histo1D(
-      {"H.gtr.th.NoCuts", "No Cuts;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
-  auto hHthTracking = dHMSGoodTrack.Histo1D(
-      {"H.gtr.th.Tracking", "Cuts: Tracking;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
-  auto hHthPID = dHMSEl.Histo1D(
-      {"H.gtr.th.PID", "Cuts: Tracking+PID;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
-  auto hHthTiming = dHMSElInTime.Histo1D(
-      {"H.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{HMS};counts", 200, -0.1, 0.1},
-      "H.gtr.th");
-  // H.gtr.ph
-  auto hHphNoCuts =
-      d_coin.Histo1D({"H.gtr.ph.NoCuts", "No Cuts;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
-  auto hHphTracking = dHMSGoodTrack.Histo1D(
-      {"H.gtr.ph.Tracking", "Cuts: Tracking;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
-  auto hHphPID = dHMSEl.Histo1D(
-      {"H.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
-  auto hHphTiming = dHMSElInTime.Histo1D(
-      {"H.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{HMS};counts", 200, -0.1, 0.1},
-      "H.gtr.ph");
-  // H.gtr.y
-  auto hHyNoCuts =
-      d_coin.Histo1D({"H.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "H.gtr.y");
-  auto hHyTracking = dHMSGoodTrack.Histo1D(
-      {"H.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "H.gtr.y");
-  auto hHyPID =
-      dHMSEl.Histo1D({"H.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "H.gtr.y");
-  auto hHyTiming = dHMSElInTime.Histo1D(
-      {"H.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "H.gtr.y");
-
-  // scalers
-  //auto total_charge        = d_sh.Max("P.BCM4B.scalerChargeCut");
-  //auto shms_el_real_scaler = d_sh.Max("P.pEL_REAL.scaler");
-  //auto hms_el_real_scaler  = d_sh.Max("P.hEL_REAL.scaler");
-  //auto time_1MHz           = d_sh.Max("P.1MHz.scalerTime");
-  //auto time_1MHz_cut       = d_sh.Max("P.1MHz.scalerTimeCut");
-
-  auto yield_all = d.Count();
-  // 5 timing cut widths worth of random backgrounds
-  auto yield_coin          = d_coin.Count();
-  auto yield_HMSGoodTrack  = dHMSGoodTrack.Count();
-  auto yield_SHMSGoodTrack = dSHMSGoodTrack.Count();
-  auto yield_COINGoodTrack = dCOINGoodTrack.Count();
-  auto yield_HMSEl         = dHMSEl.Count();
-  auto yield_SHMSEl        = dSHMSEl.Count();
-  auto yield_COINEl        = dCOINEl.Count();
-  auto yield_HMSElInTime   = dHMSElInTime.Count();
-  auto yield_HMSElRandom   = dHMSElRandom.Count();
-  auto yield_SHMSElInTime  = dSHMSElInTime.Count();
-  auto yield_SHMSElRandom  = dSHMSElRandom.Count();
-  auto yield_COINElInTime  = dCOINElInTime.Count();
-  auto yield_COINElRandom  = dCOINElRandom.Count();
-  auto yield_coin_raw      = dCOINElInTime.Count();
-  auto yield_coin_random   = dCOINElRandom.Count();
-
-
-  // -------------------------------------
-  // End lazy eval
-  // -------------------------------------
-  auto n_coin_good  = *yield_coin_raw - *yield_coin_random / 5.;
-  auto n_HMSElGood  = *yield_HMSElInTime - *yield_HMSElRandom / 5;
-  auto n_SHMSElGood = *yield_SHMSElInTime - *yield_SHMSElRandom / 5;
-  auto n_COINElGood = *yield_COINElInTime - *yield_COINElRandom / 5;
-
-  hPdpNoCuts->DrawCopy();
-  //std::cout << "  coin COUNTS : " << *(d_coin.Count()) << "\n";
-  //std::cout << " yield_HMSEl : " << *(yield_HMSEl) << "\n";
-  std::cout << " yield_COINEl : " << *(yield_COINEl) << "\n";
-  //std::cout << "  ALL COUNTS : " << *yield_all << "\n";
-  //std::cout << " GOOD COUNTS : " << n_COINElGood << "\n";
-  //
-  if( 4 != (*(yield_COINEl)) ){
-    std::exit(-1);
-  }
-  std::exit(0);
-}

tests/elastic_test.sh

-#!/bin/bash
-
-#echo "This is the elastic testing..."
-#echo " "
-#echo "There are currently 0 tests to run!"
-#which hcana
-#
-#ls -lrth
-#ls -lrth build
-#
-#git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/hallc_replay_csv.git
-#git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/online_csv.git
-#cd online_csv 
-#ln -s ../hallc_reaply_csv/PARAM
-## and the reset
-#
-#mkdir raw 
-#pushd raw
-#  wget coin.dat
-#popd
-
-
-singularity help build/Singularity.hcana.simg
-
-
-singularity exec build/Singularity.hcana.simg hcana tests/elastic_test.cxx

tests/elastic_test2.sh

-#!/bin/bash
-
-echo "This is the elastic testing..."
-echo " "
-echo "There are currently 0 tests to run!"
-which hcana
-
-ls -lrth
-ls -lrth build
-
-git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/hallc_replay_csv.git
-git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/online_csv.git
-
-cd online_csv 
-ln -s ../hallc_reaply_csv/PARAM
-# and the reset
-
-mkdir raw 
-pushd raw
-  wget coin.dat
-popd
-
-
-singularity help build/Singularity.hcana.simg
-
-singularity exec build/Singularity.hcana.simg which hcana
-
-singularity exec build/Singularity.hcana.simg hcana tests/my_root_script.cxx
-
-echo " WOOOO"

tests/my_root_script.cxx

-void my_root_script() {
-
-
-
-  std::cout << "Hello from my_root_script.cxx!\n";
-
-  std::cout << "This should be run with singularity\n";
-  double pi = 3.14;
-
-  auto pi_equals_3 = (3 == pi);
-  std::cout <<  " pi_equals_3 = " << pi_equals_3 << "\n";
-
-  if( pi_equals_3) {
-    std::cout << "what the hell?\n";
-    std::exit( 0 );
-  }
-  /* else */
-  
-  std::exit( -1 );
-}

tests/replay_elastic_data.sh

-#!/bin/bash
-
-echo "This is the elastic testing..."
-echo " "
-echo "There are currently 0 tests to run!"
-which hcana
-
-ls -lrth
-ls -lrth build
-mkdir -p ROOTfiles
-git clone https://eicweb.phy.anl.gov/whit/ci_test_data.git 
-git clone https://eicweb.phy.anl.gov/jlab/hallc/exp/CSV/hallc_replay_csv.git
-git clone https://eicweb.phy.anl.gov/jlab/hallc/exp/CSV/online_csv.git
-
-cd online_csv 
-mkdir -p logs
-ln -s ../hallc_replay_csv/PARAM
-ln -s ../hallc_replay_csv/DBASE
-ln -s ../hallc_replay_csv/CALIBRATION
-ln -s ../hallc_replay_csv/DEF-files
-ln -s ../hallc_replay_csv/MAPS
-ln -s ../hallc_replay_csv/SCRIPTS
-ln -s ../hallc_replay_csv/DATFILES
-ln -s ../ci_test_data/raw
-ln -s ../ROOTfiles
-# and the reset
-
-ls -lrth
-ls -lrth raw/
-ls -lrth ROOTfiles/
-pwd
-# run replay script
-df -h
-
-singularity exec ../build/Singularity.hcana.simg hcana -b -q "../tests/elastic_coin_replay.cxx+(6012,50000)"
-
-echo "hcana calls... the coin replay script and outputs blah.root"
-