diff --git a/CALIBRATION/shms_aero_calib/paero_calib.C b/CALIBRATION/shms_aero_calib/paero_calib.C
new file mode 100644
index 0000000000000000000000000000000000000000..6a7e864ea1440cb1d891359a7d7225b832c1b31f
--- /dev/null
+++ b/CALIBRATION/shms_aero_calib/paero_calib.C
@@ -0,0 +1,107 @@
+#include <TFile.h>
+#include <TTree.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+
+//
+// Calibrate SHMS aerogel detector by localizing SPE peaks.
+// The SPE peak positions are found from Gaussian fits to the peaks.
+// The fit ranges can be adjusted by changing fit range limits
+// flo_* and fhi_* in the body of code.
+
+// The code works on root output from hcana, the name of which is the
+// single input parameter of the program.
+
+#define NPMT 7
+#define MaxAdc 2000.
+#define NBin 1000
+
+void paero_calib(string fname) {
+
+ TFile *f = new TFile(fname.c_str());
+ TTree* tree;
+ f->GetObject("T",tree);
+
+ int nentries = tree->GetEntries();
+ cout << "nentries= " << nentries << endl;
+
+ double adc_pos[NPMT];
+ double adc_neg[NPMT];
+
+ TBranch* b_adc_pos;
+ TBranch* b_adc_neg;
+
+ tree->SetBranchAddress("P.aero.goodPosAdcPulseInt",adc_pos,&b_adc_pos);
+ tree->SetBranchAddress("P.aero.goodNegAdcPulseInt",adc_neg,&b_adc_neg);
+
+ TH1D* hpos[NPMT];
+ TH1D* hneg[NPMT];
+
+ for (int i=0; i<NPMT; i++) {
+ hpos[i] = new TH1D(Form("hpos%d",i+1), Form("ADC+ %i",i+1), NBin, 0.001, MaxAdc);
+ hneg[i] = new TH1D(Form("hneg%d",i+1), Form("ADC- %i",i+1), NBin, 0.001, MaxAdc);
+ }
+
+ for (int ientry=0; ientry<nentries; ientry++) {
+
+ tree->GetEntry(ientry);
+
+ for (int i=0; i<NPMT; i++) {
+ hpos[i]->Fill(adc_pos[i], 1.);
+ hneg[i]->Fill(adc_neg[i], 1.);
+ }
+
+ }
+
+ TCanvas* c1 = new TCanvas("adc_spec", "fADC spectra", 600, 800);
+ c1->Divide(2,NPMT);
+
+ double flo_neg[NPMT] {400., 400., 400., 400., 350., 400., 400.};
+ double fhi_neg[NPMT] {700., 650., 800., 700., 750., 650., 900.};
+
+ double flo_pos[NPMT] {250., 300., 400., 550., 250., 300., 250.};
+ double fhi_pos[NPMT] {500., 600., 700., 800., 500., 600., 550.};
+
+ float gain_pos[NPMT] {NPMT*0.};
+ float gain_neg[NPMT] {NPMT*0.};
+
+ int ip=0;
+ for (int i=0; i<NPMT; i++) {
+ c1->cd(++ip);
+ if (hpos[i]->GetSumOfWeights() > 0) {
+ hpos[i]->Fit("gaus","","",flo_pos[i],fhi_pos[i]);
+ hpos[i]->GetFunction("gaus")->SetLineColor(2);
+ hpos[i]->GetFunction("gaus")->SetLineWidth(2);
+ gain_pos[i] = hpos[i]->GetFunction("gaus")->GetParameter(1);
+ }
+ else
+ hpos[i]->Draw();
+ c1->cd(++ip);
+ if (hneg[i]->GetSumOfWeights() > 0) {
+ hneg[i]->Fit("gaus","","",flo_neg[i],fhi_neg[i]);
+ hneg[i]->GetFunction("gaus")->SetLineColor(2);
+ hneg[i]->GetFunction("gaus")->SetLineWidth(2);
+ gain_neg[i] = hneg[i]->GetFunction("gaus")->GetParameter(1);
+ }
+ else
+ hneg[i]->Draw();
+ }
+
+ ofstream of;
+ of.open("gain.r",ios::out);
+ of << "paero_neg_gain = ";
+ for (int i=0; i<NPMT; i++)
+ of << "1./" << gain_neg[i] << ", ";
+ of << endl;
+ of << "paero_pos_gain = ";
+ for (int i=0; i<NPMT; i++)
+ of << "1./" << gain_pos[i] << ", ";
+ of << endl;
+ of.close();
+
+}
diff --git a/CALIBRATION/shms_cal_calib/THcPShowerCalib.h b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
index 4200a0c24d18ddb575e305b6c2f046071f76c70e..9e1cb0694b16e6780c8af586b087c5c496f75221 100644
--- a/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
+++ b/CALIBRATION/shms_cal_calib/THcPShowerCalib.h
@@ -53,7 +53,7 @@ class THcPShowerCalib {
public:
- THcPShowerCalib(string);
+ THcPShowerCalib(string, int, int);
THcPShowerCalib();
~THcPShowerCalib();
@@ -81,6 +81,8 @@ class THcPShowerCalib {
TTree* fTree;
UInt_t fNentries;
+ UInt_t fNstart;
+ UInt_t fNstop;
// Declaration of leaves types
@@ -150,8 +152,10 @@ THcPShowerCalib::THcPShowerCalib() {};
//------------------------------------------------------------------------------
-THcPShowerCalib::THcPShowerCalib(string RunNumber) {
+THcPShowerCalib::THcPShowerCalib(string RunNumber, int nstart, int nstop) {
fRunNumber = RunNumber;
+ fNstart = nstart;
+ fNstop = nstop;
};
//------------------------------------------------------------------------------
@@ -173,7 +177,9 @@ void THcPShowerCalib::SaveRawData() {
THcPShTrack trk;
- for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ // for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ for (UInt_t ientry=TMath::Max(UInt_t(0),fNstart);
+ ientry<TMath::Min(fNstop,fNentries); ientry++) {
if (ReadShRawTrack(trk, ientry)) {
trk.SetEs(falphaC);
@@ -281,12 +287,16 @@ void THcPShowerCalib::CalcThresholds() {
// Histogram uncalibrated energy depositions, get mean and RMS from the
// histogram, establish +/-3 * RMS thresholds.
- cout << "THcPShowerCalib::CalcThresholds: FNentries = " << fNentries << endl;
+ //cout<< "THcPShowerCalib::CalcThresholds: FNentries = " << fNentries << endl;
+ cout << "THcPShowerCalib::CalcThresholds: fNstart = " << fNstart << " "
+ << " fNstop = " << fNstop << endl;
Int_t nev = 0;
THcPShTrack trk;
- for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ // for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ for (UInt_t ientry=TMath::Max(UInt_t(0),fNstart);
+ ientry<TMath::Min(fNstop,fNentries); ientry++) {
if ( ReadShRawTrack(trk, ientry)) {
@@ -425,7 +435,9 @@ void THcPShowerCalib::ComposeVMs() {
// Loop over the shower track events in the ntuples.
- for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ // for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ for (UInt_t ientry=TMath::Max(UInt_t(0),fNstart);
+ ientry<TMath::Min(fNstop,fNentries); ientry++) {
if (ReadShRawTrack(trk, ientry)) {
@@ -706,14 +718,16 @@ void THcPShowerCalib::FillHEcal() {
// Output event by event energy depositions and momenta for debug purposes.
//
- ofstream output;
- output.open("calibrated.deb",ios::out);
+ // ofstream output;
+ // output.open("calibrated.deb",ios::out);
Int_t nev = 0;
THcPShTrack trk;
- for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ // for (UInt_t ientry=0; ientry<fNentries; ientry++) {
+ for (UInt_t ientry=TMath::Max(UInt_t(0),fNstart);
+ ientry<TMath::Min(fNstop,fNentries); ientry++) {
if (ReadShRawTrack(trk, ientry)) {
// trk.Print(cout);
@@ -729,15 +743,15 @@ void THcPShowerCalib::FillHEcal() {
hESHvsEPR->Fill(trk.EPRnorm(), trk.ESHnorm());
- output << Enorm*P/1000. << " " << P/1000. << " " << delta << " "
- << trk.GetX() << " " << trk.GetY() << endl;
+ // output << Enorm*P/1000. << " " << P/1000. << " " << delta << " "
+ // << trk.GetX() << " " << trk.GetY() << endl;
nev++;
}
};
- output.close();
+ // output.close();
cout << "FillHEcal: " << nev << " events filled" << endl;
};
@@ -752,7 +766,7 @@ void THcPShowerCalib::SaveAlphas() {
//
ofstream output;
- char* fname = Form("pcal.param.%s",fRunNumber.c_str());
+ char* fname = Form("pcal.param.%s_%d-%d",fRunNumber.c_str(),fNstart,fNstop);
cout << "SaveAlphas: fname=" << fname << endl;
output.open(fname,ios::out);
diff --git a/CALIBRATION/shms_cal_calib/pcal_calib.cpp b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
index 008e7b1eea8966ea4af6baa453c6db6e73a2e9e4..0411401802208ce5e4a432e05af97317364c93ac 100644
--- a/CALIBRATION/shms_cal_calib/pcal_calib.cpp
+++ b/CALIBRATION/shms_cal_calib/pcal_calib.cpp
@@ -8,14 +8,15 @@
// A steering Root script for the SHMS calorimeter calibration.
//
-void pcal_calib(string RunNumber) {
+void pcal_calib(string RunNumber, int nstart=0, int nstop=999999999) {
// Initialize the analysis clock
clock_t t = clock();
- cout << "Calibrating run " << RunNumber << endl;
+ cout << "Calibrating run " << RunNumber << ", events "
+ << nstart << " -- " << nstop << endl;
- THcPShowerCalib theShowerCalib(RunNumber);
+ THcPShowerCalib theShowerCalib(RunNumber, nstart, nstop);
theShowerCalib.Init(); // Initialize constants and variables
theShowerCalib.CalcThresholds(); // Thresholds on the uncalibrated Edep/P
@@ -41,7 +42,7 @@ void pcal_calib(string RunNumber) {
theShowerCalib.hEuncSel->DrawCopy("same");
Canvas->cd(2);
- theShowerCalib.hESHvsEPR->Draw("colz");
+ theShowerCalib.hESHvsEPR->Draw();
// Normalized energy deposition after calibration.
@@ -56,7 +57,21 @@ void pcal_calib(string RunNumber) {
// SHMS delta(P) versus the calibrated energy deposition.
Canvas->cd(4);
- theShowerCalib.hDPvsEcal->Draw("colz");
+ theShowerCalib.hDPvsEcal->Draw();
+
+ // Save canvas in a pdf format.
+ Canvas->Print(Form("%s_%d-%d.pdf",RunNumber.c_str(),nstart,nstop));
+
+ // Save histograms in root file.
+
+ //TFile* froot=new TFile(Form("%s_%d-%d.root",RunNumber.c_str(),nstart,nstop),
+ // "RECREATE");
+ // theShowerCalib.hEunc->Write();
+ // theShowerCalib.hEuncSel->Write();
+ // theShowerCalib.hESHvsEPR->Write();
+ // theShowerCalib.hEcal->Write();
+ // theShowerCalib.hDPvsEcal->Write();
+ // froot->Close();
// Calculate the analysis rate
t = clock() - t;