#include <cassert>
#include <map>
#include <sstream>
#include <string>
#include <vector>
#include <fstream>
#include <utility>
#include <cmath>

#include <TH1I.h>
#include <TH2I.h>
#include <TSpectrum.h>
#include <TF1.h>
#include <TH2F.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TTree.h>
#include <TProfile.h>
#include <TLine.h>
#include <TBox.h>
#include <TFile.h>

#include <ICOMETLog.hxx>
#include <IOADatabase.hxx>
#include <cometEventLoop.hxx>
#include <ITFBChannelId.hxx>
#include <ITFBDigit.hxx>
#include <ITripTDigitBank.hxx>
#include <ITripTXTimeBank.hxx>
#include <IMidasTripTDigitItr.hxx>
#include <IMidasTripTDigit.hxx>
#include <IRawDataHeader.hxx>
#include <ITripTInfoBank.hxx>
#include <IGeometryDatabase.hxx>
#include <IMCMBank.hxx>
#include <IXmlBank.hxx>
#include <IXmlOdb.hxx>
#include <IOARuntimeParameters.hxx>

using namespace COMET;

class IP0dLi: public COMET::ICOMETEventLoopFunction {
public:
  IP0dLi() : fDB(COMET::IGeometryDatabase::Get()),fSaveChannelHistos(0),
       fOutputFileName("dq-p0d-li.root"){};

  virtual ~IP0dLi() {};

  void Initialize();
  void Finalize(ICOMETOutput*);
  using COMET::ICOMETEventLoopFunction::Finalize;
  bool SetOption(std::string option,std::string value);

  bool operator () (COMET::ICOMETEvent& event);

private:

  int DetLIBox(int p0dule,int layer);

  COMET::IGeometryDatabase &fDB;

  //Trigger Histogram
  //Used to track all trigger types in the data stream and their respective trigger rates
  TH1I* fTriggerType;
  TH1F* fTriggerFreq;

  //Raw ADC Map Histograms
  //Store ADC spectrum, extract the Mean and RMS of these distributions
  TH2I* fBarHighADCMap[4];
  TH2I* fBarLowADCMap[4];

  //Timing Histograms
  //Store the TDC spectrum, extract the Mean,RMS and bin with maximum TDC
  TH2F* fBarTDCMap[4];

  //Run information variables
  //Store run number, start and end time, flash number in a given amplitude, signal average
  int fRunId;
  double fendTime;
  double fstartTime;
  double fCurrentTime[4];
  int fLIFlashNum[4];
  int fLIEventNum;
  float fPrevAvg[4];

  //LI Settings
  //Store current ODB parameters. Currently hardcoded, but should be read directly from first event
  int fLIAmps[4][10];
  int fLISetNum[4];
  int fLIStartAmp;
  int fLIState;

  //Output file
  TFile* fOutputFile;

  //Output tree and variables
  TTree* fOutputTree;
  float fHighADCMean, fHighADCRMS, fLowADCMean, fLowADCRMS, fTDCMean, fTDCRMS, fTDCMaxBin;
  int fRMM,fTFB,fTripT,fChannel,fP0dule,fBar,fLayer,fPulser,fAmplitude,fFlashes;
  UInt_t fChannelId;
  double fTime;

  //Additional tree branches - ADC and TDC histograms
  TH1D *fHighADC;
  TH1D *fLowADC;
  TH1D *fTDC;

  //ADC Checker for Amplitude switching
  bool fLIADCSwitch[4];
  bool fStart;

  //Group Information
  //Will be used to produce a plot for visual checking to make sure the algorithm works
  std::vector< std::vector<int> > fOnEvents;
  std::vector< std::vector<int> > fOffEvents;
  std::vector<double> fEventTime;
  std::vector<std::vector<double> > fGroupTime;
  std::vector<std::vector<float> > fFracDiff;
  std::vector<std::vector<float> > fAvgSig;

  //prevFlash info
  int fPrevFlash[4];

  /// Cut on the fractional difference of ADC used in amplitude switching procedure.
  float fFracDiffCut;

  /// Flag whether to save Low gain ADC, High gain ADC and TDC histograms for each channel.
  /// This can be controled via runtime parameter with -O option.
  bool fSaveChannelHistos;

