Merge pull request #106 from hallc-online/vardan-kpp-work

Vardan kpp work

CALIBRATION/hms_cal_calib/THcShHit.h

+#include <iostream>
+
+// HMS calorimeter hit class for calibration.
+
+class THcShHit {
+
+  Double_t ADCpos, ADCneg;   // pedestal subtracted ADC signals.
+  Double_t Epos, Eneg;       // Energy depositions seen from pos. & neg. sides
+  UInt_t BlkNumber;
+
+ public:
+
+  THcShHit();
+  THcShHit(Double_t adc_pos, Double_t adc_neg,
+	   UInt_t blk_number);
+  ~THcShHit();
+
+  void SetADCpos(Double_t sig) {ADCpos = sig;}
+
+  void SetADCneg(Double_t sig) {ADCneg = sig;}
+
+  void SetEpos(Double_t e) {Epos = e;}
+
+  void SetEneg(Double_t e) {Eneg = e;}
+
+  void SetBlkNumber(UInt_t n) {BlkNumber = n;}
+
+  Double_t GetADCpos() {return ADCpos;}
+
+  Double_t GetADCneg() {return ADCneg;}
+
+  Double_t GetEpos() {return Epos;}
+
+  Double_t GetEneg() {return Eneg;}
+  
+  UInt_t GetBlkNumber() {return BlkNumber;}
+
+  void Print(ostream & ostrm);
+};
+
+//------------------------------------------------------------------------------
+
+THcShHit::THcShHit() {
+  ADCpos = -99999.;
+  ADCneg = -99999.;
+  Epos = -99999.;
+  Eneg = -99999.;
+  BlkNumber = 99999;
+};
+
+THcShHit::THcShHit(Double_t adc_pos, Double_t adc_neg,
+		   UInt_t blk_number) {
+  ADCpos = adc_pos;
+  ADCneg = adc_neg;
+  Epos = 0.;
+  Eneg = 0.;
+  BlkNumber = blk_number;
+};
+
+THcShHit::~THcShHit() { };
+
+//------------------------------------------------------------------------------
+
+void THcShHit::Print(ostream & ostrm) {
+
+  // Output hit data through the stream ostrm.
+
+  ostrm << ADCpos << " " << ADCneg << " " << Epos << " " << Eneg << " "
+	<< BlkNumber << endl;
+};
+
+struct pmt_hit {Double_t signal; UInt_t channel;};

CALIBRATION/hms_cal_calib/THcShTrack.h

