Index: appl-conv/trunk/src/convolute.cxx
===================================================================
--- appl-conv/trunk/src/convolute.cxx	(revision 1825)
+++ appl-conv/trunk/src/convolute.cxx	(revision 1826)
@@ -1,1845 +1,1845 @@
 #include <iostream>
 #include <iomanip>
 #include <string>
 #include <map>
 #include <vector>
 #include <cstdlib>
 
 #include <stdio.h>
 
 #include "amconfig.h"
 
 
 #include "appl_grid/appl_grid.h"
 #include "appl_grid/Directory.h"
 
 
 #include "appl_grid/appl_timer.h"
 #include "appl_grid/lumi_pdf.h"
 
 #include "TCanvas.h"
 #include "TH1D.h"
 #include "TFile.h"
 #include "TLine.h"
 #include "TPad.h"
 #include "TStyle.h"
 #include "TDirectory.h"
 
 // #include "hoppet_v1.h"
 
 #include "DrawLabel.h"
 // #include "Atlas/AtlasStyle.h"
 
 #include "Legend.h"
 
 
 // lhapdf routines
 // #include "LHAPDF/LHAPDF.h"
 // extern "C" void evolvepdf_(const double& x, const double& Q, double* xf); 
 // extern "C" double alphaspdf_(const double& Q);
 
 
 #include "LHAPDF/PDFSet.h"
 #include "LHAPDF/PDF.h"
 
 /// allow different pdfs for each target hadron
 
 bool use_photon = false;
 
 LHAPDF::PDFSet* pdfset0 = 0;
 LHAPDF::PDFSet* pdfset1 = 0;
 
 LHAPDF::PDF* pdf0 = 0;
 LHAPDF::PDF* pdf1 = 0;
 
 
 /// pdf and lphas wrappers to allow multiple simultaneous pdfs
 
 
 // extern "C" void    evolvepdf_(const double& x, const double& Q, double* xf);
 // extern "C" double  alphaspdf_(const double& Q);
 
 
 extern "C" void evolvepdf0(const double& x, const double& Q, double* xf) {
   if ( !pdf0 ) return; 
   xf += 6;
   for ( int i=-6 ; i<=6 ; i++ ) xf[i] = pdf0->xfxQ( i, x, Q );
   if ( use_photon ) xf[13] = pdf1->xfxQ( 22, x, Q );
 }
 
 
 extern "C" void _evolvepdf1(const double& x, const double& Q, double* xf) {
   if ( !pdf1 ) return; 
   xf += 6;
   for ( int i=-6 ; i<=6 ; i++ ) xf[i] = pdf0->xfxQ( i, x, Q );
   if ( use_photon ) xf[13] = pdf1->xfxQ( 22, x, Q );
 }
 
 extern "C" void (*evolvepdf1)(const double& x, const double& Q, double* xf) = _evolvepdf1;
 
 
 extern "C" double alphaspdf(const double& Q) {
   if ( pdf0 ) return pdf0->alphasQ(Q);
   return 0;
 }
 
 
 bool DBG = false;
 
 std::string watermark = "";
 
 double maxratio = 0;
 bool      quiet = false;
 bool  no_output = false;
 
 // static const int nFlavours = 5;
 
 
 std::string date() { 
   time_t _t;
   time(&_t);
   std::string a = ctime(&_t);
   std::string b = "";
   for ( unsigned i=0 ; i<a.size()-1 ; i++ ) b+=a[i];
   return b;
 }
 
 
 void ranger( TH1* h, TH1* h0 ) { 
   int ibin = 0;
   for ( int i=0 ; i<=h->GetNbinsX() ; i++ ) { 
     if ( h->GetBinContent(i) != 0 ) break;
     ibin++;
   }
   if ( ibin>1 ) { 
     h->GetXaxis()->SetRange( ibin, h->GetNbinsX() );
     h0->GetXaxis()->SetRange( ibin, h->GetNbinsX() );
   }
 }
 
 bool contains( const std::string& s, const std::string& p ) { return s.find(p)!=std::string::npos; } 
 
 
 std::string fname = "";
 
 bool make_png= false;
 
 int max_ratio = 7;
 
 void smooth(TH1D* h) { 
   for ( int i=2 ; i<h->GetNbinsX() ; i++ ) {
     double y0 = h->GetBinContent(i-1);
     double y  = h->GetBinContent(i);
     double y1 = h->GetBinContent(i+1);
 
     if ( std::fabs(y-y0)>0.2 ) { 
       y = 0.5*(y0+y1);
       h->SetBinContent( i,  y );
     }
   }
 }
 
 
 void binwidth(TH1D* h) { 
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
     double dx = h->GetBinLowEdge(i+1) - h->GetBinLowEdge(i);
     h->SetBinContent( i,  h->GetBinContent(i)/dx );
     h->SetBinError( i,  h->GetBinError(i)/dx );
   }
 }
 
 
 double getRealMinimum( TH1D* h ) { 
   double _min = 0;
   bool _set = false;
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
     double y = h->GetBinContent(i);
     if ( y!=0 ) { 
       if ( _set ) { 
 	if ( y<_min ) _min = y;
       }
       else { 
 	_min = y;
 	_set = true;
       }
     }
   }
 
   return _min;
 }
 
 
 double getRealMaximum( TH1D* h ) { 
   double _max = 0;
   bool _set = false;
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
     double y = h->GetBinContent(i);
     if ( y!=0 ) { 
       if ( _set ) { 
 	if ( y>_max ) _max = y;
       }
       else { 
 	_max = y;
 	_set = true; 
       }
     }
   }
 
   return _max;
 }
 
 
 double scale_factor = 1;
 
 TH1D* divide( const TH1D* h1, const TH1D* h2 ) {
 
   if ( h1==NULL || h2==NULL ) return NULL;
  
   TH1D* h = (TH1D*)h1->Clone();
 
   if ( DBG ) std::cout << "histograms " << h1->GetTitle() << " " << h2->GetTitle() << std::endl;
 
   
   bool unset = true;
 
   double hmin = 0;
   double hmax = 0;
   
   
   for ( int i=1 ; i<=h1->GetNbinsX() ; i++ ) { 
     double b  = h2->GetBinContent(i);
     //    double be = h2->GetBinError(i);
     double t  = h1->GetBinContent(i);
     //  double te = h1->GetBinError(i);
 
     double r  = ( b!=0 ? t/b : 0 );
     //    double re = ( b!=0 ? sqrt((r+1)*r/b) : 0 );
     double re = 0;
     
     h->SetBinContent( i, r );
     h->SetBinError( i, re ) ;
 
     if ( unset && b!=0 ) {  
       hmin=r;
       hmax=r;
       unset = false;
     }
 
     if ( b!=0 ) { 
       if ( hmin>r ) hmin=r;
       if ( hmax<r ) hmax=r; 
     }
 
     //    if ( DBG ) std::cout << "\tx=" << h->GetBinCenter(i) << "\tratio=" << r << std::endl;
   } 
 
 
   if ( DBG ) std::cout << "\tmin ratio = " << hmin << "\tmax ratio = " << hmax << std::endl;
   
   double _max = getRealMaximum( h );
   double _min = getRealMinimum( h );
 
   if ( DBG ) std::cout << "\tmin ratio = " << _min << "\tmax ratio = " << _max << std::endl;
   
   double mindelta = 0;
 
   std::string bestfile = "";
 
 
   char duff[64];
 
   for ( double delta=100 ; delta>1e-7 ; delta*=0.1 ) {
 
     if ( !quiet ) { 
 
       h->SetMaximum(1+delta);
       h->SetMinimum(1-delta);
       
       /// aha!! sutt debug code !!!
       h->DrawCopy();
       
       //    if ( watermark=="" ) watermark = label( "applgrid %s compiled %s", appl::version().c_str(), appl::compiled().c_str() );
       
       //    DrawLabel( 0.02, 0.02, label( "applgrid %s compiled %s", appl::version().c_str(), appl::compiled().c_str() )  ); 
       DrawLabel( 0.02, 0.02, watermark );
             
       sprintf(duff,"ratio-%f.pdf", delta );
       gPad->Print(duff);      
       
       if ( make_png ) {
 	sprintf(duff,"ratio-%f.png", delta );
 	gPad->Print(duff);  
       }    
     }
     
     if ( _min>(1-1*delta) && _max<(1+1*delta) ) {     
       mindelta = delta;
       if ( !quiet ) bestfile = duff;
     }
 
   }
 
   if ( maxratio && maxratio < mindelta ) mindelta = maxratio;
   h->SetMaximum( 1+mindelta );
   h->SetMinimum( 1-mindelta );
   
 
   scale_factor = 1;
 
 
   std::cout << "best agreement to within " << mindelta << "\tfile: " << bestfile << std::endl;
 
 
   if ( mindelta>1.5e-2 ) { 
     std::cerr << "WARNING!! poor convolution agreement " << fname << " " << mindelta*100 << "%" << std::endl;  
 
     if ( mindelta>0.05 ) {
       if ( std::fabs(_min-1)>std::fabs(_max-1) ) scale_factor = _min;
       else                                        scale_factor = _max;
     } 
 
   }
 
   return h;
 }
 
 
 
 template<typename T>
 T max( const std::vector<T>& v ) { 
   if ( v.empty() ) return 0;
   T _max = v.back();
   for ( int i=v.size() ; i-- ; ) if ( _max<v[i] ) _max = v[i];
   return _max;
 } 
 
 
 template<typename T>
 T min( const std::vector<T>& v ) { 
   if ( v.empty() ) return 0;
   T _min = v.back();
   for ( int i=v.size() ; i-- ; ) if ( _min>v[i] ) _min = v[i];
   return _min;
 } 
 
 
 
 template<typename T>
 void Fill( TH1D* h, const std::vector<T>& v ) { 
   for ( int i=v.size() ; i-- ; ) h->Fill( v[i] );
 }
 
 
 /// not an actual chi2, but average squared 
 /// fractional difference  
 
 double chi2( TH1D* h0, TH1D* h1 ) { 
   if ( h0==0 || h1==0 ) return 0;
   if ( h0->GetNbinsX()!=h1->GetNbinsX() ) return 0;
 
   double c2 = 0;
 
   int ibins = 0;
 
   for ( int i=0 ; i<h0->GetNbinsX() ; i++ ) { 
     double d0 = h0->GetBinContent(i+1);
     double d1 = h1->GetBinContent(i+1);
     double c = (d0-d1);
     if ( d0!=0 ) {
       ibins++;
       c2 += (c/d0)*(c/d0);
     }
   }
   
   if ( ibins>0 ) return std::sqrt(c2/ibins);
 
   return 0;
 }
 
 
 void processphoton( std::map<std::string,std::vector<int> >&  mLO,
 		    std::map<std::string,std::vector<int> >& mNLO );
 
 
 void photon_subprocesses( appl::grid& g, const std::string& basename, double rscale=1, double fscale=1 );
 
 
 /// return the head of a full path
 std::string head( std::string s ) { 
   std::string h;
   while( s.find("/")!=std::string::npos ) { 
     h+=s.substr(0,s.find("/")+1);
     s.erase(0,s.find("/")+1);
   }
   if ( h.size() ) return h;
   else            return s;
 }
 
 
 /// return the tail of a full path
 std::string tail( std::string s ) {
   while( s.find("/")!=std::string::npos ) s.erase(0,s.find("/")+1);
   return s;
 }
 
 
 
 int usage( int i=0, std::ostream& s=std::cout ){ 
   s << "Usage: appl-conv [OPTIONS] -o  input_grid.root [input_grid1.root ... input_gridN.root]\n\n";
   s << "  appl-conv performs the convolution for the grids entered on the command line and creates many output plots and a root file\n\n";
   s << "Options: \n";
   s << "    -o,  --output value\t  name for output root file (optional)\n";
   s << "    -x                 \t  do not write output root file\n\n";
   s << "    -rs, --rscale value\t  use a renormalisation scale factor of value \n";
   s << "    -fs, --fscale value\t  use a factorisation scale factor of value \n";
   s << "    -p,  --pdf    value\t  use the the particular set from the specified file\n";
   s << "    -i,  --iset   value\t  use the particular pdf value from the full set\n\n";
   s << "    -p2, --pdf2   value\t  use the the particular set from the specified file for any second input hadron\n";
   s << "    -i2, --iset2  value\t  use the particular pdf value from the full set\n\n";
   s << "         --rfile  value\t  use the reference from file value instead of from the grid\n\n";
   s << "    -d            value\t  for a pt or scale dependent x axis variable, use a dynamic scale for the convolution in each bin,\n";
   s << "    -sp               \t  plot of the subprocess contributions\n";
   s << "    -c               \t  apply the multiplicative corrections from the grid\n";
   //  s << "    -chi2           \t  report the chi2 value (no uncertainties)\n";
   s << "    -nt, --nothreads\t  disable threads\n";
   s << "\n";
   s << "    -logx           \t  use log x scale\n";
   s << "    -linx           \t  use lin x scale\n";
   s << "    -logy           \t  use log y scale\n";
   s << "    -e              \t  show uncertainties\n";
   s << "    -xrange xmin xmax\t  use xmin and xmax x-axis limits\n";
   s << "    -rp             \t  draw a plot of the cross section, and ratio\n";
   s << "    -m,  --order value\t  order of calculation for the panel\n";
   s << "    -u,          value\t  in the panel ratio plot, limit to larger than value powers of ten\n";
   s << "    -png            \t  also make png ratio plots\n";
   s << "    -q,  --quiet    \t  don't actually make any plots\n";
   s << "    -r              \t  use max ratio\n";
   s << "\n";
   s << "    -lc             \t  print out the full lumi pdf configuration\n";
   s << "    -t, --title     \t  print out the grid documentation (title) information only (do not perform the convolution)\n";
   s << "    -w, --watermark \t  watermark to be displayed\n";  
   s << "    -v, --verbose   \t  display verbose output\n";
   s << "    -f, --fixed     \t  if verbose printout selected also print out the fixed order contrbutions\n";
   s << "\n";
   s << "    -h, --help      \t  display this help\n";
   s << "\nSee " << PACKAGE_URL << " for more details\n"; 
   s << "\nReport bugs to <" << PACKAGE_BUGREPORT << ">";
   s << std::endl;
   return i;
 }
 
 
 void fonts() { 
   gStyle->SetLabelFont(42,"x");
   gStyle->SetTitleFont(42,"y");
 
   gStyle->SetLabelSize(12,"x");
   gStyle->SetTitleSize(18,"y");
  
 }
 
 
 void fonts(TH1D* h) { 
 
   h->SetLabelFont(43,"x");
   h->SetLabelSize(14,"x");
 
   h->SetLabelFont(43,"y");
   h->SetLabelSize(14,"y");
 
   h->SetTitleFont(43,"x");
   h->SetTitleSize(18,"x"); 
 
   h->SetTitleFont(43,"y");
   h->SetTitleSize(18,"y"); 
 
   h->GetXaxis()->SetTitleOffset(2.3);
   h->GetYaxis()->SetTitleOffset(2.2);
 
 }
 
 void pads(TPad* tp1, TPad* tp2 ) { 
   tp2->SetMargin(0.12, 0.03, 0.05, 0.03);
   tp1->SetMargin(0.12, 0.03, 0.17, 0.025);
 }
 
 
 double nonZeroMax( TH1D* h ) { 
   bool first = true;
   //  double hmin = 0;
   double hmax = 0;
 
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
     double y = h->GetBinContent(i);
     if ( y!=0 ) { 
       if ( first || y>hmax ) hmax = y;
       first = false;
     } 
   }
   return hmax;
 }  
 
 
 
 double nonZeroMin( TH1D* h ) { 
   bool first = true;
   //  double hmin = 0;
   double hmin = 0;
 
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) { 
     double y = h->GetBinContent(i);
     if ( y!=0 ) { 
       if ( first || y<hmin ) hmin = y;
       first = false;
     } 
   }
   return hmin;
 }  
 
 
 
 
 TCanvas* plotRatio( TH1D* _h0, TH1D* _h1, TH1D* _hratio, const std::string& options="", const std::string& poptions="" ) {  
     
   //  TH1D* hratio = (TH1D*)h0->Clone("hratio");
   //  hratio->SetDirectory(0);
 
   //  hratio->Divide( h1 );
 
   gStyle->SetPadRightMargin(0.05);
   gStyle->SetPadLeftMargin(0.1);
   
 
   
   TH1D* h0 = (TH1D*)_h0->Clone("h0");
   h0->SetDirectory(0);
   //  h0->SetTitle(";;cross section");
 
   TH1D* h1 = (TH1D*)_h1->Clone("h1");
   h1->SetDirectory(0);
   //  h1->SetTitle(";;cross section");
 
   TH1D* hratio = (TH1D*)_hratio->Clone("hratio");
   hratio->SetDirectory(0);
   hratio->GetYaxis()->SetTitle("ratio");
 
   //  TLine* tl = new TLine( hratio->GetBinLowEdge(1), 1, hratio->GetBinLowEdge( hratio->GetNbinsX()+1 ), 1 );
   TLine* tl = new TLine( hratio->GetBinLowEdge( hratio->GetXaxis()->GetFirst() ), 1, 
 			 hratio->GetBinLowEdge( hratio->GetXaxis()->GetLast()+1 ), 1 );
   tl->SetLineStyle(3);
   
 
 
   h0->GetXaxis()->SetLabelOffset(5);
   h0->GetXaxis()->SetTitleOffset(5);
 
 
   hratio->GetXaxis()->SetTitleOffset(1.5);
 
   fonts( h0 );
   fonts( h1 );
   fonts( hratio );
 
 
   //  hratio->SetMinimum(0.99);
   //  hratio->SetMaximum(1.01);
 
   std::string _options = poptions;
   std::string mname = h0->GetName();
 
   TCanvas* tc = new TCanvas( "tc", "", 550, 700 );
   //  tc->SetDirectory(0);
 
   //  tc->SetFillStyle(4000);
   tc->SetFillStyle(-4000);
 
   tc->SetFillColor(kWhite);
   tc->Draw();
   tc->cd();
 
   
 
 
   //  std::cout << "\tcreating TPads " << std::endl; 
   
   TPad* tp1 = new TPad(("xy-"+mname).c_str(), " ",  0.0, 0, 1.0, 0.413 );  
   //  tp1->SetDirectory(0);
   // tp1 = new TPad(("xy-"+mname).c_str(), " ",  0.0, 0, 1.0, 1.0 );      
   tp1->SetFillStyle(4000);
   //  tp1->SetLeftMargin(0.18);
   // tp1->SetLeftMargin(0.1);
 
 
   //  std::cout << "\tTPad name " << tp1->GetName() << " " << tp1 << std::endl; 
   
   //      TCanvas* duff = new TCanvas((mname+"duff").c_str(), "event display", (int)(mscale*4200), (int)(mscale*1400));
   
   // duff->SetFillStyle(4000);
   // duff->Draw();
   // duff->cd();
   
   TPad*  tp2 = new TPad(("rz-"+mname).c_str(), " ",  0, 0.383, 1, 1 );
   // tp2->SetDirectory(0);
 
   tp2->SetFillStyle(4000);
   //  tp2->SetRightMargin(0.05);
 
   pads( tp1, tp2 );
 
 
   tp1->Draw();  
   tp2->Draw();
   
 
   h0->SetTitle("");
 
   gStyle->SetErrorX(0);
 
 
   //  std::cout << "OPTIONS " << _options << std::endl;
   
   std::string extra = "";
 
   tp2->cd();
   h0->DrawCopy( (_options+extra).c_str() );
   h1->DrawCopy( ("same"+_options+extra).c_str() );
 
   //  h0->SetLineColor(kWhite);
   h0->DrawCopy( ("same"+_options+extra).c_str() );
 
   h0->SetLineStyle(1);
   h0->SetLineColor(kBlack);
 
 
   if ( contains( options, "logy" ) ) tp2->SetLogy(true);
   else                               tp2->SetLogy(false);
  
   if ( contains( options, "logx" ) ) { 
     tp2->SetLogx(true);
     tp1->SetLogx(true);
   }
   else { 
     tp2->SetLogx(false);
     tp1->SetLogx(false);
   }
 
   tp1->SetLogy(false);
 
 
   Legend* leg = new Legend( 0.5, 0.9, 0.7, 0.9 );
 
   h0->SetLineStyle(2);
   leg->AddEntry( h0, "Reference", "l");
   leg->AddEntry( h1, "Fast convolution", "l");
 
   //  leg->Draw( 0.2, 0.4 ); 
   
   //  std::cout << "leg" << std::endl;
 
   if ( contains( options, "logy" ) ) leg->Draw( 0.2, 0.3 ); 
   else                               leg->Draw( 0.35, 0.3 ); 
 
 
 
   tp1->cd();
 
 
   double ymax = nonZeroMax( hratio ); 
   double ymin = nonZeroMin( hratio ); 
   double maxdelta = 100;
 
   if ( max_ratio>0 ) max_ratio *= -1;
 
   for ( int ip=2 ; ip>=max_ratio ; ip-- ) { 
     double delta = std::pow( 10., ip );
     if ( ymax<(1+delta) && ymin>(1-delta) ) maxdelta = delta;
   }
 
   //  maxdelta = 0.01;
 
   if ( maxdelta!=0 ) {
     // if ( maxdelta<0.01 ) maxdelta=0.01; 
     hratio->SetMinimum(1-maxdelta);
     hratio->SetMaximum(1+maxdelta);
     hratio->DrawCopy( _options.c_str() );
     tl->Draw();
   }
 
   tc->cd();
 
   return tc;
 } 
 
 
 
 
 /// replace the annoying LHAPDF Config destructor  ...
 
 #include "LHAPDF/Config.h"
 
 LHAPDF::Config::~Config() {}
 
 
 /// why was this needed ?? to run interactively ??? 
 /// #include "TApplication.h"
 
 
 void print( TH1D* xsec, TH1D* reference ) { 
 
       std::cout << "xsection has " << xsec->GetNbinsX() << " bins" << std::endl;
       for ( int ih=0 ; ih<xsec->GetNbinsX() ; ih++ ) { 
 	std::cout << "\txsec[" << ih << "]";
 	if ( ih<10  && xsec->GetNbinsX()>9 )  std::cout << " "; 
 	if ( ih<100 && xsec->GetNbinsX()>99 ) std::cout << " ";
  
 
 	double delta = 0;
 	double ratio = 0;
 	if ( reference->GetBinContent(ih+1)!=0 ) { 
 	  ratio = xsec->GetBinContent(ih+1)/reference->GetBinContent(ih+1);
 	  delta = ratio-1;
 	}
 
 #if 0
 	double lnr = std::log( rscale*rscale );
 	
 	static const int nc = 3;
 	//TC   const int nf = 6;
 	static const int nf = 5;
 	static double beta0=(11.*nc - 2.*nf)/6.;
 
 	static double beta1=(34.*nc*nc - 3.*(nc-1./nc)*nf - 10.*nc*nf)/12.;
   
 	double fac0 = beta0*amu*lnr;
 
 	double fac1 = beta1*amu*amu*lnr;
 	
 	double fac2 = fac0*fac0;
 #endif
 
 	std::printf( " = %15.10lf    %15.10lf     (%14.10lf) \t(%14.10lf) \t(%12.10lf)",  
 		     xsec->GetBinContent(ih+1),
 		     reference->GetBinContent(ih+1),
 		     xsec->GetBinContent(ih+1)-reference->GetBinContent(ih+1),
 		     delta,	
 		     ratio );
 
 #if 0
 	if ( DBG ) std::printf( " ::  %15.10lf %15.10lf %15.10lf  ::  f0: %15.10lf  f1 %15.10lf   sf1 %15.10lf  %15.10lf  %15.10lf",  
 		     xsec->GetBinContent(ih+1)/amu,
 		     xsec->GetBinContent(ih+1)/amu/amu,
 		     xsec->GetBinContent(ih+1)/amu/amu/amu,
 		     fac0, fac1, 
 		     fac0*xsec->GetBinContent(ih+1), 
 		     fac1*xsec->GetBinContent(ih+1), 
 		     fac2*xsec->GetBinContent(ih+1)
 		     // (1+fac0)*xsec->GetBinContent(ih+1)
 		     );
 #endif
         
 	std::printf( "\n" );
      
 
 #if 0
 	std::cout << " = " << std::fixed 
 		  <<          std::setprecision(11) << std::fixed << xsec->GetBinContent(ih+1) 
 		  << "   " << std::setprecision(11) << std::fixed << reference->GetBinContent(ih+1) 
 		  << "\t(" << std::setprecision(11) << std::fixed << xsec->GetBinContent(ih+1)-reference->GetBinContent(ih+1) << ")"
 		  << "\t(" << delta << ")"
 		  << "\t(" << ratio << ")" << std::endl;
 #endif
 
       }
 
 }
 
 
 void print_contributions( TH1D* LOxsec, TH1D* NLOxsec=0, TH1D* NNLOxsec=0, TH1D* reference=0  ) { 
 
   TH1D* xsec = 0;
 
   if      ( NNLOxsec !=0 ) xsec = NNLOxsec;
   else if (  NLOxsec !=0 ) xsec =  NLOxsec;
   else                     xsec =   LOxsec;
 
   std::cout << "xsection has " << xsec->GetNbinsX() << " bins" << std::endl;
   if (   LOxsec!=0 ) std::cout << "                        LO";
   if (  NLOxsec!=0 ) std::cout << "                  NLO";
   if ( NNLOxsec!=0 ) std::cout << "                 NNLO";
   if ( reference!=0 )  std::cout << "                reference      (     ratio     )";
   
   std::cout << std::endl;
 
   for ( int ih=0 ; ih<xsec->GetNbinsX() ; ih++ ) { 
 	std::cout << "\txsec[" << ih << "]";
 	if ( ih<10  && xsec->GetNbinsX()>9 )  std::cout << " "; 
 	if ( ih<100 && xsec->GetNbinsX()>99 ) std::cout << " ";
  
 	std::printf( " = %15.10lf", LOxsec->GetBinContent(ih+1) );
 	if ( NLOxsec!=0 )  std::printf( "     %15.10lf", NLOxsec->GetBinContent(ih+1) );
 	if ( NNLOxsec!=0 ) std::printf( "     %15.10lf", NNLOxsec->GetBinContent(ih+1) );
 
 	if ( reference!=0 ) { 
 	  double r = reference->GetBinContent( ih+1 );
 	  double x = xsec->GetBinContent( ih+1 );
 	  if ( r!=0 ) printf( "     %15.10lf     (%15.10lf)", r, x/r );
 	}
 
 	std::printf( "\n" );
      	
   }
 
 }
 
 
 
 
 int main(int argc, char** argv) { 
 
   if ( argc < 2 ) return usage(-1); 
 
   std::cout << argv[0] << ": compiled " <<  __DATE__ << "\tapplgrid " << appl::version() << " compiled " << appl::compiled() << std::endl;
   
   //  SetAtlasStyle();
 
   //  fonts();
 
   //  if ( argc<2 ) return -1;
 
   gStyle->SetOptStat(0);
   gStyle->SetPadRightMargin(0.07);
   gStyle->SetPadTopMargin(0.07);
   gStyle->SetPadLeftMargin(0.12);
   gStyle->SetPadBottomMargin(0.12);
 
 
   std::vector<std::string> grids;
 
   bool prepare_subproc = false;
 
   double dynamicScale = 0;
 
   bool corrections = false;
 
   bool verbose = false;
 
   double fscale = 1;
   double rscale = 1;
 
   std::string xaxis = "";
   std::string yaxis = "";
 
   bool logx = false;
   bool logy = false;
 
   double xmin = 0;
   double xmax = 0;
 
   bool lumi_check  = false;
   bool ratio_panel = false;
   bool dochi2      = true;
   bool custompdf   = false;
   bool title_only  = false;
   bool show_errors = false;
 
   std::string timefile = "";
 
   std::string outfile = "";
 
   // setup pdf
 
   std::string _pdfname  = "";
   std::string _pdfname1 = "";
   std::string  pdfpath;
 
   if ( std::string(PDFPATH)!="" ) { 
     pdfpath = std::string(PDFPATH) + "/";
     _pdfname += pdfpath;
   }  
 
   //  std::cout << "pdfpath >" << pdfpath << "<" << std::endl; 
 
 
   _pdfname += "";
 
   int ipdf  = 0;
   int ipdf1 = 0;
  
   int user_order = -999;
 
   bool terms = false; /// fixed order terms
 
   std::string reffile = "";
 
   for ( int ia=1 ; ia<argc ; ia++ ) {  
     
     std::string arg = argv[ia];
 
     if      (   arg=="-h" || arg=="--help" ) return usage(0);
     else if ( ( arg=="-p"  || arg=="--pdf"  ) && ia<(argc-1) ) { custompdf=true; _pdfname = pdfpath + argv[++ia]; }  
     else if ( ( arg=="-p2" || arg=="--pdf2"  ) && ia<(argc-1) ) { custompdf=true; _pdfname1 = pdfpath + argv[++ia]; }  
     else if ( ( arg=="-i"  || arg=="--iset" ) && ia<(argc-1) ) { custompdf=true;  ipdf = std::atoi(argv[++ia]); }
     else if ( ( arg=="-i2" || arg=="--iset2" ) && ia<(argc-1) ) { custompdf=true;  ipdf1 = std::atoi(argv[++ia]); }
     else if (   arg=="-d" && ia<(argc-1) ) dynamicScale = std::atof(argv[++ia]);  
     else if (   arg=="-sp" ) prepare_subproc = true;  
     else if (   arg=="-c"  ) corrections = true;  
     else if (   arg=="-f" ||  arg=="--fixed"  )  terms = true;
     else if (   arg=="-v"  || arg=="--verbose" ) verbose = true;
     else if ( ( arg=="-rs" || arg=="--rscale" ) && ia<(argc-1) ) rscale = std::atof(argv[++ia]);  
     else if ( ( arg=="-fs" || arg=="--fscale" ) && ia<(argc-1) ) fscale = std::atof(argv[++ia]);  
     else if ( ( arg=="-m" ||  arg=="--order"  ) && ia<(argc-1) ) user_order = std::atoi(argv[++ia]);  
     else if ( ( arg=="-o" || arg=="--output"  ) && ia<(argc-1) ) outfile = argv[++ia];  
     else if (   arg=="--rfile" ) { 
       if ( ++ia<argc ) reffile = argv[ia]; 
       else usage(-1);
     }
     else if (   arg=="-x" ) no_output = true;  
     else if (   arg=="-q"  || arg=="--quiet"  ) quiet = true;  
     else if (   arg=="-logx" ) logx = true;
     else if (   arg=="-logy" ) logy = true;
     else if (   arg=="-u" ) { 
       if ( ++ia<argc ) max_ratio = std::atoi(argv[ia]); 
       else usage(-1);
     }
     else if (   arg=="-rp" ) ratio_panel = true;
     else if (   arg=="-e"  ) show_errors = true;
     else if (   arg=="-lc" ) lumi_check = true;
     else if (   arg=="-t"  || arg=="--title") title_only = true;
     else if (   arg=="-chi2" ) dochi2 = true;
     else if (   arg=="-png"  ) make_png = true;
     else if (   arg=="-dbg" ) DBG = true;
     else if (   arg=="-tf" ) { 
       if ( ++ia<argc ) timefile = argv[ia];
       else return usage(-1);
     }
     else if ( std::string(argv[ia])=="-xrange" ) { 
       if ( ia<(argc-1) ) xmin = std::atof(argv[++ia]);
       else return usage(-1);
       if ( ia<(argc-1) ) xmax = std::atof(argv[++ia]);
       else return usage(-1);
     }
     else if ( std::string(argv[ia])=="-nt" || std::string(argv[ia])=="--nothreads"  ) appl::grid::disable_threads(true);
     else if ( std::string(argv[ia])=="--xaxis" ) { 
       if ( ia<(argc-1) ) xaxis = argv[++ia];
       else return usage(-1);
     }
     else if ( std::string(argv[ia])=="--yaxis" ) { 
       if ( ia<(argc-1) ) yaxis = argv[++ia];
       else return usage(-1);
     }
     else if ( std::string(argv[ia])=="-w" || std::string(argv[ia])=="--watermark" ) { 
       if ( ia<(argc-1) ) watermark = argv[++ia];
       else return usage(-1);
     }
     else if ( std::string(argv[ia])=="-r" ) { 
       if ( ia<(argc-1) ) maxratio = std::atof(argv[++ia]);
       else return usage(-1);
     }
     else if ( std::string(argv[ia]).find("-")==0 ) { 
       std::cerr << "\nunknown option " << argv[ia] << std::endl;
       return usage(-1);
     }
     else grids.push_back(argv[ia]);
   }
   
   //  std::cout << "setup pdf" << std::endl; 
 
   bool pdfnotset=true;
 
   for ( unsigned ia=0 ; ia<grids.size() ; ia++ ) { 
 
     std::string gridName(grids[ia]);
     std::string basename  = tail(gridName).substr(0, gridName.find(".root") );
     std::string directory = head(gridName);
 
     fname = basename;
 
     if ( gridName.find("appl-conv.root") != std::string::npos && 
 	 gridName.find("appl-conv.root") == gridName.size()-14 ) { 
       std::cerr << "skipping convolution output file : " << gridName << "\n" ;
       continue;
     }
 
     /// check if actually a grid file ...
     { 
       TFile f(gridName.c_str());
       bool ok = false;
       if ( (TH1D*)f.Get("grid/reference") ) ok = true;
       f.Close();
       if ( !ok ) { std::cerr << "skipping non grid file : " << gridName << "\n"; continue; }
     }
 
     std::cout << "\n";
 
     try { 
       
       //    lumi_pdf::runlatex(true);
 
     appl::grid g( gridName.c_str() );
 
     use_photon = false;
-    if ( g.include_photon() ) use_photon = true;
+    //   if ( g.include_photon() ) use_photon = true;
 
     if ( title_only ) { std::cout << g.getDocumentation() << std::endl; return 0; }
 
     //    std::cout << g << std::endl;
 
 #if 1
     //    std::cout << "lumi pdf check: order " << g.nloops() << std::endl;
     
     for ( int i=0 ; i<=g.nloops() ; i++ ) { 
       //      std::cout << "printing lumi pdf:" << std::endl;
       
       const appl::appl_pdf* pdf = static_cast<const appl::grid*>(&g)->genpdf(i);
       
       if ( pdf ) { 
 	
 	const lumi_pdf* lpdf = dynamic_cast<const lumi_pdf*>( pdf );
 	
 	if ( lpdf ) {
 	  std::cout << lpdf->summary() << std::endl;
 	  if ( lumi_check ) std::cout << "lumi_pdf: " << *lpdf << std::endl;
 
 	  //	  latex( *pdf, "-pdf" );
 
 	}
 	else { 
 	  std::cout << "grid does not use lumi_pdf" << std::endl;
 	}
       }
       
     }
     //    std::cout << "lumi pdf check" << std::endl;
 #endif
 
     if ( pdfnotset ) { 
       if ( !custompdf ) { 
 	if (  g.getGeneratedPDF()!="" ) {
 	  _pdfname = g.getGeneratedPDF();
 	  ipdf     = g.getGeneratediPDF();
 	  std::cout << "reading pdf from grid: " << _pdfname << "  " << ipdf << std::endl;
 	}
 	else { 
 	  std::cerr << "no pdf specified" << std::endl;
 	  std::exit(-1);
 	}
       }
 
       std::cout << "setting up pdf set >" << _pdfname << "<" << " " << ipdf << std::endl; 
       
       /// don;t use the old lhaglue interface any longer
       /// have to fix this at some point
       /// LHAPDF::initPDFSet(_pdfname.c_str(), ipdf );
       
       pdfset0 = new LHAPDF::PDFSet( _pdfname );
       pdf0 = pdfset0->mkPDF( ipdf );
 
       if ( _pdfname1 == "" ) pdf1 = pdf0;
       else {
 	pdfset1 = new LHAPDF::PDFSet( _pdfname1 );
 	pdf1 = pdfset1->mkPDF( ipdf1 );
       }
 
       pdfnotset = false;
     }
     
 
     if ( dynamicScale ) g.setDynamicScale( dynamicScale );
     
     bool istrimmed = g.isTrimmed();
 
      
     g.untrim();
     int gusize = g.size();
 
     g.trim(); 
     int gtsize = g.size();
      
     
     
     if ( gtsize != gusize ) std::cout << "gridsize : " 
 				      << "untrimmed " << gusize 
 				      << "\ttrimmed " << gtsize 
 				      << "  (" << (gtsize*100./gusize) << "%)\tistrimmed " << istrimmed << std::endl;
     
     if ( verbose ) std::cout << g.getDocumentation() << std::endl;
 
     
     // get all the reference histograms
     
     TFile* f;
     f = new TFile(grids[ia].c_str());
     
 
     if ( outfile == "" ) {
       outfile = gridName.substr(0,gridName.find(".root"));
       outfile += std::string("_appl-conv") + ( corrections ? "_corr.root" : ".root" );
     }
    
     TFile* fout = 0;
 
     if ( !no_output ) fout = new TFile(outfile.c_str(), "recreate");
     
     
     Directory ref("reference");
     ref.push();
     
     TH1D* reference = 0;
     
     if ( reffile != "" && reffile != gridName ) {  
       TFile fref( reffile.c_str() );
       reference = (TH1D*)fref.Get("grid/reference")->Clone("reference");
       reference->SetDirectory(0);
       fref.Close();
     }
     else { 
       reference = (TH1D*)f->Get("grid/reference")->Clone("reference");
     }
 
     //  binwidth( reference );
     
     if ( !no_output ) reference->Write();
     ref.pop();
 
 
   
     // now calculate all the cross sections
   
 
 
     Directory xsDir("xSection");
     xsDir.push();
     
     int nLoops = g.nloops();
     
     if ( user_order != -999 && user_order <= g.nloops() ) nLoops = user_order;
 
     std::cout << "performing convolution at " << nLoops << " loop" << ( nLoops==1 ? "" : "s" ) << " ..." << std::endl; 
 
     TH1D* xsec = 0;
     
     if ( timefile =="" ) { 
 
       struct timeval mytimer = appl_timer_start();
 
       if ( rscale !=1 || fscale !=1 )  std::cout << "rscale: " << rscale << "\tfscale: " << fscale << "\tnloops: " << nLoops << std::endl; 
 
       xsec = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, nLoops, rscale, fscale ); 
 
       double mytime = appl_timer_stop(mytimer);
       std::cout << "done nloops: " << nLoops << "  (" << mytime << " ms)" << std::endl; 
 
     }
     else {
 
       std::vector<double> _times;
       
       for ( int it=200 ; it-- ;  ) {
 	struct timeval mytimer = appl_timer_start();
 	xsec = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, nLoops, rscale, fscale ); 
 	double mytime = appl_timer_stop(mytimer);
 	//	std::cout << "done (" << mytime << " ms)" << std::endl; 
 	_times.push_back( mytime );
       }
       
       double mx = max( _times );
       double mn = min( _times );
 
       std::cout << "times min " << mn << "\tmax " << mx << std::endl; 
       
       TDirectory* rene_cck = gDirectory;
       
       TH1D* htime = new TH1D( "ht", ";processing latency [ms];Entries", 10*int(mx*1.5), 0, int(mx*1.5) );
       
       Fill( htime, _times );
       
       TFile f( timefile.c_str(), "recreate" );
       htime->Write();
       f.Write();
       f.Close();
 
       rene_cck->cd();
 
     }
 
 
     //    std::cout << xsec << std::endl; 
 
     xsec->SetName("xsec");
     xsec->SetTitle(reference->GetTitle());
     
     xsec->SetLineColor(kBlue);
 
     xsec->GetXaxis()->SetMoreLogLabels(true);
 
     if ( !no_output ) xsec->Write();  
 
     std::vector<TH1D*> xsec_corr; 
 
     /// apply corrections if required 
 
     if ( corrections ) { 
 
       std::cout << "applying corrections" << std::endl;
       xsec_corr.resize( g.corrections().size() );
       for ( unsigned j=0 ; j<g.corrections().size() ; j++ ) { 
 	g.setApplyCorrection( j, true );
  
 	struct timeval mytimer = appl_timer_start();
 	
 	xsec_corr[j] = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, nLoops, rscale, fscale); 
 	
 	double mytime = appl_timer_stop(mytimer);
 
 	char hname[64];
 	sprintf( hname, "xsec-corr%02d-%s", j, g.correctionLabels()[j].c_str() );
 
 	xsec_corr[j]->SetName( hname );
 	xsec_corr[j]->GetXaxis()->SetMoreLogLabels(true);
 
 	if ( !no_output ) xsec_corr[j]->Write();
 
 	std::cout << "done (" << mytime << " ms)" << std::endl; 
        
       }
 
     }
 
 
     /// verbose printout of the cross section values if required 
     
     if ( dochi2 ) std::cout << "RMS difference: " << chi2( reference, xsec ) << std::endl; 
 
     
     if ( verbose && !terms ) print( xsec, reference );
 
     ranger( reference, xsec );
 
     xsec->DrawCopy();
 
     xsec->SetLineWidth(2);
 
     reference->SetLineColor(kBlack);
     reference->SetLineStyle(2);
     reference->GetXaxis()->SetMoreLogLabels(true);
     reference->DrawCopy("same");
 
     if ( show_errors ) reference->DrawCopy("samee1");
     
     
     for ( unsigned ic=0 ; ic<xsec_corr.size() ; ic++ ) { 
       xsec_corr[ic]->SetLineStyle(3);
       xsec_corr[ic]->DrawCopy("same");
     }
     
     //    if ( basename.find("pt")!=std::string::npos ) { 
     //      gPad->SetLogy(true);
     //      gPad->SetLogx(true);
     //    } 
 
     if ( !quiet ) { 
       //  gPad->Print("xsec.pdf");
       gPad->Print(("xsec"+basename+".pdf").c_str());
       gPad->Print((directory+"xsec"+basename+".pdf").c_str());
 
       if ( make_png ) { 
 	gPad->Print(("xsec"+basename+".png").c_str());
 	gPad->Print((directory+"xsec"+basename+".png").c_str());
       }    
     }
 
     //    std::cout << "head " << (directory+"xsec"+basename+".pdf") << std::endl;
 
     if ( logx ) gPad->SetLogx(true);
     else        gPad->SetLogx(false);
 
     /// for rations unset logy scale
     gPad->SetLogy(false);
 
     // now take all the ratios etc
     
     Directory ratiodir("ratio");
     ratiodir.push();
     
     // std::cout << "taking ratios" << std::endl;
     TH1D* ratio = divide( xsec, reference );
 
     if ( ratio ) {
       ratio->SetTitle("");
       ratio->SetName("ratio");
       if ( !no_output ) ratio->Write();
     }
 
     ratio->GetXaxis()->SetTitle( xsec->GetXaxis()->GetTitle() );
 
     if ( xaxis!="" ) { 
       reference->GetXaxis()->SetTitle( xaxis.c_str() );
       xsec->GetXaxis()->SetTitle( xaxis.c_str() );
       ratio->GetXaxis()->SetTitle( xaxis.c_str() );
     }
 
     if ( yaxis!="" ) { 
       reference->GetYaxis()->SetTitle( yaxis.c_str() );
       xsec->GetYaxis()->SetTitle( yaxis.c_str() );
       ratio->GetYaxis()->SetTitle( yaxis.c_str() );
     }
 
     std::string options = "";
 
     if ( contains( reference->GetXaxis()->GetTitle(), "p_{T}" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "P_{T}" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "pT" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "PT" ) ||
 	 contains( reference->GetXaxis()->GetTitle(), "pt" ) ||
 	 logy ) options += "logy";
 
     if ( contains( reference->GetXaxis()->GetTitle(), "p_{T}" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "P_{T}" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "pT" ) || 
 	 contains( reference->GetXaxis()->GetTitle(), "PT" ) ||
 	 contains( reference->GetXaxis()->GetTitle(), "pt" ) || 
 	 logx ) options += "logx";
 
     if ( !quiet ) { 
       if ( ratio_panel ) { 
 	if ( xmin!=0 || xmax!=0 ) { 
 	  reference->GetXaxis()->SetRangeUser( xmin, xmax );
 	  xsec->GetXaxis()->SetRangeUser( xmin, xmax );
 	  ratio->GetXaxis()->SetRangeUser( xmin, xmax );
 	}
 	std::string poptions = "hist";
 	if ( show_errors ) poptions = "histe1"; 
 	plotRatio( reference, xsec, ratio, options, poptions );
 
 	gPad->Print( ("ratio-panel-"+basename+".pdf").c_str() );
       }
 
       ratio->DrawCopy();
       if ( watermark=="" ) watermark = label( "applgrid %s compiled %s", appl::version().c_str(), appl::compiled().c_str() );
       DrawLabel( 0.02, 0.02, watermark );
       gPad->Print( ("ratio-"+basename+".pdf").c_str() );
       gPad->Print( (directory+"ratio-"+basename+".pdf").c_str() );
 
       if ( make_png ) gPad->Print( ("ratio-"+basename+".png").c_str() );
     }    
 
     if ( corrections ) { 
       std::cout << "\ntaking ratios (with corrections)" << std::endl;
       for ( unsigned j=0 ; j<xsec_corr.size() ; j++ ) { 
 	TH1D* ratio = divide( xsec_corr[j], reference );
 	ratio->SetName( ("ratio-" + g.correctionLabels()[j]).c_str() );
 	if ( !no_output ) ratio->Write();
 
 	if ( !quiet ) { 
 	  ratio->DrawCopy();
 
 	  gPad->Print( ("ratio-"+basename+g.correctionLabels()[j]+".pdf").c_str() );
 	  gPad->Print( (directory+"ratio-"+basename+g.correctionLabels()[j]+".pdf").c_str() );
 
 	  if ( make_png ) gPad->Print( ("ratio-"+basename+g.correctionLabels()[j]+".png").c_str() );
 	}
       }
     } 
     
     ratiodir.pop();
     
     if ( quiet && no_output && !terms ) return 0;
 
     //    std::cout << "individual components..." << std::endl;
     
     if ( user_order != -999 && user_order <= g.nloops() ) return 0;
 
     TH1D* hLOxsec    = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, 0, rscale, fscale); 
     hLOxsec->SetTitle("");
     hLOxsec->SetName("xsecLO");
 
     if ( !no_output ) hLOxsec->Write();
 
     hLOxsec->GetXaxis()->SetMoreLogLabels(true);
 
     TH1D* hNLOxsec   = 0;
     TH1D* hNNLOxsec  = 0;
 
     /// NNLO (if available)
     TH1D* hNLOOxsec  = 0;
     TH1D* hNNLOOxsec = 0;
     
     if ( g.nloops()>0 ) { 
 
       hNLOxsec  = g.convolute(evolvepdf0, evolvepdf1, alphaspdf,  1, rscale, fscale); 
       hNLOOxsec = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, -1, rscale, fscale); 
 
       hNLOxsec->SetTitle("");
       hNLOOxsec->SetTitle("");
 
       hNLOxsec->SetName("xsecNLO");
       hNLOOxsec->SetName("xsecNLOO");
 
       if ( !no_output ) { 
 	hNLOxsec->Write();
 	hNLOOxsec->Write();
       }
 
       hNLOxsec->GetXaxis()->SetMoreLogLabels(true);
       hNLOOxsec->GetXaxis()->SetMoreLogLabels(true);
       
       //      std::cout << "LOOPS " << g.nloops() << std::endl;
 
       if ( g.nloops()>1 ) { 
 
 	//	std::cout << "x-section at " << g.nloops() << ( g.nloops()>1 ? " loops" : "loop" ) << std::endl;
 
 	hNNLOxsec  = g.convolute(evolvepdf0, evolvepdf1, alphaspdf,  2, rscale, fscale); 
 	hNNLOOxsec = g.convolute(evolvepdf0, evolvepdf1, alphaspdf, -2, rscale, fscale); 
 
 	hNNLOxsec->SetTitle("");
 	hNNLOOxsec->SetTitle("");
 	
 	hNNLOxsec->SetName("xsecNNLO");
 	hNNLOOxsec->SetName("xsecNNLOO");
 
 	if ( !no_output ) {
 	  hNNLOxsec->Write();
 	  hNNLOOxsec->Write();
 	}
 
 	hNNLOxsec->GetXaxis()->SetMoreLogLabels(true);
 	hNNLOOxsec->GetXaxis()->SetMoreLogLabels(true);
       }
     }
 
     if ( verbose && terms ) print_contributions( hLOxsec, hNLOxsec, hNNLOxsec, reference );
 
     if ( quiet && no_output ) return 0;
 
 
     if ( !quiet ) { 
       hLOxsec->DrawCopy();
     
       if ( basename.find("pt")!=std::string::npos ) { 
 	gPad->SetLogy(true);
 	gPad->SetLogx(true);
       }    
       
       if ( logx ) gPad->SetLogx(true);
       if ( logy ) gPad->SetLogy(true);
       
       
       // gPad->Print("LO.pdf");
       gPad->Print(("LO"+basename+".pdf").c_str());
       gPad->Print((directory+"LO"+basename+".pdf").c_str());
       
       
       if ( hNLOxsec ) {
 	
 	reference->GetXaxis()->SetMoreLogLabels(true);
 	
 	hNLOxsec->DrawCopy();
 	
 	reference->SetLineColor(kBlack);
 	reference->SetLineStyle(2);
 	
 	reference->DrawCopy("samehist");
 	if ( show_errors ) reference->DrawCopy("samee1");
 	
 	for ( unsigned ic=0 ; ic<xsec_corr.size() ; ic++ ) { 
 	  xsec_corr[ic]->SetLineStyle(3);
 	  xsec_corr[ic]->DrawCopy("same");
 	}
 	
 	//    gPad->Print("NLO.pdf");
 	gPad->Print(("NLO"+basename+".pdf").c_str());
 	gPad->Print((directory+"NLO"+basename+".pdf").c_str());
 	
       }   
       
 
       if ( hNNLOxsec ) {
 	
 	reference->GetXaxis()->SetMoreLogLabels(true);
 	
 	hNNLOxsec->DrawCopy();
 	
 	reference->SetLineColor(kBlack);
 	reference->SetLineStyle(2);
 	
 	reference->DrawCopy("samehist");
 	if ( show_errors) reference->DrawCopy("samee1");
 	
 	for ( unsigned ic=0 ; ic<xsec_corr.size() ; ic++ ) { 
 	  xsec_corr[ic]->SetLineStyle(3);
 	  xsec_corr[ic]->DrawCopy("same");
 	}
 	
 	//    gPad->Print("NLO.pdf");
 	gPad->Print(("NNLO"+basename+".pdf").c_str());
 	gPad->Print((directory+"NNLO"+basename+".pdf").c_str());
 	
       }     
       
       if ( hNLOOxsec ) {
 	hNLOOxsec->DrawCopy();
 	// gPad->Print("NLO-only.pdf");
 	gPad->Print(("NLO-only"+basename+".pdf").c_str());
 	gPad->Print((directory+"NLO-only"+basename+".pdf").c_str());
       }
       
       if ( hNNLOOxsec ) {
 	hNNLOOxsec->DrawCopy();
 	// gPad->Print("NLO-only.pdf");
 	gPad->Print(("NNLO-only"+basename+".pdf").c_str());
 	gPad->Print((directory+"NNLO-only"+basename+".pdf").c_str());
       }
       
       
       gPad->SetLogy(false);
       gPad->SetLogx(false);
     }    
 
     xsDir.pop();
     
     
     /// write out invididual subprocesses
     
     if ( prepare_subproc && g.getGenpdf().find("photon")!=std::string::npos )   photon_subprocesses( g, basename, rscale, fscale ); 
     
     //    std::cout << "write out ..." << std::endl;
     
     if ( ! no_output ) { 
       fout->Write();
       fout->Close();
     }
 
     std::cout << "done" << std::endl;
 
 
     }    
     catch ( appl::grid::exception e ) { 
       std::cout << "caught exception " << e.what() << " for grid " << gridName << std::endl;
       continue;
     }
   }
   
 
   std::exit(0);
   return 0;
 }
 
 
 
 
 
 
 void processphoton( std::map<std::string,std::vector<int> >&  mLO,
 		    std::map<std::string,std::vector<int> >& mNLO )
 { 
   
   std::vector<int> v(2);
 
   v[0] = 2;  v[1] = 5;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "ug", v ) ); 
   v[0] = 4;  v[1] = 0;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "dg", v ) ); 
   
   v.resize(1);
   v[0] = 1;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd", v ) ); 
   v[0] = 3;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu", v ) ); 
 
   v.resize(4);
   v[0] = 10;  v[1] = 15;  v[2] = 18;  v[3] = 29;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "ug", v ) ); 
   v[0] = 3;  v[1] = 14;  v[2] = 17;  v[3] = 22;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dg", v ) ); 
 
   v.resize(3);
   v[0] = 0;  v[1] = 6;  v[2] = 23;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd", v ) ); 
   v[0] = 8;  v[1] = 13;  v[2] = 31;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu", v ) ); 
 
   v.resize(5);
   v[0] = 2;  v[1] = 4;  v[2] = 19;  v[3] = 21;  v[4] = 25;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd'", v ) ); 
   v[0] = 9;  v[1] = 12;  v[2] = 27;  v[3] = 28;  v[4] = 32;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu'", v ) ); 
 
   v.resize(8);
   v[0] = 1;  v[1] = 5;  v[2] = 7;  v[3] = 11;  v[4] = 20;  v[5] = 24;  v[6] = 26;  v[7] = 30;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "du", v ) ); 
 
   v.resize(1);
   v[0] = 16;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "gg", v ) ); 
  
 }
 
 
 
 
 void photon_subprocesses( appl::grid& g, const std::string& basename, double rscale, double fscale ) { 
   
   std::cout << "photon subprocesses ..." << std::endl;
 
   // if ( iorder==0 ) nsub = 4;
   // if ( iorder==1 ) nsub = 8;
   
   TH1D* hh = g.convolute( evolvepdf0, evolvepdf1, alphaspdf, 1, rscale, fscale);
   //  hh->Rebin(2);
 
 
   TH1D* hhf[8];
 
   std::map<std::string,std::vector<int> >  mLO;
   std::map<std::string,std::vector<int> > mNLO;
   
   processphoton( mLO, mNLO );
 
   std::string labs[8] = { "ug", "dg", "dd", "uu", "du", "dd'", "uu'", "gg" };
 
   for ( int ip=0 ; ip<8 ; ip++ ) { 
 
       std::cout << labs[ip] << " ";
       if ( labs[ip].size()<3 ) std::cout << " ";
     
       TH1D* hlo= 0;
       TH1D* hnlo= 0;
 
       std::map<std::string,std::vector<int> >::const_iterator  itr = mLO.find(labs[ip]);
       if ( itr!=mLO.end() ) { 
 	std::vector<int> vec = itr->second;
 	std::cout << "LO ";
 	for ( unsigned is=0 ; is<vec.size() ; is++ ) {
 	  std::cout << " " << vec[is];
 	  if ( hlo==0 ) hlo = g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, 0, rscale, fscale);
 	  else          hlo->Add( g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, 0, rscale, fscale) );
 	} 
 	std::cout << "\t";
       }
 
       itr = mNLO.find(labs[ip]);
       if ( itr!=mNLO.end() ) { 
 	std::vector<int> vec = itr->second;
 	std::cout << "NLO  ";
 	for ( unsigned is=0 ; is<vec.size() ; is++ ) {
 	  std::cout << " " << vec[is];
 	  if ( hnlo==0 ) hnlo = g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, -1, rscale, fscale);
 	  else           hnlo->Add( g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, -1, rscale, fscale) );
 	} 
 	std::cout << std::endl; 
       }
     
       if ( hlo!=0 && hnlo!=0 ) hnlo->Add( hlo );
 
       //      hnlo->Rebin(2);
       hnlo->SetName( labs[ip].c_str() );
       if ( !no_output ) hnlo->Write();
       
       hnlo->Divide( hh );
       hnlo->SetName( ("f"+labs[ip]).c_str() );
 
       //      smooth( hnlo );
 
       if ( !no_output ) hnlo->Write();
 
       hhf[ip] = hnlo;
   }
 
 
   if ( ! quiet ) { 
   //  hhf[0]->GetXaxis()->SetRangeUser(20, 450);
     hhf[0]->SetMaximum(1);
     hhf[0]->SetMinimum(-0.11);
     
     gStyle->SetOptStat(0);
     gStyle->SetPadTopMargin(0.05);
     gStyle->SetPadRightMargin(0.05);
     
     if ( hhf[0]->GetBinCenter(1)>100 ) hhf[0]->SetTitle(";E_{T}^{#gamma} [GeV];parton subprocess fraction");
     else                               hhf[0]->SetTitle(";|#eta^{#gamma}|;parton subprocess fraction");
     
     for ( int ip=0 ; ip<8 ; ip++ ) { 
       hhf[ip]->GetXaxis()->SetMoreLogLabels(true);
       hhf[ip]->SetLineColor( ip<4 ? ip+1 : ip+2 );
       hhf[ip]->SetLineWidth( ip<4 ? 2 : 4 );
       if ( ip==0 ) hhf[ip]->DrawCopy("lhist");
       else         hhf[ip]->DrawCopy("samelhist");
       
       // DrawLabel(0.8, 0.85-ip*0.05, labs[ip],  ip<4 ? ip+1 : ip+2, 0.038 );
     }
     
     if ( basename.find("pt")!=std::string::npos ) { 
       //  gPad->SetLogy(true);
       gPad->SetLogx(true);
     }    
     
     
     // gPad->Print("fractions.pdf"); 
     gPad->Print(("fractions"+basename+".pdf").c_str());
     // gPad->Print((directory+"fractions"+basename+".pdf").c_str());
     
     gPad->SetLogy(false);
     gPad->SetLogx(false);
   }
   
 }
 
 
 
 
 
 void process_subproc( std::map<std::string,std::vector<int> >&  mLO,
 		       std::map<std::string,std::vector<int> >& mNLO )
 { 
   
   std::vector<int> v(2);
 
   v[0] = 2;  v[1] = 5;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "ug", v ) ); 
   v[0] = 4;  v[1] = 0;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "dg", v ) ); 
   
   v.resize(1);
   v[0] = 1;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd", v ) ); 
   v[0] = 3;
   mLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu", v ) ); 
 
   v.resize(4);
   v[0] = 10;  v[1] = 15;  v[2] = 18;  v[3] = 29;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "ug", v ) ); 
   v[0] = 3;  v[1] = 14;  v[2] = 17;  v[3] = 22;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dg", v ) ); 
 
   v.resize(3);
   v[0] = 0;  v[1] = 6;  v[2] = 23;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd", v ) ); 
   v[0] = 8;  v[1] = 13;  v[2] = 31;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu", v ) ); 
 
   v.resize(5);
   v[0] = 2;  v[1] = 4;  v[2] = 19;  v[3] = 21;  v[4] = 25;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "dd'", v ) ); 
   v[0] = 9;  v[1] = 12;  v[2] = 27;  v[3] = 28;  v[4] = 32;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "uu'", v ) ); 
 
   v.resize(8);
   v[0] = 1;  v[1] = 5;  v[2] = 7;  v[3] = 11;  v[4] = 20;  v[5] = 24;  v[6] = 26;  v[7] = 30;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "du", v ) ); 
 
   v.resize(1);
   v[0] = 16;
   mNLO.insert( std::map<std::string,std::vector<int> >::value_type( "gg", v ) ); 
  
 }
 
 
 
 
 void subprocesses( appl::grid& g, const std::string& basename, double rscale, double fscale ) { 
   
   std::cout << "photon subprocesses ..." << std::endl;
 
   // if ( iorder==0 ) nsub = 4;
   // if ( iorder==1 ) nsub = 8;
   
   TH1D* hh = g.convolute( evolvepdf0, evolvepdf1,  alphaspdf, 1, rscale, fscale);
   //  hh->Rebin(2);
 
 
   TH1D* hhf[8];
 
   std::map<std::string,std::vector<int> >  mLO;
   std::map<std::string,std::vector<int> > mNLO;
   
   processphoton( mLO, mNLO );
 
   std::string labs[8] = { "ug", "dg", "dd", "uu", "du", "dd'", "uu'", "gg" };
 
   for ( int ip=0 ; ip<8 ; ip++ ) { 
 
       std::cout << labs[ip] << " ";
       if ( labs[ip].size()<3 ) std::cout << " ";
     
       TH1D* hlo= 0;
       TH1D* hnlo= 0;
 
       std::map<std::string,std::vector<int> >::const_iterator  itr = mLO.find(labs[ip]);
       if ( itr!=mLO.end() ) { 
 	std::vector<int> vec = itr->second;
 	std::cout << "LO ";
 	for ( unsigned is=0 ; is<vec.size() ; is++ ) {
 	  std::cout << " " << vec[is];
 	  if ( hlo==0 ) hlo = g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, 0, rscale, fscale);
 	  else          hlo->Add( g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, 0, rscale, fscale) );
 	} 
 	std::cout << "\t";
       }
 
       itr = mNLO.find(labs[ip]);
       if ( itr!=mNLO.end() ) { 
 	std::vector<int> vec = itr->second;
 	std::cout << "NLO  ";
 	for ( unsigned is=0 ; is<vec.size() ; is++ ) {
 	  std::cout << " " << vec[is];
 	  if ( hnlo==0 ) hnlo = g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, -1, rscale, fscale);
 	  else           hnlo->Add( g.convolute_subproc( vec[is], evolvepdf0, alphaspdf, -1, rscale, fscale) );
 	} 
 	std::cout << std::endl; 
       }
     
       if ( hlo!=0 && hnlo!=0 ) hnlo->Add( hlo );
 
       //      hnlo->Rebin(2);
       hnlo->SetName( labs[ip].c_str() );
       if ( !no_output ) hnlo->Write();
       
       hnlo->Divide( hh );
       hnlo->SetName( ("f"+labs[ip]).c_str() );
 
       //      smooth( hnlo );
 
       if ( !no_output )  hnlo->Write();
 
       hhf[ip] = hnlo;
   }
 
   if ( !quiet ) { 
     //  hhf[0]->GetXaxis()->SetRangeUser(20, 450);
     hhf[0]->SetMaximum(1);
     hhf[0]->SetMinimum(-0.11);
     
     gStyle->SetOptStat(0);
     gStyle->SetPadTopMargin(0.05);
     gStyle->SetPadRightMargin(0.05);
     
     if ( hhf[0]->GetBinCenter(1)>100 ) hhf[0]->SetTitle(";E_{T}^{#gamma} [GeV];parton subprocess fraction");
     else                               hhf[0]->SetTitle(";|#eta^{#gamma}|;parton subprocess fraction");
     
     for ( int ip=0 ; ip<8 ; ip++ ) { 
       hhf[ip]->GetXaxis()->SetMoreLogLabels(true);
       hhf[ip]->SetLineColor( ip<4 ? ip+1 : ip+2 );
       hhf[ip]->SetLineWidth( ip<4 ? 2 : 4 );
       if ( ip==0 ) hhf[ip]->DrawCopy("lhist");
       else         hhf[ip]->DrawCopy("samelhist");
       
       // DrawLabel(0.8, 0.85-ip*0.05, labs[ip],  ip<4 ? ip+1 : ip+2, 0.038 );
     }
     
     if ( basename.find("pt")!=std::string::npos ) { 
       //  gPad->SetLogy(true);
       gPad->SetLogx(true);
     }    
     
     
 
     // gPad->Print("fractions.pdf"); 
     gPad->Print(("fractions"+basename+".pdf").c_str());
     // gPad->Print((directory+"fractions"+basename+".pdf").c_str());
     
     gPad->SetLogy(false);
     gPad->SetLogx(false);
   }
 
 }