source: Sophya/trunk/SophyaPI/PIext/pawexecut.cc@ 2809

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <iostream>
#include <typeinfo>

#include "sopnamsp.h"
#include "strutil.h"
#include "strutilxx.h"
#include "histos.h"
#include "histos2.h"
#include "hisprof.h"
#include "histerr.h"
#include "ntuple.h"
#include "datatable.h"

#include "pawexecut.h"
#include "nobjmgr.h"
#include "servnobjm.h"
#include "nomgadapter.h"
#include "pistdimgapp.h"
#include "pihisto.h"

#ifdef SANS_EVOLPLANCK
#include "cvector.h"
#include "matrix.h"
#else
#include "tmatrix.h"
#include "tvector.h"
#endif

/* Reza + cmv  13/10/99 */

inline bool __vector_sort_up_r8(r_8 x1,r_8 x2) {return x1<x2;}
inline bool __vector_sort_down_r8(r_8 x1,r_8 x2) {return x1>x2;}
inline bool __vector_sort_up_r4(r_4 x1,r_4 x2) {return x1<x2;}
inline bool __vector_sort_down_r4(r_4 x1,r_4 x2) {return x1>x2;}
inline bool __vector_sort_up_i4(int_4 x1,int_4 x2) {return x1<x2;}
inline bool __vector_sort_down_i4(int_4 x1,int_4 x2) {return x1>x2;}
inline bool __vector_sort_up_u8(uint_8 x1,uint_8 x2) {return x1<x2;}
inline bool __vector_sort_down_u8(uint_8 x1,uint_8 x2) {return x1>x2;}

uint_4 PAWExecutor::autoc_counter_ = 0;

/* methode */
PAWExecutor::PAWExecutor(PIACmd *piac, PIStdImgApp* app)
: mApp(app)
{
string kw, usage;
string hgrp = "pawCmd";

kw = "reset";
usage = "Reset histograms vectors or matrix";
usage += "\n reset nameobj";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/plot";
usage = "Plot NTuple variables \"a la paw\" (alias n/pl)";
usage += "\n n/plot nameobj.x_exp [cut] [w_exp] [loop] [gratt]";
usage += "\n n/plot nameobj.y_exp%x_exp [cut] [loop] [gratt]";
usage += "\n n/plot nameobj.z_exp%y_exp%x_exp [cut] [loop] [gratt]";
usage += "\n        for default use ! , loop=i1[:i2[:di]]";
usage += "\n for 1 dimensional (1D) projection:";
usage += "\n   use graphic option \"keepbin\" to keep previous";
usage += "\n   plot binning for 1D distribution";
usage += "\n  Related commands: plot2dw plot3d";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/pl";
usage = "alias to n/plot (see n/plot)";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/proj";
usage = "Project NTuple in histogram (1D or 2D) a la paw";
usage += "\n n/proj nameproj nameobj.x_exp [cut] [w_exp] [loop] [gratt]";
usage += "\n n/proj nameproj nameobj.y_exp%x_exp [cut] [w_exp] [loop] [gratt]";
usage += "\n        for default use ! , loop=i1[:i2[:di]]";
usage += "\n for 1 dimensional (1D) projection:";
usage += "\n   no display is performed if \"nameproj\" is an existing histogram";
usage += "\n   unless a graphic option \"gratt\" is given";
usage += "\n  Related commands: projh1d projh2d projprof exptovec";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/scan";
usage = "Scan NTuple a la paw";
usage += "\n n/scan nameobj[.exp1%exp2%exp3] cut loop";
usage += "\n               [-f:filename] [list_of_variables]";
usage += "\n  loop : iev1[:iev2[:diev]] or !";
usage += "\n  cut  : cut expression or 1. or !";
usage += "\n  list_of_variables : default is all variables";
usage += "\n                    : var1 var2 var3 ... varn";
usage += "\n                    : var1 : var2   (from var1 to var2)";
usage += "\n                    :  : var2   (from first variable to var2)";
usage += "\n                    : var1 :    (from var1 to last variable)";
usage += "\n                  ex: \"v1 : v3 v7 v4 : v6 v2 v9 :\"";
usage += "\n exp1%exp2%exp3 :";
usage += "\n           if given add exp1,exp2,exp3 to the variable list";
usage += "\n -f:filename : write into \"filename\", Default is to stdout";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/read";
usage = "Read columns in an ASCII file and fill a NTuple or a DataTable";
usage += "\n n/read nt fascii [options] var_1,c_1 var_2,c_2 ... var_n,c_n ";
usage += "\n nt : NTuple or DataTable (see options) name";
usage += "\n var_i,c_i : variable name, associated column in ASCII file [0,n[";
usage += "\n where [options] are:";
usage += "\n   \"-dt\": create a DataTable of double instead of a NTuple";
usage += "\n   \"=s\": separator character is \'s\' (could be \"\t\")";
usage += "\n   \"-^abcd\": do not read lines beginning with string \"abcd\" ";
usage += "\n   \"+^abcd\": read only lines beginning with string \"abcd\" ";
usage += "\n   \"-abcd\": do not read lines which contain string \"abcd\" ";
usage += "\n   \"+abcd\": read only lines which contain string \"abcd\" ";
usage += "\n these options may be repeated (ex: \"-^abcd\" \"-^xyz\") ";
usage += "\n   - in case of \"do not read\" options are added with logical AND ";
usage += "\n   - in case of \"read only\" options are added with logical OR ";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/merge";
usage = "Merge ntuples (by adding lines)";
usage += "\n n/merge nt nt_1 nt_2 ... nt_n";
usage += "\n Merge ntuples nt_i lines into ntuple nt";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "n/merge/col";
usage = "Merge ntuples (by adding columns)";
usage += "\n n/merge/col nt nt_1[,ext] nt_2[,ext] ... nt_n[,ext]";
usage += "\n Merge ntuples nt_i columns into ntuple nt";
usage += "\n ext : character string to add at the end of the ntuple variable names";
usage += "\n   if \"ext\"=\"!\" no extension is added";
usage += "\n   if no \"ext\" is given, automatic extension \"_0\" \"_1\" \"_2\" ...";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/integ";
usage = "Integrate a 1D histogram";
usage += "\n h/integ nameh1d [norm]";
usage += "\n norm<=0 means no normalisation (def norm=1)";
usage += "\n  Related commands: h/deriv v/integ v/deriv";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "v/integ";
usage = "Integrate a TVector / vector";
usage += "\n v/integ namevec [norm]";
usage += "\n norm<=0 means no normalisation (def norm=0)";
usage += "\n  Related commands: h/integ h/deriv v/deriv";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "v/sort";
usage = "Sort a vector into itself";
usage += "\n v/sort namevec [+1/-1]";
usage += "\n \"+1\" means increasing order, \"-1\" decreasing order";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/deriv";
usage = "Derivate a 1D histogram";
usage += "\n h/deriv nameh1d";
usage += "\n  Related commands: h/integ v/integ v/deriv";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "v/deriv";
usage = "Derivate a TVector / vector";
usage += "\n v/deriv namevec [deriv_option]";
usage += "\n deriv_option -1 replace v[i] with v[i]-v[i-1]";
usage += "\n               0 replace v[i] with (v[i+1]-v[i-1])/2 (default)";
usage += "\n              +1 replace v[i] with v[i+1]-v[i]";
usage += "\n  Related commands: h/integ h/deriv v/integ";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/rebin";
usage = "Rebin a 1D histogram or profile";
usage += "\n h/rebin nameh1d nbin";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/cadd";
usage = "Add a constant to an histogram, a vector or a matrix";
usage += "\n h/cadd namehisto val";
usage += "\n  Related commands: h/cmult h/oper";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/cmult";
usage = "Multiply an histogram, a vector or a matrix by a constant";
usage += "\n h/cmult namehisto val";
usage += "\n  Related commands: h/cadd h/oper";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/oper";
usage = "Operation on histograms vectors or matrices";
usage += "\n h/oper @ h1 h2 hres";
usage += "\n        hres = h1 @ h2   with @ = (+,-,*,/)";
usage += "\n For vectors and matrices, operations are elements by elements";
usage += "\n  Related commands: h/cadd h/cmult";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/plot/2d";
usage = "Specific plot for 2D histogrammes";
usage += "\n h/plot/2d nameh2d show        : infos on 2D histogramme";
usage += "\n h/plot/2d nameh2d h    [dopt] : plot 2D histogramme";
usage += "\n h/plot/2d nameh2d px   [dopt] : plot X projection";
usage += "\n h/plot/2d nameh2d py   [dopt] : plot Y projection";
usage += "\n h/plot/2d nameh2d bx n [dopt] : plot X band number n";
usage += "\n h/plot/2d nameh2d by n [dopt] : plot Y band number n";
usage += "\n h/plot/2d nameh2d sx n [dopt] : plot X slice number n";
usage += "\n h/plot/2d nameh2d sy n [dopt] : plot Y slice number n";
usage += "\n                      n < 0 means Show Info";
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/put_vec";
usage = "Put content of vector (matrix) into content of histogram 1D or 2D";
usage += "\n h/put_vec nameh1d namevector [cont,err2]";
usage += "\n h/put_vec nameh2d namematrix [cont,err2]";
usage += "\n     cont : put into histogramme content";
usage += "\n     err2 : put into histogramme error^2";
usage += "\n  Related commands: h/get_vec";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/get_vec";
usage = "Get content of histogram 1D profile or 2D into vector (matrix)";
usage += "\n h/get_vec nameh1d namevector [cont,err2,absc] [proj]";
usage += "\n h/get_vec nameh2d namematrix [cont,err2,absc]";
usage += "\n   cont : get histogramme content";
usage += "\n   err2 : get histogramme error^2";
usage += "\n   absc : get histogramme low bin abscissa (1D only)";
usage += "\n     proj :";
usage += "\n            show  : show available projections for Histo2D";
usage += "\n            px    : get X projection";
usage += "\n            py    : get Y projection";
usage += "\n            bx n  : get X band number n";
usage += "\n            by n  : get Y band number n";
usage += "\n            sx n  : get X slice number n";
usage += "\n            sy n  : get Y slice number n";
usage += "\n  Related commands: h/put_vec";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/copy";
usage = "Copy content of object1 into object2";
usage += "\n objects are Vector,Matrix,Histo,Histo2D";
usage += "\n h/copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]";
usage += "\n   copy obj1Dfrom(i1->i2) into obj1Dto(ic1->)";
usage += "\n   copy obj2Dfrom(i1,j1->i2,j2) into obj2Dto(ic1,jc1->)";
usage += "\n Warning: elements from i1 to i2 included are copied";
usage += "\n If obj1Dto does not exist, is is created with size i2-i1+1";
usage += "\n    or obj2Dto with size i2-i1+1,j2-j1+1";
usage += "\n The adressed content of obj?Dfrom is overwritten";
usage += "\n The non-adressed content of obj?Dfrom is left unchanged";
usage += "\n  Related commands: copy";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/set/err";
usage = "Set Histo,Histo2D errors for range of bins or range of values";
usage += "\n h/set/err namehisto setvalue i1[:i2]  [j1[:j2]]";
usage += "\n     set error to setvalue for bin range i1:i2  j1:j2";
usage += "\n h/set/err namehisto setvalue v v1:v2";
usage += "\n     set error to setvalue for content values range v1:v2";
usage += "\n h/set/err namehisto setvalue e e1:e2";
usage += "\n     set error to setvalue for error values range v1:v2";
usage += "\n  Related commands: h/set/cont";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/set/cont";
usage = "Set Histo,Histo2D content for range of bins or range of values";
usage += "\n h/set/cont namehisto setvalue i1[:i2]  [j1[:j2]]";
usage += "\n     set content to setvalue for bin range i1:i2  j1:j2";
usage += "\n h/set/cont namehisto setvalue v v1:v2";
usage += "\n     set content to setvalue for content values range v1:v2";
usage += "\n h/set/cont namehisto setvalue e e1:e2";
usage += "\n     set content to setvalue for error values range v1:v2";
usage += "\n  Related commands: h/set/err";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/err";
usage = "Set Histo,Histo2D error to function of bin content value";
usage += "\n h/err namehisto expr_func";
usage += "\n  Related commands: h/set/err";  
piac->RegisterCommand(kw,usage,this,hgrp);

}

