/*ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
|                                                                       |
|  HOT CELL FINDER: analysis code                                       |
|                                                                       |
|  setup D0RunII p17.03.03                                              |
|  gmake                                                                |
|  ./FindHotCells                                                       |
|                                                                       |
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo*/

#define hotcells_cxx
#include "hotcells.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH2.h"
#include "TPad.h"
#include "TText.h"
#include <iostream>
#include <fstream>

using namespace std;

void hotcells::finder(float nsigma) {
 
  //----------------------------------------------------------------
  // OUTPUT FILE
  //----------------------------------------------------------------
 
  // Number of runs in range
  int const num_runs = last_run - first_run + 1;

  // Output file
  char outputfile[100];
  sprintf (outputfile, "HotCell-Analysis-%i-%i-sigma-%g.txt", 
	   first_run, last_run, nsigma);
  ofstream Output;
  Output.open(outputfile);

  Output << "\n---------------------------------" << endl;
  Output << "     Hot Cell Finder               " << endl;
  Output << "---------------------------------\n" << endl;

  Output << "Filelist: " << filelist << endl;

  for (int f=0; f!=filename.size(); ++f) {
    Output << filename[f] << endl;
  }

  Output << "- " << filename.size() << " files" << endl;

  Output << "\nRunlist: " << runlist << endl;

  Output << "- " << goodrun_list.size() + shortrun_list.size() 
	 << " runs of which " << goodrun_list.size() << " long enough" << endl;
  
  Output << "\nRun range: " << first_run 
	 << " - " << last_run << endl;
  Output << "Sigma threshold: " << nsigma << endl;
  Output << "Min number of events for run: " << min_events << endl;
  
  //----------------------------------------------------------------
  // OPTIONS
  //----------------------------------------------------------------
  
  // Do hot cell finding
  bool hot_cell_analysis = true;
  if (hot_cell_analysis) Output << "\nDo hot cell analysis" << endl; 

  // Make hot cell histograms
  bool hot_cell_plots = false;
  
  if (hot_cell_plots) {
    hot_cell_analysis = true;
    Output << "and hot cell histograms" << endl;
  }

  // Draw options
  bool logz = true; 
  if (logz) Output << "Draw plots with log z axis" << endl;

  // Write histogram(s) to files
  bool write = false;

  // Write only one histogram
  bool plot_one = true;
  int plot_run = 173520;
  int plot_layer = 1; // 1-17
  if (plot_one) write = true;

  if (write) {
    hot_cell_analysis = hot_cell_plots = true;
    Output << "Write plot(s) to file" << endl;
  }
  
  if (plot_one) Output << "Only Run " << plot_run << ", layer " << plot_layer << endl;
  
  //----------------------------------------------------------------
  // OUTPUT SHORT RUNS
  //----------------------------------------------------------------

  Output << "\nShort runs (not used):" << endl;

  for (map<int,int>::const_iterator it=shortrun_list.begin();
       it!=shortrun_list.end();
       ++it) {
    Output << it->first << " ";
  }
  Output << endl;

  //----------------------------------------------------------------
  // GET HISTOGRAMS
  //----------------------------------------------------------------

  // Number of calorimeter layers
  int const cal_layers = 17;
  
  TH2F *histo[num_runs][cal_layers];
  char plot_name[num_runs][cal_layers][100];

  // Mark and count existing histograms
  bool found_histo[num_runs][cal_layers];
  bool found_histo_hot[num_runs][cal_layers];
  int num_histos[num_runs];
  //bool found_run[num_runs];
  
  Output << "\nGetting cell energy histograms" << endl;

  vector<TH2F*> histo_group;
  int first_goodrun;

  // RUN LOOP: iterate through map of good runs
  for (map<int,int>::const_iterator it=goodrun_list.begin(); 
       it!=goodrun_list.end(); 
       ++it) {

    // Run number is key in map    
    int runnumber = it->first;

    // Store first run number
    if (it == goodrun_list.begin()) first_goodrun = runnumber;

    int r = runnumber - first_goodrun;

    num_histos[r] = 0;

    // LAYER LOOP
    for (int l=0; l!=cal_layers; ++l) {

      found_histo[r][l] = false;
      histo_group.clear();
      
      // Histogram name
      sprintf(plot_name[r][l],"_Cell_E_layer_%i_ieta_iphi_Run-%i",l+1,runnumber);

      // FILE LOOP
      for (int f=0; f!=filename.size(); ++f) {

	// Look for histogram
	TH2F* histo_tmp = (TH2F*)ZBfile[f]->Get(plot_name[r][l]);

	// Add to group
	if (histo_tmp) {
	  histo_tmp->Sumw2();
	  histo_group.push_back(histo_tmp);
	} 

      } // file loop

      // If histo found in any file
      if (histo_group.size()) {

	// Add histos from different files
	for (int h=0; h!=histo_group.size(); ++h) {

	  if (h == 0) histo[r][l] = (TH2F*)histo_group[0]->Clone();
	  else histo[r][l]->Add(histo_group[h],1);
	  
	}

	// Keep track of which histograms exist
	found_histo[r][l] = true;
      
	// Count number of histograms
	// for eps file closure
	++num_histos[r];
      }

    } // layers

  } // runs


  //----------------------------------------------------------------
  // CREATE HOT CELL HISTOGRAMS
  //----------------------------------------------------------------

  // Histograms showing hot cells
  TH2F *histo_hot[num_runs][cal_layers];

  // Create hot cell plots
  if (hot_cell_plots) {

    Output << "\nCreating hot cell histograms" << endl;

    char name[100];
    char title[100];
    
    // RUN LOOP
    for (map<int,int>::const_iterator it=goodrun_list.begin();
	 it!=goodrun_list.end();
	 ++it) {

      int runnumber = it->first;
      int r = runnumber - first_goodrun;
      
      for (int l=0; l!=cal_layers; ++l) {
	
	// Check if histo exists
	if (!found_histo[r][l]) continue;
	
	// Create histogram to show hot cells
	sprintf (name,"_Cell_E_layer_%i_ieta_iphi_hot_Run-%i",
		 l+1, runnumber);
	sprintf (title, "Hot Cells layer %i: Run %i",
		 l+1, runnumber);
	
	histo_hot[r][l] = 
	  new TH2F(name, title, 75,-37.5,37.5,64,0.5,64.5);
	
      }
    }
  } // if hot cell plots    

  //----------------------------------------------------------------
  // HOT CELL FILE
  //----------------------------------------------------------------
 
  // Hot cell file
  char hotcellsfile[100]; 
  ofstream HotCells;
  sprintf (hotcellsfile, "HotCell-List-%i-%i-sigma-%g.txt", first_run, last_run, nsigma);
  HotCells.open(hotcellsfile);


  //----------------------------------------------------------------
  // HOT CELL ANALYSIS
  //----------------------------------------------------------------

  // Number of hot cells found per run, total
  int HotCell[num_runs];
  int num_hot_cells = 0;

  // Calorimeter geometry (num phi varies)
  int const num_eta = 75;
  int const num_phi = 64;
  int num_cells_phi;

  // Values for s.d. calculation
  float etaring[cal_layers][num_eta];
  float bin[num_phi], deviation;
  float sum_diff_squared, sigma, threshold;
  
  if (hot_cell_analysis) {

    Output << "\nStarting hot cell analysis" << endl;

    // 1. LOOP OVER RUNS
    for (map<int,int>::const_iterator it=goodrun_list.begin();
	 it!=goodrun_list.end();
	 ++it) {

      int runnumber = it->first;
      int r = runnumber - first_goodrun;
    
      HotCell[r] = 0;
      
      // 2. LOOP OVER LAYERS
      for (int l=0; l!=cal_layers; ++l) {
	
	// Check if cell energy histogram exists
	if (!found_histo[r][l]) continue;
		
	// 3. LOOP OVER IETA
	for (int n=0; n!=num_eta; ++n) {
	  
	  int ieta = n-37;
	  bool forward = false;
	  etaring[l][n] = 0;

	  // There's no ieta = 0
	  if (ieta == 0) continue;
	  
	  // Forward calorimeter has twice the cell size in phi
	  if (abs(ieta) >= 33) forward = true;
	  if (forward) num_cells_phi = 32;
	  else num_cells_phi = 64;
	  
	  // 4. SUM ALL CELL ENERGIES IN THAT IETA
	  for (int p=0; p!=num_phi; ++p) {
	    
	    int iphi = p+1;
	    bin[p] = 0.0;
	    
	    // Forward region: phi is odd only
	    if (forward && (iphi % 2 == 0)) continue;
	    
	    // Get cell energy
	    bin[p] = histo[r][l]->GetBinContent(n+1,p+1);

	    // Add to sum	    
	    etaring[l][n] += bin[p];
	    
	  }

	  // 4. CALCULATE MEAN CELL ENERGY 	  
	  float mean = etaring[l][n]/num_cells_phi;
	  sum_diff_squared = 0;
	  
	  // 5. CALCULATE EACH CELL'S DEVIATION FROM THE MEAN
	  for (int p=0; p!=num_phi; ++p) {
	    
	    int iphi = p+1;
	    deviation = 0.0;
	    // Forward region: phi is odd only
	    if (forward && (iphi % 2 == 0)) continue;

	    // Deviation is (cell_E - mean_E)	
	    deviation = bin[p] - mean;

	    // Add squared deviations
	    sum_diff_squared += deviation*deviation;
	    
	  }
	  
	  // 6. CALCULATE THE STANDARD DEVIATION
	  sigma = sqrt(sum_diff_squared/num_cells_phi);

	  // The threshold is the mean plus some number of s.d.
	  threshold = mean + nsigma*sigma;
	  
	  // 7. COMPARE CELL TO THRESHOLD (MEAN + n*sigma)
	  for (int p=0; p!=num_phi; ++p) {
	    
	    int iphi = p+1;
	    if (forward && (iphi % 2 == 0)) continue;

	    // If cell fails test
	    if (bin[p] > threshold) {
		
	      // Write into a file, using unique number for each location
	      int key = 14 * (p + 64*n) + l;
	      HotCells << runnumber << "\t" << key << endl;

	      // Write into a hot cell histogram
	      if (hot_cell_plots) histo_hot[r][l]->Fill(ieta, iphi);
	      
	      ++HotCell[r];
	      ++num_hot_cells;
	    }
	  }
	  
	} // eta
      } // layers
      
      // Output total number of hot cells found in each run
      if (it == goodrun_list.begin()) Output << "\nRun     Hot cells\n" << endl;
    Output << runnumber << "      " << HotCell[r] << endl;
      
    } // loop over runs
    
    // Output summary
    if (goodrun_list.size()) 
      Output << "\nFound an average of " << num_hot_cells*1.0/(goodrun_list.size())*1.0 
	     << " hot cells per run" << endl;
    
  } // if hot_cell_analysis

  //----------------------------------------------------------------
  // DELETE EMPTY HOT CELL HISTOGRAMS
  //----------------------------------------------------------------  

  if (hot_cell_plots) {

    for (map<int,int>::const_iterator it=goodrun_list.begin();
         it!=goodrun_list.end();
         ++it) {

      int runnumber = it->first;
      int r = runnumber - first_goodrun;

      for (int l=0; l!=cal_layers; ++l) {
	
	found_histo_hot[r][l] = false;
	
	// Check if energy histo exists
	if (!found_histo[r][l]) continue;
	
	// Check if any hot cells found
	if ( (histo_hot[r][l]->GetEntries()) != 0 ) 
	  found_histo_hot[r][l] = true;
	
	else histo_hot[r][l]->Delete();
	
      }
    }
  }


  //----------------------------------------------------------------
  // HISTOGRAM DRAW OPTIONS
  //----------------------------------------------------------------

  // Style
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  gStyle->SetTitleBorderSize(0);  
  
  // Stats box
  gStyle->SetOptStat("");
  gStyle->SetStatX(0.99);
  gStyle->SetStatY(0.9);

  // Axis titles
  string plot_xtitle = "ieta";
  string plot_ytitle = "iphi";

  // Canvas
  TCanvas *canvas[cal_layers];
    
  // File name for plots
  char ps_name[100];
  char ps_beg[100], ps_end[100];
  
  //----------------------------------------------------------------
  // DRAW HISTOGRAMS
  //----------------------------------------------------------------

  if (write) {

    Output << "\nWriting plot(s) to file" << endl;

    // RUN LOOP
    for (map<int,int>::const_iterator it=goodrun_list.begin();
         it!=goodrun_list.end();
         ++it) {
      
      int runnumber = it->first;
      int r = runnumber - first_goodrun;

      // Exit if not plotting this run
      if (plot_one && runnumber != plot_run) continue;

      // Construct name      
      if (logz) sprintf (ps_name, "HotCell-Plots-sigma-%g-logz-Run-%i", nsigma, runnumber);
      else sprintf (ps_name, "HotCell-Plots-sigma-%g-Run-%i", nsigma, runnumber);

      if (plot_one) sprintf(ps_name, "%s-layer-%i.eps", ps_name, plot_layer);
      else sprintf(ps_name, "%s.ps", ps_name);

      sprintf (ps_beg, "%s(", ps_name);
      sprintf (ps_end, "%s)", ps_name);
      
      int canvas_count = 0;
      
      // LOOP OVER CANVASES
      for (int l=0; l!=cal_layers; ++l) {
	
	// Exit if not plotting this layer
	if (plot_one && l != plot_layer-1) continue;

	// If cell energy histogram exists
	if (!found_histo[r][l]) continue;
	  
	// Make canvas: with or without hot cell histograms
	if (hot_cell_plots) 
	  canvas[l] = new TCanvas(plot_name[r][l],plot_name[r][l], 1000, 500);
	else canvas[l] = new TCanvas(plot_name[r][l],plot_name[r][l], 500, 500);
	
	canvas[l]->SetLeftMargin(0.13);
	canvas[l]->SetRightMargin(0.13);
	canvas[l]->SetTopMargin(0.1);
	canvas[l]->SetBottomMargin(0.1); 
	
	if (hot_cell_plots) canvas[l]->Divide(2,1);

	if (logz) canvas[l]->SetLogz(true);
	
	// LOOP OVER PADS
	for (int p=0; p!=2; ++p) {
	  
	  // Don't go to 2nd pad if not plotting hot cell histo
	  if (p==1) {
	    if (!hot_cell_plots) continue;
	    if (hot_cell_plots && !found_histo_hot[r][l]) continue;
	  }
	  
	  // Enter pad
	  canvas[l]->cd(p+1);
	  
	  // First pad
	  if (p==0) {

	    // Axes
	    histo[r][l]->SetMinimum(0.000001);
	    histo[r][l]->SetXTitle(plot_xtitle.c_str());
	    histo[r][l]->SetYTitle(plot_ytitle.c_str());
	    histo[r][l]->SetTitleOffset(1.3,"Y");  
	    histo[r][l]->SetTitleOffset(1.2,"X");  
	    
	    // Draw histogram
	    histo[r][l]->Draw("colz");
	  }
	  
	  // Second pad
	  if (p==1) {

	    // Axes
	    histo_hot[r][l]->SetMinimum(0.000001);
	    histo_hot[r][l]->SetXTitle(plot_xtitle.c_str());
	    histo_hot[r][l]->SetYTitle(plot_ytitle.c_str());
	    histo_hot[r][l]->SetTitleOffset(1.3,"Y");  
	    histo_hot[r][l]->SetTitleOffset(1.2,"X");  

	    // Draw histogram
	    histo_hot[r][l]->Draw("colz");
	  }
	  
	  if (logz) canvas[l]->GetPad(p+1)->SetLogz(true);
	  canvas[l]->GetPad(p+1)->Update();
	  
	  if (p==0) ++canvas_count;
	  	  
	} // pads
	
	
	// Plot only one histogram
	if (plot_one) canvas[l]->Print(ps_name);

	// Plot all runs and layers
	else {
	  
	  // Open ps file
	  if (canvas_count == 1) canvas[l]->Print(ps_beg);	  
	  
	  // Close ps file
	  else if (canvas_count == num_histos[r]) canvas[l]->Print(ps_end);
	  	  
	  // Add canvas
	  else canvas[l]->Print(ps_name);
	    
	}
	
      } // layers
    } // runs

    Output << endl;

  } // if write
  
  
  Output.close();
  HotCells.close();
  

}