diff --git a/CALIBRATION/shms_dc_calib/README.md b/CALIBRATION/shms_dc_calib/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e16f7ed1345a6cc624e66c24770da5a82888b4bd
--- /dev/null
+++ b/CALIBRATION/shms_dc_calib/README.md
@@ -0,0 +1,88 @@
+SHMS Drift Chambers Calibration
+============================================
+This directory contains the code for calibrating the pair of SHMS drift chambers.
+
+
+
+Directory structure
+----------------------
+* hallc_replay/CALIBRATION/shms_dc_calib/run_Cal.C : steering C++ code that executes all codes in the 'scripts' directory
+* hallc_replay/CALIBRATION/shms_dc_calib/scripts : all scripts necessary to do calibration lie in this directory
+* hallc_replay/CALIBRATION/shms_dc_calib/root_files : all root files produced by the calibration are in this directory
+* hallc_replay/CALIBRATION/shms_dc_calib/data_files : all data files produced by the calibration are in this directory
+
+
+
+
+Running code
+---------------
+* First set the parameter 'p_using_tzero_per_wire = 0' in the
+ parameter file located at: hallc_replay/PARAM/SHMS/DC/pdc.param
+
+* Replay the data to produce the uncalibrated root file to be used as input in the calibration
+ * From the hallc_replay execute: ./hcana SCRIPTS/SHMS/replay_shms.C
+
+* From THIS! directory execute: root -l run_Cal.C
+
+* From the calibration results, two parameter files will be produced in:
+ * hallc_replay/PARAM/SHMS/DC/pdc_tzero_per_wire_run%d_NEW.param. %d=run_number
+ * hallc_replay/PARAM/SHMS/DC/pdriftmap_new.param
+
+* Rename the new parameter files as follows:
+ * copy: pdc_tzero_per_wire_run%d_NEW.param to pdc_tzero_per_wire.param
+ * copy: pdriftmap_new.param to pdriftmap.param
+
+* Before replaying the data again, set the parameter 'p_using_tzero_per_wire = 1' to
+ allow the source code (hcana) to read the parameter values during the replay.
+
+* Replay the data with the updated parameters to produce the new calibrated root files
+ with the corrected drift times and drift distances.
+
+* Compare the calibrated and uncalibrated root files to verify the calibration was done properly.
+
+
+
+Brief decription of code
+------------------------
+* Overview: The code determines the tzero offsets on a wire-by-wire basis. These offsets are
+ corrections by which each wire drift time spectrum must be shifted to align it with a
+ drift time of "0 ns". This time corresponds to when the electrons from the ionized gas inside
+ the chamber are in contact with the sense wire, hence a drift time of "0 ns".
+
+* Brief Description of the scripts in the '/scripts' directory
+ *** to run the scripts independenlty, open the input_RUN.txt files and write the filen_ame, run_number and number of events in this order: root_file.root run_number events
+
+ ** get_pdc_time_histo.C :
+ -- outputs root_file: 'shms_dc_time_%d.root', (%d=run_number)
+ -- contains re-binned per-plane drift time histograms
+
+ ** wire_drift_times.C (and wire_drift_times.h) - See instructions in wire_drift_times.C on how to execute independently
+ -- outputs root_file: 'shms_DC_plane_%d_wire_histos.root'
+ -- contains a 2-D histo of "drift time vs. Wire Number" and drift time spectra for all wires in the plane
+
+ ** get_wire_tzero.C
+ -- outputs root_file: 'shms_DC_plane_%d_fitted_histos.root'
+ -- contains line-fitted wire drift time histos. The extrapolation of the fit to the x-axis is defined as "tzero"
+ Also contains two "tzero vs. wirenumber" plots: one for all wires, and the other is after setting tzero=0 for
+ wires that did not have enough statistics for a good quality fit, a weighted average was calculated for the latter.
+
+ -- outputs data_file: 'shms_dc_planetzero_run%d.dat' -> contains list of "tzero" values for all wires
+ 'shms_dc_planetzero_run%d_updated.dat' -> contains list of "tzero" values for all wires that
+ had a good fit. The remaining tzeros are set to 0.
+ 'tzero_weighted_avg_run%d.dat -> contains weighted tzero values per plane
+
+ ** get_tzero_per_wire_param.C
+ -- outputs data_file: 'tzero_values_per_wire.dat'
+ -- contains list of tzero values for all wires in all planes
+
+ -- outputs param_file: /hallc_replay/PARAM/SHMS/DC/pdc_tzero_per_wire_run%d_NEW.param, where %d=run_number
+ -- contains tzero values for all wires in all planes, but the file is formatted so that the values may be read by the source code (hcana)
+
+ ** get_pdc_time_histo_tzero_corrected.C
+ -- outputs root_file: 'shms_tzero_corr_histos.root'
+ -- contains list of "t0-corrected" wire drift times, and their respective plane drift times.
+
+ ** get_LookUp_Values.C
+ -- outputs param_file: /hallc_replay/PARAM/SHMS/DC/pdriftmap_new.param
+ -- contains scaling factors calculated from the corrected plane drift times on a bin-by-in basis. These values get read by the source code
+ which will be used to scale the drift distance histograms.
diff --git a/CALIBRATION/shms_dc_calib/run_Cal.C b/CALIBRATION/shms_dc_calib/run_Cal.C
index 806e42f2c66e3a69186615c5d2143f0b6035e7dc..2f6377006371b4e766b27c5c07cce2699a569053 100644
--- a/CALIBRATION/shms_dc_calib/run_Cal.C
+++ b/CALIBRATION/shms_dc_calib/run_Cal.C
@@ -52,11 +52,6 @@
system(dir_data);
}
- if (std::system(dir_data) != 0)
- {
-
- }
-
//Create run Directories if they dont exist
diff --git a/CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo_tzero_corrected.C b/CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo_tzero_corrected.C
index e2ec3d519a90fbe023b16a229e3c95867922c882..5a24ddfbda00358f78ea2c2b47c47b5b82d10d2c 100644
--- a/CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo_tzero_corrected.C
+++ b/CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo_tzero_corrected.C
@@ -25,7 +25,7 @@ void get_pdc_time_histo_tzero_corrected()
int tot_wires[NPLANES] = {107, 107, 79, 79, 107, 107, 107, 107, 79, 79, 107, 107};
//sum over all wires in both DC
- static int wire_sum = 0;
+ int wire_sum = 0;
for (ip=0; ip<NPLANES; ip++) {
wire_sum = wire_sum + tot_wires[ip];
}
@@ -38,8 +38,9 @@ void get_pdc_time_histo_tzero_corrected()
string line;
int counter;
Double_t value;
- Double_t tzero_offset[wire_sum];
-
+ Double_t *tzero_offset;
+ tzero_offset = new Double_t[wire_sum];
+
counter = 0;
@@ -116,29 +117,31 @@ void get_pdc_time_histo_tzero_corrected()
TString file_name = Form("../root_files/run%d/shms_tzero_corr_histos.root", run_NUM);
TFile *g = new TFile(file_name, "RECREATE");
-
TH1F *h_add[NPLANES]; //t0-corrected plane drift times
-
-
+
+
+
//Loop over all planes
for (ip = 0; ip < NPLANES; ip++){
-
- //READ wire drift time root file per plane
- root_file = "../root_files/"+run+"/shms_DC_"+plane_names[ip]+Form("_%d_wire_histos.root",run_NUM);
- f[ip] = new TFile(root_file, "READ");
+ //READ wire drift time root file per plane
+ root_file = "../root_files/"+run+"/shms_DC_"+plane_names[ip]+Form("_%d_wire_histos.root",run_NUM);
+ f[ip] = new TFile(root_file, "READ");
+
h_add[ip] =new TH1F("plane_"+plane_names[ip]+"drifttime", "", nbins, -50, 350);
h_add[ip]->GetXaxis()->SetTitle("Drift Time (ns)");
h_add[ip]->GetYaxis()->SetTitle("Number of Entries / 2 ns");
- TH1F *cell_dt[tot_wires[ip]];
- TH1F *cell_dt_corr[tot_wires[ip]];
+
+ TH1F *cell_dt[107];
+ TH1F *cell_dt_corr[107];
-
+
+
//Get wire histos from root file and loop over each
// sense wire of a plane in shms Drift Chambers (DC1 or DC2)
@@ -146,8 +149,8 @@ void get_pdc_time_histo_tzero_corrected()
for (sw=1; sw<=tot_wires[ip]; sw++){
-
-
+
+
//set title of histos in root file
TString drift_time_histo = Form("wire_%d", sw);
@@ -169,7 +172,7 @@ void get_pdc_time_histo_tzero_corrected()
//************APPLY TZERO OFFSET ON A WIRE-BY-WIRE BASIS TO ALL WIRES IN ALL PLANES***************
- //INCLUDE the code 'shift.C ', which shifts a histogram
+ //INCLUDE the code 'shift.C ', which shifts a histogram
for (i=1; i<=nbins; i++) {
@@ -185,7 +188,7 @@ void get_pdc_time_histo_tzero_corrected()
//*************************************************************************************************
-
+
//write wire drift times (after tzero corrections) to file
g->cd();
@@ -195,21 +198,22 @@ void get_pdc_time_histo_tzero_corrected()
//add all cell drift times into a single plane
h_add[ip]->Add(cell_dt_corr[sw-1]);
-
-
+
} // end loop over sense wires
//Wire combined wire drift times (t0 -corrected) for each plane to file,
- g->cd();
- h_add[ip]->Write();
+ g->cd();
+ h_add[ip]->Write();
+
+
} // end loop over planes
-
+
}
diff --git a/CALIBRATION/shms_dc_calib/scripts/get_tzero_per_wire_param.C b/CALIBRATION/shms_dc_calib/scripts/get_tzero_per_wire_param.C
index 9bb118faff08bba70032ad4bd2cb7cf5cb4858b1..b6bb687d096c67bcb7ef76cc945d5e4cdf9f09f6 100644
--- a/CALIBRATION/shms_dc_calib/scripts/get_tzero_per_wire_param.C
+++ b/CALIBRATION/shms_dc_calib/scripts/get_tzero_per_wire_param.C
@@ -5,11 +5,11 @@
#define NPLANES 12
-
+using namespace std;
void get_tzero_per_wire_param()
{
-
+
//read run number from input file
int run_NUM;
Long64_t num_evts; //added
@@ -20,15 +20,14 @@ void get_tzero_per_wire_param()
TString run = Form("run%d", run_NUM);
-
TString planes[NPLANES] = {"1u1","1u2","1x1","1x2","1v1","1v2","2v2","2v1","2x2","2x1","2u2","2u1"}; //planes 0, 1, 2, 3, 4, 5 OK ::
-int fNWires[NPLANES] = {107, 107, 79, 79, 107, 107, 107, 107, 79, 79, 107, 107};
-
-int sw;
-int ip;
+ int fNWires[NPLANES] = {107, 107, 79, 79, 107, 107, 107, 107, 79, 79, 107, 107};
+ int sw;
+ int ip;
+
Int_t **wire = new Int_t*[NPLANES];
Double_t **t0 = new Double_t*[NPLANES];
//Double_t **t0_err = new Double_t*[NPLANES];
@@ -41,47 +40,55 @@ Double_t **t0 = new Double_t*[NPLANES];
}
}
-
-
+
//Loop over each plane
for(ip=0; ip<NPLANES; ip++) {
//write plane headers
- Int_t nwire;
- Double_t nt0;
- Double_t nt0_err;
- Int_t nentries;
-
+ int nwire;
+ double nt0;
+ double nt0_err;
+ int nentries;
//open and read each wire tzero file
- string line;
- ifstream input;
- input.open(Form("../data_files/run%d/shms_dc_"+planes[ip]+"tzero_run%d_updated.dat", run_NUM, run_NUM ) );
-
-
- sw = 0; //se wire counter to 0
+ TString file_name = Form("../data_files/run%d/shms_dc_", run_NUM) + planes[ip] + Form("tzero_run%d_updated.dat", run_NUM);
+
+ // cout << file_name << endl;
+
+ string line;
+ ifstream input;
+ input.open(file_name);
+
+ sw = 0;
while (getline(input, line))
{
-
- if (line != '#')
- {
- input >> nwire >> nt0 >> nt0_err >> nentries;
- // cout << nwire << " :: " << nt0 << endl;
- //cout << ip << " :: " << sw << endl;
-
- t0[ip][nwire-1] = nt0;
+ // cout << line << endl;
+ if (line[0]!='#')
+ {
+ if(!input.good()) {break;}
+ sscanf(line.c_str(), "%d %lf", &nwire, &nt0);
+
+ if(sw>0) {
- sw++;
+ t0[ip][nwire-1] = nt0;
+
+ }
+ sw++;
+
}
}
input.close();
+
+
+
+
} //end loop over planes
diff --git a/CALIBRATION/shms_dc_calib/scripts/get_wire_tzero.C b/CALIBRATION/shms_dc_calib/scripts/get_wire_tzero.C
index 4c301360ea502493d8d801a45c468a9345f6ff38..b0031a4b7a3340e45ade08032d0501cac68f0e99 100644
--- a/CALIBRATION/shms_dc_calib/scripts/get_wire_tzero.C
+++ b/CALIBRATION/shms_dc_calib/scripts/get_wire_tzero.C
@@ -6,6 +6,8 @@ a certain number of bins and this fit is extrapolated to y=0(x-axis). The extrap
#include <vector>
#include "TMath.h"
+#include <iomanip>
+#include <iostream>
#define NPLANES 12
using namespace std;
@@ -53,7 +55,7 @@ void get_wire_tzero()
ofs.open (t_zero_file);
//Set headers for subsequent columns of data
- ofs << "#WIRE " << " " << "t0" << " " << "t0_err" << " " << " entries " << endl;
+ ofs << "#WIRE " << setw(15) << "t0" << setw(15) << "t0_err" << setw(15) << "entries" << endl;
//Create root file to store fitted wire drift times histos and "t0 vs. wirenum"
TString output_root_file = "../root_files/"+run+"/shms_DC_"+plane_names[ip]+Form("_%d_fitted_histos.root", run_NUM);
@@ -239,7 +241,7 @@ void get_wire_tzero()
entries = cell_dt[sensewire-1]->GetEntries(); //number of entries (triggers) per wire
//cout << "y_int: " << y_int << " :: " << "m: " << m << " :: " << "t0: " << setprecision(6) << -y_int/m << endl;
//Write "t0" values to file
- ofs << sensewire << " " << setprecision(6) << -y_int/m << " " << t_zero_err << " " << entries << endl;
+ ofs << sensewire << setw(15) << setprecision(6) << -y_int/m << setw(15) << t_zero_err << setw(15) << entries << endl;
//Change to output root file and write fitted histos to file
g->cd();
@@ -268,13 +270,10 @@ void get_wire_tzero()
graph->Draw("AP");
t->Update();
t->Write(title); //write to a root file
+
//close dat file
ofs.close();
- //save plots
- //TString tzero_plots = "plots/"+run_NUM +"/hdc"+plane_names[ip]+Form("TESTING_tzero_v_wire_%d.eps", run);
- //t->SaveAs(tzero_plots);
-
//*****************************************************************************************//
// CALCULATE THE "t0s" WEIGHTED AVERAGE FOR WIRE DRIFT TIMES WITH ENTRIES > = 300 //
@@ -289,7 +288,7 @@ void get_wire_tzero()
//open new data file to write updated t0 values
TString t_zero_file_corr = "../data_files/" + run + "/shms_dc_"+plane_names[ip]+Form("tzero_run%d_updated.dat", run_NUM);
ofs.open(t_zero_file_corr);
- ofs << " #Wire " << " " << " t_zero " << " " << " t_zero_err " << " " << " entries " << endl;
+ ofs << "#Wire" << setw(15) << "t_zero" << setw(15) << "t_zero_err" << setw(15) << "entries" << endl;
//Initialize variables related to weighted avg
double sum_NUM; //numerator of weighted avg
@@ -298,9 +297,12 @@ void get_wire_tzero()
double weighted_AVG_err;
int counter;
- double t0_corr;
- double t0_corr_err;
- //set them to zero to start sum inside while loop
+ // double t0_corr;
+ //double t0_corr_err;
+ Double_t t0;
+ Double_t t0_err;
+
+ //set them to zero to start sum inside while loop
sum_NUM = 0.0;
sum_DEN = 0.0;
@@ -311,38 +313,40 @@ void get_wire_tzero()
//read line bt line the t_zero_file
while(getline(ifs, line)) {
- if(line!='#') {
+ if(line[0]!='#') {
- sscanf(line.c_str(), "%d %lf %lf %d", &sensewire, &t0_corr, &t0_corr_err, &entries); //assign each of the variables above a data in the t_zero_file
-
- if(sensewire<=fNWires[ip]){
-
+ sscanf(line.c_str(), "%d %lf %lf %d", &sensewire, &t0, &t0_err, &entries); //assign each of the variables above a data in the t_zero_file
+ // if(sensewire<=fNWires[ip]){
+
+ cout << sensewire << endl;
//Check if entries for each sensewire exceeds a certain number of events
- if (entries>300 && t_zero < 30)
-{
+ if (entries>300 && abs(t0)<30)
+ {
-
-
- //Calculate the weighted average of t0s
- sum_NUM = sum_NUM + t0_corr/(t0_corr_err*t0_corr_err);
- sum_DEN = sum_DEN + 1.0/(t0_corr_err*t0_corr_err);
-
- //cout << "sum_NUM : " << sum_NUM << endl;
- //cout << "sum_DEN : " << sum_DEN << endl;
-
-
+ //cout << ip << "::" << sensewire << "::" << t0 << endl;
+ //Calculate the weighted average of t0s
+ sum_NUM = sum_NUM + t0/(t0_err*t0_err);
+ sum_DEN = sum_DEN + 1.0/(t0_err*t0_err);
+
+ //cout << "sum_NUM : " << sum_NUM << endl;
+ //cout << "sum_DEN : " << sum_DEN << endl;
+
+
+
+ ofs << sensewire << setw(15) << t0 << setw(15) << t0_err << setw(15) << entries << endl;
+ //cout << "TZERO: " << t0 << endl;
+
+
- ofs << sensewire << " " << t0_corr << " " << t0_corr_err << " " << entries << endl;
- //cout << "TZERO: " << t_zero << endl;
-
-
- }
-
- else { ofs << sensewire << " " << 0.0 << " " << 0.0 << " " << entries << endl;}
}
+
+ else { ofs << sensewire << setw(15) << 0.0 << setw(15) << 0.0 << setw(15) << entries << endl;}
+
+ //}
+
} //end if statement
}
@@ -362,10 +366,7 @@ void get_wire_tzero()
ofile << " #weighted_AVG " << " " << " DC plane: " << plane_names[ip] << endl;
ofile << weighted_AVG << endl;
-
-
-
-
+
ifs.close();
@@ -402,9 +403,9 @@ void get_wire_tzero()
ltx1->DrawLatex(t1->GetUxmax()*0.75,40, Form("Weighted Average = %lf #pm %lf ns", weighted_AVG, weighted_AVG_err) );
t1->Write(title1); //write canvas to a root file
- ofs.close(); //close data file
-
-
+
+ ofs.close(); //close data file
+
diff --git a/CALIBRATION/shms_dc_calib/scripts/wire_drift_times.C b/CALIBRATION/shms_dc_calib/scripts/wire_drift_times.C
index 7cf84cb779bbee1976f6b736c38ebc2c3c589c1b..cb0a6bc1d22f510080428cadf9e4a1013dd02ed9 100644
--- a/CALIBRATION/shms_dc_calib/scripts/wire_drift_times.C
+++ b/CALIBRATION/shms_dc_calib/scripts/wire_drift_times.C
@@ -87,7 +87,7 @@ void wire_drift_times::Loop()
for (wirenum=1; wirenum<=total_wires_uv; wirenum++){
cell_dt[wirenum-1] = new TH1F(Form("wire_%d", wirenum), "", 200., -50., 350.);
- cell_dt[wirenum-1]->GetXaxis()->SetTitle("Wire Number");
+ cell_dt[wirenum-1]->GetXaxis()->SetTitle("Drift Time (ns)");
cell_dt[wirenum-1]->GetYaxis()->SetTitle("Number of Entries / 2 ns");
}
@@ -221,7 +221,7 @@ void wire_drift_times::Loop()
for (int wirenum=1; wirenum<=total_wires_x; wirenum++){
cell_dt[wirenum-1] = new TH1F(Form("wire_%d", wirenum), "", 200., -50., 350.);
- cell_dt[wirenum-1]->GetXaxis()->SetTitle("Wire Number");
+ cell_dt[wirenum-1]->GetXaxis()->SetTitle("Drift Time (ns)");
cell_dt[wirenum-1]->GetYaxis()->SetTitle("Number of Entries / 2 ns");