/* methode */
PAWExecutor::~PAWExecutor()
{
}

/* methode */
int PAWExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
{
if(kw == "reset") {
  reset(tokens); return(0);
} else if(kw == "n/plot" || kw == "n/pl") {
  n_plot(tokens); return(0);
} else if(kw == "n/proj") {
  n_proj(tokens); return(0);
} else if(kw == "n/scan") {
  n_scan(tokens); return(0);
} else if(kw == "n/read") {
  n_read(tokens); return(0);
} else if(kw == "n/merge") {
  n_merge(tokens); return(0);
} else if(kw == "n/merge/col") {
  n_merge_col(tokens); return(0);
} else if(kw == "h/integ") {
  h_integ(tokens); return(0);
} else if(kw == "v/integ") {
  v_integ(tokens); return(0);
} else if(kw == "v/sort") {
  v_sort(tokens); return(0);
} else if(kw == "h/deriv") {
  h_deriv(tokens); return(0);
} else if(kw == "v/deriv") {
  v_deriv(tokens); return(0);
} else if(kw == "h/rebin") {
  h_rebin(tokens); return(0);
} else if(kw == "h/cadd") {
  h_cadd(tokens); return(0);
} else if(kw == "h/cmult") {
  h_cmult(tokens); return(0);
} else if(kw == "h/oper") {
  h_oper(tokens); return(0);
} else if(kw == "h/plot/2d") {
  h_plot_2d(tokens); return(0);
} else if(kw == "h/put_vec") {
  h_put_vec(tokens); return(0);
} else if(kw == "h/get_vec") {
  h_get_vec(tokens); return(0);
} else if(kw == "h/copy") {
  h_copy(tokens); return(0);
} else if(kw == "h/set/err") {
  string dum = "err";
  h_set(dum,tokens); return(0);
} else if(kw == "h/set/cont") {
  string dum = "cont";
  h_set(dum,tokens); return(0);
} else if(kw == "h/err") {
  h_err(tokens); return(0);
} else return(1); 
}

/* methode */
bool PAWExecutor::IsThreadable(string const & keyw)
{
  if ( (keyw == "n/plot")  || (keyw == "n/pl") || 
       (keyw == "n/proj") )  return true;
  else return false;
}


/* methode */
void PAWExecutor::reset(vector<string>& tokens)
// Reset d'histogrammes, vecteurs et matrices
{
if(tokens.size() < 1)
  {cout<<"Usage: reset nameobj"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj == NULL)
  {cout<<"PAWExecutor::reset Error , Pas d'objet de nom "<<tokens[0]<<endl;
   return;}
string ctyp = typeid(*mobj).name();

#ifdef SANS_EVOLPLANCK
if(typeid(*mobj)==typeid(Vector))       {Vector*  ob=(Vector*) mobj; ob->Zero();}
else if(typeid(*mobj)==typeid(Matrix))  {Matrix*  ob=(Matrix*) mobj; ob->Zero();}
#else
 if(typeid(*mobj)==typeid(Vector))       {Vector*  ob=(Vector*) mobj; (*ob) = 0.; }
// ob->DataBlock().Reset(0.);}
 else if(typeid(*mobj)==typeid(Matrix))   {Matrix*  ob=(Matrix*) mobj; (*ob) = 0.; }
//ob->DataBlock().Reset(0.);}
#endif
else if(typeid(*mobj)==typeid(Histo))    {Histo* ob=(Histo*)  mobj; ob->Zero();}
else if(typeid(*mobj)==typeid(HProf))    {HProf* ob=(HProf*)  mobj; ob->Zero();}
else if(typeid(*mobj)==typeid(HistoErr)) {HistoErr* ob=(HistoErr*)  mobj; ob->Zero();}
else if(typeid(*mobj)==typeid(Histo2D))  {Histo2D* ob=(Histo2D*)mobj; ob->Zero();}
else {
  cout<<"PAWExecutor::reset Error , No reset possible on "<<ctyp<<endl;
  return;
}

return;
}

/* methode */
void PAWExecutor::n_plot(vector<string>& tokens)
// Equivalent n/plot de paw
// Plot 1D
//   n/plot nameobj.x_exp [cut] [w_exp] [gratt]
// Plot 2D (plot2dw)
//   n/plot nameobj.y_exp%x_exp [cut] [w_exp] [gratt]
// Plot 3D (plot3d)
//   n/plot nameobj.z_exp%y_exp%x_exp [cut] [gratt]
{
if(tokens.size() < 1) {
  cout
  <<"Usage: n/plot nameobj.x_exp [cut] [w_exp] [loop] [gratt] [nomh1]"<<endl
  <<"       n/plot nameobj.y_exp%x_exp [cut] [loop] [gratt]"<<endl
  <<"       n/plot nameobj.z_exp%y_exp%x_exp [cut] [loop] [gratt]"<<endl
  <<"              for default use ! , loop=i1[:i2[:di]]"<<endl;
  return;
}
string nameobj,expx,expy,expz;
int_4 nvar = decodepawstring(tokens[0],nameobj,expx,expy,expz);
string expcut = "1"; string expwt = "1."; string loop = "";
string dopt = ""; string nameproj="";
if(tokens.size()>=2) expcut = tokens[1]; if(expcut=="!") expcut="1";

NamedObjMgr omg;
Services2NObjMgr* srvo = omg.GetServiceObj();

if(nvar<=0) {
  cout<<"PAWExecutor::n_plot Error: bad coding "<<tokens[0]<<endl;
} else if(nvar==1) { // c'est un plot 1D
  if(tokens.size()>=3) expwt    = tokens[2]; if(expwt=="!") expwt="1.";
  if(tokens.size()>=4) loop     = tokens[3]; if(loop=="!") loop="";
  if(tokens.size()>=5) dopt     = tokens[4]; if(dopt=="!") dopt="";
  if(tokens.size()>=6) nameproj = tokens[5];
  if(nameproj.length()<=0 || nameproj=="!") {
    nameproj = "/autoc/paw_n_plot1D_";
    AnyDataObj* mobj = omg.GetObj(nameproj);
    if(mobj!=NULL) {
      char buff[16]; autoc_counter_++; sprintf(buff,"%d",autoc_counter_);
      nameproj += buff;
    }
    // --- Si option "keepbin" on projete avec le meme binning que l'histo 1D precedent
    if(dopt.find("keepbin")<dopt.size()) {
      PIStdImgApp* piimapp = omg.GetImgApp();
      if(piimapp) {
        PIBaseWdg* pibwdg = piimapp->CurrentBaseWdg();
        if(pibwdg) {
          int nbdrw = pibwdg->NbDrawers();
          if(nbdrw>0) {
            for(int i=nbdrw-1;i>=0;i--) {
              PIDrawer* pidwr = pibwdg->GetDrawer(i);
              PIHisto* pih = NULL;
              if( (pih = dynamic_cast<PIHisto *>(pidwr))==NULL ) continue;
              Histo* h = pih->Histogram();
              if(h==NULL) continue;
              if(h->NBins()<1) continue;
              Histo* hsame = new Histo(h->XMin(),h->XMax(),h->NBins());
              omg.AddObj(hsame,nameproj);
              break;
            }
          }
        }
      }
    }
    // --- 
  }
  srvo->ProjectH1(nameobj,expx,expwt,expcut,nameproj,dopt,loop);
} else if(nvar==2) { // c'est un plot 2D
  if(tokens.size()>=3) loop = tokens[2]; if(loop=="!") loop="";
  if(tokens.size()>=4) dopt = tokens[3];
  string err = "";
  srvo->DisplayPoints2D(nameobj,expx,expy,err,err,expcut,dopt,loop);
} else {             // c'est un plot 3D
  if(tokens.size()>=3) loop = tokens[2]; if(loop=="!") loop="";
  if(tokens.size()>=4) dopt = tokens[3];
  srvo->DisplayPoints3D(nameobj,expx,expy,expz,expcut,dopt,loop);
}

return;
}


