From 0f80f3adbf6161ee8a4aacefa1e3689862163305 Mon Sep 17 00:00:00 2001
From: Simon Zhamkochyan <szh@jlab.org>
Date: Sat, 15 Sep 2012 20:01:00 -0400
Subject: [PATCH] Added pedestal and threshold calculations for Shower counter
---
examples/DBASE/test.database | 2 +-
src/THcShower.cxx | 32 ++++++++--
src/THcShower.h | 4 +-
src/THcShowerPlane.cxx | 117 +++++++++++++++++++++++++++++++++--
src/THcShowerPlane.h | 26 +++++++-
5 files changed, 167 insertions(+), 14 deletions(-)
diff --git a/examples/DBASE/test.database b/examples/DBASE/test.database
index 66bff4a..59790ba 100644
--- a/examples/DBASE/test.database
+++ b/examples/DBASE/test.database
@@ -1,7 +1,7 @@
# ENGINE style parameter vs. run number database
50017
g_ctp_parm_filename="PARAM/general.param"
-g_decode_map_filename="MAPS/july04.map"
+g_decode_map_filename="MAPS/jan03.map"
47000-48000
g_ctp_parm_filename="PARAM/general.param"
g_decode_map_filename="MAPS/jan03.map"
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index a551db5..96d94db 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -14,6 +14,7 @@
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
+#include "THaCutList.h"
#include "THcParmList.h"
#include "VarDef.h"
#include "VarType.h"
@@ -298,19 +299,40 @@ void THcShower::DeleteArrays()
//_____________________________________________________________________________
inline
-void THcShower::ClearEvent()
+void THcShower::Clear(Option_t* opt)
{
- // Reset per-event data.
-
- fTrackProj->Clear();
+// Reset per-event data.
+ for(Int_t ip=0;ip<fNLayers;ip++) {
+ fPlanes[ip]->Clear(opt);
+ }
+ // fTrackProj->Clear();
}
//_____________________________________________________________________________
Int_t THcShower::Decode( const THaEvData& evdata )
{
-
// Get the Hall C style hitlist (fRawHitList) for this event
Int_t nhits = THcHitList::DecodeToHitList(evdata);
+
+if(gHaCuts->Result("Pedestal_event")) {
+ Int_t nexthit = 0;
+ for(Int_t ip=0;ip<fNLayers;ip++) {
+ nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
+//cout << "nexthit = " << nexthit << endl;
+ }
+ fAnalyzePedestals = 1; // Analyze pedestals first normal events
+ return(0);
+ }
+
+ if(fAnalyzePedestals) {
+ for(Int_t ip=0;ip<fNLayers;ip++) {
+ fPlanes[ip]->CalculatePedestals();
+ }
+ fAnalyzePedestals = 0; // Don't analyze pedestals next event
+ }
+
+
+
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
diff --git a/src/THcShower.h b/src/THcShower.h
index aa58aaa..9926c63 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -21,7 +21,7 @@ public:
THcShower( const char* name, const char* description = "",
THaApparatus* a = NULL );
virtual ~THcShower();
-
+ virtual void Clear( Option_t* opt="" );
virtual Int_t Decode( const THaEvData& );
virtual EStatus Init( const TDatime& run_time );
virtual Int_t CoarseProcess( TClonesArray& tracks );
@@ -38,7 +38,7 @@ public:
THcShower(); // for ROOT I/O
protected:
-
+ Int_t fAnalyzePedestals;
// Calibration
// Per-event data
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index 9389f08..9ec3dcc 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -14,6 +14,7 @@
#include "THcHitList.h"
#include "THcShower.h"
#include "TClass.h"
+#include "math.h"
#include <cstring>
#include <cstdio>
@@ -35,9 +36,10 @@ THcShowerPlane::THcShowerPlane( const char* name,
// Normal constructor with name and description
fPosADCHits = new TClonesArray("THcSignalHit",13);
fNegADCHits = new TClonesArray("THcSignalHit",13);
+#if ROOT_VERSION_CODE < ROOT_VERSION(5,32,0)
fPosADCHitsClass = fPosADCHits->GetClass();
fNegADCHitsClass = fNegADCHits->GetClass();
-
+#endif
fPosADC1 = new TClonesArray("THcSignalHit",13);
fPosADCHitsClass = fPosADC1->GetClass();
@@ -94,7 +96,8 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
strcat(parname,"_nr");
cout << " Getting value of SHOWER!!!" << parname << endl;
-// fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer();
+fNelem = 13;
+ // fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer();
//
// parname[plen]='\0';
// strcat(parname,"_spacing");
@@ -109,7 +112,7 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
// Create arrays to hold results here
-
+ InitializePedestals();
return kOK;
}
//_____________________________________________________________________________
@@ -264,8 +267,12 @@ CalADC1File = fopen("adc1_new.dat", "a");
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = nexthit;
+ //cout << "nrawhits = " << nrawhits << endl;
+ //cout << "nexthit = " << nexthit << endl;
while(ihit < nrawhits) {
THcShowerHit* hit = (THcShowerHit *) rawhits->At(ihit);
+
+//cout << "fplane = " << hit->fPlane << " Num = " << fLayerNum << endl;
if(hit->fPlane > fLayerNum) {
break;
}
@@ -286,7 +293,8 @@ THcSignalHit *sighit1 = (THcSignalHit*) fPosADC1->ConstructedAt(nPosADCHits++);
//fprintf(CalADC1File, "%d\n", hit->fADC_pos);
}
-double thresh_pos = hcal_new_threshold_pos[hit->fCounter + 13*(hit->fPlane -1) -1];
+double thresh_pos = fPosThresh[hit->fCounter -1];
+//double thresh_pos = hcal_new_threshold_pos[hit->fCounter + 13*(hit->fPlane -1) -1];
if(hit->fADC_pos > thresh_pos) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0)
THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
@@ -302,7 +310,8 @@ if(!obj->TestBit (TObject::kNotDeleted))
sighit->Set(hit->fCounter, hit->fADC_pos);
}
-double thresh_neg = hcal_new_threshold_neg[hit->fCounter + 13*(hit->fPlane -1) -1];
+double thresh_neg = fNegThresh[hit->fCounter -1];
+//double thresh_neg = hcal_new_threshold_neg[hit->fCounter + 13*(hit->fPlane -1) -1];
if(hit->fADC_neg > thresh_neg) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0)
THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
@@ -322,8 +331,106 @@ if(!obj->TestBit (TObject::kNotDeleted))
fclose(CalADC1File);
return(ihit);
}
+//_____________________________________________________________________________
+Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
+{
+ // Extract the data for this plane from hit list, accumulating into
+ // arrays for calculating pedestals.
+
+ Int_t nrawhits = rawhits->GetLast()+1;
+ // cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " " << nexthit << "/" << nrawhits << endl;
+ Int_t ihit = nexthit;
+ while(ihit < nrawhits) {
+ THcShowerHit* hit = (THcShowerHit *) rawhits->At(ihit);
+//cout << "fPlane = " << hit->fPlane << " Limit = " << fPlaneNum << endl;
+ if(hit->fPlane > fLayerNum) {
+ break;
+ }
+ Int_t element = hit->fCounter - 1; // Should check if in range
+ Int_t adcpos = hit->fADC_pos;
+ Int_t adcneg = hit->fADC_neg;
+
+ if(adcpos <= fPosPedLimit[element]) {
+ fPosPedSum[element] += adcpos;
+ fPosPedSum2[element] += adcpos*adcpos;
+ fPosPedCount[element]++;
+ if(fPosPedCount[element] == fMinPeds/5) {
+ fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
+ }
+ }
+ if(adcneg <= fNegPedLimit[element]) {
+ fNegPedSum[element] += adcneg;
+ fNegPedSum2[element] += adcneg*adcneg;
+ fNegPedCount[element]++;
+ if(fNegPedCount[element] == fMinPeds/5) {
+ fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
+ }
+ }
+ ihit++;
+ }
+
+ fNPedestalEvents++;
+
+ return(ihit);
+}
+ //_____________________________________________________________________________
+void THcShowerPlane::CalculatePedestals( )
+{
+ // Use the accumulated pedestal data to calculate pedestals
+ // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
+
+ for(Int_t i=0; i<fNelem;i++) {
+
+ // Positive tubes
+ fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
+ fPosSig[i] = sqrt((fPosPedSum2[i] - 2.*fPosPed[i]*fPosPedSum[i])/TMath::Max(1, fPosPedCount[i]) + fPosPed[i]*fPosPed[i]);
+ // fPosThresh[i] = fPosPed[i] + 15;
+ fPosThresh[i] = fPosPed[i] + TMath::Min(50., TMath::Max(10., 3.*fPosSig[i]));
+
+ // Negative tubes
+ fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
+ fNegSig[i] = sqrt((fNegPedSum2[i] - 2.*fNegPed[i]*fNegPedSum[i])/TMath::Max(1, fNegPedCount[i]) + fNegPed[i]*fNegPed[i]);
+ // fNegThresh[i] = fNegPed[i] + 15;
+ fNegThresh[i] = fNegPed[i] + TMath::Min(50., TMath::Max(10., 3.*fNegSig[i]));
+
+ // cout << i+1 << " " << 3.*fPosSig[i] << " " << 3.*fNegSig[i] << endl;
+ }
+ // cout << " " << endl;
+}
+
+//_____________________________________________________________________________
+void THcShowerPlane::InitializePedestals( )
+{
+ fNPedestalEvents = 0;
+ fMinPeds = 500; // In engine, this is set in parameter file
+ fPosPedSum = new Int_t [fNelem];
+ fPosPedSum2 = new Int_t [fNelem];
+ fPosPedLimit = new Int_t [fNelem];
+ fPosPedCount = new Int_t [fNelem];
+ fNegPedSum = new Int_t [fNelem];
+ fNegPedSum2 = new Int_t [fNelem];
+ fNegPedLimit = new Int_t [fNelem];
+ fNegPedCount = new Int_t [fNelem];
+
+ fPosSig = new Double_t [fNelem];
+ fNegSig = new Double_t [fNelem];
+ fPosPed = new Double_t [fNelem];
+ fNegPed = new Double_t [fNelem];
+ fPosThresh = new Double_t [fNelem];
+ fNegThresh = new Double_t [fNelem];
+ for(Int_t i=0;i<fNelem;i++) {
+ fPosPedSum[i] = 0;
+ fPosPedSum2[i] = 0;
+ fPosPedLimit[i] = 1000; // In engine, this are set in parameter file
+ fPosPedCount[i] = 0;
+ fNegPedSum[i] = 0;
+ fNegPedSum2[i] = 0;
+ fNegPedLimit[i] = 1000; // In engine, this are set in parameter file
+ fNegPedCount[i] = 0;
+ }
+}
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 4e44d1c..4d8c4a8 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -39,6 +39,8 @@ class THcShowerPlane : public THaSubDetector {
virtual Bool_t IsPid() { return kFALSE; }
virtual Int_t ProcessHits(TClonesArray* rawhits, Int_t nexthit);
+ virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit);
+ virtual void CalculatePedestals( );
Double_t fSpacing;
@@ -65,9 +67,31 @@ TClonesArray* fPosADC[13];
Int_t fLayerNum;
+Int_t fPlaneNum; /* Which plane am I 1-4 */
+ Int_t fNelem; /* Need since we don't inherit from
+ detector base class */
+ Int_t fNPedestalEvents; /* Number of pedestal events */
+ Int_t fMinPeds; /* Only analyze/update if num events > */
+ Int_t *fPosPedSum; /* Accumulators for pedestals */
+ Int_t *fPosPedSum2;
+ Int_t *fPosPedLimit;
+ Int_t *fPosPedCount;
+ Int_t *fNegPedSum;
+ Int_t *fNegPedSum2;
+ Int_t *fNegPedLimit;
+ Int_t *fNegPedCount;
+
+ Double_t *fPosPed;
+ Double_t *fPosSig;
+ Double_t *fPosThresh;
+ Double_t *fNegPed;
+ Double_t *fNegSig;
+ Double_t *fNegThresh;
+
+
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
-
+ virtual void InitializePedestals( );
ClassDef(THcShowerPlane,0)
};
#endif
--
GitLab