Include scripts to calibrate the SHMS drift chambers

.gitignore

@@ -4,4 +4,6 @@ ROOTfiles/*
 hcana
 raw
 .root_history
+data_files/*
+root_files/*
 REPORT_OUTPUT/*

CALIBRATION/shms_dc_calib/run_Cal.C

+//SCRIPT TO RUN OVER ALL HMS DC CALIBRATION SCRIPTS AT ONCE, AND UPDATE THE
+//NECESSARY PARAMTER FILES hdriftmap.param and  hdc.param
+void run_Cal()
+{
+      
+  //User Input Run
+  int run_NUM;
+  cout << "Enter Run Number: " << endl;
+  cin >> run_NUM;
+
+  //Create input file with run number
+  ofstream fout;
+  fout.open("scripts/input_RUN.txt");
+  fout << run_NUM << endl;
+  fout.close();
+
+
+  //Create root and data files Directories if they dont exist
+  char *dir_root = "mkdir ./root_files/";
+  char *dir_data = "mkdir ./data_files/";
+
+  if (system(dir_root || dir_data) != 0) {
+    system(dir_root);
+    system(dir_data);
+  }
+
+  //Create run Directories if they dont exist
+  char *dir0 = Form("mkdir ./root_files/run%d", run_NUM);
+  char *dir1 = Form("mkdir ./data_files/run%d", run_NUM);
+
+  if (system(dir0 || dir1) != 0) {
+    system(dir0);
+    system(dir1);
+  }
+
+  
+
+  //change directories and execute scripts
+  gSystem->cd("./scripts");
+  gSystem->Exec("root -l -q get_pdc_time_histo.C");
+  
+  //Load and Loop over Make Class events to get individual drift times
+  gROOT->LoadMacro("wire_drift_times.C");
+  gROOT->ProcessLine("wire_drift_times t"); //process line allows one to execute interactive root commands from a script, such as this one
+  gROOT->ProcessLine("t.Loop()");
+  //gROOT->ProcessLine(".q");
+  gROOT->Reset();
+  
+  
+
+gSystem->cd("./scripts");
+
+  
+  //execute code to get t0 from each wire in each plane
+  gSystem->Exec("root -l -q -b get_wire_tzero.C");
+ 
+  //execute code to update pdc parameter file
+  gSystem->Exec("root -l -q update_pdcparam.C");
+
+  //execute code to get t0 corrected drift times
+  gSystem->Exec("root -l -q get_pdc_time_histo_tzero_corrected.C");
+  
+ //execute code to update LookUp Table
+  gSystem->Exec("root -l -q get_LookUp_Values.C");
+
+  
+  
+}

CALIBRATION/shms_dc_calib/scripts/get_LookUp_Values.C

+/*This code produces a lookup table necessary to convert drift times to
+  drift distances in the shms drift chambers
+*/
+
+#define NPLANES 12
+#define TOTAL_BINS 137
+
+void get_LookUp_Values() {
+
+   
+  //Read Run Number from txt file
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;
+ 
+  //Open root file containing drift time histos
+  TFile *f = new TFile(Form("../root_files/run%d/shms_dc_t0_corrected_%d.root", run_NUM, run_NUM),"READ");
+ 
+  //Define histogram array
+  TH1F *h[NPLANES];
+ 
+  //Define the number Drift Chamber planes
+  TString plane_names[NPLANES]={"1u1", "1u2", "1x1", "1x2", "1v1", "1v2", "2v2", "2v1", "2x2", "2x1", "2u2", "2u1"};
+ 
+  //Declare bin properties
+  int bin_t0[NPLANES];
+  int bin_final[NPLANES];					/*Array to store the bin number corresponding to last bin*/
+  int bin_Content[NPLANES];          /*Array to store the content (# events) corresponding to the bin with maximum content*/
+  double binContent_TOTAL[NPLANES];     /*Array to store sum of all bin contents for each plane*/
+  double binSUM[NPLANES];
+  int bin;
+  int binx;
+  double lookup_value[NPLANES];           /*Array to store lookup values for each plane*/ 
+  
+  //Create an output file to store lookup values  
+  ofstream ofs;
+  TString lookup_table = "../../../PARAM/SHMS/DC/pdriftmap_new.param";
+  ofs.open (lookup_table);
+ 
+ 
+  //Set headers for subsequent columns of data
+  ofs << Form("; Lookup Table: RUN %d", run_NUM) << "\n";
+  ofs << "; number of bins in Carlos's time to distance lookup table" << "\n";
+  ofs << Form("pdriftbins = %d", TOTAL_BINS+1) << "\n";
+  ofs << "; number of 1st bin in Carlos's table in ns" << "\n";
+  ofs << "pdrift1stbin=0" << "\n";
+  ofs << "; bin size in ns" << "\n";
+  ofs << "pdriftbinsz=2" << "\n";
+ 
+ 
+ 
+  //Loop over each plane of shms Drift Chambers (DC1 & DC2)
+
+  for (int ip=0; ip<NPLANES; ip++){
+   
+    TString drift_time_histo = "pdc"+plane_names[ip]+"_time: t0_corr"; 
+
+    //Get drift time histograms from root file
+    h[ip] = (TH1F*)f->Get(drift_time_histo);
+
+    //Get bin corresponding to t0 = 0 ns
+    bin_t0[ip] = h[ip]->GetXaxis()->FindBin(0.0);
+   
+    //Get final bin 
+    bin_final[ip] = bin_t0[ip] + TOTAL_BINS;
+   
+   
+   
+    //Find total BIN Content over entire integration range
+    binContent_TOTAL[ip] = 0; //set counter to zero
+
+    for (bin = bin_t0[ip]; bin <= bin_final[ip]; bin ++ ) {
+     
+      bin_Content[ip] = h[ip] -> GetBinContent(bin);
+     
+      binContent_TOTAL[ip] = bin_Content[ip] + binContent_TOTAL[ip];
+     
+      //   cout << "Bin: " << bin << endl;
+      //   cout << "Content " << bin_Content[ip] << endl;
+      //   cout << "Content SUM : " << binContent_TOTAL[ip] << endl;
+    }
+   
+    TString headers = "pwc" + plane_names[ip] + "fract=";      
+    ofs << headers;	  
+   
+    //Calculate LookUp Value
+   
+    binSUM[ip] = 0.0;
+    int bin_count = 0;
+   
+    for (bin = bin_t0[ip]; bin <= bin_final[ip]; bin++) {
+     
+      bin_Content[ip] = h[ip] -> GetBinContent(bin);
+      binSUM[ip] = binSUM[ip] + bin_Content[ip];
+     
+     
+      lookup_value[ip] = binSUM[ip] / binContent_TOTAL[ip];
+      bin_count = bin_count + 1;
+     
+      if (bin_count < = 8 ) {
+	ofs << setprecision(5) << lookup_value[ip] << fixed << ",";
+      }
+     
+      else if (bin_count >8 && bin_count < 138) {
+	ofs << setprecision(5) << lookup_value[ip] << ((bin_count+1) % 10 ? "," : "\n") << fixed; 
+      }
+      else {
+	ofs << setprecision(5) << lookup_value[ip] << fixed << endl;	  
+      }
+     
+    }
+
+  }                                            
+ 
+}
+
+
+	 
+  
+
+
+
+
+
+

CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo.C

+//Script to add necessary drift time histograms/plane from original root file to new root file
+
+
+#define NPLANES 12
+
+void get_pdc_time_histo()
+{
+
+  //Read Run Number from txt file
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;
+
+  //Create RUN Directories if they dont exist
+  char *dir0 = Form("mkdir ../root_files/run%d", run_NUM);
+  char *dir1 = Form("mkdir ../data_files/run%d", run_NUM);
+
+  if (system(dir0 || dir1) != 0) {
+    system(dir0);
+    system(dir1);
+  }
+
+  //open file
+  TFile *f = new TFile(Form("../../../ROOTfiles/pdc_replay_%d.root", run_NUM), "READ");
+
+  //create new file
+  TFile *g = new TFile(Form("../root_files/run%d/shms_dc_time_%d.root", run_NUM, run_NUM), "RECREATE"); // create new file to store histo
+
+  f->cd();
+
+  //Get the tree
+  TTree *tree = (TTree*)f->Get("T");
+  TString SPECTROMETER="P";
+  TString DETECTOR="dc";
+  TString plane_names[NPLANES]={"1u1", "1u2", "1x1", "1x2", "1v1", "1v2", "2v2", "2v1", "2x2", "2x1", "2u2", "2u1"};
+  
+  //Declare Variables to Loop Over
+  Int_t Ndata[NPLANES];
+  Double_t pdc_time[NPLANES][1000];
+    
+  //Declare Histogram array to store AVG drift times per plane
+  TH1F* h[NPLANES];
+	
+  g->cd();
+	
+  //Loop over each plane
+  for(Int_t ip=0; ip<NPLANES; ip++){
+    TString base_name = SPECTROMETER+"."+DETECTOR+"."+plane_names[ip];
+    TString ndata_name = "Ndata."+base_name+".time";
+    TString drift_time = base_name+".time";
+  
+    TString drift_time_histo = "pdc"+plane_names[ip]+"_time"; 
+    TString title = "pdc"+plane_names[ip]+"_drifttime";
+     
+    //Set Branch Address
+    tree->SetBranchAddress(drift_time, &pdc_time[ip][0]);   
+    tree->SetBranchAddress(ndata_name, &Ndata[ip]);  /* Ndata represents number of triggers vs number of hits that each trigger produced.
+							A hit is refer to as when a trigger(traversing particle), ionizes the WC gas and ionized
+							electrons reach the rearest sense wire, producing a detectable signal in the O'scope */
+  
+    //Create Histograms
+    h[ip] = new TH1F(drift_time_histo, title, 200, -50, 350);  //set time to 400 ns/200 bins = 2ns/bin
+  }
+ 
+
+ 
+  //Declare number of entries in the tree
+  Long64_t nentries = tree->GetEntries(); //number of triggers (particles that passed through all 4 hodo planes)
+ 
+  //Loop over all entries
+  for(Long64_t i=0; i<nentries; i++)
+    {
+      tree->GetEntry(i);
+    
+    
+      //Loop over number of hits for each trigger in each DC plane 
+      for(ip=0; ip<NPLANES; ip++){
+       
+    
+	for(Int_t j=0; j<Ndata[ip]; j++){
+	
+	  h[ip]->Fill(pdc_time[ip][j]); 
+	}
+       
+      }
+
+    }
+
+
+
+		
+  //Write histograms to file
+  g->Write();
+
+
+
+}