/* methode */
void PAWExecutor::n_proj(vector<string>& tokens)
// Equivalent n/proj de paw
// Project NTuple in histogram a la paw
// Dans un Histo 1D
//   n/proj nameproj nameobj.x_exp [cut] [w_exp] [gratt]
// Dans un Histo 2D ou un HProf (dans ce cas nameproj doit etre cree).
//   n/proj nameproj nameobj.y_exp%x_exp [cut] [w_exp] [gratt]
{
if(tokens.size()<2)
  {cout<<"Usage: n/proj nameproj nameobj.[y_exp%]x_exp [cut] [w_exp] [loop] [gratt]"<<endl
       <<"              for default use ! , loop=i1[:i2[:di]]"<<endl; return;}
string nameproj = tokens[0];
string nameobj,expx,expy,expz;
int_4 nvar = decodepawstring(tokens[1],nameobj,expx,expy,expz);
string expcut = "1"; string expwt = "1."; string loop = ""; string dopt = "";
if(tokens.size()>=3) expcut = tokens[2]; if(expcut=="!") expcut="1";
if(tokens.size()>=4) expwt  = tokens[3]; if(expwt=="!") expwt="1.";
if(tokens.size()>=5) loop   = tokens[4]; if(loop=="!") loop="";
if(tokens.size()>=6) dopt   = tokens[5];

NamedObjMgr omg;
Services2NObjMgr* srvo = omg.GetServiceObj();

if(nvar==1) {
  // c'est une projection dans un histo 1D
  srvo->ProjectH1(nameobj,expx,expwt,expcut,nameproj,dopt,loop);
} else if(nvar==2) {
  // c'est une projection dans un histo2D
  // OU un HProf si nameproj est un HProf un deja defini
  AnyDataObj* mobj = omg.GetObj(nameproj);
  if(mobj==NULL)
    srvo->ProjectH2(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
  else if(dynamic_cast<HProf*>(mobj))
    srvo->ProjectHProf(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
  else
    srvo->ProjectH2(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
} else {
  cout<<"PAWExecutor::n_proj Error: bad coding "<<tokens[1]<<" nvar="<<nvar<<endl;
}

return;
}

/* methode */
void PAWExecutor::n_scan(vector<string>& tokens)
{
if(tokens.size()<3)
  {cerr<<"Usage: n/scan nameobj[.exp1%exp2%exp3] cut loop\n"
       <<"       [-f:filename] [list_of_variables]"<<endl; return;}

// decodage des premiers arguments
string nameobj,expr[4];
int_4 nexpr = decodepawstring(tokens[0],nameobj,expr[0],expr[1],expr[2]);
if(nexpr<4) {for(int_4 i=nexpr;i<4;i++) expr[i]="1.";}
string expcut = tokens[1]; if(expcut=="!") expcut="1";
string loop   = tokens[2]; if(loop=="!") loop="";
FILE* fout = NULL;
uint_4 ldebvar = 3;
if(tokens.size()>3) {
  if(tokens[3].find("-f:") == 0) {
    string filename = tokens[3].substr(3);
    if(filename.size()>0) {
      fout = fopen(filename.c_str(),"w");
      if(fout==NULL) 
        {cerr<<"PAWExecutor::n_scan Error, file "<<filename
             <<" not opened"<<endl; return;}
    }
    ldebvar++;
  }
}

// ntuple adaptateur
NamedObjMgr omg;
Services2NObjMgr& srvo = *omg.GetServiceObj();
NObjMgrAdapter* obja = omg.GetObjAdapter(nameobj); // Ne pas deleter
if(obja == NULL)
  {cerr<<"PAWExecutor::n_scan Error, ObjAdapter==NULL for "
       <<nameobj<<endl; return;}
bool adel = true;
NTupleInterface* objnt = obja->GetNTupleInterface(adel);
if(objnt == NULL)
  {cerr<<"PAWExecutor::n_scan Error, NTupleInterface==NULL for "
       <<nameobj<<endl; return;}

// function pour le choix
string  vardec = objnt->VarList_C("_zz6qi_");
PlotExprFunc f = srvo.LinkExprFunc(vardec,expr[0],expr[1],expr[2],expr[3],expcut);
if(!f)
  {cerr<<"PAWExecutor::n_scan Error, Creation PlotExprFunc"<<endl; 
   if(adel) delete objnt; if(fout) fclose(fout); return;}

// variables a imprimer
// "rien"          --> de 0 a ncol-1 (toutes les variables)
// : v1            --> de 0 a v1
// v1 :            --> de v1 a ncol-1
// v1 : v2         --> de v1 a v2 (si v2 apres v1)
//                     v1 et v2   (si v2 avant v1)
// v1 v2 v3 v4     --> v1 v2 v3 v4 (ordre indifferent)
// et toute combinaison : "v1 : v3 v7 v4 : v6 v2 v9 :"
//        --> v1 v2 v3  v7  v4 v5 v6  v2  v9 v10...v(ncol-1)
int_4 ncol = objnt->NbColumns();
if(ncol<=0)
  {cerr<<"PAWExecutor::n_scan Error, no columns for NTuple"<<endl;
   return;}
vector<int_4> varnum;
if(ldebvar>=tokens.size()) {   // Toutes les variables
  for(int_4 i=0;i<ncol;i++) varnum.push_back(i);
} else {   // Choix de certaines variables
  int_4 k,klast,kk; bool frlast=false;
  if(tokens[ldebvar]==":") {varnum.push_back(0); frlast=true;}
  else {k = objnt->ColumnIndex(tokens[ldebvar]); varnum.push_back(k);}
  ldebvar++;
  if(ldebvar<tokens.size()) for(uint_4 i=ldebvar;i<tokens.size();i++) {
    if(tokens[i]!=":") {  // pas un separateur
      k = klast = objnt->ColumnIndex(tokens[i]);
      if(frlast) klast = varnum[varnum.size()-1] + 1;
      if(klast>k) klast=k;
      for(kk=klast;kk<=k;kk++) varnum.push_back(kk);
      frlast=false;
    } else if(i==tokens.size()-1) {  // separateur a la fin
      k = ncol-1;
      klast = varnum[varnum.size()-1] + 1;
      if(klast>k) klast=k;
      for(kk=klast;kk<=k;kk++) varnum.push_back(kk);
    } else frlast=true;  // separateur pas a la fin
  }
}

vector<string> varname;
if(varnum.size()>0) for(int_4 i=0;i<(int)varnum.size();i++) {
  if(varnum[i]<0 || varnum[i]>=ncol)
    {cerr<<"PAWExecutor::n_scan Error, bad variable name at pos "
         <<i<<endl; if(adel) delete objnt; if(fout) fclose(fout); return;}
  string dum = objnt->ColumnName(varnum[i]);
  varname.push_back(dum);
}

// evenements a utiliser
int_8 k1=0, k2=objnt->NbLines(), dk=1;
srvo.DecodeLoopParameters(loop,k1,k2,dk);
if (k1<0) k1=0;
if (k2<0) k2=objnt->NbLines();
if (k2>(int_8)objnt->NbLines()) k2=objnt->NbLines();
if (dk<=0) dk=1;

// boucle sur les evenements et print
try  {
  int_4 i;
  if(fout) fprintf(fout,"#ev "); else printf("#ev ");
  for(i=0;i<(int)varname.size();i++)
    if(fout) fprintf(fout,"%s ",varname[i].c_str());
    else      printf(     "%s ",varname[i].c_str());
  if(nexpr>0) for(i=0;i<nexpr;i++)
    if(fout) fprintf(fout,"%s ",expr[i].c_str());
    else      printf(     "%s ",expr[i].c_str());
  if(fout) fprintf(fout,"\n"); else printf("\n");
 
  double xnt[5]={0,0,0,0,0};
  double* xn;
  for(int_8 k=k1; k<k2; k += dk)    {
    xn = objnt->GetLineD(k);
    if(f((int_8_exprf)k,xn,xnt,xnt+1,xnt+2,xnt+3) != 0) {
      if(fout) fprintf(fout,"%ld ",(long)k); else printf("%ld ",(long)k);
      for(i=0;i<(int)varnum.size();i++) {
        if(fout) fprintf(fout,"%g ",*(xn+varnum[i]));
        else      printf(     "%g ",*(xn+varnum[i]));
      }
      if(nexpr>0) for(i=0;i<nexpr;i++) {
        if(fout) fprintf(fout,"%g ",*(xnt+i));
        else      printf(     "%g ",*(xnt+i));
      }
      if(fout) fprintf(fout,"\n"); else printf("\n");
    }
  }
}   // fin du try
#ifdef SANS_EVOLPLANCK
CATCH(merr) {
  fflush(stdout); cout<<endl; cerr<<endl;
  string es = PeidaExc(merr);
  cerr<<"Services2NObjMgr::ComputeExpressions()  Exception :"<<merr<<es;
} ENDTRY;
#else
catch ( PException exc ) {
  fflush(stdout); cout<<endl; cerr<<endl;
  cerr<<"Services2NObjMgr::ComputeExpressions()  Exception :"<<exc.Msg()<<endl; 
}
#endif

if(adel) delete objnt;
if(fout) fclose(fout);
srvo.CloseDLL(); // Fermeture du fichier .so
return;
}


/* methode */
#define __LENLINE_N_READ__ 8192
void PAWExecutor::n_read(vector<string>& tokens)
{
 int lp=1;

 if(tokens.size()<3) {
   cerr<<"Usage: n/read nt fascii [options] var_1,c_1 ... var_n,c_n"<<endl;
   return;
 }

 // decodage des arguments
 string nament = tokens[0];
 string nameascii = tokens[1];
 bool create_datatable = false;
 vector<string> donotreadbeg;
 vector<string> donotreadin;
 vector<string> onlyreadbeg;
 vector<string> onlyreadin;
 vector<string> varname;
 vector<int> colnum;
 char separator = ' ';
 int numcolmaxi=-1;

 for(int i=2;i<(int_4)tokens.size();i++) {
   int lc = tokens[i].size();
   if(lc<2) continue;
   if(tokens[i].find("-dt")==0) { // create DataTable instead of NTuple
     create_datatable = true;
     continue;
   }
   const char *c = tokens[i].c_str();
   if(c[0]=='=') { // Separator
     separator = c[1];
     if(lc==3) if(c[1]=='\\' && c[2]=='t') separator = '\t';
     continue;
   }
   if(c[0]=='+') { // Selection des lignes a lire
     if(c[1]!='^') onlyreadin.push_back(&c[1]);
       else if(lc>2) onlyreadbeg.push_back(&c[2]);
     continue;
   }
   if(c[0]=='-') { // Selection des lignes commentaire
     if(c[1]!='^') donotreadin.push_back(&c[1]);
       else if(lc>2) donotreadbeg.push_back(&c[2]);
     continue;
   }
   // decodage des noms de variables et des colonnes associees
   int p = tokens[i].find(',');
   if(p<1 || p>=lc-1) continue;
   string vn = tokens[i].substr(0,p);
   string cn = tokens[i].substr(p+1,lc-p-1);
   int ic = atoi(cn.c_str());
   if(ic<0) continue;
   if( !isalpha(vn[0]) ) continue;
   if(ic>numcolmaxi) numcolmaxi = ic;
   varname.push_back(vn);
   colnum.push_back(ic);
 }

 int nvar = varname.size();
 if(nvar<=0) {
   cerr<<"n_read: no variables to be read"<<endl;
   return;
 }

 // Print what has to be done
 if(lp) {
   if(onlyreadin.size()>0) {
     cout<<"n_read Only read line containing ["<<onlyreadin.size()<<"]:";
     for(int i=0;i<(int_4)onlyreadin.size();i++) cout<<" \'"<<onlyreadin[i]<<"\'";
     cout<<endl;
   }
   if(onlyreadbeg.size()>0) {
     cout<<"n_read Only read line begining with ["<<onlyreadbeg.size()<<"]:";
     for(int i=0;i<(int_4)onlyreadbeg.size();i++) cout<<" \'"<<onlyreadbeg[i]<<"\'";
     cout<<endl;
   }
   if(donotreadin.size()>0) {
     cout<<"n_read Do not read line containing ["<<donotreadin.size()<<"]:";
     for(int i=0;i<(int_4)donotreadin.size();i++) cout<<" \'"<<donotreadin[i]<<"\'";
     cout<<endl;
   }
   if(donotreadbeg.size()>0) {
     cout<<"n_read Do not read line begining with ["<<donotreadbeg.size()<<"]:";
     for(int i=0;i<(int_4)donotreadbeg.size();i++) cout<<" \'"<<donotreadbeg[i]<<"\'";
     cout<<endl;
   }
   if(nvar>0) {
     cout<<"n_read Number of variables to be read: "<<nvar<<endl;
     for(int i=0;i<nvar;i++) cout<<" \'"<<varname[i].c_str()<<","<<colnum[i]<<"\'";
     cout<<endl;
   }
   if(separator!=' ') cout<<"n_read Separator is: \'"<<separator<<"\'"<<endl;
 }

 // Open ASCII file
 FILE * fascii = fopen(nameascii.c_str(),"r");
 if(fascii==NULL) {
   cerr<<"n_read: cannot open file "<<nameascii<<endl;
   return;
 }

 // Creation du NTuple ou de la DataTable
 NTuple *nt = NULL;
 DataTable *dt = NULL;
 if(!create_datatable) { //create NTuple
   char** ntvn = new char*[nvar];
   for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(varname[i].c_str());
   nt = new NTuple(nvar,ntvn);
   delete [] ntvn;
 } else { //create DataTable
   dt = new DataTable();
   for(int i=0;i<nvar;i++) dt->AddDoubleColumn(varname[i]);
 }
 r_8 *xnt = new r_8[nvar];

 // Read file
 char *line = new char[__LENLINE_N_READ__];
 int nline=0, nlinecom=0;
 while(fgets(line,__LENLINE_N_READ__,fascii) != NULL ) {
   nline++;
   uint_4 lc = strlen(line); if(lc<1) continue;
   // Pour enlever le \n final
   if(line[lc-1]=='\n' || line[lc-1]=='\r')
     {line[lc-1]='\0'; lc = strlen(line); if(lc<1) continue;}

   string const sline(line);
   //cout<<"\'"<<sline<<"\'  lc="<<lc<<endl;

   // Faut t'il lire cette ligne ?
   bool read_line_1 = true;
   if(onlyreadin.size()>0 || onlyreadbeg.size()>0) read_line_1 = false;
   if(onlyreadin.size()>0) {
     for(int i=0;i<(int_4)onlyreadin.size();i++) {
       uint_4 p = sline.find(onlyreadin[i].c_str());
       if(p>=lc) continue;
       read_line_1 = true;
       break;
     }
   }
   if(onlyreadbeg.size()>0) {
     for(int i=0;i<(int_4)onlyreadbeg.size();i++) {
       uint_4 p = sline.find(onlyreadbeg[i].c_str());
       if(p!=0) continue;
       read_line_1 = true;
       break;
     }
   }

   // Faut t'il ne pas lire cette ligne ?
   bool read_line_2 = true;
   if(donotreadin.size()>0) {
     for(int i=0;i<(int_4)donotreadin.size();i++) {
       uint_4 p = sline.find(donotreadin[i].c_str());
       if(p>=lc) continue;
       read_line_2 = false;
       break;
     }
   }
   if(donotreadbeg.size()>0) {
     for(int i=0;i<(int_4)donotreadbeg.size();i++) {
       uint_4 p = sline.find(donotreadbeg[i].c_str());
       if(p!=0) continue;
       read_line_2 = false;
       break;
     }
   }
   if(!read_line_2) nlinecom++;

   if(!read_line_1 || !read_line_2) continue;

   // Decodage de la ligne
   vector<string> vs;
   //FillVStringFrString(sline,vs,separator);
   SplitStringToVString(sline,vs,separator);
   int lvs = vs.size();
   //for(int i=0;i<lvs;i++) cout<<"|"<<vs[i]<<"| "; cout<<endl;
   if(lvs<numcolmaxi) continue; // Pas assez de champs decodes, mauvaise ligne

   // Remplissage du NTuple ou de la DataTable
   for(int i=0;i<nvar;i++) {
     xnt[i] = 0.;
     int ic = colnum[i];
     if(ic>=lvs) continue;
     xnt[i] = atof(vs[ic].c_str());
   }
   if(nt)      nt->Fill(xnt);
   else if(dt) dt->Fill(xnt);
   //cout<<"...xnt"; for(int i=0;i<nvar;i++) cout<<" "<<xnt[i]; cout<<endl;

 }
 int_4 nentries = (nt)? nt->NEntry(): dt->NEntry();
 cout<<"n_read: "<<nline<<" lines in file, "
     <<nlinecom<<" commentary, "
     <<nentries<<" object entries"<<endl;

 // On sauve le NTuple/DataTable si besoin, on ferme et detruit ce qu'il faut
 NamedObjMgr omg;
 if(nentries>0) {
   if(nt) omg.AddObj(nt,nament);
   if(dt) omg.AddObj(dt,nament);
 } else {
   if(nt) delete nt;
   if(dt) delete dt;
 }
 delete [] xnt;
 delete [] line;
 fclose(fascii);

 return;
}
#undef __LENLINE_N_READ__

/* methode */
void PAWExecutor::n_merge(vector<string>& tokens)
{
 if(tokens.size()<2) {
   cerr<<"Usage: n/merge nt nt_1 nt_2 ... nt_n"<<endl;
   return;
 }

 NamedObjMgr omg;

 // decodage des arguments
 string nament = tokens[0];

 // boucle sur les ntuples
 NTuple * nt = NULL;
 int nvar=0, nfill=0;
 r_8 *xnt=NULL;

 for(int i=1;i<(int_4)tokens.size();i++) {

   AnyDataObj* mobj = omg.GetObj(tokens[i]);
   if(mobj==NULL) {
     cout<<"n_merge Error: unknow object"<<tokens[i]<<endl;
     continue;
   }
   NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
   if(nt1==NULL) {
     cout<<"n_merge Error: "<<tokens[i]<<" not a NTuple"<<endl;
     continue;
   }
   if(nt1->NEntry()==0) {
     cout<<"n_merge Error: "<<tokens[i]<<" is empty"<<endl;
     continue;
   }
   if(nt1->NVar()==0) {
     cout<<"n_merge Error: "<<tokens[i]<<" has no variable"<<endl;
     continue;
   }

   // create receiving ntuple if first pass
   if(nt==NULL) {
     nvar = nt1->NVar();
     vector<string> sntvn;
     for(int i=0;i<nvar;i++) sntvn.push_back(nt1->ColumnName(i));
     char **ntvn = new char*[nvar];
     for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(sntvn[i].c_str());
     nt = new NTuple(nvar,ntvn);
     delete [] ntvn;
   }

   // filling with current ntuple
   int nvar1 = nt1->NVar();
   int n = (nvar1>nvar)? nvar1: nvar;
   r_8 *xnt1 = new r_8[n];
   for(int i=0;i<n;i++) xnt1[i]=0.;
   for(uint_4 iev=0;iev<(uint_4)nt1->NEntry();iev++) {
     nt1->GetVecD(iev,xnt1);
     nt->Fill(xnt1);
   }
   nfill++;
   delete [] xnt1;

 }

 if(xnt!=NULL) delete [] xnt;
 if(nt!=NULL) {
   cout<<"n_merge: ntuple filled with "<<nfill
       <<" ntuples, "<<nt->NEntry()<<" entries"<<endl;
   if(nt->NEntry()>0) omg.AddObj(nt,nament); else delete nt;
 }
 return;
}

/* methode */
void PAWExecutor::n_merge_col(vector<string>& tokens)
{
 if(tokens.size()<2) {
   cerr<<"Usage: n/merge/col nt nt_1[,ext] nt_2[,ext] ... nt_n[,ext]"<<endl;
   return;
 }

 NamedObjMgr omg;

 // decodage des arguments
 string nament = tokens[0];

 // boucle sur les ntuples
 uint_4 numnt=0, nentmax=0, ipos=0;
 vector<string> varname;
 vector<string> ntname;
 vector<int> posfirst;
 for(int i=1;i<(int_4)tokens.size();i++) {

   // Decode NTuple name and extension
   int_4 lc = tokens[i].size();
   string ntn = tokens[i];
   string ext = "";
   int p = ntn.find(',');
   if(p>=0 && p<lc) {
     if(p<lc-1) ext = ntn.substr(p+1,lc-p-1);
     if(p>0) ntn = ntn.substr(0,p); else ntn = "";
   }
   // choose the extentsion: given,automatic,none
   if(ext.size()<1) {
     char str[8]; sprintf(str,"_%d",numnt);
     ext = str;
   } else if(ext[0] == '!') {
     ext = "";
   }

   AnyDataObj* mobj = omg.GetObj(ntn);
   if(mobj==NULL) {
     cout<<"n_merge_col Error: unknow object "<<ntn<<" ("<<tokens[i]<<")"<<endl;
     continue;
   }
   NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
   if(nt1==NULL) {
     cout<<"n_merge_col Error: "<<ntn<<" not a NTuple"<<endl;
     continue;
   }
   if(nt1->NEntry()==0) {
     cout<<"n_merge_col Error: "<<ntn<<" is empty"<<endl;
     continue;
   }
   if(nt1->NVar()==0) {
     cout<<"n_merge_col Error: "<<ntn<<" has no variable"<<endl;
     continue;
   }

   ntname.push_back(ntn);
   posfirst.push_back(ipos);
   if(nt1->NEntry()>(int_4)nentmax) nentmax = nt1->NEntry();
   for(int i=0;i<nt1->NVar();i++) {
     string str = nt1->ColumnName(i) + ext;
     varname.push_back(str);
   }
   numnt++;
   ipos += nt1->NVar();
 }

 if(ntname.size()==0) {
   cout<<"n_merge_col Error: no ntuple found in argument list"<<endl;
   return;
 } else cout<<"Number of NTuple found "<<ntname.size()
            <<", Number of variables "<<varname.size()
            <<", Max_entries "<<nentmax<<endl;

 // create receiving ntuple
 int nvar = varname.size();
 char **ntvn = new char*[nvar];
 for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(varname[i].c_str());
 NTuple * nt = new NTuple(nvar,ntvn);
 //for(int i=0;i<(int_4)ntname.size();i++) cout<<ntname[i]<<" , p="<<posfirst[i]<<endl;
 delete [] ntvn;

 // filling with current ntuple
 r_8 *xnt = new r_8[nvar];
 for(int iev=0;iev<(int_4)nentmax;iev++) {
   for(int i=0;i<nvar;i++) xnt[i]=0.;
   for(int i=0;i<(int_4)ntname.size();i++) {
     AnyDataObj* mobj = omg.GetObj(ntname[i]);
     NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
     if(iev < nt1->NEntry())  nt1->GetVecD(iev,&xnt[posfirst[i]]);
   }
   nt->Fill(xnt);
 }

 // Ending and saving merging ntuple columns
 delete [] xnt;
 cout<<"n_merge_col: ntuple filled with "<<ntname.size()
     <<" ntuples, "<<nt->NEntry()<<" entries"<<endl;
 if(nt->NEntry()>0) omg.AddObj(nt,nament); else delete nt;

 return;
}

/* methode */
void PAWExecutor::h_integ(vector<string>& tokens)
// Pour remplacer le contenu d'un histo 1D par son integrale
{
if(tokens.size()<1)
  {cout<<"Usage: h/integ nameh1d [norm]"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_integ Error: unknow object"<<tokens[0]<<endl;
  return;}
r_8 norm = 1.;
if(tokens.size()>=2) norm = atof(tokens[1].c_str());
// attention: dynamic_cast<Histo*>(HProf)=Vrai!
Histo* h1 = dynamic_cast<Histo*>(mobj);
HProf* hp = dynamic_cast<HProf*>(mobj);
if(hp || !h1)
  {cout<<"PAWExecutor::h_integ Error: "<<tokens[0]<<" not an Histo"<<endl;
  return;}
h1->HInteg(norm);
}

/* methode */
void PAWExecutor::v_integ(vector<string>& tokens)
// Pour remplacer le contenu d'un TVector<r_8> par son integrale
//   normalisee a norm:
// Si norm <= 0 : pas de normalisation, integration seule
{
if(tokens.size()<1)
  {cout<<"Usage: v/integ namevec [norm]"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::v_integ Error: unknow object"<<tokens[0]<<endl;
  return;}
r_8 norm=-1.;
if(tokens.size()>=2) norm = atof(tokens[1].c_str());
#ifdef SANS_EVOLPLANCK
Vector* v = dynamic_cast<Vector*>(mobj);
#else
TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
#endif

if(!v)
  {cout<<"PAWExecutor::v_integ Error: "<<tokens[0]<<" not a TVector/Vector"<<endl;
  return;}

#ifdef SANS_EVOLPLANCK
uint_4 n = v->NElts();
#else
uint_4 n = v->Size();
#endif

if(n==0)
  {cout<<"PAWExecutor::v_integ Error: "<<tokens[0]<<" is an empty vector"<<endl;
  return;}
if(n>1) for(uint_4 i=1;i<n;i++) (*v)(i)+=(*v)(i-1);
if(norm<=0. || (*v)(n-1)==0.) return;
norm /= (*v)(n-1);
for(uint_4 i=0;i<n;i++) (*v)(i) *= norm;
}

/* methode */
void PAWExecutor::v_sort(vector<string>& tokens)
// Pour sort in-place d'un vecteur ascendant ou descendant
{
 if(tokens.size()<1)
   {cout<<"Usage: v/sort namevec [+1/-1]"<<endl; return;}
 NamedObjMgr omg;
 AnyDataObj* mobj = omg.GetObj(tokens[0]);
 if(mobj==NULL)
   {cout<<"PAWExecutor::v_sort Error: unknow object"<<tokens[0]<<endl;
   return;}

 bool up = true;
 if(tokens.size()>=2) if(atoi(tokens[1].c_str())<0) up=false;

 {
 TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
 if(v) {
   if(v->Size()==0) return;
   if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_r8);
     else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_r8);
   cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
   return;
 }
 }

 {
 TVector<r_4>* v = dynamic_cast<TVector<r_4>*>(mobj);
 if(v) {
   if(v->Size()==0) return;
   if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_r4);
     else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_r4);
   cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
   return;
 }
 }

 {
 TVector<int_4>* v = dynamic_cast<TVector<int_4>*>(mobj);
 if(v) {
   if(v->Size()==0) return;
   if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_i4);
     else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_i4);
   cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
   return;
 }
 }

 {
 TVector<uint_8>* v = dynamic_cast<TVector<uint_8>*>(mobj);
 if(v) {
   if(v->Size()==0) return;
   if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_u8);
     else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_u8);
   cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
   return;
 }
 }

 cout<<"PAWExecutor::v_sort Error: "<<tokens[0]
     <<" not a TVector or not a supported TVector<TYPE>"<<endl;
 return;

}

/* methode */
void PAWExecutor::h_deriv(vector<string>& tokens)
// Pour remplacer le contenu d'un histo 1D par sa derivee
{
if(tokens.size()<1)
  {cout<<"Usage: h/deriv nameh1d"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_deriv Error: unknow object"<<tokens[0]<<endl;
  return;}
// attention: dynamic_cast<Histo*>(HProf)=Vrai!
Histo* h1 = dynamic_cast<Histo*>(mobj);
HProf* hp = dynamic_cast<HProf*>(mobj);
if(hp || !h1)
  {cout<<"PAWExecutor::h_deriv Error: "<<tokens[0]<<" not an Histo"<<endl;
  return;}
h1->HDeriv();
}

/* methode */
void PAWExecutor::v_deriv(vector<string>& tokens)
// Pour remplacer le contenu d'un TVector<r_8> par sa derivee
// deriv_option = -1 replace v[i] with v[i]-v[i-1]
//              =  0 replace v[i] with (v[i+1]-v[i-1])/2 (default)
//              = +1 replace v[i] with v[i+1]-v[i]
{
if(tokens.size()<1)
  {cout<<"Usage: v/deriv namevec [deriv_option(-1,0,+1)]"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::v_deriv Error: unknow object"<<tokens[0]<<endl;
  return;}
int_4 deriv_option = 0;
if(tokens.size()>=2) deriv_option = atoi(tokens[1].c_str());

#ifdef SANS_EVOLPLANCK
Vector* v = dynamic_cast<Vector*>(mobj);
#else
TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
#endif

if(!v)
  {cout<<"PAWExecutor::v_deriv Error: "<<tokens[0]<<" not a TVector/Vector"<<endl;
  return;}


#ifdef SANS_EVOLPLANCK
uint_4 n = v->NElts();
#else
uint_4 n = v->Size();
#endif

if(n==0)
  {cout<<"PAWExecutor::v_deriv Error: "<<tokens[0]<<" is an empty vector"<<endl;
  return;}
if(n<=1) return;

#ifdef SANS_EVOLPLANCK
Vector vsave(*v);
#else
TVector<r_8> vsave(*v,false);
#endif

if(deriv_option<0) {
  for(uint_4 i=1;i<n;i++) (*v)(i) = vsave(i)-vsave(i-1);
  (*v)(0) = (*v)(1);
} else if(deriv_option>0) {
  for(uint_4 i=0;i<n-1;i++) (*v)(i) = vsave(i+1)-vsave(i);
  (*v)(n-1) = (*v)(n-2);
} else {
  for(uint_4 i=1;i<n-1;i++) (*v)(i) = (vsave(i+1)-vsave(i-1))/2.;
  (*v)(0)   = vsave(1)-vsave(0);
  (*v)(n-1) = vsave(n-1)-vsave(n-2);
}
}


