/*
 * mip_histogram_fitter.C
 * Original Author: Matt Walker
 * Update Author: J. Kevin Adkins, University of Kentucky
 * Updates: Remove old and obsolete functions and code,
 * make code more modular so it's easier to make changes
 * to fits independent of drawing tower fits, etc.
 * Fitting algorithm completely overhauled to now calculate
 * the mean of the fit over the range [0,100]. This is the value
 * we calculate the MIP relative gains with. Also, now using fit
 * to calculate RMS, and subsequently the error on the mean
 * Most recent update: Feb. 25, 2015
 */

#include <iostream>
#include <fstream>
#include <vector>
#include <sstream>
using namespace std;

#include "TFile.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TPDF.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TString.h"
#include "TLine.h"
#include "TF1Convolution.h"
#include "TROOT.h"

//void drawTower(TH1D*, Int_t, Int_t, Float_t, Float_t);
//Bool_t isBadTower(Int_t);
stringstream ss;

void mip_histogram_fitter(const Char_t* mipRootfile = "./0506/mipcombined22_0426.root", const Char_t* psName = "./0506/mips_conv.pdf", const Char_t* rootFilename = "./0506/towerMipFits_conv_f10_1.root", const Char_t* mipOutputName = "./0506/mips_conv_f10_1.gains")
{
  gROOT->Macro("loadMuDst.C");
  //StEmcGeom *mEmcGeom = StEmcGeom::instance("bemc");

  const Int_t nTowers = 30;
  Int_t towerStatus[nTowers];
  Double_t fitMean[nTowers], fitMeanError[nTowers];
  Double_t fitParams[5];

  cout << "Input Filename: " << mipRootfile << endl;
  cout << "Plot Filename: " << psName << endl;
  cout << "Gain Filename: " << mipOutputName << endl;

  gStyle->SetCanvasColor(10);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetStatColor(10);
  gStyle->SetOptTitle(0);
  gStyle->SetOptDate(0);
  gStyle->SetOptFit(111);
  gStyle->SetOptStat("e");
  gStyle->SetPalette(1);

  /* Declare pointers for histograms before output file. Instantiate after,
     this allows us to use outfile->Write() to write everything at once */
  TH1D *towerHisto[nTowers];
  //TH1F *htowerTrackDelta[nTowers];
  //TH2F *etaPhiMean, *etaPhiStatus, *etaPhiRatio;

  TFile* outfile = new TFile(rootFilename,"RECREATE");
  Float_t pi = TMath::Pi();

  // Instantiate histograms
  //etaPhiMean = new TH2F("etaPhiMean","Fit Means", 40, -1., 1., 120, -pi, pi);
  //etaPhiStatus = new TH2F("etaPhiStatus","Status Codes", 40, -1., 1., 120, -pi, pi);
  //etaPhiRatio = new TH2F("etaPhiMeanRatio","Mean Ratios", 40, -1., 1., 120, -pi, pi);
  for (Int_t i = 0; i < nTowers; ++i)
    towerHisto[i] = new TH1D(Form("towerHisto_%i",i+1),Form("Tower %i",i+1),250,-50.5,199.5);

  TFile *rootfile = new TFile(mipRootfile);
  for(Int_t y = 0; y < nTowers; ++y)
    towerHisto[y]->Add((TH1D*)rootfile->Get(Form("tower_histo_%i",y+1)));

  rootfile->Close();
  if(rootfile) delete rootfile;

  /****************** Begin MIP Fitting for each Tower ******************/
  TF1 *gaussian[nTowers];
  TF1 *landau[nTowers];
  TF1Convolution *towerFit_Initial[nTowers];

  TF1 *towerFit[nTowers];
  TString fitName;
  TString gausName;
  TString landauName;

  //TH1F *hRMS = new TH1F("hRMS","RMS",)

  for(Int_t iTow = 0; iTow < nTowers; ++iTow){
    Int_t softId = iTow+1;
    //Float_t mTowerEta = 0., mTowerPhi = 0.;
    //mEmcGeom->getEta(softId,mTowerEta);
    //mEmcGeom->getPhi(softId,mTowerPhi);
    //if(softId >= 3381 && softId <= 3540) continue;
    //
    //if (mTowerEta >= 0.1 && mTowerEta <= 0.15 && mTowerPhi <= -2.0 && mTowerPhi >= -2.2) cout << "TowerID: " << softId << endl;


    towerStatus[iTow] = 0;
    fitMean[iTow] = 0.;
    fitMeanError[iTow] = 0.;

    Double_t fitLow = towerHisto[iTow]->GetMean() - towerHisto[iTow]->GetRMS();
    Double_t fitHigh = towerHisto[iTow]->GetMean() + towerHisto[iTow]->GetRMS()*2.;
    towerHisto[iTow]->GetXaxis()->SetRangeUser(6.,100.);

    cout << "Now fitting softId: " << softId << endl;
    ss.str("");
    ss << softId;
    fitName = (TString)("towerFit_" + ss.str());
    gausName = (TString)("gausFit_" + ss.str());
    landauName = (TString)("landauFit_" + ss.str());



    //Double_t fitLow = 0;
    //Double_t fitHigh = 100;

    landau[iTow] = new TF1(landauName,"landau",fitLow,fitHigh);
    landau[iTow]->SetParameters(towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetRMS());
    //towerHisto[iTow]->Fit(landau[iTow],"RQ");
    //Double_t mpv = landau[iTow]->GetParameter(1);
    //Double_t rms_landau = landau[iTow]->GetParameter(2);
    //Double_t tail_range_low = 10;
    //Double_t tail_range_high = 40;
    gaussian[iTow] = new TF1(gausName,"gaus",fitLow,fitHigh);
    gaussian[iTow]->SetParameters(towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetMean(),towerHisto[iTow]->GetRMS());
    //towerHisto[iTow]->Fit(gaussian[iTow],"RQ");

    Double_t area_gaus_tails = gaussian[iTow]->Integral(fitLow,fitHigh);

    //cout << "area_gaus_tails " << softId << " " << area_gaus_tails << endl;

    Double_t area_tower_histograms = towerHisto[iTow]->Integral();

    //cout << "area_tower_histograms " << softId << " " << area_gaus_tails << endl;

    Double_t gaussian_scale_parameter = area_tower_histograms / area_gaus_tails;

    cout << "Tower: " << softId << " " << "Gaussian Integral: " << area_gaus_tails << " " << "Histogram Integral: " << area_tower_histograms << " " << "Scale Parameter: " << gaussian_scale_parameter << endl;

    towerFit_Initial[iTow] = new TF1Convolution(landau[iTow],gaussian[iTow],fitLow,fitHigh,true);
    //towerFit_Initial->SetRange(fitLow,fitHigh);
    towerFit_Initial[iTow]->SetNofPointsFFT(75000);

    towerFit[iTow] = new TF1(fitName, *towerFit_Initial ,fitLow,fitHigh,5);
    //towerFit[iTow]->SetNpx(90000);
    cout << "Name of Fit " << towerFit[iTow]->GetName() << endl;
    //towerFit[iTow] = new TF1(fitName,"[0]*Gaus(x,[1],[2])*Landau(x,[3],[4])",fitLow,fitHigh);

    // Fit Parameters
    towerFit[iTow]->SetParameter(0,towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin())); //Amplitude
    towerFit[iTow]->SetParameter(1,towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin())); //Landau MPV
    towerFit[iTow]->SetParameter(2,towerHisto[iTow]->GetRMS()); //Landau RMS
    towerFit[iTow]->SetParameter(3,gaussian_scale_parameter); //Gaussian Scale Parameter
    towerFit[iTow]->SetParameter(4,towerHisto[iTow]->GetRMS()); //Gaussian RMS

    //towerFit[iTow]->SetParLimits(2,0,1000);

    towerFit[iTow]->SetLineColor(kBlue);
    towerFit[iTow]->SetLineWidth(1.5);

    // If entries less than 25 we can't fit the tower, so we skip it

    if(towerHisto[iTow]->Integral(56,156) < 25){
      towerStatus[iTow] = 13;
      //int isGood = 0;
      //if (towerStatus[iTow] == 1) isGood = 1;
      //else isGood = 0;
      //etaPhiStatus->Fill(mTowerEta,mTowerPhi,towerStatus[iTow]);
      //etaPhiStatus->Fill(mTowerEta,mTowerPhi,isGood);
      continue;
    }

    towerHisto[iTow]->Fit(towerFit[iTow],"RQ");
    towerStatus[iTow]+=1;
    fitMean[iTow] = towerFit[iTow]->Mean(fitLow,fitHigh);

    // Get the fit parameters in an array
    Double_t fitParams[5] = {0.,0.,0.,0.,0.};
    fitParams[0] = towerFit[iTow]->GetParameter(0);
    fitParams[1] = towerFit[iTow]->GetParameter(1);
    fitParams[2] = towerFit[iTow]->GetParameter(2);
    fitParams[3] = towerFit[iTow]->GetParameter(3);
    fitParams[4] = towerFit[iTow]->GetParameter(4);
    //Calculate the variance of the fit
    Double_t fitVariance = towerFit[iTow]->Variance(0., 100., fitParams);
    Double_t fitSigma = sqrt(fitVariance);
    Double_t entriesInFit = towerHisto[iTow]->Integral(51,151); // Bin 51 centered over zero
    cout << "entriesInFit " << softId << " " << entriesInFit << endl;

    //Set the status codes of the towers
  //  if(fitMean[iTow] < 6) towerStatus[iTow] += 10; //Bad mean
