diff --git a/hcal_calib/THcShTrack.h b/hcal_calib/THcShTrack.h
index 6d92a49250e315fea405d80c784adba477a84467..506bc99d1fdbcc68e102d226b18ad875cb44a7de 100644
--- a/hcal_calib/THcShTrack.h
+++ b/hcal_calib/THcShTrack.h
@@ -18,7 +18,7 @@ typedef THcShHitList::iterator THcShHitIt;
class THcShTrack {
- UInt_t Nhits;
+ // UInt_t Nhits;
Double_t P; // track momentum
Double_t X; // at the calorimater face
Double_t Xp; // slope
@@ -29,12 +29,14 @@ class THcShTrack {
public:
THcShTrack();
- THcShTrack(UInt_t nh, Double_t p,
- Double_t x, Double_t xp, Double_t y, Double_t yp);
+ // 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 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,
Double_t e_pos, Double_t e_neg,
@@ -42,11 +44,11 @@ class THcShTrack {
THcShHit* GetHit(UInt_t k);
- UInt_t GetNhits() {return Nhits;};
+ UInt_t GetNhits() {return Hits.size();};
void Print();
- Bool_t CheckHitNumber();
+ // Bool_t CheckHitNumber();
void SetEs(Double_t* alpha);
@@ -77,9 +79,10 @@ class THcShTrack {
THcShTrack::THcShTrack() { };
-THcShTrack::THcShTrack(UInt_t nh, Double_t p,
+//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;
+ // Nhits = nh;
P = p;
X = x;
Xp = xp;
@@ -87,9 +90,19 @@ THcShTrack::THcShTrack(UInt_t nh, 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;
+//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) {
P = p;
X = x;
Xp = xp;
@@ -115,7 +128,8 @@ THcShHit* THcShTrack::GetHit(UInt_t k) {
void THcShTrack::Print() {
cout << "ShTrack: P=" << P << " X=" << X << " Xp=" << Xp
- << " Y=" << Y << " Yp=" << Yp << " Nhits=" << Nhits << endl;
+ << " Y=" << Y << " Yp=" << Yp << endl;
+ // << " Y=" << Y << " Yp=" << Yp << " Nhits=" << Nhits << endl;
cout << "Hits size=" << Hits.size() << endl;
for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
@@ -126,9 +140,9 @@ void THcShTrack::Print() {
// Check hit number with the size of hit collection.
//
-Bool_t THcShTrack::CheckHitNumber() {
- return (Nhits == Hits.size());
-};
+//Bool_t THcShTrack::CheckHitNumber() {
+// return (Nhits == Hits.size());
+//};
THcShTrack::~THcShTrack() {
for (THcShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
diff --git a/hcal_calib/THcShowerCalib.h b/hcal_calib/THcShowerCalib.h
index 2c44509e512214b5241eb2cd1c20f73fcf997206..29008a947ace1ab79c501dc6277f5d9b3d9a3449 100644
--- a/hcal_calib/THcShowerCalib.h
+++ b/hcal_calib/THcShowerCalib.h
@@ -29,7 +29,7 @@ class THcShowerCalib {
void ExtractData();
void Init();
- Bool_t ReadShRawTrack(THcShTrack &trk);
+ void ReadShRawTrack(THcShTrack &trk, UInt_t ientry);
void CalcThresholds();
void ComposeVMs();
void SolveAlphas();
@@ -49,7 +49,8 @@ class THcShowerCalib {
static const UInt_t fMinHitCount = 200; // Minimum number of hits for a PMT
// to be calibrated.
- ifstream fDataStream; // Output data stream
+ TTree* fTree;
+ UInt_t fNentries;
// Quantities for calculations of the calibration constants.
@@ -94,10 +95,11 @@ void THcShowerCalib::ExtractData() {
//Reset ROOT and connect tree file
gROOT->Reset();
- TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
- if (!f) {
- f = new TFile(fname);
- }
+ // TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
+ TFile *f = new TFile(fname);
+ // if (!f) {
+ // f = new TFile(fname);
+ // }
TTree* tree;
f->GetObject("T",tree);
@@ -193,6 +195,7 @@ void THcShowerCalib::ExtractData() {
} // over entries
+ f->Close();
output.close();
cout << "THcShowerCalib::ExtractData: nbytes= " << nbytes << endl;
}
@@ -203,9 +206,21 @@ void THcShowerCalib::Init() {
// Open the raw data file.
- char* FName = Form("raw_data/%d_raw.dat",fRunNumber);
- cout << "Init: FName=" << FName << endl;
- fDataStream.open(FName,ios::in);
+ // 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);
+ cout << "THcShowerCalib::Init: Root file name = " << fname << endl;
+
+ TFile *f = new TFile(fname);
+ f->GetObject("T",fTree);
+
+ fNentries = fTree->GetEntries();
+ cout << "THcShowerCalib::Init: fNentries= " << fNentries << endl;
// Histogram declarations.
@@ -262,12 +277,15 @@ void THcShowerCalib::CalcThresholds() {
// hadronic events due to the gas Cherenkov inefficiency.
Int_t nev = 0;
-
THcShTrack trk;
- while (ReadShRawTrack(trk)) {
+ // while (ReadShRawTrack(trk)) {
+ for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+
+ ReadShRawTrack(trk, ientry);
// trk.Print();
+ // getchar();
trk.SetEs(falpha0);
Double_t Enorm = trk.Enorm();
@@ -306,35 +324,96 @@ void THcShowerCalib::CalcThresholds() {
//------------------------------------------------------------------------------
-Bool_t THcShowerCalib::ReadShRawTrack(THcShTrack &trk) {
+void THcShowerCalib::ReadShRawTrack(THcShTrack &trk, UInt_t ientry) {
+
+ // Set a Shower track event from the raw data.
+
+ //Declaration of leaves types
+
+ // Calorimeter ADC signals.
+
+ Double_t H_cal_1pr_aneg_p[THcShTrack::fNrows];
+ Double_t H_cal_1pr_apos_p[THcShTrack::fNrows];
+
+ Double_t H_cal_2ta_aneg_p[THcShTrack::fNrows];
+ Double_t H_cal_2ta_apos_p[THcShTrack::fNrows];
+
+ Double_t H_cal_3ta_aneg_p[THcShTrack::fNrows];
+ Double_t H_cal_3ta_apos_p[THcShTrack::fNrows];
+
+ Double_t H_cal_4ta_aneg_p[THcShTrack::fNrows];
+ Double_t H_cal_4ta_apos_p[THcShTrack::fNrows];
+
+ // Track parameters.
+
+ Double_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.
+
+ fTree->SetBranchAddress("H.cal.1pr.aneg_p",H_cal_1pr_aneg_p);
+ fTree->SetBranchAddress("H.cal.1pr.apos_p",H_cal_1pr_apos_p);
+
+ fTree->SetBranchAddress("H.cal.2ta.aneg_p",H_cal_2ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.2ta.apos_p",H_cal_2ta_apos_p);
+
+ fTree->SetBranchAddress("H.cal.3ta.aneg_p",H_cal_3ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.3ta.apos_p",H_cal_3ta_apos_p);
- // Set a Shower track event from the read raw data.
+ fTree->SetBranchAddress("H.cal.4ta.aneg_p",H_cal_4ta_aneg_p);
+ fTree->SetBranchAddress("H.cal.4ta.apos_p",H_cal_4ta_apos_p);
- UInt_t nhit;
- Double_t p, x,xp, y,yp;
+ fTree->SetBranchAddress("H.cal.trp",&H_cal_trp);
+ fTree->SetBranchAddress("H.cal.trx",&H_cal_trx);
+ fTree->SetBranchAddress("H.cal.trxp",&H_cal_trxp);
+ fTree->SetBranchAddress("H.cal.try",&H_cal_try);
+ fTree->SetBranchAddress("H.cal.tryp",&H_cal_tryp);
- if (fDataStream >> nhit >> p >> x >> xp >> y >> yp) {
- // cout << nhit << " " << p << " " << x << " " << xp << " " << y << " "
- // << yp << endl;
+ fTree->GetEntry(ientry);
- trk.SetTrack(nhit, p, x, xp, y, yp);
+ trk.Reset(H_cal_trp, H_cal_trx, H_cal_trxp, H_cal_try, H_cal_tryp);
- for (UInt_t i=0; i<nhit; i++) {
+ // UInt_t nhit = 0;
+
+ for (UInt_t j=0; j<THcShTrack::fNrows; j++) {
+ for (UInt_t k=0; k<THcShTrack::fNcols; k++) {
Double_t adc_pos, adc_neg;
- UInt_t nb;
- fDataStream >> adc_pos >> adc_neg >> nb;
- // cout << adc_pos << " " << adc_neg << " " << nb << endl;
+ switch (k) {
+ case 0 :
+ adc_pos = H_cal_1pr_apos_p[j];
+ adc_neg = H_cal_1pr_aneg_p[j];
+ break;
+ case 1 :
+ adc_pos = H_cal_2ta_apos_p[j];
+ adc_neg = H_cal_2ta_aneg_p[j];
+ break;
+ case 2 :
+ adc_pos = H_cal_3ta_apos_p[j];
+ adc_neg = H_cal_3ta_aneg_p[j];
+ break;
+ case 3 :
+ adc_pos = H_cal_4ta_apos_p[j];
+ adc_neg = H_cal_4ta_aneg_p[j];
+ break;
+ default:
+ cout << "*** ReadShRawTrack: column number k=" << k
+ << " out of range! ***" << endl;
+ };
+
+ UInt_t nb = j+1 + k*THcShTrack::fNrows;
- trk.AddHit(adc_pos, adc_neg, 0., 0., nb);
- };
- return 1;
+ if (adc_pos>0. || adc_neg>0.) {
+ trk.AddHit(adc_pos, adc_neg, 0., 0., nb);
+ // nhit++;
+ }
+ }
}
- else {
- return 0;
- };
}
@@ -348,92 +427,92 @@ void THcShowerCalib::ComposeVMs() {
// Reset the raw data stream.
- fDataStream.clear() ;
- fDataStream.seekg(0, ios::beg) ;
+ // fDataStream.clear() ;
+ // fDataStream.seekg(0, ios::beg) ;
fNev = 0;
THcShTrack trk;
// Loop over the shower track events in the raw data stream.
- while (ReadShRawTrack(trk)) {
-
- // Check consistency of the track event.
+ // while (ReadShRawTrack(trk)) {
+ for (UInt_t ientry=0; ientry<fNentries; ientry++) {
- if (trk.CheckHitNumber()) {
+ ReadShRawTrack(trk, ientry);
- // Set energy depositions with default gains.
- // Calculate normalized to the track momentum total energy deposition,
- // check it against the thresholds.
+ // 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(falpha0);
+ Double_t Enorm = trk.Enorm();
+ if (Enorm>fLoThr && Enorm<fHiThr) {
- trk.SetEs(falpha1); // Set energies with unit gains for now.
- // trk.Print();
+ trk.SetEs(falpha1); // Set energies with unit gains for now.
+ // trk.Print();
- fe0 += trk.GetP(); // Accumulate track momenta.
+ fe0 += trk.GetP(); // Accumulate track momenta.
- vector<pmt_hit> pmt_hit_list; // Container to save PMT hits
+ vector<pmt_hit> pmt_hit_list; // Container to save PMT hits
- // Loop over hits.
+ // Loop over hits.
- for (UInt_t i=0; i<trk.GetNhits(); i++) {
+ for (UInt_t i=0; i<trk.GetNhits(); i++) {
- THcShHit* hit = trk.GetHit(i);
- // hit->Print();
+ THcShHit* hit = trk.GetHit(i);
+ // hit->Print();
- UInt_t nb = hit->GetBlkNumber();
+ UInt_t nb = hit->GetBlkNumber();
- // Fill the qe and q0 vectors (for positive side PMT).
+ // Fill the qe and q0 vectors (for positive side PMT).
- fqe[nb-1] += hit->GetEpos() * trk.GetP();
- fq0[nb-1] += hit->GetEpos();
+ fqe[nb-1] += hit->GetEpos() * trk.GetP();
+ fq0[nb-1] += hit->GetEpos();
- // Save the PMT hit.
+ // Save the PMT hit.
- pmt_hit_list.push_back( pmt_hit{hit->GetEpos(), nb} );
+ pmt_hit_list.push_back( pmt_hit{hit->GetEpos(), nb} );
- fHitCount[nb-1]++; //Accrue the hit counter.
+ fHitCount[nb-1]++; //Accrue the hit counter.
- // Do the same for the negative side PMTs.
+ // Do the same for the negative side PMTs.
- if (nb <= THcShTrack::fNnegs) {
- fqe[THcShTrack::fNblks+nb-1] += hit->GetEneg() * trk.GetP();
- fq0[THcShTrack::fNblks+nb-1] += hit->GetEneg();
+ 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} );
+ pmt_hit_list.push_back(pmt_hit{hit->GetEneg(),
+ THcShTrack::fNblks+nb} );
- fHitCount[THcShTrack::fNblks+nb-1]++;
- };
+ fHitCount[THcShTrack::fNblks+nb-1]++;
+ };
- } //over hits
+ } //over hits
- // Fill in the correlation matrix Q by retrieving the PMT 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++) {
+ 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;
+ UInt_t ic = (*i).channel;
+ Double_t is = (*i).signal;
- for (vector<pmt_hit>::iterator j=i;
- j < pmt_hit_list.end(); j++) {
+ for (vector<pmt_hit>::iterator j=i;
+ j < pmt_hit_list.end(); j++) {
- UInt_t jc = (*j).channel;
- Double_t js = (*j).signal;
+ 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;
- }
+ fQ[ic-1][jc-1] += is*js;
+ if (jc != ic) fQ[jc-1][ic-1] += is*js;
}
+ }
- fNev++;
- };
- };
- };
+ fNev++;
+
+ }; // if within thresholds
+
+ }; // over entries
// Take averages.
@@ -627,8 +706,8 @@ void THcShowerCalib::FillHEcal() {
//
// Reset.
- fDataStream.clear() ;
- fDataStream.seekg(0, ios::beg) ;
+ // fDataStream.clear() ;
+ // fDataStream.seekg(0, ios::beg) ;
ofstream output;
output.open("calibrated.d",ios::out);
@@ -637,8 +716,10 @@ void THcShowerCalib::FillHEcal() {
THcShTrack trk;
- while (ReadShRawTrack(trk)) {
+ // while (ReadShRawTrack(trk)) {
+ for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ ReadShRawTrack(trk, ientry);
// trk.Print();
trk.SetEs(falphaC); // use 'constrained' calibration constants
diff --git a/hcal_calib/db_run.dat b/hcal_calib/db_run.dat
new file mode 100755
index 0000000000000000000000000000000000000000..6dd115d7003b81a3ee817f4971d289288c794397
--- /dev/null
+++ b/hcal_calib/db_run.dat
@@ -0,0 +1,7 @@
+# Test run database
+
+# DAQ04
+--------[ 2000-01-01 01:00:00 ]
+
+#A1 2202 pedestal
+ebeam = 4.02187
diff --git a/hcal_calib/hcal_calib.cpp b/hcal_calib/hcal_calib.cpp
index c73c5c37a972078bf52711f4279b6ee6831bf0d5..4798bbb258c7902714c8232582b4237c26cef232 100644
--- a/hcal_calib/hcal_calib.cpp
+++ b/hcal_calib/hcal_calib.cpp
@@ -11,7 +11,7 @@ void hcal_calib(Int_t RunNumber) {
THcShowerCalib theShowerCalib(RunNumber);
- theShowerCalib.ExtractData(); // Extract data from the Root file
+ // theShowerCalib.ExtractData(); // Extract data from the Root file
theShowerCalib.Init(); // Initialize constants adn variables
theShowerCalib.CalcThresholds(); // Thresholds on the uncalibrated Edep/P
theShowerCalib.ComposeVMs(); // Compute vectors amd matrices for calib.
@@ -20,7 +20,7 @@ void hcal_calib(Int_t RunNumber) {
theShowerCalib.SaveAlphas(); // Save the constants
// Plot histograms
-
+
TCanvas* Canvas =
new TCanvas("Canvas", "HMS Shower Counter calibration", 1000, 667);
Canvas->Divide(2,2);
diff --git a/hcal_calib/make_cratemap.pl b/hcal_calib/make_cratemap.pl
new file mode 100755
index 0000000000000000000000000000000000000000..6ee8d5e3a2cbbe1577221886cc445dbfe63272b9
--- /dev/null
+++ b/hcal_calib/make_cratemap.pl
@@ -0,0 +1,70 @@
+#!/usr/bin/perl
+
+# Read a Hall C style MAP file and output a
+# Hall A style crate map DB file.
+#
+# 22.03.2012 (saw)
+
+%crates=();
+
+$crate = 0;
+$nsubadd = 0;
+$bsub = 0;
+$modtype = 0;
+$slot = 0;
+while(<>) {
+ $line=chomp;
+ if($line=/^\s*ROC=\s*(\d*)/i) {
+ $i++;
+ $crate = $1;
+ if(not $crates{$crate}) {
+ $slotlist={};
+ $crates{$crate} = $slotlist;
+ }
+ $modtype = 0;
+ $slot = 0;
+ } elsif ($line=/^\s*nsubadd=\s*(\d*)/i) {
+ $nsubadd = $1;
+ $modtype = 0;
+ } elsif ($line=/^\s*bsub=\s*(\d*)/i) {
+ $bsub = $1;
+ $modtype = 0;
+ } elsif ($line=/^\s*slot=\s*(\d*)/i) {
+ $slot = $1;
+ $modtype = 0;
+ } elsif ($line=/^\s*(\d*)\s*,\s*(\d*)\s*,\s*(\d*)/) {
+ if($modtype == 0) { # Slot not yet registered
+ if($nsubadd == 96) {
+ $modtype = 1877;
+ } elsif($nsubadd == 64) {
+ if($bsub == 16) {
+ $modtype = 1875;
+ } elsif($bsub == 17) {
+ $modtype = 1881;
+ }
+ }
+ if($modtype == 0) {
+ print "Unknown module Crate $crate, Slot $slot\n";
+ }
+ $crates{$crate}{$slot} = $modtype;
+ # print "$crate $slot $modtype\n";
+ }
+ }
+}
+print "# Hall C Crate map\n";
+foreach $crate (sort {$a <=> $b} keys %crates) {
+ print "==== Crate $crate type fastbus\n";
+ print "# slot model clear header mask nchan ndata\n";
+ foreach $slot (sort {$a <=> $b} keys %{ $crates{$crate}}) {
+ $modtype = $crates{$crate}{$slot};
+ if($modtype == 1877) {
+ $ndata = 256;
+ } else {
+ $ndata = 64;
+ }
+ printf " %2d %d 1 0x0 0x0 %3d %d\n"
+ ,$slot,$modtype,$nsubadd, $ndata;
+ }
+}
+
+