#include <iostream>
#include <fstream>
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"


using namespace std;

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]);
}

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

   //Negative version of langau
   //Effectively, for given x, we want to return f(-x). Then, the parameters will basically remain the same.
   //Note that this is simply an array object 

   //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;
      Double_t xneg;

      xneg=-x[0];

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

      // Range of convolution integral
      // Basically, the convolution takes place over +- 5 Gaussian sigmas
      xlow = xneg - sc * par[3];
      xupp = xneg + sc * par[3];

      // Then, convolution step size, based on range
      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      // Find integr(g(t-Tau)h(Tau)dTau

      for(i=1.0; i<=np/2; i++) {
	xx = xlow + (i-.5) * step;    // Translate i into a "Tau" value (xx)
	fland = TMath::Landau(xx,mpc,par[0]) / par[0];  // Get h(Tau) - i.e. value of Landau at xx 
	sum += fland * TMath::Gaus(xneg,xx,par[3]); // Multiply this by g(t-Tau), by getting Gaussian value at xneg when the dist is centred at xx

        xx = xupp - (i-.5) * step; // When stepping through the integral, we step from both the lower and upper limits, and meeting in the middle (symmetry of Gaussian function means this is more economical)
        fland = TMath::Landau(xx,mpc,par[0]) / par[0];
        sum += fland * TMath::Gaus(xneg,xx,par[3]);
      }

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



TF1 *langaufit(TH1F *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);
}

