Trigger/TrigAnalysis/TrigInDetUserAnalysis/Analysis/src/comparitor.cxx

  • This has all of the names of the histograms produced at the end of the process, (i.e. which I have put in ElectronSlice/EFHistograms and ElectronSlice/L2Histograms)
  • This saves the histograms in .png and .pdf formats.

Includes:

#include <sys/time.h>
#include <stdlib.h>

#include <iostream>
#include <string>
#include <vector>

#include "Resplot.h"
#include "utils.h"
#include "label.h"
#include "DrawLabel.h" 


#include "TFile.h"
#include "TH1D.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TF1.h"

Not sure:

int usage(const std::string& name, int status) { 
  std::cout << "Usage: " << name << "  test_file.root reference_file.root chain name 1 chain name 2 etc " << std::endl;
  return status;
}

Time string:

std::string stime() { 
  time_t t;
  time(&t);
  return label("%s", ctime(&t) );
}

Method for calculating the xrange required for individual histograms:

void xrange(TH1D* h, bool symmetric=false ) { 
  int ilo = 1;
  int ihi = h->GetNbinsX();
  for ( ; ilo<=ihi ; ilo++ ) if ( h->GetBinContent(ilo)>0 ) break; 
  for ( ; ihi>=ihi ; ihi-- ) if ( h->GetBinContent(ihi)>0 ) break;
 
  int delta_lo = ilo-1;
  int delta_hi = h->GetNbinsX()-ihi;

  if ( symmetric ) { 
    if ( delta_hi<delta_lo ) h->GetXaxis()->SetRange( 1+delta_hi, ihi );
    else                     h->GetXaxis()->SetRange( 1+delta_lo, h->GetNbinsX()-delta_lo );
  }
  else { 

    if ( ilo>1 ) ilo--;
    if ( ilo>1 ) ilo--;
    if ( ihi<h->GetNbinsX() ) ihi++;
    if ( ihi<h->GetNbinsX() ) ihi++;

    h->GetXaxis()->SetRange(ilo,ihi);
  }

}

Producton of a Legend class for later use:

class Legend { 

public:

  Legend() : mleg(0) { }  

  Legend(double x1, double x2, double y1, double y2) { 
    mleg = new TLegend(x1, y1, x2, y2);
    mleg->SetBorderSize(0);
    mleg->SetTextFont(42);
    mleg->SetTextSize(0.025);
    mleg->SetFillStyle(3000);
    mleg->SetFillColor(0);
    mleg->SetLineColor(0);
  } 


  Legend(const Legend& legend) : mleg(legend.mleg) { } 

  //  ~Legend() { if ( mleg!=0 ) delete mleg; } 
  ~Legend() { } 

  TLegend* legend() { return mleg; } 

  TLegend* operator->() { return mleg; }

private:

  TLegend* mleg;

}; 

(GC) Generic plotter class - better to have one of these - make sure it can be configured however you like, line styes, marker types, legends etc. Now a template so can be used for TH1D and TH2D etc.

template<typename T>
class tPlotter { 

public: 

  tPlotter(T* _htest=0, T* _href=0 , const std::string& s="") : 
    m_htest(_htest), m_href(_href), 
    m_plotfilename(s) { }

  
  tPlotter(const tPlotter& p) : 
    m_htest(p.m_htest), m_href(p.m_href), 
    m_plotfilename(p.m_plotfilename){ 
    
  }


  /// sadly, root manages all the histograms (does it really? 
  /// who can tell) so we mustn't delete anything just in case
  ~tPlotter() { } 

  std::string plotfilename() const { return m_plotfilename; }

