#include "TH1.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TF1.h"
#include "TPolyMarker.h"
#include "TRandom.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TExec.h"
#include "TPaveText.h"

#include <vector>
#include "time.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include "TFile.h"

//Make the fit function
Double_t fitfunction(Double_t *x, Double_t *par) { 
  return par[0]*(exp(0.5*(pow(((x[0]-par[1])/par[2]),2.0))))+par[3]; 
}

//run number for Clover_Energy_Spectra_Run_#.root file
//int clover_runNum = 28;
int clover_runNum = 300;


//define global things
TH1D *h1;
TFile *fin;
TDirectory *dir;
TCanvas *c1;
vector<double> Peak_Locations;
vector<double> Peak_Heights;
int x_up;
int x_down;

//Counters
int counter;
int marker_counter=0;
int fit_counter=0;

//Makers
TPolyMarker *pm[100];

//Fit Parameters
double Sigma[100];
double Sigma_Error[100];
double Mean[100];
double Mean_Error[100];
double Height[100];
double Height_Error[100];
double Area[100];
double Area_Error[100];
double Resolution[100];
double Resolution_Error[100];
double Efficiency[100];
double Efficiency_Error[100];
double Slope;
double Slope_Error;
double Intercept;
double Intercept_Error;

//Turn on and off testing
                                                                                                                                                                        bool Program_Test=false;
bool MCA_Data=false;

//Safeguard against saving under stupidity or without results
bool Dont_Save=true;

//Time stuff;
struct tm source_date;
struct tm run_date;

//Source Stuff
int N_Components;
double Half_Lives[100];  //in days

int N_Lines;
double Energy[100]; //in keV
double Energy_Unc[100]; //in keV
double CPS[100]; //CPS of the source
double CPS_Unc[100]; //in %
int Source_Comp[100]; 

double CPS_Decay_Corr[100]; // Decay corrected source intensities (CPS)
double CPS_Error_Decay_Corr[100]; //Decay corrected source intensity errors (CPS)
double Decay_Corr[100]; //factor to multiply by for decay corr

double Counts_Expected[100];
double Counts_Expected_Error[100];

//Run Stuff
double Run_Length; //in minutes
double Live_Time; //percentage

//Fit Functions
TF1 *f1[100];
TF1 *f2[100];

