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

+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.297 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.00  1.46 -1.42  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
+1764  0.659  1.080 -0.660  2.034 345.528 154.612 -0.99  1.17 -0.76  1.05  0.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 
+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.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/.rootrc

+Rint.History: .root_history

CALIBRATION/hms_cal_calib/THcShHit.h

+#include <iostream>
+
+// HMS calorimeter hit class for calibration.
+
+class THcShHit {
+
+  Double_t ADCpos, ADCneg;   // pedestal subtracted ADC signals.
+  Double_t Epos, Eneg;       // Energy depositions seen from pos. & neg. sides
+  UInt_t BlkNumber;
+
+ public:
+
+  THcShHit();
+  THcShHit(Double_t adc_pos, Double_t adc_neg,
+	   UInt_t blk_number);
+  ~THcShHit();
+
+  void SetADCpos(Double_t sig) {ADCpos = sig;}
+
+  void SetADCneg(Double_t sig) {ADCneg = sig;}
+
+  void SetEpos(Double_t e) {Epos = e;}
+
+  void SetEneg(Double_t e) {Eneg = e;}
+
+  void SetBlkNumber(UInt_t n) {BlkNumber = n;}
+
+  Double_t GetADCpos() {return ADCpos;}
+
+  Double_t GetADCneg() {return ADCneg;}
+
+  Double_t GetEpos() {return Epos;}
+
+  Double_t GetEneg() {return Eneg;}
+  
+  UInt_t GetBlkNumber() {return BlkNumber;}
+
+  void Print(ostream & ostrm);
+};
+
+//------------------------------------------------------------------------------
+
+THcShHit::THcShHit() {
+  ADCpos = -99999.;
+  ADCneg = -99999.;
+  Epos = -99999.;
+  Eneg = -99999.;
+  BlkNumber = 99999;
+};
+
+THcShHit::THcShHit(Double_t adc_pos, Double_t adc_neg,
+		   UInt_t blk_number) {
+  ADCpos = adc_pos;
+  ADCneg = adc_neg;
+  Epos = 0.;
+  Eneg = 0.;
+  BlkNumber = blk_number;
+};
+
+THcShHit::~THcShHit() { };
+
+//------------------------------------------------------------------------------
+
+void THcShHit::Print(ostream & ostrm) {
+
+  // Output hit data through the stream ostrm.
+
+  ostrm << ADCpos << " " << ADCneg << " " << Epos << " " << Eneg << " "
+	<< BlkNumber << endl;
+};
+
+struct pmt_hit {Double_t signal; UInt_t channel;};

CALIBRATION/hms_cal_calib/THcShTrack.h

