From 95db9d0937c7481eaf07eea71398754f42913ab0 Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Fri, 21 Feb 2014 21:04:29 +0400
Subject: [PATCH] Add the calibration code.
---
hcal_calib/THcShHit.h | 64 +++
hcal_calib/THcShTrack.h | 180 ++++++++
hcal_calib/THcShowerCalib.h | 696 ++++++++++++++++++++++++++++++
hcal_calib/hcal_calib.cpp | 36 ++
hcal_calib/hcal_replay.cpp | 74 ++++
hcal_calib/hcal_replay_cuts.def | 40 ++
hcal_calib/output_hcal_replay.def | 18 +
7 files changed, 1108 insertions(+)
create mode 100644 hcal_calib/THcShHit.h
create mode 100644 hcal_calib/THcShTrack.h
create mode 100644 hcal_calib/THcShowerCalib.h
create mode 100644 hcal_calib/hcal_calib.cpp
create mode 100644 hcal_calib/hcal_replay.cpp
create mode 100644 hcal_calib/hcal_replay_cuts.def
create mode 100644 hcal_calib/output_hcal_replay.def
diff --git a/hcal_calib/THcShHit.h b/hcal_calib/THcShHit.h
new file mode 100644
index 0000000..7285b63
--- /dev/null
+++ b/hcal_calib/THcShHit.h
@@ -0,0 +1,64 @@
+#include <iostream>
+
+class THcShHit {
+
+ Double_t ADCpos, ADCneg;
+ Double_t Epos, Eneg;
+ 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();
+};
+
+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() {
+ cout << "Hit: ADCpos =" << ADCpos << " ADCneg =" << ADCneg
+ << " Epos =" << Epos << " Eneg =" << Eneg
+ << " BlkNumber=" << BlkNumber << endl;
+};
+
+struct pmt_hit {Double_t signal; UInt_t channel;};
diff --git a/hcal_calib/THcShTrack.h b/hcal_calib/THcShTrack.h
new file mode 100644
index 0000000..5597364
--- /dev/null
+++ b/hcal_calib/THcShTrack.h
@@ -0,0 +1,180 @@
+#include "THcShHit.h"
+#include "TMath.h"
+
+#include <vector>
+#include <iterator>
+#include <iostream>
+
+using namespace std;
+
+// Container (collection) of hits and its iterator.
+//
+typedef vector<THcShHit*> THcShHitList;
+typedef THcShHitList::iterator THcShHitIt;
+
+class THcShTrack {
+
+ UInt_t Nhits;
+ Double_t P;
+ Double_t X;
+ Double_t Xp;
+ Double_t Y;
+ Double_t Yp;
+
+ THcShHitList Hits;
+
+ public:
+ THcShTrack();
+ THcShTrack(UInt_t nh, 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 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 Nhits;};
+ void Print();
+ Bool_t CheckHitNumber();
+ void SetEs(Double_t* alpha);
+ Double_t Enorm();
+ Double_t GetP() {return P*1000.;} //MeV
+
+ Float_t Ycor(Double_t);
+ Float_t Ycor(Double_t, Int_t);
+
+ static const Double_t fAcor = 200.;
+ static const Double_t fBcor = 8000.;
+ static const Double_t fCcor = 64.36;
+ static const Double_t fDcor = 1.66;
+
+ static const Double_t fZbl = 10;
+ static const UInt_t fNrows = 13;
+ static const UInt_t fNcols = 4;
+ static const UInt_t fNnegs = 26;
+ static const UInt_t fNpmts = 78;
+ static const UInt_t fNblks = fNrows*fNcols;
+
+};
+
+THcShTrack::THcShTrack() { };
+
+THcShTrack::THcShTrack(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;
+};
+
+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::AddHit(Double_t adc_pos, Double_t adc_neg,
+ Double_t e_pos, Double_t e_neg,
+ UInt_t blk_number) {
+ 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() {
+ cout << "ShTrack: P=" << P << " X=" << X << " Xp=" << Xp
+ << " Y=" << Y << " Yp=" << Yp << " Nhits=" << Nhits << endl;
+ cout << "Hits size=" << Hits.size() << endl;
+
+ for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+ (*iter)->Print();
+ };
+
+};
+
+Bool_t THcShTrack::CheckHitNumber() {
+ return (Nhits == Hits.size());
+};
+
+THcShTrack::~THcShTrack() {
+ for (THcShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
+ delete *i;
+ *i = 0;
+ }
+};
+
+//------------------------------------------------------------------------------
+
+void THcShTrack::SetEs(Double_t* 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() {
+
+ Double_t sum = 0;
+
+ for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+ sum += (*iter)->GetEpos();
+ sum += (*iter)->GetEneg();
+ };
+
+ return sum/P/1000.;
+}
+
+//------------------------------------------------------------------------------
+
+//Coordinate correction for single PMT modules.
+//PMT attached at right (positive) side.
+
+Float_t THcShTrack::Ycor(Double_t y) {
+ return TMath::Exp(y/fAcor)/(1. + y*y/fBcor);
+}
+
+//Coordinate correction for double PMT modules.
+//
+
+Float_t THcShTrack::Ycor(Double_t y, Int_t side) {
+ if (side!=0&&side!=1) {
+ cout << "THcShower::Ycor : wrong side " << side << endl;
+ return 0.;
+ }
+ Int_t sign = 1 - 2*side;
+ return (fCcor + sign*y)/(fCcor + sign*y/fDcor);
+}
diff --git a/hcal_calib/THcShowerCalib.h b/hcal_calib/THcShowerCalib.h
new file mode 100644
index 0000000..6f9d151
--- /dev/null
+++ b/hcal_calib/THcShowerCalib.h
@@ -0,0 +1,696 @@
+#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"
+
+using namespace std;
+
+// HMS Shower Counter calibration class
+
+class THcShowerCalib {
+
+ public:
+ THcShowerCalib(Int_t);
+ THcShowerCalib();
+ ~THcShowerCalib();
+
+ void ExtractData();
+ void Init();
+ Bool_t ReadShRawTrack(THcShTrack &trk);
+ void CalcThresholds();
+ void ComposeVMs();
+ void SolveAlphas();
+ void FillHEcal();
+ void SaveAlphas();
+ // void PlotHistos();
+
+ TH1F* hEunc;
+ TH1F* hEuncSel;
+ TH1F* hEcal;
+ TH2F* hPvsEcal;
+
+ private:
+ Int_t fRunNumber;
+ Double_t fLoThr;
+ Double_t fHiThr;
+ UInt_t fNev;
+ static const UInt_t fMinHitCount = 200;
+
+ ifstream fDataStream;
+
+ // TVectorD* fe0;
+ // TVectorD* fq0;
+ // TMatrixDSym* fQ;
+ // TVectorD* falphaU;
+ // TVectorD* falphaC;
+
+ 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];
+ Double_t falphaC[THcShTrack::fNpmts];
+ Double_t falpha0[THcShTrack::fNpmts];
+ Double_t falpha1[THcShTrack::fNpmts];
+
+ UInt_t fHitCount[THcShTrack::fNpmts];
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShowerCalib::THcShowerCalib() {};
+
+//------------------------------------------------------------------------------
+
+THcShowerCalib::THcShowerCalib(Int_t RunNumber) {
+ fRunNumber = RunNumber;
+};
+
+//------------------------------------------------------------------------------
+
+THcShowerCalib::~THcShowerCalib() {
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::ExtractData() {
+
+ char* fname = Form("Root_files/hcal_calib_%d.root",fRunNumber);
+ cout << "THcShowerCalib::ExtractData: fname= " << fname << endl;
+
+ //Reset ROOT and connect tree file
+ gROOT->Reset();
+ TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
+ if (!f) {
+ f = new TFile(fname);
+ }
+ TTree* tree;
+ f->GetObject("T",tree);
+
+ //Declaration of leaves types
+
+ // Int_t Ndata_H_cal_1pr_aneg_p;
+ 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];
+
+ Double_t H_cal_trp;
+ Double_t H_cal_trx;
+ Double_t H_cal_trxp;
+ Double_t H_cal_try;
+ Double_t H_cal_tryp;
+
+ // Set branch addresses.
+
+ // tree->SetBranchAddress("Ndata.H.cal.1pr.aneg_p",&Ndata_H_cal_1pr_aneg_p);
+ tree->SetBranchAddress("H.cal.1pr.aneg_p",H_cal_1pr_aneg_p);
+ tree->SetBranchAddress("H.cal.1pr.apos_p",H_cal_1pr_apos_p);
+
+ tree->SetBranchAddress("H.cal.2ta.aneg_p",H_cal_2ta_aneg_p);
+ tree->SetBranchAddress("H.cal.2ta.apos_p",H_cal_2ta_apos_p);
+
+ tree->SetBranchAddress("H.cal.3ta.aneg_p",H_cal_3ta_aneg_p);
+ tree->SetBranchAddress("H.cal.3ta.apos_p",H_cal_3ta_apos_p);
+
+ tree->SetBranchAddress("H.cal.4ta.aneg_p",H_cal_4ta_aneg_p);
+ tree->SetBranchAddress("H.cal.4ta.apos_p",H_cal_4ta_apos_p);
+
+ tree->SetBranchAddress("H.cal.trp",&H_cal_trp);
+ tree->SetBranchAddress("H.cal.trx",&H_cal_trx);
+ tree->SetBranchAddress("H.cal.trxp",&H_cal_trxp);
+ tree->SetBranchAddress("H.cal.try",&H_cal_try);
+ tree->SetBranchAddress("H.cal.tryp",&H_cal_tryp);
+
+ Long64_t nentries = tree->GetEntries();
+ cout << "THcShowerCalib::ExtractData: nentries= " << nentries << endl;
+
+ char* FName = Form("raw_data/%d_raw.dat",fRunNumber);
+ cout << "ExtractData: FName=" << FName << endl;
+ ofstream output;
+ output.open(FName,ios::out);
+
+ Long64_t nbytes = 0;
+ for (Long64_t i=0; i<nentries;i++) {
+ nbytes += tree->GetEntry(i);
+ Int_t nhit = 0;
+ for (Int_t j=0; j<THcShTrack::fNrows; j++) {
+ if (H_cal_1pr_apos_p[j]>0. || H_cal_1pr_aneg_p[j]>0.) nhit++;
+ if (H_cal_2ta_apos_p[j]>0. || H_cal_2ta_aneg_p[j]>0.) nhit++;
+ if (H_cal_3ta_apos_p[j]>0. || H_cal_3ta_aneg_p[j]>0.) nhit++;
+ if (H_cal_4ta_apos_p[j]>0. || H_cal_4ta_aneg_p[j]>0.) nhit++;
+ }
+ output << nhit << " " << H_cal_trp << " "
+ << H_cal_trx << " " << H_cal_trxp << " "
+ << H_cal_try << " " << H_cal_tryp << endl;
+ for (Int_t j=0; j<THcShTrack::fNrows; j++) {
+ if (H_cal_1pr_apos_p[j]>0. || H_cal_1pr_aneg_p[j]>0.)
+ output << H_cal_1pr_apos_p[j] << " " << H_cal_1pr_aneg_p[j] << " "
+ << j+1 << endl;
+ }
+ for (Int_t j=0; j<THcShTrack::fNrows; j++) {
+ if (H_cal_2ta_apos_p[j]>0. || H_cal_2ta_aneg_p[j]>0.)
+ output << H_cal_2ta_apos_p[j] << " " << H_cal_2ta_aneg_p[j] << " "
+ << j+1 + THcShTrack::fNrows << endl;
+ }
+ for (Int_t j=0; j<THcShTrack::fNrows; j++) {
+ if (H_cal_3ta_apos_p[j]>0. || H_cal_3ta_aneg_p[j]>0.)
+ output << H_cal_3ta_apos_p[j] << " " << H_cal_3ta_aneg_p[j] << " "
+ << j+1 + 2*THcShTrack::fNrows << endl;
+ }
+ for (Int_t j=0; j<THcShTrack::fNrows; j++) {
+ if (H_cal_4ta_apos_p[j]>0. || H_cal_4ta_aneg_p[j]>0.)
+ output << H_cal_4ta_apos_p[j] << " " << H_cal_4ta_aneg_p[j] << " "
+ << j+1 + 3*THcShTrack::fNrows << endl;
+ }
+
+ }
+
+ output.close();
+ cout << "THcShowerCalib::ExtractData: nbytes= " << nbytes << endl;
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::Init() {
+
+ char* FName = Form("raw_data/%d_raw.dat",fRunNumber);
+ cout << "Init: FName=" << FName << endl;
+ fDataStream.open(FName,ios::in);
+
+ 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,1.9,2.4);
+hPvsEcal = new TH2F("hPvsEcal", "P versus Edep/P ",150,0.,1.5, 100,0.7,0.9);
+
+ for (UInt_t i=0; i<THcShTrack::fNpmts; i++) fHitCount[i] = 0;
+
+ // fe0 = new TVectorD(THcShRawTrack::fNpmts);
+ // fq0 = new TVectorD(THcShRawTrack::fNpmts);
+ // fQ = new TMatrixDSym(THcShRawTrack::fNpmts,THcShRawTrack::fNpmts);
+ // falphaU = new TVectorD(THcShRawTrack::fNpmts);
+ // falphaC = new TVectorD(THcShRawTrack::fNpmts);
+
+ 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.;
+ }
+ }
+
+ 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.;
+ }
+ };
+
+ for (UInt_t ipmt=0; ipmt<THcShTrack::fNpmts; ipmt++) {
+ falpha1[ipmt] = 1.;
+ }
+
+ // for (UInt_t ipmt=0; ipmt<THcShRawTrack::fNpmts; ipmt++) {
+ // cout << "falpha0 " << ipmt << " = " << falpha0[ipmt] << endl;
+ // }
+
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::CalcThresholds() {
+
+ Int_t nev = 0;
+
+ THcShTrack trk;
+
+ while (ReadShRawTrack(trk)) {
+
+ // trk.Print();
+
+ trk.SetEs(falpha0);
+ 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;
+
+ 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;
+
+ 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_t THcShowerCalib::ReadShRawTrack(THcShTrack &trk) {
+
+ UInt_t nhit;
+ Double_t p, x,xp, y,yp;
+
+ if (fDataStream >> nhit >> p >> x >> xp >> y >> yp) {
+ // cout << nhit << " " << p << " " << x << " " << xp << " " << y << " "
+ // << yp << endl;
+
+ trk.SetTrack(nhit, p, x, xp, y, yp);
+
+ for (UInt_t i=0; i<nhit; i++) {
+
+ Double_t adc_pos, adc_neg;
+ UInt_t nb;
+
+ fDataStream >> adc_pos >> adc_neg >> nb;
+ // cout << adc_pos << " " << adc_neg << " " << nb << endl;
+
+ trk.AddHit(adc_pos, adc_neg, 0., 0., nb);
+ };
+ return 1;
+
+ }
+ else {
+ return 0;
+ };
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::ComposeVMs() {
+
+ // Reset.
+ fDataStream.clear() ;
+ fDataStream.seekg(0, ios::beg) ;
+
+ fNev = 0;
+ THcShTrack trk;
+
+ while (ReadShRawTrack(trk)) {
+
+ if (trk.CheckHitNumber()) {
+
+ trk.SetEs(falpha0);
+ Double_t Enorm = trk.Enorm();
+ if (Enorm>fLoThr && Enorm<fHiThr) {
+
+ trk.SetEs(falpha1);
+ // trk.Print();
+
+ fe0 += trk.GetP();
+
+ vector<pmt_hit> pmt_hit_list;
+
+ for (UInt_t i=0; i<trk.GetNhits(); i++) {
+
+ THcShHit* hit = trk.GetHit(i);
+ // hit->Print();
+
+ UInt_t nb = hit->GetBlkNumber();
+
+ fqe[nb-1] += hit->GetEpos() * trk.GetP();
+ fq0[nb-1] += hit->GetEpos();
+
+ pmt_hit_list.push_back( pmt_hit{hit->GetEpos(), nb} );
+
+ fHitCount[nb-1]++;
+
+ 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]++;
+ };
+
+ // cout << " fqe " << nb-1 << " =" << fqe[nb-1];
+ // if (nb <= THcShTrack::fNnegs)
+ // cout << " fqe " << THcShTrack::fNblks+nb-1 << " ="
+ // << fqe[THcShTrack::fNblks+nb-1];
+ // cout << endl;
+
+ // cout << " fq0 " << nb-1 << " =" << fq0[nb-1];
+ // if (nb <= THcShTrack::fNnegs)
+ // cout << " fq0 " << THcShTrack::fNblks+nb-1 << " ="
+ // << fq0[THcShTrack::fNblks+nb-1];
+ // cout << endl;
+ } //over hits
+
+
+ // for (vector<pmt_hit>::iterator it=pmt_hit_list.begin();
+ // it < pmt_hit_list.end(); it++) {
+ // cout << "hit: " << (*it).signal << " " << (*it).channel
+ // << endl;
+ // }
+
+ 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;
+ }
+ }
+
+ // getchar();
+
+ fNev++;
+ };
+ };
+ };
+
+ 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;
+
+ // cout << "ComposeVMs: fNev=" << fNev << endl;
+ // cout << " fe0=" << fe0 << endl;
+ // cout << " fqe:" << endl;
+ // for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+ // cout << " fqe " << i << " = " << fqe[i] << endl;
+ // cout << " fq0:" << endl;
+ // for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+ // cout << " fq0 " << i << " = " << fq0[i] << endl;
+
+ // cout << " fQ[i,i], fQ[i,1], fQ[1,i]:" << endl;
+ // for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+ // cout << i << ": " << fQ[i][i] << " " << fQ[i][1] << " " << fQ[1][i]
+ // << endl;
+
+ // cout << "ComposeVMs: vector fe0:" << endl;
+ // (*fe0).Print();
+ // cout << "ComposeVMs: matrix fQ:" << endl;
+ // (*fQ).Print();
+
+ ofstream q0out;
+ q0out.open("q0.d",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);
+ for (UInt_t i=0; i<THcShTrack::fNpmts; i++)
+ qeout << fqe[i] << " " << i << endl;
+ qeout.close();
+
+
+ ofstream Qout;
+ Qout.open("Q.d",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() {
+
+ 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;
+
+ cout << "Hit couts:" << 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;
+
+ 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++) {
+
+ 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;
+ }
+ }
+
+ 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.
+
+ cout << endl;
+ cout << "Channels with hit number less than " << fMinHitCount
+ << " will not be calibrated." << endl;
+ cout << endl;
+
+ for (UInt_t i=0; i<THcShTrack::fNpmts; i++) {
+
+ if (fHitCount[i] < fMinHitCount) {
+ cout << "Channel " << i << ", " << fHitCount[i]
+ << " hits, will not be calibrated." << endl;
+ q0[i] = 0.;
+ qe[i] = 0.;
+ for (UInt_t k=0; k<THcShTrack::fNpmts; k++) {
+ Q[i][k] = 0.;
+ Q[k][i] = 0.;
+ }
+ Q[i][i] = 1.;
+ }
+
+ }
+
+ 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;
+
+ au = lu.Solve(qe,ok);
+ cout << "au: ok=" << ok << endl;
+ // au.Print();
+
+ Double_t t1 = fe0 - au * q0;
+ // cout << "t1 =" << t1 << endl;
+
+ TVectorD Qiq0(THcShTrack::fNpmts);
+ Qiq0 = lu.Solve(q0,ok);
+ cout << "Qiq0: ok=" << ok << endl;
+ // Qiq0.Print();
+
+ Double_t t2 = q0 * Qiq0;
+ // cout << "t2 =" << t2 << endl;
+
+ ac = (t1/t2) *Qiq0 + au;
+ // cout << "ac:" << endl;
+ // ac.Print();
+
+ for (UInt_t i=0; i<THcShTrack::fNpmts; i++) {
+ falphaU[i] = au[i];
+ falphaC[i] = ac[i];
+ }
+
+}
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::FillHEcal() {
+
+ // Reset.
+ fDataStream.clear() ;
+ fDataStream.seekg(0, ios::beg) ;
+
+ ofstream output;
+ output.open("calibrated.d",ios::out);
+
+ Int_t nev = 0;
+
+ THcShTrack trk;
+
+ while (ReadShRawTrack(trk)) {
+
+ // trk.Print();
+
+ trk.SetEs(falphaC);
+ // trk.SetEs(falphaU);
+ Double_t P = trk.GetP();
+ Double_t Enorm = trk.Enorm();
+
+ hEcal->Fill(Enorm);
+ hPvsEcal->Fill(Enorm,P/1000.,1.);
+
+ output << Enorm*P/1000. << " " << P/1000. << endl;
+
+ nev++;
+ };
+
+ output.close();
+
+ cout << "FillHEcal: " << nev << " events filled" << endl;
+};
+
+//------------------------------------------------------------------------------
+
+void THcShowerCalib::SaveAlphas() {
+
+ ofstream output;
+ char* fname = Form("hcal.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;
+ 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();
+}
+
+//------------------------------------------------------------------------------
+
+/*
+void THcShowerCalib::PlotHistos() {
+
+ TCanvas* Canvas =
+ new TCanvas("Canvas", "HMS Shower Counter calibration", 1000, 667);
+ Canvas->Divide(1,2);
+
+ Canvas->cd(1);
+
+ hEunc->DrawCopy();
+
+ hEuncSel->SetFillColor(kGreen);
+ hEuncSel->DrawCopy("same");
+
+};
+*/
+
+#endif
diff --git a/hcal_calib/hcal_calib.cpp b/hcal_calib/hcal_calib.cpp
new file mode 100644
index 0000000..23c5437
--- /dev/null
+++ b/hcal_calib/hcal_calib.cpp
@@ -0,0 +1,36 @@
+#include "TCanvas.h"
+#include "THcShowerCalib.h"
+
+void hcal_calib(Int_t RunNumber) {
+
+ cout << "Calibrating run " << RunNumber << endl;
+
+ THcShowerCalib theShowerCalib(RunNumber);
+
+ theShowerCalib.ExtractData();
+ theShowerCalib.Init();
+ theShowerCalib.CalcThresholds();
+ theShowerCalib.ComposeVMs();
+ theShowerCalib.SolveAlphas();
+ theShowerCalib.FillHEcal();
+ theShowerCalib.SaveAlphas();
+
+ TCanvas* Canvas =
+ new TCanvas("Canvas", "HMS Shower Counter calibration", 1000, 667);
+ Canvas->Divide(2,2);
+
+ Canvas->cd(1);
+
+ theShowerCalib.hEunc->DrawCopy();
+
+ theShowerCalib.hEuncSel->SetFillColor(kGreen);
+ theShowerCalib.hEuncSel->DrawCopy("same");
+
+ Canvas->cd(3);
+ // theShowerCalib.hEcal->Draw();
+ theShowerCalib.hEcal->Fit("gaus");
+
+ Canvas->cd(4);
+ theShowerCalib.hPvsEcal->Draw();
+
+}
diff --git a/hcal_calib/hcal_replay.cpp b/hcal_calib/hcal_replay.cpp
new file mode 100644
index 0000000..bc6680f
--- /dev/null
+++ b/hcal_calib/hcal_replay.cpp
@@ -0,0 +1,74 @@
+void hcal_replay(Int_t RunNumber=52949, Int_t MaxEventToReplay=11000) {
+
+ char* RunFileNamePattern="daq04_%d.log.0";
+
+ gHcParms->Define("gen_run_number", "Run Number", RunNumber);
+ gHcParms->AddString("g_ctp_database_filename", "DBASE/test.database");
+
+ gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
+
+ // g_ctp_parm_filename and g_decode_map_filename should now be defined
+
+ gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
+
+ // Constants not in ENGINE PARAM files that we want to be
+ // configurable
+ gHcParms->Load("PARAM/hcana.param");
+
+ // Generate db_cratemap to correspond to map file contents
+ char command[100];
+ sprintf(command,"./make_cratemap.pl < %s > db_cratemap.dat",
+ gHcParms->GetString("g_decode_map_filename"));
+ system(command);
+
+ // Load the Hall C style detector map
+ gHcDetectorMap=new THcDetectorMap();
+ gHcDetectorMap->Load(gHcParms->GetString("g_decode_map_filename"));
+
+ // Set up the equipment to be analyzed.
+
+ THaApparatus* HMS = new THcHallCSpectrometer("H","HMS");
+ gHaApps->Add( HMS );
+
+ // Add HMS detectors
+ HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
+ HMS->AddDetector( new THcShower("cal", "Shower" ));
+ HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
+ HMS->AddDetector( new THcAerogel("aero", "Aerogel Cerenkov" ));
+ HMS->AddDetector( new THcCherenkov("cer", "Gas Cerenkov" ));
+
+ // Set up the analyzer - we use the standard one,
+ // but this could be an experiment-specific one as well.
+ // The Analyzer controls the reading of the data, executes
+ // tests/cuts, loops over Acpparatus's and PhysicsModules,
+ // and executes the output routines.
+ THaAnalyzer* analyzer = new THcAnalyzer;
+
+ // A simple event class to be output to the resulting tree.
+ // Creating your own descendant of THaEvent is one way of
+ // defining and controlling the output.
+ THaEvent* event = new THaEvent;
+
+ // Define the run(s) that we want to analyze.
+ // We just set up one, but this could be many.
+ char RunFileName[100];
+ sprintf(RunFileName,RunFileNamePattern,RunNumber);
+ THaRun* run = new THaRun(RunFileName);
+
+ // Eventually need to learn to skip over, or properly analyze
+ // the pedestal events
+ run->SetEventRange(1,MaxEventToReplay); // Physics Event number, does not
+ // include scaler or control events
+
+ // Define the analysis parameters
+ analyzer->SetCountMode( 0 ); //Mark's modification
+ analyzer->SetEvent( event );
+ analyzer->SetOutFile(Form("Root_files/hcal_calib_%05d.root",RunNumber));
+ analyzer->SetOdefFile("output_hcal_replay.def");
+ analyzer->SetCutFile("hcal_replay_cuts.def"); // optional
+
+ // File to record cuts accounting information
+ // analyzer->SetSummaryFile("summary_example.log"); // optional
+
+ analyzer->Process(run); // start the actual analysis
+}
diff --git a/hcal_calib/hcal_replay_cuts.def b/hcal_calib/hcal_replay_cuts.def
new file mode 100644
index 0000000..7bc402d
--- /dev/null
+++ b/hcal_calib/hcal_replay_cuts.def
@@ -0,0 +1,40 @@
+# Cuts for the HMS calorimeter calibration.
+#
+
+Block: RawDecode
+
+Pedestal_event g.evtyp==4
+scalar_event g.evtyp==0
+HMS_event g.evtyp==1
+SOS_event g.evtyp==2
+coin_event g.evtyp==3
+misc_event g.evtyp>=5
+hms_and_coin HMS_event||coin_event
+RawDecode_master 1
+
+Block: Decode
+Decode_master hms_and_coin
+
+Block: CoarseTracking
+CoarseTracking_master !Pedestal_event
+
+Block: CoarseReconstruct
+RawCoarseReconstruct !Pedestal_event
+
+Block: Reconstruct
+
+one_track H.dc.ntrack==1
+#one_clust H.cal.nclust==1
+one_sh_track H.cal.ntracks==1
+in_delta H.cal.trdelta>-10.&&H.cal.trdelta<10.
+good_cer H.cer.npesum>3.
+good_beta H.cal.trbeta>0.9&&H.cal.trbeta<1.1
+in_calx H.cal.trx>-60.&&H.cal.trx<60.
+in_caly H.cal.try>-30.&&H.cal.try<30.
+in_cal in_calx&&in_caly
+
+Reconstruct_master one_track&&in_delta
+#Reconstruct_master one_track&&one_sh_track&&in_delta&&good_cer
+#Reconstruct_master one_track&&one_clust&&in_delta&&good_cer
+#Reconstruct_master one_track&&one_clust&&in_delta&&good_cer&&in_cal
+#&&good_beta
diff --git a/hcal_calib/output_hcal_replay.def b/hcal_calib/output_hcal_replay.def
new file mode 100644
index 0000000..d1b7ce0
--- /dev/null
+++ b/hcal_calib/output_hcal_replay.def
@@ -0,0 +1,18 @@
+#
+# Output definition for the HMS calorimeter calibration.
+#
+
+#block H.dc.*
+#block H.hod.*
+#block H.cal.*
+#block H.aero.*
+#block H.cer.*
+#block g.evtyp
+
+block H.cal.*
+variable H.cer.npesum
+variable H.dc.ntrack
+
+#variable H.dc.delta
+#variable H.dc.P
+#variable H.hod.beta
--
GitLab