CALIBRATION/shms_dc_calib/scripts/get_pdc_time_histo_tzero_corrected.C

+//Script to add t0 correction to SHMS DC drift times
+
+#define NPLANES 12
+
+void get_pdc_time_histo_tzero_corrected()
+{
+
+	  
+  //read run number from input file
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;   
+
+  TString run = Form("run%d", run_NUM);
+
+
+	  //open file
+	  TFile *f = new TFile(Form("../../../ROOTfiles/pdc_replay_%d.root", run_NUM), "READ");
+
+	  //updates file
+	  TFile *g = new TFile(Form("../root_files/run%d/shms_dc_t0_corrected_%d.root", run_NUM, run_NUM), "UPDATE"); // create new file to store histo
+
+      f->cd();
+
+     //Get the tree
+     TTree *tree = (TTree*)f->Get("T");
+
+	TString SPECTROMETER="P";
+	TString DETECTOR="dc";
+	TString plane_names[NPLANES]={"1u1", "1u2", "1x1", "1x2", "1v1", "1v2", "2v2", "2v1", "2x2", "2x1", "2u2", "2u1"};
+
+    //Declare Variables to Loop Over
+    Int_t Ndata[NPLANES];
+    Double_t pdc_time[NPLANES][1000];
+    
+    //Declare Histogram array to store AVG drift times per plane
+    TH1F* h[NPLANES];
+	
+	g->cd();
+	
+	//Loop over each plane
+	for(Int_t ip=0; ip<NPLANES; ip++){
+		TString base_name = SPECTROMETER+"."+DETECTOR+"."+plane_names[ip];
+		TString ndata_name = "Ndata."+base_name+".time";
+		TString drift_time = base_name+".time";
+
+		TString drift_time_histo = "pdc"+plane_names[ip]+"_time: t0_corr"; 
+        TString title = "pdc"+plane_names[ip]+"_drifttime: t0-corrected";
+     
+     //Set Branch Address
+     tree->SetBranchAddress(drift_time, pdc_time[ip]);   
+     tree->SetBranchAddress(ndata_name, &Ndata[ip]);  /* Ndata represents number of triggers vs number of hits that each trigger produced.
+                                                      A hit is refer to as when a trigger(traversing particle), ionizes the WC gas and ionized
+                                                      electrons reach the rearest sense wire, producing a detectable signal in the O'scope */
+	
+     //Create Histograms
+     h[ip] = new TH1F(drift_time_histo, title, 200, -50, 350);  //set time to 400 ns/200 bins = 2ns/bin
+}
+	
+	
+	//open and read tzero data file
+    ifstream ifs;
+    ifs.open("../data_files/" + run + "/tzero.dat");
+    
+	double t_zero_offsets[NPLANES];
+
+     for (ip=0; ip < 12; ip++) {	 
+	 ifs >> t_zero_offsets[ip];  //add tzero offsets to array
+   }
+   
+    //Declare number of entries in the tree
+    Long64_t nentries = tree->GetEntries(); //number of triggers (particles that passed through all 4 hodo planes)
+
+    //Loop over all entries
+    for(Long64_t i=0; i<nentries; i++)
+    {
+	tree->GetEntry(i);
+    
+    
+    //Loop over number of hits for each trigger in each DC plane 
+    for(ip=0; ip<NPLANES; ip++){
+   
+ 
+    
+    for(Int_t j=0; j<Ndata[ip]; j++){
+	
+	h[ip]->Fill(pdc_time[ip][j] - t_zero_offsets[ip]); //add t0 offset correction 
+       }
+      
+      
+
+			
+	
+	}
+
+}
+
+
+
+		
+//Write histograms to file
+g->Write();
+
+
+
+}

CALIBRATION/shms_dc_calib/scripts/get_wire_tzero.C

