CALIBRATION/hms_cal_calib/input.dat

+-10 10	Delta range, %
+0.5 1.5	Beta range
+1.5  	Gas Cherenkov, threshold on signals in p.e.
+10	Minimum number of hits per channel required to be calibrated
+0.05 2.	Range of uncalibrated energy deposition histogram
+500	Binning of uncalibrated energy deposition histogram
+-0.1 0.1	Gaussian fit range around mean e- peak in the uncalibrated Edep histogram 
+2.	Sigma width to select events to use in calibration (green fill)
+
+; Calibration constants for file hms_ecal.root (HMS 2489-2494), 164827 events processed
+
+hcal_pos_gain_cor= 12.99,  7.13,  9.22, 10.26, 10.82, 12.49, 11.93, 12.39, 10.28, 15.14, 14.59, 12.85,  6.40,
+                   10.87, 12.23,  9.46, 13.51,  8.71,  6.30,  7.73,  7.41,  8.94, 11.40, 11.51, 12.40, 14.86,
+                   24.29, 14.78, 18.17, 21.53, 16.96, 18.97, 23.00, 19.51, 21.92, 25.86, 18.37, 22.29, 19.84,
+                   32.17, 16.75, 20.16, 18.54, 17.89, 19.25, 22.14, 17.51, 19.67, 20.51, 18.37, 18.86, 26.46,
+hcal_neg_gain_cor= 15.11, 14.11, 14.42, 10.92, 11.88, 13.32, 14.58, 18.16, 12.63, 11.47, 10.60, 11.41, 19.57,
+                   14.71, 13.20, 14.22, 14.34, 16.11, 16.93, 19.64, 17.34, 18.06, 11.49, 14.83, 13.72,  8.86,
+                    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,
+

CALIBRATION/hms_cer_calib/README

+# Calibrating the Cherenkov
+* To automatically run the calibration
+```
+root -l "run_cal.C(RunNumber,NumEvents,COIN)"
+```
+* Where COIN == 1 indicates the file was a full coincident replay
+* If options are left blank, user will be prompted for a value
+* To manually run the scripts:
+* Link your ROOT file
+* Run the script with options (listed after script is executed)
+```
+root -l ../../ROOTfiles/hms_replay_488_-1.root
+T->Process("calibration.C+", "options");
+```
+* Or, if you want to run using PROOF
+```
+root -l
+TChain ch("T");
+ch.Add("/path/to/ROOTfile");
+TProof::Open(""); \\For PROOF-lite, adjust accordingly for PROOF
+ch.SetProof();
+ch.Process("calibration.C+", "options");
+```
+### Options for the calibration script are (case/spelling important):
+* **showall** - display all calibration details (lots of windows)

CALIBRATION/hms_cer_calib/calibration.C

