Index: applgrid/trunk/src/appl_igrid.cxx
===================================================================
--- applgrid/trunk/src/appl_igrid.cxx	(revision 1752)
+++ applgrid/trunk/src/appl_igrid.cxx	(revision 1753)
@@ -1,1950 +1,1949 @@
-
 //  appl_igrid.cxx        
 
 //   grid class - all the functions needed to create and 
 //   fill the grid from an NLO calculation program. 
 //  
 //  Copyright (C) 2007 Mark Sutton (sutt@hep.ucl.ac.uk)    
 
 //   $Id: appl_igrid.cxx, v1.00 2007/10/16 17:01:39 sutt
 
 #include <stdlib.h>
 #include <iostream>
 #include <iomanip>
 #include <cmath>
 
 
 #include "appl_igrid.h"
 #include "appl_grid/appl_grid.h"
 #include "appl_grid/appl_timer.h"
 
 #include "hoppet_init.h"
 #include "threadManager.h"
 
 
 #include "TFile.h"
 #include "TObject.h"
 #include "TObjString.h"
 #include "TH3D.h"
 #include "TVectorT.h"
 #include "Splitting.h"
 
 #include "TFileString.h"
 
 
 // splitting function code
 
 
 // pdf reweighting
 // bool   igrid::m_reweight   = false;
 // bool   igrid::m_symmetrise = false;
 
 // variable tranformation parameters
 double appl::igrid::transvar = 5;
 bool   appl::igrid::threads_disabled = true;
 
 static int ithread = 0;
 
 std::string label( int i ) { 
   char lab[64];
   std::sprintf( lab, "thread-%d", i );
   return lab;
 }
 
 
 
 appl::igrid::igrid() : 
   threadManager( label(ithread++) ),
   mfy(0),   mfx(0),
   m_parent(0),
   m_Ny1(0),   m_y1min(0),   m_y1max(0),   m_deltay1(0),
   m_Ny2(0),   m_y2min(0),   m_y2max(0),   m_deltay2(0),
   m_yorder(0),   
   m_Ntau(0), m_taumin(0), m_taumax(0), m_deltatau(0),   m_tauorder(0), 
   m_Nproc(0),
   m_transform(""), 
   //  m_transvarlocal(m_transvar),
   m_transvar(transvar),
   m_reweight(false),
   m_symmetrise(false),
   m_optimised(false),
   m_weight(0),
   m_fg1(0),     m_fg2(0),
   m_fsplit1(0), m_fsplit2(0),
   m_fsplit12(0), m_fsplit22(0),
   m_alphas(0),
   m_weights(0),
   m_taufilledmin(-1),  m_taufilledmax(-1),
   m_y1filledmin(-1),   m_y1filledmax(-1),
   m_y2filledmin(-1),   m_y2filledmax(-1)
 { 
 
   //  std::cout << "igrid() (default) Ntau=" << m_Ntau << "\t" << fQ2(m_taumin) << " - " << fQ2(m_taumax) << std::endl;
 
 } 
 
 
 
 // standard constructor
 appl::igrid::igrid(int NQ2, double Q2min, double Q2max, int Q2order, 
 		   int Nx,  double xmin,  double xmax,  int xorder, 
 		   std::string transform, int Nproc, bool disflag ):
   threadManager( label(ithread++) ),
   mfy(0),   mfx(0),
   m_parent(0),
   m_Ny1(Nx),   m_Ny2( disflag ? 1 : Nx ),  m_yorder(xorder), 
   m_Ntau(NQ2), m_tauorder(Q2order), 
   m_Nproc(Nproc), 
   m_transform(transform), 
   //  m_transvarlocal(m_transvar),
   m_transvar(transvar),
   m_reweight(false),
   m_symmetrise(false), 
   m_optimised(false),
   m_weight(0),
   m_fg1(0),     m_fg2(0),  
   m_fsplit1(0), m_fsplit2(0),
   m_fsplit12(0), m_fsplit22(0),
   m_alphas(0),
   m_DISgrid(disflag),   
   m_weights(0), 
   m_taufilledmin(-1),  m_taufilledmax(-1),
   m_y1filledmin(-1),   m_y1filledmax(-1),
   m_y2filledmin(-1),   m_y2filledmax(-1)
 {
   //  std::cout << "igrid::igrid() transform=" << m_transform << std::endl;
   init_fmap();
   if ( m_fmap.find(m_transform)==m_fmap.end() ) throw exception("igrid::igrid() transform " + m_transform + " not found\n");
   mfx = m_fmap.find(m_transform)->second.mfx;
   mfy = m_fmap.find(m_transform)->second.mfy;
   
 
 
   double ymin1 = fy(xmax);
   double ymax1 = fy(xmin);
 
   m_y1min  = ymin1 ;
   m_y1max  = ymax1 ;
   m_y2min  = ymin1 ;
   m_y2max  = ymax1 ;
 
   if ( m_DISgrid ) m_y2min = m_y2max = 1;
 
   if ( m_Ny1>1 ) m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
   else           m_deltay1 = 0;
   
   if ( m_Ny2>1 ) m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
   else           m_deltay2 = 0;
   
   double taumin=ftau(Q2min);
   double taumax=ftau(Q2max);
   m_taumin   = taumin;
   m_taumax   = taumax;
   if ( m_Ntau>1 ) m_deltatau = (taumax-taumin)/(m_Ntau-1);
   else            m_deltatau = 0;
 
   if ( m_Ny1-1<m_yorder ) { 
     std::cerr << "igrid() not enough nodes for this interpolation order Ny1=" << m_Ny1 
 	      << "\tyorder=" << m_yorder << std::endl;
  
     while ( m_Ny1-1<m_yorder ) m_yorder--;
   } 
 
   if ( !m_DISgrid ) { 
     if ( m_Ny2-1<m_yorder ) { 
       std::cerr << "igrid() not enough nodes for this interpolation order Ny2=" << m_Ny2 
 		<< "\tyorder=" << m_yorder << std::endl;
       
       while ( m_Ny2-1<m_yorder ) m_yorder--;
     } 
   }
   
   if ( m_Ntau-1<m_tauorder ) { 
     std::cerr << "igrid() not enough nodes for this interpolation order Ntau=" << m_Ntau 
 	      << "\ttauorder=" << m_tauorder << std::endl;
     
     while ( m_Ntau-1<m_tauorder ) m_tauorder--;
   } 
   
   m_weight = new SparseMatrix3d*[m_Nproc];
   construct();
 
   if ( !threads_disabled ) this->start_thread();
 }
 
 
 
 
 
 // copy constructor
 appl::igrid::igrid(const appl::igrid& g) :
   threadManager( label(ithread++) ), 
   mfy(0),  mfx(0),  
   m_parent(0),
   m_Ny1(g.m_Ny1),     
   m_y1min(g.m_y1min),     m_y1max(g.m_y1max),     m_deltay1(g.m_deltay1),   
   m_Ny2(g.m_Ny2),     
   m_y2min(g.m_y2min),     m_y2max(g.m_y2max),     m_deltay2(g.m_deltay2),   
   m_yorder(g.m_yorder),   
   m_Ntau(g.m_Ntau), 
   m_taumin(g.m_taumin), m_taumax(g.m_taumax), m_deltatau(g.m_deltatau), m_tauorder(g.m_tauorder), 
   m_Nproc(g.m_Nproc),
   m_transform(g.m_transform), 
   //  m_transvarlocal(g.m_transvarlocal),
   m_transvar(g.m_transvar),
   m_reweight(g.m_reweight),
   m_symmetrise(g.m_symmetrise),
   m_optimised(g.m_optimised),
   m_weight(0),
   m_fg1(0),     m_fg2(0),
   m_fsplit1(0), m_fsplit2(0),
   m_fsplit12(0), m_fsplit22(0),
   m_alphas(0),
   m_weights(0),
   m_taufilledmin(g.m_taufilledmin),  m_taufilledmax(g.m_taufilledmax),
   m_y1filledmin(g.m_y1filledmin),    m_y1filledmax(g.m_y1filledmax),
   m_y2filledmin(g.m_y2filledmin),    m_y2filledmax(g.m_y2filledmax) 
 {
   init_fmap();
   if ( m_fmap.find(m_transform)==m_fmap.end() ) throw exception("igrid::igrid() transform " + m_transform + " not found\n");
   mfx = m_fmap.find(m_transform)->second.mfx;
   mfy = m_fmap.find(m_transform)->second.mfy;
 
 
   m_weight = new SparseMatrix3d*[m_Nproc];
 
 
   for( int ip=0 ; ip<m_Nproc ; ip++ )   m_weight[ip] = new SparseMatrix3d(*g.m_weight[ip]);
   //  construct();
 
   if ( !threads_disabled ) this->start_thread();
 }
 
 
 
 void _setlimits( int& _min, int& _max, const int _mint, const int _maxt ) {  
   if ( _mint<=_maxt ) { 
     if ( _min==-1 || _min>_mint ) _min = _mint;
     if ( _max==-1 || _max<_maxt ) _max = _maxt;
   }
 }
 
 
 
 
 // read from a file 
 appl::igrid::igrid(TFile& f, const std::string& s) :
   threadManager( label(ithread++) ), 
   mfy(0),  mfx(0),  
   m_parent(0),
   m_Ny1(0),   m_y1min(0),   m_y1max(0),   m_deltay1(0),   
   m_Ny2(0),   m_y2min(0),   m_y2max(0),   m_deltay2(0),   
   m_yorder(0),   
   m_Ntau(0), m_taumin(0), m_taumax(0), m_deltatau(0), m_tauorder(0), 
   m_Nproc(0),
   m_transform(""), 
   //  m_transvarlocal(m_transvar),
   m_transvar(transvar),
   m_reweight(false),
   m_symmetrise(false),
   m_optimised(false),
   m_weight(0), 
   m_fg1(0),     m_fg2(0),
   m_fsplit1(0), m_fsplit2(0),    
   m_fsplit12(0), m_fsplit22(0),    
   m_alphas(0), 
   m_weights(0),
   m_taufilledmin(-1),  m_taufilledmax(-1),
   m_y1filledmin(-1),   m_y1filledmax(-1),
   m_y2filledmin(-1),   m_y2filledmax(-1)
 { 
   //  std::cout << "igrid::igrid()" << std::endl;
   
   // get the name of the transform pair
   TFileString _tag = *(TFileString*)f.Get((s+"/Transform").c_str());
   m_transform = _tag[0];
 
   init_fmap();
   if ( m_fmap.find(m_transform)==m_fmap.end() ) throw exception("igrid::igrid() transform " + m_transform + " not found\n");
   mfx = m_fmap.find(m_transform)->second.mfx;
   mfy = m_fmap.find(m_transform)->second.mfy;
 
 
   // delete _transform;
   
   // retrieve setup parameters 
   TVectorT<double>* setup=(TVectorT<double>*)f.Get((s+"/Parameters").c_str());
   //  f.GetObject((s+"/Parameters").c_str(), setup);
 
   // NB: round integer variables to nearest integer 
   //     in case (unlikely) truncation error during 
   //     conversion to double when they were stored
   m_Ny1      = int((*setup)(0)+0.5);  
   m_y1min    = (*setup)(1);
   m_y1max    = (*setup)(2);
 
   m_Ny2      = int((*setup)(3)+0.5);  
   m_y2min    = (*setup)(4);
   m_y2max    = (*setup)(5);
 
   m_yorder   = int((*setup)(6)+0.5);
 
   m_Ntau     = int((*setup)(7)+0.5);
   m_taumin   = (*setup)(8);
   m_taumax   = (*setup)(9);
   m_tauorder = int((*setup)(10)+0.5);
 
   m_transvar = (*setup)(11);
 
   m_Nproc    = int((*setup)(12)+0.5);
 
   m_reweight   = ((*setup)(13)!=0 ? true : false );
   m_symmetrise = ((*setup)(14)!=0 ? true : false );
   m_optimised  = ((*setup)(15)!=0 ? true : false );
   m_DISgrid    =  ( setup->GetNoElements()>16 && (*setup)[16]!=0 ) ;
 
 
   delete setup;
 
   //  std::cout << "igrid::igrid() read setup" << std::endl;
 
   // create grids
   m_deltay1   = (m_y1max-m_y1min)/(m_Ny1-1);
   m_deltay2   = (m_y2max-m_y2min)/(m_Ny2-1);
 
   m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
   
   //  int rawsize=0;
   //  int trimsize=0;
 
   m_weight = new SparseMatrix3d*[m_Nproc];
 
   for( int ip=0 ; ip<m_Nproc ; ip++ ) {
     char _name[128];  sprintf(_name,"/weight[%i]", ip );
     // get storage histogram
     TH3D* htmp = (TH3D*)f.Get((s+_name).c_str()); 
  
     //    std::cout << "igrid::igrid() read " << _name << std::endl;
 
     // create grid
     m_weight[ip]=new SparseMatrix3d(htmp);
 
     // save some space
     m_weight[ip]->trim();
 
     // delete storage histogram
     delete htmp;
 
     // DON'T trim on reading unless the user wants it!!
     // he can trim himself if need be!! 
     // rawsize += m_weight[ip]->size();
     // m_weight[ip]->trim(); // trim the grid and do some book keeping
     // trimsize += m_weight[ip]->size();
   }
 
   /// now calculate the actual limits of this grid
   
   //  static double _ctime = 0;
 
   setlimits();
 
   if ( !threads_disabled ) this->start_thread();
 
 }
 
 
 
 
 void appl::igrid::setlimits() { 
   if ( m_weight ) { 
     //   struct timeval _tstart = appl_timer_start(); 
     for ( int i=0 ; i<m_Nproc ; i++ ) {
       const SparseMatrix3d* _weight = m_weight[i];
       if ( _weight ) { 
 	if ( _weight->empty() ) continue;
 	_setlimits( m_taufilledmin, m_taufilledmax,  _weight->xmin(),  _weight->xmax() ); 
 	_setlimits( m_y1filledmin,  m_y1filledmax,   _weight->ymin(),  _weight->ymax() ); 
 	_setlimits( m_y2filledmin,  m_y2filledmax,   _weight->zmin(),  _weight->zmax() ); 
       }
     }
   }
 }    
 
 
 
 // constructor common internals 
 void appl::igrid::construct() 
 {
   // Initialize histograms representing the weight grid
   for( int ip=0 ; ip<m_Nproc ; ip++ ) {
     m_weight[ip] = new SparseMatrix3d(m_Ntau, m_taumin,   m_taumax,    
     				      m_Ny1,  m_y1min,    m_y1max, 
     				      m_Ny2,  m_y2min,    m_y2max ); 
     
   }  
 }
 
 
 // destructor
 appl::igrid::~igrid() {
   deletepdftable();
   deleteweights();
 }
 
 
 
 
 double*** delete_array( double*** a, unsigned Nx, unsigned Ny ) { 
   if ( a ) { 
     for ( unsigned i=0 ; i<Nx ; i++ ) {
       for ( unsigned j=0 ; j<Ny ; j++ )  delete[] a[i][j];
       delete[] a[i];
     }
     delete[] a;
   }
   return 0;
 }
 
 
 // clean up internal pdf table (could be moved to local variable)
 void appl::igrid::deletepdftable() { 
 
   //  std::cout << "deleting pdf tables" << std::endl;
 
   //  std::cout << "deleting pdf and splitting tables" << std::endl;
   m_fg1       = delete_array( m_fg1,       m_Ntau, Ny1() );
   m_fsplit1   = delete_array( m_fsplit1,   m_Ntau, Ny1() );
   m_fsplit12  = delete_array( m_fsplit12,  m_Ntau, Ny1() );
   //  m_fsplit111 = delete_array( m_fsplit111, m_Ntau, Ny1() );
 
   // if grid is symmetric, then don't need to deallocate 
   // x2 pdf values and splitting functions 
   if ( isSymmetric() ) { 
     m_fg2       = 0;
     m_fsplit2   = 0;
     m_fsplit22  = 0;
     //   m_fsplit211 = 0;
   }
   else { 
     //  std::cout << "deleting pdf and splitting tables" << std::endl;
     m_fg2       = delete_array( m_fg2,       m_Ntau, Ny1() );
     m_fsplit2   = delete_array( m_fsplit2,   m_Ntau, Ny1() );
     m_fsplit22  = delete_array( m_fsplit22,  m_Ntau, Ny1() );
     //    m_fsplit211 = delete_array( m_fsplit211, m_Ntau, Ny1() );
   }
 
   //  std::cout << "delete alphas table" << std::endl;
 
   if ( m_alphas ) { 
     delete[] m_alphas;
     m_alphas=0;
   }
 }
 
 
 
 
 
 // delete internal grids
 void appl::igrid::deleteweights() { 
   if ( m_weight ) { 
     for ( int ip=0 ; ip<m_Nproc ; ip++ ) if (m_weight[ip]) delete m_weight[ip];
     delete[] m_weight;
     // m_weight=0;
   }
 }
 
 
 
 
 // write to file
 void appl::igrid::write(const std::string& _name) { 
   Directory d(_name);
   d.push();
 
   // using a TH1D to store the transform pair tag since I don't know how 
   // to write a TString to a root file
   // TH1D* _transform = new TH1D("Transform", m_transform.c_str(), 1, 0, 1);
   //  _transform->Write();
 
   // write the name of the transform pair
   TFileString("Transform",m_transform).Write();
 
 
   TVectorT<double>* setup=new TVectorT<double>(20); // a few spare
 
   (*setup)(0)  = m_Ny1;
   (*setup)(1)  = m_y1min;
   (*setup)(2)  = m_y1max;
 
   (*setup)(3)  = m_Ny2;
   (*setup)(4)  = m_y2min;
   (*setup)(5)  = m_y2max;
 
   (*setup)(6)  = m_yorder;
 
   (*setup)(7)  = m_Ntau;
   (*setup)(8)  = m_taumin;
   (*setup)(9)  = m_taumax;
   (*setup)(10) = m_tauorder;
 
   (*setup)(11) = m_transvar;
 
   (*setup)(12) = m_Nproc;
  
   (*setup)(13) = ( m_reweight   ? 1 : 0 );
   (*setup)(14) = ( m_symmetrise ? 1 : 0 );
   (*setup)(15) = ( m_optimised  ? 1 : 0 );
   (*setup)(16) = ( m_DISgrid    ? 1 : 0 );
 
   setup->Write("Parameters");
 
   int igridsize     = 0;
   int igridtrimsize = 0;
 
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) { 
 
     char hname[128];
     //    sprintf(hname,"%s[%d]", _name.c_str(), ip);
     sprintf(hname,"weight[%d]", ip);
 
     //    int oldsize = m_weight[ip]->size();
 
     igridsize += m_weight[ip]->size();
  
     // trim it so that it's quicker to copy into the TH3D
     m_weight[ip]->trim();
 
     igridtrimsize += m_weight[ip]->size();
 
     TH3D* h=m_weight[ip]->getTH3D(hname);
     h->SetDirectory(0);
     h->Write();
     delete h; // is this dengerous??? will root try to delete it later? I guess not if we SetDirectory(0)
   }
 
 #if 0
   //    std::cout << _name << " trimmed" << std::endl;
   std::cout << _name << "\tsize=" << igridsize << "\t-> " << igridtrimsize;
   if ( igridsize ) std::cout << "\t( " << igridtrimsize*100/igridsize << "% )";
   std::cout << std::endl;
 #endif
 
   d.pop();
 }
 
 
 
 
 
 void appl::igrid::fill(const double x1, const double x2, const double Q2, const double* weight) 
 {  
 
   // find preferred vertex for low end of interpolation range
   int k1=fk1(x1);
   int k2=fk2(x2);
   int k3=fkappa(Q2);
   
   double u_y1  = ( fy(x1)-gety1(k1) )/deltay1();
   double u_y2  = ( fy(x2)-gety2(k2) )/deltay2();
   double u_tau = ( ftau(Q2)-gettau(k3))/deltatau();
   
   
   double fillweight, fI_factor;
   
   // cache the interpolation coefficients so only 
   // have to calculate each one once 
   // NB: a (very naughty) hardcoded maximum interpolation order 
   //     of 16 for code simplicity.
   double _fI1[16];
   double _fI2[16];
   double _fI3[16];
   
   for( int i1=0 ; i1<=m_yorder   ; i1++ )  _fI1[i1] = fI(i1, m_yorder, u_y1);
   for( int i2=0 ; i2<=m_yorder   ; i2++ )  _fI2[i2] = fI(i2, m_yorder, u_y2);
   for( int i3=0 ; i3<=m_tauorder ; i3++ )  _fI3[i3] = fI(i3, m_tauorder, u_tau);
   
   double invwfun = 1/(weightfun(x1)*weightfun(x2));
 	 
   for( int i3=0 ; i3<=m_tauorder ; i3++ ) { // "interpolation" loop in Q2
 
     for( int i1=0 ; i1<=m_yorder ; i1++ )   { // "interpolation" loop in x1      
 
       for( int i2=0 ; i2<=m_yorder ; i2++ ) { // "interpolation" loop in x2
 	
 	fI_factor = _fI1[i1] * _fI2[i2] * _fI3[i3];
 
 	if ( m_reweight ) fI_factor *= invwfun;
 	
 	for( int ip=0 ; ip<m_Nproc ; ip++ ) {
 	  
 	  fillweight = weight[ip] * fI_factor;
 
 	  // this method only works for grids where the x1, x2 axes are identical 
 	  // ranges and binning 
 	  //  if ( m_symmetrise ) { 
 	  //	// symmetrise the grid
 	  //    if ( abs(k1+i1)<abs(k2+i2) )  (*m_weight[ip])(k3+i3, k2+i2, k1+i1) += fillweight;   
 	  //    else                          (*m_weight[ip])(k3+i3, k1+i1, k2+i2) += fillweight;  
 	  //  } 
 
 	  //	  std::cout << " (" << (*m_weight[ip])(k3+i3, k1+i1, k2+i2) << " -> ";
 
 	  (*m_weight[ip])(k3+i3, k1+i1, k2+i2) += fillweight;	  
 	 	  
 	  //	  std::cout << (*m_weight[ip])(k3+i3, k1+i1, k2+i2) << ")";
 
 	  //	  m_weight[ip]->print();	  
 	  //	  m_weight[ip]->fill_fast(k3+i3, k1+i1, k2+i2) += fillweight/wfun;	  
 	}
 
 	//	std::cout << std::endl;
       }
     }
   }
 }
 
 
 void appl::igrid::fill_phasespace(const double x1, const double x2, const double Q2, const double* weight) { 
 
   int k1=fk1(x1);
   int k2=fk2(x2);
   int k3=fkappa(Q2);
 
   //  std::cout << "appl::igrid::fill_phasespace() k: " << k1 << " " << k2 << " " << k3 << std::endl;
   //  std::cout << "\tweights: ";
 
   //  for ( int ip=0 ; ip<m_Nproc ; ip++ ) std::cout << " " << weight[ip];
   //  std::cout << std::endl;
 
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) (*m_weight[ip])(k3, k1, k2) += weight[ip];
 
   //  for ( int ip=0 ; ip<m_Nproc ; ip++ ) std::cout << "\t" << ip << "\tsparse size " << m_weight[ip]->size() << std::endl; 
 
   //  for ( int ip=0 ; ip<m_Nproc ; ip++ ) std::cout << "\t" << ip << "\tweight " << (*m_weight[ip])(k3, k1, k2) << std::endl; 
 
 } 
 
 
 
 
 
 void appl::igrid::fill_index(const int ix1, const int ix2, const int iQ2, const double* weight) { 
 
   //  int k1=ix1;
   //  int k2=ix2;
   //  int k3=iQ2;
 
   //  for ( int ip=0 ; ip<m_Nproc ; ip++ ) (*m_weight[ip])(i3, k1, k2) += weight[ip];
 
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) (*m_weight[ip])(iQ2, ix1, ix2) += weight[ip];
 
 } 
 
 
 
 
 
 
 void appl::igrid::setuppdf(double (*alphas)(const double&),
 			   NodeCache* pdf0,
 			   NodeCache* pdf1,
 			   int _nloop,
 			   double rscale_factor,
 			   double fscale_factor,
 			   double beam_scale ) 
 {
 
 
   /// if this igrid is empty, don't do anything 
 
   //  std::cout << m_taufilledmin << " " << m_taufilledmax << " " 
   //	        << m_y1filledmin  << " " << m_y1filledmax  << " " 
   //	        << m_y2filledmin  << " " << m_y2filledmax  << std::endl;
 
   if ( ( m_taufilledmin==-1 || m_taufilledmax==-1 ) ||
        (  m_y1filledmin==-1 ||  m_y1filledmax==-1 ) ||
        (  m_y2filledmin==-1 ||  m_y2filledmax==-1 ) ) return; 
   
 
   int nloop = std::fabs(_nloop);
 
   if ( pdf1==0 ) pdf1 = pdf0;
 
   bool initialise_hoppet = false;
 
 #ifdef HAVE_HOPPET
   // check if we need to use the splitting function, and if so see if we 
   // need to initialise it again, and do so if required
   if ( fscale_factor!=1 ) {
     if ( pdf1 && pdf1!=pdf0 ) inititialise_hoppet = false;
   }
 #endif
 
 
   void (*splitting)(const double& , const double&, double*, int ) = Splitting;
 
 		     
   const int n_tau = Ntau();
   const int n_y1  = Ny1();
   const int n_y2  = Ny2();
 
   static const double twopi = 2*M_PI;
 
   static const int nc = 3;
   //TC   const int nf = 6;
   static const int nf = 5;
   static double beta0=(11.*nc - 2.*nf)/(6*twopi);
   
   /// not needed here ?? 
   //  static double beta1=(34.*nc*nc - 3.*(nc-1./nc)*nf - 10.*nc*nf)/(6.*twopi);
 
   double logf = std::log( fscale_factor*fscale_factor );
   double logr = std::log( rscale_factor*rscale_factor );
 
   // pdf table
   m_fg1 = new double**[Ntau()];
   if ( isSymmetric() ) m_fg2 = m_fg1;
   else                 m_fg2 = new double**[Ntau()];
 
   const double invtwopi = 0.5/(M_PI);
 
   // splitting function table for nlo 
   // factorisation scale dependence
   if ( nloop>=1 && fscale_factor!=1 ) { 
     m_fsplit1 = new double**[Ntau()];
     if ( isSymmetric() ) m_fsplit2 = m_fsplit1;
     else                 m_fsplit2 = new double**[Ntau()];
     if ( nloop==2 ) {
       m_fsplit12 = new double**[Ntau()];
       if ( isSymmetric() ) m_fsplit22 = m_fsplit12;
       else                 m_fsplit22 = new double**[Ntau()];
     }
   }
 
 
   // alphas table
   m_alphas = new double[Ntau()];
 
 
   // allocate tables
   for ( int i=0 ; i<n_tau ; i++ ) {
     
     // pdf table
     m_fg1[i] = new double*[n_y1];
     for ( int j=0 ; j<n_y1 ; j++ ) m_fg1[i][j] = new double[13];
     if ( !isSymmetric() && !isDISgrid() ) { 
       m_fg2[i] = new double*[n_y2];
       for ( int j=0 ; j<n_y2 ; j++ ) m_fg2[i][j] = new double[13];   
     }
     
     // splitting function table
     if ( nloop>=1 && fscale_factor!=1 ) { 
       m_fsplit1[i] = new double*[n_y1];
       for ( int j=0 ; j<n_y1 ; j++ ) m_fsplit1[i][j] = new double[13];
       if ( !isSymmetric() && !isDISgrid() ) { 
 	m_fsplit2[i] = new double*[n_y2];
 	for ( int j=0 ; j<n_y2 ; j++ ) m_fsplit2[i][j] = new double[13];   
       }
       if (nloop == 2) {
 	m_fsplit12[i] = new double*[n_y1];
 	for ( int j=0 ; j<n_y1 ; j++ ) m_fsplit12[i][j] = new double[13];
 	if ( !isSymmetric() && !isDISgrid() ) { 
 	  m_fsplit22[i] = new double*[n_y2];
 	  for ( int j=0 ; j<n_y2 ; j++ ) m_fsplit22[i][j] = new double[13];   
 	}
       }
     }
   }
 
   
   bool scale_beams = false;
   if ( beam_scale!=1 ) scale_beams = true;
 
   if ( initialise_hoppet ) hoppet_init::assign( pdf0->pdf() );
   
   // set up pdf grid, splitting function grid 
   // and alpha_s grid
   //  for ( int itau=m_Ntau ; itau-- ; ) {
   for ( int itau=0 ; itau<m_Ntau ; itau++  ) {
     
     if ( itau<m_taufilledmin || itau>m_taufilledmax ) continue;
 
     double tau = gettau(itau);
     double Q2  = fQ2(tau);
     double Q   = std::sqrt(Q2); 
     
     // alpha_s table 
     m_alphas[itau] = alphas(rscale_factor*Q)*invtwopi;
 
     //    std::cout << itau << "\ttau " << tau 
     //    	    << "\tQ2 " << Q2 << "\tQ " << Q 
     //   	    << "\talphas " << m_alphas[itau] << std::endl; 
 
 #if 0
     int iymin1 = m_weight[0]->ymin();
     int iymax1 = m_weight[0]->ymax();
     
     if ( isSymmetric() ) { 
       int iymin2 = iymin1;
       int iymax2 = iymax1;
     }
     else {
       int iymin2 = m_weight[0]->zmin();
       int iymax2 = m_weight[0]->zmax();
     }  
 #endif
 
     // grid not filled for iy<iymin || iy>iymax so no need to 
     // consider outside this range
     
     // y1 tables
     for ( int iy=n_y1 ; iy-- ;  ) { 
       
       if ( iy<m_y1filledmin || iy>m_y1filledmax ) continue;
 
       double y = gety1(iy);
       double x = fx(y);
       double fun = 1;
       if ( m_reweight ) fun = weightfun(x);
       
       if ( scale_beams ) { 
 	x *= beam_scale;
 	if ( x>=1 ) { 
 	  for ( int ip=0 ; ip<13 ; ip++ ) m_fg1[itau][iy][ip]=0; 
 	  if ( nloop>=1 && fscale_factor!=1 ) { 
 	    for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit1[itau][iy][ip] = 0;
 	    if (nloop==2)
 	      for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit12[itau][iy][ip] = 0;
 	  }
 	  continue; 
 	}    
       }
 
       pdf0->evaluate(x, fscale_factor*Q, m_fg1[itau][iy]);
       
       double invx = 1/x;
       for ( int ip=0 ; ip<13 ; ip++ ) m_fg1[itau][iy][ip] *= invx;
       if ( m_reweight ) for ( int ip=0 ; ip<13 ; ip++ ) m_fg1[itau][iy][ip] *= fun;
       
       // splitting function table
       if ( nloop>=1 && fscale_factor!=1 ) { 
 
 	double PLO[13]; for ( int ip=13 ; ip-- ; ) PLO[ip] = 0;  
 	splitting( x, fscale_factor*Q, PLO, 1 ); /// PLO x pdf	
 	for ( int ip=13 ; ip-- ; ) m_fsplit1[itau][iy][ip] = PLO[ip]*invx;
 	if ( m_reweight ) for ( int ip=13 ; ip-- ; ) m_fsplit1[itau][iy][ip] *= fun;
 
 	if (nloop>1) {
 	  //	splitting( x, fscale_factor*Q, m_fsplit12[itau][iy], 1 ); placeholder for hoppet call
 
 	  double P2LO[13]; for ( int ip=13 ; ip-- ; ) P2LO[ip] = 0;  
 	  double PNLO[13]; for ( int ip=13 ; ip-- ; ) PNLO[ip] = 0;  
 
 	  splitting( x, fscale_factor*Q, P2LO, 11 ); /// PLO x PLO x pdf
 	  splitting( x, fscale_factor*Q, PNLO, 2 );  /// PNLO x pdf
 
 	  for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit12[itau][iy][ip] = 0.5*logf*(beta0*(logf - 2*logr)*PLO[ip] + logf*P2LO[ip] - 2*PNLO[ip] )*invx ;
 
 	  if ( m_reweight ) for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit12[itau][iy][ip] *= fun;
 
 	}
       }
     }
   }
 
   if ( initialise_hoppet ) hoppet_init::assign( pdf1->pdf() );
   
   if ( ( !isSymmetric() && !isDISgrid() ) || ( pdf1 && pdf1!=pdf0 ) ) {
     
     for ( int itau=0 ; itau<m_Ntau ; itau++  ) {
  
       if ( itau<m_taufilledmin || itau>m_taufilledmax ) continue;
 
       double tau = gettau(itau);
       double Q2  = fQ2(tau);
       double Q   = std::sqrt(Q2); 
 
       /// alpha_s table has already been filled  
       /// m_alphas[itau] = alphas(rscale_factor*Q)*invtwopi;
       
       // y2 tables
       //    for ( int iy=iymin2 ; iy<=iymax2 ; iy++ ) { 
       for ( int iy=n_y2 ; iy-- ;  ) { 
 
 	if ( iy<m_y2filledmin || iy>m_y2filledmax ) continue;
 		
 	double y = gety2(iy);
 	double x = fx(y);
 	double fun = 1;
 	if (m_reweight) fun = weightfun(x);
 	
 	if ( scale_beams ) { 
 	  x *= beam_scale;	
 	  if ( x>=1 ) { 
 	    for ( int ip=0 ; ip<13 ; ip++ ) m_fg2[itau][iy][ip]=0; 
 	    if ( nloop>=1 && fscale_factor!=1 ) { 
 	      for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit2[itau][iy][ip] = 0;
 	      if (nloop==2)
 		for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit22[itau][iy][ip] = 0;
 	    }
 	    continue; 
 	  }
 	}
 
 	pdf1->evaluate(x, fscale_factor*Q, m_fg2[itau][iy]);
 	
 	double invx = 1/x;
 	for ( int ip=0 ; ip<13 ; ip++ ) m_fg2[itau][iy][ip] *= invx;
 	if ( m_reweight ) for ( int ip=0 ; ip<13 ; ip++ ) m_fg2[itau][iy][ip] *= fun;      
 	
 	// splitting functions
 	if ( nloop>=1 && fscale_factor!=1 ) { 
 
 	  double PLO[13]; for ( int ip=13 ; ip-- ; ) PLO[ip] = 0;  
 	  splitting( x, fscale_factor*Q, PLO, 1 ); /// PLO x pdf
 	  for ( int ip=13 ; ip-- ; ) m_fsplit2[itau][iy][ip] = PLO[ip]*invx;
 	  if ( m_reweight ) for ( int ip=13 ; ip-- ; ) m_fsplit2[itau][iy][ip] *= fun;
 	  
 	  if (nloop > 1) {
 	    //splitting( x, fscale_factor*Q, m_fsplit22[itau][iy], 1 );
  
 	    double P2LO[13]; for ( int ip=13 ; ip-- ; ) P2LO[ip] = 0;  
 	    double PNLO[13]; for ( int ip=13 ; ip-- ; ) PNLO[ip] = 0;  
 	    
 	    splitting( x, fscale_factor*Q, P2LO, 11 ); /// PLO x PLO x pdf
 	    splitting( x, fscale_factor*Q, PNLO, 2 );  /// PNLO x pdf
 	    
 	    for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit22[itau][iy][ip] = 0.5*logf*(beta0*(logf - 2*logr)*PLO[ip] + logf*P2LO[ip] - 2*PNLO[ip] )*invx ;
 	    
 	    if ( m_reweight ) for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit22[itau][iy][ip] *= fun;
 
 	  }
 	}
       }
 
 
     } // loop over itau
 
   } // isSymmetric()
 
 }
 
 
 
 
 #if 0
 void igrid::pdfinterp(double x, double Q2, double* f)
 {
   int k1=fk1(x);
   int k3=fkappa(Q2);
 
   double u_y1  = ( fy(x)-gety1(k1) )/deltay1();
   double u_tau = ( ftau(Q2)-gettau(k3))/deltatau();
 
   //  double rint = 0.;
 
   for ( int i=0 ; i<13 ; i++ ) f[i]=0.;
 
   double fI_factor;
 
   for( int i3=0 ; i3<=m_tauorder ; i3++ ) { // interpolation loop in Q2
    
     //    double fI3 = fI(i3, m_tauorder, u_tau);
 
     for( int i1=0 ; i1<=m_yorder ; i1++ ) { // interpolation loop in x1
       
       //      double fI1 = fI(i1, m_yorder, u_y1); 
 
       fI_factor=fI(i1, m_yorder, u_y1) * fI(i3, m_tauorder, u_tau);
                     
       for ( int ip1=0 ; ip1<13 ; ip1++ ) { 
 	f[ip1] += fI_factor*m_fg1[k3+i3][std::abs(k1+i1)][ip1];     
       }
     }
   }
    
   double fun = weightfun(x); 
 
   if ( m_reweight ) for ( int ip1=0 ; ip1<13 ; ip1++ ) f[ip1] /= fun;
 
 }
 #endif
 
 
 
 
 // takes pdf as the pdf lib wrapper for the pdf set for the convolution.
 // takes genpdf as a function to form the generalised parton distribution.
 // alphas is a function for the calculation of alpha_s (surprisingly)
 // lo_order is the order of the calculation, ie for inclusive, or two jets, 
 // it's O(alpha_s) (ie lo_order=1) and for 3 jet it would be O(alpha_s^2)
 // (ie lo_order=2)
 // nloop is the number of loops, ie for the leading order component it is
 // 0, for the nlo component, interestingly it is also 0, since the leading order
 // grids are seperate from the nlo grids, if nloop=1 then we must be calculating 
 // {r,f}scale_factor ie f*mu then *scale_factor=f
 // splitting, is the splitting function 
 
 double appl::igrid::convolute(NodeCache* pdf0,
 			      NodeCache* pdf1,
 			      appl_pdf*  genpdf,
 			      double (*alphas)(const double& ), 
 			      int     lo_order,  
 			      int     _nloop, 
 			      double  rscale_factor,
 			      double  fscale_factor,
 			      double Escale) 
 { 
 
   //  m_transvar = m_transvarlocal;
 
   if ( pdf1==0 ) pdf1 = pdf0; 
 
 
   m_conv_param.pdf0 = pdf0;
   m_conv_param.pdf1 = pdf1;
   m_conv_param.alphas = alphas;
 
 
   m_conv_param.lo_order =  lo_order;
   m_conv_param._nloop   = _nloop;
   m_conv_param.Escale   =  Escale;
 
   m_conv_param.rscale_factor = rscale_factor;
   m_conv_param.fscale_factor = fscale_factor;
 
   m_conv_param.genpdf = genpdf;
 
   m_conv_param.dsigma     = 0;
   m_conv_param.dsigmaNLO  = 0;
   m_conv_param.dsigmaNNLO = 0;
 
   double dsigma     = 0; 
   double dsigmaNLO  = 0; 
   double dsigmaNNLO = 0; 
 
 
 #ifndef PDFTHREAD
 
   //char name[]="appl_grid:igrid::convolute(): ";
   //  static const double twopi = 2*M_PI;
   //  static const int nc = 3;
   //  //TC   const int nf = 6;
   //  static const int nf = 5;
   //  static double beta0=(11.*nc-2.*nf)/(6.*twopi);
   //const bool debug=false;  
 
   //  double alphas_tmp = 0.;  
   //  double _alphas = 1.;
   //  double  alphaplus1 = 0.;
   // do the convolution  
   // if (debug) std::cout<<name<<" nloop= "<<nloop<<endl;
   //  std::cout << "\torder=" << lo_order << "\tnloop=" << nloop << std::endl;
   // is the grid empty
   int size=0;
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
     if ( m_weight[ip]->empty() ) continue;
     if ( !m_weight[ip]->trimmed() )  {
       // std::cerr << "igrid::convolute() naughty, naughty!" << std::endl;
       m_weight[ip]->trim();
     }
     //    size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
     size++;
   }
 
   // grid is empty
   if ( size==0 )  return 0;
 
   int nloop = std::fabs(_nloop); 
 
   setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale);
 
 #endif
 
 
   if ( threads_disabled ) convolute_internal();
   else                    process();            
 
   dsigma     = m_conv_param.dsigma; 
   dsigmaNLO  = m_conv_param.dsigmaNLO; 
   dsigmaNNLO = m_conv_param.dsigmaNNLO; 
 
   return dsigma;
 }
 
 
 
 
 
 
 void appl::igrid::convolute_internal() { 
 
 
   //  struct timeval mytimer = appl_timer_start();
 
   int     lo_order = m_conv_param.lo_order;  
   int     _nloop    = m_conv_param._nloop;  
 
   int      nloop = std::fabs(_nloop); 
 
   double  rscale_factor = m_conv_param.rscale_factor;  
   double  fscale_factor = m_conv_param.fscale_factor;  
 
   appl_pdf* genpdf = m_conv_param.genpdf;
 
   //  static const double twopi = 2*M_PI;
   //  static const int nc = 3;
   //  //TC   const int nf = 6;
   //  static const int nf = 5;
   //  static double beta0=(11.*nc-2.*nf)/(6.*twopi);
   //  const bool debug=false;  
 
 
   double dsigma     = 0.; //, xsigma = 0.;
   double dsigmaNLO  = 0.; //, xsigma = 0.;
   double dsigmaNNLO = 0.; //, xsigma = 0.;
 
   //  std::cout << "conv_internal:  \t" << name() << " m_Nproc: " << m_Nproc << std::endl;
 
   int size=0;
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
 
     ///    std::cout << "\t" << name() << " " << ip << " :: " << m_weight[ip]->empty() << std::endl; 
     
     if ( m_weight[ip]->empty() ) continue;
 
     if ( !m_weight[ip]->trimmed() )  {
       /// std::cerr << "igrid::convolute() naughty, naughty!" << std::endl;
       m_weight[ip]->trim();
     }
     //    size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
     size++;
   }
 
   // grid is empty
   if ( size==0 )  return;
 
 
 #ifdef PDFTHREAD
 
   double (*alphas)(const double& ) = m_conv_param.alphas;   
 
   NodeCache* pdf0 = m_conv_param.pdf0;   
   NodeCache* pdf1 = m_conv_param.pdf1;   
   
   double   Escale =  m_conv_param.Escale;   
 
   //  double _alphas = 1.;
   //  double  alphaplus1 = 0.;
   // do the convolution  
   // if (debug) std::cout<<name<<" nloop= "<<nloop<<endl;
   //  std::cout << "\torder=" << lo_order << "\tnloop=" << nloop << std::endl;
   // is the grid empty
 
   //  std::cout << "conv_internal:  \t" << name() << " m_Nproc: " << m_Nproc << std::endl;
 
 
   int size=0;
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
 
     //    std::cout << "\t" << name() << " " << ip << " :: " << m_weight[ip]->empty() << std::endl; 
     
     if ( m_weight[ip]->empty() ) continue;
 
     if ( !m_weight[ip]->trimmed() )  {
       /// std::cerr << "igrid::convolute() naughty, naughty!" << std::endl;
       m_weight[ip]->trim();
     }
     size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
   }
 
   // grid is empty
   if ( size==0 )  return;
 
 
   setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale);
 
 
 #endif
 
 
   double* sig = new double[m_Nproc];  // weights from grid
 
   double* H   = new double[m_Nproc];  // generalised pdf  
 
   double* HA  = 0;  // generalised splitting functions
   double* HB  = 0;  // generalised splitting functions
 
   double* HA2 = 0;  // generalised splitting functions at NNLO
   double* HB2 = 0;  // generalised splitting functions at NNLO
   double* HAB = 0;  // generalised splitting functions at NNLO
 
   if ( nloop>=1 && fscale_factor!=1 ) { 
     HA  = new double[m_Nproc];  // generalised splitting functions
     HB  = new double[m_Nproc];  // generalised splitting functions
     if (nloop == 2 ) {
       HA2  = new double[m_Nproc];  // generalised splitting functions at NNLO
       HB2  = new double[m_Nproc];  // generalised splitting functions at NNLO
       HAB  = new double[m_Nproc];  // generalised splitting functions at NNLO
     }
   }
   
   // cross section for this igrid  
 
 
 
   
   //char name[]="appl_grid:igrid::convolute(): ";
   static const double twopi = 2*M_PI;
   static const int nc = 3;
   //TC   const int nf = 6;
   static const int nf = 5;
   static double beta0=(11.*nc - 2.*nf)/(6.*twopi);
 
   static double beta1=(34.*nc*nc - 3.*(nc-1./nc)*nf - 10.*nc*nf)/(6.*twopi);
 
 
   double alphas_tmp = 0.;  
 
   double _alphas = 1.;
   double  alphaplus1 = 0.;
   double  alphaplus2 = 0.;
 
 
   m_conv_param.dsigma     = 0;
   m_conv_param.dsigmaNLO  = 0;
   m_conv_param.dsigmaNNLO = 0;
 
 
   // loop over the grid 
   // 
 
   for ( int itau=0 ; itau<Ntau() ; itau++  ) {
     _alphas  = 1;    
 
     alphas_tmp = m_alphas[itau];
 
 
     for ( int iorder=0 ; iorder<lo_order ; iorder++ ) _alphas *= alphas_tmp;
     alphaplus1 = _alphas*alphas_tmp;
     alphaplus2 = alphaplus1*alphas_tmp;
 
 
     //    for ( int iy1=0 ; iy1<Ny1() ; iy1++ ) {            
     //      for ( int iy2=0 ; iy2<Ny2() ; iy2++ ) { 
     for ( int iy1=Ny1() ; iy1-- ;  ) {            
 
       for ( int iy2=Ny2() ; iy2-- ;  ) { 
  	// test if this element is actually filled
 	// if ( !m_weight[0]->trimmed(itau,iy1,iy2) ) continue; 
 	bool nonzero = false;
 	// basic convolution order component for either the born level
 	// or the convolution of the nlo grid with the pdf 
 	
 
 	for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
 	  if ( (sig[ip] = (*(const SparseMatrix3d*)m_weight[ip])(itau,iy1,iy2)) ) nonzero = true;
 	}
 	
 	//	for ( int ip=0 ; ip<m_Nproc ; ip++ ) std::cout << "\t" << sig[ip]; 
 	//	std::cout << std::endl;
 
 	if ( nonzero ) { 	
 	  // build the generalised pdfs from the actual pdfs
 
 	  genpdf->evaluate( m_fg1[itau][iy1],  m_fg2[itau][iy2], H );
 
 	  // do the convolution
           double xsigma=0.;
 
 	  if ( m_parent && m_parent->subproc()!=-1 ) { 
 	    int ip=m_parent->subproc();
 	    xsigma+= sig[ip]*H[ip];
 	  }
 	  else { 
 	    for ( int ip=0 ; ip<m_Nproc ; ip++ ) xsigma+= sig[ip]*H[ip];
 	  }
 
 	  /// if want NLO part only, don't add in the born term
 	  if ( _nloop >= 0 ) dsigma += _alphas*xsigma;
 
 	  // now do the convolution for the variation of factorisation and 
 	  // renormalisation scales, proportional to the leading order weights
 	  if ( nloop>=1 ) { 
 
 	    // renormalisation scale dependent bit
 	    
 	    if ( rscale_factor!=1 ) { 
 	      // nlo relative ln mu_R^2 term 
 
 	      double logr2 = std::log(rscale_factor*rscale_factor);
 	      double fac   = alphaplus1*twopi*lo_order*logr2*xsigma;
 	      
 	      dsigmaNLO += fac*beta0;
 	      
 	      if (nloop>1) { 
 		dsigmaNNLO += fac*_alphas*( beta1 + 0.5*beta0*beta0*(1+lo_order)*logr2 );
 		
 		//		  dsigmaNNLO += alphaplus2*twopi*((beta01*lo_order*log(rscale_factor*rscale_factor) + 
 		//						    (0.5*beta0*beta0*lo_order*(1+lo_order)*logr2*logr2 )*xsigma;
 	      }
 
   	    }
 
 	    
 	    // factorisation scale dependent bit
             
 	    // nlo relative ln mu_F^2 term 
 	    if ( fscale_factor!=1 ) {
 	      genpdf->evaluate( m_fg1    [itau][iy1],  m_fsplit2[itau][iy2], HA);
 	      genpdf->evaluate( m_fsplit1[itau][iy1],  m_fg2    [itau][iy2], HB);
 	      xsigma=0.; /// setting this to 0 ??? maybe this is ok because we have already used it for the rscale part
 
 	      if ( m_parent && m_parent->subproc()!=-1 ) { 
 		int ip=m_parent->subproc();
 		xsigma += sig[ip]*(HA[ip]+HB[ip]);
 	      }
 	      else { 
 		for ( int ip=0 ; ip<m_Nproc ; ip++ ) xsigma += sig[ip]*(HA[ip]+HB[ip]);
 	      }
 
 	      double logf2 = std::log(fscale_factor*fscale_factor);
 
 	      double fac = alphaplus1*logf2*xsigma;
 
 	      dsigmaNLO -= fac;
 
 	      if ( nloop>1 ) {
 
 		double logr2 = std::log(rscale_factor*rscale_factor);
 		
 		dsigmaNNLO -= fac*(1+beta0*logr2*lo_order*_alphas);
 
 		/// plus the double splitting function components ...
 		/// PDFA2 PDFB0 alphas2 +PDFA1 PDFB1 alphas2 + PDFA0 PDFB2 alphas2 
 
 		genpdf->evaluate( m_fsplit12[itau][iy1],        m_fg2[itau][iy2], HA2 );
 		genpdf->evaluate(      m_fg1[itau][iy1],   m_fsplit22[itau][iy2], HB2 );
 		genpdf->evaluate(  m_fsplit1[itau][iy1],    m_fsplit2[itau][iy2], HAB );
 
 		double fac2 = 0;
 
 	        if ( m_parent && m_parent->subproc()!=-1 ) { 
 		    int ip=m_parent->subproc();
 		    fac2 += sig[ip]*( HA2[ip] + HB2[ip] + HAB[ip] );
 		}
 		else { 
 		    for ( int ip=0 ; ip<m_Nproc ; ip++ ) fac2 += sig[ip]*( HA2[ip] + HB2[ip] + HAB[ip] );
 		}
 		 
 		xsigma += alphaplus1*logf2*xsigma*fac2;  /// can this be correct ? scaling the current
 
 		dsigmaNNLO -= xsigma;
 	      }
 
 	      //if (debug) 
               //cout <<name<<" fscale= " << fscale_factor << " dsigma= "<<dsigma << std::endl;
 	    }
 	  }
 	}  // nonzero
       }  // iy2
     }  // iy1
 
 
 
   }  // itau
 
   
   //if (debug)  std::cout << name<<"     convoluted dsigma=" << dsigma << std::endl; 
   
   delete[] sig;
   delete[] H;
   delete[] HA;
   delete[] HB;
 
   if ( HA2 ) { 
     delete[] HA2;
     delete[] HB2;
     delete[] HAB;
   }    
 
   deletepdftable();
   
   //  std::cout << "dsigma " << dsigma << std::endl;
 
   // NB!!! the return value dsigma must be scaled by Escale*Escale which 
   // is done in grid::vconvolute. It would be better here, but is reduces 
   // the number of operations if in grid. 
 
 
   m_conv_param.dsigma     = dsigma;
   m_conv_param.dsigmaNLO  = dsigmaNLO;
   m_conv_param.dsigmaNNLO = dsigmaNNLO;
 
   //  double mytime = appl_timer_stop(mytimer);
 
   //  std::printf("thread done: param: %lf internal time %lf ms\n", m_conv_param.dsigma, mytime );
 
 
 
 }
 
 
 
 
 
 
 
 /// this is the convolute routine for the amcatnlo convolution - essentially it 
 /// is the same as for the standard calculation, but the amcatnlo calculation
 /// stores weights for the NLO born contribution, and counterterms, so we need
 /// more grids than the usual two, 
 double appl::igrid::amc_convolute(NodeCache* pdf0,
 				  NodeCache* pdf1,
 				  appl_pdf*  genpdf,
 				  double (*alphas)(const double& ), 
 				  int     lo_order,  
 				  int     nloop, 
 				  double  rscale_factor,
 				  double  fscale_factor,
 				  double Escale) 
 { 
 
   m_conv_param.pdf0 = pdf0;
   m_conv_param.pdf1 = pdf1;
 
   m_conv_param.alphas = alphas;
 
   m_conv_param.lo_order =  lo_order;
 
   //  m_conv_param._nloop   =  nloop;
   //  m_conv_param.Escale   =  Escale;
 
   //  m_conv_param.rscale_factor = rscale_factor;
   //  m_conv_param.fscale_factor = fscale_factor;
 
   m_conv_param.genpdf = genpdf;
 
   m_conv_param.dsigma = 0;
 
   //  double dsigma  = 0; 
 
 
   // is the grid empty
   int size=0;
    for ( int ip=0 ; ip<m_Nproc ; ip++ ) { 
     if ( !m_weight[ip]->trimmed() )  {
       //  std::cout << "igrid::convolute() naughty, naughty!" << std::endl;
       m_weight[ip]->trim();
     }
     size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
   }
 
   // grid is empty
   if ( size==0 )  return 0;
 
   setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale);  
 
   if ( threads_disabled ) amc_convolute_internal();
   else                    process();
 
   return m_conv_param.dsigma;
 }
 
 
 
 
 
 
 
 void appl::igrid::amc_convolute_internal() { 
   
   int     lo_order = m_conv_param.lo_order;  
 
   //  int     _nloop    = m_conv_param._nloop;  
   
   //  int      nloop = std::fabs(_nloop); 
   
   //  double  rscale_factor = m_conv_param.rscale_factor;  
   //  double  fscale_factor = m_conv_param.fscale_factor;  
   
   appl_pdf* genpdf = m_conv_param.genpdf;
   
   static const double eightpisquared = 8*M_PI*M_PI;
   
   double alphas_tmp = 0.;  
   double dsigma  = 0.; //, xsigma = 0.;
   double _alphas = 1.;
 
   double* sig = new double[m_Nproc];  // weights from grid
   double* H   = new double[m_Nproc];  // generalised pdf  
   double* HA  = 0;  // generalised splitting functions
   double* HB  = 0;  // generalised splitting functions
   //  if ( nloop==1 && fscale_factor!=1 ) { 
   //    HA  = new double[m_Nproc];  // generalised splitting functions
   //    HB  = new double[m_Nproc];  // generalised splitting functions
   //  }
 
   // cross section for this igrid  
 
   // loop over the grid 
   // 
 
   for ( int itau=0 ; itau<Ntau() ; itau++  ) {
     _alphas  = 1;    
     alphas_tmp = m_alphas[itau]*eightpisquared;
     for ( int iorder=0 ; iorder<lo_order ; iorder++ ) _alphas *= alphas_tmp;
     //   alphaplus1 = _alphas*alphas_tmp;
 
     for ( int iy1=Ny1() ; iy1-- ;  ) {            
       for ( int iy2=Ny2() ; iy2-- ;  ) { 
  	// test if this element is actually filled
 	// if ( !m_weight[0]->trimmed(itau,iy1,iy2) ) continue; 
 	bool nonzero = false;
 	// basic convolution order component for either the born level
 	// or the convolution of the nlo grid with the pdf 
 	for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
 	  if ( (sig[ip] = (*(const SparseMatrix3d*)m_weight[ip])(itau,iy1,iy2)) ) nonzero = true;
 	}
 	
 	//	for ( int ip=0 ; ip<m_Nproc ; ip++ ) std::cout << "\t" << sig[ip]; 
 	//	std::cout << std::endl;
 
 	if ( nonzero ) { 	
 
 	  // build the generalised pdfs from the actual pdfs
 	  genpdf->evaluate( m_fg1[itau][iy1],  m_fg2[itau][iy2], H );
 	
 	  // do the convolution
 
           double xsigma=0.;
 	  for ( int ip=0 ; ip<m_Nproc ; ip++ ) xsigma+=sig[ip]*H[ip];
 	  dsigma += _alphas*xsigma;
 
 	}  // nonzero
       }  // iy2
     }  // iy1
   }  // itau
   
   //if (debug)  std::cout << name<<"     convoluted dsigma=" << dsigma << std::endl; 
   
   delete[] sig;
   delete[] H;
   delete[] HA;
   delete[] HB;
   
   deletepdftable();
   
   //  std::cout << "dsigma " << dsigma << std::endl;
 
   m_conv_param.dsigma = dsigma;
 
   // NB!!! the return value dsigma must be scaled by Escale*Escale which 
   // is done in grid::vconvolute. It would be better here, but is reduces 
   // the number of operations if done in grid. 
   return; 
 }
 
 
 
 
 
 
 bool appl::igrid::shrink( const std::vector<int>& keep ) {
  
   /// save the old grids
   int          Nproc = m_Nproc;
   SparseMatrix3d** w = m_weight;
 
   /// resize
   m_Nproc  = keep.size();
   m_weight = new SparseMatrix3d*[m_Nproc];
 
   /// copy across the grids we want to save
   for ( unsigned ip=0 ; ip<keep.size() ; ip++ ) {
  
     //    std::cout << "\tcopy " << keep[ip] << " -> " << ip << std::endl; 
 
     m_weight[ip] = w[keep[ip]]; /// move across the processes we want to keep
     w[keep[ip]] = 0; /// now flag as 0 
   }
 
   for ( int ip=0 ; ip<Nproc ; ip++ ) if ( w[ip]!=0 ) delete w[ip];
 
   delete[] w; 
   
   return true;
 }
 
 
 
 
 // FIXME: this has lots of legacy code and stuff commented which
 // is no longer needed or even correct, so it should be tidied.
 // the point of the algorithm here is to find the extent of the filled
 // bines + 1 on either side and create a new grid with these limits
    
 void appl::igrid::optimise(int NQ2, int Nx1, int Nx2, int extrabins) {     
 
   //  std::cout << "\tsize(untrimmed)=" << m_weight[0]->size();
 
   //  std::cout << "ymin=" << gety(0) << "\tymax=" << gety(m_Ny-1) 
   //       << "\txmin=" << fx(gety(m_Ny-1)) << "\txmax=" << fx(gety(0)) << std::endl;
 
   for ( int i=0 ; i<m_Nproc ; i++ )  m_weight[i]->trim();
 
   //  std::cout << "\tsize(trimmed)=" << m_weight[0]->size() << std::endl;
 
 #if 1
 
   // overall igrid optimisation limits
 
   int _y1setmin = Ny1();
   int _y1setmax = -1;
   
   int _y2setmin = Ny2();
   int _y2setmax = -1;
   
   int _tausetmin = Ntau(); 
   int _tausetmax = -1; 
   
   //  double oldy1min = m_y1min;
   //  double oldy1max = m_y1max;
   
   //  double oldy2min = m_y2min;
   //  double oldy2max = m_y2max;
 
   // go through all the subprocess to get the limits
   // FIXME: this is actually redundant, at the moment, it assumes all the subprocesses 
   //        occupy the same phase space, so only the first need be tested, and if
   //        they don;t all occupy the same phase space, then they should each be
   //        optimised seperately
 
   // find limits for this grid
   // initialise to original grid limits
   int y1setmin = _y1setmin;
   int y1setmax = _y1setmax;
   
   int y2setmin = _y2setmin;
   int y2setmax = _y2setmax;
   
   int tausetmin = _tausetmin; 
   int tausetmax = _tausetmax; 
   
   bool firstdebug = true;
   bool firstempty = true;
 
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) { 
     
     if ( firstdebug ) std::cout << "optimise() proc " << ip 
 				<< "\tfilled Q2 weights " << m_weight[ip]->xmax() << " - " << m_weight[ip]->xmin()
 				<< "\tx1 " << m_weight[ip]->ymax() << " - " << m_weight[ip]->ymin()
 				<< "\tx2 " << m_weight[ip]->zmax() << " - " << m_weight[ip]->zmin();
 
     // is it empty?
     if ( m_weight[ip]->xmax()-m_weight[ip]->xmin()+1 == 0 ) { 
       if ( firstdebug && firstempty ) std::cout << "\tempty" << std::endl;
       firstempty = false;
       firstdebug = false;
       continue;
     }
     
 
     firstdebug = false;
 
     //    m_weight[ip]->print();
 
     //    m_weight[ip]->yaxis().print(fx);
 
     //    std::cout << "initial ip=" << ip 
     //	       << "\tm_y2min=" << m_y2min << "\tm_y2max=" << m_y2max
     //         << "\tm_y1min=" << m_y1min << "\tm_y1max=" << m_y1max  
     //         << std::endl;  
 
     //    m_weight[ip]->print();
 
     // find actual limits
 
     // y1 optimisation
     int ymin1 = m_weight[ip]->ymin();
     int ymax1 = m_weight[ip]->ymax();
     
     if ( ymin1<y1setmin )                 y1setmin = ymin1;
     if ( ymin1<=ymax1 && y1setmax<ymax1 ) y1setmax = ymax1; 
 
     //    std::cout << "ip=" << ip << std::endl; // << "\tymin1=" << ymin1 << "\tymax1=" << ymax1 << std::endl;
     
     // y2 optimisation
     int ymin2 = m_weight[ip]->zmin();    
     int ymax2 = m_weight[ip]->zmax();
 
     if ( ymin2<y2setmin )                 y2setmin = ymin2;
     if ( ymin2<=ymax2 && y2setmax<ymax2 ) y2setmax = ymax2; 
     
     // tau optimisation
     int taumin = m_weight[ip]->xmin();    
     int taumax = m_weight[ip]->xmax();
 
     if ( taumin<tausetmin )                   tausetmin = taumin;
     if ( taumin<=taumax && taumax>tausetmax ) tausetmax = taumax;
 
   }
   
   // if grid is empty, do "nothing" ie create the grid with the same
   // limits as before but with the new required number of bins 
   if (  y1setmax==-1 || y2setmax==-1 || tausetmax==-1 ) { 
     m_Ny1  = Nx1;
     m_Ny2  = Nx2;
     m_Ntau = NQ2;
     m_deltay1  = (m_y1max-m_y1min)/(m_Ny1-1);
     m_deltay2  = (m_y2max-m_y2min)/(m_Ny2-1);
     m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
   }
   else { 
  
 
     
     // y1 optimisation
     //    double oldy1min = m_y1min;
     //    double oldy1max = m_y1max;
    
 
     if ( isOptimised() ) { 
       // add a bit on each side
       y1setmin -= extrabins;
       y1setmax += extrabins;
     }
     else { 
       // not optimised yet so filled with phase space only so add 
       // the order to the max filled grid element position also 
       y1setmin -= extrabins;
       y1setmax += m_yorder+extrabins;
     }
     
     if ( y1setmin<0 )       y1setmin = 0;
     if ( y1setmax>m_Ny1-1 ) y1setmax = m_Ny1-1;
     
     
     double _min = gety1(y1setmin); 
     double _max = gety1(y1setmax); 
     
     m_Ny1     = Nx1;
     m_y1min   = _min;
     m_y1max   = _max;
     m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
     
     
     
     // y2 optimisation
     //   double oldy2min = m_y2min;
     //   double oldy2max = m_y2max;
     
     if ( isOptimised() ) { 
       // add a bit on each side
       y2setmin -= extrabins;
       y2setmax += extrabins;
     }
     else { 
       // not optimised yet so add the order to the  
       // max filled grid element position also 
       y2setmin -= extrabins;
       y2setmax += m_yorder+extrabins;
     }
     
     if ( y2setmin<0 )      y2setmin = 0;
     if ( y2setmax>=m_Ny2 ) y2setmax = m_Ny2-1; 
 
 
     //  m_Ny2 =  y2setmax-y2setmin+1;
     _min = gety2(y2setmin); 
     _max = gety2(y2setmax); 
     
     m_Ny2     = Nx2;
     m_y2min   = _min;
     m_y2max   = _max;
     m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
     
     std::cout << "\t-> " 
 	      << y1setmin  << " - " << y1setmax << " : " 
 	      << y2setmin  << " - " << y2setmax 
 	      << std::endl;
     
     // tau optimisation
     //   double oldtaumin = m_taumin;
     //   double oldtaumax = m_taumax;
     
     // add a bit on each side
     if ( isOptimised() ) { 
       // add a bit on each side
       tausetmin--;
       tausetmax++;
     }
     else { 
       // not optimised yet so add the order to the  
       // max filled grid element position also 
       tausetmax += m_tauorder+1;
       tausetmin -= 1;
     }  
 
     if ( tausetmin<0 )         tausetmin = 0;
     if ( tausetmax>=m_Ntau-1 ) tausetmax = m_Ntau-1;
 
     
     _min   = gettau(tausetmin); 
     _max   = gettau(tausetmax);
     
     m_Ntau     = NQ2;
     m_taumin   = _min;
     m_taumax   = _max;
     m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
     
     //    std::cout << "done ip=" << ip << std::endl; 
     
   }
   
 #endif
 
   // now create the new subprocess grids with optimised limits
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) { 
     delete m_weight[ip];
     m_weight[ip] = new SparseMatrix3d(m_Ntau, m_taumin,   m_taumax,    
 				      m_Ny1,   m_y1min,   m_y1max, 
 				      m_Ny2,   m_y2min,   m_y2max ); 
   }   
   
   m_optimised = true;
 }
 
 
 
 
 
 // numerical operators
 appl::igrid& appl::igrid::operator=(const appl::igrid& g) { 
   m_Ny1     = g.m_Ny1;
   m_y1min   = g.m_y1min;
   m_y1max   = g.m_y1max;
   m_deltay1 = g.m_deltay1;
 
   m_Ny2     = g.m_Ny2;
   m_y2min   = g.m_y2min;
   m_y2max   = g.m_y2max;
   m_deltay2 = g.m_deltay2;
 
   m_yorder = g.m_yorder;
    
   m_Ntau     = g.m_Ntau;
   m_taumin   = g.m_taumin;
   m_taumax   = g.m_taumax;
   m_deltatau = g.m_deltatau;
 
   m_tauorder = g.m_tauorder;
  
   m_fg1      = 0;
   m_fg2      = 0;
   
   //  construct();
 
   //  m_weight = new SparseMatrix3d*[m_Nproc];
 
   for( int ip=0 ; ip<m_Nproc ; ip++ ) { 
     // delete old sparse matrix
     delete m_weight[ip];
     // create new
     m_weight[ip] = new SparseMatrix3d(*g.m_weight[ip]);
   }
 
   return *this;
 }
 
 
 // should really check all the limits and *everything* is the same
 appl::igrid& appl::igrid::operator+=(const igrid& g) { 
   for ( int ip=0 ; ip<m_Nproc ; ip++ ) {
     if ( m_weight[ip] && g.m_weight[ip] ) { 
       /// this is a complicated conditional:
       /// do nothing if the grid to be added is empty ...
       if ( !g.m_weight[ip]->empty() ) { 
 	if (  m_weight[ip]->empty() ) { 
 	  delete m_weight[ip];
 	  m_weight[ip] = new SparseMatrix3d( *g.m_weight[ip] );
 	}
 	else if ( m_weight[ip]->compare_axes( *g.m_weight[ip] ) ) (*m_weight[ip]) += (*g.m_weight[ip]); 
 	else { 
 	  throw exception("igrid::operator+=() grids do not match");
 	}
       }
     }
   }
   return *this;
 } 
 
 
 std::ostream& appl::igrid::header(std::ostream& s) const { 
   //  s << "interpolation orders: x=" << g.yorder() << "\tQ2=" << g.tauorder() << std::endl;  
   s << "\t x:  [ "  << std::setw(2) 
     //    << Ny1() << " ;\t" << std::setw(7) << fx(y1max())  << " - " << std::setw(7) << std::setprecision(6) << fx(y1min()) << "\t : " 
     //    << Ny2() << " ;\t" << std::setw(7) << fx(y2max())  << " - " << std::setw(7) << std::setprecision(6) << fx(y2min()) 
     << Ny1() << " :\t " << std::setw(7) << std::setprecision(6) << fx(y1max())  << " - " << std::setw(7) << std::setprecision(6) << fx(y1min()) << "\t : " 
     << Ny2() << " :\t " << std::setw(7) << std::setprecision(6) << fx(y2max())  << " - " << std::setw(7) << std::setprecision(6) << fx(y2min()) 
     << "\t : " << "\t( order=" << yorder()   << " ) ]"; 
   s << "\t Q2: [ " 
     << Ntau() << " :\t "  << std::setw(7) << std::setprecision(6) << fQ2(taumin()) << " - " << std::setw(7) << std::setprecision(6) << fQ2(taumax()) 
     << "\t( order=" << tauorder() << "  - reweight " << ( m_reweight ? "on " : "off" ) << ") ]";
   return s;
 }
 
 
 std::ostream& operator<<(std::ostream& s, const appl::igrid& g) {
   return g.header(s);
 }
 
 
 
 void appl::igrid::run_thread() {
  
   lock_proc();
   mprocessing = false;
   unlock_proc();
 
 
   while( true ) { 
 
     //    std::cout << "thread " << this << " " << mname << " waiting ... " << std::endl; 
 
     /// supend immediately
     suspend();
 
     if ( mterminate ) break;
 
     //    std::cout << "thread " << this << " " << mname << " running ... " << std::endl; 
 
     //    struct timeval mytimer = appl_timer_start();
 
     /// put the actual convoluting steps in here  
     /// once round the colbolution step, it puts itself to 
     /// sleep again
     
     if ( parent()->calculation()==grid::AMCATNLO )  amc_convolute_internal();
     else                                            convolute_internal();
 
     //    std::printf("thread done: param: %lf internal\n", m_conv_param.dsigma );
  
     /// after starting the processing step, the calling thread must 
     /// wait for all the threads to finish
 
     //    double mytime = appl_timer_stop(mytimer);
 
     //    std::printf("thread done: param: %lf internal time %lf ms\n", m_conv_param.dsigma, mytime );
 
     /// signal computation has finished
 
     //    std::cout << "thread " << this << " " << mname << " done [" << mytime << " ms" << "]" << std::endl; 
   }
 
   suspend(false);
   
 }
 
 
