// COMET software includes
#include "IMidasBank.hxx"
#include "IMidasFile.hxx"
#include "IMidasBankProxy.hxx"
#include "IMidasBankProxyRegistry.hxx"
#include "ICOMETRawEvent.hxx"

// oaRawEvent includes
#include "ITripTDigitBank.hxx"
#include "IMCMBank.hxx"
#include "IMidasTripTDigitItr.hxx"
#include "IMidasTripTDigit.hxx"

// ROOT includes
#include "TApplication.h"
#include <TSpectrum.h>
#include "TFile.h"
#include "TH1F.h"
#include "TTree.h"
#include "TH2F.h"
#include "TF1.h"
#include <TStyle.h>

#include "TString.h"
#include "TSystem.h"

#include "TClonesArray.h"

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <cstdio>
#include <iostream>
#include <fstream>

TFile* myFile    = NULL;
TH2F*  allPeds  = NULL;
TTree* TMPPCFit = NULL;
TH1F*  TimeDistribution = NULL;
TH1F*  TimeDistribution2 = NULL;

ULong64_t fTDC[300];
Short_t fC64[300];
UInt_t fTDCO[300];
Short_t fTDCOdiff[300];
ULong64_t *fMCMT0;

ULong64_t tmpTDC[4][32][64][23];
Short_t tmpC64[4][32][64][23];
UInt_t tmpTDCO[4][32][64][23];
UInt_t tmpMCMT;

Short_t tmpIDCO[10000];
Short_t tmpIDRMM[10000];
Short_t tmpIDTFB[10000];
Short_t tmpIDC64[10000];

ULong64_t *fTriggerW;
Int_t fTr;

Int_t rtc2SMRD[4][32][64];

Int_t fSMRD[300];

Float_t f1stPeakX[8192];
Float_t f1stPeakXe[8192];
Float_t f1stPeakH[8192];
Float_t f1stPeakHe[8192];
Float_t f1stPeakS[8192];
Float_t f1stPeakSe[8192];

Float_t f2ndPeakX[8192];
Float_t f2ndPeakXe[8192];
Float_t f2ndPeakH[8192];
Float_t f2ndPeakHe[8192];
Float_t f2ndPeakS[8192];
Float_t f2ndPeakSe[8192];

Float_t f3rdPeakX[8192];
Float_t f3rdPeakXe[8192];
Float_t f3rdPeakH[8192];
Float_t f3rdPeakHe[8192];
Float_t f3rdPeakS[8192];
Float_t f3rdPeakSe[8192];

Float_t fGain[8192];
Float_t fGaine[8192];
ULong64_t f1stMCMSec=0;
ULong64_t flstMCMSec=0;

ULong64_t gMCMSec;
ULong64_t gtrigSubSec;
ULong64_t gtrigTriggerW;

//____________________________________________________________________________

