{
// >> example2_2.C
// takes a z-expansion sample and a reweighting file from grwghtnp (covariance reweighting)
// generates and plots a histogram of Q2 with error bars from the weight file
// error bars are calculated from taking the standard deviation of the 
//  reweighted histogram values
//
gStyle->SetOptStat(0); //get rid of statistics box

// input files
TFile* f0 = TFile::Open("gntp.1.gst.root");
TFile* fwc = TFile::Open("wght.1.cov.root");

// get trees
TTree* t0 = f0->Get("gst");
// weight file has a covariance reweighting tree
// to find this, try fwc->ls()
TTree* twc = fwc->Get("covrwt");

// covariant reweighting parameters
int      nev = t0->GetEntries();
int      evtn;
int      ntwk;
double   q2;
TArrayD *wghtc; //< weights are stored in a TArrayD, access with wghtc.At(i)


// set data addresses to load to
t0->SetBranchAddress("Q2",&q2);
twc->SetBranchAddress("eventnum",&evtn);
twc->SetBranchAddress("n_tweaks",&ntwk);
twc->SetBranchAddress("weights",&wghtc);

twc->GetEntry(0); //< load the first event
const int n_tweaks = ntwk;

// create histograms
int    nbinsx = 8;  // number of bins in histogram
double xlow   = 0.; // x lower extreme
double xup    = .8; // x upper extreme
double xblow[nbinsx+1]; // x axis bin boundaries
double binws[nbinsx];   // bin widths
for(int i=0;i<nbinsx;i++){
 binws[i] = (xup-xlow)/nbinsx/2.;
 xblow[i] = i*(xup-xlow)/nbinsx;
}
xblow[nbinsx] = nbinsx*(xup-xlow)/nbinsx;

// nominal histogram
TH1F* hnom   = new TH1F("hnom","",nbinsx,xlow,xup);
// covariance randomized histograms
TH1F* hcov[n_tweaks];
for(int i=0;i<n_tweaks;i++)
 {
   stringstream hcname("");
   hcname << "hcov" << i; // need a unique name for each histogram, so just number them
   hcov[i]  = new TH1F(hcname.str().c_str(),"",nbinsx,xblow);
 }
hnom->Sumw2();

double wsum[n_tweaks] = {0.}; //< store the sum of the weights
// loop through events
for (int iev=0;iev<nev;iev++)
{
  t0->GetEntry(iev);
  twc->GetEntry(iev);
  
  // fill histograms
  hnom->Fill(q2);
  for(int itwk=0;itwk<n_tweaks;itwk++) {
    hcov[itwk]->Fill(q2,wghtc.fArray[itwk]);
    wsum[itwk] += wghtc.fArray[itwk];
  }
}

// calculate standard deviation for each Q2 bin
// indexing on histograms starts from 1 due to underflow bin
double sdevc[(const int)nbinsx] = {0.};
for(int i=1;i<=nbinsx+1;i++)
{
  for(itwk=0;itwk<n_tweaks;itwk++){
    sdevc[i-1] += TMath::Power(hnom->GetBinContent(i)-hcov[itwk]->GetBinContent(i),2);
  }
  sdevc[i-1] = TMath::Sqrt(sdevc[i-1]/(double)n_tweaks);
}

// create array holding central value data
double xvalc[(const int) nbinsx];
double yvalc[(const int) nbinsx];
for(int i=0;i<nbinsx;i++)
{
  xvalc[i] = (i+.5)*(xup-xlow)/nbinsx;
  yvalc[i] = hnom->GetBinContent(i+1);
}

// create a graph and plot
TCanvas *c11 = new TCanvas("c11","");
gc = new TGraphErrors(nbinsx,xvalc,yvalc,binws,sdevc);
gc->SetTitle("");
gc->GetXaxis()->SetTitle("Q^{2}[GeV^{2}]");
gc->Draw("");
}