+#define calibration_cxx
+// The class definition in calibration.h has been generated automatically
+// by the ROOT utility TTree::MakeSelector(). This class is derived
+// from the ROOT class TSelector. For more information on the TSelector
+// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
+
+
+// The following methods are defined in this file:
+//    Begin():        called every time a loop on the tree starts,
+//                    a convenient place to create your histograms.
+//    SlaveBegin():   called after Begin(), when on PROOF called only on the
+//                    slave servers.
+//    Process():      called for each event, in this function you decide what
+//                    to read and fill your histograms.
+//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
+//                    called only on the slave servers.
+//    Terminate():    called at the end of the loop on the tree,
+//                    a convenient place to draw/fit your histograms.
+//
+// To use this file, try the following session on your Tree T:
+//
+// root> T->Process("calibration.C")
+// root> T->Process("calibration.C","some options")
+// root> T->Process("calibration.C+")
+//
+
+
+#include "calibration.h"
+#include <TF1.h>
+#include <TH1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TSpectrum.h>
+#include <TPolyMarker.h>
+#include <TPaveStats.h>
+
+void calibration::Begin(TTree * /*tree*/)
+{
+   // The Begin() function is called at the start of the query.
+   // When running with PROOF Begin() is only called on the client.
+   // The tree argument is deprecated (on PROOF 0 is passed).
+
+   TString option = GetOption();
+   Info("Begin", "Starting calibration process with option: %s", option.Data());
+   Info("Begin", "To see details of calibration, use option showall");
+
+   if (option.Contains("showall")) fFullShow = kTRUE;
+}
+
+void calibration::SlaveBegin(TTree * /*tree*/)
+{
+  // The SlaveBegin() function is called after the Begin() function.
+  // When running with PROOF SlaveBegin() is called on each slave server.
+  // The tree argument is deprecated (on PROOF 0 is passed).
+
+  TString option = GetOption();
+  
+  fPulseInt = new TH1F*[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    fPulseInt[ipmt] = new TH1F(Form("PulseInt%d",ipmt+1), Form("Pulse integral PMT %d;ADC Channel (pC);Counts",ipmt+1), 1000, 0, 100);
+    GetOutputList()->Add(fPulseInt[ipmt]);
+  }
+  
+  fPulseInt_quad = new TH1F**[2];
+  fPulseInt_quad[0] = new TH1F*[4];
+  fPulseInt_quad[1] = new TH1F*[4];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    for (Int_t iquad = 0; iquad < 4; iquad++) {
+      fPulseInt_quad[ipmt][iquad] = new TH1F(Form("PulseInt_quad%d_PMT%d",iquad+1,ipmt+1),Form("Pulse Integral PMT%d quad%d; ADC Channel (pC); Counts",ipmt+1,iquad+1),1000,0,100);
+      GetOutputList()->Add(fPulseInt_quad[ipmt][iquad]);
+    }
+  }
+  
+  //Timing and Beta cut visualizations
+  fBeta_Cut = new TH1F("Beta_Cut", "Beta cut used for 'good' hits;Beta;Counts", 1000, -5, 5);
+  GetOutputList()->Add(fBeta_Cut);
+  fBeta_Full = new TH1F("Beta_Full", "Full beta for events;Beta;Counts", 1000, -5, 5);
+  GetOutputList()->Add(fBeta_Full);
+
+  fTiming_Cut = new TH1F*[2];
+  fTiming_Full = new TH1F*[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    fTiming_Cut[ipmt] = new TH1F(Form("Timing_Cut%d",ipmt+1), Form("Timing cut used for 'good' hits in PMT %d;Time (ns);Counts",ipmt+1), 10000, 50, 150);
+    GetOutputList()->Add(fTiming_Cut[ipmt]);
+    fTiming_Full[ipmt] = new TH1F(Form("Timing_Full%d",ipmt+1), Form("Full timing information for events in PMT %d;Time (ns);Counts",ipmt+1), 10000, 50, 150);
+    GetOutputList()->Add(fTiming_Full[ipmt]);
+  }
+
+  //Particle ID cut visualization
+  fCut_everything = new TH1F("Cut_everything", "Visualization of no cuts; Normalized Calorimeter Energy (GeV); Counts", 200, 0, 2.0);
+  GetOutputList()->Add(fCut_everything);
+  fCut_electron = new TH1F("Cut_electron", "Visualization of pion cut; Normalized Calorimeter Energy (GeV); Counts", 200, 0, 2.0);
+  GetOutputList()->Add(fCut_electron);
+}
+
+Bool_t calibration::Process(Long64_t entry)
+{
+   // The Process() function is called for each entry in the tree (or possibly
+   // keyed object in the case of PROOF) to be processed. The entry argument
+   // specifies which entry in the currently loaded tree is to be processed.
+   // When processing keyed objects with PROOF, the object is already loaded
+   // and is available via the fObject pointer.
+   //
+   // This function should contain the \"body\" of the analysis. It can contain
+   // simple or elaborate selection criteria, run algorithms on the data
+   // of the event and typically fill histograms.
+   //
+   // The processing can be stopped by calling Abort().
+   //
+   // Use fStatus to set the return value of TTree::Process().
+   //
+   // The return value is currently not used.
+
+   fReader.SetEntry(entry);
+   
+   //Only one track
+   if (*Ndata_H_tr_beta != 1) return kTRUE;
+
+   for (Int_t itrack = 0; itrack < *Ndata_H_tr_beta; itrack++) {
+     //Beta Cut
+     fBeta_Full->Fill(H_tr_beta[itrack]);
+     if (TMath::Abs(H_tr_beta[itrack] -1.0) > 0.2) return kTRUE;
+     fBeta_Cut->Fill(H_tr_beta[itrack]);
+
+     for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+       //Timing Cut
+       fTiming_Full[ipmt]->Fill(H_cer_goodAdcTdcDiffTime[ipmt]);
+       if (H_cer_goodAdcTdcDiffTime[ipmt] > 110 || H_cer_goodAdcTdcDiffTime[ipmt] < 90) continue;
+       fTiming_Cut[ipmt]->Fill(H_cer_goodAdcTdcDiffTime[ipmt]);
+
+       if (ipmt == 0) fCut_everything->Fill(*H_cal_etotnorm);
+       //Electron Cut
+       if (*H_cal_etotnorm < 0.5) {
+	 if (ipmt == 0) fCut_electron->Fill(*H_cal_etotnorm);
+	 fPulseInt[ipmt]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+
+	 //Quadrant Cuts
+	 if (*H_cer_xAtCer >= 0.0 && *H_cer_yAtCer >= 0.0) {
+	   fPulseInt_quad[ipmt][0]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer >= 0.0 && *H_cer_yAtCer < 0.0) {
+	   fPulseInt_quad[ipmt][1]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer < 0.0 && *H_cer_yAtCer >= 0.0) {
+	   fPulseInt_quad[ipmt][2]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+	 if (*H_cer_xAtCer < 0.0 && *H_cer_yAtCer < 0.0) {
+	   fPulseInt_quad[ipmt][3]->Fill(H_cer_goodAdcPulseInt[ipmt]);
+	 }
+       }//End of electron cut
+     }//End of PMT loop
+   }//End of tracking loop
+   
+   return kTRUE;
+}
+
+void calibration::SlaveTerminate()
+{
+   // The SlaveTerminate() function is called after all entries or objects
+   // have been processed. When running with PROOF SlaveTerminate() is called
+   // on each slave server.
+
+}
+
+void calibration::Terminate()
+{
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+
+  Info("Terminate", "'%s' showing", (fFullShow ? "full" : "minimal"));
+  Info("Terminate", "Histograms formed, now starting calibration.\n'Peak Buffer full' is a good warning!\n");
+  printf("\n");
+
+  //Have to extract the histograms from the OutputList 
+  TH1F* PulseInt[2];
+  TH1F* PulseInt_quad[2][4];
+  TH1F* Timing_Cut[2];
+  TH1F* Timing_Full[2];
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    Timing_Cut[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("Timing_Cut%d",ipmt+1)));
+    Timing_Full[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("Timing_Full%d",ipmt+1)));
+    PulseInt[ipmt] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("PulseInt%d",ipmt+1)));
+    for (Int_t iquad = 0; iquad < 4; iquad++) {
+      PulseInt_quad[ipmt][iquad] = dynamic_cast<TH1F*> (GetOutputList()->FindObject(Form("PulseInt_quad%d_PMT%d",iquad+1,ipmt+1)));
+    }
+  }
+
+  Double_t Cer_Peak[2];
+  //Begin peak Finding
+  for (Int_t ipmt = 0; ipmt < 2; ipmt++) {
+    PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,8);
+    TSpectrum *s = new TSpectrum(1);
+    s->Search(PulseInt[ipmt], 1.0, "nobackground&&nodraw", 0.001);
+    TList *functions = PulseInt[ipmt]->GetListOfFunctions();
+    TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
+    Cer_Peak[ipmt] = *pm->GetX();
+    PulseInt[ipmt]->GetXaxis()->SetRangeUser(0,100);
+  }
+
+  //Begin Fitting of SPE
+  /*TF1 *poisson =  new TF1("poisson","[1]*((pow([0],x)*exp(-[0]))/(tgamma(x+1)))",0,20);
+  poisson->SetParName(0,"#mu");
+  poisson->SetParName(1,"Amplitude");
+  poisson->SetParameters(5.0,100);
+  poisson->SetParLimits(0,0.0,10.0);
+  poisson->SetParLimits(1,0.0,9999.0);
+  poisson->SetRange(0,Cer_Peak[0]+2); PulseInt[0]->Fit("poisson","RQM");
+  TF1 *PulseIntPMT1Fit = PulseInt[0]->GetFunction("poisson");
+  poisson->SetRange(0,Cer_Peak[1]+2); PulseInt[1]->Fit("poisson","RQM");
+  TF1 *PulseIntPMT2Fit = PulseInt[1]->GetFunction("poisson");*/
+ 
+  TF1 *gaussian =  new TF1("gaussian","[1]*exp(-0.5*((x-[0])/[2])*((x-[0])/[2]))",0,20);
+  gaussian->SetParName(0,"#mu");
+  gaussian->SetParName(1,"Amplitude");
+  gaussian->SetParName(2,"#sigma");
+  gaussian->SetParameters(5.0,100,1.0);
+  gaussian->SetParLimits(0,0.0,10.0);
+  gaussian->SetParLimits(1,0.0,9999.0);
+  gaussian->SetParLimits(2,0.0,10.0);
+  gaussian->SetRange(0,Cer_Peak[0]+2); PulseInt[0]->Fit("gaussian","RQM0");
+  TF1 *PulseIntPMT1Fit = PulseInt[0]->GetFunction("gaussian");
+  gaussian->SetRange(0,Cer_Peak[1]+2); PulseInt[1]->Fit("gaussian","RQM0");
+  TF1 *PulseIntPMT2Fit = PulseInt[1]->GetFunction("gaussian");
+  
+  if (fFullShow) {
+    TCanvas *Beta = new TCanvas("Beta","Beta cuts used");
+    Beta->Divide(2,1);
+    Beta->cd(1); fBeta_Full->Draw();
+    Beta->cd(2); fBeta_Cut->Draw();
+   
+    TCanvas *Cal = new TCanvas("Cal","Calorimeter cuts used");
+    Cal->Divide(2,1);
+    Cal->cd(1); fCut_everything->Draw();
+    Cal->cd(2); fCut_electron->Draw();
+
+    TCanvas *Timing = new TCanvas("Timing","Timing cuts used");
+    Timing->Divide(2,2);
+    Timing->cd(1); Timing_Full[0]->Draw();
+    Timing->cd(2); Timing_Cut[0]->Draw();
+    Timing->cd(3); Timing_Full[1]->Draw();
+    Timing->cd(4); Timing_Cut[1]->Draw();
+    /*
+    TCanvas *PulseIntPMT1 = new TCanvas("PulseIntPMT1","Good Pulse Integral from Cherenkov PMT 1");
+    PulseIntPMT1->Divide(2,2);
+    PulseIntPMT1->cd(1); PulseInt_quad[0][0]->Draw();
+    PulseIntPMT1->cd(2); PulseInt_quad[0][1]->Draw();
+    PulseIntPMT1->cd(3); PulseInt_quad[0][2]->Draw();
+    PulseIntPMT1->cd(4); PulseInt_quad[0][3]->Draw();
+
+    TCanvas *PulseIntPMT2 = new TCanvas("PulseIntPMT2","Good Pulse Integral from Cherenkov PMT 2");
+    PulseIntPMT2->Divide(2,2);
+    PulseIntPMT2->cd(1); PulseInt_quad[1][0]->Draw();
+    PulseIntPMT2->cd(2); PulseInt_quad[1][1]->Draw();
+    PulseIntPMT2->cd(3); PulseInt_quad[1][2]->Draw();
+    PulseIntPMT2->cd(4); PulseInt_quad[1][3]->Draw();
+    */
+    gStyle->SetOptStat(111);
+    gStyle->SetOptFit(0111);
+    TCanvas *cPulseInt = new TCanvas("cPulseInt","Good Pulse Integral in both Cherenkov PMTs");
+    cPulseInt->Divide(2,1);
+    cPulseInt->cd(1); PulseInt[0]->Draw(); PulseIntPMT1Fit->Draw("same");
+    cPulseInt->cd(2); PulseInt[1]->Draw(); PulseIntPMT2Fit->Draw("same");
+  }
+
+  //Output the actual calibration information
+  cout << "Calibration constants are \nPMT#: ADC-to-NPE Value\n" << endl;
+  cout << "PMT 1:" << setw(8) << Form("%3.3f\n",PulseIntPMT1Fit->GetParameter(0));
+  cout << "PMT 2:" << setw(8) << Form("%3.3f\n",PulseIntPMT2Fit->GetParameter(0));
+
+  //Start the process of writing the calibration information to file
+  ofstream calibration;
+  calibration.open("calibration_temp.txt", ios::out);
+
+  if (!calibration.is_open()) cout << "Problem saving calibration constants, may have to update constants manually!" << endl;
+
+  else {
+    calibration << "hcer_adc_to_npe = "; 
+    calibration << Form("1./%3.3f, 1./%3.3f", PulseIntPMT1Fit->GetParameter(0), PulseIntPMT2Fit->GetParameter(0));
+    calibration.close();
+  }
+}

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_dc_calib/scripts/get_hdc_time_histo.C

