CALIBRATION/bpm_calib/Harp_files/harp_info.txt.extra.no1766.with3129.short

+1768  0.956 -1.194  1.033 -1.220 345.528 154.612 -1.00 -0.98 -0.80 -1.32 -2.0 -2.0 0.1
+1767 -0.620 -1.153 -2.499 -1.261 345.528 154.612  0.26 -1.00  1.46 -1.42  2.0  2.0 0.1
+3129  1.010 -0.420  0.400 -0.160 345.528 154.612 -1.00 -0.99 -0.40 -0.96 -1.0 -1.0 0.1 

CALIBRATION/bpm_calib/Harp_files/harp_info.txt.just1765and1766

+1766 -0.557 2.319 -2.482 2.077 345.528 154.612 0.15 1.10 1.47 0.97 2.0 2.0 0.1
+1765  0.927 1.035 0.992 1.998 345.528 154.612 -1.01 1.17 -0.76 1.05 -2.0 2.0 0.1

CALIBRATION/bpm_calib/Harp_files/harp_info.txt.no1765

+1771 0.194 -0.605 -1.548 0.328 345.528 154.612 -0.56 -0.62 0.67 -0.47 1.0 0.0 0.1 
+1770 0.838 -0.620 0.155 0.2969 345.528 154.612 -0.98 -0.61 -0.35 -0.43 -1.0 0.0 0.1 
+1769 0.683 -1.079 -0.677 -0.512 345.528 154.612 -0.99 -0.99 0.05 -1.00 0.0 -1.0 0.1
+1768 0.956 -1.194 1.033 -1.220 345.528 154.612 -1.00 -0.98 -0.80 -1.32 -2.0 -2.0 0.1
+1767  -0.620 -1.153 -2.499 -1.261 345.528 154.612 0.26 -1.0 1.46 -1.42 2.0 2.0 0.1
+1766 -0.557 2.319 -2.482 2.077 345.528 154.612 0.15 1.10 1.47 0.97 2.0 2.0 0.1
+1763 0.691 -0.645 -0.630 0.201 345.528 154.612 -1.00 -0.61 0.02 -0.47 0.0 0.0 0.1
+1762 0.825 -1.063 0.221 -0.517 345.528 154.612 -1.00 -0.99 -0.40 -0.96 -1.0 -1.0 0.1 

CALIBRATION/bpm_calib/Harp_files/harp_info.txt.no1765or1766

+1771 0.194 -0.605 -1.548 0.328 345.528 154.612 -0.56 -0.62 0.67 -0.47 1.0 0.0 0.1 
+1770 0.838 -0.620 0.155 0.2969 345.528 154.612 -0.98 -0.61 -0.35 -0.43 -1.0 0.0 0.1 
+1769 0.683 -1.079 -0.677 -0.512 345.528 154.612 -0.99 -0.99 0.05 -1.00 0.0 -1.0 0.1
+1768 0.956 -1.194 1.033 -1.220 345.528 154.612 -1.00 -0.98 -0.80 -1.32 -2.0 -2.0 0.1
+1767  -0.620 -1.153 -2.499 -1.261 345.528 154.612 0.26 -1.0 1.46 -1.42 2.0 2.0 0.1
+1763 0.691 -0.645 -0.630 0.201 345.528 154.612 -1.00 -0.61 0.02 -0.47 0.0 0.0 0.1
+1762 0.825 -1.063 0.221 -0.517 345.528 154.612 -1.00 -0.99 -0.40 -0.96 -1.0 -1.0 0.1 

CALIBRATION/bpm_calib/Harp_files/harp_info.txt.no1766

+1771  0.194  -0.605 -1.548  0.328 345.528 154.612 -0.56 -0.62  0.67 -0.47  1.0  0.0 0.10 
+1770  0.838  -0.620  0.155  0.297 345.528 154.612 -0.98 -0.61 -0.35 -0.43 -1.0  0.0 0.10 
+1769  0.683  -1.079 -0.677 -0.512 345.528 154.612 -0.99 -0.99  0.05 -1.00  0.0 -1.0 0.10
+1768  0.956  -1.194  1.033 -1.220 345.528 154.612 -1.00 -0.98 -0.80 -1.32 -2.0 -2.0 0.10
+1767 -0.620  -1.153 -2.499 -1.261 345.528 154.612  0.26 -1.00  1.46 -1.42  2.0  2.0 0.10
+1765  0.927   1.035  0.992  1.998 345.528 154.612 -1.01  1.17 -0.76  1.05 -2.0  2.0 0.10
+1763  0.691  -0.645 -0.630  0.201 345.528 154.612 -1.00 -0.61  0.02 -0.47  0.0  0.0 0.10
+1762  0.825  -1.063  0.221 -0.517 345.528 154.612 -1.00 -0.99 -0.40 -0.96 -1.0 -1.0 0.10

CALIBRATION/bpm_calib/Harp_files/harp_info.txt.orig

+1771 0.194 -0.605 -1.548 0.328 345.528 154.612 -0.56 -0.62 0.67 -0.47 1.0 0.0 0.1 
+1770 0.838 -0.620 0.155 0.2969 345.528 154.612 -0.98 -0.61 -0.35 -0.43 -1.0 0.0 0.1 
+1769 0.683 -1.079 -0.677 -0.512 345.528 154.612 -0.99 -0.99 0.05 -1.00 0.0 -1.0 0.1
+1768 0.956 -1.194 1.033 -1.220 345.528 154.612 -1.00 -0.98 -0.80 -1.32 -2.0 -2.0 0.1
+1767  -0.620 -1.153 -2.499 -1.261 345.528 154.612 0.26 -1.0 1.46 -1.42 2.0 2.0 0.1
+1766 -0.557 2.319 -2.482 2.077 345.528 154.612 0.15 1.10 1.47 0.97 2.0 2.0 0.1
+1765  0.927 1.035 0.992 1.998 345.528 154.612 -1.01 1.17 -0.76 1.05 -2.0 2.0 0.1
+1763 0.691 -0.645 -0.630 0.201 345.528 154.612 -1.00 -0.61 0.02 -0.47 0.0 0.0 0.1
+1762 0.825 -1.063 0.221 -0.517 345.528 154.612 -1.00 -0.99 -0.40 -0.96 -1.0 -1.0 0.1 

CALIBRATION/bpm_calib/README.md

+# Hall C BPM Calibration Scripts
+
+## Instructions
+1. Source the setup.sh or setup.csh as appropriate for your shell.
+
+2. It is expected that the ROOT files are already in the ~/hallc_replay/ROOTfiles directory.
+
+3. It is expected that the HARP/BPM data for the appropriate runs are in a text file called harp_info.txt, although a different filename can be specified as an argument to the ROOT macro.
+
+4.  Within ROOT, execute the bpm calibration macro as:
+
+.x bpm_calibration.C("harp_info.txt")
+
+Notes:
+
+For the Spring 2018 running, the appropriate HARP calibration file is:
+
+harp_info.txt.extra.no1766 

CALIBRATION/bpm_calib/raster_plots.C

