Shms dc calib (#141)

* added README file

* added README.md file for SHMS DC Calibrations

* Added readme file for shms dc calibration

* Added minor fix to the script to fix compilationerror in root 5.34.36

* minor fix

* get_tzero_per_wire_param.C

* fix to be compatible with root 5.34.36

* fix minor issues

* added stricter restriction on good t0 selection

CALIBRATION/shms_dc_calib/README.md

+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. 

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

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
  
  
  
-  
+    
 
 }

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
  
  

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
+
 
 
 

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");