From e5b40a347abdab4c2d41e45e79b843e769b07b39 Mon Sep 17 00:00:00 2001
From: Mark Jones <jones@jlab.org>
Date: Wed, 4 Mar 2020 09:39:28 -0500
Subject: [PATCH] Updated THcScalerEvtHandler

Dave Mack introduced two new variables to the
calculation of beam current to described correction for
the nonlinearity above 60uA in BCM1 and BCM2
seen in the Fall 18, Spring 19 and Summer 19running.
See for example summer 2019
hallcweb.jlab.org/doc-private/ShowDocument?docid=1037
Fall 18/Spring 19
hallcweb.jlab.org/doc-private/ShowDocument?docid=1042

For a BCM scaler read the difference, Dcounts, from the
previous scaler and get the time difference, Dtime.
Original calculation is:
I = (Dcounts/Dtime-fBCMOffset)/BCM_Gain
Now add an non-linearity correction
I = I + fBCM_SatQuadratic*Max(I-fBCM_SatOffset,0.0)^2

Add parameters fBCM_SatOffset and fBCM_SatQuadratic
that can be optional read-in. If parameters are
not read-in the parameters are set to zero, so
there is no correction.

Add the non-linearity correction to all BCM current
calculations.
---
 src/THcScalerEvtHandler.cxx | 43 +++++++++++++++++++++++++++++++------
 src/THcScalerEvtHandler.h   |  2 ++
 2 files changed, 38 insertions(+), 7 deletions(-)

diff --git a/src/THcScalerEvtHandler.cxx b/src/THcScalerEvtHandler.cxx
index 1fa098a..d4cf849 100644
--- a/src/THcScalerEvtHandler.cxx
+++ b/src/THcScalerEvtHandler.cxx
@@ -76,7 +76,7 @@ static const UInt_t defaultDT = 4;
 
 THcScalerEvtHandler::THcScalerEvtHandler(const char *name, const char* description)
   : THaEvtTypeHandler(name,description),
-    fBCM_Gain(0), fBCM_Offset(0), fBCM_delta_charge(0),
+    fBCM_Gain(0), fBCM_Offset(0), fBCM_SatOffset(0), fBCM_SatQuadratic(0), fBCM_delta_charge(0),
     evcount(0), evcountR(0.0), ifound(0), fNormIdx(-1),
     fNormSlot(-1),
     dvars(0),dvars_prev_read(0), dvarsFirst(0), fScalerTree(0), fUseFirstEvent(kTRUE),
@@ -104,6 +104,8 @@ THcScalerEvtHandler::~THcScalerEvtHandler()
   delete [] dvarsFirst;
   delete [] fBCM_Gain;
   delete [] fBCM_Offset;
+  delete [] fBCM_SatOffset;
+  delete [] fBCM_SatQuadratic;
   delete [] fBCM_delta_charge;
 
   for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
