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

Last change on this file since 2816 was 2816, checked in by cmv, 20 years ago

n/copy copy de ntuple avec choix de variables cmv 12/07/05

File size: 73.3 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <ctype.h>
4#include <iostream>
5#include <typeinfo>
6
7#include "sopnamsp.h"
8#include "strutil.h"
9#include "strutilxx.h"
10#include "histos.h"
11#include "histos2.h"
12#include "hisprof.h"
13#include "histerr.h"
14#include "ntuple.h"
15#include "datatable.h"
16
17#include "pawexecut.h"
18#include "nobjmgr.h"
19#include "servnobjm.h"
20#include "nomgadapter.h"
21#include "pistdimgapp.h"
22#include "pihisto.h"
23
24#ifdef SANS_EVOLPLANCK
25#include "cvector.h"
26#include "matrix.h"
27#else
28#include "tmatrix.h"
29#include "tvector.h"
30#endif
31
32/* Reza + cmv 13/10/99 */
33
34inline bool __vector_sort_up_r8(r_8 x1,r_8 x2) {return x1<x2;}
35inline bool __vector_sort_down_r8(r_8 x1,r_8 x2) {return x1>x2;}
36inline bool __vector_sort_up_r4(r_4 x1,r_4 x2) {return x1<x2;}
37inline bool __vector_sort_down_r4(r_4 x1,r_4 x2) {return x1>x2;}
38inline bool __vector_sort_up_i4(int_4 x1,int_4 x2) {return x1<x2;}
39inline bool __vector_sort_down_i4(int_4 x1,int_4 x2) {return x1>x2;}
40inline bool __vector_sort_up_u8(uint_8 x1,uint_8 x2) {return x1<x2;}
41inline bool __vector_sort_down_u8(uint_8 x1,uint_8 x2) {return x1>x2;}
42
43uint_4 PAWExecutor::autoc_counter_ = 0;
44
45/* methode */
46PAWExecutor::PAWExecutor(PIACmd *piac, PIStdImgApp* app)
47: mApp(app)
48{
49string kw, usage;
50string hgrp = "pawCmd";
51
52kw = "reset";
53usage = "Reset histograms vectors or matrix";
54usage += "\n reset nameobj";
55piac->RegisterCommand(kw,usage,this,hgrp);
56
57kw = "n/plot";
58usage = "Plot NTuple variables \"a la paw\" (alias n/pl)";
59usage += "\n n/plot nameobj.x_exp [cut] [w_exp] [loop] [gratt]";
60usage += "\n n/plot nameobj.y_exp%x_exp [cut] [loop] [gratt]";
61usage += "\n n/plot nameobj.z_exp%y_exp%x_exp [cut] [loop] [gratt]";
62usage += "\n for default use ! , loop=i1[:i2[:di]]";
63usage += "\n for 1 dimensional (1D) projection:";
64usage += "\n use graphic option \"keepbin\" to keep previous";
65usage += "\n plot binning for 1D distribution";
66usage += "\n Related commands: plot2dw plot3d";
67piac->RegisterCommand(kw,usage,this,hgrp);
68
69kw = "n/pl";
70usage = "alias to n/plot (see n/plot)";
71piac->RegisterCommand(kw,usage,this,hgrp);
72
73kw = "n/proj";
74usage = "Project NTuple in histogram (1D or 2D) a la paw";
75usage += "\n n/proj nameproj nameobj.x_exp [cut] [w_exp] [loop] [gratt]";
76usage += "\n n/proj nameproj nameobj.y_exp%x_exp [cut] [w_exp] [loop] [gratt]";
77usage += "\n for default use ! , loop=i1[:i2[:di]]";
78usage += "\n for 1 dimensional (1D) projection:";
79usage += "\n no display is performed if \"nameproj\" is an existing histogram";
80usage += "\n unless a graphic option \"gratt\" is given";
81usage += "\n Related commands: projh1d projh2d projprof exptovec";
82piac->RegisterCommand(kw,usage,this,hgrp);
83
84kw = "n/scan";
85usage = "Scan NTuple a la paw";
86usage += "\n n/scan nameobj[.exp1%exp2%exp3] cut loop";
87usage += "\n [-f:filename] [list_of_variables]";
88usage += "\n loop : iev1[:iev2[:diev]] or !";
89usage += "\n cut : cut expression or 1. or !";
90usage += "\n list_of_variables : default is all variables";
91usage += "\n : var1 var2 var3 ... varn";
92usage += "\n : var1 : var2 (from var1 to var2)";
93usage += "\n : : var2 (from first variable to var2)";
94usage += "\n : var1 : (from var1 to last variable)";
95usage += "\n ex: \"v1 : v3 v7 v4 : v6 v2 v9 :\"";
96usage += "\n exp1%exp2%exp3 :";
97usage += "\n if given add exp1,exp2,exp3 to the variable list";
98usage += "\n -f:filename : write into \"filename\", Default is to stdout";
99piac->RegisterCommand(kw,usage,this,hgrp);
100
101kw = "n/read";
102usage = "Read columns in an ASCII file and fill a NTuple or a DataTable";
103usage += "\n n/read nt fascii [options] var_1,c_1 var_2,c_2 ... var_n,c_n ";
104usage += "\n nt : NTuple or DataTable (see options) name";
105usage += "\n var_i,c_i : variable name, associated column in ASCII file [0,n[";
106usage += "\n where [options] are:";
107usage += "\n \"-dt\": create a DataTable of double instead of a NTuple";
108usage += "\n \"=s\": separator character is \'s\' (could be \"\t\")";
109usage += "\n \"-^abcd\": do not read lines beginning with string \"abcd\" ";
110usage += "\n \"+^abcd\": read only lines beginning with string \"abcd\" ";
111usage += "\n \"-abcd\": do not read lines which contain string \"abcd\" ";
112usage += "\n \"+abcd\": read only lines which contain string \"abcd\" ";
113usage += "\n these options may be repeated (ex: \"-^abcd\" \"-^xyz\") ";
114usage += "\n - in case of \"do not read\" options are added with logical AND ";
115usage += "\n - in case of \"read only\" options are added with logical OR ";
116piac->RegisterCommand(kw,usage,this,hgrp);
117
118kw = "n/merge";
119usage = "Merge ntuples (by adding lines)";
120usage += "\n n/merge nt nt_1 nt_2 ... nt_n";
121usage += "\n Merge ntuples nt_i lines into ntuple nt";
122piac->RegisterCommand(kw,usage,this,hgrp);
123
124kw = "n/merge/col";
125usage = "Merge ntuples (by adding columns)";
126usage += "\n n/merge/col nt nt_1[,ext] nt_2[,ext] ... nt_n[,ext]";
127usage += "\n Merge ntuples nt_i columns into ntuple nt";
128usage += "\n ext : character string to add at the end of the ntuple variable names";
129usage += "\n if \"ext\"=\"!\" no extension is added";
130usage += "\n if no \"ext\" is given, automatic extension \"_0\" \"_1\" \"_2\" ...";
131piac->RegisterCommand(kw,usage,this,hgrp);
132
133kw = "n/copy";
134usage = "Copy all or some variables of a ntuple into another new ntuple";
135usage += "\n n/merge ntnew nt [vname1 vname2 vname3 ...]";
136piac->RegisterCommand(kw,usage,this,hgrp);
137
138kw = "h/integ";
139usage = "Integrate a 1D histogram";
140usage += "\n h/integ nameh1d [norm]";
141usage += "\n norm<=0 means no normalisation (def norm=1)";
142usage += "\n Related commands: h/deriv v/integ v/deriv";
143piac->RegisterCommand(kw,usage,this,hgrp);
144
145kw = "v/integ";
146usage = "Integrate a TVector / vector";
147usage += "\n v/integ namevec [norm]";
148usage += "\n norm<=0 means no normalisation (def norm=0)";
149usage += "\n Related commands: h/integ h/deriv v/deriv";
150piac->RegisterCommand(kw,usage,this,hgrp);
151
152kw = "v/sort";
153usage = "Sort a vector into itself";
154usage += "\n v/sort namevec [+1/-1]";
155usage += "\n \"+1\" means increasing order, \"-1\" decreasing order";
156piac->RegisterCommand(kw,usage,this,hgrp);
157
158kw = "h/deriv";
159usage = "Derivate a 1D histogram";
160usage += "\n h/deriv nameh1d";
161usage += "\n Related commands: h/integ v/integ v/deriv";
162piac->RegisterCommand(kw,usage,this,hgrp);
163
164kw = "v/deriv";
165usage = "Derivate a TVector / vector";
166usage += "\n v/deriv namevec [deriv_option]";
167usage += "\n deriv_option -1 replace v[i] with v[i]-v[i-1]";
168usage += "\n 0 replace v[i] with (v[i+1]-v[i-1])/2 (default)";
169usage += "\n +1 replace v[i] with v[i+1]-v[i]";
170usage += "\n Related commands: h/integ h/deriv v/integ";
171piac->RegisterCommand(kw,usage,this,hgrp);
172
173kw = "h/rebin";
174usage = "Rebin a 1D histogram or profile";
175usage += "\n h/rebin nameh1d nbin";
176piac->RegisterCommand(kw,usage,this,hgrp);
177
178kw = "h/cadd";
179usage = "Add a constant to an histogram, a vector or a matrix";
180usage += "\n h/cadd namehisto val";
181usage += "\n Related commands: h/cmult h/oper";
182piac->RegisterCommand(kw,usage,this,hgrp);
183
184kw = "h/cmult";
185usage = "Multiply an histogram, a vector or a matrix by a constant";
186usage += "\n h/cmult namehisto val";
187usage += "\n Related commands: h/cadd h/oper";
188piac->RegisterCommand(kw,usage,this,hgrp);
189
190kw = "h/oper";
191usage = "Operation on histograms vectors or matrices";
192usage += "\n h/oper @ h1 h2 hres";
193usage += "\n hres = h1 @ h2 with @ = (+,-,*,/)";
194usage += "\n For vectors and matrices, operations are elements by elements";
195usage += "\n Related commands: h/cadd h/cmult";
196piac->RegisterCommand(kw,usage,this,hgrp);
197
198kw = "h/plot/2d";
199usage = "Specific plot for 2D histogrammes";
200usage += "\n h/plot/2d nameh2d show : infos on 2D histogramme";
201usage += "\n h/plot/2d nameh2d h [dopt] : plot 2D histogramme";
202usage += "\n h/plot/2d nameh2d px [dopt] : plot X projection";
203usage += "\n h/plot/2d nameh2d py [dopt] : plot Y projection";
204usage += "\n h/plot/2d nameh2d bx n [dopt] : plot X band number n";
205usage += "\n h/plot/2d nameh2d by n [dopt] : plot Y band number n";
206usage += "\n h/plot/2d nameh2d sx n [dopt] : plot X slice number n";
207usage += "\n h/plot/2d nameh2d sy n [dopt] : plot Y slice number n";
208usage += "\n n < 0 means Show Info";
209piac->RegisterCommand(kw,usage,this,hgrp);
210
211kw = "h/put_vec";
212usage = "Put content of vector (matrix) into content of histogram 1D or 2D";
213usage += "\n h/put_vec nameh1d namevector [cont,err2]";
214usage += "\n h/put_vec nameh2d namematrix [cont,err2]";
215usage += "\n cont : put into histogramme content";
216usage += "\n err2 : put into histogramme error^2";
217usage += "\n Related commands: h/get_vec";
218piac->RegisterCommand(kw,usage,this,hgrp);
219
220kw = "h/get_vec";
221usage = "Get content of histogram 1D profile or 2D into vector (matrix)";
222usage += "\n h/get_vec nameh1d namevector [cont,err2,absc] [proj]";
223usage += "\n h/get_vec nameh2d namematrix [cont,err2,absc]";
224usage += "\n cont : get histogramme content";
225usage += "\n err2 : get histogramme error^2";
226usage += "\n absc : get histogramme low bin abscissa (1D only)";
227usage += "\n proj :";
228usage += "\n show : show available projections for Histo2D";
229usage += "\n px : get X projection";
230usage += "\n py : get Y projection";
231usage += "\n bx n : get X band number n";
232usage += "\n by n : get Y band number n";
233usage += "\n sx n : get X slice number n";
234usage += "\n sy n : get Y slice number n";
235usage += "\n Related commands: h/put_vec";
236piac->RegisterCommand(kw,usage,this,hgrp);
237
238kw = "h/copy";
239usage = "Copy content of object1 into object2";
240usage += "\n objects are Vector,Matrix,Histo,Histo2D";
241usage += "\n h/copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]";
242usage += "\n copy obj1Dfrom(i1->i2) into obj1Dto(ic1->)";
243usage += "\n copy obj2Dfrom(i1,j1->i2,j2) into obj2Dto(ic1,jc1->)";
244usage += "\n Warning: elements from i1 to i2 included are copied";
245usage += "\n If obj1Dto does not exist, is is created with size i2-i1+1";
246usage += "\n or obj2Dto with size i2-i1+1,j2-j1+1";
247usage += "\n The adressed content of obj?Dfrom is overwritten";
248usage += "\n The non-adressed content of obj?Dfrom is left unchanged";
249usage += "\n Related commands: copy";
250piac->RegisterCommand(kw,usage,this,hgrp);
251
252kw = "h/set/err";
253usage = "Set Histo,Histo2D errors for range of bins or range of values";
254usage += "\n h/set/err namehisto setvalue i1[:i2] [j1[:j2]]";
255usage += "\n set error to setvalue for bin range i1:i2 j1:j2";
256usage += "\n h/set/err namehisto setvalue v v1:v2";
257usage += "\n set error to setvalue for content values range v1:v2";
258usage += "\n h/set/err namehisto setvalue e e1:e2";
259usage += "\n set error to setvalue for error values range v1:v2";
260usage += "\n Related commands: h/set/cont";
261piac->RegisterCommand(kw,usage,this,hgrp);
262
263kw = "h/set/cont";
264usage = "Set Histo,Histo2D content for range of bins or range of values";
265usage += "\n h/set/cont namehisto setvalue i1[:i2] [j1[:j2]]";
266usage += "\n set content to setvalue for bin range i1:i2 j1:j2";
267usage += "\n h/set/cont namehisto setvalue v v1:v2";
268usage += "\n set content to setvalue for content values range v1:v2";
269usage += "\n h/set/cont namehisto setvalue e e1:e2";
270usage += "\n set content to setvalue for error values range v1:v2";
271usage += "\n Related commands: h/set/err";
272piac->RegisterCommand(kw,usage,this,hgrp);
273
274kw = "h/err";
275usage = "Set Histo,Histo2D error to function of bin content value";
276usage += "\n h/err namehisto expr_func";
277usage += "\n Related commands: h/set/err";
278piac->RegisterCommand(kw,usage,this,hgrp);
279
280}
281
282/* methode */
283PAWExecutor::~PAWExecutor()
284{
285}
286
287/* methode */
288int PAWExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
289{
290if(kw == "reset") {
291 reset(tokens); return(0);
292} else if(kw == "n/plot" || kw == "n/pl") {
293 n_plot(tokens); return(0);
294} else if(kw == "n/proj") {
295 n_proj(tokens); return(0);
296} else if(kw == "n/scan") {
297 n_scan(tokens); return(0);
298} else if(kw == "n/read") {
299 n_read(tokens); return(0);
300} else if(kw == "n/merge") {
301 n_merge(tokens); return(0);
302} else if(kw == "n/merge/col") {
303 n_merge_col(tokens); return(0);
304} else if(kw == "n/copy") {
305 n_copy(tokens); return(0);
306} else if(kw == "h/integ") {
307 h_integ(tokens); return(0);
308} else if(kw == "v/integ") {
309 v_integ(tokens); return(0);
310} else if(kw == "v/sort") {
311 v_sort(tokens); return(0);
312} else if(kw == "h/deriv") {
313 h_deriv(tokens); return(0);
314} else if(kw == "v/deriv") {
315 v_deriv(tokens); return(0);
316} else if(kw == "h/rebin") {
317 h_rebin(tokens); return(0);
318} else if(kw == "h/cadd") {
319 h_cadd(tokens); return(0);
320} else if(kw == "h/cmult") {
321 h_cmult(tokens); return(0);
322} else if(kw == "h/oper") {
323 h_oper(tokens); return(0);
324} else if(kw == "h/plot/2d") {
325 h_plot_2d(tokens); return(0);
326} else if(kw == "h/put_vec") {
327 h_put_vec(tokens); return(0);
328} else if(kw == "h/get_vec") {
329 h_get_vec(tokens); return(0);
330} else if(kw == "h/copy") {
331 h_copy(tokens); return(0);
332} else if(kw == "h/set/err") {
333 string dum = "err";
334 h_set(dum,tokens); return(0);
335} else if(kw == "h/set/cont") {
336 string dum = "cont";
337 h_set(dum,tokens); return(0);
338} else if(kw == "h/err") {
339 h_err(tokens); return(0);
340} else return(1);
341}
342
343/* methode */
344bool PAWExecutor::IsThreadable(string const & keyw)
345{
346 if ( (keyw == "n/plot") || (keyw == "n/pl") ||
347 (keyw == "n/proj") ) return true;
348 else return false;
349}
350
351
352/* methode */
353void PAWExecutor::reset(vector<string>& tokens)
354// Reset d'histogrammes, vecteurs et matrices
355{
356if(tokens.size() < 1)
357 {cout<<"Usage: reset nameobj"<<endl; return;}
358NamedObjMgr omg;
359AnyDataObj* mobj = omg.GetObj(tokens[0]);
360if(mobj == NULL)
361 {cout<<"PAWExecutor::reset Error , Pas d'objet de nom "<<tokens[0]<<endl;
362 return;}
363string ctyp = typeid(*mobj).name();
364
365#ifdef SANS_EVOLPLANCK
366if(typeid(*mobj)==typeid(Vector)) {Vector* ob=(Vector*) mobj; ob->Zero();}
367else if(typeid(*mobj)==typeid(Matrix)) {Matrix* ob=(Matrix*) mobj; ob->Zero();}
368#else
369 if(typeid(*mobj)==typeid(Vector)) {Vector* ob=(Vector*) mobj; (*ob) = 0.; }
370// ob->DataBlock().Reset(0.);}
371 else if(typeid(*mobj)==typeid(Matrix)) {Matrix* ob=(Matrix*) mobj; (*ob) = 0.; }
372//ob->DataBlock().Reset(0.);}
373#endif
374else if(typeid(*mobj)==typeid(Histo)) {Histo* ob=(Histo*) mobj; ob->Zero();}
375else if(typeid(*mobj)==typeid(HProf)) {HProf* ob=(HProf*) mobj; ob->Zero();}
376else if(typeid(*mobj)==typeid(HistoErr)) {HistoErr* ob=(HistoErr*) mobj; ob->Zero();}
377else if(typeid(*mobj)==typeid(Histo2D)) {Histo2D* ob=(Histo2D*)mobj; ob->Zero();}
378else {
379 cout<<"PAWExecutor::reset Error , No reset possible on "<<ctyp<<endl;
380 return;
381}
382
383return;
384}
385
386/* methode */
387void PAWExecutor::n_plot(vector<string>& tokens)
388// Equivalent n/plot de paw
389// Plot 1D
390// n/plot nameobj.x_exp [cut] [w_exp] [gratt]
391// Plot 2D (plot2dw)
392// n/plot nameobj.y_exp%x_exp [cut] [w_exp] [gratt]
393// Plot 3D (plot3d)
394// n/plot nameobj.z_exp%y_exp%x_exp [cut] [gratt]
395{
396if(tokens.size() < 1) {
397 cout
398 <<"Usage: n/plot nameobj.x_exp [cut] [w_exp] [loop] [gratt] [nomh1]"<<endl
399 <<" n/plot nameobj.y_exp%x_exp [cut] [loop] [gratt]"<<endl
400 <<" n/plot nameobj.z_exp%y_exp%x_exp [cut] [loop] [gratt]"<<endl
401 <<" for default use ! , loop=i1[:i2[:di]]"<<endl;
402 return;
403}
404string nameobj,expx,expy,expz;
405int_4 nvar = decodepawstring(tokens[0],nameobj,expx,expy,expz);
406string expcut = "1"; string expwt = "1."; string loop = "";
407string dopt = ""; string nameproj="";
408if(tokens.size()>=2) expcut = tokens[1]; if(expcut=="!") expcut="1";
409
410NamedObjMgr omg;
411Services2NObjMgr* srvo = omg.GetServiceObj();
412
413if(nvar<=0) {
414 cout<<"PAWExecutor::n_plot Error: bad coding "<<tokens[0]<<endl;
415} else if(nvar==1) { // c'est un plot 1D
416 if(tokens.size()>=3) expwt = tokens[2]; if(expwt=="!") expwt="1.";
417 if(tokens.size()>=4) loop = tokens[3]; if(loop=="!") loop="";
418 if(tokens.size()>=5) dopt = tokens[4]; if(dopt=="!") dopt="";
419 if(tokens.size()>=6) nameproj = tokens[5];
420 if(nameproj.length()<=0 || nameproj=="!") {
421 nameproj = "/autoc/paw_n_plot1D_";
422 AnyDataObj* mobj = omg.GetObj(nameproj);
423 if(mobj!=NULL) {
424 char buff[16]; autoc_counter_++; sprintf(buff,"%d",autoc_counter_);
425 nameproj += buff;
426 }
427 // --- Si option "keepbin" on projete avec le meme binning que l'histo 1D precedent
428 if(dopt.find("keepbin")<dopt.size()) {
429 PIStdImgApp* piimapp = omg.GetImgApp();
430 if(piimapp) {
431 PIBaseWdg* pibwdg = piimapp->CurrentBaseWdg();
432 if(pibwdg) {
433 int nbdrw = pibwdg->NbDrawers();
434 if(nbdrw>0) {
435 for(int i=nbdrw-1;i>=0;i--) {
436 PIDrawer* pidwr = pibwdg->GetDrawer(i);
437 PIHisto* pih = NULL;
438 if( (pih = dynamic_cast<PIHisto *>(pidwr))==NULL ) continue;
439 Histo* h = pih->Histogram();
440 if(h==NULL) continue;
441 if(h->NBins()<1) continue;
442 Histo* hsame = new Histo(h->XMin(),h->XMax(),h->NBins());
443 omg.AddObj(hsame,nameproj);
444 break;
445 }
446 }
447 }
448 }
449 }
450 // ---
451 }
452 srvo->ProjectH1(nameobj,expx,expwt,expcut,nameproj,dopt,loop);
453} else if(nvar==2) { // c'est un plot 2D
454 if(tokens.size()>=3) loop = tokens[2]; if(loop=="!") loop="";
455 if(tokens.size()>=4) dopt = tokens[3];
456 string err = "";
457 srvo->DisplayPoints2D(nameobj,expx,expy,err,err,expcut,dopt,loop);
458} else { // c'est un plot 3D
459 if(tokens.size()>=3) loop = tokens[2]; if(loop=="!") loop="";
460 if(tokens.size()>=4) dopt = tokens[3];
461 srvo->DisplayPoints3D(nameobj,expx,expy,expz,expcut,dopt,loop);
462}
463
464return;
465}
466
467
468/* methode */
469void PAWExecutor::n_proj(vector<string>& tokens)
470// Equivalent n/proj de paw
471// Project NTuple in histogram a la paw
472// Dans un Histo 1D
473// n/proj nameproj nameobj.x_exp [cut] [w_exp] [gratt]
474// Dans un Histo 2D ou un HProf (dans ce cas nameproj doit etre cree).
475// n/proj nameproj nameobj.y_exp%x_exp [cut] [w_exp] [gratt]
476{
477if(tokens.size()<2)
478 {cout<<"Usage: n/proj nameproj nameobj.[y_exp%]x_exp [cut] [w_exp] [loop] [gratt]"<<endl
479 <<" for default use ! , loop=i1[:i2[:di]]"<<endl; return;}
480string nameproj = tokens[0];
481string nameobj,expx,expy,expz;
482int_4 nvar = decodepawstring(tokens[1],nameobj,expx,expy,expz);
483string expcut = "1"; string expwt = "1."; string loop = ""; string dopt = "";
484if(tokens.size()>=3) expcut = tokens[2]; if(expcut=="!") expcut="1";
485if(tokens.size()>=4) expwt = tokens[3]; if(expwt=="!") expwt="1.";
486if(tokens.size()>=5) loop = tokens[4]; if(loop=="!") loop="";
487if(tokens.size()>=6) dopt = tokens[5];
488
489NamedObjMgr omg;
490Services2NObjMgr* srvo = omg.GetServiceObj();
491
492if(nvar==1) {
493 // c'est une projection dans un histo 1D
494 srvo->ProjectH1(nameobj,expx,expwt,expcut,nameproj,dopt,loop);
495} else if(nvar==2) {
496 // c'est une projection dans un histo2D
497 // OU un HProf si nameproj est un HProf un deja defini
498 AnyDataObj* mobj = omg.GetObj(nameproj);
499 if(mobj==NULL)
500 srvo->ProjectH2(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
501 else if(dynamic_cast<HProf*>(mobj))
502 srvo->ProjectHProf(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
503 else
504 srvo->ProjectH2(nameobj,expx,expy,expwt,expcut,nameproj,dopt,loop);
505} else {
506 cout<<"PAWExecutor::n_proj Error: bad coding "<<tokens[1]<<" nvar="<<nvar<<endl;
507}
508
509return;
510}
511
512/* methode */
513void PAWExecutor::n_scan(vector<string>& tokens)
514{
515if(tokens.size()<3)
516 {cerr<<"Usage: n/scan nameobj[.exp1%exp2%exp3] cut loop\n"
517 <<" [-f:filename] [list_of_variables]"<<endl; return;}
518
519// decodage des premiers arguments
520string nameobj,expr[4];
521int_4 nexpr = decodepawstring(tokens[0],nameobj,expr[0],expr[1],expr[2]);
522if(nexpr<4) {for(int_4 i=nexpr;i<4;i++) expr[i]="1.";}
523string expcut = tokens[1]; if(expcut=="!") expcut="1";
524string loop = tokens[2]; if(loop=="!") loop="";
525FILE* fout = NULL;
526uint_4 ldebvar = 3;
527if(tokens.size()>3) {
528 if(tokens[3].find("-f:") == 0) {
529 string filename = tokens[3].substr(3);
530 if(filename.size()>0) {
531 fout = fopen(filename.c_str(),"w");
532 if(fout==NULL)
533 {cerr<<"PAWExecutor::n_scan Error, file "<<filename
534 <<" not opened"<<endl; return;}
535 }
536 ldebvar++;
537 }
538}
539
540// ntuple adaptateur
541NamedObjMgr omg;
542Services2NObjMgr& srvo = *omg.GetServiceObj();
543NObjMgrAdapter* obja = omg.GetObjAdapter(nameobj); // Ne pas deleter
544if(obja == NULL)
545 {cerr<<"PAWExecutor::n_scan Error, ObjAdapter==NULL for "
546 <<nameobj<<endl; return;}
547bool adel = true;
548NTupleInterface* objnt = obja->GetNTupleInterface(adel);
549if(objnt == NULL)
550 {cerr<<"PAWExecutor::n_scan Error, NTupleInterface==NULL for "
551 <<nameobj<<endl; return;}
552
553// function pour le choix
554string vardec = objnt->VarList_C("_zz6qi_");
555PlotExprFunc f = srvo.LinkExprFunc(vardec,expr[0],expr[1],expr[2],expr[3],expcut);
556if(!f)
557 {cerr<<"PAWExecutor::n_scan Error, Creation PlotExprFunc"<<endl;
558 if(adel) delete objnt; if(fout) fclose(fout); return;}
559
560// variables a imprimer
561// "rien" --> de 0 a ncol-1 (toutes les variables)
562// : v1 --> de 0 a v1
563// v1 : --> de v1 a ncol-1
564// v1 : v2 --> de v1 a v2 (si v2 apres v1)
565// v1 et v2 (si v2 avant v1)
566// v1 v2 v3 v4 --> v1 v2 v3 v4 (ordre indifferent)
567// et toute combinaison : "v1 : v3 v7 v4 : v6 v2 v9 :"
568// --> v1 v2 v3 v7 v4 v5 v6 v2 v9 v10...v(ncol-1)
569int_4 ncol = objnt->NbColumns();
570if(ncol<=0)
571 {cerr<<"PAWExecutor::n_scan Error, no columns for NTuple"<<endl;
572 return;}
573vector<int_4> varnum;
574if(ldebvar>=tokens.size()) { // Toutes les variables
575 for(int_4 i=0;i<ncol;i++) varnum.push_back(i);
576} else { // Choix de certaines variables
577 int_4 k,klast,kk; bool frlast=false;
578 if(tokens[ldebvar]==":") {varnum.push_back(0); frlast=true;}
579 else {k = objnt->ColumnIndex(tokens[ldebvar]); varnum.push_back(k);}
580 ldebvar++;
581 if(ldebvar<tokens.size()) for(uint_4 i=ldebvar;i<tokens.size();i++) {
582 if(tokens[i]!=":") { // pas un separateur
583 k = klast = objnt->ColumnIndex(tokens[i]);
584 if(frlast) klast = varnum[varnum.size()-1] + 1;
585 if(klast>k) klast=k;
586 for(kk=klast;kk<=k;kk++) varnum.push_back(kk);
587 frlast=false;
588 } else if(i==tokens.size()-1) { // separateur a la fin
589 k = ncol-1;
590 klast = varnum[varnum.size()-1] + 1;
591 if(klast>k) klast=k;
592 for(kk=klast;kk<=k;kk++) varnum.push_back(kk);
593 } else frlast=true; // separateur pas a la fin
594 }
595}
596
597vector<string> varname;
598if(varnum.size()>0) for(int_4 i=0;i<(int)varnum.size();i++) {
599 if(varnum[i]<0 || varnum[i]>=ncol)
600 {cerr<<"PAWExecutor::n_scan Error, bad variable name at pos "
601 <<i<<endl; if(adel) delete objnt; if(fout) fclose(fout); return;}
602 string dum = objnt->ColumnName(varnum[i]);
603 varname.push_back(dum);
604}
605
606// evenements a utiliser
607int_8 k1=0, k2=objnt->NbLines(), dk=1;
608srvo.DecodeLoopParameters(loop,k1,k2,dk);
609if (k1<0) k1=0;
610if (k2<0) k2=objnt->NbLines();
611if (k2>(int_8)objnt->NbLines()) k2=objnt->NbLines();
612if (dk<=0) dk=1;
613
614// boucle sur les evenements et print
615try {
616 int_4 i;
617 if(fout) fprintf(fout,"#ev "); else printf("#ev ");
618 for(i=0;i<(int)varname.size();i++)
619 if(fout) fprintf(fout,"%s ",varname[i].c_str());
620 else printf( "%s ",varname[i].c_str());
621 if(nexpr>0) for(i=0;i<nexpr;i++)
622 if(fout) fprintf(fout,"%s ",expr[i].c_str());
623 else printf( "%s ",expr[i].c_str());
624 if(fout) fprintf(fout,"\n"); else printf("\n");
625
626 double xnt[5]={0,0,0,0,0};
627 double* xn;
628 for(int_8 k=k1; k<k2; k += dk) {
629 xn = objnt->GetLineD(k);
630 if(f((int_8_exprf)k,xn,xnt,xnt+1,xnt+2,xnt+3) != 0) {
631 if(fout) fprintf(fout,"%ld ",(long)k); else printf("%ld ",(long)k);
632 for(i=0;i<(int)varnum.size();i++) {
633 if(fout) fprintf(fout,"%g ",*(xn+varnum[i]));
634 else printf( "%g ",*(xn+varnum[i]));
635 }
636 if(nexpr>0) for(i=0;i<nexpr;i++) {
637 if(fout) fprintf(fout,"%g ",*(xnt+i));
638 else printf( "%g ",*(xnt+i));
639 }
640 if(fout) fprintf(fout,"\n"); else printf("\n");
641 }
642 }
643} // fin du try
644#ifdef SANS_EVOLPLANCK
645CATCH(merr) {
646 fflush(stdout); cout<<endl; cerr<<endl;
647 string es = PeidaExc(merr);
648 cerr<<"Services2NObjMgr::ComputeExpressions() Exception :"<<merr<<es;
649} ENDTRY;
650#else
651catch ( PException exc ) {
652 fflush(stdout); cout<<endl; cerr<<endl;
653 cerr<<"Services2NObjMgr::ComputeExpressions() Exception :"<<exc.Msg()<<endl;
654}
655#endif
656
657if(adel) delete objnt;
658if(fout) fclose(fout);
659srvo.CloseDLL(); // Fermeture du fichier .so
660return;
661}
662
663
664/* methode */
665#define __LENLINE_N_READ__ 8192
666void PAWExecutor::n_read(vector<string>& tokens)
667{
668 int lp=1;
669
670 if(tokens.size()<3) {
671 cerr<<"Usage: n/read nt fascii [options] var_1,c_1 ... var_n,c_n"<<endl;
672 return;
673 }
674
675 // decodage des arguments
676 string nament = tokens[0];
677 string nameascii = tokens[1];
678 bool create_datatable = false;
679 vector<string> donotreadbeg;
680 vector<string> donotreadin;
681 vector<string> onlyreadbeg;
682 vector<string> onlyreadin;
683 vector<string> varname;
684 vector<int> colnum;
685 char separator = ' ';
686 int numcolmaxi=-1;
687
688 for(int i=2;i<(int_4)tokens.size();i++) {
689 int lc = tokens[i].size();
690 if(lc<2) continue;
691 if(tokens[i].find("-dt")==0) { // create DataTable instead of NTuple
692 create_datatable = true;
693 continue;
694 }
695 const char *c = tokens[i].c_str();
696 if(c[0]=='=') { // Separator
697 separator = c[1];
698 if(lc==3) if(c[1]=='\\' && c[2]=='t') separator = '\t';
699 continue;
700 }
701 if(c[0]=='+') { // Selection des lignes a lire
702 if(c[1]!='^') onlyreadin.push_back(&c[1]);
703 else if(lc>2) onlyreadbeg.push_back(&c[2]);
704 continue;
705 }
706 if(c[0]=='-') { // Selection des lignes commentaire
707 if(c[1]!='^') donotreadin.push_back(&c[1]);
708 else if(lc>2) donotreadbeg.push_back(&c[2]);
709 continue;
710 }
711 // decodage des noms de variables et des colonnes associees
712 int p = tokens[i].find(',');
713 if(p<1 || p>=lc-1) continue;
714 string vn = tokens[i].substr(0,p);
715 string cn = tokens[i].substr(p+1,lc-p-1);
716 int ic = atoi(cn.c_str());
717 if(ic<0) continue;
718 if( !isalpha(vn[0]) ) continue;
719 if(ic>numcolmaxi) numcolmaxi = ic;
720 varname.push_back(vn);
721 colnum.push_back(ic);
722 }
723
724 int nvar = varname.size();
725 if(nvar<=0) {
726 cerr<<"n_read: no variables to be read"<<endl;
727 return;
728 }
729
730 // Print what has to be done
731 if(lp) {
732 if(onlyreadin.size()>0) {
733 cout<<"n_read Only read line containing ["<<onlyreadin.size()<<"]:";
734 for(int i=0;i<(int_4)onlyreadin.size();i++) cout<<" \'"<<onlyreadin[i]<<"\'";
735 cout<<endl;
736 }
737 if(onlyreadbeg.size()>0) {
738 cout<<"n_read Only read line begining with ["<<onlyreadbeg.size()<<"]:";
739 for(int i=0;i<(int_4)onlyreadbeg.size();i++) cout<<" \'"<<onlyreadbeg[i]<<"\'";
740 cout<<endl;
741 }
742 if(donotreadin.size()>0) {
743 cout<<"n_read Do not read line containing ["<<donotreadin.size()<<"]:";
744 for(int i=0;i<(int_4)donotreadin.size();i++) cout<<" \'"<<donotreadin[i]<<"\'";
745 cout<<endl;
746 }
747 if(donotreadbeg.size()>0) {
748 cout<<"n_read Do not read line begining with ["<<donotreadbeg.size()<<"]:";
749 for(int i=0;i<(int_4)donotreadbeg.size();i++) cout<<" \'"<<donotreadbeg[i]<<"\'";
750 cout<<endl;
751 }
752 if(nvar>0) {
753 cout<<"n_read Number of variables to be read: "<<nvar<<endl;
754 for(int i=0;i<nvar;i++) cout<<" \'"<<varname[i].c_str()<<","<<colnum[i]<<"\'";
755 cout<<endl;
756 }
757 if(separator!=' ') cout<<"n_read Separator is: \'"<<separator<<"\'"<<endl;
758 }
759
760 // Open ASCII file
761 FILE * fascii = fopen(nameascii.c_str(),"r");
762 if(fascii==NULL) {
763 cerr<<"n_read: cannot open file "<<nameascii<<endl;
764 return;
765 }
766
767 // Creation du NTuple ou de la DataTable
768 NTuple *nt = NULL;
769 DataTable *dt = NULL;
770 if(!create_datatable) { //create NTuple
771 char** ntvn = new char*[nvar];
772 for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(varname[i].c_str());
773 nt = new NTuple(nvar,ntvn);
774 delete [] ntvn;
775 } else { //create DataTable
776 dt = new DataTable();
777 for(int i=0;i<nvar;i++) dt->AddDoubleColumn(varname[i]);
778 }
779 r_8 *xnt = new r_8[nvar];
780
781 // Read file
782 char *line = new char[__LENLINE_N_READ__];
783 int nline=0, nlinecom=0;
784 while(fgets(line,__LENLINE_N_READ__,fascii) != NULL ) {
785 nline++;
786 uint_4 lc = strlen(line); if(lc<1) continue;
787 // Pour enlever le \n final
788 if(line[lc-1]=='\n' || line[lc-1]=='\r')
789 {line[lc-1]='\0'; lc = strlen(line); if(lc<1) continue;}
790
791 string const sline(line);
792 //cout<<"\'"<<sline<<"\' lc="<<lc<<endl;
793
794 // Faut t'il lire cette ligne ?
795 bool read_line_1 = true;
796 if(onlyreadin.size()>0 || onlyreadbeg.size()>0) read_line_1 = false;
797 if(onlyreadin.size()>0) {
798 for(int i=0;i<(int_4)onlyreadin.size();i++) {
799 uint_4 p = sline.find(onlyreadin[i].c_str());
800 if(p>=lc) continue;
801 read_line_1 = true;
802 break;
803 }
804 }
805 if(onlyreadbeg.size()>0) {
806 for(int i=0;i<(int_4)onlyreadbeg.size();i++) {
807 uint_4 p = sline.find(onlyreadbeg[i].c_str());
808 if(p!=0) continue;
809 read_line_1 = true;
810 break;
811 }
812 }
813
814 // Faut t'il ne pas lire cette ligne ?
815 bool read_line_2 = true;
816 if(donotreadin.size()>0) {
817 for(int i=0;i<(int_4)donotreadin.size();i++) {
818 uint_4 p = sline.find(donotreadin[i].c_str());
819 if(p>=lc) continue;
820 read_line_2 = false;
821 break;
822 }
823 }
824 if(donotreadbeg.size()>0) {
825 for(int i=0;i<(int_4)donotreadbeg.size();i++) {
826 uint_4 p = sline.find(donotreadbeg[i].c_str());
827 if(p!=0) continue;
828 read_line_2 = false;
829 break;
830 }
831 }
832 if(!read_line_2) nlinecom++;
833
834 if(!read_line_1 || !read_line_2) continue;
835
836 // Decodage de la ligne
837 vector<string> vs;
838 //FillVStringFrString(sline,vs,separator);
839 SplitStringToVString(sline,vs,separator);
840 int lvs = vs.size();
841 //for(int i=0;i<lvs;i++) cout<<"|"<<vs[i]<<"| "; cout<<endl;
842 if(lvs<numcolmaxi) continue; // Pas assez de champs decodes, mauvaise ligne
843
844 // Remplissage du NTuple ou de la DataTable
845 for(int i=0;i<nvar;i++) {
846 xnt[i] = 0.;
847 int ic = colnum[i];
848 if(ic>=lvs) continue;
849 xnt[i] = atof(vs[ic].c_str());
850 }
851 if(nt) nt->Fill(xnt);
852 else if(dt) dt->Fill(xnt);
853 //cout<<"...xnt"; for(int i=0;i<nvar;i++) cout<<" "<<xnt[i]; cout<<endl;
854
855 }
856 int_4 nentries = (nt)? nt->NEntry(): dt->NEntry();
857 cout<<"n_read: "<<nline<<" lines in file, "
858 <<nlinecom<<" commentary, "
859 <<nentries<<" object entries"<<endl;
860
861 // On sauve le NTuple/DataTable si besoin, on ferme et detruit ce qu'il faut
862 NamedObjMgr omg;
863 if(nentries>0) {
864 if(nt) omg.AddObj(nt,nament);
865 if(dt) omg.AddObj(dt,nament);
866 } else {
867 if(nt) delete nt;
868 if(dt) delete dt;
869 }
870 delete [] xnt;
871 delete [] line;
872 fclose(fascii);
873
874 return;
875}
876#undef __LENLINE_N_READ__
877
878/* methode */
879void PAWExecutor::n_merge(vector<string>& tokens)
880{
881 if(tokens.size()<2) {
882 cerr<<"Usage: n/merge nt nt_1 nt_2 ... nt_n"<<endl;
883 return;
884 }
885
886 NamedObjMgr omg;
887
888 // decodage des arguments
889 string nament = tokens[0];
890
891 // boucle sur les ntuples
892 NTuple * nt = NULL;
893 int nvar=0, nfill=0;
894 r_8 *xnt=NULL;
895
896 for(int i=1;i<(int_4)tokens.size();i++) {
897
898 AnyDataObj* mobj = omg.GetObj(tokens[i]);
899 if(mobj==NULL) {
900 cout<<"n_merge Error: unknow object"<<tokens[i]<<endl;
901 continue;
902 }
903 NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
904 if(nt1==NULL) {
905 cout<<"n_merge Error: "<<tokens[i]<<" not a NTuple"<<endl;
906 continue;
907 }
908 if(nt1->NEntry()==0) {
909 cout<<"n_merge Error: "<<tokens[i]<<" is empty"<<endl;
910 continue;
911 }
912 if(nt1->NVar()==0) {
913 cout<<"n_merge Error: "<<tokens[i]<<" has no variable"<<endl;
914 continue;
915 }
916
917 // create receiving ntuple if first pass
918 if(nt==NULL) {
919 nvar = nt1->NVar();
920 vector<string> sntvn;
921 for(int i=0;i<nvar;i++) sntvn.push_back(nt1->ColumnName(i));
922 char **ntvn = new char*[nvar];
923 for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(sntvn[i].c_str());
924 nt = new NTuple(nvar,ntvn);
925 delete [] ntvn;
926 }
927
928 // filling with current ntuple
929 int nvar1 = nt1->NVar();
930 int n = (nvar1>nvar)? nvar1: nvar;
931 r_8 *xnt1 = new r_8[n];
932 for(int i=0;i<n;i++) xnt1[i]=0.;
933 for(uint_4 iev=0;iev<(uint_4)nt1->NEntry();iev++) {
934 nt1->GetVecD(iev,xnt1);
935 nt->Fill(xnt1);
936 }
937 nfill++;
938 delete [] xnt1;
939
940 }
941
942 if(xnt!=NULL) delete [] xnt;
943 if(nt!=NULL) {
944 cout<<"n_merge: ntuple filled with "<<nfill
945 <<" ntuples, "<<nt->NEntry()<<" entries"<<endl;
946 if(nt->NEntry()>0) omg.AddObj(nt,nament); else delete nt;
947 }
948 return;
949}
950
951/* methode */
952void PAWExecutor::n_merge_col(vector<string>& tokens)
953{
954 if(tokens.size()<2) {
955 cerr<<"Usage: n/merge/col nt nt_1[,ext] nt_2[,ext] ... nt_n[,ext]"<<endl;
956 return;
957 }
958
959 NamedObjMgr omg;
960
961 // decodage des arguments
962 string nament = tokens[0];
963
964 // boucle sur les ntuples
965 uint_4 numnt=0, nentmax=0, ipos=0;
966 vector<string> varname;
967 vector<string> ntname;
968 vector<int> posfirst;
969 for(int i=1;i<(int_4)tokens.size();i++) {
970
971 // Decode NTuple name and extension
972 int_4 lc = tokens[i].size();
973 string ntn = tokens[i];
974 string ext = "";
975 int p = ntn.find(',');
976 if(p>=0 && p<lc) {
977 if(p<lc-1) ext = ntn.substr(p+1,lc-p-1);
978 if(p>0) ntn = ntn.substr(0,p); else ntn = "";
979 }
980 // choose the extentsion: given,automatic,none
981 if(ext.size()<1) {
982 char str[8]; sprintf(str,"_%d",numnt);
983 ext = str;
984 } else if(ext[0] == '!') {
985 ext = "";
986 }
987
988 AnyDataObj* mobj = omg.GetObj(ntn);
989 if(mobj==NULL) {
990 cout<<"n_merge_col Error: unknow object "<<ntn<<" ("<<tokens[i]<<")"<<endl;
991 continue;
992 }
993 NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
994 if(nt1==NULL) {
995 cout<<"n_merge_col Error: "<<ntn<<" not a NTuple"<<endl;
996 continue;
997 }
998 if(nt1->NEntry()==0) {
999 cout<<"n_merge_col Error: "<<ntn<<" is empty"<<endl;
1000 continue;
1001 }
1002 if(nt1->NVar()==0) {
1003 cout<<"n_merge_col Error: "<<ntn<<" has no variable"<<endl;
1004 continue;
1005 }
1006
1007 ntname.push_back(ntn);
1008 posfirst.push_back(ipos);
1009 if(nt1->NEntry()>(int_4)nentmax) nentmax = nt1->NEntry();
1010 for(int i=0;i<nt1->NVar();i++) {
1011 string str = nt1->ColumnName(i) + ext;
1012 varname.push_back(str);
1013 }
1014 numnt++;
1015 ipos += nt1->NVar();
1016 }
1017
1018 if(ntname.size()==0) {
1019 cout<<"n_merge_col Error: no ntuple found in argument list"<<endl;
1020 return;
1021 } else cout<<"Number of NTuple found "<<ntname.size()
1022 <<", Number of variables "<<varname.size()
1023 <<", Max_entries "<<nentmax<<endl;
1024
1025 // create receiving ntuple
1026 int nvar = varname.size();
1027 char **ntvn = new char*[nvar];
1028 for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(varname[i].c_str());
1029 NTuple * nt = new NTuple(nvar,ntvn);
1030 //for(int i=0;i<(int_4)ntname.size();i++) cout<<ntname[i]<<" , p="<<posfirst[i]<<endl;
1031 delete [] ntvn;
1032
1033 // filling with current ntuple
1034 r_8 *xnt = new r_8[nvar];
1035 for(int iev=0;iev<(int_4)nentmax;iev++) {
1036 for(int i=0;i<nvar;i++) xnt[i]=0.;
1037 for(int i=0;i<(int_4)ntname.size();i++) {
1038 AnyDataObj* mobj = omg.GetObj(ntname[i]);
1039 NTuple* nt1 = dynamic_cast<NTuple*>(mobj);
1040 if(iev < nt1->NEntry()) nt1->GetVecD(iev,&xnt[posfirst[i]]);
1041 }
1042 nt->Fill(xnt);
1043 }
1044
1045 // Ending and saving merging ntuple columns
1046 delete [] xnt;
1047 cout<<"n_merge_col: ntuple filled with "<<ntname.size()
1048 <<" ntuples, "<<nt->NEntry()<<" entries"<<endl;
1049 if(nt->NEntry()>0) omg.AddObj(nt,nament); else delete nt;
1050
1051 return;
1052}
1053
1054/* methode */
1055void PAWExecutor::n_copy(vector<string>& tokens)
1056{
1057 if(tokens.size()<2) {
1058 cerr<<"Usage: n/copy ntnew nt [vname1 vname2 vname3 ...]"<<endl;
1059 return;
1060 }
1061
1062 NamedObjMgr omg;
1063
1064 // Le NTuple de depart
1065 AnyDataObj* mobj = omg.GetObj(tokens[1]);
1066 if(mobj==NULL) {
1067 cout<<"n_copy Error: unknow object"<<tokens[1]<<endl;
1068 return;
1069 }
1070 NTuple* nt = dynamic_cast<NTuple*>(mobj);
1071 if(nt==NULL) {
1072 cout<<"n_copy Error: "<<tokens[1]<<" not a NTuple"<<endl;
1073 return;
1074 }
1075 if(nt->NEntry()==0) {
1076 cout<<"n_copy Error: "<<tokens[1]<<" is empty"<<endl;
1077 return;
1078 }
1079 if(nt->NVar()==0) {
1080 cout<<"n_copy Error: "<<tokens[1]<<" has no variable"<<endl;
1081 return;
1082 }
1083
1084 // Les variables a copier
1085 int nvar; vector<string> varname; vector<int> ivar;
1086 if(tokens.size()==2) {
1087 nvar = nt->NVar();
1088 for(int i=0;i<nvar;i++) {
1089 ivar.push_back(i);
1090 varname.push_back(nt->NomIndex(i));
1091 }
1092 } else {
1093 nvar = tokens.size()-2;
1094 for(int i=0;i<nvar;i++) {
1095 int iv = nt->IndexNom(tokens[i+2].c_str());
1096 if(iv<0) {
1097 cout<<"n_copy Error: unkown variable "<<tokens[i+2].c_str()<<endl;
1098 return;
1099 }
1100 ivar.push_back(iv);
1101 varname.push_back(tokens[i+2]);
1102 }
1103 }
1104 cout<<"coping variables: ";
1105 for(int i=0;i<nvar;i++) cout<<" "<<varname[i];
1106 cout<<endl;
1107
1108 // Creation du nouveau NTuple
1109 char **ntvn = new char*[nvar];
1110 for(int i=0;i<nvar;i++) ntvn[i] = const_cast<char *>(varname[i].c_str());
1111 NTuple ntnew(nvar,ntvn);
1112 delete [] ntvn;
1113
1114 // Copie du NTuple
1115 r_8 *xnt = new r_8[nt->NVar()];
1116 r_8 *xntnew = new r_8[nvar];
1117 for(int iev=0;iev<nt->NEntry();iev++) {
1118 nt->GetVecD(iev,xnt);
1119 for(int i=0;i<nvar;i++) xntnew[i] = xnt[ivar[i]];
1120 ntnew.Fill(xntnew);
1121 }
1122 delete [] xnt; delete [] xntnew;
1123
1124 // Adding new NTuple
1125 omg.AddObj(ntnew,tokens[0]);
1126
1127 return;
1128}
1129
1130/* methode */
1131void PAWExecutor::h_integ(vector<string>& tokens)
1132// Pour remplacer le contenu d'un histo 1D par son integrale
1133{
1134if(tokens.size()<1)
1135 {cout<<"Usage: h/integ nameh1d [norm]"<<endl; return;}
1136NamedObjMgr omg;
1137AnyDataObj* mobj = omg.GetObj(tokens[0]);
1138if(mobj==NULL)
1139 {cout<<"PAWExecutor::h_integ Error: unknow object"<<tokens[0]<<endl;
1140 return;}
1141r_8 norm = 1.;
1142if(tokens.size()>=2) norm = atof(tokens[1].c_str());
1143// attention: dynamic_cast<Histo*>(HProf)=Vrai!
1144Histo* h1 = dynamic_cast<Histo*>(mobj);
1145HProf* hp = dynamic_cast<HProf*>(mobj);
1146if(hp || !h1)
1147 {cout<<"PAWExecutor::h_integ Error: "<<tokens[0]<<" not an Histo"<<endl;
1148 return;}
1149h1->HInteg(norm);
1150}
1151
1152/* methode */
1153void PAWExecutor::v_integ(vector<string>& tokens)
1154// Pour remplacer le contenu d'un TVector<r_8> par son integrale
1155// normalisee a norm:
1156// Si norm <= 0 : pas de normalisation, integration seule
1157{
1158if(tokens.size()<1)
1159 {cout<<"Usage: v/integ namevec [norm]"<<endl; return;}
1160NamedObjMgr omg;
1161AnyDataObj* mobj = omg.GetObj(tokens[0]);
1162if(mobj==NULL)
1163 {cout<<"PAWExecutor::v_integ Error: unknow object"<<tokens[0]<<endl;
1164 return;}
1165r_8 norm=-1.;
1166if(tokens.size()>=2) norm = atof(tokens[1].c_str());
1167#ifdef SANS_EVOLPLANCK
1168Vector* v = dynamic_cast<Vector*>(mobj);
1169#else
1170TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
1171#endif
1172
1173if(!v)
1174 {cout<<"PAWExecutor::v_integ Error: "<<tokens[0]<<" not a TVector/Vector"<<endl;
1175 return;}
1176
1177#ifdef SANS_EVOLPLANCK
1178uint_4 n = v->NElts();
1179#else
1180uint_4 n = v->Size();
1181#endif
1182
1183if(n==0)
1184 {cout<<"PAWExecutor::v_integ Error: "<<tokens[0]<<" is an empty vector"<<endl;
1185 return;}
1186if(n>1) for(uint_4 i=1;i<n;i++) (*v)(i)+=(*v)(i-1);
1187if(norm<=0. || (*v)(n-1)==0.) return;
1188norm /= (*v)(n-1);
1189for(uint_4 i=0;i<n;i++) (*v)(i) *= norm;
1190}
1191
1192/* methode */
1193void PAWExecutor::v_sort(vector<string>& tokens)
1194// Pour sort in-place d'un vecteur ascendant ou descendant
1195{
1196 if(tokens.size()<1)
1197 {cout<<"Usage: v/sort namevec [+1/-1]"<<endl; return;}
1198 NamedObjMgr omg;
1199 AnyDataObj* mobj = omg.GetObj(tokens[0]);
1200 if(mobj==NULL)
1201 {cout<<"PAWExecutor::v_sort Error: unknow object"<<tokens[0]<<endl;
1202 return;}
1203
1204 bool up = true;
1205 if(tokens.size()>=2) if(atoi(tokens[1].c_str())<0) up=false;
1206
1207 {
1208 TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
1209 if(v) {
1210 if(v->Size()==0) return;
1211 if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_r8);
1212 else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_r8);
1213 cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
1214 return;
1215 }
1216 }
1217
1218 {
1219 TVector<r_4>* v = dynamic_cast<TVector<r_4>*>(mobj);
1220 if(v) {
1221 if(v->Size()==0) return;
1222 if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_r4);
1223 else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_r4);
1224 cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
1225 return;
1226 }
1227 }
1228
1229 {
1230 TVector<int_4>* v = dynamic_cast<TVector<int_4>*>(mobj);
1231 if(v) {
1232 if(v->Size()==0) return;
1233 if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_i4);
1234 else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_i4);
1235 cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
1236 return;
1237 }
1238 }
1239
1240 {
1241 TVector<uint_8>* v = dynamic_cast<TVector<uint_8>*>(mobj);
1242 if(v) {
1243 if(v->Size()==0) return;
1244 if(up) stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_up_u8);
1245 else stable_sort(v->Data(), v->Data()+v->Size(), __vector_sort_down_u8);
1246 cout<<"median is: "<<0.5*((*v)(v->Size()/2)+(*v)((v->Size()-1)/2))<<endl;
1247 return;
1248 }
1249 }
1250
1251 cout<<"PAWExecutor::v_sort Error: "<<tokens[0]
1252 <<" not a TVector or not a supported TVector<TYPE>"<<endl;
1253 return;
1254
1255}
1256
1257/* methode */
1258void PAWExecutor::h_deriv(vector<string>& tokens)
1259// Pour remplacer le contenu d'un histo 1D par sa derivee
1260{
1261if(tokens.size()<1)
1262 {cout<<"Usage: h/deriv nameh1d"<<endl; return;}
1263NamedObjMgr omg;
1264AnyDataObj* mobj = omg.GetObj(tokens[0]);
1265if(mobj==NULL)
1266 {cout<<"PAWExecutor::h_deriv Error: unknow object"<<tokens[0]<<endl;
1267 return;}
1268// attention: dynamic_cast<Histo*>(HProf)=Vrai!
1269Histo* h1 = dynamic_cast<Histo*>(mobj);
1270HProf* hp = dynamic_cast<HProf*>(mobj);
1271if(hp || !h1)
1272 {cout<<"PAWExecutor::h_deriv Error: "<<tokens[0]<<" not an Histo"<<endl;
1273 return;}
1274h1->HDeriv();
1275}
1276
1277/* methode */
1278void PAWExecutor::v_deriv(vector<string>& tokens)
1279// Pour remplacer le contenu d'un TVector<r_8> par sa derivee
1280// deriv_option = -1 replace v[i] with v[i]-v[i-1]
1281// = 0 replace v[i] with (v[i+1]-v[i-1])/2 (default)
1282// = +1 replace v[i] with v[i+1]-v[i]
1283{
1284if(tokens.size()<1)
1285 {cout<<"Usage: v/deriv namevec [deriv_option(-1,0,+1)]"<<endl; return;}
1286NamedObjMgr omg;
1287AnyDataObj* mobj = omg.GetObj(tokens[0]);
1288if(mobj==NULL)
1289 {cout<<"PAWExecutor::v_deriv Error: unknow object"<<tokens[0]<<endl;
1290 return;}
1291int_4 deriv_option = 0;
1292if(tokens.size()>=2) deriv_option = atoi(tokens[1].c_str());
1293
1294#ifdef SANS_EVOLPLANCK
1295Vector* v = dynamic_cast<Vector*>(mobj);
1296#else
1297TVector<r_8>* v = dynamic_cast<TVector<r_8>*>(mobj);
1298#endif
1299
1300if(!v)
1301 {cout<<"PAWExecutor::v_deriv Error: "<<tokens[0]<<" not a TVector/Vector"<<endl;
1302 return;}
1303
1304
1305#ifdef SANS_EVOLPLANCK
1306uint_4 n = v->NElts();
1307#else
1308uint_4 n = v->Size();
1309#endif
1310
1311if(n==0)
1312 {cout<<"PAWExecutor::v_deriv Error: "<<tokens[0]<<" is an empty vector"<<endl;
1313 return;}
1314if(n<=1) return;
1315
1316#ifdef SANS_EVOLPLANCK
1317Vector vsave(*v);
1318#else
1319TVector<r_8> vsave(*v,false);
1320#endif
1321
1322if(deriv_option<0) {
1323 for(uint_4 i=1;i<n;i++) (*v)(i) = vsave(i)-vsave(i-1);
1324 (*v)(0) = (*v)(1);
1325} else if(deriv_option>0) {
1326 for(uint_4 i=0;i<n-1;i++) (*v)(i) = vsave(i+1)-vsave(i);
1327 (*v)(n-1) = (*v)(n-2);
1328} else {
1329 for(uint_4 i=1;i<n-1;i++) (*v)(i) = (vsave(i+1)-vsave(i-1))/2.;
1330 (*v)(0) = vsave(1)-vsave(0);
1331 (*v)(n-1) = vsave(n-1)-vsave(n-2);
1332}
1333}
1334
1335
1336/* methode */
1337void PAWExecutor::h_rebin(vector<string>& tokens)
1338// Pour re-binner un histogramme 1D
1339{
1340if(tokens.size()<2)
1341 {cout<<"Usage: h/rebin nameh1d nbin"<<endl; return;}
1342NamedObjMgr omg;
1343AnyDataObj* mobj = omg.GetObj(tokens[0]);
1344if(mobj==NULL)
1345 {cout<<"PAWExecutor::h_rebin Error: unknow object"<<tokens[0]<<endl;
1346 return;}
1347int_4 nbin = atoi(tokens[1].c_str());
1348Histo* h1 = dynamic_cast<Histo*>(mobj);
1349// Ca marche aussi pour les HProf, HRebin a ete passe virtuel
1350//HProf* hp = dynamic_cast<HProf*>(mobj); if(hp || !h1)
1351if(!h1)
1352 {cout<<"PAWExecutor::h_rebin Error: "<<tokens[0]<<" not an Histo/HProf"<<endl;
1353 return;}
1354h1->HRebin(nbin);
1355}
1356
1357/* methode */
1358void PAWExecutor::h_cadd(vector<string>& tokens)
1359// Additionne une constante a un histogramme
1360{
1361if(tokens.size()<2)
1362 {cout<<"Usage: h/cadd nameh1d val"<<endl; return;}
1363NamedObjMgr omg;
1364AnyDataObj* mobj = omg.GetObj(tokens[0]);
1365if(mobj==NULL)
1366 {cout<<"PAWExecutor::h_cadd Error: unknow object"<<tokens[0]<<endl;
1367 return;}
1368r_8 val = atof(tokens[1].c_str());
1369Histo* h1 = dynamic_cast<Histo*>(mobj);
1370HProf* hp = dynamic_cast<HProf*>(mobj);
1371Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
1372Vector* v = dynamic_cast<Vector*>(mobj);
1373Matrix* m = dynamic_cast<Matrix*>(mobj);
1374if(h1 && !hp) *h1 += val;
1375else if(h2) *h2 += val;
1376else if(v) *v += val;
1377else if(m) *m += val;
1378else cout<<"PAWExecutor::h_cadd Error: not implemented for "<<typeid(*mobj).name()<<endl;
1379}
1380
1381/* methode */
1382void PAWExecutor::h_cmult(vector<string>& tokens)
1383// Multiplie un histogramme par une constante
1384{
1385if(tokens.size()<2)
1386 {cout<<"Usage: h/cmult nameh1d val"<<endl; return;}
1387NamedObjMgr omg;
1388AnyDataObj* mobj = omg.GetObj(tokens[0]);
1389if(mobj==NULL)
1390 {cout<<"PAWExecutor::h_cmult Error: unknow object"<<tokens[0]<<endl;
1391 return;}
1392r_8 val = atof(tokens[1].c_str());
1393Histo* h1 = dynamic_cast<Histo*>(mobj);
1394HProf* hp = dynamic_cast<HProf*>(mobj);
1395Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
1396Vector* v = dynamic_cast<Vector*>(mobj);
1397Matrix* m = dynamic_cast<Matrix*>(mobj);
1398if(h1 && !hp) *h1 *= val;
1399else if(h2) *h2 *= val;
1400else if(v) *v += val;
1401else if(m) *m += val;
1402else cout<<"PAWExecutor::h_mult Error: not implemented for "<<typeid(*mobj).name()<<endl;
1403}
1404
1405/* methode */
1406void PAWExecutor::h_oper(vector<string>& tokens)
1407// Equivalent h/oper/add sub,mul,div de paw
1408// Operation entre 2 histogrammes ou 2 vecteurs ou 2 matrices
1409// h/oper @ h1 h2 hres
1410// hres = h1 @ h2 with @ = (+,-,*,/)
1411// Pour les vecteurs et les matrices, il sagit d'operations elements a elements
1412{
1413if(tokens.size()<4)
1414 {cout<<"Usage: n/oper @ h1 h2 hres with @=(+,-,*,/,@)"<<endl;
1415 return;}
1416
1417// Decode arguments
1418const char * oper = tokens[0].c_str();
1419if( oper[0]!='+' && oper[0]!='-' && oper[0]!='*' && oper[0]!='/' )
1420 {cout<<"PAWExecutor::h_oper Error: unknow operation "<<oper<<endl;
1421 return;}
1422string h1name = tokens[1];
1423string h2name = tokens[2];
1424string h3name = tokens[3];
1425
1426// Get objects
1427NamedObjMgr omg;
1428AnyDataObj* mobjh1 = omg.GetObj(h1name);
1429AnyDataObj* mobjh2 = omg.GetObj(h2name);
1430if( mobjh1==NULL || mobjh2==NULL )
1431 {cout<<"PAWExecutor::h_oper Error: unknow object(s) "<<h1name<<" or "<<h2name<<endl;
1432 return;}
1433AnyDataObj* mobjh3 = omg.GetObj(h3name);
1434
1435// Operations on HProf, only + is working
1436if( dynamic_cast<HProf*>(mobjh1) != NULL ) {
1437 if( oper[0]!='+' )
1438 { cout<<"PAWExecutor::h_oper Error: operation "<<oper
1439 <<" not implemented for HProf"<<endl; return;}
1440 if( dynamic_cast<HProf*>(mobjh2) == NULL )
1441 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
1442 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
1443 return;}
1444 HProf* h1 =(HProf*) mobjh1;
1445 HProf* h2 =(HProf*) mobjh2;
1446 if( h1->NBins() != h2->NBins() )
1447 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
1448 <<h1->NBins()<<" "<<h2->NBins()<<endl; return;}
1449 HProf* h3 = NULL;
1450 if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
1451 {h3 = new HProf(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
1452 h3 = dynamic_cast<HProf*>(mobjh3);
1453 if(h3 == NULL) // ce n'est pas un HProf
1454 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
1455 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
1456 if(h1->NBins() != h3->NBins())
1457 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
1458 <<h1->NBins()<<" "<<h3->NBins()<<endl; return;}
1459 *h3 = *h1 + *h2;
1460 h3->UpdateHisto();
1461
1462// Operations on Histo
1463} else if( dynamic_cast<Histo*>(mobjh1) != NULL ) {
1464 if( dynamic_cast<Histo*>(mobjh2) == NULL )
1465 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
1466 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
1467 return;}
1468 Histo* h1 =(Histo*) mobjh1;
1469 Histo* h2 =(Histo*) mobjh2;
1470 if( h1->NBins() != h2->NBins() )
1471 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
1472 <<h1->NBins()<<" "<<h2->NBins()<<endl; return;}
1473 Histo* h3 = NULL;
1474 if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
1475 {h3 = new Histo(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
1476 h3 = dynamic_cast<Histo*>(mobjh3);
1477 if(h3 == NULL) // ce n'est pas un Histo
1478 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
1479 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
1480 if(h1->NBins() != h3->NBins())
1481 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
1482 <<h1->NBins()<<" "<<h3->NBins()<<endl; return;}
1483 if( oper[0]=='+') *h3 = *h1 + *h2;
1484 else if( oper[0]=='-') *h3 = *h1 - *h2;
1485 else if( oper[0]=='*') *h3 = *h1 * *h2;
1486 else if( oper[0]=='/') *h3 = *h1 / *h2;
1487
1488// Operations on Histo2D
1489} else if( dynamic_cast<Histo2D*>(mobjh1) != NULL ) {
1490 if( dynamic_cast<Histo2D*>(mobjh2) == NULL )
1491 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
1492 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
1493 return;}
1494 Histo2D* h1 =(Histo2D*) mobjh1;
1495 Histo2D* h2 =(Histo2D*) mobjh2;
1496 if( h1->NBinX() != h2->NBinX() || h1->NBinY() != h2->NBinY() )
1497 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
1498 <<h1->NBinX()<<","<<h1->NBinY()<<" "
1499 <<h2->NBinX()<<","<<h2->NBinY()<<endl; return;}
1500 Histo2D* h3 = NULL;
1501 if( mobjh3 == NULL ) // l'objet n'existe pas, on le cree
1502 {h3 = new Histo2D(*h1); h3->Zero(); omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);}
1503 h3 = dynamic_cast<Histo2D*>(mobjh3);
1504 if(h3 == NULL) // ce n'est pas un Histo2D
1505 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
1506 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
1507 if( h1->NBinX() != h3->NBinX() || h1->NBinY() != h3->NBinY() )
1508 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
1509 <<h1->NBinX()<<","<<h1->NBinY()<<" "
1510 <<h3->NBinX()<<","<<h3->NBinY()<<endl; return;}
1511 if( oper[0]=='+') *h3 = *h1 + *h2;
1512 else if( oper[0]=='-') *h3 = *h1 - *h2;
1513 else if( oper[0]=='*') *h3 = *h1 * *h2;
1514 else if( oper[0]=='/') *h3 = *h1 / *h2;
1515
1516// Operations on Vector
1517} else if( dynamic_cast<Vector*>(mobjh1) != NULL ) {
1518 if( dynamic_cast<Vector*>(mobjh2) == NULL )
1519 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
1520 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
1521 return;}
1522 Vector* h1 =(Vector*) mobjh1;
1523 Vector* h2 =(Vector*) mobjh2;
1524 if( h1->NElts() != h2->NElts() )
1525 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
1526 <<h1->NElts()<<" "<<h2->NElts()<<endl; return;}
1527 Vector* h3 = NULL;
1528 if( mobjh3 == NULL ) { // l'objet n'existe pas, on le cree
1529#ifdef SANS_EVOLPLANCK
1530 h3 = new Vector(*h1);
1531#else
1532 h3 = new Vector(*h1,false);
1533#endif
1534 *h3 = 0.; omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);
1535 }
1536 h3 = dynamic_cast<Vector*>(mobjh3);
1537 if(h3 == NULL) // ce n'est pas un Vector
1538 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
1539 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
1540 if(h1->NElts() != h3->NElts())
1541 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
1542 <<h1->NElts()<<" "<<h3->NElts()<<endl; return;}
1543 if( oper[0]=='+') *h3 = *h1 + *h2;
1544 else if( oper[0]=='-') *h3 = *h1 - *h2;
1545#ifdef SANS_EVOLPLANCK
1546 else if(oper[0]=='*' || oper[0]=='/')
1547 cout<<"PAWExecutor::h_oper Error: operation "<<oper[0]
1548 <<" not implemented for Vector in Peida"<<endl;
1549#else
1550 else if( oper[0]=='*') {h3->Clone(*h1); h3->MulElt(*h2,*h3);}
1551 else if( oper[0]=='/') {h3->Clone(*h1); h3->DivElt(*h2,*h3,false,true);}
1552#endif
1553
1554// Operations on Matrix
1555} else if( dynamic_cast<Matrix*>(mobjh1) != NULL ) {
1556 if( dynamic_cast<Matrix*>(mobjh2) == NULL )
1557 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h2\n"
1558 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh2).name()<<endl;
1559 return;}
1560 Matrix* h1 =(Matrix*) mobjh1;
1561 Matrix* h2 =(Matrix*) mobjh2;
1562 if( h1->NRows() != h2->NRows() || h1->NCol() != h2->NCol() )
1563 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h2 "
1564 <<h1->NRows()<<","<<h1->NCol()<<" "
1565 <<h2->NRows()<<","<<h2->NCol()<<endl; return;}
1566 Matrix* h3 = NULL;
1567 if( mobjh3 == NULL ) { // l'objet n'existe pas, on le cree
1568#ifdef SANS_EVOLPLANCK
1569 h3 = new Matrix(*h1);
1570#else
1571 h3 = new Matrix(*h1,false);
1572#endif
1573 *h3 = 0.; omg.AddObj(h3,h3name); mobjh3 = omg.GetObj(h3name);
1574 }
1575 h3 = dynamic_cast<Matrix*>(mobjh3);
1576 if(h3 == NULL) // ce n'est pas un Matrix
1577 {cout<<"PAWExecutor::h_oper Error: type mismatch between h1 and h3\n"
1578 <<typeid(*mobjh1).name()<<" , "<<typeid(*mobjh3).name()<<endl; return;}
1579 if( h1->NRows() != h3->NRows() || h1->NCol() != h3->NCol() )
1580 {cout<<"PAWExecutor::h_oper Error: size mismatch between h1, h3 "
1581 <<h1->NRows()<<","<<h1->NCol()<<" "
1582 <<h3->NRows()<<","<<h3->NCol()<<endl; return;}
1583 if( oper[0]=='+') *h3 = *h1 + *h2;
1584 else if( oper[0]=='-') *h3 = *h1 - *h2;
1585#ifdef SANS_EVOLPLANCK
1586 else if(oper[0]=='*' || oper[0]=='/')
1587 cout<<"PAWExecutor::h_oper Error: operation "<<oper[0]
1588 <<" not implemented for Vector in Peida"<<endl;
1589#else
1590 else if( oper[0]=='*') {h3->Clone(*h1); h3->MulElt(*h2,*h3);}
1591 else if( oper[0]=='/') {h3->Clone(*h1); h3->DivElt(*h2,*h3,false,true);}
1592#endif
1593
1594// Doesn't work for other objects
1595} else {
1596 cout<<"PAWExecutor::h_oper Error: not implemented for "
1597 <<typeid(*mobjh1).name()<<endl;
1598 return;
1599}
1600
1601return;
1602}
1603
1604/* methode */
1605void PAWExecutor::h_plot_2d(vector<string>& tokens)
1606// plot for 2D histogramme: plot histo, bandx/y, slicex/y or projx/y
1607{
1608if(tokens.size()<1)
1609 {cout<<"Usage: h/plot/2d nameh2d to_plot [n/s] [dopt]"<<endl; return;}
1610string proj = "h"; if(tokens.size()>1) proj = tokens[1];
1611NamedObjMgr omg;
1612AnyDataObj* mobj = omg.GetObj(tokens[0]);
1613if(mobj==NULL)
1614 {cout<<"PAWExecutor::h_plot_2d Error: unknow object"<<tokens[0]<<endl;
1615 return;}
1616Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
1617if(!h2)
1618 {cout<<"PAWExecutor::h_plot_2d Error: "<<tokens[0]<<" not an Histo2D"<<endl;
1619 return;}
1620
1621Histo* h1p = NULL; string nametoplot = "/autoc/h_plot_2d_h1";
1622
1623string dopt = ""; if(tokens.size()>=3) dopt = tokens[2];
1624if(proj == "show") {
1625 h2->ShowProj();
1626 h2->ShowBand(2);
1627 h2->ShowSli(2);
1628 return;
1629} else if(proj == "h") {
1630 nametoplot = tokens[0];
1631} else if(proj == "px") {
1632 if((h1p=h2->HProjX())) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1633 else {h2->ShowProj(); return;}
1634} else if(proj == "py") {
1635 if((h1p=h2->HProjY())) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1636 else {h2->ShowProj(); return;}
1637} else {
1638 if(tokens.size()<3)
1639 {cout<<"Usage: h/plot/2d nameh2d bx/by/sx/sy n [dopt]"<<endl; return;}
1640 int_4 n = atoi(tokens[2].c_str());
1641 dopt = ""; if(tokens.size()>=4) dopt = tokens[3];
1642 if(proj == "bx") {
1643 if((h1p=h2->HBandX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1644 else {h2->ShowBand(); return;}
1645 } else if(proj == "by") {
1646 if((h1p=h2->HBandY(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1647 else {h2->ShowBand(); return;}
1648 } else if(proj == "sx") {
1649 if((h1p=h2->HSliX(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1650 else {h2->ShowSli(); return;}
1651 } else if(proj == "sy") {
1652 if((h1p=h2->HSliY(n))) {Histo* h1=new Histo(*h1p); omg.AddObj(h1,nametoplot);}
1653 else {h2->ShowSli(); return;}
1654 }
1655}
1656
1657omg.DisplayObj(nametoplot,dopt);
1658}
1659
1660/* methode */
1661int_4 PAWExecutor::decodepawstring(string tokens,string& nameobj
1662 ,string& xexp,string& yexp,string& zexp)
1663// Decodage general de "nameobj.xexp"
1664// "nameobj.yexp%xexp"
1665// "nameobj.zexp%yexp%xexp"
1666// Return: nombre de variables trouvees, -1 si probleme
1667{
1668nameobj = ""; xexp= ""; yexp= ""; zexp= "";
1669
1670int_4 lt = (int) tokens.length();
1671if(lt<=0) return -1;
1672
1673// decodage de la chaine de type PAW.
1674char *str = new char[lt+2];
1675strcpy(str,tokens.c_str()); strip(str,'B',' '); lt = strlen(str);
1676//cout<<"chaine1["<<lt<<"] :"<<str<<":"<<endl;
1677char *c[3] = {NULL,NULL,NULL};
1678int_4 i, np=0; bool namefound = false;
1679for(i=0;i<lt;i++) {
1680 if(!namefound && str[i]=='.') {
1681 str[i]='\0';
1682 namefound=true;
1683 c[np] = str+i+1; np++;
1684 }
1685 if( namefound && str[i]=='%') {
1686 str[i]='\0';
1687 if(np<3) {c[np] = str+i+1; np++;}
1688 }
1689}
1690//cout<<"chaine2 :"; for(i=0;i<lt;i++) cout<<str[i]; cout<<":"<<endl;
1691
1692// Remplissage du nom et des variables
1693nameobj = str;
1694if(np==1) xexp=c[0];
1695if(np==2) {yexp=c[0]; xexp=c[1];}
1696if(np==3) {zexp=c[0]; yexp=c[1]; xexp=c[2];}
1697//cout<<"pawstring str,c[0-2] "<<str<<" "<<c[0]<<" "<<c[1]<<" "<<c[2]<<endl;
1698delete [] str;
1699
1700// Comptage des variables
1701np = -1;
1702if(nameobj.length()>0)
1703 {np = 0; if(xexp.length()>0)
1704 {np++; if(yexp.length()>0)
1705 {np++; if(zexp.length()>0) np++;}}}
1706//cout<<"pawstring["<<np<<"] name="<<nameobj
1707// <<" xexp="<<xexp<<" yexp="<<yexp<<" zexp="<<zexp<<endl;
1708return np;
1709}
1710
1711/* methode */
1712void PAWExecutor::h_put_vec(vector<string>& tokens)
1713// Pour remplir un histo avec le contenu d'un vecteur ou d'une matrice
1714// L'histogramme doit deja exister. Le nombre de bins remplit
1715// depend des tailles respectives des 2 objets.
1716// tokens[2] = "cont" : on remplit les valeurs de l'histogramme (ou "!")
1717// = "err2" : on remplit les erreurs de l'histogramme
1718{
1719if(tokens.size()<2)
1720 {cout<<"Usage: h/put_vec namehisto namevector"<<endl;
1721 return;}
1722string toput = "cont"; if(tokens.size()>2) toput = tokens[2];
1723if(toput=="!") toput = "cont";
1724if(toput!="cont" && toput!="err2")
1725 {cout<<"PAWExecutor::h_put_vec Error: unknow filling "<<toput<<endl;
1726 return;}
1727string hname = tokens[0];
1728string vname = tokens[1];
1729
1730// Get objects
1731NamedObjMgr omg;
1732AnyDataObj* mobjh = omg.GetObj(hname);
1733AnyDataObj* mobjv = omg.GetObj(vname);
1734if( mobjh==NULL || mobjv==NULL )
1735 {cout<<"PAWExecutor::h_put_vec Error: unknow object(s) "<<hname<<" or "<<vname<<endl;
1736 return;}
1737
1738// Fill Histo from Vector
1739if(typeid(*mobjh) == typeid(Histo) && typeid(*mobjv) == typeid(Vector)) {
1740 Histo* h = dynamic_cast<Histo*>(mobjh);
1741 Vector* v = dynamic_cast<Vector*>(mobjv);
1742 if(toput=="cont") h->PutValue(*v);
1743 else if(toput=="err2") h->PutError2(*v);
1744
1745// Fill Histo2D from Matrix
1746} else if(typeid(*mobjh) == typeid(Histo2D) && typeid(*mobjv) == typeid(Matrix)) {
1747 Histo2D* h = dynamic_cast<Histo2D*>(mobjh);
1748 Matrix* v = dynamic_cast<Matrix*>(mobjv);
1749 if(toput=="cont") h->PutValue(*v);
1750 else if(toput=="err2") h->PutError2(*v);
1751
1752} else {
1753 cout<<"PAWExecutor::h_put_vec Error: type mismatch between histogram and vector/matrix\n"
1754 <<typeid(*mobjh).name()<<" , "<<typeid(*mobjv).name()<<endl;
1755 return;
1756}
1757
1758}
1759
1760/* methode */
1761void PAWExecutor::h_get_vec(vector<string>& tokens)
1762// Pour copier un histo dans un vecteur ou une matrice
1763// Si le vecteur (matrice) n'existe pas, il est cree.
1764// Le vecteur ou la matrice est re-dimensionne correctement.
1765// tokens[2] = "cont" : on copie les valeurs de l'histogramme (ou "!")
1766// = "err2" : on copie les erreurs de l'histogramme
1767// = "absc" : on copie les erreurs de l'histogramme (1D only)
1768// tokens[3[,4]] = show : show available projections for Histo2D
1769// px : get X projection
1770// py : get Y projection
1771// bx n : get X band number n
1772// by n : get Y band number n
1773// sx n : get X slice number n
1774// sy n : get Y slice number n
1775{
1776if(tokens.size()<2)
1777 {cout<<"Usage: h/get_vec namehisto namevector"<<endl;
1778 return;}
1779string toget = "cont"; if(tokens.size()>2) toget = tokens[2];
1780if(toget=="!") toget = "cont";
1781if(toget!="cont" && toget!="err2" && toget!="absc")
1782 {cout<<"PAWExecutor::h_get_vec Error: unknow filling "<<toget<<endl;
1783 return;}
1784string proj = "h"; if(tokens.size()>3) proj = tokens[3];
1785int_4 nproj = -1; if(tokens.size()>4) nproj = atoi(tokens[4].c_str());
1786string hname = tokens[0];
1787string vname = tokens[1];
1788
1789// Get objects
1790NamedObjMgr omg;
1791AnyDataObj* mobjh = omg.GetObj(hname);
1792AnyDataObj* mobjv = omg.GetObj(vname);
1793if( mobjh==NULL)
1794 {cout<<"PAWExecutor::h_put_vec Error: unknow object(s) "<<hname<<endl;
1795 return;}
1796
1797// Copy Histo/Histo2D/Projection to Vector/Matrix
1798Histo* h = dynamic_cast<Histo*>(mobjh);
1799Histo2D* h2 = dynamic_cast<Histo2D*>(mobjh);
1800if( h != NULL ) { // Histo ou HProf
1801 h->UpdateHisto(); // pour le cas ou c'est un HProf
1802 Vector* v = NULL;
1803 if(mobjv==NULL) // le vecteur n'existe pas
1804 {v = new Vector(1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
1805 if(typeid(*mobjv) != typeid(Vector))
1806 {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo/HProf and vector\n"
1807 <<typeid(*mobjv).name()<<endl; return;}
1808 v = dynamic_cast<Vector*>(mobjv);
1809 if(toget=="cont") h->GetValue(*v);
1810 else if(toget=="err2") h->GetError2(*v);
1811 else if(toget=="absc") h->GetAbsc(*v);
1812
1813} else if( h2 != NULL) { // Histo2D
1814
1815 if(proj == "h") { // On veut les valeurs de l'histogramme 2D
1816 Matrix* v = NULL;
1817 if(mobjv==NULL) // la matrice n'existe pas
1818 {v = new Matrix(1,1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
1819 if(typeid(*mobjv) != typeid(Matrix))
1820 {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo2D and matrix\n"
1821 <<typeid(*mobjv).name()<<endl; return;}
1822 v = dynamic_cast<Matrix*>(mobjv);
1823 if(toget=="cont") h2->GetValue(*v);
1824 else if(toget=="err2") h2->GetError2(*v);
1825 else
1826 {cout<<"PAWExecutor::h_get_vec Error: option "<<toget<<" not valid for Histo2D"<<endl;
1827 return;}
1828
1829 } else { // On veut les valeurs d'une projection de l'histo 2D
1830 h = NULL;
1831 if(proj == "show") {h2->ShowProj(); h2->ShowBand(2); h2->ShowSli(2);}
1832 else if(proj == "px") h = h2->HProjX();
1833 else if(proj == "py") h = h2->HProjY();
1834 else if(proj == "bx") h = h2->HBandX(nproj);
1835 else if(proj == "by") h = h2->HBandY(nproj);
1836 else if(proj == "sx") h = h2->HSliX(nproj);
1837 else if(proj == "sy") h = h2->HSliY(nproj);
1838 if(h==NULL)
1839 {cout<<"PAWExecutor::h_get_vec Error: unknown projection "<<proj
1840 <<" number "<<nproj<<endl; return;}
1841 Vector* v = NULL;
1842 if(mobjv==NULL) // le vecteur n'existe pas
1843 {v = new Vector(1); omg.AddObj(v,vname); mobjv = omg.GetObj(vname);}
1844 if(typeid(*mobjv) != typeid(Vector))
1845 {cout<<"PAWExecutor::h_get_vec Error: type mismatch between Histo2D "<<proj
1846 <<" and vector\n"<<typeid(*mobjv).name()<<endl; return;}
1847 v = dynamic_cast<Vector*>(mobjv);
1848 if(toget=="cont") h->GetValue(*v);
1849 else if(toget=="err2") h->GetError2(*v);
1850 else if(toget=="absc") h->GetAbsc(*v);
1851 }
1852
1853} else {
1854 cout<<"PAWExecutor::h_get_vec Error: type mismatch for histogram\n"
1855 <<typeid(*mobjh).name()<<endl;
1856 return;
1857}
1858
1859}
1860
1861/* methode */
1862void PAWExecutor::h_copy(vector<string>& tokens)
1863// Pour copier un object dans un object
1864// objects are Vector,Matrix,Histo,Histo2D
1865// h/copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]
1866// copy obj1Dfrom(i1->i2) into obj1Dto(ic1->)
1867// copy obj2Dfrom(i1,j1->i2,j2) into obj2Dto(ic1,jc1->)
1868// Attention: elements depuis i1 jusqu'a i2 COMPRIS
1869// Si obj1Dto n'existe pas, il est cree avec une taille i2-i1+1
1870// ou obj2Dto avec une taille i2-i1+1,j2-j1+1
1871// Le contenu de obj?Dfrom adresse est surecrit
1872// Le contenu de obj?Dfrom non adresse reste le meme
1873{
1874int_4 tks = tokens.size();
1875if(tks<2)
1876 {cout<<"Usage: h_copy namefrom nametocopy [i1[:i2]] [j1[:j2]] [ic1[:jc1]]"<<endl;
1877 return;}
1878
1879NamedObjMgr omg;
1880AnyDataObj* mobjv1 = omg.GetObj(tokens[0]);
1881if( mobjv1==NULL)
1882 {cout<<"PAWExecutor::h_copy Error: unknow object "<<tokens[0]<<endl;
1883 return;}
1884AnyDataObj* mobjv2 = omg.GetObj(tokens[1]);
1885
1886int_4 i1=0,i2=0, j1=0,j2=0, ic1=0,jc1=0, i,ii, j,jj, nx1,ny1, nx2,ny2;
1887
1888// Cas d'un Vector
1889if(typeid(*mobjv1) == typeid(Vector)) {
1890 Vector* v1 = dynamic_cast<Vector*>(mobjv1); nx1 = (int_4)v1->NElts();
1891 i2=nx1-1;
1892 if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
1893 if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1;
1894 if(i2<i1)
1895 {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2<<"] for NElts="
1896 <<nx1<<endl; return;}
1897 Vector* v2 = NULL;
1898 if(mobjv2==NULL)
1899 {v2 = new Vector(i2-i1+1); omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
1900 if(typeid(*mobjv2) != typeid(Vector))
1901 {cout<<"PAWExecutor::h_copy Error: type mismatch for Vector "
1902 <<typeid(*mobjv2).name()<<endl; return;}
1903 v2 = dynamic_cast<Vector*>(mobjv2); nx2 = (int_4)v2->NElts();
1904 if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d",&ic1);
1905 if(ic1<0) ic1=0;
1906 if(ic1>=nx2)
1907 {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<"] for NElts="
1908 <<nx2<<endl; return;}
1909 cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<") to "<<tokens[1]<<"("<<ic1<<"->...)"<<endl;
1910 for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++) (*v2)(ii) = (*v1)(i);
1911 return;
1912}
1913
1914// Cas d'un Histo
1915if(typeid(*mobjv1) == typeid(Histo)) {
1916 Histo* v1 = dynamic_cast<Histo*>(mobjv1); nx1 = (int_4)v1->NBins();
1917 i2=nx1-1;
1918 if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
1919 if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1;
1920 if(i2<i1)
1921 {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2<<"] for NBins="
1922 <<nx1<<endl; return;}
1923 Histo* v2 = NULL;
1924 if(mobjv2==NULL)
1925 {v2 = new Histo(0.,(r_8)(i2-i1+1),i2-i1+1);
1926 omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
1927 if(typeid(*mobjv2) != typeid(Histo))
1928 {cout<<"PAWExecutor::h_copy Error: type mismatch for Histo "
1929 <<typeid(*mobjv2).name()<<endl; return;}
1930 v2 = dynamic_cast<Histo*>(mobjv2); nx2 = (int_4)v2->NBins();
1931 if(v1->HasErrors() && !(v2->HasErrors())) v2->Errors();
1932 if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d",&ic1);
1933 if(ic1<0) ic1=0;
1934 if(ic1>=nx2)
1935 {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<"] for NBins="
1936 <<nx2<<endl; return;}
1937 cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<") to "<<tokens[1]<<"("<<ic1<<"->...)"<<endl;
1938 for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
1939 {(*v2)(ii) = (*v1)(i); if(v1->HasErrors()) v2->Error2(ii) = v1->Error2(i);}
1940 return;
1941}
1942
1943// Cas d'une Matrix
1944if(typeid(*mobjv1) == typeid(Matrix)) {
1945 Matrix* v1 = dynamic_cast<Matrix*>(mobjv1); nx1=v1->NRows(); ny1=v1->NCol();
1946 i2=nx1-1; j2=ny1-1;
1947 if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
1948 if(tks>3) if(tokens[3]!="!") sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
1949 if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1; if(j1<0) j1=0; if(j2>=ny1) j2=ny1-1;
1950 if(i2<i1 || j2<j1)
1951 {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2
1952 <<"] , ["<<j1<<":"<<j2<<"] for NRows,NCol="
1953 <<nx1<<" , "<<ny1<<endl; return;}
1954 Matrix* v2 = NULL;
1955 if(mobjv2==NULL)
1956 {v2 = new Matrix(i2-i1+1,j2-j1+1);
1957 omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
1958 if(typeid(*mobjv2) != typeid(Matrix))
1959 {cout<<"PAWExecutor::h_copy Error: type mismatch for Matrix "
1960 <<typeid(*mobjv2).name()<<endl; return;}
1961 v2 = dynamic_cast<Matrix*>(mobjv2); nx2=v2->NRows(); ny2=v2->NCol();
1962 if(tks>4) if(tokens[4]!="!") sscanf(tokens[4].c_str(),"%d:%d",&ic1,&jc1);
1963 if(ic1<0) ic1=0; if(jc1<0) jc1=0;
1964 if(ic1>=nx2 || jc1>=ny2)
1965 {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<","<<jc1
1966 <<"] for NRows,NCol="<<nx2<<","<<ny2<<endl; return;}
1967 cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<","<<j1<<":"<<j2
1968 <<") to "<<tokens[1]<<"("<<ic1<<","<<jc1<<"->...)"<<endl;
1969 for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
1970 for(j=j1,jj=jc1; j<=j2 && j<ny1 && jj<ny2; j++,jj++) (*v2)(ii,jj) = (*v1)(i,j);
1971 return;
1972}
1973
1974// Cas d'un Histo2D
1975if(typeid(*mobjv1) == typeid(Histo2D)) {
1976 Histo2D* v1 = dynamic_cast<Histo2D*>(mobjv1); nx1=v1->NBinX(); ny1=v1->NBinY();
1977 i2=nx1-1; j2=ny1-1;
1978 if(tks>2) if(tokens[2]!="!") sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
1979 if(tks>3) if(tokens[3]!="!") sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
1980 if(i1<0) i1=0; if(i2>=nx1) i2=nx1-1; if(j1<0) j1=0; if(j2>=ny1) j2=ny1-1;
1981 if(i2<i1 || j2<j1)
1982 {cout<<"PAWExecutor::h_copy Error: wrong range ["<<i1<<":"<<i2
1983 <<"] , ["<<j1<<":"<<j2<<"] for NBinX,NBinY="
1984 <<nx1<<" , "<<ny1<<endl; return;}
1985 Histo2D* v2 = NULL;
1986 if(mobjv2==NULL)
1987 {v2 = new Histo2D(0.,(r_8)(i2-i1+1),i2-i1+1,0.,(r_8)(j2-j1+1),j2-j1+1);
1988 omg.AddObj(v2,tokens[1]); mobjv2 = omg.GetObj(tokens[1]);}
1989 if(typeid(*mobjv2) != typeid(Histo2D))
1990 {cout<<"PAWExecutor::h_copy Error: type mismatch for Histo2D"
1991 <<typeid(*mobjv2).name()<<endl; return;}
1992 v2 = dynamic_cast<Histo2D*>(mobjv2); nx2=v2->NBinX(); ny2=v2->NBinY();
1993 if(v1->HasErrors() && !(v2->HasErrors())) v2->Errors();
1994 if(tks>4) sscanf(tokens[4].c_str(),"%d:%d",&ic1,&jc1);
1995 if(ic1<0) ic1=0; if(jc1<0) jc1=0;
1996 if(ic1>=nx2 || jc1>=ny2)
1997 {cout<<"PAWExecutor::h_copy Error: wrong pointer to copy ["<<ic1<<","<<jc1
1998 <<"] for NBinX,NBinY="<<nx2<<","<<ny2<<endl; return;}
1999 cout<<"copy "<<tokens[0]<<"("<<i1<<":"<<i2<<","<<j1<<":"<<j2
2000 <<") to "<<tokens[1]<<"("<<ic1<<","<<jc1<<"->...)"<<endl;
2001 for(i=i1,ii=ic1; i<=i2 && i<nx1 && ii<nx2; i++,ii++)
2002 for(j=j1,jj=jc1; j<=j2 && j<ny1 && jj<ny2; j++,jj++)
2003 {(*v2)(ii,jj) = (*v1)(i,j); if(v1->HasErrors()) v2->Error2(ii,jj) = v1->Error2(i,j);}
2004 return;
2005}
2006
2007// Cas non prevu
2008cout<<"PAWExecutor::h_copy Error: type mismatch for Vector/Matrix/Histo/Histo2D\n"
2009 <<typeid(*mobjv1).name()<<endl;
2010return;
2011}
2012
2013/* methode */
2014void PAWExecutor::h_set(string dum,vector<string>& tokens)
2015// Mise de valeurs/erreurs d'un histo 1D ou 2D a une valeur donnee
2016// dum = err : on change les valeurs des erreurs
2017// cont : on change les valeurs du contenu des bins
2018// tokens[0] : nom de l'histogramme
2019// tokens[1] : valeur de remplacement
2020// tokens[2-3] : i1:i2 j1:j2 intervalle des bins qui doivent etre remplace
2021// : v v1:v2 remplacement des bins ayant une valeur dans cet intervalle
2022// : e e1:e2 remplacement des bins ayant une erreur dans cet intervalle
2023{
2024int_4 tks = tokens.size();
2025if(tks<3)
2026 {cout<<"Usage: h/set/[err,cont] namehisto setvalue i1[:i2] [j1[:j2]]\n"
2027 <<" v v1:v2\n"
2028 <<" e e1:e2"
2029 <<endl; return;}
2030
2031// Decodage arguments
2032bool replerr = false; if(dum=="err") replerr = true;
2033r_8 setval = atof(tokens[1].c_str());
2034int_4 testcont=0; if(tokens[2]=="v") testcont=1; if(tokens[2]=="e") testcont=2;
2035int_4 i1=-1, i2=-1, j1=-1, j2=-1;
2036r_8 v1=0., v2=0.;
2037if(testcont==0) {
2038 sscanf(tokens[2].c_str(),"%d:%d",&i1,&i2);
2039 if(tks>3) sscanf(tokens[3].c_str(),"%d:%d",&j1,&j2);
2040} else {
2041 if(tks<=3)
2042 {cout<<"PAWExecutor::h_set Error: h/set/... v v1:v2, please give v1:v2"<<endl;
2043 return;}
2044 sscanf(tokens[3].c_str(),"%lf:%lf",&v1,&v2);
2045}
2046
2047if(replerr) {if(setval<0.) setval = 0.; else setval *= setval;}
2048if(testcont!=0 && v2<v1)
2049 {cout<<"PAWExecutor::h_set Error: bad value range="<<v1<<","<<v2<<endl; return;}
2050
2051// identification objet
2052NamedObjMgr omg;
2053AnyDataObj* mobj = omg.GetObj(tokens[0]);
2054if( mobj==NULL)
2055 {cout<<"PAWExecutor::h_set Error: unknow object "<<tokens[0]<<endl;
2056 return;}
2057
2058// Traitement
2059if(typeid(*mobj) == typeid(Histo)) { // Histo 1D
2060 Histo* h = dynamic_cast<Histo*>(mobj);
2061 if( (replerr || testcont==2) && !(h->HasErrors()) )
2062 {cout<<"PAWExecutor::h_set Error: histogram has no errors"<<endl; return;}
2063 if(testcont!=0) {i1=0; i2=h->NBins()-1;}
2064 if(i1<0 || i1>=h->NBins())
2065 {cout<<"PAWExecutor::h_set Error: bad bin range i1="<<i1<<endl; return;}
2066 if(i2<i1) i2=i1; if(i2>=h->NBins()) i2=h->NBins()-1;
2067 for(int_4 i=i1;i<=i2;i++) {
2068 bool change = true;
2069 if(testcont==1) {if((*h)(i)<v1 || (*h)(i)>v2) change = false;}
2070 else if(testcont==2) {if(h->Error(i)<v1 || h->Error(i)>v2) change = false;}
2071 if(!change) continue;
2072 if(replerr) h->Error2(i) = setval; else (*h)(i) = setval;
2073 }
2074
2075} else if(typeid(*mobj) == typeid(Histo2D)) { // Histo 2D
2076 Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
2077 if( (replerr || testcont==2) && !(h2->HasErrors()) )
2078 {cout<<"PAWExecutor::h_set Error: histogram has no errors"<<endl; return;}
2079 if(testcont!=0) {i1=0; i2=h2->NBinX()-1;j1=0; j2=h2->NBinY()-1;}
2080 if(i1<0 || i1>=h2->NBinX() || j1<0 || j1>=h2->NBinY())
2081 {cout<<"PAWExecutor::h_set Error: bad bin range i1,j1="<<i1<<","<<j1<<endl; return;}
2082 if(i2<i1) i2=i1; if(i2>=h2->NBinX()) i2=h2->NBinX()-1;
2083 if(j2<j1) j2=j1; if(j2>=h2->NBinY()) j2=h2->NBinY()-1;
2084 for(int_4 i=i1;i<=i2;i++) for(int_4 j=j1;j<=j2;j++) {
2085 bool change = true;
2086 if(testcont==1) {if((*h2)(i,j)<v1 || (*h2)(i,j)>v2) change = false;}
2087 else if(testcont==2) {if(h2->Error(i,j)<v1 || h2->Error(i,j)>v2) change = false;}
2088 if(!change) continue;
2089 if(replerr) h2->Error2(i,j) = setval; else (*h2)(i,j) = setval;
2090 }
2091
2092} else {
2093 cout<<"PAWExecutor::h_set Error: type mismatch for Histo/Histo2D\n"
2094 <<typeid(*mobj).name()<<endl;
2095 return;
2096}
2097
2098}
2099
2100/* methode */
2101void PAWExecutor::h_err(vector<string>& tokens)
2102// Mise des erreurs d'un histo 1D ou 2D a une valeur
2103// donnee par une fonction de la valeur du bin
2104// tokens[0] : nom de l'histogramme
2105// tokens[1] : expression de la fonction
2106{
2107if(tokens.size()<2)
2108 {cout<<"Usage: h/err namehisto expr_func"<<endl; return;}
2109
2110// identification objet
2111NamedObjMgr omg;
2112AnyDataObj* mobj = omg.GetObj(tokens[0]);
2113if( mobj==NULL)
2114 {cout<<"PAWExecutor::h_err Error: unknow object "<<tokens[0]<<endl;
2115 return;}
2116
2117// Fonction
2118FILE *fip = NULL;
2119string expfunc = ""; {for(int_4 i=1;i<(int)tokens.size();i++) expfunc += tokens[i];}
2120string fname = "func1or2_pia_dl.c";
2121string func = "func1or2_pia_dl_func";
2122if((fip = fopen(fname.c_str(), "w")) == NULL)
2123 {cout<<"PAWExecutor::h_err Error: open function file "<<fname<<endl;
2124 return;}
2125fprintf(fip,"#include <math.h>\n");
2126fprintf(fip,"double %s(double x) \n{\n",func.c_str());
2127fprintf(fip,"return(%s); \n}\n", expfunc.c_str());
2128fclose(fip);
2129DlFunctionOfX f = (DlFunctionOfX) omg.GetServiceObj()->LinkFunctionFromFile(fname,func);
2130if(!f)
2131 {cout<<"PAWExecutor::h_err Error: bad function "<<func<<","<<fname<<endl;
2132 return;}
2133
2134// Traitement
2135if(typeid(*mobj) == typeid(Histo)) { // Histo 1D
2136 Histo* h = dynamic_cast<Histo*>(mobj);
2137 if(!(h->HasErrors())) h->Errors();
2138 for(int_4 i=0;i<h->NBins();i++) {
2139 r_8 x = (*h)(i);
2140 r_8 e = f(x);
2141 h->SetErr2(i,e*e);
2142 }
2143
2144} else if(typeid(*mobj) == typeid(Histo2D)) { // Histo 2D
2145 Histo2D* h2 = dynamic_cast<Histo2D*>(mobj);
2146 if(!(h2->HasErrors())) h2->Errors();
2147 for(int_4 i=0;i<h2->NBinX();i++) for(int_4 j=0;j<h2->NBinY();j++) {
2148 r_8 x = (*h2)(i,j);
2149 r_8 e = f(x);
2150 h2->Error2(i,j) = e*e;
2151 }
2152
2153} else {
2154 cout<<"PAWExecutor::h_err Error: type mismatch for Histo/Histo2D\n"
2155 <<typeid(*mobj).name()<<endl;
2156 return;
2157}
2158
2159}
Note: See TracBrowser for help on using the repository browser.