void Book() {
  // Book some very basic histograms to show the data
  allPeds = new TH2F("allPeds","allPeds",18000,0.,18000.,1024,0.,1024.);
  TimeDistribution = new TH1F("TimeDistribution","MPPC Time Distribution in SMRD [ns]",1400,0,13999);
  TimeDistribution2 = new TH1F("TimeDistribution2","MPPC Time Distribution in SMRD [ns]",120,0,600);

  TMPPCFit =new TTree("TMPPCFit","MPPC Fit Information");
  TMPPCFit->Branch("f1stPeakH",f1stPeakH,"f1stPeakH[8192]/F");
  TMPPCFit->Branch("f1stPeakHe",f1stPeakHe,"f1stPeakHe[8192]/F");
  TMPPCFit->Branch("f1stPeakX",f1stPeakX,"f1stPeakX[8192]/F");
  TMPPCFit->Branch("f1stPeakXe",f1stPeakXe,"f1stPeakXe[8192]/F");
  TMPPCFit->Branch("f1stPeakS",f1stPeakS,"f1stPeakS[8192]/F");
  TMPPCFit->Branch("f1stPeakSe",f1stPeakSe,"f1stPeakSe[8192]/F");

  TMPPCFit->Branch("f2ndPeakH",f2ndPeakH,"f2ndPeakH[8192]/F");
  TMPPCFit->Branch("f2ndPeakHe",f2ndPeakHe,"f2ndPeakHe[8192]/F");
  TMPPCFit->Branch("f2ndPeakX",f2ndPeakX,"f2ndPeakX[8192]/F");
  TMPPCFit->Branch("f2ndPeakXe",f2ndPeakXe,"f2ndPeakXe[8192]/F");
  TMPPCFit->Branch("f2ndPeakS",f2ndPeakS,"f2ndPeakS[8192]/F");
  TMPPCFit->Branch("f2ndPeakSe",f2ndPeakSe,"f2ndPeakSe[8192]/F");

  TMPPCFit->Branch("f3rdPeakH",f3rdPeakH,"f3rdPeakH[8192]/F");
  TMPPCFit->Branch("f3rdPeakHe",f3rdPeakHe,"f3rdPeakHe[8192]/F");
  TMPPCFit->Branch("f3rdPeakX",f3rdPeakX,"f3rdPeakX[8192]/F");
  TMPPCFit->Branch("f3rdPeakXe",f3rdPeakXe,"f3rdPeakXe[8192]/F");
  TMPPCFit->Branch("f3rdPeakS",f3rdPeakS,"f3rdPeakS[8192]/F");
  TMPPCFit->Branch("f3rdPeakSe",f3rdPeakSe,"f3rdPeakSe[8192]/F");

  TMPPCFit->Branch("fGain",fGain,"fGain[8192]/F");
  TMPPCFit->Branch("fGaine",fGaine,"fGaine[8192]/F");

  TMPPCFit->Branch("f1stMCMSec",&f1stMCMSec,"f1stMCMSec/l");
  TMPPCFit->Branch("flstMCMSec",&flstMCMSec,"flstMCMSec/l");
}
//____________________________________________________________________________


Int_t ExecuteFit(TH2F* allPedsFit){

  Float_t *Xaxis;
  Float_t *Yaxis;
  char    func[200];
  int NPeaks;
  int AbbyN;
  int yanoN;
  TSpectrum *ts = new TSpectrum();
  //TSpectrum ts(10,1);

  if (((Int_t)allPedsFit->GetEntries())/4016<200*23){
    std::cerr<<"This file do not have enough entry for Fitting"<<std::endl;
    return 1;
  }

  for(int nRMM=0;nRMM<=3;nRMM++){
    for(int nTFB=0;nTFB<=31;nTFB++){
      for(int nCH=0;nCH<=63;nCH++){
        AbbyN=(100*(48*nRMM + nTFB) + nCH);
        yanoN=(2048*nRMM + 64*nTFB + nCH);
        TH1D *h1 = allPedsFit->ProjectionY("project",AbbyN+1,AbbyN+1);
        h1->SetAxisRange(100,200);
        //ts.Search(h1,1,"goff",0.013);
        ts->Search(h1,1,"goff",0.013);
        NPeaks=ts->GetNPeaks();
        if (NPeaks>3)NPeaks=3;
        Xaxis=ts->GetPositionX();
        Yaxis=ts->GetPositionY();
        sprintf(func,"gaus(0)");
        for (int i=1 ; i<NPeaks ; i++) {
          sprintf(func,"%s+gaus(%d)",func,3*i);
        }
        TF1* uTgaus = new TF1("uTgaus",func,100,200);
        for (int i=0 ; i<NPeaks ; i++) {
          uTgaus->SetParameter(3*i,Yaxis[i]);
          uTgaus->SetParameter(3*i+1,Xaxis[i]);
          uTgaus->SetParameter(3*i+2,1.5);
        }
        h1->Fit("uTgaus","0Q","",100,200);

        f1stPeakH[yanoN] = uTgaus->GetParameter(0);
        f1stPeakHe[yanoN] = uTgaus->GetParError(0);
        f1stPeakX[yanoN] = uTgaus->GetParameter(1);
        f1stPeakXe[yanoN] = uTgaus->GetParError(1);
        f1stPeakS[yanoN] = uTgaus->GetParameter(2);
        f1stPeakSe[yanoN] = uTgaus->GetParError(2);

        f2ndPeakH[yanoN] = uTgaus->GetParameter(3);
        f2ndPeakHe[yanoN] = uTgaus->GetParError(3);
        f2ndPeakX[yanoN] = uTgaus->GetParameter(4);
        f2ndPeakXe[yanoN] = uTgaus->GetParError(4);
        f2ndPeakS[yanoN] = uTgaus->GetParameter(5);
        f2ndPeakSe[yanoN] = uTgaus->GetParError(5);

        f3rdPeakH[yanoN] = uTgaus->GetParameter(6);
        f3rdPeakHe[yanoN] = uTgaus->GetParError(6);
        f3rdPeakX[yanoN] = uTgaus->GetParameter(7);
        f3rdPeakXe[yanoN] = uTgaus->GetParError(7);
        f3rdPeakS[yanoN] = uTgaus->GetParameter(8);
        f3rdPeakSe[yanoN] = uTgaus->GetParError(8);

        fGain[yanoN] = f2ndPeakX[yanoN]-f1stPeakX[yanoN];
        fGaine[yanoN] = sqrt((f2ndPeakXe[yanoN])*(f2ndPeakXe[yanoN])+(f1stPeakXe[yanoN])*(f1stPeakXe[yanoN]));

        h1->Delete();
        if(yanoN%100==0) std::cout <<"CH " <<yanoN <<" Processed!"<<std::endl << std::flush;
      };
    };
  };
  ts->Delete();
  return 0;
};