/* methode */
void PAWExecutor::h_rebin(vector<string>& tokens)
// Pour re-binner un histogramme 1D
{
if(tokens.size()<2)
  {cout<<"Usage: h/rebin nameh1d nbin"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_rebin Error: unknow object"<<tokens[0]<<endl;
  return;}
int_4 nbin = atoi(tokens[1].c_str());
Histo* h1 = dynamic_cast<Histo*>(mobj);
// Ca marche aussi pour les HProf, HRebin a ete passe virtuel
//HProf* hp = dynamic_cast<HProf*>(mobj); if(hp || !h1)
if(!h1)
  {cout<<"PAWExecutor::h_rebin Error: "<<tokens[0]<<" not an Histo/HProf"<<endl;
  return;}
h1->HRebin(nbin);
}

/* methode */
void PAWExecutor::h_cadd(vector<string>& tokens)
// Additionne une constante a un histogramme
{
if(tokens.size()<2)
  {cout<<"Usage: h/cadd nameh1d val"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_cadd Error: unknow object"<<tokens[0]<<endl;
  return;}
r_8 val = atof(tokens[1].c_str());
Histo*   h1 = dynamic_cast<Histo*>(mobj);
HProf*   hp = dynamic_cast<HProf*>(mobj);
Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
Vector*  v  = dynamic_cast<Vector*>(mobj);
Matrix*  m  = dynamic_cast<Matrix*>(mobj);
if(h1 && !hp) *h1 += val;
else if(h2)   *h2 += val;
else if(v)    *v  += val;
else if(m)    *m  += val;
else cout<<"PAWExecutor::h_cadd Error: not implemented for "<<typeid(*mobj).name()<<endl;
}

/* methode */
void PAWExecutor::h_cmult(vector<string>& tokens)
// Multiplie un histogramme par une constante
{
if(tokens.size()<2)
  {cout<<"Usage: h/cmult nameh1d val"<<endl; return;}
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_cmult Error: unknow object"<<tokens[0]<<endl;
  return;}
r_8 val = atof(tokens[1].c_str());
Histo*   h1 = dynamic_cast<Histo*>(mobj);
HProf*   hp = dynamic_cast<HProf*>(mobj);
Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
Vector*  v  = dynamic_cast<Vector*>(mobj);
Matrix*  m  = dynamic_cast<Matrix*>(mobj);
if(h1 && !hp) *h1 *= val;
else if(h2)   *h2 *= val;
else if(v)    *v  += val;
else if(m)    *m  += val;
else cout<<"PAWExecutor::h_mult Error: not implemented for "<<typeid(*mobj).name()<<endl;
}