+// this program is for bpm callibration in Hall C . 
+// Time : Fall 2017 Run.
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TF1.h"
+#include "TGraph.h"
+#include <fstream>
+#include <iostream>
+#include <TROOT.h>
+
+void
+raster_plots(Int_t run_NUM = 4251){
+  
+  
+  gStyle->SetOptStat(0);
+  gStyle->SetOptStat(0);
+  gStyle->SetTitleFontSize(0.05);
+  gStyle->SetNdivisions(505);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadTopMargin(0.10);
+  gStyle->SetPadLeftMargin(0.08);
+  gStyle->SetPadRightMargin(0.08);
+  gStyle->SetPadBottomMargin(.14);
+  gStyle->SetTitleYOffset(1.09);
+  gStyle->SetTitleXOffset(0.855);
+  gStyle->SetTitleYSize(0.03);
+  gStyle->SetTitleXSize(0.03);
+  gStyle->SetLabelFont(62,"X");
+  gStyle->SetLabelFont(62,"Y");
+  gStyle->SetTitleFont(62,"X");
+  gStyle->SetTitleFont(62,"Y");
+  gStyle->SetLabelSize(0.025,"X");
+  gStyle->SetLabelSize(0.025,"Y");
+  gStyle->SetPaperSize(23,24);
+
+  Double_t rasterAx;
+  Double_t rasterAx_ADC;
+  Double_t rasterAx_RawADC;
+  Double_t rasterBx;
+  Double_t rasterBx_ADC;
+  Double_t rasterBx_RawADC;
+  Double_t rasterAy;
+  Double_t rasterAy_ADC;
+  Double_t rasterAy_RawADC;
+  Double_t rasterBy;
+  Double_t rasterBy_ADC;
+  Double_t rasterBy_RawADC;
+  Double_t rasterA_radius;
+  Double_t rasterB_radius;
+
+  Double_t xoffA = 0.000;
+  Double_t yoffA = 0.000;
+  Double_t xoffB = 0.000;
+  Double_t yoffB = 0.000;
+  Double_t rAcut = 5.00;
+  Double_t rBcut = 5.00;
+  
+  TFile *f = new TFile(Form("/Users/brash/Dropbox/Research/analysis/hallc_replay/ROOTfiles/coin_replay_production_%d_-1.root",run_NUM),"READ"); // %d : expects integer; %f expects float 
+    
+    
+    TTree *T = (TTree*)f->Get("T");
+    Int_t totev = T->GetEntries(); 
+    //Read the branch for the BPM positions from the EPICS 
+    T->SetBranchAddress("P.rb.raster.fr_xa",&rasterAx);
+    T->SetBranchAddress("P.rb.raster.fr_xb",&rasterBx);
+    T->SetBranchAddress("P.rb.raster.fr_ya",&rasterAy);
+    T->SetBranchAddress("P.rb.raster.fr_yb",&rasterBy);
+    T->SetBranchAddress("P.rb.raster.frxa_adc",&rasterAx_ADC);
+    T->SetBranchAddress("P.rb.raster.frxb_adc",&rasterBx_ADC);
+    T->SetBranchAddress("P.rb.raster.frya_adc",&rasterAy_ADC);
+    T->SetBranchAddress("P.rb.raster.fryb_adc",&rasterBy_ADC);
+    T->SetBranchAddress("P.rb.raster.frxaRawAdc",&rasterAx_RawADC);
+    T->SetBranchAddress("P.rb.raster.frxbRawAdc",&rasterBx_RawADC);
+    T->SetBranchAddress("P.rb.raster.fryaRawAdc",&rasterAy_RawADC);
+    T->SetBranchAddress("P.rb.raster.frybRawAdc",&rasterBy_RawADC);
+    //Creating the histogram of the BPM positions and 
+    TH1F* hrasterAx =new TH1F("rasterAx","rasterAx",100,-3.00,3.00);
+    TH1F* hrasterAy =new TH1F("rasterAy","rasterAy",100,-3.00,3.00);
+    TH1F* hrasterBx =new TH1F("rasterBx","rasterBx",100,-3.00,3.00);
+    TH1F* hrasterBy =new TH1F("rasterBy","rasterBy",100,-3.00,3.00);
+    TH2F* hrasterAxAy =new TH2F("rasterAxAy","rasterAxAy",100,-3.00,3.00,100,-3.00,3.00);
+    TH2F* hrasterBxBy =new TH2F("rasterBxBy","rasterBxBy",100,-3.00,3.00,100,-3.00,3.00);
+    TH1F* hrasterA_radius =new TH1F("rasterA_radius","rasterA_radius",100,0.00,05.00);
+    TH1F* hrasterB_radius =new TH1F("rasterB_radius","rasterB_radius",100,0.00,5.55);
+    // Fill Histograms here   
+    for (Int_t iev = 0 ; iev < totev ;iev ++){
+      if (iev%1000 == 0) cout << "Event: " << iev << endl;
+      T->GetEntry(iev);
+      Double_t rA = sqrt((10.0*rasterAx-xoffA)*(10.0*rasterAx-xoffA)+(10.0*rasterAy-yoffA)*(10.0*rasterAy-yoffA));
+      Double_t rB = sqrt((10.0*rasterBx-xoffB)*(10.0*rasterBx-xoffB)+(10.0*rasterBy-yoffB)*(10.0*rasterBy-yoffB));
+      //if (ibcm1>1){
+	hrasterAx ->Fill(10.0*rasterAx);
+	hrasterBx ->Fill(10.0*rasterBx);
+	hrasterAy ->Fill(10.0*rasterAy);
+	hrasterBy ->Fill(10.0*rasterBy);
+        hrasterAxAy->Fill(10.0*rasterAx,10.0*rasterAy);
+        hrasterBxBy->Fill(10.0*rasterBx,10.0*rasterBy);
+        if (rA<rAcut) hrasterA_radius->Fill(rA);
+        if (rB<rBcut) hrasterB_radius->Fill(rB);
+      //}
+    }
+
+    TCanvas *craster = new TCanvas("craster","Raster X, Y : Hall C", 800, 900);
+    TCanvas *crasterA = new TCanvas("crasterA","RasterA X vs Y : Hall C", 800, 900);
+    TCanvas *crasterB = new TCanvas("crasterB","RasterB X vs Y : Hall C", 800, 900);
+    TCanvas *crasterRad = new TCanvas("crasterRad","Raster Radius : Hall C", 800, 900);
+    //Draw the bpm_measured vs bpm_cal
+    craster->Divide(2,2);   
+
+    craster->cd(1);
+    hrasterAx->Draw();
+    craster->cd(2);
+    hrasterAy->Draw();
+    craster->cd(3);
+    hrasterBx->Draw();
+    craster->cd(4);
+    hrasterBy->Draw();
+    crasterA->cd();
+    //crasterA->SetLogz();
+    hrasterAxAy->Draw("COLZ");
+    crasterB->cd();
+    //crasterB->SetLogz();
+    hrasterBxBy->Draw("COLZ");
+    crasterRad->cd();
+    crasterRad->Divide(1,2);
+    crasterRad->cd(1);
+    hrasterA_radius->Draw();
+    crasterRad->cd(2);
+    hrasterB_radius->Draw();
+
+}
+

CALIBRATION/bpm_calib/setup.csh

+#!/usr/bin/csh
+
+# -----------------------------------------------------------------------------
+#  Change these if this if not where hallc_replay and hcana live
+setenv hcana_dir "/home/$USER/analysis/hcana"
+setenv hallc_replay_dir "/home/$USER/analysis/hallc_replay"
+
+# -----------------------------------------------------------------------------
+#  Change if this gives you the wrong version of root, evio, etc
+#source /site/12gev_phys/production.csh 2.1
+
+# -----------------------------------------------------------------------------
+# Source setup scripts
+set curdir=`pwd`
+cd $hcana_dir
+source setup.csh
+setenv PATH "$hcana_dir/bin:$PATH"
+echo Sourced $hcana_dir/setup.csh
+
+cd $hallc_replay_dir
+source setup.csh
+echo Sourced $hallc_replay_dir/setup.csh
+
+echo cd back to $curdir
+cd $curdir
+

CALIBRATION/bpm_calib/setup.sh

