CALIBRATION/hms_cer_calib/calibration.h

+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Wed Aug 15 17:10:26 2018 by ROOT version 6.10/02
+// from TTree T/Hall A Analyzer Output DST
+// found on file: ../../ROOTfiles/coin_replay_production_3424_-1.root
+//////////////////////////////////////////////////////////
+
+#ifndef calibration_h
+#define calibration_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+#include <TTreeReader.h>
+#include <TTreeReaderValue.h>
+#include <TTreeReaderArray.h>
+
+// Headers needed by this particular selector
+
+
+class calibration : public TSelector {
+public :
+   TTreeReader     fReader;  //!the tree reader
+   TTree          *fChain = 0;   //!pointer to the analyzed TTree or TChain
+   Bool_t          fFullShow;
+   
+
+   //Declare histograms
+   TH1F          **fPulseInt; 
+   TH1F         ***fPulseInt_quad;
+   TH1F           *fCut_everything;
+   TH1F           *fCut_electron;
+   TH1F           *fBeta_Cut;
+   TH1F           *fBeta_Full;
+   TH1F          **fTiming_Cut;
+   TH1F          **fTiming_Full;
+
+   // Readers to access the data (delete the ones you do not need).
+   TTreeReaderArray<Double_t> H_cer_goodAdcPulseInt    = {fReader, "H.cer.goodAdcPulseInt"};
+   TTreeReaderArray<Double_t> H_cer_goodAdcTdcDiffTime = {fReader, "H.cer.goodAdcTdcDiffTime"};
+   TTreeReaderValue<Double_t> H_cer_xAtCer             = {fReader, "H.cer.xAtCer"};
+   TTreeReaderValue<Double_t> H_cer_yAtCer             = {fReader, "H.cer.yAtCer"};
+
+   TTreeReaderArray<Double_t> H_tr_beta                = {fReader, "H.tr.beta"};
+   TTreeReaderValue<Int_t> Ndata_H_tr_beta             = {fReader, "Ndata.H.tr.beta"};
+
+   TTreeReaderValue<Double_t> H_cal_etotnorm           = {fReader, "H.cal.etotnorm"};
+
+
+   calibration(TTree * /*tree*/ =0) {fPulseInt=0,fPulseInt_quad=0,fCut_everything=0,fCut_electron=0,fBeta_Cut=0,fBeta_Full=0,fTiming_Cut=0,fTiming_Full=0,fFullShow=kFALSE;}
+   virtual ~calibration() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(calibration,0);
+
+};
+
+#endif
+
+#ifdef calibration_cxx
+void calibration::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the reader is initialized.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   fReader.SetTree(tree);
+}
+
+Bool_t calibration::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+
+#endif // #ifdef calibration_cxx

CALIBRATION/hms_cer_calib/run_cal.C

+#include <TProof.h>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <stdio.h>
+
+void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0, Int_t coin = 0)
+{
+  if (RunNumber == 0) {
+    cout << "Enter a Run Number (-1 to exit): ";
+    cin >> RunNumber;
+    if (RunNumber <= 0) return;
+  }
+  if (NumEvents == 0) {
+    cout << "\nNumber of Events to analyze: ";
+    cin >> NumEvents;
+  }
+  if (coin == 0) {
+    cout << "\nIf this is a coincident run enter 1: ";
+    cin >> coin;
+  }
+
+  cin.ignore(numeric_limits<streamsize>::max(), '\n');
+
+  string calib_raw;
+  cout << "\nEnter options for calibration  (enter NA to skip): ";
+  getline(std::cin, calib_raw);
+  TString calib_option = calib_raw;
+
+  cout << "\n\n";
+
+  TChain ch("T");
+  if (coin == 1) ch.Add(Form("../../ROOTfiles/coin_replay_production_%d_%d.root", RunNumber, NumEvents));
+  else ch.Add(Form("../../ROOTfiles/hms_replay_production_all_%d_%d.root", RunNumber, NumEvents));
+  TProof *proof = TProof::Open("");
+  //proof->SetProgressDialog(0);  
+  ch.SetProof();
+
+  if (calib_option != "NA") {
+    //Start calibration process
+    ch.Process("calibration.C+",calib_option);
+
+    cout << "\n\nUpdate calibration constants with the better estimate (y/n)? ";
+      
+    TString user_input;
+    cin >> user_input;
+    if (user_input == "y") {
+      ifstream temp;
+      temp.open("calibration_temp.txt", ios::in);
+      if (temp.is_open()) {
+	rename("calibration_temp.txt", Form("../../PARAM/HMS/CER/hcer_calib_%d.param", RunNumber));
+      }
+
+      else cout << "Error opening calibration constants, may have to update constants manually!" << endl;
+	   
+    }
+
+    else {
+      remove("calibration_temp.txt"); 
+    }
+  }
+}

