comparitor.cxx

Go to the documentation of this file.

 
// cppcheck-suppress-file stlIfStrFind; cannot use later than C++17 for this standalone code
 
 
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <sys/time.h>
#include <sys/stat.h>
#include <sys/types.h> 
 
#include <iostream>
#include <string>
#include <vector>
#include <set>
#include <algorithm>
#include <regex>
 
#include "TrigInDetAnalysis/Efficiency1D.h"
 
#include "ReadCards.h"
 
#include "Resplot.h"
#include "utils.h"
#include "label.h"
#include "DrawLabel.h" 
 
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TH2F.h"
#include "TH1D.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TF1.h"
#include "TPave.h"
#include "TPaveStats.h"
#include "TEfficiency.h"
#include "TColor.h"
 
#include "computils.h"
 
#include "AtlasStyle.h"
#include "AtlasLabels.h"
 
#include "default_panels.h"
 
#include "TError.h"
 
 
const std::string reset  = "\033[0m";
const std::string black  = "\033[0;30m";
const std::string red    = "\033[0;31m";
const std::string green  = "\033[0;32m";
const std::string blue   = "\033[0;34m";
 
bool fulldbg = false;
 
extern bool LINEF;
extern bool LINES;
 
bool EPS = false;
 
 
 
void normaliseBinWidth( TH1* h ) {
 
  if ( !h ) return;
  
  int    nbins = h->GetNbinsX();
  double width = h->GetBinWidth(1);
 
  double eps = 1e-12;
  
  bool   variable = false;
  
  for ( int i=2 ; i<=nbins ; i++ ) { 
    if ( std::fabs( h->GetBinWidth(i)-width )>eps ) { 
    variable = true;
    break;
      }
  }
  
  if ( !variable ) return;
  
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
      double w = h->GetBinWidth(i);
      h->SetBinContent( i, h->GetBinContent(i)/w );  
      h->SetBinError(   i,   h->GetBinError(i)/w );  
  }
  
}
 
 
 
 
void testfit( TH1* htest, const std::string& s ) {
 
  if ( htest && s.find("invmass")!=std::string::npos ) { 
 
    TVirtualPad* p = gPad;
 
    TCanvas canv("canv", "", 700, 600 );
    canv.cd();
    
    std::cout << "testfit: s" << s << "\t" << htest << std::endl;
  
    // TF1 finv( "inv", "pow([0]*exp( (x-[1])*([1]-x)/([2]*[2]) ), [6]) + [3]+x*([4]+x*[5])" ); 
    TF1 finv( "inv", "pow([0]/([2]+([1]-x)*([1]-x)), [6]) + [3]+x*([4]+x*[5])" ); 
 
    finv.SetParameter(0, 10000);
    finv.SetParameter(1, 91);
    finv.SetParameter(2, 7);
    finv.SetParameter(3, 1000);
    finv.SetParameter(4, 0.01);
    finv.SetParameter(5, 0);
    finv.SetParameter(6, 1);
 
    finv.SetLineWidth(1);
 
    finv.SetParLimits(1, 89, 93 );
    finv.SetParLimits(2, 5, 50 );
    
    std::vector<double> errors(htest->GetNbinsX());
    for ( size_t i=0 ; i<errors.size() ; i++ ) {
      errors[i] = htest->GetBinError(i+1);
      htest->SetBinError(i+1, 20);
      if (htest->GetBinContent(i+1)==0 ) htest->SetBinError(i+1, 0);
    }
 
    htest->Fit( &finv, "", "", 75, 105 );
 
    static int i = 0;
 
    for ( size_t i=0 ; i<errors.size() ; i++ ) htest->SetBinError(i+1, errors[i] );
    
    htest->GetXaxis()->SetRangeUser(55,125);
    htest->DrawCopy("e1");
    
    gPad->Print( ("invmass-"+std::to_string(i++)+".pdf").c_str() );
 
    std::cout << "inv mass: " << finv.GetParameter(1) << " +- " << finv.GetParError(1) << std::endl;
    
    p->cd();
    
  }
 
}
 
 
double total_eff( TH1* hn, TH1* hd ) { 
  double n = 0;
  double d = 0;
 
  for ( int i=0 ; i<=hn->GetNbinsX()+1 ; i++ ) n += hn->GetBinContent(i);
  for ( int i=0 ; i<=hd->GetNbinsX()+1 ; i++ ) d += hd->GetBinContent(i);
 
  std::cout << "num: " << n << "\tden: " << d << "\t";
  
  if ( d!=0 ) return (n/d);
 
  return 0;
}
 
 
 
extern int colours[6];
extern int markers[6];
 
TCanvas* fitcanvas     = 0;
TCanvas* fitcanvas_rec = 0;
 
std::string eround( double x, int decfigs=4 ) { 
  char ss[1024];
  std::string s= "%" + std::to_string(decfigs+2) + "." + std::to_string(decfigs) + "lf";
  std::sprintf( ss, s.c_str(), x );
  return ss;
}
 
 
void arsaway( std::string& s, const std::string& pattern, const std::string& rep ) { 
  std::regex regexPattern(pattern);
  s = std::regex_replace( s, regexPattern, rep );
}
 
 
 
 
TH1F* trans( TH1F* h, bool t=false ) { 
 
  std::cout << "trans in" << std::endl;
 
  if ( !t ) return h;
 
  std::cout << "trans processing ..." << std::endl;
 
  h->SetDirectory(0);
 
  std::vector<double> xlim( h->GetNbinsX()+1 );
 
  for ( int i=0 ; i<h->GetNbinsX()+1 ; i++ ) { 
    xlim[i] = h->GetBinLowEdge(i+1)*(1-0.054);
  }
 
  TH1F* hn = new TH1F( (std::string(h->GetName())+"-transx").c_str(), h->GetTitle(), h->GetNbinsX(), &xlim[0] );
  hn->SetDirectory(0);
 
  hn->SetLineStyle(h->GetLineStyle());
  hn->SetLineColor(h->GetLineColor());
  hn->SetLineWidth(h->GetLineWidth());
  hn->SetMarkerStyle(h->GetMarkerStyle());
  hn->SetMarkerColor(h->GetMarkerColor());
  hn->SetMarkerSize(h->GetMarkerSize());
 
  for ( int i=0 ; i<h->GetNbinsX() ; i++ ) { 
    if ( hn->GetBinLowEdge(i)<14 ) continue; 
    hn->SetBinContent( i+1, h->GetBinContent( i+1 ) );
    hn->SetBinError( i+1, h->GetBinError( i+1 ) );
  }
 
  delete h;
 
  return hn;
 
}
 
 
 
 
void SetZeros( TH2D* h ) { 
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    for ( int j=1 ; j<=h->GetNbinsY() ; j++ ) { 
      int ibin = h->GetBin( i, j );
      if ( h->GetBinContent(ibin)==0 ) h->SetBinContent(ibin, 0.1); 
    }
  }
}
 
void Scale( TH1* h, double d=1 ) { 
  if ( d==1 ) return;
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    h->SetBinContent( i, h->GetBinContent(i)*d ); 
    h->SetBinError( i, h->GetBinError(i)*d ); 
  }
}
 
 
TH1F* Rebin( TH1F* h, double f ) { 
  
  std::cout << "\nREBIN: " << h->GetName() << " :: " << f << std::endl; 
 
  if ( int(f) == f ) { 
    TH1F* n = (TH1F*)h->Clone((std::string(h->GetName())+"-foeda").c_str()); n->SetDirectory(0);
    n->Rebin(int(f));
    return n;
  }
 
  int     fi=int(f);
  double fb=(f-fi);
 
  for ( int i=0 ; i<10 ; i++ ) { 
    fb *= 10; 
    if ( std::fabs(int(fb)-fb)<1e-6 ) break; 
  }
 
  std::vector<double> limits;
  std::vector<double> contents;
    
  for ( int i=1 ; i<=h->GetNbinsX()+1 ; i++ ) { 
    limits.push_back( h->GetBinLowEdge(i) );
    contents.push_back( h->GetBinContent(i) );
    if ( i<fb ) { 
      for ( int ib=1 ; ib<fi ; ib++ ) contents.back() += h->GetBinContent(++i);
    }
  }
  
  TH1F* n = new TH1F( "foeda", h->GetTitle(), limits.size()-1, &limits[0] ); n->SetDirectory(0);
  
  for ( size_t i=0 ; i<contents.size() ; i++ ) n->SetBinContent( i+1, contents[i] );
 
  return n;
}
 
 
 
 
TH1F* rebin_log( TH1F* h, bool low=false ) { 
 
  double limits[40] = {
    0.5, 1.5, 2.5, 3.5, 4.5, 
    5.5, 6.5, 7.5, 8.5, 9.5, 
    10.5, 11.5, 12.5, 13.5, 14.5, 
    15.5, 16.5, 17.5, 18.5, 19.5, 
    20.5, 21.5, 22.5, 23.5, 24.5, 
    25.5, 26.5, 28.5, 30.5, 32.5, 
    35.5, 38.5, 42.5, 46.5, 50.5, 
    55.5, 60.5, 66.5, 72.5, 78.5 
  };
 
  double lowlimits[24] =  { 
    0.5, 1.5, 2.5, 3.5, 4.5, 
    5.5, 6.5, 7.5, 8.5, 9.5, 
    10.5, 11.5, 12.5, 14.5, 17.5, 
    20.5, 24.5, 29.5, 35.5, 42.5, 
    50.5, 60.5, 72.5, 86.5, 
  };
 
  TH1F* n;
 
  if ( low ) n = new TH1F( "foeda", h->GetTitle(), 39, limits ); 
  else       n = new TH1F( "foeda", h->GetTitle(), 23, lowlimits );
 
  n->SetDirectory(0);
 
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    n->Fill( h->GetBinCenter(i), h->GetBinContent(i) ); 
  }
 
 
  for ( int i=1 ; i<=n->GetNbinsX() ; i++ ) { 
    n->SetBinError(i, std::sqrt( n->GetBinContent(i) ) ); 
  }
 
  return n;
}
 

 
typedef std::map<std::string,std::string> chainmap_t;
 

 
template<typename T=TH1F>
T* Get( TFile& f, const std::string& n, const std::string& dir="",
    const chainmap_t* chainmap=0, std::vector<std::string>* saved=0 ) { 
 
  std::cout << "\n\nGet: " << n << "\tdir: " << dir << "\tchainmap: " << chainmap << std::endl;
 
  std::string name;
 
  size_t pos = n.find('+');
  if ( pos!=std::string::npos ) name = n.substr( 0, pos );
  else                          name = n;
 
  std::cout << "name: " << name << "\t(file: " << f.GetName() << ")" << std::endl;
 
  T* h = (T*)f.Get( name.c_str() );
 
  bool printed = false;
 
  if ( h || dir.empty() || name.find(dir)!=std::string::npos ) { 
    std::cout << "Get() name 0 " << name << " :: " << h << std::endl;
    printed = true;
  }
  else { 
    std::string fname = dir+'/'+name;
    h = (T*)f.Get( fname.c_str() );
    //    std::cout << "Get() name 1 " << name << " :: " << h << std::endl;
  }
 
  if ( h == 0 ) { 
    if ( chainmap && chainmap->size()!=0 ) { 
      for ( chainmap_t::const_iterator itr=chainmap->begin() ; itr!=chainmap->end() ; itr++ ) { //cppcheck-suppress postfixOperator 
    if ( contains( name, itr->first ) ) { 
      
      std::cout << "\tmatch: " << itr->first << " -> " << itr->second << std::endl;
      
      name.replace( name.find(itr->first), itr->first.size(), itr->second );
      
      h = (T*)f.Get( name.c_str() );
 
      if ( h==0 ) { 
        if ( name.find(dir)==std::string::npos ) name = dir + "/" + name;
      
        h = (T*)f.Get( name.c_str() );
      }
 
      break;
    }
      }
    }
  }
 
  
  if ( ! printed ) { 
    std::cout << "Get() name 1 " << name << " :: " << h << std::endl;
  }
 
 
  if ( saved ) saved->push_back( name );
 
 
  if ( h && std::string(h->ClassName()).find("TEfficiency")==std::string::npos ) {  
    std::cout << "TH " << h->GetName() << std::endl; 
    if ( h ) h->SetDirectory(0);
  }
    
 
  return h;
}
 
 
template<typename T>
bool isTEfficiency( T* h ) {
  if ( h && std::string(h->ClassName()).find("TEfficiency")!=std::string::npos ) return true;
  return false;
}
 
 
 
class bands { 
 
public:
 
  bands() { } 
 
  bands( const bands& b ) = default;
  bands& operator= ( const bands& b ) = default;
 
  bands( const std::vector<double>& limits, const std::vector<std::string>& labels ) 
  : m_limits(limits), m_labels(labels)
  {   }
 
  void range( const std::string& chain, TH1F* h ) {
    for ( size_t i=0 ; i<m_labels.size() ; i++ ) { 
      if ( chain.find(m_labels[i])!=std::string::npos ) return range( i, h );
    }
    return range( m_labels.size(), h );
  } 
 
 
  void range( size_t i, TH1F* h ) { 
  
    double minx = h->GetBinLowEdge(1);
    double maxx = h->GetBinLowEdge(h->GetNbinsX()+1);
 
    if      ( i>=m_limits.size() )  { maxx+=1; minx=maxx; }
    else if ( (i+1)==m_limits.size() ) minx=m_limits[i];
    else { 
      minx = m_limits[i];
      maxx = m_limits[i+1];
    }
  
 
    for ( int j=1 ; j<=h->GetNbinsX() ; j++ ) { 
 
      std::cout << "range: " << j << "\tminx: " << minx << "\t" << maxx << std::endl; 
      if ( ! ( h->GetBinCenter(j)>=minx && h->GetBinCenter(j)<maxx ) ) {
    h->SetBinContent(j,0);
    h->SetBinError(j,0);
      }
    }
  }
  
  
private: 
 