+#!/usr/bin/bash
+
+# -----------------------------------------------------------------------------
+#  Change these if this if not where hallc_replay and hcana live
+export hcana_dir=/Users/brash/Dropbox/Research/analysis/hcana
+export hallc_replay_dir=/Users/brash/Dropbox/Research/analysis/hallc_replay
+
+# -----------------------------------------------------------------------------
+#  Change if this gives you the wrong version of root, evio, etc
+#source /site/12gev_phys/production.sh 2.1
+
+# -----------------------------------------------------------------------------
+# Source setup scripts
+curdir=`pwd`
+cd $hcana_dir
+source setup.sh
+export PATH=$hcana_dir/bin:$PATH
+echo Sourced $hcana_dir/setup.sh
+
+cd $hallc_replay_dir
+source setup.sh
+echo Sourced $hallc_replay_dir/setup.sh
+
+echo cd back to $curdir
+cd $curdir
+

CALIBRATION/bpm_calib/target_bpm.py

+#!/app/bin/python3
+import sys,os
+import gc
+import math
+import time
+import numpy as np
+import pylab as py
+import matplotlib.pyplot as plt
+from scipy.odr import *
+from subprocess import call
+from subprocess import check_output
+
+# Nominal beam position (note that 3H07C is offset by (+1.0, -1.0))
+NOMINAL_POS = {
+'3H07A': {'XPOS': 0.0,
+          'YPOS': 0.8},
+'3H07C': {'XPOS': -1.0,
+          'YPOS': 0.3},
+'target': {'XPOS': 0.0,
+           'YPOS': -1.8}
+}
+
+plt.ion()
+fig=plt.figure()
+# Aargh - used old, wrong value
+#AZ=-370.83 
+AZ=-320.42
+BZ=-224.97
+CZ=-129.314
+ZVEC=np.array([AZ,BZ,CZ])
+
+NOMINAL_ZVEC=np.array([AZ,CZ])
+NOMINAL_XVEC=np.array([NOMINAL_POS['3H07A']['XPOS'], 
+    NOMINAL_POS['3H07C']['XPOS']])
+NOMINAL_YVEC=np.array([NOMINAL_POS['3H07A']['YPOS'], 
+    NOMINAL_POS['3H07C']['YPOS']])
+
+#XSLOPE=np.array([-0.9843,-1.116,-0.9645])
+#XOFF=np.array([-0.05355,0.08082,-0.893])
+
+#YSLOPE=np.array([0.9687,1.166,0.8703])
+#YOFF=np.array([-0.2027,0.3751,0.466])
+
+# Fit fro DEC 2019 BPm vs HARP scan . In fit made the Harp X coordinate the same sign as the X epics
+XSLOPE=np.array([0.973,1.096,0.956])
+XOFF=np.array([-.076,-0.200,0.937])
+
+YSLOPE=np.array([0.957,1.16,0.855])
+YOFF=np.array([0.17,-0.06,-0.94])
+
+# changes to offsets from most recent survey
+#DY=np.array([0.25,0.22,0.20])
+#DX=np.array([-0.02,-0.09,-0.24])
+
+# survey numbers above were updated..
+#DY=np.array([0.21,0.21,0.29])
+#DX=np.array([0.06,0.05,-0.09])
+
+# from comparison to harp scans 8/29
+#DY=np.array([0.36,0.36,0.31])
+#DX=np.array([0.21,0.12,-0.29])
+
+# from survey, Feb. 2019 (note these are relative to 2017 survey offsets)
+#DY=np.array([-0.05,-0.23,-0.46])
+#DX=np.array([0.16,0.21,0.03])
+
+# from HARP/BPM comparison, Feb. 9 2019 !!!!WRONG!!!
+#DY=np.array([0.26,0.21,0.06])
+#DX=np.array([0.32,0.38,0.09])
+
+# from HARP/BPM comparison, Feb. 9 2019
+#DY=np.array([0.38,-0.13,-0.76])
+#DX=np.array([0.32,0.38,0.09])
+
+# For Dec 2019 comparison did a BPM vs HARP scan so do not need offset
+DY=np.array([0.,0.,0.])
+DX=np.array([0.,0.,0.])
+
+
+XOFF=XOFF+DX
+YOFF=YOFF+DY
+
+EX=np.array([0.1,0.1,0.1])
+EY=np.array([0.1,0.1,0.1])
+
+# fit line using both x and y errors
+def f(B, x):
+    '''Linear function y = m*x + b'''
+    # B is a vector of the parameters.
+    # x is an array of the current x values.
+    # x is in the same format as the x passed to Data or RealData.
+    #
+    # Return an array in the same format as y passed to Data or RealData.
+    return B[0]*x + B[1]
+
+
+ax1 = fig.add_subplot(221)
+
+ax2 = fig.add_subplot(223)
+ax2.set_xlabel('Distance from target (cm)',fontsize=15)
+
+ax3 = fig.add_subplot(122)
+
+while True:
+    time.sleep(0.5)
+    #read in epics stuff
+    ibcm=float(check_output("caget -w5 IBC3H00CRCUR4 | awk \ '{print $2}'", shell=True).strip())
+
+    AXraw=float(check_output("caget -w5 IPM3H07A.XPOS | awk \ '{print $2}'", shell=True).strip())
+    BXraw=float(check_output("caget -w5 IPM3H07B.XPOS | awk \ '{print $2}'", shell=True).strip())
+    CXraw=float(check_output("caget -w5 IPM3H07C.XPOS | awk \ '{print $2}'", shell=True).strip())
+
+    AYraw=float(check_output("caget -w5 IPM3H07A.YPOS | awk \ '{print $2}'", shell=True).strip())
+    BYraw=float(check_output("caget -w5 IPM3H07B.YPOS | awk \ '{print $2}'", shell=True).strip())
+    CYraw=float(check_output("caget -w5 IPM3H07C.YPOS | awk \ '{print $2}'", shell=True).strip())
+
+    XRAW=np.array([AXraw,BXraw,CXraw])
+    YRAW=np.array([AYraw,BYraw,CYraw])
+
+    if ibcm > 0.5:
+        XVEC=(XRAW*XSLOPE+XOFF)
+        YVEC=YRAW*YSLOPE+YOFF
+    else:
+        XVEC=XRAW*0.0
+        YVEC=YRAW*0.0
+
+    # Show X Position
+    # define 1st plot -x vs. z
+    ax1.plot(NOMINAL_ZVEC, NOMINAL_XVEC, marker='o', color='black', markersize=8,
+            fillstyle='none', linestyle='none')
+    ax1.errorbar(ZVEC, XRAW, yerr=EX, fmt='o', color='blue', markersize=4)
+    ax1.errorbar(ZVEC, XVEC, yerr=EX, fmt='o', color='red', markersize=4)
+    linear = Model(f)
+    mydata = RealData(ZVEC, XVEC, sy=EX)
+    myodr = ODR(mydata, linear, beta0=[1., 2.])
+    myoutput = myodr.run()
+    #    myoutput.pprint()
+    # get errors from covariance matrix - standard errors from ODR amplify by chi-squared
+    xtar=round(myoutput.beta[1],2)
+
+    # Plot the fit
+    x_fit = np.linspace(-400, 50, 1000)
+    y_fit = f(myoutput.beta, x_fit)
+    ax1.plot(x_fit, y_fit, color='black')
+
+    #plt.xlabel('Distance from target (cm)',fontsize=15)
+    #ax1.text(-400,2,"X Pos @ tgt: "+str(xtar))
+    ax1.set_ylim([-4,4])
+
+    # Show Y Position
+    ax2.plot(NOMINAL_ZVEC, NOMINAL_YVEC, marker='o', color='black', markersize=8,
+            fillstyle='none', linestyle='none')
+    yraw = ax2.errorbar(ZVEC, YRAW, yerr=EX, fmt='o', color='blue', markersize=4)
+    ycor = ax2.errorbar(ZVEC, YVEC, yerr=EX, fmt='o', color='red', markersize=4)
+    linear = Model(f)
+    mydata = RealData(ZVEC, YVEC, sy=EX)
+    myodr = ODR(mydata, linear, beta0=[1., 2.])
+    myoutput = myodr.run()
+    #    myoutput.pprint()
+    # get errors from covariance matrix - standard errors from ODR amplify by chi-squared
+    ytar=round(myoutput.beta[1],2)
+
+    # Plot the fit
+    x_fit = np.linspace(-400, 50, 1000)
+    y_fit = f(myoutput.beta, x_fit)
+    ax2.plot(x_fit, y_fit, color='black')
+
+    #ax2.text(-400,2,"Y Pos @ tgt: "+str(ytar));
+    ax2.set_ylim([-4,4])
+    ax1.set_ylabel('X (mm)',fontsize=15)
+    ax2.set_ylabel('Y (mm)',fontsize=15)
+    ax3.set_xlabel('Target X (mm)',fontsize=15)
+    ax3.set_ylabel('Target Y (mm)',fontsize=15)
+
+    # BPM monitoring-like Output
+    ax3.plot([NOMINAL_POS['target']['XPOS']], [NOMINAL_POS['target']['YPOS']], marker='o', color='black', markersize=8, fillstyle='none', linestyle='none')
+    ax3.errorbar(xtar,ytar, yerr=0, fmt='o', color='red', markersize=4)
+    ax3.yaxis.tick_right()
+    ax3.yaxis.set_label_position('right')
+    ax3.yaxis.grid(color='grey')
+    ax3.xaxis.grid(color='grey')
+    ax3.text(-3,1.2,"Beam position on target:\n"+8*" "+str(xtar)+" , "+str(ytar));
+    ax3.set_xlim([-4,4])
+    ax3.set_ylim([-6,2])
+    ax3.legend((yraw, ycor), ('raw', 'corrected'), loc='lower left')
+
+    print('xtar is', xtar, 'ytar', ytar)
+
+    ax3.text(-13.5, 2.5,"Nominal raw 3H07A: ({}, {}), 3H07C: ({}, {})".format(
+        NOMINAL_POS['3H07A']['XPOS'],
+        NOMINAL_POS['3H07A']['YPOS'],
+        NOMINAL_POS['3H07C']['XPOS'],
+        NOMINAL_POS['3H07C']['YPOS']),
+        size=15)
+
+    fig.canvas.draw()
+    fig.canvas.flush_events()
+    plt.pause(1)
+    ax1.clear()
+    ax2.clear()
+    ax3.clear()
+    gc.collect()