CALIBRATION/hms_hodo_calib/.gitignore

+*.pdf
+*.root
+*.so
+*.pcm

CALIBRATION/hms_hodo_calib/.rootrc

+Rint.History: .root_history

CALIBRATION/hms_hodo_calib/README.md

+# Hall C HMS Hodo Calibration Codes
+
+A set of ROOT scripts to determine the time offset parameters for the HMS hodoscope. To use these parameters a  parameter flag ptofusinginvadc is set to 0. The first script determines the time walk correction factor. The next script determines the effective propagation speeed in the paddle, the time difference between the positive and negative PMTs and then the relative time difference of all paddles compared to paddle 7 in plane S1X.
+
+
+This supercedes an older calibration code. The parameters from that code are
+used if a parameter flag ptofusinginvadc is set to 1. 
+
+The two codes have different parameters and it is possible to switch between the parameters of the two calibration codes using the  parameter flag ptofusinginvadc
+
+
+## Instructions
+
+1. Replay the data with ptofusinginvadc=0 and need to have T.hms.* and H.hod.* in the tree. 
+   Go to: PARAM/HMS/GEN/h_fadc_debug.param, and set hhodo_debug_adc = 1 --->  Use the correct hodo raw leafs variables 
+
+(In case you are calibrating HMS hodo using a coincidence run, then make sure to include T.coin.*)
+
+
+2. Determine the time walk correction parameters
+
+     a. Start "root -l" and then  .x timeWalkHistos.C+("current_dir/to/ROOT_filename.root", Run_Number, "hms") ---> If doing coincidence, then "hms"->"coin"
+
+     b. This creats the file: timeWalkHistos.root
+
+     c. Start "root -l" and then .x timeWalkCalib.C+
+
+     d. This creates the parameter file "../../PARAM/HMS/HODO/hhodo_TWcalib_runnumber.param"
+
+3.  Replay the data with ptofusinginvadc=0 and the new parameter files (the simplest is to copy phodo_TWcalib_runnumber.param to phodo_TWcalib.param).
+
+4. Determine the effective propagation speed in the paddle, the time difference between the positive and negative PMTs and then the relative time difference of all paddles compared to paddle 7 in plane S1X. The script
+puts cuts on H.cal.etracknorm, H.hgcer.npeSum and H.hod.betanotrack to select electrons. These cuts are hard coded as  etrknrm_low_cut = 0.7, npngcer_npeSum_low_cut = 0.7 , betanotrack_low_cut = 0.5 and betanotrack_hi_cut = 1.5. These may need to be modified. The event must have a track. 
+
+*** Make sure to have included H.hgcer.* and H.cal.etracknorm in the ROOT tree
+
+     a.  Start "root -l" and then  .x  fitHodoCalib.C+("current_dir/to/ROOT_filename.root",Run_Number)     
+
+     b.  This creates the parameter file "../../PARAM/HMS/HODO/hhodo_Vpcalib_runnumber.param"
+
+     c.  It also creates the root file HodoCalibPlots_runnumber.root
+
+     d.  To analyze cosmic data :  .x  fitHodoCalib.C+("current_dir/to/ROOT_filename.root",Run_Number,kTRUE) 
+
+     e. For cosmic data the speed of light is set to -30 cm/ns and the PID cut is just on P.hod.betanotrack with the default of betanotrack_low_cut = -1.2 and betanotrack_hi_cut = -.7

CALIBRATION/set_peddefault/README.md

+1) First need to prepare histograms of the goodAdcPed for each detector.
+2) the 2d histograms of goodADCPed verus the paddle number are contained
+   in HMS_PedDefault.def and SHMS_PedDefault.def
+   Move files to DEF-files/SHMS/PRODUCTION/SHMS_PedDefault.def
+   and DEF-files/HMS/PRODUCTION/HMS_PedDefault.def
+3) Include the files in your  tree/histogram def file set in replay script.
+   #include "DEF-files/SHMS/PRODUCTION/SHMS_PedDefault.def"
+   #include "DEF-files/HMS/PRODUCTION/HMS_PedDefault.def"
+4) Replay the data
+5) Start root
+   a) .L run_ped_default.C
+   b) run_shms_ped_default("entirepath/DirName/filename.root")
+   c) run_hms_ped_default("entirepath/DirName/filename.root")
+6) The SHMS does HGCER, NGCER, AERO, Preshower and Shower
+7) HSM does CER and CAL.
+8) By hand copy each set of Pedestal defaults into the detector "cuts" file
+  For example:
+phgcer_PedDefault= 2086, 2153, 2320, 1987
+
+  goes into phgcer_cuts.param or the equivelent
+

CALIBRATION/set_peddefault/run_ped_default.C

