Added checks for clustering

examples/compclustedep.C

+void compclustedep(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = emax;
+
+  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 = h414;   //edep
+  }
+
+  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.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", "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/compclustsize.C

+void compclustsize(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = mult;
+
+  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 = h413;   //mult
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster size", 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.095*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.045*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Size differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Size Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/compedeps.C

@@ -20,10 +20,10 @@ void compedeps(Int_t run)
     h1[3] = h686;   //C
     break;
   default :
-    h1[0] = h627;   //A
-    h1[1] = h628;   //B
-    h1[2] = h629;   //C
-    h1[3] = h630;   //C
+    h1[0] = h628;   //A
+    h1[1] = h629;   //B
+    h1[2] = h630;   //C
+    h1[3] = h631;   //C
   }
 
   TCanvas *c1 = new TCanvas("c1", "Shower raw Edeps", 1000, 667); 

examples/compnclusts.C

+void compnclusts(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = nclust;
+
+  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 = h412;   //hnclusters
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Cluster Map", 1000, 667); 
+
+  gPad->SetLogy();
+
+  h1->SetFillColor(kGreen);
+  h1->SetLineColor(kGreen);
+  h1->SetFillStyle(1111);
+  h1->Draw();
+
+  h->SetFillColor(kBlue);
+  h->SetLineWidth(2);
+  h->SetFillStyle(0);
+  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", "Cluster differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+
+} 

examples/compnhits.C