void ProcessTDC(Int_t i) {
  fTriggerW = 0;
  Int_t   Coinc_Bit = 0;
  Short_t Coinc_Diff;
  Int_t co,rmm,tfb;
  Int_t chan64,chan64p;
  Int_t j=0;

  //Extracting Timing
  fMCMT0    = (ULong64_t *)gtrigSubSec;
  fTr       = (Int_t )(((Long64_t )gtrigTriggerW>>48)&0xFFFF);

  // start roop for all hits. They are tagged with "co(integration window)" and channels.
  // They are numbered with "i".
  for (int k=0 ; k <= i ; k++){
    co  = tmpIDCO[k];
    rmm = tmpIDRMM[k];
    tfb = tmpIDTFB[k];
    chan64 = tmpIDC64[k];
    if (j==300){
      std::cout << "More than 300 Coinc Hits are found in an event. They will not filled to the Tree." <<std::endl;
      j++;
    }
    if (j>299){
      j++;
    }
    else if (j<300){
      //hits pairing
      if (chan64<16)
        Coinc_Bit  =  48;
      else if (chan64<32&&chan64>15)
        Coinc_Bit  =  16;
      else if (chan64<48&&chan64>31)
        Coinc_Bit  = -16;
      else if (chan64<64&&chan64>47)
        Coinc_Bit  = -48;
      chan64p  =  chan64+Coinc_Bit;

      Coinc_Diff=-255; // Time Difference between 2 pair hits.
      if(tmpTDC[rmm][tfb][chan64][co]>0 && tmpTDC[rmm][tfb][chan64p][co]>0){
        Coinc_Diff=tmpTDC[rmm][tfb][chan64][co] - tmpTDC[rmm][tfb][chan64p][co];
      }
      if (abs(Coinc_Diff)<10){
        fC64[j]  = chan64; 
        fTDC[j]  = tmpTDC[rmm][tfb][chan64][co];
        fTDCO[j]  = tmpTDCO[rmm][tfb][chan64][co];
        fTDCOdiff[j]=Coinc_Diff;
        j++;
      }
    }
  }
  for (int m=0 ; m<= j ; m++){
    if((fTr&0x1)&&fTDCO[m]>0&&fC64[m]<32) {
      Double_t abstime = ((Long64_t ) fTDC[m] - (Long64_t) gtrigSubSec)%400000000-fTDCOdiff[m]/2;
      TimeDistribution->Fill(abstime*2.5);
      TimeDistribution2->Fill(((Int_t)((abstime*2.5)-320))%581);
    }
  }
  return;
}