+#include "TFile.h"
+#include "TH1D.h"
+#include <iostream>
+void calc_ped_default(TString , TString ,TString, Double_t);
+void run_shms_ped_default(TString);
+void run_hms_ped_default(TString);
+
+//
+void run_shms_ped_default(TString file_name = "") {
+  calc_ped_default(file_name,"hgcer","p",0);
+  calc_ped_default(file_name,"ngcer","p",0);
+  calc_ped_default(file_name,"aero","p",1);
+  calc_ped_default(file_name,"aero","p",2);
+  calc_ped_default(file_name,"cal_prshwr","p",1);
+  calc_ped_default(file_name,"cal_prshwr","p",2);
+  calc_ped_default(file_name,"cal_shwr","p",0);
+}
+//
+void run_hms_ped_default(TString file_name = "") {
+  calc_ped_default(file_name,"cer","h",0);
+  calc_ped_default(file_name,"cal_hA","h",1);
+  calc_ped_default(file_name,"cal_hB","h",1);
+  calc_ped_default(file_name,"cal_hC","h",1);
+  calc_ped_default(file_name,"cal_hD","h",1);
+  calc_ped_default(file_name,"cal_hA","h",2);
+  calc_ped_default(file_name,"cal_hB","h",2);
+
+}
+
+
+void calc_ped_default(TString golden_file = "", TString detector = "",
+                           TString spect = "", Double_t polarity = -1) {
+
+  if (golden_file == "") {
+    cout << "Enter golden run root file name: " << endl;
+    cin >> golden_file;
+  }
+  if (detector == "") {
+    cout << "Enter detector prefix (hodo_1x (etc.), hgcer, aero, cal_prshwr, "
+            "cal_shwr, ngcer, cal_hA (etc.)): "
+         << endl;
+    cin >> detector;
+  }
+  if (spect == "") {
+    cout << "Enter a spectrometer (p or h): " << endl;
+    cin >> spect;
+  }
+  if (polarity == -1) {
+    cout << "Enter a detector polarity (pos = 1, neg = 2): " << endl;
+    cin >> polarity;
+  }
+
+  TString histname = Form("%s%s", spect.Data(), detector.Data());
+  if (histname.Contains("hodo") && histname.Contains("1x") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("1x") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("1y") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("1y") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("2x") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("2x") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("2y") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("2y") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("aero") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("aero") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("hgcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("ngcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("_prshwr") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("_prshwr") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("_shwr"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("hcal_h") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hcal_h") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+
+  TH2F* H1_ped_vs_pmt;
+
+  TFile* f1 = new TFile(golden_file, "READ");
+  if (f1->IsZombie()) {
+    cout << "Cannot find : " << golden_file << endl;
+    return;
+  }
+
+  f1->GetObject(histname.Data(), H1_ped_vs_pmt);
+
+  TH1D* H1_pmt;
+
+  H1_pmt = (TH1D*)H1_ped_vs_pmt->ProjectionX("H1_pmt", 1,
+                                             H1_ped_vs_pmt->GetSize() - 2);
+
+  TH1D* H1_ped[H1_pmt->GetSize() - 2];
+
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+    H1_ped[ipmt] = (TH1D*)H1_ped_vs_pmt->ProjectionY(
+        Form("H1_ped_pmt%d", ipmt + 1), ipmt + 1, ipmt + 1);
+  }
+
+  Double_t H1_ped_peak[H1_pmt->GetSize() - 2];
+  Double_t H1_ped_fitpeak[H1_pmt->GetSize() - 2];
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+
+    if (H1_ped[ipmt]->GetEntries() > 25) {
+      TSpectrum* s = new TSpectrum(1);
+      gSystem->RedirectOutput("/dev/null", "a");
+      s->Search(H1_ped[ipmt], 1.0, "nobackground&&nodraw", 0.001);
+      gSystem->RedirectOutput(0);
+      TList* functions = H1_ped[ipmt]->GetListOfFunctions();
+      TPolyMarker* pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
+      if (pm) {
+        H1_ped_peak[ipmt] = *pm->GetX();
+      } else {
+        H1_ped_peak[ipmt] = 0;
+      }
+    } else {
+      H1_ped_peak[ipmt] = 0;
+    }
+
+  }
+
+  TF1* Gaussian =
+      new TF1("Gaussian", "[0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2]))");
+  Gaussian->SetParLimits(0, 0, 1000);
+  Gaussian->SetParLimits(2, 0, 2);
+
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+    if (H1_ped_peak[ipmt] != 0) {
+    Gaussian->SetParameter(1, H1_ped_peak[ipmt]);
+    gSystem->RedirectOutput("/dev/null", "a");
+    H1_ped[ipmt]->Fit(Gaussian, "QMN");
+    gSystem->RedirectOutput(0);
+    H1_ped_peak[ipmt]=int(Gaussian->GetParameter(1)*4*4.096);
+    }
+  }
+  TString side="";
+  if (polarity ==1)  side="Pos";
+  if (polarity ==2)  side="Neg";
+  TString PedDefName = spect+detector+"_Ped"+side+"Default= ";
+  if (detector == "cal_prshwr") {
+    PedDefName = spect+"cal_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_shwr") {
+    PedDefName = spect+"cal_arr_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_hA") {
+    PedDefName = spect+"cal_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_hB" || detector == "cal_hC" || detector == "cal_hD") {
+  } else {
+  cout << PedDefName ;
+  }
+   for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << H1_ped_peak[ipmt] ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+   } 
+   cout << endl;
+   if (detector == "cal_hB" && polarity==2) {
+    for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << 0 ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+     } 
+   cout << endl;
+     for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << 0 ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+     } 
+    cout << endl;
+    
+   }
+
+}

CALIBRATION/shms_aero_calib/gain.r

+paero_neg_gain = 1./13.0732, 1./11.1519, 1./12.7171, 1./11.7946, 1./10.7654, 1./8.72566, 1./12.9509, 
+paero_pos_gain = 1./7.45058, 1./10.1029, 1./10.7762, 1./13.0162, 1./10.5675, 1./10.9487, 1./10.8312, 

CALIBRATION/shms_aero_calib/gain_4985.r

+paero_neg_gain = 1./13.5843, 1./11.1229, 1./12.5078, 1./11.5904, 1./11.019, 1./8.54602, 1./12.5383, 
+paero_pos_gain = 1./7.3637, 1./9.83655, 1./10.4292, 1./12.8879, 1./10.45, 1./10.5741, 1./10.6287, 

CALIBRATION/shms_aero_calib/paero_calib.C

@@ -21,7 +21,18 @@
 #define MaxAdc 40.   //pC
 #define NBin   1000
 
