#include <iostream>
#include "fstream.h"
#include <iomanip>
#include "TGraphErrors.h"
#include "TLatex.h"
#include <vector>

int main()
{
  int size = 11;
  vector<double> onpar(11);
  vector<double> sdev(11);
  
  //turning on the parameters
  onpar[0]=1;  // luminosity			      
  onpar[1]=1;  // Background of Z(CC)		      
  onpar[2]=1;  // Background of Z(CP)		      
  onpar[3]=1;  // Background of Z(PP)		      
  onpar[4]=1;  // Central Electron ID Efficiency      
  onpar[5]=1;  // Plug Electron ID Efficiency	      
  onpar[6]=1;  // Central Material Effect	      
  onpar[7]=1;  // Plug Material Effect		      
  onpar[8]=1;  // ZVertex Finding Efficiency in Z(PP) 
  onpar[9]=1;  // Electron Tracking Efficiency	      
  onpar[10]=1;  // No Track Events Fraction           

  //setting what sigma(uncertainty) will be deviated from the central value
  sdev[0]=1;    // luminosity			      
  sdev[1]=1;    // Background of Z(CC)		      
  sdev[2]=1;    // Background of Z(CP)		      
  sdev[3]=1;    // Background of Z(PP)		      
  sdev[4]=1;    // Central Electron ID Efficiency      
  sdev[5]=1;    // Plug Electron ID Efficiency	      
  sdev[6]=1;    // Central Material Effect	      
  sdev[7]=1;    // Plug Material Effect		      
  sdev[8]=1;    // ZVertex Finding Efficiency in Z(PP) 
  sdev[9]=1;    // Electron Tracking Efficiency	      
  sdev[10]=1;    // No Track Events Fraction           
  
  function_sys(onpar,sdev);

  return 0;
}


