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

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

#include "strutil.h"
#include "histos.h"
#include "histos.h"
#include "histos2.h"
#include "hisprof.h"
#include "ntuple.h"

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

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

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

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";
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  Related commands: plot2dw plot3d";  
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  Related commands: projh1d projh2d projprof";  
piac->RegisterCommand(kw,usage,this,hgrp);

kw = "h/integ";
usage = "Integrate a 1D histogram or profile";
usage += "\n h/integ nameh1d [norm]";
usage += "\n  Related commands: h/deriv";  
piac->RegisterCommand(kw,usage,this,hgrp);

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

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

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

kw = "h/cmult";
usage = "Multiply an histogramme by a constant";
usage += "\n h/cmult val";
usage += "\n  Related commands: h/cadd";  
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);
}

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

/* methode */
int PAWExecutor::Execute(string& kw, vector<string>& tokens)
{
if(kw == "reset") {
  reset(tokens); return(0);
} else if(kw == "n/plot") {
  n_plot(tokens); return(0);
} else if(kw == "n/proj") {
  n_proj(tokens); return(0);
} else if(kw == "h/integ") {
  h_integ(tokens); return(0);
} else if(kw == "h/deriv") {
  h_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/plot/2d") {
  h_plot_2d(tokens); return(0);
} else return(1); 
}

/* 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(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 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;
    }
  }
  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 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<=0) {
  cout<<"PAWExecutor::n_proj Error: bad coding "<<tokens[1]<<endl;
} else if(nvar==1) {
  // c'est une projection dans un histo 1D
  srvo->ProjectH1(nameobj,expx,expwt,expcut,nameproj,dopt,loop);
} else {
  // c'est une projection dans un histo2D
  // OU un HProf si nameproj est 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);
}

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;}
double norm = 1.;
if(tokens.size()>=2) norm = atof(tokens[1].c_str());
Histo* h1 = dynamic_cast<Histo*>(mobj);
// Le HProf se fait tout seul avec Histo: dynamic_cast<Histo*>(HProf)=Vrai!
//HProf* hp = dynamic_cast<HProf*>(mobj);
if(h1)      h1->HInteg(norm);
//else if(hp) hp->HInteg(norm);
else cout<<"PAWExecutor::h_integ Error: "<<tokens[0]<<" not an Histo/HProf"<<endl;
}

/* 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;}
Histo* h1 = dynamic_cast<Histo*>(mobj);
//HProf* hp = dynamic_cast<HProf*>(mobj);
if(h1)      h1->HDeriv();
//else if(hp) hp->HDeriv();
else cout<<"PAWExecutor::h_deriv Error: "<<tokens[0]<<" not an Histo/HProf"<<endl;
}

/* 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 nbin = atoi(tokens[1].c_str());
Histo* h1 = dynamic_cast<Histo*>(mobj);
//HProf* hp = dynamic_cast<HProf*>(mobj);
if(h1)      h1->HRebin(nbin);
//else if(hp) hp->HRebin(nbin);
else cout<<"PAWExecutor::h_rebin Error: "<<tokens[0]<<" not an Histo/HProf"<<endl;
}

/* 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;}
double val = atof(tokens[1].c_str());
Histo*   h1 = dynamic_cast<Histo*>(mobj);
Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
if(h1)      *h1 += val;
else if(h2) *h2 += val;
else cout<<"PAWExecutor::h_cadd Error: "<<tokens[0]<<" not an Histo/HProf/Histo2D"<<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;}
double val = atof(tokens[1].c_str());
Histo*   h1 = dynamic_cast<Histo*>(mobj);
Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
if(h1)      *h1 *= val;
else if(h2) *h2 *= val;
else cout<<"PAWExecutor::h_cmult Error: "<<tokens[0]
         <<" not an Histo/HProf/Histo2D"<<endl;
}

/* 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()<2)
  {cout<<"Usage: h/plot/2d nameh2d to_plot [n/s] [dopt]"<<endl; return;}
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";
AnyDataObj* mobjh1 = omg.GetObj(nametoplot);
// if(mobjh1!=NULL) omg.DelObj(nametoplot);  $CHECK$ Reza 12/12/99 - Ne pas faire - autoclean

string dopt = ""; if(tokens.size()>=3) dopt = tokens[2];
if(tokens[1] == "show") {
  h2->ShowProj();
  h2->ShowBand(2);
  h2->ShowSli(2);
  return;
} else if(tokens[1] == "h") {
  nametoplot = tokens[0];
} else if(tokens[1] == "px") {  
  if((h1p=h2->HProjX())) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowProj(); return;}
} else if(tokens[1] == "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 n = atoi(tokens[2].c_str());
  dopt = ""; if(tokens.size()>=4) dopt = tokens[3];
  if(tokens[1] == "bx") {
    if((h1p=h2->HBandX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowBand(); return;}
  } else if(tokens[1] == "by") {
    if((h1p=h2->HBandY(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowBand(); return;}
  } else if(tokens[1] == "sx") {
    if((h1p=h2->HSliX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
    else {h2->ShowSli(); return;}
  } else if(tokens[1] == "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 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 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 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;
}