ROOT   6.08/07 Reference Guide
Roo2DMomentMorphFunction.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * author: Max Baak (mbaak@cern.ch) *
4  *****************************************************************************/
5
6 /** \class Roo2DMomentMorphFunction
7  \ingroup Roofit
8
9 2-dimensional morph function between a list of function-numbers as a function of two input parameters (m1 and m2).
10 The matrix mrefpoints assigns a function value to each m1,m2 coordinate.
11 Morphing can be set to be linear, non-linear or a mixture of the two.
12
13 Usage example:
14 ~~~ {.cpp}
15  TMatrixD foo(9,3);
16
17  foo(0,0) = 0; // coordinate of variable m1
18  foo(0,1) = 0; // coordinate of variable m2
19  foo(0,2) = 0**2; // function value at (m1,m2) = 0,0
20  foo(1,0) = 100;
21  foo(1,1) = 0;
22  foo(1,2) = 1**2;
23  foo(2,0) = 200;
24  foo(2,1) = 0;
25  foo(2,2) = 2**2;
26
27  foo(3,0) = 0;
28  foo(3,1) = 150;
29  foo(3,2) = 4**2;
30  foo(4,0) = 100;
31  foo(4,1) = 150;
32  foo(4,2) = 5**2;
33  foo(5,0) = 200;
34  foo(5,1) = 150;
35  foo(5,2) = 8**2;
36
37  foo(6,0) = 0;
38  foo(6,1) = 300;
39  foo(6,2) = 8**2;
40  foo(7,0) = 100;
41  foo(7,1) = 300;
42  foo(7,2) = 8.5**2;
43  foo(8,0) = 200;
44  foo(8,1) = 300;
45  foo(8,2) = 9**2;
46
47  // need to provide at least 4 coordinates to Roo2DMomentMorphFunction for 2d extrapolation
48
49  foo.Print();
50
51  RooRealVar m1("m1","m1",50,0,200);
52  RooRealVar m2("m2","m2",50,0,300);
53  Roo2DMomentMorphFunction bar("bar","bar", m1, m2, foo );
54
55  bar.Print();
56 ~~~
57 */
58
59
60 #include "Riostream.h"
61
63 #include "RooAbsReal.h"
64 #include "RooAbsCategory.h"
65 #include <math.h>
66 #include "TMath.h"
67 #include "TTree.h"
68
69
70 using namespace std;
71
73
74
75 ////////////////////////////////////////////////////////////////////////////////
76 /// Constructor.
77 /// Cross-check that we have enough reference points.
78 /// \param[in] name
79 /// \param[in] title
80 /// \param[in] _m1
81 /// \param[in] _m2
82 /// \param[in] mrefpoints
83 /// \param[in] setting
84 /// \param[in] verbose
85
86 Roo2DMomentMorphFunction::Roo2DMomentMorphFunction(const char *name, const char *title,
87  RooAbsReal& _m1, RooAbsReal& _m2,
88  const TMatrixD& mrefpoints, const Setting& setting, const Bool_t& verbose ) :
89  RooAbsReal(name,title),
90  m1("m1","m1",this,_m1),
91  m2("m2","m2",this,_m2),
92  _setting(setting),
93  _verbose(verbose),
94  _npoints( mrefpoints.GetNrows() ),
95  _mref(mrefpoints)
96 {
97  if ( mrefpoints.GetNrows()<4 ) {
98  cerr << "Roo2DMomentMorphFunction::constructor(" << GetName() << ") ERROR: less than four reference points provided." << endl ;
99  assert(0);
100  }
101  // cross-check that we have enough reference points
102  if ( mrefpoints.GetNcols()!=3 ) {
103  cerr << "RooPolyMorph2D::constructor(" << GetName() << ") ERROR: no reference values provided." << endl ;
104  assert(0);
105  }
106
107  // recast matrix into more useful form
108  _frac.ResizeTo( _npoints );
109
110  initialize();
111 }
112
113
114 ////////////////////////////////////////////////////////////////////////////////
115 /// Constructor.
116 /// Cross-check that we have enough reference points.
117 /// \param[in] name
118 /// \param[in] title
119 /// \param[in] _m1
120 /// \param[in] _m2
121 /// \param[in] nrows
122 /// \param[in] dm1arr
123 /// \param[in] dm2arr
124 /// \param[in] dvalarr
125 /// \param[in] setting
126 /// \param[in] verbose
127
129  RooAbsReal& _m1, RooAbsReal& _m2,
130  const Int_t& nrows, const Double_t* dm1arr, const Double_t* dm2arr, const Double_t* dvalarr,
131  const Setting& setting, const Bool_t& verbose ) :
132  RooAbsReal(name,title),
133  m1("m1","m1",this,_m1),
134  m2("m2","m2",this,_m2),
135  _setting(setting),
136  _verbose( verbose ),
137  _npoints( nrows )
138 {
139  if ( nrows<4 ) {
140  cerr << "Roo2DMomentMorphFunction::constructor(" << GetName() << ") ERROR: less than four reference points provided." << endl ;
141  assert(0);
142  }
143
144  // recast tree into more useful form
145  _mref.ResizeTo( _npoints,3 );
147
148  for (int i=0; i<_npoints; ++i) {
149  _mref(i,0) = dm1arr[i] ;
150  _mref(i,1) = dm2arr[i] ;
151  _mref(i,2) = dvalarr[i] ;
152  }
153
154  initialize();
155 }
156
157
158 ////////////////////////////////////////////////////////////////////////////////
159 /// Copy constructor.
160
162  RooAbsReal(other,name),
163  m1("m1",this,other.m1),
164  m2("m2",this,other.m2),
165  _setting(other._setting),
166  _verbose(other._verbose),
167  _npoints(other._npoints),
168  _mref(other._mref),
169  _frac(other._frac)
170 {
171  initialize();
172 }
173
174
175 ////////////////////////////////////////////////////////////////////////////////
176 /// Destructor.
177
179 {
180 }
181
182
183 ////////////////////////////////////////////////////////////////////////////////
184
185 void
187 {
188  Double_t xmin(1e300), xmax(-1e300), ymin(1e300), ymax(-1e300);
189
190  // transformation matrix for non-linear extrapolation, needed in evaluate()
191  for (Int_t k=0; k<_npoints; ++k) {
192  if (_mref(k,0)<xmin) { xmin=_mref(k,0); _ixmin=k; }
193  if (_mref(k,0)>xmax) { xmax=_mref(k,0); _ixmax=k; }
194  if (_mref(k,1)<ymin) { ymin=_mref(k,1); _iymin=k; }
195  if (_mref(k,1)>ymax) { ymax=_mref(k,1); _iymax=k; }
196  }
197
198  // resize
199  _MSqr.ResizeTo(4,4);
200  _squareVec.ResizeTo(4,2);
201 }
202
203
204 ////////////////////////////////////////////////////////////////////////////////
205
206 Double_t
208 {
209  if (_verbose) { cout << "Now in evaluate." << endl; }
210  if (_verbose) { cout << "x = " << m1 << " ; y = " << m2 << endl; }
211
212  calculateFractions(_verbose); // verbose turned off
213
214  Double_t ret = 0.0;
215  for (Int_t i=0; i<4; ++i) { ret += ( _mref(_squareIdx[i],2) * _frac[_squareIdx[i]] ) ; }
216
217  if (_verbose) { cout << "End of evaluate : " << ret << endl; }
218
219  //return (ret>0 ? ret : 0) ;
220  return ret;
221 }
222
223
224 ////////////////////////////////////////////////////////////////////////////////
225
226 void
228 {
229  double sumfrac(0.);
230
232  // reset all fractions
233  for (Int_t i=0; i<_npoints; ++i) { _frac[i]=0.0; }
234
235  // this sets _squareVec and _squareIdx quantities and MSqr
236  (void) findSquare(m1,m2);
237
238  std::vector<double> deltavec(4,1.0);
239  deltavec[1] = m1-_squareVec(0,0) ;
240  deltavec[2] = m2-_squareVec(0,1) ;
241  deltavec[3] = deltavec[1]*deltavec[2] ;
242
243  for (Int_t i=0; i<4; ++i) {
244  double ffrac=0.;
245  for (Int_t j=0; j<4; ++j) { ffrac += _MSqr(j,i) * deltavec[j]; }
246
247  // set fractions
248  _frac[_squareIdx[i]] = ffrac;
249  if (ffrac>0) sumfrac += ffrac;
250
251  if (verbose) {
252  cout << _squareIdx[i] << " " << ffrac << " " << _squareVec(i,0) << " " << _squareVec(i,1) << endl;
253  }
254  }
255  }
256
257  if (_setting == LinearPosFractions) {
258  for (Int_t i=0; i<4; ++i) {
259  if (_frac[_squareIdx[i]]<0) { _frac[_squareIdx[i]] = 0; }
260  _frac[_squareIdx[i]] *= (1.0/sumfrac) ;
261  }
262  }
263 }
264
265
266 ////////////////////////////////////////////////////////////////////////////////
267
268 Bool_t
269 Roo2DMomentMorphFunction::findSquare(const double& x, const double& y) const
270 {
271  bool squareFound(false);
272
273  std::vector< std::pair<int,double> > idvec;
274
275  Double_t xspacing = (_mref(_ixmax,0)-_mref(_ixmin,0)) / TMath::Sqrt(_npoints) ;
276  Double_t yspacing = (_mref(_iymax,1)-_mref(_iymin,1)) / TMath::Sqrt(_npoints) ;
277
278  Double_t dx(0), dy(0), delta(0);
279  for (Int_t k=0; k<_npoints; ++k) {
280  dx = (x-_mref(k,0))/xspacing ;
281  dy = (y-_mref(k,1))/yspacing ;
282  delta = TMath::Sqrt(dx*dx+dy*dy) ;
283  idvec.push_back( std::pair<int,double>(k,delta) );
284  }
285
286  // order points closest to (x,y)
287  sort(idvec.begin(),idvec.end(),SorterL2H());
288  std::vector< std::pair<int,double> >::iterator itr = idvec.begin();
289  std::vector<int> cidx;
290  for(; itr!=idvec.end(); ++itr) {
291  cidx.push_back(itr->first);
292  }
293
294  // 1. point falls outside available ref points: pick three square points.
295  // fall-back option
296  _squareVec(0,0) = _mref(cidx[0],0);
297  _squareVec(0,1) = _mref(cidx[0],1);
298  _squareIdx[0] = cidx[0];
299  _squareVec(1,0) = _mref(cidx[1],0);
300  _squareVec(1,1) = _mref(cidx[1],1);
301  _squareIdx[1] = cidx[1];
302  _squareVec(2,0) = _mref(cidx[2],0);
303  _squareVec(2,1) = _mref(cidx[2],1);
304  _squareIdx[2] = cidx[2];
305  _squareVec(3,0) = _mref(cidx[3],0);
306  _squareVec(3,1) = _mref(cidx[3],1);
307  _squareIdx[3] = cidx[3];
308
309  // 2. try to find square enclosing square
310  if ( x>_mref(_ixmin,0) &&
311  x<_mref(_ixmax,0) &&
312  y>_mref(_iymin,1) &&
313  y<_mref(_iymax,1) )
314  {
315  for (unsigned int i=0; i<cidx.size() && !squareFound; ++i)
316  for (unsigned int j=i+1; j<cidx.size() && !squareFound; ++j)
317  for (unsigned int k=j+1; k<cidx.size() && !squareFound; ++k)
318  for (unsigned int l=k+1; l<cidx.size() && !squareFound; ++l) {
319  if ( isAcceptableSquare(_mref(cidx[i],0),_mref(cidx[i],1),_mref(cidx[j],0),_mref(cidx[j],1),_mref(cidx[k],0),_mref(cidx[k],1),_mref(cidx[l],0),_mref(cidx[l],1)) ) {
320  if ( pointInSquare(x,y,_mref(cidx[i],0),_mref(cidx[i],1),_mref(cidx[j],0),_mref(cidx[j],1),_mref(cidx[k],0),_mref(cidx[k],1),_mref(cidx[l],0),_mref(cidx[l],1)) ) {
321  _squareVec(0,0) = _mref(cidx[i],0);
322  _squareVec(0,1) = _mref(cidx[i],1);
323  _squareIdx[0] = cidx[i];
324  _squareVec(1,0) = _mref(cidx[j],0);
325  _squareVec(1,1) = _mref(cidx[j],1);
326  _squareIdx[1] = cidx[j];
327  _squareVec(2,0) = _mref(cidx[k],0);
328  _squareVec(2,1) = _mref(cidx[k],1);
329  _squareIdx[2] = cidx[k];
330  _squareVec(3,0) = _mref(cidx[l],0);
331  _squareVec(3,1) = _mref(cidx[l],1);
332  _squareIdx[3] = cidx[l];
333  squareFound=true;
334  }
335  }
336  }
337  }
338
339  // construct transformation matrix for linear extrapolation
340  TMatrixD M(4,4);
341  for (Int_t k=0; k<4; ++k) {
342  dx = _squareVec(k,0) - _squareVec(0,0) ;
343  dy = _squareVec(k,1) - _squareVec(0,1) ;
344  M(k,0) = 1.0 ;
345  M(k,1) = dx ;
346  M(k,2) = dy ;
347  M(k,3) = dx*dy ;
348  }
349
350  _MSqr = M.Invert();
351
352  return squareFound;
353 }
354
355
356 ////////////////////////////////////////////////////////////////////////////////
357 /// Returns true if p1 and p2 are on same side of line b-a.
358 /// \param[in] p1x
359 /// \param[in] p1y
360 /// \param[in] p2x
361 /// \param[in] p2y
362 /// \param[in] ax
363 /// \param[in] ay
364 /// \param[in] bx
365 /// \param[in] by
366
367 Bool_t Roo2DMomentMorphFunction::onSameSide(const double& p1x, const double& p1y, const double& p2x, const double& p2y, const double& ax, const double& ay, const double& bx, const double& by) const
368 {
369  Double_t cp1 = myCrossProduct(bx-ax, by-ay, p1x-ax, p1y-ay);
370  Double_t cp2 = myCrossProduct(bx-ax, by-ay, p2x-ax, p2y-ay);
371  if (cp1*cp2 >= 0) return true;
372  else return false;
373 }
374
375
376 ////////////////////////////////////////////////////////////////////////////////
377
378 Bool_t
379 Roo2DMomentMorphFunction::pointInSquare(const double& px, const double& py, const double& ax, const double& ay, const double& bx, const double& by, const double& cx, const double& cy, const double& dx, const double& dy) const
380 {
381  bool insquare(false);
382
383  int ntri(0);
384  if (ntri<2) ntri += static_cast<int>( pointInTriangle(px,py,ax,ay,bx,by,cx,cy) );
385  if (ntri<2) ntri += static_cast<int>( pointInTriangle(px,py,ax,ay,bx,by,dx,dy) );
386  if (ntri<2) ntri += static_cast<int>( pointInTriangle(px,py,ax,ay,cx,cy,dx,dy) );
387  if (ntri<2) ntri += static_cast<int>( pointInTriangle(px,py,bx,by,cx,cy,dx,dy) );
388
389  if (ntri>=2) insquare=true;
390  else insquare=false;
391
392  return insquare;
393 }
394
395
396 ////////////////////////////////////////////////////////////////////////////////
397
398 Bool_t
399 Roo2DMomentMorphFunction::pointInTriangle(const double& px, const double& py, const double& ax, const double& ay, const double& bx, const double& by, const double& cx, const double& cy) const
400 {
401  if (onSameSide(px,py,ax,ay,bx,by,cx,cy) && onSameSide(px,py,bx,by,ax,ay,cx,cy) && onSameSide(px,py,cx,cy,ax,ay,bx,by)) return true;
402  else return false;
403 }
404
405
406 ////////////////////////////////////////////////////////////////////////////////
407
408 Double_t
409 Roo2DMomentMorphFunction::myCrossProduct(const double& ax, const double& ay, const double& bx, const double& by) const
410 {
411  return ( ax*by - bx*ay );
412 }
413
414
415 ////////////////////////////////////////////////////////////////////////////////
416 /// Reject kinked shapes.
417 /// \param[in] ax
418 /// \param[in] ay
419 /// \param[in] bx
420 /// \param[in] by
421 /// \param[in] cx
422 /// \param[in] cy
423 /// \param[in] dx
424 /// \param[in] dy
425
426 Bool_t
427 Roo2DMomentMorphFunction::isAcceptableSquare(const double& ax, const double& ay, const double& bx, const double& by, const double& cx, const double& cy, const double& dx, const double& dy) const
428 {
429  if ( pointInTriangle(dx,dy,ax,ax,bx,by,cx,cy) ||
430  pointInTriangle(cx,cy,ax,ay,bx,by,dx,dy) ||
431  pointInTriangle(bx,by,ax,ay,cx,cy,dx,dy) ||
432  pointInTriangle(ax,ay,bx,by,cx,cy,dx,dy) ) return false;
433  else return true;
434 }
435
436
437 void
439 {
440  for( Int_t i=0; i<_npoints; i++ ){
441  cout << this << " " << i << " " << _mref(i,0) << " " << _mref(i,1) << " " << _mref(i,2) << endl;
442  }
443 }
444
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:140
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
float xmin
Definition: THbookFile.cxx:93
virtual ~Roo2DMomentMorphFunction()
Destructor.
Bool_t pointInTriangle(const double &px, const double &py, const double &ax, const double &ay, const double &bx, const double &by, const double &cx, const double &cy) const
float ymin
Definition: THbookFile.cxx:93
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
Double_t x[n]
Definition: legend1.C:17
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1396
Double_t myCrossProduct(const double &ax, const double &ay, const double &bx, const double &by) const
2-dimensional morph function between a list of function-numbers as a function of two input parameters...
float ymax
Definition: THbookFile.cxx:93
Bool_t findSquare(const double &x, const double &y) const
Bool_t isAcceptableSquare(const double &ax, const double &ay, const double &bx, const double &by, const double &cx, const double &cy, const double &dx, const double &dy) const
Reject kinked shapes.
bool verbose
TLine * l
Definition: textangle.C:4
float xmax
Definition: THbookFile.cxx:93
void calculateFractions(Bool_t verbose=kTRUE) const
#define ClassImp(name)
Definition: Rtypes.h:279
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Bool_t onSameSide(const double &p1x, const double &p1y, const double &p2x, const double &p2y, const double &ax, const double &ay, const double &bx, const double &by) const
Returns true if p1 and p2 are on same side of line b-a.
Double_t y[n]
Definition: legend1.C:17
typedef void((*Func_t)())
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
void initialize(typename Architecture_t::Matrix_t &A, EInitialization m)
Definition: Functions.h:257
Bool_t pointInSquare(const double &px, const double &py, const double &ax, const double &ay, const double &bx, const double &by, const double &cx, const double &cy, const double &dx, const double &dy) const
char name[80]
Definition: TGX11.cxx:109