From 4d1f82dcb5b278e14dfb2f9c0ca456f64cd9acbb Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Fri, 19 Jun 2015 17:12:02 +0400
Subject: [PATCH] Add SHMS calibration code and move HMS calibration code
Create hc_cal_calib directory for the calorimeter calibrations. In the
hc_cal_calib, the HMS calorimeter calibration resides in hcal_calib
directory, and the SHMS calorimeter calibration in shms_cal_calib
directory.
Change file extension for debug files from .d to .deb to avoid
confusion with make dependency files.
---
.../hcal_calib}/THcShHit.h | 0
.../hcal_calib}/THcShTrack.h | 0
.../hcal_calib}/THcShowerCalib.h | 8 +-
.../hcal_calib}/db_run.dat | 0
.../hcal_calib}/hcal_calib.cpp | 0
.../hcal_calib}/hcal_replay.cpp | 0
.../hcal_calib}/hcal_replay_cuts.def | 0
.../hcal_calib}/instructions.txt | 0
hc_cal_calib/hcal_calib/make_cratemap.pl | 1 +
.../hcal_calib}/output_hcal_replay.def | 0
hc_cal_calib/shms_cal_calib/THcPShHit.h | 59 ++
hc_cal_calib/shms_cal_calib/THcPShTrack.h | 211 ++++++
hc_cal_calib/shms_cal_calib/THcPShowerCalib.h | 645 ++++++++++++++++++
hc_cal_calib/shms_cal_calib/instructions.txt | 31 +
.../shms_cal_calib/mc_results/make_root.C | 100 +++
.../shms_cal_calib/mc_results/memo.txt | 4 +
.../new_old_diff/comp_pcal_etot.C | 127 ++++
.../shms_cal_calib/old_code/calib.kumac | 87 +++
.../shms_cal_calib/old_code/s_correct.f | 77 +++
.../shms_cal_calib/old_code/shms_calib.f | 509 ++++++++++++++
.../shms_cal_calib/old_params/11.param | 18 +
.../shms_cal_calib/old_params/5.param | 18 +
hc_cal_calib/shms_cal_calib/pcal_calib.cpp | 47 ++
hcal_calib/make_cratemap.pl | 1 -
24 files changed, 1938 insertions(+), 5 deletions(-)
rename {hcal_calib => hc_cal_calib/hcal_calib}/THcShHit.h (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/THcShTrack.h (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/THcShowerCalib.h (99%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/db_run.dat (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/hcal_calib.cpp (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/hcal_replay.cpp (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/hcal_replay_cuts.def (100%)
rename {hcal_calib => hc_cal_calib/hcal_calib}/instructions.txt (100%)
create mode 120000 hc_cal_calib/hcal_calib/make_cratemap.pl
rename {hcal_calib => hc_cal_calib/hcal_calib}/output_hcal_replay.def (100%)
create mode 100644 hc_cal_calib/shms_cal_calib/THcPShHit.h
create mode 100644 hc_cal_calib/shms_cal_calib/THcPShTrack.h
create mode 100644 hc_cal_calib/shms_cal_calib/THcPShowerCalib.h
create mode 100644 hc_cal_calib/shms_cal_calib/instructions.txt
create mode 100644 hc_cal_calib/shms_cal_calib/mc_results/make_root.C
create mode 100644 hc_cal_calib/shms_cal_calib/mc_results/memo.txt
create mode 100644 hc_cal_calib/shms_cal_calib/new_old_diff/comp_pcal_etot.C
create mode 100644 hc_cal_calib/shms_cal_calib/old_code/calib.kumac
create mode 100644 hc_cal_calib/shms_cal_calib/old_code/s_correct.f
create mode 100644 hc_cal_calib/shms_cal_calib/old_code/shms_calib.f
create mode 100644 hc_cal_calib/shms_cal_calib/old_params/11.param
create mode 100644 hc_cal_calib/shms_cal_calib/old_params/5.param
create mode 100644 hc_cal_calib/shms_cal_calib/pcal_calib.cpp
delete mode 120000 hcal_calib/make_cratemap.pl
diff --git a/hcal_calib/THcShHit.h b/hc_cal_calib/hcal_calib/THcShHit.h
similarity index 100%
rename from hcal_calib/THcShHit.h
rename to hc_cal_calib/hcal_calib/THcShHit.h
diff --git a/hcal_calib/THcShTrack.h b/hc_cal_calib/hcal_calib/THcShTrack.h
similarity index 100%
rename from hcal_calib/THcShTrack.h
rename to hc_cal_calib/hcal_calib/THcShTrack.h
diff --git a/hcal_calib/THcShowerCalib.h b/hc_cal_calib/hcal_calib/THcShowerCalib.h
similarity index 99%
rename from hcal_calib/THcShowerCalib.h
rename to hc_cal_calib/hcal_calib/THcShowerCalib.h
index c839771..b920e2b 100644
--- a/hcal_calib/THcShowerCalib.h
+++ b/hc_cal_calib/hcal_calib/THcShowerCalib.h
@@ -427,19 +427,19 @@ void THcShowerCalib::ComposeVMs() {
// Output vectors and matrixes, for debug purposes.
ofstream q0out;
- q0out.open("q0.d",ios::out);
+ 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.d",ios::out);
+ 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.d",ios::out);
+ 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;
@@ -608,7 +608,7 @@ void THcShowerCalib::FillHEcal() {
//
ofstream output;
- output.open("calibrated.d",ios::out);
+ output.open("calibrated.deb",ios::out);
Int_t nev = 0;
diff --git a/hcal_calib/db_run.dat b/hc_cal_calib/hcal_calib/db_run.dat
similarity index 100%
rename from hcal_calib/db_run.dat
rename to hc_cal_calib/hcal_calib/db_run.dat
diff --git a/hcal_calib/hcal_calib.cpp b/hc_cal_calib/hcal_calib/hcal_calib.cpp
similarity index 100%
rename from hcal_calib/hcal_calib.cpp
rename to hc_cal_calib/hcal_calib/hcal_calib.cpp
diff --git a/hcal_calib/hcal_replay.cpp b/hc_cal_calib/hcal_calib/hcal_replay.cpp
similarity index 100%
rename from hcal_calib/hcal_replay.cpp
rename to hc_cal_calib/hcal_calib/hcal_replay.cpp
diff --git a/hcal_calib/hcal_replay_cuts.def b/hc_cal_calib/hcal_calib/hcal_replay_cuts.def
similarity index 100%
rename from hcal_calib/hcal_replay_cuts.def
rename to hc_cal_calib/hcal_calib/hcal_replay_cuts.def
diff --git a/hcal_calib/instructions.txt b/hc_cal_calib/hcal_calib/instructions.txt
similarity index 100%
rename from hcal_calib/instructions.txt
rename to hc_cal_calib/hcal_calib/instructions.txt
diff --git a/hc_cal_calib/hcal_calib/make_cratemap.pl b/hc_cal_calib/hcal_calib/make_cratemap.pl
new file mode 120000
index 0000000..8596465
--- /dev/null
+++ b/hc_cal_calib/hcal_calib/make_cratemap.pl
@@ -0,0 +1 @@
+../../examples/make_cratemap.pl
\ No newline at end of file
diff --git a/hcal_calib/output_hcal_replay.def b/hc_cal_calib/hcal_calib/output_hcal_replay.def
similarity index 100%
rename from hcal_calib/output_hcal_replay.def
rename to hc_cal_calib/hcal_calib/output_hcal_replay.def
diff --git a/hc_cal_calib/shms_cal_calib/THcPShHit.h b/hc_cal_calib/shms_cal_calib/THcPShHit.h
new file mode 100644
index 0000000..06eb606
--- /dev/null
+++ b/hc_cal_calib/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/hc_cal_calib/shms_cal_calib/THcPShTrack.h b/hc_cal_calib/shms_cal_calib/THcPShTrack.h
new file mode 100644
index 0000000..e6200f5
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/THcPShTrack.h
@@ -0,0 +1,211 @@
+#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 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 x, Double_t xp, Double_t y, Double_t yp);
+ ~THcPShTrack();
+
+ void Reset(Double_t p, 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 GetP() {return P*1000.;} //MeV
+
+ Double_t GetX() {return X;}
+ Double_t GetY() {return Y;}
+
+ Float_t Ycor(Double_t, UInt_t); // coord. corection for Preshower module
+
+ // Coordinate correction constants for Preshower blocks
+ //
+ static const Double_t fAcor = 106.73;
+ static const Double_t fBcor = 2.329;
+
+ // 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 x, Double_t xp, Double_t y, Double_t yp) {
+ P = p;
+ X = x;
+ Xp = xp;
+ Y = y;
+ Yp =yp;
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShTrack::Reset(Double_t p,
+ Double_t x, Double_t xp, Double_t y, Double_t yp) {
+
+ // Reset track parameters, clear hit list.
+
+ P = p;
+ 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 << " " << 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.
+}
+
+//------------------------------------------------------------------------------
+
+// Coordinate correction for Preshower modules.
+// Fit to GEANT pion data @ 5 GeV/c (Simon).
+
+Float_t THcPShTrack::Ycor(Double_t yhit, UInt_t ncol) {
+
+ Float_t cor;
+
+ // Warn if hit does not belong to Preshower.
+ //
+ if (ncol > fNcols_pr || ncol < 1)
+ cout << "*** THcPShTrack::Ycor: wrong ncol = " << ncol << " ***" << endl;
+
+ // Check if the hit coordinate matches the fired block's column.
+ //
+ //if ((yhit < 0. && ncol == 1) || (yhit > 0. && ncol == 2)) //Vardan's MC data
+ if ((yhit < 0. && ncol == 2) || (yhit > 0. && ncol == 1)) //Simon's MC data
+ cor = 1./(1. + TMath::Power(TMath::Abs(yhit)/fAcor, fBcor));
+ else
+ cor = 1.;
+
+ // Debug output.
+ // cout << "THcShTrack::Ycor = " << cor << " yhit = " << yhit
+ // << " ncol = " << ncol << endl;
+
+ return cor;
+}
diff --git a/hc_cal_calib/shms_cal_calib/THcPShowerCalib.h b/hc_cal_calib/shms_cal_calib/THcPShowerCalib.h
new file mode 100644
index 0000000..ab3e0cf
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/THcPShowerCalib.h
@@ -0,0 +1,645 @@
+#ifndef ROOT_THcPShowerCalib
+#define ROOT_THcPShowerCalib
+
+#include "THcPShTrack.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TDecompLU.h"
+#include "TMath.h"
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+
+#include "TROOT.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#define D_CALO_FP 275. //distance from FP to the Preshower
+
+using namespace std;
+
+//
+// SHMS Calorimeter calibration class.
+//
+
+class THcPShowerCalib {
+
+ public:
+
+ THcPShowerCalib(Int_t);
+ THcPShowerCalib();
+ ~THcPShowerCalib();
+
+ void Init();
+ void 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;
+
+ private:
+
+ Int_t fRunNumber;
+ Double_t fLoThr; // Low and high thresholds on the normalized uncalibrated
+ Double_t fHiThr; // energy deposition.
+ UInt_t fNev; // Number of processed events.
+
+ static const UInt_t fMinHitCount = 5; // Minimum number of hits for a PMT
+ // to be calibrated.
+
+ TTree* fTree;
+ UInt_t fNentries;
+
+ // Quantities for calculations of the calibration constants.
+
+ Double_t fe0;
+ Double_t fqe[THcPShTrack::fNpmts];
+ Double_t fq0[THcPShTrack::fNpmts];
+ Double_t fQ[THcPShTrack::fNpmts][THcPShTrack::fNpmts];
+ Double_t falphaU[THcPShTrack::fNpmts]; // 'unconstrained' calib. constants
+ Double_t falphaC[THcPShTrack::fNpmts]; // the sought calibration constants
+ Double_t falpha0[THcPShTrack::fNpmts]; // initial gains
+ Double_t falpha1[THcPShTrack::fNpmts]; // unit gains
+
+ UInt_t fHitCount[THcPShTrack::fNpmts];
+
+};
+
+//------------------------------------------------------------------------------
+
+THcPShowerCalib::THcPShowerCalib() {};
+
+//------------------------------------------------------------------------------
+
+THcPShowerCalib::THcPShowerCalib(Int_t RunNumber) {
+ 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++) {
+ ReadShRawTrack(trk, ientry);
+ trk.SetEs(falphaC);
+ trk.Print(fout);
+ }
+
+ fout.close();
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::Init() {
+
+ //Reset ROOT and connect tree file.
+
+ gROOT->Reset();
+
+ char* fname = Form("Root_files/Pcal_calib_%d.root",fRunNumber);
+ cout << "THcPShowerCalib::Init: Root file name = " << fname << endl;
+
+ TFile *f = new TFile(fname);
+ f->GetObject("T",fTree);
+
+ fNentries = fTree->GetEntries();
+ cout << "THcPShowerCalib::Init: fNentries= " << fNentries << endl;
+
+ // Histogram declarations.
+
+ hEunc = new TH1F("hEunc", "Edep/P uncalibrated", 500, 0., 5.);
+ hEcal = new TH1F("hEcal", "Edep/P calibrated", 150, 0., 1.5);
+ hDPvsEcal = new TH2F("hDPvsEcal", "#DeltaP versus Edep/P ",
+ 350,0.,1.5, 250,-12.5,22.5);
+
+ // 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.
+
+ Int_t nev = 0;
+ THcPShTrack trk;
+
+ for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+ ReadShRawTrack(trk, ientry);
+
+ // trk.Print(cout);
+ // getchar();
+
+ trk.SetEs(falpha0); //Use initial gain constants here.
+ Double_t Enorm = trk.Enorm();
+
+ nev++;
+ // cout << "CalcThreshods: nev=" << nev << " Enorm=" << Enorm << endl;
+
+ hEunc->Fill(Enorm);
+ };
+
+ Double_t mean = hEunc->GetMean();
+ Double_t rms = hEunc->GetRMS();
+ cout << "CalcThreshods: mean=" << mean << " rms=" << rms << endl;
+
+ fLoThr = mean - 3.*rms;
+ fHiThr = mean + 3.*rms;
+ // fLoThr = 0.; // Wide open thrsholds, for
+ // fHiThr = 1.e+8; // comparison with the old code.
+
+ cout << "CalcThreshods: fLoThr=" << fLoThr << " fHiThr=" << fHiThr
+ << " nev=" << nev << endl;
+
+ Int_t nbins = hEunc->GetNbinsX();
+ Int_t nlo = hEunc->FindBin(fLoThr);
+ Int_t nhi = hEunc->FindBin(fHiThr);
+
+ cout << "CalcThresholds: nlo=" << nlo << " nhi=" << nhi
+ << " nbins=" << nbins << endl;
+
+ // Histogram of selected within the thresholds events.
+
+ hEuncSel = (TH1F*)hEunc->Clone("hEuncSel");
+
+ for (Int_t i=0; i<nlo; i++) hEuncSel->SetBinContent(i, 0.);
+ for (Int_t i=nhi; i<nbins+1; i++) hEuncSel->SetBinContent(i, 0.);
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::ReadShRawTrack(THcPShTrack &trk, UInt_t ientry) {
+
+ //
+ // Set a Shower track event from ntuple ientry.
+ //
+
+ // Declaration of leaves types
+
+ // Preshower and Shower ADC signals.
+
+ Double_t P_pr_a_p[THcPShTrack::fNrows_pr][THcPShTrack::fNcols_pr];
+ Double_t P_sh_a_p[THcPShTrack::fNrows_sh][THcPShTrack::fNcols_sh];
+
+ // Track parameters.
+
+ 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;
+
+ const Double_t adc_thr = 15.; //Low threshold on the ADC signals.
+
+ // Set branch addresses.
+
+ fTree->SetBranchAddress("P.pr.a_p",P_pr_a_p);
+ fTree->SetBranchAddress("P.sh.a_p",P_sh_a_p);
+
+ fTree->SetBranchAddress("P.tr.x", &P_tr_x);
+ fTree->SetBranchAddress("P.tr.y", &P_tr_y);
+ fTree->SetBranchAddress("P.tr.th",&P_tr_xp);
+ fTree->SetBranchAddress("P.tr.ph",&P_tr_yp);
+ fTree->SetBranchAddress("P.tr.p", &P_tr_p);
+
+ fTree->GetEntry(ientry);
+
+ // Set track coordinates and slopes at the face of Preshower.
+
+ trk.Reset(P_tr_p, 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.
+
+ for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+ for (UInt_t j=0; j<THcPShTrack::fNrows_pr; j++) {
+
+ Double_t adc = P_pr_a_p[j][k];
+
+ if (adc > adc_thr) {
+ UInt_t nb = j+1 + k*THcPShTrack::fNrows_pr;
+ trk.AddHit(adc, 0., nb);
+ }
+
+ }
+ }
+
+ // Set Shower hits.
+
+ for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+ for (UInt_t j=0; j<THcPShTrack::fNrows_sh; j++) {
+
+ Double_t adc = P_sh_a_p[j][k];
+
+ if (adc > adc_thr) {
+ UInt_t nb = THcPShTrack::fNpmts_pr + j+1 + k*THcPShTrack::fNrows_sh;
+ trk.AddHit(adc, 0., nb);
+ }
+
+ }
+ }
+
+}
+
+//------------------------------------------------------------------------------
+
+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++) {
+
+ 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 the thresholds
+
+ }; // over entries
+
+ // Take averages.
+
+ fe0 /= fNev;
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+ fqe[i] /= fNev;
+ fq0[i] /= fNev;
+ }
+
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+ for (UInt_t j=0; j<THcPShTrack::fNpmts; j++)
+ fQ[i][j] /= fNev;
+
+ // Output vectors and matrixes, for debug purposes.
+
+ ofstream q0out;
+ q0out.open("q0.deb",ios::out);
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+ q0out << setprecision(20) << fq0[i] << " " << i << endl;
+ q0out.close();
+
+ ofstream qeout;
+ qeout.open("qe.deb",ios::out);
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+ qeout << setprecision(20) << fqe[i] << " " << i << endl;
+ qeout.close();
+
+ ofstream Qout;
+ Qout.open("Q.deb",ios::out);
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++)
+ for (UInt_t j=0; j<THcPShTrack::fNpmts; j++)
+ Qout << setprecision(20) << fQ[i][j] << " " << i << " " << j << endl;
+ Qout.close();
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::SolveAlphas() {
+
+ //
+ // Solve for the sought calibration constants, by use of the Root
+ // matrix algebra package.
+ //
+
+ TMatrixD Q(THcPShTrack::fNpmts,THcPShTrack::fNpmts);
+ TVectorD q0(THcPShTrack::fNpmts);
+ TVectorD qe(THcPShTrack::fNpmts);
+ TVectorD au(THcPShTrack::fNpmts);
+ TVectorD ac(THcPShTrack::fNpmts);
+ Bool_t ok;
+
+ cout << "Solving Alphas..." << endl;
+ cout << endl;
+
+ // Print out hit numbers.
+
+ cout << "Hit counts:" << endl;
+ UInt_t j = 0;
+
+ for (UInt_t k=0; k<THcPShTrack::fNcols_pr; k++) {
+ k==0 ? cout << "Preshower:" : cout << " :";
+ for (UInt_t i=0; i<THcPShTrack::fNrows_pr; i++)
+ cout << setw(6) << fHitCount[j++] << ",";
+ cout << endl;
+ }
+
+ for (UInt_t k=0; k<THcPShTrack::fNcols_sh; k++) {
+ k==0 ? cout << "Shower :" : cout << " :";
+ for (UInt_t i=0; i<THcPShTrack::fNrows_sh; i++)
+ cout << setw(6) << fHitCount[j++] << ",";
+ cout << endl;
+ }
+
+ // Initialize the vectors and the matrix of the Root algebra package.
+
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+ q0[i] = fq0[i];
+ qe[i] = fqe[i];
+ for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+ Q[i][k] = fQ[i][k];
+ }
+ }
+
+ // Sanity check.
+
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+
+ // Check zero hit channels: the vector and matrix elements should be 0.
+
+ if (fHitCount[i] == 0) {
+
+ if (q0[i] != 0. || qe[i] != 0.) {
+
+ cout << "*** Inconsistency in chanel " << i << ": # of hits "
+ << fHitCount[i] << ", q0=" << q0[i] << ", qe=" << qe[i];
+
+ for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+ if (Q[i][k] !=0. || Q[k][i] !=0.)
+ cout << ", Q[" << i << "," << k << "]=" << Q[i][k]
+ << ", Q[" << k << "," << i << "]=" << Q[k][i];
+ }
+
+ cout << " ***" << endl;
+ }
+ }
+
+ // The hit channels: the vector elements should be non zero.
+
+ if ( (fHitCount[i] != 0) && (q0[i] == 0. || qe[i] == 0.) ) {
+ cout << "*** Inconsistency in chanel " << i << ": # of hits "
+ << fHitCount[i] << ", q0=" << q0[i] << ", qe=" << qe[i]
+ << " ***" << endl;
+ }
+
+ } //sanity check
+
+ // Low hit number channels: exclude from calculation. Assign all the
+ // correspondent elements 0, except self-correlation Q(i,i)=1.
+
+ cout << endl;
+ cout << "Channels with hit number less than " << fMinHitCount
+ << " will not be calibrated." << endl;
+ cout << endl;
+
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+
+ if (fHitCount[i] < fMinHitCount) {
+ cout << "Channel " << i << ", " << fHitCount[i]
+ << " hits, will not be calibrated." << endl;
+ q0[i] = 0.;
+ qe[i] = 0.;
+ for (UInt_t k=0; k<THcPShTrack::fNpmts; k++) {
+ Q[i][k] = 0.;
+ Q[k][i] = 0.;
+ }
+ Q[i][i] = 1.;
+ }
+
+ }
+
+ // Declare LU decomposition method for the correlation matrix Q.
+
+ TDecompLU lu(Q);
+ Double_t d1,d2;
+ lu.Det(d1,d2);
+ cout << "cond:" << lu.Condition() << endl;
+ cout << "det :" << d1*TMath::Power(2.,d2) << endl;
+ cout << "tol :" << lu.GetTol() << endl;
+
+ // Solve equation Q x au = qe for the 'unconstrained' calibration (gain)
+ // constants au.
+
+ au = lu.Solve(qe,ok);
+ cout << "au: ok=" << ok << endl;
+ // au.Print();
+
+ // Find the sought 'constrained' calibration constants next.
+
+ Double_t t1 = fe0 - au * q0; // temporary variable.
+ // cout << "t1 =" << t1 << endl;
+
+ TVectorD Qiq0(THcPShTrack::fNpmts); // an intermittent result
+ Qiq0 = lu.Solve(q0,ok);
+ cout << "Qiq0: ok=" << ok << endl;
+ // Qiq0.Print();
+
+ Double_t t2 = q0 * Qiq0; // another temporary variable
+ // cout << "t2 =" << t2 << endl;
+
+ ac = (t1/t2) *Qiq0 + au; // the sought gain constants
+ // cout << "ac:" << endl;
+ // ac.Print();
+
+ // Assign the gain arrays.
+
+ for (UInt_t i=0; i<THcPShTrack::fNpmts; i++) {
+ falphaU[i] = au[i];
+ falphaC[i] = ac[i];
+ }
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::FillHEcal() {
+
+ //
+ // Fill histogram of the normalized energy deposition, and 2-d histogram
+ // of momentum deviation versus normalized energy deposition.
+ // Output event by event energy depositions and momenta for debug purposes.
+ //
+
+ ofstream output;
+ output.open("calibrated.deb",ios::out);
+
+ Int_t nev = 0;
+
+ THcPShTrack trk;
+
+ for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+ ReadShRawTrack(trk, ientry);
+ // trk.Print(cout);
+
+ trk.SetEs(falphaC); // use the 'constrained' calibration constants
+ Double_t P = trk.GetP();
+ Double_t Enorm = trk.Enorm();
+
+ hEcal->Fill(Enorm);
+
+ Double_t delta;
+ fTree->SetBranchAddress("P.tr.tg_dp",&delta);
+ hDPvsEcal->Fill(Enorm,delta,1.);
+
+ output << Enorm*P/1000. << " " << P/1000. << " " << trk.GetX() << " "
+ << trk.GetY() << endl;
+
+ nev++;
+ };
+
+ output.close();
+
+ cout << "FillHEcal: " << nev << " events filled" << endl;
+};
+
+//------------------------------------------------------------------------------
+
+void THcPShowerCalib::SaveAlphas() {
+
+ //
+ // Output the gain constants in a format suitable for inclusion in the
+ // pcal.param file to be used in the analysis.
+ //
+
+ ofstream output;
+ char* fname = Form("pcal.param.%d",fRunNumber);
+ 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 << "shms_neg_pre_gain=" : output << "shms_pos_pre_gain=";
+ 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 << "shms_shower_gain =" : 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/hc_cal_calib/shms_cal_calib/instructions.txt b/hc_cal_calib/shms_cal_calib/instructions.txt
new file mode 100644
index 0000000..4538ba3
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/instructions.txt
@@ -0,0 +1,31 @@
+Instructions for calibration of the SHMS calorimeter.
+
+
+1. Currently the code works on the Monte Carlo simulated data. The
+simulated data are in plain text file, while the calibration code
+needs a root file as input. To prepare a root file from a text file,
+run make_root.C script in the mc_results directory. The script takes
+2 input parameters: prefix of the input file and SHMS momentum
+setting. For example, the command
+
+root [0] .x make_root.C+("5",5.)
+
+will produce a root file named 5.root from the input text file 5.r.
+
+
+2. Make Root_files directory, if not present. Copy (or link) the root
+file in the mc_results directory into the Root_files under the name
+Pcal_calib_<run_number>.root name. This file will serve as an input to
+the calibration code.
+
+
+3. Run pcal_calib.cpp Root script to calibrate a run. For instance,
+the Root command
+
+root [1] .x pcal_calib.cpp+(5)
+
+will produce pcal.param.5 text file containing the calibration
+constants, along with other text files for debug purposes. The
+graphics output will show distributions of the normalized energy
+depositions before and after the calibration, and deviation of
+momentum versus the normalized energy deposition on a scattered plot.
diff --git a/hc_cal_calib/shms_cal_calib/mc_results/make_root.C b/hc_cal_calib/shms_cal_calib/mc_results/make_root.C
new file mode 100644
index 0000000..ca2eb4f
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/mc_results/make_root.C
@@ -0,0 +1,100 @@
+#include <iostream>
+#include <fstream>
+#include <string>
+#include "TFile.h"
+#include "TTree.h"
+
+using namespace std;
+
+// Generate root file output from input file from the SHMS calorimeter
+// simulation.
+
+#define Nrows_pr 14
+#define Ncols_pr 2
+#define Nrows_sh 16
+#define Ncols_sh 14
+
+void make_root(string inp_name="5", int p_cent=5.) {
+
+ // Input and output.
+
+ string inp_file = inp_name + ".r";
+ string out_file = inp_name + ".root";
+ cout << inp_file << " ==> " << out_file << endl;
+
+ ifstream fin;
+ fin.open(inp_file.c_str(),ios::in);
+
+ TFile *file_out = new TFile(out_file.c_str(), "Recreate");
+
+ // Variables.
+
+ int nhit;
+ double x, y, xp, yp, p;
+ int ind;
+ double tg_dp;
+
+ double adc_pr[Nrows_pr][Ncols_pr];
+ double adc_sh[Nrows_sh][Ncols_sh];
+
+ // The output tree.
+
+ TTree *tr = new TTree("T", "SHMS calo. MC events");
+ string leaf = Form("P.pr.a_p[%i][%i]/D", Nrows_pr, Ncols_pr);
+ tr->Branch("P.pr.a_p", &adc_pr, leaf.c_str());
+ leaf = Form("P.sh.a_p[%i][%i]/D", Nrows_sh, Ncols_sh);
+ tr->Branch("P.sh.a_p", &adc_sh, leaf.c_str());
+ tr->Branch("P.tr.y", &x, "P.tr.y/D"); // swap X in MC and Y in hcana
+ tr->Branch("P.tr.x", &y, "P.tr.x/D"); // swap X in MC and Y in hcana
+ tr->Branch("P.tr.ph", &xp, "P.tr.ph/D"); // swap Xp in MC and Yp in hcana
+ tr->Branch("P.tr.th", &yp, "P.tr.th/D"); // swap Xp in MC and Yp in hcana
+ tr->Branch("P.tr.p", &p, "P.tr.p/D");
+ tr->Branch("P.tr.tg_dp", &tg_dp, "P.tr.tg_dp/D");
+
+ // Read and fill tree cicle.
+
+ while (fin >> nhit >> x >> y >> xp >> yp >> p >> ind) {
+
+ // cout << nhit << " " << x << " " << y << " " << xp << " " << yp << " "
+ // << p << " " << ind << endl;
+
+ tg_dp = (p/p_cent - 1.) * 100.;
+
+ for (int k=0; k<Ncols_pr; k++)
+ for (int i=0; i<Nrows_pr; i++)
+ adc_pr[i][k] = 0.;
+
+ for (int k=0; k<Ncols_sh; k++)
+ for (int i=0; i<Nrows_sh; i++)
+ adc_sh[i][k] = 0.;
+
+ for (int i=0; i<nhit; i++) {
+
+ int adc, row, col, det;
+ fin >> adc >> row >> col >> det;
+ // cout << adc << " " << row << " " << col << " " << det << endl;
+
+ if (adc > 0) {
+ switch (det) {
+ case 1 : //Preshower
+ adc_pr[row-1][col-1] = adc;
+ break;
+ case 2 : //Shower
+ adc_sh[row-1][col-1] = adc;
+ break;
+ default:
+ cout << "*** Wrong detector: " << det << " ***" << endl;
+ }
+ }
+
+ } //hits
+
+ tr->Fill();
+
+ } //events
+
+ tr->Write();
+
+ fin.close();
+ file_out->Close();
+}
diff --git a/hc_cal_calib/shms_cal_calib/mc_results/memo.txt b/hc_cal_calib/shms_cal_calib/mc_results/memo.txt
new file mode 100644
index 0000000..0c6b45b
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/mc_results/memo.txt
@@ -0,0 +1,4 @@
+From /scratch/tadevosn/shms_lg/g4.simon/mar2012.
+
+Removed degradation of quantum efficiencies of PMTs, both in Preshower
+and Shower.
diff --git a/hc_cal_calib/shms_cal_calib/new_old_diff/comp_pcal_etot.C b/hc_cal_calib/shms_cal_calib/new_old_diff/comp_pcal_etot.C
new file mode 100644
index 0000000..aa7ad0e
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/new_old_diff/comp_pcal_etot.C
@@ -0,0 +1,127 @@
+#include <TFile.h>
+#include <TROOT.h>
+#include <iostream>
+#include <TMath.h>
+#include <TLegend.h>
+
+using namespace std;
+
+void comp_pcal_etot() {
+
+ // Compare total energy depositions in the SHMS calorimeter,
+ // obtained from calibrations with old (Fortran) and new (C++?Root) codes
+ // of simulated data.
+
+ // Histograms.
+
+ TH1F* he_new = new TH1F("e_new", "Edep/P from C++/Root calibration",
+ 100, 0.011, 1.8);
+ he_new->SetLineColor(kRed);
+
+ TH1F* he_old = (TH1F*) he_new->Clone("Edep/P from Fortran calibration");
+ he_old->SetLineColor(kBlue);
+
+ TH1F* ee_diff = new TH1F("ee_diff",
+ "Event by event new - old Ecalo difference", 200,-1.,1.);
+ TH1F* ee_abs_diff = new TH1F("ee_abs_diff",
+ "Event by event absolute new - old Ecalo difference",
+ 100000,0.,0.1);
+ TH2F* evse = new TH2F("evse","Event by event new vs old Ecalo",
+ 210,-0.1,2., 210,-0.1,2.);
+ evse->SetMarkerStyle(3);
+
+ TH3F* diff_vs_xy = new TH3F("diff_vs_xy","Event by event difference vs X & Y",
+ 100, -50., 50., 120, -60., 60., 400, -0.02, 0.02);
+ // diff_vs_xy->SetMarkerStyle(3);
+
+ TH2F* diff_vs_x = new TH2F("diff_vs_x","Event by event difference vs X",
+ 100, -50., 50., 400, -0.02, 0.02);
+
+ TH2F* diff_vs_y = new TH2F("diff_vs_y","Event by event difference vs Y",
+ 120, -60., 60., 400, -0.02, 0.02);
+
+ // Read data and fill in histograms.
+
+ ifstream fin_old;
+ fin_old.open("old_calib.cal_data",ios::in);
+
+ ifstream fin_new;
+ fin_new.open("new_calibrated.d",ios::in);
+
+ double e_old, p_old, y_old;
+ double e_new, p_new, x_new, y_new;
+
+ // Loop over entries.
+
+ while ((fin_old >> e_old >> p_old >> y_old) &&
+ (fin_new >> e_new >> p_new >> x_new >> y_new)) {
+
+ double en_old = e_old/p_old;
+ double en_new = e_new/p_new;
+
+ he_new->Fill(en_new, 1.);
+ he_old->Fill(en_old, 1.);
+ ee_diff->Fill(en_new - en_old, 1.);
+ ee_abs_diff->Fill(TMath::Abs(en_new - en_old), 1.);
+ evse->Fill(en_old, en_new, 1.);
+ diff_vs_xy->Fill(x_new, y_new, en_new-en_old, 1.);
+ diff_vs_x->Fill(x_new, en_new-en_old, 1.);
+ diff_vs_y->Fill(y_new, en_new-en_old, 1.);
+
+ }
+
+ fin_new.close();
+ fin_old.close();
+
+ // Plot.
+
+ TCanvas* c1 = new TCanvas("c1","Normalized energy depositons",600,750);
+
+ c1->Divide(1,2);
+ c1->cd(1);
+ he_new->Draw();
+ he_old->Draw("same");
+
+ TLegend* leg = new TLegend(0.25, 0.70, 0.40, 0.82);
+ leg->AddEntry(he_new,"New( C++/Root) calibration","l");
+ leg->AddEntry(he_old,"Old (Fortran) calibration","l");
+ leg->Draw();
+
+ // Difference between the histograms.
+
+ c1->cd(2);
+
+ TH1F* dif = he_new->Clone("difference");
+ dif->Reset();
+ dif->Add(he_new,he_old,1.,-1.);
+
+ dif->SetTitle("new (C++/Root) -- old (Fortran) difference");
+ dif->SetFillColor(kRed);
+ dif->SetLineColor(kRed);
+ dif->SetLineWidth(1);
+ dif->SetFillStyle(1111);
+ dif->Draw();
+
+ TCanvas* c2 = new TCanvas("c2","Event by event new - old Ecalo diff.",
+ 600,750);
+
+ c2->Divide(1,2);
+ c2->cd(1);
+ gPad->SetLogy();
+ ee_diff->Draw();
+ c2->cd(2);
+ gPad->SetLogy();
+ gPad->SetLogx();
+ ee_abs_diff->Draw();
+
+ TCanvas* c3 = new TCanvas("c3","Event by event new - old vs X & Y",
+ 750, 750);
+ c3->Divide(2,2);
+ c3->cd(1);
+ diff_vs_x->Draw();
+ c3->cd(2);
+ diff_vs_y->Draw();
+ c3->cd(3);
+ diff_vs_xy->Draw();
+
+}
diff --git a/hc_cal_calib/shms_cal_calib/old_code/calib.kumac b/hc_cal_calib/shms_cal_calib/old_code/calib.kumac
new file mode 100644
index 0000000..1fc5cdd
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/old_code/calib.kumac
@@ -0,0 +1,87 @@
+macro calib P=1
+
+set *
+opt *
+ve/del *
+hi/del 0
+close 0
+
+opt date
+opt grid
+opt stat
+opt fit
+
+fn='../results/e/'[P]'.r'
+
+sh ./shms_calib < [fn]
+
+ve/read e,eb,x shms_calib.cal_data
+sigma r=e/eb
+
+1d 1 'Edep/P' 100 0.5 1.5
+**1d 1 'Edep/P' 200 0.5 1.5
+hfill r 1
+
+rms=$HINFO(1,'RMS')
+mean=$HINFO(1,'MEAN')
+lo=[mean]-3*[rms]
+hi=[mean]+3*[rms]
+**lo=[mean]-2*[rms]
+**hi=[mean]+2*[rms]
+mess lo=[lo] hi=[hi]
+hi/fit 1([lo]:[hi]) g
+
+*------------------------------------------------------------------------------
+ve/cr fpar(3) r
+ve/cr ferr(3) r
+Application COMIS Quit
+ SUBROUTINE chiq
+* NPFITS Number of points used in the fit
+* NFPAR Number of free parameters
+* FITCHI Chisquare
+* FITPAR Values of parameters
+* FITSIG Errors on parameters
+ COMMON/HCFITS/NCFITS,NPFITS,NFPAR,FITCHI,FITPAR(35),FITSIG(35)
+ + ,FITDER(35)
+ vector fpar
+ vector ferr
+ fpar(1)=fitpar(1)
+ fpar(2)=fitpar(2)
+ fpar(3)=fitpar(3)
+ ferr(1)=fitsig(1)
+ ferr(2)=fitsig(2)
+ ferr(3)=fitsig(3)
+ END
+Quit
+*------------------------------------------------------------------------------
+
+call chiq
+sig=fpar(3)
+ser=ferr(3)
+mean=fpar(2)
+sigm=[sig]/[mean]*100
+serr=[ser]/[mean]*100
+mess sigma = [sigm] +/- [serr] %
+
+lo=[mean]-3*[sig]
+hi=[mean]+3*[sig]
+1d 2 'Edep/P' 100 [lo] [hi]
+
+hfill r 2
+sum=$HINFO(2,'SUM')
+entries=$HINFO(2,'Entries')
+eff=[sum]/[entries]*100
+
+mess eff = [eff] %
+
+* Output.
+
+rms=[rms]/[mean]*100
+
+ve/cr v(5) r
+ve/inp v [P] [rms] [sigm] [serr] [eff]
+ve/write v tmp.d '(f3.0,3f6.3,f6.2)'
+sh cat tmp.d >> resolution.dat
+sh rm tmp.d
+
+return
diff --git a/hc_cal_calib/shms_cal_calib/old_code/s_correct.f b/hc_cal_calib/shms_cal_calib/old_code/s_correct.f
new file mode 100644
index 0000000..48149d9
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/old_code/s_correct.f
@@ -0,0 +1,77 @@
+ function s_correct_pr(xhit,icol)
+
+c Correct for lateral coordinate of imapct point.
+c Fit to GEANT pion data @ 5 GeV/c (Simon).
+
+ if((icol.ne.1).and.(icol.ne.2)) then
+ print*,'*** s_correct_pr: icol=',icol,' ***'
+ stop '*** wrong icol in s_correct_pr ***'
+ end if
+
+c Check hit coordinate with fired block number.
+c With Vardan's simulated data:
+c if((xhit.lt.0..and.icol.eq.1).or.(xhit.gt.0..and.icol.eq.2)) then
+c With Simon's simulated data:
+ if((xhit.lt.0..and.icol.eq.2).or.(xhit.gt.0..and.icol.eq.1)) then
+ x=xhit
+ else
+ x=0.
+ end if
+
+ if(x.lt.0.) x=-x
+
+cc s_correct_pr=1.0328+(x/73.395)**3.7036
+cc s_correct_pr=1.0328/s_correct_pr
+
+cc s_correct_pr=1.+(x/82.348)**2.537
+cc s_correct_pr=1./s_correct_pr
+
+cc s_correct_pr=1.+(x/91.51)**2.89
+cc s_correct_pr=1./s_correct_pr
+
+ s_correct_pr=1.+(x/106.73)**2.329
+ s_correct_pr=1./s_correct_pr
+
+cD print*,'s_correct_pr = ',s_correct_pr, ' x = ',x,
+cD * ' xhit = ',xhit, ' icol = ',icol
+
+ end
+
+*=======================================================================
+c
+c function s_correct_cal(xhit,iblk)
+c
+cc Correct for lateral coordinate of imapct point.
+cc GEANT muon data @ 1 GeV/c.
+c
+c real xh(14),ef(14)
+c data xh/ 0., 5., 10., 15., 20., 25., 30., 35., 40.,
+c * 45., 50., 55., 60., 65./
+c data ef/1.036,1.030,1.029,1.031,1.043,1.064,1.082,1.113,1.146,
+c * 1.192,1.271,1.360,1.497,1.697/
+c save xh,ef
+c
+c parameter (xstp=5.)
+c
+cc Check hit coordinate with fired block number.
+c if((xhit.lt.0.and.iblk.eq.1).or.(xhit.gt.0.).and.(iblk.eq.2)) then
+c x=xhit
+c else
+c x=0.
+c end if
+c
+c if(x.lt.0.) x=-x
+c
+c i=x/xstp+1
+c if(i.lt.1) i=1
+c if(i.gt.14) i=14
+c eff=ef(i)+(ef(i+1)-ef(i))/(xh(i+1)-xh(i))*(x-xh(i))
+c
+c s_correct_cal=ef(1)/eff
+c end
+
+*=======================================================================
+
+ function s_correct_sh(xh,yh,icol,irow)
+ s_correct_sh=1.
+ end
diff --git a/hc_cal_calib/shms_cal_calib/old_code/shms_calib.f b/hc_cal_calib/shms_cal_calib/old_code/shms_calib.f
new file mode 100644
index 0000000..4944a37
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/old_code/shms_calib.f
@@ -0,0 +1,509 @@
+*=======================================================================
+ program shms_calib
+*=======================================================================
+
+c
+c SHMS calorimeter calibration with electrons.
+c
+
+*
+ implicit none
+*
+
+c
+ real thr_lo,thr_hi
+ parameter (thr_lo=0.,thr_hi=1.E+8) !calibrate all the events.
+
+ print*,'========================================================='
+ print*,'Calibrating SHMS Calorimeter'
+
+ call shms_clb_det(thr_lo,thr_hi)
+
+ print*,'========================================================='
+
+ end
+
+*=======================================================================
+
+ subroutine shms_clb_det(thr_lo,thr_hi)
+ implicit none
+c
+ real thr_lo,thr_hi
+
+ integer nrow_pr,ncol_sh,nrow_sh
+cc parameter ( nrow_pr=11,ncol_sh=14,nrow_sh=11)
+ parameter ( nrow_pr=14,ncol_sh=14,nrow_sh=16)
+ integer npmts
+ parameter (npmts=2*nrow_pr+ncol_sh*nrow_sh) !shms2
+ integer npmts2
+ parameter (npmts2=npmts*npmts)
+ real*8 q0(npmts)
+ real*8 qm(npmts,npmts)
+ real*8 qe(npmts)
+ real*8 q(npmts)
+ real*8 eb
+ real*8 e0
+ real*8 ac(npmts)
+ real*8 au(npmts)
+ real*8 t
+ real*8 s
+ integer nev
+ logical*1 eod
+ integer i,j
+ integer nf(npmts)
+ integer minf
+ parameter (minf=5) ! minimum number to hit pmt to include pmt in calib.
+ integer nums(npmts)
+ integer numsel
+ real*8 q0s(npmts)
+ real*8 qes(npmts)
+ integer nsi,nsj
+ real*8 acs(npmts)
+ real*8 aus(npmts)
+ real*8 aux(npmts2)
+ integer jp
+ integer spare_id
+ parameter (spare_id=77)
+ character*40 fn
+
+ real xh
+
+ do i=1,npmts
+ q0(i)=0.
+ qe(i)=0.
+ do j=1,npmts
+ qm(i,j)=0.
+ end do
+ au(i)=0.
+ ac(i)=0.
+ nf(i)=0
+ end do
+ e0=0.
+c
+ nev=0
+ eod=.false.
+ do while(.not.eod)
+ call s_get_data(eb,q,xh,eod,thr_lo,thr_hi)
+ if(.not.eod) then
+ do i=1,npmts
+ if(q(i).gt.0.) then
+ q0(i)=q0(i)+q(i)
+ qe(i)=qe(i)+eb*q(i)
+ if(i.eq.6) then
+ccD print*,'q0(6)=',q0(6)
+ccD print*,'qe(6)=',qe(6)
+ccD print*,'q(6)=',q(6)
+ end if
+ do j=1,npmts
+ qm(i,j)=qm(i,j)+q(i)*q(j)
+ccD if(i.eq.6.and.j.eq.6) print*,'qm(6,6)=',qm(6,6)
+ end do
+ nf(i)=nf(i)+1
+ end if
+ end do
+ e0=e0+eb
+ nev=nev+1
+ if(nev/1000*1000.eq.nev) write(*,'(e10.3,i7)') e0,nev
+ end if
+ end do
+
+ccD print*,'q0(6)=',q0(6)
+ccD print*,'qe(6)=',qe(6)
+ccD print*,'qm(6,6)=',qm(6,6)
+ccD pause
+
+ do i=1,npmts
+ q0(i)=q0(i)/nev
+ qe(i)=qe(i)/nev
+ do j=1,npmts
+ qm(i,j)=qm(i,j)/nev
+ end do
+ end do
+ e0=e0/nev
+
+c Debug outputs.
+
+ open(spare_id, file='qm.d')
+ do j=1,npmts
+ do i=1,npmts
+ write(spare_id,*) qm(i,j),i,j
+ end do
+ end do
+ close(spare_id)
+
+ open(spare_id, file='q0.d')
+ do i=1,npmts
+ write(spare_id,*) q0(i),i
+ end do
+ close(spare_id)
+
+ open(spare_id, file='qe.d')
+ do i=1,npmts
+ write(spare_id,*) qe(i),i
+ end do
+ close(spare_id)
+
+ numsel=0
+ do i=1,npmts
+ if(nf(i).ge.minf) then
+ numsel=numsel+1
+ nums(numsel)=i
+ccD print*,nums(numsel),numsel,nf(i)
+ else
+ write(*,*) ' PMT ',i,' only ',nf(i),
+ * ' events. Will not to be calibrated. Gain is set to 0.'
+ end if
+ end do
+ccD print*,'numsel =',numsel
+
+ write(*,'(''Total '',i6,'' events processed'')') nev
+ write(*,*) ' PMT with less than', minf,
+ * ' events are not included in calibration.'
+ write(*,*)
+ write(*,11) 'shms_neg_pre_hits=',(nf(i),i= 1, nrow_pr)
+ write(*,11) 'shms_pos_pre_hits=',(nf(i),i=nrow_pr+1,2*nrow_pr)
+ write(*,11) 'shms_shower_hits =',
+ * (nf(i),i=2*nrow_pr+1,2*nrow_pr+nrow_sh)
+ do j=2,ncol_sh
+ write(*,11) ' ',
+ * (nf(i),i=2*nrow_pr+(j-1)*nrow_sh+1,
+ * 2*nrow_pr+j*nrow_sh)
+ end do
+c
+ do i=1,numsel
+ nsi=nums(i)
+ q0s(i)=q0(nsi)
+ qes(i)=qe(nsi)
+ do j=1,numsel
+ nsj=nums(j)
+ jp=j+(i-1)*numsel
+ aux(jp)=qm(nsj,nsi)
+ccD write(65,'(e12.5)') aux(jp)
+ end do
+ end do
+
+ call calib(e0,q0s,qes,aux,numsel,numsel*numsel,aus,acs)
+
+ do i=1,numsel
+ nsi=nums(i)
+ au(nsi)=aus(i)
+ ac(nsi)=acs(i)
+ end do
+
+ccD write(*,'(2e10.3,i5)') (ac(i),au(i),i,i=1,npmts)
+
+ if(e0.lt.10.) then
+ write(fn,'(i1,a6)') int(e0),'.param'
+ else
+ write(fn,'(i2,a6)') int(e0),'.param'
+ end if
+
+ open(spare_id,file=fn)
+
+cc write(spare_id,'(''; Calibration constants for SHMS momentum'',f5.1,
+cc * '' GeV/c, '',i6,'' events processed'')') e0,nev
+ write(spare_id,20) e0,nev
+ write(spare_id,*)
+
+ write(spare_id,10)
+ & 'shms_neg_pre_gain=',(ac(i)*1.D+3,i= 1, nrow_pr)
+ write(spare_id,10)
+ & 'shms_pos_pre_gain=',(ac(i)*1.D+3,i=nrow_pr+1,2*nrow_pr)
+ write(spare_id,10) 'shms_shower_gain =',
+ & (ac(i)*1.D+3,i=2*nrow_pr+1,2*nrow_pr+nrow_sh)
+ do j=2,ncol_sh
+ write(spare_id,10) ' ',
+ & (ac(i)*1.D+3,i=2*nrow_pr+(j-1)*nrow_sh+1,
+ * 2*nrow_pr+j*nrow_sh)
+ end do
+ close(spare_id)
+
+ write(*,20) e0,nev
+ write(*,*)
+ write(*,10)
+ & 'shms_neg_pre_gain=',(ac(i)*1.D+3,i= 1, nrow_pr)
+ write(*,10)
+ & 'shms_pos_pre_gain=',(ac(i)*1.D+3,i=nrow_pr+1,2*nrow_pr)
+ write(*,10) 'shms_shower_gain =',
+ & (ac(i)*1.D+3,i=2*nrow_pr+1,2*nrow_pr+nrow_sh)
+ do j=2,ncol_sh
+ write(*,10) ' ',
+ & (ac(i)*1.D+3,i=2*nrow_pr+(j-1)*nrow_sh+1,
+ * 2*nrow_pr+j*nrow_sh)
+ end do
+
+ 20 format('; Calibration constants for SHMS momentum',f5.1,
+ * ' GeV/c, ',i6,' events processed')
+ 10 format(a18,30(f6.3,','))
+ 11 format(a18,30(i5,','))
+
+ open(spare_id,file='shms_calib.cal_data')
+ write(*,*) 'Creating file s_cal_calib.cal_data, channel ',
+ * spare_id
+
+ rewind(5)
+ eod=.false.
+ nev=0
+ do while(.not.eod)
+ call s_get_data(eb,q,xh,eod,0.,1.E+8)
+ if(.not.eod) then
+ s=0.
+* t=0.
+ do i=1,npmts
+ s=s+q(i)*ac(i)
+* t=t+q(i)*au(i)
+ end do
+ write(spare_id,*) s,eb,xh
+ end if
+ end do
+
+ close(5)
+ close(spare_id)
+
+ end
+
+*=======================================================================
+
+ subroutine calib(e0,q0,qe,aux,npmts,npmts2,au,ac)
+ implicit none
+ integer npmts,npmts2
+ real*8 e0
+ real*8 q0(npmts)
+ real*8 qe(npmts)
+ real*8 aux(npmts2)
+ real*8 ac(npmts)
+ real*8 au(npmts)
+ real*8 qm(npmts,npmts)
+ real*8 t
+ real*8 s
+ integer ifail
+ integer i,j
+ integer jp
+c
+ do i=1,npmts
+ do j=1,npmts
+ jp=j+(i-1)*npmts
+ qm(j,i)=aux(jp)
+cD write(66,'(e12.5,2i4)') qm(j,i),j,i
+ end do
+ end do
+c
+ print*,'Calib: npmts =',npmts
+ print*,' '
+c
+ print*,'Inversing the Matrix...'
+ call smxinv(qm,npmts,ifail)
+ if(ifail.ne.0) then
+ stop '*** Singular Matrix ***'
+ else
+ print*,' ...done.'
+ end if
+c
+ do i=1,npmts
+ au(i)=0.
+ do j=1,npmts
+ au(i)=au(i)+qm(i,j)*qe(j)
+ end do
+ end do
+c
+ s=0.
+ do i=1,npmts
+ t=0.
+ do j=1,npmts
+ t=t+qm(i,j)*q0(j)
+ end do
+ s=s+q0(i)*t
+ end do
+c
+ t=0.
+ do i=1,npmts
+ t=t+au(i)*q0(i)
+ end do
+ s=(e0-t)/s
+c
+ do i=1,npmts
+ t=0.
+ do j=1,npmts
+ t=t+qm(i,j)*q0(j)
+ end do
+ ac(i)=s*t+au(i)
+ end do
+c
+ end
+
+*-----------------------------------------------------------------------
+
+ subroutine s_get_data(eb,q,xh,eod,thr_lo,thr_hi)
+ implicit none
+c
+ real*8 eb
+ integer nrow_pr,ncol_sh,nrow_sh
+cc parameter (nrow_pr=11,ncol_sh=14,nrow_sh=11)
+ parameter (nrow_pr=14,ncol_sh=14,nrow_sh=16)
+ integer*4 num_pmts
+ parameter (num_pmts=2*nrow_pr+ncol_sh*nrow_sh) !shms2
+ real*8 q(num_pmts)
+ logical*1 eod
+
+ integer*4 nhit
+ real*4 qdc
+ integer*4 nh
+ integer*4 nb
+ integer irow,icol,idet
+c
+ real*4 xh,yh, xph,yph
+ real*4 s_correct_pr,s_correct_sh
+ real*4 thr_lo,thr_hi
+ logical*1 good_ev
+ real*4 qnet
+
+ real thr
+ parameter (thr=15.)
+
+ good_ev=.false.
+ do while(.not.good_ev)
+
+ eb=0.d0
+ do nb=1,num_pmts
+ q(nb)=0.d0
+ end do
+ qnet=0.
+ eod=.true.
+
+ read(*,*,end=5) nhit,xh,yh,xph,yph,eb !x,y,xp,yp @ FP
+cc read(*,*,end=5) nhit,yh,xh,eb
+ do nh=1,nhit
+
+ read(*,*,end=5) qdc,irow,icol,idet
+ if(qdc.gt.thr) then
+
+ if(idet.eq.1) then !preshower
+ nb=(icol-1)*nrow_pr+irow
+ q(nb)=qdc*s_correct_pr(xh+275.*xph,icol) !275cm - PrSwr dist.
+ccD if(icol.eq.1.and.irow.eq.6) then
+ccD print*,'get_data: qdc, q, nb=',qdc,q(nb),nb
+ccD print*,' xh, cor=',xh,s_correct_pr(xh,icol)
+ccD end if
+ if(nb.eq.18) print*,'nb, xh, icol, cor = ',
+ * nb, xh, icol, s_correct_pr(xh,icol)
+
+ elseif(idet.eq.2) then !shower
+ nb=2*nrow_pr+(icol-1)*nrow_sh+irow
+ q(nb)=qdc*s_correct_sh(xh,yh,icol,irow)
+ else
+ stop '*** wrong detector identifier ***'
+ end if
+
+ qnet=qnet+q(nb)
+ccD print*,qdc,idet,irow,icol,nb
+
+ end if !qdc>thr
+
+ enddo !1,nhit
+ eod=.false.
+
+ qnet=qnet/(eb*1000.)
+ good_ev=(qnet.gt.thr_lo).and.(qnet.lt.thr_hi)
+
+cD write(99,*) qnet
+
+ end do !.not.good_ev
+
+ 5 continue
+
+ end
+
+*-----------------------------------------------------------------------
+
+ SUBROUTINE SMXINV (A,NDIM,IFAIL)
+C
+C CERN PROGLIB# F107 SMXINV .VERSION KERNFOR 1.0 720503
+C ORIG. 03/05/72 CL
+C
+*** REAL*8 A(*),RI(100)
+*** INTEGER*4 INDEX(100)
+cc REAL*8 A(*),RI(200)
+cc INTEGER*4 INDEX(200)
+ REAL*8 A(*),RI(300)
+ INTEGER*4 INDEX(300)
+C
+ DATA TOL / 1.D-14/
+C
+ IFAIL=0
+ N=NDIM
+ NP1=N+1
+ DO 10 I=1,N
+ 10 INDEX(I)=1
+C
+ DO 80 I=1,N
+C
+C-- FIND PIVOT
+ PIVOT=0.0D0
+ JJ=1
+ DO 20 J=1,N
+ IF (INDEX(J).EQ.0) GO TO 19
+ ELM=DABS (A(JJ))
+ IF (ELM.LE.PIVOT) GO TO 19
+ PIVOT=ELM
+ K=J
+ KK=JJ
+ 19 JJ=JJ+NP1
+ 20 CONTINUE
+ IF (PIVOT/DABS(A(1)).LT.TOL) GO TO 100
+ INDEX(K)=0
+ PIVOT=-A(KK)
+C
+C-- ELIMINATION
+ KJ=K
+ NP=N
+C
+ DO 70 J=1,N
+ IF (J-K) 34,30,34
+C
+ 30 A(KJ)=1.0D0/PIVOT
+ RI(J)=0.0D0
+ NP=1
+ GO TO 70
+C
+ 34 ELM=-A(KJ)
+ 40 RI(J)=ELM/PIVOT
+ IF (ELM.EQ.0.0D0) GO TO 50
+C
+ JL=J
+ DO 45 L=1,J
+ A(JL)=A(JL)+ELM*RI(L)
+ 45 JL=JL+N
+C
+ 50 A(KJ)=RI(J)
+C
+ 70 KJ=KJ+NP
+C
+ 80 CONTINUE
+C
+C-- CHANGE THE SIGN AND PROVISIONAL FILL-UP
+ IJ0=1
+ JI0=1
+ DO 95 I=1,N
+ IJ=IJ0
+ JI=JI0
+C
+ DO 90 J=1,I
+ A(IJ)=-A(IJ)
+ A(JI)=A(IJ)
+ IJ=IJ+N
+ JI=JI+1
+ 90 CONTINUE
+C
+ IJ0=IJ0+1
+ JI0=JI0+N
+ 95 CONTINUE
+ RETURN
+C
+C-- FAILURE RETURN
+ 100 IFAIL=1
+ RETURN
+ END
+
+*=======================================================================
+ include 's_correct.f'
+
diff --git a/hc_cal_calib/shms_cal_calib/old_params/11.param b/hc_cal_calib/shms_cal_calib/old_params/11.param
new file mode 100644
index 0000000..2bfbad5
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/old_params/11.param
@@ -0,0 +1,18 @@
+; Calibration constants for SHMS momentum 11.6 GeV/c, 40000 events processed
+
+shms_neg_pre_gain= 0.806, 1.085, 0.738, 0.989, 0.633, 1.412, 0.749, 0.779, 0.655, 1.033, 0.604, 0.799, 5.938, 0.000,
+shms_pos_pre_gain= 1.517, 0.936, 1.248, 0.644, 1.047, 1.013, 0.689, 1.348, 1.356, 0.691, 1.587, 0.555, 2.527, 0.000,
+shms_shower_gain = 0.000, 0.000, 0.000,-2.352, 0.000, 0.000, 0.000,-0.101, 0.000, 0.000, 1.047, 0.000, 0.000, 0.000, 0.000, 0.000,
+ 0.565,-0.823, 2.383,-1.089,25.026, 0.000, 1.715, 3.599, 3.860, 1.577,-0.397, 0.370, 3.988,-6.333, 0.000, 0.000,
+ 0.261, 0.555, 0.780, 1.034,-0.634, 1.067, 8.527, 0.696, 0.847, 0.703, 0.639, 0.895, 0.617, 0.779,-1.255, 0.000,
+ 0.660, 0.553, 1.057, 0.861, 0.861, 0.228, 1.213, 0.830, 1.287, 1.173, 0.554, 0.538, 0.556, 1.052, 0.139, 0.000,
+ 0.862, 0.640, 0.901, 0.698, 0.490, 0.504, 0.624, 0.941, 0.813, 1.241, 0.687, 1.382, 1.052, 0.737, 1.545, 4.147,
+ 1.125, 1.001, 0.618, 0.566, 0.523, 0.546, 0.584, 0.708, 0.543, 0.531, 0.789, 0.721, 0.920, 0.689, 1.756, 3.531,
+ 0.796, 0.570, 0.536, 1.228, 0.479, 0.462, 0.668, 1.015, 0.916, 0.529, 1.088, 0.557, 1.107, 0.515, 0.137,-0.291,
+ 0.739, 0.649, 0.485, 0.963, 0.474, 1.245, 0.536, 0.520, 0.487, 0.568, 0.645, 0.810, 0.684, 1.302, 0.251,-0.176,
+ 0.556, 0.498, 0.796, 0.868, 0.926, 0.536, 1.224, 0.774, 0.677, 1.304, 1.345, 0.992, 0.636, 0.910,-1.690, 0.000,
+ 0.700, 1.322, 0.523, 0.607, 0.709, 0.971, 0.563, 0.673, 0.871, 0.600, 1.123, 0.828, 0.500, 1.326, 0.589, 0.000,
+ 0.476, 0.666, 0.624, 0.663, 0.619, 1.907, 1.371, 1.130, 0.740, 0.631, 0.909, 0.768, 0.873, 0.803, 0.323, 0.000,
+ 3.263, 0.467, 0.907, 1.066,-0.290, 0.701,-1.482, 0.606, 0.880, 0.431, 0.575, 0.988, 0.906, 0.540, 0.723,-1.256,
+ 1.098, 1.428, 0.000,-1.427,-1.641, 8.426, 8.996, 0.000, 4.226, 8.884, 2.032,-1.224, 0.819,-4.338, 0.000, 0.000,
+ 0.000, 0.732, 0.000, 0.000, 0.000, 0.000, 0.000, 4.147, 8.043,26.418, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
diff --git a/hc_cal_calib/shms_cal_calib/old_params/5.param b/hc_cal_calib/shms_cal_calib/old_params/5.param
new file mode 100644
index 0000000..7ec5453
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/old_params/5.param
@@ -0,0 +1,18 @@
+; Calibration constants for SHMS momentum 5.3 GeV/c, 40000 events processed
+
+shms_neg_pre_gain= 0.717, 1.025, 0.703, 0.941, 0.595, 1.329, 0.704, 0.728, 0.621, 0.973, 0.569, 0.745, 3.579, 6.302,
+shms_pos_pre_gain= 1.405, 0.889, 1.184, 0.612, 0.987, 0.958, 0.663, 1.272, 1.274, 0.652, 1.503, 0.513, 0.863, 0.000,
+shms_shower_gain = 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, 3.664,-0.585, 1.641,-1.201, 0.000, 6.724, 0.904, 0.000, 0.000, 4.384, 1.765, 0.000,-1.670, 0.000, 0.000,
+ 1.697, 0.555, 0.768, 2.119,-0.278, 1.736, 0.000,-0.820, 2.365, 0.555, 0.646, 0.819, 0.542, 1.208, 0.932, 0.000,
+ 0.424, 0.554, 1.075, 0.880, 0.859, 0.817, 0.387, 0.339, 1.224, 1.197, 0.564, 0.549, 0.574, 0.964, 1.815, 0.000,
+ 0.897, 0.656, 0.927, 0.712, 0.498, 0.509, 0.646, 0.949, 0.836, 1.260, 0.713, 1.407, 1.072, 0.777, 1.772, 0.000,
+ 1.198, 1.026, 0.633, 0.580, 0.538, 0.558, 0.598, 0.725, 0.559, 0.551, 0.817, 0.743, 0.951, 0.718, 0.373, 0.000,
+ 0.822, 0.581, 0.550, 1.258, 0.493, 0.478, 0.686, 1.036, 0.942, 0.550, 1.119, 0.578, 1.144, 0.544, 0.114, 0.000,
+ 0.834, 0.664, 0.495, 0.984, 0.486, 1.280, 0.553, 0.535, 0.505, 0.583, 0.665, 0.833, 0.710, 1.335, 0.894, 1.903,
+ 0.661, 0.509, 0.816, 0.884, 0.946, 0.556, 1.256, 0.801, 0.701, 1.340, 1.389, 1.015, 0.656, 0.952, 0.067, 0.000,
+ 0.799, 1.332, 0.530, 0.618, 0.726, 1.009, 0.578, 0.691, 0.904, 0.615, 1.154, 0.847, 0.518, 1.318, 3.123, 0.000,
+ 0.629, 0.682, 0.643, 0.674, 0.650, 1.405, 0.865, 0.837, 0.731, 0.645, 0.923, 0.786, 0.900, 0.818, 0.000, 0.000,
+ 0.035, 0.496, 0.836,-0.061, 1.129, 0.564, 0.000, 0.118, 0.505, 0.851, 0.616, 0.984, 0.946, 0.402,-0.605, 0.000,
+ 0.000, 4.184, 0.000,-0.322, 0.000, 0.000, 3.362, 0.000, 0.000,-2.433,-0.359, 1.658,-1.390, 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,
diff --git a/hc_cal_calib/shms_cal_calib/pcal_calib.cpp b/hc_cal_calib/shms_cal_calib/pcal_calib.cpp
new file mode 100644
index 0000000..292bc0b
--- /dev/null
+++ b/hc_cal_calib/shms_cal_calib/pcal_calib.cpp
@@ -0,0 +1,47 @@
+#include "TCanvas.h"
+#include "THcPShowerCalib.h"
+
+//
+// A steering Root script for the SHMS calorimeter calibration.
+//
+
+void pcal_calib(Int_t RunNumber) {
+
+ cout << "Calibrating run " << RunNumber << endl;
+
+ THcPShowerCalib theShowerCalib(RunNumber);
+
+ theShowerCalib.Init(); // Initialize constants and variables
+ theShowerCalib.CalcThresholds(); // Thresholds on the uncalibrated Edep/P
+ theShowerCalib.ComposeVMs(); // Compute vectors amd matrices for calib.
+ theShowerCalib.SolveAlphas(); // Solve for the calibration constants
+ theShowerCalib.SaveAlphas(); // Save the constants
+ theShowerCalib.SaveRawData(); // Save raw data into file for debug purposes
+ theShowerCalib.FillHEcal(); // Fill histograms
+
+ // Plot histograms
+
+ TCanvas* Canvas =
+ new TCanvas("Canvas", "PHMS Shower Counter calibration", 1000, 667);
+ Canvas->Divide(2,2);
+
+ Canvas->cd(1);
+
+ // Normalized uncalibrated energy deposition.
+
+ theShowerCalib.hEunc->DrawCopy();
+
+ theShowerCalib.hEuncSel->SetFillColor(kGreen);
+ theShowerCalib.hEuncSel->DrawCopy("same");
+
+ // Normalized energy deposition after calibration.
+
+ Canvas->cd(3);
+ theShowerCalib.hEcal->Fit("gaus");
+
+ // SHMS delta(P) versus the calibrated energy deposition.
+
+ Canvas->cd(4);
+ theShowerCalib.hDPvsEcal->Draw();
+
+}
diff --git a/hcal_calib/make_cratemap.pl b/hcal_calib/make_cratemap.pl
deleted file mode 120000
index 8d3d2ef..0000000
--- a/hcal_calib/make_cratemap.pl
+++ /dev/null
@@ -1 +0,0 @@
-../examples/make_cratemap.pl
\ No newline at end of file
--
GitLab