---

rootEnergyHist.cxx

Go to the documentation of this file.

//  class rootEnergyHist
//  Author:  Theodore Hierath
//  This class contains a histogram and graphing class for
//  energies.
//  The bins are numbered from 0 to number of bins - 1;

#include "rootEnergyHist.h"
#include <iostream>
#include <iomanip>
#include <fstream>

rootEnergyHist::rootEnergyHist(int bins, 
                               double min_energy, 
                               double max_energy) :
                               emin(min_energy),
                               emax(max_energy),
                               num_bins(bins)
{
   range = log10(emax/emin);
   currentType = linear;
   flux_hist = new rootHist(bins);
   flux_sigma_hist = new rootHist(bins);
   raw_hist = new rootHist(bins);
   graphTitle = "No Title";
   xlabel = "";
   ylabel = "";
   fluxMode = false;
   use_flux_min = false;
   use_flux_max = false;
}

rootEnergyHist::~rootEnergyHist(void)
{
   delete flux_hist;
   delete flux_sigma_hist;
   delete raw_hist;
}

rootEnergyHist::rootEnergyHist(const rootEnergyHist& oldHist) :
                               emin(oldHist.emin),
                               emax(oldHist.emax),
                               num_bins(oldHist.num_bins)
{
   range = oldHist.range;
   currentType = oldHist.currentType;
   flux_hist = new rootHist(num_bins);
   flux_sigma_hist = new rootHist(num_bins);
   raw_hist = new rootHist(num_bins);

   for(int i = 0; i < num_bins; i++)
   {  
      double oldCount = (oldHist.raw_hist)->retrieveBin(i);
      raw_hist->updateBin(i, oldCount );

      oldCount = (oldHist.flux_hist)->retrieveBin(i);
      flux_hist->updateBin(i, oldCount );

      oldCount = (oldHist.flux_sigma_hist)->retrieveBin(i);
      flux_sigma_hist->updateBin(i, oldCount );
   }

   graphTitle = oldHist.graphTitle;
   xlabel = oldHist.xlabel;
   ylabel = oldHist.ylabel;
   fluxMode = oldHist.fluxMode;
   use_flux_min = oldHist.use_flux_min;
   use_flux_max = oldHist.use_flux_max;
}


void rootEnergyHist::setTitle(std::string title)
{
   graphTitle = title;
}

void rootEnergyHist::setXLabel(std::string label)
{
   xlabel = label;
}

void rootEnergyHist::setYLabel(std::string label)
{
   ylabel = label;
}

void rootEnergyHist::setGraphType(const char *graph_type)
{
   if(0 == strcmp(graph_type,"linear"))
      currentType = linear;
   else if(0 == strcmp(graph_type,"semilogx"))
      currentType = semilogx;
   else if(0 == strcmp(graph_type,"semilogy"))
      currentType = semilogy;
   else if(0 == strcmp(graph_type,"log"))
      currentType = loglog;
   else
      std::cerr << "ERROR:  Invalid Graph Type" << std::endl;
}

void rootEnergyHist::store(double energy)
{ 
   if(energy < emax && energy > emin)
   {
      if(currentType == loglog || currentType == semilogx)
      {
         int currentBin = int(num_bins * (log10(energy/emin) / range));
         double contents = raw_hist->retrieveBin(currentBin);
         raw_hist->updateBin(currentBin,++contents);
      }
      else
      {
         int currentBin = int(floor(num_bins * (energy - emin) / (emax - emin)));
         double contents = raw_hist->retrieveBin(currentBin);
         raw_hist->updateBin(currentBin,++contents);
      }
   }
}

double rootEnergyHist::retrieveCount(int binNumber)
{
   return raw_hist->retrieveBin(binNumber);
}

double rootEnergyHist::retrieveFlux(int binNumber)
{
   return flux_hist->retrieveBin(binNumber);
}

double rootEnergyHist::retrieveFluxUncertainty(int binNumber)
{
   return flux_sigma_hist->retrieveBin(binNumber);
}

double rootEnergyHist::retrieveEnergy(int binNumber)
{
   return emin*pow(10,((binNumber + 0.5)/num_bins) * range );
}

double rootEnergyHist::retrieveRange(void)
{
   return range;
}


void rootEnergyHist::setFluxMode(void)
{
   fluxMode = true;
}

void rootEnergyHist::setFluxMin(double f)
{
   use_flux_min = true;
   flux_min = f;
}

void rootEnergyHist::setFluxMax(double f)
{
   use_flux_max = true;
   flux_max = f;
}

// This function must be called before the drawing function
void rootEnergyHist::apply(double scale_factor)
{
   for(int i = 0; i < num_bins; ++i)
   {
      double current_flux = flux_hist->retrieveBin(i);
      double current_count = raw_hist->retrieveBin(i);
      current_flux += (current_count * scale_factor);
      flux_hist->updateBin(i,current_flux);

      // calculate uncertainty by addition under quadrature
      current_flux = sqrt(pow(flux_sigma_hist->retrieveBin(i),2)+current_count) * scale_factor;
      flux_sigma_hist->updateBin(i,current_flux);
      
      // clear temporary array
      raw_hist->updateBin(i,0);
   }
}