+#include "THcShHit.h"
+#include "TMath.h"
+
+#include <vector>
+#include <iterator>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+// Track class for the HMS calorimeter calibration.
+// Comprises the spectrometer track parameters and calorimeter hits.
+//
+
+// Container (collection) of hits and its iterator.
+//
+typedef vector<THcShHit*> THcShHitList;
+typedef THcShHitList::iterator THcShHitIt;
+
+class THcShTrack {
+
+  Double_t P;   // track momentum
+  Double_t Dp;  // track momentum deviation, %
+  Double_t X;   // at the calorimater face
+  Double_t Xp;  // slope
+  Double_t Y;   // at the calorimater face
+  Double_t Yp;  // slope
+
+  THcShHitList Hits;
+
+ public:
+  THcShTrack();
+  THcShTrack(Double_t p, Double_t dp,
+	     Double_t x, Double_t xp, Double_t y, Double_t yp);
+  ~THcShTrack();
+
+  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_pos, Double_t adc_neg,
+	      Double_t e_pos, Double_t e_neg,
+	      UInt_t blk_number);
+
+  THcShHit* 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 ETAnorm();
+
+  Double_t GetP() {return P*1000.;}      //MeV
+
+  Double_t GetDp() {return Dp;}
+
+  Double_t GetX() {return X;}
+  Double_t GetY() {return Y;}
+
+  Float_t Ycor(Double_t);         // coord. corection for single PMT module
+  Float_t Ycor(Double_t, Int_t);  // coord. correction for double PMT module
+
+  // Coordinate correction constants from hcana.param.
+  //
+  static constexpr Double_t fAcor = 200.;
+  static constexpr Double_t fBcor = 8000.;
+  static constexpr Double_t fCcor = 64.36;
+  static constexpr Double_t fDcor = 1.66;
+
+  // Calorimeter geometry constants.
+  //
+  static constexpr Double_t fZbl = 10;   //cm, block transverse size
+  static constexpr UInt_t fNrows = 13;
+  static constexpr UInt_t fNcols =  4;
+  static constexpr UInt_t fNnegs = 26;  // number of blocks with neg. side PMTs.
+  static constexpr UInt_t fNpmts = 78;  // total number of PMTs.
+  static constexpr UInt_t fNblks = fNrows*fNcols;
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShTrack::THcShTrack() { };
+
+THcShTrack::THcShTrack(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 THcShTrack::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 THcShTrack::AddHit(Double_t adc_pos, Double_t adc_neg,
+			Double_t e_pos, Double_t e_neg,
+			UInt_t blk_number) {
+
+  // Add a hit to the hit list.
+
+  THcShHit* hit = new THcShHit(adc_pos, adc_neg, blk_number);
+  hit->SetEpos(e_pos);
+  hit->SetEneg(e_neg);
+  Hits.push_back(hit);
+};
+
+//------------------------------------------------------------------------------
+
+THcShHit* THcShTrack::GetHit(UInt_t k) {
+  THcShHitIt it = Hits.begin();
+  for (UInt_t i=0; i<k; i++) it++;
+  return *it;
+}
+
+void THcShTrack::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 (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    (*iter)->Print(ostrm);
+  };
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShTrack::~THcShTrack() {
+  for (THcShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
+    delete *i;
+    *i = 0;
+  }
+};
+
+//------------------------------------------------------------------------------
+
+void THcShTrack::SetEs(Double_t* alpha) {
+
+  // Set hit energy depositions seen from postive and negative sides,
+  // by use of calibration (gain) constants alpha.
+  
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+  
+    Double_t adc_pos = (*iter)->GetADCpos();
+    Double_t adc_neg = (*iter)->GetADCneg();
+    UInt_t nblk = (*iter)->GetBlkNumber();
+
+    Int_t ncol=(nblk-1)/fNrows+1;
+    Double_t xh=X+Xp*(ncol-0.5)*fZbl;
+    Double_t yh=Y+Yp*(ncol-0.5)*fZbl;
+    if (nblk <= fNnegs) {
+      (*iter)->SetEpos(adc_pos*Ycor(yh,0)*alpha[nblk-1]);
+      (*iter)->SetEneg(adc_neg*Ycor(yh,1)*alpha[fNblks+nblk-1]);
+    }
+    else {
+      (*iter)->SetEpos(adc_pos*Ycor(yh)*alpha[nblk-1]);
+      (*iter)->SetEneg(0.);
+    };
+
+  };
+
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcShTrack::Enorm() {
+
+  // Normalized to the track momentum energy depostion in the calorimeter.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    sum += (*iter)->GetEpos();
+    sum += (*iter)->GetEneg();
+  };
+
+  return sum/P/1000.;
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcShTrack::EPRnorm() {
+
+  // Normalized to the track momentum energy depostion in Preshower.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    UInt_t nblk = (*iter)->GetBlkNumber();
+    Int_t ncol=(nblk-1)/fNrows+1;
+    if (ncol==1) {
+      sum += (*iter)->GetEpos();
+      sum += (*iter)->GetEneg();
+    }
+  }
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+//------------------------------------------------------------------------------
+
+Double_t THcShTrack::ETAnorm() {
+
+  // Normalized to the track momentum energy depostion in Preshower.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    UInt_t nblk = (*iter)->GetBlkNumber();
+    Int_t ncol=(nblk-1)/fNrows+1;
+    if (ncol!=1) {
+      sum += (*iter)->GetEpos();
+      sum += (*iter)->GetEneg();
+    }
+  }
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+
+//------------------------------------------------------------------------------
+
+//Coordinate correction for single PMT modules.
+//PMT attached at right (positive) side.
+
+Float_t THcShTrack::Ycor(Double_t y) {
+  return TMath::Exp(y/fAcor)/(1. + y*y/fBcor);
+}
+
+//Coordinate correction for double PMT modules.
+//
+
+Float_t THcShTrack::Ycor(Double_t y, Int_t side) {
+  if (side!=0&&side!=1) {
+    cout << "THcShower::Ycor : wrong side " << side << endl;
+    return 0.;
+  }
+  Int_t sign = 1 - 2*side;
+  return (fCcor + sign*y)/(fCcor + sign*y/fDcor);
+}

CALIBRATION/hms_cal_calib/hcal_calib.cpp

+#include "TCanvas.h"
+#include <TStyle.h>
+#include <TH1.h>
+#include <TF1.h>
+#include "THcShowerCalib.h"
+
+//
+// A steering Root script for the HMS calorimeter calibration.
+//
+
+void hcal_calib(string RunNumber) {
+
+  // Initialize the analysis clock
+  clock_t t = clock();
+  
+ cout << "Calibrating run " << RunNumber << endl;
+
+ THcShowerCalib 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", "HMS 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");
+
+ // E(Tail) vs E(Preshower) plot.
+
+ Canvas->cd(2);
+ theShowerCalib.hETAvsEPR->Draw();
+
+ // Normalized energy deposition after calibration.
+
+ Canvas->cd(3);
+ gStyle->SetOptFit();
+
+ theShowerCalib.hEcal->Fit("gaus","","",0.95,1.05);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineColor(2);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineWidth(1);
+ theShowerCalib.hEcal->GetFunction("gaus")->SetLineStyle(1);
+
+ // HMS delta(P) versus the calibrated energy deposition.
+
+ Canvas->cd(4);
+ theShowerCalib.hDPvsEcal->Draw();
+
+ // Calculate the analysis rate
+ t = clock() - t;
+ printf ("The analysis took %.1f seconds \n", ((float) t) / CLOCKS_PER_SEC);
+}

CALIBRATION/hms_cal_calib/howto.txt

+How to calibrate HMS calorimeter with real electrons.
+
+The calibration scripts reside in
+hallc_replay/CALIBRATION/hms_cal_calib directory. They consist of 3
+header files called THcShHit.h,THcShTrack.h, THcShowerCalib.h, and a
+steering script hcal_calib.cpp.
+
+The scripts work on root files from hcana analysis and make use of
+quantities pertained to tracking, gas Cherenkov, and TOF from
+hodoscopes. Hence it is convenient to calibrate on root files from
+full HMS analysis. The root files are assumed in a linked ROOTfiles
+directory, under names hms_replay_<run_number>.root.
+
+Once your hcana, hallc_replay and Root are set up, you can compile and
+run hcal_calib under hcana, by issuing command
+
+.x hcal_calib.cpp+("run number").
+
+For instance, for calibration on hms_replay_303.root file in
+ROOTfiles, the correct commad would be
+
+.x pcal_calib.cpp+("303") .
+
+Upon calibration, a canvas with representative plots will pop up. The
+calibration constants will be written in output file
+hcal.param.<run_number>, in a format suitable for plugging them into
+your hcal.param file for subsequent use.
+
+If you want to modify selection cuts used in calibration (cuts on
+delta, beta, gas Cherenkov signals), you can find them at the
+beginning of THcShowerCalib.h file, in the #define directives.

CALIBRATION/shms_cal_calib/ROOTfiles

+/net/cdaq/cdaql3data/cdaq/hallc-online/ROOTfiles
\ No newline at end of file

CALIBRATION/shms_cal_calib/THcPShHit.h

+#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;};

CALIBRATION/shms_cal_calib/THcPShTrack.h

+#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;}
+
+  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 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 == 2) || (yhit > 0. && ncol == 1))
+    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;
+}

