diff --git a/examples/PARAM/52949/hcal.param.vt.52949 b/examples/PARAM/52949/hcal.param.vt.52949
index 8e5f719ae9da73b3c944f4cd480a6fa016ccdaac..3970fb3c358fdc38f4f883f2d8ab5276b27a2c1b 100644
--- a/examples/PARAM/52949/hcal.param.vt.52949
+++ b/examples/PARAM/52949/hcal.param.vt.52949
@@ -4,7 +4,7 @@ hcal_slop = 7.5
;Turn on HMS cal. fiducial volume cut. 0="no cut"
;Default hcal_fv_test=0
-hcal_fv_test = 0
+hcal_fv_test = 1
hcal_pos_cal_const =0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001
diff --git a/examples/PARAM/52949/hcana.param b/examples/PARAM/52949/hcana.param
index 61d60c67416cd2b92a6e69416e20b1281f30689b..9d0dd731686ec55d14980c3a739ba254c91b5aeb 100644
--- a/examples/PARAM/52949/hcana.param
+++ b/examples/PARAM/52949/hcana.param
@@ -7,6 +7,8 @@ hhodo_num_planes = 4
hcal_num_layers = 4
+hcal_fv_delta = 5.
+
haero_num_pairs = 8
# Names of planes so that parameter names can be constructed
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 1d10da96accea46954445bf1200c9f289e65e91b..f96bc72110ededc4d15a7e597dd08e3cf49a8328 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -165,17 +165,19 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
DBRequest list[]={
{"cal_num_neg_columns", &fNegCols, kInt},
{"cal_slop", &fSlop, kDouble},
- {"cal_fv_test", &fvTest, kDouble},
+ {"cal_fv_test", &fvTest, kInt},
+ {"cal_fv_delta", &fvDelta, kDouble},
{"dbg_clusters_cal", &fdbg_clusters_cal, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
}
- cout << "Number of neg. columns = " << fNegCols << endl;
- cout << "Slop parameter = " << fSlop << endl;
- cout << "Fiducial volum test flag = " << fvTest << endl;
- cout << "Cluster debug flag = " << fdbg_clusters_cal << endl;
+ cout << "Number of neg. columns = " << fNegCols << endl;
+ cout << "Slop parameter = " << fSlop << endl;
+ cout << "Fiducial volume test flag = " << fvTest << endl;
+ cout << "Fiducial volume excl. width = " << fvDelta << endl;
+ cout << "Cluster debug flag = " << fdbg_clusters_cal << endl;
BlockThick = new Double_t [fNLayers];
fNBlocks = new Int_t [fNLayers];
@@ -220,6 +222,18 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
}
+ // Fiducial volume limits, based on Plane 1 positions.
+ //
+
+ fvXmin = XPos[0][0] + fvDelta;
+ fvXmax = XPos[0][fNBlocks[0]-1] + BlockThick[0] - fvDelta;
+ fvYmin = YPos[0] + fvDelta;
+ fvYmax = YPos[1] - fvDelta;
+
+ cout << "Fiducial volume limits:" << endl;
+ cout << " Xmin = " << fvXmin << " Xmax = " << fvXmax << endl;
+ cout << " Ymin = " << fvYmin << " Ymax = " << fvYmax << endl;
+
//Calibration related parameters (from hcal.param).
fNtotBlocks=0; //total number of blocks
@@ -716,68 +730,77 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
<< " Phi = " << Track->GetPhi()
<< endl;
- Double_t xtrk = kBig;
- Double_t ytrk = kBig;
+ Double_t XTrFront = kBig;
+ Double_t YTrFront = kBig;
Double_t pathl = kBig;
// Track interception with face of the calorimeter. The coordinates are
// in the calorimeter's local system.
- CalcTrackIntercept(Track, pathl, xtrk, ytrk);
+ fOrigin = fPlanes[0]->GetOrigin();
+
+ CalcTrackIntercept(Track, pathl, XTrFront, YTrFront);
// Transform coordiantes to the spectrometer's coordinate system.
- xtrk += GetOrigin().X();
- ytrk += GetOrigin().Y();
+ XTrFront += GetOrigin().X();
+ YTrFront += GetOrigin().Y();
- cout << "Track at Calorimeter:"
- << " X = " << xtrk
- << " Y = " << ytrk
+ cout << "Track at the front of Calorimeter:"
+ << " X = " << XTrFront
+ << " Y = " << YTrFront
<< " Pathl = " << pathl
<< endl;
- // Match a cluster to the track.
+ Bool_t inFidVol = true;
- Int_t mclust = -1; // The match cluster #, initialize with a bogus value.
- Double_t deltaX = kBig; // Track to cluster distance
+ if (fvTest) {
- // hcal_zmin= hcal_1pr_zpos
- // hcal_zmax= hcal_4ta_zpos
- // hcal_fv_xmin=hcal_xmin+5.
- // hcal_fv_xmax=hcal_xmax-5.
- // hcal_fv_ymin=hcal_ymin+5.
- // hcal_fv_ymax=hcal_ymax-5.
- // hcal_fv_zmin=hcal_zmin
- // hcal_fv_zmax=hcal_zmax
+ fOrigin = fPlanes[fNLayers-1]->GetOrigin();
- // dz_f=hcal_zmin-hz_fp(nt)
- // dz_b=hcal_zmax-hz_fp(nt)
+ Double_t XTrBack = kBig;
+ Double_t YTrBack = kBig;
- // xf=hx_fp(nt)+hxp_fp(nt)*dz_f
- // xb=hx_fp(nt)+hxp_fp(nt)*dz_b
+ CalcTrackIntercept(Track, pathl, XTrBack, YTrBack);
- // yf=hy_fp(nt)+hyp_fp(nt)*dz_f
- // yb=hy_fp(nt)+hyp_fp(nt)*dz_b
+ XTrBack += GetOrigin().X();
+ YTrBack += GetOrigin().Y();
- // track_in_fv = (xf.le.hcal_fv_xmax .and. xf.ge.hcal_fv_xmin .and.
- // & xb.le.hcal_fv_xmax .and. xb.ge.hcal_fv_xmin .and.
- // & yf.le.hcal_fv_ymax .and. yf.ge.hcal_fv_ymin .and.
- // & yb.le.hcal_fv_ymax .and. yb.ge.hcal_fv_ymin)
+ cout << "Track at the back of Calorimeter:"
+ << " X = " << XTrBack
+ << " Y = " << YTrBack
+ << " Pathl = " << pathl
+ << endl;
- for (Int_t i=0; i<fNclust; i++) {
+ inFidVol = (XTrFront < fvXmax) && (XTrFront > fvXmin) &&
+ (YTrFront < fvYmax) && (YTrFront > fvYmin) &&
+ (XTrBack < fvXmax) && (XTrBack > fvXmin) &&
+ (YTrBack < fvYmax) && (YTrBack > fvYmin);
- THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+ cout << "Fiducial volume test: inFidVol = " << inFidVol << endl;
+ }
- Double_t dx = TMath::Abs( (*cluster).clX() - xtrk );
+ // Match a cluster to the track.
- if (dx <= (0.5*BlockThick[0] + fSlop)) {
+ Int_t mclust = -1; // The match cluster #, initialize with a bogus value.
+ Double_t deltaX = kBig; // Track to cluster distance
+
+ if (inFidVol) {
- if (dx <= deltaX) {
- mclust = i;
- deltaX = dx;
+ for (Int_t i=0; i<fNclust; i++) {
+
+ THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+
+ Double_t dx = TMath::Abs( (*cluster).clX() - XTrFront );
+
+ if (dx <= (0.5*BlockThick[0] + fSlop)) {
+
+ if (dx <= deltaX) {
+ mclust = i;
+ deltaX = dx;
+ }
}
}
-
}
cout << "MatchCluster: mclust= " << mclust << " delatX= " << deltaX << endl;
diff --git a/src/THcShower.h b/src/THcShower.h
index 1cde6022556c75ae1afb7c17fdee3a6ba41bceca..9bfe9abf42c0d109f2118c8dd66f6a0e1409d1e0 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -36,6 +36,22 @@ public:
Int_t GetNBlocks(Int_t NLayer) const { return fNBlocks[NLayer];}
+ Double_t GetXPos(Int_t NLayer, Int_t NRaw) const {
+ return XPos[NLayer][NRaw];
+ }
+
+ Double_t GetYPos(Int_t NLayer, Int_t Side) const {
+
+ //Side = 0 for postive (left) side
+ //Side = 1 for negative (right) side
+
+ return YPos[2*NLayer+(1-Side)];
+ }
+
+ Double_t GetZPos(Int_t NLayer) const {return fNLayerZPos[NLayer];}
+
+ Double_t GetBlockThick(Int_t NLayer) {return BlockThick[NLayer];}
+
// friend class THaScCalib; not needed for now.
Int_t GetPedLimit(Int_t NBlock, Int_t NLayer, Int_t Side) {
@@ -110,6 +126,12 @@ public:
Int_t fNegCols; // # of columns with neg. side PMTs only.
Double_t fSlop; // Track to cluster vertical slop distance.
Int_t fvTest; // fiducial volume test flag for tracking
+ Double_t fvDelta; // Exclusion band width for fiducial volume
+
+ Double_t fvXmin; // Fiducial volume limits
+ Double_t fvXmax;
+ Double_t fvYmin;
+ Double_t fvYmax;
Int_t fdbg_clusters_cal; // Shower debug flag
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index d579a85ba14255d58fd50dc36403e6e470b62c4c..66cec1e4e161bd60916fbff9f4add6c90e69fe8b 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -16,6 +16,8 @@
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
+#include "THaTrack.h"
+#include "THaTrackProj.h"
#include <cstring>
#include <cstdio>
@@ -115,6 +117,30 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
// cout << "THcShowerPlane::ReadDatabase: fLayerNum=" << fLayerNum
// << " fNelem=" << fNelem << endl;
+ // Origin of the plane.
+
+ Double_t BlockThick = fParent->GetBlockThick(fLayerNum-1);
+
+ Double_t xOrig = (fParent->GetXPos(fLayerNum-1,0) +
+ fParent->GetXPos(fLayerNum-1,fNelem-1))/2 +
+ BlockThick/2;
+
+ Double_t yOrig = (fParent->GetYPos(fLayerNum-1,0) +
+ fParent->GetYPos(fLayerNum-1,1))/2;
+
+ Double_t zOrig = fParent->GetZPos(fLayerNum-1);
+
+ fOrigin.SetXYZ(xOrig, yOrig, zOrig);
+
+ cout << "Origin of Layer " << fLayerNum << ": "
+ << fOrigin.X() << " " << fOrigin.Y() << " " << fOrigin.Z() << endl;
+
+ // Detector axes. Assume no rotation.
+ //
+ // DefineAxes(0.); not for subdetector
+
+ //
+
fA_Pos = new Double_t[fNelem];
fA_Neg = new Double_t[fNelem];
fA_Pos_p = new Double_t[fNelem];
@@ -216,6 +242,22 @@ Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks )
cout << "THcShowerPlane::CoarseProcess called ---------------------" << endl;
+ Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+
+ cout << " Number of reconstructed tracks = " << Ntracks << endl;
+
+ /*
+ for (Int_t i=0; i<Ntracks; i++) {
+
+ THaTrack* theTrack = static_cast<THaTrack*>( tracks[i] );
+
+ Double_t pathl;
+ Double_t xtrk;
+ Double_t ytrk;
+ CalcTrackIntercept(theTrack, pathl, xtrk, ytrk);
+ }
+ */
+
return 0;
}
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 9f2467c0d669937c685729fa8cbbe0df8bf5843e..6c0aa5124fc13cfa0cc8c62fc38c44bcec20179d 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -46,6 +46,10 @@ class THcShowerPlane : public THaSubDetector {
TClonesArray* fParentHitList;
+ TVector3 GetOrigin() {
+ return fOrigin;
+ }
+
Double_t GetEplane() {
return fEplane;
};