#include <unistd.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>

#include <fstream>


// find signal amplitude of each PMT at +- 5,4,3,2,1ns
// using a simple landau fit to the data from
// the projection of the 2d histogram APMvsDTRawXXX
// where XXX is the PMT number.
// The histograms are generated by the plugin TOFmon
// The root file is called tofmon.root the location of
// the file is selected by the valiable "choice" either in
// ./localdir/runXXXXX/ or int genroot/.....
//
// write data to file: tof_amps_run%d.dat
//


int DEBUG = 0;

float findpeak(TH1D *h){

  float locmax = 0;
  TH1D *hist = (TH1D*) h->Rebin(30,"hist");
  int nbins = hist->GetSize();
  for (int n=1;n<nbins-5;n++){
    hist->Fit("pol2","Q","QR",hist->GetBinCenter(n),hist->GetBinCenter(n+4));

    TF1 *f1 = hist->GetFunction("pol2");
    float quad = f1->GetParameter(2);
    if (quad > 0){
      hist->GetXaxis()->SetRange(n+4,nbins-2);
      int maxbin = hist->GetMaximumBin();
      locmax = hist->GetBinCenter(maxbin);
      //cout<<"Maximum Found at "<<locmax<<endl;
      break;
    }

  }
  return locmax;
}


void findamp(int RunNumber, int choice){
 char ofnam[128];
  if (!choice) {
    sprintf(ofnam,"localdir/run%d/tof_amps_run%d.dat",RunNumber,RunNumber);
  } else {
    sprintf(ofnam,"genroot/hists/%06d/tof_amps_run%d.dat",RunNumber,RunNumber);
  }
  
  ofstream OUTF(ofnam);

  char Rfnam[128]; // root file name
  if (!choice) {
    sprintf(Rfnam,"localdir/run%d/tofmon.root",RunNumber);
  } else {
    sprintf(Rfnam,"genroot/hists/%06d/tofmon_%06d.root",RunNumber,RunNumber);
  }
  TFile *RF = new TFile(Rfnam,"READ");
  TCanvas *c1 = new TCanvas("c1","PLOTS",900,600);

  for (int plane=0;plane<2;plane++){

    for (int paddle=0;paddle<44;paddle++){

      if ((paddle==21) || (paddle==22)){
	continue;
      }
      
      int idx[2];
      idx[0] = paddle;
      idx[1] = paddle + 44;

      for (int k=0;k<2;k++){

	char hnam[128];
	sprintf(hnam,"TOFmon/AMPvsDTRaw%d%d",plane,idx[k]);
	TH2F *h2d = (TH2F*)RF->Get(hnam);
        if (h2d == NULL){
          cout<<"ERROR not histogram:"<<hnam<<endl;
          return;
        }
	TH1D *h1[8] ;
	h1[0] = h2d->ProjectionY("h10",37,39); // projection -5ns
	h1[1] = h2d->ProjectionY("h11",40,41); // projection -4ns
	h1[2] = h2d->ProjectionY("h12",43,44); // projection -3ns
	h1[3] = h2d->ProjectionY("h13",45,46); // projection -2ns
	h1[4] = h2d->ProjectionY("h14",55,56); // projection  2ns
	h1[5] = h2d->ProjectionY("h15",57,58); // projection  3ns
	h1[6] = h2d->ProjectionY("h16",60,61); // projection  4ns
	h1[7] = h2d->ProjectionY("h17",62,64); // projection  5ns
	
	for (int n=0;n<8;n++) {

	  h1[n]->GetXaxis()->SetRangeUser(100., 4000.);
	  float max = findpeak(h1[n]);

	  //int mbin = h1[n]->GetMaximumBin();
	  //float max = h1[n]->GetBinCenter(mbin);
	  float width = max*0.15;
	  float integr = 5000.;
	  
	  float loli = max - 1.5*width;
	  float hili = max + 7*width;
	  
	  h1[n]->Fit("landau","Q","QR",loli,hili);
	  TF1 *f1 = h1[n]->GetFunction("landau");
	  
	  float MPV = f1->GetParameter(1);
	  float SIG = f1->GetParameter(2);

	  //cout<<MPV<<endl;
	  if ((MPV < 450.) || (paddle==20)){
	    max = h1[n]->GetBinCenter( h1[n]->GetMaximumBin());
	    width = max*0.15; 
	    loli = max - 1.5*width;
	    hili = max + 7*width;
	    h1[n]->Fit("landau","Q","QR",loli,hili);
	    f1 = h1[n]->GetFunction("landau");
	    MPV = f1->GetParameter(1);
	    SIG = f1->GetParameter(2);
	    //cout<<MPV<<endl;
	  }
	  if (MPV<300.){
	    max = 600.;
	    width = max*0.15; 
	    loli = max - 1.5*width;
	    hili = max + 7*width;
	    h1[n]->Fit("landau","Q","QR",loli,hili);
	    f1 = h1[n]->GetFunction("landau");
	    MPV = f1->GetParameter(1);
	    SIG = f1->GetParameter(2);
	    //cout<<MPV<<endl;
	  }
	  //cout<<MPV<<endl;

	  width = SIG; 
	  loli = MPV - 1.5*width;
	  hili = MPV + 7*width;
	  h1[n]->Fit("landau","Q","QR",loli,hili);
	  f1 = h1[n]->GetFunction("landau");
	  MPV = f1->GetParameter(1);
	  SIG = f1->GetParameter(2);
	  float dMPV = f1->GetParError(1);
	  
	  cout<<plane<<" "<<paddle<<" PMT "<<idx[k]<<"  "<<n<<" "<<MPV<<"  "<<dMPV<<"  "<<SIG<<endl;
	  OUTF<<plane<<" "<<paddle<<"     "<<idx[k]<<"  "<<n<<" "<<MPV<<"  "<<dMPV<<"  "<<SIG<<endl;
	  h1[n]->Draw();
	  c1->Modified();
	  if (DEBUG) {
	    sleep(1);
	  }
	  if (DEBUG>1){
	    h1[n]->GetXaxis()->SetTitle("ADC amplitude [arb.]");
	    c1->Modified();
	    c1->SaveAs("landau_fit.pdf");
	    getchar();
	  }
	}

      }
    }

  }
  OUTF.close();
  RF->Close();
}