diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 15172b1d71a6692245305a5028dc8fc0f7f5c225..e1dbf46bdbec3bd8569e010e7ea950008f75bea1 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -1169,6 +1169,13 @@ Int_t THcShower::FineProcess( TClonesArray& tracks )
}
} //over tracks
+ if (Ntracks>0) {
+ for(UInt_t ip=0;ip<fNLayers;ip++) {
+ fPlanes[ip]-> AccumulateStat(tracks);
+ }
+ if(fHasArray) fArray->AccumulateStat(tracks);
+ }
+
//Debug output.
if (fdbg_tracks_cal) {
diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
index e64a577bae1ecb2424de56fc8bc7399bb6816c73..007f5df9b5f154b5d07d120dbc9f7ccfcd685cb4 100644
--- a/src/THcShowerArray.cxx
+++ b/src/THcShowerArray.cxx
@@ -17,6 +17,8 @@
#include "math.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
+#include "THcCherenkov.h" //for efficiency calculations
+#include "THcHallCSpectrometer.h"
#include <cstring>
#include <cstdio>
@@ -129,6 +131,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
fPedSampHigh=9;
fDataSampLow=23;
fDataSampHigh=49;
+ fStatCerMin=1.;
+ fStatSlop=3.;
+ fStatMaxChi2=10.;
DBRequest list[]={
{"cal_arr_nrows", &fNRows, kInt},
{"cal_arr_ncolumns", &fNColumns, kInt},
@@ -148,6 +153,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
{"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
{"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1},
{"cal_debug_adc", &fDebugAdc, kInt, 0, 1},
+ {"stat_cermin", &fStatCerMin, kDouble, 0, 1},
+ {"stat_slop_array", &fStatSlop, kDouble, 0, 1},
+ {"stat_maxchisq", &fStatMaxChi2, kDouble, 0, 1},
{0}
};
@@ -342,6 +350,13 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
//fE = new Double_t[fNelem];
+ // Numbers of tracks and hits , for efficiency calculations.
+
+ fStatNumTrk = vector<Int_t> (fNelem, 0);
+ fStatNumHit = vector<Int_t> (fNelem, 0);
+ fTotStatNumTrk = 0;
+ fTotStatNumHit = 0;
+
#ifdef HITPIC
hitpic = new char*[fNRows];
for(Int_t row=0;row<fNRows;row++) {
@@ -370,6 +385,14 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
+ // Register counters for efficiency calculations in gHcParms so that the
+ // variables can be used in end of run reports.
+
+ gHcParms->Define(Form("%sstat_trksum_array", GetParent()->GetPrefix()),
+ "Number of tracks in calo. array", fTotStatNumTrk);
+ gHcParms->Define(Form("%sstat_hitsum_array", GetParent()->GetPrefix()),
+ "Number of hits in calo. array", fTotStatNumHit);
+
// Register variables in global list
if (fDebugAdc) {
RVarDef vars[] = {
@@ -1163,3 +1186,99 @@ Double_t THcShowerArray::clMaxEnergyBlock(THcShowerCluster* cluster) {
}
return max_block;
}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::AccumulateStat(TClonesArray& tracks )
+{
+ // Accumumate statistics for efficiency calculations.
+ //
+ // Choose electron events in gas Cherenkov with good Chisq of the best track.
+ // Project best track to Array,
+ // calculate module number for the track,
+ // accrue number of tracks for the module,
+ // accrue number of hits for the module, if it is hit.
+ // Accrue total numbers of tracks and hits for Array.
+
+ THaTrack* BestTrack = static_cast<THaTrack*>( tracks[0]);
+ if (BestTrack->GetChi2()/BestTrack->GetNDoF() > fStatMaxChi2) return 0;
+
+ THcHallCSpectrometer *app=dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+
+ THaDetector* detc;
+ if (GetParent()->GetPrefix()[0] == 'P')
+ detc = app->GetDetector("hgcer");
+ else
+ detc = app->GetDetector("cer");
+
+ THcCherenkov* hgcer = dynamic_cast<THcCherenkov*>(detc);
+ if (!hgcer) {
+ cout << "****** THcShowerArray::AccumulateStat: HGCer not found! ******"
+ << endl;
+ return 0;
+ }
+
+ if (hgcer->GetCerNPE() < fStatCerMin) return 0;
+
+ Double_t XTrk = kBig;
+ Double_t YTrk = kBig;
+ Double_t pathl = kBig;
+
+ // Track interception with Array. The coordinates are in the calorimeter's
+ // local system.
+
+ fOrigin = GetOrigin();
+ THcShower* fParent = (THcShower*) GetParent();
+ fParent->CalcTrackIntercept(BestTrack, pathl, XTrk, YTrk);
+
+ // Transform coordiantes to the spectrometer's coordinate system.
+ XTrk += GetOrigin().X();
+ YTrk += GetOrigin().Y();
+
+ for (Int_t i=0; i<fNelem; i++) {
+
+ Int_t row = i%fNRows;
+ Int_t col =i/fNRows;
+
+ if (TMath::Abs(XTrk - fXPos[row][col]) < fStatSlop &&
+ TMath::Abs(YTrk - fYPos[row][col]) < fStatSlop) {
+
+ fStatNumTrk.at(i)++;
+ fTotStatNumTrk++;
+
+ if (fGoodAdcPulseInt.at(i) > 0.) {
+ fStatNumHit.at(i)++;
+ fTotStatNumHit++;
+ }
+
+ }
+
+ }
+
+ if ( ((THcShower*) GetParent())->fdbg_tracks_cal ) {
+ cout << "---------------------------------------------------------------\n";
+ cout << "THcShowerArray::AccumulateStat:" << endl;
+ cout << " Chi2/NDF = " <<BestTrack->GetChi2()/BestTrack->GetNDoF() <<endl;
+ cout << " HGCER Npe = " << hgcer->GetCerNPE() << endl;
+ cout << " XTrk, YTrk = " << XTrk << " " << YTrk << endl;
+ for (Int_t i=0; i<fNelem; i++) {
+
+ Int_t row = i%fNRows;
+ Int_t col =i/fNRows;
+
+ if (TMath::Abs(XTrk - fXPos[row][col]) < fStatSlop &&
+ TMath::Abs(YTrk - fYPos[row][col]) < fStatSlop) {
+
+ cout << " Module " << i << ", X=" << fXPos[i/fNRows][i%fNRows]
+ << ", Y=" << fYPos[i/fNRows][i%fNRows] << " matches track" << endl;
+
+ if (fGoodAdcPulseInt.at(i) > 0.)
+ cout << " PulseIntegral = " << fGoodAdcPulseInt.at(i) << endl;
+ }
+ }
+
+ cout << "---------------------------------------------------------------\n";
+ // getchar();
+ }
+
+ return 1;
+}
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
index 8e45937cb1e1e15c7c38c007cc4a7f3ea289638c..0ac6f1a87c440fdfc10ec01d0b45638fe4d5705a 100644
--- a/src/THcShowerArray.h
+++ b/src/THcShowerArray.h
@@ -89,6 +89,8 @@ public:
Double_t fvYmin();
Double_t clMaxEnergyBlock(THcShowerCluster* cluster);
+ Int_t AccumulateStat(TClonesArray& tracks);
+
protected:
#ifdef HITPIC
@@ -192,6 +194,15 @@ protected:
TClonesArray* frAdcPulseInt;
TClonesArray* frAdcPulseAmp;
+ //Quatitites for efficiency calculations.
+
+ Double_t fStatCerMin;
+ Double_t fStatSlop;
+ Double_t fStatMaxChi2;
+ vector<Int_t> fStatNumTrk;
+ vector<Int_t> fStatNumHit;
+ Int_t fTotStatNumTrk;
+ Int_t fTotStatNumHit;
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index 4b5d7cce17a9e5fbf9e7900c98d144fe9b30774f..4df9594a93574e97c203baac002c2e3de7dc068e 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -17,6 +17,8 @@ One plane of shower blocks with side readout
#include "math.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
+#include "THcCherenkov.h" //for efficiency calculations
+#include "THcHallCSpectrometer.h"
#include <cstring>
#include <cstdio>
@@ -158,6 +160,9 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
fDataSampHigh=49;
fAdcNegThreshold=0.;
fAdcPosThreshold=0.;
+ fStatCerMin=1.;
+ fStatSlop=2.;
+ fStatMaxChi2=10.;
DBRequest list[]={
{"cal_AdcNegThreshold", &fAdcNegThreshold, kDouble, 0, 1},
{"cal_AdcPosThreshold", &fAdcPosThreshold, kDouble, 0, 1},
@@ -166,6 +171,9 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
{"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
{"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1},
{"cal_debug_adc", &fDebugAdc, kInt, 0, 1},
+ {"stat_cermin", &fStatCerMin, kDouble, 0, 1},
+ {"stat_slop", &fStatSlop, kDouble, 0, 1},
+ {"stat_maxchisq", &fStatMaxChi2, kDouble, 0, 1},
{0}
};
@@ -243,6 +251,13 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
// fEneg = new Double_t[fNelem];
// fEmean= new Double_t[fNelem];
+ // Numbers of tracks and hits , for efficiency calculations.
+
+ fStatNumTrk = vector<Int_t> (fNelem, 0);
+ fStatNumHit = vector<Int_t> (fNelem, 0);
+ fTotStatNumTrk = 0;
+ fTotStatNumHit = 0;
+
// Debug output.
if (fParent->fdbg_init_cal) {
@@ -304,6 +319,19 @@ Int_t THcShowerPlane::DefineVariables( EMode mode )
DefineVarsFromList( vars, mode);
} //end debug statement
+ // Register counters for efficiency calculations in gHcParms so that the
+ // variables can be used in end of run reports.
+
+ gHcParms->Define(Form("%sstat_trksum%d", GetParent()->GetPrefix(), fLayerNum),
+ Form("Number of tracks in calo. layer %d",fLayerNum), fTotStatNumTrk);
+ gHcParms->Define(Form("%sstat_hitsum%d", GetParent()->GetPrefix(), fLayerNum),
+ Form("Number of hits in calo. layer %d", fLayerNum), fTotStatNumHit);
+
+ cout << "THcShowerPlane::DefineVariables: registered counters "
+ << Form("%sstat_trksum%d",GetParent()->GetPrefix(),fLayerNum) << " and "
+ << Form("%sstat_hitsum%d",GetParent()->GetPrefix(),fLayerNum) << endl;
+ // getchar();
+
RVarDef vars[] = {
{"posAdcErrorFlag", "List of positive raw ADC Error Flags", "frPosAdcErrorFlag.THcSignalHit.GetData()"},
{"negAdcErrorFlag", "List of negative raw ADC Error Flags ", "frNegAdcErrorFlag.THcSignalHit.GetData()"},
@@ -891,3 +919,92 @@ void THcShowerPlane::InitializePedestals( )
fNegPedCount[i] = 0;
}
}
+
+//_____________________________________________________________________________
+Int_t THcShowerPlane::AccumulateStat(TClonesArray& tracks )
+{
+ // Accumumate statistics for efficiency calculations.
+ //
+ // Choose electron events in gas Cherenkov with good Chisq of the best track.
+ // Project best track to the plane,
+ // calculate row number for the track,
+ // accrue number of tracks for the row,
+ // accrue number of hits for the row, if row is hit.
+ // Accrue total numbers of tracks and hits for plane.
+
+ THaTrack* BestTrack = static_cast<THaTrack*>( tracks[0]);
+ if (BestTrack->GetChi2()/BestTrack->GetNDoF() > fStatMaxChi2) return 0;
+
+ THcHallCSpectrometer *app=dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+
+ THaDetector* detc;
+ if (GetParent()->GetPrefix()[0] == 'P')
+ detc = app->GetDetector("hgcer");
+ else
+ detc = app->GetDetector("cer");
+
+ THcCherenkov* hgcer = dynamic_cast<THcCherenkov*>(detc);
+ if (!hgcer) {
+ cout << "****** THcShowerPlane::AccumulateStat: HGCer not found! ******"
+ << endl;
+ return 0;
+ }
+
+ if (hgcer->GetCerNPE() < fStatCerMin) return 0;
+
+ Double_t XTrk = kBig;
+ Double_t YTrk = kBig;
+ Double_t pathl = kBig;
+
+ // Track interception with plane. The coordinates are in the calorimeter's
+ // local system.
+
+ fOrigin = GetOrigin();
+ THcShower* fParent = (THcShower*) GetParent();
+ fParent->CalcTrackIntercept(BestTrack, pathl, XTrk, YTrk);
+
+ // Transform coordiantes to the spectrometer's coordinate system.
+ XTrk += GetOrigin().X();
+ YTrk += GetOrigin().Y();
+
+ for (Int_t i=0; i<fNelem; i++) {
+
+ if (TMath::Abs(XTrk - fParent->GetXPos(fLayerNum-1,i)) < fStatSlop &&
+ YTrk > fParent->GetYPos(fLayerNum-1,1) &&
+ YTrk < fParent->GetYPos(fLayerNum-1,0) ) {
+
+ fStatNumTrk.at(i)++;
+ fTotStatNumTrk++;
+
+ if (fGoodPosAdcPulseInt.at(i) > 0. || fGoodNegAdcPulseInt.at(i) > 0.) {
+ fStatNumHit.at(i)++;
+ fTotStatNumHit++;
+ }
+
+ }
+
+ }
+
+ if ( ((THcShower*) GetParent())->fdbg_tracks_cal ) {
+ cout << "---------------------------------------------------------------\n";
+ cout << "THcShowerPlane::AccumulateStat:" << endl;
+ cout << " Chi2/NDF = " <<BestTrack->GetChi2()/BestTrack->GetNDoF() << endl;
+ cout << " HGCER Npe = " << hgcer->GetCerNPE() << endl;
+ cout << " XTrk, YTrk = " << XTrk << " " << YTrk << endl;
+ for (Int_t i=0; i<fNelem; i++) {
+ if (TMath::Abs(XTrk - fParent->GetXPos(fLayerNum-1,i)) < fStatSlop) {
+
+ cout << " Module " << i << ", X=" << fParent->GetXPos(fLayerNum-1,i)
+ << " matches track" << endl;
+
+ if (fGoodPosAdcPulseInt.at(i) > 0. || fGoodNegAdcPulseInt.at(i) > 0.)
+ cout << " PulseIntegrals = " << fGoodPosAdcPulseInt.at(i) << " "
+ << fGoodNegAdcPulseInt.at(i) << endl;
+ }
+ }
+ cout << "---------------------------------------------------------------\n";
+ // getchar();
+ }
+
+ return 1;
+}
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 6a029e13ee808848d1b85d63eb58db573929c8a1..cc9f94ff62b8165e37383936b5b92a859a1df849 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -109,6 +109,8 @@ public:
return fNegPed[i];
};
+ Int_t AccumulateStat(TClonesArray& tracks);
+
protected:
// Flash ADC parameters
@@ -218,6 +220,17 @@ protected:
virtual void FillADC_SampleIntegral( );
virtual void FillADC_SampIntDynPed( );
virtual void FillADC_Standard( );
+
+ //Quatitites for efficiency calculations.
+
+ Double_t fStatCerMin;
+ Double_t fStatSlop;
+ Double_t fStatMaxChi2;
+ vector<Int_t> fStatNumTrk;
+ vector<Int_t> fStatNumHit;
+ Int_t fTotStatNumTrk;
+ Int_t fTotStatNumHit;
+
ClassDef(THcShowerPlane,0); // Calorimeter bars in a plane
};
#endif