diff --git a/CALIBRATION/hms_cal_calib/THcShHit.h b/CALIBRATION/hms_cal_calib/THcShHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..20afcced1f71e1fe0c234c87a60eff0db32298c5
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/THcShHit.h
@@ -0,0 +1,72 @@
+#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;};
diff --git a/CALIBRATION/hms_cal_calib/THcShTrack.h b/CALIBRATION/hms_cal_calib/THcShTrack.h
new file mode 100644
index 0000000000000000000000000000000000000000..ffd11b71c393acabca51c05003dc4ad87b7c18c6
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/THcShTrack.h
@@ -0,0 +1,230 @@
+#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);
+}
diff --git a/CALIBRATION/hms_cal_calib/THcShowerCalib.h b/CALIBRATION/hms_cal_calib/THcShowerCalib.h
new file mode 100644
index 0000000000000000000000000000000000000000..654f0c493e24a4e59b85fa1866ab3e054ca9c316
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/THcShowerCalib.h
@@ -0,0 +1,756 @@
+#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
diff --git a/CALIBRATION/hms_cal_calib/hcal_calib.cpp b/CALIBRATION/hms_cal_calib/hcal_calib.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b07474614024c1bd1ac5d737fd2705c3782f5dfd
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/hcal_calib.cpp
@@ -0,0 +1,61 @@
+#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);
+}