void L_cer_sp_langaus(TString runlist = "run_list.txt",TString path = "/fs/work1/g2p/jie/Rootfiles/",Int_t mom = 1792){

  //open list of runs to use

  ifstream file_list(runlist.Data(),ios::in);

  //check to make sure the file opened
    if(!file_list){
    cout << "Whoops, can't find the file list!" << endl;
    cout << "Exiting the program." << endl;
    exit(1);
  }

    Int_t runnum;
    Int_t nsplit;


    TChain *T = new TChain("T");

    file_list >> runnum;

      while(!file_list.eof()){

      file_list >> nsplit;

    cout << "runnum = " << runnum << endl;
     cout << "nsplit = " << nsplit << endl;
 

         if(runnum == 0){
 cout << "Run Number 0 does not exist. Don't forget to specify a run number!" <<endl;
       } 
	 /*
        else{ 
	  T->Add(Form("/g2p_L_%d.root",runnum));
	  cout << "file 0 for run  " << runnum << endl;
	  T->Add(Form("/w/halla-1/g2p/melissa/replay/Rootfiles/g2p_L_%d_1.root",runnum));
	  cout << "file 1 for run  " << runnum << endl;
	  T->Add(Form("/w/halla-1/g2p/melissa/replay/Rootfiles/g2p_L_%d_2.root",runnum));
	  cout << "file 2 for run  " << runnum << endl;  
	}

	 */

          else{ 
	 cout << "test = "<< nsplit << endl;
	 for(Int_t i=0;i<=nsplit;i++){
	cout << "i=" << i << endl;
	if(i==0){
	  TString input_filename = path + "g2p_L_";
	  input_filename += runnum;
	  input_filename += ".root";
	  cout << "Input File = " << input_filename << endl;
	  T->Add(input_filename);
	}
	else{
          TString input_filename = path + "g2p_L_";
	  input_filename += runnum;
	  input_filename += "_";
	  input_filename += i;
	  input_filename += ".root";
	  cout << "Input File = " << input_filename << endl;
	  T->Add(input_filename);
	}

      } // end of for loop  
	

	  }
	/*     TString input_filename = path + "g2p_L_";
      input_filename += runnum;
      input_filename += ".root";
      cout << "Input File = " << input_filename << endl;
      T->Add(input_filename); */

	 //            } // end of else loop


       file_list >> runnum;
     

      } //end of while loop
      file_list.close();

 
      Int_t num = T->GetEntries();
      cout << "number of entries = " << num << endl;


   
      TString s1,s2,s3;

      TCanvas *c1 = new TCanvas("c1","",1200,1000);
      c1->Divide(2,5);




 /* for(int i=0;i<10;i++){

   s1 = Form("L.cer.t[%d]>>h%d(200,0,3500)",i,i);
   s2 = Form("L.cer.t[%d]>0",i);
   s3 = Form("LHRS Cherenkov TDC %d",i); // need to figure out a way to set each histogram title individually...

   c1->cd(i+1);
   T->Draw(s1,s2);

   gPad->SetLogy();



   } */



 //define fit
 /*     
p1 = new TF1("m1", "gaus",65,180);
p2 = new TF1("m2", "gaus",65,210);
p3 = new TF1("m3", "gaus",60,140);
p4 = new TF1("m4", "gaus",65,215);
p5 = new TF1("m5", "gaus",65,160);
p6 = new TF1("m6", "gaus",65,190);
 */
//p1 = new TF1("m1", "gaus",70,200);

 Double_t mean_sp[10] = {0};
 Double_t mp_sp[10] = {0};
 Double_t fit_err[10] = {0};
 Double_t mp_width[10] = {0};
 Double_t mirror[10] = {1,2,3,4,5,6,7,8,9,10};

 //*******************************************
 // Setting fit range and start values
 Double_t fr[2],fr1[2];
   Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
   // fr[0]=0.3*hSNR->GetMean();
   //fr[1]=3.0*hSNR->GetMean();
   fr[0]=60;
   fr[1]=400;
   fr1[0]=40;
   fr1[1]=800;
   pllo[0]=0.1; pllo[1]=100.0; pllo[2]=1.0; pllo[3]=0.0;
   plhi[0]=100.0; plhi[1]=170.0; plhi[2]=100000000.0; plhi[3]=50.0;
   sv[0]=10; sv[1]=130.0; sv[2]=100000.0; sv[3]=1.0;

   Double_t chisqr;
   Int_t    ndf;
   Double_t SNRPeak, SNRFWHM;
   //********************************************
 //Making TDC Cuts

 gStyle->SetOptFit(1111);
 gStyle->SetOptStat(1111);
 gStyle->SetLabelSize(0.03,"x");
 gStyle->SetLabelSize(0.03,"y");

 c1->cd(1);
 TH1F *h0=new TH1F("h0","L.cer.a_p[0]",190,20,400);
 TCut cut0="(L.cer.t[0]>1630&&L.cer.t[0]<1665)||(L.cer.t[0]>1695&&L.cer.t[0]<1705)&&abs(L.cer.trx[0]+0.7391)>0.3";
 T->Project("h0","L.cer.a_p[0]",cut0);
 TF1 *fitsnr0 = langaufit(h0,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m1","L.cer.a_p[0]>>h0(100,15,1000)","(L.cer.t[0]>1630&&L.cer.t[0]<1665)||(L.cer.t[0]>1695&&L.cer.t[0]<1705)&&abs(L.cer.trx[0]+0.7391)>0.35","R");
 h0->Draw();   fitsnr0->Draw("lsame");

 mean_sp[0] = fp[1];
 fit_err[0] = fp[3];
 mp_sp[0]=SNRPeak;
 mp_width[0]=SNRFWHM;
 h0->SetTitle("LHRS Cherenkov Mirror[0]");

 //&&abs(L.cer.try[0]+0.0445)>0.2

 
 c1->cd(2);
 TH1F *h1=new TH1F("h1","L.cer.a_p[1]",190,20,800);
 TCut cut1="(L.cer.t[1]>1630&&L.cer.t[1]<1660)||(L.cer.t[1]<1710&&L.cer.t[1]>1700)&&abs(L.cer.trx[0]+0.7000)>0.3";
 T->Project("h1","L.cer.a_p[1]",cut1);
 TF1 *fitsnr1 = langaufit(h1,fr1,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m4","L.cer.a_p[1]>>h1(100,15,1000)","(L.cer.t[1]>1630&&L.cer.t[1]<1660)||(L.cer.t[1]<1710&&L.cer.t[1]>1700)&&L.cer.t[1]>100&&abs(L.cer.trx[0]+0.7000)>0.35","R");
 h1->Draw();   fitsnr1->Draw("lsame");
 mean_sp[1] = fp[1];
 fit_err[1] = fp[3];
 mp_sp[1]=SNRPeak;
 mp_width[1]=SNRFWHM;
 h1->SetTitle("LHRS Cherenkov Mirror[1]");
 //&&abs(L.cer.try[0]-0.0256)>0.2


 c1->cd(3);
 TH1F *h2=new TH1F("h2","L.cer.a_p[2]",190,20,400);
 TCut cut2="(L.cer.t[2]>1630&&L.cer.t[2]<1660)||(L.cer.t[2]<1710&&L.cer.t[2]>1700)&&abs(L.cer.trx[0]+0.5153)>0.3";
 T->Project("h2","L.cer.a_p[2]",cut2);
 TF1 *fitsnr2 = langaufit(h2,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m3","L.cer.a_p[2]>>h2(100,15,1000)","(L.cer.t[2]>1630&&L.cer.t[2]<1660)||(L.cer.t[2]<1710&&L.cer.t[2]>1700)&&L.cer.t[2]>100&&abs(L.cer.trx[0]+0.5153)>0.35","R");
 h2->Draw();   fitsnr2->Draw("lsame");
 mean_sp[2] = fp[1];
 fit_err[2] = fp[3];
 mp_sp[2]=SNRPeak;
 mp_width[2]=SNRFWHM;
 // mean_sp[2] = p3->GetParameter(1);
 // fit_err[2] = p3->GetParError(1);
 h2->SetTitle("LHRS Cherenkov Mirror[2]");
 //&&abs(L.cer.try[0]+0.0244)>0.2

 c1->cd(4);
 TH1F *h3=new TH1F("h3","L.cer.a_p[3]",190,20,800);
 TCut cut3="(L.cer.t[3]>1630&&L.cer.t[3]<1660)||(L.cer.t[3]<1700&&L.cer.t[3]>1690)&&abs(L.cer.trx[0]+0.4264)>0.3";
 T->Project("h3","L.cer.a_p[3]",cut3);
 TF1 *fitsnr3 = langaufit(h3,fr1,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 //T->Fit("m4","L.cer.a_p[3]>>h3(100,0,1000)","L.cer.t[3]>1600&&L.cer.t[3]<1730&&abs(L.cer.trx[0]+0.4264)>0.35","R");
 h3->Draw();   fitsnr3->Draw("lsame");
 mean_sp[3] = fp[1];
 fit_err[3] = fp[3];
 mp_sp[3]=SNRPeak;
 mp_width[3]=SNRFWHM;
 // mean_sp[3] = p4->GetParameter(1);
 // fit_err[3] = p4->GetParError(1);
 h3->SetTitle("LHRS Cherenkov Mirror[3]");
 //&&abs(L.cer.try[0]-0.0435)>0.2


 c1->cd(5);
 TH1F *h4=new TH1F("h4","L.cer.a_p[4]",190,20,400);
 TCut cut4="(L.cer.t[4]>1630&&L.cer.t[4]<1660)||(L.cer.t[4]<1710&&L.cer.t[4]>1700)&&abs(L.cer.trx[0]+0.1463)>0.3";
 T->Project("h4","L.cer.a_p[4]",cut4);
 TF1 *fitsnr4 = langaufit(h4,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m1","L.cer.a_p[4]>>h4(120,0,1000)","L.cer.t[4]>1600&&L.cer.t[4]<1750&&abs(L.cer.trx[0]+0.1463)>0.35","R");
 h4->Draw();   fitsnr4->Draw("lsame");
 mean_sp[4] = fp[1];
 fit_err[4] = fp[3];
 mp_sp[4]=SNRPeak;
 mp_width[4]=SNRFWHM;
 // mean_sp[4] = p1->GetParameter(1);
 // fit_err[4] = p1->GetParError(1);
 h4->SetTitle("LHRS Cherenkov Mirror[4]");
 //&&abs(L.cer.try[0]+0.0110)>0.2


 c1->cd(6);
 TH1F *h5=new TH1F("h5","L.cer.a_p[5]",190,20,400);
 TCut cut5="(L.cer.t[5]>1630&&L.cer.t[5]<1660)||(L.cer.t[5]<1710&&L.cer.t[5]>1700)&&abs(L.cer.trx[0]+0.0703)>0.3";
 T->Project("h5","L.cer.a_p[5]",cut5);
 TF1 *fitsnr5 = langaufit(h5,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m5","L.cer.a_p[5]>>h5(120,0,1000)","L.cer.t[5]>1600&&L.cer.t[5]<1750&&abs(L.cer.trx[0]+0.0703)>0.35","R");
 h5->Draw();   fitsnr5->Draw("lsame");
 // mean_sp[5] = p5->GetParameter(1);
 //fit_err[5] = p5->GetParError(1);
 mean_sp[5] = fp[1];
 fit_err[5] = fp[3];
 mp_sp[5]=SNRPeak;
 mp_width[5]=SNRFWHM;
 h5->SetTitle("LHRS Cherenkov Mirror[5]");
 //&&abs(L.cer.try[0]-0.0594)>0.2


 c1->cd(7);
 TH1F *h6=new TH1F("h6","L.cer.a_p[6]",190,20,400);
 TCut cut6="(L.cer.t[6]>1630&&L.cer.t[6]<1660)||(L.cer.t[6]<1710&&L.cer.t[6]>1700)&&abs(L.cer.trx[0]-0.2405)>0.3";
 T->Project("h6","L.cer.a_p[6]",cut6);
 TF1 *fitsnr6 = langaufit(h6,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);

 // T->Fit("m1","L.cer.a_p[6]>>h6(120,0,1000)","L.cer.t[6]>1600&&L.cer.t[6]<1750&&abs(L.cer.trx[0]-0.2405)>0.35","R");
 h6->Draw();   fitsnr6->Draw("lsame");
 mean_sp[6] = fp[1];
 fit_err[6] = fp[3];
 mp_sp[6]=SNRPeak;
 mp_width[6]=SNRFWHM;
 // mean_sp[6] = p1->GetParameter(1);
 // fit_err[6] = p1->GetParError(1);
 h6->SetTitle("LHRS Cherenkov Mirror[6]");
 //&&abs(L.cer.try[0]+0.0039)>0.2


 c1->cd(8);
 TH1F *h7=new TH1F("h7","L.cer.a_p[7]",190,20,400);
 TCut cut7="(L.cer.t[7]>1630&&L.cer.t[7]<1660)||(L.cer.t[7]<1710&&L.cer.t[7]>1700)&&abs(L.cer.trx[0]-0.2937)>0.3";
 T->Project("h7","L.cer.a_p[7]",cut7);
 TF1 *fitsnr7 = langaufit(h7,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 //T->Fit("m1","L.cer.a_p[7]>>h7(120,0,1000)","L.cer.t[7]>1600&&L.cer.t[7]<1750&&abs(L.cer.trx[0]-0.2937)>0.35","R");
 h7->Draw();   fitsnr7->Draw("lsame");
 mean_sp[7] = fp[1];
 fit_err[7] = fp[3];
 mp_sp[7]=SNRPeak;
 mp_width[7]=SNRFWHM;
 //mean_sp[7] = p1->GetParameter(1);
 //fit_err[7] = p1->GetParError(1);
 h7->SetTitle("LHRS Cherenkov Mirror[7]");
 //&&abs(L.cer.try[0]-0.0448)>0.2


 c1->cd(9);
 TH1F *h8=new TH1F("h8","L.cer.a_p[8]",190,20,400);
 TCut cut8="(L.cer.t[8]>1630&&L.cer.t[8]<1660)||(L.cer.t[8]<1710&&L.cer.t[8]>1700)&&abs(L.cer.trx[0]-0.5471)>0.3";
 T->Project("h8","L.cer.a_p[8]",cut8);
 TF1 *fitsnr8 = langaufit(h8,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);
 // T->Fit("m1","L.cer.a_p[8]>>h8(120,0,1000)","L.cer.t[8]>1600&&L.cer.t[8]<1750&&abs(L.cer.trx[0]-0.5471)>0.35","R");
 h8->Draw();    fitsnr8->Draw("lsame");
 mean_sp[8] = fp[1];
 fit_err[8] = fp[3];
 mp_sp[8]=SNRPeak;
 mp_width[8]=SNRFWHM;
 // mean_sp[8] = p1->GetParameter(1);
 // fit_err[8] = p1->GetParError(1);
 h8->SetTitle("LHRS Cherenkov Mirror[8]");
 //&&abs(L.cer.try[0]-0.0001)>0.2


 c1->cd(10);
 TH1F *h9=new TH1F("h9","L.cer.a_p[9]",190,20,400);
 TCut cut9="(L.cer.t[9]>1630&&L.cer.t[9]<1660)||(L.cer.t[9]<1710&&L.cer.t[9]>1700)&&abs(L.cer.trx[0]-0.5722)>0.3";
 T->Project("h9","L.cer.a_p[9]",cut9);
 TF1 *fitsnr9 = langaufit(h9,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
 langaupro(fp,SNRPeak,SNRFWHM);

 // T->Fit("m6","L.cer.a_p[9]>>h9(120,0,1000)","L.cer.t[9]>1600&&L.cer.t[9]<1750&&abs(L.cer.trx[0]-0.5722)>0.35","R");
 h9->Draw();   fitsnr9->Draw("lsame");
 mean_sp[9] = fp[1];
 fit_err[9] = fp[3];
 mp_sp[9]=SNRPeak;
 mp_width[9]=SNRFWHM;
 // mean_sp[9] = p6->GetParameter(1);
 //fit_err[9] = p6->GetParError(1);
 h9->SetTitle("LHRS Cherenkov Mirror[9]");
 //&&abs(L.cer.try[0]-0.0619)>0.2



 TString sp_filename =  "/adaqfs/home/adaq/g2p/jie/analysis/HRScalib/cer/lhrs/myoutput/";
 sp_filename += "sp_origin_";
 sp_filename += mom;
 sp_filename += ".txt";
 

 ofstream outfile(sp_filename.Data(),ios::app);
 cout << sp_filename << endl;

 Double_t calib[10] = {0};

 for(int i=0;i<10;i++){

   calib[i] = 100.0/mean_sp[i];


   //need to add more to output file, error, difference from mean value of histo, mean value??
   outfile << mirror[i] << "\t"  << mean_sp[i]<< "\t" << fit_err[i] << "\t" << calib[i]  << endl; 
   cout << mirror[i] << "\t"  << mean_sp[i]<< "\t" << fit_err[i] << "\t" << calib[i]  << "\t" << mp_sp[i]<<"\t"<<mp_width[i]<<endl; 
 }

 outfile.close();


 /*
 TCanvas *c3 = new TCanvas("c3","Peak Location");
c3->SetGrid();
gr1 = new TGraph(10,mirror,mean_sp);
gr1->SetMarkerStyle(8);
gr1->Draw("AP");


TLine* l2 = new TLine(0,100,11,100);
l2->SetLineStyle(1);
l2->SetLineColor(2);

gr1->GetXaxis()->SetRangeUser(0.00,11.00);
gr1->GetYaxis()->SetRangeUser(0.00,200.00);
gr1->Draw("AP");
gr1->SetTitle("LHRS Cherenkov SPE  Peak Alignment: Channel[#] vs ADC Channel");
l2->Draw();TLatex *tex3 = new TLatex(8,200,"");
tex3->SetTextSize(0.04);
tex3->SetTextColor(2);
tex3->Draw("same");

c3->Update(); 

 */



}