Add HMS calorimeter calibration.

    Works now with Root6 and KPP data.

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 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 const Double_t fAcor = 200.;
+  ////  static const Double_t fBcor = 8000.;
+  ////  static const Double_t fCcor = 64.36;
+  ////  static const Double_t fDcor = 1.66;
+  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 const Double_t fZbl = 10;   //cm, block transverse size
+  ////  static const UInt_t fNrows = 13;
+  ////  static const UInt_t fNcols =  4;
+  ////static const UInt_t fNnegs = 26;  // number of blocks with neg. side PMTs.
+  ////  static const UInt_t fNpmts = 78;   // total number of PMTs.
+  ////  static const UInt_t fNblks = fNrows*fNcols;
+  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.;
+}
+
+//------------------------------------------------------------------------------
+
+//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/THcShowerCalib.h

+#ifndef ROOT_THcShowerCalib
+#define ROOT_THcShowerCalib
+
+#include "THcShTrack.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 338.69    //distance from FP to the calorimeter face
+
+//Whole calorimeter
+#define XMIN -65.4
+#define XMAX  54.6
+#define YMIN -30.
+#define YMAX  30.
+
+#define DELTA_MIN -10   //HMS nominal acceptance
+#define DELTA_MAX  10
+
+#define PR_ADC_THR 0
+#define SH_ADC_THR 0
+
+#define BETA_MIN 0.5
+#define BETA_MAX 1.5
+
+using namespace std;
+
+//
+// HMS Shower Counter calibration class.
+//
+
+class THcShowerCalib {
+
+ public:
+  ////  THcShowerCalib(Int_t);
+  THcShowerCalib(string);
+  THcShowerCalib();
+  ~THcShowerCalib();
+
+  void Init();
+  bool ReadShRawTrack(THcShTrack &trk, UInt_t ientry);
+  void CalcThresholds();
+  void ComposeVMs();
+  void SolveAlphas();
+  void FillHEcal();
+  void SaveAlphas();
+  void SaveRawData();
+
+  TH1F* hEunc;
+  TH1F* hEuncSel;
+  TH1F* hEcal;
+  TH2F* hDPvsEcal;
+
+ 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 = 100;   // 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[THcShTrack::fNpmts];
+  Double_t fq0[THcShTrack::fNpmts];
+  Double_t fQ[THcShTrack::fNpmts][THcShTrack::fNpmts];
+  Double_t falphaU[THcShTrack::fNpmts];   // 'unconstrained' calib. constants
+  Double_t falphaC[THcShTrack::fNpmts];   // the sought calibration constants
+  Double_t falpha0[THcShTrack::fNpmts];   // initial gains
+  Double_t falpha1[THcShTrack::fNpmts];   // unit gains
+
+  UInt_t fHitCount[THcShTrack::fNpmts];
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShowerCalib::THcShowerCalib() {};
+
+//------------------------------------------------------------------------------
+
+////THcShowerCalib::THcShowerCalib(Int_t RunNumber) {
+THcShowerCalib::THcShowerCalib(string RunNumber) {
+  fRunNumber = RunNumber;
+};
+
+//------------------------------------------------------------------------------
+
+THcShowerCalib::~THcShowerCalib() {
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::SaveRawData() {
+
+  // Output raw data into file for debug purposes. To be called after
+  // calibration constants are determined.
+
+  cout << "SaveRawData: Output raw data into hcal_calib.raw_data." << endl;
+
+  ofstream fout;
+  fout.open("hcal_calib.raw_data",ios::out);
+
+  THcShTrack trk;
+
+  for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+    if (ReadShRawTrack(trk, ientry)) {
+      trk.SetEs(falphaC);
+      trk.Print(fout);
+    }
+
+  }
+
+  fout.close();
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::Init() {
+
+  //Reset ROOT and connect tree file.
+
+  gROOT->Reset();
+
+  ////  char* fname = Form("Root_files/hcal_calib_%d.root",fRunNumber);
+  char* fname = Form("ROOTfiles/hms_replay_%s.root",fRunNumber.c_str());
+  cout << "THcShowerCalib::Init: Root file name = " << fname << endl;
+
+  TFile *f = new TFile(fname);
+  f->GetObject("T",fTree);
+
+  fNentries = fTree->GetEntries();
+  cout << "THcShowerCalib::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 ",
+		       150,0.,1.5, 250,-12.5,12.5);
+
+  // Initialize qumulative quantities.
+  
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++) fHitCount[i] = 0;
+
+  fe0 = 0.;
+
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++) {
+    fqe[i] = 0.;
+    fq0[i] = 0.;
+    falphaU[i] = 0.;
+    falphaC[i] = 0.;
+    for (UInt_t j=0; j<THcShTrack::fNpmts; j++) {
+      fQ[i][j] = 0.;
+    }
+  }
+
+  // Initial gains (0.5 for the 2 first columns, 1 for others).
+
+  for (UInt_t iblk=0; iblk<THcShTrack::fNblks; iblk++) {
+    if (iblk < THcShTrack::fNnegs) {
+      falpha0[iblk] = 0.5;
+      falpha0[THcShTrack::fNblks+iblk] = 0.5;
+    }
+    else {
+      falpha0[iblk] = 1.;
+    }
+  };
+
+  // Unit gains.
+
+  for (UInt_t ipmt=0; ipmt<THcShTrack::fNpmts; ipmt++) {
+    falpha1[ipmt] = 1.;
+  }
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::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.
+
+  Int_t nev = 0;
+  THcShTrack 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++;
+      //    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 selected wthin 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 THcShowerCalib::ReadShRawTrack(THcShTrack &trk, UInt_t ientry) {
+
+  //
+  // Set a Shower track event from ntuple ientry.
+  //
+
+  // Declaration of leaves types
+
+  // Calorimeter ADC signals.
+
+  Double_t        H_cal_1pr_aneg_p[THcShTrack::fNrows];
+  Double_t        H_cal_1pr_apos_p[THcShTrack::fNrows];
+
+  Double_t        H_cal_2ta_aneg_p[THcShTrack::fNrows];
+  Double_t        H_cal_2ta_apos_p[THcShTrack::fNrows];
+
+  Double_t        H_cal_3ta_aneg_p[THcShTrack::fNrows];
+  Double_t        H_cal_3ta_apos_p[THcShTrack::fNrows];
+
+  Double_t        H_cal_4ta_aneg_p[THcShTrack::fNrows];
+  Double_t        H_cal_4ta_apos_p[THcShTrack::fNrows];
+
+  // Track parameters.
+
+  double          H_tr_n;
+  Double_t        H_tr_p;
+  Double_t        H_tr_x;   //X FP
+  Double_t        H_tr_xp;
+  Double_t        H_tr_y;   //Y FP
+  Double_t        H_tr_yp;
+
+  Double_t        H_tr_tg_dp;
+
+  Double_t        H_cer_npe[2];
+  Double_t        H_tr_beta;
+ 
+  // Set branch addresses.
+
+  fTree->SetBranchAddress("H.cal.1pr.aneg_p",H_cal_1pr_aneg_p);
+  fTree->SetBranchAddress("H.cal.1pr.apos_p",H_cal_1pr_apos_p);
+
+  fTree->SetBranchAddress("H.cal.2ta.aneg_p",H_cal_2ta_aneg_p);
+  fTree->SetBranchAddress("H.cal.2ta.apos_p",H_cal_2ta_apos_p);
+
+  fTree->SetBranchAddress("H.cal.3ta.aneg_p",H_cal_3ta_aneg_p);
+  fTree->SetBranchAddress("H.cal.3ta.apos_p",H_cal_3ta_apos_p);
+
+  fTree->SetBranchAddress("H.cal.4ta.aneg_p",H_cal_4ta_aneg_p);
+  fTree->SetBranchAddress("H.cal.4ta.apos_p",H_cal_4ta_apos_p);
+
+  fTree->SetBranchAddress("H.tr.n", &H_tr_n);
+  fTree->SetBranchAddress("H.tr.x",&H_tr_x);
+  fTree->SetBranchAddress("H.tr.y",&H_tr_y);
+  fTree->SetBranchAddress("H.tr.th",&H_tr_xp);
+  fTree->SetBranchAddress("H.tr.ph",&H_tr_yp);
+  fTree->SetBranchAddress("H.tr.p",&H_tr_p);
+
+  fTree->SetBranchAddress("H.tr.tg_dp", &H_tr_tg_dp);
+ 
+  fTree->SetBranchAddress("H.cer.npe", H_cer_npe);
+  fTree->SetBranchAddress("H.tr.beta", &H_tr_beta);
+
+  fTree->GetEntry(ientry);
+
+  if (ientry%100000 == 0) cout << "   ReadShRawTrack: " << ientry << endl;
+
+  if (H_tr_n != 1) return 0;
+
+  bool good_trk =   H_tr_tg_dp > DELTA_MIN &&
+		    H_tr_tg_dp < DELTA_MAX &&
+		    H_tr_x + H_tr_xp*D_CALO_FP > XMIN &&
+		    H_tr_x + H_tr_xp*D_CALO_FP < XMAX &&
+                    H_tr_y + H_tr_yp*D_CALO_FP > YMIN &&
+		    H_tr_y + H_tr_yp*D_CALO_FP < YMAX ;
+  if (!good_trk) return 0;
+
+  bool good_cer = H_cer_npe[0] > 0.125 ||
+                  H_cer_npe[1] > 0.125 ;
+  if(!good_cer) return 0;
+
+  bool good_beta = H_tr_beta > BETA_MIN &&
+                   H_tr_beta < BETA_MAX ;
+  if(!good_beta) return 0;
+
+  trk.Reset(H_tr_p, H_tr_tg_dp, H_tr_x+D_CALO_FP*H_tr_xp, H_tr_xp,
+	    H_tr_y+D_CALO_FP*H_tr_yp, H_tr_yp);
+
+  for (UInt_t j=0; j<THcShTrack::fNrows; j++) {
+    for (UInt_t k=0; k<THcShTrack::fNcols; k++) {
+
+      Double_t adc_pos, adc_neg;
+
+      switch (k) {
+      case 0 : 
+	adc_pos = H_cal_1pr_apos_p[j];
+	adc_neg = H_cal_1pr_aneg_p[j];
+	break;
+      case 1 : 
+	adc_pos = H_cal_2ta_apos_p[j];
+	adc_neg = H_cal_2ta_aneg_p[j];
+	break;
+      case 2 : 
+	adc_pos = H_cal_3ta_apos_p[j];
+	adc_neg = H_cal_3ta_aneg_p[j];
+	break;
+      case 3 : 
+	adc_pos = H_cal_4ta_apos_p[j];
+	adc_neg = H_cal_4ta_aneg_p[j];
+	break;
+      default:
+	cout << "*** ReadShRawTrack: column number k=" << k
+	     << " out of range! ***" << endl;
+      };
+
+      UInt_t nb = j+1 + k*THcShTrack::fNrows;
+
+      if (adc_pos>0. || adc_neg>0.) {
+	trk.AddHit(adc_pos, adc_neg, 0., 0., nb);
+      }
+
+    }
+  }
+
+  return 1;
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::ComposeVMs() {
+
+  //
+  // Fill in vectors and matrixes for the gain constant calculations.
+  //
+
+  fNev = 0;
+  THcShTrack 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++) {
+
+	  THcShHit* hit = trk.GetHit(i);
+	  // hit->Print(cout);
+	  
+	  UInt_t nb = hit->GetBlkNumber();
+
+	  // Fill the qe and q0 vectors (for positive side PMT).
+
+	  fqe[nb-1] += hit->GetEpos() * trk.GetP();
+	  fq0[nb-1] += hit->GetEpos();
+
+	  // Save the PMT hit.
+
+	  pmt_hit_list.push_back( pmt_hit{hit->GetEpos(), nb} );
+
+	  fHitCount[nb-1]++;   //Accrue the hit counter.
+
+	  // Do same for the negative side PMTs.
+
+	  if (nb <= THcShTrack::fNnegs) {
+	    fqe[THcShTrack::fNblks+nb-1] += hit->GetEneg() * trk.GetP();
+	    fq0[THcShTrack::fNblks+nb-1] += hit->GetEneg();
+
+	    pmt_hit_list.push_back(pmt_hit{hit->GetEneg(),
+		  THcShTrack::fNblks+nb} );
+
+	    fHitCount[THcShTrack::fNblks+nb-1]++;
+	  };
+
+	}      //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 the thresholds
+
+    };   // success in reading
+
+  };     // over entries
+
+  // Take averages.
+
+  fe0 /= fNev;
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++) {
+    fqe[i] /= fNev;
+    fq0[i] /= fNev;
+  }
+
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+    for (UInt_t j=0; j<THcShTrack::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<THcShTrack::fNpmts; i++)
+    q0out << fq0[i] << " " << i << endl;
+  q0out.close();
+
+  ofstream qeout;
+  qeout.open("qe.deb",ios::out);
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+    qeout << fqe[i] << " " << i << endl;
+  qeout.close();
+
+  ofstream Qout;
+  Qout.open("Q.deb",ios::out);
+  for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+    for (UInt_t j=0; j<THcShTrack::fNpmts; j++)
+      Qout << fQ[i][j] << " " << i << " " << j << endl;
+  Qout.close();
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::SolveAlphas() {
+
+  //
+  // Solve for the sought calibration constants, by use of the Root
+  // matrix algebra package.
+  //
+
+  TMatrixD Q(THcShTrack::fNpmts,THcShTrack::fNpmts);
+  TVectorD q0(THcShTrack::fNpmts);
+  TVectorD qe(THcShTrack::fNpmts);
+  TVectorD au(THcShTrack::fNpmts);
+  TVectorD ac(THcShTrack::fNpmts);
+  Bool_t ok;
+
+  cout << "Solving Alphas..." << endl;
+  cout << endl;
+
+  // Print out hit numbers.
+
+  cout << "Hit counts:" << endl;
+  UInt_t j = 0;
+  cout << "Positives:";
+  for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+    cout << setw(6) << fHitCount[j++] << ",";
+  cout << endl;
+  for (Int_t k=0; k<3; k++) {
+    cout << "          ";
+    for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+      cout << setw(6) << fHitCount[j++] << ",";
+    cout << endl;
+  }
+  cout << "Negatives:";
+  for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+    cout << setw(6) << fHitCount[j++] << ",";
+  cout << endl;
+  cout << "          ";
+  for (UInt_t i=0; i<THcShTrack::fNrows; 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<THcShTrack::fNpmts; i++) {
+    q0[i] = fq0[i];
+    qe[i] = fqe[i];
+    for (UInt_t k=0; k<THcShTrack::fNpmts; k++) {
+      Q[i][k] = fQ[i][k];
+    }
+  }
+
+  // Sanity check.
+
+  for (UInt_t i=0; i<THcShTrack::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<THcShTrack::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<THcShTrack::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<THcShTrack::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(THcShTrack::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<THcShTrack::fNpmts; i++) {
+    falphaU[i] = au[i];
+    falphaC[i] = ac[i];
+  }
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::FillHEcal() {
+
+  //
+  // Fill histogram of the normalized energy deposition, 2-d histogram
+  // of momentum deviation versus normalized energy deposition.
+  //
+
+  ofstream output;
+  output.open("calibrated.deb",ios::out);
+
+  Int_t nev = 0;
+
+  THcShTrack 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);
+
+      ////      Double_t delta;
+      ////      fTree->SetBranchAddress("H.tr.tg_dp",&delta);
+      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 THcShowerCalib::SaveAlphas() {
+
+  //
+  // Output the gain constants in a format suitable for inclusion in the
+  // hcal.param file to be used in the analysis.
+  //
+
+  ofstream output;
+  ////  char* fname = Form("hcal.param.%d",fRunNumber);
+  char* fname = Form("hcal.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;
+  output << "hcal_pos_gain_cor=";
+  for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+    output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+  output << endl;
+  for (Int_t k=0; k<3; k++) {
+    output << "                  ";
+    for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+      output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+    output << endl;
+  }
+  output << "hcal_neg_gain_cor=";
+  for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+    output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+  output << endl;
+  output << "                  ";
+  for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+    output << fixed << setw(6) << setprecision(3) << falphaC[j++] << ",";
+  output << endl;
+  for (Int_t k=0; k<2; k++) {
+    output << "                  ";
+    for (UInt_t i=0; i<THcShTrack::fNrows; i++)
+      output << fixed << setw(6) << setprecision(3) << 0. << ",";
+    output << endl;
+  }
+
+  output.close();
+}
+
+#endif

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");
+
+ // 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);
+}