Spline.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 */
10 //____________________________________________________________________________
11 
12 #include <cassert>
13 #include <iomanip>
14 #include <cfloat>
15 
16 #include "libxml/parser.h"
17 #include "libxml/xmlmemory.h"
18 
19 #include <TFile.h>
20 #include <TNtupleD.h>
21 #include <TTree.h>
22 #include <TSQLServer.h>
23 #include <TGraph.h>
24 #include <TMath.h>
25 #include <TH2F.h>
26 #include <TROOT.h>
27 
32 
33 using std::endl;
34 using std::setw;
35 using std::setprecision;
36 using std::setfill;
37 using std::ios;
38 
39 using namespace genie;
40 
42 
43 //___________________________________________________________________________
45 {
46  ostream & operator << (ostream & stream, const Spline & spl)
47  {
48  spl.Print(stream);
49  return stream;
50  }
51 }
52 //___________________________________________________________________________
54 {
55  this->InitSpline();
56 }
57 //___________________________________________________________________________
58 Spline::Spline(string filename, string xtag, string ytag, bool is_xml) :
59 TObject()
60 {
61  string fmt = (is_xml) ? "XML" : "ASCII";
62 
63  LOG("Spline", pDEBUG)
64  << "Constructing spline from data in " << fmt << " file: " << filename;
65 
66  this->InitSpline();
67 
68  if(is_xml)
69  this->LoadFromXmlFile(filename, xtag, ytag);
70  else
71  this->LoadFromAsciiFile(filename);
72 }
73 //___________________________________________________________________________
74 Spline::Spline(TNtupleD * ntuple, string var, string cut) :
75 TObject()
76 {
77  LOG("Spline", pDEBUG) << "Constructing spline from data in a TNtuple";
78 
79  this->InitSpline();
80  this->LoadFromNtuple(ntuple, var, cut);
81 }
82 //___________________________________________________________________________
83 Spline::Spline(TTree * tree, string var, string cut) :
84 TObject()
85 {
86  LOG("Spline", pDEBUG) << "Constructing spline from data in a TTree";
87 
88  this->InitSpline();
89  this->LoadFromTree(tree, var, cut);
90 }
91 //___________________________________________________________________________
92 Spline::Spline(TSQLServer * db, string query) :
93 TObject()
94 {
95  LOG("Spline", pDEBUG) << "Constructing spline from data in a MySQL server";
96 
97  this->InitSpline();
98  this->LoadFromDBase(db, query);
99 }
100 //___________________________________________________________________________
101 Spline::Spline(int nentries, double x[], double y[]) :
102 TObject()
103 {
104  LOG("Spline", pDEBUG)
105  << "Constructing spline from the arrays passed to the ctor";
106 
107  this->InitSpline();
108  this->BuildSpline(nentries, x, y);
109 }
110 //___________________________________________________________________________
111 Spline::Spline(int nentries, float x[], float y[]) :
112 TObject()
113 {
114  LOG("Spline", pDEBUG)
115  << "Constructing spline from the arrays passed to the ctor";
116 
117  double * dblx = new double[nentries];
118  double * dbly = new double[nentries];
119 
120  for(int i = 0; i < nentries; i++) {
121  dblx[i] = double( x[i] );
122  dbly[i] = double( y[i] );
123  }
124 
125  this->InitSpline();
126  this->BuildSpline(nentries, dblx, dbly);
127 
128  delete [] dblx;
129  delete [] dbly;
130 }
131 //___________________________________________________________________________
133  TObject(), fInterpolator(0)
134 {
135  LOG("Spline", pDEBUG) << "Spline copy constructor";
136 
137  this->LoadFromTSpline3( *spline.GetAsTSpline(), spline.NKnots() );
138 }
139 //___________________________________________________________________________
140 Spline::Spline(const TSpline3 & spline, int nknots) :
141  TObject(), fInterpolator(0)
142 {
143  LOG("Spline", pDEBUG)
144  << "Constructing spline from the input TSpline3 object";
145 
146  this->LoadFromTSpline3( spline, nknots );
147 }
148 //___________________________________________________________________________
150 {
151  if(fInterpolator) delete fInterpolator;
152 }
153 //___________________________________________________________________________
154 bool Spline::LoadFromXmlFile(string filename, string xtag, string ytag)
155 {
156  LOG("Spline", pDEBUG) << "Retrieving data from file: " << filename;
157 
158  xmlDocPtr xml_doc = xmlParseFile(filename.c_str());
159 
160  if(xml_doc==NULL) {
161  LOG("Spline", pERROR)
162  << "XML file could not be parsed! [filename: " << filename << "]";
163  return false;
164  }
165 
166  xmlNodePtr xmlCur = xmlDocGetRootElement(xml_doc);
167 
168  if(xmlCur==NULL) {
169  LOG("Spline", pERROR)
170  << "XML doc. has null root element! [filename: " << filename << "]";
171  return false;
172  }
173  if( xmlStrcmp(xmlCur->name, (const xmlChar *) "spline") ) {
174  LOG("Spline", pERROR)
175  << "XML doc. has invalid root element! [filename: " << filename << "]";
176  return false;
177  }
178 
179  string name = utils::str::TrimSpaces(
180  utils::xml::GetAttribute(xmlCur, "name"));
181  string snkn = utils::str::TrimSpaces(
182  utils::xml::GetAttribute(xmlCur, "nknots"));
183  int nknots = atoi( snkn.c_str() );
184 
185  LOG("Spline", pINFO)
186  << "Parsing XML spline: " << name << ", nknots = " << nknots;
187 
188  int iknot = 0;
189  double * vx = new double[nknots];
190  double * vy = new double[nknots];
191 
192  xmlNodePtr xmlSplChild = xmlCur->xmlChildrenNode; // <knots>'s
193 
194  // loop over all xml tree nodes that are children of the <spline> node
195  while (xmlSplChild != NULL) {
196  LOG("Spline", pDEBUG)
197  << "Got <spline> children node: " << xmlSplChild->name;
198 
199  // enter everytime you find a <knot> tag
200  if( (!xmlStrcmp(xmlSplChild->name, (const xmlChar *) "knot")) ) {
201 
202  xmlNodePtr xmlKnotChild = xmlSplChild->xmlChildrenNode;
203 
204  // loop over all xml tree nodes that are children of this <knot>
205  while (xmlKnotChild != NULL) {
206  LOG("Spline", pDEBUG)
207  << "Got <knot> children node: " << xmlKnotChild->name;
208 
209  // enter everytime you find a <E> or a <xsec> tag
210  const xmlChar * tag = xmlKnotChild->name;
211  bool is_xtag = ! xmlStrcmp(tag,(const xmlChar *) xtag.c_str());
212  bool is_ytag = ! xmlStrcmp(tag,(const xmlChar *) ytag.c_str());
213  if (is_xtag || is_ytag) {
214  xmlNodePtr xmlValTagChild = xmlKnotChild->xmlChildrenNode;
215  string val = utils::xml::TrimSpaces(
216  xmlNodeListGetString(xml_doc, xmlValTagChild, 1));
217 
218  if (is_xtag) vx[iknot] = atof(val.c_str());
219  if (is_ytag) vy[iknot] = atof(val.c_str());
220 
221  xmlFree(xmlValTagChild);
222  LOG("Spline", pDEBUG) << "tag: " << tag << ", value: " << val;
223  }//if current tag is <E>,<xsec>
224 
225  xmlKnotChild = xmlKnotChild->next;
226  }//[end of] loop over tags within <knot>...</knot> tags
227  xmlFree(xmlKnotChild);
228  iknot++;
229  } // if <knot>
230  xmlSplChild = xmlSplChild->next;
231  } //[end of] loop over tags within <spline>...</spline> tags
232  xmlFree(xmlSplChild);
233 
234  this->BuildSpline(nknots, vx, vy);
235  delete [] vx;
236  delete [] vy;
237 
238  return true;
239 }
240 //___________________________________________________________________________
242 {
243  LOG("Spline", pDEBUG) << "Retrieving data from file: " << filename;
244 
245  TNtupleD nt("ntuple","","x:y");
246  nt.ReadFile(filename.c_str());
247 
248  this->LoadFromNtuple(&nt, "x:y");
249  return true;
250 }
251 //___________________________________________________________________________
252 bool Spline::LoadFromNtuple(TNtupleD * nt, string var, string cut)
253 {
254  if(!nt) return false;
255 
256  TTree * tree = dynamic_cast<TTree *> (nt);
257  return this->LoadFromTree(tree,var,cut);
258 }
259 //___________________________________________________________________________
260 bool Spline::LoadFromTree(TTree * tree, string var, string cut)
261 {
262  LOG("Spline", pDEBUG) << "Retrieving data from tree: " << tree->GetName();
263 
264  if(!cut.size()) tree->Draw(var.c_str(), "", "GOFF");
265  else tree->Draw(var.c_str(), cut.c_str(), "GOFF");
266 
267  TH2F * hst = (TH2F*)gROOT->FindObject("htemp");
268  if(hst) { hst->SetDirectory(0); delete hst; }
269 
270  // Now, take into account that the data retrieved from the ntuple would
271  // not be sorted in x and the resulting spline will be bogus...
272  // Sort the x array, use the x-sorting to re-arrange the y array and only
273  // then build the spline..
274 
275  int n = tree->GetSelectedRows();
276 
277  int * idx = new int[n];
278  double * x = new double[n];
279  double * y = new double[n];
280 
281  TMath::Sort(n,tree->GetV1(),idx,false);
282 
283  for(int i=0; i<n; i++) {
284  int ii = idx[i];
285  x[i] = (tree->GetV1())[ii];
286  y[i] = (tree->GetV2())[ii];
287  }
288 
289  this->BuildSpline(n,x,y);
290 
291  delete [] idx;
292  delete [] x;
293  delete [] y;
294 
295  return true;
296 }
297 //___________________________________________________________________________
298 bool Spline::LoadFromDBase(TSQLServer * /*db*/, string /*query*/)
299 {
300  LOG("Spline", pDEBUG) << "Retrieving data from data-base: ";
301  return false;
302 }
303 //___________________________________________________________________________
304 bool Spline::LoadFromTSpline3(const TSpline3 & spline, int nknots)
305 {
306  int nentries = nknots;
307 
308  double * x = new double[nentries];
309  double * y = new double[nentries];
310  double cx = 0, cy = 0;
311 
312  for(int i = 0; i < nentries; i++) {
313  spline.GetKnot(i, cx, cy);
314  x[i] = cx;
315  y[i] = cy;
316  }
317  this->BuildSpline(nentries, x, y);
318 
319  delete [] x;
320  delete [] y;
321 
322  return true;
323 }
324 //___________________________________________________________________________
325 void Spline::GetKnot(int iknot, double & x, double & y) const
326 {
327  if(!fInterpolator) {
328  LOG("Spline", pWARN) << "Spline has not been built yet!";
329  return;
330  }
331  fInterpolator->GetKnot(iknot,x,y);
332 }
333 //___________________________________________________________________________
334 double Spline::GetKnotX(int iknot) const
335 {
336  if(!fInterpolator) {
337  LOG("Spline", pWARN) << "Spline has not been built yet!";
338  return 0;
339  }
340  double x,y;
341  fInterpolator->GetKnot(iknot,x,y);
342  return x;
343 }
344 //___________________________________________________________________________
345 double Spline::GetKnotY(int iknot) const
346 {
347  if(!fInterpolator) {
348  LOG("Spline", pWARN) << "Spline has not been built yet!";
349  return 0;
350  }
351  double x,y;
352  fInterpolator->GetKnot(iknot,x,y);
353  return y;
354 }
355 //___________________________________________________________________________
356 bool Spline::IsWithinValidRange(double x) const
357 {
358  bool is_in_range = (fXMin <= x && x <= fXMax);
359  return is_in_range;
360 }
361 //___________________________________________________________________________
362 double Spline::Evaluate(double x) const
363 {
364  LOG("Spline", pDEBUG) << "Evaluating spline at point x = " << x;
365  assert(!TMath::IsNaN(x));
366 
367  double y = 0;
368  if( this->IsWithinValidRange(x) ) {
369 
370  // we can interpolate within the range of spline knots - be careful with
371  // strange cubic spline behaviour when close to knots with y=0
372  bool is0p = this->ClosestKnotValueIsZero(x, "+");
373  bool is0n = this->ClosestKnotValueIsZero(x, "-");
374 
375  if(!is0p && !is0n) {
376  // both knots (on the left and right are non-zero) - just interpolate
377  LOG("Spline", pDEBUG) << "Point is between non-zero knots";
378  y = fInterpolator->Eval(x);
379  } else {
380  // at least one of the neighboring knots has y=0
381  if(is0p && is0n) {
382  // both neighboring knots have y=0
383  LOG("Spline", pDEBUG) << "Point is between zero knots";
384  y=0;
385  } else {
386  // just 1 neighboring knot has y=0 - do a linear interpolation
387  LOG("Spline", pDEBUG)
388  << "Point has zero" << (is0n ? " left " : " right ") << "knot";
389  double xpknot=0, ypknot=0, xnknot=0, ynknot=0;
390  this->FindClosestKnot(x, xnknot, ynknot, "-");
391  this->FindClosestKnot(x, xpknot, ypknot, "+");
392  if(is0n) y = ypknot * (x-xnknot)/(xpknot-xnknot);
393  else y = ynknot * (x-xnknot)/(xpknot-xnknot);
394  }
395  }
396 
397  } else {
398  LOG("Spline", pDEBUG) << "x = " << x
399  << " is not within spline range [" << fXMin << ", " << fXMax << "]";
400  }
401 
402  if(y<0 && !fYCanBeNegative) {
403  LOG("Spline", pINFO) << "Negative y (" << y << ")";
404  LOG("Spline", pINFO) << "x = " << x;
405  LOG("Spline", pINFO) << "spline range [" << fXMin << ", " << fXMax << "]";
406  }
407 
408  LOG("Spline", pDEBUG) << "Spline(x = " << x << ") = " << y;
409 
410  return y;
411 }
412 //___________________________________________________________________________
414  string filename, string xtag, string ytag, string name) const
415 {
416  ofstream outxml(filename.c_str());
417  if(!outxml.is_open()) {
418  LOG("Spline", pERROR) << "Couldn't create file = " << filename;
419  return;
420  }
421  string spline_name = (name.size()>0 ? name : fName);
422 
423  this->SaveAsXml(outxml, xtag, ytag, spline_name);
424 
425  outxml << endl;
426  outxml.close();
427 }
428 //___________________________________________________________________________
430  ofstream & ofs, string xtag, string ytag, string name) const
431 {
432  string spline_name = (name.size()>0 ? name : fName);
433 
434  // create a spline tag with the number of knots as an attribute
435  int nknots = this->NKnots();
436  string padding = " ";
437  ofs << padding << "<spline name=\"" << spline_name
438  << "\" nknots=\"" << nknots << "\">" << endl;
439 
440  // start printing the knots
441  double x=0, y=0;
442  for(int iknot = 0; iknot < nknots; iknot++)
443  {
444  fInterpolator->GetKnot(iknot, x, y);
445 
446  ofs << std::fixed << setprecision(5);
447  ofs << "\t<knot>"
448  << " <" << xtag << "> " << setfill(' ')
449  << setw(10) << x << " </" << xtag << ">";
450  ofs << std::scientific << setprecision(10);
451  ofs << " <" << ytag << "> " << setfill(' ')
452  << setw(10) << y << " </" << ytag << ">"
453  << " </knot>" << endl;
454  }
455  ofs << padding << "</spline>" << endl;
456 }
457 //___________________________________________________________________________
458 void Spline::SaveAsText(string filename, string format) const
459 {
460  ofstream outtxt(filename.c_str());
461  if(!outtxt.is_open()) {
462  LOG("Spline", pERROR) << "Couldn't create file = " << filename;
463  return;
464  }
465  int nknots = this->NKnots();
466  outtxt << nknots << endl;
467 
468  double x=0, y=0;
469  for(int iknot = 0; iknot < nknots; iknot++) {
470  fInterpolator->GetKnot(iknot, x, y);
471  char line[1024];
472  sprintf(line,format.c_str(),x,y);
473  outtxt << line << endl;
474  }
475  outtxt << endl;
476  outtxt.close();
477 }
478 //___________________________________________________________________________
479 void Spline::SaveAsROOT(string filename, string name, bool recreate) const
480 {
481  string spline_name = (name.size()>0 ? name : fName);
482 
483  string opt = ( (recreate) ? "RECREATE" : "UPDATE" );
484 
485  TFile f(filename.c_str(), opt.c_str());
486  if(fInterpolator) fInterpolator->Write(spline_name.c_str());
487  f.Close();
488 }
489 //___________________________________________________________________________
491  int np, bool scale_with_x, bool in_log, double fx, double fy) const
492 {
493  double xmin = fXMin;
494  double xmax = fXMax;
495 
496  np = TMath::Max(np,2);
497 
498  bool use_log = in_log && (fXMin>0) && (fXMax>0);
499 
500  if(use_log) {
501  xmin = TMath::Log10(xmin);
502  xmax = TMath::Log10(xmax);
503  }
504 
505  double dx = (xmax - xmin) / (np-1);
506 
507  double * x = new double[np];
508  double * y = new double[np];
509 
510  for(int i=0; i<np; i++) {
511  x[i] = ( (use_log) ? TMath::Power(10, xmin+i*dx) : xmin + i*dx );
512  y[i] = this->Evaluate( x[i] );
513 
514  // scale with x if needed
515  if (scale_with_x) y[i] /= x[i];
516 
517  // apply x,y scaling
518  y[i] *= fy;
519  x[i] *= fx;
520  }
521 
522  TGraph * graph = new TGraph(np, x, y);
523  delete[] x;
524  delete[] y;
525  return graph;
526 }
527 //___________________________________________________________________________
529  double x, double & xknot, double & yknot, Option_t * opt) const
530 {
531  string option(opt);
532 
533  bool pos = (option.find("+") != string::npos);
534  bool neg = (option.find("-") != string::npos);
535 
536  if(!pos && !neg) return;
537 
538  int iknot = fInterpolator->FindX(x);
539 
540  double xp=0, yp=0, xn=0, yn=0;
541  fInterpolator->GetKnot(iknot, xn,yn);
542  fInterpolator->GetKnot(iknot+1,xp,yp);
543 
544  bool p = (TMath::Abs(x-xp) < TMath::Abs(x-xn));
545 
546  if(pos&&neg) {
547  if(p) { xknot = xp; yknot = yp; }
548  else { xknot = xn; yknot = yn; }
549  } else {
550  if(pos) { xknot = xp; yknot = yp; }
551  if(neg) { xknot = xn; yknot = yn; }
552  }
553 }
554 //___________________________________________________________________________
555 bool Spline::ClosestKnotValueIsZero(double x, Option_t * opt) const
556 {
557  double xknot = 0, yknot = 0;
558  this->FindClosestKnot(x, xknot, yknot, opt);
559  if(utils::math::AreEqual(yknot,0)) return true;
560  return false;
561 }
562 //___________________________________________________________________________
563 void Spline::Print(ostream & stream) const
564 {
565  int nknots = this->NKnots();
566  double xmin = this->XMin();
567  double xmax = this->XMax();
568 
569  stream << endl;
570  stream << "** Spline: " << this->Name() << endl;
571  stream << "Has " << nknots
572  << " knots in the [" << xmin << ", " << xmax << "] range" << endl;
573  double x,y;
574  for(int i=0; i<nknots; i++) {
575  this->GetKnot(i,x,y);
576  stream << "-- knot : " << i
577  << " -> (x = " << x << ", y = " << y << ")" << endl;
578  }
579 }
580 //___________________________________________________________________________
581 void Spline::Add(const Spline & spl, double c)
582 {
583  // continue only if the input spline is defined at a wider x-range
584  double xmin = this->XMin();
585  double xmax = this->XMax();
586  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
587  if(!ok) {
588  LOG("Spline", pERROR) << "** Can't add splines. X-range mismatch!";
589  return;
590  }
591 
592  int nknots = this->NKnots();
593  double * x = new double[nknots];
594  double * y = new double[nknots];
595 
596  for(int i=0; i<nknots; i++) {
597  this->GetKnot(i,x[i],y[i]);
598  y[i] += (c * spl.Evaluate(x[i]));
599  }
600  this->ResetSpline();
601  this->BuildSpline(nknots,x,y);
602  delete [] x;
603  delete [] y;
604 }
605 //___________________________________________________________________________
606 void Spline::Multiply(const Spline & spl, double c)
607 {
608  // continue only if the input spline is defined at a wider x-range
609  double xmin = this->XMin();
610  double xmax = this->XMax();
611  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
612  if(!ok) {
613  LOG("Spline", pERROR) << "** Can't multiply splines. X-range mismatch!";
614  return;
615  }
616 
617  int nknots = this->NKnots();
618  double * x = new double[nknots];
619  double * y = new double[nknots];
620 
621  for(int i=0; i<nknots; i++) {
622  this->GetKnot(i,x[i],y[i]);
623  y[i] *= (c * spl.Evaluate(x[i]));
624  }
625  this->ResetSpline();
626  this->BuildSpline(nknots,x,y);
627  delete [] x;
628  delete [] y;
629 }
630 //___________________________________________________________________________
631 void Spline::Divide(const Spline & spl, double c)
632 {
633  // continue only if the input spline is defined at a wider x-range
634  double xmin = this->XMin();
635  double xmax = this->XMax();
636  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
637  if(!ok) {
638  LOG("Spline", pERROR) << "** Can't divide splines. X-range mismatch!";
639  return;
640  }
641 
642  int nknots = this->NKnots();
643  double * x = new double[nknots];
644  double * y = new double[nknots];
645 
646  for(int i=0; i<nknots; i++) {
647  this->GetKnot(i,x[i],y[i]);
648  double denom = c * spl.Evaluate(x[i]);
649  bool denom_is_zero = TMath::Abs(denom) < DBL_EPSILON;
650  if(denom_is_zero) {
651  LOG("Spline", pERROR) << "** Refusing to divide spline knot by 0";
652  delete [] x;
653  delete [] y;
654  return;
655  }
656  y[i] /= denom;
657  }
658  this->ResetSpline();
659  this->BuildSpline(nknots,x,y);
660  delete [] x;
661  delete [] y;
662 }
663 //___________________________________________________________________________
664 void Spline::Add(double a)
665 {
666  int nknots = this->NKnots();
667  double * x = new double[nknots];
668  double * y = new double[nknots];
669 
670  for(int i=0; i<nknots; i++) {
671  this->GetKnot(i,x[i],y[i]);
672  y[i]+=a;
673  }
674  this->ResetSpline();
675  this->BuildSpline(nknots,x,y);
676  delete [] x;
677  delete [] y;
678 }
679 //___________________________________________________________________________
680 void Spline::Multiply(double a)
681 {
682  int nknots = this->NKnots();
683  double * x = new double[nknots];
684  double * y = new double[nknots];
685 
686  for(int i=0; i<nknots; i++) {
687  this->GetKnot(i,x[i],y[i]);
688  y[i]*=a;
689  }
690  this->ResetSpline();
691  this->BuildSpline(nknots,x,y);
692  delete [] x;
693  delete [] y;
694 }
695 //___________________________________________________________________________
696 void Spline::Divide(double a)
697 {
698  bool a_is_zero = TMath::Abs(a) < DBL_EPSILON;
699  if(a_is_zero==0) {
700  LOG("Spline", pERROR) << "** Refusing to divide spline by 0";
701  return;
702  }
703  int nknots = this->NKnots();
704  double * x = new double[nknots];
705  double * y = new double[nknots];
706 
707  for(int i=0; i<nknots; i++) {
708  this->GetKnot(i,x[i],y[i]);
709  y[i]/=a;
710  }
711  this->ResetSpline();
712  this->BuildSpline(nknots,x,y);
713  delete [] x;
714  delete [] y;
715 }
716 //___________________________________________________________________________
718 {
719  LOG("Spline", pDEBUG) << "Initializing spline...";
720 
721  fName = "genie-spline";
722  fXMin = 0.0;
723  fXMax = 0.0;
724  fYMax = 0.0;
725 
726  fInterpolator = 0;
727 
728  fYCanBeNegative = false;
729 
730  LOG("Spline", pDEBUG) << "...done initializing spline";
731 }
732 //___________________________________________________________________________
734 {
735  if(fInterpolator) delete fInterpolator;
736  this->InitSpline();
737 }
738 //___________________________________________________________________________
739 void Spline::BuildSpline(int nentries, double x[], double y[])
740 {
741  LOG("Spline", pDEBUG) << "Building spline...";
742 
743  double xmin = x[0]; // minimum x in spline
744  double xmax = x[nentries-1]; // maximum x in spline
745 
746  fNKnots = nentries;
747  fXMin = xmin;
748  fXMax = xmax;
749  fYMax = y[ TMath::LocMax(nentries, y) ]; // maximum y in spline
750 
751  if(fInterpolator) delete fInterpolator;
752 
753  fInterpolator = new TSpline3("spl3", x, y, nentries, "0");
754 
755  LOG("Spline", pDEBUG) << "...done building spline";
756 }
757 //___________________________________________________________________________
const XML_Char * name
Definition: expat.h:151
TGraph * GetAsTGraph(int np=500, bool xscaling=false, bool inlog=false, double fx=1., double fy=1.) const
Definition: Spline.cxx:490
std::map< std::string, double > xmax
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
TTree * ntuple
Definition: macro.C:6
bool AreEqual(double x1, double x2)
Definition: MathUtils.cxx:243
string TrimSpaces(xmlChar *xmls)
const char * p
Definition: xmltok.h:285
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
void Print(ostream &stream) const
Definition: Spline.cxx:563
bool LoadFromAsciiFile(string filename)
Definition: Spline.cxx:241
void BuildSpline(int nentries, double x[], double y[])
Definition: Spline.cxx:739
string filename
Definition: shutoffs.py:106
bool fYCanBeNegative
Definition: Spline.h:131
double Evaluate(double x) const
Definition: Spline.cxx:362
bool ClosestKnotValueIsZero(double x, Option_t *opt="-+") const
Definition: Spline.cxx:555
bool IsWithinValidRange(double x) const
Definition: Spline.cxx:356
TSpline3 * GetAsTSpline(void) const
Definition: Spline.h:97
bool LoadFromXmlFile(string filename, string xtag, string ytag)
Definition: Spline.cxx:154
double GetKnotX(int iknot) const
Definition: Spline.cxx:334
double dx[NP][NC]
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ResetSpline(void)
Definition: Spline.cxx:733
const Int_t nknots
Definition: testXsec.C:24
Long64_t nentries
const double a
ClassImp(Spline) namespace genie
Definition: Spline.cxx:41
double fXMax
Definition: Spline.h:128
void SaveAsROOT(string filename, string name="", bool recreate=false) const
Definition: Spline.cxx:479
#define pINFO
Definition: Messenger.h:63
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
void Add(const Spline &spl, double c=1)
Definition: Spline.cxx:581
void Multiply(const Spline &spl, double c=1)
Definition: Spline.cxx:606
void GetKnot(int iknot, double &x, double &y) const
Definition: Spline.cxx:325
#define pWARN
Definition: Messenger.h:61
double fYMax
Definition: Spline.h:129
string TrimSpaces(string input)
Definition: StringUtils.cxx:24
double GetKnotY(int iknot) const
Definition: Spline.cxx:345
double XMax(void) const
Definition: Spline.h:78
string fName
Definition: Spline.h:125
bool LoadFromTSpline3(const TSpline3 &spline, int nknots)
Definition: Spline.cxx:304
void SaveAsText(string filename, string format="%10.6f\t%10.6f") const
Definition: Spline.cxx:458
void SaveAsXml(string filename, string xtag, string ytag, string name="") const
Definition: Spline.cxx:413
const Cut cut
Definition: exporter_fd.C:30
int fNKnots
Definition: Spline.h:126
TSpline3 * fInterpolator
Definition: Spline.h:130
bool LoadFromNtuple(TNtupleD *nt, string xy, string cut="")
Definition: Spline.cxx:252
bool LoadFromDBase(TSQLServer *db, string query)
Definition: Spline.cxx:298
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
string Name(void) const
Definition: Spline.h:84
int NKnots(void) const
Definition: Spline.h:73
void InitSpline(void)
Definition: Spline.cxx:717
assert(nhit_max >=nhit_nbins)
bool LoadFromTree(TTree *tr, string xy, string cut="")
Definition: Spline.cxx:260
virtual ~Spline()
Definition: Spline.cxx:149
double fXMin
Definition: Spline.h:127
void FindClosestKnot(double x, double &xknot, double &yknot, Option_t *opt="-+") const
Definition: Spline.cxx:528
double XMin(void) const
Definition: Spline.h:77
TNtuple * nt
Definition: drawXsec.C:2
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
void Divide(const Spline &spl, double c=1)
Definition: Spline.cxx:631
#define pDEBUG
Definition: Messenger.h:64