  /// Name of the output root file.
  std::string fOutputFileName;
};

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

void IP0dLi::Initialize() {

  std::cout << "\n===============================================" << std::endl;
  std::cout << "      Starting production of the LI tree       " << std::endl;
  std::cout << "===============================================" << std::endl;

  //Bit to trigger start of the run code (ODB reads)
  fStart=true;

  //Total event counter
  fLIEventNum=0;
  fLIStartAmp=0;

  //Initialize all box specific variables
  for(int i=0;i<4;i++){
    fOnEvents.push_back(std::vector<int>(0));
    fOffEvents.push_back(std::vector<int>(0));
    fGroupTime.push_back(std::vector<double>(0));
    fFracDiff.push_back(std::vector<float>(0));
    fAvgSig.push_back(std::vector<float>(0));

    fLIFlashNum[i]=0;
    fLISetNum[i]=0;

    fPrevAvg[i]=0;
    fCurrentTime[i]=0;

    fLIADCSwitch[i]=true;
    fPrevFlash[i]=0;
  }

  //Initialize trigger type histogram
  fTriggerType = new TH1I("Triggers","Triggers",1000,0,1000);
  fTriggerType->GetXaxis()->SetTitle("Trigger type");
  fTriggerType->GetXaxis()->SetTitle("N");

  fTriggerFreq = new TH1F("Trigger Frequency","Trigger Frequency",1000,0,1000);
  fTriggerType->GetXaxis()->SetTitle("Trigger type");
  fTriggerType->GetXaxis()->SetTitle("Frequency (Hz)");

  //Hardcoded ODB amplitude values, need to add code to read this from the ODB instead (add in number of flashes and width as well)
  for(int box=0;box<4;box++){
    fLIAmps[box][0]=1000;
    fLIAmps[box][1]=500;
    fLIAmps[box][2]=900;
    fLIAmps[box][3]=400;
    fLIAmps[box][4]=800;
    fLIAmps[box][5]=300;
    fLIAmps[box][6]=700;
    fLIAmps[box][7]=200;
    fLIAmps[box][8]=600;
    fLIAmps[box][9]=100;
  }

  std::cout << "----> Amplitude settings summary (#setting <--> amplitude) : " << std::endl;
  for(int box = 0; box < 4; box++){
    std::cout << "--> Box: " << box+1 << std::endl;
    for(int set = 0; set < 10; set++){
      std::cout << set << " <--> " << fLIAmps[box][set] << std::endl;
    }
  }

  //Get the value of the cut on fractional difference from the parameters file
  //fFracDiffCut = COMET::IOARuntimeParameters::Get().GetParameterD("p0dLI.FracDiffCut");
  fFracDiffCut = 0.04; // 4% relative change in the amplitude

  fHighADC = new TH1D("fHighADC","fHighADC;HADC",1024,0,1024);
  fLowADC = new TH1D("fLowADC","fLowADC;HADC",1024,0,1024);
  fTDC = new TH1D("fTDC","fTDC;TDC",800,0,800);

  std::cout << "-----------------------------------------------" << std::endl;
}

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

