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..4d9a9d73e717f4a08c89133fbe3f0d9bac182b92
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/THcShTrack.h
@@ -0,0 +1,261 @@
+#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);
+}
diff --git a/CALIBRATION/hms_cal_calib/THcShowerCalib.h b/CALIBRATION/hms_cal_calib/THcShowerCalib.h
new file mode 100644
index 0000000000000000000000000000000000000000..523ea33c18433a774d1921d1d6e9777c14004ced
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/THcShowerCalib.h
@@ -0,0 +1,791 @@
+#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 BETA_MIN 0.5
+#define BETA_MAX 1.5
+
+#define CER_MIN 1.5
+
+#define MIN_HIT_COUNT 100 // Minimum number of hits for a PMT to be calibrated.
+
+using namespace std;
+
+//
+// HMS Shower Counter calibration class.
+//
+
+class THcShowerCalib {
+
+ public:
+ 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;
+ TH2F* hETAvsEPR;
+
+ 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
+
+ // 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;
+
+ Double_t H_cal_nclust;
+
+ TBranch* b_H_cal_1pr_aneg_p;
+ TBranch* b_H_cal_1pr_apos_p;
+
+ TBranch* b_H_cal_2ta_aneg_p;
+ TBranch* b_H_cal_2ta_apos_p;
+
+ TBranch* b_H_cal_3ta_aneg_p;
+ TBranch* b_H_cal_3ta_apos_p;
+
+ TBranch* b_H_cal_4ta_aneg_p;
+ TBranch* b_H_cal_4ta_apos_p;
+
+ TBranch* b_H_tr_n;
+ TBranch* b_H_tr_x;
+ TBranch* b_H_tr_y;
+ TBranch* b_H_tr_xp;
+ TBranch* b_H_tr_yp;
+ TBranch* b_H_tr_p;
+
+ TBranch* b_H_tr_tg_dp;
+
+ TBranch* b_H_cer_npe;
+ TBranch* b_H_tr_beta;
+
+ TBranch* b_H_cal_nclust;
+
+ // 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(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("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;
+
+ // Set branch addresses.
+
+ fTree->SetBranchAddress("H.cal.1pr.aneg_p",H_cal_1pr_aneg_p,
+ &b_H_cal_1pr_aneg_p);
+ fTree->SetBranchAddress("H.cal.1pr.apos_p",H_cal_1pr_apos_p,
+ &b_H_cal_1pr_apos_p);
+
+ fTree->SetBranchAddress("H.cal.2ta.aneg_p",H_cal_2ta_aneg_p,
+ &b_H_cal_2ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.2ta.apos_p",H_cal_2ta_apos_p,
+ &b_H_cal_2ta_apos_p);
+
+ fTree->SetBranchAddress("H.cal.3ta.aneg_p",H_cal_3ta_aneg_p,
+ &b_H_cal_3ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.3ta.apos_p",H_cal_3ta_apos_p,
+ &b_H_cal_3ta_apos_p);
+
+ fTree->SetBranchAddress("H.cal.4ta.aneg_p",H_cal_4ta_aneg_p,
+ &b_H_cal_4ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.4ta.apos_p",H_cal_4ta_apos_p,
+ &b_H_cal_4ta_apos_p);
+
+ fTree->SetBranchAddress("H.tr.n", &H_tr_n,&b_H_tr_n);
+ fTree->SetBranchAddress("H.tr.x",&H_tr_x,&b_H_tr_x);
+ fTree->SetBranchAddress("H.tr.y",&H_tr_y,&b_H_tr_y);
+ fTree->SetBranchAddress("H.tr.th",&H_tr_xp,&b_H_tr_xp);
+ fTree->SetBranchAddress("H.tr.ph",&H_tr_yp,&b_H_tr_yp);
+ fTree->SetBranchAddress("H.tr.p",&H_tr_p,&b_H_tr_p);
+
+ fTree->SetBranchAddress("H.tr.tg_dp", &H_tr_tg_dp,&b_H_tr_tg_dp);
+
+ fTree->SetBranchAddress("H.cer.npe", H_cer_npe,&b_H_cer_npe);
+ fTree->SetBranchAddress("H.tr.beta", &H_tr_beta,&b_H_tr_beta);
+
+ fTree->SetBranchAddress("H.cal.nclust", &H_cal_nclust,&b_H_cal_nclust);
+
+ // 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);
+ hETAvsEPR = new TH2F("hETAvsEPR", "E_{TA} versus E_{PR}",
+ 300,0.,1.5, 300,0.,1.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.
+ //
+
+ fTree->GetEntry(ientry);
+
+ if (ientry%100000 == 0) cout << " ReadShRawTrack: " << ientry << endl;
+
+ if (H_tr_n != 1) return 0;
+
+ if (H_cal_nclust != 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] > CER_MIN ||
+ H_cer_npe[1] > CER_MIN ;
+ 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 " << MIN_HIT_COUNT
+ << " will not be calibrated." << endl;
+ cout << endl;
+
+ for (UInt_t i=0; i<THcShTrack::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<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);
+
+ hDPvsEcal->Fill(Enorm,delta,1.);
+
+ hETAvsEPR->Fill(trk.EPRnorm(), trk.ETAnorm());
+
+ 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.%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..6594d8c0d4a77ee89ea814974de4d1a2b5b1a7f0
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/hcal_calib.cpp
@@ -0,0 +1,66 @@
+#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);
+}
diff --git a/CALIBRATION/hms_cal_calib/howto.txt b/CALIBRATION/hms_cal_calib/howto.txt
new file mode 100644
index 0000000000000000000000000000000000000000..764d780cda0d1da05adae45a1cdc0a473c2c4860
--- /dev/null
+++ b/CALIBRATION/hms_cal_calib/howto.txt
@@ -0,0 +1,31 @@
+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.
diff --git a/CALIBRATION/shms_cal_calib/ROOTfiles b/CALIBRATION/shms_cal_calib/ROOTfiles
new file mode 120000
index 0000000000000000000000000000000000000000..431cd2c92c0de8a3bae4131fc34816f3af0c523b
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/ROOTfiles
@@ -0,0 +1 @@
+/net/cdaq/cdaql3data/cdaq/hallc-online/ROOTfiles
\ No newline at end of file
diff --git a/CALIBRATION/shms_cal_calib/THcPShHit.h b/CALIBRATION/shms_cal_calib/THcPShHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..06eb6067b3296317a6c86c05082e77bfca39b69a
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShHit.h
@@ -0,0 +1,59 @@
+#include <iostream>
+
+// SHMS calorimeter hit class for calibration.
+
+class THcPShHit {
+
+ Double_t ADC; // pedestal subtracted ADC signal.
+ Double_t Edep; // Energy deposition.
+ UInt_t BlkNumber; // Block number.
+
+ public:
+
+ THcPShHit();
+ THcPShHit(Double_t adc, UInt_t blk_number);
+ ~THcPShHit();
+
+ void SetADC(Double_t sig) {ADC = sig;}
+
+ void SetEdep(Double_t e) {Edep = e;}
+
+ void SetBlkNumber(UInt_t n) {BlkNumber = n;}
+
+ Double_t GetADC() {return ADC;}
+
+ Double_t GetEdep() {return Edep;}
+
+ UInt_t GetBlkNumber() {return BlkNumber;}
+
+ void Print(ostream & ostrm);
+};
+
+//------------------------------------------------------------------------------
+
+THcPShHit::THcPShHit() {
+ ADC = -99999.;
+ Edep = -99999.;
+ BlkNumber = 99999;
+};
+
+THcPShHit::THcPShHit(Double_t adc, UInt_t blk_number) {
+ ADC = adc;
+ Edep = 0.;
+ BlkNumber = blk_number;
+};
+
+THcPShHit::~THcPShHit() { };
+
+//------------------------------------------------------------------------------
+
+void THcPShHit::Print(ostream & ostrm) {
+
+ // Output hit data through the stream ostrm.
+
+ ostrm << ADC << " " << Edep << " " << BlkNumber << endl;
+};
+
+//------------------------------------------------------------------------------
+
+struct pmt_hit {Double_t signal; UInt_t channel;};
diff --git a/CALIBRATION/shms_cal_calib/THcPShTrack.h b/CALIBRATION/shms_cal_calib/THcPShTrack.h
new file mode 100644
index 0000000000000000000000000000000000000000..1592e5cd57e2fb5d23a4521538fff0d8546bdefd
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShTrack.h
@@ -0,0 +1,251 @@
+#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;
+}
diff --git a/CALIBRATION/shms_cal_calib/THcPShowerCalib.h b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
new file mode 100644
index 0000000000000000000000000000000000000000..bb4613517a37970d16b98cdb7330c35907dacf6c
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
@@ -0,0 +1,737 @@
+#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
diff --git a/CALIBRATION/shms_cal_calib/howto.txt b/CALIBRATION/shms_cal_calib/howto.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1335b2ba03a026d42e22d4b6e3660e2fb0f21710
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/howto.txt
@@ -0,0 +1,31 @@
+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.
diff --git a/CALIBRATION/shms_cal_calib/pcal_calib.cpp b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..df917354358bacf3a139690d53b0338acbc6bdf6
--- /dev/null
+++ b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
@@ -0,0 +1,64 @@
+#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);
+}
diff --git a/PARAM/HMS/CAL/hcal.param b/PARAM/HMS/CAL/hcal.param
index 78b490f157d7f69d8739ecb0d51102995784cefe..88cc67b47ef420a4279d9361e6117aa6895d401c 100644
--- a/PARAM/HMS/CAL/hcal.param
+++ b/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,
+