TMultiGraph.cxx

Go to the documentation of this file.
00001 // @(#)root/hist:$Id: TMultiGraph.cxx 36530 2010-11-08 10:24:42Z couet $
00002 // Author: Rene Brun   12/10/2000
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #include "TROOT.h"
00013 #include "TMultiGraph.h"
00014 #include "TGraph.h"
00015 #include "TH1.h"
00016 #include "TVirtualPad.h"
00017 #include "Riostream.h"
00018 #include "TVirtualFitter.h"
00019 #include "TPluginManager.h"
00020 #include "TClass.h"
00021 #include "TMath.h"
00022 #include "TSystem.h"
00023 #include <stdlib.h>
00024 
00025 #include "HFitInterface.h"
00026 #include "Fit/DataRange.h"
00027 #include "Math/MinimizerOptions.h"
00028 
00029 #include <ctype.h>
00030 
00031 extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
00032 
00033 ClassImp(TMultiGraph)
00034 
00035 
00036 //______________________________________________________________________________
00037 /* Begin_Html
00038 <center><h2>TMultiGraph class</h2></center>
00039 A TMultiGraph is a collection of TGraph (or derived) objects
00040 Use <tt>TMultiGraph::Add</tt> to add a new graph to the list.
00041 The TMultiGraph owns the objects in the list.
00042 Drawing options are the same as for TGraph.
00043 <p>
00044 Example:
00045 <pre>
00046      TGraph *gr1 = new TGraph(...
00047      TGraphErrors *gr2 = new TGraphErrors(...
00048      TMultiGraph *mg = new TMultiGraph();
00049      mg->Add(gr1,"lp");
00050      mg->Add(gr2,"cp");
00051      mg->Draw("a");
00052 </pre>
00053 The drawing option for each TGraph may be specified as an optional
00054 second argument of the Add function.
00055 If a draw option is specified, it will be used to draw the graph,
00056 otherwise the graph will be drawn with the option specified in
00057 <tt>TMultiGraph::Draw</tt>.
00058 End_Html */
00059 
00060 
00061 //______________________________________________________________________________
00062 TMultiGraph::TMultiGraph(): TNamed()
00063 {
00064    // TMultiGraph default constructor
00065 
00066    fGraphs    = 0;
00067    fFunctions = 0;
00068    fHistogram = 0;
00069    fMaximum   = -1111;
00070    fMinimum   = -1111;
00071 }
00072 
00073 
00074 //______________________________________________________________________________
00075 TMultiGraph::TMultiGraph(const char *name, const char *title)
00076        : TNamed(name,title)
00077 {
00078    // constructor with name and title
00079 
00080    fGraphs    = 0;
00081    fFunctions = 0;
00082    fHistogram = 0;
00083    fMaximum   = -1111;
00084    fMinimum   = -1111;
00085 }
00086 
00087 
00088 //______________________________________________________________________________
00089 TMultiGraph::TMultiGraph(const TMultiGraph& mg) :
00090   TNamed (mg),
00091   fGraphs(mg.fGraphs),
00092   fFunctions(mg.fFunctions),
00093   fHistogram(mg.fHistogram),
00094   fMaximum(mg.fMaximum),
00095   fMinimum(mg.fMinimum)
00096 {
00097    //copy constructor
00098 }
00099 
00100 
00101 //______________________________________________________________________________
00102 TMultiGraph& TMultiGraph::operator=(const TMultiGraph& mg)
00103 {
00104    //assignement operator
00105    if(this!=&mg) {
00106       TNamed::operator=(mg);
00107       fGraphs=mg.fGraphs;
00108       fFunctions=mg.fFunctions;
00109       fHistogram=mg.fHistogram;
00110       fMaximum=mg.fMaximum;
00111       fMinimum=mg.fMinimum;
00112    }
00113    return *this;
00114 }
00115 
00116 
00117 //______________________________________________________________________________
00118 TMultiGraph::~TMultiGraph()
00119 {
00120    // TMultiGraph destructor
00121 
00122    if (!fGraphs) return;
00123    TGraph *g;
00124    TIter   next(fGraphs);
00125    while ((g = (TGraph*) next())) {
00126       g->ResetBit(kMustCleanup);
00127    }
00128    fGraphs->Delete();
00129    delete fGraphs;
00130    fGraphs = 0;
00131    delete fHistogram;
00132    fHistogram = 0;
00133    if (fFunctions) {
00134       fFunctions->SetBit(kInvalidObject);
00135       //special logic to support the case where the same object is
00136       //added multiple times in fFunctions.
00137       //This case happens when the same object is added with different
00138       //drawing modes
00139       TObject *obj;
00140       while ((obj  = fFunctions->First())) {
00141          while(fFunctions->Remove(obj)) { }
00142          delete obj;
00143       }
00144       delete fFunctions;
00145    }
00146 }
00147 
00148 
00149 //______________________________________________________________________________
00150 void TMultiGraph::Add(TGraph *graph, Option_t *chopt)
00151 {
00152    // add a new graph to the list of graphs
00153    // note that the graph is now owned by the TMultigraph.
00154    // Deleting the TMultiGraph object will automatically delete the graphs.
00155    // You should not delete the graphs when the TMultigraph is still active.
00156 
00157    if (!fGraphs) fGraphs = new TList();
00158    graph->SetBit(kMustCleanup);
00159    fGraphs->Add(graph,chopt);
00160 }
00161 
00162 
00163 //______________________________________________________________________________
00164 void TMultiGraph::Add(TMultiGraph *multigraph, Option_t *chopt)
00165 {
00166    // add all the graphs in "multigraph" to the list of graphs.
00167 
00168    TList *graphlist = multigraph->GetListOfGraphs();
00169    if (!graphlist) return;
00170 
00171    if (!fGraphs) fGraphs = new TList();
00172 
00173    TGraph *gr;
00174    gr = (TGraph*)graphlist->First();
00175    fGraphs->Add(gr,chopt);
00176    for(Int_t i = 1; i < graphlist->GetSize(); i++){
00177       gr = (TGraph*)graphlist->After(gr);
00178       fGraphs->Add(gr,chopt);
00179    }
00180 }
00181 
00182 
00183 //______________________________________________________________________________
00184 void TMultiGraph::Browse(TBrowser *)
00185 {
00186    // Browse multigraph.
00187 
00188    Draw("alp");
00189    gPad->Update();
00190 }
00191 
00192 
00193 //______________________________________________________________________________
00194 Int_t TMultiGraph::DistancetoPrimitive(Int_t px, Int_t py)
00195 {
00196    // Compute distance from point px,py to each graph
00197 
00198    // Are we on the axis?
00199    const Int_t kMaxDiff = 10;
00200    Int_t distance = 9999;
00201    if (fHistogram) {
00202       distance = fHistogram->DistancetoPrimitive(px,py);
00203       if (distance <= 0) return distance;
00204    }
00205 
00206    // Loop on the list of graphs
00207    if (!fGraphs) return distance;
00208    TGraph *g;
00209    TIter   next(fGraphs);
00210    while ((g = (TGraph*) next())) {
00211       Int_t dist = g->DistancetoPrimitive(px,py);
00212       if (dist <= 0) return 0;
00213       if (dist < kMaxDiff) {gPad->SetSelected(g); return dist;}
00214    }
00215    return distance;
00216 }
00217 
00218 
00219 //______________________________________________________________________________
00220 void TMultiGraph::Draw(Option_t *option)
00221 {
00222    // Draw this multigraph with its current attributes.
00223    //
00224    //   Options to draw a graph are described in TGraph::PainGraph
00225    //
00226    //  The drawing option for each TGraph may be specified as an optional
00227    //  second argument of the Add function. You can use GetGraphDrawOption
00228    //  to return this option.
00229    //  If a draw option is specified, it will be used to draw the graph,
00230    //  otherwise the graph will be drawn with the option specified in
00231    //  TMultiGraph::Draw. Use GetDrawOption to return the option specified
00232    //  when drawin the TMultiGraph.
00233 
00234    AppendPad(option);
00235 }
00236 
00237 
00238 //______________________________________________________________________________
00239 TFitResultPtr TMultiGraph::Fit(const char *fname, Option_t *option, Option_t *, Axis_t xmin, Axis_t xmax)
00240 {
00241    // Fit this graph with function with name fname.
00242    //
00243    //  interface to TF1::Fit(TF1 *f1...
00244 
00245    char *linear;
00246    linear= (char*)strstr(fname, "++");
00247    TF1 *f1=0;
00248    if (linear)
00249       f1=new TF1(fname, fname, xmin, xmax);
00250    else {
00251       f1 = (TF1*)gROOT->GetFunction(fname);
00252       if (!f1) { Printf("Unknown function: %s",fname); return -1; }
00253    }
00254 
00255    return Fit(f1,option,"",xmin,xmax);
00256 }
00257 
00258 
00259 //______________________________________________________________________________
00260 TFitResultPtr TMultiGraph::Fit(TF1 *f1, Option_t *option, Option_t *goption, Axis_t rxmin, Axis_t rxmax)
00261 {
00262    // Fit this multigraph with function f1.
00263    //
00264    //   In this function all graphs of the multigraph are fitted simultaneously
00265    //
00266    //   f1 is an already predefined function created by TF1.
00267    //   Predefined functions such as gaus, expo and poln are automatically
00268    //   created by ROOT.
00269    //
00270    //   The list of fit options is given in parameter option.
00271    //      option = "W"  Set all errors to 1
00272    //             = "U" Use a User specified fitting algorithm (via SetFCN)
00273    //             = "Q" Quiet mode (minimum printing)
00274    //             = "V" Verbose mode (default is between Q and V)
00275    //             = "B" Use this option when you want to fix one or more parameters
00276    //                   and the fitting function is like "gaus","expo","poln","landau".
00277    //             = "R" Use the Range specified in the function range
00278    //             = "N" Do not store the graphics function, do not draw
00279    //             = "0" Do not plot the result of the fit. By default the fitted function
00280    //                   is drawn unless the option"N" above is specified.
00281    //             = "+" Add this new fitted function to the list of fitted functions
00282    //                   (by default, any previous function is deleted)
00283    //             = "C" In case of linear fitting, not calculate the chisquare
00284    //                    (saves time)
00285    //             = "F" If fitting a polN, switch to minuit fitter
00286    //             = "ROB" In case of linear fitting, compute the LTS regression
00287    //                     coefficients (robust(resistant) regression), using
00288    //                     the default fraction of good points
00289    //               "ROB=0.x" - compute the LTS regression coefficients, using
00290    //                           0.x as a fraction of good points
00291    //
00292    //   When the fit is drawn (by default), the parameter goption may be used
00293    //   to specify a list of graphics options. See TGraph::Paint for a complete
00294    //   list of these options.
00295    //
00296    //   In order to use the Range option, one must first create a function
00297    //   with the expression to be fitted. For example, if your graph
00298    //   has a defined range between -4 and 4 and you want to fit a gaussian
00299    //   only in the interval 1 to 3, you can do:
00300    //        TF1 *f1 = new TF1("f1","gaus",1,3);
00301    //        graph->Fit("f1","R");
00302    //
00303    //
00304    //   who is calling this function
00305    //   ============================
00306    //   Note that this function is called when calling TGraphErrors::Fit
00307    //   or TGraphAsymmErrors::Fit ot TGraphBentErrors::Fit
00308    //   see the discussion below on the errors calulation.
00309    //
00310    //   Setting initial conditions
00311    //   ==========================
00312    //   Parameters must be initialized before invoking the Fit function.
00313    //   The setting of the parameter initial values is automatic for the
00314    //   predefined functions : poln, expo, gaus, landau. One can however disable
00315    //   this automatic computation by specifying the option "B".
00316    //   You can specify boundary limits for some or all parameters via
00317    //        f1->SetParLimits(p_number, parmin, parmax);
00318    //   if parmin>=parmax, the parameter is fixed
00319    //   Note that you are not forced to fix the limits for all parameters.
00320    //   For example, if you fit a function with 6 parameters, you can do:
00321    //     func->SetParameters(0,3.1,1.e-6,0.1,-8,100);
00322    //     func->SetParLimits(4,-10,-4);
00323    //     func->SetParLimits(5, 1,1);
00324    //   With this setup, parameters 0->3 can vary freely
00325    //   Parameter 4 has boundaries [-10,-4] with initial value -8
00326    //   Parameter 5 is fixed to 100.
00327    //
00328    //  Fit range
00329    //  =========
00330    //  The fit range can be specified in two ways:
00331    //    - specify rxmax > rxmin (default is rxmin=rxmax=0)
00332    //    - specify the option "R". In this case, the function will be taken
00333    //      instead of the full graph range.
00334    //
00335    //  Changing the fitting function
00336    //  =============================
00337    //   By default a chi2 fitting function is used for fitting the TGraphs's.
00338    //   The function is implemented in FitUtil::EvaluateChi2.
00339    //   In case of TGraphErrors an effective chi2 is used 
00340    //   (see TGraphErrors fit in TGraph::Fit) and is implemented in 
00341    //   FitUtil::EvaluateChi2Effective
00342    //   To specify a User defined fitting function, specify option "U" and
00343    //   call the following functions:
00344    //     TVirtualFitter::Fitter(mygraph)->SetFCN(MyFittingFunction)
00345    //   where MyFittingFunction is of type:
00346    //   extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
00347    //
00348    //  Access to the fit result
00349    //  ========================
00350    //  The function returns a TFitResultPtr which can hold a  pointer to a TFitResult object.
00351    //  By default the TFitResultPtr contains only the status of the fit and it converts
00352    //  automatically to an integer. If the option "S" is instead used, TFitResultPtr contains
00353    //  the TFitResult and behaves as a smart pointer to it. For example one can do:
00354    //     TFitResultPtr r = graph->Fit("myFunc","S");
00355    //     TMatrixDSym cov = r->GetCovarianceMatrix();  //  to access the covariance matrix
00356    //     Double_t par0   = r->Parameter(0); // retrieve the value for the parameter 0
00357    //     Double_t err0   = r->ParError(0); // retrieve the error for the parameter 0
00358    //     r->Print("V");     // print full information of fit including covariance matrix
00359    //     r->Write();        // store the result in a file
00360    //
00361    //   The fit parameters, error and chi2 (but not covariance matrix) can be retrieved also
00362    //   from the fitted function.
00363    //
00364    //
00365    //   Associated functions
00366    //   ====================
00367    //  One or more object (typically a TF1*) can be added to the list
00368    //  of functions (fFunctions) associated to each graph.
00369    //  When TGraph::Fit is invoked, the fitted function is added to this list.
00370    //  Given a graph gr, one can retrieve an associated function
00371    //  with:  TF1 *myfunc = gr->GetFunction("myfunc");
00372    //
00373    //  If the graph is made persistent, the list of
00374    //  associated functions is also persistent. Given a pointer (see above)
00375    //  to an associated function myfunc, one can retrieve the function/fit
00376    //  parameters with calls such as:
00377    //    Double_t chi2 = myfunc->GetChisquare();
00378    //    Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
00379    //    Double_t err0 = myfunc->GetParError(0);  //error on first parameter
00380    //
00381    //   Fit Statistics
00382    //   ==============
00383    //  You can change the statistics box to display the fit parameters with
00384    //  the TStyle::SetOptFit(mode) method. This mode has four digits.
00385    //  mode = pcev  (default = 0111)
00386    //    v = 1;  print name/values of parameters
00387    //    e = 1;  print errors (if e=1, v must be 1)
00388    //    c = 1;  print Chisquare/Number of degress of freedom
00389    //    p = 1;  print Probability
00390    //
00391    //  For example: gStyle->SetOptFit(1011);
00392    //  prints the fit probability, parameter names/values, and errors.
00393    //  You can change the position of the statistics box with these lines
00394    //  (where g is a pointer to the TGraph):
00395    //
00396    //  Root > TPaveStats *st = (TPaveStats*)g->GetListOfFunctions()->FindObject("stats")
00397    //  Root > st->SetX1NDC(newx1); //new x start position
00398    //  Root > st->SetX2NDC(newx2); //new x end position
00399 
00400    // internal multigraph fitting methods
00401    Foption_t fitOption;
00402    ROOT::Fit::FitOptionsMake(option,fitOption);
00403 
00404    // create range and minimizer options with default values 
00405    ROOT::Fit::DataRange range(rxmin,rxmax); 
00406    ROOT::Math::MinimizerOptions minOption; 
00407    return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range); 
00408 
00409 }
00410 
00411 //______________________________________________________________________________
00412 void TMultiGraph::FitPanel()
00413 {
00414 //   -*-*-*-*-*Display a panel with all histogram fit options*-*-*-*-*-*
00415 //             ==============================================
00416 //
00417 //      See class TFitPanel for example
00418 
00419    if (!gPad)
00420       gROOT->MakeDefCanvas();
00421 
00422    if (!gPad) {
00423       Error("FitPanel", "Unable to create a default canvas");
00424       return;
00425    }
00426 
00427    // use plugin manager to create instance of TFitEditor
00428    TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
00429    if (handler && handler->LoadPlugin() != -1) {
00430       if (handler->ExecPlugin(2, gPad, this) == 0)
00431          Error("FitPanel", "Unable to crate the FitPanel");
00432    }
00433    else
00434          Error("FitPanel", "Unable to find the FitPanel plug-in");
00435 }
00436 
00437 //______________________________________________________________________________
00438 Option_t *TMultiGraph::GetGraphDrawOption(const TGraph *gr) const
00439 {
00440    // Return the draw option for the TGraph gr in this TMultiGraph
00441    // The return option is the one specified when calling TMultiGraph::Add(gr,option).
00442 
00443    if (!fGraphs || !gr) return "";
00444    TListIter next(fGraphs);
00445    TObject *obj;
00446    while ((obj = next())) {
00447       if (obj == (TObject*)gr) return next.GetOption();
00448    }
00449    return "";
00450 }
00451 
00452 
00453 //______________________________________________________________________________
00454 void TMultiGraph::InitGaus(Double_t xmin, Double_t xmax)
00455 {
00456    // Compute Initial values of parameters for a gaussian.
00457 
00458    Double_t allcha, sumx, sumx2, x, val, rms, mean;
00459    Int_t bin;
00460    const Double_t sqrtpi = 2.506628;
00461 
00462    // Compute mean value and RMS of the graph in the given range
00463    Int_t np = 0;
00464    allcha = sumx = sumx2 = 0;
00465    TGraph *g;
00466    TIter next(fGraphs);
00467    Double_t *px, *py;
00468    Int_t npp; //number of points in each graph
00469    while ((g = (TGraph*) next())) {
00470       px=g->GetX();
00471       py=g->GetY();
00472       npp=g->GetN();
00473       for (bin=0; bin<npp; bin++){
00474          x=px[bin];
00475          if (x<xmin || x>xmax) continue;
00476          np++;
00477          val=py[bin];
00478          sumx+=val*x;
00479          sumx2+=val*x*x;
00480          allcha+=val;
00481       }
00482    }
00483    if (np == 0 || allcha == 0) return;
00484    mean = sumx/allcha;
00485    rms  = TMath::Sqrt(sumx2/allcha - mean*mean);
00486 
00487    Double_t binwidx = TMath::Abs((xmax-xmin)/np);
00488    if (rms == 0) rms = 1;
00489    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
00490    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
00491    f1->SetParameter(0,binwidx*allcha/(sqrtpi*rms));
00492    f1->SetParameter(1,mean);
00493    f1->SetParameter(2,rms);
00494    f1->SetParLimits(2,0,10*rms);
00495 }
00496 
00497 
00498 //______________________________________________________________________________
00499 void TMultiGraph::InitExpo(Double_t xmin, Double_t xmax)
00500 {
00501    // Compute Initial values of parameters for an exponential.
00502 
00503    Double_t constant, slope;
00504    Int_t ifail;
00505 
00506    LeastSquareLinearFit(-1, constant, slope, ifail, xmin, xmax);
00507 
00508    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
00509    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
00510    f1->SetParameter(0,constant);
00511    f1->SetParameter(1,slope);
00512 }
00513 
00514 
00515 //______________________________________________________________________________
00516 void TMultiGraph::InitPolynom(Double_t xmin, Double_t xmax)
00517 {
00518    // Compute Initial values of parameters for a polynom.
00519 
00520    Double_t fitpar[25];
00521 
00522    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
00523    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
00524    Int_t npar   = f1->GetNpar();
00525 
00526    LeastSquareFit(npar, fitpar, xmin, xmax);
00527 
00528    for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
00529 }
00530 
00531 
00532 //______________________________________________________________________________
00533 void TMultiGraph::LeastSquareFit(Int_t m, Double_t *a, Double_t xmin, Double_t xmax)
00534 {
00535    // Least squares lpolynomial fitting without weights.
00536    //
00537    //  m     number of parameters
00538    //  a     array of parameters
00539    //  first 1st point number to fit (default =0)
00540    //  last  last point number to fit (default=fNpoints-1)
00541    //
00542    //   based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
00543 
00544    const Double_t zero = 0.;
00545    const Double_t one = 1.;
00546    const Int_t idim = 20;
00547 
00548    Double_t  b[400]        /* was [20][20] */;
00549    Int_t i, k, l, ifail, bin;
00550    Double_t power;
00551    Double_t da[20], xk, yk;
00552 
00553 
00554    //count the total number of points to fit
00555    TGraph *g;
00556    TIter next(fGraphs);
00557    Double_t *px, *py;
00558    Int_t n=0;
00559    Int_t npp;
00560    while ((g = (TGraph*) next())) {
00561       px=g->GetX();
00562       npp=g->GetN();
00563       for (bin=0; bin<npp; bin++){
00564          xk=px[bin];
00565          if (xk < xmin || xk > xmax) continue;
00566          n++;
00567       }
00568    }
00569    if (m <= 2) {
00570       LeastSquareLinearFit(n, a[0], a[1], ifail, xmin, xmax);
00571       return;
00572    }
00573    if (m > idim || m > n) return;
00574    da[0] = zero;
00575    for (l = 2; l <= m; ++l) {
00576       b[l-1]           = zero;
00577       b[m + l*20 - 21] = zero;
00578       da[l-1]          = zero;
00579    }
00580    Int_t np = 0;
00581 
00582    next.Reset();
00583    while ((g = (TGraph*) next())) {
00584       px=g->GetX();
00585       py=g->GetY();
00586       npp=g->GetN();
00587 
00588       for (k = 0; k <= npp; ++k) {
00589          xk     = px[k];
00590          if (xk < xmin || xk > xmax) continue;
00591          np++;
00592          yk     = py[k];
00593          power  = one;
00594          da[0] += yk;
00595          for (l = 2; l <= m; ++l) {
00596             power   *= xk;
00597             b[l-1]  += power;
00598             da[l-1] += power*yk;
00599          }
00600          for (l = 2; l <= m; ++l) {
00601             power            *= xk;
00602             b[m + l*20 - 21] += power;
00603          }
00604       }
00605    }
00606    b[0]  = Double_t(np);
00607    for (i = 3; i <= m; ++i) {
00608       for (k = i; k <= m; ++k) {
00609          b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
00610       }
00611    }
00612    H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
00613 
00614    if (ifail < 0) {
00615       //a[0] = fY[0];
00616       py=((TGraph *)fGraphs->First())->GetY();
00617       a[0]=py[0];
00618       for (i=1; i<m; ++i) a[i] = 0;
00619       return;
00620    }
00621    for (i=0; i<m; ++i) a[i] = da[i];
00622 }
00623 
00624 
00625 //______________________________________________________________________________
00626 void TMultiGraph::LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail, Double_t xmin, Double_t xmax)
00627 {
00628    // Least square linear fit without weights.
00629    //
00630    //  Fit a straight line (a0 + a1*x) to the data in this graph.
00631    //  ndata:  number of points to fit
00632    //  first:  first point number to fit
00633    //  last:   last point to fit O(ndata should be last-first
00634    //  ifail:  return parameter indicating the status of the fit (ifail=0, fit is OK)
00635    //
00636    //   extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
00637 
00638    Double_t xbar, ybar, x2bar;
00639    Int_t i;
00640    Double_t xybar;
00641    Double_t fn, xk, yk;
00642    Double_t det;
00643 
00644    ifail = -2;
00645    xbar  = ybar = x2bar = xybar = 0;
00646    Int_t np = 0;
00647    TGraph *g;
00648    TIter next(fGraphs);
00649    Double_t *px, *py;
00650    Int_t npp;
00651    while ((g = (TGraph*) next())) {
00652       px=g->GetX();
00653       py=g->GetY();
00654       npp=g->GetN();
00655       for (i = 0; i < npp; ++i) {
00656          xk = px[i];
00657          if (xk < xmin || xk > xmax) continue;
00658          np++;
00659          yk = py[i];
00660          if (ndata < 0) {
00661             if (yk <= 0) yk = 1e-9;
00662             yk = TMath::Log(yk);
00663          }
00664          xbar  += xk;
00665          ybar  += yk;
00666          x2bar += xk*xk;
00667          xybar += xk*yk;
00668       }
00669    }
00670    fn    = Double_t(np);
00671    det   = fn*x2bar - xbar*xbar;
00672    ifail = -1;
00673    if (det <= 0) {
00674       if (fn > 0) a0 = ybar/fn;
00675       else        a0 = 0;
00676       a1 = 0;
00677       return;
00678    }
00679    ifail = 0;
00680    a0 = (x2bar*ybar - xbar*xybar) / det;
00681    a1 = (fn*xybar - xbar*ybar) / det;
00682 }
00683 
00684 
00685 //______________________________________________________________________________
00686 Int_t TMultiGraph::IsInside(Double_t x, Double_t y) const
00687 {
00688    // Return 1 if the point (x,y) is inside one of the graphs 0 otherwise.
00689 
00690    Int_t in = 0;
00691    if (!fGraphs) return in;
00692    TGraph *g;
00693    TIter next(fGraphs);
00694    while ((g = (TGraph*) next())) {
00695       in = g->IsInside(x, y);
00696       if (in) return in;
00697    }
00698    return in;
00699 }
00700 
00701 
00702 //______________________________________________________________________________
00703 TH1F *TMultiGraph::GetHistogram() const
00704 {
00705    // Returns a pointer to the histogram used to draw the axis
00706    // Takes into account the two following cases.
00707    //    1- option 'A' was specified in TMultiGraph::Draw. Return fHistogram
00708    //    2- user had called TPad::DrawFrame. return pointer to hframe histogram
00709 
00710    if (fHistogram) return fHistogram;
00711    if (!gPad) return 0;
00712    gPad->Modified();
00713    gPad->Update();
00714    if (fHistogram) return fHistogram;
00715    TH1F *h1 = (TH1F*)gPad->FindObject("hframe");
00716    return h1;
00717 }
00718 
00719 
00720 //______________________________________________________________________________
00721 TF1 *TMultiGraph::GetFunction(const char *name) const
00722 {
00723    // Return pointer to function with name.
00724    //
00725    // Functions such as TGraph::Fit store the fitted function in the list of
00726    // functions of this graph.
00727 
00728    if (!fFunctions) return 0;
00729    return (TF1*)fFunctions->FindObject(name);
00730 }
00731 
00732 //______________________________________________________________________________
00733 TList *TMultiGraph::GetListOfFunctions()
00734 {
00735    // Return pointer to list of functions
00736    // if pointer is null create the list
00737 
00738    if (!fFunctions) fFunctions = new TList();
00739    return fFunctions;
00740 }
00741 
00742 
00743 //______________________________________________________________________________
00744 TAxis *TMultiGraph::GetXaxis() const
00745 {
00746    // Get x axis of the graph.
00747 
00748    if (!gPad) return 0;
00749    TH1 *h = GetHistogram();
00750    if (!h) return 0;
00751    return h->GetXaxis();
00752 }
00753 
00754 
00755 //______________________________________________________________________________
00756 TAxis *TMultiGraph::GetYaxis() const
00757 {
00758    // Get y axis of the graph.
00759 
00760    if (!gPad) return 0;
00761    TH1 *h = GetHistogram();
00762    if (!h) return 0;
00763    return h->GetYaxis();
00764 }
00765 
00766 
00767 //______________________________________________________________________________
00768 void TMultiGraph::Paint(Option_t *option)
00769 {
00770    // paint all the graphs of this multigraph
00771 
00772    if (!fGraphs) return;
00773    if (fGraphs->GetSize() == 0) return;
00774 
00775    char *l;
00776    static char chopt[33];
00777    Int_t nch = strlen(option);
00778    Int_t i;
00779    for (i=0;i<nch;i++) chopt[i] = toupper(option[i]);
00780    chopt[nch] = 0;
00781    TGraph *g;
00782 
00783    l = strstr(chopt,"A");
00784    if (l) {
00785       *l = ' ';
00786       TIter   next(fGraphs);
00787       Int_t npt = 100;
00788       Double_t maximum, minimum, rwxmin, rwxmax, rwymin, rwymax, uxmin, uxmax, dx, dy;
00789       rwxmin    = gPad->GetUxmin();
00790       rwxmax    = gPad->GetUxmax();
00791       rwymin    = gPad->GetUymin();
00792       rwymax    = gPad->GetUymax();
00793       char *xtitle = 0;
00794       char *ytitle = 0;
00795       Int_t firstx = 0;
00796       Int_t lastx  = 0;
00797 
00798       if (fHistogram) {
00799          //cleanup in case of a previous unzoom
00800          if (fHistogram->GetMinimum() >= fHistogram->GetMaximum()) {
00801             nch = strlen(fHistogram->GetXaxis()->GetTitle());
00802             firstx = fHistogram->GetXaxis()->GetFirst();
00803             lastx  = fHistogram->GetXaxis()->GetLast();
00804             if (nch) {
00805                xtitle = new char[nch+1];
00806                strlcpy(xtitle,fHistogram->GetXaxis()->GetTitle(),nch+1);
00807             }
00808             nch = strlen(fHistogram->GetYaxis()->GetTitle());
00809             if (nch) {
00810                ytitle = new char[nch+1];
00811                strlcpy(ytitle,fHistogram->GetYaxis()->GetTitle(),nch+1);
00812             }
00813             delete fHistogram;
00814             fHistogram = 0;
00815          }
00816       }
00817       if (fHistogram) {
00818          minimum = fHistogram->GetYaxis()->GetXmin();
00819          maximum = fHistogram->GetYaxis()->GetXmax();
00820          uxmin   = gPad->PadtoX(rwxmin);
00821          uxmax   = gPad->PadtoX(rwxmax);
00822       } else {
00823          g = (TGraph*) next();
00824          if (g) g->ComputeRange(rwxmin, rwymin, rwxmax, rwymax);
00825          while ((g = (TGraph*) next())) {
00826             Double_t rx1,ry1,rx2,ry2;
00827             g->ComputeRange(rx1, ry1, rx2, ry2);
00828             if (rx1 < rwxmin) rwxmin = rx1;
00829             if (ry1 < rwymin) rwymin = ry1;
00830             if (rx2 > rwxmax) rwxmax = rx2;
00831             if (ry2 > rwymax) rwymax = ry2;
00832             if (g->GetN() > npt) npt = g->GetN();
00833          }
00834          if (rwxmin == rwxmax) rwxmax += 1.;
00835          if (rwymin == rwymax) rwymax += 1.;
00836          dx = 0.05*(rwxmax-rwxmin);
00837          dy = 0.05*(rwymax-rwymin);
00838          uxmin    = rwxmin - dx;
00839          uxmax    = rwxmax + dx;
00840          if (gPad->GetLogy()) {
00841             if (rwymin <= 0) rwymin = 0.001*rwymax;
00842             minimum = rwymin/(1+0.5*TMath::Log10(rwymax/rwymin));
00843             maximum = rwymax*(1+0.2*TMath::Log10(rwymax/rwymin));
00844          } else {
00845             minimum  = rwymin - dy;
00846             maximum  = rwymax + dy;
00847          }
00848          if (minimum < 0 && rwymin >= 0) minimum = 0;
00849          if (maximum > 0 && rwymax <= 0) maximum = 0;
00850       }
00851 
00852       if (fMinimum != -1111) rwymin = minimum = fMinimum;
00853       if (fMaximum != -1111) rwymax = maximum = fMaximum;
00854       if (uxmin < 0 && rwxmin >= 0) {
00855          if (gPad->GetLogx()) uxmin = 0.9*rwxmin;
00856          //else                 uxmin = 0;
00857       }
00858       if (uxmax > 0 && rwxmax <= 0) {
00859          if (gPad->GetLogx()) uxmax = 1.1*rwxmax;
00860          //else                 uxmax = 0;
00861       }
00862       if (minimum < 0 && rwymin >= 0) {
00863          if(gPad->GetLogy()) minimum = 0.9*rwymin;
00864          //else                minimum = 0;
00865       }
00866       if (maximum > 0 && rwymax <= 0) {
00867          if(gPad->GetLogy()) maximum = 1.1*rwymax;
00868          //else                maximum = 0;
00869       }
00870       if (minimum <= 0 && gPad->GetLogy()) minimum = 0.001*maximum;
00871       if (uxmin <= 0 && gPad->GetLogx()) {
00872          if (uxmax > 1000) uxmin = 1;
00873          else              uxmin = 0.001*uxmax;
00874       }
00875       rwymin = minimum;
00876       rwymax = maximum;
00877       if (fHistogram) {
00878          fHistogram->GetYaxis()->SetLimits(rwymin,rwymax);
00879       }
00880 
00881       // Create a temporary histogram to draw the axis
00882       if (!fHistogram) {
00883          // the graph is created with at least as many channels as there are points
00884          // to permit zooming on the full range
00885          rwxmin = uxmin;
00886          rwxmax = uxmax;
00887          fHistogram = new TH1F(GetName(),GetTitle(),npt,rwxmin,rwxmax);
00888          if (!fHistogram) return;
00889          fHistogram->SetMinimum(rwymin);
00890          fHistogram->SetBit(TH1::kNoStats);
00891          fHistogram->SetMaximum(rwymax);
00892          fHistogram->GetYaxis()->SetLimits(rwymin,rwymax);
00893          fHistogram->SetDirectory(0);
00894          if (xtitle) {fHistogram->GetXaxis()->SetTitle(xtitle); delete [] xtitle;}
00895          if (ytitle) {fHistogram->GetYaxis()->SetTitle(ytitle); delete [] ytitle;}
00896          if (firstx != lastx) fHistogram->GetXaxis()->SetRange(firstx,lastx);
00897       }
00898       fHistogram->Paint("0");
00899    }
00900 
00901    TGraph *gfit = 0;
00902    if (fGraphs) {
00903       TObjOptLink *lnk = (TObjOptLink*)fGraphs->FirstLink();
00904       TObject *obj = 0;
00905 
00906       while (lnk) {
00907          obj = lnk->GetObject();
00908          if (strlen(lnk->GetOption())) obj->Paint(lnk->GetOption());
00909          else                          obj->Paint(chopt);
00910          lnk = (TObjOptLink*)lnk->Next();
00911       }
00912       gfit = (TGraph*)obj; // pick one TGraph in the list to paint the fit parameters.
00913    }
00914 
00915    TObject *f;
00916    TF1 *fit = 0;
00917    if (fFunctions) {
00918       TIter   next(fFunctions);
00919       while ((f = (TObject*) next())) {
00920          if (f->InheritsFrom(TF1::Class())) {
00921             if (f->TestBit(TF1::kNotDraw) == 0) f->Paint("lsame");
00922             fit = (TF1*)f;
00923          } else  {
00924             f->Paint();
00925          }
00926       }
00927    }
00928 
00929    if (fit) gfit->PaintStats(fit);
00930 }
00931 
00932 
00933 //______________________________________________________________________________
00934 void TMultiGraph::Print(Option_t *option) const
00935 {
00936    // Print the list of graphs
00937 
00938    TGraph *g;
00939    if (fGraphs) {
00940       TIter   next(fGraphs);
00941       while ((g = (TGraph*) next())) {
00942          g->Print(option);
00943       }
00944    }
00945 }
00946 
00947 
00948 //______________________________________________________________________________
00949 void TMultiGraph::RecursiveRemove(TObject *obj)
00950 {
00951    // Recursively remove this object from a list. Typically implemented
00952    // by classes that can contain mulitple references to a same object.
00953 
00954    if (!fGraphs) return;
00955    TObject *objr = fGraphs->Remove(obj);
00956    if (!objr) return;
00957    delete fHistogram; fHistogram = 0;
00958    if (gPad) gPad->Modified();
00959 }
00960 
00961 
00962 //______________________________________________________________________________
00963 void TMultiGraph::SavePrimitive(ostream &out, Option_t *option /*= ""*/)
00964 {
00965    // Save primitive as a C++ statement(s) on output stream out
00966 
00967    char quote = '"';
00968    out<<"   "<<endl;
00969    if (gROOT->ClassSaved(TMultiGraph::Class())) {
00970       out<<"   ";
00971    } else {
00972       out<<"   TMultiGraph *";
00973    }
00974    out<<"multigraph = new TMultiGraph();"<<endl;
00975    out<<"   multigraph->SetName("<<quote<<GetName()<<quote<<");"<<endl;
00976    out<<"   multigraph->SetTitle("<<quote<<GetTitle()<<quote<<");"<<endl;
00977 
00978    if (fGraphs) {
00979       TObjOptLink *lnk = (TObjOptLink*)fGraphs->FirstLink();
00980       TObject *g;
00981 
00982       while (lnk) {
00983          g = lnk->GetObject();
00984          g->SavePrimitive(out, Form("multigraph%s",lnk->GetOption()));
00985          lnk = (TObjOptLink*)lnk->Next();
00986       }
00987    }
00988    const char *l = strstr(option,"th2poly");
00989    if (l) {
00990       out<<"   "<<l+7<<"->AddBin(multigraph);"<<endl;
00991    } else {
00992       out<<"   multigraph->Draw(" <<quote<<option<<quote<<");"<<endl;
00993    }
00994    TAxis *xaxis = GetXaxis();
00995    TAxis *yaxis = GetYaxis();
00996 
00997    if (xaxis) xaxis->SaveAttributes(out, "multigraph","->GetXaxis()");
00998    if (yaxis) yaxis->SaveAttributes(out, "multigraph","->GetYaxis()");
00999 }
01000 
01001 
01002 //______________________________________________________________________________
01003 void TMultiGraph::SetMaximum(Double_t maximum)
01004 {
01005    // Set multigraph maximum.
01006 
01007    fMaximum = maximum;
01008    if (fHistogram)  fHistogram->SetMaximum(maximum);
01009 }
01010 
01011 
01012 //______________________________________________________________________________
01013 void TMultiGraph::SetMinimum(Double_t minimum)
01014 {
01015    // Set multigraph minimum.
01016 
01017    fMinimum = minimum;
01018    if (fHistogram) fHistogram->SetMinimum(minimum);
01019 }

Generated on Tue Jul 5 14:24:21 2011 for ROOT_528-00b_version by  doxygen 1.5.1