From 14caa016e52b78f391f0d6c50a6299d9d8c86867 Mon Sep 17 00:00:00 2001
From: Ryan Ambrose <george.ryan.ambrose@gmail.com>
Date: Wed, 4 Apr 2018 15:15:51 -0600
Subject: [PATCH] Generalized calibration to first determine appropriate timing
cuts (hence new "preprocess" files), removing the need of arbitrary
hard-coded timing values.
---
CALIBRATION/shms_hgcer_calib/calibration.C | 94 ++++++----
CALIBRATION/shms_hgcer_calib/calibration.h | 16 +-
CALIBRATION/shms_hgcer_calib/efficiencies.C | 10 +-
CALIBRATION/shms_hgcer_calib/preprocess.C | 191 ++++++++++++++++++++
CALIBRATION/shms_hgcer_calib/preprocess.h | 99 ++++++++++
CALIBRATION/shms_hgcer_calib/run_cal.C | 36 +++-
6 files changed, 401 insertions(+), 45 deletions(-)
create mode 100644 CALIBRATION/shms_hgcer_calib/preprocess.C
create mode 100644 CALIBRATION/shms_hgcer_calib/preprocess.h
diff --git a/CALIBRATION/shms_hgcer_calib/calibration.C b/CALIBRATION/shms_hgcer_calib/calibration.C
index 3c7c23e0..a627ceec 100644
--- a/CALIBRATION/shms_hgcer_calib/calibration.C
+++ b/CALIBRATION/shms_hgcer_calib/calibration.C
@@ -49,8 +49,9 @@ void calibration::Begin(TTree * /*tree*/)
printf("\n\n");
TString option = GetOption();
+ TString report_option = option(0,option.Length()-79);
Info("Begin", "Script will fail unless 'calibration.C+' is used");
- Info("Begin", "Starting calibration process with option: %s", option.Data());
+ Info("Begin", "Starting calibration process with option: %s", report_option.Data());
Info("Begin", "To load all branches, use option readall (warning, very slow)");
Info("Begin", "To see details of calibration, use option showall");
Info("Begin", "Default calibration is the HGC, for NGC use option NGC");
@@ -65,7 +66,7 @@ void calibration::Begin(TTree * /*tree*/)
if (option.Contains("showall")) fFullShow = kTRUE;
if (option.Contains("trackfired")) fTrack = kTRUE;
if (option.Contains("pions") || option.Contains("pion")) fPions = kTRUE;
- if (option.Contains("cut") || fPions || option.Contains("cuts")) fCut = kTRUE;
+ if (option.Contains("cut") || fPions || option.Contains("cuts")) fCut = kTRUE;
}
void calibration::SlaveBegin(TTree * /*tree*/)
@@ -76,7 +77,25 @@ void calibration::SlaveBegin(TTree * /*tree*/)
printf("\n\n");
TString option = GetOption();
-
+
+ TString timing_mean_1 = option(option.Length()-79,option.Length()-69);
+ TString timing_std_1 = option(option.Length()-67,option.Length()-60);
+ TString timing_mean_2 = option(option.Length()-59,option.Length()-49);
+ TString timing_std_2 = option(option.Length()-47,option.Length()-40);
+ TString timing_mean_3 = option(option.Length()-39,option.Length()-29);
+ TString timing_std_3 = option(option.Length()-27,option.Length()-20);
+ TString timing_mean_4 = option(option.Length()-19,option.Length()-9);
+ TString timing_std_4 = option(option.Length()-7,option.Length()-0);
+
+ timing_mean[0] = timing_mean_1.Atof();
+ timing_std[0] = timing_std_1.Atof();
+ timing_mean[1] = timing_mean_2.Atof();
+ timing_std[1] = timing_std_2.Atof();
+ timing_mean[2] = timing_mean_3.Atof();
+ timing_std[2] = timing_std_3.Atof();
+ timing_mean[3] = timing_mean_4.Atof();
+ timing_std[3] = timing_std_4.Atof();
+
//Check option
if (option.Contains("readall")) fFullRead = kTRUE;
if (option.Contains("NGC")) fNGC = kTRUE;
@@ -100,8 +119,8 @@ void calibration::SlaveBegin(TTree * /*tree*/)
if (fNGC) //Set up histograms for NGC
{
ADC_min = -10;
- ADC_max = 250;
- bins = 2*(abs(ADC_min) + abs(ADC_max));
+ ADC_max = 200;
+ bins = 12*(abs(ADC_min) + abs(ADC_max));
}
if (!fNGC) //Set up histograms for HGC
@@ -148,7 +167,7 @@ void calibration::SlaveBegin(TTree * /*tree*/)
GetOutputList()->Add(fCut_electron);
fCut_pion = new TH2F("Cut_pion", "Visualization of pion cut; Calorimeter Energy (GeV); Pre-Shower Energy (GeV)", 250, 0, 1.0, 250, 0, 1.0);
GetOutputList()->Add(fCut_pion);
-
+
printf("\n\n");
}
@@ -172,10 +191,17 @@ Bool_t calibration::Process(Long64_t entry)
//
// The return value is currently not used.
-
//Output to verify script is working, and store the total number of events
if (entry % 100000 == 0) printf("Processing Entry number %lld\n",entry);
+ if (entry == 1)
+ {
+ cout << timing_mean[0] << " " << timing_std[0] << endl;
+ cout << timing_mean[1] << " " << timing_std[1] << endl;
+ cout << timing_mean[2] << " " << timing_std[2] << endl;
+ cout << timing_mean[3] << " " << timing_std[3] << endl;
+ }
+
//Define quantities to loop over
Int_t fpmts;
fpmts = fNGC ? fngc_pmts : fhgc_pmts; //Note HGC & NGC have the same # of PMTS
@@ -183,7 +209,7 @@ Bool_t calibration::Process(Long64_t entry)
//Get the entry to loop over
if (fFullRead) fChain->GetTree()->GetEntry(entry);
else b_Ndata_P_tr_p->GetEntry(entry);
-
+
//Require only one good track reconstruction for the event
if (Ndata_P_tr_p != 1) return kTRUE;
@@ -200,11 +226,10 @@ Bool_t calibration::Process(Long64_t entry)
for (Int_t ipmt = 0; ipmt < fpmts; ipmt++)
{
//Perform a loose timing cut
- if (!fFullRead) fNGC ? b_P_ngcer_goodAdcPulseTime->GetEntry(entry) : b_P_hgcer_goodAdcTdcDiffTime->GetEntry(entry);
- fTiming_Full->Fill(fNGC ? P_ngcer_goodAdcPulseTime[ipmt] : P_hgcer_goodAdcTdcDiffTime[ipmt]);
- if (fNGC ? P_ngcer_goodAdcPulseTime[ipmt] < 50 || P_ngcer_goodAdcPulseTime[ipmt] > 125 :
- P_hgcer_goodAdcTdcDiffTime[ipmt] > -17.0 || P_hgcer_goodAdcTdcDiffTime[ipmt] < -30.0) continue;
- fTiming_Cut->Fill(fNGC ? P_ngcer_goodAdcPulseTime[ipmt] : P_hgcer_goodAdcTdcDiffTime[ipmt]);
+ if (!fFullRead) fNGC ? b_P_ngcer_goodAdcTdcDiffTime->GetEntry(entry) : b_P_hgcer_goodAdcTdcDiffTime->GetEntry(entry);
+ fTiming_Full->Fill(fNGC ? P_ngcer_goodAdcTdcDiffTime[ipmt] : P_hgcer_goodAdcTdcDiffTime[ipmt]);
+ if (fNGC ? P_ngcer_goodAdcTdcDiffTime[ipmt] > -10.0 || P_ngcer_goodAdcTdcDiffTime[ipmt] < -35.0 : TMath::Abs(P_hgcer_goodAdcTdcDiffTime[ipmt] - timing_mean[ipmt]) > 3*timing_std[ipmt]) continue;
+ fTiming_Cut->Fill(fNGC ? P_ngcer_goodAdcTdcDiffTime[ipmt] : P_hgcer_goodAdcTdcDiffTime[ipmt]);
//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);
@@ -491,17 +516,17 @@ void calibration::Terminate()
}
//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++)
- {
- fNGC ? PulseInt_quad[iquad][ipmt]->Rebin(20) : PulseInt_quad[iquad][ipmt]->Rebin(20);
- }
- fNGC ? PulseInt[ipmt]->Rebin(20) : PulseInt[ipmt]->Rebin(20);
- }
-*/
+ if (fTrack) {
+ for (Int_t ipmt=0; ipmt < (fNGC ? fngc_pmts : fhgc_pmts); ipmt++)
+ {
+ for (Int_t iquad=0; iquad<4; iquad++)
+ {
+ fNGC ? PulseInt_quad[iquad][ipmt]->Rebin(4) : PulseInt_quad[iquad][ipmt]->Rebin(4);
+ }
+ fNGC ? PulseInt[ipmt]->Rebin(4) : PulseInt[ipmt]->Rebin(4);
+ }
+ }
+
//Canvases to display cut information
if (fFullShow)
{
@@ -598,7 +623,7 @@ 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.0, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 2.5, "nobackground&&nodraw", 0.001);
+ fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 2.5, "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();
@@ -811,25 +836,24 @@ 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(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);
+ fNGC ? PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(0,30) : PulseInt_quad[ipmt][ipmt]->GetXaxis()->SetRangeUser(0,30);
+ fFullShow ? s->Search(PulseInt_quad[iquad][ipmt], 1.0, "nobackground", 0.001) : s->Search(PulseInt_quad[iquad][ipmt], 1.5, "nobackground&&nodraw", 0.001);
TList *functions = PulseInt_quad[iquad][ipmt]->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
Double_t *xpeaks = pm->GetX();
+ PulseInt_quad[iquad][ipmt]->GetXaxis()->SetRangeUser(-1,200);
//Use the peak to fit the SPE with a Gaussian to determine the mean
- Gauss1->SetRange(xpeaks[1]-3, xpeaks[1]+3);
- Gauss1->SetParameter(1, xpeaks[1]);
+ Gauss1->SetRange(xpeaks[0]-3, xpeaks[0]+3);
+ Gauss1->SetParameter(1, xpeaks[0]);
Gauss1->SetParameter(2, 10.);
Gauss1->SetParLimits(0, 0., 2000.);
- Gauss1->SetParLimits(1, xpeaks[1]-3, xpeaks[1]+3);
+ Gauss1->SetParLimits(1, xpeaks[0]-3, xpeaks[0]+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
- //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);
+ if (xpeaks[0] != 0.0 && PulseInt_quad[iquad][ipmt]->GetBinContent(PulseInt_quad[iquad][ipmt]->GetXaxis()->FindBin(xpeaks[0])) > 10 && ipmt != iquad) mean[iquad] = Gauss1->GetParameter(1);
}
Double_t xscale = 0.0;
@@ -866,7 +890,7 @@ void calibration::Terminate()
if (fFullShow) final_spectra_ipmt->cd(1);
//Find the location of the SPE and subtract from 1.0 to determine accuracy of calibration
- Gauss1->SetRange(0.50, 2.0);
+ Gauss1->SetRange(0.50, 1.50);
Gauss1->SetParameter(0, 0.05);
Gauss1->SetParameter(1, 1.0);
Gauss1->SetParameter(2, 0.3);
@@ -898,7 +922,7 @@ void calibration::Terminate()
if (fFullShow) final_spectra_mk2_ipmt->cd(1);
//Find the location of the SPE and subtract from 1.0 to determine accuracy of calibration
- Gauss1->SetRange(0.50,2.0);
+ Gauss1->SetRange(0.50, 1.50);
Gauss1->SetParameter(0, 0.05);
Gauss1->SetParameter(1, 1.0);
Gauss1->SetParameter(2, 0.3);
diff --git a/CALIBRATION/shms_hgcer_calib/calibration.h b/CALIBRATION/shms_hgcer_calib/calibration.h
index f73c1ef2..2ae4093c 100644
--- a/CALIBRATION/shms_hgcer_calib/calibration.h
+++ b/CALIBRATION/shms_hgcer_calib/calibration.h
@@ -14,6 +14,7 @@
#include <TSelector.h>
#include <TH1.h>
#include <TH2.h>
+#include <TNtuple.h>
const Int_t fhgc_pmts = 4;
const Int_t fngc_pmts = 4;
@@ -67,8 +68,11 @@ public :
TCanvas *final_spectra_combined;
TCanvas *final_spectra_combined_mk2;
TCanvas *scaled_poisson;
- TCanvas *scaled_total;
-
+ TCanvas *scaled_total;
+
+ // Declaration of preprocessing quantities
+ Double_t timing_mean[4];
+ Double_t timing_std[4];
// Declaration of leaf types
Int_t Ndata_P_aero_goodNegAdcPed;
@@ -623,6 +627,8 @@ public :
Double_t P_ngcer_goodAdcPulseIntRaw[4]; //[Ndata.P.ngcer.goodAdcPulseIntRaw]
Int_t Ndata_P_ngcer_goodAdcPulseTime;
Double_t P_ngcer_goodAdcPulseTime[4]; //[Ndata.P.ngcer.goodAdcPulseTime]
+ Int_t Ndata_P_ngcer_goodAdcTdcDiffTime;
+ Double_t P_ngcer_goodAdcTdcDiffTime[4];
Int_t Ndata_P_ngcer_npe;
Double_t P_ngcer_npe[4]; //[Ndata.P.ngcer.npe]
Int_t Ndata_P_ngcer_numGoodAdcHits;
@@ -1452,6 +1458,8 @@ public :
TBranch *b_P_ngcer_goodAdcPulseIntRaw; //!
TBranch *b_Ndata_P_ngcer_goodAdcPulseTime; //!
TBranch *b_P_ngcer_goodAdcPulseTime; //!
+ TBranch *b_Ndata_P_ngcer_goodAdcTdcDiffTime;
+ TBranch *b_P_ngcer_goodAdcTdcDiffTime;
TBranch *b_Ndata_P_ngcer_npe; //!
TBranch *b_P_ngcer_npe; //!
TBranch *b_Ndata_P_ngcer_numGoodAdcHits; //!
@@ -1727,7 +1735,7 @@ public :
TBranch *b_Event_Branch_fEvtHdr_fTargetPol; //!
TBranch *b_Event_Branch_fEvtHdr_fRun; //!
- calibration(TTree * /*tree*/ =0) : fChain(0) {fPulseInt = 0, fPulseInt_quad = 0, fCut_everything = 0, fCut_electron = 0, fCut_pion = 0, fBeta_Cut = 0, fBeta_Full = 0, fTiming_Cut = 0, fTiming_Full = 0, fFullRead = kFALSE, fFullShow = kFALSE, fNGC = kFALSE, fTrack = kFALSE, fCut = kFALSE, fPions = kFALSE;}
+ calibration(TTree * /*tree*/ =0) : fChain(0) {fPulseInt = 0, fPulseInt_quad = 0, fCut_everything = 0, fCut_electron = 0, fCut_pion = 0, fBeta_Cut = 0, fBeta_Full = 0, fTiming_Cut = 0, fTiming_Full = 0, fFullRead = kFALSE, fFullShow = kFALSE, fNGC = kFALSE, fTrack = kFALSE, fCut = kFALSE, fPions = kFALSE, timing_mean[0] = 0, timing_mean[1] = 0, timing_mean[2] = 0, timing_mean[3] = 0, timing_std[0] = 0, timing_std[1] = 0, timing_std[2] = 0, timing_std[3] = 0;}
virtual ~calibration() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
@@ -2316,6 +2324,8 @@ void calibration::Init(TTree *tree)
fChain->SetBranchAddress("P.ngcer.goodAdcPulseIntRaw", P_ngcer_goodAdcPulseIntRaw, &b_P_ngcer_goodAdcPulseIntRaw);
fChain->SetBranchAddress("Ndata.P.ngcer.goodAdcPulseTime", &Ndata_P_ngcer_goodAdcPulseTime, &b_Ndata_P_ngcer_goodAdcPulseTime);
fChain->SetBranchAddress("P.ngcer.goodAdcPulseTime", P_ngcer_goodAdcPulseTime, &b_P_ngcer_goodAdcPulseTime);
+ fChain->SetBranchAddress("Ndata.P.ngcer.goodAdcTdcDiffTime", &Ndata_P_ngcer_goodAdcTdcDiffTime, &b_Ndata_P_ngcer_goodAdcTdcDiffTime);
+ fChain->SetBranchAddress("P.ngcer.goodAdcTdcDiffTime", P_ngcer_goodAdcTdcDiffTime, &b_P_ngcer_goodAdcTdcDiffTime);
fChain->SetBranchAddress("Ndata.P.ngcer.npe", &Ndata_P_ngcer_npe, &b_Ndata_P_ngcer_npe);
fChain->SetBranchAddress("P.ngcer.npe", P_ngcer_npe, &b_P_ngcer_npe);
fChain->SetBranchAddress("Ndata.P.ngcer.numGoodAdcHits", &Ndata_P_ngcer_numGoodAdcHits, &b_Ndata_P_ngcer_numGoodAdcHits);
diff --git a/CALIBRATION/shms_hgcer_calib/efficiencies.C b/CALIBRATION/shms_hgcer_calib/efficiencies.C
index 3d299a0e..543909e8 100644
--- a/CALIBRATION/shms_hgcer_calib/efficiencies.C
+++ b/CALIBRATION/shms_hgcer_calib/efficiencies.C
@@ -262,8 +262,8 @@ Bool_t efficiencies::Process(Long64_t entry)
Float_t esemimajor_axis = 0.28;
Float_t esemiminor_axis = 0.04;*/
Float_t eangle = 3.0*3.14159/4.0;
- Float_t ex_center = 0.375;
- Float_t ey_center = 0.360;
+ Float_t ex_center = 0.490;
+ Float_t ey_center = 0.400;
Float_t esemimajor_axis = 0.38;
Float_t esemiminor_axis = 0.05;
if (pow((P_cal_fly_earray/p - ex_center)*cos(eangle) + (P_cal_pr_eplane/p - ey_center)*sin(eangle),2)/pow(esemimajor_axis,2) +
@@ -290,9 +290,9 @@ Bool_t efficiencies::Process(Long64_t entry)
Float_t piy_center = 0.03;
Float_t pisemimajor_axis = 0.1;
Float_t pisemiminor_axis = 0.02;*/
- Float_t pix_center = 0.00;
- Float_t piy_center = 0.00;
- Float_t pisemimajor_axis = 0.01;
+ Float_t pix_center = 0.30;
+ Float_t piy_center = 0.05;
+ Float_t pisemimajor_axis = 0.08;
Float_t pisemiminor_axis = 0.02;
if (pow((P_cal_fly_earray/p - pix_center)*cos(piangle) + (P_cal_pr_eplane/p - piy_center)*sin(piangle),2)/pow(pisemimajor_axis,2) +
pow((P_cal_fly_earray/p - pix_center)*sin(piangle) - (P_cal_pr_eplane/p - piy_center)*cos(piangle),2)/pow(pisemiminor_axis,2) < 1)
diff --git a/CALIBRATION/shms_hgcer_calib/preprocess.C b/CALIBRATION/shms_hgcer_calib/preprocess.C
new file mode 100644
index 00000000..a6bbff4d
--- /dev/null
+++ b/CALIBRATION/shms_hgcer_calib/preprocess.C
@@ -0,0 +1,191 @@
+#define preprocess_cxx
+// The class definition in preprocess.h has been generated automatically
+// by the ROOT utility TTree::MakeSelector(). This class is derived
+// from the ROOT class TSelector. For more information on the TSelector
+// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
+
+
+// The following methods are defined in this file:
+// Begin(): called every time a loop on the tree starts,
+// a convenient place to create your histograms.
+// SlaveBegin(): called after Begin(), when on PROOF called only on the
+// slave servers.
+// Process(): called for each event, in this function you decide what
+// to read and fill your histograms.
+// SlaveTerminate: called at the end of the loop on the tree, when on PROOF
+// called only on the slave servers.
+// Terminate(): called at the end of the loop on the tree,
+// a convenient place to draw/fit your histograms.
+//
+// To use this file, try the following session on your Tree T:
+//
+// root> T->Process("preprocess.C")
+// root> T->Process("preprocess.C","some options")
+// root> T->Process("preprocess.C+")
+//
+
+
+#include "preprocess.h"
+#include <TH1.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TSpectrum.h>
+#include <TPolyMarker.h>
+#include <TStyle.h>
+#include <TNtuple.h>
+
+void preprocess::Begin(TTree * /*tree*/)
+{
+ // The Begin() function is called at the start of the query.
+ // When running with PROOF Begin() is only called on the client.
+ // The tree argument is deprecated (on PROOF 0 is passed).
+
+ printf("\n\n");
+
+ Info("Begin", "Starting preprocessing script to: determine correct fADC timing cuts");
+
+ printf("\n");
+ TString option = GetOption();
+ if (option.Contains("NGC")) fNGC = kTRUE;
+ if (option.Contains("showall")) fFullShow = kTRUE;
+}
+
+void preprocess::SlaveBegin(TTree * /*tree*/)
+{
+ // The SlaveBegin() function is called after the Begin() function.
+ // When running with PROOF SlaveBegin() is called on each slave server.
+ // The tree argument is deprecated (on PROOF 0 is passed).
+
+ TString option = GetOption();
+
+ fTiming = new TH1F*[4];
+ for (Int_t ipmt = 0; ipmt < 4; ipmt++)
+ {
+ fTiming[ipmt] = new TH1F(Form("Timing_PMT%d",ipmt+1),Form("ADC TDC Diff Time PMT%d; Time (ns); Counts",ipmt+1), 400, -50, 0);
+ GetOutputList()->Add(fTiming[ipmt]);
+ }
+
+}
+
+Bool_t preprocess::Process(Long64_t entry)
+{
+ // The Process() function is called for each entry in the tree (or possibly
+ // keyed object in the case of PROOF) to be processed. The entry argument
+ // specifies which entry in the currently loaded tree is to be processed.
+ // When processing keyed objects with PROOF, the object is already loaded
+ // and is available via the fObject pointer.
+ //
+ // This function should contain the \"body\" of the analysis. It can contain
+ // simple or elaborate selection criteria, run algorithms on the data
+ // of the event and typically fill histograms.
+ //
+ // The processing can be stopped by calling Abort().
+ //
+ // Use fStatus to set the return value of TTree::Process().
+ //
+ // The return value is currently not used.
+
+ fReader.SetEntry(entry);
+
+ //Output to verify script is working, and store the total number of events
+ if (entry % 100000 == 0) printf("Processing Entry number %lld\n",entry);
+
+ //Require only one good track reconstruction for the event
+ if (Ndata_P_tr_p[0] != 1) return kTRUE;
+
+ //Redundant, but useful if multiple tracks will be allowed
+ for (Int_t itrack = 0; itrack < Ndata_P_tr_p[0]; itrack++)
+ {
+ //Require cut on particle velocity
+ if (TMath::Abs(P_tr_beta[itrack] -1.0) > 0.2) return kTRUE;
+
+ //Fill the histograms
+ for (Int_t ipmt = 0; ipmt < 4; ipmt++)
+ {
+ fTiming[ipmt]->Fill(P_hgcer_goodAdcTdcDiffTime[ipmt]);
+ }
+ }
+
+ return kTRUE;
+}
+
+void preprocess::SlaveTerminate()
+{
+ // The SlaveTerminate() function is called after all entries or objects
+ // have been processed. When running with PROOF SlaveTerminate() is called
+ // on each slave server.
+
+}
+
+void preprocess::Terminate()
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+ printf("\n\n");
+ Info("Terminate", "calibrating '%s'", (fNGC ? "NGC" : "HGC"));
+ Info("Terminate", "'%s' showing", (fFullShow ? "full" : "minimal"));
+ Info("Terminate","Histograms formed, now determining timing peaks");
+
+ //Create temporary file to store fit information
+ TFile ofile("Timing_Cuts.root","RECREATE");
+
+ //Create arrays to store Gaussian fit information
+ Float_t mean[4];
+ Float_t sigma[4];
+
+ //Create Ntuple to store timing data
+ TNtuple timing_data("timing_data","Storage of Timing Information","Mean:Std");
+
+ //Single Gaussian to find mean and spread of peak
+ TF1 *Gauss1 = new TF1("Gauss1",pregauss,-10,200,3);
+ Gauss1->SetParNames("Amplitude","Mean","Std. Dev.");
+
+ //Have to extract histograms from OutputList
+ TH1F *Timing[4];
+ for (Int_t ipmt = 0; ipmt < 4; ipmt++)
+ {
+ Timing[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("Timing_PMT%d",ipmt+1)));
+ }
+
+ //Visualize the timing information
+ if (fFullShow)
+ {
+ cTiming = new TCanvas("Timing", "Timing Information");
+ cTiming->Divide(2,2);
+ }
+
+ TSpectrum *s = new TSpectrum(2);
+
+ //Main loop to find timing peaks
+ for (Int_t ipmt = 0; ipmt < 4; ipmt++)
+ {
+ //Use peak finding algorithm from TSpectrum
+
+ if (fFullShow) cTiming->cd(ipmt+1);
+
+ //Perform search for largest peak and save into array xpeaks
+ fFullShow ? s->Search(Timing[ipmt],2.0,"nobackground",0.9) : s->Search(Timing[ipmt],2.0,"nobackground&&nodraw",0.9);
+ TList *functions = Timing[ipmt]->GetListOfFunctions();
+ TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
+ Double_t *xpeaks = pm->GetX();
+
+ //Use the peak to fit with a Gaussian to determine range of cut
+ Gauss1->SetRange(xpeaks[0]-3, xpeaks[0]+3);
+ Gauss1->SetParameter(0, 10000.);
+ Gauss1->SetParameter(1, xpeaks[0]);
+ Gauss1->SetParameter(2, 0.4);
+ //Gauss1->SetParLimits(0, 0., 10000.);
+ Gauss1->SetParLimits(1, xpeaks[0]-3, xpeaks[0]+3);
+ Gauss1->SetParLimits(2, 0.2, 1.);
+ fFullShow ? Timing[ipmt]->Fit("Gauss1","RQ") : Timing[ipmt]->Fit("Gauss1","RQN");
+ mean[ipmt] = Gauss1->GetParameter(1);
+ sigma[ipmt] = Gauss1->GetParameter(2);
+ timing_data.Fill(mean[ipmt],sigma[ipmt]);
+ }
+
+ cout << "\nTiming cut results are (mean & sigma):\n" << Form("PMT 1: %f, %f\nPMT 2: %f, %f\nPMT 3: %f, %f\nPMT 4: %f, %f\n\n",mean[0],sigma[0],mean[1],sigma[1],mean[2],sigma[2],mean[3],sigma[3]);
+
+ timing_data.Write();
+ ofile.Close();
+}
diff --git a/CALIBRATION/shms_hgcer_calib/preprocess.h b/CALIBRATION/shms_hgcer_calib/preprocess.h
new file mode 100644
index 00000000..db4376a0
--- /dev/null
+++ b/CALIBRATION/shms_hgcer_calib/preprocess.h
@@ -0,0 +1,99 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar 2 15:04:26 2018 by ROOT version 6.10/02
+// from TTree T/Hall A Analyzer Output DST
+// found on file: /home/rambrose/hallc_replay/ROOTfiles/shms_replay_production_all_2732_-1.root
+//////////////////////////////////////////////////////////
+
+#ifndef preprocess_h
+#define preprocess_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+#include <TTreeReader.h>
+#include <TTreeReaderValue.h>
+#include <TTreeReaderArray.h>
+
+// Headers needed by this particular selector
+
+
+class preprocess : public TSelector {
+public :
+ TTreeReader fReader; //!the tree reader
+ TTree *fChain = 0; //!pointer to the analyzed TTree or TChain
+
+ Bool_t fNGC;
+ Bool_t fFullShow;
+
+ // Declaration of Histograms
+ TH1F **fTiming;
+
+ // Declaration of Canvases
+ TCanvas *cTiming;
+
+ // Readers to access the data (delete the ones you do not need).
+ TTreeReaderArray<Double_t> P_hgcer_goodAdcTdcDiffTime = {fReader, "P.hgcer.goodAdcTdcDiffTime"};
+ TTreeReaderArray<Int_t> Ndata_P_tr_p = {fReader, "Ndata.P.tr.p"};
+ TTreeReaderArray<Double_t> P_tr_beta = {fReader, "P.tr.beta"};
+
+
+ preprocess(TTree * /*tree*/ =0) {fNGC = kFALSE, fFullShow = kFALSE; }
+ virtual ~preprocess() { }
+ virtual Int_t Version() const { return 2; }
+ virtual void Begin(TTree *tree);
+ virtual void SlaveBegin(TTree *tree);
+ virtual void Init(TTree *tree);
+ virtual Bool_t Notify();
+ virtual Bool_t Process(Long64_t entry);
+ virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+ virtual void SetOption(const char *option) { fOption = option; }
+ virtual void SetObject(TObject *obj) { fObject = obj; }
+ virtual void SetInputList(TList *input) { fInput = input; }
+ virtual TList *GetOutputList() const { return fOutput; }
+ virtual void SlaveTerminate();
+ virtual void Terminate();
+
+ ClassDef(preprocess,0);
+
+};
+
+#endif
+
+#ifdef preprocess_cxx
+void preprocess::Init(TTree *tree)
+{
+ // The Init() function is called when the selector needs to initialize
+ // a new tree or chain. Typically here the reader is initialized.
+ // It is normally not necessary to make changes to the generated
+ // code, but the routine can be extended by the user if needed.
+ // Init() will be called many times when running on PROOF
+ // (once per file to be processed).
+
+ fReader.SetTree(tree);
+}
+
+Bool_t preprocess::Notify()
+{
+ // The Notify() function is called when a new file is opened. This
+ // can be either for a new TTree in a TChain or when when a new TTree
+ // is started when using PROOF. It is normally not necessary to make changes
+ // to the generated code, but the routine can be extended by the
+ // user if needed. The return value is currently not used.
+
+ return kTRUE;
+}
+
+//Gaussian distribution is used to find the mean of the SPE and determine spacing between subsequent NPE
+Double_t pregauss(Double_t *x, Double_t *par)
+{
+ Double_t result1 = par[0]*exp((-0.5)*(pow((x[0] - par[1]),2)/pow((par[2]),2)));
+ Double_t result2 = par[3]*exp((-0.5)*(pow((x[0] - par[4]),2)/pow((par[5]),2)));
+ Double_t result3 = par[6]*exp((-0.5)*(pow((x[0] - par[7]),2)/pow((par[8]),2)));
+ Double_t result4 = par[9]*exp((-0.5)*(pow((x[0] - par[10]),2)/pow((par[11]),2)));
+ Double_t result5 = par[12]*exp((-0.5)*(pow((x[0] - par[13]),2)/pow((par[14]),2)));
+ return result1 + result2 + result3 + result4 + result5;
+}
+
+#endif // #ifdef preprocess_cxx
diff --git a/CALIBRATION/shms_hgcer_calib/run_cal.C b/CALIBRATION/shms_hgcer_calib/run_cal.C
index e3d32b47..868a8ca6 100644
--- a/CALIBRATION/shms_hgcer_calib/run_cal.C
+++ b/CALIBRATION/shms_hgcer_calib/run_cal.C
@@ -4,7 +4,7 @@
#include <string>
#include <stdio.h>
-void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0)
+void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0, Int_t coin = 0)
{
if (RunNumber == 0)
{
@@ -17,6 +17,11 @@ void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0)
cout << "\nNumber of Events to analyze: ";
cin >> NumEvents;
}
+ if (coin == 0)
+ {
+ cout << "\nIf this is a coincident run enter 1: ";
+ cin >> coin;
+ }
cin.ignore(numeric_limits<streamsize>::max(), '\n');
@@ -30,14 +35,41 @@ void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0)
getline(std::cin, eff_raw);
TString eff_option = eff_raw;
+ cout << "\n\n";
+
TChain ch("T");
- ch.Add(Form("../../ROOTfiles/shms_replay_production_%d_%d.root", RunNumber, NumEvents));
+ if (coin == 1) ch.Add(Form("../../ROOTfiles/coin_replay_production_%d_%d.root", RunNumber, NumEvents));
+ else ch.Add(Form("../../ROOTfiles/shms_replay_production_all_%d_%d.root", RunNumber, NumEvents));
TProof *proof = TProof::Open("");
proof->SetProgressDialog(0);
ch.SetProof();
if (calib_option != "NA")
{
+ //Perform a "preprocess" of replay to extact timing information
+ ch.Process("preprocess.C+",calib_option);
+ if (calib_option.Contains("showall"))
+ {
+ TString preprocessing_input;
+ cout << "\n\nPlease verify if timing cuts are correct (y/n): ";
+ cin >> preprocessing_input;
+ if (preprocessing_input != "y") return;
+ }
+
+ //Obtain cut information from preprocessing
+ TFile preprocess_file("Timing_Cuts.root");
+ TNtuple *input_data = new TNtuple("input_data","Storage for Client Timing Information","Mean:Std");
+ preprocess_file.GetObject("timing_data",input_data);
+ float *timing_row_content;
+ for (Int_t irow = 0; irow < input_data->GetEntries(); irow++)
+ {
+ input_data->GetEntry(irow);
+ timing_row_content = input_data->GetArgs();
+ calib_option.Append(Form(" %f %f",timing_row_content[0],timing_row_content[1]));
+ }
+ preprocess_file.Close();
+
+ //Start calibration process
ch.Process("calibration.C+",calib_option);
cout << "\n\nUpdate calibration constants with the better estimate (y/n)? ";
--
GitLab