@@ -12,17 +12,8 @@ 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/test_%d.root", run_NUM), "READ");
+TFile *f = new TFile(Form("../../../ROOTfiles/hms_replay_%d.root", run_NUM), "READ");
 
 //create new file
  TFile *g = new TFile(Form("../root_files/run%d/hms_dc_time_%d.root", run_NUM, run_NUM), "RECREATE"); // create new file to store histo
@@ -38,58 +29,58 @@ TString plane_names[NPLANES]={"1x1", "1y1", "1u1", "1v1", "1y2", "1x2", "2x1", "
 //Declare Variables to Loop Over
 Int_t Ndata[NPLANES];
 Double_t hdc_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 = "hdc"+plane_names[ip]+"_time"; 
+
+  TString drift_time_histo = "hdc"+plane_names[ip]+"_time";
   TString title = "hdc"+plane_names[ip]+"_drifttime";
-     
+
   //Set Branch Address
-  tree->SetBranchAddress(drift_time, &hdc_time[ip][0]);   
+  tree->SetBranchAddress(drift_time, &hdc_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 
+
+
+     //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(hdc_time[ip][j]); 
+
+	 h[ip]->Fill(hdc_time[ip][j]);
        }
-       
+
      }
 
    }
 
 
 
-		
+
  //Write histograms to file
  g->Write();
 

CALIBRATION/hms_hodo_calib/.gitignore

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

CALIBRATION/hms_hodo_calib/.rootrc

+Rint.History: .root_history

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,