+#include "THcShHit.h"
+#include "TMath.h"
+
+#include <vector>
+#include <iterator>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+// Track class for the HMS calorimeter calibration.
+// Comprises the spectrometer track parameters and calorimeter hits.
+//
+
+// Container (collection) of hits and its iterator.
+//
+typedef vector<THcShHit*> THcShHitList;
+typedef THcShHitList::iterator THcShHitIt;
+
+class THcShTrack {
+
+  Double_t P;   // track momentum
+  Double_t Dp;  // track momentum deviation, %
+  Double_t X;   // at the calorimater face
+  Double_t Xp;  // slope
+  Double_t Y;   // at the calorimater face
+  Double_t Yp;  // slope
+
+  THcShHitList Hits;
+
+ public:
+  THcShTrack();
+  THcShTrack(Double_t p, Double_t dp,
+	     Double_t x, Double_t xp, Double_t y, Double_t yp);
+  ~THcShTrack();
+
+  void Reset(Double_t p, Double_t dp,
+	     Double_t x, Double_t xp, Double_t y, Double_t yp);
+
+  void AddHit(Double_t adc_pos, Double_t adc_neg,
+	      Double_t e_pos, Double_t e_neg,
+	      UInt_t blk_number);
+
+  THcShHit* GetHit(UInt_t k);
+
+  UInt_t GetNhits() {return Hits.size();};
+
+  void Print(ostream & ostrm);
+
+  void SetEs(Double_t* alpha);
+  void SetEsNoCor(Double_t* alpha);
+
+  Double_t Enorm();
+  Double_t EPRnorm();
+  Double_t ETAnorm();
+
+  Double_t GetP() {return P*1000.;}      //MeV
+
+  Double_t GetDp() {return Dp;}
+
+  Double_t GetX() {return X;}
+  Double_t GetY() {return Y;}
+
+  Float_t Ycor(Double_t);         // coord. corection for single PMT module
+  Float_t Ycor(Double_t, Int_t);  // coord. correction for double PMT module
+
+  // Coordinate correction constants from hcana.param.
+  //
+  static constexpr Double_t fAcor = 200.;
+  static constexpr Double_t fBcor = 8000.;
+  static constexpr Double_t fCcor = 64.36;
+  static constexpr Double_t fDcor = 1.66;
+
+  // Calorimeter geometry constants.
+  //
+  static constexpr Double_t fZbl = 10;   //cm, block transverse size
+  static constexpr UInt_t fNrows = 13;
+  static constexpr UInt_t fNcols =  4;
+  static constexpr UInt_t fNnegs = 26;  // number of blocks with neg. side PMTs.
+  static constexpr UInt_t fNpmts = 78;  // total number of PMTs.
+  static constexpr UInt_t fNblks = fNrows*fNcols;
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShTrack::THcShTrack() { };
+
+THcShTrack::THcShTrack(Double_t p, Double_t dp,
+		       Double_t x, Double_t xp, Double_t y, Double_t yp) {
+  P = p;
+  Dp = dp;
+  X = x;
+  Xp = xp;
+  Y = y;
+  Yp =yp;
+};
+
+//------------------------------------------------------------------------------
+
+void THcShTrack::Reset(Double_t p, Double_t dp,
+		       Double_t x, Double_t xp, Double_t y, Double_t yp) {
+
+  // Reset track parameters, clear hit list.
+
+  P = p;
+  Dp = dp;
+  X = x;
+  Xp = xp;
+  Y = y;
+  Yp =yp;
+  Hits.clear();
+};
+
+//------------------------------------------------------------------------------
+
+void THcShTrack::AddHit(Double_t adc_pos, Double_t adc_neg,
+			Double_t e_pos, Double_t e_neg,
+			UInt_t blk_number) {
+
+  // Add a hit to the hit list.
+
+  THcShHit* hit = new THcShHit(adc_pos, adc_neg, blk_number);
+  hit->SetEpos(e_pos);
+  hit->SetEneg(e_neg);
+  Hits.push_back(hit);
+};
+
+//------------------------------------------------------------------------------
+
+THcShHit* THcShTrack::GetHit(UInt_t k) {
+  THcShHitIt it = Hits.begin();
+  for (UInt_t i=0; i<k; i++) it++;
+  return *it;
+}
+
+void THcShTrack::Print(ostream & ostrm) {
+
+  // Output the track parameters and hit list through the stream ostrm.
+
+  ostrm << P << " " << Dp << " " << X << " " << Xp << " " << Y << " " << Yp
+	<< " " << Hits.size() << endl;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    (*iter)->Print(ostrm);
+  };
+
+};
+
+//------------------------------------------------------------------------------
+
+THcShTrack::~THcShTrack() {
+  for (THcShHitIt i = Hits.begin(); i != Hits.end(); ++i) {
+    delete *i;
+    *i = 0;
+  }
+};
+
+//------------------------------------------------------------------------------
+
+void THcShTrack::SetEs(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) {
+      //      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.);
+    };
+
+  };
+
+}
+
+//------------------------------------------------------------------------------
+
+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.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    sum += (*iter)->GetEpos();
+    sum += (*iter)->GetEneg();
+  };
+
+  return sum/P/1000.;
+}
+
+//------------------------------------------------------------------------------
+
+Double_t THcShTrack::EPRnorm() {
+
+  // Normalized to the track momentum energy depostion in Preshower.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    UInt_t nblk = (*iter)->GetBlkNumber();
+    Int_t ncol=(nblk-1)/fNrows+1;
+    if (ncol==1) {
+      sum += (*iter)->GetEpos();
+      sum += (*iter)->GetEneg();
+    }
+  }
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+//------------------------------------------------------------------------------
+
+Double_t THcShTrack::ETAnorm() {
+
+  // Normalized to the track momentum energy depostion in Preshower.
+
+  Double_t sum = 0;
+
+  for (THcShHitIt iter = Hits.begin(); iter != Hits.end(); iter++) {
+    UInt_t nblk = (*iter)->GetBlkNumber();
+    Int_t ncol=(nblk-1)/fNrows+1;
+    if (ncol!=1) {
+      sum += (*iter)->GetEpos();
+      sum += (*iter)->GetEneg();
+    }
+  }
+
+  return sum/P/1000.;         //Momentum in MeV.
+}
+
+//------------------------------------------------------------------------------
+
+//Coordinate correction for single PMT modules.
+//PMT attached at right (positive) side.
+
+Float_t THcShTrack::Ycor(Double_t y) {
+  return TMath::Exp(y/fAcor)/(1. + y*y/fBcor);
+}
+
+//Coordinate correction for double PMT modules.
+//
+
+Float_t THcShTrack::Ycor(Double_t y, Int_t side) {
+  if (side!=0&&side!=1) {
+    cout << "THcShower::Ycor : wrong side " << side << endl;
+    return 0.;
+  }
+  Int_t sign = 1 - 2*side;
+  return (fCcor + sign*y)/(fCcor + sign*y/fDcor);
+}

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 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 file " << Prefix << ".root, events "
+       << nstart << " -- " << nstop << endl;
+
+  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.
+ theShowerCalib.SolveAlphas();     // Solve for the calibration constants
+ 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(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");
+
+ // Uncalibrated e/p at calorimeter
+ Canvas->cd(3);
+ //  theShowerCalib.hCaloPosNormU->GetZaxis()->SetRangeUser(0.5,1.5);
+ theShowerCalib.hCaloPosNormU->Draw("COLZ");
+ theShowerCalib.hCaloPosNormU->SetStats(0);
+ theShowerCalib.hCaloPosNormU->GetZaxis()->SetRangeUser(0.5,1.5);
+
+ // 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(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 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 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+("Prefix"),
+
+where Prefix is the prefix part in the name of root file.
+
+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.<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.