From 6df1766e953d2abb03b3a56426e89cfd7128acc1 Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Sat, 15 Mar 2014 10:50:41 +0400
Subject: [PATCH] Clean up and comment HMS calorimeter calibration code.
---
hcal_calib/THcShHit.h | 12 ++++--
hcal_calib/THcShTrack.h | 46 +++++++++-------------
hcal_calib/THcShowerCalib.h | 76 +++++++++++++++----------------------
hcal_calib/hcal_calib.cpp | 8 ++--
4 files changed, 59 insertions(+), 83 deletions(-)
diff --git a/hcal_calib/THcShHit.h b/hcal_calib/THcShHit.h
index 5807edf..20afcce 100644
--- a/hcal_calib/THcShHit.h
+++ b/hcal_calib/THcShHit.h
@@ -5,7 +5,7 @@
class THcShHit {
Double_t ADCpos, ADCneg; // pedestal subtracted ADC signals.
- Double_t Epos, Eneg;
+ Double_t Epos, Eneg; // Energy depositions seen from pos. & neg. sides
UInt_t BlkNumber;
public:
@@ -38,6 +38,8 @@ class THcShHit {
void Print(ostream & ostrm);
};
+//------------------------------------------------------------------------------
+
THcShHit::THcShHit() {
ADCpos = -99999.;
ADCneg = -99999.;
@@ -57,10 +59,12 @@ THcShHit::THcShHit(Double_t adc_pos, Double_t adc_neg,
THcShHit::~THcShHit() { };
+//------------------------------------------------------------------------------
+
void THcShHit::Print(ostream & ostrm) {
- // ostrm << "Hit: ADCpos =" << ADCpos << " ADCneg =" << ADCneg
- // << " Epos =" << Epos << " Eneg =" << Eneg
- // << " BlkNumber=" << BlkNumber << endl;
+
+ // Output hit data through the stream ostrm.
+
ostrm << ADCpos << " " << ADCneg << " " << Epos << " " << Eneg << " "
<< BlkNumber << endl;
};
diff --git a/hcal_calib/THcShTrack.h b/hcal_calib/THcShTrack.h
index defd5be..d8c74d9 100644
--- a/hcal_calib/THcShTrack.h
+++ b/hcal_calib/THcShTrack.h
@@ -19,7 +19,6 @@ typedef THcShHitList::iterator THcShHitIt;
class THcShTrack {
- // UInt_t Nhits;
Double_t P; // track momentum
Double_t X; // at the calorimater face
Double_t Xp; // slope
@@ -30,13 +29,9 @@ class THcShTrack {
public:
THcShTrack();
- // THcShTrack(UInt_t nh, Double_t p,
- // Double_t x, Double_t xp, Double_t y, Double_t yp);
THcShTrack(Double_t p, Double_t x, Double_t xp, Double_t y, Double_t yp);
~THcShTrack();
- // void SetTrack(UInt_t nh, Double_t p,
- // Double_t x, Double_t xp, Double_t y, Double_t yp);
void Reset(Double_t p, Double_t x, Double_t xp, Double_t y, Double_t yp);
void AddHit(Double_t adc_pos, Double_t adc_neg,
@@ -49,8 +44,6 @@ class THcShTrack {
void Print(ostream & ostrm);
- // Bool_t CheckHitNumber();
-
void SetEs(Double_t* alpha);
Double_t Enorm();
@@ -78,12 +71,12 @@ class THcShTrack {
};
+//------------------------------------------------------------------------------
+
THcShTrack::THcShTrack() { };
-//THcShTrack::THcShTrack(UInt_t nh, Double_t p,
THcShTrack::THcShTrack(Double_t p,
Double_t x, Double_t xp, Double_t y, Double_t yp) {
- // Nhits = nh;
P = p;
X = x;
Xp = xp;
@@ -91,19 +84,13 @@ THcShTrack::THcShTrack(Double_t p,
Yp =yp;
};
-//void THcShTrack::SetTrack(UInt_t nh, Double_t p,
-// Double_t x, Double_t xp, Double_t y, Double_t yp) {
-// Nhits = nh;
-// P = p;
-// X = x;
-// Xp = xp;
-// Y = y;
-// Yp =yp;
-// // Hits.clear();
-//};
+//------------------------------------------------------------------------------
void THcShTrack::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;
@@ -112,15 +99,22 @@ void THcShTrack::Reset(Double_t p,
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++;
@@ -128,9 +122,9 @@ THcShHit* THcShTrack::GetHit(UInt_t k) {
}
void THcShTrack::Print(ostream & ostrm) {
- // ostrm << "ShTrack: P=" << P << " X=" << X << " Xp=" << Xp
- // << " Y=" << Y << " Yp=" << Yp << endl;
- // ostrm << "Hits size=" << Hits.size() << endl;
+
+ // Output the track parameters and hit list through the stream ostrm.
+
ostrm << P << " " << X << " " << Xp << " " << Y << " " << Yp << " "
<< Hits.size() << endl;
@@ -140,11 +134,7 @@ void THcShTrack::Print(ostream & ostrm) {
};
-// Check hit number with the size of hit collection.
-//
-//Bool_t THcShTrack::CheckHitNumber() {
-// return (Nhits == Hits.size());
-//};
+//------------------------------------------------------------------------------
THcShTrack::~THcShTrack() {
for (THcShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
@@ -186,7 +176,7 @@ void THcShTrack::SetEs(Double_t* alpha) {
Double_t THcShTrack::Enorm() {
- // Normalized to track momentum energy depostion in the calorimeter.
+ // Normalized to the track momentum energy depostion in the calorimeter.
Double_t sum = 0;
diff --git a/hcal_calib/THcShowerCalib.h b/hcal_calib/THcShowerCalib.h
index 9c5e03a..0d92f5d 100644
--- a/hcal_calib/THcShowerCalib.h
+++ b/hcal_calib/THcShowerCalib.h
@@ -18,7 +18,9 @@
using namespace std;
-// HMS Shower Counter calibration class
+//
+// HMS Shower Counter calibration class.
+//
class THcShowerCalib {
@@ -39,7 +41,6 @@ class THcShowerCalib {
TH1F* hEunc;
TH1F* hEuncSel;
TH1F* hEcal;
- // TH2F* hPvsEcal;
TH2F* hDPvsEcal;
private:
@@ -87,7 +88,8 @@ THcShowerCalib::~THcShowerCalib() {
void THcShowerCalib::SaveRawData() {
- // Output raw data into file for debug purposes.
+ // 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;
@@ -110,13 +112,8 @@ void THcShowerCalib::SaveRawData() {
void THcShowerCalib::Init() {
- // Open the raw data file.
+ //Reset ROOT and connect tree file.
- // char* FName = Form("raw_data/%d_raw.dat",fRunNumber);
- // cout << "Init: FName=" << FName << endl;
- // // fDataStream.open(FName,ios::in);
-
- //Reset ROOT and connect tree file
gROOT->Reset();
char* fname = Form("Root_files/hcal_calib_%d.root",fRunNumber);
@@ -132,7 +129,6 @@ void THcShowerCalib::Init() {
hEunc = new TH1F("hEunc", "Edep/P uncalibrated", 500, 0., 5.);
hEcal = new TH1F("hEcal", "Edep/P calibrated", 150, 0., 1.5);
- // hPvsEcal = new TH2F("hPvsEcal", "P versus Edep/P ",150,0.,1.5, 100,0.7,0.9);
hDPvsEcal = new TH2F("hDPvsEcal", "#DeltaP versus Edep/P ",
150,0.,1.5, 250,-12.5,12.5);
@@ -152,7 +148,7 @@ void THcShowerCalib::Init() {
}
}
- // Initial gains (0.5 for the 2 first columns, 1 for others),
+ // 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) {
@@ -170,24 +166,22 @@ void THcShowerCalib::Init() {
falpha1[ipmt] = 1.;
}
- // for (UInt_t ipmt=0; ipmt<THcShRawTrack::fNpmts; ipmt++) {
- // cout << "falpha0 " << ipmt << " = " << falpha0[ipmt] << endl;
- // }
-
};
//------------------------------------------------------------------------------
void THcShowerCalib::CalcThresholds() {
- // Calculate +/-3 sigma thresholds on the uncalibrated total energy
+ // 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;
- // while (ReadShRawTrack(trk)) {
for (UInt_t ientry=0; ientry<fNentries; ientry++) {
ReadShRawTrack(trk, ientry);
@@ -195,7 +189,7 @@ void THcShowerCalib::CalcThresholds() {
// trk.Print(cout);
// getchar();
- trk.SetEs(falpha0);
+ trk.SetEs(falpha0); //Use initial gain constants here.
Double_t Enorm = trk.Enorm();
nev++;
@@ -234,9 +228,11 @@ void THcShowerCalib::CalcThresholds() {
void THcShowerCalib::ReadShRawTrack(THcShTrack &trk, UInt_t ientry) {
- // Set a Shower track event from the raw data.
+ //
+ // Set a Shower track event from ntuple ientry.
+ //
- //Declaration of leaves types
+ // Declaration of leaves types
// Calorimeter ADC signals.
@@ -284,8 +280,6 @@ void THcShowerCalib::ReadShRawTrack(THcShTrack &trk, UInt_t ientry) {
trk.Reset(H_cal_trp, H_cal_trx, H_cal_trxp, H_cal_try, H_cal_tryp);
- // UInt_t nhit = 0;
-
for (UInt_t j=0; j<THcShTrack::fNrows; j++) {
for (UInt_t k=0; k<THcShTrack::fNcols; k++) {
@@ -317,7 +311,6 @@ void THcShowerCalib::ReadShRawTrack(THcShTrack &trk, UInt_t ientry) {
if (adc_pos>0. || adc_neg>0.) {
trk.AddHit(adc_pos, adc_neg, 0., 0., nb);
- // nhit++;
}
}
@@ -333,17 +326,11 @@ void THcShowerCalib::ComposeVMs() {
// Fill in vectors and matrixes for the gain constant calculations.
//
- // Reset the raw data stream.
-
- // fDataStream.clear() ;
- // fDataStream.seekg(0, ios::beg) ;
-
fNev = 0;
THcShTrack trk;
- // Loop over the shower track events in the raw data stream.
+ // Loop over the shower track events in the ntuples.
- // while (ReadShRawTrack(trk)) {
for (UInt_t ientry=0; ientry<fNentries; ientry++) {
ReadShRawTrack(trk, ientry);
@@ -383,7 +370,7 @@ void THcShowerCalib::ComposeVMs() {
fHitCount[nb-1]++; //Accrue the hit counter.
- // Do the same for the negative side PMTs.
+ // Do same for the negative side PMTs.
if (nb <= THcShTrack::fNnegs) {
fqe[THcShTrack::fNblks+nb-1] += hit->GetEneg() * trk.GetP();
@@ -418,7 +405,7 @@ void THcShowerCalib::ComposeVMs() {
fNev++;
- }; // if within thresholds
+ }; // if within the thresholds
}; // over entries
@@ -462,7 +449,7 @@ void THcShowerCalib::ComposeVMs() {
void THcShowerCalib::SolveAlphas() {
//
- // Solve the sought calibration constants, by use of the Root
+ // Solve for the sought calibration constants, by use of the Root
// matrix algebra package.
//
@@ -476,6 +463,8 @@ void THcShowerCalib::SolveAlphas() {
cout << "Solving Alphas..." << endl;
cout << endl;
+ // Print out hit numbers.
+
cout << "Hit counts:" << endl;
UInt_t j = 0;
cout << "Positives:";
@@ -497,6 +486,8 @@ void THcShowerCalib::SolveAlphas() {
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];
@@ -509,7 +500,7 @@ void THcShowerCalib::SolveAlphas() {
for (UInt_t i=0; i<THcShTrack::fNpmts; i++) {
- // Check zero hit channels: the vector and matrix elements should equal 0.
+ // Check zero hit channels: the vector and matrix elements should be 0.
if (fHitCount[i] == 0) {
@@ -539,7 +530,7 @@ void THcShowerCalib::SolveAlphas() {
} //sanity check
// Low hit number channels: exclude from calculation. Assign all the
- // correspondent elements 0, except swelf-correalation Q(i,i)=1.
+ // correspondent elements 0, except self-correlation Q(i,i)=1.
cout << endl;
cout << "Channels with hit number less than " << fMinHitCount
@@ -583,12 +574,12 @@ void THcShowerCalib::SolveAlphas() {
Double_t t1 = fe0 - au * q0; // temporary variable.
// cout << "t1 =" << t1 << endl;
- TVectorD Qiq0(THcShTrack::fNpmts); // intermittend result
+ 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 tempoary variable
+ Double_t t2 = q0 * Qiq0; // another temporary variable
// cout << "t2 =" << t2 << endl;
ac = (t1/t2) *Qiq0 + au; // the sought gain constants
@@ -610,13 +601,9 @@ void THcShowerCalib::FillHEcal() {
//
// Fill histogram of the normalized energy deposition, 2-d histogram
- // of momentum versus normalized energy deposition.
+ // of momentum deviation versus normalized energy deposition.
//
- // Reset.
- // fDataStream.clear() ;
- // fDataStream.seekg(0, ios::beg) ;
-
ofstream output;
output.open("calibrated.d",ios::out);
@@ -624,14 +611,12 @@ void THcShowerCalib::FillHEcal() {
THcShTrack trk;
- // while (ReadShRawTrack(trk)) {
for (UInt_t ientry=0; ientry<fNentries; ientry++) {
ReadShRawTrack(trk, ientry);
// trk.Print(cout);
- trk.SetEs(falphaC); // use 'constrained' calibration constants
- // trk.SetEs(falphaU);
+ trk.SetEs(falphaC); // use the 'constrained' calibration constants
Double_t P = trk.GetP();
Double_t Enorm = trk.Enorm();
@@ -639,7 +624,6 @@ void THcShowerCalib::FillHEcal() {
Double_t delta;
fTree->SetBranchAddress("H.cal.trdelta",&delta);
- // hPvsEcal->Fill(Enorm,P/1000.,1.);
hDPvsEcal->Fill(Enorm,delta,1.);
output << Enorm*P/1000. << " " << P/1000. << endl;
diff --git a/hcal_calib/hcal_calib.cpp b/hcal_calib/hcal_calib.cpp
index 032364c..3ef151e 100644
--- a/hcal_calib/hcal_calib.cpp
+++ b/hcal_calib/hcal_calib.cpp
@@ -11,12 +11,12 @@ void hcal_calib(Int_t RunNumber) {
THcShowerCalib theShowerCalib(RunNumber);
- theShowerCalib.Init(); // Initialize constants adn variables
+ 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.SaveRawData(); // Save raw data into file for debug purposes
theShowerCalib.FillHEcal(); // Fill histograms
// Plot histograms
@@ -37,13 +37,11 @@ void hcal_calib(Int_t RunNumber) {
// Normalized energy deposition after calibration.
Canvas->cd(3);
- // theShowerCalib.hEcal->Draw();
theShowerCalib.hEcal->Fit("gaus");
- // Momentum versus the calibrated energy deposition.
+ // HMS delta(P) versus the calibrated energy deposition.
Canvas->cd(4);
- // theShowerCalib.hPvsEcal->Draw();
theShowerCalib.hDPvsEcal->Draw();
}
--
GitLab