CALIBRATION/shms_cal_calib/THcPShowerCalib.h

+#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"
+
+#include <time.h>
+
+#define D_CALO_FP 275.    //distance from FP to the Preshower
+
+//Whole calorimeter fid. limits
+#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 0
+#define SH_ADC_THR 0
+
+#define BETA_MIN 0.5
+#define BETA_MAX 1.5
+
+#define HGCER_NPE_MIN 2
+
+#define MIN_HIT_COUNT 5   //Minimum number of hits for a PMT to be calibrated.
+
+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:
+
+  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.
+
+  TTree* fTree;
+  UInt_t fNentries;
+
+  // Declaration of leaves types
+
+  // Preshower and Shower ADC signals.
+
+  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];
+  Double_t        P_tr_beta;
+
+  Double_t        P_cal_nclust;
+
+  TBranch* b_P_tr_p;
+  TBranch* b_P_pr_apos_p;
+  TBranch* b_P_pr_aneg_p;
+  TBranch* b_P_sh_a_p;
+  TBranch* b_P_tr_n;
+  TBranch* b_P_tr_x;
+  TBranch* b_P_tr_y;
+  TBranch* b_P_tr_xp;
+  TBranch* b_P_tr_yp;
+  TBranch* b_P_tr_tg_dp;
+  TBranch* b_P_hgcer_npe;
+  TBranch* b_P_tr_beta;
+
+  TBranch* b_P_cal_nclust;
+
+  // 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(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("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;
+
+  fTree->SetBranchAddress("P.cal.pr.apos_p",P_pr_apos_p,&b_P_pr_apos_p);
+  fTree->SetBranchAddress("P.cal.pr.aneg_p",P_pr_aneg_p,&b_P_pr_aneg_p);
+  fTree->SetBranchAddress("P.cal.fly.a_p",  P_sh_a_p,   &b_P_sh_a_p);
+
+  fTree->SetBranchAddress("P.tr.n", &P_tr_n,&b_P_tr_n);
+  fTree->SetBranchAddress("P.tr.x", &P_tr_x,&b_P_tr_x);
+  fTree->SetBranchAddress("P.tr.y", &P_tr_y,&b_P_tr_y);
+  fTree->SetBranchAddress("P.tr.th",&P_tr_xp,&b_P_tr_xp);
+  fTree->SetBranchAddress("P.tr.ph",&P_tr_yp,&b_P_tr_yp);
+  fTree->SetBranchAddress("P.tr.p", &P_tr_p,&b_P_tr_p);
+
+  fTree->SetBranchAddress("P.tr.tg_dp", &P_tr_tg_dp,&b_P_tr_tg_dp);
+ 
+  fTree->SetBranchAddress("P.hgcer.npe", P_hgcer_npe,&b_P_hgcer_npe);
+
+  fTree->SetBranchAddress("P.tr.beta", &P_tr_beta,&b_P_tr_beta);
+
+  fTree->SetBranchAddress("P.cal.nclust", &P_cal_nclust,&b_P_cal_nclust);
+
+  // Histogram declarations.
+
+  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;
+
+  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.
+  //
+
+  fTree->GetEntry(ientry);
+
+  if (ientry%100000 == 0) cout << "   ReadShRawTrack: " << ientry << endl;
+
+  // Request single electron track in calorimeter's fid. volume.
+  //
+
+  if (P_tr_n != 1) return 0;
+
+  if (P_cal_nclust != 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] > HGCER_NPE_MIN ||
+		    P_hgcer_npe[1] > HGCER_NPE_MIN ||
+		    P_hgcer_npe[2] > HGCER_NPE_MIN ||
+		    P_hgcer_npe[3] > HGCER_NPE_MIN  ;
+  if(!good_hgcer) return 0;
+
+  bool good_beta = P_tr_beta > BETA_MIN &&
+                   P_tr_beta < BETA_MAX ;
+  if(!good_beta) 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);
+    }
+  }
+
+  // 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);
+      }
+    }
+  }
+
+  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 " << MIN_HIT_COUNT 
+       << " will not be calibrated." << endl;
+  cout << endl;
+
+  for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+
+    if (fHitCount[i] < MIN_HIT_COUNT) {
+      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;
+
+  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.%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

CALIBRATION/shms_cal_calib/howto.txt

+How to calibrate SHMS calorimeter with real electrons.
+
+The calibration scripts reside in
+hallc_replay/CALIBRATION/shms_cal_calib directory. They consist of 3
+header files called THcPShHit.h,THcPShTrack.h, THcPShowerCalib.h, and
+a steering script pcal_calib.cpp.
+
+The scripts work on root files from hcana analysis and make use of
+quantities pertained to tracking, heavy gas Cherenkov, and TOF from
+hodoscopes. Hence it is convenient to calibrate on root files from
+full SHMS analysis. The root files are assumed in a linked ROOTfiles
+directory, under names shms_replay_<run_number>.root.
+
+Once your hcana, hallc_replay and Root are set up, you can compile and
+run pcal_calib under hcana, by issuing command
+
+.x pcal_calib.cpp+("run number").
+
+For instance, for calibration on shms_replay_464_50000.root file in
+ROOTfiles, the correct commad would be
+
+.x pcal_calib.cpp+("464_50000") .
+
+Upon calibration, a canvas with representative plots will pop up. The
+calibration constants will be written in output file
+pcal.param.<run_number>, in a format suitable for plugging them into
+your pcal.param file for subsequent use.
+
+If you want to modify selection cuts used in calibration (cuts on
+delta, beta, gas Cherenkov signals), you can find them at the
+beginning of THcPShowerCalib.h file, in the #define directives.

CALIBRATION/shms_cal_calib/pcal_calib.cpp

+#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(string RunNumber) {
+
+  // Initialize the analysis clock
+  clock_t t = clock();
+ 
+  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 debuging
+  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();
+
+  // Calculate the analysis rate
+  t = clock() - t;
+  printf ("The analysis took %.1f seconds \n", ((float) t) / CLOCKS_PER_SEC);
+}