+
+
+/*Script to extract reference time "t0"  for each sense wire in a given HMS Wire Chamber Plane with COSMIC RUNS.  
+20% (MAX BIN CONTENT) is calculated per wire, and the corresponding bin is fitted linearly about +/- 
+a certain number of bins and this fit is extrapolated to y=0(x-axis). The extrapolated value is take to be t0*/
+
+#include <vector>
+#include <TMath>
+
+#define NPLANES 12
+
+void get_wire_tzero()
+{
+  using namespace std;
+
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;   
+
+  //check if tzero_weighted_avg text file exists (if it does, DELETE IT, otherwise new values will be appended to it, in addition to pre-existing tzero values)
+  std::ifstream stream(Form("../data_files/run%d/tzero_weighted_avg_run%d.txt",run_NUM, run_NUM));
+ if (stream.good())
+  {
+  gSystem->Exec(Form("rm ../data_files/run%d/tzero_weighted_avg_run%d.txt",run_NUM, run_NUM));
+  }
+
+ TString run = Form("run%d", run_NUM);
+
+ //Declare plane names to loop over
+ TString plane_names[NPLANES]={"1u1", "1u2", "1x1", "1x2", "1v1", "1v2", "2v2", "2v1", "2x2", "2x1", "2u2", "2u1"};
+
+ //Declare a root file array to store individual DC cell drift times
+ TString root_file;
+ TFile *f[NPLANES];
+   
+ int total_wires;  //integer to store total sense wires for a plane chosen by the user
+        
+ //Loop over all planes
+ for (int ip = 0; ip < NPLANES; ip++){
+
+   //READ root file
+   root_file = "../root_files/"+run+"/shms_DC_"+plane_names[ip]+Form("_%d.root",run_NUM);
+   f[ip] = new TFile(root_file, "READ");
+
+   //Create a file output file stream object to write t0 values to data file
+   ofstream ofs;
+   TString t_zero_file = "../data_files/" + run + "/hdc_"+plane_names[ip]+Form("tzero_run%d.dat", run_NUM);
+   ofs.open (t_zero_file);
+
+   //Set headers for subsequent columns of data
+   ofs << "#WIRE " << "   "  << "t0" << "   " << "t0_err" << "   " << " entries " << endl;
+
+   //Create root file to store fitted wire drift times histos and "t0 vs. wirenum"
+   TString output_root_file = "../root_files/"+run+"/shmsDC_"+plane_names[ip]+Form("run%d_fitted_histos.root", run_NUM);
+   TFile *g = new TFile(output_root_file,"RECREATE");
+
+   f[ip]->cd();  //change to file containing the wire drift times histos
+ 
+   int total_wires;  //integer to store total sense wires for a plane chosen by the user
+   
+   //Set variables depending on which plane is being studied
+   if(ip == 0 || ip == 1 || ip == 4 || ip == 5 || ip == 6 || ip == 7 || ip == 10 || ip == 11) {
+     TH1F *cell_dt[107]; //declare array of histos to store drift times     
+     total_wires=107; 
+
+     //Declare bin properties for given sense wires in a plane
+
+     int bin_max[107];                    /*Array to store the bin number corresponding to the drift time distribution peak*/
+     int bin_maxContent[107];             /*Array to store the content (# events) corresponding to the bin with maximum content*/
+     double time_max[107];                /*Array to store the x-axis(drift time (ns)) corresponding to bin_max*/
+     double twenty_perc_maxContent[107];  /*Array to store 20% of maximum bin content (peak)*/						     
+     double ref_time[107];               /*Array to store reference times for each sense wire*/
+
+   }
+
+   else if(ip == 2 || ip == 3 || ip == 8 || ip == 9) {
+     TH1F *cell_dt[79];
+     total_wires=79;      
+   
+     int bin_max[79];                                 
+     int bin_maxContent[79];                           
+     double time_max[79];                               
+     double twenty_perc_maxContent[79];                
+     double ref_time[79];          
+
+   }	   
+   
+ 	
+   /*Get wire histos from root file and loop over each 
+     sense wire of a plane in shms Drift Chambers (DC1 or DC2)*/
+ 
+   for (int sensewire=1; sensewire<=total_wires; sensewire++){
+
+     //Get title of histos in root file
+     TString drift_time_histo = Form("wire_%d", sensewire); 
+ 
+     //Get drift time histograms from root file
+     cell_dt[sensewire-1] = (TH1F*)f[ip]->Get(drift_time_histo);
+
+     
+     //Get bin with Maximum Content
+     bin_max[sensewire-1] = cell_dt[sensewire-1]->GetMaximumBin();
+
+     //Get content of bin_max
+     bin_maxContent[sensewire-1] = cell_dt[sensewire-1]->GetBinContent(bin_max[sensewire-1]);
+     
+     //Get time (ns) [x-axis] corresponding to bin_max 
+     time_max[sensewire-1] = cell_dt[sensewire-1]->GetXaxis()->GetBinCenter(bin_max[sensewire-1]);
+     
+     //Calculate 20% of max content
+     twenty_perc_maxContent[sensewire-1] = bin_maxContent[sensewire-1] * 0.20;
+     
+
+   }
+   
+   
+   
+   //****************************************************//
+   //Determine  which bin has around 20% max_BinContent *//
+   //****************************************************//
+   
+   
+   //Declarations
+   int content_bin;      //stores content for each bin
+   int counts;           //a counter used to count the number of bins that have >20% max bin content for a plane 
+   int bin;              //store bin number
+   int j;                //jth bin, used to loop over n bins 
+   
+   //Declare vector arrays 
+   vector<int> content;               //stores bin content
+   vector <int> bin_num;           //stores bin number
+   
+   
+   //Loop over each wire 
+   for(sensewire=1; sensewire<=total_wires; sensewire++) {
+     
+     //Loop over each bin for individual wire drift time histo
+     for(bin=0; bin < bin_max[sensewire-1]; bin++) {
+       
+       content_bin = cell_dt[sensewire-1]->GetBinContent(bin);              //get bin content for all bins in a wire
+       
+       content.push_back(content_bin);                                      //add bin content to array
+       bin_num.push_back(bin);                                              //add bin number to array
+       
+       
+       // check if 2 bin contents have been stored and examine if these contents exceed or not 20% of peak
+       if (content.size() == 2) {
+	 
+	 //initialize counter to count how many bin contents >= 20%
+	 counts = 0;
+	 
+	 // Loop over 2 bin contents stored in array content
+	 for (j=0; j<2; j++){
+	   
+	   if(content[j] > =  twenty_perc_maxContent[sensewire-1]){
+	     counts = counts+1;
+             
+	     if(counts >= 2) { goto stop;}
+	     
+	     
+	   }
+	   
+	   content.clear();
+	   bin_num.clear();
+	   
+	 }
+	 
+       }
+     }
+     
+     //Print the time(ns) and BIN NUM corresponding to 20% of MAX content 
+     //if 2/2 elements exceeds 20% of Max content (for each plane)
+     
+   stop:
+     ref_time[sensewire-1] = cell_dt[sensewire-1] ->GetXaxis() -> GetBinCenter(bin_num[0]); //Get time corresponding ~20% Max BIN CONTENT  
+     
+     //cout << " ******* " << "Wire " << sensewire << " ******* " << endl;
+     //cout << "time (20% of Max BIN): " << ref_time[sensewire-1] << " ns" << endl;
+     //cout << "BIN: " << bin_num[0] << endl;
+     
+     
+     //*********************************************************//
+     //*******Extract the "t0" Using a Fitting Procedure********//
+     //*********************************************************//
+     
+     //Declarations
+     int time_init;           //start fit value 
+     int time_final;          //end fit value
+     int t_zero;
+     int entries;             //entries for each wire
+     
+     double m;                //slope
+     double y_int;            //y-intercept
+     double m_err;
+     double y_int_err;
+     double t_zero_err;
+     
+     //Get time corresponding to bin (fit range) 
+     time_init = cell_dt[sensewire-1] -> GetXaxis() -> GetBinCenter(bin_num[0]-5); //choose bin range over which to fit
+     time_final = cell_dt[sensewire-1] -> GetXaxis() -> GetBinCenter(bin_num[0]+5); 
+     
+     //Create Fit Function
+     TF1* tZero_fit = new TF1("tZero_fit", "[0]*x + [1]", time_init, time_final);
+     
+     //Set Parameter Names and Values
+     tZero_fit->SetParName(0, "slope");
+     tZero_fit->SetParName(1, "y-int");
+     tZero_fit->SetParameter(0, 1.0);
+     tZero_fit->SetParameter(1, 1.0);
+     
+     //Fit Function in specified range
+     cell_dt[sensewire-1]->Fit("tZero_fit", "QR");
+     
+     //Get Parameters and their errors
+     m = tZero_fit->GetParameter(0);
+     y_int = tZero_fit->GetParameter(1);
+     m_err = tZero_fit->GetParError(0);
+     y_int_err = tZero_fit->GetParError(1);
+     
+     //Calculate error on t0 using error propagation method of expanding partial derivatives
+     t_zero = - y_int/m;
+     t_zero_err = sqrt(y_int_err*y_int_err/(m*m) + y_int*y_int*m_err*m_err/(m*m*m*m) );
+     entries = cell_dt[sensewire-1]->GetEntries();  //number of entries (triggers) per wire
+     
+     //Write "t0" values to file
+     ofs << sensewire << "          " << t_zero << "          " << t_zero_err << "          " << entries << endl;
+     
+     //Change to output root file and write fitted histos to file
+     g->cd();
+     cell_dt[sensewire-1]->Write();
+     
+   }
+   
+   // Make Plot of t0 versus Wire Number 
+   
+   TCanvas *t = new TCanvas("t", "", 2000,500);
+   t->SetGrid();
+
+   
+   TGraphErrors *graph = new TGraphErrors(t_zero_file, "%lg %lg %lg");
+   graph->SetName("graph");
+   TString title = "DC"+plane_names[ip]+": t0 versus sensewire";
+   graph->SetTitle(title);
+   graph->SetMarkerStyle(20);
+   graph->SetMarkerColor(1);
+   graph->GetXaxis()->SetLimits(0., total_wires);
+   graph->GetXaxis()->SetTitle("Wire Number");
+   graph->GetXaxis()->CenterTitle();
+   graph->GetYaxis()->SetTitle("t-Zero (ns)");
+   graph->GetYaxis()->CenterTitle();
+   graph->GetYaxis()->SetRangeUser(-50.0, 50.0);
+   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   //
+   //*****************************************************************************************//
+   
+   
+   //open t0 dat file
+   ifstream ifs;
+   ifs.open (t_zero_file);
+   string line;
+   
+   //open new data file to write updated t0 values
+   TString t_zero_file_corr = "../data_files/" + run + "/hdc_"+plane_names[ip]+Form("tzero_run%d_updated.txt", run_NUM);
+   ofs.open(t_zero_file_corr);
+   ofs << " #Wire " << "     " << " t_zero " << "     " << " t_zero_err " << "     " << " entries " << endl; 
+   
+   //Initialize variables related to weighted avg
+   double sum_NUM;  //numerator of weighted avg
+   double sum_DEN;   //denominator of weighted avg
+   double weighted_AVG;
+   double weighted_AVG_err; 
+  
+   //set them to zero to start sum inside while loop 
+   sum_NUM = 0.0;
+   sum_DEN = 0.0;
+   
+   weighted_AVG;
+   weighted_AVG_err; 
+   
+   //read line bt line the t_zero_file
+   while(getline(ifs, line)) {
+     if(!line.length()|| line[0] == '#')
+       continue;
+     //	sensewire = 0, t_zero = 0.0, t_zero_err = 0.0, entries = 0 ; //set values to zero
+     
+     sscanf(line.c_str(), "%d %d %lf %d", &sensewire, &t_zero, &t_zero_err, &entries); //assign each of the variables above a data in the t_zero_file
+     
+     //Check if entries for each sensewire exceeds a certain number of events
+     
+     if (entries>300 && t_zero < 30) {
+	
+       //Calculate the weighted average of t0s
+       sum_NUM = sum_NUM + t_zero/(t_zero_err*t_zero_err);
+       sum_DEN = sum_DEN + 1.0/(t_zero_err*t_zero_err);      
+       
+       //cout << "sum_NUM : " << sum_NUM << endl;  
+       //cout << "sum_DEN : " << sum_DEN << endl;  
+       
+    
+
+
+       ofs << sensewire << "        " << t_zero << "        " << t_zero_err << "        " << entries << endl;
+
+       
+       
+     }
+     
+   }
+   
+   
+   
+   weighted_AVG = sum_NUM / sum_DEN;
+   weighted_AVG_err = sqrt( 1.0 / sum_DEN );
+   
+   
+   
+   //open new data file to write weighted average of updated t_zero values
+   
+   TString t_zero_AVG = Form("../data_files/run%d/tzero_weighted_avg_run%d.txt", run_NUM, run_NUM);
+   
+   ofstream ofile;
+   ofile.open(t_zero_AVG, std::ofstream::out | std::ofstream::app); //open file in and output and append mode
+   
+   ofile << " #weighted_AVG " << "     " << " DC plane: " <<  plane_names[ip] << endl; 
+   ofile << weighted_AVG << endl;
+   
+  
+   
+   
+   
+   ifs.close();
+
+   // Make Plot of t0 versus Wire Number for entries > 300 events
+
+   TCanvas *t1 = new TCanvas("t1", "", 2000,500);
+   t1->SetGrid();
+
+   //TString mygraph = "hdc"+plane_names[ip]+Form("_t_zero_run%d.txt", run);
+   TGraphErrors *graph1 = new TGraphErrors(t_zero_file_corr, "%lg %lg %lg");
+   graph1->SetName("graph1");
+   TString title1 = "hdc"+plane_names[ip]+": t0 versus sensewire_corrected";
+   graph1->SetTitle(title1);
+   graph1->SetMarkerStyle(20);
+   graph1->SetMarkerColor(1);
+   //graph1->GetXaxis()->SetLimits(0., total_wires);
+   graph1->GetXaxis()->SetTitle("Wire Number");
+   graph1->GetXaxis()->CenterTitle();
+   graph1->GetYaxis()->SetTitle("t-Zero (ns)");
+   graph1->GetYaxis()->CenterTitle();
+   graph1->GetYaxis()->SetRangeUser(-50.0, 50.0);
+   graph1->Draw("AP");
+   t1->Update();
+
+   // Draw TLine
+   TLine *wght_avg = new TLine(t1->GetUxmin(), weighted_AVG, t1->GetUxmax(), weighted_AVG);
+   wght_avg->SetLineColor(kRed);
+   wght_avg->SetLineWidth(2);
+   wght_avg->SetLineStyle(2);
+   wght_avg->Draw();
+   
+   //Add text to canvas
+   TLatex* ltx1 =  new TLatex();
+   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
+
+   
+
+
+
+
+
+
+
+
+ }
+ 
+ 
+
+
+
+
+}