void IP0dLi::Finalize(ICOMETOutput* output)
{
  std::cout << "===============================================" << std::endl;
  std::cout << "                 Finalizing...                 " << std::endl;
  std::cout << "===============================================" << std::endl;

  //Check to make On and off are the same and delete leftover histograms
  for(int i=0;i<4;i++){
    if(fOnEvents[i].size()!=fOffEvents[i].size()){
      COMETLog("Box: " << i+1 << ", On=" << fOnEvents[i].size() << ", Off=" << fOffEvents[i].size());
    }
    if(fLIFlashNum[i]>=1){
      delete fBarHighADCMap[i];
      delete fBarLowADCMap[i];
      delete fBarTDCMap[i];
    }
  }

  //gSystem->ChangeDirectory(fOutputFile);
  fOutputFile->cd();

  //Create Frequency histogram and write Trigger histograms
  if ( fEventTime.size() > 0 ) {
    int mysize = fEventTime.size()-1;
    double runtime = fEventTime.at(mysize)-fEventTime.at(0);
    
    for(int i=0;i<1000;i++){
      double bincontent = fTriggerType->GetBinContent(i+1);
      fTriggerFreq->Fill(i,bincontent/runtime);
    }
  }

  //fTriggerType->Write();

  fTriggerFreq->GetXaxis()->SetTitle("Trigger type");
  fTriggerFreq->GetYaxis()->SetTitle("Trigger frequency (Hz)");
  fTriggerFreq->Write();

  //Code to produce the AvgGain plot and frac diff plot
  const int length[4] = {
    static_cast<int>(fOnEvents[0].size()),
    static_cast<int>(fOnEvents[1].size()),
    static_cast<int>(fOnEvents[2].size()),
    static_cast<int>(fOnEvents[3].size())
  };

  TLegend *leg = new TLegend(0.9,0.9,0.999,0.999);

  TGraph *AvgSigGraph[4];
  TGraph *FracDiffGraph[4];

  TMultiGraph *allAvgSigGraph = new TMultiGraph();
  TMultiGraph *allFracDiffGraph = new TMultiGraph();

  for(int i=0;i<4;i++){
    AvgSigGraph[i] = new TGraph();
    FracDiffGraph[i] = new TGraph();

    AvgSigGraph[i]->GetXaxis()->SetTitle("Number of flashes");
    AvgSigGraph[i]->GetYaxis()->SetTitle("Amplitude setting");

    FracDiffGraph[i]->GetXaxis()->SetTitle("Number of flashes");
    FracDiffGraph[i]->GetYaxis()->SetTitle("Fractional difference ((ADC-PrevADC)/PrevADC)");

    for(uint j=1;j<fAvgSig[i].size();j++){
      AvgSigGraph[i]->SetPoint(j,j,(fAvgSig[i][j]));
    }

    for(uint j=1;j<fFracDiff[i].size();j++){
      FracDiffGraph[i]->SetPoint(j,j,(fFracDiff[i][j]));
    }

    AvgSigGraph[i]->SetLineColor(i+1);
    AvgSigGraph[i]->SetMarkerColor(i+1);
    AvgSigGraph[i]->SetTitle(Form("Box-%d",i));

    FracDiffGraph[i]->SetLineColor(i+1);
    FracDiffGraph[i]->SetMarkerColor(i+1);
    FracDiffGraph[i]->SetTitle(Form("Box-%d",i));

    AvgSigGraph[i]->SetMaximum(1024);
    AvgSigGraph[i]->SetMinimum(0);

    allAvgSigGraph->Add(AvgSigGraph[i],"PL");
    allFracDiffGraph->Add(FracDiffGraph[i],"PL");

    leg->AddEntry(AvgSigGraph[i],Form("Box-%d",i+1),"L");
  }

  //Add in boxes
  TLine *mylinesstart[4][200];
  TLine *mylinesend[4][200];
  for(int i=0;i<4;i++){
    COMETLog("Box: " << i+1);

    int start=0;
    int end=0;
    int total=0;
    for(int j=0;j<length[i];j++){
      if(j==0){
  start=0;
  end=fOnEvents[i][j];
  total=end+fOffEvents[i][j];
      }
      else{
  start=total;
  end=total+fOnEvents[i][j];
  total=end+fOffEvents[i][j];
      }

      COMETLog("Aplitude changes (in #flashes): start=" << start << ", end=" << end << ", On=" <<fOnEvents[i][j] << ", Off=" << fOffEvents[i][j]);

      mylinesstart[i][j] = new TLine(start,0,start,1024);
      mylinesend[i][j] = new TLine(end,0,end,1024);

      mylinesstart[i][j]->SetLineColor(kGreen);
      mylinesend[i][j]->SetLineColor(kRed);
    }
  }

  TCanvas *c2 = new TCanvas("AllAvgSig","AllAvgSig",10,10,1000,700);
  //allAvgSigGraph->GetXaxis()->SetTitle("Number of flashes");
  //allAvgSigGraph->GetYaxis()->SetTitle("Amplitude setting");
  allAvgSigGraph->Draw("AP");
  leg->Draw("SAME");

  if(c2)
    c2->Write();

  TCanvas *c3 = new TCanvas("AllFracDiff","AllFracDiff",10,10,1000,700);
  //allFracDiffGraph->GetXaxis()->SetTitle("Number of flashes");
  //allFracDiffGraph->GetYaxis()->SetTitle("Fractional difference ((ADC-PrevADC)/PrevADC)");
  allFracDiffGraph->Draw("AP");
  leg->Draw("SAME");

  if(c3)
    c3->Write();

  TCanvas *c5 = new TCanvas("AvgSig","AvgSig",10,10,1000,700);
  c5->Divide(2,2);
  for(int i=0;i<4;i++){
    c5->cd(i+1);
    AvgSigGraph[i]->Draw("ALP");

    /*
    for(int j=0;j<length[i];j++){
      mylinesstart[i][j]->Draw("SAME");
      mylinesend[i][j]->Draw("SAME");
    }
    */
  }

  if(c5)
    c5->Write();

  TCanvas *c6 = new TCanvas("FracDiff","FracDiff",10,10,1000,700);
  c6->Divide(2,2);
  for(int i=0;i<4;i++){
    c6->cd(i+1);
    FracDiffGraph[i]->Draw("ALP");
  }

  if(c6)
    c6->Write();

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

  if(fOutputTree){
    fOutputTree->Write();
    fOutputTree->Print();
  }

  fOutputFile->Write();
  fOutputFile->Close();
}

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

