Skip to content
Snippets Groups Projects
Commit 6df1766e authored by Vardan Tadevosyan's avatar Vardan Tadevosyan Committed by Stephen A. Wood
Browse files

Clean up and comment HMS calorimeter calibration code.

parent d721dced
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
hcal_calib/THcShHit.h
+ 8
− 4
View file @ 6df1766e
......@@ -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;
};
......
hcal_calib/THcShTrack.h
+ 18
− 28
View file @ 6df1766e
......@@ -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;
......
hcal_calib/THcShowerCalib.h
+ 30
− 46
View file @ 6df1766e
......@@ -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;
......
hcal_calib/hcal_calib.cpp
+ 3
− 5
View file @ 6df1766e
......@@ -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();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment