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

Subversion Admin
ViewVC Help
Powered by ViewVC 1.0.9