/*****************************************************************************************
This program plots resolution data and the fits for all the threshold values.
It computes the errors for all the data points. You can turn errorbars off for
all the values of threshold except 90% threshold value, if you need.
Input: Needs files Reso_0.XX.data where XX is 70,85,90,95,98 as input
Output: Resolution for all values of threshold and the root file.

****************************************************************************************/



main(TString path,Int_t ComputeErrors = 0){
        #define Err_Quadrature 1
        #define MulitplyBy 1
	#define Constrain 1

        TCanvas *c = new TCanvas("reso", "",10,10,1000,750);
	leg1=new TLegend(0.45,0.65,.91,0.91);
     c->SetFillColor(10);
     c->SetFillStyle(4000);        
 //    c->SetBorderSize(2);
     c->SetGridx();
     c->SetGridy();
     //c->SetFrameFillColor(10);
     //c->SetFrameFillStyle(0);   

        TFile *f=new TFile("../RootFiles/Reso_Thresh_ALLPLOTS.root","RECREATE");   
        TString GIF_File("../GIFs/Reso_Thresh_ALLPLOTS.gif");
        gStyle->SetOptFit();
//--------------------------------------------------------------------------------------------------
        Int_t n=9,ngr,i,j,dof;
	Static Int_t NGrafs = 5;
        Double_t y[NGrafs][n],ey[NGrafs][n], x[NGrafs][n],deltaMean[NGrafs][n],deltaSigma[NGrafs][n];
	ifstream in;
        Double_t Mean1[NGrafs][n],deltaMean1[NGrafs][n],Sigma1[NGrafs][n],deltaSigma1[NGrafs][n],ChiSq1[NGrafs][n],Dof1[NGrafs][n];
        Double_t Mean2[NGrafs][n],deltaMean2[NGrafs][n],Sigma2[NGrafs][n],deltaSigma2[NGrafs][n],ChiSq2[NGrafs][n],Dof2[NGrafs][n];
        Double_t ex[NGrafs][n],E[NGrafs][n],Mean[NGrafs][n],Sigma[NGrafs][n];
	Double_t intercept,chi,slope;
        Char_t q,str[NGrafs][120];
        Double_t temp1, temp2, dm1, dm2, ds1, ds2;
	Double_t ErrZero[] = {0,0,0,0,0,0,0,0,0};
	TString ip_file[] = {"Reso_0.70.data","Reso_0.85.data","Reso_0.90.data","Reso_0.95.data","Reso_0.98.data" };
 //	TString str[NGrafs];
      for(ngr= 0; ngr < NGrafs; ngr++) {
         in.open(path + ip_file[ngr],ios::in);
         for(i = 0; i < n; i++)
          in >> E[ngr][i] >>  Mean1[ngr][i] >> q >> deltaMean1[ngr][i] >> q >> ChiSq1[ngr][i] >> q >> Dof1[ngr][i] >> q >> 
		Sigma1[ngr][i] >> q >> deltaSigma1[ngr][i];

         for(i = 0; i < n; i++)
          in >> E[ngr][i] >>  Mean2[ngr][i] >> q >> deltaMean2[ngr][i] >> q >> ChiSq2[ngr][i] >> q >> Dof2[ngr][i] >> q >> 
		Sigma2[ngr][i] >> q >> deltaSigma2[ngr][i];
         in.close();
	} // for ngr

       for(ngr= 0; ngr < NGrafs; ngr++) {
         for(i = 0; i < n ; i++) {

             Sigma[ngr][i]=Sigma2[ngr][i];
             Mean[ngr][i]=Mean2[ngr][i];
             dm1 = Mean1[ngr][i] - Mean2[ngr][i];
             dm2 = deltaMean2[ngr][i];
             dm2 = dm2 * TMath::Sqrt(ChiSq2[ngr][i]/Dof2[ngr][i]);
             ds1 = Sigma1[ngr][i] - Sigma2[ngr][i];
             ds2 = deltaSigma2[ngr][i];
             ds2 = ds2 * TMath::Sqrt(ChiSq2[ngr][i]/Dof2[ngr][i]);

             deltaMean[ngr][i] = TMath::Sqrt(dm1 * dm1  + dm2 * dm2);
             deltaSigma[ngr][i] = TMath::Sqrt(ds1 * ds1  + ds2 * ds2);

             x[ngr][i] = 1.0 / TMath::Sqrt(Mean[ngr][i]);         // x - axis = 1/Mean

             ex[ngr][i] = deltaMean[ngr][i] / 2 * Mean[ngr][i] * TMath::Sqrt(Mean[ngr][i]); // x axis err = dE / 2E(E)^1/2
	     y[ngr][i] = Sigma[ngr][i] * 100.0 / Mean[ngr][i] ;		// y - axis = Sigma / Mean (%)
	  
	     temp1 = deltaSigma[ngr][i] / Sigma[ngr][i];   
	     temp2 = deltaMean[ngr][i] / Mean[ngr][i];   
	     ey[ngr][i] = y[ngr][i]*TMath::Sqrt(temp1 * temp1 + temp2 * temp2); // y axis err = S/E [ (dS/S)^2 + (dE/E)^2 ]	
//	     ey[ngr][i] = 0.239* ey[ngr][i];

         }
	} // for ngr 

       TGraph *gr[NGrafs];
       for(ngr= 0; ngr < NGrafs; ngr++) {
	   if(ComputeErrors == 1)
             gr[ngr]=new TGraphErrors(n,Mean[ngr],y[ngr],deltaMean[ngr],ey[ngr]);
	   else {
	   if(ngr == 2) 
             gr[ngr]=new TGraphErrors(n,Mean[ngr],y[ngr],deltaMean[ngr],ey[ngr]);
	   else
             gr[ngr]=new TGraphErrors(n,Mean[ngr],y[ngr],ErrZero,ErrZero);
	   }
           gr[ngr]->SetLineColor(ngr+1);
           gr[ngr]->SetMarkerColor(ngr+1);
           gr[ngr]->SetMarkerStyle(20+ngr);
           gr[ngr]->SetMarkerSize(1.0);
	} // for ngr

      TF1 *f1[NGrafs];
      for(ngr= 0; ngr < NGrafs; ngr++) {
           f1[ngr] = new TF1("f1","([1]*sqrt(1/x)+[0])",5,100);
           f1[ngr]->SetMarkerStyle(20+ngr);
           f1[ngr]->SetMarkerColor(ngr+1);
           f1[ngr]->SetMarkerSize(1.0);
           f1[ngr]->SetLineColor(ngr+1);
           gr[ngr]->Fit("f1");
         intercept=f1[ngr]->GetParameter(0);
         slope=f1[ngr]->GetParameter(1);   
         chi=f1[ngr]->GetChisquare();
         dof=f1[ngr]->GetNDF();
         //fline->SetLineColor(ngr+1);
         sprintf(str[ngr],"#frac{#sigma}{#LTE#GT} = #frac{%5.2f%s}{#sqrt{ #LTE#GT }} + %5.2f%s (GEM- %s) ",slope,"%",intercept,"%",(const char *)ip_file[ngr]);
         leg1->AddEntry(f1[ngr],str[ngr],"L");
       } // for ngr

   TMultiGraph *mg = new TMultiGraph();
   mg->Add(gr[0]);
   mg->Add(gr[1]);
   mg->Add(gr[2]);
   mg->Add(gr[3]);
   mg->Add(gr[4]);
//   mg->Add(bmark);
   mg->Draw("AP");
   leg1->Draw();
   mg->GetXaxis()->SetTitle("#LTE#GT (GeV)");
   mg->GetXaxis()->CenterTitle(kFALSE);
   mg->GetXaxis()->SetTitleSize(0.03);
   mg->GetXaxis()->SetLabelSize(0.02);
   mg->GetXaxis()->SetNdivisions(520);
   mg->GetYaxis()->SetTitle("#frac{#sigma}{#LTE#GT} (%)");
   mg->GetYaxis()->CenterTitle(kFALSE);
   mg->GetYaxis()->SetLabelSize(0.02);
   mg->GetYaxis()->SetTitleSize(0.03);


//--------------------------------------------------------------------------------------------------------------------------------



 	c->cd();
        c->Print(*GIF_file);
        f->Write();

}