#include <iostream>
#include <fstream>
#include <string>
#include "TTree.h"
#include "TTreeReader.h"
#include "TBranch.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TBranch.h"
#include "TCut.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TString.h"
#include "TMultiGraph.h"

using namespace std;

int main(int argc, char* argv[]) {
  
  string runlistFileName  = argv[1]; //when I execute the macro with ./ I have to write the file's name which contains the list of the scans' file
  //string runlistFileName = argv[1];
  cout << "runlistFileName = " << runlistFileName << endl;

  ifstream runlistFile(runlistFileName.data());
  string out_file_name = "2023-02-28_scan.root";
  TFile *out_file = new TFile(out_file_name.c_str(), "recreate");

  if (!runlistFile) {
      cout << "No run list file found!" << endl;
      exit(1);
   }

  cout << "List of file taken from: " << argv[1] << endl;

  Int_t file_counter = 0;
  
  while(!runlistFile.eof())
    { //loop over input file
     
      TString fileName;
      runlistFile >> fileName;
      cout << "Input file: " << fileName << endl;
      ifstream current_infile(fileName);
      //current_infile = fileName;
      //current_infile.open(fileName);
   
      //line counters
      Int_t Npts = 0;
      string line;
      while(getline(current_infile,line))
	{
	  Npts++;
	}
      cout << "Entries: " << Npts<< endl;
      current_infile.close();
          current_infile.open(fileName);
      Double_t B[Npts];   //magnetic field 
      Double_t eB[Npts];    
      Double_t i_DF[Npts];  //current of our faraday cup
      Double_t ei_DF[Npts];
      Double_t i_MiB[Npts];   //current of Milan's faraday cup 
      Double_t ei_MiB[Npts];
      Double_t dummy;
      
      for(Int_t i=0; i< Npts; i++)
	{
	  current_infile >> B[i] >> eB[i] >> i_DF[i] >> ei_DF[i] >> dummy >> dummy >> i_MiB[i] >> ei_MiB[i];
	  cout << B[i] << " " << eB[i] << " " << i_DF[i] << " " << ei_DF[i] << endl;
	}

      TGraphErrors *g = new TGraphErrors(Npts, B, i_DF, eB, ei_DF);
      TString g_name = Form("g_%02i", file_counter);
      g->SetTitle(g_name);
      g->Draw("AP*");
      g->Write();
      file_counter++;

    //current_infile.clear();
    
    //current_infile.close(fileName);
  }
    
  
  
  //end of while loop


 
    
    
    
  //tempo / corrente vector used for global graph
  vector<Float_t> tempo;
  vector<Float_t> corrente;
  vector<Float_t> campo;
  Float_t offset = 0;
   
  while(!runlistFile.eof()) { //loop over input file

    //open current file
    TString fileName;
    runlistFile >> fileName;
    cout << "Input file: " << fileName << endl;
    TFile *myFile = TFile::Open(fileName);

    //Define TreeReader for current file
    TTreeReader myReader("run", myFile);
    TTreeReaderValue<Float_t> time(myReader, "time");
    TTreeReaderValue<Float_t> I_b(myReader, "I_b");
    TTreeReaderValue<Float_t> field(myReader, "field");
    TTreeReaderValue<Float_t> I_discharge_r(myReader, "I_discharge_r");
    TTreeReaderValue<Float_t> I_filament_r(myReader, "I_filament_r");
    
    //define offset for time realignmetn of each file
    //loop over current file events
    while(myReader.Next()) {

   //   cout << *time << " " << *I_b << endl;
      //if(*field < 6420 || *field > 6470 || *I_discharge_r < 10) continue;

      if(*field < 6430 || *field > 6460) continue;
       tempo.push_back(*time/3600 + offset);//divided by 3600 to move to hour
       if(*I_b>1.5e-6) *I_b =1e-9; //taglio su massima corrente elimina spike
      if(*I_b<1e-8) *I_b =1e-9;
      corrente.push_back(*I_b);
      
    }

    myReader.Delete();
    // time.Delete();
    //I_b.Delete();
    //new offset is the last time
    offset = tempo[tempo.size() - 1];
    
    cout << "Al file: " << fileName << " numero di entries nell' array: " << tempo.size() << endl;
    cout << "Al file: " << fileName << " valore di tempo all'ultima entries: " << offset << endl;
  }//end loop on file


  TCanvas *can = new TCanvas();
  can->Divide(1,3);can->cd(1);
  
  //plot graph 
  TGraph *g = new TGraph(tempo.size(), &tempo[0], &corrente[0]);
 // TGraph *g = new TGraph(tempo.size(), time, I_b);
 g->GetXaxis()->SetTitle("Time [h]");
 g->GetYaxis()->SetTitle("Current [A]");
 // g->GetYaxis()->SetRangeUser(-1e-8,1.2e-6);
 g->SetMarkerStyle(20);
 g->SetMarkerSize(0.1);
 g->Draw("AP");
 can->cd(2);
  
  
 //find max of each scan

 //get vector entries
  // cout << "numero di entries corrente: " << corrente.size() << endl;
  // cout << "numero di entries tempo: " << tempo.size() << endl;

 //vectors to store max current and corresponding time
 vector<Float_t> max_i;
 vector<Float_t> max_t;

 Float_t temp_i_max = -1e-5;
 Float_t temp_t_max = 0;
 max_i.push_back(1e-7);
 Int_t counter = 0; //contatore per eliminare massimi casuali dovuti al riallineamento del campo. 
 for(Int_t i=1; i<corrente.size(); i++) {
   
   // cout << "tempo " << i << " " << tempo[i] << endl;
   if(corrente[i]>temp_i_max) {//qui controllo se corrente aumenta

     temp_i_max = corrente[i];
     temp_t_max = tempo[i];
     // cout << "evento " << i << " tempo " << temp_t_max  << " corrente " << temp_i_max << endl;
     counter++;

   }
   
   if((tempo[i+1] - tempo[i])>0.5/3600) { //divided by 3600 to consider hour
     
     // cout << "salto al picco successivo e riempo vettore dei massimi " << endl;
     //  if(abs(temp_i_max - max_i[i]) < (0.95 * max_i[i-1]) ) continue;
     if(counter <2) continue; //se non ho almeno 3 conteggi attorno al massimo non lo salvo
     if(temp_i_max <100e-9) continue;
     max_i.push_back(temp_i_max);
     // cout << " tempo salvato " << temp_t_max << " corrente salvata " <<  temp_i_max << endl;
     max_t.push_back(temp_t_max);
     temp_i_max = -1e-5; //ripristino valore di controllo per massimo della corrente
     counter = 0;
   }

 }
 //closeloop over current value for maximum search
 
 
 //print max value for check
 for(Int_t i=0; i<max_i.size()-1; i++) {//check questo -1
 
  cout << "massimo " << i << " tempo " << max_t[i] << " altezza " << max_i[i] << endl;
 
 }

 
 TGraph *g_max = new TGraph(max_i.size(), &max_t[0], &max_i[0]);
 // TGraph *g = new TGraph(tempo.size(), time, I_b);
 g_max->GetXaxis()->SetTitle("Time [h]");
 g_max->GetYaxis()->SetTitle("Current [A]");
 // g_max->GetYaxis()->SetRangeUser(-1e-8,1.2e-6);
 g_max->SetMarkerStyle(20);
 g_max->SetMarkerSize(0.4);
 g_max->Draw("AP");
 can->cd(3);

 vector<Float_t> integral;
 integral.push_back(0.);
 
 for(Int_t i=0; i<max_i.size()-1; i++) { //qui aggiungo fattore 1.20 per t< 7h

   if(max_t[i] < 7) {
   integral.push_back(integral[i]+(max_t[i+1]-max_t[i])*max_i[i+1]*3600*6.2e18*1.2);//multiplied by 3600 to return to s, multiplied by 6.2e18 to convert to ions
   cout << "integrale " << integral[i] << endl;
   }
   if(max_t[i]>7) {

     integral.push_back(integral[i]+(max_t[i+1]-max_t[i])*max_i[i+1]*3600*6.2e18);//multiplied by 3600 to return to s, multiplied by 6.2e18 to convert to ions
   cout << "integrale " << integral[i] << endl;


   }
   
 }

 TGraph *g_integral = new TGraph(integral.size(), &max_t[0], &integral[0]);
 // TGraph *g = new TGraph(tempo.size(), time, I_b);
 g_integral->GetXaxis()->SetTitle("Time [h]");
 g_integral->GetYaxis()->SetTitle("Integrated Ho atoms");
 // g_max->GetYaxis()->SetRangeUser(-1e-8,1.2e-6);
 g_integral->SetMarkerStyle(20);
 g_integral->SetMarkerSize(0.4);
 g_integral->Draw("AP");
 //can->cd(3);

 TCanvas *can1 = new TCanvas();
 g->Draw("AP");
  TCanvas *can2 = new TCanvas();
 g_max->Draw("AP");
  TCanvas *can3 = new TCanvas();
 g_integral->Draw("APF");
 
     
    
//return 0; 

}//close macro