{
gROOT->Reset();

  // Edit for appropriate filename
  TFile* fh = new TFile("run00568.root","read");
  TTree* sourceTree = (TTree*)fh->Get("Trigger");

  // Sourced values are left as integers, calculated to floats
  Int_t TDCS[8], ADCS[12];

  // Define data sources from the loaded tree
  sourceTree->SetBranchAddress("TDCS",TDCS);
  sourceTree->SetBranchAddress("ADCS",ADCS);

  gStyle->SetOptStat(11);
  gStyle->SetOptFit(1);
  TCut selectCut;
  char cutString[512];
  Double_t lowEdge, highEdge;
  Double_t spread1 = 440; //total spread of ADC1 divisions
  Double_t spread2 = 1000; //total spread of ADC2 divisions 
  const Int_t numPts = 14;
  const Double_t peakEnergy = 2; // MeV
  Double_t energy1[numPts], energy1_e[numPts], timeRes1[numPts], timeRes1_e[numPts], energy2[numPts], energy2_e[numPts], timeRes2[numPts], timeRes2_e[numPts];
  Double_t corrTimeRes1[numPts], corrTimeRes1_e[numPts], corrTimeRes2[numPts], corrTimeRes2_e[numPts];
  Double_t pedBin1, peakBin1, binEnergy1, pedBin2, peakBin2, binEnergy2;
  Int_t maxBin1, maxBin2;
  

  // Find the pedestal and max of /////// ADC1 ///////
  cout << "Finding energy calibration for ADC1..." << endl;
  sourceTree->Draw("ADCS[1] >> pedADC1","(ADCS[1])<50");
  pedADC1->Fit("gaus","Q", "", 0, 50);
  pedBin1 = pedADC1->GetFunction("gaus")->GetParameter(1);
  sourceTree->Draw("ADCS[1] >> peakADC1","(ADCS[1])>50");
  maxBin1 = peakADC1->GetMaximumBin();
  peakBin1 = peakADC1->GetBinLowEdge(maxBin1);
  peakADC1->Fit("gaus","Q","",peakBin1-100,peakBin1+100);
  peakBin1 = peakADC1->GetFunction("gaus")->GetParameter(1);

  binEnergy1 = peakEnergy / (peakBin1 - pedBin1);
  cout << "Found pedestal at " << pedBin1 << " and peak at " << peakBin1 << endl;
  cout << "Calculated " << binEnergy1 << " MeV/bin." << endl;

  // Find the pedestal and max of ////// ADC2 ///////
  cout << "Finding energy calibration for ADC2..." << endl;
  sourceTree->Draw("ADCS[2] >> pedADC2","(ADCS[2])<280");
  pedADC2->Fit("gaus","Q", "", 200, 280);
  pedBin2 = pedADC2->GetFunction("gaus")->GetParameter(1);
  sourceTree->Draw("ADCS[2] >> peakADC2","(ADCS[2])>280");
  maxBin2 = peakADC2->GetMaximumBin();
  peakBin2 = peakADC2->GetBinLowEdge(maxBin2);
  peakADC2->Fit("gaus","Q","",peakBin2-100,peakBin2+100);
  peakBin2 = peakADC2->GetFunction("gaus")->GetParameter(1);

  binEnergy2 = peakEnergy / (peakBin2 - pedBin2);
  cout << "Found pedestal at " << pedBin2 << " and peak at " << peakBin2 << endl;
  cout << "Calculated " << binEnergy2 << " MeV/bin." << endl;
   
  ///////////////////////////////////
  //Corrections made - find Energy and Time Resolution data points
  ///////////////////////////////////

  for (Int_t i=0; i<numPts; i++) {
    cout << "Finding ADC1 :"  << i << endl; 
    lowEdge = (spread1 - pedBin1)*(i)/(numPts) + pedBin1 + spread1/40;
    highEdge = (spread1 - pedBin1)*(i+1)/(numPts) + pedBin1 + spread1/40;
    sprintf(cutString,"(ADCS[1])>%f && (ADCS[1])<%f",lowEdge,highEdge);
    selectCut = cutString;
    sourceTree->Draw("ADCS[1] >> sampleHist",selectCut);    
    energy1[i] = ((sampleHist->GetMean() - pedBin1) * binEnergy1)/1000;
    energy1_e[i] = 0;
    sourceTree->Draw("TDCS[0] >> timeHist", selectCut);
    timeHist->Fit("gaus","Q", "", 0, 4000);
    timeRes1[i] = timeHist->GetFunction("gaus")->GetParameter(2)* 50;
    timeRes1_e[i] = timeHist->GetFunction("gaus")->GetParError(2)* 50;

    //cout << energy1[i] <<","<< timeRes1[i] <<","<< timeRes1_e[i] <<","<< corrTimeRes1[i] <<","<< corrTimeRes1_e[i] << endl;
  }

  for (Int_t i=0; i<numPts; i++) {
    cout << "Finding ADC2 :"  << i << endl; 
    lowEdge = (spread2 - pedBin2)*(i)/(numPts) + pedBin2 + spread2/40;
    highEdge = (spread2 - pedBin2)*(i+1)/(numPts) + pedBin2 + spread2/40;
    sprintf(cutString,"(ADCS[2])>%f && (ADCS[2])<%f",lowEdge,highEdge);
    selectCut = cutString;
    sourceTree->Draw("ADCS[2] >> sampleHist",selectCut);    
    energy2[i] = ((sampleHist->GetMean() - pedBin2) * binEnergy2)/1000;
    energy2_e[i] = 0;
    sourceTree->Draw("TDCS[1] >> timeHist", selectCut);
    timeHist->Fit("gaus","Q", "", 0, 4000);
    timeRes2[i] = timeHist->GetFunction("gaus")->GetParameter(2)* 50;
    timeRes2_e[i] = timeHist->GetFunction("gaus")->GetParError(2)* 50;

    //cout << energy2[i] <<","<< timeRes2[i] <<","<< timeRes2_e[i] <<","<< corrTimeRes2[i] <<","<< corrTimeRes2_e[i] << endl;
  }

  ///////////////////////////////////
  //Graph the calculated Time Resolution vs. Energy
  ///////////////////////////////////

  gROOT->LoadMacro("TimeResEnergy.cxx");

  //start segment

  //TF1* t2func = new TF1("TimeResEnergy",0,10,3);
  TF1* t2func = new TF1("t2func",TimeResEnergy,0,2,3);
  t2func->SetParName(0,"k");
  t2func->SetParName(1,"E_{0}");
  t2func->SetParName(2,"T");
  t2func->SetParameters(68,0.0001,1300);

  TCanvas* timeCan2 = new TCanvas("timeCan2","Time Resolution - ADC2");
  timeCan2->cd();
  TGraphErrors* timeGraph2 = new TGraphErrors(numPts,energy2,timeRes2,energy2_e,timeRes2_e);
  timeGraph2->SetTitle("Uncorrected Time Resolution vs. Energy - ADC2");
  timeGraph2->SetMarkerStyle(21);
  timeGraph2->SetMarkerColor(4);
  timeGraph2->GetXaxis()->SetTitle("E (GeV)");
  timeGraph2->GetYaxis()->SetTitle("#deltaT (pS)");
  timeGraph2->GetYaxis()->SetRangeUser(700, 5000);
  timeGraph2->Fit(t2func, "V", "",0,0.1);
  timeGraph2->Draw("AP");
  //timeGraph2->GetFunction("t2func")->Draw("same");
  timeCan2->Draw();

  TCanvas* timeCan1 = new TCanvas("timeCan1","Time Resolution - ADC1");
  timeCan1->cd();
  //TF1* t1func = new TF1("TimeResEnergy",0,10,3);
  TF1* t1func = new TF1("t1func",TimeResEnergy,0,2,3);
  t1func->SetParName(0,"k");
  t1func->SetParName(1,"E_{0}");
  t1func->SetParName(2,"T");
  t1func->SetParameters(68,2e-5,1000);
  TGraphErrors* timeGraph1 = new TGraphErrors(numPts,energy1,timeRes1,energy1_e,timeRes1_e);
  timeGraph1->SetTitle("Uncorrected Time Resolution vs. Energy - ADC1");
  timeGraph1->SetMarkerStyle(21);
  timeGraph1->SetMarkerColor(4);
  timeGraph1->GetXaxis()->SetTitle("E (GeV)");
  timeGraph1->GetYaxis()->SetTitle("#deltaT (pS)");
  timeGraph1->GetYaxis()->SetRangeUser(700, 5000);
  timeGraph1->Fit(t1func, "V", "",0,0.1);
  timeGraph1->Draw("AP");
  timeCan1->Draw();
 


}