CALIBRATION/hms_cal_calib/THcShTrack.h

@@ -48,6 +48,7 @@ class THcShTrack {
   void Print(ostream & ostrm);
 
   void SetEs(Double_t* alpha);
+  void SetEsNoCor(Double_t* alpha);
 
   Double_t Enorm();
   Double_t EPRnorm();
@@ -172,10 +173,13 @@ void THcShTrack::SetEs(Double_t* alpha) {
     Double_t xh=X+Xp*(ncol-0.5)*fZbl;
     Double_t yh=Y+Yp*(ncol-0.5)*fZbl;
     if (nblk <= fNnegs) {
+      //      cout << adc_pos <<"\t"<< Ycor(yh,0) <<"\t"<< alpha[nblk-1] << endl;
+      //      cout << adc_neg <<"\t"<< Ycor(yh,1) <<"\t"<< alpha[fNblks+nblk-1] << endl;
       (*iter)->SetEpos(adc_pos*Ycor(yh,0)*alpha[nblk-1]);
       (*iter)->SetEneg(adc_neg*Ycor(yh,1)*alpha[fNblks+nblk-1]);
     }
     else {
+      //      cout << adc_pos <<"\t"<< Ycor(yh) <<"\t"<< alpha[nblk-1] << endl;
       (*iter)->SetEpos(adc_pos*Ycor(yh)*alpha[nblk-1]);
       (*iter)->SetEneg(0.);
     };
@@ -186,6 +190,35 @@ void THcShTrack::SetEs(Double_t* alpha) {
 
 //------------------------------------------------------------------------------
 
+void THcShTrack::SetEsNoCor(Double_t* alpha) {
+
+  // Set hit energy depositions seen from postive and negative sides,
+  // by use of calibration (gain) constants alpha.
+  
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+  
+    Double_t adc_pos = (*iter)->GetADCpos();
+    Double_t adc_neg = (*iter)->GetADCneg();
+    UInt_t nblk = (*iter)->GetBlkNumber();
+
+    Int_t ncol=(nblk-1)/fNrows+1;
+    Double_t xh=X+Xp*(ncol-0.5)*fZbl;
+    Double_t yh=Y+Yp*(ncol-0.5)*fZbl;
+    if (nblk <= fNnegs) {
+      (*iter)->SetEpos(adc_pos*alpha[nblk-1]);
+      (*iter)->SetEneg(adc_neg*alpha[fNblks+nblk-1]);
+    }
+    else {
+      (*iter)->SetEpos(adc_pos*alpha[nblk-1]);
+      (*iter)->SetEneg(0.);
+    };
+
+  };
+
+}
+
+//------------------------------------------------------------------------------
+
 Double_t THcShTrack::Enorm() {
 
   // Normalized to the track momentum energy depostion in the calorimeter.

CALIBRATION/hms_cal_calib/hcal_calib.cpp

 #include "TCanvas.h"
+#include "TLine.h"
 #include <TStyle.h>
 #include <TH1.h>
+#include <TEllipse.h>
 #include <TF1.h>
+#include <TPaveText.h>
 #include "THcShowerCalib.h"
 
 //
 // A steering Root script for the HMS calorimeter calibration.
 //
 
-void hcal_calib(string RunNumber) {
+void hcal_calib(string Prefix, int nstop=-1, int nstart=0) {
+
+  bool DRAW = 1;  //flag to draw extra plots
 
   // Initialize the analysis clock
   clock_t t = clock();
   
- cout << "Calibrating run " << RunNumber << endl;
+  cout << "Calibrating file " << Prefix << ".root, events "
+       << nstart << " -- " << nstop << endl;
 
- THcShowerCalib theShowerCalib(RunNumber);
+  THcShowerCalib theShowerCalib(Prefix, nstart, nstop);
 
+ theShowerCalib.ReadThresholds();  // Read in threshold param-s and intial gains
  theShowerCalib.Init();            // Initialize constants and variables
  theShowerCalib.CalcThresholds();  // Thresholds on the uncalibrated Edep/P
  theShowerCalib.ComposeVMs();      // Compute vectors amd matrices for calib.
@@ -24,43 +31,267 @@ void hcal_calib(string RunNumber) {
  theShowerCalib.SaveAlphas();      // Save the constants
  // theShowerCalib.SaveRawData();  // Save raw data into file for debug purposes
  theShowerCalib.FillHEcal();       // Fill histograms
+ // theShowerCalib.FillHEcalNoCor();  // Removes Y correction (pulseInt * gain * yCor) = E
+ theShowerCalib.FillCutBranch();       // Fill diagnostic histos
+ theShowerCalib.FillHitsGains();       // Plot hits and gains
 
  // Plot histograms
-
+ //___________Canvas 1_______________
  TCanvas* Canvas =
-   new TCanvas("Canvas", "HMS Shower Counter calibration", 1000, 667);
- Canvas->Divide(2,2);
-
- Canvas->cd(1);
+ new TCanvas("Canvas", "HMS Shower Counter calibration", 1000, 667);
+ Canvas->Divide(3,2);
 
  // Normalized uncalibrated energy deposition.
-
+ Canvas->cd(1);
  theShowerCalib.hEunc->DrawCopy();
-  
  theShowerCalib.hEuncSel->SetFillColor(kGreen);
  theShowerCalib.hEuncSel->DrawCopy("same");
 
  // E(Tail) vs E(Preshower) plot.
-
  Canvas->cd(2);
  theShowerCalib.hETAvsEPR->Draw("colz");
 
- // Normalized energy deposition after calibration.
-
+ // Uncalibrated e/p at calorimeter
  Canvas->cd(3);
- gStyle->SetOptFit();
+ //  theShowerCalib.hCaloPosNormU->GetZaxis()->SetRangeUser(0.5,1.5);
+ theShowerCalib.hCaloPosNormU->Draw("COLZ");
+ theShowerCalib.hCaloPosNormU->SetStats(0);
+ theShowerCalib.hCaloPosNormU->GetZaxis()->SetRangeUser(0.5,1.5);
 
- theShowerCalib.hEcal->Fit("gaus","","",0.95,1.05);
+ // Normalized energy deposition after calibration.
+ Canvas->cd(4);
+ gStyle->SetOptFit();
+ //______wph Fit over range centered by max bin_________
+ //**commented out next line to improve location of fit mean**
+ // theShowerCalib.hEcal->Fit("gaus","O","",0.5,1.5); 
+ Double_t maxBin= theShowerCalib.hEcal->GetMaximumBin();
+ Double_t maxValue= theShowerCalib.hEcal->GetBinCenter(maxBin);
+ theShowerCalib.hEcal->Fit("gaus","O","",maxValue-.1, maxValue+.2);
+ TF1 *fit = theShowerCalib.hEcal->GetFunction("gaus");
+ Double_t gmean  = fit->GetParameter(1);
+ Double_t gsigma = fit->GetParameter(2);
+ double gLoThr = gmean - 1.*gsigma;
+ double gHiThr = gmean + 2.*gsigma;
+ cout << "gLoThr=" << gLoThr << "  gHiThr=" << gHiThr << endl;
+ theShowerCalib.hEcal->Fit("gaus","","",gLoThr,gHiThr);
  theShowerCalib.hEcal->GetFunction("gaus")->SetLineColor(2);
  theShowerCalib.hEcal->GetFunction("gaus")->SetLineWidth(1);
  theShowerCalib.hEcal->GetFunction("gaus")->SetLineStyle(1);
+ // theShowerCalib.hEcalNoCor->SetLineColor(kGreen);
+ // theShowerCalib.hEcalNoCor->Draw("same");
 
  // HMS delta(P) versus the calibrated energy deposition.
-
- Canvas->cd(4);
+ Canvas->cd(5);
  theShowerCalib.hDPvsEcal->Draw("colz");
 
+ 
+ // Canvas->cd(5);
+ // theShowerCalib.hCaloPosWt->Draw("COLZ");
+ Canvas->cd(6);
+ theShowerCalib.hCaloPosNorm->GetZaxis()->SetRangeUser(0.05,1.2);
+ theShowerCalib.hCaloPosNorm->SetTitle("Normalized E/p at Calorimeter");
+  theShowerCalib.hCaloPosNorm->Draw("COLZ");
+ 
+ // Save canvas in a pdf format.
+ Canvas->Print(Form("%s_%d_%d.pdf",Prefix.c_str(),nstart,nstop));
+
+// A bunch of diagnostic plots added Oct. 2019
+//wph
+ if(DRAW==1){ 
+ TCanvas* Canvas2 =
+   new TCanvas("Canvas2", "Hits, gains, and Projections", 1000, 667);
+
+ TCanvas* Canvas3 =
+   new TCanvas("Canvas3", "Branches for cuts", 1000, 667);
+
+ TCanvas* Canvas4 =
+   new TCanvas("Canvas4", "Pulse Integrals, x/y_calo vs E/p", 1000, 667);
+
+ //___________Canvas 2_______________
+ Canvas2->Divide(3,2);
+ // Projection at calorimeter
+ Canvas2->cd(1);
+ theShowerCalib.hCaloPos->Draw("COLZ");
+
+ // Projection at dipole
+ Canvas2->cd(2);
+ TEllipse *el1=new TEllipse( 0, 2.8, 46.507, 46.507);
+ el1->SetLineColor(kRed);
+ theShowerCalib.hExitPos->Draw("COLZ");
+ el1->Draw("same");
+ theShowerCalib.hExitPos->Draw("COLZ SAME");
+
+ // Hits and gains for each block by layer
+ Canvas2->cd(3);
+ theShowerCalib.pr1->SetStats(0);
+ theShowerCalib.pr1->Draw("TEXT COLZ");
+ theShowerCalib.pr1a->Draw("TEXT SAME");
+ Canvas2->cd(4);
+ theShowerCalib.ta2->Draw("TEXT COLZ");
+ theShowerCalib.ta2->SetStats(0);
+ theShowerCalib.ta2a->Draw("TEXT SAME");
+ Canvas2->cd(5);
+ theShowerCalib.ta3->Draw("TEXT COLZ");
+ theShowerCalib.ta3->SetStats(0);
+ theShowerCalib.ta3a->Draw("TEXT SAME");
+ Canvas2->cd(6);
+ theShowerCalib.ta4->SetStats(0);
+ theShowerCalib.ta4->Draw("TEXT COLZ");
+ theShowerCalib.ta4a->Draw("TEXT SAME");
+ // Set range for hits
+ Double_t maxHit=1.3*theShowerCalib.pr1->GetMaximum();
+ theShowerCalib.pr1->GetZaxis()->SetRangeUser(0,maxHit);
+ theShowerCalib.ta2->GetZaxis()->SetRangeUser(0,maxHit);
+ theShowerCalib.ta3->GetZaxis()->SetRangeUser(0,maxHit);
+ theShowerCalib.ta4->GetZaxis()->SetRangeUser(0,maxHit);
+ //gains font color
+ theShowerCalib.pr1a->SetMarkerColor(kRed);
+ theShowerCalib.ta2a->SetMarkerColor(kRed);
+ theShowerCalib.ta3a->SetMarkerColor(kRed);
+ theShowerCalib.ta4a->SetMarkerColor(kRed);
+
+ //___________Canvas 3_______________
+ Canvas3->Divide(3,2);
+ //Cherenkov Cut
+ Canvas3->cd(1);
+ theShowerCalib.hCer->Draw();
+ gPad->SetLogy();
+ Double_t top=theShowerCalib.hCer->GetMaximum();
+ Double_t locx=theShowerCalib.GetCerMin();
+ TLine *cerl=new TLine(locx,0,locx,top);
+ cerl->SetLineColor(kMagenta);
+ cerl->Draw("same");
+
+ // H.gtr.p
+ Canvas3->cd(2);
+ theShowerCalib.hP->Draw();
+
+ //Delta Cut 
+ Canvas3->cd(3);
+ theShowerCalib.hDelta->Draw();
+ top=theShowerCalib.hDelta->GetMaximum();
+ locx=theShowerCalib.GetDeltaMin();
+ TLine *dplmin=new TLine(locx,0,locx,top);
+ dplmin->SetLineColor(kMagenta);
+ dplmin->Draw("same");
+ locx=theShowerCalib.GetDeltaMax();
+ TLine *dplmax=new TLine(locx,0,locx,top);
+ dplmax->SetLineColor(kMagenta);
+ dplmax->Draw("same");
+
+ //Beta Cut 
+ Canvas3->cd(4);
+ theShowerCalib.hBeta->Draw();
+ gPad->SetLogy();
+ top=theShowerCalib.hBeta->GetMaximum();
+ locx=theShowerCalib.GetBetaMin();
+ TLine *betalmin=new TLine(locx,0,locx,top);
+ betalmin->SetLineColor(kMagenta);
+ betalmin->Draw("same");
+ locx=theShowerCalib.GetBetaMax();
+ TLine *betalmax=new TLine(locx,0,locx,top);
+ betalmax->SetLineColor(kMagenta);
+ betalmax->Draw("same");
+ // ncluster and ntrack
+ Canvas3->cd(5);
+ theShowerCalib.hClusTrk->Draw("BOX");
+ //gPad->SetLogy();
+ // theShowerCalib.hNclust->Draw();
+ // Canvas3->cd(6);
+ // gPad->SetLogy();
+ // theShowerCalib.hNtrack->Draw();
+
+ Canvas3->cd(6);
+  TPaveText *pt=new TPaveText(0.1,0.1,.9,.9);
+  pt->AddText("Percentage of Events passing cuts");
+  Double_t val= theShowerCalib.GetRatio()*100;
+  pt->AddText(Form("All Cuts:                    %.2f%%",val));
+  //__________________________________________________________
+  Double_t fb=theShowerCalib.GetCerMin();
+  Double_t bmin=theShowerCalib.hCer->GetXaxis()->FindBin(fb);
+  fb=theShowerCalib.hCer->GetXaxis()->GetXmax();
+  Double_t bmax=theShowerCalib.hCer->GetXaxis()->FindBin(fb);
+  Double_t intg=theShowerCalib.hCer->Integral(bmin,bmax);
+  Double_t nentries=theShowerCalib.hCer->GetEntries();
+  val=intg/nentries*100;
+  pt->AddText(Form("Cherenkov:                   %.2f%%",val));
+  //__________________________________________________________
+  fb=theShowerCalib.GetBetaMin();
+  bmin=theShowerCalib.hBeta->GetXaxis()->FindBin(fb);
+  fb=theShowerCalib.GetBetaMax();
+  bmax=theShowerCalib.hBeta->GetXaxis()->FindBin(fb);
+  intg=theShowerCalib.hBeta->Integral(bmin,bmax);
+  nentries=theShowerCalib.hBeta->GetEntries();
+  val=intg/nentries*100;
+  pt->AddText(Form("Beta:                        %.2f%%",val));
+  //__________________________________________________________
+  fb=theShowerCalib.GetDeltaMin();
+  bmin=theShowerCalib.hDelta->GetXaxis()->FindBin(fb);
+  fb=theShowerCalib.GetDeltaMax();
+  bmax=theShowerCalib.hDelta->GetXaxis()->FindBin(fb);
+  intg=theShowerCalib.hDelta->Integral(bmin,bmax);
+  nentries=theShowerCalib.hDelta->GetEntries();
+  val=intg/nentries*100;
+  pt->AddText(Form("Delta:                       %.2f%%",val));
+  //__________________________________________________________
+  Double_t bx=theShowerCalib.hClusTrk->GetXaxis()->FindBin(1);
+  Double_t by=theShowerCalib.hClusTrk->GetYaxis()->FindBin(1);
+  Double_t tot=theShowerCalib.hClusTrk->GetBinContent(bx,by);
+  nentries=theShowerCalib.hClusTrk->GetEntries();
+  val=tot/nentries*100;
+  pt->AddText(Form("nClusters==1 && nTracks==1:  %.2f%%",val));
+  //__________________________________________________________
+  pt->Draw();
+
+ //___________Canvas 4_______________
+ Canvas4->Divide(2,2);
+
+ //Pulse Intregral vs Block Number
+ Canvas4->cd(1);
+ theShowerCalib.pmtList->Draw("COLZ");
+
+ // All 78 pulse integrals
+ Canvas4->cd(2);
+ Float_t adcmax=0;
+ Float_t temp=0;
+ for(UInt_t i=1; i<78; i++){
+   temp=theShowerCalib.hAdc[i]->GetMaximum();
+   if (temp>adcmax)adcmax=temp;}
+ theShowerCalib.hAdc[0]->Draw();
+ theShowerCalib.hAdc[0]->GetYaxis()->SetRangeUser(0,adcmax);
+ for(UInt_t i=1; i<78; i++){theShowerCalib.hAdc[i]->Draw("same");}
+
+ // Y_calo vs E/p
+ Canvas4->cd(3);
+ theShowerCalib.yCalVsEp->Draw("COLZ");
+
+ // X_calo vs E/p
+ Canvas4->cd(4);
+ theShowerCalib.xCalVsEp->Draw("COLZ");
+ //  Lines representing block edges
+ Int_t nRows=13;
+ Double_t thkRow=10;
+ Double_t startRow=-70.4;
+ TLine *t[nRows+1];
+ for (Int_t i=0;i<=nRows;i++)
+ {
+   Double_t y= startRow+i*thkRow;
+   t[i]=new TLine(0.1,y,1.4,y);
+   t[i]->SetLineColor(kRed);
+   t[i]->Draw("same");
+ }
+ 
+ // Save canvases
+ Canvas2->Print(Form("PDFs/hits_%s_%d_%d.pdf",Prefix.c_str(),nstart,nstop));
+ Canvas3->Print(Form("PDFs/cuts_%s_%d_%d.pdf",Prefix.c_str(),nstart,nstop));
+ Canvas3->Print(Form("PDFs/pInt_%s_%d_%d.pdf",Prefix.c_str(),nstart,nstop));
+ }
+
  // Calculate the analysis rate
+ //   TFile *oFile=new TFile("hcal_out.root","RECREATE");
+ //  theShowerCalib.hEcal->Write();
+ //  oFile->Close();
  t = clock() - t;
  printf ("The analysis took %.1f seconds \n", ((float) t) / CLOCKS_PER_SEC);
+
 }

CALIBRATION/hms_cal_calib/howto.txt

 How to calibrate HMS calorimeter with real electrons.
 
-The calibration scripts reside in
-hallc_replay/CALIBRATION/hms_cal_calib directory. They consist of 3
-header files called THcShHit.h,THcShTrack.h, THcShowerCalib.h, and a
-steering script hcal_calib.cpp.
+The calibration code resides in hallc_replay/CALIBRATION/hms_cal_calib
+directory. It consists of 3 header files called
+THcShHit.h,THcShTrack.h, THcShowerCalib.h, and a steering script
+hcal_calib.cpp. There is also an input file input.dat containing
+threshold parameters and initial gain constants necessary for
+calibration.
 
 The scripts work on root files from hcana analysis and make use of
 quantities pertained to tracking, gas Cherenkov, and TOF from
 hodoscopes. Hence it is convenient to calibrate on root files from
-full HMS analysis. The root files are assumed in a linked ROOTfiles
-directory, under names hms_replay_<run_number>.root.
+full HMS analysis. The root files are assumed to be in a linked
+ROOTfiles directory.
+
+Th input file contains thresholds on the quantities used in the
+calibration, such as HMS: delta range, the hodoscope beta range,
+threshold on the gas Cherenkov signals, required minimum number of
+hits in the calorimeter channels in order to be calibrated, initial
+gain constants, and parameters for an auxiliary histogram of energy
+deposition calculated with initial gain constants. The gas Cherenkov
+signals and beta are used to select good electron events for
+calibration.  Additional clean up of the electron sample from hadron
+contamination (if any) left after PID with gas Cherenkov and beta is
+done by a rough localization of the electron peak in the auxiliary
+histogram. Note that code does not iterate the gain constants. The
+user is free to modify cuts in the input file as needed. Caution must
+be exercised not to alter the format of the file.
+
+The steering script hcal_calib.cpp takes 3 parameters: prefix of the
+root file to be calibrated, last (-1 by default) and the first (0 by
+default) events to be used for calibration. The -1 value for last
+event number defaults to the number of entries in the root file.
 
 Once your hcana, hallc_replay and Root are set up, you can compile and
 run hcal_calib under hcana, by issuing command
 
-.x hcal_calib.cpp+("run number").
+.x hcal_calib.cpp+("Prefix"),
 
-For instance, for calibration on hms_replay_303.root file in
-ROOTfiles, the correct commad would be
+where Prefix is the prefix part in the name of root file.
 
-.x hcal_calib.cpp+("303") .
+For instance, for calibration with hms_replay_303.root file in
+ROOTfiles, the correct command would be
+
+.x hcal_calib.cpp+("hms_replay_303") .
 
 Upon calibration, a canvas with representative plots will pop up. The
 calibration constants will be written in output file
-hcal.param.<run_number>, in a format suitable for plugging them into
-your hcal.param file for subsequent use.
-
-If you want to modify selection cuts used in calibration (cuts on
-delta, beta, gas Cherenkov signals), you can find them at the
-beginning of THcShowerCalib.h file, in the #define directives.
+hcal.param.<Prefix>, in a format suitable for plugging them into your
+hcal.param file for subsequent use. The representative canvas is saved
+in a pdf file <Prefix>.pdf.

CALIBRATION/hms_cal_calib/input.dat

+-10 10	Delta range, %
+0.5 1.5	Beta range
+1.5  	Gas Cherenkov, threshold on signals in p.e.
+10	Minimum number of hits per channel required to be calibrated
+0.05 2.	Range of uncalibrated energy deposition histogram
+500	Binning of uncalibrated energy deposition histogram
+-0.1 0.1	Gaussian fit range around mean e- peak in the uncalibrated Edep histogram 
+2.	Sigma width to select events to use in calibration (green fill)
+
+; Calibration constants for file hms_ecal.root (HMS 2489-2494), 164827 events processed
+
+hcal_pos_gain_cor= 12.99,  7.13,  9.22, 10.26, 10.82, 12.49, 11.93, 12.39, 10.28, 15.14, 14.59, 12.85,  6.40,
+                   10.87, 12.23,  9.46, 13.51,  8.71,  6.30,  7.73,  7.41,  8.94, 11.40, 11.51, 12.40, 14.86,
+                   24.29, 14.78, 18.17, 21.53, 16.96, 18.97, 23.00, 19.51, 21.92, 25.86, 18.37, 22.29, 19.84,
+                   32.17, 16.75, 20.16, 18.54, 17.89, 19.25, 22.14, 17.51, 19.67, 20.51, 18.37, 18.86, 26.46,
+hcal_neg_gain_cor= 15.11, 14.11, 14.42, 10.92, 11.88, 13.32, 14.58, 18.16, 12.63, 11.47, 10.60, 11.41, 19.57,
+                   14.71, 13.20, 14.22, 14.34, 16.11, 16.93, 19.64, 17.34, 18.06, 11.49, 14.83, 13.72,  8.86,
+                    0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+                    0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+

CALIBRATION/hms_cer_calib/README

+# Calibrating the Cherenkov
+* To automatically run the calibration
+```
+root -l "run_cal.C(RunNumber,NumEvents,COIN)"
+```
+* Where COIN == 1 indicates the file was a full coincident replay
+* If options are left blank, user will be prompted for a value
+* To manually run the scripts:
+* Link your ROOT file
+* Run the script with options (listed after script is executed)
+```
+root -l ../../ROOTfiles/hms_replay_488_-1.root
+T->Process("calibration.C+", "options");
+```
+* Or, if you want to run using PROOF
+```
+root -l
+TChain ch("T");
+ch.Add("/path/to/ROOTfile");
+TProof::Open(""); \\For PROOF-lite, adjust accordingly for PROOF
+ch.SetProof();
+ch.Process("calibration.C+", "options");
+```
+### Options for the calibration script are (case/spelling important):
+* **showall** - display all calibration details (lots of windows)

CALIBRATION/hms_cer_calib/calibration.C

+#define calibration_cxx
+// The class definition in calibration.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("calibration.C")
+// root> T->Process("calibration.C","some options")
+// root> T->Process("calibration.C+")
+//
+
+
+#include "calibration.h"
+#include <TF1.h>
+#include <TH1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TSpectrum.h>
+#include <TPolyMarker.h>
+#include <TPaveStats.h>
+
+void calibration::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).
+
+   TString option = GetOption();
+   Info("Begin", "Starting calibration process with option: %s", option.Data());
+   Info("Begin", "To see details of calibration, use option showall");
+
+   if (option.Contains("showall")) fFullShow = kTRUE;
+}
+
+void calibration::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();
+  
+  fPulseInt = new TH1F*[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    fPulseInt[ipmt] = new TH1F(Form("PulseInt%d",ipmt+1), Form("Pulse integral PMT %d;ADC Channel (pC);Counts",ipmt+1), 1000, 0, 100);
+    GetOutputList()->Add(fPulseInt[ipmt]);
+  }
+  
+  fPulseInt_quad = new TH1F**[2];
+  fPulseInt_quad[0] = new TH1F*[4];
+  fPulseInt_quad[1] = new TH1F*[4];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    for (Int_t iquad = 0; iquad < 4; iquad++) {
+      fPulseInt_quad[ipmt][iquad] = 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),1000,0,100);
+      GetOutputList()->Add(fPulseInt_quad[ipmt][iquad]);
+    }
+  }
+  
+  //Timing and Beta cut visualizations
+  fBeta_Cut = new TH1F("Beta_Cut", "Beta cut used for 'good' hits;Beta;Counts", 1000, -5, 5);
+  GetOutputList()->Add(fBeta_Cut);
+  fBeta_Full = new TH1F("Beta_Full", "Full beta for events;Beta;Counts", 1000, -5, 5);
+  GetOutputList()->Add(fBeta_Full);
+
+  fTiming_Cut = new TH1F*[2];
+  fTiming_Full = new TH1F*[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    fTiming_Cut[ipmt] = new TH1F(Form("Timing_Cut%d",ipmt+1), Form("Timing cut used for 'good' hits in PMT %d;Time (ns);Counts",ipmt+1), 10000, 50, 150);
+    GetOutputList()->Add(fTiming_Cut[ipmt]);
+    fTiming_Full[ipmt] = new TH1F(Form("Timing_Full%d",ipmt+1), Form("Full timing information for events in PMT %d;Time (ns);Counts",ipmt+1), 10000, 50, 150);
+    GetOutputList()->Add(fTiming_Full[ipmt]);
+  }
+
+  //Particle ID cut visualization
+  fCut_everything = new TH1F("Cut_everything", "Visualization of no cuts; Normalized Calorimeter Energy (GeV); Counts", 200, 0, 2.0);
+  GetOutputList()->Add(fCut_everything);
+  fCut_electron = new TH1F("Cut_electron", "Visualization of pion cut; Normalized Calorimeter Energy (GeV); Counts", 200, 0, 2.0);
+  GetOutputList()->Add(fCut_electron);
+}
+
+Bool_t calibration::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);
+   
+   //Only one track
+   if (*Ndata_H_tr_beta != 1) return kTRUE;
+
+   for (Int_t itrack = 0; itrack < *Ndata_H_tr_beta; itrack++) {
+     //Beta Cut
+     fBeta_Full->Fill(H_tr_beta[itrack]);
+     if (TMath::Abs(H_tr_beta[itrack] -1.0) > 0.2) return kTRUE;
+     fBeta_Cut->Fill(H_tr_beta[itrack]);
+
+     for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+       //Timing Cut
+       fTiming_Full[ipmt]->Fill(H_cer_goodAdcTdcDiffTime[ipmt]);
+       if (H_cer_goodAdcTdcDiffTime[ipmt] > 110 || H_cer_goodAdcTdcDiffTime[ipmt] < 90) continue;
+       fTiming_Cut[ipmt]->Fill(H_cer_goodAdcTdcDiffTime[ipmt]);
+
+       if (ipmt == 0) fCut_everything->Fill(*H_cal_etotnorm);
+       //Electron Cut
+       if (*H_cal_etotnorm < 0.5) {
+	 if (ipmt == 0) fCut_electron->Fill(*H_cal_etotnorm);
+	 fPulseInt[ipmt]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+
+	 //Quadrant Cuts
+	 if (*H_cer_xAtCer >= 0.0 && *H_cer_yAtCer >= 0.0) {
+	   fPulseInt_quad[ipmt][0]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer >= 0.0 && *H_cer_yAtCer < 0.0) {
+	   fPulseInt_quad[ipmt][1]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer < 0.0 && *H_cer_yAtCer >= 0.0) {
+	   fPulseInt_quad[ipmt][2]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer < 0.0 && *H_cer_yAtCer < 0.0) {
+	   fPulseInt_quad[ipmt][3]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+       }//End of electron cut
+     }//End of PMT loop
+   }//End of tracking loop
+   
+   return kTRUE;
+}
+
+void calibration::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 calibration::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.
+
+  Info("Terminate", "'%s' showing", (fFullShow ? "full" : "minimal"));
+  Info("Terminate", "Histograms formed, now starting calibration.\n'Peak Buffer full' is a good warning!\n");
+  printf("\n");
+
+  //Have to extract the histograms from the OutputList 
+  TH1F* PulseInt[2];
+  TH1F* PulseInt_quad[2][4];
+  TH1F* Timing_Cut[2];
+  TH1F* Timing_Full[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    Timing_Cut[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("Timing_Cut%d",ipmt+1)));
+    Timing_Full[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("Timing_Full%d",ipmt+1)));
+    PulseInt[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("PulseInt%d",ipmt+1)));
+    for (Int_t iquad = 0; iquad < 4; iquad++) {
+      PulseInt_quad[ipmt][iquad] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("PulseInt_quad%d_PMT%d",iquad+1,ipmt+1)));
+    }
+  }
+
+  Double_t Cer_Peak[2];
+  //Begin peak Finding
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,8);
+    TSpectrum *s = new TSpectrum(1);
+    s->Search(PulseInt[ipmt], 1.0, "nobackground&&nodraw", 0.001);
+    TList *functions = PulseInt[ipmt]->GetListOfFunctions();
+    TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
+    Cer_Peak[ipmt] = *pm->GetX();
+    PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,100);
+  }
+
+  //Begin Fitting of SPE
+  /*TF1 *poisson =  new TF1("poisson","[1]*((pow([0],x)*exp(-[0]))/(tgamma(x+1)))",0,20);
+  poisson->SetParName(0,"#mu");
+  poisson->SetParName(1,"Amplitude");
+  poisson->SetParameters(5.0,100);
+  poisson->SetParLimits(0,0.0,10.0);
+  poisson->SetParLimits(1,0.0,9999.0);
+  poisson->SetRange(0,Cer_Peak[0]+2); PulseInt[0]->Fit("poisson","RQM");
+  TF1 *PulseIntPMT1Fit = PulseInt[0]->GetFunction("poisson");
+  poisson->SetRange(0,Cer_Peak[1]+2); PulseInt[1]->Fit("poisson","RQM");
+  TF1 *PulseIntPMT2Fit = PulseInt[1]->GetFunction("poisson");*/
+ 
+  TF1 *gaussian =  new TF1("gaussian","[1]*exp(-0.5*((x-[0])/[2])*((x-[0])/[2]))",0,20);
+  gaussian->SetParName(0,"#mu");
+  gaussian->SetParName(1,"Amplitude");
+  gaussian->SetParName(2,"#sigma");
+  gaussian->SetParameters(5.0,100,1.0);
+  gaussian->SetParLimits(0,0.0,10.0);
+  gaussian->SetParLimits(1,0.0,9999.0);
+  gaussian->SetParLimits(2,0.0,10.0);
+  gaussian->SetRange(0,Cer_Peak[0]+2); PulseInt[0]->Fit("gaussian","RQM0");
+  TF1 *PulseIntPMT1Fit = PulseInt[0]->GetFunction("gaussian");
+  gaussian->SetRange(0,Cer_Peak[1]+2); PulseInt[1]->Fit("gaussian","RQM0");
+  TF1 *PulseIntPMT2Fit = PulseInt[1]->GetFunction("gaussian");
+  
+  if (fFullShow) {
+    TCanvas *Beta = new TCanvas("Beta","Beta cuts used");
+    Beta->Divide(2,1);
+    Beta->cd(1); fBeta_Full->Draw();
+    Beta->cd(2); fBeta_Cut->Draw();
+   
+    TCanvas *Cal = new TCanvas("Cal","Calorimeter cuts used");
+    Cal->Divide(2,1);
+    Cal->cd(1); fCut_everything->Draw();
+    Cal->cd(2); fCut_electron->Draw();
+
+    TCanvas *Timing = new TCanvas("Timing","Timing cuts used");
+    Timing->Divide(2,2);
+    Timing->cd(1); Timing_Full[0]->Draw();
+    Timing->cd(2); Timing_Cut[0]->Draw();
+    Timing->cd(3); Timing_Full[1]->Draw();
+    Timing->cd(4); Timing_Cut[1]->Draw();
+    /*
+    TCanvas *PulseIntPMT1 = new TCanvas("PulseIntPMT1","Good Pulse Integral from Cherenkov PMT 1");
+    PulseIntPMT1->Divide(2,2);
+    PulseIntPMT1->cd(1); PulseInt_quad[0][0]->Draw();
+    PulseIntPMT1->cd(2); PulseInt_quad[0][1]->Draw();
+    PulseIntPMT1->cd(3); PulseInt_quad[0][2]->Draw();
+    PulseIntPMT1->cd(4); PulseInt_quad[0][3]->Draw();
+
+    TCanvas *PulseIntPMT2 = new TCanvas("PulseIntPMT2","Good Pulse Integral from Cherenkov PMT 2");
+    PulseIntPMT2->Divide(2,2);
+    PulseIntPMT2->cd(1); PulseInt_quad[1][0]->Draw();
+    PulseIntPMT2->cd(2); PulseInt_quad[1][1]->Draw();
+    PulseIntPMT2->cd(3); PulseInt_quad[1][2]->Draw();
+    PulseIntPMT2->cd(4); PulseInt_quad[1][3]->Draw();
+    */
+    gStyle->SetOptStat(111);
+    gStyle->SetOptFit(0111);
+    TCanvas *cPulseInt = new TCanvas("cPulseInt","Good Pulse Integral in both Cherenkov PMTs");
+    cPulseInt->Divide(2,1);
+    cPulseInt->cd(1); PulseInt[0]->Draw(); PulseIntPMT1Fit->Draw("same");
+    cPulseInt->cd(2); PulseInt[1]->Draw(); PulseIntPMT2Fit->Draw("same");
+  }
+
+  //Output the actual calibration information
+  cout << "Calibration constants are \nPMT#: ADC-to-NPE Value\n" << endl;
+  cout << "PMT 1:" << setw(8) << Form("%3.3f\n",PulseIntPMT1Fit->GetParameter(0));
+  cout << "PMT 2:" << setw(8) << Form("%3.3f\n",PulseIntPMT2Fit->GetParameter(0));
+
+  //Start the process of writing the calibration information to file
+  ofstream calibration;
+  calibration.open("calibration_temp.txt", ios::out);
+
+  if (!calibration.is_open()) cout << "Problem saving calibration constants, may have to update constants manually!" << endl;
+
+  else {
+    calibration << "hcer_adc_to_npe = "; 
+    calibration << Form("1./%3.3f, 1./%3.3f", PulseIntPMT1Fit->GetParameter(0), PulseIntPMT2Fit->GetParameter(0));
+    calibration.close();
+  }
+}