diff --git a/examples/compclustedep.C b/examples/compclustedep.C
new file mode 100644
index 0000000000000000000000000000000000000000..c251ae614d23f1440cbde67db14b2fb3e9192027
--- /dev/null
+++ b/examples/compclustedep.C
@@ -0,0 +1,57 @@
+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();
+}
diff --git a/examples/compclustsize.C b/examples/compclustsize.C
new file mode 100644
index 0000000000000000000000000000000000000000..4707a0d3a2d17a11a1b71a184628f541d35187e0
--- /dev/null
+++ b/examples/compclustsize.C
@@ -0,0 +1,57 @@
+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();
+}
diff --git a/examples/compedeps.C b/examples/compedeps.C
index 5998740ed353b461180aeef97d72e59bb0659aa2..dcb8f0b3d7beeda77dd232f0917a4e8b51808221 100644
--- a/examples/compedeps.C
+++ b/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);
diff --git a/examples/compnclusts.C b/examples/compnclusts.C
new file mode 100644
index 0000000000000000000000000000000000000000..2f16142e2e3bf1d91cc3f007933a70a104bebe60
--- /dev/null
+++ b/examples/compnclusts.C
@@ -0,0 +1,59 @@
+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();
+
+}
diff --git a/examples/compnhits.C b/examples/compnhits.C
new file mode 100644
index 0000000000000000000000000000000000000000..b30f829f4da4d6d4e7a68b1914e5463e0f755543
--- /dev/null
+++ b/examples/compnhits.C
@@ -0,0 +1,57 @@
+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();
+}
diff --git a/examples/comphits.C b/examples/comprawhits.C
similarity index 97%
rename from examples/comphits.C
rename to examples/comprawhits.C
index 6489c463a1d83f40e113dcc6f9173b58bdcab562..6b8890f10760a786c6e8f6abaec06c8f4ef97037 100644
--- a/examples/comphits.C
+++ b/examples/comprawhits.C
@@ -1,4 +1,4 @@
-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));
diff --git a/examples/hodtest_mkj.C b/examples/hodtest_mkj.C
new file mode 100644
index 0000000000000000000000000000000000000000..64c47cd3b12581be54588eb4ad5268f50d10ed13
--- /dev/null
+++ b/examples/hodtest_mkj.C
@@ -0,0 +1,82 @@
+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
+}
diff --git a/src/THcRawShowerHit.cxx b/src/THcRawShowerHit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..75c3a322642ddda42996eced2697a45ed91a6b34
--- /dev/null
+++ b/src/THcRawShowerHit.cxx
@@ -0,0 +1,75 @@
+///////////////////////////////////////////////////////////////////////////////
+// //
+// 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)
+
+
diff --git a/src/THcRawShowerHit.h b/src/THcRawShowerHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..1671ddcbfcc7645315c2e99e7440c1b40927d067
--- /dev/null
+++ b/src/THcRawShowerHit.h
@@ -0,0 +1,36 @@
+#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
+
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index db9b4ab11e929a4aeba7d97e7ff30bd6e7846171..b49a109136009ef94aa54603477bf3e7cf2588a3 100644
--- a/src/THcShower.cxx
+++ b/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;
diff --git a/src/THcShower.h b/src/THcShower.h
index dbab5b596f715f803dd0bdd12509db2e94377d5a..24f4c5601147cb1a26a37d3a08c61b9162093de0 100644
--- a/src/THcShower.h
+++ b/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
diff --git a/src/THcShowerCluster.h b/src/THcShowerCluster.h
new file mode 100644
index 0000000000000000000000000000000000000000..513e71c9860d0d47b8e3b5abb02099c05db0d40e
--- /dev/null
+++ b/src/THcShowerCluster.h
@@ -0,0 +1,200 @@
+#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
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index a6d642f3b8aa2506c3baa8fd85422826049c3d2d..39d94b259619039839a72aa7d59ca6818bc9de9f 100644
--- a/src/THcShowerPlane.cxx
+++ b/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;
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 149cb2014e7e6466d3618ef63a009270c799be7f..c38218b475cb98421605fcb60948d385bef25bda 100644
--- a/src/THcShowerPlane.h
+++ b/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