  std::vector<double>      m_limits;  
  std::vector<std::string> m_labels;  
 
};
 
 
TH1F* d0rebin( TH1F* h ) { 
 
  return h;
  
  std::vector<double> limits;
 
  for ( int i=1 ; i<h->GetNbinsX() ; i++ ) { 
    limits.push_back( h->GetBinLowEdge(i) );
  }
 
  std::vector<double> alimits;
  std::vector<double> nlimits;
 
  double x = 0; 
 
  //  for ( int i=0 ; x<20 ; i++ ) { 
  while ( x<20 ) { 
  
    alimits.push_back(x);
 
#if 0
    if      ( x<0.5 ) x+=0.2;
    else if ( x<1.0 ) x+=0.2;
    else if ( x<1.5 ) x+=0.3;
    else if ( x<2.0 ) x+=0.4;
    else if ( x<2.5 ) x+=0.5;
    else              x+=0.6;
#else
    if      ( x<0.5 ) x+=0.1;
    else if ( x<1.0 ) x+=0.2;
    else if ( x<1.5 ) x+=0.2;
    else if ( x<2.0 ) x+=0.2;
    else if ( x<2.5 ) x+=0.2;
    else if ( x<3.0 ) x+=0.3;
    else              x+=0.3;
#endif
 
  }
 
  for ( size_t i=alimits.size() ; i-- ; ) {
    if ( alimits[i]<=15 ) nlimits.push_back( -alimits[i] );
  } 
  
  for ( size_t i=1 ; i<alimits.size() ; i++ ) { 
    if ( alimits[i]<=15 )  nlimits.push_back( alimits[i] );
  }
 
  std::cout << "limits: " << limits.size() << " " << nlimits.size() << std::endl;
 
  TH1F* hnew = new TH1F( "h", h->GetTitle(), nlimits.size()-1, &nlimits[0] ); 
 
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    double d = h->GetBinCenter(i);
    int j = hnew->FindBin( d );
    hnew->SetBinContent( j, hnew->GetBinContent(j)+h->GetBinContent(i) ); 
  }
  
  return hnew;
}
 
 
 
 
int usage(const std::string& name, int status, const std::string& err_msg="" ) { 
  if ( err_msg != "" ) std::cerr << err_msg << "\n" << std::endl;
  std::ostream& s = std::cout;
  s << "Usage: " << name << "\t [OPTIONS]  test.root reference.root    chain1 chain2 chain2 ...\n\n"; 
  s << "\t" << " plots comparison histograms"; 
  s << " - compiled on " << __DATE__ << " at " << __TIME__ << "\n\n";
  s << "Options: \n";
  s << "    -c,  --config       value \t configure which histograms to plot from config file,\n\n";
  s << "    -t,  --tag          value \t appends tag 'value' to the end of output plot names, \n";
  s << "    -k,  --key          value \t prepends key 'value' to the front of output plots name, \n";
  s << "    -t,  --tag          value \t post pend tag 'value' to the end of output plots name, \n";
  s << "    -d,  --dir          value \t creates output files into directory, \"value\" \n\n";
  s << "         --ncols        value \t creates panels with \"value\" columns\n\n";
 
  s << "    -e,  --efficiencies       \t make test efficiencies with respect to reference \n";
  s << "    -es, --effscale     value \t scale efficiencies to value\n";
  s << "    -er, --effscaleref  value \t scale reference efficiencies to value\n";
  s << "    -nb  --nobayes            \t do not calculate Basyesian efficiency uncertaintiesr\n\n";
 
  s << "    -r,  --refit               \t refit all test resplots\n";
  s << "    -rr, --refitref            \t also refit all reference resplots\n";
  s << "         --oldrms              \t use fast rms95 when refitting resplots\n\n";
 
  s << "    -as, --atlasstyle          \t use ATLAS style\n";
  s << "    -l,  --labels        values\t use specified labels for key\n";
  s << "         --taglabels     values\t use specified additional labels \n";
  s << "    -al, --atlaslable    value \t set value for atlas label\n";
  s << "    -sx, --swapaxtitles  exp pattern\t swap the expression \"exp \" in the axis titles with \"pattern\"\n";
  s << "    -ac, --addchains           \t if possible, add chain names histogram labels \n\n";   
 
  s << "    -m,  --mapfile            \t remap file for reference histograms \n\n";   
 
  s << "    -rc, --refchains values ..\t allow different reference chains for comparison\n";
  s << "    -s,  --swap  pattern regex \t swap \"pattern\" in the chain names by \"regex\"\n";
  s << "         --swapr pattern regex \t swap \"pattern\" in the ref chain names by \"regex\"\n";
  s << "         --swapt pattern regex \t swap \"pattern\" in the test chain names by \"regex\"\n";
  s << "    -nr, --noref              \t do not plot reference histograms\n";
  s << "         --normref            \t normalise the reference counting histograms to test histograms\n";
  s << "    -us, --usechainref        \t use the histograms from chain definied in the \"Chain\" histogram as reference\n\n";
  s << "         --lumi               \t scale pile-up mu values down by 5.4%\n";
  s << "         --lumitest           \t scale pile-up mu values down by 5.4% for test histograms\n";
  s << "         --lumiref            \t ecale pile-up mu values down by 5.4% for reference histograms\n";
  s << "    -rb, --refbands           \t draw the reference as an error band, rather than a histogram\n\n";
  
  s << "    -ns, --nostats            \t do not show stats for mean and rms\n";
  s << "    -nm, --nomeans            \t do not show stats for the means\n";
  s << "         --means              \t show stats for the means\n";
  s << "         --chi2               \t show the chi2 with respect to the reference\n\n";
 
  s << "    -np, --noplots            \t do not actually make any plot\n";
  s << "    -q,  --quiet              \t make the plots but do not print them out\n\n";
 
  s << "         --unscalepix         \t do not scale the number of pixels by 0.5 (unscaled by default)\n";
  s << "         --scalepix           \t do scale the number of pixels by 0.5 (unscaled by default)\n";
  s << "         --yrange     min max \t use specified y axis range\n";  
  s << "    -xo, --xoffset      value \t relative x offset for the key\n"; 
  s << "    -yp, --ypos         value \t relative yposition for the key\n"; 
  s << "    -xe, --xerror       value \t size of the x error tick marks\n"; 
  s << "    -nw, --nowatermark        \t do not plot the release watermark\n\n"; 
 
  s << "    -C,  --Cfiles             \t write C files also\n"; 
  s << "         --eps                \t create eps files\n"; 
  s << "         --nopng              \t do not print png files\n"; 
  s << "         --nopdf              \t do not print pdf files\n"; 
  s << "         --deleteref          \t delete unused reference histograms\n\n";
 
  s << "         --printconfig        \t print the configuration being used in the form useable as an input file\n\n";
 
  s << "         --run          value \t print the run number on the plot\n\n";
 
  s << "    -h,  --help              \t this help\n";
  //  s << "\nSee " << PACKAGE_URL << " for more details\n"; 
  //  s << "\nReport bugs to <" << PACKAGE_BUGREPORT << ">";
  s << std::endl;
  return status;
}
 
 
void binwidth( TH1F* h ) { 
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    double w = h->GetBinLowEdge(i+1) - h->GetBinLowEdge(i);
    h->SetBinContent( i, h->GetBinContent(i)/w );
    h->SetBinError( i, h->GetBinError(i)/w );
  }
}
 
void ascale( TH1F* h, double s_ ) { 
  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
    h->SetBinContent( i, h->GetBinContent(i)*s_ );
    h->SetBinError( i, h->GetBinError(i)*s_ );
  }
}
 
 
 
// replace from a string
std::string fullreplace( std::string s, const std::string& s2, const std::string& s3) {
  if ( s2=="" || s2==s3 ) return s;  
  std::string::size_type pos;
  while ( (pos=s.find(s2)) != std::string::npos )  s.replace(pos, s2.size(), s3);
  return s;
} 
 
void print_pad( const std::string& s ) { 
  std::cout << green << "Printing " << s << reset << std::endl;
  gPad->Print( s.c_str() );
  if ( EPS ) gPad->Print( (s+".eps").c_str() );
}
 
 
TH1F* makeplot( TObject* hreft ) {
 
  TH1F* href = 0; // (TH1F*)hreft->Clone();
      
  if ( isTEfficiency( hreft ) ) { 
    
    TH1F* hn = (TH1F*)dynamic_cast<TEfficiency*>(hreft)->GetPassedHistogram();
    TH1F* hd = (TH1F*)dynamic_cast<TEfficiency*>(hreft)->GetTotalHistogram();
    // TH1F* hn = (TH1F*)((TEfficiency*)hreft)->GetPassedHistogram();
    // TH1F* hd = (TH1F*)((TEfficiency*)hreft)->GetTotalHistogram();
 
    //    hd->DrawCopy();
    //    gPad->SetLogy(true);
    //    gPad->Print("den.pdf");
 
    //    hn->DrawCopy("same");
    //    gPad->SetLogy(true);
    //    gPad->Print("num.pdf");
 
    href = (TH1F*)hn->Clone("ars"); href->SetDirectory(0);
    for ( int ie=1 ; ie<=href->GetNbinsX() ; ie++ ) {
      if ( hd->GetBinContent(ie)==0 ) { href->SetBinContent( ie, 0 ); href->SetBinError( ie, 0 ); }
      else {
    double e = hn->GetBinContent(ie)/hd->GetBinContent(ie);
    double ee = e*(1-e)/hn->GetBinContent(ie);
    //if (e>1) e=1;
    href->SetBinContent( ie, e );
    href->SetBinError( ie, ee ); 
      }
    }
  }
 
  // if ( href) href->DrawCopy();
  // gPad->SetLogy(false);
  // gPad->Print("plot.pdf");
  
  return href;
}
 
 
 

void zero( TH2* h ) { 
  int Nx = h->GetNbinsX();
  int Ny = h->GetNbinsY();
  for ( int i=1 ; i<=Nx ; i++ ) { 
    for ( int j=1 ; j<=Ny ; j++ ) { 
      if ( h->GetBinContent( i, j )==0 ) { 
    h->SetBinContent( i, j, 1e-10 );
    h->SetBinError( i, j, 1e-10 );
      }
    }
  }
}
 

void normy( TH2* h ) { 
  int Nx = h->GetNbinsX();
  int Ny = h->GetNbinsY();
  for ( int i=1 ; i<=Nx ; i++ ) { 
    double n = 0;
    for ( int j=1 ; j<=Ny ; j++ ) n += h->GetBinContent(i,j);
    if ( n>0 ) for ( int j=1 ; j<=Ny ; j++ ) h->SetBinContent( i, j, h->GetBinContent(i,j)/n );
  }
}
 

void autox( TH2* h, bool t=false ) {
 
  int Nx = h->GetNbinsX();
  int Ny = h->GetNbinsY();
 
  double xhi = h->GetXaxis()->GetBinLowEdge(1);
  double xlo = h->GetXaxis()->GetBinLowEdge( h->GetXaxis()->GetNbins() );
 
  for ( int i=1 ; i<=Nx ; i++ ) { 
    double x = h->GetXaxis()->GetBinLowEdge(i);
    for ( int j=1 ; j<=Ny ; j++ ) { 
      double n = h->GetBinContent(i,j);
      if ( n!=0 ) { 
    if ( x<xlo ) { 
      if ( i>1 ) xlo = h->GetXaxis()->GetBinLowEdge(i-1);
      else  xlo = x;
    }
    if ( x>xhi ) { 
      if ( i<h->GetXaxis()->GetNbins() ) xhi = h->GetXaxis()->GetBinLowEdge(i+1);
      else  xhi = x;
    }
      }
    }
  }
  
  if ( t ) { 
    xlo -= 5.4;
    xhi -= 5.4;
  }
 
  h->GetXaxis()->SetRangeUser( xlo, xhi );
}
 

void autoy( TH2* h ) {
 
  int Nx = h->GetNbinsX();
  int Ny = h->GetNbinsY();
 
  double yhi = h->GetYaxis()->GetBinLowEdge(1);
  double ylo = h->GetYaxis()->GetBinLowEdge( h->GetYaxis()->GetNbins() );
 
  for ( int i=1 ; i<=Ny ; i++ ) { 
    double y = h->GetYaxis()->GetBinLowEdge(i);
    for ( int j=1 ; j<=Nx ; j++ ) { 
      double n = h->GetBinContent(j,i);
      if ( n!=0 ) { 
    if ( y<ylo ) { 
      if ( i>1 ) ylo = h->GetYaxis()->GetBinLowEdge(i-1);
      else  ylo = y;
    }
    if ( y>yhi ) { 
      if ( i<h->GetYaxis()->GetNbins() ) yhi = h->GetYaxis()->GetBinLowEdge(i+1);
      else  yhi = y;
    }
      }
    }
  }
  
  h->GetYaxis()->SetRangeUser( ylo, yhi );
}
 
 
 
double chi2( TH1* h0, TH1* h1 ) {
 
  std::cout << "chi2" << std::endl;
 
  double c2 = 0;
 
  std::cout << "h: " << h0 << " " << h1 << std::endl; 
 
  if ( h0==0 || h1==0 ) {
    std::cerr << "histogram not found: h0: " << h0 << "\th1: " << h1 << std::endl;
    return 999;
  }
  
  std::cout << "h: " << h0->GetName() << std::endl; 
  std::cout << "h: " << h1->GetName() << std::endl; 
  
  std::cout << "h0->GetNbinsX() " << h0->GetNbinsX() << std::endl;
  std::cout << "h0->GetNbinsX() " << h1->GetNbinsX() << std::endl;
 
  for ( int i=0 ; i<h0->GetNbinsX() ; i++ ) {
 
    double d0 = h0->GetBinContent(i+1);
    double d1 = h1->GetBinContent(i+1);
 
    double e0 = h0->GetBinError(i+1);
    double e1 = h1->GetBinError(i+1);
 
    double  e2 = e0*e0+e1*e1;
    
    if ( e2>0 ) c2 += (d0-d1)*(d0-d1)/e2; 
    
  }
  
  return c2;
}
 