void Event(COMET::ICOMETRawEvent* re) {
  Int_t i=0;
  //Int_t lowadc; 
  Int_t highadc;
  Int_t chan64;
  Int_t tfb ;
  UInt_t timeo ;
  //ULong64_t timet0;
  //UInt_t timeti ;
  Int_t cycle  ;
  //Int_t rmm ;
  Int_t chanID;

  ULong64_t time ;
  //Int_t trip ;

  Int_t start_cycle ;
  Int_t co;
  //Int_t chan64p;
  // Loop over all banks of type TTripTHitBank
  COMET::IHandle<COMET::ITripTDigitBank> triptBank;
  COMET::IHandle<COMET::IMCMBank> mcmBank;
  const Int_t mcmTick(4);

  //MCM Trigger Time Inplementation Start
  while (mcmBank = re->GetMidasBank<COMET::IMCMBank>("OMCM",mcmBank)){
    gMCMSec = mcmBank->GetUnixTimeTrig();
    gtrigSubSec = mcmBank->GetUnixTimeSSecTrig()*mcmTick;
    gtrigTriggerW = (ULong64_t)mcmBank->GetTriggerWord();
    if (f1stMCMSec==0)f1stMCMSec = gMCMSec;
    flstMCMSec = gMCMSec;
  }
  //MCM Trigger Time Inplementation End

  while ( triptBank = re->GetMidasBank<COMET::ITripTDigitBank>("",triptBank) ) {
    const Char_t *name = triptBank->GetName();
    Int_t rmm = (name[2]>'9')?(name[2]-'A'+10):(name[2]-'0');
    // RMM is third digit in the bank name, 0123456789ABCDEFG... 
    /*
       Int_t detector = (name[0]=='P') ? 0 :      // P0D
       (name[0]=='I') ? 1 :      // INGRID
       (name[0]=='E') ? 2 :      // ECAL
       (name[0]=='S') ? 3 : 4;   // SMRD : 4=Unknown
     */

    // printf("Bank name %s rmm %d detector %d\n",name,rmm,detector);

    if (name[0]=='S')
    {
      // Create an iterator over digits
      COMET::IMidasTripTDigitItr itr(triptBank->GetMidasTripTDigitItr());

      while ( ! itr.EOD() ) {
        COMET::IMidasTripTDigit digit(itr.Get());

        //lowadc  =  digit.GetLowGainADC();
        highadc =  digit.GetHighGainADC();
        chan64  =  digit.GetChannelNum() + 16*digit.GetTripTNum();
        tfb     =  digit.GetTFBNum();
        timeo   =  digit.GetTimeOffset();
        //timeti  =  digit.GetTimeOffsetTI();
        cycle   =  digit.GetIntegrationNum();
        rmm     =  digit.GetRMMNum();
        //timet0  =  (ULong64_t) digit.GetTimeOffsetT0();
        time    =  (ULong64_t) digit.GetTime();
        //trip    =  digit.GetTripTNum();
        start_cycle   =  digit.GetTFBStartCycle();
        co      =  cycle - start_cycle;
        chanID  =  static_cast<Int_t>(100.*(48.*rmm + tfb) + chan64);
        if(co<0) co += 23;

        if ( (((Long64_t )gtrigTriggerW>>49)&0x1) ==0x1 )//Pedestal Trigger
        {
          allPeds->Fill(chanID,highadc);
        }

        if(highadc>100&&timeo!=16777201){
          tmpC64[rmm][tfb][chan64][co]   = chan64;
          tmpTDCO[rmm][tfb][chan64][co]  = timeo;
          tmpTDC[rmm][tfb][chan64][co]   = time;
          if (i>9999){
            std::cout << "More than 10000 Hits are found in an event. They will not filled to the Tree." <<std::endl;
            i++;
          }
          else if (i<=9999){
            tmpIDCO[i]  = co;
            tmpIDRMM[i] = rmm;
            tmpIDTFB[i] = tfb;
            tmpIDC64[i] = chan64;
          }
          i++;
        }
      }   // End of loop over digits in this bank
    }     // End of loop over banks of digits in this event
  }
  ProcessTDC(i);
  return;
}
//____________________________________________________________________________