-void paero_calib(string fname) {
+void paero_calib(string rootname = "") {
+
+  if (rootname=="") {
+  cout << "Please enter root file name...\n";
+  cin >> rootname;
+  }
+  
+  string fname = "../../ROOTfiles/" + rootname + ".root";
+  TString outputpdf;
+  outputpdf = "OUTPUT/" + rootname + ".pdf";
+  TString outputpng;
+  outputpng = "OUTPUT/" + rootname + ".png";
 
   TFile *f = new TFile(fname.c_str());
   TTree* tree;
@@ -55,17 +66,22 @@ void paero_calib(string fname) {
       hpos[i]->Fill(adc_pos[i], 1.);
       hneg[i]->Fill(adc_neg[i], 1.);
     }
+    
+    cout << "Entry " << ientry << "/" << nentries << '\n';
 
   }
 
   TCanvas* c1 = new TCanvas("adc_spec", "fADC spectra", 600, 800);
   c1->Divide(2,NPMT);
 
-  double flo_neg[NPMT] { 8., 8.,10., 8., 7., 5., 8.};
-  double fhi_neg[NPMT] {16.,13.,18.,14.,15.,10.,18.};
 
-  double flo_pos[NPMT] { 5., 6., 8.,10., 5., 6., 5.};
-  double fhi_pos[NPMT] {10.,12.,14.,17.,10.,12.,11.};
+  //ADC                {-1, -2, -3, -4, -5, -6, -7}
+  double flo_neg[NPMT] { 9., 7.,8., 8., 6., 5., 8.};
+  double fhi_neg[NPMT] {17.,14.,18.,15.,15.,12.,16.};
+
+  //ADC                {+1, +2, +3, +4, +5, +6, +7}
+  double flo_pos[NPMT] { 5., 6., 7.,9., 6., 7., 7.};
+  double fhi_pos[NPMT] {9.,14.,14.,17.,14.,14.,14.};
 
   float gain_pos[NPMT] {NPMT*0.};
   float gain_neg[NPMT] {NPMT*0.};
@@ -77,7 +93,7 @@ void paero_calib(string fname) {
       hpos[i]->Fit("gaus","","",flo_pos[i],fhi_pos[i]);
       hpos[i]->GetFunction("gaus")->SetLineColor(2);
       hpos[i]->GetFunction("gaus")->SetLineWidth(2);
-      gain_pos[i] = hpos[i]->GetFunction("gaus")->GetParameter(1);
+      gain_pos[i] = hpos[i]->GetFunction("gaus")->GetParameter(1);      
     }
     else
       hpos[i]->Draw();
@@ -86,11 +102,14 @@ void paero_calib(string fname) {
       hneg[i]->Fit("gaus","","",flo_neg[i],fhi_neg[i]);
       hneg[i]->GetFunction("gaus")->SetLineColor(2);
       hneg[i]->GetFunction("gaus")->SetLineWidth(2);
-      gain_neg[i] = hneg[i]->GetFunction("gaus")->GetParameter(1);
+      gain_neg[i] = hneg[i]->GetFunction("gaus")->GetParameter(1);    
     }
     else
       hneg[i]->Draw();
   }
+  
+  c1->Print(outputpdf);
+  c1->Print(outputpng);
 
   ofstream of;
   of.open("gain.r",ios::out);

CALIBRATION/shms_cal_calib/THcPShTrack.h

@@ -166,6 +166,7 @@ void THcPShTrack::SetEs(Double_t* alpha) {
       UInt_t ncol = 1;
       if (nblk > fNrows_pr) ncol = 2;
       (*iter)->SetEdep(adc*Ycor(Y,ncol)*alpha[nblk-1]);
+      //(*iter)->SetEdep(adc*alpha[nblk-1]);
     }
     else
       //Shower block, no coordinate correction.

CALIBRATION/shms_cal_calib/howto.txt

 How to calibrate SHMS calorimeter with real electrons.
 
-The calibration scripts reside in
-hallc_replay/CALIBRATION/shms_cal_calib directory. They consist of 3
+The calibration code resides in
+hallc_replay/CALIBRATION/shms_cal_calib directory. It consists of 3
 header files called THcPShHit.h,THcPShTrack.h, THcPShowerCalib.h, and
-a steering script pcal_calib.cpp.
+a steering script pcal_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, heavy gas Cherenkov, and TOF from
+quantities pertained to tracking, gas Cherenkovs, and TOF from
 hodoscopes. Hence it is convenient to calibrate on root files from
-full SHMS analysis. The root files are assumed in a linked ROOTfiles
-directory, under names shms_replay_<run_number>.root.
+full SHMS 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: SHMS delta range, the hodoscope beta range,
+thresholds on heavy and noble 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-s 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 pcal_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 pcal_calib under hcana, by issuing command
 
-.x pcal_calib.cpp+("run number").
+.x pcal_calib.cpp+("Prefix"),
 
-For instance, for calibration on shms_replay_464_50000.root file in
-ROOTfiles, the correct commad would be
+where Prefix is the prefix part in the name of root file.
 
-.x pcal_calib.cpp+("464_50000") .
+For instance, for calibration with
+shms_coin_replay_production_all_4313.root file in ROOTfiles, the
+correct command would be
+
+.x pcal_calib.cpp+("shms_coin_replay_production_all_4313") .
 
 Upon calibration, a canvas with representative plots will pop up. The
 calibration constants will be written in output file
-pcal.param.<run_number>, in a format suitable for plugging them into
-your pcal.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 THcPShowerCalib.h file, in the #define directives.
+pcal.param.<Prefix>, in a format suitable for plugging them into your
+pcal.param file for subsequent use. The representative canvas is saved
+in a pdf file <Prefix>.pdf, and histograms in <Prefix>.root file.

CALIBRATION/shms_cal_calib/input.dat

+-15. 27.	Delta range, %
+0.5 1.5	Beta range
+0.  	Heavy Gas Cherenkov, threshold on signals in p.e.
+10.  	Noble Gas Cherenkov, threshold on signals in p.e.
+5.	Minimum number of hits per channel required to be calibrated
+0. 3.	Range of uncalibrated energy deposition histogram
+500	Binning of uncalibrated energy deposition histogram
+-.1 1.25	Gaussian fit range around e- peak in the uncalibrated Edep histogram
+2.5	Gaussian fit range around mean e- peak in the uncalibrated Edep histogram 
+
+; Calibration constants for file shms_defocused.root (SHMS 9643, 9644, 9645, 9646, 9648, 9649, 9650, 9652 and 9654) 306065 events processed
+
+pcal_neg_gain_cor = 23.98, 26.50, 27.96, 17.52, 22.82, 16.46, 22.16, 29.95, 30.88, 43.38, 21.45, 22.89, 21.67, 55.90,
+pcal_pos_gain_cor = 11.72, 53.90, 23.71, 24.34, 28.44, 30.40, 23.79, 26.81, 28.26, 18.71, 26.45, 29.51, 22.84,  0.00,
+pcal_arr_gain_cor =  0.00,  0.00,  0.00,  8.94, 26.83, 22.93, 56.80, 38.60, 28.83, 62.52, 44.36,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00, 27.21, 29.68, 35.31, 27.88, 36.87, 33.00, 30.82, 35.13, 30.02, 34.52, 26.17,  0.00,  0.00,  0.00,  0.00,
+                     0.00, 45.97, 27.32, 37.10, 29.30, 31.68, 28.79, 31.26, 30.81, 26.81, 30.30, 33.84, 27.26,  0.00,  0.00,  0.00,
+                     0.00, 37.45, 27.89, 37.59, 25.33, 27.25, 34.88, 29.72, 32.25, 29.53, 31.04, 26.59, 25.45, 37.25, 18.66,  0.00,
+                     0.00, 39.42, 32.60, 29.71, 30.02, 30.66, 27.36, 31.91, 31.44, 27.97, 29.05, 32.15, 31.40, 39.06, 39.77,  0.00,
+                     0.00, 40.89, 28.40, 30.11, 31.55, 33.34, 26.43, 28.06, 32.35, 28.88, 32.93, 33.22, 36.02, 25.97, 33.14,  0.00,
+                     0.00, 33.96, 28.59, 29.46, 29.55, 34.78, 30.58, 39.43, 41.19, 40.14, 35.08, 34.17, 28.21, 28.97, 28.91,  0.00,
+                     18.83, 29.93, 28.87, 28.48, 27.35, 29.80, 32.57, 35.65, 43.36, 44.82, 49.44, 30.55, 33.73, 32.49, 37.63,  8.54,
+                     24.72, 29.79, 28.03, 28.83, 30.09, 30.23, 30.97, 32.67, 31.17, 34.81, 27.39, 30.92, 32.68, 32.03, 36.40,  0.00,
+                     18.14, 26.43, 28.73, 29.75, 27.81, 28.17, 28.83, 32.40, 30.13, 34.18, 31.61, 33.15, 32.58, 32.02, 29.99,  0.00,
+                     9.24, 25.62, 26.91, 24.97, 26.43, 31.29, 30.08, 27.68, 31.54, 33.71, 31.36, 32.03, 29.98, 33.83, 34.00,  0.00,
+                     0.00, 26.08, 23.04, 26.06, 29.22, 30.54, 28.60, 27.62, 26.95, 27.16, 27.36, 30.72, 28.00, 20.02,  0.00,  0.00,
+                     0.00, 31.38, 26.62, 32.41, 57.55, 28.71, 23.28, 26.16, 27.63, 22.19, 22.46, 36.94,  0.00,  0.00,  0.00,  0.00,
+                     0.00,  0.00,  0.00,  0.00, 42.28,  5.81, 45.36, 25.30, 21.53, 46.54, 16.55,  0.00,  0.00,  0.00,  0.00,  0.00,
+

