#include <cstdlib>
#include <vector>
#include <iostream>
#include <map>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TStopwatch.h"
//
#include "TROOT.h"
#include "TRint.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TVector3.h"
#include "TMath.h"
#include "TFile.h"
#include "TSystem.h"
#include "TGStatusBar.h"
#include "TSystem.h"
#include "TXMLEngine.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TList.h"
#include "TMatrixT.h"
#include "TMatrixD.h"
#include "TLatex.h"
#include "TLorentzRotation.h"

#include <cmath>
//#include <iostream>
#include <string>
#include <cstdlib>
#include<iomanip>
#include <fstream>
#include <map>
#include <sstream>
#include <stdexcept>
#include <vector>


using namespace std;
#define SQR(a) ((a)*(a))



int main()
{   

//
int Per_val = 9;

  int sit_limit = 9;

  int BPol_val = 9;

  float CL_cut = 0.0;

  float ebin_width = -9.0;
  // in MeV                                                                                                    
  float lowest_energy = -9.0; 
  //in MeV                                                                                                                 

  int NewEBin_val = 9;
  int NewEBin = -1;
  int NewEBin_car = -1;

  int NewPBin_val = 9;
  int NewPBin = 9;
  int NewPBin_car = 9; 

  

//
char name[500],name2[500],name3[500],name11[500],name12[500];
Float_t q_err_nt;

int NNevent_nt[1000];
float q_err_i, q_err_j;
int count_common;


int prefix;
  Int_t event, run;
  
  Float_t kin_beamE, beamE;
 Float_t sum_q[8]={0};
  
  

  Float_t z;
  
  Float_t pip_px;
  Float_t pip_py;
  Float_t pip_pz;
  Float_t pim_px;
  Float_t pim_py;
  Float_t pim_pz;
  Float_t kin_pip_px;
  Float_t kin_pip_py;
  Float_t kin_pip_pz;
  Float_t kin_pim_px;
  Float_t kin_pim_py;
  Float_t kin_pim_pz;
  
  Float_t p_px;
  Float_t p_py;
  Float_t p_pz;
  Float_t kin_p_px;
  Float_t kin_p_py;
  Float_t kin_p_pz;
  Float_t pip_E;
  Float_t pim_E;
  Float_t p_E;
  Float_t kin_pip_E;
  Float_t kin_pim_E;
  Float_t kin_p_E;
  

 
  Float_t bec_pi0_prob;
  Float_t bec_nop_prob;
Float_t Qvalue, Qerror;

int beampol;
Float_t RealPhi,tetastar,cosprocm;
RealPhi = 0;
tetastar = 0;

for (int iii=Per_val; iii<=Per_val; iii++) //period
for (int jjj=NewEBin_val; jjj<=NewEBin_val; jjj++) //binenergy
for (int kkk=BPol_val; kkk<=BPol_val; kkk++) //beampol
{
{
{

Float_t correl_error_sqr[20]={0};


//TChain *ttree;

//ttree = new TChain("finaltree");

//input data
sprintf(name2,"/d/grid13/akbar/crosskaon/qvalue_100MeV/ROOTFILES/output_Sit05_Per%02d_NewEBin%03d_BPol%02d_NewPBin01.root",iii,jjj,kkk);
//(name,"outQvalue_Sit%02d_Per%02d_NewEBin%02d_BPol%02d_NewPBin%02d.txt",5,Per_val,NewEBin_val,BPol_val,NewPBin_val);
//ttree->Add(name2);
cerr<<"get the file from :"<<name2<<endl;

TFile *f = new TFile(name2); 


    
   

TTree *tree = (TTree*)gROOT->FindObject("finaltree");


tree->SetBranchAddress( "NNarray_nt", &NNevent_nt);
tree->SetBranchAddress( "Qfit_error_nt", &q_err_nt);
 tree->SetBranchAddress("run_number_nt",&run);
  tree->SetBranchAddress("eventnum_nt",&event);
  tree->SetBranchAddress("z_vertex_nt",&z);
 
  tree->SetBranchAddress("pip_pz_kinf_nt",&kin_pip_pz);
  tree->SetBranchAddress("pip_px_kinf_nt",&kin_pip_px);
  tree->SetBranchAddress("pip_py_kinf_nt",&kin_pip_py);
  tree->SetBranchAddress("pim_pz_kinf_nt",&kin_pim_pz);
  tree->SetBranchAddress("pim_px_kinf_nt",&kin_pim_px);
  tree->SetBranchAddress("pim_py_kinf_nt",&kin_pim_py);
  tree->SetBranchAddress("pro_pz_kinf_nt",&kin_p_pz);
  tree->SetBranchAddress("pro_px_kinf_nt",&kin_p_px);
  tree->SetBranchAddress("pro_py_kinf_nt",&kin_p_py);
  tree->SetBranchAddress("pip_E_kinf_nt",&kin_pip_E);
  tree->SetBranchAddress("pim_E_kinf_nt",&kin_pim_E);
  tree->SetBranchAddress("pro_E_kinf_nt",&kin_p_E);
  tree->SetBranchAddress("pho_energy_kinf_nt",&kin_beamE);
  
  tree->SetBranchAddress("pip_pz_nt",&pip_pz);
  tree->SetBranchAddress("pip_px_nt",&pip_px);
  tree->SetBranchAddress("pip_py_nt",&pip_py);
  tree->SetBranchAddress("pim_pz_nt",&pim_pz);
  tree->SetBranchAddress("pim_px_nt",&pim_px);
  tree->SetBranchAddress("pim_py_nt",&pim_py);
  tree->SetBranchAddress("pro_pz_nt",&p_pz);
  tree->SetBranchAddress("pro_px_nt",&p_px);
  tree->SetBranchAddress("pro_py_nt",&p_py);
  tree->SetBranchAddress("pip_E_nt",&pip_E);
  tree->SetBranchAddress("pim_E_nt",&pim_E);
  tree->SetBranchAddress("pro_E-nt",&p_E);
  tree->SetBranchAddress("pho_energy_nt",&beamE);
    
    tree->SetBranchAddress("cl_3pi_nt",&bec_pi0_prob);
    tree->SetBranchAddress("cl_2pi_nt",&bec_nop_prob);

tree->SetBranchAddress("realphi_nt",&RealPhi);
    tree->SetBranchAddress("cosrealtheta_nt",&tetastar);
tree->SetBranchAddress("coscmstheta_pro_nt",&cosprocm);

tree->SetBranchAddress("Qval_nt",&Qvalue);
tree->SetBranchAddress("bpol_nt",&beampol);

//
Int_t nentries = (Int_t)tree->GetEntries();






Float_t counting;
counting = 0;

sprintf(name3,"/d/grid13/akbar/crosskaon/qerror/sigma100MeV/output_Sit05_Per%02d_NewEBin%03d_BPol%02d_NewPBin01.root",iii,jjj,kkk);
TFile hh (name3, "RECREATE", "Cross Sections from histos");

TH1D *histo;
TH1D *histo_data;
 
 float ener = lowest_energy/1000.0 + (ebin_width/1000.0 *(NewEBin_val-1));
  sprintf(name11,"histo_energy_%f",ener);  
  sprintf(name12,"cross_energy_%f",ener); 
  histo = new TH1D(name11,name11, 20,-1,1);
  histo_data = new TH1D(name12,name12, 20,-1,1);
  histo->Sumw2();
  histo_data->Sumw2();



                                         for (long i=0; i<nentries; i++)
					    {
                                              
					      //cerr<<i+1<<" "<<nentries<<" "<<endl;
                                               
                                              if(i%20 == 0) cerr << " QVALUEINFO: eventnumber = " << i << " OF  "<<nentries<<endl;
                                                tree->GetEntry(i);

					
					           q_err_i = q_err_nt;
                                                   double Wbin_i = TMath::Sqrt((2*kin_beamE*938.272)+(938.27*938.27));
                                                   //int binenergy_i = int ((Wbin_i-(lowest_energy))/ebin_width) ;
						   int binenergy_i = int ((kin_beamE-(lowest_energy))/ebin_width) ;
                                                   int binteta_i= int ((-1*cosprocm+1)/0.1);
                                                 float Qvalue_i = Qvalue;
                                                   histo_data->Fill(-1*cosprocm,Qvalue_i);
				                   int NN_event_i[1000];
                                                  // sum_q[binteta_i]=sum_q[binteta_i]+Qvalue_i*Qvalue_i;
				                   for(int count=0; count<1000; count++)
						     {
					               NN_event_i[count] = NNevent_nt[count];
					             }
						   
                                                for(long j=0; j<nentries; j++)
						     {
						       
					                   tree->GetEntry(j);							   
							       q_err_j = q_err_nt;
                                                               double Wbin_j = TMath::Sqrt((2*kin_beamE*938.272)+(938.27*938.27));
                                                               //int binenergy_j = int ((Wbin_j-(lowest_energy))/ebin_width) ;
							       int binenergy_j = int ((kin_beamE-(lowest_energy))/ebin_width) ;
                                                               int binteta_j= int ((-1*cosprocm+1)/0.1);
                                                         if ((binteta_i == binteta_j) /*&& (binenergy_i == binenergy_j) && binenergy_i>=0 */)
                                                           {
					                        int NN_event_j[1000];
					                        for(int count=0; count<1000; count++)
								  {
						                     NN_event_j[count] = NNevent_nt[count];
						        	  }
					                        count_common = 0;
					                        for(int a=0; a<1000; a++)
								  {
						                    for(int b=0; b<1000; b++)
								      {
						                        if (NN_event_i[a] == NN_event_j[b])
									  {
							                    count_common++;
									  }
								      }
								  }
                                                                
      					                         correl_error_sqr[binteta_i] = correl_error_sqr[binteta_i] + q_err_i*q_err_j*(count_common/1000.0);          
                                                            }
                                                                                                      
							     //cerr<<i<<" "<<j<<" "<<count_common<<" "<<Qvalue<<" "<<correl_error_sqr<<endl;
						    }

				     }
					





for (int l=1; l<=20; l++)
 {
  double a = histo_data->GetBinContent(l);
  histo->SetBinContent(l,a);
  double b = TMath::Sqrt(correl_error_sqr[l-1]);
  histo->SetBinError(l,b);
  
  
 }
 


histo->Write();
histo_data->Write();




cerr<<"FINISH : "<<nentries<<endl;
}
}
}



}