Merge branch 'shower4' into develop

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
@@ -15,3 +17,6 @@ hdc_plane_names = "1x1 1y1 1u1 1v1 1y2 1x2 2x1 2y1 2u1 2v1 2y2 2x2"
 hcal_layer_names = "1pr 2ta 3ta 4ta"
 
 hhodo_plane_names = "1x 1y 2x 2y"
+
+hpcentral=0.811
+

examples/comptrackedep.C

+void comptrackedep(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = tre;
+
+  TFile* f1 = new TFile(Form("%d_hbk.root",run));
+  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+  TH1F* h1;
+  switch (run) {
+  case 50017 :
+  //    h1 = h212;   //A+
+    break;
+  default :
+    h1 = h420;   //edep
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower 1 track cluster Edep", 1000, 667); 
+
+  //  gPad->SetLogy();
+
+  h1->SetFillColor(kGreen);
+  h1->SetLineColor(kGreen);
+  h1->Draw();
+
+  h->SetLineColor(kBlue);
+  h->SetFillStyle(0);
+  h->SetLineWidth(2);
+  h->Draw("same");
+
+  TLatex l;
+  l.SetTextSize(0.04);
+  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+  Float_t xmin = h1->GetXaxis()->GetXmin();
+  Float_t xmax = h1->GetXaxis()->GetXmax();
+  Float_t xt = xmin + 0.67*(xmax-xmin);
+
+  l.SetTextColor(kGreen);
+  l.DrawLatex(xt,0.65*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.75*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Edep differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Edep Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/comptrackepr.C

+void comptrackepr(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = trepr;
+
+  TFile* f1 = new TFile(Form("%d_hbk.root",run));
+  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+  TH1F* h1;
+  switch (run) {
+  case 50017 :
+  //    h1 = h212;   //A+
+    break;
+  default :
+    h1 = h421;   //epr
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster Edep", 1000, 667); 
+
+  //  gPad->SetLogy();
+
+  h1->SetFillColor(kGreen);
+  h1->SetLineColor(kGreen);
+  h1->Draw();
+
+  h->SetLineColor(kBlue);
+  h->SetFillStyle(0);
+  h->SetLineWidth(2);
+  h->Draw("same");
+
+  TLatex l;
+  l.SetTextSize(0.04);
+  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+  Float_t xmin = h1->GetXaxis()->GetXmin();
+  Float_t xmax = h1->GetXaxis()->GetXmax();
+  Float_t xt = xmin + 0.75*(xmax-xmin);
+
+  l.SetTextColor(kGreen);
+  l.DrawLatex(xt,0.65*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.75*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Epr differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Epr Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/comptrackx.C

+void comptrackx(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = trx;
+
+  TFile* f1 = new TFile(Form("%d_hbk.root",run));
+  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+  TH1F* h1;
+  switch (run) {
+  case 50017 :
+  //    h1 = h212;   //A+
+    break;
+  default :
+    h1 = h417;   //X
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster X", 1000, 667); 
+
+  //  gPad->SetLogy();
+
+  h1->SetFillColor(kGreen);
+  h1->SetLineColor(kGreen);
+  h1->Draw();
+
+  h->SetLineColor(kBlue);
+  h->SetFillStyle(0);
+  h->SetLineWidth(2);
+  h->Draw("same");
+
+  TLatex l;
+  l.SetTextSize(0.04);
+  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+  Float_t xmin = h1->GetXaxis()->GetXmin();
+  Float_t xmax = h1->GetXaxis()->GetXmax();
+  Float_t xt = xmin + 0.67*(xmax-xmin);
+
+  l.SetTextColor(kGreen);
+  l.DrawLatex(xt,0.65*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.75*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Track X differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Hcal track X Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/comptrackxy.C

+void comptrackxy(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH2F* h = trx_vs_try;
+
+  TFile* f1 = new TFile(Form("%d_hbk.root",run));
+  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+  TH2F* h1;
+  switch (run) {
+  case 50017 :
+  //    h1 = h212;   //A+
+    break;
+  default :
+    h1 = h419;   //Y vs Y
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster X", 1000, 667); 
+
+  //  gPad->SetLogy();
+
+  h1->SetMarkerColor(kGreen);
+  //  h1->SetLineColor(kGreen);
+  h1->Draw();
+
+  h->SetMarkerColor(kBlue);
+  //  h->SetFillStyle(0);
+  //  h->SetLineWidth(2);
+  h->Draw("same");
+
+  /*
+  TLatex l;
+  l.SetTextSize(0.04);
+  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+  Float_t xmin = h1->GetXaxis()->GetXmin();
+  Float_t xmax = h1->GetXaxis()->GetXmax();
+  Float_t xt = xmin + 0.67*(xmax-xmin);
+
+  l.SetTextColor(kGreen);
+  l.DrawLatex(xt,0.095*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.045*maxy,"hcana");
+  */
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "X differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  //  dif->SetTitle("X Difference");
+  //  dif->SetFillColor(kRed);
+  //  dif->SetLineColor(kRed);
+  //  dif->SetLineWidth(1);
+  //  dif->SetFillStyle(1111);
+  dif->Draw("LEGO2");
+} 

examples/comptracky.C

+void comptracky(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = try;
+
+  TFile* f1 = new TFile(Form("%d_hbk.root",run));
+  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+  TH1F* h1;
+  switch (run) {
+  case 50017 :
+  //    h1 = h212;   //A+
+    break;
+  default :
+    h1 = h418;   //Y
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower track Y", 1000, 667); 
+
+  //  gPad->SetLogy();
+
+  h1->SetFillColor(kGreen);
+  h1->SetLineColor(kGreen);
+  h1->Draw();
+
+  h->SetLineColor(kBlue);
+  h->SetFillStyle(0);
+  h->SetLineWidth(2);
+  h->Draw("same");
+
+  TLatex l;
+  l.SetTextSize(0.04);
+  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+  Float_t xmin = h1->GetXaxis()->GetXmin();
+  Float_t xmax = h1->GetXaxis()->GetXmax();
+  Float_t xt = xmin + 0.67*(xmax-xmin);
+
+  l.SetTextColor(kGreen);
+  l.DrawLatex(xt,0.65*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.755*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Hcal Y differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Y Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/output_52949.def

@@ -63,6 +63,17 @@ TH1F eprmax 'HMS PrSh Max Clust. Edep' H.cal.eprmax 100 -0.1 0.811
 # X coordinate the largest cluster.
 TH1F xmax 'HMS Cal Max Clust. X' H.cal.xmax 150 -75. 75.
 
+# Track coordinates at the front of calorimeter
+TH2F trx_vs_try 'HMS Track X vs Y at calor.' H.cal.try H.cal.trx 80 -40. 40. 150 -75. 75. H.dc.ntrack==1
+TH1F trx 'HMS Track X at calorimeter' H.cal.trx 150 -75. 75. H.dc.ntrack==1
+TH1F try 'HMS Track Y at calorimeter' H.cal.try  80 -40. 40. H.dc.ntrack==1
+
+# Energy of the cluster associated to the main (currently first) track
+TH1F tre 'HMS Cal Track Clust. Edep' H.cal.tre 200 -0.1 1.622 H.dc.ntrack==1
+
+# Preshower Energy of the cluster associated to the main (currently 1st) track
+TH1F trepr 'HMS PrSh Track Clust. Edep' H.cal.trepr 100 -0.1 0.811 H.dc.ntrack==1
+
 #Calorimeter ADC channels
 TH1F hcaladc_A1p 'HMS Cal ADC A1p - PED' H.cal.1pr.apos_p[0] 150 50 500
 TH1F hcaladc_A2p 'HMS Cal ADC A2p - PED' H.cal.1pr.apos_p[1] 150 50 500

src/THcShower.cxx

@@ -74,7 +74,9 @@ void THcShower::Setup(const char* name, const char* description)
   vector<string> layer_names = vsplit(layernamelist);
 
   if(layer_names.size() != (UInt_t) fNLayers) {
-    cout << "ERROR: Number of layers " << fNLayers << " doesn't agree with number of layer names " << layer_names.size() << endl;
+    cout << "ERROR: Number of layers " << fNLayers 
+	 << " doesn't agree with number of layer names "
+	 << layer_names.size() << endl;
     // Should quit.  Is there an official way to quit?
   }
 
@@ -165,17 +167,25 @@ 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_decoded_cal", &fdbg_decoded_cal, kInt},
+      {"dbg_sparsified_cal", &fdbg_sparsified_cal, kInt},
       {"dbg_clusters_cal", &fdbg_clusters_cal, kInt},
+      {"dbg_tracks_cal", &fdbg_tracks_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 << "Decode debug flag           = " << fdbg_decoded_cal  << endl;
+  cout << "Sparsify debug flag         = " << fdbg_sparsified_cal  << endl;
+  cout << "Cluster debug flag          = " << fdbg_clusters_cal  << endl;
+  cout << "Tracking debug flag         = " << fdbg_tracks_cal  << endl;
 
   BlockThick = new Double_t [fNLayers];
   fNBlocks = new Int_t [fNLayers];
@@ -220,6 +230,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
@@ -430,8 +452,10 @@ Int_t THcShower::DefineVariables( EMode mode )
     { "mult",   "Multiplicity of largest cluster",    "fMult" },
  //   { "nblk",   "Numbers of blocks in main cluster",  "fNblk" },
  //   { "eblk",   "Energies of blocks in main cluster", "fEblk" },
- //   { "trx",    "track x-position in det plane",      "fTRX" },
- //   { "try",    "track y-position in det plane",      "fTRY" },
+    { "trx",    "track x-position in det plane (1st track)",      "fTRX" },
+    { "try",    "track y-position in det plane (1st track)",      "fTRY" },
+    { "tre",    "track energy in the calorimeter (1st track)",    "fTRE" },
+    { "trepr",  "track energy in the preshower (1st track)",      "fTREpr" },
     { 0 }
   };
   return DefineVarsFromList( vars, mode );
@@ -486,7 +510,11 @@ void THcShower::Clear(Option_t* opt)
   fMult = 0;
   fE = 0.;
   fEpr = 0.;
-  fX = -75.;   //out of acceptance
+  fX = -75.;    //out of acceptance
+  fTRX = -75.;  //out of acceptance
+  fTRY = -40.;  //out of acceptance
+  fTRE = -1.;
+  fTREpr = -1.;
 }
 
 //_____________________________________________________________________________
@@ -502,7 +530,9 @@ Int_t THcShower::Decode( const THaEvData& evdata )
     Int_t nexthit = 0;
     for(Int_t ip=0;ip<fNLayers;ip++) {
       nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
-      //cout << "nexthit = " << nexthit << endl;
+      if (fdbg_decoded_cal) {
+	cout << "THcShower::Decode: nexthit = " << nexthit << endl;
+      }
     }
     fAnalyzePedestals = 1;	// Analyze pedestals first normal events
     return(0);
@@ -521,13 +551,16 @@ Int_t THcShower::Decode( const THaEvData& evdata )
   }
 
   //   fRawHitList is TClones array of THcRawShowerHit objects
-  //  cout << "THcShower::Decode: Shower raw hit list:" << endl;
-  //  for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
-  //    THcRawShowerHit* hit = (THcRawShowerHit *) fRawHitList->At(ihit);
-  //    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
-  //	 << hit->fADC_pos << " " << hit->fADC_neg << " "  << endl;
+
+  //  if (fdbg_decoded_cal) {
+  //    cout << "THcShower::Decode: Shower raw hit list:" << endl;
+  //    for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+  //      THcRawShowerHit* hit = (THcRawShowerHit *) fRawHitList->At(ihit);
+  //      cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+  //	   << hit->fADC_pos << " " << hit->fADC_neg << " "  << endl;
+  //    }
+  //    cout << endl;
   //  }
-  //  cout << endl;
 
   return nhits;
 }
@@ -565,7 +598,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
   // Clustering of hits.
   //
 
-  THcShowerHitList HitList;                    //list of unclusterd hits
+  THcShowerHitList HitList;                    //list of unclustered hits
 
   for(Int_t j=0; j < fNLayers; j++) {
 
@@ -672,30 +705,54 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
     fX = (*MaxCluster).clX();
   }
 
- if (fdbg_clusters_cal)
-    cout << fEpr << " " << fE << " " << fX << " PrSh" << endl;
+  if (fdbg_clusters_cal)
+    cout << "Max.cluster: fEpr = " << fEpr << " fE = " << fE << " fX = " << fX
+	 << endl;
 
   // Track-to-cluster  matching.
   //
 
   Int_t Ntracks = tracks.GetLast()+1;   // Number of reconstructed tracks
 
-  cout << "Number of reconstructed tracks = " << Ntracks << endl;
+  if (fdbg_tracks_cal) {
+    cout << endl;
+    cout << "Number of reconstructed tracks = " << Ntracks << endl;
+  }
 
   for (Int_t i=0; i<Ntracks; i++) {
 
     THaTrack* theTrack = static_cast<THaTrack*>( tracks[i] );
 
-    //    cout << "   Track " << i << ": "
-    //	 << "  X = " << theTrack->GetX()
-    //	 << "  Y = " << theTrack->GetY()
-    //	 << "  Theta = " << theTrack->GetTheta()
-    //	 << "  Phi = " << theTrack->GetPhi()
-    //	 << endl;
+    if (fdbg_tracks_cal) {
+      cout << "   Track " << i << ": "
+	   << "  X = " << theTrack->GetX()
+	   << "  Y = " << theTrack->GetY()
+	   << "  Theta = " << theTrack->GetTheta()
+	   << "  Phi = " << theTrack->GetPhi()
+	   << endl;
+    }
+
+    Double_t Xtr;
+    Double_t Ytr;
+
+    Int_t mclust = MatchCluster(theTrack, ClusterList, Xtr, Ytr);
+
+    if (i==0) {
+      fTRX = Xtr;
+      fTRY = Ytr;
+      if (mclust >= 0) {
+	THcShowerCluster* cluster = (*ClusterList).ListedCluster(mclust);
+	fTRE = (*cluster).clE();
+	fTREpr = (*cluster).clEpr();
+      }
+    }
 
-    MatchCluster(theTrack);
   }
 
+  if (fdbg_tracks_cal)
+    cout << "1st track cluster:  fTRX = " << fTRX << "  fTRY = " << fTRY 
+	 << endl;
+
   if (fdbg_clusters_cal)
   cout << "THcShower::CoarseProcess return ---------------------------" <<endl;
 
@@ -704,32 +761,102 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
 
 //-----------------------------------------------------------------------------
 
-void THcShower::MatchCluster(THaTrack* Track)
+Int_t THcShower::MatchCluster(THaTrack* Track,
+			      THcShowerClusterList* ClusterList,
+			      Double_t& XTrFront, Double_t& YTrFront)
 {
   // Match a cluster to a given track.
 
-  cout << "Track at DC:"
-       << "  X = " << Track->GetX()
-       << "  Y = " << Track->GetY()
-       << "  Theta = " << Track->GetTheta()
-       << "  Phi = " << Track->GetPhi()
-       << endl;
+  if (fdbg_tracks_cal) {
+    cout << "THcShower::MatchCluster: Track at DC:"
+	 << "  X = " << Track->GetX()
+	 << "  Y = " << Track->GetY()
+	 << "  Theta = " << Track->GetTheta()
+	 << "  Phi = " << Track->GetPhi()
+	 << endl;
+  }
+
+  XTrFront = kBig;
+  YTrFront = kBig;
+  Double_t pathl = kBig;
+
+  // Track interception with face of the calorimeter. The coordinates are
+  // in the calorimeter's local system.
+
+  fOrigin = fPlanes[0]->GetOrigin();
+
+  CalcTrackIntercept(Track, pathl, XTrFront, YTrFront);
+
+  // Transform coordiantes to the spectrometer's coordinate system.
 
-  Double_t xc = 1.e8;
-  Double_t yc = 1.e8;
-  Double_t pathl = 1.e8;
+  XTrFront += GetOrigin().X();
+  YTrFront += GetOrigin().Y();
 
-  CalcTrackIntercept(Track, pathl, xc, yc);
+  if (fdbg_tracks_cal)
+    cout << "  Track at the front of Calorimeter:"
+	 << "  X = " << XTrFront
+	 << "  Y = " << YTrFront
+	 << "  Pathl = " << pathl
+	 << endl;
+
+  Bool_t inFidVol = true;
+
+  if (fvTest) {
+
+    fOrigin = fPlanes[fNLayers-1]->GetOrigin();
+
+    Double_t XTrBack = kBig;
+    Double_t YTrBack = kBig;
+
+    CalcTrackIntercept(Track, pathl, XTrBack, YTrBack);
+
+    XTrBack += GetOrigin().X();
+    YTrBack += GetOrigin().Y();
+
+    if (fdbg_tracks_cal)
+      cout << "  Track at the back of Calorimeter:"
+	   << "  X = " << XTrBack
+	   << "  Y = " << YTrBack
+	   << "  Pathl = " << pathl
+	   << endl;
+
+    inFidVol = (XTrFront < fvXmax) && (XTrFront > fvXmin) &&
+               (YTrFront < fvYmax) && (YTrFront > fvYmin) &&
+               (XTrBack < fvXmax) && (XTrBack > fvXmin) &&
+               (YTrBack < fvYmax) && (YTrBack > fvYmin);
+
+    if (fdbg_tracks_cal) 
+      cout << "  Fiducial volume test: inFidVol = " << inFidVol << endl;
+  }
+
+  // Match a cluster to the track.
+
+  Int_t mclust = -1;    // The match cluster #, initialize with a bogus value.
+  Double_t deltaX = kBig;   // Track to cluster distance
+
+  if (inFidVol) {
+
+    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;
+	}
+      }
+    }
+  }
 
-  //  Double_t dx = 0.;
-  //  Int_t pad = -1;
-  //  new ( (*fTrackProj)[0]) THaTrackProj(xc,yc,pathl,dx,pad,this);
+  if (fdbg_tracks_cal)
+    cout << "MatchCluster: mclust= " << mclust << "  delatX= " << deltaX 
+	 << endl;
 
-  cout << "Track at Calorimeter:"
-       << "  X = " << xc
-       << "  Y = " << yc
-       << "  Pathl = " << pathl
-       << endl;
+  return mclust;
 }
 
 

src/THcShower.h

@@ -11,6 +11,7 @@
 #include "THaNonTrackingDetector.h"
 #include "THcHitList.h"
 #include "THcShowerPlane.h"
+#include "THcShowerCluster.h"
 
 class THaScCalib;
 
@@ -35,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) {
@@ -90,10 +107,12 @@ public:
   Double_t fX;               // x-position (cm) of the largest cluster
   Double_t fZ;               // z-position (cm) of the largest cluster
   Int_t fMult;               // # of hits in the largest cluster
-  //  Int_t fNblk;      // Number of blocks in main cluster
-  //  Double_t* fEblk;  // Energies of blocks in main cluster
-  //  Double_t fTRX;    // track x-position in det plane"
-  //  Double_t fTRY;    // track y-position in det plane",
+  //  Int_t fNblk;           // Number of blocks in main cluster
+  //  Double_t* fEblk;       // Energies of blocks in main cluster
+  Double_t fTRX;             // track x-position in det plane (1st track)
+  Double_t fTRY;             // track y-position in det plane (1st track)
+  Double_t fTRE;             // Energy (MeV) of the cluster associated to track
+  Double_t fTREpr;           // Preshower Energy (MeV) of the track's cluster
 
   // Potential Hall C parameters.  Mostly here for demonstration
 
@@ -106,10 +125,21 @@ public:
   Double_t** XPos;		// [fNLayers] X,Y,Z positions of blocks
   Double_t* YPos;
   Double_t* ZPos;
-  Int_t fNegCols; //number of columns with PMTTs on the negative side only.
-  Double_t fSlop;               //Track to cluster vertical slop distance.
-  Int_t fvTest;                 //fiducial volume test flag
-  Int_t fdbg_clusters_cal;      // Shower debug flag
+  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_decoded_cal;      // Shower debug flags
+  Int_t fdbg_sparsified_cal;
+  Int_t fdbg_clusters_cal;
+  Int_t fdbg_tracks_cal;
 
   THcShowerPlane** fPlanes;     // [fNLayers] Shower Plane objects
 
@@ -122,7 +152,9 @@ public:
 
   void Setup(const char* name, const char* description);
 
-  void MatchCluster(THaTrack*);
+  Int_t MatchCluster(THaTrack*, THcShowerClusterList*, Double_t&, Double_t&);
+
+  friend class THcShowerPlane;
 
   ClassDef(THcShower,0)         // Generic class
 };

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];
@@ -214,8 +240,26 @@ Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks )
  
   //  HitCount();
 
+  /*
+    if (THcShower::fdbg_tracks_cal)
   cout << "THcShowerPlane::CoarseProcess called ---------------------" << endl;
 
+  Int_t Ntracks = tracks.GetLast()+1;   // Number of reconstructed tracks
+
+  if (THcShower::fdbg_tracks_cal)
+  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;
 }
 
@@ -231,6 +275,12 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
   // Assumes that the hit list is sorted by layer, so we stop when the
   // plane doesn't agree and return the index for the next hit.
 
+  THcShower* fParent;
+  fParent = (THcShower*) GetParent();
+
+  if (fParent->fdbg_decoded_cal)
+    cout << "THcShowerPlane::ProcessHits called ----" << endl;
+
   Int_t nPosADCHits=0;
   Int_t nNegADCHits=0;
   fPosADCHits->Clear();
@@ -251,12 +301,19 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
   Int_t nrawhits = rawhits->GetLast()+1;
 
   Int_t ihit = nexthit;
-  //cout << "nrawhits =  " << nrawhits << endl;
-  //cout << "nexthit =  " << nexthit << endl;
+
+  if (fParent->fdbg_decoded_cal)
+ {
+    cout << "   nrawhits =  " << nrawhits << endl;
+    cout << "   nexthit =  " << nexthit << endl;
+  }
+
   while(ihit < nrawhits) {
     THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
 
-    //cout << "fplane =  " << hit->fPlane << " Num = " << fLayerNum << endl;
+    if (fParent->fdbg_decoded_cal)
+      cout << "   fplane =  " << hit->fPlane << " Num = " << fLayerNum << endl;
+
     if(hit->fPlane > fLayerNum) {
       break;
     }
@@ -268,13 +325,11 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
     fA_Pos_p[hit->fCounter-1] = hit->fADC_pos - fPosPed[hit->fCounter -1];
     fA_Neg_p[hit->fCounter-1] = hit->fADC_neg - fNegPed[hit->fCounter -1];
 
-    THcShower* fParent;
-    fParent = (THcShower*) GetParent();
-
     Double_t thresh_pos = fPosThresh[hit->fCounter -1];
     if(hit->fADC_pos >  thresh_pos) {
 
-      THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+      THcSignalHit *sighit =
+	(THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
       sighit->Set(hit->fCounter, hit->fADC_pos);
 
       fEpos[hit->fCounter-1] += fA_Pos_p[hit->fCounter-1]*
@@ -284,7 +339,8 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
     Double_t thresh_neg = fNegThresh[hit->fCounter -1];
     if(hit->fADC_neg >  thresh_neg) {
 
-      THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+      THcSignalHit *sighit = 
+	(THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
       sighit->Set(hit->fCounter, hit->fADC_neg);
 
       fEneg[hit->fCounter-1] += fA_Neg_p[hit->fCounter-1]*
@@ -296,6 +352,10 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
 
     ihit++;
   }
+
+  if (fParent->fdbg_decoded_cal) 
+    cout << "THcShowerPlane::ProcessHits return ----" << endl;
+
   return(ihit);
 }
 
@@ -306,7 +366,10 @@ Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
   // arrays for calculating pedestals.
 
   Int_t nrawhits = rawhits->GetLast()+1;
-  //  cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " " << nexthit << "/" << nrawhits << endl;
+
+  //  cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " " 
+  //  << nexthit << "/" << nrawhits << endl;
+
   Int_t ihit = nexthit;
   while(ihit < nrawhits) {
     THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
@@ -362,8 +425,8 @@ void THcShowerPlane::CalculatePedestals( )
 		      - fNegPed[i]*fNegPed[i]);
     fNegThresh[i] = fNegPed[i] + TMath::Min(50., TMath::Max(10., 3.*fNegSig[i]));
 
-    cout << "Ped&Thr: " << fPosPed[i] << " " << fPosThresh[i] << " " <<
-      fNegPed[i] << " " << fNegThresh[i] << " " << i+1 << endl;
+    //    cout << "Ped&Thr: " << fPosPed[i] << " " << fPosThresh[i] << " " <<
+    //      fNegPed[i] << " " << fNegThresh[i] << " " << i+1 << endl;
   }
   //  cout << " " << endl;
   

src/THcShowerPlane.h

@@ -46,6 +46,10 @@ class THcShowerPlane : public THaSubDetector {
 
   TClonesArray* fParentHitList;
 
+  TVector3 GetOrigin() {
+    return fOrigin;
+  }
+
   Double_t GetEplane() {
     return fEplane;
   };