  //  void Draw(const std::string& chain, unsigned int i, Legend& leg ) { 
  void Draw( unsigned int i, Legend& leg ) { 
    if ( href() && htest() ) {
      gStyle->SetOptStat(0);
      href()->SetLineColor(i+1);
      href()->SetLineStyle(2);
      href()->SetMarkerColor(href()->GetLineColor());
      href()->SetMarkerStyle(24+i%4);
      htest()->SetLineColor(i+1);
      htest()->SetLineStyle(1);
      htest()->SetMarkerColor(htest()->GetLineColor());
      htest()->SetMarkerStyle(20+i%4);

      std::cout << "\n\nhist " << href()->GetName() << std::endl; 

      if(i==0) { 
   href()->Draw("ep");
   // href()->Draw("lhist");
   //   xrange(href());
      }
      if ( std::string(href()->GetName()).find("_vs_")!=std::string::npos  ||
      std::string(href()->GetName()).find("_eff")!=std::string::npos ) {
   href()->Draw("hist same][");
      }
      else {
   href()->Draw("hist same");
      }
      htest()->Draw("ep same");
      // href()->Draw("lhistsame");
      // htest()->Draw("lhistsame");
      leg->AddEntry(href(),(m_plotfilename+" reference").c_str(),"lp");
      leg->AddEntry(htest(),m_plotfilename.c_str(),"lp");
      leg->Draw();
    }
  }


  /// print the output 
  void Print(const std::string& s="") {
    std::cout<<gPad<<std::endl;
    if ( s!="" ) gPad->Print(s.c_str());
    else         gPad->Print(m_plotfilename.c_str());
  } 

  T* htest() { return m_htest; }
  T* href()  { return m_href; }

private:

  T* m_htest;
  T* m_href;


  std::string m_plotfilename;

};

typeDef = Create a synonym name for this type. i.e. our class 'tPlotter' (in the case of TH1D) may be referred to simply as Plotter, for clarity.

typedef tPlotter<TH1D> Plotter;

Create a class known as 'plots', which is a vector of 'Plotter' (tPlotter) class types.

class Plots : public std::vector<Plotter> { 
  
public:

  Plots(const std::string& s="") : m_name(s) { }
  
  void Max(double scale=1.1) {  
    
    if ( size()<1 ) return;

    double tmax = 0;
    for ( unsigned i=0 ; i<size() ; i++ ) { 
      if ( i==0 ) tmax = at(i).href()->GetMaximum();
      if ( tmax<at(i).href()->GetMaximum() )  tmax = at(i).href()->GetMaximum(); 
      if ( tmax<at(i).htest()->GetMaximum() ) tmax = at(i).htest()->GetMaximum(); 
    }
    
    for ( unsigned i=0 ; i<size() ; i++ ) {
      at(i).href()->SetMaximum(scale*tmax);
      at(i).htest()->SetMaximum(scale*tmax);
    } 

  }


  void xrange(bool symmetric=false) { 
    for ( unsigned i=0 ; i<size() ; i++ ) { 
      ::xrange( at(i).href(), symmetric );
      ::xrange( at(i).htest(), symmetric );
    }
  }
   
  void Draw( Legend& leg ) {  
    //   Max();
    for ( unsigned i=0 ; i<size() ; i++ ) at(i).Draw( i, leg);
    
    DrawLabel(0.1, 0.02, "built on "+stime(), kBlack, 0.03 );
  }

private:
  
  std::string m_name;

};


Start of the main program code.