CALIBRATION/shms_cal_calib/pcal_calib_7109.param

+; Calibration constants for file shms_replay_production_all_7109_500000, 26606 events processed
+
+pcal_neg_gain_cor =  0.00, 30.14, 25.74, 15.50, 19.70, 15.53, 26.19, 30.16, 17.18, 25.95, 20.34, 18.94, 19.31,  0.00,
+pcal_pos_gain_cor =  0.00, 15.52, 21.25, 21.53, 21.22, 21.46, 16.55, 15.88, 23.98, 18.74, 27.22, 27.62,  0.00,  0.00,
+pcal_arr_gain_cor =  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,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00,  0.00,  0.00, 25.01, 48.26, 56.36, 29.61,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00,  0.00, 27.66, 42.45, 29.87, 24.86, 71.85, 51.58,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00, 24.33, 26.15, 23.41, 50.80, 44.25, 37.29, 28.90, 28.76, 30.54,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00,  0.00, 25.66, 59.30, 15.68, 36.50, 34.78, 41.79, 60.85, 28.60, 44.03, 49.46, 76.67,  9.74,  0.00,  0.00,
+                     0.00,  0.00, 41.49, 37.45, 20.81, 42.18, 49.36, 37.50, 30.75, 30.73, 26.24, 21.50, 31.72, 29.17,  0.00,  0.00,
+                     0.00,  1.87, 13.88, 11.45, 14.07, 10.49, 15.18, 17.80, 12.16, 15.86, 22.85, 18.77, 27.38, 15.60, 36.10,  0.00,
+                     0.00, 40.56, 27.36, 27.88, 44.07, 37.18, 21.70, 40.66, 52.96, 65.24, 21.19, 55.98, 26.12, 23.60,  0.00,  0.00,
+                     0.00, 37.00, 16.72, 37.56, 21.27, 50.57, 39.43, 52.75, 28.40, 26.12, 50.23, 37.45, 50.94, 31.13,  0.00,  0.00,
+                     0.00,  0.00, 59.09, 33.38, 27.97, 37.41, 19.72, 15.23,  0.00,  0.00,  0.00,  0.00,  0.00, 31.92,  0.00,  0.00,
+                     0.00,  0.00,  0.00,  0.00, 42.05, 23.08, 18.68,  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,  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,