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.h"
#include <iomanip>
#include <iostream>

#define NPLANES 12
  using namespace std;
void get_wire_tzero()
{
 

  int run_NUM;
  Long64_t num_evts;        //added
  string input_file;   //added

  TString f0 = "input_RUN.txt";
  ifstream infile(f0);
  infile >> input_file >> run_NUM >> num_evts;   //read the input file (added)   

  //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.dat",run_NUM, run_NUM));
 if (stream.good())
  {
  gSystem->Exec(Form("rm ../data_files/run%d/tzero_weighted_avg_run%d.dat",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"};
 int fNWires[NPLANES] = {107, 107, 79, 79, 107, 107, 107, 107, 79, 79, 107, 107};

 //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_wire_histos.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 + "/shms_dc_"+plane_names[ip]+Form("tzero_run%d.dat", run_NUM);
   ofs.open (t_zero_file);

   //Set headers for subsequent columns of data
   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);
   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

   //INITIALIZE VARIABLES
   int total_wires;
   int sensewire;
   TH1F *cell_dt[107];
   Int_t *bin_max;
   Int_t *bin_maxContent;
   Double_t *time_max;
   Double_t *twenty_perc_maxContent;
   Double_t *ref_time;



     //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;

   //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) {
     //declare array of histos to store drift times     
     total_wires=107; 

     //Declare bin properties for given sense wires in a plane

     bin_max = new Int_t[total_wires];                    /*Array to store the bin number corresponding to the drift time distribution peak*/
     bin_maxContent= new Int_t[total_wires];             /*Array to store the content (# events) corresponding to the bin with maximum content*/
     time_max= new Double_t[total_wires];                /*Array to store the x-axis(drift time (ns)) corresponding to bin_max*/
     twenty_perc_maxContent= new Double_t[total_wires];  /*Array to store 20% of maximum bin content (peak)*/						     
     ref_time= new Double_t[total_wires];               /*Array to store reference times for each sense wire*/

   }
   
   else if(ip == 2 || ip == 3 || ip == 8 || ip == 9) {
     total_wires=79;      
     bin_max = new Int_t[total_wires];                               
     bin_maxContent= new Int_t[total_wires]; 
     time_max= new Double_t[total_wires];
     twenty_perc_maxContent= new Double_t[total_wires];
     ref_time= new Double_t[total_wires];
   }	   
   
 	
   /*Get wire histos from root file and loop over each 
     sense wire of a plane in shms Drift Chambers (DC1 or DC2)*/
 
   for (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********//
     //*********************************************************//
     

     
     //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
     //cout << "y_int: " << y_int << " :: " << "m: " << m << " :: " << "t0: " << setprecision(6) << -y_int/m << endl;
     //Write "t0" values to file
     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();
     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();
   
   //*****************************************************************************************//
   //        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 + "/shms_dc_"+plane_names[ip]+Form("tzero_run%d_updated.dat", run_NUM);
   ofs.open(t_zero_file_corr);
   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
   double sum_DEN;   //denominator of weighted avg
   double weighted_AVG;
   double weighted_AVG_err; 
  
    int counter;
    // 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;
   
   weighted_AVG =0.0 ;
   weighted_AVG_err= 0.0; 
   
   counter = 0;
   //read line bt line the t_zero_file
   while(getline(ifs, line)) {

     if(line[0]!='#') {
       
       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 && abs(t0)<30) 
       {

	 //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;
	 
	 

       }
     
     else { ofs << sensewire << setw(15) << 0.0 << setw(15) << 0.0 << setw(15) << entries << endl;}
     
     //}
     
     } //end if statement
   }
   
   
   
   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.dat", 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 = "shms_dc"+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









 }
 
 




}