[root] / trunk / hist / hist / src / TMultiGraph.cxx Repository:
ViewVC logotype

Annotation of /trunk/hist/hist/src/TMultiGraph.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 36493 - (view) (download) (as text)

1 : brun 24702 // @(#)root/hist:$Id$
2 : brun 754 // Author: Rene Brun 12/10/2000
3 :    
4 :     /*************************************************************************
5 :     * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 :     * All rights reserved. *
7 :     * *
8 :     * For the licensing terms see $ROOTSYS/LICENSE. *
9 :     * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 :     *************************************************************************/
11 :    
12 :     #include "TROOT.h"
13 :     #include "TMultiGraph.h"
14 :     #include "TGraph.h"
15 :     #include "TH1.h"
16 :     #include "TVirtualPad.h"
17 : rdm 3748 #include "Riostream.h"
18 : brun 11226 #include "TVirtualFitter.h"
19 : moneta 25487 #include "TPluginManager.h"
20 : pcanal 14336 #include "TClass.h"
21 : brun 17338 #include "TMath.h"
22 : moneta 25487 #include "TSystem.h"
23 : rdm 22419 #include <stdlib.h>
24 : brun 754
25 : moneta 31207 #include "HFitInterface.h"
26 :     #include "Fit/DataRange.h"
27 :     #include "Math/MinimizerOptions.h"
28 :    
29 : brun 754 #include <ctype.h>
30 :    
31 : brun 11226 extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
32 : brun 754
33 :     ClassImp(TMultiGraph)
34 :    
35 : couet 13304
36 : brun 754 //______________________________________________________________________________
37 : couet 20297 /* Begin_Html
38 :     <center><h2>TMultiGraph class</h2></center>
39 :     A TMultiGraph is a collection of TGraph (or derived) objects
40 :     Use <tt>TMultiGraph::Add</tt> to add a new graph to the list.
41 :     The TMultiGraph owns the objects in the list.
42 :     Drawing options are the same as for TGraph.
43 :     <p>
44 :     Example:
45 : couet 20433 <pre>
46 : couet 20297 TGraph *gr1 = new TGraph(...
47 :     TGraphErrors *gr2 = new TGraphErrors(...
48 :     TMultiGraph *mg = new TMultiGraph();
49 :     mg->Add(gr1,"lp");
50 :     mg->Add(gr2,"cp");
51 :     mg->Draw("a");
52 :     </pre>
53 :     The drawing option for each TGraph may be specified as an optional
54 :     second argument of the Add function.
55 :     If a draw option is specified, it will be used to draw the graph,
56 :     otherwise the graph will be drawn with the option specified in
57 :     <tt>TMultiGraph::Draw</tt>.
58 :     End_Html */
59 : brun 754
60 : couet 13304
61 : brun 754 //______________________________________________________________________________
62 :     TMultiGraph::TMultiGraph(): TNamed()
63 :     {
64 : couet 13304 // TMultiGraph default constructor
65 : brun 754
66 :     fGraphs = 0;
67 : brun 11226 fFunctions = 0;
68 : brun 754 fHistogram = 0;
69 :     fMaximum = -1111;
70 :     fMinimum = -1111;
71 :     }
72 :    
73 : couet 13304
74 : brun 754 //______________________________________________________________________________
75 :     TMultiGraph::TMultiGraph(const char *name, const char *title)
76 :     : TNamed(name,title)
77 :     {
78 : couet 13304 // constructor with name and title
79 :    
80 : brun 754 fGraphs = 0;
81 : brun 11226 fFunctions = 0;
82 : brun 754 fHistogram = 0;
83 :     fMaximum = -1111;
84 :     fMinimum = -1111;
85 :     }
86 :    
87 : couet 20433
88 : brun 15134 //______________________________________________________________________________
89 :     TMultiGraph::TMultiGraph(const TMultiGraph& mg) :
90 :     TNamed (mg),
91 :     fGraphs(mg.fGraphs),
92 :     fFunctions(mg.fFunctions),
93 :     fHistogram(mg.fHistogram),
94 :     fMaximum(mg.fMaximum),
95 :     fMinimum(mg.fMinimum)
96 : couet 20433 {
97 : brun 15171 //copy constructor
98 :     }
99 : couet 13304
100 : couet 20433
101 : brun 754 //______________________________________________________________________________
102 : brun 15134 TMultiGraph& TMultiGraph::operator=(const TMultiGraph& mg)
103 :     {
104 : brun 15740 //assignement operator
105 : brun 15171 if(this!=&mg) {
106 :     TNamed::operator=(mg);
107 :     fGraphs=mg.fGraphs;
108 :     fFunctions=mg.fFunctions;
109 :     fHistogram=mg.fHistogram;
110 :     fMaximum=mg.fMaximum;
111 :     fMinimum=mg.fMinimum;
112 : couet 20433 }
113 : brun 15171 return *this;
114 : brun 15134 }
115 :    
116 : couet 20433
117 : brun 15134 //______________________________________________________________________________
118 : brun 754 TMultiGraph::~TMultiGraph()
119 :     {
120 : couet 13304 // TMultiGraph destructor
121 : brun 754
122 :     if (!fGraphs) return;
123 : brun 5552 TGraph *g;
124 :     TIter next(fGraphs);
125 :     while ((g = (TGraph*) next())) {
126 : couet 13304 g->ResetBit(kMustCleanup);
127 : brun 5552 }
128 : brun 754 fGraphs->Delete();
129 :     delete fGraphs;
130 :     fGraphs = 0;
131 :     delete fHistogram;
132 :     fHistogram = 0;
133 : brun 11226 if (fFunctions) {
134 :     fFunctions->SetBit(kInvalidObject);
135 :     //special logic to support the case where the same object is
136 :     //added multiple times in fFunctions.
137 :     //This case happens when the same object is added with different
138 :     //drawing modes
139 :     TObject *obj;
140 :     while ((obj = fFunctions->First())) {
141 : rdm 22546 while(fFunctions->Remove(obj)) { }
142 : brun 11226 delete obj;
143 :     }
144 :     delete fFunctions;
145 :     }
146 : brun 754 }
147 :    
148 : couet 13304
149 : brun 754 //______________________________________________________________________________
150 : brun 4037 void TMultiGraph::Add(TGraph *graph, Option_t *chopt)
151 : brun 754 {
152 :     // add a new graph to the list of graphs
153 : brun 6639 // note that the graph is now owned by the TMultigraph.
154 :     // Deleting the TMultiGraph object will automatically delete the graphs.
155 :     // You should not delete the graphs when the TMultigraph is still active.
156 : rdm 3742
157 : brun 754 if (!fGraphs) fGraphs = new TList();
158 : brun 5552 graph->SetBit(kMustCleanup);
159 : brun 4037 fGraphs->Add(graph,chopt);
160 : brun 754 }
161 :    
162 : couet 13304
163 : brun 754 //______________________________________________________________________________
164 : couet 16436 void TMultiGraph::Add(TMultiGraph *multigraph, Option_t *chopt)
165 :     {
166 :     // add all the graphs in "multigraph" to the list of graphs.
167 :    
168 :     TList *graphlist = multigraph->GetListOfGraphs();
169 :     if (!graphlist) return;
170 :    
171 :     if (!fGraphs) fGraphs = new TList();
172 : couet 20433
173 : couet 16436 TGraph *gr;
174 :     gr = (TGraph*)graphlist->First();
175 :     fGraphs->Add(gr,chopt);
176 : couet 22379 for(Int_t i = 1; i < graphlist->GetSize(); i++){
177 : couet 16436 gr = (TGraph*)graphlist->After(gr);
178 :     fGraphs->Add(gr,chopt);
179 :     }
180 :     }
181 :    
182 :    
183 :     //______________________________________________________________________________
184 : brun 754 void TMultiGraph::Browse(TBrowser *)
185 :     {
186 : couet 13304 // Browse multigraph.
187 :    
188 :     Draw("alp");
189 :     gPad->Update();
190 : brun 754 }
191 :    
192 : couet 13304
193 : brun 754 //______________________________________________________________________________
194 :     Int_t TMultiGraph::DistancetoPrimitive(Int_t px, Int_t py)
195 :     {
196 : couet 13304 // Compute distance from point px,py to each graph
197 : brun 754
198 : couet 13304 // Are we on the axis?
199 : brun 754 const Int_t kMaxDiff = 10;
200 :     Int_t distance = 9999;
201 :     if (fHistogram) {
202 :     distance = fHistogram->DistancetoPrimitive(px,py);
203 :     if (distance <= 0) return distance;
204 :     }
205 :    
206 : couet 13304 // Loop on the list of graphs
207 : brun 754 if (!fGraphs) return distance;
208 :     TGraph *g;
209 :     TIter next(fGraphs);
210 :     while ((g = (TGraph*) next())) {
211 :     Int_t dist = g->DistancetoPrimitive(px,py);
212 : brun 4915 if (dist <= 0) return 0;
213 : brun 754 if (dist < kMaxDiff) {gPad->SetSelected(g); return dist;}
214 :     }
215 :     return distance;
216 :     }
217 :    
218 : couet 13304
219 : brun 754 //______________________________________________________________________________
220 :     void TMultiGraph::Draw(Option_t *option)
221 :     {
222 : couet 13304 // Draw this multigraph with its current attributes.
223 :     //
224 :     // Options to draw a graph are described in TGraph::PainGraph
225 :     //
226 :     // The drawing option for each TGraph may be specified as an optional
227 :     // second argument of the Add function. You can use GetGraphDrawOption
228 :     // to return this option.
229 :     // If a draw option is specified, it will be used to draw the graph,
230 :     // otherwise the graph will be drawn with the option specified in
231 :     // TMultiGraph::Draw. Use GetDrawOption to return the option specified
232 :     // when drawin the TMultiGraph.
233 : brun 754
234 : couet 13304 AppendPad(option);
235 : brun 754 }
236 :    
237 : couet 13304
238 : brun 754 //______________________________________________________________________________
239 : moneta 31491 TFitResultPtr TMultiGraph::Fit(const char *fname, Option_t *option, Option_t *, Axis_t xmin, Axis_t xmax)
240 : brun 11226 {
241 : couet 13304 // Fit this graph with function with name fname.
242 :     //
243 :     // interface to TF1::Fit(TF1 *f1...
244 : brun 11226
245 :     char *linear;
246 : brun 11231 linear= (char*)strstr(fname, "++");
247 : brun 11226 TF1 *f1=0;
248 :     if (linear)
249 : brun 12643 f1=new TF1(fname, fname, xmin, xmax);
250 : brun 11226 else {
251 :     f1 = (TF1*)gROOT->GetFunction(fname);
252 :     if (!f1) { Printf("Unknown function: %s",fname); return -1; }
253 :     }
254 :    
255 :     return Fit(f1,option,"",xmin,xmax);
256 :     }
257 :    
258 : couet 13304
259 : brun 11226 //______________________________________________________________________________
260 : moneta 31491 TFitResultPtr TMultiGraph::Fit(TF1 *f1, Option_t *option, Option_t *goption, Axis_t rxmin, Axis_t rxmax)
261 : brun 11226 {
262 : couet 13304 // Fit this multigraph with function f1.
263 :     //
264 :     // In this function all graphs of the multigraph are fitted simultaneously
265 :     //
266 :     // f1 is an already predefined function created by TF1.
267 :     // Predefined functions such as gaus, expo and poln are automatically
268 :     // created by ROOT.
269 :     //
270 :     // The list of fit options is given in parameter option.
271 :     // option = "W" Set all errors to 1
272 :     // = "U" Use a User specified fitting algorithm (via SetFCN)
273 :     // = "Q" Quiet mode (minimum printing)
274 :     // = "V" Verbose mode (default is between Q and V)
275 :     // = "B" Use this option when you want to fix one or more parameters
276 :     // and the fitting function is like "gaus","expo","poln","landau".
277 :     // = "R" Use the Range specified in the function range
278 :     // = "N" Do not store the graphics function, do not draw
279 :     // = "0" Do not plot the result of the fit. By default the fitted function
280 :     // is drawn unless the option"N" above is specified.
281 :     // = "+" Add this new fitted function to the list of fitted functions
282 :     // (by default, any previous function is deleted)
283 :     // = "C" In case of linear fitting, not calculate the chisquare
284 :     // (saves time)
285 :     // = "F" If fitting a polN, switch to minuit fitter
286 :     // = "ROB" In case of linear fitting, compute the LTS regression
287 : couet 20433 // coefficients (robust(resistant) regression), using
288 : couet 13304 // the default fraction of good points
289 :     // "ROB=0.x" - compute the LTS regression coefficients, using
290 :     // 0.x as a fraction of good points
291 :     //
292 :     // When the fit is drawn (by default), the parameter goption may be used
293 :     // to specify a list of graphics options. See TGraph::Paint for a complete
294 :     // list of these options.
295 :     //
296 :     // In order to use the Range option, one must first create a function
297 :     // with the expression to be fitted. For example, if your graph
298 :     // has a defined range between -4 and 4 and you want to fit a gaussian
299 :     // only in the interval 1 to 3, you can do:
300 :     // TF1 *f1 = new TF1("f1","gaus",1,3);
301 :     // graph->Fit("f1","R");
302 :     //
303 :     //
304 :     // who is calling this function
305 :     // ============================
306 :     // Note that this function is called when calling TGraphErrors::Fit
307 :     // or TGraphAsymmErrors::Fit ot TGraphBentErrors::Fit
308 :     // see the discussion below on the errors calulation.
309 :     //
310 :     // Setting initial conditions
311 :     // ==========================
312 :     // Parameters must be initialized before invoking the Fit function.
313 :     // The setting of the parameter initial values is automatic for the
314 :     // predefined functions : poln, expo, gaus, landau. One can however disable
315 :     // this automatic computation by specifying the option "B".
316 :     // You can specify boundary limits for some or all parameters via
317 :     // f1->SetParLimits(p_number, parmin, parmax);
318 :     // if parmin>=parmax, the parameter is fixed
319 :     // Note that you are not forced to fix the limits for all parameters.
320 :     // For example, if you fit a function with 6 parameters, you can do:
321 :     // func->SetParameters(0,3.1,1.e-6,0.1,-8,100);
322 :     // func->SetParLimits(4,-10,-4);
323 :     // func->SetParLimits(5, 1,1);
324 :     // With this setup, parameters 0->3 can vary freely
325 :     // Parameter 4 has boundaries [-10,-4] with initial value -8
326 :     // Parameter 5 is fixed to 100.
327 :     //
328 :     // Fit range
329 :     // =========
330 :     // The fit range can be specified in two ways:
331 :     // - specify rxmax > rxmin (default is rxmin=rxmax=0)
332 :     // - specify the option "R". In this case, the function will be taken
333 :     // instead of the full graph range.
334 :     //
335 : moneta 34699 // Changing the fitting function
336 :     // =============================
337 :     // By default a chi2 fitting function is used for fitting the TGraphs's.
338 :     // The function is implemented in FitUtil::EvaluateChi2.
339 :     // In case of TGraphErrors an effective chi2 is used
340 :     // (see TGraphErrors fit in TGraph::Fit) and is implemented in
341 :     // FitUtil::EvaluateChi2Effective
342 :     // To specify a User defined fitting function, specify option "U" and
343 :     // call the following functions:
344 :     // TVirtualFitter::Fitter(mygraph)->SetFCN(MyFittingFunction)
345 :     // where MyFittingFunction is of type:
346 :     // extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
347 : couet 13304 //
348 : moneta 34699 // Access to the fit result
349 :     // ========================
350 :     // The function returns a TFitResultPtr which can hold a pointer to a TFitResult object.
351 :     // By default the TFitResultPtr contains only the status of the fit and it converts
352 :     // automatically to an integer. If the option "S" is instead used, TFitResultPtr contains
353 :     // the TFitResult and behaves as a smart pointer to it. For example one can do:
354 :     // TFitResultPtr r = graph->Fit("myFunc","S");
355 :     // TMatrixDSym cov = r->GetCovarianceMatrix(); // to access the covariance matrix
356 :     // Double_t par0 = r->Parameter(0); // retrieve the value for the parameter 0
357 :     // Double_t err0 = r->ParError(0); // retrieve the error for the parameter 0
358 :     // r->Print("V"); // print full information of fit including covariance matrix
359 :     // r->Write(); // store the result in a file
360 : couet 13304 //
361 : moneta 34699 // The fit parameters, error and chi2 (but not covariance matrix) can be retrieved also
362 :     // from the fitted function.
363 : couet 13304 //
364 :     //
365 :     // Associated functions
366 :     // ====================
367 :     // One or more object (typically a TF1*) can be added to the list
368 :     // of functions (fFunctions) associated to each graph.
369 :     // When TGraph::Fit is invoked, the fitted function is added to this list.
370 :     // Given a graph gr, one can retrieve an associated function
371 :     // with: TF1 *myfunc = gr->GetFunction("myfunc");
372 :     //
373 :     // If the graph is made persistent, the list of
374 :     // associated functions is also persistent. Given a pointer (see above)
375 :     // to an associated function myfunc, one can retrieve the function/fit
376 :     // parameters with calls such as:
377 :     // Double_t chi2 = myfunc->GetChisquare();
378 :     // Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
379 :     // Double_t err0 = myfunc->GetParError(0); //error on first parameter
380 :     //
381 :     // Fit Statistics
382 :     // ==============
383 :     // You can change the statistics box to display the fit parameters with
384 :     // the TStyle::SetOptFit(mode) method. This mode has four digits.
385 :     // mode = pcev (default = 0111)
386 :     // v = 1; print name/values of parameters
387 :     // e = 1; print errors (if e=1, v must be 1)
388 :     // c = 1; print Chisquare/Number of degress of freedom
389 :     // p = 1; print Probability
390 :     //
391 :     // For example: gStyle->SetOptFit(1011);
392 :     // prints the fit probability, parameter names/values, and errors.
393 :     // You can change the position of the statistics box with these lines
394 :     // (where g is a pointer to the TGraph):
395 :     //
396 :     // Root > TPaveStats *st = (TPaveStats*)g->GetListOfFunctions()->FindObject("stats")
397 :     // Root > st->SetX1NDC(newx1); //new x start position
398 :     // Root > st->SetX2NDC(newx2); //new x end position
399 : brun 11226
400 : moneta 31207 // internal multigraph fitting methods
401 :     Foption_t fitOption;
402 :     ROOT::Fit::FitOptionsMake(option,fitOption);
403 : rdm 11750
404 : moneta 31207 // create range and minimizer options with default values
405 :     ROOT::Fit::DataRange range(rxmin,rxmax);
406 :     ROOT::Math::MinimizerOptions minOption;
407 :     return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range);
408 :    
409 : moneta 25487 }
410 : rdm 11750
411 : moneta 25487 //______________________________________________________________________________
412 :     void TMultiGraph::FitPanel()
413 :     {
414 :     // -*-*-*-*-*Display a panel with all histogram fit options*-*-*-*-*-*
415 :     // ==============================================
416 :     //
417 :     // See class TFitPanel for example
418 : rdm 11750
419 : moneta 25487 if (!gPad)
420 :     gROOT->MakeDefCanvas();
421 : brun 11226
422 : moneta 25487 if (!gPad) {
423 :     Error("FitPanel", "Unable to create a default canvas");
424 :     return;
425 : brun 12136 }
426 : brun 11226
427 : moneta 25487 // use plugin manager to create instance of TFitEditor
428 :     TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
429 :     if (handler && handler->LoadPlugin() != -1) {
430 :     if (handler->ExecPlugin(2, gPad, this) == 0)
431 :     Error("FitPanel", "Unable to crate the FitPanel");
432 : brun 11226 }
433 : couet 26029 else
434 : moneta 25487 Error("FitPanel", "Unable to find the FitPanel plug-in");
435 : brun 11226 }
436 :    
437 :     //______________________________________________________________________________
438 : brun 11584 Option_t *TMultiGraph::GetGraphDrawOption(const TGraph *gr) const
439 :     {
440 : couet 13304 // Return the draw option for the TGraph gr in this TMultiGraph
441 :     // The return option is the one specified when calling TMultiGraph::Add(gr,option).
442 : rdm 11750
443 : brun 11584 if (!fGraphs || !gr) return "";
444 :     TListIter next(fGraphs);
445 :     TObject *obj;
446 :     while ((obj = next())) {
447 : brun 11600 if (obj == (TObject*)gr) return next.GetOption();
448 : brun 11584 }
449 :     return "";
450 :     }
451 :    
452 : couet 13304
453 : brun 11584 //______________________________________________________________________________
454 : brun 11226 void TMultiGraph::InitGaus(Double_t xmin, Double_t xmax)
455 :     {
456 : couet 13304 // Compute Initial values of parameters for a gaussian.
457 : brun 11226
458 :     Double_t allcha, sumx, sumx2, x, val, rms, mean;
459 :     Int_t bin;
460 :     const Double_t sqrtpi = 2.506628;
461 :    
462 : couet 13304 // Compute mean value and RMS of the graph in the given range
463 : brun 11226 Int_t np = 0;
464 :     allcha = sumx = sumx2 = 0;
465 :     TGraph *g;
466 :     TIter next(fGraphs);
467 :     Double_t *px, *py;
468 :     Int_t npp; //number of points in each graph
469 :     while ((g = (TGraph*) next())) {
470 :     px=g->GetX();
471 :     py=g->GetY();
472 :     npp=g->GetN();
473 :     for (bin=0; bin<npp; bin++){
474 : brun 12643 x=px[bin];
475 :     if (x<xmin || x>xmax) continue;
476 :     np++;
477 :     val=py[bin];
478 :     sumx+=val*x;
479 :     sumx2+=val*x*x;
480 :     allcha+=val;
481 : brun 11226 }
482 :     }
483 :     if (np == 0 || allcha == 0) return;
484 :     mean = sumx/allcha;
485 :     rms = TMath::Sqrt(sumx2/allcha - mean*mean);
486 :    
487 :     Double_t binwidx = TMath::Abs((xmax-xmin)/np);
488 :     if (rms == 0) rms = 1;
489 :     TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
490 :     TF1 *f1 = (TF1*)grFitter->GetUserFunc();
491 :     f1->SetParameter(0,binwidx*allcha/(sqrtpi*rms));
492 :     f1->SetParameter(1,mean);
493 :     f1->SetParameter(2,rms);
494 :     f1->SetParLimits(2,0,10*rms);
495 :     }
496 :    
497 : couet 13304
498 : brun 11226 //______________________________________________________________________________
499 :     void TMultiGraph::InitExpo(Double_t xmin, Double_t xmax)
500 :     {
501 : couet 13304 // Compute Initial values of parameters for an exponential.
502 : brun 11226
503 :     Double_t constant, slope;
504 :     Int_t ifail;
505 :    
506 :     LeastSquareLinearFit(-1, constant, slope, ifail, xmin, xmax);
507 :    
508 :     TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
509 :     TF1 *f1 = (TF1*)grFitter->GetUserFunc();
510 :     f1->SetParameter(0,constant);
511 :     f1->SetParameter(1,slope);
512 :     }
513 :    
514 : couet 13304
515 : brun 11226 //______________________________________________________________________________
516 :     void TMultiGraph::InitPolynom(Double_t xmin, Double_t xmax)
517 :     {
518 : couet 13304 // Compute Initial values of parameters for a polynom.
519 : brun 11226
520 :     Double_t fitpar[25];
521 :    
522 :     TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
523 :     TF1 *f1 = (TF1*)grFitter->GetUserFunc();
524 :     Int_t npar = f1->GetNpar();
525 :    
526 :     LeastSquareFit(npar, fitpar, xmin, xmax);
527 :    
528 :     for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
529 :     }
530 :    
531 : couet 13304
532 : brun 11226 //______________________________________________________________________________
533 :     void TMultiGraph::LeastSquareFit(Int_t m, Double_t *a, Double_t xmin, Double_t xmax)
534 :     {
535 : couet 13304 // Least squares lpolynomial fitting without weights.
536 :     //
537 :     // m number of parameters
538 :     // a array of parameters
539 :     // first 1st point number to fit (default =0)
540 :     // last last point number to fit (default=fNpoints-1)
541 :     //
542 :     // based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
543 : brun 11226
544 : couet 13304 const Double_t zero = 0.;
545 :     const Double_t one = 1.;
546 :     const Int_t idim = 20;
547 : brun 11226
548 : couet 13304 Double_t b[400] /* was [20][20] */;
549 :     Int_t i, k, l, ifail, bin;
550 :     Double_t power;
551 :     Double_t da[20], xk, yk;
552 : brun 11226
553 :    
554 : couet 13304 //count the total number of points to fit
555 :     TGraph *g;
556 :     TIter next(fGraphs);
557 :     Double_t *px, *py;
558 :     Int_t n=0;
559 :     Int_t npp;
560 :     while ((g = (TGraph*) next())) {
561 :     px=g->GetX();
562 :     npp=g->GetN();
563 :     for (bin=0; bin<npp; bin++){
564 :     xk=px[bin];
565 :     if (xk < xmin || xk > xmax) continue;
566 :     n++;
567 :     }
568 :     }
569 :     if (m <= 2) {
570 :     LeastSquareLinearFit(n, a[0], a[1], ifail, xmin, xmax);
571 :     return;
572 :     }
573 :     if (m > idim || m > n) return;
574 :     da[0] = zero;
575 :     for (l = 2; l <= m; ++l) {
576 :     b[l-1] = zero;
577 :     b[m + l*20 - 21] = zero;
578 :     da[l-1] = zero;
579 :     }
580 :     Int_t np = 0;
581 : brun 11226
582 : couet 13304 next.Reset();
583 :     while ((g = (TGraph*) next())) {
584 :     px=g->GetX();
585 :     py=g->GetY();
586 :     npp=g->GetN();
587 :    
588 :     for (k = 0; k <= npp; ++k) {
589 :     xk = px[k];
590 :     if (xk < xmin || xk > xmax) continue;
591 :     np++;
592 :     yk = py[k];
593 :     power = one;
594 :     da[0] += yk;
595 :     for (l = 2; l <= m; ++l) {
596 :     power *= xk;
597 : brun 12643 b[l-1] += power;
598 :     da[l-1] += power*yk;
599 : couet 13304 }
600 :     for (l = 2; l <= m; ++l) {
601 :     power *= xk;
602 :     b[m + l*20 - 21] += power;
603 :     }
604 :     }
605 :     }
606 :     b[0] = Double_t(np);
607 :     for (i = 3; i <= m; ++i) {
608 :     for (k = i; k <= m; ++k) {
609 :     b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
610 :     }
611 :     }
612 :     H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
613 : brun 11226
614 : couet 13304 if (ifail < 0) {
615 :     //a[0] = fY[0];
616 :     py=((TGraph *)fGraphs->First())->GetY();
617 :     a[0]=py[0];
618 :     for (i=1; i<m; ++i) a[i] = 0;
619 :     return;
620 :     }
621 :     for (i=0; i<m; ++i) a[i] = da[i];
622 : brun 11226 }
623 :    
624 : couet 13304
625 : brun 11226 //______________________________________________________________________________
626 :     void TMultiGraph::LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail, Double_t xmin, Double_t xmax)
627 :     {
628 : couet 13304 // Least square linear fit without weights.
629 :     //
630 :     // Fit a straight line (a0 + a1*x) to the data in this graph.
631 :     // ndata: number of points to fit
632 :     // first: first point number to fit
633 :     // last: last point to fit O(ndata should be last-first
634 :     // ifail: return parameter indicating the status of the fit (ifail=0, fit is OK)
635 :     //
636 :     // extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
637 : brun 11226
638 : couet 13304 Double_t xbar, ybar, x2bar;
639 :     Int_t i;
640 :     Double_t xybar;
641 :     Double_t fn, xk, yk;
642 :     Double_t det;
643 : brun 11226
644 : couet 13304 ifail = -2;
645 :     xbar = ybar = x2bar = xybar = 0;
646 :     Int_t np = 0;
647 :     TGraph *g;
648 :     TIter next(fGraphs);
649 :     Double_t *px, *py;
650 :     Int_t npp;
651 :     while ((g = (TGraph*) next())) {
652 :     px=g->GetX();
653 :     py=g->GetY();
654 :     npp=g->GetN();
655 :     for (i = 0; i < npp; ++i) {
656 :     xk = px[i];
657 :     if (xk < xmin || xk > xmax) continue;
658 :     np++;
659 :     yk = py[i];
660 :     if (ndata < 0) {
661 :     if (yk <= 0) yk = 1e-9;
662 :     yk = TMath::Log(yk);
663 :     }
664 :     xbar += xk;
665 :     ybar += yk;
666 :     x2bar += xk*xk;
667 :     xybar += xk*yk;
668 :     }
669 :     }
670 :     fn = Double_t(np);
671 :     det = fn*x2bar - xbar*xbar;
672 :     ifail = -1;
673 :     if (det <= 0) {
674 :     if (fn > 0) a0 = ybar/fn;
675 :     else a0 = 0;
676 :     a1 = 0;
677 :     return;
678 :     }
679 :     ifail = 0;
680 :     a0 = (x2bar*ybar - xbar*xybar) / det;
681 :     a1 = (fn*xybar - xbar*ybar) / det;
682 :     }
683 : brun 11226
684 :    
685 :     //______________________________________________________________________________
686 : couet 36493 Int_t TMultiGraph::IsInside(Double_t x, Double_t y) const
687 :     {
688 :     // Return 1 if the point (x,y) is inside one of the graphs 0 otherwise.
689 :    
690 :     Int_t in = 0;
691 :     if (!fGraphs) return in;
692 :     TGraph *g;
693 :     TIter next(fGraphs);
694 :     while ((g = (TGraph*) next())) {
695 :     in = g->IsInside(x, y);
696 :     if (in) return in;
697 :     }
698 :     return in;
699 :     }
700 :    
701 :    
702 :     //______________________________________________________________________________
703 : brun 1205 TH1F *TMultiGraph::GetHistogram() const
704 : brun 754 {
705 : couet 13304 // Returns a pointer to the histogram used to draw the axis
706 :     // Takes into account the two following cases.
707 :     // 1- option 'A' was specified in TMultiGraph::Draw. Return fHistogram
708 :     // 2- user had called TPad::DrawFrame. return pointer to hframe histogram
709 : brun 754
710 :     if (fHistogram) return fHistogram;
711 :     if (!gPad) return 0;
712 :     gPad->Modified();
713 :     gPad->Update();
714 :     if (fHistogram) return fHistogram;
715 :     TH1F *h1 = (TH1F*)gPad->FindObject("hframe");
716 :     return h1;
717 :     }
718 :    
719 : couet 13304
720 : brun 754 //______________________________________________________________________________
721 : brun 11226 TF1 *TMultiGraph::GetFunction(const char *name) const
722 :     {
723 : couet 13304 // Return pointer to function with name.
724 :     //
725 :     // Functions such as TGraph::Fit store the fitted function in the list of
726 :     // functions of this graph.
727 : brun 11226
728 :     if (!fFunctions) return 0;
729 :     return (TF1*)fFunctions->FindObject(name);
730 :     }
731 :    
732 : moneta 25487 //______________________________________________________________________________
733 : couet 26029 TList *TMultiGraph::GetListOfFunctions()
734 : moneta 25487 {
735 :     // Return pointer to list of functions
736 :     // if pointer is null create the list
737 : couet 13304
738 : moneta 25487 if (!fFunctions) fFunctions = new TList();
739 : couet 26029 return fFunctions;
740 : moneta 25487 }
741 :    
742 :    
743 : brun 11226 //______________________________________________________________________________
744 : brun 1205 TAxis *TMultiGraph::GetXaxis() const
745 : brun 754 {
746 :     // Get x axis of the graph.
747 :    
748 :     if (!gPad) return 0;
749 : brun 4083 TH1 *h = GetHistogram();
750 :     if (!h) return 0;
751 :     return h->GetXaxis();
752 : brun 754 }
753 :    
754 : couet 13304
755 : brun 754 //______________________________________________________________________________
756 : brun 1205 TAxis *TMultiGraph::GetYaxis() const
757 : brun 754 {
758 :     // Get y axis of the graph.
759 :    
760 :     if (!gPad) return 0;
761 : brun 4083 TH1 *h = GetHistogram();
762 :     if (!h) return 0;
763 :     return h->GetYaxis();
764 : brun 754 }
765 :    
766 : couet 13304
767 : brun 754 //______________________________________________________________________________
768 :     void TMultiGraph::Paint(Option_t *option)
769 :     {
770 : couet 13304 // paint all the graphs of this multigraph
771 : brun 754
772 : couet 33513 if (!fGraphs) return;
773 : couet 13304 if (fGraphs->GetSize() == 0) return;
774 : rdm 11750
775 : couet 13304 char *l;
776 :     static char chopt[33];
777 :     Int_t nch = strlen(option);
778 :     Int_t i;
779 :     for (i=0;i<nch;i++) chopt[i] = toupper(option[i]);
780 :     chopt[nch] = 0;
781 :     TGraph *g;
782 : rdm 3742
783 : couet 13304 l = strstr(chopt,"A");
784 :     if (l) {
785 :     *l = ' ';
786 :     TIter next(fGraphs);
787 :     Int_t npt = 100;
788 :     Double_t maximum, minimum, rwxmin, rwxmax, rwymin, rwymax, uxmin, uxmax, dx, dy;
789 :     rwxmin = gPad->GetUxmin();
790 :     rwxmax = gPad->GetUxmax();
791 :     rwymin = gPad->GetUymin();
792 :     rwymax = gPad->GetUymax();
793 :     char *xtitle = 0;
794 :     char *ytitle = 0;
795 :     Int_t firstx = 0;
796 :     Int_t lastx = 0;
797 : rdm 11750
798 : couet 13304 if (fHistogram) {
799 :     //cleanup in case of a previous unzoom
800 :     if (fHistogram->GetMinimum() >= fHistogram->GetMaximum()) {
801 : brun 23506 nch = strlen(fHistogram->GetXaxis()->GetTitle());
802 : couet 13304 firstx = fHistogram->GetXaxis()->GetFirst();
803 :     lastx = fHistogram->GetXaxis()->GetLast();
804 :     if (nch) {
805 :     xtitle = new char[nch+1];
806 : brun 35498 strlcpy(xtitle,fHistogram->GetXaxis()->GetTitle(),nch+1);
807 : couet 13304 }
808 :     nch = strlen(fHistogram->GetYaxis()->GetTitle());
809 :     if (nch) {
810 :     ytitle = new char[nch+1];
811 : brun 35498 strlcpy(ytitle,fHistogram->GetYaxis()->GetTitle(),nch+1);
812 : couet 13304 }
813 :     delete fHistogram;
814 :     fHistogram = 0;
815 :     }
816 :     }
817 :     if (fHistogram) {
818 :     minimum = fHistogram->GetYaxis()->GetXmin();
819 :     maximum = fHistogram->GetYaxis()->GetXmax();
820 :     uxmin = gPad->PadtoX(rwxmin);
821 :     uxmax = gPad->PadtoX(rwxmax);
822 :     } else {
823 : couet 21573 g = (TGraph*) next();
824 : brun 32499 if (g) g->ComputeRange(rwxmin, rwymin, rwxmax, rwymax);
825 : couet 13304 while ((g = (TGraph*) next())) {
826 : couet 20433 Double_t rx1,ry1,rx2,ry2;
827 :     g->ComputeRange(rx1, ry1, rx2, ry2);
828 :     if (rx1 < rwxmin) rwxmin = rx1;
829 :     if (ry1 < rwymin) rwymin = ry1;
830 :     if (rx2 > rwxmax) rwxmax = rx2;
831 :     if (ry2 > rwymax) rwymax = ry2;
832 : couet 13304 if (g->GetN() > npt) npt = g->GetN();
833 :     }
834 :     if (rwxmin == rwxmax) rwxmax += 1.;
835 :     if (rwymin == rwymax) rwymax += 1.;
836 :     dx = 0.05*(rwxmax-rwxmin);
837 :     dy = 0.05*(rwymax-rwymin);
838 :     uxmin = rwxmin - dx;
839 :     uxmax = rwxmax + dx;
840 :     if (gPad->GetLogy()) {
841 :     if (rwymin <= 0) rwymin = 0.001*rwymax;
842 :     minimum = rwymin/(1+0.5*TMath::Log10(rwymax/rwymin));
843 :     maximum = rwymax*(1+0.2*TMath::Log10(rwymax/rwymin));
844 :     } else {
845 :     minimum = rwymin - dy;
846 :     maximum = rwymax + dy;
847 :     }
848 :     if (minimum < 0 && rwymin >= 0) minimum = 0;
849 :     if (maximum > 0 && rwymax <= 0) maximum = 0;
850 :     }
851 : brun 754
852 : couet 13304 if (fMinimum != -1111) rwymin = minimum = fMinimum;
853 :     if (fMaximum != -1111) rwymax = maximum = fMaximum;
854 :     if (uxmin < 0 && rwxmin >= 0) {
855 :     if (gPad->GetLogx()) uxmin = 0.9*rwxmin;
856 :     //else uxmin = 0;
857 :     }
858 :     if (uxmax > 0 && rwxmax <= 0) {
859 :     if (gPad->GetLogx()) uxmax = 1.1*rwxmax;
860 :     //else uxmax = 0;
861 :     }
862 :     if (minimum < 0 && rwymin >= 0) {
863 :     if(gPad->GetLogy()) minimum = 0.9*rwymin;
864 :     //else minimum = 0;
865 :     }
866 :     if (maximum > 0 && rwymax <= 0) {
867 :     if(gPad->GetLogy()) maximum = 1.1*rwymax;
868 :     //else maximum = 0;
869 :     }
870 :     if (minimum <= 0 && gPad->GetLogy()) minimum = 0.001*maximum;
871 :     if (uxmin <= 0 && gPad->GetLogx()) {
872 :     if (uxmax > 1000) uxmin = 1;
873 :     else uxmin = 0.001*uxmax;
874 :     }
875 :     rwymin = minimum;
876 :     rwymax = maximum;
877 :     if (fHistogram) {
878 :     fHistogram->GetYaxis()->SetLimits(rwymin,rwymax);
879 :     }
880 : brun 3575
881 : couet 13304 // Create a temporary histogram to draw the axis
882 :     if (!fHistogram) {
883 :     // the graph is created with at least as many channels as there are points
884 :     // to permit zooming on the full range
885 :     rwxmin = uxmin;
886 :     rwxmax = uxmax;
887 :     fHistogram = new TH1F(GetName(),GetTitle(),npt,rwxmin,rwxmax);
888 :     if (!fHistogram) return;
889 :     fHistogram->SetMinimum(rwymin);
890 :     fHistogram->SetBit(TH1::kNoStats);
891 :     fHistogram->SetMaximum(rwymax);
892 :     fHistogram->GetYaxis()->SetLimits(rwymin,rwymax);
893 :     fHistogram->SetDirectory(0);
894 :     if (xtitle) {fHistogram->GetXaxis()->SetTitle(xtitle); delete [] xtitle;}
895 :     if (ytitle) {fHistogram->GetYaxis()->SetTitle(ytitle); delete [] ytitle;}
896 :     if (firstx != lastx) fHistogram->GetXaxis()->SetRange(firstx,lastx);
897 :     }
898 :     fHistogram->Paint("0");
899 : brun 3575 }
900 : rdm 3742
901 : couet 21962 TGraph *gfit = 0;
902 : brun 754 if (fGraphs) {
903 : brun 4037 TObjOptLink *lnk = (TObjOptLink*)fGraphs->FirstLink();
904 : couet 21962 TObject *obj = 0;
905 : brun 4037
906 :     while (lnk) {
907 :     obj = lnk->GetObject();
908 :     if (strlen(lnk->GetOption())) obj->Paint(lnk->GetOption());
909 :     else obj->Paint(chopt);
910 :     lnk = (TObjOptLink*)lnk->Next();
911 :     }
912 : couet 21962 gfit = (TGraph*)obj; // pick one TGraph in the list to paint the fit parameters.
913 : brun 754 }
914 : brun 11226
915 :     TObject *f;
916 : couet 21962 TF1 *fit = 0;
917 : brun 11226 if (fFunctions) {
918 : couet 13304 TIter next(fFunctions);
919 :     while ((f = (TObject*) next())) {
920 :     if (f->InheritsFrom(TF1::Class())) {
921 :     if (f->TestBit(TF1::kNotDraw) == 0) f->Paint("lsame");
922 : couet 21962 fit = (TF1*)f;
923 : couet 13304 } else {
924 :     f->Paint();
925 :     }
926 : brun 11226 }
927 :     }
928 : couet 21962
929 : brun 24084 if (fit) gfit->PaintStats(fit);
930 : couet 13304 }
931 : brun 11226
932 :    
933 : brun 754 //______________________________________________________________________________
934 : brun 1205 void TMultiGraph::Print(Option_t *option) const
935 : brun 754 {
936 : couet 13304 // Print the list of graphs
937 : brun 754
938 :     TGraph *g;
939 :     if (fGraphs) {
940 : couet 13304 TIter next(fGraphs);
941 :     while ((g = (TGraph*) next())) {
942 :     g->Print(option);
943 :     }
944 : brun 754 }
945 :     }
946 :    
947 : couet 13304
948 : brun 754 //______________________________________________________________________________
949 : brun 5552 void TMultiGraph::RecursiveRemove(TObject *obj)
950 :     {
951 :     // Recursively remove this object from a list. Typically implemented
952 :     // by classes that can contain mulitple references to a same object.
953 :    
954 :     if (!fGraphs) return;
955 :     TObject *objr = fGraphs->Remove(obj);
956 :     if (!objr) return;
957 :     delete fHistogram; fHistogram = 0;
958 :     if (gPad) gPad->Modified();
959 :     }
960 :    
961 : couet 13304
962 : brun 5552 //______________________________________________________________________________
963 : brun 15672 void TMultiGraph::SavePrimitive(ostream &out, Option_t *option /*= ""*/)
964 : brun 754 {
965 : couet 13304 // Save primitive as a C++ statement(s) on output stream out
966 : brun 754
967 :     char quote = '"';
968 :     out<<" "<<endl;
969 :     if (gROOT->ClassSaved(TMultiGraph::Class())) {
970 : couet 13304 out<<" ";
971 : brun 754 } else {
972 : couet 13304 out<<" TMultiGraph *";
973 : brun 754 }
974 :     out<<"multigraph = new TMultiGraph();"<<endl;
975 :     out<<" multigraph->SetName("<<quote<<GetName()<<quote<<");"<<endl;
976 :     out<<" multigraph->SetTitle("<<quote<<GetTitle()<<quote<<");"<<endl;
977 :    
978 :     if (fGraphs) {
979 : brun 10847 TObjOptLink *lnk = (TObjOptLink*)fGraphs->FirstLink();
980 :     TObject *g;
981 :    
982 :     while (lnk) {
983 :     g = lnk->GetObject();
984 : couet 20297 g->SavePrimitive(out, Form("multigraph%s",lnk->GetOption()));
985 : brun 10847 lnk = (TObjOptLink*)lnk->Next();
986 : couet 13304 }
987 : brun 754 }
988 : couet 20297 out<<" multigraph->Draw(" <<quote<<option<<quote<<");"<<endl;
989 : couet 14017
990 :     TAxis *xaxis = GetXaxis();
991 :     TAxis *yaxis = GetYaxis();
992 : couet 20433
993 :     if (xaxis) xaxis->SaveAttributes(out, "multigraph","->GetXaxis()");
994 :     if (yaxis) yaxis->SaveAttributes(out, "multigraph","->GetYaxis()");
995 : brun 754 }
996 :    
997 : couet 13304
998 : brun 754 //______________________________________________________________________________
999 :     void TMultiGraph::SetMaximum(Double_t maximum)
1000 :     {
1001 : couet 13304 // Set multigraph maximum.
1002 :    
1003 : brun 754 fMaximum = maximum;
1004 :     if (fHistogram) fHistogram->SetMaximum(maximum);
1005 :     }
1006 :    
1007 : couet 13304
1008 : brun 754 //______________________________________________________________________________
1009 :     void TMultiGraph::SetMinimum(Double_t minimum)
1010 :     {
1011 : couet 13304 // Set multigraph minimum.
1012 :    
1013 : brun 754 fMinimum = minimum;
1014 :     if (fHistogram) fHistogram->SetMinimum(minimum);
1015 :     }

Subversion Admin
ViewVC Help
Powered by ViewVC 1.0.9