//-----------------------------------------------------------------------
// Fit ADC spectra with landau convoluted with Gaussian as a function of
// the distance from the pmt for both left and right side
// then fit these data points do determine attenuation length
// use adchists_run%d.root to get the 2d histograms x-position(TDC) vs. ADCIntegral
//
//
//-----------------------------------------------------------------------
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TText.h"
#include "TMinuit.h"

#include <iostream>
#include <fstream>

using namespace std;

int DEBUG = 99;
Double_t SNRPeak, SNRFWHM;
int RunNumber;

double XPos[2][100];
double YPos[2][100];
double dXPos[2][100];
double dYPos[2][100];
int NResults[2];
float Ratio;
float AllRatios[88];
float AllAmplitudes[88];

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
  
  //calculate chisquare
  Double_t chisq = 0;
  Double_t delta;
  Double_t amp;
  //first Left NResults[0],XPos[0],YPos[0],dXPos[0],dYPos[0]);
  const Int_t nbins = NResults[0] + NResults[1];

  for (int i=0; i<NResults[0]; i++) {
    
    if (XPos[0][i]<126.) {
      amp = par[0]*(TMath::Exp(-XPos[0][i]/par[1]));
    } else {
      amp = par[0]*(TMath::Exp(-(XPos[0][i]-126.)/par[2])) * (TMath::Exp(-126./par[1]));
    }
    
    delta = (YPos[0][i] - amp)/dYPos[0][i];
    chisq += delta*delta;
  }
  
  for (int i=0; i<NResults[1]; i++) {
    
    if (XPos[1][i]<126.) {
      amp = par[0]*(TMath::Exp(-XPos[0][i]/par[2]));
    } else {
      amp = par[0]*(TMath::Exp(-(XPos[0][i]-126.)/par[1])) * (TMath::Exp(-126./par[2]));
    }
    
    delta = (Ratio*YPos[1][i] - amp)/dYPos[1][i];
    chisq += delta*delta;
  }
  
  f = chisq;
}

Double_t langaufun(Double_t *x, Double_t *par) {

   //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation), 
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 100.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0]; 

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);

         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");
   
   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}


Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {

   // Seaches for the location (x value) at the maximum of the 
   // Landau-Gaussian convolute and its full width at half-maximum.
   //
   // The search is probably not very efficient, but it's a first try.

   Double_t p,x,fy,fxr,fxl;
   Double_t step;
   Double_t l,lold;
   Int_t i = 0;
   Int_t MAXCALLS = 10000;


   // Search for maximum

   p = params[1] - 0.1 * params[0];
   step = 0.05 * params[0];
   lold = -2.0;
   l    = -1.0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = langaufun(&x,params);
 
      if (l < lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-1);

   maxx = x;

   fy = l/2;


   // Search for right x location of fy

   p = maxx + params[0];
   step = params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-2);

   fxr = x;


   // Search for left x location of fy

   p = maxx - 0.5 * params[0];
   step = -params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;

   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-3);


   fxl = x;

   FWHM = fxr - fxl;
   return (0);
}