void rootEnergyHist::reset(void)
{
   for(int i = 0; i < num_bins; ++i)
   {
      flux_hist->updateBin(i,0);
      flux_sigma_hist->updateBin(i,0);
      raw_hist->updateBin(i,0);
   }

   graphTitle = "No Title";
   xlabel = "";
   ylabel = "";
   fluxMode = false;
   use_flux_min = false;
   use_flux_max = false;
   range = log10(emax/emin);
   currentType = linear;
}


// scale_factor = number to multiply each bin by
// mode = "normal", "begin", or "end"
// current_plot = current_plot number (numbering starts at 0)
// total_plots = the total number of plots to go on the graph
void rootEnergyHist::draw(double scale_factor, std::string mode, int current_plot, int total_plots)
{
   char *window_title = "Graph Window";

   if(current_plot >= total_plots)
   {
      std::cerr << "Error:  Invalid plot number" << std::endl;
      exit(0);
   }

   if(current_plot == 0)
   {
      std::ofstream out_file;

      if(mode != "end")
      {
         out_file.open("graph.cxx");

         if(false == out_file.is_open())
         {
            std::cerr << "Unable to open temporary file for writing." << std::endl;
            exit(0);
         }

         out_file << 
            "{\n"
            "   gROOT->Reset();\n";
      }
      else
      {
         out_file.open("graph.cxx", std::ios::app);   
      }

      out_file <<         
         "   char *energy_window_title = \"" << window_title << "\";\n"
         "   char *energy_graph_title = \"" << "Particle Flux vs. Kinetic Energy" << "\";\n"
         "   char *energy_y_label = \"" << ylabel << "\";\n"
         "   char *energy_x_label = \"" << xlabel << "\";\n"
         
         "   double energy_min= "<< emin << ", energy_max=" << emax << ";\n"
         "   int num_bins = "<< num_bins << ";\n"
         "   double log_energy_range = " << retrieveRange() << ";\n"
         "   double energy[num_bins];\n"
         "   double e_energy[num_bins];\n"
         
         "   for(int i = 0; i < num_bins; i++) {\n"
         "      energy[i] = energy_min*pow(10,((i+0.5)/num_bins) * log_energy_range);\n"
         "      e_energy[i] = 0;\n"
         "   }\n"
         "   //                 name, title, wtopx, wtopy, ww, wh\n"
         "   c1 = new TCanvas(\"c1\",energy_window_title, 200, 200, 700, 500);\n"
         "   c1->SetGrid();\n"
         "   c1->GetFrame()->SetFillColor(21);\n"
         "   c1->GetFrame()->SetBorderSize(12);\n";
           
      if(currentType == loglog)        out_file << "   c1->SetLogx(1); c1->SetLogy(1);\n";
      else if(currentType == semilogx) out_file << "   c1->SetLogx(1);\n";
      else if(currentType == semilogy) out_file << "   c1->SetLogy(1);\n";
      else if(currentType == linear)   ;
      else
         std::cerr << "Error:  invalid graph type" << std::endl;
      
      out_file << 
         "   leg = new TLegend(0.73,0.83,0.99,0.99);\n"
         "   double scale_factor" << current_plot << " = " << scale_factor << ";\n"
         "   double count" << current_plot << "[] = {\n";
      {for(int i = 0; i < num_bins; ++i)
         out_file << std::setw(7) << retrieveFlux(i) << (i%5==4? ",\n" : ",");
      }

      out_file << 
         "   };\n"
         "   double e_count" << current_plot << "[num_bins] = {\n";
      {for(int i = 0; i < num_bins; ++i)
         out_file << std::setw(7) << retrieveFluxUncertainty(i) << (i%5==4? ",\n" : ",");
      }

      out_file <<
         "   };\n"
         "   for(i = 0; i < num_bins; i++) {\n"
         "      e_count" << current_plot << "[i] *= scale_factor" << current_plot;
      if(fluxMode == true)
         out_file << "/(energy[i]*1000);\n";
      else
         out_file << ";\n";

      out_file <<
         "      count" << current_plot << "[i] *=scale_factor" << current_plot;

      if(fluxMode == true)
         out_file << "/(energy[i]*1000);\n";
      else
         out_file << ";\n";

      out_file <<
         "   }\n"
         "   gr" << current_plot << " = new TGraphErrors(num_bins,energy,count" << current_plot 
                 << ",e_energy,e_count" << current_plot << ");\n"
         "   gr" << current_plot << "->SetTitle(energy_graph_title);\n"
         "   gr" << current_plot << "->SetMarkerColor(2);\n"
         "   gr" << current_plot << "->SetMarkerStyle(21);\n"
         "   leg->AddEntry(gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n";
     
      out_file.close();
   }   
   else if(current_plot <= total_plots - 1)
   {
      std::ofstream out_file("graph.cxx", std::ios::app);
      out_file <<
         "   c1->cd();\n"
         "   double scale_factor" << current_plot << " = " << scale_factor << ";\n"
         "   double count" << current_plot << "[] = {\n";
      {for(int i = 0; i < num_bins; i++)
         out_file << std::setw(5) << retrieveFlux(i) << (i%5==4? ",\n" : ",");
      }

      out_file << 
         "   };\n"
         "   double e_count" << current_plot << "[num_bins] = {\n";
      {for(int i = 0; i < num_bins; ++i)
         out_file << std::setw(7) << retrieveFluxUncertainty(i) << (i%5==4? ",\n" : ",");
      }

      out_file <<
         "   };\n"
         "   for(i = 0; i < num_bins; i++) {\n"
         "      e_count" << current_plot << "[i] *= scale_factor" << current_plot;
      if(fluxMode == true)
         out_file << "/(energy[i]*1000);\n";
      else
         out_file << ";\n";

      out_file <<
         "      count" << current_plot << "[i] *=scale_factor" << current_plot;

      if(fluxMode == true)
         out_file << "/(energy[i]*1000);\n";
      else
         out_file << ";\n";

      out_file <<
         "   }\n"
         "   gr" << current_plot << " = new TGraphErrors(num_bins,energy,count" << current_plot 
                 << ",e_energy,e_count" << current_plot << ");\n"
         "   gr" << current_plot << "->SetMarkerColor(" << 2+current_plot << ");\n"
         "   gr" << current_plot << "->SetMarkerStyle(" << 21 << ");\n"
         "   leg->AddEntry(gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n";
      out_file.close();
   }
   
   if(current_plot >= total_plots - 1)
   {
      std::ofstream out_file("graph.cxx", std::ios::app);

      out_file << 
         "   double max_count = 0;\n"
         "   double min_count = 1e12;\n"
         "   for(int i = 0; i < num_bins; i++)\n"
         "   {\n";
      {for(int plot = 0; plot < total_plots; plot++) {
         out_file << 
            "      if(count" << plot << "[i] > max_count)\n"
            "         max_count = count" << plot << "[i];\n"
            "      if(count" << plot << "[i] < min_count && count" << plot << "[i] > 0)\n"
            "         min_count = count" << plot << "[i];\n";
      }}
      out_file <<
         "   }\n"
         "   double energy_limits[] = {energy_min,energy_max};\n"
         "   double count_limits[] = {";

      if(use_flux_min)
         out_file << flux_min << ",";
      else
         out_file << "min_count,";

      if(use_flux_max)
         out_file << flux_max;
      else
         out_file << "max_count";

      out_file <<
         "};\n"
         "   int bin_limits = 2;\n"
         "   graph0 = new TGraph(bin_limits,energy_limits,count_limits);\n"
         "   graph0->SetTitle(energy_graph_title);\n"
         "   graph0->Draw(\"AP\");\n"
         "   TAxis *ax = graph0->GetXaxis();\n"
         "   TAxis *ay = graph0->GetYaxis();\n"
         "   ay->SetTitle(energy_y_label); ay->CenterTitle(1);\n"
         "   ax->SetLimits(energy_min, energy_max);\n"
         "   ax->SetTitle(energy_x_label); ax->CenterTitle(1); \n"
         "   ax->SetTitleOffset(1.2);\n";

      {for(int plot = 0; plot < total_plots; plot++) {
         out_file <<
            "   gr" << plot << "->Draw(\"P\");\n";
      }}

      out_file << 
         "   leg->Draw();\n"
         "   c1->Modified();\n"
         "   c1->Update();\n";

      if(mode != "begin")
         out_file << "\n}\n";

      out_file.close();

      if(mode != "begin")
         system("root -l graph.cxx");
   }
}


// If the output is going to be used for another source, the histogram should have 
// a large enough number of bins to minimize rounding errors caused by bin size.
void rootEnergyHist::writeFile(double scale_factor, std::ostream& out_file)
{
   out_file << "%Energy (in MeV)  vs. Flux (in particles/m2-s-sr-MeV)" << std::endl;

   for(int i = 0; i < num_bins; i++)
   {
      // Use center of bin to determine the energies
      double energy = retrieveEnergy(i) * 1000; // GeV->MeV
      double flux = double(retrieveFlux(i)) * scale_factor / energy; // E*Flux -> Flux

      out_file << std::setiosflags(std::ios::right|std::ios::scientific) 
               << std::setw(14) << std::setprecision(5) << energy 
               << std::setiosflags(std::ios::right|std::ios::scientific)
               << std::setw(14) << std::setprecision(5) << flux 
               << std::endl;
   }

}

---