diff --git a/src/THcHodoEff.cxx b/src/THcHodoEff.cxx
index 8a793926caade67ce1c541b77584f2d0d322747b..acb0b0a59c731de86dfc8e3cee6744601f876530 100644
--- a/src/THcHodoEff.cxx
+++ b/src/THcHodoEff.cxx
@@ -9,6 +9,17 @@
// For now trying to emulate work done in h_scin_eff/h_scin_eff_shutdown //
///////////////////////////////////////////////////////////////////////////////
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
#include "THcHodoEff.h"
#include "THaApparatus.h"
#include "THcHodoHit.h"
@@ -107,6 +118,120 @@ THaAnalysisObject::EStatus THcHodoEff::Init( const TDatime& run_time )
return fStatus = kOK;
}
+//_____________________________________________________________________________
+Int_t THcHodoEff::ReadDatabase( const TDatime& date )
+{
+ // Read database. Gets variable needed for efficiency calculation
+ // Get # of planes and their z positions here.
+
+ fNPlanes = fHod->GetNPlanes();
+ fPlanes = new THcScintillatorPlane* [fNPlanes];
+ fPosZ = new Double_t[fNPlanes];
+ fSpacing = new Double_t[fNPlanes];
+ fCenterFirst = new Double_t[fNPlanes];
+ fNCounters = new Int_t[fNPlanes];
+ fHodoSlop = new Double_t[fNPlanes];
+ fHodoPosEffi = new Int_t[100];
+ fHodoNegEffi = new Int_t[100];
+ fHodoOrEffi = new Int_t[100];
+ fHodoAndEffi = new Int_t[100];
+ fStatTrk = new Int_t[100];
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fPlanes[ip] = fHod->GetPlane(ip);
+ fPosZ[ip] = fPlanes[ip]->GetZpos() + 0.5*fPlanes[ip]->GetDzpos();
+ fSpacing[ip] = fPlanes[ip]->GetSpacing();
+ fCenterFirst[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ fNCounters[ip] = fPlanes[ip]->GetNelem();
+ }
+
+ char prefix[2];
+ prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
+ prefix[1] = '\0';
+
+ DBRequest list[]={
+ {"stat_slop", &fStatSlop, kDouble},
+ {"stat_maxchisq",&fMaxChisq, kDouble},
+ {"hodo_slop", fHodoSlop, kDouble, fNPlanes},
+ {0}
+ };
+ // fMaxShTrk = 0.05; // For cut on fraction of momentum seen in shower
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ cout << "\n\nTHcHodoEff::ReadDatabase nplanes=" << fHod->GetNPlanes() << endl;
+ // Setup statistics arrays
+ // Better method to put this in?
+ // These all need to be cleared in Begin
+ fHitPlane = new Int_t[fNPlanes];
+ fStatTrkDel.resize(fNPlanes);
+ fStatAndHitDel.resize(fNPlanes);
+ fStatPosHit.resize(fNPlanes);
+ fStatNegHit.resize(fNPlanes);
+ fStatAndHit.resize(fNPlanes);
+ fStatOrHit.resize(fNPlanes);
+ fBothGood.resize(fNPlanes);
+ fPosGood.resize(fNPlanes);
+ fNegGood.resize(fNPlanes);
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+
+ cout << "Plane = " << ip + 1 << " counters = " << fNCounters[ip] << endl;
+
+ fStatTrkDel[ip].resize(fNCounters[ip]);
+ fStatAndHitDel[ip].resize(fNCounters[ip]);
+ fStatPosHit[ip].resize(fNCounters[ip]);
+ fStatNegHit[ip].resize(fNCounters[ip]);
+ fStatAndHit[ip].resize(fNCounters[ip]);
+ fStatOrHit[ip].resize(fNCounters[ip]);
+ fBothGood[ip].resize(fNCounters[ip]);
+ fPosGood[ip].resize(fNCounters[ip]);
+ fNegGood[ip].resize(fNCounters[ip]);
+ for(Int_t ic=0;ic<fNCounters[ip];ic++) {
+ fStatTrkDel[ip][ic].resize(20); // Max this settable
+ fStatAndHitDel[ip][ic].resize(20); // Max this settable
+
+ fHodoPosEffi[fHod->GetScinIndex(ip,ic)] = 0;
+ fHodoNegEffi[fHod->GetScinIndex(ip,ic)] = 0;
+ fHodoOrEffi[fHod->GetScinIndex(ip,ic)] = 0;
+ fHodoAndEffi[fHod->GetScinIndex(ip,ic)] = 0;
+ fStatTrk[fHod->GetScinIndex(ip,ic)] = 0;
+ }
+ }
+
+ // Int_t fHodPaddles = fNCounters[0];
+ // gHcParms->Define(Form("%shodo_pos_hits[%d][%d]",fPrefix,fNPlanes,fHodPaddles),
+ // "Golden track's pos pmt hit",*&fStatPosHit);
+
+ Int_t fTotalPaddles = 100;
+ gHcParms->Define(Form("%shodo_pos_eff[%d]", prefix,fTotalPaddles), "Hodo positive effi",*fHodoPosEffi);
+ gHcParms->Define(Form("%shodo_neg_eff[%d]", prefix,fTotalPaddles), "Hodo negative effi",*fHodoNegEffi);
+ gHcParms->Define(Form("%shodo_or_eff[%d]", prefix,fTotalPaddles), "Hodo or effi", *fHodoOrEffi);
+ gHcParms->Define(Form("%shodo_and_eff[%d]", prefix,fTotalPaddles), "Hodo and effi", *fHodoAndEffi);
+ gHcParms->Define(Form("%shodo_gold_hits[%d]",prefix,fTotalPaddles), "Hodo golden hits", *fStatTrk);
+
+ return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoEff::DefineVariables( EMode mode )
+{
+
+ if( mode == kDefine && fIsSetup ) return kOK;
+ fIsSetup = ( mode == kDefine );
+
+ fEffiTest = 0;
+ gHcParms->Define(Form("hodoeffi"),"Testing effi",fEffiTest);
+
+ const RVarDef vars[] = {
+ // Move these into THcHallCSpectrometer using track fTracks
+ // {"effitestvar", "efficiency test var", "fEffiTest"},
+ // {"goldhodposhit", "pos pmt hit in hodo", "fStatPosHit"},
+ { 0 }
+ };
+ return DefineVarsFromList( vars, mode );
+ // return kOK;
+}
+
//_____________________________________________________________________________
Int_t THcHodoEff::Process( const THaEvData& evdata )
{
@@ -116,6 +241,7 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
if( !IsOK() ) return -1;
+
// Project the golden track to each
// plane. Need to get track at Focal Plane, not tgt.
//
@@ -127,7 +253,7 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
if(!theTrack) return 0;
Int_t trackIndex = theTrack->GetTrkNum()-1;
-
+
// May make these member variables
Double_t hitPos[fNPlanes];
Double_t hitDistance[fNPlanes];
@@ -140,19 +266,25 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
// Should really have plane object self identify as X or Y
if(ip%2 == 0) { // X Plane
hitPos[ip] = theTrack->GetX() + theTrack->GetTheta()*fPosZ[ip];
- } else { // Y Plane
- hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip];
- }
- // Which counter does track pass through?
- hitCounter[ip] = TMath::Min(TMath::Max(TMath::Nint((hitPos[ip] - fCenterFirst[ip])/fSpacing[ip]) + 1,1),fNCounters[ip]);
- // How far from paddle center is track?
- if(ip%2 == 0) { // X Plane
+ hitCounter[ip] = TMath::Max(
+ TMath::Min(
+ TMath::Nint((hitPos[ip]-fCenterFirst[ip])/
+ fSpacing[ip]+1),
+ TMath::Nint( fNCounters[ip] )),1);
hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) +
fCenterFirst[ip]);
} else { // Y Plane
- hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) -
- fCenterFirst[ip]);
+ hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip];
+ hitCounter[ip] = TMath::Max(
+ TMath::Min(
+ TMath::Nint((fCenterFirst[ip]-hitPos[ip])/
+ fSpacing[ip]+1),
+ TMath::Nint( fNCounters[ip] )),1);
+ hitDistance[ip] = hitPos[ip] -(fCenterFirst[ip] -
+ fSpacing[ip]*(hitCounter[ip]-1));
}
+
+
}
// Fill dpos histograms and set checkHit for each plane.
@@ -186,27 +318,29 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
Int_t hitcounter = hitCounter[ip];
Double_t dist = hitDistance[ip];
if(TMath::Abs(dist) <= fStatSlop &&
- theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
- // && hsshtrk >= 0.05
- Double_t delta = theTrack->GetDp();
- Int_t idel = TMath::Floor(delta+10.0);
- // Should
- if(idel >=0 && idel < 20) {
- fStatTrkDel[ip][hitcounter][idel]++;
+ theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq &&
+ theTrack->GetEnergy() >= 0.05 )
+ {
+ fStatTrk[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ Double_t delta = theTrack->GetDp();
+ Int_t idel = TMath::Floor(delta+10.0);
+ // Should
+ // if(idel >=0 && idel < 20) {
+ // fStatTrkDel[ip][hitcounter][idel]++;
+ // }
+ // lookat[ip] = TRUE;
}
- // lookat[ip] = TRUE;
- }
fHitPlane[ip] = 0;
- }
+ }
// Is there a hit on or adjacent to paddle that track
// passes through?
-
+
// May collapse this loop into last
-
+
// record the hits as a "didhit" if track is near center of
// scintillator, the chisqared of the track is good and it is the
// first "didhit" in that plane.
-
+
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t hitcounter = hitCounter[ip];
Double_t dist = hitDistance[ip];
@@ -219,13 +353,14 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
Bool_t onTrack, goodScinTime, goodTdcNeg, goodTdcPos;
fHod->GetFlags(trackIndex,ip,ihit,
onTrack, goodScinTime, goodTdcNeg, goodTdcPos);
+
if(TMath::Abs(dist) <= fStatSlop &&
TMath::Abs(hitcounter-counter) <= checkHit[ip] &&
fHitPlane[ip] == 0 &&
- theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
- // && hsshtrk >= 0.05
+ theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq &&
+ theTrack->GetEnergy() >= 0.05 ) {
fHitPlane[ip]++;
-
+
// Need to find out hgood_tdc_pos(igoldentrack,ihit) and neg
if(goodTdcPos) {
if(goodTdcNeg) { // Both fired
@@ -233,131 +368,78 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
fStatNegHit[ip][hitcounter]++;
fStatAndHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
+
+ fHodoPosEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ fHodoNegEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ fHodoAndEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+
Double_t delta = theTrack->GetDp();
- Int_t idel = TMath::Floor(delta+10.0);
- if(idel >=0 && idel < 20) {
- fStatAndHitDel[ip][hitcounter][idel]++;
- }
+ // Int_t idel = TMath::Floor(delta+10.0);
+ // if(idel >=0 && idel < 20) {
+ // fStatAndHitDel[ip][hitcounter][idel]++;
+ // }
} else {
fStatPosHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
+ fHodoPosEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
}
- } else if(goodTdcNeg) {
+ } else if (goodTdcNeg) {
fStatNegHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
+ fHodoNegEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
+ fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
}
- }
- // Increment pos/neg/both fired. Track independent, so
- // no chisquared cut, but note that only scintillators on the
- // track are examined.
- if(goodTdcPos) {
- if(goodTdcNeg) {
- fBothGood[ip][hitcounter]++;
- } else {
- fPosGood[ip][hitcounter]++;
+
+ // Increment pos/neg/both fired. Track independent, so
+ // no chisquared cut, but note that only scintillators on the
+ // track are examined.
+ if(goodTdcPos) {
+ if(goodTdcNeg) {
+ fBothGood[ip][hitcounter]++;
+ } else {
+ fPosGood[ip][hitcounter]++;
+ }
+ } else if (goodTdcNeg) {
+ fNegGood[ip][hitcounter]++;
+ }
+ // Determine if one or both PMTs had a good tdc
+
+ if(goodTdcPos && goodTdcNeg) {
+ goodTdcBothSides[ip] = kTRUE;
+ }
+ if(goodTdcPos || goodTdcNeg) {
+ goodTdcOneSide[ip] = kTRUE;
}
- } else if (goodTdcNeg) {
- fNegGood[ip][hitcounter]++;
}
- // Determine if one or both PMTs had a good tdc
- if(goodTdcPos && goodTdcNeg) {
- goodTdcBothSides[ip] = kTRUE;
- }
- if(goodTdcPos || goodTdcNeg) {
- goodTdcOneSide[ip] = kTRUE;
- }
+ /*
+ For each plane, see of other 3 fired. This means that they were enough
+ to form a 3/4 trigger, and so the fraction of times this plane fired is
+ the plane trigger efficiency. NOTE: we only require a TDC hit, not a
+ TDC hit within the SCIN 3/4 trigger window, so high rates will make
+ this seem better than it is. Also, make sure we're not near the edge
+ of the hodoscope (at the last plane), using the same hhodo_slop param.
+ as for h_tof.f
+ NOTE ALSO: to make this check simpler, we are assuming that all planes
+ have identical active areas. y_scin = y_cent + y_offset, so shift track
+ position by offset for comparing to edges.
+ */
+
+ // Need to add calculation and cuts on
+ // xatback and yatback in order to set the
+ // htrig_hododidflag, htrig_hodoshouldflag and otherthreehit flags
+ //
+
+ ++fNevt;
}
}
- /*
- For each plane, see of other 3 fired. This means that they were enough
- to form a 3/4 trigger, and so the fraction of times this plane fired is
- the plane trigger efficiency. NOTE: we only require a TDC hit, not a
- TDC hit within the SCIN 3/4 trigger window, so high rates will make
- this seem better than it is. Also, make sure we're not near the edge
- of the hodoscope (at the last plane), using the same hhodo_slop param. as for h_tof.f
- NOTE ALSO: to make this check simpler, we are assuming that all planes
- have identical active areas. y_scin = y_cent + y_offset, so shift track
- position by offset for comparing to edges.
- */
-
- // Need to add calculation and cuts on
- // xatback and yatback in order to set the
- // htrig_hododidflag, htrig_hodoshouldflag and otherthreehit flags
- //
-
- ++fNevt;
-
return 0;
}
//_____________________________________________________________________________
-Int_t THcHodoEff::ReadDatabase( const TDatime& date )
-{
- // Read database. Gets variable needed for efficiency calculation
- // Get # of planes and their z positions here.
-
- fNPlanes = fHod->GetNPlanes();
- fPlanes = new THcScintillatorPlane* [fNPlanes];
- fPosZ = new Double_t[fNPlanes];
- fSpacing = new Double_t[fNPlanes];
- fCenterFirst = new Double_t[fNPlanes];
- fNCounters = new Int_t[fNPlanes];
- fHodoSlop = new Double_t[fNPlanes];
-
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- fPlanes[ip] = fHod->GetPlane(ip);
- fPosZ[ip] = fPlanes[ip]->GetZpos() + 0.5*fPlanes[ip]->GetDzpos();
- fSpacing[ip] = fPlanes[ip]->GetSpacing();\
- fCenterFirst[ip] = fPlanes[ip]->GetPosCenter(0);
- fNCounters[ip] = fPlanes[ip]->GetNelem();
- }
-
- char prefix[2];
- prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
- prefix[1] = '\0';
-
- DBRequest list[]={
- {"stat_slop", &fStatSlop, kDouble},
- {"stat_maxchisq",&fMaxChisq, kDouble},
- {"hodo_slop", fHodoSlop, kDouble, fNPlanes},
- {0}
- };
- // fMaxShTrk = 0.05; // For cut on fraction of momentum seen in shower
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- cout << "THcHodoEff::ReadDatabase nplanes=" << fHod->GetNPlanes() << endl;
- // Setup statistics arrays
- // Better method to put this in?
- // These all need to be cleared in Begin
- fHitPlane = new Int_t[fNPlanes];
- fStatTrkDel.resize(fNPlanes);
- fStatAndHitDel.resize(fNPlanes);
- fStatPosHit.resize(fNPlanes);
- fStatNegHit.resize(fNPlanes);
- fStatAndHit.resize(fNPlanes);
- fStatOrHit.resize(fNPlanes);
- fBothGood.resize(fNPlanes);
- fPosGood.resize(fNPlanes);
- fNegGood.resize(fNPlanes);
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- fStatTrkDel[ip].resize(fNCounters[ip]);
- fStatAndHitDel[ip].resize(fNCounters[ip]);
- fStatPosHit[ip].resize(fNCounters[ip]);
- fStatNegHit[ip].resize(fNCounters[ip]);
- fStatAndHit[ip].resize(fNCounters[ip]);
- fStatOrHit[ip].resize(fNCounters[ip]);
- fBothGood[ip].resize(fNCounters[ip]);
- fPosGood[ip].resize(fNCounters[ip]);
- fNegGood[ip].resize(fNCounters[ip]);
- for(Int_t ic=0;ic<fNCounters[ip];ic++) {
- fStatTrkDel[ip][ic].resize(20); // Max this settable
- fStatAndHitDel[ip][ic].resize(20); // Max this settable
- }
- }
- return kOK;
-}
-//_____________________________________________________________________________
ClassImp(THcHodoEff)
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcHodoEff.h b/src/THcHodoEff.h
index eaed7f07f95f7e23365abfbde41428daffd59251..ec7e2d7c7f5f6ad08738d5b049c88673b49b0cf6 100644
--- a/src/THcHodoEff.h
+++ b/src/THcHodoEff.h
@@ -7,6 +7,17 @@
// //
///////////////////////////////////////////////////////////////////////////////
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
#include "THaPhysicsModule.h"
#include "THcHodoscope.h"
#include "THaSpectrometer.h"
@@ -27,6 +38,8 @@ public:
protected:
virtual Int_t ReadDatabase( const TDatime& date);
+ virtual Int_t DefineVariables( EMode mode = kDefine );
+ /* Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle); */
// Data needed for efficiency calculation for one Hodoscope paddle
@@ -41,13 +54,19 @@ protected:
// Information about the hodoscopes that we get from the
// THcHodoscope object
+ Int_t fEffiTest;
Int_t fNPlanes;
THcScintillatorPlane** fPlanes;
Double_t* fPosZ;
Double_t* fSpacing;
Double_t* fCenterFirst;
Int_t* fNCounters;
- // Int_t fMaxNcounters;
+ // Int_t* fHodoPlnContHit;
+ Int_t* fHodoPosEffi;
+ Int_t* fHodoNegEffi;
+ Int_t* fHodoOrEffi;
+ Int_t* fHodoAndEffi;
+ Int_t* fStatTrk;
Double_t fStatSlop;
Double_t fMaxChisq;
Double_t* fHodoSlop;