/*
    if((softId)%20==0){
      if (abs(fitMean[iTow]-towerHisto[iTow]->GetMean()) > 15)
	towerStatus[iTow]+=8;
    }
    else{
      if(abs(fitMean[iTow]-towerHisto[iTow]->GetMean()) > 10)
	towerStatus[iTow]+=8;
    }
*/
    //Check for hard coded bad towers

    //if(isBadTower(softId)) towerStatus[iTow] = 7;
    //if(softId == 4) towerStatus[iTow] = 1;
    //if(softId >= 3381 && softId <= 3540) towerStatus[iTow] = 9;

    //Calculate error on mean for good statuses
    if (towerStatus[iTow] == 1)fitMeanError[iTow] = fitSigma/sqrt(entriesInFit);


    // Fill eta/phi for good statuses
    /*
    if (towerStatus[iTow] == 1){
      etaPhiMean->Fill(mTowerEta,mTowerPhi,fitMean[iTow]);
      etaPhiRatio->Fill(mTowerEta,mTowerPhi,fitMean[iTow]/towerHisto[iTow]->GetMean());
    }*/
    //etaPhiStatus->Fill(mTowerEta,mTowerPhi,towerStatus[iTow]);
  }// Towers loop
  /******************** End MIP Fitting for each Tower ********************/

  /******************* Begin Drawing the Fitted Spectra *******************/
  TPDF *ps = new TPDF(psName);
  TCanvas *canvas = new TCanvas("can","can1",100,100,600.,800.);
  Int_t pad;

  cout << endl << "Begin Tower Drawing Loop" << endl;
  for (Int_t iTow = 0; iTow < nTowers; ++iTow){
    gStyle->SetOptFit(1);
    Int_t softId = iTow + 1;
    //Float_t mTowerEta = 0., mTowerPhi = 0.;
    //mEmcGeom->getEta(softId,mTowerEta);
    //mEmcGeom->getPhi(softId,mTowerPhi);
    //if(softId >= 3381 && softId <= 3540) continue;


    cout << "Drawing tower " << softId << endl;
    if(iTow%20 == 0){
      canvas->Update();
      ps->NewPage();
      canvas->Clear();
      canvas->Divide(4,5);
      pad = 1;
    }

    canvas->cd(pad);
    if(towerStatus[iTow] > 1)
      towerHisto[iTow]->SetLineColor(kRed);

    // Draw the tower onto the pad
    //drawTower(towerHisto[iTow],softId,towerStatus[iTow],mTowerEta,mTowerPhi);
    towerHisto[iTow]->GetXaxis()->SetTitle("ADC");
    towerHisto[iTow]->GetXaxis()->CenterTitle();
    towerHisto[iTow]->Draw();

    Double_t lineMax = towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()) + 30;
    if(towerStatus[iTow] == 1){
      TLine *avgValueLine = new TLine(fitMean[iTow],0.,fitMean[iTow],lineMax);
      avgValueLine->SetLineColor(kRed);
      avgValueLine->SetLineWidth(1.25);
      avgValueLine->Draw("same");
    }
    pad++;
  }//End tower draw loop
