//Created by Laszo Varga & Jan Glorius
//last update: 20.11.2020

//run in root:
//root -l
//.L script.c+
//run()
////////////////////////

// INPUT FILES
#define INPUT1 "/data.local3/test_data_2020/Si2_run002.lmd.root"

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

//Energy thresholds for the silicon:
#define Si_Thrs_LOW 100
#define Si_Thrs_HIGH 8000

////////////////////////

#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 "TH2.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 "TStyle.h"
#include "TString.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TSystem.h"
#include "TProfile.h"
#include "TVirtualFitter.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TColor.h"

using namespace std;

void loop(TChain *fChain){

#include "/data.local3/test_data_2020/e127.h"

//creating histos

//counting of variables starts from 0 --> eg.: first Si strip is 0.

TH2D *h_pos_Si_xy_raw=new TH2D("h_pos_si_xy_raw", "h_pos_si_xy_raw",16,-0.5,15.5,16,-0.5,15.5);
TH2D *h_pos_Si_xy_noDC=new TH2D("h_pos_si_xy_noDC", "h_pos_si_xy_noDC",16,-0.5,15.5,16,-0.5,15.5);
 
TH2D *h_E_Si_all=new TH2D("h_E_Si_all", "E vs. strip Nr.; strip nr.; energy ",32,-0.5,31.5,8192,0.5,8191.5);
 
TH1D *h_E_Si_x[16];
TH1D *h_E_Si_y[16];
TH2D *h_E_ev_x[16];
TH2D *h_E_ev_y[16];
 
for(int a=0;a<16;a++){
	h_E_Si_x[a]=new TH1D(Form("h_E_Si_x%d",a),Form("h_E_Si_x%d",a),8192,0.5,8191.5);
	h_E_Si_y[a]=new TH1D(Form("h_E_Si_y%d",a),Form("h_E_Si_y%d",a),8192,0.5,8191.5);
	h_E_ev_x[a]=new TH2D(Form("h_E_ev_x%d",a),Form("h_E_ev_x%d;ev;E",a),1024,0,55e6,4096,0.5,8191.5);
	h_E_ev_y[a]=new TH2D(Form("h_E_ev_y%d",a),Form("h_E_ev_y%d;ev;E",a),1024,0,55e6,4096,0.5,8191.5);
}
 
////////////////////////////////////////////////////////////////////////

Long64_t nentries = fChain->GetEntries();
Long64_t nbytes = 0;

int dc_y_max=-999;
int dc_x_max=-999;
int dc_Ey_max=0;
int dc_Ex_max=0;
int dc_ty_max=0;
int dc_tx_max=0;

for (Long64_t i=0; i<nentries;i++){
//for (Long64_t i=0; i<100000;i++){
  nbytes += fChain->GetEntry(i);
	// event countdown
	if ((float(i)/100000.)==int(i/100000)){cout << "event: " << i << " \tof " << nentries << endl;}
	
	if(trigger==1){

		for (int i_x=0;i_x<16;i_x++){
			for (int i_y=0;i_y<16;i_y++){
					
			  if(int(E_Si_x[i_x])> 2000 && int(E_Si_y[i_y])> 2000)
			    {
			      h_pos_Si_xy_raw->Fill(i_x,i_y);
			    }

			  if(int(E_Si_x[i_x])>dc_Ex_max && int(E_Si_y[i_y])>dc_Ey_max){
			    dc_y_max = i_y;
			    dc_x_max = i_x;
			    dc_Ey_max = E_Si_y[i_y];
			    dc_Ex_max = E_Si_x[i_x];
			    dc_ty_max = t_Si_y[i_y];
			    dc_tx_max = t_Si_x[i_x];
			  }					
			}

			//x-loop is missused for y also
			if(t_Si_x[i_x] > 0){
			  h_E_Si_all->Fill(i_x,E_Si_x[i_x]);
			  h_E_Si_x[i_x]->Fill(E_Si_x[i_x]);
			}
			if(t_Si_y[i_x] > 0){
			  h_E_Si_all->Fill(i_x + 16,E_Si_y[i_x]);
			  h_E_Si_y[i_x]->Fill(E_Si_y[i_x]);
			}
			
			
		}
	
		if(dc_x_max!=-999 && dc_y_max!=-999){
			h_pos_Si_xy_noDC->Fill(dc_x_max,dc_y_max);
			h_E_ev_x[dc_x_max]->Fill(i,dc_Ex_max);
			h_E_ev_y[dc_y_max]->Fill(i,dc_Ey_max);
		}
		
		dc_y_max=-999;
		dc_x_max=-999;
		dc_Ex_max=0;
		dc_Ey_max=0;
		dc_tx_max=0;
		dc_ty_max=0;
	}	
 }//event loop

 TFile *graphfile = TFile::Open(OUTPUT, "RECREATE");

//graphfile->mkdir("Si_map");
//graphfile->cd("Si_map");
 h_pos_Si_xy_raw->Write();
 h_pos_Si_xy_noDC->Write();
 h_E_Si_all->Write();
 
//graphfile->mkdir("Si_E_histos");
//graphfile->cd("Si_E_histos");
 for(int a=0;a<16;a++){
   h_E_Si_x[a]->Write();
   h_E_Si_y[a]->Write();
   h_E_ev_x[a]->Write();
   h_E_ev_y[a]->Write();
 }

graphfile->Close();
cout << "\033[0;32m" << OUTPUT << " is created!\033[0m" << endl;

} // loop

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

//sprintf(filename,OUTPUT);
//	TFile *outfile = new TFile(filename,"RECREATE");

  loop(fChain);

//  outfile->Write();

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

int main(){
 run();

  return(0);
}