RootTools.cxx

#include "RootTools.h"





//---------------------------------------------------------------------------------------------------------
Double_t RootTools::getSumOfYVals(TGraph* gr){
  Double_t sum = 0;
  for(int i=0; i<gr->GetN(); i++){
    sum += gr->GetY()[i];
  }
  return sum;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::printArray(int n, double* array, TString delimiter, TString start ,TString end){
  std::cout << start.Data();
  for(int i=0; i<n; i++){
    std::cout << array[i];
    if(i<n-1) std::cout << delimiter.Data();
  }  
  std::cout << end.Data();
}





//---------------------------------------------------------------------------------------------------------
void RootTools::printYVals(TGraph* gr, TString delimiter, TString start, TString end){
  printArray(gr->GetN(), gr->GetY(), delimiter, start, end);
}





//---------------------------------------------------------------------------------------------------------
void RootTools::printXVals(TGraph* gr, TString delimiter, TString start, TString end){
  printArray(gr->GetN(), gr->GetX(), delimiter, start, end);
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeDerivativeTGraph(TGraph* gr){
  TGraph* grDer = new TGraph();
  Int_t nDer = gr->GetN()-1;
  grDer->Set(nDer);
  for(int samp=0; samp<nDer; samp++){
    Double_t dy = gr->GetY()[samp+1]-gr->GetY()[samp];
    Double_t dx = gr->GetX()[samp+1]-gr->GetX()[samp];
    grDer->SetPoint(samp, gr->GetX()[samp], dy/dx);
  }
  return grDer;
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeUnwrappedCorrelationGraph(TGraph* gr){
  /* Unwraps a circularly correlated array so the n > 2 points get -ve values */
  Int_t n = gr->GetN();
  Double_t dt = gr->GetX()[1] - gr->GetX()[0];
  TGraph* grUnwrapped = new TGraph();
  grUnwrapped->Set(n);
  Int_t samp2=0;
  for(Int_t samp=n/2; samp<n; samp++){
    grUnwrapped->SetPoint(samp2, (samp - n)*dt, gr->GetY()[samp]);
    samp2++;
  }
  for(Int_t samp=0; samp<n/2; samp++){
    grUnwrapped->SetPoint(samp2, samp*dt, gr->GetY()[samp]);
    samp2++;
  }

  return grUnwrapped;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::getMaxMin(TGraph* gr, Double_t& maxY, Double_t& minY){

  Double_t maxX=0, minX=0;
  RootTools::getMaxMin(gr, maxY, maxX, minY, minX);
}




//---------------------------------------------------------------------------------------------------------
void RootTools::getMaxMin(TGraph* gr, Double_t& maxY, Double_t& maxX, 
              Double_t& minY, Double_t& minX){

  maxY=gr->GetY()[0];
  maxX=gr->GetX()[0];
  minY=gr->GetY()[0];
  minX=gr->GetX()[0];

  for(int i=0; i<gr->GetN(); i++){
    if(gr->GetY()[i] > maxY){
      maxY = gr->GetY()[i];
      maxX = gr->GetX()[i];
    }
    if(gr->GetY()[i] < minY){
      minY = gr->GetY()[i];
      minX = gr->GetX()[i];
    }
  } 
}





//---------------------------------------------------------------------------------------------------------
void RootTools::getMaxMinWithinLimits(TGraph* gr, Double_t& maxY, Double_t& maxX, 
                      Double_t& minY, Double_t& minX, 
                      Double_t lowerLimit, Double_t upperLimit){

  // Macro values from <cfloat> header
  Double_t minPoss = -DBL_MAX;
  Double_t maxPoss = DBL_MAX;
  maxY=minPoss;
  maxX=minPoss;
  minY=maxPoss;
  minX=maxPoss;

  for(int i=0; i<gr->GetN(); i++){
    if(gr->GetX()[i] >= lowerLimit && gr->GetX()[i] <=upperLimit){
      if(gr->GetY()[i] > maxY){
    maxY = gr->GetY()[i];
    maxX = gr->GetX()[i];
      }
      if(gr->GetY()[i] < minY){
    minY = gr->GetY()[i];
    minX = gr->GetX()[i];
      }
    }
  }

  if(maxY==minPoss){
    std::cerr << "Warning in " << __PRETTY_FUNCTION__ << "!" << std::endl;
    std::cerr << "Unable to find value in " << gr->GetName() << " > " << minPoss << std::endl;
  }
  if(minY==maxPoss){
    std::cerr << "Warning in " << __PRETTY_FUNCTION__ << "!" << std::endl;
    std::cerr << "Unable to find value in " << gr->GetName() << " < " << maxPoss << std::endl;
  }
}





//---------------------------------------------------------------------------------------------------------
void RootTools::subtractOffset(TGraph* gr, Double_t offset){
  for(int i=0; i<gr->GetN(); i++){
    gr->GetY()[i] -= offset;
  }
}





//---------------------------------------------------------------------------------------------------------
Double_t RootTools::getDeltaAngleDeg(Double_t angle1, Double_t angle2){
  
  Double_t deltaAngle = angle1 - angle2;
  Int_t loopCount = 0;
  const Int_t maxLoopCount = 5;
  if(deltaAngle > 180){
    while(deltaAngle >= 180){      
      deltaAngle -= 360;
      loopCount++;
      if(loopCount >= maxLoopCount){
    deltaAngle = -9999;
    break;
      }
    }
  }
  else{
    while(deltaAngle < -180){
      deltaAngle += 360;
      loopCount++;
      if(loopCount >= maxLoopCount){
    deltaAngle = -9999;
    break;
      }
    }
  }
  
  return deltaAngle;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::normalize(TGraph* gr){
  double mean, rms;
  normalize(gr, mean, rms);
}





//---------------------------------------------------------------------------------------------------------
void RootTools::normalize(TGraph* gr, Double_t& mean, Double_t& rms){
  RootTools::getMeanAndRms(gr, mean, rms);
  if(rms>0){ // Don't make any NaNs
    for(int i=0; i<gr->GetN(); i++){    
      gr->GetY()[i] -= mean;
      gr->GetY()[i] /= rms;
    }
  }
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeNormalized(TGraph* gr){
  /* Copies the TGraph and normalizes that */
  TGraph* grCopy = (TGraph*) gr->Clone();
  normalize(grCopy);
  return grCopy;
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeNormalized(TGraph* gr, Double_t& mean, Double_t& rms){
  /* Copies the TGraph and normalizes that */
  TGraph* grCopy = (TGraph*) gr->Clone();
  normalize(grCopy, mean, rms);
  return grCopy;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::getMeanAndRms(TGraph* gr, Double_t& mean, Double_t& rms){
  Double_t sum = 0;
  Double_t square = 0;
  for(int i=0; i<gr->GetN(); i++){
    sum += gr->GetY()[i];
    square += gr->GetY()[i]*gr->GetY()[i];
  }
  mean = sum/gr->GetN();
  rms = TMath::Sqrt(square/gr->GetN() - mean*mean);
}





//---------------------------------------------------------------------------------------------------------
Int_t RootTools::getIndexOfMaximum(Int_t len, Double_t* arr){
  Double_t max=-DBL_MAX;
  Int_t maxIndex = -1;
  for(Int_t i=0; i < len; ++i){
    if(arr[i] > max){
      max = arr[i];
      maxIndex = i;
    }
  }
  if(maxIndex==-1){
    std::cerr << "Warning in " << __PRETTY_FUNCTION__ << "!" << std::endl;
    std::cerr << "Could not find an index > " << max << std::endl;
  }
  
  return maxIndex;
}





//---------------------------------------------------------------------------------------------------------
std::vector<Int_t> RootTools::getIndicesOfNans(TGraph* gr){
  std::vector<Int_t> nanIndices;
  for(Int_t i=0; i<gr->GetN(); i++){
    if(TMath::IsNaN(gr->GetY()[i])){
      nanIndices.push_back(i);
    }
  }
  return nanIndices;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::printTGraphInfo(TGraph* gr){
  std::cout << "******************************************************************************" << std::endl;
  std::cout << "RootTools::printTGraphValues(TGraph* gr = " << gr << "):" << std::endl;
  std::cout << "Name: " << gr->GetName() << std::endl;
  std::cout << "Title: " << gr->GetTitle() << std::endl;
  std::cout << "Num points: " << gr->GetN() << std::endl;
  std::cout << "Xvals: ";
  for(Int_t i=0; i<gr->GetN(); i++){
    std::cout << gr->GetX()[i];
    if(i<gr->GetN()-1){
      std::cout << ", ";
    }
  }
  std::cout << std::endl;
  std::cout << "Yvals: ";
  for(Int_t i=0; i<gr->GetN(); i++){
    std::cout << gr->GetY()[i];
    if(i<gr->GetN()-1){
      std::cout << ", ";
    }
  }
  std::cout << std::endl;
  std::cout << "******************************************************************************" << std::endl;
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeSortedTGraph(TTree* tree, TString drawText, TString cutString, Double_t wrapValue){
  
  // Draw
  const Int_t nEntries = tree->Draw(drawText, cutString, "goff"); // "goff" means graphics off

  // Sort
  std::vector<Int_t> sortedIndices(nEntries);
  TMath::Sort(nEntries, tree->GetV2(), &sortedIndices.front(), kFALSE);
  std::vector<Double_t> newX(nEntries);
  std::vector<Double_t> newY(nEntries);
  
  for(int i=0; i<nEntries; i++){
    newX.at(i) = tree->GetV2()[sortedIndices.at(i)];
    newY.at(i) = tree->GetV1()[sortedIndices.at(i)];
  }

  // Unwrap data here so interpolation works smoothly
  // will unwrap when getting entries from graph
  if(wrapValue != 0){
    for(int i=1; i<nEntries; i++){
      // If y[i] >> y[i-1] => then we went below zero to the wrap value
      // can only modify y[i], so we should subtract wrapValue enough times 
      // that the graph becomes smooth
      while (newY.at(i) - newY.at(i-1) > wrapValue/2){
    newY.at(i) -= wrapValue;
      }

      // If y[i] << y[i-1] => then we went add zero to the wrap value
      // can only modify y[i], so we should add wrapValue enough times 
      // that the graph becomes smooth
      while (newY.at(i) - newY.at(i-1) < -wrapValue/2){
    newY.at(i) += wrapValue;
      }    

    }   
  }

  // sorted TGraph
  TGraph* gr(new TGraph(nEntries,&newX.front(), &newY.front()));
  TString title = cutString.Length() == 0 ? drawText : drawText + ", " + cutString;
  gr->SetTitle(title);
  TObjArray* tokens = drawText.Tokenize(":");
  gr->GetXaxis()->SetTitle(((TObjString*) tokens->At(1))->String());
  gr->GetYaxis()->SetTitle(((TObjString*) tokens->At(0))->String());
  delete tokens;
  
  return gr;
}





//---------------------------------------------------------------------------------------------------------
TH1D* RootTools::plotsZaxisDist(TH2* h2, TString hName, Int_t nBins, Double_t xMin, Double_t xMax){
  TH1D* h = new TH1D(hName, hName, nBins, xMin, xMax);
  for(int xBin=1; xBin<=h2->GetNbinsX(); xBin++){
    for(int yBin=1; yBin<=h2->GetNbinsY(); yBin++){
      Double_t val = h2->GetBinContent(xBin, yBin);
      h->Fill(val);
    }
  }
  return h;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::zeroPadTGraph(TGraph* gr, Int_t newLen, Double_t dt){
  Int_t oldLen = gr->GetN();
  gr->Set(newLen);
  if(dt >= 0){
    Double_t x0 = gr->GetX()[oldLen-1];
    for(Int_t samp=oldLen; samp<newLen; samp++){
      gr->GetX()[samp] = x0 + samp*dt;
    }
  }
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::makeLinearlyInterpolatedGraph(TGraph* grIn, Double_t dt){

  Int_t nIn = grIn->GetN();
  Int_t newPoints = Int_t((grIn->GetX()[nIn-1] - grIn->GetX()[0])/dt) + 1;
  std::vector<Double_t> newTimes(newPoints);
  std::vector<Double_t> newVolts(newPoints);

  Double_t time = grIn->GetX()[0];
  Double_t lastTime = grIn->GetX()[nIn-1];

  Int_t sampOut = 0; // Will iterate through 
  Int_t sampIn = 0;

  Double_t y1 = grIn->GetY()[sampIn];
  Double_t x1 = grIn->GetX()[sampIn];
  Double_t x2 = grIn->GetX()[sampIn+1];
  Double_t m = (grIn->GetY()[sampIn+1]-y1)/(grIn->GetX()[sampIn+1]-x1);
  
  while(time < lastTime){
    if(time > x2){

      Int_t countLoop = 0;
      while(time > x2){
    if(sampIn < (nIn - 1)){ // if last point is equal don't need to increase?
      sampIn++;
    }
    x2 = grIn->GetX()[sampIn+1];
    countLoop++;
    if(countLoop==1000){
      std::cerr << "You can't do logic very well" << std::endl;
      std::cerr << time << "\t" << x1 << "\t" << x2 << "\t" << nIn << "\t" << sampIn << std::endl;
      exit(0);
    }
      }
      x1 = grIn->GetX()[sampIn];
      y1 = grIn->GetY()[sampIn];

      m = (grIn->GetY()[sampIn+1]-y1)/(x2-x1);

    }

    // std::cout << time << "\t" << x1 << "\t" << x2;

    // if(time < x1 || time > x2 ){
    //   std::cout << "!!!!!!!!!!!!!!!!!!!!!!!";
    // }
    // std::cout << std::endl;

    Double_t newY = y1 + m * (time - x1);

    newTimes.at(sampOut) = time;
    newVolts.at(sampOut) = newY;
    sampOut++;
    time += dt;
  }

  TGraph* grOut = new TGraph(newTimes.size(), &newTimes[0], &newVolts[0]);
  // std::cout << grOut->GetN() - newPoints << " " << grIn->GetN() << " " << sampIn << " " 
  //        << grIn->GetX()[0] << " " << grIn->GetX()[nIn-1] << " " << dt << std::endl;
        



  return grOut;
}





//---------------------------------------------------------------------------------------------------------
TGraph* RootTools::interpolateWithStartTime(TGraph* grIn, Double_t startTime, Double_t dt, Int_t nSamp){  

  std::vector<Double_t> newTimes = std::vector<Double_t>(nSamp, 0);
  std::vector<Double_t> newVolts = std::vector<Double_t>(nSamp, 0);
  Double_t thisStartTime = grIn->GetX()[0];
  Double_t lastTime = grIn->GetX()[grIn->GetN()-1];


  // Quantizes the start and end times so data poInt_ts lie at Int_teger multiples of nominal sampling 
  // startTime = correlationDeltaT*TMath::Nint(startTime/correlationDeltaT + 0.5);
  // lastTime = correlationDeltaT*TMath::Nint(lastTime/correlationDeltaT - 0.5);
  startTime = dt*TMath::Nint(startTime/dt + 0.5);
  lastTime = dt*TMath::Nint(lastTime/dt - 0.5);  

   //ROOT Int_terpolator object constructor takes std::vector objects
  std::vector<Double_t> tVec(grIn->GetX(), grIn->GetX() + grIn->GetN());
  std::vector<Double_t> vVec(grIn->GetY(), grIn->GetY() + grIn->GetN());
  
  // This is ROOT's Int_terpolator object
  ROOT::Math::Interpolator chanInterp(tVec,vVec,ROOT::Math::Interpolation::kAKIMA);
  
  // Put new data Int_to arrays
  Double_t time = startTime;
  for(Int_t samp = 0; samp < nSamp; samp++){    
    newTimes.at(samp) = time;
    if(time >= thisStartTime && time <= lastTime){
      newVolts.at(samp) = chanInterp.Eval(time);
    }
    else{
      newVolts.at(samp) = 0;
    }
    time += dt;
  }
  return new TGraph(nSamp, &newTimes[0], &newVolts[0]);  

}




//---------------------------------------------------------------------------------------------------------
TCanvas* RootTools::drawArrayOfHistosPrettily(TH1D* hs[], Int_t numHists, TCanvas* can, Double_t* colWeights){

  if(can==NULL){
    can = new TCanvas();
  }
  can->cd();

  Double_t minCol = 0;
  Double_t maxCol = 255;
  if(colWeights!=NULL){
    maxCol = TMath::MaxElement(numHists, colWeights);
    minCol = TMath::MinElement(numHists, colWeights);
  }

  Int_t numDrawn = 0;
  Int_t first = 0;
  Double_t max = -1;
  for(Int_t histInd=0; histInd < numHists; histInd++){
    if(hs[histInd] != NULL){
      if(numDrawn==0){
    first = histInd;
      }
      TString opt = numDrawn == 0 ? "" : "same";
      if(colWeights!=NULL){
    hs[histInd]->SetLineColor(gStyle->GetColorPalette(Int_t(250.*(colWeights[histInd] - minCol)/(maxCol-minCol))));
    // std::cout << colWeights[histInd] << "\t" << maxCol << "\t" << minCol << "\t" 
    //    << ((colWeights[histInd] - minCol)/(maxCol-minCol)) << "\t" 
    //    << Int_t(250.*(colWeights[histInd] - minCol)/(maxCol-minCol))
    //    << std::endl;
      }
      else{
    hs[histInd]->SetLineColor(gStyle->GetColorPalette(histInd*Int_t(255./numHists)));
      }
      hs[histInd]->Draw(opt);
      if(hs[histInd]->GetMaximum() > max){
    max = hs[histInd]->GetMaximum()*1.1;
    hs[first]->SetMaximum(max);
      }
      numDrawn++;
    }
  }
  return can;
}




//---------------------------------------------------------------------------------------------------------
void RootTools::offsetTGraphXAxes(Int_t numGrs, TGraph* grs[], Double_t offsets[]){
  for(Int_t grInd=0; grInd < numGrs; grInd++){
    for(Int_t samp=0; samp < grs[grInd]->GetN(); samp++){
      grs[grInd]->GetX()[samp] += offsets[grInd];
    }
  }
}





//---------------------------------------------------------------------------------------------------------
void RootTools::multiplyTGraphYAxes(Int_t numGrs, TGraph* grs[], Double_t factors[]){
  for(Int_t grInd=0; grInd < numGrs; grInd++){
    for(Int_t samp=0; samp < grs[grInd]->GetN(); samp++){
      grs[grInd]->GetY()[samp] *= factors[grInd];
    }
  }
}





//---------------------------------------------------------------------------------------------------------
TCanvas* RootTools::drawArrayOfTGraphsPrettily(TGraph* grs[], Int_t numGrs, 
                           TString drawOpt, TCanvas* can, Double_t* colWeights){

  if(can==NULL){
    can = new TCanvas();
  }
  can->cd();

  Double_t colFactor = 1;
  if(colWeights!=NULL){
    Double_t maxCol = TMath::MaxElement(numGrs, colWeights);
    Double_t minCol = TMath::MinElement(numGrs, colWeights);
    colFactor = 254./(maxCol-minCol);
  }

  Double_t max = -DBL_MAX;
  Double_t min = DBL_MAX;
  Int_t numDrawn = 0;
  Int_t first = 0;
  for(Int_t grInd=0; grInd < numGrs; grInd++){
    if(grs[grInd]!=NULL){
      if(grs[grInd]->GetN() > 0){
    TString opt = numDrawn == 0 ? "a" + drawOpt : drawOpt + "same";
    if(numDrawn==0){
      first = grInd;
    }

    if(colWeights!=NULL){
      grs[grInd]->SetLineColor(gStyle->GetColorPalette(Color_t(colWeights[grInd]*colFactor)));
      grs[grInd]->SetMarkerColor(gStyle->GetColorPalette(Color_t(colWeights[grInd]*colFactor)));
      grs[grInd]->SetMarkerColor(gStyle->GetColorPalette(Color_t(colWeights[grInd]*colFactor)));
    }
    else{
      grs[grInd]->SetLineColor(gStyle->GetColorPalette(grInd*Int_t(254./(numGrs-1))));
      grs[grInd]->SetMarkerColor(gStyle->GetColorPalette(grInd*Int_t(254./(numGrs-1))));
      grs[grInd]->SetMarkerColor(gStyle->GetColorPalette(grInd*Int_t(254./(numGrs-1))));
    }

    grs[grInd]->Draw(opt);

    Double_t thisMax = TMath::MaxElement(grs[grInd]->GetN(), grs[grInd]->GetY());
    if(thisMax > max){
      Double_t fact = thisMax > 0 ? 1.1 : 0.9;
      max = thisMax*fact;
      grs[first]->SetMaximum(max);
    }

    Double_t thisMin = TMath::MinElement(grs[grInd]->GetN(), grs[grInd]->GetY());
    if(thisMin < min){
      Double_t fact = thisMin < 0 ? 1.1 : 0.9;
      min = thisMin*fact;
      grs[first]->SetMinimum(min);
    }
    numDrawn++;
      }
    }
  }

  return can;

}





//---------------------------------------------------------------------------------------------------------
TLegend* RootTools::makeLegend(TGraph* grs[], Int_t numGrs, 
                   TString titles[], TString opt, 
                   Double_t minX, Double_t minY,
                   Double_t maxX, Double_t maxY){

  TLegend* l = new TLegend(minX, minY, maxX, maxY);
  for(Int_t grInd=0; grInd < numGrs; grInd++){
    if(grs[grInd]!=NULL){
      l->AddEntry(grs[grInd], titles[grInd], opt);
    }
  }
  return l;
}





//---------------------------------------------------------------------------------------------------------
TLegend* RootTools::makeLegend(TH1D* hs[], Int_t numHists, 
                   TString titles[], TString opt, 
                   Double_t minX, Double_t minY,
                   Double_t maxX, Double_t maxY){

  TLegend* l = new TLegend(minX, minY, maxX, maxY);
  for(Int_t histInd=0; histInd < numHists; histInd++){
    if(hs[histInd]!=NULL){
      l->AddEntry(hs[histInd], titles[histInd], opt);
    }
  }
  return l;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::writeTGraph(TGraph* gr, TString name){
  // Taking laziness to a new level... replace two lines with a one line function.
  gr->SetName(name);
  gr->Write();
}





//---------------------------------------------------------------------------------------------------------
void RootTools::getLocalMaxToMin(TGraph* gr, 
                 Double_t& maxY, Double_t& maxX, 
                 Double_t& minY, Double_t& minX){
  getLocalMaxToMinWithinLimits(gr, maxY, maxX, minY, minX, gr->GetX()[0], gr->GetX()[gr->GetN()-1]);

}





//---------------------------------------------------------------------------------------------------------
void RootTools::getLocalMaxToMinWithinLimits(TGraph* gr, 
                         Double_t& maxY, Double_t& maxX, 
                         Double_t& minY, Double_t& minX,
                         Double_t lowerLimit, Double_t upperLimit){


  // assumes unique value at local minima and maxima, e.g. 0.1, 0.2, 0.1 for maxima.
  // if it goes 0.1, 0.2, 0.2, 0.1 then this will need to get more complicated
  std::vector<Int_t> extremaSamps;

  for(Int_t samp=1; samp < gr->GetN()-1; samp++){
    if(gr->GetX()[samp] >= lowerLimit && gr->GetX()[samp] <=upperLimit){
      Double_t y0 = gr->GetY()[samp-1];
      Double_t y1 = gr->GetY()[samp];
      Double_t y2 = gr->GetY()[samp+1];
    

      if(y0 > y1 && y1 < y2){
    // Is a local minimum
    extremaSamps.push_back(samp);
      }
      else if(y0 < y1 && y1 > y2){
    // Is a local maximum
    extremaSamps.push_back(samp);
      }
    }
  }
  
  Double_t maxOffset = -DBL_MAX;
  Int_t maxSampInd = -1;
  for(UInt_t sampInd = 0; sampInd < extremaSamps.size()-1; sampInd++){
    Double_t offset = TMath::Abs(gr->GetY()[extremaSamps.at(sampInd)] - gr->GetY()[extremaSamps.at(sampInd+1)]);
    if(offset > maxOffset){
      maxOffset = offset;
      maxSampInd = sampInd;
    }
  }

  if(gr->GetY()[extremaSamps.at(maxSampInd)] > gr->GetY()[extremaSamps.at(maxSampInd+1)]){
    maxY = gr->GetY()[extremaSamps.at(maxSampInd)];
    maxX = gr->GetX()[extremaSamps.at(maxSampInd)];

    minY = gr->GetY()[extremaSamps.at(maxSampInd+1)];
    minX = gr->GetX()[extremaSamps.at(maxSampInd+1)];
  }
  else{
    maxY = gr->GetY()[extremaSamps.at(maxSampInd+1)];
    maxX = gr->GetX()[extremaSamps.at(maxSampInd+1)];

    minY = gr->GetY()[extremaSamps.at(maxSampInd)];
    minX = gr->GetX()[extremaSamps.at(maxSampInd)];
  }
}



//---------------------------------------------------------------------------------------------------------
void RootTools::saveCanvas(TCanvas* c, TString fileName){

  std::cout << "Saving this canvas as an .eps, .png, and .C file..." << std::endl;
  TString fName = fileName + ".eps";
  c->SaveAs(fName);
  fName = fileName + ".png";
  c->SaveAs(fName);
  fName = fileName + ".C";
  c->SaveAs(fName);  
  std::cout << "...Complete!" << std::endl;
  c->Update();
}





//---------------------------------------------------------------------------------------------------------
Double_t RootTools::getFullWidthHalfMax(TH1D* h){
  
  Int_t maxBin = RootTools::getPeakBinOfHistogram(h);
  Int_t n = h->GetNbinsX();
  Double_t max = h->GetBinContent(maxBin);


  Double_t halfMaxPos = 0;
  for(Int_t bin=maxBin; bin<=n; bin++){
    Double_t val = h->GetBinContent(bin);
    if(val < max/2){
      halfMaxPos = h->GetBinLowEdge(bin);
      break;
    }
  }

  Double_t halfMaxNeg = 0;
  for(Int_t bin=maxBin; bin>=1; bin--){
    Double_t val = h->GetBinContent(bin);
    if(val < max/2){
      halfMaxNeg = h->GetBinLowEdge(bin);
      break;
    }
  }

  return halfMaxPos - halfMaxNeg;
}





//---------------------------------------------------------------------------------------------------------
Int_t RootTools::getPeakBinOfHistogram(TH1D* h){

  Int_t n = h->GetNbinsX();
  Double_t max = -DBL_MAX;
  Int_t maxBin = 0;
  for(Int_t bin=1; bin<n; bin++){
    Double_t val = h->GetBinContent(bin);
    if(val > max){
      max = val;
      maxBin = bin;
    }
  }
  return maxBin;
}





//---------------------------------------------------------------------------------------------------------
Double_t RootTools::getPeakBinOfHistogram(TH2D* hist, Int_t& binx, Int_t& biny){

  Int_t nx = hist->GetNbinsX();
  Int_t ny = hist->GetNbinsY();  
  Double_t histPeak = -DBL_MAX;
  for(Int_t by = 1; by<=ny; by++){
    for(Int_t bx = 1; bx<=nx; bx++){
      Double_t val = hist->GetBinContent(bx, by);
      if(val > histPeak){
    histPeak = val;
    binx = bx;
    biny = by;
      }
    }
  }

  return histPeak;
  
}





//---------------------------------------------------------------------------------------------------------
Double_t RootTools::getLowBinEdgeOfHistogramPeak(TH1D* h){

  Int_t peakBin = getPeakBinOfHistogram(h);
  Double_t lowBinEdgeOfPeak = h->GetXaxis()->GetBinLowEdge(peakBin);
  return lowBinEdgeOfPeak;
}






//---------------------------------------------------------------------------------------------------------
void RootTools::makeZaxisScaleEqualAboutZero(TH2D* h){

  Double_t min = h->GetMinimum();
  Double_t max = h->GetMaximum();

  if(TMath::Abs(min) > max){
    h->SetMaximum(-min);
  }
  else{
    h->SetMinimum(-max);
  }
}





//---------------------------------------------------------------------------------------------------------
Int_t RootTools::getBit(UInt_t bitIndex, UInt_t bitMask){
  return ((bitMask >> bitIndex) & 1);
}





//---------------------------------------------------------------------------------------------------------
Int_t RootTools::getNumBitsSet(Int_t numBitsToCheck, UInt_t bitMask){
  Int_t numBitsHigh = 0;
  for(Int_t bitInd=0; bitInd < numBitsToCheck; bitInd++){
    if(getBit(bitInd, bitMask) > 0){
      numBitsHigh++;
    }    
  }
  return numBitsHigh;
}





//---------------------------------------------------------------------------------------------------------
void RootTools::setWhiteZeroColorScale(){
  const int NRGBs = 3, NCont = 999;
  gStyle->SetNumberContours(NCont);
  Double_t stops[NRGBs] = { 0.00, 0.50, 1.00};
  Double_t red[NRGBs]   = { 0.00, 1.00, 1.00};
  Double_t green[NRGBs] = { 0.00, 1.00, 0.00};
  Double_t blue[NRGBs]  = { 1.00, 1.00, 0.00};
  TColor color;
  color.InitializeColors();
  color.CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); 
}






//---------------------------------------------------------------------------------------------------------
void RootTools::draw2D(TH2D* hist, TString opt){
  hist->Draw(opt);
  Double_t max = TMath::Abs(hist->GetMaximum());
  Double_t min = TMath::Abs(hist->GetMinimum());
  Double_t limit = min > max ? min : max;
  hist->SetMaximum(limit);
  hist->SetMinimum(-limit);
}





//---------------------------------------------------------------------------------------------------------
TChain* RootTools::getHeadChain(Int_t firstRun, Int_t lastRun, RawAnitaHeader*& headPtr){
  TChain* c = new TChain("headTree");
  for(Int_t run=firstRun; run <= lastRun; run++){
    c->Add(TString::Format("~/UCL/ANITA/flight1415/root/run%d/headFile%d.root", run, run));
  }
  c->SetBranchAddress("header", &headPtr);
  return c;
}





//---------------------------------------------------------------------------------------------------------
TChain* RootTools::getAdu5PatChain(Int_t firstRun, Int_t lastRun, Adu5Pat*& pat){
  TChain* c = new TChain("adu5PatTree");
  for(Int_t run=firstRun; run <= lastRun; run++){
    c->Add(TString::Format("~/UCL/ANITA/flight1415/root/run%d/gpsFile%d.root", run, run));
  }
  c->SetBranchAddress("pat", &pat);
  return c;
}






//---------------------------------------------------------------------------------------------------------
TString RootTools::getAntName(AnitaPol::AnitaPol_t pol, Int_t antInd){
  // Assumes antInd starts at 0. Returns the name w/ +1
  Int_t phi = (antInd % NUM_PHI) + 1;
  
  Int_t ring = antInd / NUM_PHI;

  TString antName = TString::Format("%d", phi);
  
  antName += AnitaRing::ringAsChar(AnitaRing::AnitaRing_t(ring));
  antName += AnitaPol::polAsChar(pol);  
  
  return antName;
}






//---------------------------------------------------------------------------------------------------------
Int_t RootTools::getColorFracThroughPalette(Int_t index, Int_t maxVal){
  return gStyle->GetColorPalette(index*Int_t(255./maxVal));  
}






//---------------------------------------------------------------------------------------------------------
TCanvas* RootTools::drawHistsWithStatsBoxes(Int_t numHists, TH1D* hs[], TString drawOpt, TString statsOption){

  gStyle->SetOptStat(statsOption);
  TCanvas* theCan = new TCanvas();

  Double_t boxSize = 0.2; // In normalized coordinates
  const Int_t maxNumBoxes = 1./boxSize;
  
  Int_t divisor = (numHists/maxNumBoxes) + 1; 

  // This is dumb and will get messy If I have more than 2*maxnumBoxes
  boxSize/=divisor;
  
  for(Int_t histInd=0; histInd < numHists; histInd++){

    TString thisDrawOpt = histInd == 0 ? drawOpt : drawOpt+"sames";
    hs[histInd]->Draw(thisDrawOpt);

    theCan->Update();
    TPaveStats* sbox = (TPaveStats*) hs[histInd]->FindObject("stats");
    sbox->SetTextColor(hs[histInd]->GetLineColor());
    sbox->SetX1NDC(0.8);
    sbox->SetY1NDC(1 - (histInd+1)*boxSize);
    sbox->SetX2NDC(0.98);
    sbox->SetY2NDC(0.98 - histInd*boxSize);
  }

  theCan->Update();
  return theCan; 
}