Skip to content
Snippets Groups Projects
Commit 082b78d7 authored by Ryan-Ambrose's avatar Ryan-Ambrose Committed by Eric Pooser
Browse files

Devel (#288) 

* Calibration script now uses normalized energy for particle ID and can update calibration constants directly. Both scripts use PROOF-lite by default. Added a run_cal script for ease of use. README now reflects updated scripts.

* Calibration/running script cleanup

* Updating scripts to be compatible with new PRODUCTION replay

* Altered scripts to be compatible with new SI units and changed axis labels.

* Updated NGC & HGC calibration to align with new SI units.
parent e25a8a78
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CALIBRATION/shms_hgcer_calib/calibration.C
+ 46
− 43
View file @ 082b78d7
......@@ -99,16 +99,16 @@ void calibration::SlaveBegin(TTree * /*tree*/)
if (fNGC) //Set up histograms for NGC
{
ADC_min = 0;
ADC_max = 12000;
bins = (abs(ADC_min) + abs(ADC_max));
ADC_min = -10;
ADC_max = 250;
bins = 2*(abs(ADC_min) + abs(ADC_max));
}
if (!fNGC) //Set up histograms for HGC
{
ADC_min = -4000;
ADC_max = 12000;
bins = abs(ADC_min) + abs(ADC_max);
ADC_min = -10;
ADC_max = 200;
bins = 2*(abs(ADC_min) + abs(ADC_max));
}
fPulseInt = new TH1F*[4];
......@@ -119,11 +119,11 @@ void calibration::SlaveBegin(TTree * /*tree*/)
fPulseInt_quad[3] = new TH1F*[4];
for (Int_t ipmt = 0; ipmt < 4; ipmt++)
{
fPulseInt[ipmt] = new TH1F(Form("PulseInt_PMT%d",ipmt+1),Form("Pulse Integral PMT%d; ADC Channel; Counts",ipmt+1), bins, ADC_min, ADC_max);
fPulseInt[ipmt] = new TH1F(Form("PulseInt_PMT%d",ipmt+1),Form("Pulse Integral PMT%d; ADC Channel (pC); Counts",ipmt+1), bins, ADC_min, ADC_max);
GetOutputList()->Add(fPulseInt[ipmt]);
for (Int_t iquad = 0; iquad < 4; iquad++)
{
fPulseInt_quad[iquad][ipmt] = new TH1F(Form("PulseInt_quad%d_PMT%d",iquad+1,ipmt+1),Form("Pulse Integral PMT%d quad%d; ADC Channel; Counts",ipmt+1,iquad+1),bins,ADC_min,ADC_max);
fPulseInt_quad[iquad][ipmt] = new TH1F(Form("PulseInt_quad%d_PMT%d",iquad+1,ipmt+1),Form("Pulse Integral PMT%d quad%d; ADC Channel (pC); Counts",ipmt+1,iquad+1),bins,ADC_min,ADC_max);
GetOutputList()->Add(fPulseInt_quad[iquad][ipmt]);
}
}
......@@ -185,8 +185,8 @@ Bool_t calibration::Process(Long64_t entry)
{
//Perform a loose timing cut
if (!fFullRead) fNGC ? b_P_ngcer_goodAdcPulseTime->GetEntry(entry) : b_P_hgcer_goodAdcPulseTime->GetEntry(entry);
if (fNGC ? P_ngcer_goodAdcPulseTime[ipmt] < 1000 || P_ngcer_goodAdcPulseTime[ipmt] > 2000 :
P_hgcer_goodAdcPulseTime[ipmt] < 1000 || P_hgcer_goodAdcPulseTime[ipmt] > 2000) continue;
if (fNGC ? P_ngcer_goodAdcPulseTime[ipmt] < 50 || P_ngcer_goodAdcPulseTime[ipmt] > 125 :
P_hgcer_goodAdcPulseTime[ipmt] < 70 || P_hgcer_goodAdcPulseTime[ipmt] > 135) continue;
//Cuts to remove entries corresponding to a PMT not registering a hit
if (!fFullRead) fNGC ? b_P_ngcer_goodAdcPulseInt->GetEntry(entry) : b_P_hgcer_goodAdcPulseInt->GetEntry(entry);
......@@ -467,7 +467,9 @@ void calibration::Terminate()
}
}
//Rebin the histograms into something more sensible, add functionality to bin HGC & NGC independently
//Rebin the histograms, add functionality to bin HGC & NGC independently
//Not needed since unit conversion into SI, but available in the future
/*
for (Int_t ipmt=0; ipmt < (fNGC ? fngc_pmts : fhgc_pmts); ipmt++)
{
for (Int_t iquad=0; iquad<4; iquad++)
......@@ -476,7 +478,7 @@ void calibration::Terminate()
}
fNGC ? PulseInt[ipmt]->Rebin(20) : PulseInt[ipmt]->Rebin(20);
}
*/
//Show the particle cuts performed in the histogram forming
if (fCut)
......@@ -492,7 +494,7 @@ void calibration::Terminate()
gStyle->SetOptFit(111);
//Single Gaussian to find mean of SPE
TF1 *Gauss1 = new TF1("Gauss1",gauss,-500,7000,3);
TF1 *Gauss1 = new TF1("Gauss1",gauss,-10,200,3);
Gauss1->SetParNames("Amplitude","Mean","Std. Dev.");
//Sum of three Gaussians to determine NPE spacing
......@@ -551,23 +553,23 @@ void calibration::Terminate()
if (fFullShow) quad_cuts_ipmt->cd(ipad);
//Perform search for the SPE and save the peak into the array xpeaks
fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 2.5, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 2.5, "nobackground&&nodraw", 0.001);
fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 2.0, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 2.5, "nobackground&&nodraw", 0.001);
TList *functions = PulseInt_quad[iquad][ipmt]->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
Double_t *xpeaks = pm->GetX();
if (xpeaks[1] < xpeaks[0]) xpeaks[1] = xpeaks[0];
//Use the peak to fit the SPE with a Gaussian to determine the mean
Gauss1->SetRange(xpeaks[1]-150, xpeaks[1]+150);
Gauss1->SetRange(xpeaks[1]-5, xpeaks[1]+5);
Gauss1->SetParameter(1, xpeaks[1]);
Gauss1->SetParameter(2, 200.);
Gauss1->SetParameter(2, 10.);
Gauss1->SetParLimits(0, 0., 2000.);
Gauss1->SetParLimits(1, xpeaks[1]-150, xpeaks[1]+150);
Gauss1->SetParLimits(2, 10., 500.);
Gauss1->SetParLimits(1, xpeaks[1]-10, xpeaks[1]+10);
Gauss1->SetParLimits(2, 0.5, 10.);
fFullShow ? PulseInt_quad[iquad][ipmt]->Fit("Gauss1","RQ") : PulseInt_quad[iquad][ipmt]->Fit("Gauss1","RQN");
//Store the mean of the SPE in the mean array provided it is not zero and passes a loose statistical cut. Note that indexing by ipad-1 is for convienience
if (xpeaks[1] != 0.0 && PulseInt_quad[iquad][ipmt]->GetBinContent(PulseInt_quad[iquad][ipmt]->GetXaxis()->FindBin(xpeaks[1])) > 90) mean[ipad-1] = Gauss1->GetParameter(1);
if (xpeaks[1] > 2.0 && PulseInt_quad[iquad][ipmt]->GetBinContent(PulseInt_quad[iquad][ipmt]->GetXaxis()->FindBin(xpeaks[1])) > 90) mean[ipad-1] = Gauss1->GetParameter(1);
ipad++;
}
......@@ -583,23 +585,23 @@ void calibration::Terminate()
if (num_peaks != 0.0) xscale = xscale/num_peaks;
//Perform check if the statistics were too low to get a good estimate of the SPE mean
if (xscale < 10.0)
if (xscale < 1.0)
{
//Repeat the exact same procedure for the SPE of each quadrant, except now its for the full PMT spectra
if(fFullShow) low_stats_ipmt = new TCanvas(Form("low_stats_%d",ipmt),Form("Low stats analysis for PMT%d",ipmt+1));
if(fFullShow) low_stats_ipmt->cd(1);
PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,1000);
PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,20);
fFullShow ? s->Search(PulseInt[ipmt], 3.5, "nobackground", 0.001) : s->Search(PulseInt[ipmt], 2.0, "nobackground&&nodraw", 0.001);
TList *functions = PulseInt[ipmt]->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
Double_t *xpeaks = pm->GetX();
Gauss1->SetRange(xpeaks[1]-100, xpeaks[1]+100);
Gauss1->SetRange(xpeaks[1]-5, xpeaks[1]+5);
Gauss1->SetParameter(1, xpeaks[1]);
Gauss1->SetParameter(2, 200.);
Gauss1->SetParameter(2, 10.);
Gauss1->SetParLimits(0, 0., 2000.);
Gauss1->SetParLimits(1, xpeaks[1]-50, xpeaks[1]+50);
Gauss1->SetParLimits(2, 10., 500.);
PulseInt[ipmt]->GetXaxis()->SetRangeUser(-500,7000);
Gauss1->SetParLimits(1, xpeaks[1]-5, xpeaks[1]+5);
Gauss1->SetParLimits(2, 0.5, 20.);
PulseInt[ipmt]->GetXaxis()->SetRangeUser(-10,200);
fFullShow ? PulseInt[ipmt]->Fit("Gauss1","RQ") : PulseInt[ipmt]->Fit("Gauss1","RQN");
xscale = Gauss1->GetParameter(1);
}
......@@ -608,7 +610,7 @@ void calibration::Terminate()
nbins = (PulseInt[ipmt]->GetXaxis()->GetNbins());
//With the scale of ADC to NPE create a histogram that has the conversion applied
fscaled[ipmt] = new TH1F(Form("fscaled_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 360 : 280, 0, fNGC ? 25 : 20);
fscaled[ipmt] = new TH1F(Form("fscaled_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 210 : 210, -1, fNGC ? 20 : 20);
//Fill this histogram bin by bin
for (Int_t ibin=0; ibin<nbins; ibin++)
......@@ -632,7 +634,7 @@ void calibration::Terminate()
fFullShow ? fscaled[ipmt]->Fit("Poisson","RQ") : fscaled[ipmt]->Fit("Poisson","RQN");
//Make and fill histogram with the background removed
fscaled_nobackground[ipmt] = new TH1F(Form("fscaled_nobackground_pmt%d", ipmt+1), Form("NPE spectra background removed for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 360 : 280, 0, fNGC ? 25 : 20);
fscaled_nobackground[ipmt] = new TH1F(Form("fscaled_nobackground_pmt%d", ipmt+1), Form("NPE spectra background removed for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 210 : 210, -1, fNGC ? 20 : 20);
for (Int_t ibin=0; ibin<nbins; ibin++)
{
......@@ -679,7 +681,7 @@ void calibration::Terminate()
Double_t xscale_mk2 = xscale * Gauss3->GetParameter(1);
//Take this new xscale and repeat the exact same procedure as before
fscaled_mk2[ipmt] = new TH1F(Form("fhgc_scaled_mk2_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 360 : 280, 0, fNGC ? 25 : 20);
fscaled_mk2[ipmt] = new TH1F(Form("fhgc_scaled_mk2_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 210 : 210, -1, fNGC ? 20 : 20);
//Fill this histogram bin by bin
for (Int_t ibin=0; ibin<nbins; ibin++)
......@@ -703,7 +705,7 @@ void calibration::Terminate()
fFullShow ? fscaled_mk2[ipmt]->Fit("Poisson","RQ"):fscaled_mk2[ipmt]->Fit("Poisson","RQN");
//Make and fill histogram with the background removed
fscaled_mk2_nobackground[ipmt] = new TH1F(Form("fscaled_mk2_nobackground_pmt%d", ipmt+1), Form("NPE spectra background removed for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 360 : 280, 0, fNGC ? 25 : 20);
fscaled_mk2_nobackground[ipmt] = new TH1F(Form("fscaled_mk2_nobackground_pmt%d", ipmt+1), Form("NPE spectra background removed for PMT%d; NPE; Normalized Counts",ipmt+1), fNGC ? 210 : 210, -1, fNGC ? 20 : 20);
for (Int_t ibin=0; ibin<nbins; ibin++)
{
......@@ -764,24 +766,25 @@ void calibration::Terminate()
//Perform search for the SPE and save the peak into the array xpeaks
if (fFullShow) quad_cuts_ipmt->cd(iquad+1);
fNGC ? PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(150,2000) : PulseInt_quad[ipmt][ipmt]->GetXaxis()->SetRangeUser(150,600);
fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 1.5, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 1.5, "nobackground&&nodraw", 0.001);
//fNGC ? PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(150,2000) : PulseInt_quad[ipmt][ipmt]->GetXaxis()->SetRangeUser(5,20);
fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 2.0, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 2.0, "nobackground&&nodraw", 0.001);
TList *functions = PulseInt_quad[iquad][ipmt]->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
Double_t *xpeaks = pm->GetX();
//Use the peak to fit the SPE with a Gaussian to determine the mean
Gauss1->SetRange(xpeaks[0]-150, xpeaks[0]+400);
Gauss1->SetParameter(1, xpeaks[0]);
Gauss1->SetParameter(2, 200.);
Gauss1->SetRange(xpeaks[1]-3, xpeaks[1]+3);
Gauss1->SetParameter(1, xpeaks[1]);
Gauss1->SetParameter(2, 10.);
Gauss1->SetParLimits(0, 0., 2000.);
Gauss1->SetParLimits(1, xpeaks[0]-100, xpeaks[0]+200);
Gauss1->SetParLimits(2, 10., 500.);
PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(-500,12000);
Gauss1->SetParLimits(1, xpeaks[1]-3, xpeaks[1]+3);
Gauss1->SetParLimits(2, 0.5, 10.);
//PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(5,20);
fFullShow ? PulseInt_quad[iquad][ipmt]->Fit("Gauss1","RQ") : PulseInt_quad[iquad][ipmt]->Fit("Gauss1","RQN");
//Store the mean of the SPE in the mean array provided it is not zero, passes a loose statistical cut, and is above a minimum channel number
if (xpeaks[0] != 0.0 && Gauss1->GetParameter(1) > 300.) mean[iquad] = Gauss1->GetParameter(1);
//Added condition that iquad != ipmt since spectrum is quite different
if (xpeaks[1] != 0.0 && PulseInt_quad[iquad][ipmt]->GetBinContent(PulseInt_quad[iquad][ipmt]->GetXaxis()->FindBin(xpeaks[1])) > 50 && iquad != ipmt) mean[iquad] = Gauss1->GetParameter(1);
}
Double_t xscale = 0.0;
......@@ -799,7 +802,7 @@ void calibration::Terminate()
Int_t nbins;
nbins = PulseInt_quad[ipmt][ipmt]->GetXaxis()->GetNbins();
fscaled[ipmt] = new TH1F(Form("fscaled_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), 300, 0, fNGC ? 30 : 20);
fscaled[ipmt] = new TH1F(Form("fscaled_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), 210, -1, fNGC ? 20 : 20);
//Fill this histogram bin by bin
for (Int_t ibin=0; ibin<nbins; ibin++)
......@@ -831,7 +834,7 @@ void calibration::Terminate()
pmt_calib[ipmt] = abs(1.0 - Gauss1->GetParameter(1));
//Scale full ADC spectra according to the mean of the SPE. This requires filling a new histogram with the same number of bins but scaled min/max
fscaled_mk2[ipmt] = new TH1F(Form("fscaled_mk2_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), 300, 0, fNGC ? 30 : 20);
fscaled_mk2[ipmt] = new TH1F(Form("fscaled_mk2_PMT%d", ipmt+1), Form("Scaled ADC spectra for PMT%d; NPE; Normalized Counts",ipmt+1), 210, -1, fNGC ? 20 : 20);
//Fill this histogram bin by bin
for (Int_t ibin=0; ibin<nbins; ibin++)
......@@ -935,7 +938,7 @@ void calibration::Terminate()
cout << "Calibration constants are (where the '*' indicates the better value)\nPMT#: First Guess Second Guess\n" << endl;
for (Int_t i=0; i<4; i++)
{
cout << Form("PMT%d:", i+1) << setw(8) << Form("%3.0f", calibration_mk1[i]) << (pmt_calib[i] < pmt_calib_mk2[i] ? "*" : " ") << setw(13) << Form("%3.0f", calibration_mk2[i]) << (pmt_calib[i] > pmt_calib_mk2[i] ? "*\n" : "\n");
cout << Form("PMT%d:", i+1) << setw(8) << Form("%3.3f", calibration_mk1[i]) << (pmt_calib[i] < pmt_calib_mk2[i] ? "*" : " ") << setw(13) << Form("%3.3f", calibration_mk2[i]) << (pmt_calib[i] > pmt_calib_mk2[i] ? "*\n" : "\n");
}
printf("\n");
......@@ -954,7 +957,7 @@ void calibration::Terminate()
calibration << "phgcer_adc_to_npe = ";
for (Int_t ipmt = 0; ipmt < 4; ipmt++)
{
calibration << Form("1./%3.0f. ", (pmt_calib[ipmt] < pmt_calib_mk2[ipmt]) ? calibration_mk1[ipmt] : calibration_mk2[ipmt]);
calibration << Form("1./%3.3f. ", (pmt_calib[ipmt] < pmt_calib_mk2[ipmt]) ? calibration_mk1[ipmt] : calibration_mk2[ipmt]);
}
calibration.close();
......
CALIBRATION/shms_hgcer_calib/efficiencies.C
+ 7
− 7
View file @ 082b78d7
......@@ -171,20 +171,20 @@ void efficiencies::SlaveBegin(TTree * /*tree*/)
fBeta_Full = new TH1F("Beta_Full", "Full beta for events;Beta;Counts", 1000, -5, 5);
GetOutputList()->Add(fBeta_Full);
fTiming_Cut = new TH1F("Timing_Cut", "Timing cut used for 'good' hits;Time (ns);Counts", 10000, 0, 5000);
fTiming_Cut = new TH1F("Timing_Cut", "Timing cut used for 'good' hits;Time (ns);Counts", 10000, 1, 200);
GetOutputList()->Add(fTiming_Cut);
fTiming_Full = new TH1F("Timing_Full", "Full timing information for events;Time (ns);Counts", 10000, 1, 5000);
fTiming_Full = new TH1F("Timing_Full", "Full timing information for events;Time (ns);Counts", 10000, 1, 200);
GetOutputList()->Add(fTiming_Full);
//Histograms examining particle ID cuts
fFly_Pr_Full = new TH2F("Fly_Pr_Full", "Particle ID from calorimeter & preshower;Calorimeter.Fly.Earray;Calorimeter.Pr.Eplane", 200, 0.0, 1.0, 200, 0.0, 1.0);
fFly_Pr_Full = new TH2F("Fly_Pr_Full", "Particle ID from calorimeter & preshower;Calorimeter (GeV);Pre-Shower (GeV)", 200, 0.0, 1.0, 200, 0.0, 1.0);
GetOutputList()->Add(fFly_Pr_Full);
fFly_Pr_eCut = new TH2F("Fly_Pr_eCut", "calorimeter & preshower electrons;Calorimeter.Fly.Earray;Calorimeter.Pr.Eplane", 200, 0.0, 1.0, 200, 0.0, 1.0);
fFly_Pr_eCut = new TH2F("Fly_Pr_eCut", "calorimeter & preshower electrons;Calorimeter (GeV);Pre-Shower (GeV)", 200, 0.0, 1.0, 200, 0.0, 1.0);
GetOutputList()->Add(fFly_Pr_eCut);
fFly_Pr_piCut = new TH2F("Fly_Pr_piCut", "calorimeter & preshower pions;Calorimeter.Fly.Earray;Calorimeter.Pr.Eplane", 200, 0.0, 1.0, 200, 0.0, 1.0);
fFly_Pr_piCut = new TH2F("Fly_Pr_piCut", "calorimeter & preshower pions;Calorimeter (GeV);Pre-Shower (GeV)", 200, 0.0, 1.0, 200, 0.0, 1.0);
GetOutputList()->Add(fFly_Pr_piCut);
printf("\n\n");
......@@ -234,8 +234,8 @@ Bool_t efficiencies::Process(Long64_t entry)
//Require the signal passes a timing cut
fNGC ? b_P_ngcer_goodAdcPulseTime->GetEntry(entry) : b_P_hgcer_goodAdcPulseTime->GetEntry(entry);
fTiming_Full->Fill(fNGC ? P_ngcer_goodAdcPulseTime[ipmt] : P_hgcer_goodAdcPulseTime[ipmt]);
if (fNGC ? P_ngcer_goodAdcPulseTime[ipmt] < 1000 || P_ngcer_goodAdcPulseTime[ipmt] > 2000 :
P_hgcer_goodAdcPulseTime[ipmt] < 1000 || P_hgcer_goodAdcPulseTime[ipmt] > 2000) continue;
if (fNGC ? P_ngcer_goodAdcPulseTime[ipmt] < 50 || P_ngcer_goodAdcPulseTime[ipmt] > 125 :
P_hgcer_goodAdcPulseTime[ipmt] < 70 || P_hgcer_goodAdcPulseTime[ipmt] > 135) continue;
fTiming_Cut->Fill(fNGC ? P_ngcer_goodAdcPulseTime[ipmt] : P_hgcer_goodAdcPulseTime[ipmt]);
//Require the signal passes a tracking cut, with a threshold NPE cut as well
......
CALIBRATION/shms_hgcer_calib/run_cal.C
+ 1
− 1
View file @ 082b78d7
......@@ -31,7 +31,7 @@ void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0)
TString eff_option = eff_raw;
TChain ch("T");
ch.Add(Form("../../ROOTfiles/shms_replay_%d_%d.root", RunNumber, NumEvents));
ch.Add(Form("../../ROOTfiles/shms_replay_production_%d_%d.root", RunNumber, NumEvents));
TProof *proof = TProof::Open("");
proof->SetProgressDialog(0);
ch.SetProof();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment