//
// anal1.C:   Analyze FDC test stand data. The root data files contain data from
//            f1TDCs and FADC250
//            the code analyze1.C does the looping over the root file and fills the
//            histogram while this file initializes all the historgrams and calls 
//            the analysis function.
//
// hisotrams: fadc[i]          : raw FADC histograms of sum over the readout window
//            fadcNP[i]        : Sum minus the calculated pedestal from the first 6 
//                               bins of the readout window
//            histT[i]         : i=0,...,15 time difference in units of [ns]
//            histT[i]         : i=16,..,31 hit position along paddle in units of [cm]
//            hist2d           : hit position top vs. bottom scintillator paddles
//            qadc[32]         : raw ADC spectra from the V792 CAEN 32-channel ADC
//            DriftTimes[320]   : Drift times of the wires 
//         

// include files */
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
using namespace std;

// root include files
#include "TFile.h"
#include "TH1.h"
#include "TH1F.h" 
#include "TH2F.h"
#include "TH3F.h"
#include "TF1.h"
#include "TMath.h"


#include "anal1.h"

int main(int argc, char **argv) {


  //gROOT->Reset();
  //gROOT->ProcessLine(".L analyze1.C");
  

  TH1F *fadc[32];
  TH1F *fadcNP[32];
  Int_t adc_chan = 2001;
  Double_t adc_low = 20000.;
  Double_t adc_hig = 148000.;
  Double_t adc_l1 = -2000.;
  Double_t adc_h1 = 6000.;
  
  Int_t Flag1 = 0;
  Int_t Flag2 = 0;  

  TList * fadc_list= new TList();
  TList * fadcNP_list= new TList();
  char hnam[128];
  char htit[128];
  for (Int_t k=0;k<32;k++){
    
    sprintf(hnam,"fadc%d",k);
    sprintf(htit,"FADC channel %d",k);
    fadc[k] = new TH1F(hnam,"htit", adc_chan, adc_low, adc_hig);
    fadc_list->Add(fadc[k]);

    sprintf(hnam,"fadcNP%d",k);
    sprintf(htit,"FADC pedestal subtr. channel %d",k);
    fadcNP[k] = new TH1F(hnam, htit, adc_chan, adc_l1, adc_h1);
    fadcNP_list->Add(fadcNP[k]);
  }

  // create TDC histograms

  TH1F *histT[87];   
  TList * histT_list= new TList();
  for (Int_t k=0;k<16;k++){
    sprintf(hnam,"histT%d",k);
    if (k<8) {
      sprintf(htit,"TOP %d Tdiff [ns]",k);
    }else{
      sprintf(htit,"BOT %d Tdiff [ns]",k);
    }
    histT[k] = new TH1F(hnam,htit,100,-20., 20.);
    histT[k]->GetXaxis()->SetTitle("Time difference [ns]");
    histT_list->Add(histT[k]);
  }
  
  
  for (Int_t k=16;k<32;k++){
    sprintf(hnam,"histT%d",k);
    if (k<8) {
      sprintf(htit,"TOP %d dist [cm]",k);
    }else{
      sprintf(htit,"BOT %d dist [cm]",k);
    }
    histT[k] = new TH1F(hnam,htit,50,-150., 150.);
    histT[k]->GetXaxis()->SetTitle("Impact position [cm]");
    histT_list->Add(histT[k]);
 }
  
  histT[32] = new TH1F("histT32","TOF TOP to BOT",100,-15., 15.);
  histT[33] = new TH1F("histT33", "Cathode wire position X", 3600, -900., 900.);
  histT[34] = new TH1F("histT34", "Cathode wire position Y", 2000, -500., 500.);
  histT[35] = new TH1F("histT35", "Cathode wire fit position X", 200, -100., 100.);
  histT[36] = new TH1F("histT36", "Cathode wire fit position Y", 200, -100., 100.);
  histT[37] = new TH1F("histT37", "Cathode wire fit position Y", 2000, -500., 500.);
  histT[38] = new TH1F("histT38", "Cathode wire position X", 3600, -900., 900.);
  histT[39] = new TH1F("histT39", "Proj. of cath-y vs. y_IU to diagonal", 2000, 0., 500.);

  histT[40] = new TH1F("histT40", "Chi2 IU chamber fit x", 1000, 0., 100.);
  histT[41] = new TH1F("histT41", "Chi2 IU chamber fit y", 1000, 0., 100.);
  histT[42] = new TH1F("histT42", "IU chamber fit slope x", 100, -0.5, 0.5);
  histT[43] = new TH1F("histT43", "IU chamber fit slope y", 100, -0.5, 0.5);

  histT[44] = new TH1F("histT44", "Chi2 IU chamber fit x with slope cut", 1000, 0., 100.);
  histT[45] = new TH1F("histT45", "Chi2 IU chamber fit y with slope cut", 1000, 0., 100.);

  histT[46] = new TH1F("histT46", "Proj. of cath-y vs. y_IU to diagonal TOP1", 2000, 0., 500.);
  histT[47] = new TH1F("histT47", "Proj. of cath-y vs. y_IU to diagonal TOP2", 2000, 0., 500.);
  histT[48] = new TH1F("histT48", "Proj. of cath-y vs. y_IU to diagonal TOP3", 2000, 0., 500.);
  histT[49] = new TH1F("histT49", "Proj. of cath-y vs. y_IU to diagonal TOP4", 2000, 0., 500.);
  histT[50] = new TH1F("histT50", "Proj. of cath-y vs. y_IU to diagonal TOP5", 2000, 0., 500.);
  histT[51] = new TH1F("histT51", "Proj. of cath-y vs. y_IU to diagonal TOP6", 2000, 0., 500.);

  histT[52] = new TH1F("histT52", "cath-y - y_IU", 1600, -400., 0.);
  histT[53] = new TH1F("histT53", "FDC-x - x_IU", 600, -30., 30.);
  histT[54] = new TH1F("histT54", "FDC-x - x_IU chamber1", 600, -30., 30.);

  histT[55] = new TH1F("histT55", "Distance from FDC wire 1", 160, -20., 20.);
  histT[56] = new TH1F("histT56", "Distance from FDC wire 2", 160, -20., 20.);
  histT[57] = new TH1F("histT57", "Distance from FDC wire 3", 160, -20., 20.);
  histT[58] = new TH1F("histT58", "Distance from FDC wire 4", 160, -20., 20.);
  histT[59] = new TH1F("histT59", "Distance from FDC wire 5", 160, -20., 20.);
  histT[60] = new TH1F("histT60", "Distance from FDC wire 6", 160, -20., 20.);
  histT[61] = new TH1F("histT61", "Distance from FDC wire 7", 160, -20., 20.);
  histT[62] = new TH1F("histT62", "Distance from FDC wire 8", 160, -20., 20.);
  histT[63] = new TH1F("histT63", "Distance from FDC wire 9", 160, -20., 20.);
  histT[64] = new TH1F("histT64", "Distance from FDC wire 10", 160, -20., 20.);
  histT[65] = new TH1F("histT65", "Distance from FDC wire 11", 160, -20., 20.);
  histT[66] = new TH1F("histT66", "Distance from FDC wire 12", 160, -20., 20.);
  histT[67] = new TH1F("histT67", "Distance from FDC wire 13", 160, -20., 20.);
  histT[68] = new TH1F("histT68", "Distance from FDC wire 14", 160, -20., 20.);
  histT[69] = new TH1F("histT69", "Distance from FDC wire 15", 160, -20., 20.);
  histT[70] = new TH1F("histT70", "Distance from FDC wire 16", 100, -20., 20.);

  histT[71] = new TH1F("histT71", "FDC-IU wire 1", 160, -20., 20.);
  histT[72] = new TH1F("histT72", "FDC-IU wire 2", 160, -20., 20.);
  histT[73] = new TH1F("histT73", "FDC-IU wire 3", 160, -20., 20.);
  histT[74] = new TH1F("histT74", "FDC-IU wire 4", 160, -20., 20.);
  histT[75] = new TH1F("histT75", "FDC-IU wire 5", 160, -20., 20.);
  histT[76] = new TH1F("histT76", "FDC-IU wire 6", 160, -20., 20.);
  histT[77] = new TH1F("histT77", "FDC-IU wire 7", 160, -20., 20.);
  histT[78] = new TH1F("histT78", "FDC-IU wire 8", 160, -20., 20.);
  histT[79] = new TH1F("histT79", "FDC-IU wire 9", 160, -20., 20.);
  histT[80] = new TH1F("histT80", "FDC-IU wire 10", 160, -20., 20.);
  histT[81] = new TH1F("histT81", "FDC-IU wire 11", 160, -20., 20.);
  histT[82] = new TH1F("histT82", "FDC-IU wire 12", 160, -20., 20.);
  histT[83] = new TH1F("histT83", "FDC-IU wire 13", 160, -20., 20.);
  histT[84] = new TH1F("histT84", "FDC-IU wire 14", 160, -20., 20.);
  histT[85] = new TH1F("histT85", "FDC-IU wire 15", 160, -20., 20.);
  histT[86] = new TH1F("histT86", "FDC-IU wire 16", 160, -20., 20.);

  for ( Int_t k=32;k<87;k++){
    histT_list->Add(histT[k]);
  }

  histT[32]->GetXaxis()->SetTitle("Flight time [ns]");
  histT[33]->GetXaxis()->SetTitle("Cathode strip wire position [mm]");
  
  TH2F *hist2d[66]; 
  TList * hist2d_list= new TList();
  hist2d[0]= new TH2F("hist2d0","top vs bot",8, -0.5, 7.5, 8, -0.5,7.5);
  hist2d[1]= new TH2F("hist2d1","Cathode ypos vs xpos", 200, -100.5, 99.5, 200, -100.5, 99.5);
  hist2d[2]= new TH2F("hist2d2","Cathode wireypos vs wirexpos", 3600, -800., 800.,3600, -800.,800.);
  hist2d[3]= new TH2F("hist2d3","All Chambers", 800, 200., 1000., 351, -0.5, 350.5);
  hist2d[4]= new TH2F("hist2d4","FDCWire vs. Cathode pos x", 1000, -600., 400., 20, -1.5, 18.5);
  hist2d[5]= new TH2F("hist2d5","IU pos y vs. x", 1800, -100., 800., 1800, -100., 800.);
  hist2d[6]= new TH2F("hist2d6","FDCWire vs. Cathode pos y", 800, -400., 400., 20, -1.5, 18.5);
  hist2d[7]= new TH2F("hist2d7","Cathode wireypos vs wirexpos", 3600, -800., 800., 3600, -800., 800.);
  hist2d[8]= new TH2F("hist2d8","QADCs", 1000, 0., 4096., 64, -0.5, 63.5);
  hist2d[9]= new TH2F("hist2d9","QADCs calib", 1000, 0., 4096., 64, -0.5, 63.5);


  hist2d[10]= new TH2F("hist2d10","FDCwire vs. IU_x", 801, 0., 800.,20, -1.5, 18.5);
  hist2d[11]= new TH2F("hist2d11","Cathode x-pos vs. IU_x", 1801, -600., 600., 401, -200., 200.);
  hist2d[12]= new TH2F("hist2d12","Cathode y-pos vs. IU_y", 801, 0., 800., 801, -200., 200.);
  hist2d[13]= new TH2F("hist2d13","FDCwire vs. IU_y", 801, 0., 800.,20, -1.5, 18.5);

  hist2d[14]= new TH2F("hist2d14","slopey vs. IU_y_chath_y projection", 2001, 0., 500., 100, -0.5, 0.5);

  hist2d[15]= new TH2F("hist2d15","Trigger Paddle vs. histpos BOT", 300, -150., 150., 8, -1.5, 6.5);
  hist2d[16]= new TH2F("hist2d16","Trigger Paddle vs. histpos TOP", 300, -150., 150., 8, -1.5, 6.5);

  hist2d[17]= new TH2F("hist2d17","FDC_x vs. IU_x", 600, 0., 400., 600, 0., 400. );
  hist2d[18]= new TH2F("hist2d18","IU_slopex vs. IU_slopey", 101, -0.5, 0.5, 101, -0.5, 0.5 );

  hist2d[19]= new TH2F("hist2d19","IU pos y vs. x FDCw1 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[20]= new TH2F("hist2d20","IU pos y vs. x FDCw2 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[21]= new TH2F("hist2d21","IU pos y vs. x FDCw3 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[22]= new TH2F("hist2d22","IU pos y vs. x FDCw4 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[23]= new TH2F("hist2d23","IU pos y vs. x FDCw5 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[24]= new TH2F("hist2d24","IU pos y vs. x FDCw6 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[25]= new TH2F("hist2d25","IU pos y vs. x FDCw7 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[26]= new TH2F("hist2d26","IU pos y vs. x FDCw8 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[27]= new TH2F("hist2d27","IU pos y vs. x FDCw9 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[28]= new TH2F("hist2d28","IU pos y vs. x FDCw10 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[29]= new TH2F("hist2d29","IU pos y vs. x FDCw11 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[30]= new TH2F("hist2d30","IU pos y vs. x FDCw12 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[31]= new TH2F("hist2d31","IU pos y vs. x FDCw13 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[32]= new TH2F("hist2d32","IU pos y vs. x FDCw14 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[33]= new TH2F("hist2d33","IU pos y vs. x FDCw15 P1", 301, -100., 800., 301, -100., 800.);
  hist2d[34]= new TH2F("hist2d34","IU pos y vs. x FDCw16 P1", 301, -100., 800., 301, -100., 800.);

  hist2d[35]= new TH2F("hist2d35","FDC1 vs FDC2", 18, -0.5, 17.5, 18, -0.5, 17.5);
 
  hist2d[36]= new TH2F("hist2d36","FDC_IU_diff vs trackangle", 150, -0.5, 0.5, 60, -20., 20.);
  hist2d[37]= new TH2F("hist2d37","FDC_IU_diff vs trackazimuth", 150, -0.5, 0.5, 60, -20., 20.);


  hist2d[38]= new TH2F("hist2d38","IU pos y vs. x FDCw1 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[39]= new TH2F("hist2d39","IU pos y vs. x FDCw2 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[40]= new TH2F("hist2d40","IU pos y vs. x FDCw3 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[41]= new TH2F("hist2d41","IU pos y vs. x FDCw4 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[42]= new TH2F("hist2d42","IU pos y vs. x FDCw5 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[43]= new TH2F("hist2d43","IU pos y vs. x FDCw6 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[44]= new TH2F("hist2d44","IU pos y vs. x FDCw7 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[45]= new TH2F("hist2d45","IU pos y vs. x FDCw8 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[46]= new TH2F("hist2d46","IU pos y vs. x FDCw9 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[47]= new TH2F("hist2d47","IU pos y vs. x FDCw10 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[48]= new TH2F("hist2d48","IU pos y vs. x FDCw11 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[49]= new TH2F("hist2d49","IU pos y vs. x FDCw12 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[50]= new TH2F("hist2d50","IU pos y vs. x FDCw13 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[51]= new TH2F("hist2d51","IU pos y vs. x FDCw14 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[52]= new TH2F("hist2d52","IU pos y vs. x FDCw15 P2", 301, -100., 800., 301, -100., 800.);
  hist2d[53]= new TH2F("hist2d53","IU pos y vs. x FDCw16 P2", 301, -100., 800., 301, -100., 800.);

  hist2d[54]= new TH2F("hist2d54","QADC1 vs QADC2", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[55]= new TH2F("hist2d55","QADC2 vs QADC3", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[56]= new TH2F("hist2d56","QADC3 vs QADC4", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[57]= new TH2F("hist2d57","QADC4 vs QADC5", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[58]= new TH2F("hist2d58","QADC5 vs QADC6", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[59]= new TH2F("hist2d59","QADC6 vs QADC7", 1500, 0., 1500., 1500, 0., 1500.);

  hist2d[60]= new TH2F("hist2d60","QADC33 vs QADC34", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[61]= new TH2F("hist2d61","QADC34 vs QADC35", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[62]= new TH2F("hist2d62","QADC35 vs QADC36", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[63]= new TH2F("hist2d63","QADC36 vs QADC37", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[64]= new TH2F("hist2d64","QADC37 vs QADC38", 1500, 0., 1500., 1500, 0., 1500.);
  hist2d[65]= new TH2F("hist2d65","QADC38 vs QADC39", 1500, 0., 1500., 1500, 0., 1500.);


  for ( Int_t k=0;k<66;k++){
    hist2d_list->Add(hist2d[k]);
  }
  hist2d[2]->GetXaxis()->SetTitle("Cathode wire pos [mm]");
  hist2d[2]->GetYaxis()->SetTitle("Cathode wire length [arb.units]");
  hist2d[3]->GetXaxis()->SetTitle("All chambers drift time");
  hist2d[3]->GetYaxis()->SetTitle("All chambers TDC channel");
  hist2d[4]->GetYaxis()->SetTitle("FDC wire hit");
  hist2d[4]->GetXaxis()->SetTitle("Cathode x position [mm]");
  hist2d[5]->GetYaxis()->SetTitle("y position [mm]");
  hist2d[5]->GetXaxis()->SetTitle("x position [mm]");
  hist2d[18]->GetXaxis()->SetTitle("IU-chamber yslope");
  hist2d[18]->GetYaxis()->SetTitle("IU-chamber xslope");
 

  // histograms for CAEN ADC values

  TH1F *qadc[128];
  TList * qadc_list= new TList();
  for (Int_t k=0;k<128;k++){
    char hnam[128];
    char htit[128];
    
    sprintf(htit,"CAEN V792 ADC channel %d",k);
    sprintf(hnam,"qadc%d",k);  
    qadc[k] = new TH1F(hnam, htit, 4097, 0., 4096.);
    qadc_list->Add(qadc[k]);
  }
  
  // histograms for wire drift times
  TH1F *DriftTimes[350];
  TList * DriftTimes_list= new TList();
 for (Int_t k=0;k<350;k++){
    char hnam[128];
    char htit[128];
    
    sprintf(htit,"Drift Time wire number %d",k);
    sprintf(hnam,"DriftTimes%d",k);  
    DriftTimes[k] = new TH1F(hnam, htit, 1000, -2000., 2000.);
    DriftTimes_list->Add(DriftTimes[k]);

  }
   

  Int_t run=0;
  cout << "Run number : " ;
  std::cin >> run;
 
  char file[256];  
  char datadir[] = ".";
  char file1[] = "fdcTest_00000";
  char file2[] = "fdcTest_0000";
  char file3[] = "fdcTest_000";
  char file4[] = "fdcTest_00";
  char file5[] = "fdcTest_0";
  char *fp;
  if (run < 10) 
    fp = file1;
  else if (run < 100) 
    fp = file2;
  else if (run < 1000) 
    fp = file3;
  else if (run < 10000) 
    fp = file4;
  else 
    fp = file5;

  Int_t fext = 0;
  cout << "File extention: ";
  std::cin >> fext; 

  sprintf(file,"%s/%s%d.root.00%d",datadir,fp,run,fext);

  cout << endl 
       << "----------------------------------------------------------------------" << endl;
  cout << "analyzing : " << file << endl;
  cout << "----------------------------------------------------------------------" << endl;

 
  analyze1(run, file, fadc, fadcNP, histT, hist2d, qadc, DriftTimes, Flag1, Flag2);

  cout<<"End of run "<<run<<", save histograms ......"<<endl;
  flush(cout);

  char fnam[128];
  sprintf(fnam,"histograms_run%d_00%d.root",run,fext);
  TFile *fout = new TFile(fnam,"RECREATE");
  fadc_list->Write("fadc_list", TObject::kSingleKey);
  fadcNP_list->Write("fadcNP_list", TObject::kSingleKey);
  histT_list->Write("histT_list", TObject::kSingleKey);
  hist2d_list->Write("hist2d_list", TObject::kSingleKey);
  qadc_list->Write("qadc_list", TObject::kSingleKey);
  DriftTimes_list->Write("DriftTimes_list", TObject::kSingleKey);
  fout->Close();
}