void function_sys(vector<double> onpar, vector<double> sdev)
{
  cout.setf(ios::right);
  cout.setf(ios::fixed);
  cout.precision(4);

  //double a[10];
  //for(int i=0; i<10; i++) a[i] = onpar[i];
  
  //dsigma/dy = (N(CC)-B(CC) + N(CP)-B(CP) + N(PP)-B(PP))/(L1*(A(CC)*E(CC) + A(CP)+E(CP)) + L2*(A(PP)*E(PP)));
  
  //independent systematic parameters :
  //B(CC), B(CP), B(PP), Luminosity, zvertex finding efficiency in Z(PP) - global factor - related in data
  //Electron ID efficiency - Z(PP) efficiency has a correlation with tracking
  //However, systematic uncertainty from Z(PP) ID efficiency is pretty small compared to BKG and tracking
  //Ignore the correlation
  //Material effect

  const int nbin=60;
  const int anbin=30;

  //parameter : 1(luminosity), 2(bkg,CC), 3(bkg,CP), 4(bkg,PP), 5(ID,cen),6(ID,plg), 7(tracking), 8(material,cen), 9(material,plg), 10(zvtx,pp)
  //parameter : 8 and 9 is relatively small compared to other sources

  ifstream in1("./file_list/adsigmady_table.dat");
  ifstream in2("./file_list/accp_eff_adsigmady.dat");
  ifstream in3("./file_list/nevents_dndy_abs.dat");

  double par[11][30];
  double fpar[11][30];

  double ybin[30]={0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,
		   1.05,1.15,1.25,1.35,1.45,1.55,1.65,1.75,1.85,1.95,
		   2.05,2.15,2.25,2.35,2.45,2.55,2.65,2.75,2.85,2.95};
  
  double dsb[30]={69.41,70.94,70.86,69.67,67.62,67.79,66.06,66.19,64.30,63.87,
		  61.90,61.26,58.19,55.58,53.10,49.90,46.46,40.82,37.01,33.20,
		  28.05,22.43,19.08,14.87,9.43,6.20,3.31,1.51,1.01,0.00};
  
  double statber[30]={0.73,0.74,0.73,0.72,0.70,0.70,0.69,0.69,0.68,0.68,
		      0.66,0.65,0.65,0.64,0.63,0.62,0.60,0.58,0.56,0.55,
		      0.53,0.51,0.51,0.52,0.48,0.48,0.47,0.41,0.60,0.00};
  
  double dsc[30]={69.41,70.94,70.86,69.68,67.62,67.79,66.07,66.19,64.31,63.88,
		  61.91,61.27,58.21,55.59,53.11,49.91,46.48,40.83,37.03,33.21,
		  28.05,22.43,19.07,14.85,9.40,6.16,3.27,1.48,0.96,0.00};
  
  double statcer[30]={0.73,0.74,0.73,0.72,0.70,0.70,0.69,0.69,0.68,0.68,
		      0.66,0.65,0.65,0.64,0.63,0.62,0.60,0.58,0.56,0.55,
		      0.53,0.51,0.51,0.52,0.48,0.48,0.46,0.41,0.57,0.00};
  
  double ccae[30]={60.993,57.381,52.278,46.649,41.351,35.915,29.499,21.522,12.955,5.7620,
		   1.3341,0.1061,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		   0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  
  double cpae[30]={3.9277,7.2482,13.347,20.632,27.654,33.987,40.332,47.533,56.115,64.045,
		   69.263,69.620,63.057,53.897,43.990,33.623,23.771,14.688,7.4906,2.9393,
		   0.7691,0.1558,0.0227,0.0025,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  
  double ppae[30]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		   0.1755,1.7418,6.4640,13.965,22.973,31.755,39.807,46.372,50.943,52.341,
		   50.055,43.791,35.987,27.748,20.469,13.377,7.6139,4.4123,1.4205,0.0000};  

  double ncc[30]={8510.0,8140.0,7420.0,6603.0,5679.0,4906.0,3936.0,2943.0,1684.0,764.0,
		  155.0,12.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
		  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};   
  
  double ncp[30]={564.0,1093.0,1956.0,2868.0,3754.0,4693.0,5425.0,6335.0,7346.0,8301.0,
		  8736.0,8696.0,7513.0,6120.0,4812.0,3505.0,2313.0,1276.0,610.0,223.0,
		  51.0,13.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
  
  double npp[30]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
		  24.0,219.0,750.0,1611.0,2495.0,3223.0,3805.0,3930.0,3900.0,3601.0,2894.0,
		  2046.0,1428.0,851.0,396.0,171.0,54.0,14.0,3.0,0.0};  

  double dev[anbin];
  double mdev[anbin];
  double pdev[anbin];

  double devf[anbin];
  double mdevf[anbin];
  double pdevf[anbin];

  //luminosity effect : 6% uncertainty in luminosity factor - global effect in total lum (L1 and L2)
  double lum_dev[anbin];
  
  for(int i=0; i<anbin; i++){
    
    if((ccae[i]+cpae[i]+ppae[i])!=0.0){
      lum_dev[i] = -1.0*dsc[i]*(ccae[i]+cpae[i]+ppae[i])*0.06*sdev[0]/(ccae[i]+cpae[i]+ppae[i]);
    }else{
      lum_dev[i] = 0.0;
    }
    
    par[0][i]=lum_dev[i];
    
  }

  //Background estimation : Z(CC) and Z(CP), and Z(PP) independent each other
  //Z(CP) and Z(PP) has three different components - 2track, 1track in fid, 1track in nofid
  //dijet and gammajet has 100% correlation in Z mass fit

  double bkgcc_dev[anbin];
  double bkgcp_dev[anbin];
  double bkgpp_dev[anbin];
 
  double accb[30] = {20.3874,19.2681,17.1894,16.0301,12.8321,12.2324,10.0338,7.0357,4.0775,2.1587,
                     0.4797,0.0800,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
                     0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};

  double acpb[30] = {8.9285,15.1512,26.8970,41.4853,55.1167,68.1634,75.1223,84.6674,97.6891,103.9387,
                     112.0996,117.1921,113.3532,100.8637,89.9579,77.3566,61.8287,42.6334,26.2479,
                     12.7574,4.3160,1.1950,0.1639,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};

  double appb[30] = {0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.2940,
                     6.8594,19.0954,53.7532,83.0881,104.1627,119.2991,152.0450,138.8262,130.6140,
                     82.9790,80.7053,48.6119,24.7114,9.7034,4.9723,3.6007,0.6007,0.1296,0.0266};
  
  double atcc[30] = {1.4251,1.3548,1.2238,1.1507,0.9485,0.9105,0.7706,0.5779,0.3827,0.2475,0.1021,
                     0.0403,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
                     0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};

  double atcp[30] = {1.9620,2.5232,4.2838,6.0612,7.4026,9.5898,10.3810,11.6099,12.8124,12.5361,14.7719,
                     15.0440,16.6066,13.4362,12.9912,11.8384,10.1300,6.8303,4.8285,2.9390,1.7060,0.8312,
                     0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};

  double atpp[30] = {0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.2442,
                     2.0957,3.0175,6.9488,9.3999,9.7043,10.6854,12.8717,12.5126,11.5310,7.8861,
                     7.3927,5.0873,3.1862,1.2439,0.9247,1.2323,0.3085,0.1283,0.0000};  

  double sumdev[anbin];

  for(int i=0; i<anbin; i++){
        
    //bin by bin has 100% correlation -> global shift in total rate
    
    if(ncc[i]<(accb[i]+atcc[i])) atcc[i]=0.0;
    if(ncp[i]<(acpb[i]+atcp[i])) atcp[i]=0.0;
    if(npp[i]<(appb[i]+atpp[i])) atpp[i]=0.0;

    if((ccae[i]+cpae[i]+ppae[i])!=0.0) bkgcc_dev[i] = -1.0*atcc[i]*sdev[1]/(ccae[i]+cpae[i]+ppae[i]);
    else bkgcc_dev[i] = 0.0;
       
    if((ccae[i]+cpae[i]+ppae[i])!=0.0) bkgcp_dev[i] = -1.0*atcp[i]*sdev[2]/(ccae[i]+cpae[i]+ppae[i]);
    else bkgcp_dev[i] = 0.0;
    
    if((ccae[i]+cpae[i]+ppae[i])!=0.0) bkgpp_dev[i] = -1.0*atpp[i]*sdev[3]/(ccae[i]+cpae[i]+ppae[i]);
    else bkgpp_dev[i] = 0.0;
    
    par[1][i]=bkgcc_dev[i];
    par[2][i]=bkgcp_dev[i];
    par[3][i]=bkgpp_dev[i];
    
    sumdev[i] = TMath::Sqrt(bkgcc_dev[i]*bkgcc_dev[i] + bkgcp_dev[i]*bkgcp_dev[i] + bkgpp_dev[i]*bkgpp_dev[i]);
    
  }
  
  //central electron efficiencies in Z(CC) and Z(CP) are 100 % correlated 
  //plug electron efficiencies in Z(CP) and Z(PP) are 100 % correlated
  //all bin by bin has 100 % correlation - global scale factor applied
  
  double cceff[anbin]    = {0.9066,0.9066,0.9035,0.9013,0.9006,0.9017,0.9038,0.9044,0.9089,0.9075,
			    0.9057,0.8958,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
			    0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  double cceffer[anbin]  = {0.0037,0.0037,0.0037,0.0037,0.0037,0.0037,0.0037,0.0037,0.0037,0.0037,
			    0.0037,0.0037,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
			    0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  double cpeff[anbin]    = {0.7135,0.7170,0.7283,0.7337,0.7298,0.7291,0.7263,0.7221,0.7216,0.7130,
			    0.7060,0.6975,0.6853,0.6706,0.6484,0.6259,0.6045,0.5879,0.5805,0.5664,
			    0.5571,0.5650,0.5308,1.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  double cpeffcer[anbin] = {0.0029,0.0029,0.0030,0.0030,0.0030,0.0030,0.0030,0.0029,0.0029,0.0029,
			    0.0029,0.0028,0.0028,0.0027,0.0026,0.0025,0.0024,0.0024,0.0023,0.0023,
			    0.0022,0.0023,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  double cpeffper[anbin] = {0.0034,0.0034,0.0035,0.0035,0.0035,0.0035,0.0035,0.0035,0.0035,0.0034,
			    0.0034,0.0034,0.0033,0.0032,0.0031,0.0030,0.0029,0.0028,0.0028,0.0027,
			    0.0027,0.0027,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000};
  double ppeff[anbin]    = {0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
			    0.8071,0.7847,0.7831,0.7777,0.7800,0.7783,0.7758,0.7667,0.7456,0.7324,
			    0.7210,0.7128,0.7076,0.7023,0.6972,0.6978,0.6531,0.6954,0.5359,0.0000};
  double ppeffer[anbin]  = {0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
			    0.0059,0.0045,0.0036,0.0027,0.0021,0.0018,0.0016,0.0014,0.0013,0.0013,
			    0.0013,0.0013,0.0014,0.0015,0.0016,0.0018,0.0018,0.0025,0.0028,0.0000};
  
  double ccfer[anbin];
  double cpcfer[anbin];
  double cppfer[anbin];
  double ppfer[anbin];
  
  double cpfer[anbin];
  
  double id_dev[anbin];
  double cid_dev[anbin];
  double pid_dev[anbin];
  
  for(int i=0; i<anbin; i++){
    
    if(cceff[i]!=0.0) ccfer[i]=cceffer[i]*sdev[4]/cceff[i];
    else ccfer[i]=0.0;
    
    if(cpeff[i]!=0.0) cpcfer[i]=cpeffcer[i]*sdev[4]/cpeff[i];
    else cpcfer[i]=0.0;
    
    if(cpeff[i]!=0.0) cppfer[i]=cpeffper[i]*sdev[5]/cpeff[i];
    else cppfer[i]=0.0;
    
    if(ppeff[i]!=0.0) ppfer[i]=ppeffer[i]*sdev[5]/ppeff[i];
    else ppfer[i]=0.0;
    
    if(cpeff[i]!=0.0) cpfer[i]=TMath::Sqrt(cpeffcer[i]*cpeffcer[i] + cpeffper[i]*cpeffper[i])/cpeff[i];
    else cpfer[i]=0.0;
    
    if((ccae[i]+cpae[i]+ppae[i])!=0.0){
      cid_dev[i] = -dsc[i]*(ccae[i]*ccfer[i] + cpae[i]*cpcfer[i])/(ccae[i]+cpae[i]+ppae[i]);
      pid_dev[i] = -dsc[i]*(cpae[i]*cppfer[i] + ppae[i]*ppfer[i])/(ccae[i]+cpae[i]+ppae[i]); 
    }else{
      cid_dev[i] = 0.0;
      pid_dev[i] = 0.0;
    }
    
    id_dev[i] = TMath::Sqrt( cid_dev[i]*cid_dev[i] + pid_dev[i]*pid_dev[i]);
    par[4][i] = cid_dev[i];
    par[5][i] = pid_dev[i];
    
  }
  
  
  //material effect : material_estimator(type, y) returned the acceptance ratio (new/standard)
  //estimated the acceptance ratio in global
  double cm_dev[anbin];
  double pm_dev[anbin];
  
  for(int i=0; i<anbin; i++){
    cm_dev[i] = sdev[6]*dsc[i]*(1.0/(material_estimator(1,ybin[i]))-1.0);
    pm_dev[i] = sdev[7]*dsc[i]*(1.0/(material_estimator(2,ybin[i]))-1.0);
    
    par[6][i] = cm_dev[i];
    par[7][i] = pm_dev[i];
  }
  
  //ZVertex finding efficiency    
  double vtxscl = 0.98168;
  double vtxscler = -0.004697;
  
  double vtx_dev[anbin];
  for(int i=0; i<30; i++){
    if((ccae[i]+cpae[i]+ppae[i])!=0.0) vtx_dev[i] = -1.0*dsc[i]*(ppae[i]*(vtxscler*sdev[8]/vtxscl))/(ccae[i]+cpae[i]+ppae[i]);
    else vtx_dev[i] = 0.0;
    
    par[8][i] = vtx_dev[i];
  }
  
  //Tracking efficiency
  
  double pnxeff[30]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.6998,
		     0.8300,0.8851,0.8954,0.8973,0.8888,0.8799,0.8674,0.8569,0.8408,0.8213,
		     0.7949,0.7646,0.7350,0.6911,0.6427,0.6250,0.6257,0.6142,0.0000};
  
  double pnxerr[30]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0077,
		     0.0040,0.0027,0.0023,0.0023,0.0026,0.0027,0.0030,0.0031,0.0034,0.0037,
		     0.0042,0.0047,0.0053,0.0058,0.0067,0.0072,0.0085,0.0140,0.0000};
  
  double tpnxerr[30]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0036,
		     0.0015,0.0010,0.0004,0.0002,0.0002,0.0002,0.0002,0.0004,0.0005,0.0005,
		     0.0007,0.0009,0.0012,0.0015,0.0021,0.0027,0.0038,0.0064,0.0000};

  double npnxerr[30]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,-0.0075,
		     -0.0041,-0.0026,-0.0025,-0.0025,-0.0028,-0.0030,-0.0033,-0.0034,-0.0036,-0.0040,
		     -0.0045,-0.0050,-0.0056,-0.0061,-0.0069,-0.0073,-0.0081,-0.0124,0.0000};

  double trk_dev[anbin];
  double pnxfer[anbin];
  double ntrk_dev[anbin];
  double npnxfer[anbin];
  for(int i=0; i<anbin; i++){
    
    if(pnxeff[i]!=0.0){
      pnxfer[i] = tpnxerr[i]*sdev[9]/pnxeff[i];
      npnxfer[i] = npnxerr[i]*sdev[10]/pnxeff[i];
    }else{
      pnxfer[i] = 0.0;
      npnxfer[i] = 0.0;
    }

    if((ccae[i]+cpae[i]+ppae[i])!=0.0){
      trk_dev[i] = -1.0*dsc[i]*(ppae[i]*pnxfer[i])/(ccae[i]+cpae[i]+ppae[i]);
      ntrk_dev[i] = -1.0*dsc[i]*(ppae[i]*npnxfer[i])/(ccae[i]+cpae[i]+ppae[i]);
    }else{
      trk_dev[i] = 0.0;
      ntrk_dev[i] = 0.0;
    }
    par[9][i] = trk_dev[i];      
    par[10][i] = ntrk_dev[i];      
  }

  double fpstat[anbin];
  double fmstat[anbin];
  double pstat[anbin];
  double mstat[anbin];
  
  for(int i=0; i<anbin; i++){
    if(dsc[i]!=0.0) fpstat[i] = statcer[i]/dsc[i];
    else fpstat[i] = 0.0;

    pstat[i] = statcer[i];
    mstat[i] = -1.0*statcer[i];

    fmstat[i] = -1.0*fpstat[i];

    for(int j=0; j<11; j++){
      if(dsc[i]!=0.0) fpar[j][i] = par[j][i]/dsc[i];
      else fpar[j][i] = 0.0;

      if(par[j][i]==0) par[j][i]=0.0;
    }
  }  

  cout.precision(1);  
  cout<<"Systematics are calculated for :"<<endl;
  if(onpar[0]==1) cout<<"Luminosity factor            with "<<sdev[0]<<" sigma change"<<endl;
  if(onpar[1]==1) cout<<"Background for Z(CC)         with "<<sdev[1]<<" sigma change"<<endl;
  if(onpar[2]==1) cout<<"Background for Z(CP)         with "<<sdev[2]<<" sigma change"<<endl;
  if(onpar[3]==1) cout<<"Background for Z(PP)         with "<<sdev[3]<<" sigma change"<<endl;
  if(onpar[4]==1) cout<<"Central electron ID effi.    with "<<sdev[4]<<" sigma change"<<endl;
  if(onpar[5]==1) cout<<"Plug electron ID effi.       with "<<sdev[5]<<" sigma change"<<endl;
  if(onpar[6]==1) cout<<"Central extra material       with "<<sdev[6]<<" sigma change"<<endl;
  if(onpar[7]==1) cout<<"Plug extra material          with "<<sdev[7]<<" sigma change"<<endl;
  if(onpar[8]==1) cout<<"ZVertex finding effi.        with "<<sdev[8]<<" sigma change"<<endl;
  if(onpar[9]==1) cout<<"Plug Electron Tracking effi. with "<<sdev[9]<<" sigma change"<<endl;
  if(onpar[10]==1) cout<<"No Track Events fraction     with "<<sdev[10]<<" sigma change"<<endl;
  cout.precision(4);  
  double total[anbin];
  double tot_sum=0.0;
  cout<<"---------------------------------------------------------------------------------------------------------------------------------------------------------"<<endl;
  cout<<"  y bin   |";
  cout<<"   sigma   |";
  cout<<"   stat.  |";
  if(onpar[0]==1) cout<<"  lum    |";
  if(onpar[1]==1) cout<<"  B(CC)  |";
  if(onpar[2]==1) cout<<"  B(CP)  |";
  if(onpar[3]==1) cout<<"  B(PP)  |";
  if(onpar[4]==1) cout<<"  CID    |";
  if(onpar[5]==1) cout<<"  PID    |";
  if(onpar[6]==1) cout<<"  CMat   |";
  if(onpar[7]==1) cout<<"  PMat   |";
  if(onpar[8]==1) cout<<"  ZVtx   |";
  if(onpar[9]==1) cout<<"  Trkeff |";
  if(onpar[10]==1) cout<<"  NoTrk  |";
  cout<<" Tot errors "<<endl;
  cout<<"---------------------------------------------------------------------------------------------------------------------------------------------------------"<<endl;
  for(int i=0; i<anbin; i++){
    total[i] = TMath::Sqrt(TMath::Power(onpar[0]*par[0][i],2)+TMath::Power(onpar[1]*par[1][i],2)+
			   TMath::Power(onpar[2]*par[2][i],2)+TMath::Power(onpar[3]*par[3][i],2)+  
			   TMath::Power(onpar[4]*par[4][i],2)+TMath::Power(onpar[5]*par[5][i],2)+ 
			   TMath::Power(onpar[6]*par[6][i],2)+TMath::Power(onpar[7]*par[7][i],2)+ 
			   TMath::Power(onpar[8]*par[8][i],2)+TMath::Power(onpar[9]*par[9][i],2)+
			   TMath::Power(onpar[10]*par[10][i],2));
    
    cout<<"  "<<ybin[i]<<"  | ";
    if(dsc[i]>=10) cout<<" " <<dsc[i]<<"  | ";
    else cout<<"  " <<dsc[i]<<"  | ";
    cout<<" " <<statcer[i]<<"  | ";
    for(int t=0; t<11; t++) {
      if(onpar[t]==1){
	if(par[t][i]>=0) cout <<" "<<par[t][i]<<" | ";
	else cout<<par[t][i]<<" | ";
      }
    }
    if(total[i]>=0) cout <<" "<<total[i]<<endl;
    else cout<<total[i]<<endl;

    pdev[i] = total[i];
    mdev[i] = -1.0*total[i];
    
    if(dsc[i]!=0.0){
      pdevf[i] = total[i]/dsc[i];    
      mdevf[i] = -1.0*total[i]/dsc[i];    
    }else{
      pdevf[i] = 0.0;
      mdevf[i] = 0.0;
    }

    tot_sum+=total[i];

  }
  cout<<"---------------------------------------------------------------------------------------------------------------------------------------------------------"<<endl;
  cout<<"Total systematic uncertainty = "<<tot_sum*0.1*2.0<<endl;
  cout<<"---------------------------------------------------------------------------------------------------------------------------------------------------------"<<endl;
  gStyle->SetFillColor(10);
  
  TLatex lat_comp;
  lat_comp.SetTextSize(0.06);
  lat_comp.SetTextColor(1);

  TLine *l1 = new TLine(0.0,0.0,3.0,0.0);
  l1->SetLineWidth(2);

  TCanvas * Zee_comp1 = new TCanvas("Zee_comp1","Systematic Deviation",1200,900);
  Zee_comp1->Divide(4,3);
  Zee_comp1->SetFillColor(10);
  Zee_comp1->SetFillStyle(4000);
  
  TGraph *syscomp[11];
    
  for(int t=0; t<11; t++){
    
    syscomp[t] = new TGraph(anbin,ybin,par[t]);
    
    Zee_comp1->cd(t+1);
    
    TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
    padmr->Draw();
    padmr->cd();
    if(t==0){
      padmr->Range(0.0,-5.0,3.0,5.0);
      padmr->DrawFrame(0.0,-5.0,3.0,5.0);
    }else{
      padmr->Range(0.0,-1.0,3.0,1.0);
      padmr->DrawFrame(0.0,-1.0,3.0,1.0);
    }
    padmr->GetFrame()->SetFillColor(10);
    padmr->SetLeftMargin(0.18);
    padmr->SetTopMargin(0.1);
    padmr->SetBottomMargin(0.15);
    padmr->SetGrid();
    
    TH1 *hframe = new TH1F("hframe","",1000,0.0,3.0);
    if(t==0){
      hframe->SetMinimum(-5.0);
      hframe->SetMaximum(5.0);
    }else{
      hframe->SetMinimum(-1.0);
      hframe->SetMaximum(1.0);
    }
    
    hframe->SetDirectory(0);
    hframe->SetStats(0);
    hframe->GetXaxis()->CenterTitle(1);
    hframe->GetYaxis()->CenterTitle(1);
    hframe->GetXaxis()->SetTitle("Boson Rapidity");
    hframe->GetYaxis()->SetTitle("#Delta #sigma_{i}");
    hframe->GetXaxis()->SetTitleSize(0.05);
    hframe->GetYaxis()->SetTitleSize(0.07);
    hframe->GetYaxis()->SetTitleOffset(1.0);
    hframe->Draw("a");
    
    if(t==0) lat_comp.DrawLatex(1.0,3.5,"Luminosity");
    else if(t==1) lat_comp.DrawLatex(1.0,0.7,"CC Bkg");
    else if(t==2) lat_comp.DrawLatex(1.0,0.7,"CP Bkg");
    else if(t==3) lat_comp.DrawLatex(1.0,0.7,"PP Bkg");
    else if(t==4) lat_comp.DrawLatex(1.0,0.7,"Cen ID");
    else if(t==5) lat_comp.DrawLatex(1.0,0.7,"Plg ID");
    else if(t==6) lat_comp.DrawLatex(1.0,0.7,"Cen Mat");
    else if(t==7) lat_comp.DrawLatex(1.0,0.7,"Plg Mat");
    else if(t==8) lat_comp.DrawLatex(1.0,0.7,"ZVtx");
    else if(t==9) lat_comp.DrawLatex(1.0,0.7,"Trk Efficiency");
    else if(t==10) lat_comp.DrawLatex(1.0,0.7,"No Trk Fraction");
  
    if(onpar[t]==1){  
      syscomp[t]->SetLineColor(4);
      syscomp[t]->SetLineWidth(2);
      syscomp[t]->SetMarkerColor(4);
      syscomp[t]->Draw("CP");
      l1->Draw();  
    }else lat_comp.DrawLatex(0.7,0.0,"No Contribution");
  
  }

  Zee_comp1->cd(12);
  TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
  padmr->Draw();
  padmr->cd();


  TCanvas * Zee_comp2 = new TCanvas("Zee_comp2","Systematic Fraction",1200,900);
  Zee_comp2->Divide(4,3);
  Zee_comp2->SetFillColor(10);
  Zee_comp2->SetFillStyle(4000);
  
  TGraph *syscompf[11];
    
  for(int t=0; t<11; t++){
    
    syscompf[t] = new TGraph(anbin,ybin,fpar[t]);
    
    Zee_comp2->cd(t+1);
    
    TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
    padmr->Draw();
    padmr->cd();
    if(t==0 || t==3){
      padmr->Range(0.0,-0.1,3.0,0.1);
      padmr->DrawFrame(0.0,-0.1,3.0,0.1);
    }else{
      padmr->Range(0.0,-0.015,3.0,0.015);
      padmr->DrawFrame(0.0,-0.015,3.0,0.015);
    }
    padmr->GetFrame()->SetFillColor(10);
    padmr->SetLeftMargin(0.18);
    padmr->SetTopMargin(0.1);
    padmr->SetBottomMargin(0.15);
    padmr->SetGrid();
    
    TH1 *hframe = new TH1F("hframe","",1000,0.0,3.0);
    if(t==0 || t==3){
      hframe->SetMinimum(-0.1);
      hframe->SetMaximum(0.1);
    }else{
      hframe->SetMinimum(-0.02);
      hframe->SetMaximum(0.02);
    }
    
    hframe->SetDirectory(0);
    hframe->SetStats(0);
    hframe->GetXaxis()->CenterTitle(1);
    hframe->GetYaxis()->CenterTitle(1);
    hframe->GetXaxis()->SetTitle("Boson Rapidity");
    hframe->GetYaxis()->SetTitle("#Delta #sigma_{i}/#sigma_{i}");
    hframe->GetXaxis()->SetTitleSize(0.05);
    hframe->GetYaxis()->SetTitleSize(0.07);
    hframe->GetYaxis()->SetTitleOffset(1.1);
    hframe->Draw("a");
    
    if(t==0) lat_comp.DrawLatex(1.0,0.07,"Luminosity");
    else if(t==1) lat_comp.DrawLatex(1.0,0.014,"CC Bkg");
    else if(t==2) lat_comp.DrawLatex(1.0,0.014,"CP Bkg");
    else if(t==3) lat_comp.DrawLatex(1.0,0.070,"PP Bkg");
    else if(t==4) lat_comp.DrawLatex(1.0,0.014,"Cen ID");
    else if(t==5) lat_comp.DrawLatex(1.0,0.014,"Plg ID");
    else if(t==6) lat_comp.DrawLatex(1.0,0.014,"Cen Mat");
    else if(t==7) lat_comp.DrawLatex(1.0,0.014,"Plg Mat");
    else if(t==8) lat_comp.DrawLatex(1.0,0.014,"ZVtx");
    else if(t==9) lat_comp.DrawLatex(1.0,0.014,"Trk Efficiency");
    else if(t==10) lat_comp.DrawLatex(1.0,0.014,"No Trk Fraction");
  
    if(onpar[t]==1){  
      syscompf[t]->SetLineColor(4);
      syscompf[t]->SetLineWidth(2);
      syscompf[t]->SetMarkerColor(4);
      syscompf[t]->Draw("CP");
      l1->Draw();  
    }else lat_comp.DrawLatex(0.7,0.0,"No Contribution");
  
  }

  Zee_comp2->cd(12);
  TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
  padmr->Draw();
  padmr->cd();


  TCanvas * Zee = new TCanvas("Zee","Zee anal",550,800);
  Zee->Divide(1,2);
  Zee->SetFillColor(10);
  Zee->SetFillStyle(4000);

  TGraph * psysgrp = new TGraph(anbin,ybin,pdev);
  TGraph * msysgrp = new TGraph(anbin,ybin,mdev);
  TGraph * pstsgrp = new TGraph(anbin,ybin,pstat);
  TGraph * mstsgrp = new TGraph(anbin,ybin,mstat);

  TGraph * psysgrpf = new TGraph(anbin,ybin,pdevf);
  TGraph * msysgrpf = new TGraph(anbin,ybin,mdevf);
  TGraph * pstsgrpf = new TGraph(anbin,ybin,fpstat);
  TGraph * mstsgrpf = new TGraph(anbin,ybin,fmstat);

  Zee->cd(1);
  TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
  padmr->Draw();
  padmr->cd();
  if(onpar[0]==1){
    padmr->Range(0.0,-5.0,3.0,5.0);
    padmr->DrawFrame(0.0,-5.0,3.0,5.0);
  }else{
    padmr->Range(0.0,-1.0,3.0,1.0);
    padmr->DrawFrame(0.0,-1.0,3.0,1.0);
  }
  padmr->GetFrame()->SetFillColor(10);
  padmr->SetLeftMargin(0.18);
  padmr->SetRightMargin(0.15);
  padmr->SetTopMargin(0.1);
  padmr->SetBottomMargin(0.15);
  padmr->SetGrid();
  
  TH1 *hframe = new TH1F("hframe","",1000,0.0,3.0);
  if(onpar[0]==1){
    hframe->SetMinimum(-5.0);
    hframe->SetMaximum(5.0);
  }else{
    hframe->SetMinimum(-1.0);
    hframe->SetMaximum(1.0);
  }
  hframe->SetDirectory(0);
  hframe->SetStats(0);
  hframe->GetXaxis()->CenterTitle(1);
  hframe->GetYaxis()->CenterTitle(1);
  hframe->GetXaxis()->SetTitle("Boson Rapidity");
  hframe->GetYaxis()->SetTitle("#Delta #sigma_{i}");
  hframe->GetXaxis()->SetTitleSize(0.05);
  hframe->GetYaxis()->SetTitleSize(0.07);
  hframe->GetYaxis()->SetTitleOffset(0.8);
  hframe->Draw("a");
  
  psysgrp->SetLineColor(4);
  psysgrp->SetLineWidth(2);
  psysgrp->SetMarkerColor(4);
  psysgrp->Draw("CP");

  msysgrp->SetLineColor(4);
  msysgrp->SetLineWidth(2);
  msysgrp->SetMarkerColor(4);
  msysgrp->Draw("CP");

  pstsgrp->SetLineColor(2);
  pstsgrp->SetLineStyle(2);
  pstsgrp->SetLineWidth(2);
  pstsgrp->SetMarkerColor(2);
  pstsgrp->Draw("CP");

  mstsgrp->SetLineColor(2);
  mstsgrp->SetLineStyle(2);
  mstsgrp->SetLineWidth(2);
  mstsgrp->SetMarkerColor(2);
  mstsgrp->Draw("CP");

  TLine *l1 = new TLine(0.0,0.0,3.0,0.0);
  l1->SetLineWidth(2);
  l1->Draw();  

  TLatex lat;
  lat.SetTextSize(0.03);
  lat.SetTextColor(1);
  if(onpar[0]==1){
    if(onpar[0]==1) lat.DrawLatex(3.05,4.5,"Luminosity");
    if(onpar[1]==1) lat.DrawLatex(3.05,4.0,"CC Bkg");
    if(onpar[2]==1) lat.DrawLatex(3.05,3.5,"CP Bkg");
    if(onpar[3]==1) lat.DrawLatex(3.05,3.0,"PP Bkg");
    if(onpar[4]==1) lat.DrawLatex(3.05,2.5,"Cen ID");
    if(onpar[5]==1) lat.DrawLatex(3.05,2.0,"Plg ID");
    if(onpar[6]==1) lat.DrawLatex(3.05,1.0,"Cen Mat");
    if(onpar[7]==1) lat.DrawLatex(3.05,0.5,"Plg Mat");
    if(onpar[8]==1) lat.DrawLatex(3.05,0.0,"ZVtx");
    if(onpar[9]==1) lat.DrawLatex(3.05,1.5,"Tracking");
    if(onpar[10]==1) lat.DrawLatex(3.05,1.5,"Tracking");
  }else{
    if(onpar[1]==1) lat.DrawLatex(3.05,0.8,"CC Bkg");
    if(onpar[2]==1) lat.DrawLatex(3.05,0.7,"CP Bkg");
    if(onpar[3]==1) lat.DrawLatex(3.05,0.6,"PP Bkg");
    if(onpar[4]==1) lat.DrawLatex(3.05,0.5,"Cen ID");
    if(onpar[5]==1) lat.DrawLatex(3.05,0.4,"Plg ID");
    if(onpar[6]==1) lat.DrawLatex(3.05,0.2,"Cen Mat");
    if(onpar[7]==1) lat.DrawLatex(3.05,0.1,"Plg Mat");
    if(onpar[8]==1) lat.DrawLatex(3.05,0.0,"ZVtx");
    if(onpar[9]==1) lat.DrawLatex(3.05,0.3,"Tracking");
    if(onpar[10]==1) lat.DrawLatex(3.05,0.3,"Tracking");
  }

  TLegend * tld = new TLegend(0.5,0.8,0.84,0.95);
  tld->SetBorderSize(0);
  tld->SetFillColor(10);
  tld->AddEntry(psysgrp,"Systematic Error","l");
  tld->AddEntry(pstsgrp,"Statistic Error","l");
  tld->Draw();

  Zee->cd(2);
  TPad * padmr = new TPad("padm","y",0.0,0.0,1.0,1.0,10);
  padmr->Draw();
  padmr->cd();
  if(onpar[0]==1){
    padmr->Range(0.0,-0.20,3.0,0.20);
    padmr->DrawFrame(0.0,-0.20,3.0,0.20);
  }else{
    padmr->Range(0.0,-0.01,3.0,0.01);
    padmr->DrawFrame(0.0,-0.01,3.0,0.01);
  }
  padmr->GetFrame()->SetFillColor(10);
  padmr->SetLeftMargin(0.18);
  padmr->SetRightMargin(0.15);
  padmr->SetTopMargin(0.1);
  padmr->SetBottomMargin(0.15);
  padmr->SetGrid();
  
  TH1 *hframe = new TH1F("hframe","",1000,0.0,3.0);
  if(onpar[0]==1 || onpar[3]==1 || onpar[6]==1){
    hframe->SetMinimum(-0.20);
    hframe->SetMaximum(0.20);
  }else{
    hframe->SetMinimum(-0.01);
    hframe->SetMaximum(0.01);
  }
  hframe->SetDirectory(0);
  hframe->SetStats(0);
  hframe->GetXaxis()->CenterTitle(1);
  hframe->GetYaxis()->CenterTitle(1);
  hframe->GetXaxis()->SetTitle("Boson Rapidity");
  hframe->GetYaxis()->SetTitle("#Delta #sigma_{i}/#sigma_{i}");
  hframe->GetXaxis()->SetTitleSize(0.05);
  hframe->GetYaxis()->SetTitleSize(0.07);
  hframe->GetYaxis()->SetTitleOffset(0.8);
  hframe->Draw("a");
  
  psysgrpf->SetLineColor(4);
  psysgrpf->SetLineWidth(2);
  psysgrpf->SetMarkerColor(4);
  psysgrpf->Draw("CP");

  msysgrpf->SetLineColor(4);
  msysgrpf->SetLineWidth(2);
  msysgrpf->SetMarkerColor(4);
  msysgrpf->Draw("CP");

  pstsgrpf->SetLineColor(2);
  pstsgrpf->SetLineStyle(2);
  pstsgrpf->SetLineWidth(2);
  pstsgrpf->SetMarkerColor(2);
  pstsgrpf->Draw("CP");

  mstsgrpf->SetLineColor(2);
  mstsgrpf->SetLineStyle(2);
  mstsgrpf->SetLineWidth(2);
  mstsgrpf->SetMarkerColor(2);
  mstsgrpf->Draw("CP");

  l1->Draw();

  if(onpar[0]==1 || onpar[3]==1 || onpar[6]==1){
    if(onpar[0]==1) lat.DrawLatex(3.05,0.20-0.20*0.1,"Luminosity");
    if(onpar[1]==1) lat.DrawLatex(3.05,0.20-0.20*0.2,"CC Bkg");
    if(onpar[2]==1) lat.DrawLatex(3.05,0.20-0.20*0.3,"CP Bkg");
    if(onpar[3]==1) lat.DrawLatex(3.05,0.20-0.20*0.4,"PP Bkg");
    if(onpar[4]==1) lat.DrawLatex(3.05,0.20-0.20*0.5,"Cen ID");
    if(onpar[5]==1) lat.DrawLatex(3.05,0.20-0.20*0.6,"Plg ID");
    if(onpar[6]==1) lat.DrawLatex(3.05,0.20-0.20*0.8,"Cen Mat");
    if(onpar[7]==1) lat.DrawLatex(3.05,0.20-0.20*0.9,"Plg Mat");
    if(onpar[8]==1) lat.DrawLatex(3.05,0.20-0.20*1.0,"ZVtx");
    if(onpar[9]==1) lat.DrawLatex(3.05,0.20-0.20*0.7,"Tracking");
    if(onpar[10]==1) lat.DrawLatex(3.05,0.20-0.20*0.7,"Tracking");
  }else{
    if(onpar[1]==1) lat.DrawLatex(3.05,0.01-0.01*0.2,"CC Bkg");
    if(onpar[2]==1) lat.DrawLatex(3.05,0.01-0.01*0.3,"CP Bkg");
    if(onpar[3]==1) lat.DrawLatex(3.05,0.01-0.01*0.4,"PP Bkg");
    if(onpar[4]==1) lat.DrawLatex(3.05,0.01-0.01*0.5,"Cen ID");
    if(onpar[5]==1) lat.DrawLatex(3.05,0.01-0.01*0.6,"Plg ID");
    if(onpar[6]==1) lat.DrawLatex(3.05,0.01-0.01*0.8,"Cen Mat");
    if(onpar[7]==1) lat.DrawLatex(3.05,0.01-0.01*0.9,"Plg Mat");
    if(onpar[8]==1) lat.DrawLatex(3.05,0.01-0.01*1.0,"ZVtx");
    if(onpar[9]==1) lat.DrawLatex(3.05,0.01-0.01*0.7,"Tracking");
    if(onpar[10]==1) lat.DrawLatex(3.05,0.01-0.01*0.7,"Tracking");
  }

}


double material_estimator(int type, double x){
  
  //central material : extra 1% X0 more : total acceptance ratio
  Double_t cenpar[2] = {0.9953,0.001725};
  Double_t cenparer[2] = {0.002587,0.002075};
  Double_t cenfun = cenpar[0] + cenpar[1]*TMath::Abs(x);
  
  //central material : extra 1/6 X0 more : total acceptance ratio
  Double_t plgpar[2] = {0.9966,0.00103};
  Double_t plgparer[2] = {0.002609,0.002212};
  Double_t plgfun = plgpar[0] + plgpar[1]*TMath::Abs(x);

  if(type==1) return cenfun;
  else if(type==2) return plgfun;
}