void ProcessFile(const char *FileName) {
  // Start by listing the available Access Classes.
  COMETLog("Access Classes available to interpret the MIDAS banks");
  COMET::IMidasBankProxyRegistry::Instance().Print();

  // Make sure file exists
  FileStat_t fs;
  if ( gSystem->GetPathInfo(FileName,fs) ) {
    std::cerr << "Cannot find file: " << FileName << std::endl;
    return;
  }

  COMET::IMidasFile mf;
  mf.Open(FileName);

  // Loop over events in this file
  while ( COMET::ICOMETRawEvent* re = mf.ReadRawEvent() ) 
  {
    //
    //re->PromoteMidasBanks(false);
    Event(re);
    delete re;

  }  // End loop over events
  ExecuteFit(allPeds);
  TMPPCFit->Fill();
}

//____________________________________________________________________________
/*
   void Ly2s(){
   Int_t RMM,TFB,CH;
   Int_t Nsmrd, Nyano;
   char Trash[200];
   Float_t POS_X,POS_Y,POS_Z;

   Int_t ch;
   Float_t peakh1,peakx1,sigma1;
   Float_t peakh2,peakx2,sigma2;
   Float_t peakh3,peakx3,sigma3;
   Float_t peakh4,peakx4,sigma4;

   Float_t mean,rms;





   std::ifstream datalist("yano_SMRD_pos.lst");
   if (!datalist){
   std::cout << "NO FILE!"<<std::endl;
   };
   datalist>>Trash;
//cout<<"Start y2s Exec"<<endl;
while( datalist >> Nyano >> Nsmrd
>> RMM>>TFB>>CH
>>POS_X
>>POS_Y
>>POS_Z
){
rtc2SMRD[RMM][TFB][CH]=Nsmrd;
}
//cout<<"End y2s Exec"<<endl;
datalist.close();


std::ifstream datalist2("FitAllCH.dat");
if (!datalist2){
std::cout << "NO FILE! for Gain"<<std::endl;
};
datalist2>>Trash;
while( datalist2 >> ch
>> peakh1 >> peakx1 >> sigma1
>> peakh2 >> peakx2 >> sigma2
>> peakh3 >> peakx3 >> sigma3
>> peakh4 >> peakx4 >> sigma4
){
RMM=ch/2048;
TFB=(ch%2048)/64;
CH=ch%64;
Mped[RMM][TFB][CH]=peakx1;
Mgain[RMM][TFB][CH]=peakx2-peakx1;
}
datalist2.close();

std::ifstream datalist3("TDC_Calib.dat");
if (!datalist3){
std::cout << "NO FILE! for TDC Calib"<<std::endl;
};
datalist3>>Trash;
while( datalist3 >> RMM >> TFB >> mean >>rms){
MtdcC[RMM][TFB]=(Int_t )mean+1964;
}
datalist3.close();
}
 */
//____________________________________________________________________________


int main(int argc, char **argv) {
  // Read from File
  //Ly2s();

  // Get the filename from argv[1] before TApplication modifies the
  // argument list!!!
  if (argc != 3) {
    std::cerr << "usage:  " << argv[0] << " <infilename> <outfilename>" << std::endl;
  }
  char InFileName[1000];
  strncpy(InFileName,argv[1],1000);

  char OutFileName[1000];
  strncpy(OutFileName,argv[2],1000);

  printf("Infile: %s, Outfile %s\n",InFileName, OutFileName);

  // Start the root system
  TApplication *app = new TApplication("TCRTestRead", &argc, argv);
  gROOT->SetBatch();
  if (app) ;    // This just removes 'unused variable app' compiler warning

  // Open an output root file.
  myFile = new TFile(OutFileName,"RECREATE","Demonstration ROOT file");
  myFile->SetCompressionLevel(2);

  // Book the histograms
  Book();

  // Analyse the data
  ProcessFile(InFileName);

  // Close root histogram output file
  myFile->Write();
  myFile->Close();
  delete myFile;
  myFile = NULL;

  return 0;
}