Skip to content
Snippets Groups Projects
Commit bd9f164f authored by Mark K Jones's avatar Mark K Jones Committed by GitHub
Browse files

Merge pull request #155 from vardant/develop 

Add calibration code for the SHMS aerogel detector
parents cb1620e2 6403cc21
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CALIBRATION/shms_aero_calib/paero_calib.C 0 → 100644
+ 107
− 0
View file @ bd9f164f
#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();
}
CALIBRATION/shms_cal_calib/THcPShowerCalib.h
+ 27
− 13
View file @ bd9f164f
...@@ -53,7 +53,7 @@ class THcPShowerCalib { ...@@ -53,7 +53,7 @@ class THcPShowerCalib {
public: public:
THcPShowerCalib(string); THcPShowerCalib(string, int, int);
THcPShowerCalib(); THcPShowerCalib();
~THcPShowerCalib(); ~THcPShowerCalib();
...@@ -81,6 +81,8 @@ class THcPShowerCalib { ...@@ -81,6 +81,8 @@ class THcPShowerCalib {
TTree* fTree; TTree* fTree;
UInt_t fNentries; UInt_t fNentries;
UInt_t fNstart;
UInt_t fNstop;
// Declaration of leaves types // Declaration of leaves types
...@@ -150,8 +152,10 @@ THcPShowerCalib::THcPShowerCalib() {}; ...@@ -150,8 +152,10 @@ THcPShowerCalib::THcPShowerCalib() {};
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
THcPShowerCalib::THcPShowerCalib(string RunNumber) { THcPShowerCalib::THcPShowerCalib(string RunNumber, int nstart, int nstop) {
fRunNumber = RunNumber; fRunNumber = RunNumber;
fNstart = nstart;
fNstop = nstop;
}; };
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
...@@ -173,7 +177,9 @@ void THcPShowerCalib::SaveRawData() { ...@@ -173,7 +177,9 @@ void THcPShowerCalib::SaveRawData() {
THcPShTrack trk; 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)) { if (ReadShRawTrack(trk, ientry)) {
trk.SetEs(falphaC); trk.SetEs(falphaC);
...@@ -281,12 +287,16 @@ void THcPShowerCalib::CalcThresholds() { ...@@ -281,12 +287,16 @@ void THcPShowerCalib::CalcThresholds() {
// Histogram uncalibrated energy depositions, get mean and RMS from the // Histogram uncalibrated energy depositions, get mean and RMS from the
// histogram, establish +/-3 * RMS thresholds. // 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; Int_t nev = 0;
THcPShTrack trk; 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)) { if ( ReadShRawTrack(trk, ientry)) {
...@@ -425,7 +435,9 @@ void THcPShowerCalib::ComposeVMs() { ...@@ -425,7 +435,9 @@ void THcPShowerCalib::ComposeVMs() {
// Loop over the shower track events in the ntuples. // 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)) { if (ReadShRawTrack(trk, ientry)) {
...@@ -706,14 +718,16 @@ void THcPShowerCalib::FillHEcal() { ...@@ -706,14 +718,16 @@ void THcPShowerCalib::FillHEcal() {
// Output event by event energy depositions and momenta for debug purposes. // Output event by event energy depositions and momenta for debug purposes.
// //
ofstream output; // ofstream output;
output.open("calibrated.deb",ios::out); // output.open("calibrated.deb",ios::out);
Int_t nev = 0; Int_t nev = 0;
THcPShTrack trk; 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)) { if (ReadShRawTrack(trk, ientry)) {
// trk.Print(cout); // trk.Print(cout);
...@@ -729,15 +743,15 @@ void THcPShowerCalib::FillHEcal() { ...@@ -729,15 +743,15 @@ void THcPShowerCalib::FillHEcal() {
hESHvsEPR->Fill(trk.EPRnorm(), trk.ESHnorm()); hESHvsEPR->Fill(trk.EPRnorm(), trk.ESHnorm());
output << Enorm*P/1000. << " " << P/1000. << " " << delta << " " // output << Enorm*P/1000. << " " << P/1000. << " " << delta << " "
<< trk.GetX() << " " << trk.GetY() << endl; // << trk.GetX() << " " << trk.GetY() << endl;
nev++; nev++;
} }
}; };
output.close(); // output.close();
cout << "FillHEcal: " << nev << " events filled" << endl; cout << "FillHEcal: " << nev << " events filled" << endl;
}; };
...@@ -752,7 +766,7 @@ void THcPShowerCalib::SaveAlphas() { ...@@ -752,7 +766,7 @@ void THcPShowerCalib::SaveAlphas() {
// //
ofstream output; 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; cout << "SaveAlphas: fname=" << fname << endl;
output.open(fname,ios::out); output.open(fname,ios::out);
......
CALIBRATION/shms_cal_calib/pcal_calib.cpp
+ 20
− 5
View file @ bd9f164f
...@@ -8,14 +8,15 @@ ...@@ -8,14 +8,15 @@
// A steering Root script for the SHMS calorimeter calibration. // 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 // Initialize the analysis clock
clock_t t = 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.Init(); // Initialize constants and variables
theShowerCalib.CalcThresholds(); // Thresholds on the uncalibrated Edep/P theShowerCalib.CalcThresholds(); // Thresholds on the uncalibrated Edep/P
...@@ -41,7 +42,7 @@ void pcal_calib(string RunNumber) { ...@@ -41,7 +42,7 @@ void pcal_calib(string RunNumber) {
theShowerCalib.hEuncSel->DrawCopy("same"); theShowerCalib.hEuncSel->DrawCopy("same");
Canvas->cd(2); Canvas->cd(2);
theShowerCalib.hESHvsEPR->Draw("colz"); theShowerCalib.hESHvsEPR->Draw();
// Normalized energy deposition after calibration. // Normalized energy deposition after calibration.
...@@ -56,7 +57,21 @@ void pcal_calib(string RunNumber) { ...@@ -56,7 +57,21 @@ void pcal_calib(string RunNumber) {
// SHMS delta(P) versus the calibrated energy deposition. // SHMS delta(P) versus the calibrated energy deposition.
Canvas->cd(4); 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 // Calculate the analysis rate
t = clock() - t; t = clock() - t;
......
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