/* methode */
void PAWExecutor::h_oper(vector<string>& tokens)
// Equivalent h/oper/add sub,mul,div de paw
// Operation entre 2 histogrammes ou 2 vecteurs ou 2 matrices
// h/oper @ h1 h2 hres
// hres = h1 @ h2    with @ = (+,-,*,/)
// Pour les vecteurs et les matrices, il sagit d'operations elements a elements
{
if(tokens.size()<4)
  {cout<<"Usage: n/oper @ h1 h2 hres with @=(+,-,*,/,@)"<<endl;
   return;}

// Decode arguments
const char * oper = tokens[0].c_str();
if( oper[0]!='+' && oper[0]!='-' && oper[0]!='*' && oper[0]!='/' )
  {cout<<"PAWExecutor::h_oper Error: unknow operation "<<oper<<endl;
  return;}
string h1name = tokens[1];
string h2name = tokens[2];
string h3name = tokens[3];

// Get objects
NamedObjMgr omg;
AnyDataObj* mobjh1 = omg.GetObj(h1name);
AnyDataObj* mobjh2 = omg.GetObj(h2name);
if( mobjh1==NULL || mobjh2==NULL )
  {cout<<"PAWExecutor::h_oper Error: unknow object(s) "<<h1name<<" or "<<h2name<<endl;
  return;}
AnyDataObj* mobjh3 = omg.GetObj(h3name);

// Operations on HProf, only + is working
if( dynamic_cast<HProf*>(mobjh1) != NULL ) {
  if( oper[0]!='+' )
    { cout<<"PAWExecutor::h_oper Error: operation "<<oper
          <<" not implemented for HProf"<<endl; return;}
  if( dynamic_cast<HProf*>(mobjh2) == NULL )
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
    return;}
  HProf* h1 =(HProf*) mobjh1;
  HProf* h2 =(HProf*) mobjh2;
  if( h1->NBins() != h2->NBins() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
         <<h1->NBins()<<" "<<h2->NBins()<<endl; return;}
  HProf* h3 = NULL;
  if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
    {h3 = new HProf(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
  h3 = dynamic_cast<HProf*>(mobjh3);
  if(h3 == NULL)  // ce n'est pas un HProf
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
  if(h1->NBins() != h3->NBins())
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
         <<h1->NBins()<<" "<<h3->NBins()<<endl; return;}
  *h3 = *h1 + *h2;
  h3->UpdateHisto();

// Operations on Histo
} else if( dynamic_cast<Histo*>(mobjh1) != NULL ) {
  if( dynamic_cast<Histo*>(mobjh2) == NULL )
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
    return;}
  Histo* h1 =(Histo*) mobjh1;
  Histo* h2 =(Histo*) mobjh2;
  if( h1->NBins() != h2->NBins() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
         <<h1->NBins()<<" "<<h2->NBins()<<endl; return;}
  Histo* h3 = NULL;
  if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
    {h3 = new Histo(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
  h3 = dynamic_cast<Histo*>(mobjh3);
  if(h3 == NULL)  // ce n'est pas un Histo
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
  if(h1->NBins() != h3->NBins())
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
         <<h1->NBins()<<" "<<h3->NBins()<<endl; return;}
  if( oper[0]=='+')      *h3 = *h1 + *h2;
  else if( oper[0]=='-') *h3 = *h1 - *h2;
  else if( oper[0]=='*') *h3 = *h1 * *h2;
  else if( oper[0]=='/') *h3 = *h1 / *h2;

// Operations on Histo2D
} else if( dynamic_cast<Histo2D*>(mobjh1) != NULL ) {
  if( dynamic_cast<Histo2D*>(mobjh2) == NULL )
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
    return;}
  Histo2D* h1 =(Histo2D*) mobjh1;
  Histo2D* h2 =(Histo2D*) mobjh2;
  if( h1->NBinX() != h2->NBinX() || h1->NBinY() != h2->NBinY() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
         <<h1->NBinX()<<","<<h1->NBinY()<<"   "
         <<h2->NBinX()<<","<<h2->NBinY()<<endl; return;}
  Histo2D* h3 = NULL;
  if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
    {h3 = new Histo2D(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
  h3 = dynamic_cast<Histo2D*>(mobjh3);
  if(h3 == NULL)  // ce n'est pas un Histo2D
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
  if( h1->NBinX() != h3->NBinX() || h1->NBinY() != h3->NBinY() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
         <<h1->NBinX()<<","<<h1->NBinY()<<"   "
         <<h3->NBinX()<<","<<h3->NBinY()<<endl; return;}
  if( oper[0]=='+')      *h3 = *h1 + *h2;
  else if( oper[0]=='-') *h3 = *h1 - *h2;
  else if( oper[0]=='*') *h3 = *h1 * *h2;
  else if( oper[0]=='/') *h3 = *h1 / *h2;

// Operations on Vector
} else if( dynamic_cast<Vector*>(mobjh1) != NULL ) {
  if( dynamic_cast<Vector*>(mobjh2) == NULL )
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
    return;}
  Vector* h1 =(Vector*) mobjh1;
  Vector* h2 =(Vector*) mobjh2;
  if( h1->NElts() != h2->NElts() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
         <<h1->NElts()<<" "<<h2->NElts()<<endl; return;}
  Vector* h3 = NULL;
  if( mobjh3 == NULL ) {  // l'objet n'existe pas, on le cree
#ifdef SANS_EVOLPLANCK
    h3 = new Vector(*h1);
#else
    h3 = new Vector(*h1,false);
#endif
    *h3 = 0.; omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);
  }
  h3 = dynamic_cast<Vector*>(mobjh3);
  if(h3 == NULL)  // ce n'est pas un Vector
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
  if(h1->NElts() != h3->NElts())
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
         <<h1->NElts()<<" "<<h3->NElts()<<endl; return;}
  if( oper[0]=='+')      *h3 = *h1 + *h2;
  else if( oper[0]=='-') *h3 = *h1 - *h2;
#ifdef SANS_EVOLPLANCK
  else if(oper[0]=='*' || oper[0]=='/')
    cout<<"PAWExecutor::h_oper Error: operation "<<oper[0]
        <<" not implemented for Vector in Peida"<<endl;
#else
  else if( oper[0]=='*') {h3->Clone(*h1); h3->MulElt(*h2,*h3);}
  else if( oper[0]=='/') {h3->Clone(*h1); h3->DivElt(*h2,*h3,false,true);}
#endif

// Operations on Matrix
} else if( dynamic_cast<Matrix*>(mobjh1) != NULL ) {
  if( dynamic_cast<Matrix*>(mobjh2) == NULL )
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
    return;}
  Matrix* h1 =(Matrix*) mobjh1;
  Matrix* h2 =(Matrix*) mobjh2;
  if( h1->NRows() != h2->NRows() || h1->NCol() != h2->NCol() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
         <<h1->NRows()<<","<<h1->NCol()<<"   "
         <<h2->NRows()<<","<<h2->NCol()<<endl; return;}
  Matrix* h3 = NULL;
  if( mobjh3 == NULL ) {  // l'objet n'existe pas, on le cree
#ifdef SANS_EVOLPLANCK
    h3 = new Matrix(*h1);
#else
    h3 = new Matrix(*h1,false);
#endif
    *h3 = 0.; omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);
  }
  h3 = dynamic_cast<Matrix*>(mobjh3);
  if(h3 == NULL)  // ce n'est pas un Matrix
    {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
         <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
  if( h1->NRows() != h3->NRows() || h1->NCol() != h3->NCol() )
    {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
         <<h1->NRows()<<","<<h1->NCol()<<"   "
         <<h3->NRows()<<","<<h3->NCol()<<endl; return;}
  if( oper[0]=='+')      *h3 = *h1 + *h2;
  else if( oper[0]=='-') *h3 = *h1 - *h2;
#ifdef SANS_EVOLPLANCK
  else if(oper[0]=='*' || oper[0]=='/')
    cout<<"PAWExecutor::h_oper Error: operation "<<oper[0]
        <<" not implemented for Vector in Peida"<<endl;
#else
  else if( oper[0]=='*') {h3->Clone(*h1); h3->MulElt(*h2,*h3);}
  else if( oper[0]=='/') {h3->Clone(*h1); h3->DivElt(*h2,*h3,false,true);}
#endif

// Doesn't work for other objects
} else {
  cout<<"PAWExecutor::h_oper Error: not implemented for "
      <<typeid(*mobjh1).name()<<endl;
  return;
}

return;
}