CALIBRATION/shms_dc_calib/scripts/input_RUN.txt

+407

CALIBRATION/shms_dc_calib/scripts/update_pdcparam.C

+//This scirpt will produce an updated version of hdc.param file, with
+//the necessary t-zero corrections
+#define time_shift 1280.0
+
+void update_pdcparam()
+{
+	  
+  //read run number from input file
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;   
+
+  TString run = Form("run%d", run_NUM);
+
+    int lin_NUM = 0;
+	string t_zero[12];
+	double tzero[12];
+	string line;
+   //open t_zero file
+   ifstream ifs;
+   ifs.open("../data_files/"+ run +"/tzero_weighted_avg_" + run + ".txt");
+   
+
+while (getline(ifs, line))
+  {
+    
+    istringstream ss(line);
+    char id;
+    
+    if ( ss >> t_zero)
+      {
+		
+	if (id != '#') //skip comments   
+	  {
+	    //count lines  
+	    lin_NUM = lin_NUM + 1;
+	    cout << lin_NUM << endl;    
+	       t_zero[lin_NUM-1] = line;
+	       tzero[lin_NUM-1] = atof(t_zero[lin_NUM-1].c_str()); // convert string to double
+	    cout << tzero[lin_NUM-1] << endl;
+	  }       
+	
+      }     
+    
+  }
+ifs.close();
+
+//Update pdc.param parameter file 
+TString pdc_param = "../../../PARAM/SHMS/DC/pdc_tracking_new.param";
+ofstream ofs(pdc_param);
+
+// ofs << ";---------------------------------------------------------------------" << endl;
+// ofs <<"; SHMS_TRACKING"<< endl;
+// ofs <<"; CTP parameter file containing all tracking parameters for the SHMS "<< endl;
+// ofs <<";----------------------------------------------------------------------"<< endl;
+// ofs <<"; sigma of wire chamber resolution for each plane "<< endl;
+// ofs <<" pdc_sigma = 0.020 "<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"             0.020"<< endl;
+// ofs <<"  pdc_tdc_min_win = -25000,-25000,-25000,-25000,-25000,-25000 "<< endl;
+// ofs <<"  		    -25000,-25000,-25000,-25000,-25000,-25000 "<< endl;
+// ofs <<"  pdc_tdc_max_win = 25000,25000,25000,25000,25000,25000 "<< endl;
+// ofs <<"                    25000,25000,25000,25000,25000,25000 "<< endl;
+// ofs <<"; hms drift chamber tdc's time per channel "<< endl;
+// ofs <<"        pdc_tdc_time_per_channel = -0.10 "<< endl;
+// ofs <<"; hms zero time for drift chambers	!DECREASING this number moves the hdtime plots to LOWER time. "<< endl;
+// ofs <<"pdc_plane_time_zero =  ";
+
+
+//*****************************************************************
+//output all t_0 corrected values to pdc.param    
+for (int i=0; i<12; i++) { 
+{
+if (i < = 5){
+ofs  <<  time_shift - tzero[i]   << ","; 
+}
+if (i ==6) {ofs << "\n" << time_shift - tzero[6]  << ",";}
+else if (i>6 && i <11) {
+ofs << time_shift - tzero[i] << ","; 
+}
+if (i==11){ ofs << time_shift - tzero[i] << endl;}
+}
+}
+//*****************************************************************
+// ofs << "\n";
+// ofs <<"; Dave Abbott's wire velocity correction "<< endl;
+// ofs <<"pdc_wire_velocity = 12.0 "<< endl;
+// ofs <<"pdc_central_time = 7,9,3,4,6,5 "<< endl;
+// ofs << "                   7,5,3,4,6,6" << endl;
+ofs.close();
+
+//create a t_zero data file copy in another directory that will also use these values
+TString tzero_dat = "../data_files/" + run + "/tzero.dat";
+ofstream ofs(tzero_dat);
+
+for (int i=0; i<12; i++) 
+{
+ofs  <<  tzero[i]  << endl; 
+}
+ 
+}