TF1*  langaus(TH1D *hist) {

   // Fitting histogram
   printf("Fitting...\n");

   // Setting fit range and start values
   Double_t FitRange[2];
   Double_t StartValues[4], LimitLow[4], LimitHigh[4], fp[4], fpe[4];
   FitRange[0] = 1000.;//0.5*hist->GetMean();
   FitRange[1] = 2.0*hist->GetMean();
   hist->GetXaxis()->SetRangeUser(1000.,20000.);
   int maxbin = hist->GetMaximumBin();
   double max  = hist->GetBinCenter(maxbin);
   FitRange[0] = max - 2000.;
   FitRange[1] = max + 6000.;
   if (FitRange[1]>23000.){
     FitRange[1]=23000.;
   }

   cout<<"Fitrange:  "<<FitRange[0]<<"  "<<FitRange[1]<<endl;

   StartValues[0] = 100; //Width
   StartValues[1] = max; //MP
   StartValues[2] = hist->GetEntries(); //Area
   StartValues[3] = 10.0; //GSigma

   LimitHigh[0] = 1500.0;
   LimitLow[0]  = 0.5;
   LimitHigh[1] = FitRange[1];
   LimitLow[1]  = FitRange[0];
   LimitHigh[2] = hist->GetEntries()*500.;
   LimitLow[2]  = hist->GetEntries()*0.01;
   LimitHigh[3] = 3000.0;
   LimitLow[3]  = 0.00001;

   Double_t chisqr;
   Int_t    ndf;
   TF1 *fitfunc = langaufit(hist,FitRange,StartValues,
			    LimitLow,LimitHigh,fp,fpe,&chisqr,&ndf);
   
   double w = fitfunc->GetParameter(0);
   double mp = fitfunc->GetParameter(1);
   if (FitRange[0]>(mp-w)){
     FitRange[0] = mp-w*1.5;
     if (FitRange[0]<0.){
       FitRange[0] = 10.;
     } 
     cout<<"Fitrange:  "<<FitRange[0]<<"  "<<FitRange[1]<<endl;
     fitfunc = langaufit(hist,FitRange,StartValues,
			 LimitLow,LimitHigh,fp,fpe,&chisqr,&ndf);
   } 
   
   langaupro(fp,SNRPeak,SNRFWHM);

   return fitfunc;
}