+void compnhits(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = nhits;
+
+  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 = h411;   //hnhits
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Sparsified Hit Map", 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.095*maxy,"Engine");
+  l.SetTextColor(kBlue);
+  l.DrawLatex(xt,0.045*maxy,"hcana");
+    
+  // Difference between the histograms.
+
+  TCanvas *c2 = new TCanvas("c2", "Hit differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/comphits.C → examples/comprawhits.C

-void comphits(Int_t run)
+void comprawhits(Int_t run)
 {
   //  TFile* f = new TFile("hodtest.root");
   TFile* f = new TFile(Form("hodtest_%d.root",run));

examples/hodtest_mkj.C

+void hodtest_mkj(Int_t RunNumber=50017, Int_t MaxEventToReplay=5000) {
+
+  //
+  //  Steering script to test hodoscope decoding
+  //
+  
+  if (RunNumber == 50017) {
+    char* RunFileNamePattern="daq04_%d.log.0";
+  } else {
+    char* RunFileNamePattern="/cache/mss/hallc/daq04/raw/daq04_%d.log.0";
+  }
+  gHcParms->Define("gen_run_number", "Run Number", RunNumber);
+  gHcParms->AddString("g_ctp_database_filename", "DBASE/test.database");
+  
+  gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
+
+  // g_ctp_parm_filename and g_decode_map_filename should now be defined
+
+
+  gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
+
+  // Constants not in ENGINE PARAM files that we want to be
+  // configurable
+  gHcParms->Load("PARAM/hcana.param");
+
+  // Generate db_cratemap to correspond to map file contents
+  char command[100];
+  sprintf(command,"./make_cratemap.pl < %s > db_cratemap.dat",gHcParms->GetString("g_decode_map_filename"));
+  system(command);
+
+  // Load the Hall C style detector map
+  gHcDetectorMap=new THcDetectorMap();
+  gHcDetectorMap->Load(gHcParms->GetString("g_decode_map_filename"));
+
+  // Set up the equipment to be analyzed.
+
+  THaApparatus* HMS = new THcHallCSpectrometer("H","HMS");
+  gHaApps->Add( HMS );
+
+  // Add hodoscope
+  HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
+  HMS->AddDetector( new THcShower("cal", "Shower" ));
+  HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
+  THcAerogel* aerogel = new THcAerogel("aero", "Aerogel Cerenkov" );
+  HMS->AddDetector( aerogel );
+
+  // Set up the analyzer - we use the standard one,
+  // but this could be an experiment-specific one as well.
+  // The Analyzer controls the reading of the data, executes
+  // tests/cuts, loops over Acpparatus's and PhysicsModules,
+  // and executes the output routines.
+  THaAnalyzer* analyzer = new THcAnalyzer;
+  
+
+  // A simple event class to be output to the resulting tree.
+  // Creating your own descendant of THaEvent is one way of
+  // defining and controlling the output.
+  THaEvent* event = new THaEvent;
+  
+  // Define the run(s) that we want to analyze.
+  // We just set up one, but this could be many.
+  char RunFileName[100];
+  sprintf(RunFileName,RunFileNamePattern,RunNumber);
+  THaRun* run = new THaRun(RunFileName);
+
+  // Eventually need to learn to skip over, or properly analyze
+  // the pedestal events
+  run->SetEventRange(1,MaxEventToReplay);//  Physics Event number, does not
+                                // include scaler or control events
+
+  // Define the analysis parameters
+  analyzer->SetCountMode( 0 );   //Mark's modification
+  analyzer->SetEvent( event );
+  analyzer->SetOutFile(Form("hodtest_%05d.root",RunNumber));
+  analyzer->SetOdefFile(Form("output_%d.def",RunNumber));
+  analyzer->SetCutFile("hodtest_cuts_mkj.def");        // optional
+  
+  // File to record cuts accounting information
+  //  analyzer->SetSummaryFile("summary_example.log"); // optional
+  
+  analyzer->Process(run);     // start the actual analysis
+}

src/THcRawShowerHit.cxx

+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// THcRawShowerHit                                                           //
+//                                                                           //
+// Class representing a single raw hit for a hodoscope paddle                //
+//                                                                           //
+// Contains plane, counter and pos/neg adc and tdc values                    //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcRawShowerHit.h"
+
+using namespace std;
+
+
+void THcRawShowerHit::SetData(Int_t signal, Int_t data) {
+  if(signal==0) {
+    fADC_pos = data;
+  } else if (signal==1) {
+    fADC_neg = data;
+  } 
+}
+
+Int_t THcRawShowerHit::GetData(Int_t signal) {
+  if(signal==0) {
+    return(fADC_pos);
+  } else if (signal==1) {
+    return(fADC_neg);
+  } 
+
+  return(-1); // Actually should throw exception
+}
+
+#if 0
+Int_t THcRawShowerHit::Compare(const TObject* obj) const
+{
+  // Compare to sort by plane and counter
+
+  const THcRawShowerHit* hit = dynamic_cast<const THcRawShowerHit*>(obj);
+
+  if(!hit) return -1;
+  Int_t p1 = fPlane;
+  Int_t p2 = hit->fPlane;
+  if(p1 < p2) return -1;
+  else if(p1 > p2) return 1;
+  else {
+    Int_t c1 = fCounter;
+    Int_t c2 = hit->fCounter;
+    if(c1 < c2) return -1;
+    else if (c1 == c2) return 0;
+    else return 1;
+  }
+}
+#endif
+//_____________________________________________________________________________
+THcRawShowerHit& THcRawShowerHit::operator=( const THcRawShowerHit& rhs )
+{
+  // Assignment operator.
+
+  THcRawHit::operator=(rhs);
+  if ( this != &rhs ) {
+    fPlane = rhs.fPlane;
+    fCounter = rhs.fCounter;
+    fADC_pos = rhs.fADC_pos;
+    fADC_neg = rhs.fADC_neg;
+
+  }
+  return *this;
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+ClassImp(THcRawShowerHit)
+
+ 

src/THcRawShowerHit.h

+#ifndef ROOT_THcRawShowerHit
+#define ROOT_THcRawShowerHit
+
+#include "THcRawHit.h"
+
+class THcRawShowerHit : public THcRawHit {
+
+ public:
+
+ THcRawShowerHit(Int_t plane=0, Int_t counter=0) : THcRawHit(plane, counter), 
+    fADC_pos(-1), fADC_neg(-1){
+  }
+  THcRawShowerHit& operator=( const THcRawShowerHit& );
+  virtual ~THcRawShowerHit() {}
+
+  virtual void Clear( Option_t* opt="" )
+    { fADC_pos = -1; fADC_neg = -1; }
+
+  void SetData(Int_t signal, Int_t data);
+  Int_t GetData(Int_t signal);
+
+  //  virtual Bool_t  IsSortable () const {return kTRUE; }
+  //  virtual Int_t   Compare(const TObject* obj) const;
+
+  Int_t fADC_pos;
+  Int_t fADC_neg;
+
+ protected:
+
+ private:
+
+  ClassDef(THcRawShowerHit, 0);	// Shower hit class
+};  
+
+#endif
+ 

src/THcShower.cxx

@@ -383,25 +383,26 @@ Int_t THcShower::DefineVariables( EMode mode )
 
   // Register variables in global list
 
- //  RVarDef vars[] = {
- //   { "nhit",   "Number of hits",                     "fNhits" },
+  RVarDef vars[] = {
+    { "nhits",  "Number of hits",                     "fNhits" },
  //   { "a",      "Raw ADC amplitudes",                 "fA" },
  //   { "a_p",    "Ped-subtracted ADC amplitudes",      "fA_p" },
  //   { "a_c",    "Calibrated ADC amplitudes",          "fA_c" },
  //   { "asum_p", "Sum of ped-subtracted ADCs",         "fAsum_p" },
  //   { "asum_c", "Sum of calibrated ADCs",             "fAsum_c" },
- //   { "nclust", "Number of clusters",                 "fNclust" },
- //   { "e",      "Energy (MeV) of largest cluster",    "fE" },
+    { "nclust", "Number of clusters",                 "fNclust" },
+    { "emax",   "Energy (MeV) of largest cluster",    "fE" },
  //   { "x",      "x-position (cm) of largest cluster", "fX" },
  //   { "y",      "y-position (cm) of largest cluster", "fY" },
- //   { "mult",   "Multiplicity of largest cluster",    "fMult" },
+    { "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" },
- //   { 0 }
- // };
- //  return DefineVarsFromList( vars, mode );
+    { 0 }
+  };
+  return DefineVarsFromList( vars, mode );
+
   return kOK;
 }
 
@@ -438,16 +439,27 @@ void THcShower::DeleteArrays()
 inline 
 void THcShower::Clear(Option_t* opt)
 {
+
 //   Reset per-event data.
+
   for(Int_t ip=0;ip<fNLayers;ip++) {
     fPlanes[ip]->Clear(opt);
   }
+
  // fTrackProj->Clear();
+
+  fNhits = 0;
+  fNclust = 0;
+  fMult = 0;
+  fE = 0.;
 }
 
 //_____________________________________________________________________________
 Int_t THcShower::Decode( const THaEvData& evdata )
 {
+
+  Clear();
+
   // Get the Hall C style hitlist (fRawHitList) for this event
   Int_t nhits = THcHitList::DecodeToHitList(evdata);
 
@@ -525,10 +537,25 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
 
     for (Int_t i=0; i<fNBlocks[j]; i++) {
 
-      Float_t Edep = fPlanes[j]->GetEmean(i);
-      if (Edep > 0.) {                                    //hit
-	Float_t y = YPos[j][i] + BlockThick[j]/2.;        //top + thick/2
-	Float_t z = fNLayerZPos[j] + BlockThick[j]/2.;    //front + thick/2
+      //May be should be done this way.
+      //
+      //      Float_t Edep = fPlanes[j]->GetEmean(i);
+      //      if (Edep > 0.) {                                    //hit
+      //	Float_t y = YPos[j][i] + BlockThick[j]/2.;        //top + thick/2
+      //	Float_t z = fNLayerZPos[j] + BlockThick[j]/2.;    //front + thick/2
+      //      	THcShowerHit* hit = new THcShowerHit(i,j,y,z,Edep);
+
+      //ENGINE way.
+      //
+      //      if (fPlanes[j]->GetApos(i)> fPlanes[j]->GetPosThr(i) ||
+      //	  fPlanes[j]->GetAneg(i)> fPlanes[j]->GetNegThr(i)) {    //hit
+      if (fPlanes[j]->GetApos(i) - fPlanes[j]->GetPosPed(i) >
+	  fPlanes[j]->GetPosThr(i) - fPlanes[j]->GetPosPed(i) ||
+	  fPlanes[j]->GetAneg(i) - fPlanes[j]->GetNegPed(i) >
+	  fPlanes[j]->GetNegThr(i) - fPlanes[j]->GetNegPed(i)) {    //hit
+	Float_t Edep = fPlanes[j]->GetEmean(i);
+      	Float_t y = YPos[j][i] + BlockThick[j]/2.;        //top + thick/2
+      	Float_t z = fNLayerZPos[j] + BlockThick[j]/2.;    //front + thick/2
 	THcShowerHit* hit = new THcShowerHit(i,j,y,z,Edep);
 
 	HitList.push_back(hit);
@@ -542,8 +569,9 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
 
   //Print out hits before clustering.
   //
-  cout << "Total hits:     " << HitList.size() << endl;
-  for (unsigned int i=0; i!=HitList.size(); i++) {
+  fNhits = HitList.size();
+  cout << "Total hits:     " << fNhits << endl;
+  for (unsigned int i=0; i!=fNhits; i++) {
     cout << "unclustered hit " << i << ": ";
     (*(HitList.begin()+i))->show();
   }
@@ -553,9 +581,10 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
 
   //Print out the cluster list.
   //
-  cout << "Cluster_list size: " << (*ClusterList).NbClusters() << endl;
+  fNclust = (*ClusterList).NbClusters();
+  cout << "Cluster_list size: " << fNclust << endl;
 
-  for (unsigned int i=0; i!=(*ClusterList).NbClusters(); i++) {
+  for (unsigned int i=0; i!=fNclust; i++) {
 
     THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
 
@@ -573,6 +602,39 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
     }
 
   }
+
+  
+  //The largest cluster.
+  //
+
+  if (fNclust != 0) {
+
+    THcShowerCluster* MaxCluster;
+
+    for (unsigned int i=0; i!=fNclust; i++) {
+
+      THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+
+      Int_t size = (*cluster).clSize();
+
+      if (fMult < size) {
+	fMult = size;
+	MaxCluster = cluster;
+      }
+    }
+
+    fE = (*MaxCluster).clE();
+  }
+
+  /*
+    cout << "Cluster #" << i 
+         <<":  E=" << (*cluster).clE() 
+         << "  Epr=" << (*cluster).clEpr()
+         << "  Y=" << (*cluster).clY()
+         << "  Z=" << (*cluster).clZ()
+         << "  size=" << (*cluster).clSize()
+         << endl;
+  */
    
   cout << "THcShower::CoarseProcess return ---------------------------" <<endl;
   return 0;

src/THcShower.h

@@ -63,6 +63,10 @@ public:
     return fShMinPeds;
   }
 
+  Int_t GetMult() {
+    return fMult;
+  }
+
   THcShower();  // for ROOT I/O
 
 protected:
@@ -81,6 +85,21 @@ protected:
 
   // Per-event data
 
+  Int_t fNhits;    // Number of hits
+  Float_t** fA;    // Raw ADC amplitudes
+  Float_t** fA_p;  // Ped-subtracted ADC amplitudes
+  Float_t** fA_c;  // Calibrated ADC amplitudes
+  Float_t fAsum_p; // Sum of ped-subtracted ADCs
+  Float_t fAsum_c; // Sum of calibrated ADCs
+  Int_t fNclust;   // Number of clusters
+  Float_t fE;      // Energy (MeV) of largest cluster
+  Float_t fX;      // x-position (cm) of largest cluster
+  Float_t fY;      // y-position (cm) of largest cluster
+  Int_t fMult;     // Multiplicity of largest cluster
+  Int_t fNblk;     // Number of blocks in main cluster
+  Float_t* fEblk;  // Energies of blocks in main cluster
+  Float_t fTRX;    // track x-position in det plane"
+  Float_t fTRY;    // track y-position in det plane",      
 
   // Potential Hall C parameters.  Mostly here for demonstration
 

src/THcShowerCluster.h

+#ifndef ROOT_THcShowerCluster
+#define ROOT_THcShowerCluster
+
+//HMS calorimeter hit cluster, version 3.
+
+#include "THcShowerHit.h"
+
+//HMS calorimeter hit cluster
+//
+class THcShowerCluster : THcShowerHitList {
+
+ public:
+
+  THcShowerCluster() {
+    //    cout << "Dummy THcShowerCluster object created" << endl;
+  }
+
+  ~THcShowerCluster() {
+    for (THcShowerHitIt i = THcShowerHitList::begin();
+    	 i != THcShowerHitList::end(); i++) {
+      delete *i;
+      *i = 0;
+    }
+    //    purge(THcShowerHitList::);
+    //    THcShowerHitList::clear();
+    //    cout << "THcShowerCluster object destructed" << endl;
+  }
+
+  // Add a hit to the cluster hit list
+  //
+  void grow(THcShowerHit* hit) {
+    THcShowerHitList::push_back(hit);
+  }
+
+  //Pointer to the hit #i in the cluster hit list
+  //
+  THcShowerHit* ClusteredHit(unsigned int i) {
+    return * (THcShowerHitList::begin()+i);
+  }
+
+  //Print out a hit in the cluster
+  //
+  void showHit(unsigned int num) {
+    (*(THcShowerHitList::begin()+num))->show();
+  }
+
+  //Y coordinate of cluster's center of gravity.
+  //
+  float clY() {
+    float y_sum=0.;
+    float Etot=0.;
+    for (THcShowerHitIt it=THcShowerHitList::begin();
+	 it!=THcShowerHitList::end(); it++) {
+      y_sum += (*it)->hitY() * (*it)->hitE();
+      Etot += (*it)->hitE();
+    }
+    //    cout << "y_sum=" << y_sum << "  Etot=" << Etot << endl;
+    return y_sum/Etot;
+  }
+
+  //Z coordinate for a cluster, calculated as a weighted by energy average.
+  //
+  float clZ() {
+    float z_sum=0.;
+    float Etot=0.;
+    for (THcShowerHitIt it=THcShowerHitList::begin();
+	 it!=THcShowerHitList::end(); it++) {
+      z_sum += (*it)->hitZ() * (*it)->hitE();
+      Etot += (*it)->hitE();
+    }
+    //    cout << "z_sum=" << z_sum << "  Etot=" << Etot << endl;
+    return z_sum/Etot;
+  }
+
+  //Energy depostion in a cluster
+  //
+  float clE() {
+    //    cout << "In ECl:" << endl;
+    float Etot=0.;
+    for (THcShowerHitIt it=THcShowerHitList::begin();
+	 it!=THcShowerHitList::end(); it++) {
+      Etot += (*it)->hitE();
+    }
+    return Etot;
+  }
+
+  //Energy deposition in the Preshower (1st layer) for a cluster
+  //
+  float clEpr() {
+    float Epr=0.;
+    for (THcShowerHitIt it=THcShowerHitList::begin();
+	 it!=THcShowerHitList::end(); it++) {
+      if ((*it)->hitColumn() == 0) Epr += (*it)->hitE();
+    }
+    return Epr;
+  }
+
+  //Cluster size.
+  //
+  unsigned int clSize() {
+    return THcShowerHitList::size();
+  }
+
+};
+
+//-----------------------------------------------------------------------------
+
+//Alias for cluster container and for its iterator
+//
+typedef vector<THcShowerCluster*> THcShClusterList;
+typedef THcShClusterList::iterator THcShClusterIt;
+
+//List of clusters
+
+class THcShowerClusterList : private THcShClusterList {
+
+ public:
+
+  THcShowerClusterList() {
+    //    cout << "Dummy THcShowerClusterList object created" << endl;
+  }
+
+  ~THcShowerClusterList() {
+    for (THcShClusterIt i = THcShClusterList::begin();
+    	 i != THcShClusterList::end(); i++) {
+      delete *i;
+      *i = 0;
+    }
+    //    purge(THcShClusterList);
+    //    THcShClusterList::clear();
+    //    cout << "THcShowerClusterList object destroyed" << endl;
+  }
+
+  //Put a cluster in the cluster list
+  //
+  void grow(THcShowerCluster* cluster) {
+    THcShClusterList::push_back(cluster);
+  }
+
+  //Pointer to the cluster #i in the cluster list
+  //
+  THcShowerCluster* ListedCluster(unsigned int i) {
+    return *(THcShClusterList::begin()+i);
+  }
+
+  //Cluster list size.
+  //
+  unsigned int NbClusters() {
+    return THcShClusterList::size();
+  }
+
+//_____________________________________________________________________________
+
+void ClusterHits(THcShowerHitList HitList) {
+
+//Cluster hits from the HitList. The resultant hit clusters are saved
+//in the ClusterList.
+
+  while (HitList.size() != 0) {
+
+    THcShowerCluster* cluster = new THcShowerCluster;
+
+    (*cluster).grow(*(HitList.end()-1)); //move the last hit from the hit list
+    HitList.erase(HitList.end()-1);      //into the 1st cluster
+    bool clustered = true;
+
+    while (clustered) {                   //while a hit is clustered
+
+      clustered = false;
+
+      for (THcShowerHitIt i=HitList.begin(); i!=HitList.end(); i++) {
+
+	for (unsigned int k=0; k!=(*cluster).clSize(); k++) {
+
+	  if ((**i).isNeighbour((*cluster).ClusteredHit(k))) {
+
+	    (*cluster).grow(*i);          //If hit i is neighbouring a hit
+	    HitList.erase(i);             //in the cluster, then move it
+	                                  //into cluster.
+	    clustered = true;
+	  }
+
+	  if (clustered) break;
+	}                                 //k
+
+	if (clustered) break;
+      }                                   //i
+
+    }                                     //while clustered
+
+    //    (*ClusterList).grow(cluster);  //Put the cluster in the cluster list
+    grow(cluster);        //Put the cluster in the cluster list
+
+  }                                       //while hit_list not exhausted
+
+}
+
+};
+
+#endif

src/THcShowerPlane.cxx

@@ -362,7 +362,8 @@ void THcShowerPlane::CalculatePedestals( )
     //  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 << "Ped&Thr: " << fPosPed[i] << " " << fPosThresh[i] << " " <<
+      fNegPed[i] << " " << fNegThresh[i] << " " << i+1 << endl;
   }
   //  cout << " " << endl;
   

src/THcShowerPlane.h

@@ -54,6 +54,38 @@ class THcShowerPlane : public THaSubDetector {
     return fEmean[i];
   };
 
+  Float_t GetAposP(Int_t i) {
+    return fA_Pos_p[i];
+  };
+
+  Float_t GetAnegP(Int_t i) {
+    return fA_Neg_p[i];
+  };
+
+  Float_t GetApos(Int_t i) {
+    return fA_Pos[i];
+  };
+
+  Float_t GetAneg(Int_t i) {
+    return fA_Neg[i];
+  };
+
+  Float_t GetPosThr(Int_t i) {
+    return fPosThresh[i];
+  };
+
+  Float_t GetNegThr(Int_t i) {
+    return fNegThresh[i];
+  };
+
+  Float_t GetPosPed(Int_t i) {
+    return fPosPed[i];
+  };
+
+  Float_t GetNegPed(Int_t i) {
+    return fNegPed[i];
+  };
+
  protected:
 
   Float_t*   fA_Pos;         // [fNelem] Array of ADC amplitudes of blocks