@@ -152,11 +154,15 @@ Int_t THcScalerEvtHandler::ReadDatabase(const TDatime& date )
   if(fNumBCMs > 0) {
     fBCM_Gain = new Double_t[fNumBCMs];
     fBCM_Offset = new Double_t[fNumBCMs];
+    fBCM_SatOffset = new Double_t[fNumBCMs];
+    fBCM_SatQuadratic = new Double_t[fNumBCMs];
     fBCM_delta_charge= new Double_t[fNumBCMs];
     string bcm_namelist;
     DBRequest list2[]={
       {"BCM_Gain",      fBCM_Gain,         kDouble, (UInt_t) fNumBCMs},
       {"BCM_Offset",     fBCM_Offset,       kDouble,(UInt_t) fNumBCMs},
+      {"BCM_SatQuadratic",     fBCM_SatQuadratic,       kDouble,(UInt_t) fNumBCMs,1},
+      {"BCM_SatOffset",     fBCM_SatOffset,       kDouble,(UInt_t) fNumBCMs,1},
       {"BCM_Names",     &bcm_namelist,       kString},
       {"BCM_Current_threshold",     &fbcm_Current_Threshold,       kDouble,0, 1},
       {"BCM_Current_threshold_index",     &fbcm_Current_Threshold_Index,       kInt,0,1},
@@ -164,6 +170,10 @@ Int_t THcScalerEvtHandler::ReadDatabase(const TDatime& date )
     };
     fbcm_Current_Threshold = 0.0;
     fbcm_Current_Threshold_Index = 0;
+    for(Int_t i=0;i<fNumBCMs;i++) {
+      fBCM_SatOffset[i]=0.;
+      fBCM_SatQuadratic[i]=0.;
+    }
     gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
     vector<string> bcm_names = vsplit(bcm_namelist);
     for(Int_t i=0;i<fNumBCMs;i++) {
@@ -412,12 +422,17 @@ Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
 	      }
 	    if (scalerloc[ivar]->ikind == ICURRENT) {
               dvars[ivar]=0.;
-		if (bcm_ind != -1) dvars[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+	      if (bcm_ind != -1) {
+                 dvars[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		 dvars[ivar]=dvars[ivar]+fBCM_SatOffset[bcm_ind]*TMath::Max(dvars[ivar]-fBCM_SatOffset[i],0.0);
+	      }
          	if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvars[ivar];
 	    }
 	    if (scalerloc[ivar]->ikind == ICHARGE) {
 	      if (bcm_ind != -1) {
-		fBCM_delta_charge[bcm_ind]=fDeltaTime*((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		Double_t cur_temp=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		cur_temp=cur_temp+fBCM_SatQuadratic[bcm_ind]*TMath::Power(TMath::Max(cur_temp-fBCM_SatOffset[bcm_ind],0.0),2.0);
+		fBCM_delta_charge[bcm_ind]=fDeltaTime*cur_temp;
 		dvars[ivar]+=fBCM_delta_charge[bcm_ind];
 	      }
 	    }
@@ -452,12 +467,17 @@ Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
 		}
 	    if (scalerloc[ivar]->ikind == ICURRENT) {
 	        dvarsFirst[ivar]=0.0;
-                if (bcm_ind != -1) dvarsFirst[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+                if (bcm_ind != -1) {
+                 dvarsFirst[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		 dvarsFirst[ivar]=dvarsFirst[ivar]+fBCM_SatQuadratic[bcm_ind]*TMath::Power(TMath::Max(dvars[ivar]-fBCM_SatOffset[i],0.0),2.);
+		}
          	if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvarsFirst[ivar];
 	    }
 	    if (scalerloc[ivar]->ikind == ICHARGE) {
 	      if (bcm_ind != -1) {
-		fBCM_delta_charge[bcm_ind]=fDeltaTime*((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		Double_t cur_temp=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		cur_temp=cur_temp+fBCM_SatQuadratic[bcm_ind]*TMath::Power(TMath::Max(cur_temp-fBCM_SatOffset[bcm_ind],0.0),2.);
+		fBCM_delta_charge[bcm_ind]=fDeltaTime*cur_temp;
                dvarsFirst[ivar]+=fBCM_delta_charge[bcm_ind];
 	      }
 	    }
@@ -519,7 +539,12 @@ Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
 		  diff = scaldata - scal_prev_read[nscal-1];
 		}
 		dvars[ivar]=0.;
- 		if (fDeltaTime>0) dvars[ivar]=(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+ 		if (fDeltaTime>0) {
+		Double_t cur_temp=(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		cur_temp=cur_temp+fBCM_SatQuadratic[bcm_ind]*TMath::Power(TMath::Max(cur_temp-fBCM_SatOffset[bcm_ind],0.0),2.);
+		  
+		dvars[ivar]=cur_temp;
+		}
 	      }
 	      if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvars[ivar];
 	    }
@@ -533,7 +558,11 @@ Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
 		  diff = scaldata - scal_prev_read[nscal-1];
 		}
 		fBCM_delta_charge[bcm_ind]=0;
-		if (fDeltaTime>0)  fBCM_delta_charge[bcm_ind]=fDeltaTime*(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		if (fDeltaTime>0)  {
+		  Double_t cur_temp=(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
+		  cur_temp=cur_temp+fBCM_SatQuadratic[bcm_ind]*TMath::Power(TMath::Max(cur_temp-fBCM_SatOffset[bcm_ind],0.0),2.);
+		fBCM_delta_charge[bcm_ind]=fDeltaTime*cur_temp;
+		}
                  dvars[ivar]+=fBCM_delta_charge[bcm_ind];
 	      }
 	    }
diff --git a/src/THcScalerEvtHandler.h b/src/THcScalerEvtHandler.h
index c4749ff..1e99e4a 100644
--- a/src/THcScalerEvtHandler.h
+++ b/src/THcScalerEvtHandler.h
@@ -59,6 +59,8 @@ private:
    Int_t fNumBCMs;
    Double_t *fBCM_Gain;
    Double_t *fBCM_Offset;
+   Double_t *fBCM_SatOffset;
+   Double_t *fBCM_SatQuadratic;
    Double_t *fBCM_delta_charge;
    Double_t fTotalTime;
    Double_t fDeltaTime;
-- 
GitLab