PARAM/HMS/CAL/hcal.param

@@ -53,7 +53,6 @@ hcal_neg_ped_limit =1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,
                     1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000
 
 
-
 hcal_pos_gain_cor= 0.217, 0.133, 0.135, 0.238, 0.168, 0.147, 0.146, 0.281, 0.154, 0.269, 0.224, 0.186, 0.186,
 0.216, 0.151, 0.127, 0.231, 0.187, 0.145, 0.174, 0.202, 0.175, 0.199, 0.290, 0.243, 0.140,
 0.552, 0.334, 0.238, 0.435, 0.302, 0.308, 0.378, 0.418, 0.344, 0.427, 0.313, 0.444, 0.284,
@@ -61,4 +60,5 @@ hcal_pos_gain_cor= 0.217, 0.133, 0.135, 0.238, 0.168, 0.147, 0.146, 0.281, 0.154
 hcal_neg_gain_cor= 0.434, 0.195, 0.196, 0.215, 0.138, 0.189, 0.294, 0.433, 0.214, 0.300, 0.162, 0.151, 0.202,
 0.265, 0.364, 0.235, 0.317, 0.279, 0.268, 0.276, 0.229, 0.257, 0.279, 0.296, 0.223, 0.228,
 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
-0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
\ No newline at end of file
+0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+