Index: applgrid/trunk/src/appl_grid.cxx
===================================================================
--- applgrid/trunk/src/appl_grid.cxx	(revision 1752)
+++ applgrid/trunk/src/appl_grid.cxx	(revision 1753)
@@ -1,2954 +1,2953 @@
 
 //  appl_grid.cxx        
 
 //   grid class - all the functions needed to create and 
 //   fill the grid from an NLO calculation program. 
 //  
 //  Copyright (C) 2007 Mark Sutton (sutt@hep.ucl.ac.uk)    
 
 // $Id: appl_grid.cxx, v1.00 2007/10/16 17:01:39 sutt $
 
 #include <cstdlib>
 #include <stdio.h>
 #include <sys/stat.h>
 
 #include <vector>
 #include <map>
 #include <set>
 #include <iostream>
 #include <fstream>
 #include <iomanip>
 #include <cmath>
 
 #include "appl_grid/appl_pdf.h"
 #include "appl_grid/appl_timer.h"
 #include "appl_grid/Directory.h"
 #include "appl_grid/appl_grid.h"
 
 #include "appl_grid/generic_pdf.h"
 #include "appl_grid/lumi_pdf.h"
 
 #include "appl_igrid.h"
 #include "Cache.h"
 
 
 #include "TFileString.h"
 #include "TFileVector.h"
 
 #include "TFile.h"
 #include "TObjString.h"
 #include "TVectorT.h"
 
 
 #include "amconfig.h"
 
 #include "Splitting.h"
 
 
 /// in case of non c++11 compliant compilation, define
 /// some c++11 type replacemant functions
 
 #if __cplusplus <= 199711L
 // #warning Not using C++11 compilation - using substitute functions
 
 namespace std {
 
 int stoi( const std::string s ) { return atoi(s.c_str()); }
 
 std::string to_string( int i ) {
     char a[16];
     sprintf( a, "%d", i );
     return a;
 }
 
 
 std::string to_string( double f ) {
     char a[178];
     sprintf( a, "%lf", f );
     return a;
 }
 
 }
 
 #endif
 
 
 
 template<typename T>
 std::ostream& operator<<(std::ostream& s, const std::vector<T>& v) {
   for ( unsigned i=0 ; i<v.size() ; i++ ) s << "\t" << v[i];
   return s;
 }
 
 
 static bool appl_first = true;
 
 
 /// this is a compatability flag for persistent versions 
 /// of the grid
 /// NB: ONLY change the major version if the persistent 
 ///     class changes in a non-backwards compatible way.
 
 // const std::string appl::grid::m_version = "version-3.3";
 const std::string appl::grid::m_version = PACKAGE_VERSION;
 
 std::string appl::grid::appl_version() const { return PACKAGE_VERSION; }
 
 #include "hoppet_init.h"
 
 #ifdef HAVE_HOPPET
 
 #include "hoppet_v1.h"
 
 // include hoppet splitting function code
 
 static appl::hoppet_init* hoppet = 0;
 
 void Splitting(const double& x, const double& Q, double* xf, int nLoops) {
   //  static const int nLoops    = 1;
   static const int nFlavours = 5;
   hoppetevalsplit_( x, Q, nLoops, nFlavours, xf); 
   return;
 }
 
 #else
 
 void Splitting(const double& , const double& , double* , int ) {
   throw appl::grid::exception( "hoppet library not included - cannot call splitting function"  ); 
   return; // technically, don't need this - should throw an exception
 }
 
 #endif
 
 
 #define appl_line() std::cout << __FILE__ << " " << __LINE__ << std::endl;
 
 
 
 /// helper function
 static bool contains(const std::string& s, const std::string& reg ) { 
   return s.find(reg)!=std::string::npos;
 }
 
 /// make sure pdf std::map is initialised
 // bool pdf_ready = appl::appl_pdf::create_map(); 
 
 std::string appl::compiled() {  return __DATE__; }
 std::string appl::version()  { return VERSION; }
 
 
 /// current time
 std::string appl::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;
 }
 
 
 /// simple test if a file exists
 bool appl::file_exists(const std::string& s) {
   struct stat sb;
   if ( stat( s.c_str(), &sb)==0 ) return true; // && S_ISREG(sb.st_mode ))
   else return false;
 }
 
 
 
 
 appl::grid::grid(int NQ2, double Q2min, double Q2max, int Q2order, 
 		 int Nx,  double xmin,  double xmax,  int xorder,
 		 int Nobs,  double obsmin, double obsmax, 
 		 std::string genpdfname,
 		 int leading_order, int _nloops, 
 		 std::string transform ) :
   m_leading_order(leading_order), m_order(_nloops+1), 
   m_run(0), m_optimised(false), m_trimmed(false), m_normalised(false), m_symmetrise(false), 
   m_transform(transform), m_genpdfname(genpdfname), m_cmsScale(0), m_dynamicScale(0),
   m_applyCorrections(false),
   m_documentation(""),
   m_type(STANDARD),
   m_read(false),
   m_subproc(-1),
   m_bin(-1),
   m_genwithpdf("")
 {
   // Initialize histogram that saves the correspondence obsvalue<->obsbin
   m_obs_bins=new TH1D("referenceInternal","Bin-Info for Observable", Nobs, obsmin, obsmax);
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2(); /// grrr root is so rubbish - not scaling errors properly
 
   m_obs_bins_combined = m_obs_bins;
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf
   //  if      ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   //  else if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
   findgenpdf( m_genpdfname );
 
   construct(Nobs, NQ2, Q2min, Q2max, Q2order, Nx, xmin, xmax, xorder, m_order, m_transform); 
 }
 
 
 
 
 appl::grid::grid(int Nobs, const double* obsbins, 
 		 int NQ2,  double Q2min, double Q2max, int Q2order, 
 		 int Nx,   double xmin, double xmax,   int xorder, 
 		 std::string genpdfname, 
 		 int leading_order, int _nloops, 
 		 std::string transform ) :
   m_leading_order(leading_order), m_order(_nloops+1), 
   m_run(0), m_optimised(false), m_trimmed(false),  m_normalised(false), m_symmetrise(false),
   m_transform(transform), m_genpdfname(genpdfname), m_cmsScale(0), m_dynamicScale(0),
   m_applyCorrections(false),
   m_documentation(""),
   m_type(STANDARD),
   m_read(false),
   m_subproc(-1),
   m_bin(-1),
   m_genwithpdf("")
 {
   
   // Initialize histogram that saves the correspondence obsvalue<->obsbin
   m_obs_bins=new TH1D("referenceInternal","Bin-Info for Observable", Nobs, obsbins);
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2(); /// grrr root is so rubbish - not scaling errors properly
 
   m_obs_bins_combined = m_obs_bins;
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf
   //  if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
 
   findgenpdf( m_genpdfname );
 
   construct(Nobs, NQ2, Q2min, Q2max, Q2order, Nx, xmin, xmax, xorder, m_order, m_transform );
 }
 
 
 
 
 appl::grid::grid(const std::vector<double>& obs, 
 		 int NQ2, double Q2min, double Q2max, int Q2order,
 		 int Nx,  double xmin,  double xmax,  int xorder, 
 		 std::string genpdfname,
 		 int leading_order, int _nloops, 
 		 std::string transform )  :
   m_leading_order(leading_order), m_order(_nloops+1), 
   m_run(0), m_optimised(false), m_trimmed(false), m_normalised(false), m_symmetrise(false),  
   m_transform(transform), m_genpdfname(genpdfname), m_cmsScale(0), m_dynamicScale(0),
   m_applyCorrections(false),
   m_documentation(""),
   m_type(STANDARD),
   m_read(false),
   m_subproc(-1),
   m_bin(-1),
   m_genwithpdf("")
 {
   
   if ( obs.size()==0 ) { 
     std::cerr << "grid::not enough bins in observable" << std::endl;
     std::exit(0);
   } 
   
   //  double* obsbins = new double[obs.size()];  
   //  for ( unsigned i=0 ; i<obs.size() ; i++ ) obsbins[i] = obs[i];
   //  int Nobs = obs.size()-1;
 
   // Initialize histogram that saves the correspondence obsvalue<->obsbin
   m_obs_bins=new TH1D("referenceInternal","Bin-Info for Observable", obs.size()-1, &obs[0] );
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2(); /// grrr root is so rubbish - not scaling errors properly
   //  delete[] obsbins;
 
   m_obs_bins_combined = m_obs_bins;
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf
   //   if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
   findgenpdf( m_genpdfname );
 
   construct( obs.size()-1, NQ2, Q2min, Q2max, Q2order, Nx, xmin, xmax, xorder, m_order, m_transform); 
 }
 
 
 
 appl::grid::grid(const std::vector<double>& obs, 
 		 std::string genpdfname,
 		 int leading_order, int _nloops, 
 		 std::string transform )  :
   m_leading_order(leading_order), m_order(_nloops+1), 
   m_run(0), m_optimised(false), m_trimmed(false), m_normalised(false), m_symmetrise(false),  
   m_transform(transform), m_genpdfname(genpdfname), m_cmsScale(0), m_dynamicScale(0),
   m_applyCorrections(false),  
   m_documentation(""),
   m_type(STANDARD),
   m_read(false),
   m_subproc(-1),
   m_bin(-1),
   m_genwithpdf("")
 { 
 
   if ( obs.size()==0 ) { 
     std::cerr << "grid::not enough bins in observable" << std::endl;
     std::exit(0);
   } 
   
   //  double* obsbins = new double[obs.size()];  
   //  for ( unsigned i=0 ; i<obs.size() ; i++ ) obsbins[i] = obs[i];
   //  int Nobs = obs.size()-1;
 
   // Initialize histogram that saves the correspondence obsvalue<->obsbin
   m_obs_bins=new TH1D("referenceInternal","Bin-Info for Observable", obs.size()-1, &obs[0] );
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2(); /// grrr root is so rubbish - not scaling errors properly
   //  delete[] obsbins;
 
   m_obs_bins_combined = m_obs_bins;
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf
   //  if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
   findgenpdf( m_genpdfname );
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) m_grids[iorder] = new igrid*[obs.size()-1];
 
 }
 
 
 
 appl::grid::grid(const std::string& filename, const std::string& dirname)  :
   m_leading_order(0),  m_order(0),
   m_optimised(false),  m_trimmed(false), 
   m_normalised(false),
   m_symmetrise(false), m_transform(""), 
   m_dynamicScale(0),
   m_applyCorrections(false),
   m_documentation(""),
   m_type(STANDARD),
   m_read(false),
   m_subproc(-1),
   m_bin(-1),
   m_genwithpdf("")
 {
 
   if ( appl_first ) { 
     appl_first = false;
     std::cout << "appl::grid() version " << VERSION << " compiled " << __DATE__ << std::endl;
   }
   
   m_obs_bins_combined = m_obs_bins = 0;
 
   struct timeval tstart = appl_timer_start();
   
   struct stat _fileinfo;
   if ( stat(filename.c_str(),&_fileinfo) )   {    
     throw exception( std::string("grid::grid() cannot open file ") + filename  ); 
   }
 
   std::cout << "appl::grid() reading grid from file " << filename;
   
   TFile* gridfilep = TFile::Open(filename.c_str());
 
   if (gridfilep==0 ) {
     std::cout << std::endl;
     throw exception( std::string("grid::grid() cannot open file: ") + filename );
   }
 
   if (gridfilep->IsZombie()) {
     std::cout << std::endl;
     delete gridfilep;
     throw exception( std::string("grid::grid() cannot open file: zombie ") + filename );
   }
 
   // TFile gridfile(filename.c_str());
   
   //  gDirectory->cd(dirname.c_str());
 
   //  std::cout << "pwd=" << gDirectory->GetName() << std::endl;
 
   //  gDirectory->cd(dirname.c_str());
 
   //  std::cout << "pwd=" << gDirectory->GetName() << std::endl;
 
   //  Directory d(dirname);
   //  d.push();
 
   TFileString* _tagsp = (TFileString*)gridfilep->Get((dirname+"/Tags").c_str());  
 
   if ( _tagsp==0 ) { 
     std::cout << std::endl;
     throw exception( std::string("grid::grid() cannot get tags: ") + filename );
   }
 
   TFileString _tags = *_tagsp;
   // TFileString _tags = *(TFileString*)gridfile.Get("Tags");  
 
   m_transform  = _tags[0];
   m_genpdfname = _tags[1];
   std::string _version = _tags[2];
 
   //  std::cout << "tags:: transform: " << m_transform << "\tpdfname: " << m_genpdfname << "\tdoc: " << m_documentation << std::endl;  
 
   if ( _tags.size()>3 ) m_documentation = _tags[3];
 
   std::cout << "\t:  applgrid version " << _version;
   if ( _version != m_version ) std::cout  << " (transformed to " << m_version << ")";
 
   /// encoded pdfset used for the generation 
   if ( _tags.size()>4 ) { 
     std::string tmp = _tags[4];
     if ( tmp.find(":")!=std::string::npos ) {  
       m_genwithpdf  = tmp.substr(0,tmp.find(":"));
       m_genwithipdf = std::stoi(tmp.substr(tmp.find(":")+1,std::string::npos));
     }
     else { 
       m_genwithpdf  = tmp;
       m_genwithipdf = 0;
     }
     std::cout << "\t:  grid generated with " << m_genwithpdf << " : set " << m_genwithipdf;
   }
 
 
   // check it has the correct version
   // if ( _version > m_version ) { 
   //      throw exception(cerr << "incorrect version " << _version << " expected " << m_version ); 
   // }
   //  m_version = _version;
   
   std::cout << std::endl;
 
   if ( getDocumentation()!="" ) std::cout << getDocumentation() << std::endl; 
   
   bool added = false;
 
   if ( m_genpdfname=="basic" ) { 
     added = true;
     m_genpdfname = "basic.config"; 
     if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
     findgenpdf( m_genpdfname );
   }
 
   //  std::cout << "Tags=" << _tags << std::endl;
 
   //  std::cout << "grid::grid() use transform " << m_transform << std::endl;
 
   // read state information
   // hmmm, have to use TVectorT<double> since TVector<int> 
   // apparently has no constructor (???)
   TVectorT<double>* setup=(TVectorT<double>*)gridfilep->Get((dirname+"/State").c_str());
  
   m_run        = (*setup)(0);
   m_optimised  = ( (*setup)(1)!=0 ? true : false );
   m_symmetrise = ( (*setup)(2)!=0 ? true : false );  
 
   m_leading_order = int((*setup)(3)+0.5);  
   m_order         = int((*setup)(4)+0.5);  
 
   if ( setup->GetNoElements()>5 ) m_cmsScale = (*setup)(5);
   else                            m_cmsScale = 0;
  
   if ( setup->GetNoElements()>6 ) m_normalised = ( (*setup)(6)!=0 ? true : false );
   else                            m_normalised = true;
 
   if ( setup->GetNoElements()>7 ) m_applyCorrections = ( (*setup)(7)!=0 ? true : false );
   else                            m_applyCorrections = false;
 
   int n_userdata = 0;
   if ( setup->GetNoElements()>10 ) n_userdata = int((*setup)(10)+0.5);
 
   //  std::vector<double> _ckmsum;
   std::vector<std::vector<double> > _ckm2;
   std::vector<std::vector<double> > _ckm;
 
 
 
   /// check whether we need to read in the ckm matrices
 
   if ( setup->GetNoElements()>8 && (*setup)(8)!=0 ) {
 
     std::cout << "grid::grid() read ckm matrices" << std::endl;
     
     /// try 13x13 squared ckm matrix 
 
     TVectorT<double>* ckm2flat=(TVectorT<double>*)gridfilep->Get((dirname+"/CKM2").c_str());
 
     if ( ckm2flat ) { 
       if ( ckm2flat->GetNrows()>0 ) { 
 	//	_ckm2 = std::vector<std::vector<double> >(14, std::vector<double>(14) );
 	_ckm2 = std::vector<std::vector<double> >(13, std::vector<double>(13) );
       
 	for ( int ic=0 ; ic<13 ; ic++ ) { 
 	  for ( int id=0 ; id<13 ; id++ ) _ckm2[ic][id] = (*ckm2flat)(ic*13+id); 
 	}
 
 	//	for ( int ic=0 ; ic<14 ; ic++ ) _ckm2[13][ic] = 0; 
 	//	for ( int ic=0 ; ic<14 ; ic++ ) _ckm2[ic][13] = 0; 
 
       }  
     }
 
     /// now try usual 3x3 matrix
 
     TVectorT<double>* ckmflat=(TVectorT<double>*)gridfilep->Get((dirname+"/CKM").c_str());
 
     if ( ckmflat ) { 
       if ( ckmflat->GetNrows()>0 ) { 
 	_ckm = std::vector<std::vector<double> >(3, std::vector<double>(3) );
 	
 	for ( int ic=0 ; ic<3 ; ic++ ) { 
 	  for ( int id=0 ; id<3 ; id++ ) _ckm[ic][id] = (*ckmflat)(ic*3+id); 
 	}
       }  
     }
 
   }
 
   if ( setup->GetNoElements()>9 ) m_type = (CALCULATION)int( (*setup)(9)+0.5 );
   else                            m_type = STANDARD;
 
   //  std::cout << "appl::grid() reading grid calculation type: " << _calculation(m_type) << std::endl;
 
   //  std::cout << "appl::grid() normalised: " << getNormalised() << std::endl;
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf (or lumi_pdf for amcatnlo)
   //  if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
 
   //  std::cout << "appl::grid() requested pdf combination " << m_genpdfname << std::endl;
   
 
   if ( !added && contains(m_genpdfname, ".config") ) { 
     /// decode the pdf combination if appropriate
 
     // find out if we have one combination per order or one overall
     std::vector<std::string> namevec = parse( m_genpdfname, ":" );
 
     std::string label = dirname+"/Combinations";
 
     for ( unsigned i=0 ; i<namevec.size() ; i++ ) { 
 
       /// again have to use TVectorT<double> because TVectorT<int> has no constructor!!!
       /// I ask you!! what's the point of a template if it doesn't actually instantiate
       /// it's pathetic!
 
       TVectorT<double>* _combinations = (TVectorT<double>*)gridfilep->Get( label.c_str() );
 
       label += "N"; /// add an N for each order, N-LO, NN-LO etc
 
       if ( _combinations==0 ) throw exception( std::string( "grid::grid() cannot read pdf combination " ) + namevec[i] );
 
       std::vector<int> combinations(_combinations->GetNoElements());
 
       for ( unsigned ic=0 ; ic<combinations.size() ; ic++ )  combinations[ic] = int((*_combinations)(ic)); 
     
       //      addpdf(m_genpdfname, combinations);
       addpdf( namevec[i], combinations);
 
     }
   }
   else { 
     /// of just create the generic from the file
     if ( contains(m_genpdfname, ".dat") ) addpdf(m_genpdfname);
   }
   
   /// retrieve the pdf routine 
   findgenpdf( m_genpdfname );
 
   //  std::cout << "grid::grid() read " << m_genpdfname << " " << m_genpdf[0]->getckmsum().size() << std::endl; 
 
   // set the ckm matrices 
   if      ( _ckm.size()>0 )  setckm( _ckm );
   else if ( _ckm2.size()>0 ) setckm2( _ckm2 );
 
   delete setup;
 
   //  std::cout << "grid::grid() read setup" << std::endl;
 
   // Read observable bins information
   //  gridfile.GetObject("obs_bins", m_obs_bins );
 
   m_obs_bins = (TH1D*)gridfilep->Get((dirname+"/reference_internal").c_str());
   if ( m_obs_bins ) { 
     m_obs_bins_combined = (TH1D*)gridfilep->Get((dirname+"/reference").c_str());
     m_obs_bins_combined->SetDirectory(0);
     if ( run() ) m_obs_bins_combined->Scale(run());
     /// for ( int i=0 ; i<m_obs_bins_combined->GetNbinsX ; i++ ) m_obs_bins_combined->SetBinContent( i+1, m_obs_bins_combined->GetBinContent( i+1, run(i) ) );
   }
   else { 
     m_obs_bins = (TH1D*)gridfilep->Get((dirname+"/reference").c_str());
     m_obs_bins_combined = m_obs_bins;
   }
 
   m_obs_bins->SetDirectory(0);
   if ( run() ) m_obs_bins->Scale(run());
   /// for ( int i=0 ; i<m_obs_bins->GetNbinsX ; i++ ) m_obs_bins->SetBinContent( i+1, m_obs_bins->GetBinContent( i+1, run(i) ) );
 
   m_obs_bins->SetName("referenceInternal");
   if ( m_normalised && m_optimised ) m_read = true;
 
 
   ///  std::cout << "grid::grid() read obs bins" << std::endl;
 
   /// also calculate maximum, values for Q2 and y=ln(1/x)
   /// for hoppet initialisation
 
   double _Q2max = 0;
   double _xmin  = 2; // always bigger than 1
 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     //    std::cout << "grid::grid() iorder=" << iorder << std::endl;
     m_grids[iorder] = new igrid*[Nobs_internal()];  
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
       char name[128];  sprintf(name, (dirname+"/weight[alpha-%d][%03d]").c_str(), iorder, iobs);
       //   std::cout << "grid::grid() reading " << name << "\tiobs=" << iobs << std::endl;
 
       m_grids[iorder][iobs] = new igrid(*gridfilep, name);
       m_grids[iorder][iobs]->setparent( this ); 
 
       if (  !m_grids[iorder][iobs]->empty() ) { 
 
 	double x1 = m_grids[iorder][iobs]->x1filledmin();
 	double x2 = m_grids[iorder][iobs]->x1filledmin();	
 	double Q2 = m_grids[iorder][iobs]->Q2filledmax();
 	
 	
 	if (  Q2>_Q2max ) _Q2max = Q2;
 	if (  x1<_xmin )   _xmin = x1;
 	if (  x2<_xmin )   _xmin = x2;
       }
 
       //    _size += m_grids[iorder][iobs]->size();
       //      std::cout << "grid::grid() done" << std::endl;
     }
   }
 
   if ( _Q2max>0 )           m_Qmax = std::sqrt(_Q2max);
   if ( _xmin>0 && _xmin<2 ) m_ymax = std::log(1/_xmin);
 
   //  std::cout << "Qmax " << m_Qmax << std::endl; 
   //  std::cout << "ymax " << m_ymax << std::endl;
 
   //  d.pop();
 
   /// bin-by-bin correction labels                                       
   TFileString* correctionLabels = (TFileString*)gridfilep->Get((dirname+"/CorrectionLabels").c_str());  
   if ( correctionLabels ) { 
     for ( unsigned i=0 ; i<correctionLabels->size() ; i++ ) {
       m_correctionLabels.push_back( (*correctionLabels)[i] ); // copy the correction label
     }
   }
 
   /// bin-by-bin correction values
   TFileVector* corrections = (TFileVector*)gridfilep->Get((dirname+"/Corrections").c_str());  
   if ( corrections ) { 
     for ( unsigned i=0 ; i<corrections->size() ; i++ ) {
       m_corrections.push_back( (*corrections)[i] ); // copy the correction histograms
       m_applyCorrection.push_back(false);
     }
   }
 
 
   /// now read in vector of bins to be combined if present
   TVectorT<double>* _combine = (TVectorT<double>*)gridfilep->Get((dirname+"/CombineBins").c_str());
 
   if ( _combine!=0 ) { 
     //    std::cout << "read in " << _combine->GetNrows() << " entries in combine" << std::endl;
     m_combine = std::vector<int>(_combine->GetNrows(),0);
     for ( unsigned i=_combine->GetNrows() ; i-- ; ) m_combine[i] = int((*_combine)(i));
 
     if ( m_combine.size() ) combineReference();
   }
 
   //  std::cout << "grid::grid() read from file Nobs = " << Nobs_internal() << std::endl;
 
   //  std::cout << "read grid" << std::endl;
 
   /// read in the user data
   if ( n_userdata ) { 
     TVectorT<double>* userdata=(TVectorT<double>*)gridfilep->Get((dirname+"/UserData").c_str());
     m_userdata.clear();
     for ( int i=0 ; i<n_userdata ; i++ ) m_userdata.push_back( (*userdata)(i) );
   }
 
   gridfilep->Close();
   delete gridfilep;
 
   double tstop = appl_timer_stop( tstart );
 
   unsigned usize = size();
 
   struct timeval tstart2 =  appl_timer_start();
 
   trim();
 
   double tstop2 = appl_timer_stop( tstart2 );
  
   std::cout << "appl::grid() read grid, size ";
   if ( usize>1024*10 ) std::cout << usize/1024./1024. << " MB";
   else                 std::cout << usize/1024.      << " kB";
   std::cout << "\tin " << tstop << " ms";
 
   std::cout << "\ttrim in " << tstop2 << " ms" << std::endl;
 
   std::cout << "appl::grid::weights " << run() << std::endl;
 
 }
 
 
 appl::grid::grid(const grid& g) : 
   m_obs_bins(new TH1D(*g.m_obs_bins)), 
   m_leading_order(g.m_leading_order), m_order(g.m_order), 
   m_run(g.m_run), m_optimised(g.m_optimised), m_trimmed(g.m_trimmed), 
   m_normalised(g.m_normalised),
   m_symmetrise(g.m_symmetrise),
   m_transform(g.m_transform),
   m_genpdfname(g.m_genpdfname), 
   m_cmsScale(g.m_cmsScale),
   m_dynamicScale(g.m_dynamicScale),
   m_applyCorrections(g.m_applyCorrections), 
   m_documentation(g.m_documentation),
   m_ckmsum(g.m_ckmsum), /// need a deep copy of the contents
   m_ckm2(g.m_ckm2),     /// need a deep copy of the contents
   m_ckm(g.m_ckm),       /// need a deep copy of the contents
   m_type(g.m_type),
   m_read(g.m_read),
   m_bin(-1),
   m_genwithpdf("")
 {
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2();
 
   m_obs_bins_combined = m_obs_bins;
   if ( g.m_obs_bins_combined != g.m_obs_bins) { 
     m_obs_bins_combined = new TH1D(*g.m_obs_bins_combined); 
     m_obs_bins_combined->SetDirectory(0);
   }
 
   /// check to see if we require a generic pdf from a text file, and 
   /// and if so, create the required generic pdf
   //  if ( m_genpdfname.find(".dat")!=std::string::npos ) addpdf(m_genpdfname);
   if ( contains(m_genpdfname, ".dat") ||  contains(m_genpdfname, ".config") ) addpdf(m_genpdfname);
   findgenpdf( m_genpdfname );
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     m_grids[iorder] = new igrid*[Nobs_internal()];
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  { 
       m_grids[iorder][iobs] = new igrid(*g.m_grids[iorder][iobs]);
       m_grids[iorder][iobs]->setparent( this ); 
     }
   }
 } 
 
 
 // Initialize histogram that saves the correspondence obsvalue<->obsbin
 
 // constructor common internals 
 void appl::grid::construct(int Nobs, 
 			   int NQ2,  double Q2min, double Q2max, int Q2order,
 			   int Nx,   double xmin,  double xmax,  int xorder, 
 			   int order, 
 			   std::string transform ) { 
   
   //  std::cout << "appl::grid::construct() m_order " << m_order << "\tNobs " << Nobs << std::endl; 
 
   if ( m_order!=order ) std::cerr << "appl::grid::construct() order mismatch" << std::endl;
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) m_grids[iorder] = 0;
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     m_grids[iorder] = new igrid*[Nobs];
     for ( int iobs=0 ; iobs<Nobs ; iobs++ ) {
       m_grids[iorder][iobs] = new igrid(NQ2, Q2min, Q2max, Q2order, 
 					Nx, xmin, xmax, xorder, 
 					transform, m_genpdf[iorder]->Nproc());
       m_grids[iorder][iobs]->setparent( this ); 
     }
   }
   //  std::cout << "appl::grid::construct() return" << std::endl; 
 }
 
 
 
 
   // number of subprocesses 
 int appl::grid::subProcesses(int i) const { 
   if ( i<0 || i>=m_order ) { 
     std::stringstream s;
     s << "grid::subProcess(int i) " << i << " out of range [0-" << m_order-1 << "]";
     throw exception( s.str() );
   }
   return m_grids[i][0]->SubProcesses();     
 }  
 
 
 /// access the transform functions for the appl::igrid so that the 
 /// igrid can be hidden 
 double appl::grid::transformvar()         { return igrid::transformvar(); }
 double appl::grid::transformvar(double v) { return igrid::transformvar(v); }
 
 
 
 // add a single grid
 void appl::grid::add_igrid(int bin, int order, igrid* g) { 
 
   if ( !(order>=0 && order<m_order) ) { 
     std::cerr << "grid::add_igrid() order out of range " << order << std::endl; 
     return;
   } 
 
   if ( !(bin>=0 && bin<Nobs_internal() ) ) {
     //    std::cerr << "grid::add_igrid() observable bin out of range " << bin << std::endl; 
     return;
   }
 
   m_grids[order][bin] = g;
   m_grids[order][bin]->setparent(this);
 
   if ( g->transform()!=m_transform ) { 
     std::cerr << "grid::add_igrid() transform " << m_transform 
 	      << " does not match that from added igrid, " << g->transform() << std::endl;
     m_transform = g->transform();
   }
 
 }
 
 
 
 
 appl::grid::~grid() {
 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {  
     if( m_grids[iorder] ) { 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) { 
 	delete m_grids[iorder][iobs];
       }
       delete[] m_grids[iorder];
       m_grids[iorder] = 0;
     }
   }
 
 
   if (m_obs_bins_combined) {
     if ( m_obs_bins_combined!=m_obs_bins) delete m_obs_bins_combined;
   }
 
   if (m_obs_bins) delete m_obs_bins;
   m_obs_bins=0;
   m_obs_bins_combined = 0;
 
 #ifdef HAVE_HOPPET
   if ( hoppet ) delete hoppet; 
   hoppet=0; 
 #endif
 
 }
 
 
 
 
 // algebraic operators
 
 appl::grid& appl::grid::operator=(const appl::grid& g) { 
   // clear out the old...
   if ( m_obs_bins_combined!=m_obs_bins ) delete m_obs_bins_combined;
   delete m_obs_bins;
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  delete m_grids[iorder][iobs];
     delete m_grids[iorder];
   }
   
   // copy the new
   m_obs_bins = new TH1D(*g.m_obs_bins);
   m_obs_bins->SetDirectory(0);
   m_obs_bins->Sumw2();
 
   m_obs_bins_combined = m_obs_bins;
 
   // copy the state
   m_leading_order = g.m_leading_order;
   m_order         = g.m_order;
 
   m_run       = g.m_run;
   m_optimised = g.m_optimised;
   m_trimmed   = g.m_trimmed;
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     m_grids[iorder] = new igrid*[Nobs_internal()];
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  { 
       m_grids[iorder][iobs] = new igrid(*g.m_grids[iorder][iobs]);
       m_grids[iorder][iobs]->setparent( this ); 
     }
   }
   return *this;
 } 
   
 
 
 
 double integral( TH1D* h ) { 
   double d = 0;
   for ( int i=0 ; i<h->GetNbinsX() ; i++ ) d += h->GetBinContent(i+1);
   return d;
 }
 
 
 appl::grid& appl::grid::operator*=(const double& d) { 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) (*m_grids[iorder][iobs]) *= d; 
   }
   getReference()->Scale( d );
   if(getReference()!=getReference_internal()) getReference_internal()->Scale( d );
   combineReference(true);
 
   return *this;
 }
 
 
 void hprint( TH1D* h ) { 
   for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) std::cout << h->GetBinContent(i) << " ";
   std::cout << std::endl;
 }
 
 
 
 /// this is messed up - how can we apply "per bin" normalisations
 /// it we have specified more actual bins and are combining them 
 /// a-posteriori ? 
 /// Fixme: at some point *remove* the bin combining functionality
 ///        it is far too much of a pain to maintain
 
 appl::grid& appl::grid::operator*=(const std::vector<double>& v) {
 
   /// don't do anything if the scale vector is not the correct size
   if ( v.size() != size_t(Nobs_internal()) ) return *this;
 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) (*m_grids[iorder][iobs])*=v[iobs]; 
   }
 
 #if 0   
   for ( int i=0 ; i<getReference()->GetNbinsX() ; i++ ) {
     int ih=i+1;
     getReference()->SetBinContent( ih,  getReference()->GetBinContent(ih)*v[i] );
     getReference()->SetBinError(   ih,  getReference()->GetBinError(ih)*v[i] );
   }
 #endif
 
   for ( int i=0 ; i<getReference_internal()->GetNbinsX() ; i++ ) {
     int ih=i+1;
     //   std::cout << "v[" << i << "] = " << v[i] << std::endl;
     getReference_internal()->SetBinContent( ih,  getReference_internal()->GetBinContent(ih)*v[i] );
     getReference_internal()->SetBinError(   ih,  getReference_internal()->GetBinError(ih)*v[i] );
   }
 
   combineReference(true);
 
   return *this;
 }
 
 
 
 double appl::grid::fx(double x) const { 
   if ( m_order>0 && Nobs_internal()>0 ) return m_grids[0][0]->fx(x);
   else return 0;
 }
 
 double appl::grid::fy(double x) const { 
   if ( m_order>0 && Nobs_internal()>0 ) return m_grids[0][0]->fy(x);
   else return 0;
 }
 
 
 appl::grid& appl::grid::operator+=(const appl::grid& g) {
   m_run      += g.m_run;
   m_optimised = g.m_optimised;
   m_trimmed   = g.m_trimmed;
   if ( Nobs_internal()!=g.Nobs_internal() ) throw exception("grid::operator+ Nobs bin mismatch");
   if ( m_order!=g.m_order )                 throw exception("grid::operator+ different order grids");
   if ( m_leading_order!=g.m_leading_order ) throw exception("grid::operator+ different order processes in grids");
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) { 
       (*m_grids[iorder][iobs]) += (*g.m_grids[iorder][iobs]); 
     }
   }
 
   /// grrr use root TH1::Add() even though I don't like it. 
   //  getReference()->Add( g.getReference() );
   getReference_internal()->Add( g.getReference_internal() );
 
   combineReference(true);
 
   return *this;
 }
 
 
 
 
 /// check grids match
 bool appl::grid::operator==(const appl::grid& g) const {
   
   bool match = true; 
 
   if ( Nobs_internal()!=g.Nobs_internal() )    match = false;
   if ( m_order!=g.m_order )  match = false;
   if ( m_leading_order!=g.m_leading_order ) match = false;
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     //    for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) match &= ( (*m_grids[iorder][iobs]) == (*g.m_grids[iorder][iobs]) ); 
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) match &= ( m_grids[iorder][iobs]->compare_axes( *g.m_grids[iorder][iobs] ) ); 
   }
   return match;
 }
 
 
 
 // fill the appropriate igrid with these weights
 void appl::grid::fill(const double x1, const double x2, const double Q2, 
 		      const double obs, 
 		      const double* weight, const int iorder)  {  
   int iobs = m_obs_bins->FindBin(obs)-1;
 
   //  std::cout << "grid::fill() iobs " << iobs << "\tobs=" << obs << std::endl;
 
   if ( iobs<0 || iobs>=Nobs_internal() ) {
     //    cerr << "grid::fill() iobs out of range " << iobs << "\tobs=" << obs << std::endl;
     //    cerr << "obs=" << obs << "\tobsmin=" << obsmin() << "\tobsmax=" << obsmax() << std::endl;
     return;
   }
   
   //  std::cout << "iobs=" << iobs << "\tobs=" << obs;
   //  for ( int i=0 ; i<subProcesses(iorder) ; i++ ) std::cout << "\t" << weight[i];
   //  std::cout << std::endl;
 
   //  std::cout << "\tiobs=" << iobs << std::endl;
   if ( m_symmetrise && x2<x1 )  m_grids[iorder][iobs]->fill(x2, x1, Q2, weight);
   else                          m_grids[iorder][iobs]->fill(x1, x2, Q2, weight);
 }
 
 
 // fast fill pre-optimisation don't perform the interpolation and so on
 void appl::grid::fill_phasespace(const double x1, const double x2, const double Q2, 
 				 const double obs, 
 				 const double* weight, const int iorder) {
   int iobs = m_obs_bins->FindBin(obs)-1;
 
   //  std::cout << "grid::fill_phasespace() iobs " << iobs << "\tobs=" << obs << std::endl;
 
   if ( iobs<0 || iobs>=Nobs_internal() ) {
     // std::cerr << "grid::fill() iobs out of range " << iobs << "\tobs=" << obs << std::endl;
     // std::cerr << "obs=" << obs << "\tobsmin=" << obsmin() << "\tobsmax=" << obsmax() << std::endl;
     return;
   }
   if ( m_symmetrise && x2<x1 )  m_grids[iorder][iobs]->fill_phasespace(x2, x1, Q2, weight);
   else                          m_grids[iorder][iobs]->fill_phasespace(x1, x2, Q2, weight);
 }
 
 
 
 
 // fast fill pre-optimisation don't perform the interpolation and so on
 void appl::grid::fill_index(const int ix1, const int ix2, const int iQ2, 
 			    const int iobs, 
 			    const double* weight, const int iorder) {
 
   if ( iobs<0 || iobs>=Nobs_internal() ) {
     //  cerr << "grid::fill() iobs out of range " << iobs << "\tobs=" << obs << std::endl;
     //  cerr << "obs=" << obs << "\tobsmin=" << obsmin() << "\tobsmax=" << obsmax() << std::endl;
     return;
   }
   if ( m_symmetrise && ix2<ix1 )  m_grids[iorder][iobs]->fill_index(ix2, ix1, iQ2, weight);
   else                            m_grids[iorder][iobs]->fill_index(ix1, ix2, iQ2, weight);
 }
 
 
 void appl::grid::trim(int iorder) {
   if ( iorder>=0 ) { 
     if ( iorder<m_order ) for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) m_grids[iorder][iobs]->trim(); 
   }
   else { 
     m_trimmed = true;
     for( int iord=0 ; iord<m_order ; iord++ ) {
       for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) m_grids[iord][iobs]->trim(); 
     }
   }
 }
 
 void appl::grid::untrim(int iorder) {
   if ( iorder>=0 ) { 
     if ( iorder<m_order ) for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) m_grids[iorder][iobs]->untrim(); 
   }
   else { 
     m_trimmed = false;
     for( int iord=0 ; iord<m_order ; iord++ ) {
       for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) m_grids[iord][iobs]->untrim(); 
     }
   }
 }
 
 std::ostream& appl::grid::print(std::ostream& s) const {
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {     
       s << iobs << "\t" 
 	<< std::setprecision(5) << std::setw(6) << getReference()->GetBinLowEdge(iobs+1) << "\t- " 
 	<< std::setprecision(5) << std::setw(6) << getReference()->GetBinLowEdge(iobs+2) << "\t"; 
       m_grids[iorder][iobs]->print(s);       
     }
   }
   return s;
 }
 
 
 
 
  /// get the required pdf combinations from those registered   
 void appl::grid::findgenpdf( std::string s ) { 
     std::vector<std::string> names = parse( s, ":" );
     if ( names.size()==unsigned(m_order) ) for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i] = appl_pdf::getpdf( names[i] );
     else  if ( names.size()==1 )           for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i] = appl_pdf::getpdf( names[0] );
     else  { 
       std::stringstream s_;
       s_ << "requested " << m_order << " pdf combination but given " << names.size();
       throw exception( s_.str() );
     }
 }
 
 
 void appl::grid::addpdf( const std::string& s, const std::vector<int>& combinations ) {
 
   //  std::cout << "addpdf() in " << std::endl;
 
     /// parse names, if they contain .dat, then create the new generic pdfs
     /// they will be added to the pdf std::map automatically 
     std::vector<std::string> names = parse( s, ":" );
 
     unsigned imax = unsigned(m_order); 
 
     /// check to see whether we have a different pdf for each order
     if ( names.size()!=imax ) { 
       if ( names.size()==1 ) imax = 1;
       else { 
 	std::stringstream s_;
 	s_ << "requested " << m_order << " pdf combination but given " << names.size();
 	throw exception( s_.str() );
       }
     }
 
     //    std::cout << "imax " << imax << std::endl; 
 
     /// loop through all the required pdfs checking whether they exist already,
     /// if not (from thrown exception) then create it, otherwise, don't need to 
     /// do anything 
     for ( unsigned i=0 ; i<imax ; i++ ) { 
 
       //      std::cout << "\ti " << i<< std::endl; 
 
       if ( names[i].find(".dat")!=std::string::npos ) {
 
 	if ( appl_pdf::getpdf(names[i])==0 ) { 
  	  std::cout << "appl::grid::addpdf() creating new generic_pdf " << names[i] << std::endl;
 	  new generic_pdf(names[i]);
 	}
 
 	// 	try {
 	// 	  appl_pdf::getpdf(names[i]); // , false);
 	// 	}
 	// 	catch ( appl_pdf::exception e ) { 
 	// 	  std::cout << "creating new generic_pdf " << names[i] << std::endl;
 	// 	  new generic_pdf(names[i]);
 	// 	}
 
       }
       else if ( names[i].find(".config")!=std::string::npos ) { 
 
 	if ( appl_pdf::getpdf(names[i])==0 ) { 
 	  std::cout << "appl::grid::addpdf() creating new lumi_pdf " << names[i] << std::endl;
 	  lumi_pdf* lp = new lumi_pdf(names[i], combinations);
 	  //  latex( *lp );
 	  lp->removeDuplicates();
 	  //	  latex( *lp, "duff" );
 	}
 
 	// 	try {
 	// 	  appl_pdf::getpdf(names[i]); // , false);
 	// 	}
 	// 	catch ( appl_pdf::exception e ) { 
 	// 	  std::cout << "creating new lumi_pdf " << names[i] << std::endl;
 	// 	  new lumi_pdf(names[i], combinations);
 	// 	  //	  std::cout << "created" << names[i] << std::endl;
 	
 	// 	}
 
       }
 
     }
 
 }
 
 
 
 void appl::grid::setckm2( const std::vector<std::vector<double> >& ckm2 ) { 
   for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i]->setckm2(ckm2);
 }
 
 
 void appl::grid::setckm( const std::vector<std::vector<double> >& ckm ) { 
   for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i]->setckm(ckm);
 }
 
 
 void appl::grid::setckm( const std::vector<double>& ckm ) {
   std::vector<std::vector<double> > _ckm(3, std::vector<double>(3,0) );
   for ( unsigned i=0 ; i<ckm.size() && i<9 ; i++ ) _ckm[i/3][i%3] = ckm[i]; 
   for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i]->setckm(_ckm);
 }
 
 void appl::grid::setckm( const double* ckm ) { 
   std::vector<std::vector<double> > _ckm(3, std::vector<double>(3,0) );
   for ( unsigned i=0 ; i<9 ; i++ ) _ckm[i/3][i%3] = ckm[i]; 
   for ( int i=0 ; i<m_order ; i++ ) m_genpdf[i]->setckm(_ckm);
 }
 
 
 const std::vector<std::vector<double> >& appl::grid::getckm()  const { return m_genpdf[0]->getckm(); }  
 
 const std::vector<std::vector<double> >& appl::grid::getckm2() const { return m_genpdf[0]->getckm2(); }  
 
 
 // void appl::grid::setuppdf(void (*pdf)(const double&, const double&, double* ) )  {  }
 // void grid::pdfinterp(double x, double Q2, double* f) {  }
 
 
 
 // set the rewight flag of the internal grids
 bool appl::grid::reweight(bool t) { 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) m_grids[iorder][iobs]->reweight(t);       
   }
   return t;
 }
 
 
 bool exists( const std::string& filename ) { 
   struct stat sb;
   if ( stat( filename.c_str(), &sb)==0 ) return true; // && S_ISREG(sb.st_mode ))
   else return false;
 }
 
 
 // dump to file
 void appl::grid::Write(const std::string& filename, 
 		       const std::string& dirname, 
 		       const std::string& pdfname) 
 { 
 
   std::cout << "appl::grid::Write() " << filename << "\tdirname " << dirname << "\tpdfname " << pdfname << std::endl; 
 
   if ( exists( filename ) ) { 
     std::string filename_save = filename + "-save";
     if ( std::rename( filename.c_str(), filename_save.c_str() ) ) std::cerr << "could not rename grid file " << filename << std::endl;
   } 
 
   if ( pdfname!="" ) shrink( pdfname, m_genpdf[0]->getckmcharge() );
 
   //  std::cout << "grid::Write() writing to file " << filename << std::endl;
   TFile rootfile(filename.c_str(),"recreate");
 
   //  std::cout << "pwd=" << gDirectory->GetName() << std::endl;
 
   Directory d(dirname);
   d.push();
   
   //  std::cout << "pwd=" << gDirectory->GetName() << std::endl;
 
   // write the name of the transform pair and the
   // generalised pdf
   TFileString _tags("Tags");
   _tags.add(m_transform);
   _tags.add(m_genpdfname);
   _tags.add(m_version);
   if ( m_documentation!="" ) _tags.add( m_documentation );
   if ( m_genwithpdf!="" )    _tags.add( m_genwithpdf + ":"+std::to_string(m_genwithipdf) );
   _tags.Write();
 
   std::cout << "genwith pdf " << m_genwithpdf << " : " << m_genwithipdf << std::endl; 
 
   //  TH1D* _transform = new TH1D("Transform", m_transform.c_str(), 1, 0, 1);  
   //  _transform->Write();
 
   //  TH1D* _genpdfname = new TH1D("GenPDF", m_genpdfname.c_str(), 1, 0, 1);
   //  _genpdfname->Write();
 
 
   //  std::cout << "state std::vector=" << std::endl;
 
   // state information
   TVectorT<double>* setup=new TVectorT<double>(12); // add a few extra just in case 
   (*setup)(0) = m_run;
   (*setup)(1) = ( m_optimised  ? 1 : 0 );
   (*setup)(2) = ( m_symmetrise ? 1 : 0 );
   (*setup)(3) =   m_leading_order ;
   (*setup)(4) =   m_order ;
   (*setup)(5) =   m_cmsScale ;
   (*setup)(6) = ( m_normalised ? 1 : 0 );
   (*setup)(7) = ( m_applyCorrections ? 1 : 0 );
 
   if ( m_genpdf[0]->getckmsum().size()==0 ) (*setup)(8) = 0;
   else                                      (*setup)(8) = 1;
 
   (*setup)(9) = (int)m_type;
 
   (*setup)(10) = m_userdata.size();
 
   setup->Write("State");
   
   if ( (*setup)(8) == 1 ) { 
     
     /// no longer write out squared ckm matrix - just use the 3x3
     TVectorT<double>* ckmflat = new TVectorT<double>(9);
     const std::vector<std::vector<double> >& _ckm = m_genpdf[0]->getckm();
 
     for ( int ic=0 ; ic<3 ; ic++ ) { 
       for ( int id=0 ; id<3 ; id++ ) (*ckmflat)(ic*3+id) = _ckm[ic][id];
     }
 
     ckmflat->Write("CKM");
 
   }
 
   /// encode the pdf combination if appropriate
   
   if ( contains( m_genpdfname, ".config" ) ) { 
     
     // find out if we have one combination per order or one overall
     
     std::vector<std::string> namevec = parse( m_genpdfname, ":" );
     
     /// now write out all the ones we need
     
     std::string label = "Combinations";
     
     for ( unsigned i=0 ; i<namevec.size() && i<unsigned(m_order) ; i++ ) {  
 
       //      dynamic_cast<lumi_pdf*>(m_genpdf[i])->restoreDuplicates();
 
       lumi_pdf* lpdf = dynamic_cast<lumi_pdf*>(m_genpdf[i]);
 
       std::cout << "lumi_pdf:   " << lpdf->name() << "\t processes " << lpdf->Nproc() << "\n";
 
       //      std::cout << "lumi pdf: " << *lpdf << std::endl;
 
       lumi_pdf lpdf_tmp(*lpdf);
 
       lpdf_tmp.restoreDuplicates();
 
       //      std::cout << "lumi pdf: " << *lpdf     << std::endl;
       //      std::cout << "lumi pdf: " <<  lpdf_tmp << std::endl;
 
       //      std::vector<int>   combinations = dynamic_cast<lumi_pdf*>(m_genpdf[i])->serialise();
       std::vector<int>   combinations = lpdf_tmp.serialise();
       TVectorT<double>* _combinations = new TVectorT<double>(combinations.size());
       for ( unsigned ic=0 ; ic<combinations.size() ; ic++ ) { 
 	if ( combinations[ic]<0 ) (*_combinations)(ic) = double(combinations[ic]-0.5);
 	else                      (*_combinations)(ic) = double(combinations[ic]+0.5);
       }
       
       _combinations->Write( label.c_str() );
       
       std::cout << "writing " << m_genpdf[i]->name() << std::endl;
       
       label += "N";  /// add an N for each order, N-LO, NN-LO etc !!!! ARGH!!! Need to be more general !!!
     }
   } 
   
   //  std::cout << "grids Nobs = " << Nobs_internal() << std::endl;
   
   untrim();
   int untrim_size = size();
   trim();
   int trim_size = size();
   std::cout <<"grid::Write()"
 	    << "size(untrimmed)=" << untrim_size/1024 << " kB"
 	    << "\tsize(trimmed)=" <<   trim_size/1024 << " kB"
 	    << " (" << 0.1*int(trim_size*1000./untrim_size) << "%)" << std::endl;
   
 
   // internal grids
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
       char name[128];  sprintf(name, "weight[alpha-%d][%03d]", iorder, iobs);
       // std::cout << "writing grid " << name << std::endl;
       //   _size += m_grids[iorder][iobs]->size();
       m_grids[iorder][iobs]->write(name);
       //   trim_size += m_grids[iorder][iobs]->size();
     }
   }
  
   
   //  d.pop();
 
   //  std::cout << "reference" << std::endl;
 
   TH1D* reference          = 0;
   TH1D* reference_internal = 0;
 
 
   if ( m_obs_bins_combined != m_obs_bins )  { 
     reference = (TH1D*)m_obs_bins_combined->Clone("reference");
     reference->SetDirectory(0);
 
     reference_internal = (TH1D*)m_obs_bins->Clone("reference_internal");
     reference_internal->SetDirectory(0);
   }
   else { 
     reference = (TH1D*)m_obs_bins->Clone("reference");
     reference->SetDirectory(0);
   }
 
   /// grrr, what is this? This whole normalisation issue is a real pain - 
   /// we need to address this - but how to do it in a backwards compatible way?
   /// will need to check the grid version and do things differently probably
 
   //  std::cout << "normalised() " << getNormalised() << "\tread " << m_read << std::endl; 
   
   //  if ( !getNormalised() || m_read )  if ( run() ) reference->Scale(1/double(run()));
   if ( run() ) { 
     reference->Scale(1/double(run()));
     if ( reference_internal ) reference_internal->Scale(1/double(run()));
   }
 
   // if ( run() ) reference->Scale(1/double(run()));
 
 #if 0
 
   double hmin = reference->GetBinContent(1); 
   double hmax = reference->GetBinContent(1); 
   for ( int i=1 ; i<=reference->GetNbinsX() ; i++ ) { 
     std::cout << "reference " << i << " " << reference->GetBinContent(i) << std::endl;
     if ( hmin<reference->GetBinContent(i) ) hmin = reference->GetBinContent(i); 
     if ( hmax>reference->GetBinContent(i) ) hmax = reference->GetBinContent(i); 
   }
 
   reference->SetMinimum(hmin);
   reference->SetMaximum(hmax);
 
 #endif
 
 
 
 
   reference->Write();
 
   delete reference;
   if ( reference_internal ) { 
     reference_internal->Write();
     delete reference_internal;
   }
 
   //  std::cout << "corrections" << std::endl;
 
   /// correction histograms
 
   if ( m_corrections.size()>0 ) {
 
     /// Fixme: should add the labels to the actual corrections rather than save separately
     /// write labels
     TFileVector* corrections = new TFileVector("Corrections");
     for ( unsigned i=0 ; i<m_corrections.size() ; i++ )  corrections->add( m_corrections[i] );    
     corrections->Write("Corrections");
 
     /// write actual corrections
     TFileString correctionLabels("CorrectionLabels");
     for ( unsigned i=0 ; i<m_correctionLabels.size() ; i++ )  correctionLabels.add( m_correctionLabels[i] );
     correctionLabels.Write("CorrectionLabels");
 
   }
 
   /// now write out vector of bins to be combined if this has been set
   if ( m_combine.size()>0 ) { 
     TVectorT<double>* _combine = new TVectorT<double>(m_combine.size()); 
     for ( unsigned i=m_combine.size() ; i-- ; ) (*_combine)(i) = m_combine[i]+0.5; /// NB: add 0.5 to prevent root double -> int rounding errors
     _combine->Write( "CombineBins" );
   }
 
 
   /// now write out the user data
 
   if ( m_userdata.size() ) { 
     TVectorT<double>* userdata=new TVectorT<double>(m_userdata.size()); // add a few extra just in case 
     for ( unsigned i=0 ; i<m_userdata.size() ; i++ ) (*userdata)(i) = m_userdata[i];
     userdata->Write("UserData");
   }
 
   d.pop();
   rootfile.Close();
 
 }
 
 
 void appl::grid::disable_threads(bool b) { appl::igrid::disable_threads(b); }
 
 
 bool duff() { 
   appl::grid::disable_threads(true);
   return true;
 }
 
 
 // takes pdf as the pdf lib wrapper for the pdf set for the convolution.
 // type specifies which sort of partons should be included:
 
 std::vector<double> appl::grid::vconvolute(void (*pdf)(const double& , const double&, double* ), 
 					   double (*alphas)(const double& ), 
 					   int     _nloops, 
 					   double  rscale_factor,
 					   double  fscale_factor,
 					   double Escale )
 { 
   return vconvolute( pdf, 0, alphas, _nloops, rscale_factor, fscale_factor, Escale );
 }
 
 
 
 std::vector<double> appl::grid::vconvolute(void (*pdf1)(const double& , const double&, double* ), 
 					   void (*pdf2)(const double& , const double&, double* ), 
 					   double (*alphas)(const double& ), 
 					   int     _nloops, 
 					   double  rscale_factor,
 					   double  fscale_factor,
 					   double Escale )
 { 
-
+ 
   NodeCache cache1( pdf1 );
   NodeCache cache2;
 
  
   cache1.reset();
 
 
   NodeCache* _pdf1 = &cache1;
   NodeCache* _pdf2 = 0;
   
   if ( pdf2!=0 && pdf1!=pdf2 ) { 
     cache2 = NodeCache( pdf2 );
     cache2.reset();
   
     _pdf2    = &cache2;
   }
 
 
   //  struct timeval _ctimer = appl_timer_start();
   
   double Escale2 = 1;
  
   if ( Escale!=1 ) Escale2 = Escale*Escale;
   
   std::vector<double> hvec;
   hvec.clear(); /// explicit clear just in case
 
   double invNruns = 1;
   if ( (!m_normalised) && run() ) invNruns /= double(run());
 
   //  TH1D* h = new TH1D(*m_obs_bins);
   //  h->SetName("xsec");
 
 
 #ifdef HAVE_HOPPET
   // check if we need to use the splitting function, and if so see if we 
   // need to initialise it again, and do so if required
 
 
 
 #if 0
   for ( int ig=0 ; ig<m_order ; ig++ ) { 
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  { 	
       if ( m_grids[ig][iobs] ) { 
 	std::cout << "\tgrid limits " << iobs 
 		  << "\t x1min " << m_grids[ig][iobs]->getx1min() << " (" << std::log(1/m_grids[ig][iobs]->getx1min()) << ")"
 		  << "\t x1max " << m_grids[ig][iobs]->getx1max() << " (" << std::log(1/m_grids[ig][iobs]->getx1max()) << ")"
 		  << "\t x2min " << m_grids[ig][iobs]->getx2min() << " (" << std::log(1/m_grids[ig][iobs]->getx2min()) << ")"
 		  << "\t x2max " << m_grids[ig][iobs]->getx2max() << " (" << std::log(1/m_grids[ig][iobs]->getx2max()) << ")"
 		  << "\t Q2min " << m_grids[ig][iobs]->getQ2min() 
 		  << "\t Q2max " << m_grids[ig][iobs]->getQ2max() << std::endl;
       }
     }    
   }
 #endif
 
 
 
   if ( fscale_factor!=1 || m_dynamicScale ) {
 
     if ( pdf2==0 || pdf1==pdf2 ) { 
       
       /// fixed and shared between 
       /// all instances of the grid
       static double Qmax = 26000;  /// hmmm is this correct ???
       static double ymax = 12; 
 
       //      double _xmin = xmin();
       // double _ymax = std::log(1/_xmin);
 
       double _ymax = ymax;
 
       bool restart_hoppet = false;
 
       /// Qmax limit needs to be recalulated
       if ( m_cmsScale>Qmax ) { 
 	if ( hoppet ) restart_hoppet = true;
 	Qmax = m_cmsScale;
       }
 
       /// ymax limit needs to be recalculated
       if ( _ymax>ymax ) { 
 	if ( hoppet ) restart_hoppet = true;
 
 	double dymax = (_ymax - ymax);
 
 	/// NB: this 0.1 is just dy from hoppet init - would be better to caluclate the values and set this in hoppet init
 	double iy = dymax/0.1; 
 
 	std::cout << "appl::grid changing hoppet yrange " << ymax << "\t" << iy << "\t" << (0.1+int(iy)*0.1) << std::endl;
 
 	/// make sure that we are on a dy integer multiple
 	if ( int(iy)!=iy )  ymax += 0.1+int(iy)*0.1;
 	else                ymax +=     int(iy)*0.1;
       }
 
       if ( hoppet && restart_hoppet ) { 
 	std::cout << "appl::grid deleting old hoppet" << std::endl;
 	delete hoppet;
 	hoppet = 0;
       }  
 
       if ( hoppet == 0 ) {  
 	std::cout << "appl::grid initialising hoppet" << std::endl;
 	//	double Qmax = 15000;
 	//	double ymax
 	hoppet = new appl::hoppet_init( Qmax, ymax );	
       } 
 
       bool newpdf = hoppet->compareCache( pdf1 );
       
       if ( newpdf ) hoppet->fillCache( pdf1 );
 
     }
 
   }
 #endif
 
 
   std::string label;
 
   if ( _nloops>=m_order ) { 
     std::cerr << "too many loops for grid: " << _nloops << " > " << nloops() << " " << std::endl;   
     return hvec;
   } 
   
   /// NB: in the following code, the return values is stored in a class variable that 
   ///     must be accessed when each convolution has finished since the convolution 
   ///     itself might not be happening in this thread
 
 
   if ( m_type==STANDARD ) { 
 
     static bool first = true;
  
     if ( first && m_dynamicScale ) std::cout << "** grid::vconvolute() emulating dynamic scale **" << std::endl;
 
     //    std::cout << "standard convolution" << std::endl;
 
     if      ( _nloops==0  ) label = "lo        ";
     else if ( _nloops==1  ) label = "nlo       ";
     else if ( _nloops==2  ) label = "nnlo      ";
     else if ( _nloops==-1 ) label = "nlo only  ";
     else if ( _nloops==-2 ) label = "nnlo only ";
     
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {  
 
       /// here we see whether we need to emulate a dynamic scale by simply 
       /// changing the renormalisation and factorisation scale terms to give 
       /// what a dynamic scale would be in this bin
  
       if ( m_bin!=-1 && iobs!=m_bin ) continue;
 
       double dynamic_factor = 1;
 
 
       if ( m_dynamicScale ) {
 	double var = m_obs_bins->GetBinCenter(iobs+1);
 	dynamic_factor = var/m_dynamicScale;
 	//	if ( first ) std::cout << "grid::vconvolute() bin " << iobs << "\tscale " << var << "\tdynamicScale " << m_dynamicScale << "\t scale factor " << dynamic_factor << std::endl;
       } 
 
       /// now do the convolution proper
 
       //      std::cout << iobs << "\tnloops " << _nloops << std::endl;
 
 
       if ( _nloops==0 ) {
 	/// leading order cross section
 	if ( subproc()==-1 ) {  
 	  m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	}
 	else {
 
  
 	  /// fixme: for the subproceses, this is technically incorrect - the "LO" contribution 
 	  ///        includes the scale dependent "NLO" terms that are proportional to the LO 
 	  ///        coefficient functions - it could use the strict "LO" part without the scale 
 	  ///        dependent parts using order "0" rather than order "1" but see the comment 
 	  ///        for the "NLO only" contribution. The reason for this is that the subprocesses
 	  ///        can be different for LO and NLO, so this NLO part with "LO coefficients", 
 	  ///        can have different subprocesses from the actual NLO part, so the correct 
 	  ///        LO/NLO separation is only guaranteed for the full convolution, and not by 
 	  ///        subprocess
 	  m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 1, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	}
 
       }
       else if ( _nloops==1 ) {
 
 
 	// next to leading order cross section
 	// std::cout << "convolute() nloop=1" << std::endl;
 	// leading order contribution and scale dependent born dependent terms
 	m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 1, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	// std::cout << "dsigma_lo=" << dsigma_lo << std::endl;
 	// next to leading order contribution
 	//      double dsigma_nlo = m_grids[1][iobs]->convolute(pdf, m_genpdf, alphas, m_leading_order+1, 0);
 	// GPS: the NLO piece must use the same rscale_factor and fscale_factor as
 	//      the LO piece -- that's the convention that defines how NLO calculations
 	//      are done.
 	m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	// std::cout << "dsigma_nlo=" << dsigma_nlo << std::endl;
 
 
       }
       else if ( _nloops==-1 ) {
 	label = "nlo only";
 
-
 	// nlo contribution only 
 	if ( subproc()==-1 ) { 
 	  m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, -1, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	  m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1,  0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	}
 	else { 
 	  /// fixme: this is technically incorrect - the "LO" component contains the 
 	  ///        scale dependent NLO contribution dependent on the LO coefficient
 	  ///        functions. This part is difficult to include for individual 
 	  ///        subprocesses, since different subprocesses can be present 
 	  ///        at LO and NLO, so speifying subprocess X at NLO does not 
 	  ///        neccessarily correspond to subprocess X at LO, so adding the 
 	  ///        subprocesses - so these terms are only strict LO And NLO when 
 	  ///        *not* specifying subprocess
 	  m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1,  0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale);
 	}
 
       } 
       else if ( _nloops==2 ) {
 
 
 
 	// FIXME: not implemented completely yet - no scale variation
 	//      return hvec;
 	label = "nnlo    ";
 	// next to next to leading order contribution 
 	// NB: NO scale dependendent parts so only  muR=muF=mu
-	m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 2);
+	m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 2,   dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale  );
 	// next to leading order contribution      
-	m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 1);
+	m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 1, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale  );
 	// next to next to leading order contribution
