#define e121_analysis_cxx

// INPUT FILES
#define INPUT1 "/data.local1/E121/root-files/part1/e121_run_0253.root"
//#define INPUT2 ""
//#define INPUT3 ""
//#define INPUT4 ""
//#define INPUT5 ""
//#define INPUT6 ""
//#define INPUT7 ""
//#define INPUT8 ""
//#define INPUT9 ""
//#define INPUT10 ""

//OUTPUT FILE
#define OUTPUT "run253.ana.root"

#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <cmath>
#include <string>
#include <cstdlib>
#include <cstdio>
#include <ctime>
#include <fstream>
#include <iostream>
#include <stdint.h>
#include "TROOT.h"
#include "TAttText.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TCut.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1.h"
#include "THistPainter.h"
#include "TKey.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TMultiGraph.h"
#include "TNtuple.h"
#include "TPave.h"
#include "TPaveText.h"
#include "TPoint.h"
#include "TRandom.h"
#include "TRint.h"
#include "TString.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TSystem.h"
#include "TProfile.h"
#include "TVirtualFitter.h"
#include "TLegend.h"
#include "TColor.h"
#include "TBranch.h"
#include "TList.h"

using namespace std;

void loop(TChain *fChain)
{
  if (fChain == 0) return;
  
  //#include "e121_raw.h" // use this for raw mapping of channels
#include "e121_mapped.h"  // use this for mapped channel names 

  ///////////////////////////
  // histogram definitions //
  ///////////////////////////
  // energy
  TList *e_histo= new TList();
    
  TH1D *h_e_pad_1[7];
  TH1D *h_e_pad_2[7];
  TH1D *h_e_pad_3[7];
  TH1D *h_e_pad_4[7];
  TH1D *h_e_pad_5[7];
  TH1D *h_e_pad_6[7];
  TH1D *h_e_pad_n[7];
  TH1D *h_e_pad_n_lo[7];
  TH1D *h_e_pad_n_hi[7];

  
  char chnumber[300];
  for (int i=0;i<7;i++)
    {
      sprintf(chnumber,"e_pad_1_strip_%d",i);
      h_e_pad_1[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_1[i]);

      sprintf(chnumber,"e_pad_2_strip_%d",i);
      h_e_pad_2[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_2[i]);
      
      sprintf(chnumber,"e_pad_3_strip_%d",i);
      h_e_pad_3[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_3[i]);
      
      sprintf(chnumber,"e_pad_4_strip_%d",i);
      h_e_pad_4[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_4[i]);
      
      sprintf(chnumber,"e_pad_5_strip_%d",i);
      h_e_pad_5[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_5[i]);
      
      sprintf(chnumber,"e_pad_6_strip_%d",i);
      h_e_pad_6[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_6[i]);
      
      sprintf(chnumber,"e_pad_%d_n-side",i);
      h_e_pad_n[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_n[i]);

      sprintf(chnumber,"e_pad_%d_n-side_lo",i);
      h_e_pad_n_lo[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_n_lo[i]);

      sprintf(chnumber,"e_pad_%d_n-side_hi",i);
      h_e_pad_n_hi[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
      e_histo->Add(h_e_pad_n_hi[i]);
    }

  TH1D *h_e_dssd_left=new TH1D("e_dssd_left","e_dssd_left",4096,0.5,4096.5);
  e_histo->Add(h_e_dssd_left);
  TH1D *h_e_dssd_right=new TH1D("e_dssd_right","e_dssd_right",4096,0.5,4096.5);
  e_histo->Add(h_e_dssd_right);
  TH1D *h_e_dssd_top=new TH1D("e_dssd_top","e_dssd_top",4096,0.5,4096.5);
  e_histo->Add(h_e_dssd_top);
  TH1D *h_e_dssd_bottom=new TH1D("e_dssd_bottom","e_dssd_bottom",4096,0.5,4096.5);
  e_histo->Add(h_e_dssd_bottom);

  TH1D *h_e_csi[2];
  h_e_csi[0]=new TH1D("e_csi_0","e_csi_0",4096,0.5,4096.5);
  e_histo->Add(h_e_csi[0]);
  h_e_csi[1]=new TH1D("e_csi_1","e_csi_1",4096,0.5,4096.5);
  e_histo->Add(h_e_csi[1]);
  
  // position
   
  
  // time
    
  
  // scaler
  TList *sc_histo= new TList();
    
  Int_t sc_bins = 500;
  TGraph *g_sc_si=new TGraph();
  g_sc_si->SetNameTitle("g_sc_si", "rate Si OR");
  sc_histo->Add(g_sc_si);
  TGraph *g_sc_mwpc=new TGraph();
  g_sc_mwpc->SetNameTitle("g_sc_mwpc", "rate MWPC anode");
  sc_histo->Add(g_sc_mwpc);
  TGraph *g_sc_target=new TGraph();
  g_sc_target->SetNameTitle("g_sc_target", "target density");
  sc_histo->Add(g_sc_target);
  TGraph *g_sc_trafo=new TGraph();
  g_sc_trafo->SetNameTitle("g_sc_trafo", "ESR ion current");
  sc_histo->Add(g_sc_trafo);
  TGraph *g_sc_inhibit=new TGraph();
  g_sc_inhibit->SetNameTitle("g_sc_inhibit", "target inhibit");
  sc_histo->Add(g_sc_inhibit);
  TGraph *g_sc_clock=new TGraph();
  g_sc_clock->SetNameTitle("g_sc_clock", "1.5 MHz clock");
  sc_histo->Add(g_sc_clock);
  
  
  ULong64_t sc_si_sum = 0;
  ULong64_t sc_mwpc_sum = 0;
  ULong64_t sc_clock_sum = 0;
  ULong64_t sc_target_sum = 0;
  ULong64_t sc_trafo_sum = 0;
  ULong64_t sc_inhibit_sum = 0;
  
  Int_t sc_interval = 1; //interval time [sec] for scaler plotting
  Double_t clock_ref = 1.5e6; // clock reference rate = 1.5 MHz
  Double_t t_diff = 0;
  Double_t t_elapsed = 0;
  
  // get number of events to process
  Long64_t nentries = fChain->GetEntries();
  Long64_t nbytes = 0, nb = 0;

  //////////////////////
  // EVENT LOOP START //
  //////////////////////
  for (Long64_t i=0; i<nentries;i++)
  //  for (Long64_t i=0; i<10000;i++)
    {
      nb = fChain->GetEntry(i);   nbytes += nb;

      // event countdown
      if ((float(i)/100000.)==int(i/100000)){cout << "event: " << i << " \tof " << nentries << endl;}

      //////////////////////
      //scaler processing //
      //////////////////////
      sc_si_sum += sc_silicon_or;
      sc_mwpc_sum += sc_mwpc_anode;
      sc_clock_sum += sc_clock;
      sc_target_sum += sc_target;
      sc_trafo_sum += sc_trafo;
      sc_inhibit_sum += sc_inhibit;
      t_elapsed += (double)(sc_clock/clock_ref);
      t_diff = (double)(sc_clock_sum/clock_ref); 

      // fill the scaler sums to histo for very interval 
      if( !((int)t_diff % sc_interval) && (int)t_diff > 0 )
      //if( true )
	{

	  if ( (int)t_diff <= 2*sc_interval ) //exclude large time intervals >> crap data
	    {

//	      cout << "time: " << t_elapsed << " trafo:" << sc_trafo_sum;
//	      cout << " diff: " << t_diff << " modulo: " << ((int)t_diff % sc_interval) << endl;
	      
	      g_sc_si->SetPoint(g_sc_si->GetN()          ,t_elapsed, sc_si_sum);
	      g_sc_mwpc->SetPoint(g_sc_mwpc->GetN()      ,t_elapsed, sc_mwpc_sum);
	      g_sc_target->SetPoint(g_sc_target->GetN()  ,t_elapsed, sc_target_sum);
	      g_sc_trafo->SetPoint(g_sc_trafo->GetN()    ,t_elapsed, sc_trafo_sum);
	      g_sc_inhibit->SetPoint(g_sc_inhibit->GetN(),t_elapsed, sc_inhibit_sum);
	      g_sc_clock->SetPoint(g_sc_clock->GetN()    ,t_elapsed, sc_clock_sum);
	    }
	  
	  sc_si_sum = 0;
	  sc_mwpc_sum = 0;
	  sc_clock_sum = 0;
	  sc_target_sum = 0;
	  sc_trafo_sum = 0;
	  sc_inhibit_sum = 0;
	}
	
      ////////////////////
      // energy spectra //
      ////////////////////
      if(TRIGGER==1) 
	{
	  for (int i=0;i<7;i++)
	    {
	      h_e_pad_1[i]->Fill(e_pad_1[i]);
	      h_e_pad_2[i]->Fill(e_pad_2[i]);
	      h_e_pad_3[i]->Fill(e_pad_3[i]);
	      h_e_pad_4[i]->Fill(e_pad_4[i]);
	      h_e_pad_5[i]->Fill(e_pad_5[i]);
	      h_e_pad_6[i]->Fill(e_pad_6[i]);
	      if ( i > 0 )
		{
		  h_e_pad_n[i]->Fill(e_pad_n[i]);
		  if ( t_elapsed >= 30  && t_elapsed <= 100 ) h_e_pad_n_hi[i]->Fill(e_pad_n[i]);
		  if ( t_elapsed >= 250 && t_elapsed <= 350 ) h_e_pad_n_lo[i]->Fill(e_pad_n[i]);
		}
	    }
	}


      
    }
  ////////////////////
  // EVENT LOOP END //
  ////////////////////

  // write output file
  
  TFile *outfile = TFile::Open(OUTPUT, "RECREATE");

  e_histo->Write("energy", TObject::kSingleKey);
  sc_histo->Write("scaler", TObject::kSingleKey);
  
  outfile->Close();
  
  cout << "\033[0;32m" << OUTPUT << " is created!\033[0m" << endl;
   
  // direct plotting (optional)

    TCanvas *c_test = new TCanvas("test","test", 1000, 1200);
    c_test->Divide(1,5); 
    c_test->cd(1);
    g_sc_si->SetFillColor(15);
    g_sc_si->Draw("AB");
    c_test->cd(2);
    g_sc_mwpc->SetFillColor(15);
    g_sc_mwpc->Draw("AB");
    c_test->cd(3);
    g_sc_target->SetFillColor(15);
    g_sc_target->Draw("AB");
    c_test->cd(4);
    g_sc_trafo->SetFillColor(15);
    g_sc_trafo->Draw("AB");
    c_test->cd(5);
    g_sc_inhibit->SetFillColor(15);
    g_sc_inhibit->Draw("AB");

    TCanvas *c_pad = new TCanvas("pad_1_6","strips of pads 1 to 6", 2000, 1400);
    c_pad->Divide(6,7);
    for (int i=0;i<7;i++)
      {
	c_pad->cd(i*6+1);
	h_e_pad_1[i]->Draw();
	c_pad->cd(i*6+2);
	h_e_pad_2[i]->Draw();
	c_pad->cd(i*6+3);
	h_e_pad_3[i]->Draw();
	c_pad->cd(i*6+4);
	h_e_pad_4[i]->Draw();
	c_pad->cd(i*6+5);
	h_e_pad_5[i]->Draw();
	c_pad->cd(i*6+6);
	h_e_pad_6[i]->Draw();
      }

    TCanvas *c_pad_n = new TCanvas("pad_n","n-side of pads 1 to 6", 2000, 1200);
    c_pad_n->Divide(2,3);
    for (int i=1;i<7;i++)
      {
	c_pad_n->cd(i);
	h_e_pad_n[i]->Draw();
	h_e_pad_n_lo[i]->SetLineColor(kRed);
	h_e_pad_n_lo[i]->Draw("SAME");
	h_e_pad_n_hi[i]->SetLineColor(kBlack);
	h_e_pad_n_hi[i]->Draw("SAME");
      }
}

void run(){

  const char    *command = new char[1000];
  char filename[100];
  TChain *fChain = new TChain("h101");

  sprintf(filename,INPUT1);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);

#ifdef INPUT2 // optional second input file
  sprintf(filename,INPUT2);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT3 // optional third input file
  sprintf(filename,INPUT3);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT4 // optional third input file
  sprintf(filename,INPUT4);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT5 // optional third input file
  sprintf(filename,INPUT5);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT6 // optional third input file
  sprintf(filename,INPUT6);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT7 // optional third input file
  sprintf(filename,INPUT7);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT8 // optional third input file
  sprintf(filename,INPUT8);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT9 // optional third input file
  sprintf(filename,INPUT9);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif
#ifdef INPUT10 // optional third input file
  sprintf(filename,INPUT10);
  cout<<"\033[0;37m//loading run: "<<filename << "\033[0m" <<endl;
  fChain->Add(filename);
#endif

  loop(fChain);

  command = "rm  *.so";
  gSystem->Exec(command);
  command = "rm  *.d";
  gSystem->Exec(command);
  command = "rm  *.pcm";
  gSystem->Exec(command);

  return;
}

int main(){
  run();

  return(0);
}