diff --git a/CALIBRATION/shms_cal_calib/THcPShHit.h b/CALIBRATION/shms_cal_calib/THcPShHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..06eb6067b3296317a6c86c05082e77bfca39b69a
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShHit.h
@@ -0,0 +1,59 @@
+#include <iostream>
+
+// SHMS calorimeter hit class for calibration.
+
+class THcPShHit {
+
+  Double_t ADC;        // pedestal subtracted ADC signal.
+  Double_t Edep;       // Energy deposition.
+  UInt_t BlkNumber;    // Block number.
+
+ public:
+
+  THcPShHit();
+  THcPShHit(Double_t adc, UInt_t blk_number);
+  ~THcPShHit();
+
+  void SetADC(Double_t sig) {ADC = sig;}
+
+  void SetEdep(Double_t e) {Edep = e;}
+
+  void SetBlkNumber(UInt_t n) {BlkNumber = n;}
+
+  Double_t GetADC() {return ADC;}
+
+  Double_t GetEdep() {return Edep;}
+
+  UInt_t GetBlkNumber() {return BlkNumber;}
+
+  void Print(ostream & ostrm);
+};
+
+//------------------------------------------------------------------------------
+
+THcPShHit::THcPShHit() {
+  ADC = -99999.;
+  Edep = -99999.;
+  BlkNumber = 99999;
+};
+
+THcPShHit::THcPShHit(Double_t adc, UInt_t blk_number) {
+  ADC = adc;
+  Edep = 0.;
+  BlkNumber = blk_number;
+};
+
+THcPShHit::~THcPShHit() { };
+
+//------------------------------------------------------------------------------
+
+void THcPShHit::Print(ostream & ostrm) {
+
+  // Output hit data through the stream ostrm.
+
+  ostrm << ADC << " " << Edep << " " << BlkNumber << endl;
+};
+
+//------------------------------------------------------------------------------
+
+struct pmt_hit {Double_t signal; UInt_t channel;};
diff --git a/CALIBRATION/shms_cal_calib/THcPShTrack.h b/CALIBRATION/shms_cal_calib/THcPShTrack.h
new file mode 100644
index 0000000000000000000000000000000000000000..35c63724251449f950e84afa2d786579ea16f272
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShTrack.h
@@ -0,0 +1,254 @@
+#include "THcPShHit.h"
+#include "TMath.h"
+
+#include <vector>
+#include <iterator>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+// Track class for the SHMS calorimeter calibration.
+// Comprises the spectrometer track parameters and calorimeter hits.
+//
+
+// Container (collection) of hits and its iterator.
+//
+typedef vector<THcPShHit*> THcPShHitList;
+typedef THcPShHitList::iterator THcPShHitIt;
+
+class THcPShTrack {
+
+  Double_t P;   // track momentum
+  Double_t Dp;  // track momentum deviation, %/
+  Double_t X;   // at the Preshower face
+  Double_t Xp;  // slope
+  Double_t Y;   // at the Preshower face
+  Double_t Yp;  // slope
+
+  THcPShHitList Hits;
+
+ public:
+
+  THcPShTrack();
+  THcPShTrack(Double_t p, Double_t dp, Double_t x, Double_t xp,
+	      Double_t y, Double_t yp);
+  ~THcPShTrack();
+
+  void Reset(Double_t p, Double_t dp, Double_t x, Double_t xp,
+	     Double_t y, Double_t yp);
+
+  void AddHit(Double_t adc, Double_t edep, UInt_t blk_number);
+
+  THcPShHit* GetHit(UInt_t k);
+
+  UInt_t GetNhits() {return Hits.size();};
+
+  void Print(ostream & ostrm);
+
+  void SetEs(Double_t* alpha);
+
+  Double_t Enorm();
+  Double_t EPRnorm();
+  Double_t ESHnorm();
+
+  Double_t GetP() {return P*1000.;}     //MeV
+
+  Double_t GetDp() {return Dp;}     //MeV
+
+  Double_t GetX() {return X;}
+  Double_t GetY() {return Y;}
+
+  Float_t Ycor(Double_t, UInt_t);       // coord. corection for Preshower module
+
+  // Coordinate correction constants for Preshower blocks
+  //
+  //  static const Double_t fAcor = 106.73;
+  //  static const Double_t fBcor = 2.329;
+  static constexpr Double_t fAcor = 106.73;
+  static constexpr Double_t fBcor = 2.329;
+
+  // Calorimeter geometry constants.
+  //
+  static const UInt_t fNrows_pr = 14;    //Row number for Preshower
+  static const UInt_t fNrows_sh = 16;    //Row number for Shower
+  static const UInt_t fNcols_pr =  2;    //2 columns in Preshower
+  static const UInt_t fNcols_sh = 14;    //14 columnsin Shower
+  static const UInt_t fNpmts_pr = fNrows_pr*fNcols_pr;
+  static const UInt_t fNpmts = fNpmts_pr + fNrows_sh*fNcols_sh;;
+
+};
+
+//------------------------------------------------------------------------------
+
+THcPShTrack::THcPShTrack() { };
+
+THcPShTrack::THcPShTrack(Double_t p, Double_t dp,
+		       Double_t x, Double_t xp, Double_t y, Double_t yp) {
+  P = p;
+  Dp = dp;
+  X = x;
+  Xp = xp;
+  Y = y;
+  Yp =yp;
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShTrack::Reset(Double_t p, Double_t dp,
+		       Double_t x, Double_t xp, Double_t y, Double_t yp) {
+
+  // Reset track parameters, clear hit list.
+
+  P = p;
+  Dp = dp;
+  X = x;
+  Xp = xp;
+  Y = y;
+  Yp =yp;
+  Hits.clear();
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShTrack::AddHit(Double_t adc, Double_t edep, UInt_t blk_number) {
+
+  // Add a hit to the hit list.
+
+  THcPShHit* hit = new THcPShHit(adc, blk_number);
+  hit->SetEdep(edep);
+  Hits.push_back(hit);
+};
+
+//------------------------------------------------------------------------------
+
+THcPShHit* THcPShTrack::GetHit(UInt_t k) {
+  THcPShHitIt it = Hits.begin();
+  for (UInt_t i=0; i<k; i++) it++;
+  return *it;
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShTrack::Print(ostream & ostrm) {
+
+  // Output the track parameters and hit list through the stream ostrm.
+
+  ostrm << P << " " << Dp << " " << X << " " << Xp << " " << Y << " " << Yp
+	<< " " << Hits.size() << endl;
+
+  for (THcPShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    (*iter)->Print(ostrm);
+  };
+
+};
+
+//------------------------------------------------------------------------------
+
+THcPShTrack::~THcPShTrack() {
+  for (THcPShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
+    delete *i;
+    *i = 0;
+  }
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShTrack::SetEs(Double_t* alpha) {
+
+  // Set hit energy depositions by use of calibration (gain) constants alpha.
+  
+  for (THcPShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+  
+    Double_t adc = (*iter)->GetADC();
+    UInt_t nblk = (*iter)->GetBlkNumber();
+
+    if(nblk <= fNrows_pr*fNcols_pr) {
+      //Preshower block, correct for Y coordinate
+      UInt_t ncol = 1;
+      if (nblk > fNrows_pr) ncol = 2;
+      (*iter)->SetEdep(adc*Ycor(Y,ncol)*alpha[nblk-1]);
+    }
+    else
+      //Shower block, no coordinate correction.
+      (*iter)->SetEdep(adc*alpha[nblk-1]);
+
+  };
+
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcPShTrack::Enorm() {
+
+  // Normalized to the track momentum energy depostion in the calorimeter.
+
+  Double_t sum = 0;
+
+  for (THcPShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    sum += (*iter)->GetEdep();
+  };
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcPShTrack::EPRnorm() {
+
+  // Normalized to the track momentum energy depostion in Preshower.
+
+  Double_t sum = 0;
+
+  for (THcPShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    if ((*iter)->GetBlkNumber() <= fNpmts_pr)
+      sum += (*iter)->GetEdep();
+  };
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcPShTrack::ESHnorm() {
+
+  // Normalized to the track momentum energy depostion in Shower.
+
+  Double_t sum = 0;
+
+  for (THcPShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    if ((*iter)->GetBlkNumber() > fNpmts_pr)
+      sum += (*iter)->GetEdep();
+  };
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+
+//------------------------------------------------------------------------------
+
+// Coordinate correction for Preshower modules.
+// Fit to GEANT pion data @ 5 GeV/c (Simon).
+
+Float_t THcPShTrack::Ycor(Double_t yhit, UInt_t ncol) {
+
+  Float_t cor;
+
+  // Warn if hit does not belong to Preshower.
+  //
+  if (ncol > fNcols_pr || ncol < 1)
+    cout << "*** THcPShTrack::Ycor: wrong ncol = " << ncol << " ***" << endl;
+
+  // Check if the hit coordinate matches the fired block's column.
+  //
+  //if ((yhit < 0. && ncol == 1) || (yhit > 0. && ncol == 2)) //Vardan's MC data
+  if ((yhit < 0. && ncol == 2) || (yhit > 0. && ncol == 1))   //Simon's MC data
+    cor = 1./(1. + TMath::Power(TMath::Abs(yhit)/fAcor, fBcor));
+  else
+    cor = 1.;
+
+  // Debug output.
+  //  cout << "THcShTrack::Ycor = " << cor << "  yhit = " << yhit
+  //       << "  ncol = " << ncol << endl;
+
+  return cor;
+}
diff --git a/CALIBRATION/shms_cal_calib/THcPShowerCalib.h b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
new file mode 100644
index 0000000000000000000000000000000000000000..46168526feda23d8f353f8f468b99f5d17004233
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
@@ -0,0 +1,757 @@
+#ifndef ROOT_THcPShowerCalib
+#define ROOT_THcPShowerCalib
+
+#include "THcPShTrack.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TDecompLU.h"
+#include "TMath.h"
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+
+#include "TROOT.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#define D_CALO_FP 275.    //distance from FP to the Preshower
+
+//Whole calorimeter
+#define XMIN -60.
+#define XMAX  60.
+#define YMIN -58.
+#define YMAX  58.
+
+#define DELTA_MIN -10   //SHMS nominal acceptance
+#define DELTA_MAX  22
+
+#define PR_ADC_THR 250
+#define SH_ADC_THR 250
+
+//#define MAX_TRACKS 10
+
+using namespace std;
+
+//
+// SHMS Calorimeter calibration class.
+//
+
+class THcPShowerCalib {
+
+ public:
+
+  THcPShowerCalib(string);
+  THcPShowerCalib();
+  ~THcPShowerCalib();
+
+  void Init();
+  bool ReadShRawTrack(THcPShTrack &trk, UInt_t ientry);
+  void CalcThresholds();
+  void ComposeVMs();
+  void SolveAlphas();
+  void FillHEcal();
+  void SaveAlphas();
+  void SaveRawData();
+
+  TH1F* hEunc;
+  TH1F* hEuncSel;
+  TH1F* hEcal;
+  TH2F* hDPvsEcal;
+  TH2F* hESHvsEPR;
+
+ private:
+
+  ////  Int_t fRunNumber;
+  string fRunNumber;
+  Double_t fLoThr;     // Low and high thresholds on the normalized uncalibrated
+  Double_t fHiThr;     // energy deposition.
+  UInt_t fNev;         // Number of processed events.
+
+  static const UInt_t fMinHitCount = 5;     // Minimum number of hits for a PMT
+                                            // to be calibrated.
+
+  TTree* fTree;
+  UInt_t fNentries;
+
+  // Quantities for calculations of the calibration constants.
+
+  Double_t fe0;
+  Double_t fqe[THcPShTrack::fNpmts];
+  Double_t fq0[THcPShTrack::fNpmts];
+  Double_t fQ[THcPShTrack::fNpmts][THcPShTrack::fNpmts];
+  Double_t falphaU[THcPShTrack::fNpmts];   // 'unconstrained' calib. constants
+  Double_t falphaC[THcPShTrack::fNpmts];   // the sought calibration constants
+  Double_t falpha0[THcPShTrack::fNpmts];   // initial gains
+  Double_t falpha1[THcPShTrack::fNpmts];   // unit gains
+
+  UInt_t fHitCount[THcPShTrack::fNpmts];
+
+};
+
+//------------------------------------------------------------------------------
+
+THcPShowerCalib::THcPShowerCalib() {};
+
+//------------------------------------------------------------------------------
+
+////THcPShowerCalib::THcPShowerCalib(Int_t RunNumber) {
+THcPShowerCalib::THcPShowerCalib(string RunNumber) {
+  fRunNumber = RunNumber;
+};
+
+//------------------------------------------------------------------------------
+
+THcPShowerCalib::~THcPShowerCalib() {
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::SaveRawData() {
+
+  // Output raw data into file for debug purposes. To be called after
+  // calibration constants are determined.
+
+  cout << "SaveRawData: Output raw data into Pcal_calib.raw_data." << endl;
+
+  ofstream fout;
+  fout.open("Pcal_calib.raw_data",ios::out);
+
+  THcPShTrack trk;
+
+  for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+    if (ReadShRawTrack(trk, ientry)) {
+      trk.SetEs(falphaC);
+      trk.Print(fout);
+    }
+
+  }
+
+  fout.close();
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::Init() {
+
+  //Reset ROOT and connect tree file.
+
+  gROOT->Reset();
+
+  //  char* fname = Form("Root_files/Pcal_calib_%d.root",fRunNumber);
+  char* fname = Form("ROOTfiles/shms_replay_%s.root",fRunNumber.c_str());
+  cout << "THcPShowerCalib::Init: Root file name = " << fname << endl;
+
+  TFile *f = new TFile(fname);
+  f->GetObject("T",fTree);
+
+  fNentries = fTree->GetEntries();
+  cout << "THcPShowerCalib::Init: fNentries= " << fNentries << endl;
+
+  // Histogram declarations.
+
+  //  hEunc = new TH1F("hEunc", "Edep/P uncalibrated", 500, 0., 5.);
+  //  hEcal = new TH1F("hEcal", "Edep/P calibrated", 150, 0., 1.5);
+  //  hDPvsEcal = new TH2F("hDPvsEcal", "#DeltaP versus Edep/P ",
+  //		       350,0.,1.5, 250,-12.5,22.5);
+  hEunc = new TH1F("hEunc", "Edep/P uncalibrated", 500, 0., 10.);
+  hEcal = new TH1F("hEcal", "Edep/P calibrated", 200, 0., 2.);
+  hDPvsEcal = new TH2F("hDPvsEcal", "#DeltaP versus Edep/P ",
+		       400,0.,2., 250,-12.5,22.5);
+  hESHvsEPR = new TH2F("hESHvsEPR", "E_{SH} versus E_{PR}",
+		       300,0.,1.5, 300,0.,1.5);
+
+
+  // Initialize cumulative quantities.
+  
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) fHitCount[i] = 0;
+
+  fe0 = 0.;
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+    fqe[i] = 0.;
+    fq0[i] = 0.;
+    falphaU[i] = 0.;
+    falphaC[i] = 0.;
+    for (UInt_t j=0; j<THcPShTrack::fNpmts; j++) {
+      fQ[i][j] = 0.;
+    }
+  }
+
+  // Initial gains, 1 for all.
+
+  for (UInt_t ipmt=0; ipmt<THcPShTrack::fNpmts; ipmt++) {
+    falpha0[ipmt] = 1.;
+  };
+
+  // Unit gains.
+
+  for (UInt_t ipmt=0; ipmt<THcPShTrack::fNpmts; ipmt++) {
+    falpha1[ipmt] = 1.;
+  }
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::CalcThresholds() {
+
+  // Calculate +/-3 RMS thresholds on the uncalibrated total energy
+  // depositions. These thresholds are used mainly to exclude potential
+  // hadronic events due to the gas Cherenkov inefficiency.
+
+  // Histogram uncalibrated energy depositions, get mean and RMS from the
+  // histogram, establish +/-3 * RMS thresholds.
+
+  cout << "THcPShowerCalib::CalcThresholds: FNentries = " << fNentries << endl;
+
+  Int_t nev = 0;
+  THcPShTrack trk;
+
+  for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+    if ( ReadShRawTrack(trk, ientry)) {
+
+	//    trk.Print(cout);
+	//    getchar();
+
+	trk.SetEs(falpha0);             //Use initial gain constants here.
+	Double_t Enorm = trk.Enorm();
+
+	hEunc->Fill(Enorm);
+
+	nev++;
+    //    if (nev%100000 == 0)
+    //      cout << "CalcThreshods: nev=" << nev << "  Enorm=" << Enorm << endl;
+    }
+
+  };
+
+  Double_t mean = hEunc->GetMean();
+  Double_t rms = hEunc->GetRMS();
+  cout << "CalcThreshods: mean=" << mean << "  rms=" << rms << endl;
+
+  fLoThr = mean - 3.*rms;
+  fHiThr = mean + 3.*rms;
+  //  fLoThr = 0.;              // Wide open thrsholds, for
+  //  fHiThr = 1.e+8;           // comparison with the old code.
+
+  cout << "CalcThreshods: fLoThr=" << fLoThr << "  fHiThr=" << fHiThr 
+       << "  nev=" << nev << endl;
+
+  Int_t nbins = hEunc->GetNbinsX();
+  Int_t nlo = hEunc->FindBin(fLoThr);
+  Int_t nhi = hEunc->FindBin(fHiThr);
+
+  cout << "CalcThresholds: nlo=" << nlo << "  nhi=" << nhi 
+       << "  nbins=" << nbins << endl;
+
+  // Histogram of selected within the thresholds events.
+  
+  hEuncSel = (TH1F*)hEunc->Clone("hEuncSel");
+  
+  for (Int_t i=0; i<nlo; i++) hEuncSel->SetBinContent(i, 0.);
+  for (Int_t i=nhi; i<nbins+1; i++) hEuncSel->SetBinContent(i, 0.);
+
+};
+
+//------------------------------------------------------------------------------
+
+bool THcPShowerCalib::ReadShRawTrack(THcPShTrack &trk, UInt_t ientry) {
+
+  //
+  // Set a Shower track event from ntuple ientry.
+  //
+
+  // Declaration of leaves types
+
+  // Preshower and Shower ADC signals.
+
+  ////Double_t        P_pr_a_p[THcPShTrack::fNrows_pr][THcPShTrack::fNcols_pr];
+  ////Double_t        P_sh_a_p[THcPShTrack::fNrows_sh][THcPShTrack::fNcols_sh];
+  Double_t        P_pr_apos_p[THcPShTrack::fNrows_pr];
+  Double_t        P_pr_aneg_p[THcPShTrack::fNrows_pr];
+  Double_t        P_sh_a_p[THcPShTrack::fNcols_sh][THcPShTrack::fNrows_sh];
+
+  // Track parameters.
+
+  double          P_tr_n;
+  Double_t        P_tr_p;
+  Double_t        P_tr_x;   //X FP
+  Double_t        P_tr_xp;
+  Double_t        P_tr_y;   //Y FP
+  Double_t        P_tr_yp;
+
+  Double_t        P_tr_tg_dp;
+
+  Double_t        P_hgcer_npe[4];
+ 
+  ////  const Double_t adc_thr = 15.;   //Low threshold on the ADC signals.
+
+  // Set branch addresses.
+
+  fTree->SetBranchAddress("P.cal.pr.apos_p",P_pr_apos_p);
+  fTree->SetBranchAddress("P.cal.pr.aneg_p",P_pr_aneg_p);
+  fTree->SetBranchAddress("P.cal.fly.a_p",P_sh_a_p);
+
+  fTree->SetBranchAddress("P.tr.n", &P_tr_n);
+  fTree->SetBranchAddress("P.tr.x", &P_tr_x);
+  fTree->SetBranchAddress("P.tr.y", &P_tr_y);
+  fTree->SetBranchAddress("P.tr.th",&P_tr_xp);
+  fTree->SetBranchAddress("P.tr.ph",&P_tr_yp);
+  fTree->SetBranchAddress("P.tr.p", &P_tr_p);
+
+  fTree->SetBranchAddress("P.tr.tg_dp", &P_tr_tg_dp);
+ 
+  fTree->SetBranchAddress("P.hgcer.npe", P_hgcer_npe);
+
+  fTree->GetEntry(ientry);
+
+  if (ientry%100000 == 0) cout << "   ReadShRawTrack: " << ientry << endl;
+
+  ////
+  //  cout << "   ReadShRawTrack: tracks " << P_tr_n << endl;
+  //  cout << "      delta = " << P_tr_tg_dp << endl;
+  //  cout << "      x, y = " << P_tr_x + P_tr_xp*D_CALO_FP << " " 
+  //       << P_tr_y + P_tr_yp*D_CALO_FP << endl;
+  //  cout << "      Cer: " << P_hgcer_npe[0] << " " <<  P_hgcer_npe[1] << " "
+  //       << P_hgcer_npe[2] << " " << P_hgcer_npe[3] << endl;
+
+  if (P_tr_n != 1) return 0;
+
+  bool good_trk =   P_tr_tg_dp > DELTA_MIN &&
+		    P_tr_tg_dp < DELTA_MAX &&
+		    P_tr_x + P_tr_xp*D_CALO_FP > XMIN &&
+		    P_tr_x + P_tr_xp*D_CALO_FP < XMAX &&
+                    P_tr_y + P_tr_yp*D_CALO_FP > YMIN &&
+		    P_tr_y + P_tr_yp*D_CALO_FP < YMAX ;
+  if (!good_trk) return 0;
+
+  bool good_hgcer = P_hgcer_npe[0] > 2 ||
+		    P_hgcer_npe[1] > 2 ||
+		    P_hgcer_npe[2] > 2 ||
+		    P_hgcer_npe[3] > 2  ;
+  if(!good_hgcer) return 0;
+
+  //  cout << "    Track is good." << endl << endl;
+  //  getchar();
+
+  // Set track coordinates and slopes at the face of Preshower.
+
+  trk.Reset(P_tr_p, P_tr_tg_dp, P_tr_x+D_CALO_FP*P_tr_xp, P_tr_xp,
+  	    P_tr_y+D_CALO_FP*P_tr_yp, P_tr_yp);
+
+  // Set Preshower hits.
+
+  UInt_t nb = 0;
+
+  for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+    for (UInt_t j=0; j<THcPShTrack::fNrows_pr; j++) {
+      nb++;
+      Double_t adc;
+      switch (k) {
+      case 0: adc = P_pr_aneg_p[j]; break;
+      case 1: adc = P_pr_apos_p[j]; break;
+      default : cout << "*** Wrong PreShower column! ***" << endl;
+      }
+      if (adc > PR_ADC_THR) trk.AddHit(adc, 0., nb);
+    }
+  }
+
+  ////  for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+  ////    for (UInt_t j=0; j<THcPShTrack::fNrows_pr; j++) {
+  ////
+  ////      Double_t adc = P_pr_a_p[j][k];
+  ////
+  ////      if (adc > adc_thr) {
+  ////    if (adc > PR_ADC_THR) {
+  ////	UInt_t nb = j+1 + k*THcPShTrack::fNrows_pr;
+  ////	trk.AddHit(adc, 0., nb);
+  ////      }
+  ////
+  ////    }
+  ////  }
+
+  // Set Shower hits.
+
+  for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+    for (UInt_t j=0; j<THcPShTrack::fNrows_sh; j++) {
+      nb++;
+      Double_t adc = P_sh_a_p[k][j];
+      if (adc > SH_ADC_THR) {
+	trk.AddHit(adc, 0., nb);
+      }
+    }
+  }
+
+  ////  for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+  ////    for (UInt_t j=0; j<THcPShTrack::fNrows_sh; j++) {
+  ////
+  ////      Double_t adc = P_sh_a_p[j][k];
+  ////
+  ////      if (adc > SH_ADC_THR) {
+  ////	UInt_t nb = THcPShTrack::fNpmts_pr + j+1 + k*THcPShTrack::fNrows_sh;
+  ////	trk.AddHit(adc, 0., nb);
+  ////      }
+  ////
+  ////    }
+  ////  }
+
+  return 1;
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::ComposeVMs() {
+
+  //
+  // Fill in vectors and matrixes for the gain constant calculations.
+  //
+
+  fNev = 0;
+  THcPShTrack trk;
+
+  // Loop over the shower track events in the ntuples.
+
+  for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+    if (ReadShRawTrack(trk, ientry)) {
+
+      // Set energy depositions with default gains.
+      // Calculate normalized to the track momentum total energy deposition,
+      // check it against the thresholds.
+
+      trk.SetEs(falpha0);
+      Double_t Enorm = trk.Enorm();
+      if (Enorm>fLoThr && Enorm<fHiThr) {
+
+	trk.SetEs(falpha1);   // Set energies with unit gains for now.
+	// trk.Print(cout);
+
+	fe0 += trk.GetP();    // Accumulate track momenta.
+
+	vector<pmt_hit> pmt_hit_list;     // Container to save PMT hits
+
+	// Loop over hits.
+
+	for (UInt_t i=0; i<trk.GetNhits(); i++) {
+
+	  THcPShHit* hit = trk.GetHit(i);
+	  // hit->Print(cout);
+
+	  UInt_t nb = hit->GetBlkNumber();
+
+	  // Fill the qe and q0 vectors.
+
+	  fqe[nb-1] += hit->GetEdep() * trk.GetP();
+	  fq0[nb-1] += hit->GetEdep();
+
+	  // Save the PMT hit.
+
+	  pmt_hit_list.push_back( pmt_hit{hit->GetEdep(), nb} );
+
+	  fHitCount[nb-1]++;   //Accrue the hit counter.
+
+	}      //over hits
+
+	// Fill in the correlation matrix Q by retrieving the PMT hits.
+
+	for (vector<pmt_hit>::iterator i=pmt_hit_list.begin();
+	     i < pmt_hit_list.end(); i++) {
+
+	  UInt_t ic = (*i).channel;
+	  Double_t is = (*i).signal;
+
+	  for (vector<pmt_hit>::iterator j=i;
+	       j < pmt_hit_list.end(); j++) {
+
+	    UInt_t jc = (*j).channel;
+	    Double_t js = (*j).signal;
+
+	    fQ[ic-1][jc-1] += is*js;
+	    if (jc != ic) fQ[jc-1][ic-1] += is*js;
+	  }
+	}
+
+	fNev++;
+
+      };   // if within enorm thresholds
+
+    };   // success in reading
+
+  };     // over entries
+
+  // Take averages.
+
+  fe0 /= fNev;
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+    fqe[i] /= fNev;
+    fq0[i] /= fNev;
+  }
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+    for (UInt_t j=0; j<THcPShTrack::fNpmts; j++)
+      fQ[i][j] /= fNev;
+
+  // Output vectors and matrixes, for debug purposes.
+
+  ofstream q0out;
+  q0out.open("q0.deb",ios::out);
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+    q0out << setprecision(20) << fq0[i] << " " << i << endl;
+  q0out.close();
+
+  ofstream qeout;
+  qeout.open("qe.deb",ios::out);
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+    qeout << setprecision(20) << fqe[i] << " " << i << endl;
+  qeout.close();
+
+  ofstream Qout;
+  Qout.open("Q.deb",ios::out);
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+    for (UInt_t j=0; j<THcPShTrack::fNpmts; j++)
+      Qout << setprecision(20) << fQ[i][j] << " " << i << " " << j << endl;
+  Qout.close();
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::SolveAlphas() {
+
+  //
+  // Solve for the sought calibration constants, by use of the Root
+  // matrix algebra package.
+  //
+
+  TMatrixD Q(THcPShTrack::fNpmts,THcPShTrack::fNpmts);
+  TVectorD q0(THcPShTrack::fNpmts);
+  TVectorD qe(THcPShTrack::fNpmts);
+  TVectorD au(THcPShTrack::fNpmts);
+  TVectorD ac(THcPShTrack::fNpmts);
+  Bool_t ok;
+
+  cout << "Solving Alphas..." << endl;
+  cout << endl;
+
+  // Print out hit numbers.
+
+  cout << "Hit counts:" << endl;
+  UInt_t j = 0;
+  
+  for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+    k==0 ? cout << "Preshower:" : cout << "        :";
+    for (UInt_t i=0; i<THcPShTrack::fNrows_pr; i++)
+      cout << setw(6) << fHitCount[j++] << ",";
+    cout << endl;
+  }
+
+  for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+    k==0 ? cout << "Shower   :" : cout << "        :";
+    for (UInt_t i=0; i<THcPShTrack::fNrows_sh; i++)
+      cout << setw(6) << fHitCount[j++] << ",";
+    cout << endl;
+  }
+
+  // Initialize the vectors and the matrix of the Root algebra package.
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+    q0[i] = fq0[i];
+    qe[i] = fqe[i];
+    for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+      Q[i][k] = fQ[i][k];
+    }
+  }
+
+  // Sanity check.
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+
+    // Check zero hit channels: the vector and matrix elements should be 0.
+
+    if (fHitCount[i] == 0) {
+
+      if (q0[i] != 0. || qe[i] != 0.) {
+
+	cout << "*** Inconsistency in chanel " << i << ": # of hits  "
+	     << fHitCount[i] << ", q0=" << q0[i] << ", qe=" << qe[i];
+
+	for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+	  if (Q[i][k] !=0. || Q[k][i] !=0.)
+	    cout << ", Q[" << i << "," << k << "]=" << Q[i][k]
+		 << ", Q[" << k << "," << i << "]=" << Q[k][i];
+	}
+
+	cout << " ***" << endl;
+      }
+    }
+
+    // The hit channels: the vector elements should be non zero.
+
+    if ( (fHitCount[i] != 0) && (q0[i] == 0. || qe[i] == 0.) ) {
+      cout << "*** Inconsistency in chanel " << i << ": # of hits  "
+	   << fHitCount[i] << ", q0=" << q0[i] << ", qe=" << qe[i]
+	   << " ***" << endl;
+    }
+
+  } //sanity check
+
+  // Low hit number channels: exclude from calculation. Assign all the
+  // correspondent elements 0, except self-correlation Q(i,i)=1.
+
+  cout << endl;
+  cout << "Channels with hit number less than " << fMinHitCount 
+       << " will not be calibrated." << endl;
+  cout << endl;
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+
+    if (fHitCount[i] < fMinHitCount) {
+      cout << "Channel " << i << ", " << fHitCount[i]
+	   << " hits, will not be calibrated." << endl;
+      q0[i] = 0.;
+      qe[i] = 0.;
+      for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+	Q[i][k] = 0.;
+	Q[k][i] = 0.;
+      }
+      Q[i][i] = 1.;
+    }
+
+  }
+
+  // Declare LU decomposition method for the correlation matrix Q.
+
+  TDecompLU lu(Q);
+  Double_t d1,d2;
+  lu.Det(d1,d2);
+  cout << "cond:" << lu.Condition() << endl;
+  cout << "det :" << d1*TMath::Power(2.,d2) << endl;
+  cout << "tol :" << lu.GetTol() << endl;
+
+  // Solve equation Q x au = qe for the 'unconstrained' calibration (gain)
+  // constants au.
+
+  au = lu.Solve(qe,ok);
+  cout << "au: ok=" << ok << endl;
+  //  au.Print();
+
+  // Find the sought 'constrained' calibration constants next.
+
+  Double_t t1 = fe0 - au * q0;         // temporary variable.
+  //  cout << "t1 =" << t1 << endl;
+
+  TVectorD Qiq0(THcPShTrack::fNpmts);   // an intermittent result
+  Qiq0 = lu.Solve(q0,ok);
+  cout << "Qiq0: ok=" << ok << endl;
+  //  Qiq0.Print();
+
+  Double_t t2 = q0 * Qiq0;             // another temporary variable
+  //  cout << "t2 =" << t2 << endl;
+
+  ac = (t1/t2) *Qiq0 + au;             // the sought gain constants
+  // cout << "ac:" << endl;
+  //  ac.Print();
+
+  // Assign the gain arrays.
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+    falphaU[i] = au[i];
+    falphaC[i] = ac[i];
+  }
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::FillHEcal() {
+
+  //
+  // Fill histogram of the normalized energy deposition, and 2-d histogram
+  // of momentum deviation versus normalized energy deposition.
+  // Output event by event energy depositions and momenta for debug purposes.
+  //
+
+  ofstream output;
+  output.open("calibrated.deb",ios::out);
+
+  Int_t nev = 0;
+
+  THcPShTrack trk;
+
+  //  Double_t delta;
+  //  fTree->SetBranchAddress("P.tr.tg_dp",&delta);
+
+  for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+    if (ReadShRawTrack(trk, ientry)) {
+      //    trk.Print(cout);
+
+      trk.SetEs(falphaC);       // use the 'constrained' calibration constants
+      Double_t P = trk.GetP();
+      Double_t delta = trk.GetDp();
+      Double_t Enorm = trk.Enorm();
+
+      hEcal->Fill(Enorm);
+
+      hDPvsEcal->Fill(Enorm,delta,1.);
+
+      hESHvsEPR->Fill(trk.EPRnorm(), trk.ESHnorm());
+
+      output << Enorm*P/1000. << " " << P/1000. << " " << delta << " "
+	     << trk.GetX() << " " << trk.GetY() << endl;
+
+      nev++;
+    }
+
+  };
+
+  output.close();
+
+  cout << "FillHEcal: " << nev << " events filled" << endl;
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::SaveAlphas() {
+
+  //
+  // Output the gain constants in a format suitable for inclusion in the
+  // pcal.param file to be used in the analysis.
+  //
+
+  ofstream output;
+  ////  char* fname = Form("pcal.param.%d",fRunNumber);
+  char* fname = Form("pcal.param.%s",fRunNumber.c_str());
+  cout << "SaveAlphas: fname=" << fname << endl;
+
+  output.open(fname,ios::out);
+
+  output << "; Calibration constants for run " << fRunNumber 
+	 << ", " << fNev << " events processed" << endl;
+  output << endl;
+
+  UInt_t j = 0;
+
+  for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+    k==0 ? output << "pcal_neg_gain_cor =" : output << "pcal_pos_gain_cor =";
+    for (UInt_t i=0; i<THcPShTrack::fNrows_pr; i++)
+      output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+    output << endl;
+  }
+
+  for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+    k==0 ? output << "pcal_arr_gain_cor =" : output << "                  ";
+    for (UInt_t i=0; i<THcPShTrack::fNrows_sh; i++)
+      output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+    output << endl;
+  }
+
+  output.close();
+}
+
+#endif
diff --git a/CALIBRATION/shms_cal_calib/pcal_calib.cpp b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..75015706e5a4ef230d1450c29e8ed32caadfccb4
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
@@ -0,0 +1,59 @@
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TF1.h>
+#include "THcPShowerCalib.h"
+
+//
+// A steering Root script for the SHMS calorimeter calibration.
+//
+
+////void pcal_calib(Int_t RunNumber) {
+void pcal_calib(string RunNumber) {
+ 
+ cout << "Calibrating run " << RunNumber << endl;
+
+ THcPShowerCalib theShowerCalib(RunNumber);
+
+ theShowerCalib.Init();            // Initialize constants and variables
+ theShowerCalib.CalcThresholds();  // Thresholds on the uncalibrated Edep/P
+ theShowerCalib.ComposeVMs();      // Compute vectors amd matrices for calib.
+ theShowerCalib.SolveAlphas();     // Solve for the calibration constants
+ theShowerCalib.SaveAlphas();      // Save the constants
+ //theShowerCalib.SaveRawData();   // Save raw data into file for debug purposes
+ theShowerCalib.FillHEcal();       // Fill histograms
+
+ // Plot histograms
+
+ TCanvas* Canvas =
+   new TCanvas("Canvas", "PHMS Shower Counter calibration", 1000, 667);
+ Canvas->Divide(2,2);
+
+ Canvas->cd(1);
+
+ // Normalized uncalibrated energy deposition.
+
+ theShowerCalib.hEunc->DrawCopy();
+  
+ theShowerCalib.hEuncSel->SetFillColor(kGreen);
+ theShowerCalib.hEuncSel->DrawCopy("same");
+
+ Canvas->cd(2);
+ theShowerCalib.hESHvsEPR->Draw();
+
+ // Normalized energy deposition after calibration.
+
+ Canvas->cd(3);
+ gStyle->SetOptFit();
+
+ theShowerCalib.hEcal->Fit("gaus","","",0.9,1.1);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineColor(2);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineWidth(1);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineStyle(1);
+
+ // SHMS delta(P) versus the calibrated energy deposition.
+
+ Canvas->cd(4);
+ theShowerCalib.hDPvsEcal->Draw();
+
+}