-	m_grids[2][iobs]->convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+2, 0);
+	m_grids[2][iobs]->convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+2, 0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale  );
       }
       else if ( _nloops==-2 ) {
 	// next to next to leading order contribution
-	m_grids[2][iobs]->convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+2, 0);
+	m_grids[2][iobs]->convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+2, 0, dynamic_factor*rscale_factor, dynamic_factor*fscale_factor, Escale  );
       }
       else { 
 	std::stringstream s_;
 	s_ << "invalid value for nloops " << _nloops; 
 	throw grid::exception( s_.str() );
       }
 
       //      double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
       //      hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 
     }
 
 
     /// wait on convolution results and combine the values from the different igrids
    
     //    std::cout << "vconvolute:: nloops " << _nloops << std::endl; 
 
     if ( _nloops==0 ) { 
       /// LO only 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma  = 0;
 	if ( m_grids[0][iobs]->ready() ) dsigma  = m_grids[0][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 	//      printf("iobs %d   xsec raw: %lf   %lf \n",  iobs, dsigma, invNruns*Escale2*dsigma/deltaobs );     
       }    
     }
     else if ( _nloops==1 || ( _nloops==-1 && subproc()==-1 ) ) { 
       /// NLO only or overall NLO only part
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigLO  = 0;
 	double dsigNLO = 0;
 	if ( m_grids[0][iobs]->ready() ) { 
 	  dsigLO  = m_grids[0][iobs]->xsec(); 
 	  dsigNLO = m_grids[0][iobs]->xsecNLO(); 
 	} 
 	if ( m_grids[1][iobs]->ready() ) dsigNLO += m_grids[1][iobs]->xsec(); 
 	double dsigma = dsigLO + dsigNLO;
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 	//      printf("iobs %d   xsec raw: %lf   %lf \n",  iobs, dsigma, invNruns*Escale2*dsigma/deltaobs );     
       }    
     }
     else if ( _nloops==-1 && subproc()!=-1 ) { 
       /// NLO only suboprocess only
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[1][iobs]->ready() ) dsigma = m_grids[1][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
     else if ( _nloops==2 ) { 
       /// NNLO only suboprocess only
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigLO  = 0;
 	double dsigNLO = 0;
 	double dsigNNLO = 0;
 	if ( m_grids[0][iobs]->ready() ) dsigLO   = m_grids[0][iobs]->xsec(); 
 	if ( m_grids[1][iobs]->ready() ) dsigNLO  = m_grids[1][iobs]->xsec(); 
 	if ( m_grids[2][iobs]->ready() ) dsigNNLO = m_grids[2][iobs]->xsec(); 
 	double dsigma = dsigLO + dsigNLO  + dsigNNLO;
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 	//      printf("iobs %d   xsec raw: %lf   %lf \n",  iobs, dsigma, invNruns*Escale2*dsigma/deltaobs );     
       }    
     }
     else if ( _nloops==-2 ) {
       /// NNLO only
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[2][iobs]->ready() ) dsigma = m_grids[2][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
 
     first = false;
     
 
   }
   else if ( m_type==AMCATNLO ) {  
 
     //    std::cout << "amc@NLO convolution " << _nloops << std::endl;
 
     if      ( _nloops==0  )  label = "lo      ";
     else if ( _nloops==1  )  label = "nlo     ";
     else if ( _nloops==-1 )  label = "nlo only";
     else if ( _nloops==-2 )  label = "nlo_w0  ";
     else if ( _nloops==-3 )  label = "nlo_wR  ";
     else if ( _nloops==-4 )  label = "nlo_wF  ";
     
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {  
 
       if ( _nloops==0 ) {
 	  /// this is the amcatnlo LO calculation (without FKS shower)
 	  /// work out how to call from the igrid - maybe just implement additional 
 	  m_grids[3][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[3], alphas, m_leading_order,   0, rscale_factor, fscale_factor,  Escale );
        }
       else if ( _nloops==1 || _nloops==-1 ) {
 	  /// this is the amcatnlo NLO calculation (without FKS shower)
 	  /// Next-to-leading order contribution
 
 	  // Scale independent contribution
 	  m_grids[0][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order+1, 0,  rscale_factor, fscale_factor,  Escale );
 
 	  // Renormalization scale dependent contribution
 	  if ( rscale_factor!=1 ) { 
 	    m_grids[1][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 0,  rscale_factor, fscale_factor,  Escale );
 	  }
 
 	  // Factorization scale dependent contribution
 	  if ( fscale_factor!=1 ) { 
 	    m_grids[2][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+1, 0,  rscale_factor, fscale_factor,  Escale );
 	  }
       
 	  /// Add the LO contribution if we want full NLO 
 	  /// rather than specific NLO contribution only  
 	  if ( _nloops==1 ) { 
 	    /// this is the amcatnlo NLO calculation (without FKS shower)
 	    /// work out how to call from the igrid - maybe just implement additional 
 	    /// convolution routines and call them here
 	    m_grids[3][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[3], alphas, m_leading_order,   0,  rscale_factor, fscale_factor,  Escale );	
 	  }
       }
       else if ( _nloops==-2 ) { 
         /// Only the convolution from the W0 grid
 	m_grids[0][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order+1, 0,  rscale_factor, fscale_factor,  Escale );
       }
       else if ( _nloops==-3 ) {
 	/// Only the convolution from the WR grid
 	m_grids[1][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 0, rscale_factor, fscale_factor,  Escale );
       }
       else if ( _nloops==-4 ) { 
         /// Only the convolution from the WF grid
 	m_grids[2][iobs]->amc_convolute( _pdf1, _pdf2, m_genpdf[2], alphas, m_leading_order+1, 0,  rscale_factor, fscale_factor,  Escale );
       }
       else { 
 	std::stringstream s_;
 	s_ << "invalid value for nloops " << _nloops; 
 	throw grid::exception( s_.str() );
       }
 
       //  double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
       //  hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
     }
 
     /// wait on convolutions
 
     if ( _nloops==0 ) { 
       /// LO only 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma  = 0;
 	if ( m_grids[3][iobs]->ready() ) dsigma  = m_grids[3][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 	//      printf("iobs %d   xsec raw: %lf   %lf \n",  iobs, dsigma, invNruns*Escale2*dsigma/deltaobs );     
       }    
     }
     else if ( _nloops==1 || _nloops==-1 ) { 
       /// NLO only or overall NLO only part
 	
       double logF2 = 0;
       double logR2 = 0;
    
       if ( rscale_factor!=1 ) logR2 = std::log(rscale_factor*rscale_factor);
       if ( fscale_factor!=1 ) logF2 = std::log(fscale_factor*fscale_factor);
      
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 
 	double dsig0 = 0;
 	double dsigF = 0;
 	double dsigR = 0;
 	double dsigB = 0;
 
 	if ( m_grids[0][iobs]->ready() ) dsig0  = m_grids[0][iobs]->xsec(); 
 	double dsigma = dsig0;
 
 	if ( _nloops==1 ) { 
 	  if ( m_grids[3][iobs]->ready() ) dsigB  = m_grids[3][iobs]->xsec(); 
 	  dsigma += dsigB;
 	}
 	
 	if ( rscale_factor!=1 ) { 
 	  if ( m_grids[1][iobs]->ready() ) dsigR  = m_grids[1][iobs]->xsec(); 
 	  dsigma += dsigR*logR2;
 	}
 
 	if ( fscale_factor!=1 ) { 
 	  if ( m_grids[2][iobs]->ready() ) dsigF  = m_grids[2][iobs]->xsec(); 
 	  dsigma += dsigF*logF2;
 	}
 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
 	//      printf("iobs %d   xsec raw: %lf   %lf \n",  iobs, dsigma, invNruns*Escale2*dsigma/deltaobs );     
       }         
     }
     else if ( _nloops==-2) { 
       /// NLO only suboprocess only
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[0][iobs]->ready() ) dsigma = m_grids[0][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
     else if ( _nloops==-3) { 
       /// NLO only suboprocess only
 
       double logR2 = 0;
       if ( rscale_factor!=1 ) logR2 = std::log(rscale_factor*rscale_factor);
 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[1][iobs]->ready() ) dsigma = m_grids[1][iobs]->xsec(); 
 	dsigma  *= logR2;
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
     else if ( _nloops==-4) { 
       /// NLO only suboprocess only
 
       double logF2 = 0;
       if ( fscale_factor!=1 ) logF2 = std::log(fscale_factor*fscale_factor);
 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[2][iobs]->ready() ) dsigma = m_grids[2][iobs]->xsec(); 
 	dsigma  *= logF2;
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
 
   }
   else if ( m_type == SHERPA ) { 
     
     //    std::cout << "sherpa convolution" << std::endl;
 
     if      ( _nloops==0  ) label = "lo      ";
     else if ( _nloops==1  ) label = "nlo     ";
     else if ( _nloops==-1 ) label = "nlo only";
 
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {  
     
       if ( _nloops==0 ) {
 	label = "lo      ";
 	// leading order cross section
 	m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order, 0, 1, 1, Escale);
       }
       else if ( _nloops==1 ) { 
 	label = "nlo     ";
 	// next to leading order cross section
 	// leading order contribution and scale dependent born dependent terms
 	// will eventually add the other nlo terms ...
 	m_grids[0][iobs]->convolute( _pdf1, _pdf2, m_genpdf[0], alphas, m_leading_order,   0, rscale_factor, fscale_factor, Escale);
 	m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 0, rscale_factor, fscale_factor, Escale);
       }
       else if ( _nloops==-1 ) { 
 	label = "nlo only";
 	// next to leading order cross section
 	// leading order contribution and scale dependent born dependent terms
 	m_grids[1][iobs]->convolute( _pdf1, _pdf2, m_genpdf[1], alphas, m_leading_order+1, 0, rscale_factor, fscale_factor, Escale);
       }
 
     }
 
     /// wait on convolutions
 
     if    ( _nloops==0 ) { 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[0][iobs]->ready() ) dsigma = m_grids[0][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
     else if ( _nloops==1 ) { 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigLO = 0;
 	double dsigNLO = 0;
 	if ( m_grids[0][iobs]->ready() ) dsigLO  = m_grids[0][iobs]->xsec(); 
 	if ( m_grids[1][iobs]->ready() ) dsigNLO = m_grids[1][iobs]->xsec(); 
 	double dsigma  = dsigLO + dsigNLO;
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
     else if ( _nloops==-1 ) { 
       for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) {
 	double dsigma = 0;
 	if ( m_grids[1][iobs]->ready() ) dsigma = m_grids[1][iobs]->xsec(); 
 	double deltaobs = m_obs_bins->GetBinLowEdge(iobs+2)-m_obs_bins->GetBinLowEdge(iobs+1);      
 	hvec.push_back( invNruns*Escale2*dsigma/deltaobs );
       }
     }
 
 
   }
 
 
   /// now combine bins if required ...
 
   std::vector<bool> applied(m_corrections.size(),false);
 
   /// apply corrrections on the *uncombined* bins
   if ( getApplyCorrections() ) applyCorrections(hvec, applied);
   else { 
     for ( unsigned i=m_corrections.size() ; i-- ; ) if ( getApplyCorrection(i) ) applied[i] = applyCorrection(i,hvec);
   }
 
 
   /// combine bins if required
   if ( m_combine.size()!=0 ) { 
     combineBins( hvec );
 
     /// apply additional corrrections on the *combined* bins
     /// NB: Only apply those that have not already been applied
     if ( getApplyCorrections() ) applyCorrections(hvec, applied);
     else { 
       for ( unsigned i=m_corrections.size() ; i-- ; ) if ( getApplyCorrection(i) && !applied[i] ) applied[i] = applyCorrection(i,hvec);
     }
   }
 
 
   /// check all corrections have been applied correctly
   if ( getApplyCorrections() ) { 
     unsigned appliedcorrections = 0;
     for ( unsigned i=applied.size() ; i-- ; ) if ( applied[i] ) appliedcorrections++;
     if ( appliedcorrections!=applied.size() ) throw grid::exception( "correction vector size does not match data "  ); 
   }
   else { 
     unsigned Ncorrections = 0;
     unsigned appliedcorrections = 0;
     for ( unsigned i=m_corrections.size() ; i-- ; ) { 
       if ( getApplyCorrection(i) ) { 
 	Ncorrections++;
 	if ( applied[i] ) appliedcorrections++;
       }
     }
     if ( appliedcorrections!=Ncorrections ) throw grid::exception( "correction vector size does not match data "  ); 
   }
 
   //  double _ctime = appl_timer_stop(_ctimer);
   //  std::cout << "grid::convolute() " << label << " convolution time=" << _ctime << " ms" << std::endl;
 
   
   cache1.stats();
   if ( cache2.ncalls() ) cache2.stats();
   
   return hvec;
 }
 
 
 
 
 TH1D* appl::grid::convolute(void (*pdf)(const double& , const double&, double* ), 
 			    double (*alphas)(const double& ), 
 			    int     _nloops, 
 			    double  rscale_factor,
 			    double  fscale_factor,
 			    double Escale )
 {
   return convolute( pdf, 0, alphas, _nloops, rscale_factor, fscale_factor, Escale );
 }
 
 
 
 
 /// a dirty hack to tell the sub grid it should only 
 /// use a single subprocess
 
 std::vector<double> appl::grid::vconvolute_subproc(int subproc,
 						   void (*pdf)(const double& , const double&, double* ), 
 						   double (*alphas)(const double& ), 
 						   int     _nloops, 
 						   double  rscale_factor, double Escale ) 
 { 
   /// set the subprocess index - this is tested by the 
   /// igrid convolution
   m_subproc = subproc;
   std::vector<double> xsec = vconvolute( pdf, 0, alphas, _nloops, rscale_factor, rscale_factor, Escale ); 
   m_subproc = -1;
 
   return xsec;
 
 }
 
 
 
 double appl::grid::vconvolute_bin(int bin,
 				  void (*pdf)(const double& , const double&, double* ), 
 				  double (*alphas)(const double&) ) {
   /// set the subprocess index - this is tested by the 
   /// igrid convolution
   m_bin = bin;
   std::vector<double> xsec = vconvolute( pdf, alphas );
   m_bin = -1;
   return xsec[bin];
 }
 
 
 
 
 
 TH1D* appl::grid::convolute(void (*pdf1)(const double& , const double&, double* ), 
 			    void (*pdf2)(const double& , const double&, double* ), 
 			    double (*alphas)(const double& ), 
 			    int     _nloops, 
 			    double  rscale_factor,
 			    double  fscale_factor,
 			    double Escale ) 
 {
   //  int nbins = m_obs_bins->GetNbinsX();
 
   TH1D* h = 0;
   if ( m_combine.size() ) { 
     
     //    nbins = m_combine.size();
 
     std::vector<double> limits(m_combine.size()+1);
     
     int i=0;
     limits[0] = m_obs_bins->GetBinLowEdge(i+1);
     for ( unsigned ib=0 ; ib<m_combine.size() ; ib++) { 
       i += m_combine[ib]; 
       limits[ib+1] = m_obs_bins->GetBinLowEdge(i+1);
     }
 
     h = new TH1D("xsec", "xsec", m_combine.size(), &limits[0] );
     h->SetDirectory(0);
   }
   else { 
     h = new TH1D(*m_obs_bins);
     h->SetName("xsec");
     h->SetDirectory(0);
   }
 
   std::vector<double> dvec = vconvolute( pdf1, pdf2, alphas, _nloops, rscale_factor, fscale_factor, Escale );
  
   for ( unsigned i=0 ; i<dvec.size() ; i++ ) { 
       h->SetBinContent( i+1, dvec[i] );
       h->SetBinError( i+1, 0 );
   }
   
   return h;
   
 }
 
 
 
 
 
 TH1D* appl::grid::convolute_subproc(int subproc,
 				    void (*pdf)(const double& , const double&, double* ), 
 				    double (*alphas)(const double& ), 
 				    int     _nloops, 
 				    double  rscale_factor, double Escale ) {
   
     TH1D* h = new TH1D(*m_obs_bins);
     h->SetName("xsec");
     
     std::vector<double> dvec = vconvolute_subproc( subproc, pdf, alphas, _nloops, rscale_factor, Escale );
     
     for ( unsigned i=0 ; i<dvec.size() ; i++ ) { 
       h->SetBinContent( i+1, dvec[i] );
       h->SetBinError( i+1, 0 );
     }
   
     return h;
   
 }
 
 
 
 
 
 void appl::grid::optimise(bool force, int extrabins ) {
   if ( !force && m_optimised ) return;
   m_optimised = true;
   m_read = false;
   std::cout << "grid::optimise() " << std::endl;
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  { 
       m_grids[iorder][iobs]->optimise(extrabins);
     }
   }
   m_obs_bins->Reset();
 }
 
 
 
 void appl::grid::optimise(int NQ2, int Nx) {  optimise(NQ2, Nx, Nx);  }
 
 
 
 void appl::grid::optimise(int NQ2, int Nx1, int Nx2) {
   m_optimised = true;
   m_read = false;
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     for ( int iobs=0 ; iobs<Nobs_internal() ; iobs++ )  { 
       std::cout << "grid::optimise() bin " << iobs << "\t";
       m_grids[iorder][iobs]->optimise(NQ2, Nx1, Nx2);
     }
   }
   m_obs_bins->Reset();
 }
 
 
 
 
 // redefine the limits by hand
 void appl::grid::redefine(int iobs, int iorder,
 			  int NQ2, double Q2min, double Q2max, 
 			  int Nx,  double  xmin, double  xmax ) 
 { 
   
   if  ( iorder>=m_order ) { 
     std::cerr << "grid does not extend to this order" << std::endl;
     return;
   }
   
   if ( iobs<0 || iobs>=Nobs_internal() ) { 
     std::cerr << "observable bin out of range" << std::endl;
     return;
   }
   
   if ( iorder==0 ) { 
     std::cout << "grid::redefine() iobs=" << iobs 
 	      << "NQ2="  << NQ2 << "\tQmin=" << std::sqrt(Q2min) << "\tQmax=" << std::sqrt(Q2max) 
 	      << "\tNx=" << Nx  << "\txmin=" <<            xmin  << "\txmax=" <<            xmax << std::endl; 
   }
   
   igrid* oldgrid = m_grids[iorder][iobs];
   
   //  m_grids[iorder][iobs]->redefine(NQ2, Q2min, Q2max, Nx, xmin, xmax);
 
   m_grids[iorder][iobs] = new igrid(NQ2, Q2min, Q2max, oldgrid->tauorder(),
 				    Nx,  xmin,  xmax,  oldgrid->yorder(), 
 				    oldgrid->transform(), m_genpdf[iorder]->Nproc());
 
   m_grids[iorder][iobs]->setparent( this ); 
 
   delete oldgrid;
 }
   
 
 
 void appl::grid::setBinRange(int ilower, int iupper, double xScaleFactor) { 
   if ( ilower>=0 && iupper <Nobs_internal() ) {  
     double lower = getReference()->GetBinLowEdge(ilower+1);
     double upper = getReference()->GetBinLowEdge(iupper+2); 
     setRange( lower, upper, xScaleFactor );
   }
 }
 
 
 
 /// Fixme: need to fix this so that if the m_combine vector has values 
 ///        then the range of this vector is also modified
 ///        and then the combined reference recalculated:
 ///        at the moment, if the range is changed, the combine vector
 ///        needs to be recalculated and reset by hand
 
 void appl::grid::setRange(double lower, double upper, double xScaleFactor) { 
   
   std::cout << "grid::SetRange() " << lower << " " << upper << std::endl; 
 
   std::vector<bool>   used;
   std::vector<double> limits;
   std::vector<double> contents;
   std::vector<double> errors;
 
   m_combine = std::vector<int>(0);
 
   /// get the occupied bins
   //  int Nbins = 0;
   double last = 0;
   for ( int i=1 ; i<=m_obs_bins->GetNbinsX() ; i++ ) { 
     double bin = m_obs_bins->GetBinCenter(i);
     if ( bin>lower && bin<upper ) { 
       limits.push_back( m_obs_bins->GetBinLowEdge(i) );
       contents.push_back( m_obs_bins->GetBinContent(i) );
       errors.push_back( m_obs_bins->GetBinError(i) );
       //  std::cout << "keep bin " << i << " : " << m_obs_bins->GetBinLowEdge(i) << " - " << m_obs_bins->GetBinLowEdge(i+1) << std::endl;
       last = m_obs_bins->GetBinLowEdge(i+1);
       used.push_back(true);
     }
     else { 
       // std::cout << "skip bin " << i << " : " << m_obs_bins->GetBinLowEdge(i) << " - " << m_obs_bins->GetBinLowEdge(i+1) << std::endl;
       used.push_back(false);
     }
   }
 
   int firstbin = 0;
   int lastbin  = 0;
 
   for ( unsigned i=0 ; i<used.size() ; i++ ) { 
     if ( used[i] ) { 
       firstbin = i;
       break;
     }
   }
   
   for ( unsigned i=used.size() ; i-- ; ) { 
     if ( used[i] ) { 
       lastbin = i;
       break;
     }
   }
 
   std::cout << "grid::SetRange() bins chosen " << firstbin << " - " << lastbin << std::endl;
 
   std::cout << "before setrange: " << m_obs_bins->GetNbinsX() << std::endl;
 
 
   /// copy the range of the reference histogram
   if ( limits.size()>0 ) limits.push_back( last );
   else { 
     throw grid::exception( "new range does not include any bins"  ); 
   }
 
   if ( xScaleFactor!=1 ) { 
     for ( unsigned i=0 ; i<limits.size(); i++ ) limits[i] *= xScaleFactor;
   }
 
   /// delete combined reference if it exists
 
   if ( m_obs_bins_combined != m_obs_bins ) delete m_obs_bins_combined;
   
 
   /// save bins somewhere so can overwrite them 
   TH1D* h = m_obs_bins;
   
   m_obs_bins = new TH1D( h->GetTitle(), h->GetName(), limits.size()-1, &(limits[0]) );
   
   for ( int i=0 ; i<m_obs_bins->GetNbinsX() ; i++ ) { 
     m_obs_bins->SetBinContent( i+1, contents[i] );
     m_obs_bins->SetBinError( i+1, errors[i] );
   }
 
   std::cout << "after combine:  " << m_obs_bins->GetNbinsX() << std::endl;
 
   /// copy the igrids for the observable bins in the range 
 
   igrid** grids[appl::MAXGRIDS];
 
   /// save old grids
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) grids[iorder] = m_grids[iorder];
   
   int _Nobs = m_obs_bins->GetNbinsX();
 
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     m_grids[iorder] = new igrid*[_Nobs];
     int iobs = 0;
     for ( int igrid=0 ; igrid<h->GetNbinsX() ; igrid++ ) {
       if ( used[igrid] ) m_grids[iorder][iobs++] = grids[iorder][igrid];
       else               delete grids[iorder][igrid];                           
     }
   }
   
   m_obs_bins_combined = m_obs_bins; 
 
   delete h;
 
 }
 
 
 /// methods to handle the documentation
 void appl::grid::setDocumentation(const std::string& s) { m_documentation = s; }
 void appl::grid::addDocumentation(const std::string& s) {   
   if ( m_documentation.size() ) m_documentation += s;
   else                          setDocumentation(s);    
 }
 
 
 
 
 
 
 /// methods to handle bin-by-bin corrections
 
 /// add a correction as a std::vector
 void appl::grid::addCorrection( std::vector<double>& v, const std::string& label, bool ) {
   //  std::cout << "addCorrections(vector) " << v.size() << " " << m_obs_bins->GetNbinsX() << std::endl;
   if ( v.size()==unsigned(m_obs_bins->GetNbinsX()) || v.size()==unsigned(m_obs_bins_combined->GetNbinsX()) ) {
     m_corrections.push_back(v);
     m_correctionLabels.push_back(label);
     m_applyCorrection.push_back(false);
     //  std::cout << "appl::grid::addCorrection(vector) now " << m_corrections.size() << " corrections" << std::endl;
   }
 }
 
 
 /// add a correction by histogram
 void appl::grid::addCorrection(TH1D* h, const std::string& label, double scale, bool ) {
   
   TH1D* hobs = 0;
   
   if      ( h->GetNbinsX()==m_obs_bins->GetNbinsX() )          hobs = m_obs_bins;
   else if ( h->GetNbinsX()==m_obs_bins_combined->GetNbinsX() ) hobs = m_obs_bins_combined;
 
 
   if ( hobs ) { 
     for ( int i=1 ; i<=h->GetNbinsX()+1 ; i++ ) { 
       if ( std::fabs(h->GetBinLowEdge(i+1)*scale-hobs->GetBinLowEdge(i+1))>1e-10 ) { 
 	std::cerr << "bins " << h->GetBinLowEdge(i+1) << " " << hobs->GetBinLowEdge(i+1) << std::endl; 
 	std::cerr << "grid::addCorrection(TH1D* h): bin mismatch, not adding correction" << std::endl;
 	return;
       }
     }
 
     std::vector<double> v(h->GetNbinsX());
     for ( int i=0 ; i<h->GetNbinsX() ; i++ ) v[i] = h->GetBinContent(i+1);
     if ( label=="" ) addCorrection(v, h->GetName());
     else             addCorrection(v, label);
   }
   else { 
     std::cerr << "grid::addCorrection(TH1D* h): bin mismatch, not adding correction" << std::endl;
   }
   
 }
 
 
 
 // find the number of words used for storage
 int appl::grid::size() const { 
     int _size = 0;
     for( int iorder=0 ; iorder<m_order ; iorder++ ) {
       for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) _size += m_grids[iorder][iobs]->size();
     }
     return _size;
 }
 
 
 /// apply corrections to a std::vector
 void appl::grid::applyCorrections(std::vector<double>& v, std::vector<bool>& applied) {
  
   if ( applied.size()!=m_corrections.size() ) throw grid::exception( "wrong number of corrections expected" ); 
  
   for ( unsigned i=m_corrections.size() ; i-- ; ) { 
  
     std::vector<double>& correction = m_corrections[i];
     
     if ( applied[i] || v.size()!=correction.size() ) continue; /// correction applied already or wrong size     
 
     for ( unsigned j=v.size() ; j-- ; ) v[j] *= correction[j];
     applied[i] = true;
   }
   //  std::cout << "grid::applyCorrections(vector) done" << std::endl;
 }
 
 
 
 
 /// apply correction to a std::vector
 bool appl::grid::applyCorrection(unsigned i, std::vector<double>& v) {
 
   if ( i>=m_corrections.size() ) throw grid::exception( "correction index out of range"  ); 
  
   std::vector<double>& correction = m_corrections[i];
 
   if ( v.size()!=correction.size() ) return false; /// wrong size
 
   for ( unsigned k=v.size() ; k-- ; ) v[k] *= correction[k];
   return true;
 }
 
 
 
 
 
 /// do a deep comparison of all the different sub processes - if any 
 /// are the same, then remove them
 void appl::grid::shrink(const std::string& name, int ckmcharge) { 
 
   std::cout << "appl::grid::shrink()" << std::endl;
 
   std::string label[3] = { "LO", "NLO", "NNLO" };
 
   std::string genpdfname="";
 
   std::vector<std::vector<int> > keep(m_order);
 
   bool found = false;
 
 
   /// loop over orders 
   for( int iorder=0 ; iorder<m_order ; iorder++ ) {
 
     std::cout << "appl::grid::shrink() order " << iorder << std::endl;
     
     std::vector< std::vector<int> > pdf_combinations;
     pdf_combinations.reserve( Nobs_internal() );
     
     //    std::cout << "appl::grid::shrink() observable bins " << Nobs_internal() << std::endl;
 
     /// ... for each observable bin ...
  
     for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) { 
       
       //      std::cout << "shrink() order: " << iorder << "\t obs: " << iobs;
       
       igrid* ig = m_grids[iorder][iobs];
     
       std::map< int, std::vector<int> > same;
 
       std::set<int> duplicates;
 
       std::set<int> empty;
 
       for ( int i=0 ; i<ig->SubProcesses() ; i++ ) { 
 
 	if ( duplicates.find(i)!=duplicates.end() ) continue;
 
 	//	std::pair< std::map< int, std::vector<int> >, bool> itr = 
 	//	  same.insert( std::map< int, std::vector<int> >::value_type( i, std::vector<int>() ) );
 	//	
 	//  if ( itr->second==false ) std::cerr << "OH SHIT!!!" << std::endl;
 
 	int isize = ig->weightgrid(i)->size();
 
 	if ( isize==0 ) { 
 	  empty.insert(i); 
 	  continue;
 	}
 
 	//	std::cout << i << " : ";
 
 	keep[iorder].push_back( i );
 
 	std::vector<int> vec; 
 	
 	found = true;
 
 	for ( int j=i+1 ; j<ig->SubProcesses() ; j++ ) {
 
 	  if ( duplicates.find(j)!=duplicates.end() ) continue;
 	  if ( empty.find(j)!=empty.end() ) continue;
 
 	  int jsize = ig->weightgrid(j)->size();
 	  
 	  if ( jsize==0 ) { 
 	    empty.insert(j); 
 	    continue;
 	  }
 	  
 	  // std::cout << "\t" << j << ":" << jsize;
 	   
 	  if ( (*ig->weightgrid(i)) == (*ig->weightgrid(j)) ) { 
 	    //  std::cout << "!"; 
 	    vec.push_back(j);
 	    duplicates.insert(j);
 	  }
 
 	}
 	//	std::cout << std::endl;
 
 	same.insert( std::map< int, std::vector<int> >::value_type( i, vec ) );
 
       }
 
 
       if ( same.empty() ) continue;
 
       //      for ( int ik=0 ; ik<m_order ; ik++ ) std::cout << m_genpdf[iorder]->name() << std::endl;
 
       lumi_pdf*  _pdf = 0;
       if ( m_genpdf[iorder]->name().find(".config")!=std::string::npos ){ 
 	_pdf = dynamic_cast<lumi_pdf*>(m_genpdf[iorder]);
 	if ( _pdf ) std::cout << "read lumi pdf: " << *_pdf << std::endl; 
       }
       else { 
 	std::cerr << __FUNCTION__ << "\tpdf not found:" << m_genpdf[iorder]->name() << std::endl;
 	return;
       }
 
       std::vector<combination> combinations;
 
       int i=0;
 
       std::map< int, std::vector<int> >::iterator itr  = same.begin();
       std::map< int, std::vector<int> >::iterator iend = same.end();
 
       while ( itr!=iend ) { 
 
 	std::vector<int>& v = itr->second;
 	//	std::cout << "\t" << i++ << "\tproc: " << itr->first << ":" << sizeitr->second << "\t" << itr->second << " ( " << (*_pdf)[itr->first] << " )" << std::endl;
 	//	std::cout << "\t" << i << "\tproc: " << itr->first << ":" << sizeitr->second << "\t" << (*_pdf)[itr->first] << std::endl;
 	//	for ( unsigned iv=0 ; iv<v.size() ; iv++ ) std::cout << "\t\t\t" << v[iv] << "\t" << (*_pdf)[v[iv]] << std::endl;
 
 	std::vector<int> c(2);
 	c[0] = i;
 	c[1] = 0;
 
 	const combination& comb = (*_pdf)[itr->first]; 
 
 	for ( unsigned ic=0 ; ic<comb.size() ; ic++ ) { 
 	  c[1]++;
 	  c.push_back( comb[ic].first );
 	  c.push_back( comb[ic].second );
 	}
 
 	for ( unsigned iv=0 ; iv<v.size() ; iv++ ) { 
 
 	  const combination& comb = (*_pdf)[v[iv]]; 
 	  
 	  for ( unsigned ic=0 ; ic<comb.size() ; ic++ ) { 
 	    c[1]++;
 	    c.push_back( comb[ic].first );
 	    c.push_back( comb[ic].second );
 	  }
 	}
 
 	combinations.push_back( combination( c ) );
 
 	i++;
 	itr++;
       }
 
       //      std::set< int>::iterator eitr  = empty.begin();
       //      std::set< int>::iterator eend  = empty.end();
 
       //      std::cout << "\tempty: " << empty.size() << " sub-processes";
       //      while ( eitr!=eend ) std::cout << " " << (*eitr++);
       //      std::cout << std::endl; 
      
       lumi_pdf newpdf( "newpdf.config", combinations, 0 );
 
       pdf_combinations.push_back( newpdf.serialise() );
       
       //      std::cout << newpdf << std::endl;
 
       if ( found ) break;
 
     }
 
     bool common = true;
     for ( unsigned ipdf=1 ; ipdf<pdf_combinations.size() ; ipdf++ ) { 
     
       if ( pdf_combinations[ipdf].size()>0 ) {
 	if ( pdf_combinations[ipdf]!=pdf_combinations[ipdf-1] ) { 
 	  std::cout << "pdfs " << ipdf << " and " << ipdf-1 << " don't match" << std::endl; 
 	  // << lumi_pdf("duff.config", pdf_combinations[ipdf]) 
 	  common = false;
 	}  
       } 
     }      
     
     /// create the actual lumi_pdfs from these combinations
 
     if ( common && pdf_combinations.size()>0 ) { 
       pdf_combinations[0].push_back(ckmcharge);
 
       lumi_pdf lpdf( name+label[iorder]+".config",   pdf_combinations[0]);
 
       std::cout << lpdf.name() << std::endl;
 
       lpdf.write( lpdf.name() );
 
       if ( genpdfname.size() ) genpdfname += ":";
       genpdfname += name+label[iorder]+".config";
     }  
   }
 
   std::cout << "appl::grid::shrink() genpdfname " << genpdfname << std::endl;
 
   //  if ( addpdf(m_genpdfname) )  findgenpdf( m_genpdfname );
 
   /// horray!! here we have the optimised pdf combinations written out, so now we 
   /// to delete the duplicated (or empty) grids ...
 
   /// loop over the igrids, telling each grid which processes to keep
   
   for( int iobs=0 ; iobs<Nobs_internal() ; iobs++ ) { 
     for( int iorder=0 ; iorder<m_order ; iorder++ ) {
       //      std::cout << "appl::grid::shrink()  obs " << iobs << "\torder " << iorder << std::endl;       
       m_grids[iorder][iobs]->shrink( keep[iorder] );
     }
   }
 
   /// ... and need to replace the genpdf with these new ones
   
   m_genpdfname = genpdfname;
   addpdf( m_genpdfname );
   findgenpdf( genpdfname );
 
   std::cout << "appl::grid::shrink()  new " 
 	    << "\t" << m_genpdf[0]->name() << ":" << m_genpdf[0]->size() 
 	    << "\t" << m_genpdf[1]->name() << ":" << m_genpdf[1]->size() 
 	    <<  std::endl;
 
 }
 
 
 void appl::grid::combineReference(bool force) { 
   //  std::cout << "grid::combineReference() " << m_obs_bins->GetNbinsX();
 
   if ( m_combine.empty() ) return;
 
   if ( force ) { 
     if ( m_obs_bins_combined ) { 
       if ( m_obs_bins_combined!=m_obs_bins ) delete m_obs_bins_combined;
       m_obs_bins_combined = 0;
     }
   }
 
   //  std::cout << "\tcombine " << m_obs_bins_combined << " " << m_obs_bins << std::endl;; 
 
   if ( m_obs_bins_combined && m_obs_bins_combined!=m_obs_bins) return; 
 
   std::vector<double> hvec(  m_obs_bins->GetNbinsX(), 0 );
   std::vector<double> hvece( m_obs_bins->GetNbinsX(), 0 );
 
 
   for ( int i=m_obs_bins->GetNbinsX() ; i-- ; )  { 
     hvec[i]  = m_obs_bins->GetBinContent( i+1 );
     hvece[i] = m_obs_bins->GetBinError( i+1 );
   }
 
   combineBins( hvec );
   combineBins( hvece, 2 );
 
   std::vector<double> limits(m_combine.size()+1);
   
   int i=0;
   limits[0] = m_obs_bins->GetBinLowEdge(i+1);
   for ( unsigned ib=0 ; ib<m_combine.size() ; ib++) { 
     i += m_combine[ib]; 
     limits[ib+1] = m_obs_bins->GetBinLowEdge(i+1);
   }
   
   /// need to make this a class variable set to 0 so we don't 
   /// recalculate this every time, only if m_combine changes 
   TH1D* h = new TH1D("reference", "xsec", m_combine.size(), &limits[0] );
   h->SetDirectory(0);
 
   for ( unsigned i=0 ; i<hvec.size() ; i++ ) { 
     h->SetBinContent( i+1, hvec[i] );
     h->SetBinError( i+1, hvece[i] );
   }
 
   m_obs_bins_combined = h;
 
   //  std::cout << " -> " << m_obs_bins->GetNbinsX() << std::endl;
   
 }
 
 
 
 void appl::grid::combineBins(std::vector<double>& hvec, int power ) const { 
   
   /// now combine bins if required ...
   
   if ( m_combine.size()!=0 ) { 
     /// need to go through, scaling by bin width, adding and then dividing by bin width again
     /// in the TH1D* version, will need to recalculate the bin limits to create the new histogram
     
     std::vector<double> _hvec(m_combine.size(),0);
     
     unsigned nbins = 0;
 
     unsigned i=0;
 
     for ( unsigned ic=0 ; ic<m_combine.size() ; ic++ ) { 
 
       nbins += m_combine[ic];
 
       if ( nbins>hvec.size() ) throw grid::exception( "too many bins specified for rebinning"  ); 
 
       double sigma = 0;
       double width = 0;
 
       for ( int ib=0 ; ib<m_combine[ic] && i<hvec.size() ; ib++, i++ ) { 
 	double deltaobs = m_obs_bins->GetBinLowEdge(i+2)-m_obs_bins->GetBinLowEdge(i+1);
 	if ( power==1 ) sigma +=  hvec[i]*deltaobs;
 	if ( power==2 ) sigma += (hvec[i]*deltaobs*hvec[i]*deltaobs);
 	width += deltaobs;
       }
 
       if ( power==1 ) _hvec[ic] = sigma/width;
       if ( power==2 ) _hvec[ic] = std::sqrt(sigma)/width;
     }
 
     hvec = _hvec;
   }
 }
 
 
 std::ostream& operator<<(std::ostream& s, const appl::grid& g) {
   
   s << "==================================================" << std::endl;
   
   //  s << "appl::grid version " << g.version() << "\t(" << g.subProcesses(0) << " initial states, " << g.Nobs_internal() << " observable bins)" << std::endl;
 
   std::string basis[5] = {  "-LO, ",  "-NLO, ",  "-NNLO, ", "-Xtra0", "-Xtra1" };  
   std::string order[appl::MAXGRIDS];
   for ( int i=0 ; i<appl::MAXGRIDS ; i++ ) { 
     if ( i<5) order[i] = basis[i];
     else      order[i] = "-Unknown";
   }
 
   s << "appl::grid version " << g.version() << "\t( "; 
 
   for ( int i=0 ; i<g.nloops()+1 ; i++ ) s << g.subProcesses(i) << order[i];
   
   s << "initial states, " << g.Nobs_internal() << " observable bins )" << std::endl;
   
   if ( g.isOptimised() ) s << "Optimised grid" << std::endl;
   
   if ( g.isSymmetric() ) s << "Symmetrised in x1, x2" << std::endl;
   else                   s << "Unsymmetrised in x1, x2" << std::endl;
   
   if ( g.getNormalised() ) s << "Normalised " << std::endl;
   
   s << "leading order of processes  "  << g.leadingOrder() << std::endl;
   s << "number of loops for grid    "  << g.nloops() << std::endl;   
   s << "x->y coordinate transform:  "  << g.getTransform() << std::endl;
 
   s << "genpdf in use: " << g.getGenpdf() << std::endl;
   s << "--------------------------------------------------" << std::endl;
   s << "Observable binning: [ " << g.Nobs_internal() 
     << " bins : " << g.obsmin() << ",  " << g.obsmax() << " ]" << std::endl;
 
   //  for( int iorder=0 ; iorder<g.nloops()+1 ; iorder++ ) {
   for( int iorder=0 ; iorder<g.ninternalgrids() ; iorder++ ) {
     s << "order: " << iorder << "\n";
     for( int iobs=0 ; iobs<g.Nobs_internal() ; iobs++ ) {
       s << "  " 
 	<< iobs << "\t" 
 	<< std::setprecision(5) << std::setw(5) << g.getReference()->GetBinLowEdge(iobs+1) << "\t- " 
 	<< std::setprecision(5) << std::setw(5) << g.getReference()->GetBinLowEdge(iobs+2) << "\t"; 
       s << "   " << *(g.weightgrid(iorder,iobs)) << std::endl;
     }
   }
 
   s << std::endl;
   
   return s;
 }
 
 
 void appl::grid::replaceBin( int iobs, grid& g ) { 
 
   std::cout << "replace bin " << iobs << std::endl;
  
   for ( int iorder=0 ; iorder<m_order ; iorder++ ) { 
     add_igrid( iobs, iorder, g.m_grids[iorder][iobs] );
   }
 
 #if 0
   int iorder = 2;
   //  for ( int iorder=0 ; iorder<m_order ; iorder++ ) 
   { 
     std::cout << "order: " << iorder << std::endl;
     
     std::cout <<  m_grids[iorder][iobs] << " " << g.m_grids[iorder][iobs] << std::endl;
 
     igrid* ig = m_grids[iorder][iobs];
 
     igrid* hig = new igrid( *ig );
 
     //    m_grids[iorder][iobs] = 
     //    m_grids[iorder][iobs]->setparent( this ); 
   }
 #endif
 
   std::cout << "fixing reference: " << iobs << std::endl;
   getReference_internal()->SetBinContent(iobs+1, g.getReference_internal()->GetBinContent(iobs+1) );
   getReference_internal()->SetBinError(iobs+1, g.getReference_internal()->GetBinError(iobs+1) );
 
   /// grrr, is this correct - why, why, why, why why o why can't we 
   /// just use the correct binning in the first place !!!
   combineReference(true); 
 }
 
 
 
 static bool  no_threads = duff();