int main(int argc, char** argv) { 

  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadTopMargin(0.1);
  
  std::cout << argc << std::endl;

  if ( argc<4 ) return usage(argv[0], -1);

Creation of vectors to hold chain names. Creation of vectors to hold histogram titles and names.

  TFile ftest(argv[1]);
  TFile fref(argv[2]);
  std::vector<std::string> chains;
  for(int i=3; i<argc; i++){
    chains.push_back(argv[i]);
  }
  std::vector<std::string> chainnames;
  chainnames.resize(chains.size());
  chainnames.clear();

  const int Nhistos = 14;
  std::string histos[Nhistos] = { 
    /// efficiencies
    "pT_eff",       
    "eta_eff", 
    "phi_eff", 
    "d0_eff", 
    /// standard residuals
    "ipT_res", 
    "eta_res", 
    "phi_res", 
    "z0_res",
    /// residuals vs track parameters
    "ript_vs_pt/sigma",
    "reta_vs_pt/sigma",
    "reta_vs_eta/sigma",
    "rzed_vs_eta/sigma",
    "rzed_vs_pt/sigma",
    "rzed_vs_zed/sigma"
  };
  
  std::string histonames[Nhistos] = { 
    "Efficiency P_{T}", 
    "Efficiency #eta", 
    "Efficiency #phi", 
    "Efficiency d0", 
    
    "Residual 1/P_{T}", 
    "Residual #eta", 
    "Residual #phi", 
    "Residual z0",
    
    "Residual 1/P_{T} vs P_{T}",
    "Residual #eta vs P_{T}",
    "Residual #eta vs #eta",
    "Residual z vs #eta",
    "Residual z vs P_{T}",
    "Residual z vs z"
  }; 


  std::cout<<"size of chains " << chains.size()<<std::endl;
  
  if ( chains.empty() ) return 0;

Loop over histograms, set up Plots object, setup canvas, setup Legend object, setup vectors for Mean and RMS. Work out x-axis labels.

  for ( int i=0 ; i<Nhistos ; i++ ) {

    std::cout << "main() processing histos[" << i <<  "] " << histos[i] << std::endl;

    //    std::vector<Plotter> plots;
    Plots plots;
    plots.clear();
  
    TCanvas* c1 = new TCanvas(label("canvas-%d",i).c_str(),"histogram",800,600);
    c1->cd();

    Legend legend;
    // efficiencies or residuals?
    if ( histos[i].find("_eff")!=std::string::npos ) legend = Legend( 0.2, 0.3, 0.4-chains.size()*0.07, 0.40 );
    else                                             legend = Legend( 0.30, 0.4, 0.9-chains.size()*0.07, 0.90 );
   
    std::vector<std::string> Mean;
    std::vector<std::string> RMS;
  
    Mean.clear();
    RMS.clear();

    std::string plotname; 

    std::string xaxis = "";
    std::string yaxis = "";

    /// work out axis labels
    /// x axis
    if ( histos[i].find("vs_pt")!=std::string::npos )  xaxis = "true P_{T} [GeV]"; 
    if ( histos[i].find("vs_eta")!=std::string::npos ) xaxis = "true #eta"; 
    if ( histos[i].find("vs_zed")!=std::string::npos ) xaxis = "true z_{0} [mm]"; 

    if ( histos[i].find("ript")!=std::string::npos ) yaxis = "rms_{95} 1/P_{T} residual [GeV^{-1}]"; 
    if ( histos[i].find("reta")!=std::string::npos ) yaxis = "rms_{95} #eta residual"; 
    if ( histos[i].find("rzed")!=std::string::npos ) yaxis = "rms_{95} z_{0} residual [mm]"; 

Loop over chains. Grab histograms from test and reference chains and put them into htest and href. Set their titles and axis labels. Define names for each plot to save (EF/L2 + histogram name). Put Plotter object for test and href onto Plots vector object.

    
    for ( unsigned int j=0; j<chains.size(); j++)  {

      std::cout << "main() processing chain[" << j << "] " << chains[j] << std::endl;
            
      TH1D* htest = (TH1D*)ftest.Get((chains[j]+"/"+histos[i]).c_str()) ;
      TH1D* href  = (TH1D*)fref.Get((chains[j]+"/"+histos[i]).c_str()) ;
      
      if ( htest==0 ) { 
   std::cerr << "main() no test histogram : " << (chains[j]+"/"+histos[i]) << std::endl;
   return -1;
      }
      
      if ( href==0 ) { 
   std::cerr << "main() no ref histogram : " << (chains[j]+"/"+histos[i]) << std::endl;
   return -1;
      }
      
      if ( xaxis!="" ) {
   htest->GetYaxis()->SetTitleOffset(1.5);
   href->GetYaxis()->SetTitleOffset(1.5);
   htest->GetXaxis()->SetTitle(xaxis.c_str());
   htest->GetYaxis()->SetTitle(yaxis.c_str());
   href->GetXaxis()->SetTitle(xaxis.c_str());
   href->GetYaxis()->SetTitle(yaxis.c_str());
      }

      if ( histos[i].find("vs_pt")!=std::string::npos ) { 
   htest->GetXaxis()->SetRangeUser(0,150);
   href->GetXaxis()->SetRangeUser(0,150);
      }

      //axis labels
      if(chains[j].find("EF")!=std::string::npos){
   htest->SetTitle(("Event Filter "+histonames[i]).c_str());
        href->SetTitle(("Event Filter "+histonames[i]).c_str());
   plotname = "EF_";
      }
      else if(chains[j].find("L2")!=std::string::npos){
   htest->SetTitle(("Level 2 "+histonames[i]).c_str());
   href->SetTitle(("Level 2 "+histonames[i]).c_str());
   plotname = "L2_";
      }

      plotname += histos[i]; 

      /// replace the "/" in the filename so we don't try to 
      /// make the plots in subdirectories  
      replace(plotname, "/", "_"); 

      //change efficiencies
      if(histos[i].find("_eff")!=std::string::npos) {
   if(histos[i]=="pT_eff") { 
     //     href->SetMinimum(60);
     href->GetXaxis()->SetRangeUser(0,150);
   }
   href->SetMinimum(70);
   href->SetMaximum(105);
      }

      plots.push_back( Plotter( htest, href, chains[j] ) );
      
      //      plots.back().Draw( j, legend );

TF1 is a generic 1-dimensional function defined between an upper and lower limit. FitNull95 is defined in Resplot.h and Resplot.cxx . FILL IN WHAT IT DOES.

      if(histos[i].find("_res")!=std::string::npos)  {
   
   TF1* d95 = Resplot::FitNull95( htest );
   
   double   mean_95 = d95->GetParameter(1);
   double  dmean_95 = d95->GetParError(1);
   double    rms_95 = d95->GetParameter(2);
   double   drms_95 = d95->GetParError(2);

   std::cout <<  histos[i] 
        << "\tmean: " << mean_95 << " +- " << dmean_95 
        << "\trms: "  <<  rms_95 << " +- " << drms_95 << std::endl; 
          
   if(histos[i]=="ipT_res") {
     Mean.push_back(label("mean_{95} = %2.3lf #pm %2.3lf GeV^{-1}", mean_95, dmean_95));
     RMS.push_back(label( "rms_{95}   = %2.3lf #pm %2.3lf GeV^{-1}", rms_95,  drms_95));
   }
   else { 
     Mean.push_back(label("mean_{95} = %2.5lf #pm %2.5lf ", mean_95, dmean_95));
     RMS.push_back(label( "rms_{95}   = %2.5lf #pm %2.5lf ", rms_95,  drms_95));
   }
   
   //   href->Rebin(2);
   //   htest->Rebin(2);

      }
      
    }     

Set the maximum (I assume height and width?) of each plot to allow for room for the Legend. then draw the Legend.

    if ( histos[i].find("_res")!=std::string::npos) {
      /// use 20 times the max from any of the plots as the upper limit so there's room for the key  
      plots.Max(20);  

      if ( histos[i].find("ipT")!=std::string::npos) plots.xrange(false);   // sets the xrange to something sensible (and symmetric)
      else                                           plots.xrange(true);  // sets the xrange to something sensible
    }

    if ( histos[i].find("_vs_")!=std::string::npos) { 
      /// use 1.5 times the acyual maximum so room for the key
      plots.Max(1.5); 
    }

    plots.Draw( legend );

Decide whether to set the y-axis to a log axis.

    if ( histos[i].find("_res")!=std::string::npos) { 
      c1->SetLogy();
      for ( unsigned j=0 ; j<chains.size() ; j++ ) { 
   DrawLabel( 0.63, (0.6-j*0.035), Mean[j], j+1 );
   DrawLabel( 0.63, (0.6-0.035*chains.size()-j*0.035), RMS[j], j+1 );
      }
    }
    else  c1->SetLogy(false);
    //    if ( histos[i].find("rzed")!=std::string::npos ) c1->SetLogy(true);

Plot the histograms into .png and .pdf files to be saved. Delete the canvas and end the code.

    plots.back().Print( plotname+".pdf" );
    plots.back().Print( plotname+".png" );
    
    //    std::cout << "delete c1 " << c1 << std::endl;
    delete c1;
    //  std::cout << "deleted " << std::endl;
    
  }

  
  return 0;
}
Edit | Attach | Watch | Print version | History: r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r1 - 2013-03-05 - JacobKempster
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback