#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string>
#include <string.h>
#include <exception>


#include "TTree.h"
#include "TFile.h"
#include "TString.h"

#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TChain.h"
#include "TEllipse.h"
#include "TLine.h"
#include "TMarker.h"
#include "TVector3.h"
#include "TRotation.h"
#include "TMath.h"
#include "TROOT.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TLatex.h" 
#include "TPolyLine.h" 


#define DEUTERON_MASS 1.875614
#define PROTON_MASS 0.93827
#define RAD2DEG 57.295779
#define DEG2RAD 0.017453293

using namespace std;

void HallA_style() {

  gROOT->SetStyle("Plain");
  gStyle->SetPaperSize(TStyle::kUSLetter);
  gStyle->SetPaperSize(18,22);
  gStyle->SetOptFit(1111);
  gStyle->SetPalette(1);
  gStyle->SetNdivisions(508);
  
  gStyle->SetCanvasColor(10);
  gStyle->SetPadTopMargin(.05);
  gStyle->SetPadLeftMargin(.15);
  gStyle->SetPadRightMargin(.1);
  gStyle->SetPadBottomMargin(.10);
  gStyle->SetTitleYOffset(1.7);
  gStyle->SetTitleXOffset(1.2);
  gStyle->SetLabelFont(42,"X");
  gStyle->SetLabelFont(42,"Y");

  gStyle->SetTitleX(0.2);
  gStyle->SetTitleW(0.7);
  gStyle->SetTitleFillColor(18);

  // prepare gStyle to be useful
  //   1 = solid
  //   2 = long dash (30 10)
  //   3 = dotted (4 8)
  //   4 = dash-dot (15 12 4 12)
  //   5 = short dash ( 15 15 )
  //   6 = dash-dot-dot 
  gStyle->SetLineStyleString(1,"[]");
  gStyle->SetLineStyleString(2,"[30 10]");
  gStyle->SetLineStyleString(3,"[4 8]");
  gStyle->SetLineStyleString(4,"[15 12 4 12]");
  gStyle->SetLineStyleString(5,"[15 15]");
  gStyle->SetLineStyleString(6,"[15 12 4 12 4 12]");
  gStyle->SetLabelSize(0.04,"X");
  gStyle->SetLabelSize(0.04,"Y");
  gStyle->SetNdivisions(508,"Y");
  gStyle->SetOptDate(0);
  gStyle->SetDateY(.98);
  gStyle->SetStripDecimals(kFALSE);


  gStyle->SetOptStat(1111);

}



//void BeamChargeAsymmetry(int runnumber = 2268)
//void BeamChargeAsymmetry(int runnumber = 1338)
//void BeamChargeAsymmetry(int runnumber = 1339)
//void BeamChargeAsymmetry(int runnumber = 2300)
//void BeamChargeAsymmetry(int runnumber = 2814)
//void BeamChargeAsymmetry(int runnumber = 2618)
void BeamChargeAsymmetry(int runnumber = 3155)
{

  HallA_style();


//####################################################################################
// // Calibration of BCMs.  Runs 1338, 1339 E05-102 Apr 2009
// // Parameters determined by Yawei.

  double calib_u1 = 2067.7;  double calib_u3 = 6370.86 ;  double calib_u10 = 19394.0; 
  double calib_d1 = 2172.32; double calib_d3 = 6712.81;   double calib_d10 = 21036.1; 
  double off_u1 = 396.056;   double off_u3 = 490.003;     double off_u10 = 930.143;  
  double off_d1 = 170.152;   double off_d3 = 211.956;     double off_d10 = 548.635;

//####################################################################################


  //--- Module variables
  Double_t DL_t1[100];
  Double_t DL_t2[100];
  Double_t DL_t3[100];
  Double_t DL_t5[100];
  Double_t DL_t8[100];
  Double_t DL_bbretime[100];
  Int_t Ndata_DL_t3;
  Int_t Ndata_DL_t1;
  Int_t Ndata_DL_t2;
  Int_t Ndata_DL_t5;
  Int_t Ndata_DL_t8;
  Int_t Ndata_bbretime;
  Int_t Ndata_DL_edtmbb;
  Double_t DL_evtypebits;
  Double_t g0hel_L_helicity;

  Double_t N_bcmu3;
  Double_t N_bcmu30;
  Double_t N_clocktime;
  Double_t N_PS1, N_PS2, N_PS3, N_PS4, N_PS5, N_PS6, N_PS7, N_PS8; 


  Int_t Count_bbretime = 0;
  Int_t Count_t3 = 0;


  //--- Scaler variables
  Double_t bcmu3count;
  Double_t bcmu3rate;
  Double_t G1_bcmu3count;
  Double_t G1_bcmu3rate;
  Double_t G2_bcmu3count;
  Double_t G2_bcmu3rate;
  Double_t G3_bcmu3count;
  Double_t G3_bcmu3rate;
  Double_t G4_bcmu3count;
  Double_t G4_bcmu3rate;

  Double_t clock1024rate;
  Double_t clock1024count;
  Double_t IBC0L02Current;
  Double_t hac_bcm_average;

  Double_t evbbite_BBretimerate;
  Double_t evbbite_BBretimecount;


  TChain *t1 =new TChain("T");


  //TString filename = Form("./eed_BigBiteScalar_2300.root");
  //TString filename = Form("./eed_BigBiteScalar_1338.root");
  TString filename = Form("./eed_BigBiteScalar_%d.root", runnumber);
  t1->Add(filename.Data());

  assert(t1);

//############################ VARIABLES ############################



  t1->SetBranchAddress( "DL.evtypebits", &DL_evtypebits);
  t1->SetBranchAddress( "DL.t1", &DL_t1);
  t1->SetBranchAddress( "DL.t2", &DL_t2);
  t1->SetBranchAddress( "DL.t3", &DL_t3);
  t1->SetBranchAddress( "DL.t5", &DL_t5);
  t1->SetBranchAddress( "Ndata.DL.t1", &Ndata_DL_t1);
  t1->SetBranchAddress( "Ndata.DL.t2", &Ndata_DL_t2);
  t1->SetBranchAddress( "Ndata.DL.t3", &Ndata_DL_t3);
  t1->SetBranchAddress( "Ndata.DL.t5", &Ndata_DL_t5);
  t1->SetBranchAddress( "Ndata.DL.edtmbb", &Ndata_DL_edtmbb);
  t1->SetBranchAddress( "adchel.L.helicity", &g0hel_L_helicity);

  t1->SetBranchAddress( "N.bcmu3", &N_bcmu3);
  t1->SetBranchAddress( "N.bcmu30", &N_bcmu30);
  t1->SetBranchAddress( "N.clocktime", &N_clocktime);
  t1->SetBranchAddress( "N.PS1", &N_PS1);
  t1->SetBranchAddress( "N.PS2", &N_PS2);
  t1->SetBranchAddress( "N.PS3", &N_PS3);
  t1->SetBranchAddress( "N.PS4", &N_PS4);
  t1->SetBranchAddress( "N.PS5", &N_PS5);
  t1->SetBranchAddress( "N.PS6", &N_PS6);
  t1->SetBranchAddress( "N.PS7", &N_PS7);
  t1->SetBranchAddress( "N.PS8", &N_PS8);


  t1->SetBranchAddress( "evleft_bcmu3cnt", &bcmu3count);
  t1->SetBranchAddress( "evleft_bcmu3", &bcmu3rate);
  t1->SetBranchAddress( "evleft_G1_bcmu3cnt", &G1_bcmu3count);
  t1->SetBranchAddress( "evleft_G1_bcmu3", &G1_bcmu3rate);
  t1->SetBranchAddress( "evleft_G2_bcmu3cnt", &G2_bcmu3count);
  t1->SetBranchAddress( "evleft_G2_bcmu3", &G2_bcmu3rate);
  t1->SetBranchAddress( "evleft_G3_bcmu3cnt", &G3_bcmu3count);
  t1->SetBranchAddress( "evleft_G3_bcmu3", &G3_bcmu3rate);
  t1->SetBranchAddress( "evleft_G4_bcmu3cnt", &G4_bcmu3count);
  t1->SetBranchAddress( "evleft_G4_bcmu3", &G4_bcmu3rate);


  t1->SetBranchAddress( "evleft_lclk", &clock1024rate);
  t1->SetBranchAddress( "evleft_lclkcount", &clock1024count);
  //t1->SetBranchAddress( "evleft_clk104k", &clock1024rate);
  //t1->SetBranchAddress( "evleft_clk104kcnt", &clock1024count);

  t1->SetBranchAddress( "IBC0L02Current", &IBC0L02Current);
  t1->SetBranchAddress( "hac_bcm_average", &hac_bcm_average);

  t1->SetBranchAddress( "evbbite_BBretime", &evbbite_BBretimerate);
  t1->SetBranchAddress( "evbbite_BBretimecnt", &evbbite_BBretimecount);
  



//############################ Histograms ############################

  TH1F *hBeamHelicityPlus = new TH1F("hBeamHelicityPlus","hBeamHelicityPlus",10,-2,2);
  TH1F *hBeamHelicityMinus = new TH1F("hBeamHelicityMinus","hBeamHelicityMinus",10,-2,2);

  TGraph *grBCMu3cnt_vs_time = new TGraph();
  TH2F *hBCMu3cnt_vs_time = new TH2F("hBCMu3cnt_vs_time", "hBCMu3cnt_vs_time", 100, 0, 2000, 100, 0, 2E8);

  TGraph *grNA_BCMu3cnt_vs_time = new TGraph();
  TH2F *hNA_BCMu3cnt_vs_time = new TH2F("hNA_BCMu3cnt_vs_time", "hNA_BCMu3cnt_vs_time", 100, 0, 2000, 100, 0, 2E8);

  TH1F *hTimeDiff = new TH1F("hTimeDiff","hTimeDiff",100,-0.01,1);

  TH1F *hCurrentDistribution = new TH1F("hCurrentDistribution","hCurrentDistribution",200,-4,4);

  TH1F *hCollectChargeDifference = new TH1F("hCollectChargeDifference","hCollectChargeDifference",200,-10.,10);
  TH1F *hCollectChargeDifferenceScalers = new TH1F("hCollectChargeDifferenceScalers","hCollectChargeDifferenceScalers",200,-10.,10);

  TGraph *grCharge_vs_time = new TGraph();
  TGraph *grCollectedCharge_vs_time = new TGraph();
  TGraph *grCollectedChargePlus_vs_time = new TGraph();
  TGraph *grCollectedChargeMinus_vs_time = new TGraph();

  TGraph *grCollectedChargeScalers_vs_time = new TGraph();
  TGraph *grCollectedChargePlusScalers_vs_time = new TGraph();
  TGraph *grCollectedChargeMinusScalers_vs_time = new TGraph();


  TGraph *grCollectedChargeAsymmetry_vs_time = new TGraph();
  TGraph *grCollectedChargeAsymmetryScalers_vs_time = new TGraph();

  TGraph *grCurrent_vs_time = new TGraph();
  TGraph *grEpicsCurrent_vs_time = new TGraph();
  TGraph *grEpicsCurrentHacMean_vs_time = new TGraph();
  

//############################ Analysis ############################

  int num = t1->GetEntries();
  //num = 30000;
  cout<<"Number of entries : "<<num<<endl;

  double TotalCharge = 0.0;

  double time = 0.0;
  double oldtime = 0.0;
  double CollectedCharge = 0.0;
  double CollectedChargePlus = 0.0;
  double CollectedChargeMinus = 0.0;
  double ChargeAsymmetry = 0.0;
  double ChargeAsymmetryError = 0.0;

  double ScalersCollectedCharge = 0.0;
  double ScalersCollectedChargePlus = 0.0;
  double ScalersCollectedChargeMinus = 0.0;
  double ScalersChargeAsymmetry = 0.0;
  double ScalersChargeAsymmetryError = 0.0;


  for (int i = 0; i<num; i++)
  {
    t1->GetEntry(i);
    if ((i%1000==0) && i>0) cout<<i<<endl; 

//################ Primary Cuts ################################

    if (fabs(g0hel_L_helicity)<1) continue; // use only well defined helicity states

    //############### ENABLE THIS CUT ###############

    if (Ndata_DL_edtmbb!=0  ) continue; // we are interested only in real events;


                                                  
//##############################################################  

    time = clock1024count/1024; // time in [s]
    //double time = clock1024count/103700; // time in [s]
    double timediff = time - oldtime;
    //cout<<"Time: "<<time<<", Oldtime: "<<oldtime<<", TimeDiff: "<<timediff<<endl;
    oldtime = time;

    TotalCharge  =  (bcmu3count - time*off_u3)/calib_u3;
    double TotalChargeNP  =  (G1_bcmu3count/(G1_bcmu3count+0.000001))*(G1_bcmu3count - time*off_u3)/calib_u3;  // (T,h) = (N,P)
    double TotalChargeNN  =  (G2_bcmu3count/(G2_bcmu3count+0.000001))*(G2_bcmu3count - time*off_u3)/calib_u3;  // (T,h) = (N,N)

    double TotalChargePP  =  (G3_bcmu3count/(G3_bcmu3count+0.000001))*(G3_bcmu3count - time*off_u3)/calib_u3;  // (T,h) = (P,P)
    double TotalChargePN  =  (G4_bcmu3count/(G4_bcmu3count+0.000001))*(G4_bcmu3count - time*off_u3)/calib_u3;  // (T,h) = (P,N)


    ScalersCollectedChargePlus  = TotalChargeNP + TotalChargePP; // these two charges are XOR
    ScalersCollectedChargeMinus  = TotalChargeNN + TotalChargePN; // these two charges are XOR

    ScalersCollectedCharge = ScalersCollectedChargePlus + ScalersCollectedChargeMinus;
    
    double Current  =  (bcmu3rate - off_u3)/calib_u3;
    

    // First lets determine beam helicity;
    if (g0hel_L_helicity > 0 ){
      hBeamHelicityPlus->Fill(g0hel_L_helicity);
      CollectedChargePlus += Current*timediff;
    }
    else if (g0hel_L_helicity < 0 )
    { 
      hBeamHelicityMinus->Fill(g0hel_L_helicity);
      CollectedChargeMinus += Current*timediff;
    }


    CollectedCharge = CollectedChargePlus + CollectedChargeMinus;

    hCollectChargeDifference->Fill(TotalCharge - CollectedCharge);
    hCollectChargeDifferenceScalers->Fill(TotalCharge - ScalersCollectedCharge);

    grBCMu3cnt_vs_time->SetPoint(i, time, bcmu3count);
    grNA_BCMu3cnt_vs_time->SetPoint(i, N_clocktime, N_bcmu30);

    hBCMu3cnt_vs_time->Fill(time, bcmu3count);
    hNA_BCMu3cnt_vs_time->Fill(N_clocktime, N_bcmu30);

    hTimeDiff->Fill(timediff);
    hCurrentDistribution->Fill(Current -IBC0L02Current );

    grCharge_vs_time->SetPoint(i, time, TotalCharge);
    grCurrent_vs_time->SetPoint(i, time, Current);
    grEpicsCurrent_vs_time->SetPoint(i, time, IBC0L02Current);
    grEpicsCurrentHacMean_vs_time->SetPoint(i, time, hac_bcm_average);

    grCollectedCharge_vs_time->SetPoint(i, time, CollectedCharge);
    grCollectedChargePlus_vs_time->SetPoint(i, time, CollectedChargePlus);
    grCollectedChargeMinus_vs_time->SetPoint(i, time, CollectedChargeMinus);

    grCollectedChargeScalers_vs_time->SetPoint(i, time, ScalersCollectedCharge);
    grCollectedChargePlusScalers_vs_time->SetPoint(i, time, ScalersCollectedChargePlus);
    grCollectedChargeMinusScalers_vs_time->SetPoint(i, time, ScalersCollectedChargeMinus);
 
    
    if (CollectedChargePlus >0.0 || CollectedChargeMinus >0.0){
      ChargeAsymmetry = 
      (CollectedChargePlus - CollectedChargeMinus)/(CollectedChargePlus + CollectedChargeMinus);
      ChargeAsymmetryError = 0.0;

      ScalersChargeAsymmetry = 
      (ScalersCollectedChargePlus - ScalersCollectedChargeMinus)/(ScalersCollectedChargePlus + ScalersCollectedChargeMinus);
      ScalersChargeAsymmetryError = 0.0;

    }
    else { 
     ChargeAsymmetry = 0.0;
     ChargeAsymmetryError = 0.0;
     ScalersChargeAsymmetry = 0.0;
     ScalersChargeAsymmetryError = 0.0;
    }
    grCollectedChargeAsymmetry_vs_time->SetPoint(i, time, ChargeAsymmetry );
    grCollectedChargeAsymmetryScalers_vs_time->SetPoint(i, time, ScalersChargeAsymmetry );

  }
  





  //t1->Delete();
  delete t1; t1=NULL;


 
  //--- Canvas C1
  TCanvas *c1 = new TCanvas("c1","Beam Plots",800,800);
  c1->Divide(1,1);
  c1->cd(1);
  hBeamHelicityPlus->Draw();
  hBeamHelicityPlus->GetXaxis()->SetTitle("g0.L.helicity [1]");
  hBeamHelicityPlus->SetTitle("Beam Helicity");
  hBeamHelicityMinus->Draw("same");

  double BeamHelicityAsymmetry = ( hBeamHelicityPlus->GetEntries() - hBeamHelicityMinus->GetEntries() )/
                                 ( hBeamHelicityPlus->GetEntries() + hBeamHelicityMinus->GetEntries() ) ;
  double BeamHelicityAsymmetryError = 1.0 /
				sqrt( hBeamHelicityPlus->GetEntries() + hBeamHelicityMinus->GetEntries() + 0.00001 ) ;
  
  TString BeamAsymmetryString = Form(" Helicity Asymmetry = %f +/- %f",  BeamHelicityAsymmetry, BeamHelicityAsymmetryError);

  TLatex TextBeamAsymmetry;
  TextBeamAsymmetry.SetTextAlign(12);
  TextBeamAsymmetry.SetTextSize(0.03);
  TextBeamAsymmetry.SetTextColor(kRed);
  TextBeamAsymmetry.DrawLatex(-1.7,(int)hBeamHelicityPlus->GetMaximum()/2,BeamAsymmetryString);

  TString FigureFileName = Form("./Figures/figure_BeamPlots_%d.png",runnumber);
  c1->Print(FigureFileName.Data());


  //--- Canvas C2
  TCanvas *c2 = new TCanvas("c2","BCM",800,800);
  c2->Divide(2,2);
  c2->cd(1);
  gPad->SetGrid();
  grBCMu3cnt_vs_time->SetMarkerStyle(0);
  grBCMu3cnt_vs_time->SetMarkerColor(kBlue);
  grBCMu3cnt_vs_time->SetMarkerSize(1);
  grBCMu3cnt_vs_time->GetXaxis()->SetTitle("Time [s]");
  grBCMu3cnt_vs_time->GetYaxis()->SetTitle("BCM-u3");
  grBCMu3cnt_vs_time->Draw("AP");
  grNA_BCMu3cnt_vs_time->SetMarkerStyle(23);
  grNA_BCMu3cnt_vs_time->SetMarkerColor(kRed);
  grNA_BCMu3cnt_vs_time->SetMarkerSize(1);
  grNA_BCMu3cnt_vs_time->Draw("P");

  c2->cd(2);
  gPad->SetGrid();
  hBCMu3cnt_vs_time->GetXaxis()->SetTitle("Time [s]");
  hBCMu3cnt_vs_time->GetYaxis()->SetTitle("BCM-u3 (roc-11 scaler)");
  hBCMu3cnt_vs_time->Draw("colz");

  c2->cd(3);
  gPad->SetGrid();
  hNA_BCMu3cnt_vs_time->GetXaxis()->SetTitle("Time [s]");
  hNA_BCMu3cnt_vs_time->GetYaxis()->SetTitle("N.BCM-u3 (ev. 140)");
  hNA_BCMu3cnt_vs_time->Draw("colz");

  c2->cd(4);
  gPad->SetGrid();
  hTimeDiff->GetXaxis()->SetTitle("Time Difference [s]");
  hTimeDiff->Draw("");
  gPad->SetLogy();

  FigureFileName = Form("./Figures/figure_BCM_%d.png",runnumber);
  c2->Print(FigureFileName.Data());



  //--- Canvas C3
  TCanvas *c3 = new TCanvas("c3","Charge",800,800);
  c3->Divide(3,2);
  c3->cd(1);
  gPad->SetGrid();
  grCurrent_vs_time->SetMarkerStyle(0);
  grCurrent_vs_time->SetMarkerColor(kBlue);
  grCurrent_vs_time->SetMarkerSize(1);
  grCurrent_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCurrent_vs_time->GetYaxis()->SetTitle("I [uA]");
  grCurrent_vs_time->Draw("AP");

  grEpicsCurrent_vs_time->SetMarkerStyle(0);
  grEpicsCurrent_vs_time->SetMarkerColor(kRed);
  grEpicsCurrent_vs_time->SetMarkerSize(1);
  grEpicsCurrent_vs_time->Draw("P");

  grEpicsCurrentHacMean_vs_time->SetMarkerStyle(0);
  grEpicsCurrentHacMean_vs_time->SetMarkerColor(kGreen);
  grEpicsCurrentHacMean_vs_time->SetMarkerSize(1);
  grEpicsCurrentHacMean_vs_time->Draw("P");

  c3->cd(2);
  gPad->SetGrid();
  hCurrentDistribution->Draw();
  hCurrentDistribution->GetYaxis()->SetTitle("I_{BCM} - I_{IBC0L02} [A]");

  c3->cd(3);
  gPad->SetGrid();
  grCharge_vs_time->SetMarkerStyle(22);
  grCharge_vs_time->SetMarkerColor(kCyan);
  grCharge_vs_time->SetMarkerSize(1);
  grCharge_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCharge_vs_time->GetYaxis()->SetTitle("Q [uAs]");
  grCharge_vs_time->Draw("AP");
  grCollectedCharge_vs_time->SetMarkerStyle(0);
  grCollectedCharge_vs_time->SetMarkerColor(kRed);
  grCollectedCharge_vs_time->SetMarkerSize(1);
  grCollectedCharge_vs_time->Draw("P");

  grCollectedChargeScalers_vs_time->SetMarkerStyle(0);
  grCollectedChargeScalers_vs_time->SetMarkerColor(kBlue);
  grCollectedChargeScalers_vs_time->SetMarkerSize(1);
  grCollectedChargeScalers_vs_time->Draw("P");

  TString TotalChargeString = Form(" Charge (ungated BCM) = %f uAs", TotalCharge);
  TLatex TextTotalCharge;
  TextTotalCharge.SetTextAlign(12);
  TextTotalCharge.SetTextSize(0.04);
  TextTotalCharge.SetTextColor(kCyan);
  TextTotalCharge.DrawLatex( 10.0,TotalCharge*0.9, TotalChargeString);

  TString ChargeString = Form(" Calculated Charge Sum (+/-) = %f uAs",  CollectedCharge);
  TLatex TextCollectedCharge;
  TextCollectedCharge.SetTextAlign(12);
  TextCollectedCharge.SetTextSize(0.04);
  TextCollectedCharge.SetTextColor(kRed);
  TextCollectedCharge.DrawLatex( 10.0,CollectedCharge*0.7, ChargeString);

  TString ChargeScalersString = Form(" Charge Sum (Gated Scalers) = %f uAs", ScalersCollectedCharge );
  TLatex TextCollectedChargeScalers;
  TextCollectedChargeScalers.SetTextAlign(12);
  TextCollectedChargeScalers.SetTextSize(0.04);
  TextCollectedChargeScalers.SetTextColor(kBlue);
  TextCollectedChargeScalers.DrawLatex( 10.0,ScalersCollectedCharge*0.5, ChargeScalersString);


  c3->cd(4);
  gPad->SetGrid();
  grCollectedChargePlus_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCollectedChargePlus_vs_time->GetYaxis()->SetTitle("Q+ [uAs]");
  grCollectedChargePlus_vs_time->SetMarkerStyle(0);
  grCollectedChargePlus_vs_time->SetMarkerColor(kRed);
  grCollectedChargePlus_vs_time->SetMarkerSize(1);
  grCollectedChargePlus_vs_time->Draw("AP");

  grCollectedChargePlusScalers_vs_time->SetMarkerStyle(0);
  grCollectedChargePlusScalers_vs_time->SetMarkerColor(kBlue);
  grCollectedChargePlusScalers_vs_time->SetMarkerSize(1);
  grCollectedChargePlusScalers_vs_time->Draw("P");

  TString ChargePlusString = Form(" Collected (h+) Charge = %f uAs",  CollectedChargePlus);
  TLatex TextCollectedChargePlus;
  TextCollectedChargePlus.SetTextAlign(12);
  TextCollectedChargePlus.SetTextSize(0.04);
  TextCollectedChargePlus.SetTextColor(kRed);
  TextCollectedChargePlus.DrawLatex( 10.0,CollectedChargePlus*0.9, ChargePlusString);

  TString ChargePlusScalersString = Form(" Gated scalers (h+) Charge = %f uAs",  ScalersCollectedChargePlus);
  TLatex TextCollectedChargePlusScalers;
  TextCollectedChargePlusScalers.SetTextAlign(12);
  TextCollectedChargePlusScalers.SetTextSize(0.04);
  TextCollectedChargePlusScalers.SetTextColor(kBlue);
  TextCollectedChargePlusScalers.DrawLatex( 10.0,ScalersCollectedChargePlus*0.7, ChargePlusScalersString);

  c3->cd(5);
  gPad->SetGrid();
  grCollectedChargeMinus_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCollectedChargeMinus_vs_time->GetYaxis()->SetTitle("Q- [uAs]");
  grCollectedChargeMinus_vs_time->SetMarkerStyle(0);
  grCollectedChargeMinus_vs_time->SetMarkerColor(kRed);
  grCollectedChargeMinus_vs_time->SetMarkerSize(1);
  grCollectedChargeMinus_vs_time->Draw("AP");

  grCollectedChargeMinusScalers_vs_time->SetMarkerStyle(0);
  grCollectedChargeMinusScalers_vs_time->SetMarkerColor(kBlue);
  grCollectedChargeMinusScalers_vs_time->SetMarkerSize(1);
  grCollectedChargeMinusScalers_vs_time->Draw("P");

  TString ChargeMinusString = Form(" Collected (h-) Charge = %f uAs",  CollectedChargeMinus);
  TLatex TextCollectedChargeMinus;
  TextCollectedChargeMinus.SetTextAlign(12);
  TextCollectedChargeMinus.SetTextSize(0.04);
  TextCollectedChargeMinus.SetTextColor(kRed);
  TextCollectedChargeMinus.DrawLatex( 10.0,CollectedChargeMinus*0.9, ChargeMinusString);

  TString ChargeMinusScalersString = Form(" Gated scalers (h-) Charge = %f uAs",  ScalersCollectedChargeMinus);
  TLatex TextCollectedChargeMinusScalers;
  TextCollectedChargeMinusScalers.SetTextAlign(12);
  TextCollectedChargeMinusScalers.SetTextSize(0.04);
  TextCollectedChargeMinusScalers.SetTextColor(kBlue);
  TextCollectedChargeMinusScalers.DrawLatex( 10.0,ScalersCollectedChargeMinus*0.7, ChargeMinusScalersString);
/*
  c3->cd(6);
  gPad->SetGrid();
  hCollectChargeDifference->SetLineColor(kRed);
  hCollectChargeDifference->Draw();
  hCollectChargeDifferenceScalers->SetLineColor(kBlue);
  hCollectChargeDifferenceScalers->Draw("same");
*/


  FigureFileName = Form("./Figures/figure_Charge_%d.png",runnumber);
  c3->Print(FigureFileName.Data());


  //--- Canvas C3a
  TCanvas *c3a = new TCanvas("c3a","Charge",800,800);
  c3a->Divide(1,1);
  c3a->cd(1);
  gPad->SetGrid();
  grCurrent_vs_time->Draw("AP");
  grEpicsCurrent_vs_time->Draw("P");
  grEpicsCurrentHacMean_vs_time->Draw("P");

  FigureFileName = Form("./Figures/figure_Charge2_%d.png",runnumber);
  c3a->Print(FigureFileName.Data());

  //--- Canvas C4
  TCanvas *c4 = new TCanvas("c4","Charge Ratio",1200,800);
  c4->Divide(2,1);
  c4->cd(1);
  gPad->SetGrid();
  grCollectedChargeAsymmetry_vs_time->SetMarkerStyle(0);
  grCollectedChargeAsymmetry_vs_time->SetMarkerColor(kRed);
  grCollectedChargeAsymmetry_vs_time->SetMarkerSize(1);
  grCollectedChargeAsymmetry_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCollectedChargeAsymmetry_vs_time->GetYaxis()->SetTitle("Beam Charge Asymmetry [1]");
  grCollectedChargeAsymmetry_vs_time->Draw("AP");

  TString ChargeAsymmetryString = Form(" Final Charge Asymmetry = %f +/- ???",  ChargeAsymmetry, ChargeAsymmetryError);

  TLatex TextChargeAsymmetry;
  TextChargeAsymmetry.SetTextAlign(12);
  TextChargeAsymmetry.SetTextSize(0.03);
  TextChargeAsymmetry.SetTextColor(kRed);
  TextChargeAsymmetry.DrawLatex( 10.0,-0.1,ChargeAsymmetryString);

  c4->cd(2);
  gPad->SetGrid();
  grCollectedChargeAsymmetryScalers_vs_time->SetMarkerStyle(0);
  grCollectedChargeAsymmetryScalers_vs_time->SetMarkerColor(kRed);
  grCollectedChargeAsymmetryScalers_vs_time->SetMarkerSize(1);
  grCollectedChargeAsymmetryScalers_vs_time->GetXaxis()->SetTitle("Time [s]");
  grCollectedChargeAsymmetryScalers_vs_time->GetYaxis()->SetTitle("Beam Charge Asymmetry (Gated Scalers) [1]");
  grCollectedChargeAsymmetryScalers_vs_time->Draw("AP");

  TString ChargeAsymmetryScalersString = Form(" Final Charge Asymmetry (Gated Scalers) = %f +/- ???",  ScalersChargeAsymmetry, ScalersChargeAsymmetryError);

  TLatex TextChargeAsymmetryScalers;
  TextChargeAsymmetryScalers.SetTextAlign(12);
  TextChargeAsymmetryScalers.SetTextSize(0.03);
  TextChargeAsymmetryScalers.SetTextColor(kRed);
  TextChargeAsymmetryScalers.DrawLatex( 10.0,-0.1,ChargeAsymmetryScalersString);

  FigureFileName = Form("./Figures/figure_ChargeRatio_%d.png",runnumber);
  c4->Print(FigureFileName.Data());


  cout<<"######################################################################\n"<<endl;
  cout<<" RunNumber: "<<runnumber<<endl;
  cout<<" NumberOfEvents: "<<num<<endl;
  cout<<" Calibration Constants: calib_u3 = "<<calib_u3<<", off_u3 = "<<off_u3<<endl;
  cout<<" Time of run: "<<time<<" sec."<<endl;
  cout<<TotalChargeString.Data()<<endl;
  cout<<ChargeString.Data()<<endl;
  cout<<ChargePlusString.Data()<<endl;
  cout<<ChargeMinusString.Data()<<endl;
  cout<<BeamAsymmetryString.Data()<<endl;
  cout<<ChargeAsymmetryString.Data()<<endl;
  cout<<" ------------------------------------------------"<<endl;
  cout<<ChargeScalersString.Data()<<endl;
  cout<<ChargePlusScalersString.Data()<<endl;
  cout<<ChargeMinusScalersString.Data()<<endl;
  cout<<ChargeAsymmetryScalersString.Data()<<endl;
  cout<<"\n######################################################################"<<endl;

  //--- Now lets write this information to summaryfile
  ofstream outfile;
  TString OutFileName = Form("./summaryfiles/summaryfiles_BeamChargeAsymmetry_%d.dat",runnumber);
  outfile.open (OutFileName.Data());
  outfile<<"######################################################################\n"<<endl;
  outfile<<" RunNumber: "<<runnumber<<endl;
  outfile<<" NumberOfEvents: "<<num<<endl;
  outfile<<" Calibration Constants: calib_u3 = "<<calib_u3<<", off_u3 = "<<off_u3<<endl;
  outfile<<" Time of run: "<<time<<" sec."<<endl;
  outfile<<TotalChargeString.Data()<<endl;
  outfile<<ChargeString.Data()<<endl;
  outfile<<ChargePlusString.Data()<<endl;
  outfile<<ChargeMinusString.Data()<<endl;
  outfile<<BeamAsymmetryString.Data()<<endl;
  outfile<<ChargeAsymmetryString.Data()<<endl;
  outfile<<" ------------------------------------------------"<<endl;
  outfile<<ChargeScalersString.Data()<<endl;
  outfile<<ChargePlusScalersString.Data()<<endl;
  outfile<<ChargeMinusScalersString.Data()<<endl;
  outfile<<ChargeAsymmetryScalersString.Data()<<endl;
  outfile<<"\n######################################################################"<<endl;
  outfile.close();



}