//Executive Function
void Position_Selector() {
  TObject *selected = gPad->GetSelected();
  if(!selected) return;
  if (!selected->InheritsFrom(TAxis::Class())) {
    gPad->SetUniqueID(0);
    return;
    
  }
  gPad->GetCanvas()->FeedbackMode(kTRUE);
  
  int event = gPad->GetEvent();
  int eventx = gPad->GetEventX();

  //cout << "event = " << event <<", eventx = " << eventx <<endl;
  
  if(event ==1) {
    x_down = eventx;
  }
   if(event == 11) {
    x_up = eventx;
  }

  if(event == 11 && (x_up == x_down)) {
    // if(event == 5 && (x_up == x_down)) {
    double px =0;
    double py=0;
    //look before and after the peak for the height (only for asthetics)
    //for(int k= gPad->PixeltoX(eventx)-25; k<gPad->PixeltoX(eventx)+25; k++) {
    for(int k= gPad->PixeltoX(eventx)-100; k<gPad->PixeltoX(eventx)+100; k++) {
      Int_t bin = h1->GetXaxis()->FindBin(k);
      //Double_t yp = h1->GetBinContent(bin);
      
       if(h1->GetBinContent(bin) > py) {
	py = h1->GetBinContent(k);
	px = k;
      }
    }
    py = 1.03*py;
    pm[marker_counter] = new TPolyMarker(1,&px,&py);
    pm[marker_counter]->SetMarkerStyle(20);
    pm[marker_counter]->SetMarkerColor(2);
    pm[marker_counter]->Draw();
    gPad->Update();

    //  Peak_Locations.push_back(gPad->PixeltoX(eventx));
    Peak_Locations.push_back(px);
    Peak_Heights.push_back(py);
    cout<< marker_counter <<"  "<<px<<endl;
    marker_counter++;
    
  }
  
  //if you press "s" you can view the peaks selected
  if(event == 24 && eventx == 115) {
    
    if(Dont_Save) {
      cout<<"I refuse to save! Check the efficency and energy calibration histograms"<<endl;
    }
    else {
     // printf("Saving to ''Clover_%d_Results.txt''",counter);
        printf("Saving to ''DSSD-B_%d_Results.txt''",counter);
      stringstream ss;
      ss.str();
     // string s1 = "Clover_";
    //  string s2 = "_Results.txt";
      
        string s1 = "DSSD-B_";
        string s2 = "_Results.txt";

        
      ss<<s1<<counter<<s2;
      ofstream out;
      out.open(ss.str().c_str());
      out<<Slope<<"  "<<Slope_Error<<"\n";
      out<<Intercept<<"  "<<Intercept_Error<<"\n";
      
      for(int n =0; n<N_Lines; n++) {
	out<<Energy[n]<<"   "<<Mean[n]<<"  "<<Mean_Error[n]<<"   "<<Area[n]<<"  "<<Area_Error[n]<<"   "<<Efficiency[n]<<"  "<<Efficiency_Error[n]<<"\n";
      }
      
      cout<<endl;
      for(int n = 0; n<Peak_Locations.size(); n++) {
	cout<<n<<"  "<<Peak_Locations[n]<<endl;
      }
    }
  }

  //if you press "f" you can fit all the peaks selected
  if(event == 24 && eventx == 102) {
    fit_counter=0;

    if(Peak_Locations.size() != N_Lines) {
      cout<<"Warning!!!!"<<endl;
      cout<<"The number of peaks selected from the spectrum doesn't match the source definition!"<<endl;
      cout<<"The results will be worthless.  This will not Save!"<<endl;
      Dont_Save=true;
    }
    else {
      Dont_Save=false;
    }
    
    
   
    //bubble sort so they in order of lowest to highest
    for(int m=0; m<Peak_Locations.size(); m++) {
      for(int n = 0; n<Peak_Locations.size()-1; n++) {
    	if(Peak_Locations[n+1]<Peak_Locations[n]) {
    	  cout<<"Ordering the locations appropriately"<<endl;
    	  double temp = Peak_Locations[n];
    	  Peak_Locations[n]=Peak_Locations[n+1];
    	  Peak_Locations[n+1]=temp;
	  double temp2 = Peak_Heights[n];
    	  Peak_Heights[n]=Peak_Heights[n+1];
    	  Peak_Heights[n+1]=temp2;
    	}
      }
    }

    for(int i=0; i<Peak_Locations.size(); i++)
      cout << "i " << i << " pos " << Peak_Locations[i] << " height " << Peak_Heights[i] << endl;
    getc(stdin);
    
    //Fit the histogram
    for(int n = 0; n<Peak_Locations.size(); n++) {
      f1[n] = new TF1(Form("Fit_%d",n),"gaus(0)+pol1(3)",(Peak_Locations[n]-220),(Peak_Locations[n]+220)); //200

      f1[n]->SetParameter(0,Peak_Heights[n]);
      f1[n]->SetParameter(1,Peak_Locations[n]);
      f1[n]->SetParameter(2,120); //5.0 //100
      f1[n]->SetParameter(3,Peak_Heights[n]*0.05);
      f1[n]->SetParameter(4,0.0);
      // f1[n]->SetParameter(5,1.0);

      // cout << " INITIAL FIT " << n << endl;
      h1->Fit(Form("Fit_%d",n),"SELRMNQ");
      // getc(stdin);

      int fit_low = 210;//20
      int fit_high = 220;//20

      // if(Energy[n]==42.8) {
      // 	fit_high= 5;
      // }          
      // if(Energy[n]==591.8) {
      // 	fit_low = 20;
      // }
      // if(Energy[n]==996.3) {
      // 	fit_high=9;
      // }
      // if(Energy[n]==1004.7) {
      // 	fit_low=10;
      // }
      // if(Energy[n]==123.1) {
      // 	fit_high=50;
      // 	fit_low=50;
      // }
      // if(Energy[n]==723.3) {
      // 	fit_high=25;
      // 	fit_low=25;
      // }
      
      f2[n] = new TF1(Form("Fit2_%d",n),"gaus(0)+pol1(3)",(Peak_Locations[n]-1.7*fit_low),(Peak_Locations[n]+1.7*fit_high));

      f2[n]->SetParameter(0,f1[n]->GetParameter(0));
      f2[n]->SetParameter(1,f1[n]->GetParameter(1));
      // f2[n]->SetParLimits(1,(Peak_Locations[n]-0.5*fit_low),(Peak_Locations[n]+0.5*fit_high));
      f2[n]->SetParameter(2,f1[n]->GetParameter(2));
      f2[n]->SetParameter(3,f1[n]->GetParameter(3));
      f2[n]->SetParameter(4,f1[n]->GetParameter(4));
      f2[n]->SetParameter(5,f1[n]->GetParameter(5));
		 
      h1->Fit(Form("Fit2_%d",n),"RMQ+");
      f2[n]->SetLineColor(2);
      f2[n]->Draw("Same");
      fit_counter++;
      gPad->Modified();
      gPad->Update();
      
      //Get the Parameters of the fit
      Sigma[n] = fabs(f2[n]->GetParameter(2));
      Mean[n] = f2[n]->GetParameter(1);
      Height[n] = f2[n]->GetParameter(0);
      
      //Get the associated Errors
      Sigma_Error[n] = f2[n]->GetParError(2);
      Mean_Error[n] = f2[n]->GetParError(1);
      Height_Error[n] = f2[n]->GetParError(0);
      
      //Calculate the Area
      Area[n] = pow((2.0*3.14152),0.5)*Height[n]*Sigma[n];
      Area_Error[n] = pow((2.0*3.14152),0.5)*pow(((pow(Sigma_Error[n],2.0)*pow(Height[n],2.0))+(pow(Height_Error[n],2.0)*pow(Sigma[n],2.0))),0.5);
      //Calculate % Resolution
      Resolution[n] = 100*2.35482*Sigma[n]/Mean[n];
      Resolution_Error[n] = 100*2.35482*pow(((pow(Sigma_Error[n],2.0)*pow((1.0/Mean[n]),2.0))+(pow(Mean_Error[n],2.0)*pow(-1.0*Sigma[n]/(pow(Mean[n],2.0)),2.0))),0.5);
      
      //Calculate the % Efficiency
      Efficiency[n] = 100.0*Area[n]/Counts_Expected[n];
      Efficiency_Error[n] = 100.0*pow(((pow(Area_Error[n],2.0)*pow((1.0/Counts_Expected[n]),2.0))+(pow(Counts_Expected_Error[n],2.0)*pow(-1.0*Area[n]/(pow(Counts_Expected[n],2.0)),2.0))),0.5);
      
      cout<<"Peak: "<<n<<endl;
      cout<<"Area: "<<Area[n]<<" +/- "<<Area_Error[n]<<endl;
      cout<<"Resolution: "<<Resolution[n]<<" +/- "<<Resolution_Error[n]<<endl;
      cout<<"Efficiency: "<<Efficiency[n]<<" +/- "<<Efficiency_Error[n]<<endl;
    }
    
    for(int n = 0; n<Peak_Locations.size(); n++) {
      cout<<n<<"  "<<Peak_Locations[n]<<endl;
    }

    h1->GetXaxis()->SetRangeUser(Peak_Locations[0]-200, Peak_Locations[Peak_Locations.size()-1]+200);
    gPad->Modified();
    gPad->Update();
    
    //Make energy calibrations
    TGraph *gr1 = new TGraphErrors(N_Lines, Mean, Energy, Mean_Error, Energy_Unc);
    gr1->SetName("En_Cal");
    gr1->SetTitle("Energy vs Channel");
    gr1->GetXaxis()->SetTitle("Channels (Arb. Units)");
    gr1->GetYaxis()->SetTitle("Energy (keV)");
    gr1->SetMarkerStyle(20);

    //Make efficiency calibrations
    TGraph *gr2 = new TGraphErrors(N_Lines, Energy, Efficiency, Energy_Unc, Efficiency_Error);
    gr2->SetName("Eff_Cal");
    gr2->SetTitle("Efficiency vs Energy");
    gr2->GetXaxis()->SetTitle("Energy (keV)");
    gr2->GetYaxis()->SetTitle("Efficiency (%)");
    gr2->SetMarkerStyle(20);

    //Make resolutions plot
    TGraph *gr3 = new TGraphErrors(N_Lines, Energy, Resolution, Energy_Unc, Resolution_Error);
    gr3->SetName("Res");
    gr3->SetTitle("Resolution vs Energy");
    gr3->GetXaxis()->SetTitle("Energy (keV)");
    gr3->GetYaxis()->SetTitle("Resolution (%)");
    gr3->SetMarkerStyle(20);
    
    
    
    TCanvas *c2 = new TCanvas("Results","Results",640,480);
    c2->Divide(3,1);
    
    c2->cd(1);
    gr1->Draw("AP");
    //TF1 *en_fit = new TF1("en_fit","pol1(0)",(Peak_Locations[3]),(Peak_Locations[Peak_Locations.size()-1]+200));

     TF1 *en_fit = new TF1("en_fit","pol1(0)",(Peak_Locations[0]-200),(Peak_Locations[Peak_Locations.size()-1]+200));
    en_fit->SetParameter(1,0.3);
    en_fit->SetParameter(0,0.2);
    gr1->Fit("en_fit","RMQ");
    gr1->Fit("en_fit","RMQ");
    gr1->Fit("en_fit","RMQ");

    Slope = en_fit->GetParameter(1);
    Intercept = en_fit->GetParameter(0);

    Slope_Error = en_fit->GetParError(1);
    Intercept_Error = en_fit->GetParError(0);

    cout<<"Energy Calibration"<<endl;
    cout<<"Slope: "<<Slope<<" +/- "<<Slope_Error<<endl;
    cout<<"Intercept: "<<Intercept<<" +/- "<<Intercept_Error<<endl;
    
    c2->cd(1);
    TPaveText *pten = new TPaveText(200,1400,3500,1600);
    pten->SetFillColor(0);
    pten->SetTextAlign(12);
    pten->AddText(Form("Slope: %3.5f +/- %3.5f",Slope, Slope_Error));
    pten->AddText(Form("Intercept: %3.5f +/- %3.5f",Intercept, Intercept_Error));
    pten->Draw("same");
    c2->Modified();
    c2->Update();
    
    c2->cd(2);
    gr2->Draw("AP");
    TPaveText *pteff = new TPaveText(400,2.5,1600,6.5);
    pteff->SetFillColor(0);
    pteff->SetTextAlign(32);
    pteff->AddText("Energy (keV)    Efficiency (%)");
    // for(i=3; i<N_Lines; i++) {
    for(int i=0; i<N_Lines; i++) {
      pteff->AddText(Form("%4.1f     %3.4f +/- %3.4f",Energy[i],Efficiency[i],Efficiency_Error[i]));      
    }
    pteff->Draw("same");
    c2->Modified();
    c2->Update();
    
    
    c2->cd(3);
    gr3->Draw("AP");

    TPaveText *ptres = new TPaveText(400,1,1600,5);
    ptres->SetFillColor(0);
    ptres->SetTextAlign(32);
    ptres->AddText("Energy (keV)    Resolution (%)");
    // for(i=3; i<N_Lines; i++) {
    for(int i=0; i<N_Lines; i++) {
      ptres->AddText(Form("%4.1f     %3.4f +/- %3.4f",Energy[i],Resolution[i],Resolution_Error[i]));      
    }
    ptres->Draw("same");
    c2->Modified();
    c2->Update();
    
    c2->Modified();
    c2->Update();      

  }

  //if you press "v" you can view the peaks selected
  if(event == 24 && eventx == 118) {
    
    //Bubble Sort
    for(int m=0; m<Peak_Locations.size(); m++) {
      for(int n = 0; n<Peak_Locations.size()-1; n++) {
	if(Peak_Locations[n+1]<Peak_Locations[n]) {
	  double temp = Peak_Locations[n];
	  Peak_Locations[n]=Peak_Locations[n+1];
	  Peak_Locations[n+1]=temp;
	}
      }
    }
    
    cout<<"In the terminal press ''n'' for next or ''q'' quit"<<endl;
    for(int n = 0; n<Peak_Locations.size()+1; n++) {
      
      if(n<Peak_Locations.size()) {
	h1->GetXaxis()->SetRangeUser(Peak_Locations[n]-200, Peak_Locations[n]+200);
	gPad->Modified();
	gPad->Update();
	
	cout<<"Peak Number: "<<n<<" Located at: "<<Peak_Locations[n]<<endl;
	
	string text;
	cin>>text;
	
	if(text=="n") continue;
	if(text=="q") {
	  h1->GetXaxis()->SetRangeUser(Peak_Locations[0]-200, Peak_Locations[Peak_Locations.size()-1]+200);
	  gPad->Modified();
	  gPad->Update();
	  return;
	}      
      }
      else {
	h1->GetXaxis()->SetRangeUser(Peak_Locations[0]-200, Peak_Locations[Peak_Locations.size()-1]+200);
	gPad->Modified();
	gPad->Update();
	return;
      }
    }
  }
  
  //pressing "r" removes the last peak
  if(event == 24 && eventx == 114) {
    if(Peak_Locations.size() > 0) {
      cout<<"Removing Last Point"<<endl;
      Peak_Locations.pop_back();
      
      marker_counter--;
      cout<<"Marker "<<marker_counter<<" delete" <<endl;
      TList *li = gPad->GetListOfPrimitives();
      li->Remove(pm[marker_counter]);
      //   pm[marker_counter]->Delete();
      gPad->Modified();
      gPad->Update();
    
    }
    else {
      cout<<"Cant remove points because there is nothing selected!"<<endl;
    }
  }
}


