#include <iostream>
#include <stdio.h>
#include <vector>

#include "TROOT.h"
#include "TStyle.h"
#include "TString.h"
#include "TFile.h"
#include "TH1I.h"
#include "TH1F.h"
#include "TNtuple.h"
#include "TBranch.h"
#include "TH1F.h"
#include "TF1.h"
#include "TGraphErrors.h"


#define BIGSEPARATOR std::cout << "********************************************************************************\n\n";

// method to get unique index (idx) and relative name (hGainName)
void hGainNameCreator( Int_t iLayer, Int_t iChip, Int_t iChannel, Int_t& idx, TString &hGainName );

// ntuple name 
const TString ntupleName = "ntuple";
// active layers histo name
const TString activeLayersHistoName = "activelayers";
/*
 * ranges for some indices
 */
const Int_t nLayers = 16;
const Int_t nChips = 3;
const Int_t nChannels = 128;
const Int_t nThresholds = 8;
const Int_t nADC = 8;
/*
 * arrays to feed the gain graphs
 */
Float_t meanArray[nLayers][nChips][nChannels][nThresholds]={0.};
Float_t errMeanArray[nLayers][nChips][nChannels][nThresholds]={0.};
Float_t VthDACArray[nLayers][nChips][nChannels][nThresholds]={0.};
Float_t DisVtnDACArray[nLayers][nChips][nChannels]={0.};
/* NOTA BENE : error for threshold dac set to 1. */
const Float_t errVthDACArray[nLayers][nChips][nChannels][nThresholds]={1.}; 

/* NOTA BENE: physical limits for Vin values */
const Float_t minVinDAC = 0.;
const Float_t maxVinDAC = 255.;

