void NDpairStudy()
{
  gROOT->ProcessLine(".L ND_pair.C");
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  gStyle->SetPadGridX(1);
  gStyle->SetPadGridY(1);  

  TGaxis::SetMaxDigits(3);

  const int bins = 20;
  float NT[bins];
  float NTbinWidth[bins];
  float minNT[bins];
  float maxNT[bins];

  float NT_ave[bins];
  float NTbinWidthPM[bins];
  float *upstream_up;
  float current_upstream[bins];
  float current_upstream_err[bins];

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

//Binning
NT[i] = 0.5+float(i);
NTbinWidth[i] = 1.;
minNT[i] = NT[i] - 0.5;
maxNT[i] = NT[i] + 0.5;


NT_ave[i] = minNT[i]+0.5;
NTbinWidthPM[i] = NTbinWidth[i]/2.0;

    cout << "i = "<<i<<" minT = "<<minNT[i]<<", maxT ="<<maxNT[i]<<endl;
    
    upstream_up = ND_pair(minNT[i],maxNT[i]);

    current_upstream[i] = upstream_up[0];
    current_upstream_err[i] = upstream_up[1];
     

  }

TCanvas *c2 = new TCanvas("c2","Neutron Detector Asymmetry Study",900,900);
  c2->Divide(1,1);
c2->cd(1);
gStyle->SetOptStat(0);
  gPad->SetLogy(0);
  char graphname[256];
  sprintf(graphname,"Beam current stability defined by the pair spectrometer ");
  hr1 = new TH2F("hr1",graphname,500,-1,92,500,0,3000);

  gStyle->SetTitleX(0.18);
  char yTitle[256];
  char xTitle[256]; 
  sprintf(xTitle,"time,min");  
  sprintf(yTitle,"e+ #");  
  hr1->GetYaxis()->SetLabelColor(1);
  hr1->GetYaxis()->SetTitleColor(1);
  hr1->SetYTitle(yTitle);
  hr1->GetYaxis()->SetTitleOffset(1.3);
  hr1->GetXaxis()->SetLabelColor(1);
  hr1->GetXaxis()->SetTitleColor(1);
  hr1->SetXTitle(xTitle);
  hr1->Draw();

// TGraphErrors *gr1 = new TGraphErrors(bins,maxNE, asym_TDC1);
TGraphErrors *gr1 = new TGraphErrors(bins, NT_ave, current_upstream, NTbinWidthPM, current_upstream_err);
  gr1->SetMarkerColor(kBlue);
  gr1->SetMarkerStyle(20);   
  gr1->SetMarkerSize(1);
  //  gr1->SetFillColor(0);
  gr1->Draw("P,same");
c2->Update();


}