source: trunk/MagicSoft/Mars/mtemp/mmpi/macros/SCNexSignificanceVSAlphaCut.C@ 6232

// This program produces a ps file with 
// two graphs in the same TPad, but 
// with different axis. 

// The aim is to plot Nex and significance vs 
// cut in alpha parameter.

// As input, this macro needs a root file that is 
// generated by the SupercutsONOFF programs.

// It will read the histograms where Nex and significance 
// vs alpha are stored and will plot this info in a single 
// ps file. 

// The user must write

// 1) Name of the root file with the info (Path included!!)

// 2) Name of the output ps file (with path included!!)

// Depending on the values to be plotted the axis scale
// as well as the axis offsets for labels might have to be 
// modified. If I have time I will modify the program to 
// make it automatic...




#include <string.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <math.h>


gROOT -> Reset();

void SCNexSignificanceVSAlphaCut() 
{

  
  
  char* SCRootFile = {"/.magic/magicserv01/scratch/Daniel/SuperCuts/Mrk421/2004_04_22/4slices_3520_nc/E800_1200_Opt/RootFileDynCuts.root"};

  char* output_ps_file = {"/.magic/magicserv01/scratch/Daniel/SuperCuts/Mrk421/2004_04_22/4slices_3520_nc/E800_1200_Opt/NexSignificanceVSCutinAlpha.ps"};

  
 
  char* ps_title = 
    {"Excess events (green) and significance (red) vs cut in alpha"};




  // Name of the histograms (contained in root file) that contain the number of 
  // excess events and significance vs alpha cut

  TString ExcessEventsHistoName = ("NexVSAlphaTrainThetaRange0_1570mRad");
  TString SignificanceHistoName = ("SignificanceVSAlphaTrainThetaRange0_1570mRad");


  // Axis titles

  char* xaxis_title = "Cut in alpha (degrees)";
  char* yaxis_title = "Excess events";
  char* yaxis_2_title = "Significance";

  // Axis for plot

  const Double_t pad_xlow = 0;
  const Double_t pad_ylow = 0;
  const Double_t pad_xup = 30;
  const Double_t pad_yup = 1000;

  
  
  
  // Axis labels offsets for plot
 
  const Double_t xaxis_label_xoffset_pad = 8;
  const Double_t xaxis_label_yoffset_pad = -300;
  const Double_t yaxis_label_xoffset_pad = -5;
  const Double_t yaxis_label_yoffset_pad = 1200;
  const Double_t yaxis_2_label_xoffset_pad = 35;
  const Double_t yaxis_2_label_yoffset_pad = 1200;

 

  // Vector of pointer to histograms that will be plotted
  // Colors and marker types are also defined 


  const Int_t n_graphs = 2; 
  TGraph* graph_vector [n_graphs];
  Int_t graph_colors_vector[n_graphs] = {3,2};
  Int_t polimarker_style_vector[n_graphs] = {21,25};
  char* additional_info_vector[n_graphs] = 
    {"Excess events", "Significance"};



  // TMP VARIABLES

  Int_t tmp_int = 0;

  TH1F* histo1;
  TH1F* histo2;


 
  Double_t y_axis_2_min_value = 0.0;
  Double_t y_axis_2_max_value = 0.0;
  
  Double_t graph_y_minmax_vector [n_graphs][2];
  // It contains y values min and max for graphs that 
  // are plotted.
  
  Double_t graph_2_scale_factor = 0.0;

  // ******************************************************************
  // Data is read from the histograms contained in the input root file 
  // and stored in vectors 
  // ******************************************************************
  
  TFile rootfile (SCRootFile, "READ");

  histo1 = (TH1F*) rootfile.Get(ExcessEventsHistoName);
  histo2 = (TH1F*) rootfile.Get(SignificanceHistoName);

  

  // check that dimension is the same 

  if (histo1 -> GetNbinsX() != histo2 -> GetNbinsX())
    {
      cout << "Dimension of histogram " << ExcessEventsHistoName 
           << " is not equal to dimension of histogram "<< SignificanceHistoName << endl
           << "Aborting ..." << endl;
      
      Exit(1);
    }

  
  tmp_int = histo1 -> GetNbinsX();
  
  const Int_t vectors_dimension = tmp_int;


  // Vectors that will be used to plot graphs are defined and filled

  Double_t alpha_cut_vector[vectors_dimension] = {0.0};
  Double_t excess_events_vector[vectors_dimension] = {0.0};
  Double_t significance_vector[vectors_dimension] = {0.0};


  for (Int_t i = 0 ; i < vectors_dimension; i++)
    {
      alpha_cut_vector[i] =  histo1-> GetBinCenter(i+1) + ((histo1->GetBinWidth(i+1))/2);
      excess_events_vector[i] = histo1-> GetBinContent(i+1);
      significance_vector[i] = histo2-> GetBinContent(i+1);
    }

  
  

   // Dynamic memory from is released and root file closed 

  delete histo1;
  delete histo2;

  rootfile.Close();


    // Information retrieved from histos is displayed in terminal
  
  cout << endl
       << "Cut in alpha, excess events and significance values retrieved from root file: "
       << endl;

  for (Int_t i = 0 ; i < vectors_dimension; i++)
    {
      cout << alpha_cut_vector[i] << "        " 
           << excess_events_vector[i] << "        " 
           << significance_vector[i] << endl;
    }
  





  //******************************************************
  // Graph 2 is scaled so that it can be plotted in the 
  // same TPad... kind of a trick in root.

  // First, values min and max are needed for the 
  // graphs to be plotted.

  graph_y_minmax_vector [0][0] = pad_ylow;
  graph_y_minmax_vector [0][1] = find_max (excess_events_vector, 
                                           vectors_dimension);
    
cout << " min set to : " <<  graph_y_minmax_vector [0][0] << endl;
cout << " max set to : " <<  graph_y_minmax_vector [0][1] << endl;

  graph_y_minmax_vector [1][0] = 0;
                                        
  graph_y_minmax_vector [1][1] = find_max (significance_vector, 
                                           vectors_dimension);
  //TEST
  /*
  cout << graph_y_minmax_vector [0][1] << "  "  
       << graph_y_minmax_vector [0][0] << endl;

   cout << graph_y_minmax_vector [1][1] << "  "  
       << graph_y_minmax_vector [1][0] << endl;
   */


   
  graph_2_scale_factor = 
    (graph_y_minmax_vector [0][1] - 
     graph_y_minmax_vector [0][0])/(graph_y_minmax_vector [1][1] 
                                    - graph_y_minmax_vector [1][0]);

  //  cout << "scale factor :  " << graph_2_scale_factor << endl;

  for (Int_t i = 0 ; i < vectors_dimension; i++)
    { 
      
      significance_vector[i] = ((graph_2_scale_factor * significance_vector[i])
                           + graph_y_minmax_vector [0][0]);
    }
  
  // ***********  END OF SCALING   ********************************

  // Graphs are initialized and filled with vectors

   
  // Graph i initialized and filled 
  graph_vector[0] = new TGraph(vectors_dimension,
                               alpha_cut_vector, 
                               excess_events_vector);
      

  graph_vector[1] = new TGraph(vectors_dimension,
                               alpha_cut_vector, 
                               significance_vector);
  
  // Some properties of graphs are set

  for (Int_t i = 0; i < n_graphs; i++)
    {
      graph_vector[i]->SetTitle();
      graph_vector[i]->SetFillColor(19);
      graph_vector[i]->SetLineColor(graph_colors_vector[i]);
      graph_vector[i]->SetLineWidth(1);
      graph_vector[i]->SetMarkerColor(graph_colors_vector[i]);
      graph_vector[i]->SetMarkerStyle(polimarker_style_vector[i]);
      graph_vector[i]->SetMarkerSize(1);
    }


// Graphgrams are plotted


 TCanvas *Window = new TCanvas("Plot","Plot",600,600); 
  
 gStyle->SetOptStat(0);   
 // So that statistic information is not printed

  Window -> SetFillColor(10);
  gStyle -> SetFrameFillColor(10);
  gStyle -> SetPadLeftMargin (0.15);
  gStyle -> SetPadRightMargin (0.15);
  gStyle -> SetPadTopMargin (0.05);
  gStyle -> SetPadBottomMargin (0.15);

  TPaveText *title = new TPaveText (.05,.95,.95,.99);
  title -> SetFillColor (38);  
  title -> AddText (ps_title);
  title -> Draw();

 
  
  
  /* TPad(const char *name, const char *title, 
     Double_t xlow, Double_t ylow, 
     Double_t xup, Double_t yup, Color_t color, 
     Short_t bordersize,Short_t bordermode) 
     : TVirtualPad(name,title,xlow,ylow,xup,yup,color,
     bordersize,bordermode) 
  */

  TPad *pad = new TPad ("plot", "plot",0.0, 0.0, 1.0, 0.92, 
                        10, 0, 0); 
  // I do not want any border effect.
  pad -> Draw();
   
   
  // Axis and axis labels for pad1 are created and drawn.

  pad -> cd();
  pad -> SetGrid(1,1); // Grid for x axis and for y axis.
  pad -> DrawFrame(pad_xlow,pad_ylow,pad_xup,pad_yup); 
  // dimensions of the pad are set.
  
  
   
  // Labels of the graph are set

  TText *t = new TText();
  t->SetTextFont(62);
  t->SetTextColor(1);
  t->SetTextSize(0.05);
  t->SetTextAlign(11);

  // X axis
  t->SetTextAngle(0);
  t->DrawText(xaxis_label_xoffset_pad,
              xaxis_label_yoffset_pad,
              xaxis_title);
  // Y axis
  t->SetTextAngle(90); 
  t->DrawText(yaxis_2_label_xoffset_pad,
              yaxis_2_label_yoffset_pad, 
              yaxis_2_title);

 
  

  // Graphs are plotted in pad


  for (Int_t j = 0; j < n_graphs; j++)
    {  
      
      pad -> cd();
      graph_vector[j]-> Draw("LP");
      // "A" option is removed in order to 
      //control de dimensionsn of the pad, 
      // by means of the " DrawFrame()" function.
      
    }
  
  // **************************************************
  // Second Axis is drawn with its label
  
  // First of all we must look for units at
  // bottom and top of axis. This is found 

  
  y_axis_2_min_value = ((pad_ylow - graph_y_minmax_vector [0][0])/
                        graph_2_scale_factor);

  y_axis_2_max_value = ((pad_yup - graph_y_minmax_vector [0][0])/
                        graph_2_scale_factor);

  

  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
                            gPad->GetUxmax(), gPad->GetUymax(),
                            y_axis_2_min_value,
                            y_axis_2_max_value,510,"+L");

   axis->SetLineColor(1);
   axis->SetTextColor(1);
   axis->Draw();

    // Y axis 2
  t->SetTextAngle(90); 
  t->DrawText(yaxis_label_xoffset_pad,
              yaxis_label_yoffset_pad, 
              yaxis_title);

  // Graph information is produced and drawn.

  // Spectrum names
  /*
  char info_string[150] = 0;
  char info_string_tmp [20] = 0;
  Double_t xpos = 0;
  Double_t ypos = 0;
  TText *tt = new TText();
  tt->SetTextFont(62);
  tt->SetTextSize(0.04);
  tt->SetTextAlign(12);
  tt->SetTextAngle(0);
  
  for (Int_t i = 0; i < n_graphs; i++)
    {
      if (i < 1)
        {xpos = 0.15; ypos =  0.91;} 
      else 
        {xpos = 0.15; ypos = 0.85;}
      
      //xpos = 0.02 + i*0.12 + i*(i-1)*0.15; //for three curves
      //xpos = 0.1 + i*0.4; // for two curves
      strcpy(info_string, additional_info_vector[i]);
      tt->SetTextColor(graph_colors_vector[i]);
      Window -> cd();
      tt->DrawText(xpos,ypos,info_string);
    }
  */


  Window -> SaveAs(output_ps_file);

}




// FUNCTION DEFINITIONS





Double_t find_value (Double_t vector[], Int_t dimension, Int_t min_max)
{
  Double_t value = vector[0];
  if (min_max < 0.5) // we look for min
    {
      for (Int_t i = 1; i < dimension; i++)
        {
          if (value > vector[i])
            {
              value = vector[i];
            }
        }
    }

  else // we look for max
    {
      for (Int_t i = 1; i < dimension; i++)
        {
          if (value < vector[i])
            {
              value = vector[i];
            }
        }
    }

  return value;
}
  


Double_t find_min (Double_t vector[], Int_t dimension)
{
  Double_t min = 0.0;
  min = find_value (vector, dimension, 0);
  return min;
}



Double_t find_max (Double_t vector[], Int_t dimension)
{
  Double_t max = 0.0;
  max = find_value (vector, dimension, 1);
  return max;
}