// vector to list active layers
std::vector<Int_t> activeLayers;

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

  // erf data file is needed as the one and the only input parameter
  if ( argc != 2 ) {
    std::cerr << "Usage: " << argv[0] << " <ErfDataFile>\n";
    exit ( 1 );
  }

  // cleaning stuff and setting proper style;
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(1111111);
  
  // this is the file name you passed
  TString *ErfDataFileName = new TString(argv[1]);
  std::cout << "\nYour ErfDataFile is " << ErfDataFileName->Data() << "\n\n";

  // trying to open input file
  TFile *ErfDataFile = new TFile(ErfDataFileName->Data(),"READ");
  if ( ErfDataFile->IsZombie() ) {
    std::cerr << "Unable to open ErfDataFile " << ErfDataFileName->Data() << "\n";
    exit( 2 );
  }

  ErfDataFile->cd();

  // looking for activeLayers histo
  TH1I* activeLayersHisto = NULL;
  activeLayersHisto = (TH1I*)ErfDataFile->Get(activeLayersHistoName.Data());
  if ( activeLayersHisto == NULL ) {
    std::cerr << "Unable to find Histo " << activeLayersHistoName.Data() << "\n";
    exit( 3 );
  }

  // looking for ntuple
  TNtuple *ntuple = NULL;
  ntuple = (TNtuple*)ErfDataFile->Get(ntupleName.Data());
  if ( ntuple == NULL ) {
    std::cerr << "Unable to find ntuple " << ntupleName.Data() << "\n";
    exit( 4 );
  }

  // everything is ok... let's go!
  BIGSEPARATOR
  std::cout << "Ntuple and activeLayers histo have been found.\n";
  std::cout << "Analysis will start now...\n";

  // Extracting the root of the file-name
  // it will be used as base for all the ps and
  // root file which will be created
  TString globalName =  ((TObjString*)ErfDataFileName->Tokenize(".")->At(0))->GetString();


  // booking TFile for output
  TString outputFileName("");
  outputFileName.Form("%s-Gain.root",globalName.Data());
  std::cout << "Output file name will be: " << outputFileName.Data() << "\n\n";
  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
  if ( outputFile ->IsZombie() ) {
    std::cerr << "Unable to open outputFile " << outputFileName.Data() << "\n";
    exit( 2 );
  }

  /*
   * A loop over all active channels will pick up 
   * the mean information of erf fit
   * We will plot it as a function of input the threshold
   */

   /*
    * variables to be read
    */

   Float_t fLayer(-1.), fChip(-1.), fChannel(-1.), fADC(-1.),
         Vth0DAC(-1.), Vth1DAC(-1.), Vth2DAC(-1.), Vth3DAC(-1.),
         Vth4DAC(-1.), Vth5DAC(-1.), Vth6DAC(-1.), Vth7DAC(-1.),
         DisVtnDAC(-1.);

   Float_t  chi2Prob(-1.),
             mean(-1.), errMean(-1.),
             sigma(-1.), errSigma(-1.);

   // some useful indices
   Int_t iLayer(-1), iChip(-1), iChannel(-1), iADC(-1);

   // information on active layers
   for ( Int_t iLayer = 0; iLayer < 16; iLayer++ ) {
     if ( activeLayersHisto->GetBinContent( iLayer + 1 ) > 0 ) {
        activeLayers.push_back( iLayer );
     }
   }
   std::cout << "A total of " << activeLayers.size() << " active layers have been found.\n";
   std::cout << "They are: \n\t";
   for ( iLayer = 0; iLayer < activeLayers.size(); iLayer++ ) {
     std::cout << activeLayers[iLayer] << " ";
   }
   std::cout << "\n\n";
 
   /*
    * branches to be filled
    */

   ntuple->SetBranchAddress("fLayer",&fLayer);

   ntuple->SetBranchAddress("fChip",&fChip);

   ntuple->SetBranchAddress("fChannel",&fChannel);

   ntuple->SetBranchAddress("fADC",&fADC);

   ntuple->SetBranchAddress("Vth0DAC",&Vth0DAC);

   ntuple->SetBranchAddress("Vth1DAC",&Vth1DAC);

   ntuple->SetBranchAddress("Vth2DAC",&Vth2DAC);

   ntuple->SetBranchAddress("Vth3DAC",&Vth3DAC);

   ntuple->SetBranchAddress("Vth4DAC",&Vth4DAC);

   ntuple->SetBranchAddress("Vth5DAC",&Vth5DAC);

   ntuple->SetBranchAddress("Vth6DAC",&Vth6DAC);

   ntuple->SetBranchAddress("Vth7DAC",&Vth7DAC);

   ntuple->SetBranchAddress("DisVtnDAC",&DisVtnDAC);

   ntuple->SetBranchAddress("chi2Prob",&chi2Prob);

   ntuple->SetBranchAddress("mean",&mean);

   ntuple->SetBranchAddress("errMean",&errMean);

   ntuple->SetBranchAddress("sigma",&sigma);

   ntuple->SetBranchAddress("errSigma",&errSigma);

   // this is the graph for the gain plot
   TString hGainName;
   Int_t idx;
   TH1F **hGainPos = (TH1F**) (new TH1F[nLayers*nChips*nChannels]);
   TH1F **hGainNeg = (TH1F**) (new TH1F[nLayers*nChips*nChannels]);

   /*
    * Loop over the ntuple!
    */
         
   Int_t nentries = Int_t(ntuple->GetEntries());
   for (Int_t i=0;i<8;i++) {
      ntuple->GetEntry(i);
      iLayer   = (Int_t)(fLayer);
      iChip    = (Int_t)(fChip);
      iChannel = (Int_t)(fChannel);
      iADC     = (Int_t)(fADC);
      if ( activeLayersHisto->GetBinContent( iLayer + 1 ) > 0 ) {
        #ifdef  _DEBUG_ 
        std::cout << "\nEvent # " << i;
        std::cout << " fLayer # " << fLayer;
        std::cout << " fChip # " << fChip;
        std::cout << " fChannel # " << fChannel << "\n";
        #endif
        Float_t VthMap[nThresholds] = { Vth0DAC, Vth1DAC, Vth2DAC, Vth3DAC,
                              Vth4DAC, Vth5DAC, Vth6DAC, Vth7DAC };
        meanArray[iLayer][iChip-1][iChannel][iADC] = mean;
        errMeanArray[iLayer][iChip-1][iChannel][iADC] = errMean;
        VthDACArray[iLayer][iChip-1][iChannel][iADC] = VthMap[iADC]; 
        DisVtnDACArray[iLayer][iChip-1][iChannel] = DisVtnDAC; 
        #ifdef  _DEBUG_ 
        std::cout << "\tRead from ntuple: \n";
        std::cout << "\t   fADC # " << fADC;
        std::cout << " mean # " << mean;
        std::cout << " errMean # " << errMean;
        std::cout << " VthDAC # " <<  VthMap[iADC] << "\n";
	std::cout << "\tWrote in TGraph: \n";
        std::cout << "\t   fADC # " << fADC;
        std::cout << " mean # " << meanArray[iLayer][iChip-1][iChannel][iADC];
        std::cout << " errMean # " << errMeanArray[iLayer][iChip-1][iChannel][iADC] ;
        std::cout << " VthDAC # " <<  VthDACArray[iLayer][iChip-1][iChannel][iADC] << "\n";
        #endif

      }
   }

   // moving to output file
   outputFile->cd();

   // booking and filling graph
   for ( iLayer = 0; iLayer < 16; iLayer++ ) {
     if ( activeLayersHisto->GetBinContent( iLayer + 1 ) > 0 ) {
       for ( iChip = 1; iChip <= 3; iChip++ ) {
         for ( iChannel = 0; iChannel < 128; iChannel++ ) {
           hGainNameCreator( iLayer, iChip, iChannel, idx, hGainName );
           TString hGainPosName("");
           TString hGainNegName(""); 
           hGainPosName += hGainName.Data();
           hGainPosName += "_pos";
           hGainNegName += hGainName.Data();
           hGainNegName += "_neg";
           hGainPos[idx] = new TH1F( hGainPosName.Data(), "(Gain)^-1",256,-0.5,254.5); 
           hGainPos[idx]->GetXaxis()->SetTitle("Vth (Threshold DAC)");
           hGainPos[idx]->GetYaxis()->SetTitle("Vin (Pulser    DAC)");
           hGainNeg[idx] = new TH1F( hGainNegName.Data(), "(Gain)^-1",256,-0.5,254.5); 
           hGainNeg[idx]->GetXaxis()->SetTitle("Vth (Threshold DAC)");
           hGainNeg[idx]->GetYaxis()->SetTitle("Vin (Pulser    DAC)");
           /* 
            * NOTA BENE: histo will be filled only if the fitted mean
            * lies in the proper range
            */
           Float_t currDisVtnDAC = DisVtnDACArray[iLayer][iChip-1][iChannel];
           Int_t currDisVtnDACIndex = Int_t(currDisVtnDAC);
           for ( iADC = 0; iADC < nADC; iADC++ ) {
             Float_t currMean = meanArray[iLayer][iChip-1][iChannel][iADC];
             Float_t currErrMean = errMeanArray[iLayer][iChip-1][iChannel][iADC];
             Float_t currVthDAC = VthDACArray[iLayer][iChip-1][iChannel][iADC];
             Int_t currVthDACIndex = Int_t(currVthDAC);
             if ( currMean >= minVinDAC && currMean <= maxVinDAC ) {
               if ( DisVtnDAC == 0 ) {
                 hGainPos[idx]->SetBinContent( currVthDACIndex+1, currMean );
                 hGainPos[idx]->SetBinError( currVthDACIndex+1, currErrMean );
               } else {
                 hGainNeg[idx]->SetBinContent( currDisVtnDACIndex-currVthDACIndex+1, currMean );
                 hGainNeg[idx]->SetBinError( currDisVtnDACIndex-currVthDACIndex+1, currErrMean );
               }
             }
           }   
         }
       }
     }
   }


   // saving histos
   Int_t nPosGains = 0;
   Int_t nNegGains = 0;
   for ( iLayer = 0; iLayer < 16; iLayer++ ) {
     if ( activeLayersHisto->GetBinContent( iLayer + 1 ) > 0 ) {
       for ( iChip = 1; iChip <= 3; iChip++ ) {
         for ( iChannel = 0; iChannel < 128; iChannel++ ) {
           if ( iChip == 1 && iChannel == 0 ) {
             Float_t currDisVtnDAC = DisVtnDACArray[iLayer][iChip-1][iChannel];
             hGainNameCreator( iLayer, iChip, iChannel, idx, hGainName );
             std::cout << "idx: " << idx << "\n";
             gStyle->SetOptFit(11111);
             gStyle->SetOptStat(11111);
             TF1 *f1 = new TF1("f1","pol1",0.,500.);
             if ( currDisVtnDAC == 0. ) {
               hGainPos[idx]->Write();
               hGainPos[idx]->Fit("f1");
               nPosGains++;
             } else {
               hGainNeg[idx]->Write();
               hGainNeg[idx]->Fit("f1");
               nNegGains++;
             }
           }
         }
       }
     }
   }

   activeLayersHisto->Write();
   outputFile->Close();
   ErfDataFile->cd();
   ErfDataFile->Close();

   Int_t nGains = nPosGains+nNegGains;
   
  std::cout << "\nA total of " << nPosGains << " positive-polarity channels have been analyzed.\n\n";
  std::cout << "\nA total of " << nNegGains << " negative-polarity channels have been analyzed.\n\n";
  std::cout << "\nA total of " << nGains << " channels have been analyzed.\n\n";
  BIGSEPARATOR
  std::cout << "Output file name is: " << outputFileName.Data() << "\n\n";
  return ( 0 );
}

void hGainNameCreator( Int_t iLayer, Int_t iChip, Int_t iChannel, Int_t &idx, TString &hGainName ) {
  
      idx = iLayer*3*128+(iChip-1)*128+iChannel;
      hGainName.Form("hGain_strip%03d_chip%02d_lay%02d", iChannel,iChip,iLayer);
}