/* methode */
void PAWExecutor::h_plot_2d(vector<string>& tokens)
// plot for 2D histogramme: plot histo, bandx/y, slicex/y or projx/y
{
if(tokens.size()<1)
  {cout<<"Usage: h/plot/2d nameh2d to_plot [n/s] [dopt]"<<endl; return;}
string proj = "h"; if(tokens.size()>1) proj = tokens[1];
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if(mobj==NULL)
  {cout<<"PAWExecutor::h_plot_2d Error: unknow object"<<tokens[0]<<endl;
  return;}
Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
if(!h2)
  {cout<<"PAWExecutor::h_plot_2d Error: "<<tokens[0]<<" not an Histo2D"<<endl;
  return;}

Histo* h1p = NULL; string nametoplot = "/autoc/h_plot_2d_h1";

string dopt = ""; if(tokens.size()>=3) dopt = tokens[2];
if(proj == "show") {
  h2->ShowProj();
  h2->ShowBand(2);
  h2->ShowSli(2);
  return;
} else if(proj == "h") {
  nametoplot = tokens[0];
} else if(proj == "px") {  
  if((h1p=h2->HProjX())) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowProj(); return;}
} else if(proj == "py") {
  if((h1p=h2->HProjY())) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowProj(); return;}
} else {
  if(tokens.size()<3)
    {cout<<"Usage: h/plot/2d nameh2d bx/by/sx/sy n [dopt]"<<endl; return;}
  int_4 n = atoi(tokens[2].c_str());
  dopt = ""; if(tokens.size()>=4) dopt = tokens[3];
  if(proj == "bx") {
    if((h1p=h2->HBandX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowBand(); return;}
  } else if(proj == "by") {
    if((h1p=h2->HBandY(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowBand(); return;}
  } else if(proj == "sx") {
    if((h1p=h2->HSliX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowSli(); return;}
  } else if(proj == "sy") {
    if((h1p=h2->HSliY(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowSli(); return;}
  }
}

omg.DisplayObj(nametoplot,dopt);
}

/* methode */
int_4 PAWExecutor::decodepawstring(string tokens,string& nameobj
                                ,string& xexp,string& yexp,string& zexp)
// Decodage general de "nameobj.xexp"
//                     "nameobj.yexp%xexp"
//                     "nameobj.zexp%yexp%xexp"
// Return: nombre de variables trouvees, -1 si probleme
{
nameobj = ""; xexp= ""; yexp= ""; zexp= "";

int_4 lt = (int) tokens.length();
if(lt<=0) return -1;

// decodage de la chaine de type PAW.
char *str = new char[lt+2];
strcpy(str,tokens.c_str()); strip(str,'B',' '); lt = strlen(str);
//cout<<"chaine1["<<lt<<"] :"<<str<<":"<<endl;
char *c[3] = {NULL,NULL,NULL};
int_4 i, np=0; bool namefound = false;
for(i=0;i<lt;i++) {
  if(!namefound && str[i]=='.') {
    str[i]='\0';
    namefound=true;
    c[np] = str+i+1; np++;
  }
  if( namefound && str[i]=='%') {
    str[i]='\0';
    if(np<3) {c[np] = str+i+1; np++;}
  }
}
//cout<<"chaine2 :"; for(i=0;i<lt;i++) cout<<str[i]; cout<<":"<<endl;

// Remplissage du nom et des variables
nameobj = str;
if(np==1)  xexp=c[0];
if(np==2) {yexp=c[0]; xexp=c[1];}
if(np==3) {zexp=c[0]; yexp=c[1]; xexp=c[2];}
//cout<<"pawstring str,c[0-2] "<<str<<" "<<c[0]<<" "<<c[1]<<" "<<c[2]<<endl;
delete [] str;

// Comptage des variables
np = -1;
if(nameobj.length()>0) 
  {np = 0; if(xexp.length()>0)
             {np++; if(yexp.length()>0)
                      {np++; if(zexp.length()>0) np++;}}}
//cout<<"pawstring["<<np<<"] name="<<nameobj
//    <<" xexp="<<xexp<<" yexp="<<yexp<<" zexp="<<zexp<<endl;
return np;
}

/* methode */
void PAWExecutor::h_put_vec(vector<string>& tokens)
// Pour remplir un histo avec le contenu d'un vecteur ou d'une matrice
// L'histogramme doit deja exister. Le nombre de bins remplit
// depend des tailles respectives des 2 objets.
// tokens[2] = "cont" : on remplit les valeurs de l'histogramme (ou "!")
//           = "err2" : on remplit les erreurs de l'histogramme
{
if(tokens.size()<2)
  {cout<<"Usage: h/put_vec  namehisto namevector"<<endl;
   return;}
string toput = "cont"; if(tokens.size()>2) toput = tokens[2];
if(toput=="!") toput = "cont";
if(toput!="cont" && toput!="err2")
  {cout<<"PAWExecutor::h_put_vec Error: unknow filling "<<toput<<endl;
   return;}
string hname = tokens[0];
string vname = tokens[1];

// Get objects
NamedObjMgr omg;
AnyDataObj* mobjh = omg.GetObj(hname);
AnyDataObj* mobjv = omg.GetObj(vname);
if( mobjh==NULL || mobjv==NULL )
  {cout<<"PAWExecutor::h_put_vec Error: unknow object(s) "<<hname<<" or "<<vname<<endl;
  return;}

// Fill Histo from Vector
if(typeid(*mobjh) == typeid(Histo) && typeid(*mobjv) == typeid(Vector)) {
  Histo*  h = dynamic_cast<Histo*>(mobjh);
  Vector* v = dynamic_cast<Vector*>(mobjv);
  if(toput=="cont")      h->PutValue(*v);
  else if(toput=="err2") h->PutError2(*v);

// Fill Histo2D from Matrix
} else if(typeid(*mobjh) == typeid(Histo2D) && typeid(*mobjv) == typeid(Matrix)) {
  Histo2D* h = dynamic_cast<Histo2D*>(mobjh);
  Matrix*  v = dynamic_cast<Matrix*>(mobjv);
  if(toput=="cont")      h->PutValue(*v);
  else if(toput=="err2") h->PutError2(*v);

} else {
  cout<<"PAWExecutor::h_put_vec Error: type mismatch between histogram and vector/matrix\n"
      <<typeid(*mobjh).name()<<" , "<<typeid(*mobjv).name()<<endl;
  return;
}

}

/* methode */
void PAWExecutor::h_get_vec(vector<string>& tokens)
// Pour copier un histo dans un vecteur ou une matrice
// Si le vecteur (matrice) n'existe pas, il est cree.
// Le vecteur ou la matrice est re-dimensionne correctement.
// tokens[2] = "cont" : on copie les valeurs de l'histogramme (ou "!")
//           = "err2" : on copie les erreurs de l'histogramme
//           = "absc" : on copie les erreurs de l'histogramme (1D only)
// tokens[3[,4]] = show : show available projections for Histo2D
//                 px   : get X projection
//                 py   : get Y projection
//                 bx n : get X band number n
//                 by n : get Y band number n
//                 sx n : get X slice number n
//                 sy n : get Y slice number n
{
if(tokens.size()<2)
  {cout<<"Usage: h/get_vec namehisto namevector"<<endl;
   return;}
string toget = "cont"; if(tokens.size()>2) toget = tokens[2];
if(toget=="!") toget = "cont";
if(toget!="cont" && toget!="err2" && toget!="absc")
  {cout<<"PAWExecutor::h_get_vec Error: unknow filling "<<toget<<endl;
   return;}
string proj = "h";  if(tokens.size()>3) proj = tokens[3];
int_4 nproj = -1;   if(tokens.size()>4) nproj = atoi(tokens[4].c_str());
string hname = tokens[0];
string vname = tokens[1];

// Get objects
NamedObjMgr omg;
AnyDataObj* mobjh = omg.GetObj(hname);
AnyDataObj* mobjv = omg.GetObj(vname);
if( mobjh==NULL)
  {cout<<"PAWExecutor::h_put_vec Error: unknow object(s) "<<hname<<endl;
  return;}

// Copy Histo/Histo2D/Projection to Vector/Matrix
Histo*   h  = dynamic_cast<Histo*>(mobjh);
Histo2D* h2 = dynamic_cast<Histo2D*>(mobjh);
if( h != NULL ) {  // Histo ou HProf
  h->UpdateHisto();  // pour le cas ou c'est un HProf
  Vector* v = NULL;
  if(mobjv==NULL)   // le vecteur n'existe pas
    {v = new Vector(1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
  if(typeid(*mobjv) != typeid(Vector))
    {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo/HProf and vector\n"
         <<typeid(*mobjv).name()<<endl; return;}
  v = dynamic_cast<Vector*>(mobjv);
  if(toget=="cont")      h->GetValue(*v);
  else if(toget=="err2") h->GetError2(*v);
  else if(toget=="absc") h->GetAbsc(*v);

} else if( h2 != NULL) {  // Histo2D

  if(proj == "h") {   // On veut les valeurs de l'histogramme 2D
    Matrix*  v = NULL;
    if(mobjv==NULL)   // la matrice n'existe pas
      {v = new Matrix(1,1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
    if(typeid(*mobjv) != typeid(Matrix))
      {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo2D and matrix\n"
           <<typeid(*mobjv).name()<<endl; return;}
    v = dynamic_cast<Matrix*>(mobjv);
    if(toget=="cont")      h2->GetValue(*v);
    else if(toget=="err2") h2->GetError2(*v);
    else
      {cout<<"PAWExecutor::h_get_vec Error: option "<<toget<<" not valid for Histo2D"<<endl;
      return;}

  } else {   // On veut les valeurs d'une projection de l'histo 2D
    h = NULL;
    if(proj == "show") {h2->ShowProj(); h2->ShowBand(2); h2->ShowSli(2);}
    else if(proj == "px") h = h2->HProjX();
    else if(proj == "py") h = h2->HProjY();
    else if(proj == "bx") h = h2->HBandX(nproj);
    else if(proj == "by") h = h2->HBandY(nproj);
    else if(proj == "sx") h = h2->HSliX(nproj);
    else if(proj == "sy") h = h2->HSliY(nproj);
    if(h==NULL)
      {cout<<"PAWExecutor::h_get_vec Error: unknown projection "<<proj
           <<" number "<<nproj<<endl; return;}
    Vector* v = NULL;
    if(mobjv==NULL)   // le vecteur n'existe pas
      {v = new Vector(1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
    if(typeid(*mobjv) != typeid(Vector))
      {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo2D "<<proj
       <<" and vector\n"<<typeid(*mobjv).name()<<endl; return;}
    v = dynamic_cast<Vector*>(mobjv);
    if(toget=="cont")      h->GetValue(*v);
    else if(toget=="err2") h->GetError2(*v);
    else if(toget=="absc") h->GetAbsc(*v);
  }

} else {
  cout<<"PAWExecutor::h_get_vec Error: type mismatch for histogram\n"
      <<typeid(*mobjh).name()<<endl;
  return;
}

}

/* methode */
void PAWExecutor::h_copy(vector<string>& tokens)
// Pour copier un object dans un object
// objects are Vector,Matrix,Histo,Histo2D
//   h/copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]
//   copy obj1Dfrom(i1->i2) into obj1Dto(ic1->)
//   copy obj2Dfrom(i1,j1->i2,j2) into obj2Dto(ic1,jc1->)
// Attention: elements depuis i1 jusqu'a i2 COMPRIS
// Si obj1Dto n'existe pas, il est cree avec une taille i2-i1+1
//    ou obj2Dto avec une taille i2-i1+1,j2-j1+1
// Le contenu de obj?Dfrom adresse est surecrit
// Le contenu de obj?Dfrom non adresse reste le meme
{
int_4 tks = tokens.size();
if(tks<2)
  {cout<<"Usage: h_copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]"<<endl;
   return;}

NamedObjMgr omg;
AnyDataObj* mobjv1 = omg.GetObj(tokens[0]);
if( mobjv1==NULL)
  {cout<<"PAWExecutor::h_copy Error: unknow object "<<tokens[0]<<endl;
  return;}
AnyDataObj* mobjv2 = omg.GetObj(tokens[1]);

int_4 i1=0,i2=0, j1=0,j2=0, ic1=0,jc1=0, i,ii, j,jj, nx1,ny1, nx2,ny2;

// Cas d'un Vector
if(typeid(*mobjv1) == typeid(Vector)) {
  Vector* v1 = dynamic_cast<Vector*>(mobjv1); nx1 = (int_4)v1->NElts();
  i2=nx1-1;
  if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
  if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1;
  if(i2<i1)
    {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2<<"] for NElts="
         <<nx1<<endl; return;}
  Vector* v2 = NULL;
  if(mobjv2==NULL)
    {v2 = new Vector(i2-i1+1); omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
  if(typeid(*mobjv2) != typeid(Vector))
    {cout<<"PAWExecutor::h_copy Error: type mismatch for Vector "
         <<typeid(*mobjv2).name()<<endl; return;}
  v2 = dynamic_cast<Vector*>(mobjv2); nx2 = (int_4)v2->NElts();
  if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d",&ic1);
  if(ic1<0) ic1=0;
  if(ic1>=nx2)
    {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<"] for NElts="
         <<nx2<<endl; return;}
  cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<") to "<<tokens[1]<<"("<<ic1<<"->...)"<<endl;
  for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++) (*v2)(ii) = (*v1)(i);
  return;
}

// Cas d'un Histo
if(typeid(*mobjv1) == typeid(Histo)) {
  Histo* v1 = dynamic_cast<Histo*>(mobjv1); nx1 = (int_4)v1->NBins();
  i2=nx1-1;
  if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
  if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1;
  if(i2<i1)
    {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2<<"] for NBins="
         <<nx1<<endl; return;}
  Histo* v2 = NULL;
  if(mobjv2==NULL)
    {v2 = new Histo(0.,(r_8)(i2-i1+1),i2-i1+1);
     omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
  if(typeid(*mobjv2) != typeid(Histo))
    {cout<<"PAWExecutor::h_copy Error: type mismatch for Histo "
         <<typeid(*mobjv2).name()<<endl; return;}
  v2 = dynamic_cast<Histo*>(mobjv2); nx2 = (int_4)v2->NBins();
  if(v1->HasErrors() && !(v2->HasErrors())) v2->Errors();
  if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d",&ic1);
  if(ic1<0) ic1=0;
  if(ic1>=nx2)
    {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<"] for NBins="
         <<nx2<<endl; return;}
  cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<") to "<<tokens[1]<<"("<<ic1<<"->...)"<<endl;
  for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
    {(*v2)(ii) = (*v1)(i); if(v1->HasErrors()) v2->Error2(ii) = v1->Error2(i);}
  return;
}

// Cas d'une Matrix
if(typeid(*mobjv1) == typeid(Matrix)) {
  Matrix* v1 = dynamic_cast<Matrix*>(mobjv1); nx1=v1->NRows(); ny1=v1->NCol();
  i2=nx1-1; j2=ny1-1;
  if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
  if(tks>3) if(tokens[3]!="!") sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
  if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1; if(j1<0) j1=0; if(j2>=ny1) j2=ny1-1;
  if(i2<i1 || j2<j1)
    {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2
         <<"] , ["<<j1<<":"<<j2<<"] for NRows,NCol="
         <<nx1<<" , "<<ny1<<endl; return;}
  Matrix* v2 = NULL;
  if(mobjv2==NULL)
    {v2 = new Matrix(i2-i1+1,j2-j1+1);
     omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
  if(typeid(*mobjv2) != typeid(Matrix))
    {cout<<"PAWExecutor::h_copy Error: type mismatch for Matrix "
         <<typeid(*mobjv2).name()<<endl; return;}
  v2 = dynamic_cast<Matrix*>(mobjv2); nx2=v2->NRows(); ny2=v2->NCol();
  if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d:%d",&ic1,&jc1);
  if(ic1<0) ic1=0; if(jc1<0) jc1=0;
  if(ic1>=nx2 || jc1>=ny2)
    {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<","<<jc1
         <<"] for NRows,NCol="<<nx2<<","<<ny2<<endl; return;}
  cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<","<<j1<<":"<<j2
      <<") to "<<tokens[1]<<"("<<ic1<<","<<jc1<<"->...)"<<endl;
  for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
    for(j=j1,jj=jc1; j<=j2 && j<ny1 && jj<ny2; j++,jj++) (*v2)(ii,jj) = (*v1)(i,j);
  return;
}

// Cas d'un Histo2D
if(typeid(*mobjv1) == typeid(Histo2D)) {
  Histo2D* v1 = dynamic_cast<Histo2D*>(mobjv1); nx1=v1->NBinX(); ny1=v1->NBinY();
  i2=nx1-1; j2=ny1-1;
  if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
  if(tks>3) if(tokens[3]!="!") sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
  if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1; if(j1<0) j1=0; if(j2>=ny1) j2=ny1-1;
  if(i2<i1 || j2<j1)
    {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2
         <<"] , ["<<j1<<":"<<j2<<"] for NBinX,NBinY="
         <<nx1<<" , "<<ny1<<endl; return;}
  Histo2D* v2 = NULL;
  if(mobjv2==NULL)
    {v2 = new Histo2D(0.,(r_8)(i2-i1+1),i2-i1+1,0.,(r_8)(j2-j1+1),j2-j1+1);
     omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
  if(typeid(*mobjv2) != typeid(Histo2D))
    {cout<<"PAWExecutor::h_copy Error: type mismatch for Histo2D"
         <<typeid(*mobjv2).name()<<endl; return;}
  v2 = dynamic_cast<Histo2D*>(mobjv2); nx2=v2->NBinX(); ny2=v2->NBinY();
  if(v1->HasErrors() && !(v2->HasErrors())) v2->Errors();
  if(tks>4) sscanf(tokens[4].c_str(),"%d:%d",&ic1,&jc1);
  if(ic1<0) ic1=0; if(jc1<0) jc1=0;
  if(ic1>=nx2 || jc1>=ny2)
    {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<","<<jc1
         <<"] for NBinX,NBinY="<<nx2<<","<<ny2<<endl; return;}
  cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<","<<j1<<":"<<j2
      <<") to "<<tokens[1]<<"("<<ic1<<","<<jc1<<"->...)"<<endl;
  for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
    for(j=j1,jj=jc1; j<=j2 && j<ny1 && jj<ny2; j++,jj++)
      {(*v2)(ii,jj) = (*v1)(i,j); if(v1->HasErrors()) v2->Error2(ii,jj) = v1->Error2(i,j);}
  return;
}

// Cas non prevu
cout<<"PAWExecutor::h_copy Error: type mismatch for Vector/Matrix/Histo/Histo2D\n"
    <<typeid(*mobjv1).name()<<endl;
return;
}

/* methode */
void PAWExecutor::h_set(string dum,vector<string>& tokens)
// Mise de valeurs/erreurs d'un histo 1D ou 2D a une valeur donnee
// dum = err  : on change les valeurs des erreurs
//       cont : on change les valeurs du contenu des bins
// tokens[0] : nom de l'histogramme
// tokens[1] : valeur de remplacement
// tokens[2-3] : i1:i2 j1:j2  intervalle des bins qui doivent etre remplace
//             : v v1:v2  remplacement des bins ayant une valeur dans cet intervalle
//             : e e1:e2  remplacement des bins ayant une erreur dans cet intervalle
{
int_4 tks = tokens.size();
if(tks<3)
  {cout<<"Usage: h/set/[err,cont] namehisto setvalue i1[:i2]  [j1[:j2]]\n"
       <<"                                           v v1:v2\n"
       <<"                                           e e1:e2"
       <<endl; return;}

// Decodage arguments
bool replerr = false; if(dum=="err") replerr = true;
r_8 setval = atof(tokens[1].c_str());
int_4 testcont=0; if(tokens[2]=="v") testcont=1; if(tokens[2]=="e") testcont=2;
int_4 i1=-1, i2=-1, j1=-1, j2=-1;
r_8 v1=0., v2=0.;
if(testcont==0) {
  sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
  if(tks>3) sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
} else {
  if(tks<=3)
    {cout<<"PAWExecutor::h_set Error: h/set/... v v1:v2, please give v1:v2"<<endl;
    return;}
  sscanf(tokens[3].c_str(),"%lf:%lf",&v1,&v2);
}

if(replerr) {if(setval<0.) setval = 0.; else setval *= setval;}
if(testcont!=0 && v2<v1)
  {cout<<"PAWExecutor::h_set Error: bad value range="<<v1<<","<<v2<<endl; return;}

// identification objet
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if( mobj==NULL)
  {cout<<"PAWExecutor::h_set Error: unknow object "<<tokens[0]<<endl;
  return;}

// Traitement
if(typeid(*mobj) == typeid(Histo)) {  // Histo 1D
  Histo* h  = dynamic_cast<Histo*>(mobj);
  if( (replerr || testcont==2) && !(h->HasErrors()) )
    {cout<<"PAWExecutor::h_set Error: histogram has no errors"<<endl; return;}
  if(testcont!=0) {i1=0; i2=h->NBins()-1;}
  if(i1<0 || i1>=h->NBins())
    {cout<<"PAWExecutor::h_set Error: bad bin range i1="<<i1<<endl; return;}
  if(i2<i1) i2=i1; if(i2>=h->NBins()) i2=h->NBins()-1;
  for(int_4 i=i1;i<=i2;i++) {
    bool change = true;
    if(testcont==1) {if((*h)(i)<v1 || (*h)(i)>v2) change = false;}
    else if(testcont==2) {if(h->Error(i)<v1 || h->Error(i)>v2) change = false;}
    if(!change) continue;
    if(replerr) h->Error2(i) = setval; else (*h)(i) = setval;
  }

} else if(typeid(*mobj) == typeid(Histo2D)) {  // Histo 2D
  Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
  if( (replerr || testcont==2) && !(h2->HasErrors()) )
    {cout<<"PAWExecutor::h_set Error: histogram has no errors"<<endl; return;}
  if(testcont!=0) {i1=0; i2=h2->NBinX()-1;j1=0; j2=h2->NBinY()-1;}
  if(i1<0 || i1>=h2->NBinX() || j1<0 || j1>=h2->NBinY())
    {cout<<"PAWExecutor::h_set Error: bad bin range i1,j1="<<i1<<","<<j1<<endl; return;}
  if(i2<i1) i2=i1; if(i2>=h2->NBinX()) i2=h2->NBinX()-1;
  if(j2<j1) j2=j1; if(j2>=h2->NBinY()) j2=h2->NBinY()-1;
  for(int_4 i=i1;i<=i2;i++) for(int_4 j=j1;j<=j2;j++) {
    bool change = true;
    if(testcont==1) {if((*h2)(i,j)<v1 || (*h2)(i,j)>v2) change = false;}
    else if(testcont==2) {if(h2->Error(i,j)<v1 || h2->Error(i,j)>v2) change = false;}
    if(!change) continue;
    if(replerr) h2->Error2(i,j) = setval; else (*h2)(i,j) = setval;
  }

} else {
  cout<<"PAWExecutor::h_set Error: type mismatch for Histo/Histo2D\n"
    <<typeid(*mobj).name()<<endl;
  return;
}

}

/* methode */
void PAWExecutor::h_err(vector<string>& tokens)
// Mise des erreurs d'un histo 1D ou 2D a une valeur
// donnee par une fonction de la valeur du bin
// tokens[0] : nom de l'histogramme
// tokens[1] : expression de la fonction
{
if(tokens.size()<2)
  {cout<<"Usage: h/err namehisto expr_func"<<endl; return;}

// identification objet
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(tokens[0]);
if( mobj==NULL)
  {cout<<"PAWExecutor::h_err Error: unknow object "<<tokens[0]<<endl;
  return;}

// Fonction
FILE *fip = NULL;
string expfunc = ""; {for(int_4 i=1;i<(int)tokens.size();i++) expfunc += tokens[i];}
string fname = "func1or2_pia_dl.c";
string func = "func1or2_pia_dl_func";
if((fip = fopen(fname.c_str(), "w")) == NULL)
  {cout<<"PAWExecutor::h_err Error: open function file "<<fname<<endl;
  return;}
fprintf(fip,"#include <math.h>\n");
fprintf(fip,"double %s(double x) \n{\n",func.c_str());
fprintf(fip,"return(%s); \n}\n", expfunc.c_str());
fclose(fip);
DlFunctionOfX f = (DlFunctionOfX) omg.GetServiceObj()->LinkFunctionFromFile(fname,func);
if(!f) 
  {cout<<"PAWExecutor::h_err Error: bad function "<<func<<","<<fname<<endl;
  return;}

// Traitement
if(typeid(*mobj) == typeid(Histo)) {  // Histo 1D
  Histo* h  = dynamic_cast<Histo*>(mobj);
  if(!(h->HasErrors())) h->Errors();
  for(int_4 i=0;i<h->NBins();i++) {
    r_8 x = (*h)(i);
    r_8 e = f(x);
    h->SetErr2(i,e*e); 
  }

} else if(typeid(*mobj) == typeid(Histo2D)) {  // Histo 2D
  Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
  if(!(h2->HasErrors())) h2->Errors();
  for(int_4 i=0;i<h2->NBinX();i++) for(int_4 j=0;j<h2->NBinY();j++) {
    r_8 x = (*h2)(i,j);
    r_8 e = f(x);
    h2->Error2(i,j) = e*e; 
  }

} else {
  cout<<"PAWExecutor::h_err Error: type mismatch for Histo/Histo2D\n"
    <<typeid(*mobj).name()<<endl;
  return;
}

}