/*
  canvas->Update();
  ps->NewPage();
  canvas->Clear();
  canvas->Divide(1,2);
  canvas->cd(1);
  etaPhiMean->GetXaxis()->SetTitle("#eta");
  etaPhiMean->GetXaxis()->CenterTitle();
  etaPhiMean->GetYaxis()->SetTitle("#phi");
  etaPhiMean->GetYaxis()->CenterTitle();
  etaPhiMean->GetZaxis()->SetRangeUser(0,45);
  etaPhiMean->Draw("colz");
  canvas->cd(2);
  //etaPhiStatus->GetZaxis()->SetRangeUser(0,20);
  cout << endl << etaPhiStatus->GetNbinsX();
  for (int ix = 1; ix <= etaPhiStatus->GetNbinsX(); ix++) {
      for (int iy = 1; iy <= etaPhiStatus->GetNbinsY(); iy++) {
          if (etaPhiStatus->GetBinContent(ix, iy) != 1) etaPhiStatus->SetBinContent(ix, iy, 0);
      }
  }
  etaPhiStatus->GetZaxis()->SetRangeUser(0,2);
  etaPhiStatus->GetXaxis()->SetTitle("#eta");
  etaPhiStatus->GetXaxis()->CenterTitle();
  etaPhiStatus->GetYaxis()->SetTitle("#phi");
  etaPhiStatus->GetYaxis()->CenterTitle();
  etaPhiStatus->Draw("colz");

  canvas->Update();
  ps->NewPage();
  canvas->Clear();
  canvas->Divide(1,2);
  canvas->cd(1);
  etaPhiRatio->GetXaxis()->SetTitle("#eta");
  etaPhiRatio->GetXaxis()->CenterTitle();
  etaPhiRatio->GetYaxis()->SetTitle("#phi");
  etaPhiRatio->GetYaxis()->CenterTitle();
  etaPhiRatio->Draw("colz");
  canvas->cd(2);
  TH2F *ratioClone = (TH2F*)etaPhiRatio->Clone();
  ratioClone->GetZaxis()->SetRangeUser(0.7,1.);
  ratioClone->GetXaxis()->SetTitle("#eta");
  ratioClone->GetXaxis()->CenterTitle();
  ratioClone->GetYaxis()->SetTitle("#phi");
  ratioClone->GetYaxis()->CenterTitle();
  ratioClone->Draw("colz");
  ps->Close();
*/

  cout << "Tower Drawing Loop Complete" << endl;


  /******************** End Drawing the Fitted Spectra ********************/

  /******************** Begin Output of MIP Fit Data ********************/
  ofstream mipOutput(mipOutputName);
  cout << "Now writing to Gains file: " << mipOutputName << endl;

  Int_t nGood = 0;
  Int_t nZero = 0;
  Int_t nBad = 0;
/*
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerFit[iTow]->GetName() << endl;
  cout << " " << endl;
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetMean() << endl;
  cout << " " << endl;
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetRMS() << endl;
  cout << " " << endl;
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()) << endl;
  cout << " " << endl;
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin()) << endl;
  cout << " " << endl;
  for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << fitMean[iTow]/towerHisto[iTow]->GetMean() << endl;
*/

  for(Int_t iTow = 0; iTow < nTowers; ++iTow){
    //cout << fitMeanError[iTow] << endl;
  //cout << "Inside of Tower loop to write to Gains file" << endl;


    Int_t softId = iTow + 1;
    Double_t mipAdc = fitMean[iTow];
    Double_t mipError = fitMeanError[iTow];
  //  Double_t mipRMS = towerFit[iTow]->GetParameter(2);
    Double_t histMean = towerHisto[iTow]->GetMean();
    Double_t histRMS = towerHisto[iTow]->GetRMS();
    Double_t histIP1 = towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin());
    Double_t histIP3 = towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin());
    Double_t Ratio = fitMean[iTow]/towerHisto[iTow]->GetMean();
    //cout << "Writing to Gains file" << endl;
    mipOutput << softId << " " << mipAdc << " " << mipError << " " << towerStatus[iTow] <<  endl;

    if(towerStatus[iTow]==1)nGood++;
    if(towerStatus[iTow]==13)nZero++;
    if(towerStatus[iTow]==7)nBad++;
  }
  cout << "Finished writing to Gains file: " << mipOutputName << endl;
  mipOutput.close();


  /********************* End Output of MIP Fit Data *********************/

  outfile->Write();
  outfile->Close();

  cout << "Finished tower fits" << endl << endl;
  cout << "Number of Good Statuses = " << nGood << endl;
  cout << "Number of Bad Statuses = " << nBad << endl;
  cout << "Number of Empty Towers = " << nZero << endl;


  // Delete geometry pointer
  //if (mEmcGeom) delete mEmcGeom;
}
  // Load StEmcGeom for tower eta & phi

 // End of main routine
/*
void drawTower(TH1D* histo, Int_t id, Int_t status, Float_t towerEta, Float_t towerPhi)
{
  // Set histogram options
  histo->GetXaxis()->SetTitle("ADC");
  histo->GetXaxis()->CenterTitle();
  histo->Draw(); // Used Sumw2() in histogram maker to calculate errors

  // Write softId and status code on plot
  Char_t towerTitle[100];
  Char_t etaTitle[100];
  Char_t phiTitle[100];
  sprintf(etaTitle,"#eta:%1.2f",towerEta);
  sprintf(phiTitle,"#phi:%1.2f",towerPhi);
  TLatex etaLatex;
  TLatex phiLatex;
  etaLatex.SetTextSize(0.1);
  phiLatex.SetTextSize(0.1);
  sprintf(towerTitle,"%i",id);
  TLatex titleLatex;
  titleLatex.SetTextSize(0.15);
  //if(status!=1) titleLatex.SetTextColor(kRed);  // Non - colorblind friendly color
  if(status!=1) titleLatex.SetTextColor(kOrange);
  phiLatex.DrawTextNDC(0.55,0.2,phiTitle);
  etaLatex.DrawTextNDC(0.55,0.33,etaTitle);
  titleLatex.DrawTextNDC(0.55,0.45,towerTitle);
  if(status!=1){
    Char_t towerCode[100];
    sprintf(towerCode,"%i",status);
    TLatex statusLatex;
    statusLatex.SetTextSize(0.15);
    //statusLatex.SetTextColor(kRed); // Non - colorblind friendly color
    statusLatex.SetTextColor(kOrange);
    statusLatex.DrawTextNDC(0.47,0.78,towerCode);
  }
}

Bool_t isBadTower(Int_t id) {
  //Empty towers = {31,113,134,166,257,426,509,541,793,817,832,844,846,956,989,1005,1012,1027,1057,1081,1122,1158,1159,1180,1189,1190,1202,1203,1204,1207,1214,1222,1224,1242
                // 1243,1244,1329,1341,1354,1369,1388,1434,1438,1443,1448,1533,1534,1575,1593,1594,1595,1597,1654,1668,1674,1679,1705,1807,1866,2019,2104,2111,2168,2214,2281
                // 2282,2283,2386,2392}
  //Hot towers = {365,375,533,939,945,1211,1223,1826,2032,2234,2300}

  Int_t badTowers[] = {34,106,139,160,245,266,267,286,287,365,368,375,504,514,533,561,562,615,616,633,637,638,649,650,653,657,671,673,674,740,760,
                      813,814,837,857,875,899,901,939,945,953,954,993,1026,1044,1045,1046,1048,1080,1100,1198,1199,1200,1207,1211,1219,1220,1223,
                      1237,1280,1303,1320,1337,1341,1348,1353,1354,1377,1380,1407,1447,1448,1530,1574,1575,1597,1612,1620,1654,1753,
                      1765,1766,1826,1856,1877,1878,2032,2040,2073,2077,2092,2093,2097,2219,2234,2300,2305,2529,2569,
                      2589,2590,2613,2700,2834,2890,2969,3005,3061,3070,3071,3381,3384,3387,3391,3392,3393,3394,3395,
                      3396,3397,3398,3403,3407,3408,3411,3413,3416,3419,3426,3430,3432,3436,3437,3442,3445,3447,3448,
                      3451,3452,3455,3456,3458,3460,3464,3468,3474,3476,3478,3481,3483,3484,3485,3486,3488,3489,3496,
                      3497,3498,3499,3500,3501,3502,3503,3504,3507,3509,3510,3515,3516,3519,3523,3526,3528,3531,3532,3534,3535,3538,3539,3678,3679,3718,3738,3794,4017,4018,4019,4057,4177,4223,
                      4339,4560,4657,4677,4678,4684};

  for (int i = 0; i < sizeof(badTowers); i++)
  {
    if (badTowers[i] == id) return true;
  }

  return false;
}*/