Added more checks on clustering

examples/compclustepr.C

+void compclustepr(Int_t run)
+{
+  TFile* f = new TFile(Form("hodtest_%d.root",run));
+  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+  TH1F* h = eprmax;
+
+  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 = h415;   //edep
+  }
+
+  TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster Eprsh", 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", "Eprsh differences", 1000, 667); 
+
+  TH1F* dif = h->Clone();
+
+  dif->Add(h,h1,1.,-1.);
+
+  dif->SetTitle("Eprsh Difference");
+  dif->SetFillColor(kRed);
+  dif->SetLineColor(kRed);
+  dif->SetLineWidth(1);
+  dif->SetFillStyle(1111);
+  dif->Draw();
+} 

examples/compclustx.C

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

src/THcShower.cxx

@@ -173,15 +173,15 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
   BlockThick = new Double_t [fNLayers];
   fNBlocks = new Int_t [fNLayers];
   fNLayerZPos = new Double_t [fNLayers];
-  XPos = new Double_t [2*fNLayers];
+  YPos = new Double_t [2*fNLayers];
 
   for(Int_t i=0;i<fNLayers;i++) {
     DBRequest list[]={
       {Form("cal_%s_thick",fLayerNames[i]), &BlockThick[i], kDouble},
       {Form("cal_%s_nr",fLayerNames[i]), &fNBlocks[i], kInt},
       {Form("cal_%s_zpos",fLayerNames[i]), &fNLayerZPos[i], kDouble},
-      {Form("cal_%s_left",fLayerNames[i]), &XPos[2*i], kDouble},
-      {Form("cal_%s_right",fLayerNames[i]), &XPos[2*i+1], kDouble},
+      {Form("cal_%s_right",fLayerNames[i]), &YPos[2*i], kDouble},
+      {Form("cal_%s_left",fLayerNames[i]), &YPos[2*i+1], kDouble},
       {0}
     };
     gHcParms->LoadParmValues((DBRequest*)&list, prefix);
@@ -189,11 +189,11 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
 
   //Caution! Z positions (fronts) are off in hcal.param! Correct later on.
 
-  YPos = new Double_t* [fNLayers];
+  XPos = new Double_t* [fNLayers];
   for(Int_t i=0;i<fNLayers;i++) {
-    YPos[i] = new Double_t [fNBlocks[i]];
+    XPos[i] = new Double_t [fNBlocks[i]];
     DBRequest list[]={
-      {Form("cal_%s_top",fLayerNames[i]),YPos[i], kDouble, fNBlocks[i]},
+      {Form("cal_%s_top",fLayerNames[i]),XPos[i], kDouble, fNBlocks[i]},
       {0}
     };
     gHcParms->LoadParmValues((DBRequest*)&list, prefix);
@@ -204,10 +204,10 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
     cout << "    Block thickness: " << BlockThick[i] << endl;
     cout << "    NBlocks        : " << fNBlocks[i] << endl;
     cout << "    Z Position     : " << fNLayerZPos[i] << endl;
-    cout << "    X Positions    : " << XPos[2*i] << ", " << XPos[2*i+1] << endl;
-    cout << "    Y Positions    :";
+    cout << "    Y Positions    : " << YPos[2*i] << ", " << YPos[2*i+1] <<endl;
+    cout << "    X Positions    :";
     for(Int_t j=0; j<fNBlocks[i]; j++) {
-      cout << " " << YPos[i][j];
+      cout << " " << XPos[i][j];
     }
     cout << endl;
 
@@ -392,8 +392,9 @@ Int_t THcShower::DefineVariables( EMode mode )
  //   { "asum_c", "Sum of calibrated ADCs",             "fAsum_c" },
     { "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" },
+    { "eprmax",   "Preshower Energy (MeV) of largest cluster",    "fEpr" },
+    { "xmax",      "x-position (cm) of largest cluster", "fX" },
+ //   { "z",      "z-position (cm) of largest cluster", "fZ" },
     { "mult",   "Multiplicity of largest cluster",    "fMult" },
  //   { "nblk",   "Numbers of blocks in main cluster",  "fNblk" },
  //   { "eblk",   "Energies of blocks in main cluster", "fEblk" },
@@ -430,7 +431,7 @@ void THcShower::DeleteArrays()
   delete [] fNBlocks;  fNBlocks = NULL;
   delete [] fNLayerZPos;  fNLayerZPos = NULL;
   delete [] XPos;  XPos = NULL;
-  delete [] YPos;  YPos = NULL;
+  delete [] ZPos;  ZPos = NULL;
   //delete [] fSpacing;  fSpacing = NULL;
   //delete [] fCenter;   fCenter = NULL; // This 2D. What is correct way to delete?
 }
@@ -452,6 +453,8 @@ void THcShower::Clear(Option_t* opt)
   fNclust = 0;
   fMult = 0;
   fE = 0.;
+  fEpr = 0.;
+  fX = -75.;   //out of acceptance
 }
 
 //_____________________________________________________________________________
@@ -541,9 +544,9 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
       //
       //      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 x = 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);
+      //      	THcShowerHit* hit = new THcShowerHit(i,j,x,z,Edep);
 
       //ENGINE way.
       //
