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

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

n/plot keepbin option do not impose anymore the "same" option cmv 29/04/2005

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