bool IP0dLi::SetOption(std::string option,std::string value="") {
  if (option == "save-channel-histos")
    fSaveChannelHistos = true;
  else if (option == "outputfile")
    fOutputFileName = value;

  return true;
}

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

//Determine which pulser box
int IP0dLi::DetLIBox(int p0dule,int layer){

  int boxnumber=-1;

  if(p0dule<20 && layer==0) boxnumber=3;
  if(p0dule<20 && layer==1) boxnumber=1;
  if(p0dule>=20 && layer==0) boxnumber=2;
  if(p0dule>=20 && layer==1) boxnumber=0;

  return boxnumber;
}

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

bool IP0dLi::operator () (COMET::ICOMETEvent& event) {
  COMETVerbose( "TriggerBits = " << event.GetHeader().GetTriggerBits() );

  //P0D electronic addressing to physical addressing (channel,layer,p0dule)
  const COMET::IP0DChannelMap& P0DChanMap = fDB.P0D();

  //Get Run Context
  fRunId = event.GetContext().GetRun();
  const int eventId = event.GetContext().GetEvent();

  //Fill trigger histogram (all triggers)
  fTriggerType->Fill(event.GetHeader().GetTriggerBits());

  //Grab the raw event data banks
  IHandle<ICOMETRawEvent> rawEvt(event.GetRaw());
  assert(rawEvt);

  //Access the P0DLI current setting value in the ODB (Written at the beginning of each file, trigger =0)
  if(event.GetHeader().GetTriggerBits()==0){

    COMETLog("Reading ODB dump at the begining of the file...");
    //std::cout << "----> Reading ODB dump at the begining of the file..." << std::endl;

    COMET::IHandle<COMET::IXmlBank> xmlBank(rawEvt->GetMidasBank<COMET::IXmlBank>("XMLD"));
    if(!xmlBank)
      COMETLog("Cannot find an XMLD bank in the first record.");
    else{
      const COMET::IXmlOdb* xmlOdb = xmlBank->GetXmlOdb();

      if(!xmlOdb)
  COMETLog("Cannot parse the XMLD bank in the first record.");
      else{
  fLIState = xmlOdb->odbReadInt("/STATE/P0D/P0DLI/CurrentSetting",-100);

  if(fLIState==-100)
    COMETLog("Cannot find the CurrentSetting in XML bank of the first record.");
  else {
    COMETLog("LI setting according to ODB = " << fLIState);

    for(int box=0;box<4;box++){
      if((fLIStartAmp+fLISetNum[box])%10!=fLIState && fLIFlashNum[box]>1){
        COMETLog("Warning LI setting according to code and ODB do not match!! For box "<< box <<".");
        COMETLog("ODB= " << fLIState << " Code= " << (fLIStartAmp+fLISetNum[box])%10);
      }
    }
  }
      }
    }
  }// end of if(event.GetHeader().GetTriggerBits()==0){

  //Setup output file and LI amplitude state at the begining of the fill
  if(fStart==true){
    COMETLog("LI setting according to code = " << fLIState);

    fOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE");
    fOutputTree = new TTree("LI","LI");

    fOutputTree->Branch("RunNumber",&fRunId,"RunNumber/I");
    fOutputTree->Branch("HighADCMean",&fHighADCMean,"HighADCMean/F");
    fOutputTree->Branch("HighADCRMS",&fHighADCRMS,"HighADCRMS/F");
    fOutputTree->Branch("LowADCMean",&fLowADCMean,"LowADCMean/F");
    fOutputTree->Branch("LowADCRMS",&fLowADCRMS,"LowADCRMS/F");
    fOutputTree->Branch("TDCMean",&fTDCMean,"TDCMean/F");
    fOutputTree->Branch("TDCRMS",&fTDCRMS,"TDCRMS/F");
    fOutputTree->Branch("TDCMaxBin",&fTDCMaxBin,"TDCMaxBin/F");
    fOutputTree->Branch("RMM",&fRMM,"RMM/I");
    fOutputTree->Branch("TFB",&fTFB,"TFB/I");
    fOutputTree->Branch("TripT",&fTripT,"TripT/I");
    fOutputTree->Branch("Channel",&fChannel,"Channel/I");
    fOutputTree->Branch("ChannelId",&fChannelId,"ChannelId/i");
    fOutputTree->Branch("P0dule",&fP0dule,"P0dule/I");
    fOutputTree->Branch("Bar",&fBar,"Bar/I");
    fOutputTree->Branch("Layer",&fLayer,"Layer/I");
    fOutputTree->Branch("Pulser",&fPulser,"Pulser/I");
    fOutputTree->Branch("Amplitude",&fAmplitude,"Amplitude/I");
    fOutputTree->Branch("Flashes",&fFlashes,"Flashes/I");
    fOutputTree->Branch("Time",&fTime,"Time/D");

    if(fSaveChannelHistos){
      fOutputTree->Branch("fHighADC","TH1D",&fHighADC,32000,0);
      fOutputTree->Branch("fLowADC","TH1D",&fLowADC,32000,0);
      fOutputTree->Branch("fTDC","TH1D",&fTDC,32000,0);
    }

    fLIStartAmp = fLIState;
    fStart=false;
  }//end of if(fStart==true)

  //Sum of the signal (ADC) on each of the boxes
  float summation[4]= {0,0,0,0};

  //Only look at trigger type 8 (P0DLED)
  if(event.GetHeader().GetTriggerBits()==8){

    fLIEventNum++;

    fEventTime.push_back(event.GetContext().GetTimeStamp());

    //Increment box flash counter by 1 and check for amplitude change conditions
    for(int box=0;box<4;box++){
      fLIFlashNum[box]++;

      //Setup histograms for each amplitude
      if(fLIADCSwitch[box]==true){
  fCurrentTime[box]=event.GetContext().GetTimeStamp();
  fLIFlashNum[box]=1;
  //Output the current amplitude setting
  COMETLog("Amplitude according to code " << fLIAmps[box][(fLIStartAmp+fLISetNum[box])%10]<< " : corresponds to ODB array element "
     <<(fLIStartAmp+fLISetNum[box])%10
     << " for pulser box " <<box);

  //Vector containing the timestamp of the first event in the amplitude setting
  fGroupTime[box].push_back(event.GetContext().GetTimeStamp());

  //Raw high gain ADC signal for all 10400 channels
  fBarHighADCMap[box] = new TH2I(Form("HighADCMap-Box-%d",box),
          Form("HighADCMap-Box-%d",box),
          6002, -2, 6000, 1050, 0, 1050);
  fBarHighADCMap[box]->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)");
  fBarHighADCMap[box]->GetYaxis()->SetTitle("High Gain ADC");

  //Raw low gain ADC signal for all 10400 channels
  fBarLowADCMap[box] = new TH2I(Form("LowADCMap-Box-%d",box),
         Form("LowADCMap-Box-%d",box),
         6002, -2, 6000, 1050, 0, 1050);
  fBarLowADCMap[box]->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)");
  fBarLowADCMap[box]->GetYaxis()->SetTitle("Low Gain ADC");

  fBarTDCMap[box] = new TH2F(Form("TDCMap-Box-%d",box),Form("TDC Map-Box-%d",box), 6002, -2, 6000, 800, 0, 800);
  fBarTDCMap[box]->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)");
  fBarTDCMap[box]->GetYaxis()->SetTitle("TDC (2.5ns per bin)");

  //Boolean switch to indicate when to change amplitudes
  fLIADCSwitch[box] = false;
      }
    }
    //Timing information from the MCM bank
    int trigtime=0;
    //int timesincebusy_mcm=0;
    //int evttime=0;

    IHandle<IMCMBank> mcmInfo;
    while(mcmInfo = rawEvt->GetMidasBank<COMET::IMCMBank>("", mcmInfo)){
      std::stringstream rr;
      rr << mcmInfo->GetName();
      std::string name(rr.str());
      std::string p0d="P";
      if(name[0]!=p0d) continue;
      trigtime = mcmInfo->GetUnixTimeSSecTrig();
      //timesincebusy_mcm = mcmInfo->GetTimeSinceBusy();
      //evttime = mcmInfo->GetEventTime();
    }

    //Data bank that holds all raw ADC/TDC information for all channels
    IHandle<ITripTDigitBank> triptBank;
    while (triptBank = rawEvt->GetMidasBank<ITripTDigitBank>("",triptBank) ) {

      if (triptBank->GetSubDetector() != 1){
  // Not a P0D Bank, ignore it
  continue;
      }

      IMidasTripTDigitItr itr(triptBank->GetMidasTripTDigitItr());
      while (!itr.EOD() ) {
  IMidasTripTDigit digit(itr.Get());

  ITFBChannelId tfbChannelId(digit.GetChannelId());

  const int capacitor = digit.GetIntegrationNum();
  const int time_tfb = digit.GetTime();

  COMET::IGeometryId geometryId;
  if (!fDB.GetGeometryId(geometryId, tfbChannelId)) {
    continue;
  }

  int p0dule, layer, bar;
  if(!P0DChanMap.GeomIDToLocation(p0dule, layer, bar, geometryId)) {
    continue;
  }

  int position = p0dule*140+50+bar;
  int timeRelToTrig_tfb = (time_tfb)%400000000-trigtime*4;

  int box = DetLIBox(p0dule,layer);

  //Fill raw ADC maps
  if(capacitor==2){
    fBarHighADCMap[box]->Fill(position, digit.GetHighGainADC());
    fBarLowADCMap[box]->Fill(position, digit.GetLowGainADC());

    //Fill timing histograms
    fBarTDCMap[box]->Fill(position,timeRelToTrig_tfb);
    summation[box]+=digit.GetLowGainADC();
  }
      }
    }

    /*
    for(int i=0;i<4;i++){
      std::cout << summation[i] << std::endl;
      std::cout << fPrevAvg[i] << std::endl;
    }
    */

    //int modifier=500;

    //Check each box for amplitude change condition
    for(int box=0;box<4;box++){

      //Need 2 events to calculate a proper difference
      //std::cout << "Flashes " << fLIFlashNum[box] << " for box " << box << std::endl;
      if(fLIEventNum<2) continue;

      int channels=0;
      int p0dules=0;
      int layer=0;
      int startp0d=0;
      int endp0d=0;

      if(box==0||box==1){
  channels=134;
  p0dules=20;
  layer=1;
      }
      if(box==2||box==3){
  channels=126;
  p0dules=20;
  layer=0;
      }
      if(box==0||box==2){
  startp0d=20;
  endp0d=40;
      }
      if(box==1||box==3){
  startp0d=0;
  endp0d=20;
      }

      float fracdiff = ((summation[box]/(channels*p0dules)) - fPrevAvg[box])/fPrevAvg[box];
      fFracDiff[box].push_back(fracdiff);

      //      std::cout << "Fracdiff = " << fracdiff << " for box " << box << std::endl;

      // Readout histograms and input values into tree
      // 10% value will miss the shift from amplitude 100 to off and off to 100
      //     which will add ~1second * trigger rate *2 events to the amplitude 100 histogram
      //     need to test how this affects the data.

      //if( (TMath::Abs(fracdiff) > fFracDiffCut) && (fLIFlashNum[box] > 100) ){
      if( (TMath::Abs(fracdiff) > fFracDiffCut)){ // This should work when processing data subrun-by-subrun
  COMETLog("Entering end of amplitude for box " << box << " and event " << eventId << " with flashes " << fLIFlashNum[box] );

  fAmplitude = fLIAmps[box][(fLIStartAmp+fLISetNum[box])%10];
  fLayer = layer;
  fPulser = box+1;
  fFlashes = fLIFlashNum[box];
  fTime = fCurrentTime[box];

  for(int p0d=startp0d;p0d<endp0d;p0d++){
    for(int chan=0;chan<channels;chan++){
      int bin = 53+140*p0d+chan;

      fHighADCMean = fBarHighADCMap[box]->ProjectionY("t_py",bin,bin)->GetMean();
      gDirectory->Delete("t_py");

      fHighADCRMS = fBarHighADCMap[box]->ProjectionY("t_py",bin,bin)->GetRMS();
      gDirectory->Delete("t_py");

      fLowADCMean = fBarLowADCMap[box]->ProjectionY("t_py",bin,bin)->GetMean();
      gDirectory->Delete("t_py");

      fLowADCRMS = fBarLowADCMap[box]->ProjectionY("t_py",bin,bin)->GetRMS();
      gDirectory->Delete("t_py");

      fTDCMean = fBarTDCMap[box]->ProjectionY("t_py",bin,bin)->GetMean();
      gDirectory->Delete("t_py");

      fTDCRMS = fBarTDCMap[box]->ProjectionY("t_py",bin,bin)->GetRMS();
      gDirectory->Delete("t_py");

      fTDCMaxBin = fBarTDCMap[box]->ProjectionY("t_py",bin,bin)->GetMaximumBin();
      gDirectory->Delete("t_py");

      //////////
      // Histograms for every channel
      if(fSaveChannelHistos){
        fHighADC = fBarHighADCMap[box]->ProjectionY("fHighADC",bin,bin);
        fLowADC = fBarLowADCMap[box]->ProjectionY("fLowADC",bin,bin);
        fTDC = fBarTDCMap[box]->ProjectionY("fTDC",bin,bin);
      }

      /////////
      //Physical geometry

      fP0dule = p0d;
      fBar = chan;

      //Electronic geometry
      COMET::IGeometryId mygeom = COMET::GeomId::P0D::Bar(p0d,layer,chan);
      ITFBChannelId mychanid;
      if(!fDB.GetChannelId(mychanid,mygeom)) continue;

      fRMM = mychanid.GetRMM();
      fTFB = mychanid.GetTFB();
      fTripT = mychanid.GetTripChip();
      fChannel = mychanid.GetChannel();
      fChannelId = mychanid.AsUInt();

      //std::cout << "fChannelId=" << fChannelId << std::endl;

      //Fill the output tree once we have the all channels read
      fOutputTree->Fill();
    }
  }
  delete fBarHighADCMap[box];
  delete fBarLowADCMap[box];
  delete fBarTDCMap[box];

  fLISetNum[box]+= 1;

  if(fOnEvents[box].size()> fOffEvents[box].size()) fOffEvents[box].push_back(0);
  fOnEvents[box].push_back(fLIFlashNum[box]);
  fLIADCSwitch[box]=true;
  fPrevFlash[box]=fLIFlashNum[box];
  fLIFlashNum[box]=-1;

      }
    }

    //Save current average for the next event
    for(int i=0;i<4;i++){
      int channels;
      if(i<2){
  channels=134*20;
      }
      else{
  channels=126*20;
      }

      fPrevAvg[i] = summation[i]/channels;
      fAvgSig[i].push_back(fPrevAvg[i]);
    }
  }
  return false;
}

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

int main(int argc, char **argv) {
  IP0dLi userCode;
  cometEventLoop(argc,argv,userCode);
}