void ifdbg( const std::string& /* s */ ) {
  //  std::cout << "if dbg: " << s << std::endl;
}
  
 
 
 
int main(int argc, char** argv) { 
 
  gErrorIgnoreLevel = kError;
 
  std::cout << black;
  std::cout << "\n---------------------------------\n";
  std::cout << "\n comparitor is off ...\n";
 
  if ( argc<4 ) return usage(argv[0], -1, "too few arguments");

  
  std::string tag = "";
  std::string key = "";
 
  std::string dir = "";
 
  std::string ftestname = "";
  std::string frefname  = "";
 
  TFile* ftest_ = 0;
  TFile* fref_  = 0;
 
  bool   effset = false;
  double effmin =  90;
  double effmax = 102;
 
  std::string defreflabel = "";

 
  std::vector<std::string> usrlabels;
  bool uselabels    = false;
  bool addinglabels = false;
 
  std::vector<std::string> taglabels;
  std::vector<std::string> taglabels2;
  bool addingtags  = false;
  bool addingtags2 = false;
 
  bool lumiref_trans  = false;
  bool lumitest_trans = false;
 
  bool make_efficiencies = true;
  
  bool make_ref_efficiencies = false;
  bool refit_resplots        = false;
  bool refitref_resplots     = false;
  bool bayes      = true;
  bool nostats     = false;
  bool nomeans     = false;
  bool means       = false;
  bool noref       = false;
  bool atlasstyle  = false;
  bool deleteref   = false;
  bool nowatermark = false;
  bool noplots     = false;
  bool nopng       = false;
  bool nopdf       = false;
  bool Cfile       = false;
  bool notitle     = true;
  bool dochi2      = false;
  bool normref     = false;
  bool scalepix    = false;
  bool oldrms      = false;
  bool addchains   = false;
  bool usechainref = false;
  bool quiet       = false;
 
  bool RANGEMAP    = false;
  bool ALLRANGEMAP = false;
 
  double xerror    = 0;
 
  std::string atlaslabel_tmp = "Internal";
 
 
  double scale_eff     = -1;
  double scale_eff_ref = -1;
  
  std::string configfile = "";
 
  double xoffset = 0;
 
  
  double ypos_in = 0;
 
  std::string pattern = "";
  std::string regex   = "";
 
  std::string patternr = "";
  std::string regexr   = "";
 
  std::string patternt = "";
  std::string regext   = "";
 
  std::string basedir = "";
 
  std::string xpattern = "";
  std::string xregex   = "";
 
  std::vector<std::string> chains;
  std::vector<std::string> refchains;
 
  bool addingrefchains = false;
 
  std::string mapfile = "";
 
  std::vector<std::string> chainfiles;
 
  int ncols = 2;
 
  std::string runnumber = "";
 
  bool refbands = false;
  
  for(int i=1; i<argc; i++){
    std::string arg  = argv[i];
 
    ifdbg( "\ncnt: " + arg );
 
    
    if ( arg.find('-')!=0 && addinglabels ) {
      std::string label = arg;
      fullreplace( label, "__", " " );
      replace( label, '#', ' ' );
      usrlabels.push_back( label );
      continue;
    }
    else addinglabels = false;
 
    if ( arg.find('-')!=0 && addingrefchains ) { 
      refchains.push_back( arg );
      continue;
    }
    else addingrefchains = false;
 
    if ( arg.find('-')!=0 && addingtags ) { 
      taglabels.push_back( fullreplace( arg, "__", " " ) );
      std::cout << "\tadding tag label: " << taglabels.back() << std::endl;
      continue;
    }
    else addingtags = false;
 
    if ( arg.find('-')!=0 && addingtags2 ) { 
      taglabels2.push_back( fullreplace( arg, "__", " " ) );
      std::cout << "\tadding tag label: " << taglabels2.back() << std::endl;
      continue;
    }
    else addingtags2 = false;
 
    ifdbg( "starting if cascade: "+arg );
 
    if ( arg=="-h" || arg=="--help" ) { 
      return usage(argv[0], 0);
    }
    else if ( arg=="-c" || arg=="--config" ) { 
      if ( ++i<argc ) configfile=argv[i];
      else return usage(argv[0], -1, "no config file provided");
      ifdbg("-c");
    }
    else if ( arg=="--jl" ) { 
      JLflag = true;
      ifdbg("--jl");
    }
    else if ( arg=="--eps" ) { 
      EPS = true;
      ifdbg("--eps");
    }
    else if ( arg=="-t" || arg=="--tag" ) { 
      if ( ++i<argc ) tag=std::string("-")+argv[i];
      else return usage(argv[0], -1, "no tag provided");
      ifdbg("--tag");
    }
    else if ( arg=="-l" || arg=="--labels" ) { 
      addinglabels = true;
      ifdbg("--labels");
    }
    else if ( arg=="-el" ) { 
      LINES = true;
      ifdbg("-el");
    }
    else if ( arg=="-k" || arg=="--key" ) { 
      if ( ++i<argc ) key=argv[i];
      else return usage(argv[0], -1, "no key provided");
      ifdbg("--key");
    }
    else if ( arg=="--run" ) { 
      if ( ++i<argc ) runnumber=argv[i];
      else return usage(argv[0], -1, "no run number provided");
      ifdbg("--run");
    }
    else if ( arg=="-m" || arg=="--mapfile" ) { 
      if ( ++i<argc ) mapfile=argv[i];
      else return usage(argv[0], -1, "no mapfile provided");
      ifdbg("--mapfile");
    }
    else if ( arg=="-d" || arg=="--dir" ) { 
      if ( ++i<argc ) dir=argv[i];
      else return usage(argv[0], -1, "no directory provided");
      ifdbg("--dir");
    }
    else if ( arg=="--lumi" ) { 
      lumitest_trans = true;
      lumiref_trans = true;
      ifdbg("--lumi");
    }
    else if ( arg=="--lumitest" ) { 
      lumitest_trans = true;
      ifdbg("--lumitest");
    }
    else if ( arg=="--lumiref" ) { 
      lumiref_trans = true;
      ifdbg("--lumiref");
    }
    else if ( arg=="-rb" || arg=="--refbands" ) { 
      refbands = true;
      ifdbg("--refbands");
    }
    else if ( arg=="-b" || arg=="--bdir" ) { 
      if ( ++i<argc ) basedir=argv[i];
      else return usage(argv[0], -1, "no directory provided");
      ifdbg("--bdir");
    }
    else if ( arg=="--taglabels" ) { 
      addingtags = true;
      ifdbg("--taglabels");
    }
    else if ( arg=="--taglabels2" ) { 
      addingtags2 = true;
      ifdbg("--taglabels2");
    }
    else if ( arg=="--unscalepix" ) { 
      scalepix = false;
      ifdbg("--unscalepix");
    }
    else if ( arg=="--scalepix" ) { 
      scalepix = true;
      ifdbg("--scalepix");
    }
    else if ( arg=="-ac" || arg=="--addchains" ) { 
      addchains = true;
      ifdbg("--sddchains");
    }
    else if ( arg=="-yrange" ) { 
      effset = true;
      if ( ++i<argc ) effmin=std::atof(argv[i]);
      else return usage(argv[0], -1, "no range specified");
      if ( ++i<argc ) effmax=std::atof(argv[i]);
      else return usage(argv[0], -1, "no upper y limit specified");
      ifdbg("-yrange");
    }
    else if ( arg=="-e" || arg=="--efficiencies" ) { 
      make_ref_efficiencies = true;
      ifdbg("--efficiencies");
    }
    else if ( arg=="-r" || arg=="--refit" ) { 
      refit_resplots = true;
      ifdbg("--refit");
    }
    else if ( arg=="-rr" || arg=="--refitref" ) { 
      refitref_resplots = true;
      ifdbg("--refitref");
    }
    else if ( arg=="--oldrms" ) { 
      oldrms = true;
      ifdbg("--oldrms");
    }
    else if ( arg=="-nw" || arg=="--nowatermark" ) {
      nowatermark = true;
      Plots::setwatermark(!nowatermark);
      ifdbg("--nowatermark");
    }
    else if ( arg=="--chi2" ) { 
      dochi2 = true;
      ifdbg("--chi2");
    }
    else if ( arg=="-ns" || arg=="--nostats" ) { 
      nostats = true;
      ifdbg("--nostats");
    }
    else if ( arg=="-nm" || arg=="--nomeans" ) { 
      nomeans = true;
      ifdbg("--nomeans");
    }
    else if ( arg=="--means" ) { 
      means = true;
      ifdbg("--means");
    }
    else if ( arg=="-nt" || arg=="--notitle" ) { 
      notitle = true;
      ifdbg("--notitle");
    }
    else if ( arg=="-nr" || arg=="--noref" ) { 
      Plotter::setplotref(false);
      noref = true;
      ifdbg("--noref");
    }
    else if ( arg=="--normref" ) { 
      normref = true;
      ifdbg("--normref");
    }
    else if ( arg=="-rc" || arg=="--refchains" ) { 
      addingrefchains = true;
      ifdbg("--refchains");
    }
    else if ( arg=="-uc" || arg=="--usechainref" ) { 
      usechainref = true;
      ifdbg("--usechainref");
    }
    else if ( arg=="-nb" || arg=="--nobayes" ) { 
      bayes = false;
      ifdbg("--nobayes");
    }
    else if ( arg=="-es" || arg=="--effscale" ) { 
      if ( ++i<argc ) scale_eff=std::atof(argv[i]);
      else return usage(argv[0], -1, "no efficiency scale provided");
      ifdbg("--effscale");
    } 
    else if ( arg=="-er" || arg=="--effscaleref" ) { 
      if ( ++i<argc ) scale_eff_ref=std::atof(argv[i]);
      else return usage(argv[0], -1, "no efficiency scale for the reference histograms provided");
      ifdbg("--effscaleref");
    } 
    else if (               arg=="--ncols" ) { 
      if ( ++i<argc ) ncols=std::atoi(argv[i]);
      else return usage(argv[0], -1, "no number of columns provided");
      ifdbg("--ncols");
    } 
    else if ( arg=="-np" || arg=="--noplots" ) { 
      noplots = true;
      ifdbg("--noplots");
    }
    else if ( arg=="-C" || arg=="--Cfiles" ) { 
      Cfile = true;
      ifdbg("--Cfiles");
    }
    else if ( arg=="--deleteref" ) { 
      deleteref = true;
      ifdbg("--deleteref");
    }
    else if (              arg=="--nopng" ) { 
      nopng = true;
      ifdbg("--nopng");
    }
    else if (              arg=="--nopdf" ) { 
      nopdf = true;
      ifdbg("--nopdf");
    }
    else if ( arg=="-as" || arg=="--atlasstyle" ) { 
      atlasstyle = true;
      //      std::cout << "atlas style" << std::endl;      
      ifdbg("--atlasstyle");
    }
    else if ( arg=="--tp" ) { 
      make_efficiencies = false;
      ifdbg("--tp");
    }
    else if ( arg=="-q" || arg=="--quiet" ) { 
      quiet = true;
      ifdbg("--quiet");     
    } 
    else if ( arg=="-al" || arg=="--atlaslabel" ) { 
      if ( ++i<argc ) atlaslabel_tmp=argv[i];
      else return usage(argv[0], -1, "no label provided");
      ifdbg("--atlaslabel");
    }
    else if ( arg=="-xo" || arg=="--xoffset" ) { 
      if ( ++i<argc ) xoffset=std::atof(argv[i]);
      else return usage(argv[0], -1, "no xoffset provided");
      ifdbg("--xoffset");       
    }
    else if ( arg=="-yp" || arg=="--ypos" ) { 
      if ( ++i<argc ) ypos_in=std::atof(argv[i]);
      else return usage(argv[0], -1, "no y position provided");
      ifdbg("--ypos");      
    }
    else if ( arg=="-xe" || arg=="--xerror" ) { 
      if ( ++i<argc ) xerror=std::atof(argv[i]);
      else return usage(argv[0], -1, "no x error provided");
      ifdbg("--xerror");        
    }
    else if ( arg=="-s" || arg=="--swap" ) { 
      if ( ++i<argc ) pattern=argv[i];
      else return usage(argv[0], -1, "no patterns provided");
      if ( ++i<argc ) regex=argv[i];
      else return usage(argv[0], -1, "no target pattern provided");
      ifdbg("--swap");      
    }
    else if (             arg=="--swapt" ) { 
      if ( ++i<argc ) patternt=argv[i];
      else return usage(argv[0], -1, "no patterns provided");
      if ( ++i<argc ) regext=argv[i];
      else return usage(argv[0], -1, "no target pattern provided");
      ifdbg("--swapt");     
    }
    else if (             arg=="--swapr" ) { 
      if ( ++i<argc ) patternr=argv[i];
      else return usage(argv[0], -1, "no patterns provided");
      if ( ++i<argc ) regexr=argv[i];
      else return usage(argv[0], -1, "no target pattern provided");
      ifdbg("--swapr");     
    }
    else if ( arg=="-sx" || arg=="--swapaxtitles" ) { 
      if ( ++i<argc ) xregex=argv[i];
      else return usage(argv[0], -1, "no target pattern provided");
      if ( ++i<argc ) xpattern=argv[i];
      else return usage(argv[0], -1, "no patterns provided");
      ifdbg("--swapxtitles");       
    }
    else if ( arg.find('-')==0 ) {
      std::cerr << "unknown option: " << arg << "\n" << std::endl;
      return usage(argv[0], -4);
    }
    else {
      ifdbg("arg: + "+arg);     
      if      ( ftestname=="" ) ftestname = arg;
      else if (  frefname=="" ) frefname = arg;
      else { 
    std::string file = "";
    
    if ( arg.find(":file:")!=std::string::npos ) { 
      file = arg.substr( 0, arg.find(":file:") );
      chainfiles.push_back( file );
      arg  = arg.substr(arg.find(":file:")+6,arg.size());
    }
 
    std::string chain = arg;
 
    replace ( chain, ':', '_' );
    replace ( chain, ';', '_' );
    chains.push_back(chain);
 
    std::cout << "file: " << file << "\tchain: " << chain << std::endl;
 
      }
    }
  }
 
  if ( ftestname.empty() )  { 
    std::cerr << "main(): test file not specified " << std::endl;
    return -1;
  }
 
  if ( !exists(ftestname) ) { 
    std::cerr << "main(): test file " << ftestname << " does not exist" << std::endl;
    return -1;
  }
 
  std::vector<TFile*> chainTFiles;
 
  if ( chainfiles.size()==0 )  ftest_ = TFile::Open( ftestname.c_str() );
  else {
    noref=true;
    chainTFiles.resize(chainfiles.size());
    for ( size_t i=0 ; i<chainfiles.size() ; i++ ) { 
      chainTFiles[i] = TFile::Open( chainfiles[i].c_str() );
      if ( chainTFiles[i] == 0 ) { 
    std::cerr << "\tfile: " << chainfiles[i] << " could not be opened" << std::endl;
    return -1;
      }
      else std::cout << "\tchainfiles: " << chainfiles[i] << " " << chainTFiles[i] << std::endl;
    }
  }
 
 
  if ( noref==false ) { 
    if ( frefname.empty() )  { 
      std::cerr << "main(): ref file not specified " << std::endl;
      Plotter::setplotref(false);
      noref = true;
      //       return -1;
    }
    
    if ( frefname==ftestname )    fref_ = ftest_;
    else if ( exists(frefname) )  fref_ = TFile::Open( frefname.c_str() );
    else { 
      std::cerr << "main(): ref file " << frefname << " does not exist" << std::endl;
      //  return -1;
      Plotter::setplotref(false);
      noref=true;
    }
  }
  else fref_ = ftest_;
  
  if ( chainfiles.size()==0 ) { 
    if ( ftest_==0 ) { 
      std::cerr << "could not open test file " << ftestname << std::endl;
      return -1;
    }
 
    if ( noref==false && fref_==0 ) { 
      std::cerr << "could not open files " << std::endl;
      noref=true;
      fref_=ftest_;
      defreflabel = "failed to open reference file";
    }
    
  }
  
  if ( scale_eff     == -1 ) scale_eff     = 100;
  if ( scale_eff_ref == -1 ) scale_eff_ref = scale_eff;
  
 
  bool noreftmp = noref;
 
  if ( chains.size()==0 ) return usage(argv[0], -1, "no chains specified");
 
 
  if ( basedir.size()>0 ) {
    if ( basedir[basedir.size()-1]!='/' ) basedir += "/";
    for ( size_t ic=chains.size() ; ic-- ; ) chains[ic] = basedir+chains[ic];
  }
 
  if ( refchains.size()>0 && refchains.size()!=chains.size() ) return usage(argv[0], -1, "not enough chains specified");
  
  if ( refchains.size()==0 ) refchains = chains;
  
  std::vector<std::string> chainref(chains.size(),"");
  std::vector<std::string> chain_name(chains.size(),"");
 
  std::vector<std::string> refchain(chainref.size(),"");
    
 
 
  std::cout << argv[0] << " options:" << std::endl;
  std::cout << "\tATLAS style:                 " << ( atlasstyle ? "true" : "false" ) << std::endl; 
  std::cout << "\tBayesian uncertainties:      " << ( bayes ? "true" : "false" ) << std::endl; 
  std::cout << "\trefit resplot uncertainties: " << ( refit_resplots ? "true" : "false" ) << std::endl; 
  std::cout << "\tsuppress mean and rms stats: " << ( nostats ? "true" : "false" ) << std::endl;  
  if ( !nostats ) std::cout << "\tsuppress meanstats:          " << ( nomeans ? "true" : "false" ) << std::endl;  
  std::cout << "\tsuppress png output:         " << ( nopng ? "true" : "false" ) << std::endl;  
  std::cout << "\tsuppress pdf output:         " << ( nopdf ? "true" : "false" ) << std::endl;  
  std::cout << "\tsuppress reference output:   " << ( noref ? "true" : "false" ) << std::endl;  
  std::cout << "\tuse chain references:        " << ( usechainref ? "true" : "false" ) << std::endl;
  std::cout << "\tpanel ncols:                 " << ncols << std::endl;
 
  if ( usrlabels.size()>0 )  std::cout << "\tlabels:                      " << usrlabels.size() << std::endl;        
  if ( taglabels.size()>0 )  std::cout << "\textra text:                  " << taglabels  << std::endl;  
  if ( taglabels2.size()>0 ) std::cout << "\textra text:                  " << taglabels2 << std::endl;  
  
 
  for ( size_t il=0 ; il<usrlabels.size() ; il++ ) { 
    std::cout << "usr label[" << il << "] : " << usrlabels[il] << std::endl;
  } 
 
  std::cout << "atlas style : " << atlasstyle << std::endl;
 
  if ( atlasstyle ) { 
    SetAtlasStyle();
    gStyle = AtlasStyle();
    gStyle->cd();
  }
  else { 
    gROOT->SetStyle("Plain");
    gStyle->cd();
    gStyle->SetLineScalePS(1);
    xoffset += 0.02;
  }
 
  gStyle->SetErrorX(xerror);
 
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadBottomMargin(0.15);
 
  gStyle->SetPadRightMargin(0.02);
  gStyle->SetPadTopMargin(0.05);
 
  std::cout << "Chains: " << std::endl;
  for ( unsigned ic=0 ; ic<chains.size() ; ic++ ) std::cout << "\t" << chains[ic] << std::endl;    
 
  if ( usrlabels.size()>0 ) std::cout << "labels: " << usrlabels << std::endl;
 
  if ( usrlabels.size()>0 && usrlabels.size()==chains.size() ) uselabels = true;

 
  TTree*   dataTree = 0;
  TString* releaseData = new TString("");
  std::vector<std::string> release_data;
 
 
  if ( !nowatermark && ftest_ ) { 
 
    dataTree = (TTree*)ftest_->Get("dataTree");
  
    if ( dataTree ) { 
      dataTree->SetBranchAddress( "ReleaseMetaData", &releaseData);
      
      for (unsigned int i=0; i<dataTree->GetEntries() ; i++ ) {
    dataTree->GetEntry(i);      
    release_data.push_back( releaseData->Data() );
    std::cout << "main() release data: " << release_data.back() << " : " << *releaseData << std::endl;
      }
    }
 
    std::cout << "release data size " << release_data.size() << std::endl;
 
    if ( release_data.size()>0 ) { 
      if ( release_data.size()>1 ) std::cerr << "main() more than one release - using only the first" << std::endl;  
      
      //    std::cout << "release: " << chop(release_data[0], " " ) << std::endl;
 
      std::string rd = release_data[0];
      
      while ( rd.size()>0 ) {
    
    std::string nightly = chop(rd, " " );
 
    if ( contains(nightly,"private" ) ) { 
      for ( int ic=0 ; ic<4 ; ic++ ) chop(release_data[0], " " );
      release += "  (" + release_data[0]+")";
      break;
    }
    else {
 
      if ( contains(nightly,"nightly") ) {
        chop( nightly,"/" );        
        release += "  ( " + chop( nightly,"/"); 
        release += " " + chop(nightly, "/" ) + " )";
        break;
      }
      
    }
      }
    }
  }
    
  std::cout << "release: " << release << std::endl;
  
 
  // Make output directory                                                                                                                           
  if (!dir.empty()) {
    
    gDirectory->pwd();
    
    std::cout << "trying to make directory" << std::endl;
    dir += '/';
    if ( !quiet && !noplots && !exists(dir) ) { 
      if ( mkdir( dir.c_str(), 0777 ) ) std::cerr << "main() couldn't create directory " << dir << std::endl;
      else                              std::cout << "main() output will be sent to directory " << dir << std::endl; 
    }
  }
 
  if ( ftest_==0 && chainTFiles.size()>0 ) fref_ = ftest_ = chainTFiles[0]; 
 
  TFile* ftest = ftest_;
  TFile* fref  = fref_;
 
  std::string testrun = findrun( ftest );
 
  std::string  rawrefrun = "";
  std::string  refrun = "";
 
  if ( fref_ ) refrun = rawrefrun = findrun( fref  );
 
  std::cout << "testrun: " << testrun << "\nrefrun:  " << refrun << std::endl;
 
  if ( !testrun.empty() && refrun != testrun ) { 
    if ( pattern.empty() ) { 
      pattern = testrun;
      regex   =  refrun;
    }
 
  }
  
  if ( !refrun.empty() ) { 
 
    std::string newtag = "Reference: ";
 
    std::cout << "refrun: " << refrun << std::endl;
 
    size_t pos;
    while ( (pos=refrun.find('_'))!=std::string::npos ) refrun.replace( pos, 1, " " );
    newtag += refrun;
 
    std::string rawrun = refrun.erase( refrun.find("run"), 4 );
 
    if ( contains(frefname, rawrun) ) { 
      
      std::string release = frefname;
 
      release.erase( 0, release.find(rawrun) ); 
   
      if    ( contains(release,"HIST") ) release.erase( 0, release.find("HIST")+5 ); 
      while ( contains(release,".") ) release.erase( release.find('.'), release.size() ); 
      while ( contains(release,"-") ) release.erase( release.find('-'), release.size() ); 
      while ( contains(release,"_p") ) release.erase( release.find("_p"), release.size() ); 
      while ( contains(release,"_t") ) release.erase( release.find("_t"), release.size() ); 
 
      newtag += " ";
      newtag += release;
 
    } 
 
    taglabels.push_back( newtag );
    
    std::cout << "tag labels: " << taglabels << std::endl;
 
  } 
  
  int NeventTest = 1;
  int NeventRef  = 1;
 
 
  std::vector<std::string> savedhistos;
 
  
  if ( !nowatermark ) { 
 
    
    TH1D* htestev = (TH1D*)ftest->Get("event") ;
    TH1D* hrefev  = (TH1D*)fref->Get("event") ;
 
 
    
    std::cout << "htestev " << htestev << " " << hrefev  << std::endl;
    
    if ( htestev ) NeventTest = htestev->GetEntries(); 
    if ( hrefev )  NeventRef  = hrefev->GetEntries();
  
    savedhistos.push_back("event");
  }
  else { 
    NeventTest = 1; 
    NeventRef  = 1; 
  }
 
 
  if ( NeventTest>1 ) std::cout << "Nevents Test: " << NeventTest << std::endl;
  if ( NeventRef>1 )  std::cout << "Nevents Ref:  " << NeventRef << std::endl;
 
  chainmap_t* chainmap = nullptr;
 
  if ( mapfile.empty() ) mapfile = configfile;
 
 
  std::cout << "mapfile: " << mapfile << std::endl;
 
  if ( !mapfile.empty() ) {
 
    ReadCards m( mapfile );
 
    if ( m.isTagDefined( "ChainMap" ) )  {
 
      std::vector<std::string> chains = m.GetStringVector( "ChainMap" );
 
      for ( size_t i=0 ; i<chains.size() ; i+=2 ) std::cout << "map: " << chains[i] << " -> " << chains[i+1] << std::endl;
      
      chainmap = new chainmap_t();
      
      for ( size_t i=0 ; i<chains.size() ; i+=2 ) { 
    chainmap->insert( chainmap_t::value_type( chains[i], chains[i+1] ) );
      }
      
      std::cout << "\nusing chain map:" << std::endl;
      
      for ( chainmap_t::iterator itr=chainmap->begin() ; itr!=chainmap->end() ; itr++ ) { //cppcheck-suppress postfixOperator 
    std::cout << "\t" << itr->first << "\t" << itr->second << std::endl;
      }
    }
  }

 
  std::cout << "\ncreating chain and reference information ..." << std::endl;
 
  for ( size_t j=0; j<chains.size(); j++)  {
 
    if ( !regex.empty() )  chains[j] = fullreplace( chains[j], pattern, regex );
    if ( !regext.empty() ) chains[j] = fullreplace( chains[j], patternt, regext );
    

    std::cout << "chain:        " << chains[j] << "\taddchains: " << addchains << std::endl;
 
    if ( addchains && ( contains(chains[j],"Shifter") || ( !contains(chains[j],"HLT_") && !contains(chains[j], "Fullscan" ) ) ) ) { 
 
      TFile* fftest = ftest;
  
      if ( chainfiles.size()>0 && chainfiles.size()>j ) fftest = chainTFiles[j];
 
      TH1F* hchain = Get( *fftest, chains[j]+"/Chain", testrun );
      
      if ( hchain ) { 
    
    std::string name = hchain->GetTitle();
    
    if ( usechainref && !contains(chains[j],"Purity") ) { 
    
      chainref[j] = name;
      
      std::cout << "new chainref:   " << chainref[j] << std::endl;
 
      std::string::size_type pos = chainref[j].find(":for");
      if ( pos!=std::string::npos ) chainref[j].replace( pos, 4, "_for" );
      std::replace( chainref[j].begin(), chainref[j].end(), ':', '/');      
      std::string newchain = dirname( dirname( chains[j]) ) + "/Expert/" + chainref[j];
      chainref[j] = newchain;
 
      std::cout << "final chainref: " << chainref[j] << std::endl;
 
    }
 
    std::cout << "chainref:     " << chainref[j] << std::endl;
    
    while ( contains( name, "HLT_" ) ) name = name.erase( name.find("HLT_"), 4 );
 
    std::cout << "name:         " << name << std::endl;
 
    if ( contains( name, ":" ) )  chain_name[j] = name.substr( 0, name.find(':') ) + " : ";
    else                          chain_name[j] = name;
 
    if  ( chain_name[j] == " : "  )  chain_name[j] = "";
    
    std::cout << "chain_name: " << chain_name[j] << std::endl;
      }
    }
  }
 
  std::cout << "chainref size: " << chainref.size() << "  " << refchains.size() << std::endl;
 
  for ( size_t j=0 ; j<chainref.size() ; j++ ) { 
    std::cout << "chainref: " << chainref[j] << " :: " << refchains[j] << std::endl;

    if ( chainref[j]!="" ) refchain[j] = fullreplace( chainref[j], pattern, regex );
    else                   refchain[j] = fullreplace( refchains[j], pattern, regex );
 
    if ( !patternr.empty() ) { 
      if ( chainref[j]!="" ) refchain[j] = fullreplace( chainref[j], patternr, regexr );
      else                   refchain[j] = fullreplace( refchains[j], patternr, regexr );
    }
 
    std::cout << "refchain: " << refchain[j] << std::endl;
  }
 
  std::cout << "done chains" << std::endl;
  
 
  std::vector<Panel> panels;

 
  bands bnd;
 
 
  std::vector<int>         ccolours;
  std::vector<int>         cstyles;
 
  std::vector<std::string> ctags;
  std::vector<std::string> ctaglabels;
 
 
  std::cout << "\n" << argv[0] << "\tconfigfile: " << configfile << std::endl;
 
  bool use_file_config = false;
 
  if ( !configfile.empty() ) { 
 
    if ( exists(configfile) ) { 
 
      std::cout << argv[0] << ":\treading configuration file " << configfile << std::endl; 
    
      ReadCards rc(configfile);

      
      std::cout << "searching for panels" << std::endl;

 
      if ( rc.isTagDefined( "panels" ) ) { 
    
    std::cout << argv[0] << ":\treading histogram panel configuration from file " << configfile << std::endl; 
 
    use_file_config = true;
 
    std::vector<std::string> panel_config = rc.GetStringVector( "panels" );
    
    std::vector<string> panel_columns;
 
    if ( rc.isTagDefined("panel_columns") ) { 
      panel_columns = rc.GetStringVector( "panel_columns" );
      if ( panel_columns.size()%2 ) return usage( argv[0], -1, "incorrect panel settings" );  
    }
 
 
    for ( size_t ipanel=panel_config.size() ; ipanel-- ;  ) { 
      
      std::vector<std::string> raw_input = rc.GetStringVector( panel_config[ipanel] );
 
      int tncols = ncols;
      
      if ( panel_columns.size() ) {
        std::vector<string>::iterator itr = find( panel_columns.begin(), panel_columns.end(), panel_config[ipanel] );
        if ( itr!=panel_columns.end() ) tncols = std::atoi( (++itr)->c_str() );
      }
  
      Panel p( panel_config[ipanel], tncols );
      
      if ( raw_input.empty() ) {
        std::cerr << "no plots provided" << std::endl;
        return 1;
      }
 
      for ( size_t iraw=0 ; iraw+6<=raw_input.size() ; iraw += 6 ) p.push_back( raw_input.begin()+iraw, raw_input.begin()+iraw+6 );
          
      panels.push_back( p ); 
 
    }
 
      }
      else if ( rc.isTagDefined( "histos" ) ) { 
    
    std::cout << argv[0] << ":\treading histogram configuration from file " << configfile << std::endl; 
    
    use_file_config = true;
    
    std::vector<std::string> raw_input = rc.GetStringVector( "histos" );
    
    for ( size_t iraw=0 ; iraw+6<=raw_input.size() ; iraw += 6) {
      HistDetails h( raw_input.begin()+iraw, raw_input.begin()+iraw+6 );
      Panel p( h.name(), 1 ); 
      p.push_back( h );
      panels.push_back( p );
    }
          
      }
      
      if ( rc.isTagDefined( "Bands" ) && rc.isTagDefined( "Labels" ) ) { 
    bnd = bands( rc.GetVector( "Bands"), rc.GetStringVector( "Labels" ) );
      }
 
 
      if ( rc.isTagDefined("Colours") ) { 
 
    std::vector<int> ccolours = rc.GetIntVector("Colours");
 
    std::cout << "remap colours: " << ccolours << std::endl;
 
    int new_colours[6];
 
    for ( size_t i=6 ; i-- ; ) new_colours[i] = colours[i];
 
    for ( size_t i=0 ; i<ccolours.size() && i<6 ; i++ ) {
      if ( ccolours[i]<6 ) new_colours[i] = colours[ccolours[i]];
    }
 
    for ( size_t i=6 ; i-- ; ) colours[i] = new_colours[i];
 
      }
 
 
      if ( rc.isTagDefined("Styles") ) { 
      
    std::vector<int> cstyles = rc.GetIntVector("Styles");
 
    std::cout << "remap markers: " << cstyles << std::endl;
 
    int new_markers[6];
 
    for ( size_t i=6 ; i-- ; ) new_markers[i] = markers[i];
     
    for ( size_t i=0 ; i<cstyles.size() && i<6 ; i++ ) {
      if ( cstyles[i]<6 ) new_markers[i] = markers[cstyles[i]];
    }
    
    for ( size_t i=6 ; i-- ; ) markers[i] = new_markers[i];
 
      }
 
      if ( fulldbg ) { 
    for ( int i=6 ; i-- ; ) std::cout << "\tcolours[" << i << "] = " << colours[i] << std::endl;
    for ( int i=6 ; i-- ; ) std::cout << "\tmarkers[" << i << "] = " << markers[i] << std::endl;
      }
 
      if ( rc.isTagDefined("Tags") )      ctags      = rc.GetStringVector("Tags");
      if ( rc.isTagDefined("TagLabels") ) ctaglabels = rc.GetStringVector("TagLabels");
      if ( rc.isTagDefined("TagLabels") ) usrlabels  = rc.GetStringVector("TagLabels");
    
 
      if ( rc.isTagDefined("Styles") )  { 
        for ( size_t is=0 ; is<cstyles.size() && is<6 ; is++ ) { 
          markers[is] = cstyles[is];
        }
      }
 
 
      uselabels = true;
 
 
      std::cout << argv[0] << "\tuserlabels :" << usrlabels << ":" << std::endl;
 
      if ( rc.isTagDefined("RANGEMAP") )    RANGEMAP    = true;
      if ( rc.isTagDefined("ALLRANGEMAP") ) ALLRANGEMAP = true;
 
      std::cout << "Extra: " << rc.isTagDefined("Extra") << std::endl;
      
      if ( rc.isTagDefined("Extra") ) taglabels.push_back( fullreplace( rc.GetString("Extra"), "__", " " ) );
 
 
    }
    else { 
      std::cerr << argv[0] << ":\t config file not found: " << configfile << std::endl;
      return -1;
    }
 
  }
 
  if ( !use_file_config ) {
 
    std::cout << "using default panels" << std::endl;

 
    // for ( size_t iraw=0 ; iraw<histos_default.size() ; iraw++ ) { 
    //    Panel p( histos_default[iraw][0], 1 );
    //    p.push_back( histos_default[iraw] );
    //    panels.push_back( p );
    // }
 
    // default panel efficiencies plotted from 0 to 100, 
    // so scale efficiencies no matter what
 
    scale_eff     = 1;
    scale_eff_ref = 1;

 
    std::vector<panel_type> inpanels = { eff_panel, res_panel, diff_panel };
 
    std::string pnames[3] = { "eff", "res", "diff" };
 
    for ( size_t ip=0 ; ip<inpanels.size() ; ip++ ) { 
      Panel p( pnames[ip]+"_panel", 2 );
      for ( size_t ih=0 ; ih<inpanels[ip].size() ; ih++ ) p.push_back( inpanels[ip][ih] );
      panels.push_back( p );
    }
 
  }
 
 
  std::cout << "taglabels" << std::endl;
 
  for ( size_t it=0 ; it<taglabels.size()  ; it++ ) std::cout << taglabels[it]  << std::endl;
  for ( size_t it=0 ; it<taglabels2.size() ; it++ ) std::cout << taglabels2[it] << std::endl;
 
 
  std::cout << "\npanels: " << panels.size() << std::endl; 
 
  if ( panels.size()==0 ) return usage(argv[0], -1, "no panels to plot");
 
  for ( size_t ip=0 ; ip<panels.size() ; ip++ ) std::cout << panels[ip] << std::endl;
 

  if ( true ) { 
    gStyle->SetPadRightMargin(0.05);
    gStyle->SetPadTopMargin(0.05);

    const Int_t Number = 3;
    Double_t Red[Number]    = { 0.00, 0.00, 1.00 };
    Double_t Green[Number]  = { 0.00, 5.00, 1.00 };
    Double_t Blue[Number]   = { 0.00, 0.50, 0.00 };
    Double_t Length[Number] = { 0.00, 0.50, 1.00 };
    Int_t nb = 50;
    TColor::CreateGradientColorTable( Number, Length, Red, Green, Blue, nb );
  }
  else gStyle->SetPalette(1);
 
  if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
  double rightmargin = gStyle->GetPadRightMargin(); 
  gStyle->SetPadRightMargin(0.1); 

 
  gStyle->SetPadRightMargin(rightmargin); 
    
  // gStyle->SetLineScalePS(1);
 
  for ( size_t ipanel=0 ; ipanel<panels.size() ; ipanel++ ) { 
 
    Panel& panel = panels[ipanel]; 
 
    std::cout << blue << "\n\n---------------------------------------------\n" << reset;
 
    std::cout << panel << "\n" << std::endl;
 
    int ncolsp = panel.ncols();
    int nrowsp = panel.nrows();
    
    double extraw = 1;
 
    std::cout << "\nncols: " << ncolsp << "\tnrows: " << nrowsp << std::endl;
 
    bool multipanel = ( panel.size() > 1 );
 
    if ( panel.size()==0 ) { 
      std::cout << "panel empty: " << panel.name() << std::endl;
      continue;
    }
 
    gStyle->cd();
    
    if ( panel.size()>4 ) gStyle->SetLineScalePS(0.5);
 
    if ( multipanel ) extraw = 1.05;
 
    TCanvas* tc = new TCanvas( "tc", "", extraw*ncolsp*800, nrowsp*600 );
 
    tc->cd();
 
    const std::string& atlaslabel = atlaslabel_tmp;
 
    if ( multipanel ) { 
 
      gStyle->SetLineScalePS(1);



 
      tc->Divide( ncolsp, nrowsp, 0.0001, 0.0003 );
      //  atlaslabel = "     " + atlaslabel_tmp;
    }

 
    std::string plotname = ""; 
 
    for ( size_t i=0 ; i<panel.size() ; i++ ) {
 
      HistDetails histo = panel[i];
 
      bool drawmeans     = false;
      bool drawresiduals = true;
      
      //      if ( contains(histo.detail(), "+mean" ) || means )   drawmeans = true;
      if ( contains(histo.detail(), "+mean" )  )   drawmeans = true;
      if ( contains(histo.detail(), "-residual") ) drawresiduals = false;
 
      std::string xaxis = histo.xtitle();
      std::string yaxis = histo.ytitle();
      
      if ( !xregex.empty() ) { 
        size_t pos = xaxis.find(xregex);
    if ( pos!=std::string::npos ) xaxis.replace( pos, xregex.size(), xpattern );  
    pos = yaxis.find(xregex);
    if ( pos!=std::string::npos ) yaxis.replace( pos, xregex.size(), xpattern );  
      }
 
      const AxisInfo& xinfo = histo.xaxis(); 
      const AxisInfo& yinfo = histo.yaxis(); 
 
      // AxisInfo xinfo = histo.xaxis(); 
      // AxisInfo yinfo = histo.yaxis(); 
 
      std::string hname  = histo.name();
      std::string detail = histo.detail();
 
      bool translate_x = false;
 
      if ( ( lumiref_trans || lumitest_trans ) && contains(histo.xtitle(), "Pile-up" ) ) translate_x = true;
 
#if 0
      if ( translate_x ) { 
     xinfo.lo( xinfo.lo()*(1-0.054) );
     xinfo.lo(10);
     xinfo.hi( xinfo.hi()*(1-0.054) );
      }
#endif
      
      double rebin = 1;
 
#ifdef D0REBIN
      bool d0rebin_flag = false;
 
      if ( contains( detail, "+d0rebin" ) ) d0rebin_flag = true; 
#endif
      
      if ( contains( detail, "+Rebin" ) ) { 
    rebin = std::atof( detail.substr( detail.find("+Rebin")+6, detail.size() ).c_str() ); 
      } 
      
      if ( contains( detail, "+rebin" ) ) { 
    rebin = std::atof( detail.substr( detail.find("+rebin")+6, detail.size() ).c_str() ); 
      } 
      
      int  labelcolour = kBlack;
 
      int npanel = nrowsp*(i/nrowsp) + i%nrowsp + 1 ;
      
      std::cout << "panel: panel: " << panel.name() << "\tsubpanel: " << npanel << std::endl;
 
      if ( multipanel ) tc->cd( npanel );
 
      SetAtlasStyle();
 
      noreftmp = noref;
 
      Plotter::setplotref(!noreftmp);
      
      if ( fulldbg )  std::cout << __LINE__ << std::endl;
      
      std::cout << "main() processing histo[" << i <<  "] " << (i<10 ? " " : "" ) << histo.name() << "\t" << histo.xaxis() << std::endl;
      
      Plots plots_eff( "", yinfo.trim() );
      plots_eff.clear();
      
      Plots plots( "", yinfo.trim() );
      plots.clear();
      
      std::string noreflabel=defreflabel;
                  
      double xpos  = 0.18;
      double ypos  = 0.91;
      
      if ( contains(histo.name(),"eff") || contains(histo.name(),"Eff_") ) ypos = 0.19;
 
      //      ypos = 0.19;

 
      if ( atlasstyle ) { 
    xpos  = 0.18;
    if ( ypos>0.5 ) ypos = 0.85; 
    else            ypos = 0.18;
      }
 
      if ( ypos_in!=0 ) ypos = ypos_in;
      
      double xpos_original = xpos;
 
      xpos += xoffset;
      
      if ( xinfo.offset() != 0 ) { 
 
    std::cout << "HA ! xinfo.offset:  " << xinfo.offset() << std::endl; 
 
    xpos = xinfo.offset(); 
 
      }
      
 
      if ( yinfo.offset() != 0 ) { 
 
    std::cout << "HA ! yinfo.offset:  " << yinfo.offset() << std::endl; 
 
    ypos = yinfo.offset(); 
 
      }
      
 
      
 
      size_t Nrows = chains.size();
 
      if ( ALLRANGEMAP || (RANGEMAP && xaxis.find("p_{T}")!=std::string::npos && ccolours.size() ) ) {
    Nrows = ( Nrows < ccolours.size() ? Nrows : ccolours.size() );
      }
 
      int Nlines = Nrows + taglabels.size();
      
      std::vector<double> ypositions;
      std::vector<double> ypositions2;
      
      double deltay = (Nrows*0.055-0.005)/Nrows;
 
      double ylo = ypos;
      double yhi = ypos-0.01;
      
      if ( ypos>0.5 ) ylo -= Nlines*deltay;
      else            yhi += Nlines*deltay;
      
      ypositions.reserve(Nlines);
      for ( int ilines=0 ; ilines<Nlines ; ilines++ ) { 
        ypositions.push_back( yhi - deltay*(ilines+0.5) );
        ypositions2.push_back( yhi - deltay*(ilines-0.5) );
      }
 
      // specify different legends for efficiencies or residuals?
      
 
      if ( JLflag ) { 
    ylo -= 0.02;
      }
      else {
    if ( !contains(histo.name(),"eff") && !contains(histo.name(),"Eff_") ) ylo -= 0.02;
      }
 
      Legend legend(     xpos, xpos+0.1, ylo, ylo+Nrows*0.06-0.005 );
      Legend legend_eff( xpos, xpos+0.1, ylo, ylo+Nrows*0.06-0.005 );
            
      
      std::vector<std::string> Mean;
      std::vector<std::string> MeanLF;
      std::vector<std::string> RMS;
      
      Mean.clear();
      MeanLF.clear();
      RMS.clear();
      
      std::vector<std::string> Chi2;
      std::vector<std::string> MeanRef;
      std::vector<std::string> RMSRef;
      
      //    int colours[6] = { 1, 2, 4, 6, 7, 8 };
      
      Chi2.clear();
      MeanRef.clear();
      RMSRef.clear();
      
      int mean_power = 0;
      int rms_power  = 0;
      bool power_set = false;
      
      LINEF = true;
      
      //    bool uselogx = xinfo.log();
      //    bool uselogy = yinfo.log();
      
      for ( unsigned int j=0; j<chains.size(); j++)  {
    
      TFile* fftest = ftest;
      TFile* ffref  = fref;
 
      std::string cc = ftest->GetName();
      std::string cr = "";
 
      if ( fref ) cr = fref->GetName();
 
      if ( chainfiles.size()>0 && chainfiles.size()>j ) ffref = fftest = chainTFiles[j];

    
    std::cout << "chain:    " << chains[j] << "\taddchains: " << addchains << std::endl;
    std::cout << "chainref: " << chainref[j] << std::endl;
 
    noreftmp = noref;
    Plotter::setplotref(!noreftmp);
    
    TH1F* htest = 0;
    TH1F* href  = 0;
    
    TH1F* htestnum = 0;
    TH1F* htestden = 0;
    TH1F* hrefnum  = 0;
 
    TH2D* h2test = 0;
 
    
    TGraphAsymmErrors* tgtest = 0;
    
    std::cout << "refchain.size() " << refchain.size() << std::endl;
    
    std::cout << "refchain: " << refchain[j] << std::endl;
    
    
    gPad->SetRightMargin(0.03);
 
    bool do2D = false;
    
    labelcolour = kBlack;
 
 
    if ( contains(histo.name(),"/2d") ) { 
 
      do2D = true;
      labelcolour = kWhite;
      
      gPad->SetRightMargin(0.13);
 
      std::cout << "\n\nsee! it is a 2d histograms !!! " << histo.name() << "\n\n" << std::endl; 
      
      h2test = Get<TH2D>( *fftest, chains[j]+"/"+histo.name(), testrun, 0, &savedhistos );
 
      std::cout << "\t2d: " << h2test << std::endl;
      
      DrawLabel( 0.5, 0.5, "a test" );
      
      h2test->GetYaxis()->SetTitleOffset(1.55); 
      h2test->GetXaxis()->SetTitleOffset(1.5); 
      h2test->GetXaxis()->SetTitle(xaxis.c_str());
      h2test->GetYaxis()->SetTitle(yaxis.c_str());
          
      const AxisInfo& xinfo = histo.xaxis(); 
      const AxisInfo& yinfo = histo.yaxis(); 
      
      std::cout << xinfo << std::endl;
      std::cout << yinfo << std::endl;
 
      if ( yinfo.rangeset() ) { 
        h2test->GetYaxis()->SetRangeUser( yinfo.lo(), yinfo.hi() );
      }
 
      if ( yinfo.normset() ) normy( h2test );
 
      if ( xinfo.rangeset() ) { 
        h2test->GetXaxis()->SetRangeUser( xinfo.lo(), xinfo.hi() );
      }
 
      if ( xinfo.autoset() ) autox( h2test ); 
      if ( yinfo.autoset() ) autoy( h2test ); 
 
      SetZeros( h2test );
         
      h2test->DrawCopy("colz");
      
      if ( histo.detail().find("logz")!=std::string::npos ) gPad->SetLogz(true);
      else gPad->SetLogz(false);
     
          
    }
    else if ( refit_resplots && ( contains(histo.name(),"/sigma") || contains(histo.name(),"/mean") ) ) { 
 
      Resplot::setoldrms95(oldrms);
      Resplot::setscalerms95(true);
        
      std::cout << red << "fetch 2d: " << histo.name() << reset << std::endl;
      
      bool bsigma = false;
      if ( contains(histo.name(),"/sigma") ) bsigma = true;
 
      bool bmean = false;
      if ( contains(histo.name(),"/mean") ) bmean = true;
 
      std::string tmp_  = histo.name();
      std::string base;
        
      if ( bsigma ) base = chop( tmp_, "/sigma" );
      if ( bmean )  base = chop( tmp_, "/mean" );
 
      std::cout << "refitting:  " << histo.name() << "\tbase: " << base << "\tbsigma: " << bsigma << "\tbmean: " << bmean << std::endl;
        
      TH2D* htest2d = Get<TH2D>( *fftest, chains[j]+"/"+base+"/2d", testrun, 0, &savedhistos );
 
      if ( htest2d==0 ) htest2d = Get<TH2D>( *fftest, chains[j]+"/"+base+"/2D", testrun, 0, &savedhistos );
      
 
      TH2D* href2d = 0;
 
      if ( ffref )     href2d = Get<TH2D>( *ffref,  chains[j]+"/"+base+"/2d", testrun, chainmap  );
 
      if ( href2d==0 ) href2d = Get<TH2D>( *ffref,  chains[j]+"/"+base+"/2D", testrun, chainmap  );
      
      if ( htest2d==0 ) continue;
      if ( !noreftmp && href2d==0 )  noreftmp = true;
 
      std::cout << "href2d : " << href2d << std::endl; 
 
      Plotter::setplotref(!noreftmp);
      
      std::cout << "plotter" << std::endl;
      
        
      Resplot rtest("tmp", htest2d );
 
      std::cout << "Resplot rtest\t" << rtest.finalised() << std::endl;
 
      if ( rtest.finalised() ) { 
        std::cout << "refitting ..." << std::endl;
        if ( contains(histo.name(),"npix") || contains(histo.name(),"nsct") || contains(histo.name(),"nsi") || contains(histo.name(),"nbl") ) rtest.Refit(Resplot::FitNull);
        else  rtest.Refit(Resplot::FitNull95);
      }
      else {
        std::cout << "refitting (finalising) ..." << std::endl;
        if ( contains(histo.name(),"npix") || contains(histo.name(),"nsct") || contains(histo.name(),"nsi") || contains(histo.name(),"nbl") ) rtest.Finalise(Resplot::FitNull);
        else  rtest.Finalise(Resplot::FitNull95);
        
        rtest.Mean()->DrawCopy();
        gPad->Print("Duff.pdf");
      }
 
      std::cout << "bsigma: " << bsigma << std::endl;
      std::cout << "bmean:  " << bmean  << std::endl;
      
      if ( bsigma ) { htest = (TH1F*)rtest.Sigma()->Clone("rtest_sigma"); htest->SetDirectory(0); }
      if ( bmean  ) { htest = (TH1F*)rtest.Mean()->Clone("rtest_mean");   htest->SetDirectory(0); }
 
      std::cout << reset << std::endl;
      
      if ( htest==0 ) { 
        std::cerr << red << "missing test histogram: " << (refchain[j]+" / "+histo.name()) << " " << htest
              << "(test)" << reset <<  std::endl; 
        continue;
      }
      else {
        std::cout << green << "htest: " << htest->GetName() << reset << std::endl;
      }
 
      std::cout << "\nhisto.name(): " << histo.name() << std::endl;
 
      std::cout << reset << std::endl;
      
      if ( true && histo.name().find("d0_vs_phi")!=std::string::npos ) {
        
        TVirtualPad* old = gPad;
        
        std::cout << "old canvas: " << gPad << std::endl;
        
        std::cout << "dbg histo name: " << histo.name() << std::endl;       
 
        bool first = false;
 
        TCanvas* tc = 0;
 
        if ( histo.name().find("_rec")!=std::string::npos ) { 
          //if ( fitcanvas_rec == 0 ) { 
        fitcanvas_rec = new TCanvas( "fit_rec", "", 700, 600 );
        first = true;
        // }
          tc = fitcanvas_rec;
        }
 
        if ( histo.name().find("_rec")==std::string::npos ) { 
          //if ( fitcanvas == 0 ) { 
        fitcanvas = new TCanvas( "fit", "", 700, 600 );
        first = true;
        //}
          tc = fitcanvas;
        }
 
        TVirtualPad* tp = gPad; 
 
        if ( tc ) tc->cd();

        //      TH1D* hf = (TH1D*)rtest.Mean()->Clone("cck");
        //      TH1D* hf = rtest.Mean(); 
 
        rtest.Mean()->SetDirectory(0);
        
        TH1* hf = htest; 
 
        if ( histo.name().find("_rec")!=std::string::npos ) hf->SetTitle(";Trigger #phi;Trigger d_{0} [mm]");
        else                                                 hf->SetTitle(";Offline #phi;Offline d_{0} [mm]");
 
        TF1* tf = new TF1( "sinus", "[0]*sin(x+[1])+[2]" );
 
        tf->SetLineWidth(1);
        
        hf->GetYaxis()->SetTitleOffset( hf->GetYaxis()->GetTitleOffset()*1.1 );
 
        { 
          TCanvas* tfc = new TCanvas( "tmp", "", 700, 600 ); 
          tfc->cd();
          hf->Fit( tf );      
          delete tfc;
        }       
 
        double phi0 = tf->GetParameter(1);
        //      double x = 0.0001*int(tf->GetParameter(0)*cos(phi0)*10000);
        //      double y = 0.0001*int(tf->GetParameter(0)*sin(phi0)*10000);
        double x = tf->GetParameter(0)*cos(phi0);
        double y = tf->GetParameter(0)*sin(phi0);
        double off = tf->GetParameter(2);
 
        std::cout << "\tx = " << x << std::endl;
        std::cout << "\ty = " << y << std::endl;
        std::cout << "\toffset = " << off << std::endl;
 
        std::string cckstr = chains[j];
 
        cckstr = cckstr.substr( cckstr.find("HLT_IDTrack") );
 
        std::string vstr = chains[j];
 
        vstr.resize(vstr.find("_boffperf"));
        vstr.resize(vstr.find("_L1"));
 
        hf->SetMinimum( -0.03 );
        hf->SetMaximum(  0.03 );
 
        tc->cd();
 
        if ( first ) hf->DrawCopy();
        
        hf->DrawCopy("same");
        hf->DrawCopy("samee");
 
        
        DrawLabel( 0.2, 0.89, cckstr );
        DrawLabel( 0.2, 0.85, vstr );
 
        DrawLabel( 0.2, 0.81, "x = " + eround(x) + " mm" );
        DrawLabel( 0.2, 0.77, "y = " + eround(y) + " mm" );
        DrawLabel( 0.2, 0.73, "offset = " + eround(off) + " mm" );
 
        std::string sht = chains[j]+"-"+histo.name()+"-fit.pdf";
        
        size_t p = sht.find("/");
 
        while ( p!=std::string::npos ) { 
          sht.erase( p, 1 );
          p = sht.find("/");
        } 
          
 
        gPad->Print( sht.c_str() );
        gPad->Print( "sht.pdf" );
 
        std::cout << "CNT !!!" << std::endl;
 
        delete tf;
 
        tp->cd();
 
        gPad = old; 
        
      }
      

 
      TH1F* hreft  = 0;
 
      std::cout << "Resplot noreftmp " << noreftmp << " ... " << std::endl;
 
      if ( !noreftmp ) {
        if ( refitref_resplots ) { 
          
          std::cout << "Resplot rref:  " << href2d << std::endl;
 
          Resplot rref("tmp", href2d );
 
          std::cout << "Resplot rref:  " << href2d << " (2)" << std::endl;
 
          if ( rref.finalised() ) {
        std::cout << "refitting (2) ..." << std::endl;
        if ( contains(histo.name(),"npix") || contains(histo.name(),"nsct") ) rref.Refit(Resplot::FitNull);
        else  rref.Refit(Resplot::FitNull95);
          }
          else {
        std::cout << "refitting (2 - finalising) ..." << std::endl;
        if ( contains(histo.name(),"npix") || contains(histo.name(),"nsct") ) rref.Finalise(Resplot::FitNull);
        else  rref.Finalise(Resplot::FitNull95);
          }
 
          std::cout << "Resplot rref:  " << href2d << " (3)" << std::endl;
                  
          if ( bsigma ) { hreft = (TH1F*)rref.Sigma()->Clone("rref_sigma"); hreft->SetDirectory(0); }
          if ( bmean )  { hreft = (TH1F*)rref.Mean()->Clone("rref_mean");   hreft->SetDirectory(0); }
 
          std::cout << "Resplot rref:  " << href2d << " (4)" << std::endl;
        
        }
        else { 
 
          std::cout << "Resplot rref:  " << href2d << " (5)" << reset << std::endl;
          
          hreft = Get( *ffref, refchain[j]+"/"+histo.name(), rawrefrun, chainmap );
          if ( hreft==0 ) {
        std::cerr << "ERROR: could not find " << (refchain[j]+"/"+histo.name()) << std::endl;
          }
 
 
        }
      }
          
      std::cout << "Resplot rref:  " << href2d << " (6)" << std::endl;
 
 
      if ( !noreftmp && hreft==0 ) { 
 
        std::cerr << "missing ref histogram: " << (refchain[j]+" / "+histo.name()) << " " << htest << "(ref)" << std::endl; 
 
        noreftmp = true;
        Plotter::setplotref(!noreftmp);
        noreflabel="reference not found";
        //      std::exit(-1);
      }
          
      if ( !noreftmp ) { 
 
        std::cout << "Resplot rref:  " << href2d << " (7)" << std::endl;
 
        href = (TH1F*)hreft->Clone();
        href->SetDirectory(0);
 
        std::cout << "Resplot href:  " << href << " (7a)" << std::endl;
 
      }
                  
      std::cout << "Resplot rref:  " << href2d << " (8)" << std::endl;
 

      //        std::cout << "\tget " << (refchain[j]+"/"+histos[i]) << "\t" << href << std::endl;
        
      savedhistos.push_back( refchain[j]+"/"+histo.name() );
 
      std::cout << "Resplot rref:  " << href2d << " (9)" << std::endl;
       
    }
    else { 

        
      std::string reghist = histo.name();
 
      std::cout << "hist:  " << (chains[j]+"/"+reghist) << "\tftest " << ftest << std::endl;
      
 
      
      htest = Get( *fftest, chains[j]+"/"+reghist, testrun, 0, &savedhistos );
        
      std::cout << "\nhist: " << htest << std::endl;
 
      std::cout << "fftest: " << fftest->GetName() << std::endl;

      // if ( contains( fftest->GetName(), "-mc" ) ) translate_x = false;
 
      if ( lumitest_trans && translate_x && !contains( cc, "-mc" ) ) htest = trans( htest, translate_x );
      
      std::cout << xaxis << std::endl;
 
      if ( htest==0 ) { 
        std::cerr << "missing test histogram: " << (chains[j]+" / "+reghist) << " " << htest<< std::endl; 
        continue;
      }
 
      testfit( htest, htest->GetName() );
      
      TH1F* hreft = 0;
 
      if ( ffref ) hreft = Get( *ffref, refchain[j]+"/"+reghist, rawrefrun, chainmap );
      else noreftmp = true;
      
      if ( hreft && lumiref_trans && translate_x && !contains( cr, "-mc" ) ) hreft = trans( hreft, translate_x );
 
      std::cout << "hreft: " << hreft << std::endl;
 
      if ( std::string(htest->ClassName()).find("TH2")!=std::string::npos ) { 
        std::cout << "Class TH2: " << htest->GetName() << std::endl;
        continue;
      }
 
 
      if ( std::string(htest->ClassName()).find("TH1")!=std::string::npos ) { 
        std::cout << "Class TH1: " << htest->GetName() << std::endl; 
      }
      else if ( std::string(htest->ClassName()).find("TProfile")!=std::string::npos ) {  
        std::cout << "Class TProf: " << htest->GetName() << std::endl; 
      }
      else if ( std::string(htest->ClassName()).find("TEfficiency")!=std::string::npos ) {  
        std::cout << "Class TEff: " << htest->GetName() << std::endl; 
      }
 
 
      std::cout << "Resplot hreft:  " << hreft << " (10)" << " " << noreftmp << std::endl;
 
      if ( !noreftmp && hreft==0 ) { 
        std::cerr << "missing ref histogram: " << (refchain[j]+" / "+reghist)  
              << " " << hreft << std::endl;
        
        noreftmp = true;
        Plotter::setplotref(false);
        noreflabel="reference not found";
        //      std::exit(-1);
 
      }
 
      std::cout << "Resplot hreft:  " << hreft << " (10)" << std::endl;
              
      if ( hreft!=0 ) {
 
        std::cout << "Resplot hreft:  " << hreft->ClassName() << " (11)" << std::endl;
 
        if ( std::string(hreft->ClassName()).find("TH1")!=std::string::npos ) {
          href = (TH1F*)hreft->Clone();
          href->SetDirectory(0);
        }
        else if ( isTEfficiency(hreft) ) {
          href = makeplot( hreft );       
        }
      }
      else { 
        noreftmp = true;
        href = 0;
      }
 
      std::cout << " \tget " << (chains[j]+"/"+reghist)    << "\thtest " << htest << std::endl;
      std::cout << " \tget " << (refchain[j]+"/"+reghist)  << "\thref  " << href << std::endl;
 
      if ( htest==0 ) continue;
 
      if ( !noreftmp && href==0 ) { 
        noreftmp = true;
        Plotter::setplotref(!noreftmp);
      }
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
    
 
      if ( rebin!=1 ) {
        std::cout << "rebin: " << hname << "\t" << rebin << std::endl;
        if ( htest ) htest->Rebin(rebin);
        if ( href  ) href->Rebin(rebin);
        for ( int ip=0 ; ip<10 ; ip++ ) std::cout << std::endl; 
      }
 
#if 0
      if ( !contains( histo.name(), "rdz_vs_zed" ) && contains( histo.name(), "1d") ) {
        std::cout << "Rebinning histogram: " << histo.name() << std::endl;
        if (        htest->GetNbinsX()>500 ) htest->Rebin(10);
        if ( href && href->GetNbinsX()>500 ) href->Rebin(10);
      }
#endif
 
 
      if ( make_efficiencies && histo.name().find("zed_eff")!=std::string::npos ) { 
        if (        htest->GetNbinsX()>100 ) htest->Rebin(5);
        if ( href && href->GetNbinsX()>100 ) href->Rebin(5);
      }
 
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      if ( scalepix && std::string(htest->GetName()).find("npix")!=std::string::npos ) Scale(htest,0.5);
      if ( scalepix && href && std::string(htest->GetName()).find("npix")!=std::string::npos ) Scale(href,0.5);
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      if ( notitle ) { 
        htest->SetTitle("");
        if( href ) href->SetTitle("");
      } 
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    }
 
    std::cout << "done else" << std::endl;
 
    if ( do2D ) continue;
    
    std::cout << "done else:\t" << htest << std::endl;
    std::cout << "done else:\t" << htest->GetName() << std::endl;
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    if ( !do2D && make_ref_efficiencies ) { 
 
      std::cout << "make ref efficiencies" << std::endl; 
 
      if ( htest && href ) { 

        //    std::cout << "contains _eff " << contains( std::string(htest->GetName()), "eff" ) << std::endl;
 
        if ( contains( std::string(htest->GetName()), "eff" ) && !contains( std::string(htest->GetName()), "_d" ) ){
 
          std::string effhist = histo.name();
 
 
          
          htestnum = Get( *fftest, chains[j]+"/"+effhist+"_n", testrun, 0, &savedhistos );
 
 
          
          TH1F* hrefnumt = Get( *ffref, refchain[j]+"/"+effhist+"_n", rawrefrun, chainmap, &savedhistos );
 
 
          
          if ( !noreftmp && hrefnumt!=0 ) { 
        hrefnum = (TH1F*)hrefnumt->Clone();
        hrefnum->SetDirectory(0);
          }
 
        }
      }
 
    }      
 
    std::cout << "done effs ? " << bayes << "\thtest: " << htest << " :: " << htest->GetName() << std::endl;
 
    if ( bayes ) { 
 
      std::cout << "make ref efficiencies " << htest->GetName() << std::endl; 
 
      if ( make_efficiencies && htest && contains( std::string(htest->GetName()), "eff" ) && !contains( std::string(htest->GetName()), "_d" ) ) {
 
        std::string effhist = histo.name();
        
        
        if ( rebin!=1 ) std::cout << effhist << "\trebin: " << rebin << std::endl; 
 
        
        htestnum = Get( *fftest, chains[j]+"/"+effhist+"_n", testrun, 0, &savedhistos ) ;
        htestden = Get( *fftest, chains[j]+"/"+effhist+"_d", testrun, 0, &savedhistos ) ;
 
        if ( htestnum==0 && htestden==0 ) {
          TEfficiency* eff = 0;
          
          eff = Get<TEfficiency>( *fftest, chains[j]+"/"+effhist, testrun, 0, &savedhistos ) ;
 
          if ( eff ) {
        htestnum = (TH1F*)eff->GetPassedHistogram();
        htestden = (TH1F*)eff->GetTotalHistogram();
          }
          
        }
 
        std::cout << "eff: htestnum: " << htestnum << "\thtestden: " << htestden << std::endl; 
        
 
        double ars = total_eff( htestnum, htestden );
 
        std::cout << "heff: " << effhist << "\t" << ars << std::endl;
 
 
        std::cout << "1: Bayesian error calculation " << htestnum << " " << htestden << "\tscale " << scale_eff << std::endl;
 
        if ( htestnum && htestden ) { 
 
#ifdef D0REBIN
          if ( d0rebin_flag ) {
        htestnum = d0rebin( htestnum );
        htestden = d0rebin( htestden );
          }
          else if ( rebin!=1 ) { 
        htestnum = Rebin(htestnum, rebin );
        htestden = Rebin(htestden, rebin );
          }
#endif        
          std::cout << "test histogram name: : " << htestnum->GetName() << "\txaxis: " << xaxis << "\t" << std::endl;
          
          if ( make_efficiencies && std::string(htestnum->GetName()).find("ntrax_eff")!=std::string::npos ) {
 
        bool low = true;
 
        htestnum = rebin_log( htestnum, low );
        htestden = rebin_log( htestden, low );
          }
 
#if 0
          if ( contains( htest->GetName(), "_vs_lb" )  ) { 
        std::cout << "rebin " << histo.name() << std::endl;
        htestnum->Rebin(3);
        htestden->Rebin(3);
          }
 
          if ( make_efficiencies && contains( htest->GetName(), "eta_eff" )  ) { 
        std::cout << "rebin " << histo.name() << std::endl;
        htestnum->Rebin(2);
        htestden->Rebin(2);
          }
 
#endif  
 
          //    if ( RANGEMAP && (effhist.find("pT")!=std::string::npos || effhist.find("pt")!=std::string::npos ) ) { 
          if ( RANGEMAP && (effhist.find("ET")!=std::string::npos ) ) { 
        std::cout << "\trange: " << j << " " << htest << std::endl;
        bnd.range( chains[j], htestnum );
        bnd.range( chains[j], htestden );
          }
 
          if ( lumitest_trans && translate_x && !contains( cc, "-mc" ) ) { 
        htestnum = trans( htestnum, translate_x );
        htestden = trans( htestden, translate_x );
          }    
 
          Efficiency1D e( htestnum, htestden, "", scale_eff );
 
          tgtest = e.Bayes(scale_eff);
 
          htest = e.Hist();
          
          htest->SetName( (std::string(htestnum->GetName())+"_eff").c_str() );
 
          std::cout << "effhist: " << effhist << std::endl;
 
        }
    
 
        std::cout << "recalculating reference efficiencies ..." << std::endl; 
 
        if ( href ) { 
 
          std::cout << "doin ..." << std::endl; 
 
 
          TH1F* hrefnum = Get( *ffref, refchain[j]+"/"+histo.name()+"_n", rawrefrun, chainmap );
          TH1F* hrefden = Get( *ffref, refchain[j]+"/"+histo.name()+"_d", rawrefrun, chainmap );
 
 
          if ( hrefnum==0 && hrefden==0 ) {
        TEfficiency* eff = 0;
        
        eff = Get<TEfficiency>( *ffref, refchain[j]+"/"+histo.name(), rawrefrun, chainmap ) ;
        
        if ( eff ) {
          hrefnum = (TH1F*)eff->GetPassedHistogram();
          hrefden = (TH1F*)eff->GetTotalHistogram();
        }
        
          }
 
          
          
          std::cout << "2. Bayesian error calculation " << htestnum << " " << htestden << "\tscale " << scale_eff     << std::endl;
          std::cout << "3. Bayesian error calculation " << hrefnum  << " " << hrefden  << "\tscale " << scale_eff_ref << std::endl;
 
        
          if ( hrefnum && hrefden ) { 
 
        if ( lumiref_trans && translate_x && !contains( cr, "-mc" ) ) { 
          hrefnum = trans( hrefnum, translate_x );
          hrefden = trans( hrefden, translate_x );
        }
 
#ifdef D0REBIN
        if ( d0rebin_flag ) {
          hrefnum = d0rebin( hrefnum );
          hrefden = d0rebin( hrefden );
        }
        else if ( rebin!=1 ) { 
          hrefnum = Rebin(hrefnum, rebin );
          hrefden = Rebin(hrefden, rebin );
        }
#endif
        
        std::cout << hrefnum << " " << hrefden << std::endl;
 
        Efficiency1D e( hrefnum, hrefden, "", scale_eff_ref );
        // tgref = e.Bayes(scale_eff);
 
        href = e.Hist();
          
          }
        }
      }
 
      std::cout << "dbg: htest:     " << htest << std::endl; 
      std::cout << "dbg: htest name " << htest->GetName() << std::endl; 
 
      if ( true && histo.name().find("invmass")!=std::string::npos ) {
        testfit( htest, histo.name() );
      }
 
      
    }
 
      
    std::cout << "chains[j]    " << chains[j] << std::endl;
    std::cout << "histo.name() " << histo.name() << std::endl;
    std::cout << "href         " << href << "\t (TEfficiency " << isTEfficiency(href) << ")" << std::endl;
 
    if ( isTEfficiency(href) ) {
      //      TH1F* hrefeff = makeplot<TEfficiency>(href);
      TH1F* hrefeff = makeplot(href);
      href = hrefeff;
      std::cout << "efficiency" << std::endl;
    }
      
    std::cout << "href         " << href << "\t (TEfficiency " << isTEfficiency(href) << ")" << std::endl;
 
    
    if ( htest==0 ) { 
      std::cout << "       no test histogram :   " << (chains[j]+"/"+histo.name()) << std::endl;
      continue;
    }
 
 
    if ( !noreftmp && href==0 ) { 
      std::cout << "       no ref histogram :    " << (chains[j]+"/"+histo.name()) << std::endl;
      noreftmp = true;
      Plotter::setplotref(!noreftmp);
      noreflabel="reference not found";
      //      continue;
    }
 
    std::cout << "htest: " << isTEfficiency(htest) << std::endl;
 
    if ( isTEfficiency(htest) ) htest = makeplot(htest);
    
    htest->GetYaxis()->SetTitleOffset(1.55); 
    htest->GetXaxis()->SetTitleOffset(1.5); 
    htest->GetXaxis()->SetTitle(xaxis.c_str());
    htest->GetYaxis()->SetTitle(yaxis.c_str());
 
    std::cout << "htest: " << htest << "\t" << __LINE__ << std::endl;
 
    
    if ( !noreftmp ) { 
 
      href->GetYaxis()->SetTitleOffset(1.5);
      href->GetXaxis()->SetTitleOffset(1.5);
      std::cout << "xaxis: " << href->GetXaxis()->GetTitle() << std::endl;
 
    }   
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
#if 0
    if ( !contains(histo.name(),"vtx_eff") && contains(histo.name(),"ntracks") ) {
 
      double xm  = htest->GetMean();
 
      if ( xm>=10 ) { 
        double lxm = std::log10(xm);
        int newbins = int(0.5+xm/std::pow(10,int(lxm)))*pow(10,int(lxm));
        int nrebin   = int( (newbins+5)/10 );
        
        if ( nrebin>1 ) { 
          std::cout << "rebin: " << htest->GetName() << "\tbins: " << nrebin << std::endl;
          htest->Rebin(nrebin);
          htest->Sumw2();
          if ( !noreftmp ) { 
        href->Rebin(nrebin);
        href->Sumw2();
          }
        }
      }
    }
#endif
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 

    if ( plotname == "" ) { 
 
      if ( key!="" ) { 
        htest->SetTitle("");
        if ( href ) href->SetTitle("");
        plotname = key+"_";
      }
      else if (fcontains(chains[j],"HLT_")) {
        htest->SetTitle("");
        if ( href ) href->SetTitle("");
        plotname = "HLT_";
      }
      else if (fcontains(chains[j],"EF_")) {
        htest->SetTitle("");
        if ( href ) href->SetTitle("");
        plotname = "EF_";
      }
      else if (fcontains(chains[j],"L2_")) {
        htest->SetTitle("");
        if ( href ) href->SetTitle("");
        plotname = "L2_";
      }
      else if (contains(chains[j],"FTK") && ! contains(chains[j],"HLT_") ) { 
        htest->SetTitle(("FTK "+ histo.name()).c_str());
        if ( href ) href->SetTitle(("FTK "+ histo.name()).c_str());
        plotname = "FTK_";
      }
    
      plotname += histo.name(); 

      replace(plotname, '/', '_'); 
 
    }
 
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    bool residual = false;
      
    if ( contains(histo.name(),"_res") ||  contains(histo.name(),"residual_") || contains(histo.name(),"1d") ) residual = true; 
 
    std::string collection = basename( chains[j] );
 
    std::string actual_chain = basename( dirname( chains[j] ) );
 
    if ( collection.find("_InDet")!=std::string::npos ) collection.erase( 0, collection.find("_InDet")+1 );
    if ( actual_chain.find("_InDet")!=std::string::npos ) actual_chain.erase( actual_chain.find("_InDet") );
    
 
    std::cout << "raw:              " << chains[j]    << std::endl;
    std::cout << "track collection: " << collection   << std::endl;
    std::cout << "actual chain:     " << actual_chain << std::endl;
 
    std::regex rx("_HLT_[^_]*RoI.*");
    std::regex rx1("_HLT_[^_]*_RoI.*");
 
    std::cout << "\n\n\nactual_chain: " << actual_chain << std::endl;
 
    arsaway( actual_chain, "_L1.*", "" );
    arsaway( actual_chain, "L1.*", "" );
 
        actual_chain = std::regex_replace( actual_chain, std::regex( "L1.*:" ), "" );
 
    if ( actual_chain.find("jvt")!=std::string::npos )  actual_chain.erase( 
                                    actual_chain.find("jvt"), 
                                    actual_chain.find("HLT_IDTrack")-actual_chain.find("jvt") );
 
    std::cout << "            : " << actual_chain << "\n\n\n" << std::endl;
 
    
        actual_chain = std::regex_replace( actual_chain, std::regex( "_HLT_IDTrack.*" ), "" );
 
    
        collection = std::regex_replace(  collection, std::regex(".*HLT_IDTrack_"), "IDTrack " );
        collection = std::regex_replace(  collection, std::regex("IDTrack "), "" );
    collection = std::regex_replace(  std::regex_replace( collection, rx, "" ), rx1, "" );
 
    if ( actual_chain.find("HLT_IDTrack_")!=std::string::npos ) actual_chain.erase( actual_chain.find("HLT_IDTrack_"), 12 );
    if ( actual_chain.find("_IDTrack_")!=std::string::npos    ) actual_chain.erase( actual_chain.find("_IDTrack_"), 9 );
    if ( actual_chain.find("IDTrack")!=std::string::npos )   actual_chain.erase( actual_chain.find("IDTrack"), 7 );
    if ( actual_chain.find("_idperf")!=std::string::npos )   actual_chain.erase( actual_chain.find("_idperf"), 7 );
    if ( actual_chain.find("_bperf")!=std::string::npos  )   actual_chain.erase( actual_chain.find("_bperf"), 6 );
    if ( actual_chain.find("_boffperf")!=std::string::npos ) actual_chain.erase( actual_chain.find("_boffperf"), 9 );
    if ( actual_chain.find("_HLT_")!=std::string::npos )     actual_chain.replace( actual_chain.find("_HLT_"), 5, " " );
    if ( actual_chain.find("HLT_")!=std::string::npos  )     actual_chain.erase( actual_chain.find("HLT_"), 4 );
 
 
 
    if ( collection.find("_IDTrkNoCut")!=std::string::npos  )    collection.erase( collection.find("_IDTrkNoCut"), 11 );
    if ( collection.find("xAODCnv")!=std::string::npos      )    collection.erase( collection.find("xAODCnv"), 7 );
    if ( collection.find("HLT_IDTrack_")!=std::string::npos )    collection.erase( collection.find("HLT_IDTrack_"), 12 );
    if ( collection.find("HLT_IDTrack")!=std::string::npos  )    collection.erase( collection.find("HLT_IDTrack"), 11 );
    if ( collection.find("Tracking")!=std::string::npos     )    collection.replace( collection.find("Tracking"), 8, "Trk" );    
    if ( collection.find("InDetTrigTrk_")!=std::string::npos   ) collection.erase( collection.find("InDetTrigTrk_"), 13 );    
    if ( collection.find("HLT_xAODTracks_")!=std::string::npos ) collection.erase( collection.find("HLT_xAODTracks_"), 15 );    
    if ( collection.find("_HLT_")!=std::string::npos  )  collection.replace( collection.find("_HLT_"), 5, " " );    
    if ( collection.find("HLT_")!=std::string::npos   )  collection.erase( collection.find("HLT_"), 4 );    
    if ( collection.find("_extra")!=std::string::npos )  collection.resize( collection.find("_extra") );
 
    if ( actual_chain.size()>25 ) {
      size_t pos = actual_chain.find_last_of("_");
      while ( pos!=std::string::npos && actual_chain.size()>25 ) {
        actual_chain.resize(pos);
        pos = actual_chain.find_last_of("_");
      }
    }
 
    std::cout << blue << "actual chain: " << actual_chain << "\tcollection:   " << collection   << reset << std::endl; 
    
    std::string c = actual_chain + " : " + collection;
 
    std::cout << "track collection: " << collection   << "   <-" << std::endl;
    std::cout << "actual chain:     " << actual_chain << "   <-" << std::endl;
       
    replace( c, "_In", " :  " );
    
    c = "  " + c;
 
    std::cout << "use label: " << c << "\tchains size " << chains.size() << "\t" << usrlabels.size() << std::endl;
 
    std::cout << "chains.size() " << chains.size() << "  " << j <<  std::endl; 
    std::cout << "chains.size() " << chains.size() << "  " << j <<  std::endl; 
 
    std::cout << "chains[j] : " << j << " " << chains[j] << std::endl; 
 
    std::cout << "userlabels.size() " << usrlabels.size() << std::endl; 
 
    if ( usrlabels.size() < j+1 ) { 
      if ( usrlabels.size()!=0 ) std::cerr << "userlabels not large enough - not using userlabels" << std::endl;
    }
    else c = usrlabels[ j ]; 
 
    std::cout << "use label: c: " << c << std::endl;

 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
            
    std::cout << blue << "xaxis: " << xaxis << reset << std::endl;
 
    if ( ALLRANGEMAP || xaxis.find("p_{T}")!=std::string::npos || xaxis.find("E_{T}")!=std::string::npos ) {
 
      if ( RANGEMAP && xaxis.find("p_{T}")!=std::string::npos ) {
        bnd.range( chains[j], htest );
        if ( href ) bnd.range( chains[j], href );
      }
 
          if ( RANGEMAP || ALLRANGEMAP ) { 
 
        LINEF = false;
      
        std::cout << "\n\n\nctags " << ctags.size() << "\n\n" << std::endl;
 
        for ( size_t ic=0 ; ic<ctags.size() ; ic++ ) { 
          
          std::cout << "\n\nic: " << ic << " " << ctags[ic] << " " << ccolours[ic] << "\n\n" << std::endl;
 
          if ( chains[j].find(ctags[ic])!=std::string::npos ) {
        std::cout << "\ttag " << ctags[ic] << "  \tcolour: " << ccolours[ic] << "\tstyle: " << cstyles[ic] << std::endl;
        htest->SetLineColor( ccolours[ic] );
        htest->SetMarkerColor( ccolours[ic] );
        htest->SetMarkerStyle( cstyles[ic] );
        c = ctaglabels[ic];
        if ( tgtest ) { 
          tgtest->SetLineColor(htest->GetMarkerColor());
          tgtest->SetMarkerStyle(htest->GetMarkerStyle());
          tgtest->SetMarkerColor(htest->GetMarkerColor());
        }
        break;
          }
        }
 
      }
 
      std::cout << "test: " << chains[j] << "chains colour: " << htest->GetMarkerColor() << std::endl;
 
    } 
 
    
    std::cout << "movin' on ..." << std::endl;
 
    std::cout << "chain: " << chains[j] << " \t marker colour: " << htest->GetMarkerColor() << std::endl;
 
    std::cout << "Plotter marker : " << htest->GetMarkerColor() << " " << htest->GetMarkerStyle() << std::endl;
 
        // do we really want ...  chain_name+usrlabels[j] ?
    if ( uselabels )  plots.push_back( Plotter( htest, href, " " + chain_name[j] + c, tgtest ) );
    else {
      std::cout << "using label: " << c << std::endl;
      plots.push_back( Plotter( htest, href, c, tgtest ) );
    }
 
    std::cout << "AxisInfo: " << xaxis << std::endl;
 
    plots.back().xaxis( &(histo.xaxis()) );
    plots.back().yaxis( &(histo.yaxis()) );
    
    std::cout << "c: " << c << "\t" << ftest->GetName() << std::endl;
 
    if ( refbands ) { 
      std::cout << "refbands: " << refbands << std::endl;
      plots.back().mc( true ); 
      translate_x = false;
    }
 
    if ( ALLRANGEMAP || ( RANGEMAP && xaxis.find("p_{T}")!=std::string::npos )  ) plots.back().max_entries( ccolours.size() );
    
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    if ( make_ref_efficiencies ) { 
    
      if ( htestnum && hrefnum ) { 
        Efficiency1D e( htestnum, hrefnum, "", scale_eff );
 
        TH1* h = e.Hist();
 
        double range = h->GetMaximum()-h->GetMinimum();
 
        if ( range<0.2*scale_eff ) {
 
          double fmax = int( (h->GetMaximum() + 20)*0.1 )*0.1*scale_eff;
          double fmin = int( (h->GetMinimum() - 10)*0.1 )*0.1*scale_eff;
        
          if ( fmax>1*scale_eff ) fmax = 1.02*scale_eff;
          if ( fmin<0 )           fmin = 0;
        
          h->SetMinimum(fmin);
          h->SetMaximum(fmax);
       
        }
 
        plots_eff.push_back( Plotter( e.Hist(), 0, c ) );
      
      }     
    }
      
 
    std::cout << "means: " << drawmeans << " " <<  means << std::endl;
 
    std::cout << "chi2: " << std::endl;
 
    std::cout << "htest " << htest << std::endl; 
    std::cout << "href  " << href  << std::endl; 
    
    std::cout << "chi2: " << chi2( htest, href ) << "\t histo name: " << histo.name() << std::endl;
    
    testfit( htest, histo.name() );
 
 
    if ( href ) Chi2.push_back( label( "chi2 = %5.2lf / %2.0lf", chi2( htest, href ), double(htest->GetNbinsX()) ) );
 
    std::cout << "Chi2.push_back()" << std::endl; 
 
    if ( drawmeans || means ) {
      
          double   mean_95 = htest->GetMean();
          double  dmean_95 = htest->GetMeanError();
          double    rms_95 = htest->GetRMS();
          double   drms_95 = htest->GetRMSError();
          
          Mean.push_back(label("     mean = %4.2lf #pm %4.2lf", mean_95, dmean_95) );
          MeanLF.push_back(label("mean = %6.4lf #pm %6.4lf", mean_95, dmean_95) );
          RMS.push_back(label( "     rms   = %4.2lf #pm %4.2lf", rms_95,  drms_95 ) );
          
        }
    
    std::cout << "residual: " << residual << " " << drawresiduals << std::endl;
 
    if ( residual && drawresiduals ) {
    
 
      xpos = xpos_original;
 
      std::cout << "calculating resolutions : " << histo.name() << " " << htest->GetName() << std::endl;
 
      TF1* d95 = Resplot::FitNull95(reinterpret_cast<TH1D*>(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 <<  "\t\t" << histo.name()
            << "\tmean:     " << mean_95 << " +- " << dmean_95 
            << "\trms:     "  <<  rms_95 << " +- " << drms_95 << std::endl; 

 
      if ( !power_set ) { 
        for ( int ip=-2 ; ip<9 ; ip++ ) { 
          if ( std::fabs(mean_95) >= std::pow( 10., double(-ip) ) ) { 
        mean_power = ip;
        break;
          }
        }
      
        for ( int ip=-2 ; ip<9 ; ip++ ) { 
          if ( std::fabs(rms_95)  >= std::pow( 10., double(-ip) ) ) { 
        rms_power = ip;
        break;
          }
        }
      }
 
      power_set = true;
 
      std::cout <<  "\t\t" << histo.name() 
            << "\tmean:     " << mean_95 << " +- " << dmean_95 << " : pow " << mean_power
            << "\trms:     "  <<  rms_95 << " +- " << drms_95  << " : pow " << rms_power << std::endl; 
 
 
      if ( mean_power == 0 ) { 
        Mean.push_back(label("mean_{95}  = %4.2lf #pm %4.2lf", mean_95, dmean_95) );
      }
      else { 
        Mean.push_back(label("mean_{95}  = ( %4.2lf #pm %4.2lf ) #times 10^{%d}", 
                 mean_95*std::pow(10.,double(mean_power)), dmean_95*std::pow(10,double(mean_power)), -mean_power ) );
      
      }
 
      if ( rms_power == 0 ) { 
        RMS.push_back(label( "rms_{95}   = %4.2lf #pm %4.2lf", rms_95,  drms_95 ) ); 
      }
      else { 
        RMS.push_back(label( "rms_{95}   = ( %4.2lf #pm %4.2lf ) #times 10^{%d}", 
                 rms_95*std::pow(10.,double(rms_power)),  drms_95*std::pow(10,double(rms_power)), -rms_power ) );
      }
 
      if ( href ) { 
        TF1* d95ref = Resplot::FitNull95(reinterpret_cast<TH1D*>(href));
 
        double   mean_95ref = d95ref->GetParameter(1);
        double  dmean_95ref = d95ref->GetParError(1);
        double    rms_95ref = d95ref->GetParameter(2);
        double   drms_95ref = d95ref->GetParError(2);
      
        std::cout <<  "\t\t" << histo.name()
              << "\tmean ref: " << mean_95ref << " +- " << dmean_95ref << " : pow " << mean_power
              << "\trms ref: "  <<  rms_95ref << " +- " << drms_95ref  << " : pow " << rms_power << std::endl;
      
        if ( mean_power == 0 ) { 
          MeanRef.push_back(label("mean_{95} ref = %4.2lf #pm %4.2lf", mean_95ref, dmean_95ref) );
        }
        else { 
          MeanRef.push_back(label("mean_{95} ref = ( %4.2lf #pm %4.2lf ) #times 10^{%d}", 
                      mean_95ref*std::pow(10,double(mean_power)), dmean_95ref*std::pow(10,double(mean_power)), -mean_power ) );
        
        }
      
      
        if ( rms_power == 0 ) { 
          RMSRef.push_back(label( "rms_{95}  ref = %4.2lf #pm %4.2lf", rms_95ref,  drms_95ref ) ); 
        }
        else { 
          RMSRef.push_back(label( "rms_{95}  ref = ( %4.2lf #pm %4.2lf ) #times 10^{%d}", 
                      rms_95ref*std::pow(10,double(rms_power)),  drms_95ref*std::pow(10,double(rms_power)), -rms_power ) );
        }
      }
    
      htest->Sumw2();
      if ( href ) href->Sumw2();
    }
 
    if ( yinfo.normset() ) { 
      Norm( htest );
      if ( href ) Norm( href  );
      if ( xinfo.rangeset() ) { 
        Norm( htest, 1, xinfo.lo(), xinfo.hi() );
        if ( href ) Norm( href, 1, xinfo.lo(), xinfo.hi() );
      } 
    }
    else if ( yinfo.refnormset() ) { 
      if ( href ) { 
        Norm( href, Entries(htest) );
        if ( xinfo.rangeset() ) { 
          Norm( href, 1, xinfo.lo(), xinfo.hi() );
        } 
      }
    }
 
    if ( yinfo.binwidth() ) { 
      binwidth( htest );
      if ( href ) binwidth( href );
    }
 
    if ( !noreftmp && normref && 
         !contains( histo.name(), "mean") && !contains( histo.name(), "sigma" ) && 
         !contains( histo.name(), "Eff")  && !contains( histo.name(), "eff") &&
         !contains( histo.name(), "Res")  && !contains( histo.name(), "vs") && !contains( histo.name(), "_lb") ) { 
      Norm( href, Entries( htest ) );
    }
 
      }
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      if ( !noplots ) { 

 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
    
    std::cout << "duff: " << xinfo << std::endl;

    plots.sortx( xinfo );
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    double  yminset = 0;
    double  ymaxset = 0;
 
    double rmin = 0;
    double rmax = 0;
      
    if ( xinfo.rangeset() ) { 
      rmin = plots.realmin( plots.lo(), plots.hi() );
      rmax = plots.realmax( plots.lo(), plots.hi() );
      // if ( yinfo.normset() ) { 
      //  Norm( htest, 1, rmin, rmax ); 
      //  if ( href ) Norm( href, 1, rmin, rmax );
      // } 
    }
    else {
      rmin = plots.realmin();
      rmax = plots.realmax();
    }
        
    if ( yinfo.autoset() ) { 
    
      int csize = chains.size() + taglabels.size() + ( atlasstyle ? 1 : 0 );
 
      if ( yinfo.log() && rmin>0 && rmax>0 ) { 

        double delta = std::log10(rmax)-std::log10(rmin);

        //    ymaxset =  rmax*std::pow(10,delta*0.15*csize); 
 
        yminset =  rmin*std::pow(10,-delta*0.1);
 
        double newdelta = std::log10(rmax) - std::log10(yminset) + 0.05*delta;
 
        if ( csize<10 ) ymaxset =  rmin*std::pow(10,newdelta/(1-0.07*csize));
        else            ymaxset =  rmin*std::pow(10,newdelta*2);
 
        if ( yminset!=yminset ) { 
          std::cerr << " range error " << delta << " " << yminset << " " << ymaxset << "\t(" << rmin << " " << rmax << ")" << std::endl;
          continue;
        }
 
      }
      else { 

 
        if ( ypos>0.5 ) { 
          double delta = rmax-rmin;
        
          yminset = rmin-0.1*delta;
        
          if ( rmin>=0 && yminset<=0 ) yminset = 0;
        
          double newdelta = rmax - yminset + 0.05*delta;
        
          if ( csize<10 ) ymaxset = yminset + newdelta/(1-0.09*csize);
          else            ymaxset = yminset + newdelta*2;
        }
        else { 
          double delta = rmax-rmin;
  
          ymaxset = rmax+0.1*delta;
                
          double newdelta = ymaxset - rmin - 0.05*delta;
        
          if ( csize<10 ) yminset = ymaxset - newdelta/(1-0.09*csize);
          else            yminset = ymaxset - newdelta*2;
 
          if ( rmin>=0 && yminset<=0 ) yminset = 0;
 
        }
 
      }
    
    }
    else {  
      if ( yinfo.rangeset() ) { 
        yminset = yinfo.lo();
        ymaxset = yinfo.hi();
      }
    }
      
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    // more useful debugging ...
    //      std::cout <<  "yauto: " << yinfo.autoset() << "\tyrange " << yinfo.rangeset() << std::endl;
 
    //      std::cout << "yminset " << yminset << "\tymaxset " << ymaxset << std::endl;  
 
    if ( yinfo.autoset() && yinfo.rangeset() ) {
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      if ( yminset>yinfo.lo() ) yminset = yinfo.lo();
      if ( ymaxset<yinfo.hi() ) ymaxset = yinfo.hi();
    }
      
 
    if ( make_efficiencies && contains(histo.name(),"_eff") ) {
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
      
      if ( effset ) { 
        ymaxset = effmax;
        yminset = effmin;
      }
    }
      
    if ( ymaxset!=0 || yminset!=0 ) {
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      plots.Max( ymaxset );
      plots.Min( yminset );
    }
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    if ( yminset!=0 || ymaxset!=0 ) { 
      if ( yminset>0 ) plots.SetLogy(yinfo.log());
      else             plots.SetLogy(false);
    }   
    else plots.SetLogy(yinfo.log());
 
    //      plots.SetLogy(false);
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
            
    //    plots.limits();
      
      
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
    
    plots.Draw( legend );
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
    if ( atlasstyle ) ATLASLabel( xpos, ypositions[0]+deltay, atlaslabel, labelcolour, ncolsp, nrowsp );
 
    if ( fulldbg ) std::cout << __LINE__ << std::endl;
          
    for ( unsigned it=0 ; it<taglabels.size() ; it++ ) { 
      DrawLabel( xpos, ypositions[it], taglabels[it], labelcolour, 0.04 );
    }
 
    for ( unsigned it=0 ; it<taglabels2.size() ; it++ ) { 
      DrawLabel( xpos+0.35, ypositions2[it], taglabels2[it], labelcolour, 0.04 );
    }
 
      }
    
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      std::cout << "means: " << means << std::endl;
      
      //      if ( means ) for ( unsigned  j=0 ; j<MeanLF.size() ; j++ ) DrawLabel( 0.5, (0.65-(j+chains.size()-1)*0.035), MeanLF[j],  colours[j%6] ); // DrawLabel( 0.75, 0.85-j*0.035, Chi2[j], colours[j%6] );
      if ( means ) { 
    if ( runnumber!="" ) DrawLabel( 0.2, 0.89, "run "+runnumber );
    for ( unsigned  j=0 ; j<MeanLF.size() ; j++ ) DrawLabel( 0.2, (0.85-(j+chains.size()-1)*0.035), MeanLF[j],  colours[j%6] );
      }
                                                                  
 
      //      double ytop = 0.67;
      double ytop = 0.5;
      
      if ( ( !nostats || !nomeans ) && !noplots ) { 
    if ( dochi2 ) for ( unsigned  j=0 ; j<Chi2.size() ; j++ ) DrawLabel( 0.75, 0.85-j*0.035, Chi2[j], colours[j%6] );
    if ( ( (contains(histo.name(),"_res") || 
        contains(histo.name(),"1d")   ||
        histo.name()=="pT"            || 
        contains(histo.name(),"residual_") ||
        contains(histo.name(),"vs_pt") ) && !contains(histo.name(),"sigma") ) || drawmeans ) { 
    
      if ( contains(histo.name(),"_res") || contains(histo.name(),"residual_") || contains(histo.name(),"1d") || drawresiduals ) {
        for ( unsigned j=0 ; j<chains.size() ; j++ ) { 
          if ( !noreftmp ) { 
        if ( j<MeanRef.size() ) {
          if ( !nomeans ) DrawLabel( xpos_original, (ytop-j*0.035), MeanRef[j], colours[j%6] );
          DrawLabel( xpos_original, (ytop-0.035*chains.size()-j*0.035)-0.01, RMSRef[j],  colours[j%6] );
        }
          }
          if ( j<Mean.size() ) {
        if ( !nomeans ) DrawLabel( 0.62, (ytop-j*0.035), Mean[j],  colours[j%6] );
        DrawLabel( 0.62, (ytop-0.035*chains.size()-j*0.035)-0.01, RMS[j],  colours[j%6] );
          }
        }
      }
    }
      }
 
      if ( xinfo.log() ) gPad->SetLogx(true);
      else               gPad->SetLogx(false);
 
      if ( yinfo.log() ) gPad->SetLogy(true);
      else               gPad->SetLogy(false);
 
      if ( fulldbg ) std::cout << __LINE__ << std::endl;
 
      if ( !noplots ) { 
      
    if ( make_ref_efficiencies ) { 
    
      plots_eff.SetXaxisTitle( plots.GetXaxisTitle() ); 
      plots_eff.SetYaxisTitle( plots.GetYaxisTitle() ); 
    
      plots_eff.Draw( legend_eff );
    }
 
    if ( !noreflabel.empty() ) DrawLabel(0.1, 0.06, noreflabel, kRed, 0.03 );
 
      } // no plots
    
    } 
 

 
    if ( !noplots ) { 
 
      if ( !quiet ) { 
 
    tc->cd();
 
    std::string useplotname;
 
        if ( panel.size()>1 ) {
      useplotname = panel.name();
      replace( useplotname, '/', '_' );
    }   
    else { 
      useplotname = plotname;
    }
 
    useplotname.erase( std::remove( useplotname.begin(), useplotname.end(), '+' ), useplotname.end() );
 
    std::string printbase = dir + useplotname + tag;
 
    tc->Update();
 
    if ( !nopdf ) print_pad( printbase+".pdf" );
    if ( !nopng ) print_pad( printbase+".png" );
    if ( Cfile  ) print_pad( printbase+".C" );
    
    std::cout << std::endl;
      }
      
    }
 
    if ( tc ) delete tc;
 
  } 
 
  
  if ( fulldbg ) std::cout << __LINE__ << std::endl;
  
  
  
  bool files_duplicated = ( fref_==ftest_ );
  
  if ( deleteref && !files_duplicated ) {
    
    if ( fref_ ) { 
      
      
      std::cout << "main() cleaning up reference file" << std::endl;
      
      TFile* newout = new TFile(".newout.root","recreate"); 
      newout->cd();
      
      copyReleaseInfo( fref, newout );
      
      TDirectory* base = gDirectory;
      
      for ( unsigned i=0 ; i<savedhistos.size() ; i++ ) { 
    
    //  std::cout << i << " " << savedhistos[i] << std::endl;
    
    std::vector<std::string> dirs = AxisInfo::split( savedhistos[i], "/" );
    
    for ( unsigned j=0 ; j<dirs.size()-1 ; j++ ) { 
      std::cout << "\t" << dirs[j] << std::endl;
      TDirectory* renedir = gDirectory->GetDirectory( dirs[j].c_str() );
      if ( renedir==0 ) gDirectory->mkdir( dirs[j].c_str() );
      gDirectory->cd( dirs[j].c_str() );
    }
    
 
    TH1* href  = Get( *fref, savedhistos[i].c_str(), "", chainmap );
    if ( !noreftmp && href ) {
      std::cout << i << " " << savedhistos[i] << " 0x" << href << std::endl;
      href->Write( dirs.back().c_str() );
    }
      
            
    base->cd();
      }
    
      newout->Close();
    
    }
  }
 
 
  if ( fulldbg ) std::cout << __LINE__ << std::endl;
      
 
  if ( fref_ && !files_duplicated ) fref_->Close();
  if ( ftest_ ) ftest_->Close();
      
 
  if ( deleteref && !noref ) { 
    std::cout << "ref " << frefname << "\ttest " << ftestname << std::endl; 
    if ( frefname != ftestname && !files_duplicated ) { 
      std::string cmd = std::string("mv ") + frefname + " " + frefname + ".bak";
      std::system( cmd.c_str() );
      
      cmd = std::string("mv .newout.root ") + std::string(frefname);
      std::system( cmd.c_str() );
    }
    else {
      std::cerr << "reference file and test file are the same - not replacing" << std::endl;
    }
  }  
 
  std::cout << "deleting at " << __LINE__ << std::endl;
 
  if ( fref_ && !files_duplicated ) delete fref_;
  if ( ftest_ ) delete ftest_;
 
  return 0;
}