#include <TH1F.h>
#include <TH2.h>
#include <TF1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <TLegend.h>
#include <TPaveText.h>
#include <TPostScript.h>
#include <TLatex.h>
#include <TFile.h>


#define theCut 5

#define binN 20000
#define binDown 0
#define binUp 20

#define  OutputDirectory "profiles/run"

void EnergyplotsNS(int _lumi, Int_t ieta=-13, TString _identifier="_aug_zb", bool _correct=0) {
  gROOT->SetStyle("Plain");
  // Prepare histograms
  TH1F* average = new TH1F("average","average",200,0,10);
  TH1F* average_fine = new TH1F("average_fine","average_fine",binN,binDown,binUp);
  TH1F* chan[65], *scale_ch[65];
  int _no_of_bad=0;
  
  float calibNS=1;
  
  if (ieta==-12)
    calibNS=1.11961;
  if (ieta==-13)
    calibNS=1.25654;
  if (ieta==-14)
    calibNS=0.921234;
  
//read list of bad channels
//   ifstream fileNames("/home/penning/interphi/badCellsFile.list");
//   char file[100];
  std::multimap<int, int> _badCellMap;
//   while(fileNames >> file)
//     {      
//       if(file[0] == '*') continue;
//       std::cout << "bad cell file: " << file << std::endl;
      
  fstream inExcl("excl.list",ios::in);
  while (inExcl.good())
    {
      int _ieta,_iphi ;
      inExcl >>  _ieta >> _iphi ;
      _badCellMap.insert(make_pair(_ieta,_iphi));
    }
  //    }
  
  assert(_badCellMap.size()>0);
  // Read information about stores
  TString inN="calib/calib";
  inN+=_identifier;
  inN+="/calib";
  inN+=_identifier;
  inN+=_lumi;
  inN+="_";
  inN+=ieta;
  inN+=".txt";
  float _calibs[65];
  _calibs[0]=-99;
  int _phi, _Bin, _N;
  float _calib, _calib_error;
  cout<<inN<<endl;
  fstream in(inN,ios::in);
  assert(in);
  string line;
  while (in.good()) {
    in>>_phi>>_Bin>>_N>>_calib>>_calib_error;
    _calibs[_phi]=_calib;
    if (in.eof()) break;
  }
  in.close();
  // Loop over the 64 iphis
  for (Int_t i=1; i<=64; i++) {
    TString nn("chani");
    TString cc("scale");
    nn+=_lumi;
    cc+=_lumi;
    nn+=ieta;
    cc+=ieta;
    nn+=i;
    chan[i]=new TH1F(nn,nn,binN,binDown,binUp);
    TString n(OutputDirectory);
    n+=_identifier;
    n+="/profile";
    n+=ieta;
    n+="_";
    n+=i;
    n+=".txt";
    // loop over data
    Int_t cnt=0;
    Float_t E;
    int _lumiBin;
    fstream rein(n,ios::in);
    assert(rein);
    bool _badSet=0;
    while (rein.good()) {
      rein>>E>>_lumiBin;
      if(_lumiBin!=_lumi) continue;
      if (rein.eof()) break;
      cnt++;
      if(_correct){
        chan[i]->Fill(E*_calibs[i]*calibNS);
        //        chan[iP]->Scale(_calibs[iP]);
        //        cout<<"correcting... "<<_calibs[iP]<<endl;
      }
      else
        chan[i]->Fill(E);		// Fill chan[i] always!  
      bool _isBad=0;
      typedef multimap<int, int>::const_iterator iter;
      pair<iter, iter> range=_badCellMap.equal_range(ieta);
      for(iter it=range.first; it!= range.second; ++it){
        if(it->second==i){
          _isBad=1;
          if (!_badSet)
            _no_of_bad++;
          _badSet=1;
        }
      }
      if (_isBad) continue;  
      average->Fill(E*calibNS);	// average and average_fine will be filled only if the channels are good!
      average_fine->Fill(E*calibNS);
    }
    rein.close();
    cout<<"Read "<<cnt<<" events for iphi "<<i<<endl;
  }
  average->Draw();

  // Prepare intercalibration
  Int_t lowBin=((Float_t)(theCut-binDown)/(binUp-binDown)*binN);
  Float_t avgExpectation=0;
  for (Int_t i=lowBin; i<=binN+1; i++) {
    avgExpectation+=average->GetBinContent(i);
  }

  float _divide=64-_no_of_bad;
  avgExpectation/=_divide;
  
  cout<<" _divide"<<_divide<<endl;
  TH1F *avfine = (TH1F*) average_fine->Clone();
  avfine->SetStats(0);
  avfine->GetXaxis()->SetRangeUser(binDown,7);
  avfine->Scale(1/_divide);

  // Find out the scale here
  for(Int_t sc = 1; sc <= 64; sc++ ) {
    //avfine->Sumw2();
    //chan[sc]->Sumw2();
    if(sc == 59 || sc == 60) continue;
    scale_ch[sc] = (TH1F *) avfine->Clone();
    scale_ch[sc]->Reset();
    scale_ch[sc]->SetLineColor(kGreen);
    scale_ch[sc]->Divide(avfine,chan[sc],1.,1.,"B");
    scale_ch[sc]->Rebin(400);
    scale_ch[sc]->SetTitle(Form("%d/%d",ieta,sc));
  }  


  TString psName="calib/calib";
  psName+=_identifier;
  psName+="/EProfile";
  psName+=_identifier;
  psName+="_lumi";
  psName+=_lumi;
  psName+="_";
  psName+=("ieta");
  psName+=ieta;
  TString cvname = psName;
  psName+=("NS.eps");
  avfine->SetLineColor(kRed); 
  TCanvas *cv1 = new TCanvas("c1", "c1", 6400,6400);
  cv1->Divide(8,8);

  avfine->Rebin(400);
  for (Int_t iP=1; iP<=64; iP++) {
    cv1->cd(iP);
    gPad->SetLogy(1); 
    chan[iP]->SetStats(0);
    chan[iP]->Rebin(400);
    chan[iP]->GetXaxis()->SetRangeUser(binDown,7);
    avfine->SetStats(0);
    avfine->GetXaxis()->SetRangeUser(binDown,7);
    avfine->Draw("e1");
    chan[iP]->Draw("e1same");
    ostringstream mystream;
    mystream.clear();
    mystream << iP;
    string s=mystream.str();
    TLegend* leg = new TLegend(0.5,0.75,0.88,0.99);
    leg->SetBorderSize(1);
    leg->SetFillColor(0);
    leg->AddEntry(avfine,"Avg","LP");
    leg->AddEntry(chan[1],"E","LP");
    leg->Draw();
    TPaveText *pt = new TPaveText(0.55,0.45,0.7,0.7,"brNDC");
    pt->SetTextSize(0.25);
    pt->SetFillColor(0);
    pt->SetBorderSize(0);
    pt->AddText(s.c_str());
    pt->Draw();
    typedef multimap<int, int>::const_iterator iter;
    pair<iter, iter> range=_badCellMap.equal_range(ieta);
    for(iter it=range.first; it!= range.second; ++it){
      if(it->second==iP){
        TPaveText *pt = new TPaveText(0.2,0.2,0.3,0.3,"brNDC");
        pt->SetTextSize(0.08);
        pt->SetFillColor(0);
        pt->SetBorderSize(0);
        TString _addText="excl.";
        pt->AddText(_addText);
        pt->Draw();
      }
    }

  }
  cv1->SaveAs(psName);
  cvname+="NS.root";
  TFile *outFile = new TFile(cvname,"RECREATE");
  avfine->Write();
  outFile->Close();
}