@@ -554,13 +557,13 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
 	  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 x = XPos[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);
+	THcShowerHit* hit = new THcShowerHit(i,j,x,z,Edep);
 
 	HitList.push_back(hit);
 
-	//cout << "Hit: Edep = " << Edep << " Y = " << y << " Z = " << z <<
+	//cout << "Hit: Edep = " << Edep << " X = " << x << " Z = " << z <<
 	//	  " Block " << i << " Layer " << j << endl;
       };
 
@@ -591,7 +594,7 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
     cout << "Cluster #" << i 
          <<":  E=" << (*cluster).clE() 
          << "  Epr=" << (*cluster).clEpr()
-         << "  Y=" << (*cluster).clY()
+         << "  X=" << (*cluster).clX()
          << "  Z=" << (*cluster).clZ()
          << "  size=" << (*cluster).clSize()
          << endl;
@@ -624,13 +627,17 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
     }
 
     fE = (*MaxCluster).clE();
+    fEpr = (*MaxCluster).clEpr();
+    fX = (*MaxCluster).clX();
   }
 
+  cout << fEpr << " " << fE << " " << fX << " PrSh" << endl;
+
   /*
     cout << "Cluster #" << i 
          <<":  E=" << (*cluster).clE() 
          << "  Epr=" << (*cluster).clEpr()
-         << "  Y=" << (*cluster).clY()
+         << "  X=" << (*cluster).clX()
          << "  Z=" << (*cluster).clZ()
          << "  size=" << (*cluster).clSize()
          << endl;

src/THcShower.h

@@ -93,8 +93,9 @@ protected:
   Float_t fAsum_c; // Sum of calibrated ADCs
   Int_t fNclust;   // Number of clusters
   Float_t fE;      // Energy (MeV) of largest cluster
+  Float_t fEpr;    // Preshower Energy (MeV) of largest cluster
   Float_t fX;      // x-position (cm) of largest cluster
-  Float_t fY;      // y-position (cm) of largest cluster
+  Float_t fZ;      // z-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
@@ -110,8 +111,9 @@ protected:
   Double_t* BlockThick;		// Thickness of shower counter blocks, blocks
   Int_t* fNBlocks;              // Number of shower counter blocks per layer
   Int_t fNtotBlocks;            // Total number of shower counter blocks
-  Double_t** YPos;		//X,Y positions of shower counter blocks
-  Double_t* XPos;
+  Double_t** XPos;		//X,Y,Z positions of shower counter 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

src/THcShowerCluster.h

@@ -44,18 +44,18 @@ class THcShowerCluster : THcShowerHitList {
     (*(THcShowerHitList::begin()+num))->show();
   }
 
-  //Y coordinate of cluster's center of gravity.
+  //X coordinate of cluster's center of gravity.
   //
-  float clY() {
-    float y_sum=0.;
+  float clX() {
+    float x_sum=0.;
     float Etot=0.;
     for (THcShowerHitIt it=THcShowerHitList::begin();
 	 it!=THcShowerHitList::end(); it++) {
-      y_sum += (*it)->hitY() * (*it)->hitE();
+      x_sum += (*it)->hitX() * (*it)->hitE();
       Etot += (*it)->hitE();
     }
-    //    cout << "y_sum=" << y_sum << "  Etot=" << Etot << endl;
-    return y_sum/Etot;
+    //    cout << "x_sum=" << x_sum << "  Etot=" << Etot << endl;
+    return (Etot != 0. ? x_sum/Etot : -75.);
   }
 
   //Z coordinate for a cluster, calculated as a weighted by energy average.

src/THcShowerHit.h

@@ -13,22 +13,22 @@ class THcShowerHit {           //HMS calorimeter hit class
 
  private:
   unsigned int fCol, fRow; //hit colomn and row
-  float fY, fZ;            //hit Y (vert.) and Z (along spect.axis) coordinates
+  float fX, fZ;            //hit X (vert.) and Z (along spect.axis) coordinates
   float fE;                 //hit energy deposition
 
  public:
 
   THcShowerHit() {             //default constructor
   fCol=fRow=0;
-  fY=fZ=0.;
+  fX=fZ=0.;
   fE=0.;
   }
 
-  THcShowerHit(unsigned int hRow, unsigned int hCol, float hY, float hZ,
+  THcShowerHit(unsigned int hRow, unsigned int hCol, float hX, float hZ,
 	       float hE) {
     fRow=hRow;
     fCol=hCol;
-    fY=hY;
+    fX=hX;
     fZ=hZ;
     fE=hE;
   }
@@ -45,8 +45,8 @@ class THcShowerHit {           //HMS calorimeter hit class
     return fRow;
   }
 
-  float hitY() {
-    return fY;
+  float hitX() {
+    return fX;
   }
 
   float hitZ() {
@@ -67,7 +67,7 @@ class THcShowerHit {           //HMS calorimeter hit class
   //
   void show() {
     cout << "row=" << fRow << "  column=" << fCol << 
-      "  y=" << fY << "  z=" << fZ << "  E=" << fE << endl;
+      "  x=" << fX << "  z=" << fZ << "  E=" << fE << endl;
   }
 
 };