//Main function
void Clover_Analyzer() {
  
  /* cout<<"What Clover do you want to analyze"<<endl; */
  cout<<"Which DSSD strip (0-39) do you want to analyze"<<endl;
  cin>>counter;
  cout<<"Detector Number: "<<counter<<" Will be Analyzed"<<endl;
  TExec *ex1;
  
  
  /* if(Program_Test) { */
  /*     h1 = new TH1D(Form("Clover_Spectrum_%d",counter), Form("Clover_Spectrum_%d",counter), 31768, 0, 31768); */
    
  /*   //Make a fake histogram for testing purposes */
    
  /*   for(int j=0; j<32768; j++) { */
  /*     //  h1->Fill(j,gRandom->Uniform(1.0,5.0)); */
  /*   } */
  /*   for(int j=0; j<2000; j++) { */
  /*     h1->Fill(gRandom->Gaus(1500.0, 5.0)); */
  /*   } */
  /*   for(int j=0; j<1000; j++) { */
  /*     h1->Fill(gRandom->Gaus(5500.0,5.0)); */
  /*   } */
  /*   for(int j=0; j<1500; j++) { */
  /*     h1->Fill(gRandom->Gaus(10000.0, 5.0)); */
  /*   } */
  /*   for(int j=0; j<100000; j++) { */
  /*     h1->Fill(gRandom->Exp(2000.0)); */
  /*   } */
  /* } */
  /* else { */
  /*   if(MCA_Data) { */
  /*     h1 = new TH1D(Form("Clover_Spectrum_%d",counter), Form("Clover_Spectrum_%d",counter), 8192, 0, 8192); */
      
  /*     ifstream mca_data; */
  /*     mca_data.open("MCA_Clover_9_Ch1.txt"); */
  /*     int temp; */
  /*     for(int i=0; i<8192; i++) { */
  /* 	mca_data>>temp; */
  /* 	h1->Fill(i,temp); */
  /*     }       */
  /*   } */

  /*   else{ */
      //Get the histogram from the file
      // TFile *fin = new TFile(Form("Clover_Energy_Spectra_Run_%d.root", clover_runNum));
      // h1 = (TH1D*)fin->Get(Form("Clover_Spectrum_%d",counter));
      //  TFile *fin = new TFile(Form("DSSDFB_HiGain_Run_300.root"));
      TFile *fin =  new TFile( Form("Test-for-Chain-003.root") );
      //   cout<<"file opened"<<fin<<endl;
      // std::vector< TString > folder(0);
      cout <<"enetr"<<endl;
      //folder[0] = "SSSD/HiERaw";
      //cout<<"******"<<folder[0]<<endl;
      TDirectory* dir = fin ->GetDirectory("SSSD");
      cout<<"*****"<<dir<<endl;
      if(dir ==0)
	{
	  cout<<"Fatal error"<<dir<<"doesnt exsit"<<endl;
	  return;
	}
      
      dir->cd();
      TDirectory* dirSSSD = dir->GetDirectory("HiERaw/");
      dirSSSD->cd();
      
      TH1D *h1;
      TString name = Form("SSSD_Hi_%2d_ERaw",counter);
      dirSSSD->GetObject(name.Data(),h1);
      cout<<"******* H! "<<h1<<endl;
      //  h1->Draw();
      
      // TH1D* h1 = (TH1D*)fin->FindObject(Form("SSSD/HiERaw/SSSD_Hi_%2d_ERaw",counter));      
	//   h1= (TH1D*)fin->Get(Form("SSSD/HiERaw/SSSD_Hi_%2d_ERaw",counter));
      //  h1 = (TH1D*)fin->Get(Form("DSSD_HiGain_Front%d",counter));
      //	h1 = (TH1D*)fin->Get(Form("DSSD_HiGain_Back%d",counter));

	
	// h1 = (TH1D*)fin->Get(Form("Clover_Spectrum_%d",counter));
    


  // cout<<h1<<endl;
  
  //Read in the data from the source file
  cout<<"Reading in Source Data"<<endl;
  ifstream source_data;
  
  if(Program_Test) {
    // source_data.open(".SRM");
    source_data.open(".Program_Test");
  }
  else {
    //source_data.open(".Ra");
    //  source_data.open(".SRM");
     source_data.open(".228Th");
    //source_data.open(".Background");
  }
  
  //Read in the number of compenents and what the half life is (in days) for each
  source_data>>N_Components;
  for(int n = 0; n<N_Components; n++) {
    source_data>>Half_Lives[n];
  }
  
  //Read in the date of the source  Year Month Day Hour
  source_data>>source_date.tm_year>>source_date.tm_mon>>source_date.tm_mday>>source_date.tm_hour>>source_date.tm_min;
  source_date.tm_sec =0;
  
  //Read in the decay energies and intensities  Energy (kev) Unc (kev)  Intensity (CPS)  Unc (%)
  source_data>>N_Lines;
  for(int n = 0; n<N_Lines; n++) {
    source_data>>Energy[n]>>Energy_Unc[n]>>CPS[n]>>CPS_Unc[n]>>Source_Comp[n];
  }
  
  source_data.close();
  cout<<"Source found! "<<N_Components<<" Components with "<<N_Lines<<" Lines Read In Succesfully"<<endl;
  
  //Read in the data from the run file
  cout<<"Reading in Run Data (.Run_Data)"<<endl;
  
  ifstream run_data;
  run_data.open(".Run_Data");
  
  run_data>>Run_Length; //in minutes
  run_data>>Live_Time; //percentage
  
  //Read in the date of the run  Year Month Day Hour
  run_data>>run_date.tm_year>>run_date.tm_mon>>run_date.tm_mday>>run_date.tm_hour>>run_date.tm_min;
  run_date.tm_sec =0;

  cout.precision(10);
  //calucate the time difference (in seconds) between the source date and the run date
  double time_diff = difftime(mktime(&run_date),mktime(&source_date));
  cout<<"It has been "<<time_diff<<" seconds between the run and the source production date"<<endl;
  
  //Decay correct each component
  for(int n=0; n<N_Components; n++) {
    Decay_Corr[n] = exp(-1.0*time_diff*log(2)/(Half_Lives[n]*24*3600.0));
    cout<<"Component: "<<n<<" Decay Correction: "<<Decay_Corr[n]<<endl;
  }
  
  //correct for the decay 
  for(int n=0; n<N_Lines; n++) {
    CPS_Decay_Corr[n] = CPS[n]*Decay_Corr[Source_Comp[n]];
    CPS_Error_Decay_Corr[n] = (CPS_Unc[n]/100.0)*CPS[n]*Decay_Corr[Source_Comp[n]];
    cout<<CPS_Decay_Corr[n]<<" +/- "<<CPS_Error_Decay_Corr[n]<<endl;
  }

  //correct for and live time then multiply by the run time
  for(int n=0; n<N_Lines; n++) {
    Counts_Expected[n] = CPS_Decay_Corr[n]*(Live_Time/100.0)*Run_Length*60.0;
    Counts_Expected_Error[n] = CPS_Error_Decay_Corr[n]*(Live_Time/100.0)*Run_Length*60.0;
    //   cout<<"Line "<<n<<" at "<<Energy[n]<<" keV: "<<Counts_Expected[n]<<" +/- "<<Counts_Expected_Error[n]<<" Counts Expected"<<endl;
    }
  
  int duplicate_counter = 0;
  
  //check for any duplicate energies (i.e. xrays from mixed sources)
  for(int n=0; n<N_Lines-1; n++) {
    if(Energy[n] == Energy[n+1]) {
      Counts_Expected[n] = Counts_Expected[n]+Counts_Expected[n+1];
      Counts_Expected_Error[n] = pow((pow(Counts_Expected_Error[n],2.0)+pow(Counts_Expected_Error[n+1],2.0)),0.5); 
      duplicate_counter++;
      for(int m = n+1; m<N_Lines; m++) {
	Energy[m] = Energy[m+1];
	Counts_Expected[m] = Counts_Expected[m+1];
	Counts_Expected_Error[m] = Counts_Expected_Error[m+1];
      }
    }    
  }
  
  cout<<"Combined "<<duplicate_counter<<" Lines"<<endl;
  //correct the number of lines
  N_Lines = N_Lines-duplicate_counter;

  //output the results 
  for(int n=0; n<N_Lines; n++) {
    cout<<"Line "<<n<<" at "<<Energy[n]<<" keV: "<<Counts_Expected[n]<<" +/- "<<Counts_Expected_Error[n]<<" Counts Expected"<<endl;
  } 
  //End of calculations


  //Draw the histogram;  
  c1 = new TCanvas("Clover_Spectrum","Clover_Spectrum",1024,768);      
  //cout<<"Made Canvas "<<h1<<endl;
  c1->cd();
  // h1->Draw("hist");
  // cout<<"Draw Hist "<<h1<<endl;
  h1->Draw("hist");
  h1->GetXaxis()->SetTitle("Energy (Channels)");
  h1->GetYaxis()->SetTitle("Counts / Channel");
  h1->Rebin(2);
  gPad->AddExec(Form("Exec_%d",counter),"Position_Selector()");
  
  
}