void atten1(int R){

  RunNumber = R;
  char fnam[128];
  if (R) {
    sprintf(fnam,"localdir/tofdata_run%d.root", RunNumber);
  } else {
    sprintf(fnam,"localdir/tofdata_all.root");
  }
  TFile *f = new TFile(fnam,"READ");
  //f->ls();

  TCanvas *c2 = new TCanvas("c2","Attenuation length",800.,600.);
  float a1[176]   ; // amplitude double attenuation fit
  float d1[176]   ; // first attenuation length 
  float a2[176]   ; // second amplitude not used 
  float d2[176]   ; // second attenuation length
  float Norm[176] ; // normalization factor
  float SingleAtten[176]; // attenuation length for single fit
  float ADC100S[176]; // attenuation at 100cm for single attenuation length fit
  float Chi2_D[176];  // chi2 of double attenuation length fit
  float Chi2_S[176];  // chi2 of single attenuation length fit
  for (int k=0;k<176;k++){
    a1[k]   = 0.;
    d1[k]   = 0.;
    a2[k]   = 0.;
    d2[k]   = 0.;
    Norm[k] = 0.; 
    SingleAtten[k] = 0.;
    ADC100S[k] = 0.;
    Chi2_D[176] = 0.;  
    Chi2_S[176] = 0.;  
  }

  for (int k=0;k<88;k++){
    AllAmplitudes[k] = 0;
    AllRatios[k] = 0;
  }

  double AttenuationData[176];
  double GainVal[176][5];

  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 + Plane*88;
      idx[1] = Paddle + Plane*88 + 44;
      
      

      for (int n=0;n<2;n++){
	
	char hnam[128];
	sprintf(hnam,"adc_hist%d",idx[n]);
	
	TH2D *h = (TH2D*)f->Get(hnam);
	if (h == NULL){
	  cout<<"ERROR not histogram:"<<hnam<<endl;
	  return;
	}
	int ev = h->GetEntries();
	if (!ev){
	  continue;
	}

	int NBins = h->GetNbinsY();

	int NFits = 0;

	for ( int k=1; k<NBins-2; k++) {

	  if  ((k==21) || (k==22)){
	    continue;
	  }
	  
	  TH1D *hp = h->ProjectionX("hp",k,k);
	  int nev = hp->GetEntries();

	  if (nev<150){
	    continue;
	  }

	  TH1D *hnew = (TH1D*)hp->Rebin(2,"hnew");
	  char newnam[128];
	  sprintf(newnam,"Energy Response Plane%d Paddle%d End%d xbin=%d",
		  Plane, Paddle, n, k);
	  hnew->SetTitle(newnam);

	  TF1* TheFitFunction = langaus(hnew);
	  
	  XPos[n][NFits] = h->GetYaxis()->GetBinCenter(k);
	  dXPos[n][NFits] = 0.;
	  YPos[n][NFits] = SNRPeak;
	  dYPos[n][NFits++] = SNRFWHM/12.;
	  
	  if (DEBUG>1) {
	    hnew->Draw();
	    TheFitFunction->Draw("lsame");
	    gPad->SetGrid();
	    gPad->Update();
	    char tn[128];
	    sprintf(tn,"plots/landau_fit_%d_%d_run%d.pdf",idx[n],k,RunNumber);
	    gPad->SaveAs(tn);
	    if (DEBUG==99){
	      //getchar();
	    }
	  }
	}
	NResults[n] = NFits;
      }
      
      if ((NResults[0]==0) || (NResults[1]==0)) {
	continue;
      }

      cout<<"Now get attenuation lengths"<<endl;
      TGraphErrors *grf[3];
      grf[0] = new TGraphErrors(NResults[0],XPos[0],YPos[0],dXPos[0],dYPos[0]);
      grf[1] = new TGraphErrors(NResults[1],XPos[1],YPos[1],dXPos[1],dYPos[1]);
      
      grf[0]->SetMarkerColor(4);
      grf[0]->SetMarkerStyle(21);
      
      grf[1]->SetMarkerColor(6);
      grf[1]->SetMarkerStyle(21);
      
      char grftit[128];
      sprintf(grftit,"Paddle %d, Plane %d",Paddle+1, Plane);
      TMultiGraph *mgr = new TMultiGraph("mgr",grftit);
      mgr->Add(grf[0]);
      mgr->Add(grf[1]);
      
      c2->Clear();
      c2->SetGrid();
      mgr->Draw("AP");
      gPad->SetGrid();
      //fitf->Draw();
      gPad->Update();
      //getchar();
      
      float cnt[2] = {0., 0.};
      float sums[2] = {0.,0.};

      for (int n=0; n<2; n++) {
	for (int k=0;k<NResults[n];k++){
	  if  ((XPos[n][k]>70) && (XPos[n][k]<180.)){
	    sums[n] += YPos[n][k];
	    cnt[n] +=1.;
	  } 
	}
	if (cnt[n]>0){
	  sums[n] /= cnt[n];
	}
      }

      Ratio = sums[0]/sums[1];
      AllRatios[Plane*44+Paddle] = Ratio;

      TF1 *f11,*f2,*f12;

      f11 = new TF1("f11","([0]*exp(-x/[1]) + [0]*exp(-x/[2]))", 0., 250.);
      f11->SetParameter(0,5000.);
      f11->SetParameter(1,50.);
      f11->SetParameter(2,500.);
      f11->SetParLimits(0, 1000., 50000.);
      f11->SetParLimits(1, 30., 150.);
      f11->SetParLimits(2, 150., 1000.);

      f12 = new TF1("f12","([0]*exp(-x/[1]) + [0]*exp(-x/[2]))", 0., 250.);
      f12->SetParameter(0,5000.);
      f12->SetParameter(1,50.);
      f12->SetParameter(2,500.);
      f12->SetParLimits(0, 1000., 50000.);
      f11->SetParLimits(1, 30., 150.);
      f12->SetParLimits(2, 150., 1000.);
      
      f2 = new TF1("f2","[0]*exp(-x/[1])", 0., 250.);
      f2->SetParameter(0,5000.);
      f2->SetParameter(1,500.);

      
      grf[0]->Fit("f11", "R", "0",0.,250. );
      grf[0]->GetFunction("f11")->SetLineColor(4);
      double NDF = (double)grf[0]->GetFunction("f11")->GetNDF();
      double chi2 = grf[0]->GetFunction("f11")->GetChisquare();
      Chi2_D[idx[0]] = chi2/NDF; 
      //f11->SetLineColor(4);
      //f11->Draw("same");
      mgr->GetXaxis()->SetTitle("Distance from PMT [cm]");
      mgr->GetYaxis()->SetTitle("Landau Peak [arb.]");
      gPad->Update();
      
      float att[2] = {0.,0.};
      a1[idx[0]] = f11->GetParameter(0);
      d1[idx[0]] = f11->GetParameter(1);
      d2[idx[0]] = f11->GetParameter(2);	
      grf[0]->Fit("f2", "R", "0",0.,250.);
      NDF = (double)grf[0]->GetFunction("f2")->GetNDF();
      chi2 = grf[0]->GetFunction("f2")->GetChisquare();
      Chi2_S[idx[0]] = chi2/NDF; 
      //f2->SetLineColor(7);
      //f2->Draw("same");
      att[0] = f2->GetParameter(1);
      a2[idx[0]] = f2->GetParameter(0);
      SingleAtten[idx[0]] = att[0];
      ADC100S[idx[0]] = f2->GetParameter(0)*exp(-100./att[0]);
      

      grf[1]->Fit("f12", "R", "0",0.,250. );
      grf[1]->GetFunction("f12")->SetLineColor(6);
      NDF = (double)grf[1]->GetFunction("f12")->GetNDF();
      chi2 = grf[1]->GetFunction("f12")->GetChisquare();
      Chi2_D[idx[1]] = chi2/NDF; 
      //grf[1]->GetFunction("f12")->Draw("same");
      //f12->SetLineColor(6);
      //f12->Draw("same");
      a1[idx[1]] = f12->GetParameter(0);
      d1[idx[1]] = f12->GetParameter(1);
      d2[idx[1]] = f12->GetParameter(2);	
      a2[idx[1]] = 0.0; 
      grf[1]->Fit("f2", "R", "0",0.,250.);
      NDF = (double)grf[1]->GetFunction("f2")->GetNDF();
      chi2 = grf[1]->GetFunction("f2")->GetChisquare();
      Chi2_S[idx[1]] = chi2/NDF; 
      //f2->SetLineColor(7);
      //f2->Draw("same");
      att[1] = f2->GetParameter(1);
      a2[idx[1]] = f2->GetParameter(0);
      SingleAtten[idx[1]] = att[1];
      ADC100S[idx[1]] = f2->GetParameter(0)*exp(-100./att[1]);
    
      grf[0]->Fit("f12", "R", "0",0.,250. );
      grf[0]->GetFunction("f12")->SetLineColor(4);
      grf[1]->Fit("f12", "R", "0",0.,250. );
      grf[1]->GetFunction("f12")->SetLineColor(6);

      TText *t0 ;
      t0 = new TText(0.4,0.85,"Fit: a1(exp(-x/d1)+exp(-x/d2))");
      
      char lin[128];
      sprintf(lin,"a1:   %6.1f | %6.1f ",a1[idx[0]],a1[idx[1]]);
      TText *t2 = new TText(0.5,0.75,lin);
      sprintf(lin,"d1:   %6.1f | %6.1f ",d1[idx[0]],d1[idx[1]]);
      TText *t3 = new TText(0.5,0.72,lin);
      sprintf(lin,"a2:   %6.1f | %6.1f ",a2[idx[0]],a2[idx[1]]);
      TText *t4 = new TText(0.5,0.63,lin);
      sprintf(lin,"d2:   %6.1f | %6.1f ",d2[idx[0]],d2[idx[1]]);
      TText *t5 = new TText(0.5,0.69,lin);
      
      sprintf(lin,"Ratio blue/pink:   %.2f",Ratio);
      TText *t6 = new TText(0.6, 0.5,lin);

      sprintf(lin,"Single Atten.:   %6.1f | %6.1f ",att[0],att[1]);
      TText *t7 = new TText(0.5, 0.6,lin);
      
      t0->SetNDC();
      t0->SetTextSize(0.038);
      t0->Draw();
      t2->SetNDC();
      t2->SetTextSize(0.03);
      t2->Draw();
      t3->SetNDC();
      t3->SetTextSize(0.03);
      t3->Draw();

      t4->SetNDC();
      t4->SetTextSize(0.03);
      t4->Draw();
    
      t5->SetNDC();
      t5->SetTextSize(0.03);
      t5->Draw();
      t6->SetNDC();
      t6->SetTextSize(0.03);
      t6->Draw();
      t7->SetNDC();
      t7->SetTextSize(0.03);
      t7->Draw();
      
      gPad->Update();
      c2->Update();
      sprintf(grftit,"plots/paddle%d_plane%d_atten_run%d.pdf",Paddle+1, Plane, RunNumber);
      c2->SaveAs(grftit);
      if (DEBUG == 99){
	//getchar();
      }

      if (DEBUG=101) {
	TMinuit *gMinuit = new TMinuit(3);  //initialize TMinuit with a maximum of 5 params
	gMinuit->SetFCN(fcn);
	
	Double_t arglist[10];
	Int_t ierflg = 0;
	
	arglist[0] = 1;
	gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
	
	// Set starting values and step sizes for parameters
	static Double_t vstart[4] = {8000, 400. , 400. , 0.01};
	static Double_t step[4] = {0.1 , 0.1 , 0.01 , 0.001};
	gMinuit->mnparm(0, "amp", vstart[0], step[0], 0, 0, ierflg);
	gMinuit->mnparm(1, "d1", vstart[1], step[1],  0, 0, ierflg);
	gMinuit->mnparm(2, "d2", vstart[2], step[2],  0, 0, ierflg);

	
	// Now ready for minimization step
	arglist[0] = 500;
	arglist[1] = 1.;
	gMinuit->mnexcm("MIGRAD", arglist , 2, ierflg);
	
	// Print results
	Double_t amin,edm,errdef;
	Int_t nvpar,nparx,icstat;
	gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
	gMinuit->mnprin(3,amin);

	for (int s=0; s<NResults[1]; s++){
	  YPos[1][s] *= Ratio;
	}
	grf[2] = new TGraphErrors(NResults[1],XPos[1],YPos[1],dXPos[1],dYPos[1]);
	
	grf[2]->SetMarkerColor(6);
	grf[2]->SetMarkerStyle(21);
      
	char grftit[128];
	sprintf(grftit,"Paddle %d, Plane %d",Paddle+1, Plane);
	TMultiGraph *mgr1 = new TMultiGraph("mgr1",grftit);
	mgr1->Add(grf[0]);
	mgr1->Add(grf[2]);
	
	c2->Clear();
	c2->SetGrid();
	mgr1->Draw("AP");
	gPad->SetGrid();
	mgr1->GetXaxis()->SetTitle("Distance from PMT [cm]");
	gPad->Update();
	
	double TheAmp, dTheAmp;
	gMinuit->GetParameter(0, TheAmp, dTheAmp);
	double Atten1, dAtten1;
	gMinuit->GetParameter(1, Atten1, dAtten1);
	double Atten2, dAtten2;
	gMinuit->GetParameter(2, Atten2, dAtten2);

	TF1 *TheFit1 = new TF1("TheFit1", "[0]*TMath::Exp(-x/[1])", 0., 126.); 
	TF1 *TheFit2 = new TF1("TheFit2", "[0]*TMath::Exp(-(x-126)/[2])*TMath::Exp(-126/[1])", 126., 252.); 

	TheFit1->SetParameter(0, TheAmp);
	TheFit2->SetParameter(0, TheAmp);
	TheFit1->SetParameter(1, Atten1);
	TheFit2->SetParameter(1, Atten1);
	TheFit2->SetParameter(2, Atten2);
	TheFit1->SetLineColor(1);
	TheFit2->SetLineColor(2);

	TheFit1->Draw("same");
	TheFit2->Draw("same");
	gPad->Update();

	int idxx = Plane*88 + Paddle;
	AttenuationData[idxx] = Atten1;
	AttenuationData[idxx+44] = Atten2;
	AllAmplitudes[Plane*44+Paddle] = TheAmp;

	sprintf(lin,"Atten1 = %6.1f", Atten1);
	TText *l1 = new TText(0.6, 0.8,lin);
	sprintf(lin,"Atten2 = %6.1f", Atten2);
	TText *l2 = new TText(0.6, 0.75,lin);
	sprintf(lin,"Amp    = %6.1f", TheAmp);
	TText *l3 = new TText(0.6, 0.7,lin);
	l1->SetNDC();
	l1->SetTextSize(0.03);
	l1->Draw();
	l2->SetNDC();
	l2->SetTextSize(0.03);
	l2->Draw();
	l3->SetNDC();
	l3->SetTextSize(0.03);
	l3->Draw();
	gPad->Update();
	sprintf(grftit,"plots/paddle%d_plane%d_atten_run%d_MINUIT.pdf",Paddle+1, Plane, RunNumber);
	c2->SaveAs(grftit);
	
	//getchar();
      }
    }
  }
  f->Close();

  TH1F *ParA1 = new TH1F("ParA1"," Parameter a1 in attenuation fit", 20., 4000.,10000.);
  TH1F *ParA2 = new TH1F("ParA2"," MPV at x=100cm from Paddle End", 178, -0.5, 177.5);
  TH1F *ParD1 = new TH1F("ParD1"," Parameter d1 in attenuation fit", 20., 0.,200.);
  TH1F *ParD2 = new TH1F("ParD2"," Parameter d2 in attenuation fit", 20., 0.,1000.);
  ParA2->GetXaxis()->SetTitle("PMT #");
  ParA2->GetYaxis()->SetTitle("MPV in ADC units");

  TH1F *ParSingleD = new TH1F("ParSingleD"," Single attenuation parameter fit", 30., 0.,1000.);
  TH1F *ParA3 = new TH1F("ParA3"," MPV at x=100cm from Paddle End Single Att.", 178, -0.5, 177.5);
  TH1F *ParSingleAmp = new TH1F("ParSingleAmp","Amplitude Single atten. fit", 30., 4000.,13000.);

  TH1F *c2_2 = new TH1F("c2_2"," Chi2 double atten. fit", 30, 0.,2.);
  TH1F *c2_1 = new TH1F("c2_1"," Chi2 single atten. fit", 30, 0.,2.);

  char of[128];

  sprintf(of,"calibration/attenuation_values_full_run%d.dat",RunNumber);
  ofstream OUTF2;
  OUTF2.open(of);

  sprintf(of,"calibration/attenuation_values_single_run%d.dat",RunNumber);
  ofstream OUTF3;
  OUTF3.open(of);

  for (int k=0;k<176;k++){
    double y =0;
    if (SingleAtten[k]>0.) {
      ParA1->Fill(a1[k]);
      y = a1[k]*(exp(-100./d1[k]) + exp(-100./d2[k]));
      ParA2->Fill((float)k,y);
      ParD1->Fill(d1[k]);
      ParD2->Fill(d2[k]);
      ParSingleD->Fill(SingleAtten[k]);
      ParA3->Fill((float)k,ADC100S[k]);
 
      ParSingleAmp->Fill(a2[k]);
      c2_2->Fill(Chi2_D[k]);
      c2_1->Fill(Chi2_S[k]);

   }

    OUTF2<<k<<"   "<<a1[k]<<"   "<<d1[k]<<"   "<<d2[k]<<"  "<<y<<endl; 
    OUTF3<<k<<"   "<<SingleAtten[k]<<"   "<<ADC100S[k]<<endl; 

  }

  OUTF2.close();
  OUTF3.close();

  TCanvas *c3 = new TCanvas("c3","Attenuation length Parameters",800.,900.);

  c3->Clear();
  c3->Divide(2,3);
  c3->cd(1);
  gPad->SetGrid();
  ParA1->Draw();
  c3->cd(2);
  gPad->SetGrid();
  ParA2->Draw();
  c3->cd(3);
  gPad->SetGrid();
  ParD1->Draw();
  c3->cd(4);
  gPad->SetGrid();
  ParD2->Draw();
  c3->cd(5);
  gPad->SetGrid();
  ParSingleD->Draw();
  c3->cd(6);
  gPad->SetGrid();
  ParA3->Draw();
  gPad->Update();
  sprintf(fnam,"calibration/attenuation_run%d.pdf", RunNumber);
  c3->SaveAs(fnam);
  sprintf(fnam,"calibration/attenuation_run%d.gif", RunNumber);
  c3->SaveAs(fnam);

  sprintf(fnam,"calibration/attenuation_run%d.root", RunNumber);
  TFile *fnew = new TFile(fnam,"RECREATE");
  fnew->cd();
  ParA1->Write();
  ParA2->Write();
  ParD1->Write();
  ParD2->Write();
  ParSingleD->Write();
  ParA3->Write();
  ParSingleAmp->Write();
  c2_2->Write();
  c2_1->Write();
  
  TH1F *a1MIN = new TH1F("a1MIN", "Attenuation length Atten1 from Minuit", 20, 100., 300.);
  a1MIN->GetXaxis()->SetTitle("attenuation legnth left side [cm]");
  TH1F *a2MIN = new TH1F("a2MIN", "Attenuation length Atten2 from Minuit", 20, 100., 300.);
  a2MIN->GetXaxis()->SetTitle("attenuation legnth right side [cm]");
  TH2F *a1vsa2 = new TH2F("a1vsa2","Minuit fit a1 vs a2", 20, 140., 340.,20, 140., 340.); 
  a1vsa2->GetXaxis()->SetTitle("attenuation legnth right side [cm]");
  a1vsa2->GetYaxis()->SetTitle("attenuation legnth left side [cm]");

  for (int Plane=0; Plane<2; Plane++){
    for (int Paddle=0; Paddle<44; Paddle++){
      int k1 = Plane*88 + Paddle;
      int k2 = k1+44;
      a1MIN->Fill(AttenuationData[k1]);
      a2MIN->Fill(AttenuationData[k2]);
      a1vsa2->Fill(AttenuationData[k1], AttenuationData[k2]);
    }
  }
  TCanvas *c4 = new TCanvas("c4","Attenuation MINUIT",800.,900.);
  c4->Clear();
  c4->Divide(2,2);
  c4->cd(1);
  a1MIN->Draw();
  gPad->SetGrid();
  c4->cd(2);
  a1vsa2->Draw("box");
  gPad->SetGrid();
  c4->cd(3);
  a2MIN->Draw();
  gPad->SetGrid();
  gPad->Update();

  sprintf(fnam,"calibration/attenuation_run%d_MINUIT.pdf", RunNumber);
  c4->SaveAs(fnam);
  sprintf(fnam,"calibration/attenuation_run%d_MINUIT.gif", RunNumber);
  c4->SaveAs(fnam);
  
  a1MIN->Write();
  a2MIN->Write();
  a1vsa2->Write();


  fnew->Close();


  ofstream OUTF1;
  OUTF1.open("calibration/minuit_results_attenuation.dat");
  for (int k=0;k<88;k++){

    OUTF1 << k << "   " << AllRatios[k] << "   " << AllAmplitudes[k]<<endl;

  }  
  OUTF1.close();
  
}