// conatins some basic histogramming utilities
// CJJ with bug fixes from DW and KC
// June 2004

// averages numHists 1-d histos and returns the average. the 
// errors ineach bin on the average are the std deviations (bin by bin)
// of the input hists
TH1F* averageHist(TH1F** hin, Int_t numHists)  {  
  TH1F* h = (TH1F*)hin[0]->Clone();
  h->SetName("have");
  h->SetTitle("Average Histogram");
  h->Reset();
  Int_t i;
  Int_t nbinsx = h->GetNbinsX();
  Int_t k;
  Double_t mean, stdDev;
  Double_t*x = new Double_t[numHists];

  for( i=1 ; i<=nbinsx ; i++ ) {
    for( k=0 ; k<numHists ; k++ ) {
      x[k] = hin[k]->GetBinContent(i); 
    }
    meanStdDev(x, numHists, &mean, &stdDev);
    h->SetBinContent(i, mean);
    h->SetBinError(i, stdDev);
  }
  delete x;
  return h;
}  

// as averageHist but for 2d histograms
TH2F* average2d(TH2F** hin, Int_t numHists) {
  TH2F* h = (TH2F*)hin[0]->Clone();
  h->SetName("have");
  h->SetTitle("Average Histogram");
  h->Reset();
  Int_t i,j, bin;
  Int_t nbinsx = h->GetNbinsX();
  Int_t nbinsy = h->GetNbinsY();
  Int_t k;
  Double_t mean, stdDev;
  Double_t*x = new Double_t[numHists];

  for( i=1 ; i<=nbinsx ; i++ ) {
    for( j=1 ; j<=nbinsy ; j++ ) {
      bin = h->GetBin(i,j);
      for( k=0 ; k<numHists ; k++ ) x[k] = hin[k]->GetBinContent(bin);
      meanStdDev(x, numHists, &mean, &stdDev);
      h->SetBinContent(bin, mean);
      h->SetBinError(bin, stdDev);
    }
  }
  delete x;
  return h;
}

//calc mean asn std dev of an array of doubles...
void meanStdDev(Double_t* x, Int_t numData, Double_t* mean, Double_t* stdDev) {
  int i;
  if( numData == 0 ) {
    *mean = 0;
    *stdDev = 1e32;
  }
  Double_t sum = 0;
  for( i=0 ; i<numData ; i++ ) sum += x[i];
  *mean = sum/numData;

  if( numData==1 ) *stdDev = 1e32;;
  sum = 0;
  for( i=0 ; i<numData ; i++ ) {
    sum += pow((x[i]- (*mean)), 2);
  }
  *stdDev = sqrt(sum/(numData-1));
}

// for a given input histogram, returns the cumulative
// histogram CJJ
TH1D* integrateHist(TH1D* hin) {
  char htitle[200];
  sprintf(htitle,"Cumulative of %s",hin->GetTitle());
  TH1D* hint = new TH1D("hint",htitle,hin->GetNbinsX(), 
			hin->GetXaxis()->GetXmin(),
			hin->GetXaxis()->GetXmax());
  int i;
  Double_t sum = 0.0;
  Double_t scaleFactor = hin->Integral();
  for( i=1 ; i<=hin->GetNbinsX() ; i++ ) {
    sum += hin->GetBinContent(i);
    hint->SetBinContent(i,sum/scaleFactor);
  }
  return hint;
}