CALIBRATION/shms_dc_calib/scripts/wire_drift_times.C

+#define wire_drift_times_cxx
+#include "wire_drift_times.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#define NPLANES 12
+
+void wire_drift_times::Loop()
+{
+  //   In a ROOT session, you can do:
+  //      Root > .L wire_drift_times.C
+  //      Root > wire_drift_times t
+  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
+  //      Root > t.Show();       // Show values of entry 12
+  //      Root > t.Show(16);     // Read and show values of entry 16
+  //      Root > t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  //Read Run Number from txt file
+  int run_NUM;
+  TString f0 = "input_RUN.txt";
+  ifstream infile(f0);
+  infile >> run_NUM;
+
+  TString run = Form("run%d", run_NUM);
+  //Declare plane names to loop over
+  TString plane_names[NPLANES]={"1u1", "1u2", "1x1", "1x2", "1v1", "1v2", "2v2", "2v1", "2x2", "2x1", "2u2", "2u1"};
+
+  //Declare a root file array to store individual DC cell drift times
+  TString root_file[NPLANES];
+  TFile *g[NPLANES];
+   
+  //integer to store total sense wires for a plane chosen by the user
+  static const Int_t total_wires_x  = 79;  
+  static const Int_t total_wires_uv = 107;
+     
+  Long64_t nbytes = 0, nb = 0;
+   
+  //Loop over all planes
+  for (int ip = 0; ip < NPLANES; ip++){
+
+    //Initialize a root file array to store individual DC cell drift times
+    root_file[ip] = "../root_files/" + run + "/shms_DC_"+plane_names[ip]+Form("_%d.root", run_NUM);
+    g[ip] = new TFile(root_file[ip], "RECREATE");
+    g[ip]->cd();
+
+    /*========================PLANES 1X1,1X2,2X1,2X2=====================================*/
+	
+    //If specific planes are encountered, treat them as follows:
+
+    /*PLANE 1U1, 1V1, 2U1, 2V1*/				
+    //If specific planes are encountered, treat them as follows:
+    if(ip == 0 || ip == 1 || ip == 4 || ip == 5 || ip == 6 || ip == 7 || ip == 10 || ip == 11) {
+
+      TH1F *cell_dt[total_wires_uv];    
+      TH2F *wire_vs_dt = new TH2F("wire_vs_dt", "", 200., -50., 350., 107., 0.,107.);
+    
+      //Initialize wire drift time histograms
+      for (int wirenum=1; wirenum<=total_wires_uv; wirenum++){
+	cell_dt[wirenum-1] = new TH1F(Form("wire_%d", wirenum), "", 200., -50., 350.);
+      }
+	
+      //Loop over all entries (triggers or events)   
+      for (Long64_t jentry=0; jentry<nentries; jentry++) {
+	Long64_t ientry = LoadTree(jentry);
+	if (ientry < 0) break;
+	nb = fChain->GetEntry(jentry);   nbytes += nb;
+	// if (Cut(ientry) < 0) continue;
+ 
+	if (ip == 0) {
+	  for (int i=0; i< Ndata_P_dc_1u1_wirenum; i++){
+	    wirenum = int(P_dc_1u1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1u1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1u1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1u1_time[i]);
+	    wire_vs_dt->Fill(P_dc_1u1_time[i], P_dc_1u1_wirenum[i]);
+
+	  }
+	}
+							
+	if (ip == 1) {
+	  for (int i=0; i< Ndata_P_dc_1u2_wirenum; i++){
+	    wirenum = int(P_dc_1u2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1u2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1u2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1u2_time[i]);
+	    wire_vs_dt->Fill(P_dc_1u2_time[i], P_dc_1u2_wirenum[i]);
+
+	  }
+	}						
+			
+	if (ip == 4) {
+	  for (int i=0; i< Ndata_P_dc_1v1_wirenum; i++){
+	    wirenum = int(P_dc_1v1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1v1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1v1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1v1_time[i]);
+	    wire_vs_dt->Fill(P_dc_1v1_time[i], P_dc_1v1_wirenum[i]);
+
+	  }
+	}		
+			
+	if (ip == 5) {
+	  for (int i=0; i< Ndata_P_dc_1v2_wirenum; i++){
+	    wirenum = int(P_dc_1v2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1v2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1v2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1v2_time[i]);
+	    wire_vs_dt->Fill(P_dc_1v2_time[i], P_dc_1v2_wirenum[i]);
+
+	  }
+	}
+
+	if (ip == 6) {
+	  for (int i=0; i< Ndata_P_dc_2v2_wirenum; i++){
+	    wirenum = int(P_dc_2v2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2v2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2v2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2v2_time[i]);
+	    wire_vs_dt->Fill(P_dc_2v2_time[i], P_dc_2v2_wirenum[i]);
+
+	  }
+	}
+
+	if (ip == 7) {
+	  for (int i=0; i< Ndata_P_dc_2v1_wirenum; i++){
+	    wirenum = int(P_dc_2v1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2v1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2v1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2v1_time[i]);
+	    wire_vs_dt->Fill(P_dc_2v1_time[i], P_dc_2v1_wirenum[i]);
+
+	  }
+	}
+
+	if (ip == 10) {
+	  for (int i=0; i< Ndata_P_dc_2u2_wirenum; i++){
+	    wirenum = int(P_dc_2u2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2u2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2u2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2u2_time[i]);
+	    wire_vs_dt->Fill(P_dc_2u2_time[i], P_dc_2u2_wirenum[i]);
+
+	  }
+	}
+
+	if (ip == 11) {
+	  for (int i=0; i< Ndata_P_dc_2u1_wirenum; i++){
+	    wirenum = int(P_dc_2u1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2u1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2u1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2u1_time[i]);
+	    wire_vs_dt->Fill(P_dc_2u1_time[i], P_dc_2u1_wirenum[i]);
+
+	  }
+	}						
+					
+					
+					
+      }
+    }
+	
+    if(ip == 2 || ip == 3 || ip == 8 || ip == 9) {
+
+      TH1F *cell_dt[total_wires_x];    
+      TH2F *wire_vs_dt = new TH2F("wire_vs_dt", "", 200., -50., 350., 79., 0., 79.);
+    
+      //Initialize wire drift time histograms
+      for (int wirenum=1; wirenum<=total_wires_x; wirenum++){
+	cell_dt[wirenum-1] = new TH1F(Form("wire_%d", wirenum), "", 200., -50., 350.);
+      }
+	
+      //Loop over all entries (triggers or events)   
+      for (Long64_t jentry=0; jentry<nentries; jentry++) {
+	Long64_t ientry = LoadTree(jentry);
+	if (ientry < 0) break;
+	nb = fChain->GetEntry(jentry);   nbytes += nb;
+	// if (Cut(ientry) < 0) continue;
+ 
+	if (ip == 2) {
+	  for (int i=0; i< Ndata_P_dc_1x1_wirenum; i++){
+	    wirenum = int(P_dc_1x1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1x1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1x1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1x1_time[i]);
+	    wire_vs_dt->Fill(P_dc_1x1_time[i], P_dc_1x1_wirenum[i]);
+
+	  }
+	}
+							
+	if (ip == 3) {
+	  for (int i=0; i< Ndata_P_dc_1x2_wirenum; i++){
+	    wirenum = int(P_dc_1x2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_1x2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_1x2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_1x2_time[i]);
+	    wire_vs_dt->Fill(P_dc_1x2_time[i], P_dc_1x2_wirenum[i]);
+
+	  }
+	}						
+			
+	if (ip == 8) {
+	  for (int i=0; i< Ndata_P_dc_2x2_wirenum; i++){
+	    wirenum = int(P_dc_2x2_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2x2_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2x2_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2x2_time[i]);
+	    wire_vs_dt->Fill(P_dc_2x2_time[i], P_dc_2x2_wirenum[i]);
+
+	  }
+	}		
+			
+	if (ip == 9) {
+	  for (int i=0; i< Ndata_P_dc_2x1_wirenum; i++){
+	    wirenum = int(P_dc_2x1_wirenum[i]);
+	    //cout << " wire num: " << P_dc_2x1_wirenum[i] << endl;
+	    //cout << "Time: " << P_dc_2x1_time[i] << endl;
+
+	    //Fill the Histograms
+	    cell_dt[wirenum-1]->Fill(P_dc_2x1_time[i]);
+	    wire_vs_dt->Fill(P_dc_2x1_time[i], P_dc_2x1_wirenum[i]);
+
+	  }
+	}						
+					
+					
+					
+      }
+    }
+					
+
+
+    //Write wire drift time histos to file
+    g[ip]->Write();
+    cout << "EVERYTHING OK in plane:" << ip << endl;	
+
+  }
+
+
+  //     cout << "\r          \r" << (float)sensewire / total_wires * 100.0 << "%" << flush; 
+
+}