PWG1/AliTreeDraw.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 /*
  19
  20 Origin: marian.ivanov@cern.ch
  21 Frequenlty used function for visualization
  22 marian.ivanov@cern.ch
  23 */
  24
  25 #if !defined(__CINT__) || defined(__MAKECINT__)
  26 #include <stdio.h>
  27 #include <string.h>
  28 //ROOT includes
  29 #include "TROOT.h"
  30 #include "Rtypes.h"
  31 #include "TFile.h"
  32 #include "TTree.h"
  33 #include "TChain.h"
  34 #include "TCut.h"
  35 #include "TString.h"
  36 #include "TBenchmark.h"
  37 #include "TStopwatch.h"
  38 #include "TParticle.h"
  39 #include "TSystem.h"
  40 #include "TTimer.h"
  41 #include "TVector3.h"
  42 #include "TH1F.h"
  43 #include "TH2F.h"
  44 #include "TCanvas.h"
  45 #include "TPad.h"
  46 #include "TF1.h"
  47 #include "TView.h"
  48 #include "TView3D.h"
  49 #include "TPolyLine3D.h"
  50 #include "TPolyMarker3D.h"
  51 #include "TObjString.h"
  52
  53
  54 //ALIROOT includes
  55 #include "AliTrackPointArray.h"
  56 #include "AliTreeDraw.h"
  57
  58 #endif
  59
  60 //
  61 //     Class for visualization and some statistacal analysis using tree
  62 //     To be used in comparisons
  63 //                and calib viewers based on tree
  64
  65
  66 ClassImp(AliTreeDraw)
  67
  68
  69 AliTreeDraw::AliTreeDraw():
  70   fTree(0),
  71   fRes(0),
  72   fMean(0),
  73   fPoints(0){
  74   //
  75   // default constructor
  76   //
  77 }
  78
  79 void  AliTreeDraw::ClearHisto(){
  80   //
  81   //
  82   delete fRes;
  83   delete fMean;
  84   fRes=0;
  85   fMean=0;
  86 }
  87
  88
  89
  90 TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection,
  91                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
  92 {
  93   //
  94   Double_t* bins = CreateLogBins(nbins, minx, maxx);
  95   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
  96   char cut[1000];
  97   sprintf(cut,"%s&&%s",selection,quality);
  98   char expression[1000];
  99   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 100   fTree->Draw(expression, cut, "groff");
 101   TH1F* hMean=0;
 102   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 103   AliLabelAxes(hRes, chx, chy);
 104   //
 105   delete hRes2;
 106   delete[] bins;
 107   ClearHisto();
 108   fRes  = hRes;
 109   fMean = hMean;
 110   return hRes;
 111 }
 112
 113
 114
 115 TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
 116                                     const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
 117 {
 118   //
 119   //
 120   //
 121   Double_t* bins = CreateLogBins(nbins, minx, maxx);
 122   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
 123   char cut[1000];
 124   sprintf(cut,"%s&&%s",selection,quality);
 125   char expression[1000];
 126   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 127   fTree->Draw(expression, cut, "groff");
 128   TH1F* hMean=0;
 129   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 130   AliLabelAxes(hRes, chx, chy);
 131   //
 132   delete hRes2;
 133   delete[] bins;
 134   ClearHisto();
 135   fRes  = hRes;
 136   fMean = hMean;
 137   return hRes;
 138 }
 139
 140 ///////////////////////////////////////////////////////////////////////////////////
 141 ///////////////////////////////////////////////////////////////////////////////////
 142 TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality,
 143                               Int_t nbins, Float_t min, Float_t max)
 144 {
 145   //
 146   //
 147   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
 148   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
 149   char inputGen[1000];
 150   sprintf(inputGen,"%s>>hGen", variable);
 151   fTree->Draw(inputGen, selection, "groff");
 152   char selectionRec[256];
 153   sprintf(selectionRec, "%s && %s", selection, quality);
 154   char inputRec[1000];
 155   sprintf(inputRec,"%s>>hRec", variable);
 156   fTree->Draw(inputRec, selectionRec, "groff");
 157   //
 158   TH1F* hEff = CreateEffHisto(hGen, hRec);
 159   AliLabelAxes(hEff, variable, "#epsilon [%]");
 160   fRes = hEff;
 161   delete hRec;
 162   delete hGen;
 163   return hEff;
 164 }
 165
 166
 167
 168 ///////////////////////////////////////////////////////////////////////////////////
 169 ///////////////////////////////////////////////////////////////////////////////////
 170 TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality,
 171                               Int_t nbins, Float_t min, Float_t max)
 172 {
 173   //
 174   //
 175   Double_t* bins = CreateLogBins(nbins, min, max);
 176   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
 177   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
 178   char inputGen[1000];
 179   sprintf(inputGen,"%s>>hGen", variable);
 180   fTree->Draw(inputGen, selection, "groff");
 181   char selectionRec[256];
 182   sprintf(selectionRec, "%s && %s", selection, quality);
 183   char inputRec[1000];
 184   sprintf(inputRec,"%s>>hRec", variable);
 185   fTree->Draw(inputRec, selectionRec, "groff");
 186   //
 187   TH1F* hEff = CreateEffHisto(hGen, hRec);
 188   AliLabelAxes(hEff, variable, "#epsilon [%]");
 189   fRes = hEff;
 190   delete hRec;
 191   delete hGen;
 192   delete[] bins;
 193   return hEff;
 194 }
 195
 196
 197 ///////////////////////////////////////////////////////////////////////////////////
 198 ///////////////////////////////////////////////////////////////////////////////////
 199
 200 Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
 201 {
 202   Double_t* bins = new Double_t[nBins+1];
 203   bins[0] = xMin;
 204   Double_t factor = pow(xMax/xMin, 1./nBins);
 205   for (Int_t i = 1; i <= nBins; i++)
 206     bins[i] = factor * bins[i-1];
 207   return bins;
 208 }
 209
 210
 211
 212
 213 void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
 214 {
 215   //
 216   histo->GetXaxis()->SetTitle(xAxisTitle);
 217   histo->GetYaxis()->SetTitle(yAxisTitle);
 218 }
 219
 220
 221 TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
 222 {
 223   //
 224   Int_t nBins = hGen->GetNbinsX();
 225   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 226   hEff->SetTitle("");
 227   hEff->SetStats(kFALSE);
 228   hEff->SetMinimum(0.);
 229   hEff->SetMaximum(110.);
 230   //
 231   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 232     Double_t nGen = hGen->GetBinContent(iBin);
 233     Double_t nRec = hRec->GetBinContent(iBin);
 234     if (nGen > 0) {
 235       Double_t eff = nRec/nGen;
 236       hEff->SetBinContent(iBin, 100. * eff);
 237       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
 238       //      if (error == 0) error = sqrt(0.1/nGen);
 239       //
 240       if (error == 0) error = 0.0001;
 241       hEff->SetBinError(iBin, 100. * error);
 242     } else {
 243       hEff->SetBinContent(iBin, 100. * 0.5);
 244       hEff->SetBinError(iBin, 100. * 0.5);
 245     }
 246   }
 247   return hEff;
 248 }
 249
 250
 251
 252 TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
 253                      Bool_t overflowBinFits)
 254 {
 255   TVirtualPad* currentPad = gPad;
 256   TAxis* axis = hRes2->GetXaxis();
 257   Int_t nBins = axis->GetNbins();
 258   TH1F* hRes, *hMean;
 259   if (axis->GetXbins()->GetSize()){
 260     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 261     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 262   }
 263   else{
 264     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 265     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 266
 267   }
 268   hRes->SetStats(false);
 269   hRes->SetOption("E");
 270   hRes->SetMinimum(0.);
 271   //
 272   hMean->SetStats(false);
 273   hMean->SetOption("E");
 274
 275   // create the fit function
 276   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 277
 278   fitFunc->SetLineWidth(2);
 279   fitFunc->SetFillStyle(0);
 280   // create canvas for fits
 281   TCanvas* canBinFits = NULL;
 282   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 283   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 284   Int_t ny = (nPads-1) / nx + 1;
 285   if (drawBinFits) {
 286     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 287     if (canBinFits) delete canBinFits;
 288     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 289     canBinFits->Divide(nx, ny);
 290   }
 291
 292   // loop over x bins and fit projection
 293   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 294   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 295     if (drawBinFits) canBinFits->cd(bin + dBin);
 296     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
 297     //
 298     if (hBin->GetEntries() > 5) {
 299       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 300       hBin->Fit(fitFunc,"s");
 301       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 302
 303       if (sigma > 0.){
 304         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 305         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 306       }
 307       else{
 308         hRes->SetBinContent(bin, 0.);
 309         hMean->SetBinContent(bin,0);
 310       }
 311       hRes->SetBinError(bin, fitFunc->GetParError(2));
 312       hMean->SetBinError(bin, fitFunc->GetParError(1));
 313
 314       //
 315       //
 316
 317     } else {
 318       hRes->SetBinContent(bin, 0.);
 319       hRes->SetBinError(bin, 0.);
 320       hMean->SetBinContent(bin, 0.);
 321       hMean->SetBinError(bin, 0.);
 322     }
 323
 324
 325     if (drawBinFits) {
 326       char name[256];
 327       if (bin == 0) {
 328         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 329       } else if (bin == nBins+1) {
 330         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 331       } else {
 332         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 333                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 334       }
 335       canBinFits->cd(bin + dBin);
 336       hBin->SetTitle(name);
 337       hBin->SetStats(kTRUE);
 338       hBin->DrawCopy("E");
 339       canBinFits->Update();
 340       canBinFits->Modified();
 341       canBinFits->Update();
 342     }
 343
 344     delete hBin;
 345   }
 346
 347   delete fitFunc;
 348   currentPad->cd();
 349   *phMean = hMean;
 350   return hRes;
 351 }
 352
 353 TH1F* AliTreeDraw::CreateResHistoI(TH2F* hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits)
 354 {
 355   TVirtualPad* currentPad = gPad;
 356   TAxis* axis = hRes2->GetXaxis();
 357   Int_t nBins = axis->GetNbins();
 358   Bool_t overflowBinFits = kFALSE;
 359   TH1F* hRes, *hMean;
 360   if (axis->GetXbins()->GetSize()){
 361     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 362     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 363   }
 364   else{
 365     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 366     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 367
 368   }
 369   hRes->SetStats(false);
 370   hRes->SetOption("E");
 371   hRes->SetMinimum(0.);
 372   //
 373   hMean->SetStats(false);
 374   hMean->SetOption("E");
 375
 376   // create the fit function
 377   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 378
 379   fitFunc->SetLineWidth(2);
 380   fitFunc->SetFillStyle(0);
 381   // create canvas for fits
 382   TCanvas* canBinFits = NULL;
 383   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 384   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 385   Int_t ny = (nPads-1) / nx + 1;
 386   if (drawBinFits) {
 387     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 388     if (canBinFits) delete canBinFits;
 389     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 390     canBinFits->Divide(nx, ny);
 391   }
 392
 393   // loop over x bins and fit projection
 394   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 395   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 396     if (drawBinFits) canBinFits->cd(bin + dBin);
 397     Int_t bin0=TMath::Max(bin-integ,0);
 398     Int_t bin1=TMath::Min(bin+integ,nBins);
 399     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
 400     //
 401     if (hBin->GetEntries() > 5) {
 402       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 403       hBin->Fit(fitFunc,"s");
 404       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 405
 406       if (sigma > 0.){
 407         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 408         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 409       }
 410       else{
 411         hRes->SetBinContent(bin, 0.);
 412         hMean->SetBinContent(bin,0);
 413       }
 414       hRes->SetBinError(bin, fitFunc->GetParError(2));
 415       hMean->SetBinError(bin, fitFunc->GetParError(1));
 416
 417       //
 418       //
 419
 420     } else {
 421       hRes->SetBinContent(bin, 0.);
 422       hRes->SetBinError(bin, 0.);
 423       hMean->SetBinContent(bin, 0.);
 424       hMean->SetBinError(bin, 0.);
 425     }
 426
 427
 428     if (drawBinFits) {
 429       char name[256];
 430       if (bin == 0) {
 431         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 432       } else if (bin == nBins+1) {
 433         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 434       } else {
 435         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 436                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 437       }
 438       canBinFits->cd(bin + dBin);
 439       hBin->SetTitle(name);
 440       hBin->SetStats(kTRUE);
 441       hBin->DrawCopy("E");
 442       canBinFits->Update();
 443       canBinFits->Modified();
 444       canBinFits->Update();
 445     }
 446
 447     delete hBin;
 448   }
 449
 450   delete fitFunc;
 451   currentPad->cd();
 452   *phMean = hMean;
 453   return hRes;
 454 }
 455
 456 TH1F* AliTreeDraw::CreateResHistoII(TH2F* hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits, Int_t cut)
 457 {
 458   TVirtualPad* currentPad = gPad;
 459   TAxis* axis = hRes2->GetXaxis();
 460   Int_t nBins = axis->GetNbins();
 461   Bool_t overflowBinFits = kFALSE;
 462   TH1F* hRes, *hMean;
 463   if (axis->GetXbins()->GetSize()){
 464     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 465     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 466   }
 467   else{
 468     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 469     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 470
 471   }
 472   hRes->SetStats(false);
 473   hRes->SetOption("E");
 474   hRes->SetMinimum(0.);
 475   //
 476   hMean->SetStats(false);
 477   hMean->SetOption("E");
 478
 479   // create the fit function
 480   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 481
 482   fitFunc->SetLineWidth(2);
 483   fitFunc->SetFillStyle(0);
 484   // create canvas for fits
 485   TCanvas* canBinFits = NULL;
 486   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 487   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 488   Int_t ny = (nPads-1) / nx + 1;
 489   if (drawBinFits) {
 490     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 491     if (canBinFits) delete canBinFits;
 492     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 493     canBinFits->Divide(nx, ny);
 494   }
 495
 496   // loop over x bins and fit projection
 497   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 498   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 499     if (drawBinFits) canBinFits->cd(bin + dBin);
 500     Int_t bin0=TMath::Max(bin-integ,0);
 501     Int_t bin1=TMath::Min(bin+integ,nBins);
 502     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
 503     //
 504     if (hBin->GetEntries() > cut) {
 505       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 506       hBin->Fit(fitFunc,"s");
 507       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 508
 509       if (sigma > 0.){
 510         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 511         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 512       }
 513       else{
 514         hRes->SetBinContent(bin, 0.);
 515         hMean->SetBinContent(bin,0);
 516       }
 517       hRes->SetBinError(bin, fitFunc->GetParError(2));
 518       hMean->SetBinError(bin, fitFunc->GetParError(1));
 519
 520       //
 521       //
 522
 523     } else {
 524       hRes->SetBinContent(bin, 0.);
 525       hRes->SetBinError(bin, 0.);
 526       hMean->SetBinContent(bin, 0.);
 527       hMean->SetBinError(bin, 0.);
 528     }
 529
 530
 531     if (drawBinFits) {
 532       char name[256];
 533       if (bin == 0) {
 534         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 535       } else if (bin == nBins+1) {
 536         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 537       } else {
 538         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 539                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 540       }
 541       canBinFits->cd(bin + dBin);
 542       hBin->SetTitle(name);
 543       hBin->SetStats(kTRUE);
 544       hBin->DrawCopy("E");
 545       canBinFits->Update();
 546       canBinFits->Modified();
 547       canBinFits->Update();
 548     }
 549
 550     delete hBin;
 551   }
 552
 553   delete fitFunc;
 554   currentPad->cd();
 555   *phMean = hMean;
 556   return hRes;
 557 }
 558
 559
 560
 561
 562 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
 563                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
 564   //
 565   // load selected points from tree
 566   //
 567    if (!fPoints) fPoints = new TObjArray;
 568    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
 569    TBranch * br = tpoints->GetBranch(chpoints);
 570    if (!br) return;
 571    AliTrackPointArray * points = new AliTrackPointArray;
 572    br->SetAddress(&points);
 573
 574    Int_t npoints = fTree->Draw(label,selection);
 575    Float_t xyz[30000];
 576    rmin*=rmin;
 577    for (Int_t ii=0;ii<npoints;ii++){
 578      Int_t index = (Int_t)fTree->GetV1()[ii];
 579      tpoints->GetEntryWithIndex(index,index);
 580      if (points->GetNPoints()<2) continue;
 581      Int_t goodpoints=0;
 582      for (Int_t i=0;i<points->GetNPoints(); i++){
 583        xyz[goodpoints*3]   = points->GetX()[i];
 584        xyz[goodpoints*3+1] = points->GetY()[i];
 585        xyz[goodpoints*3+2] = points->GetZ()[i];
 586        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
 587      }
 588      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
 589      marker->SetMarkerColor(color);
 590      marker->SetMarkerStyle(1);
 591      fPoints->AddLast(marker);
 592    }
 593 }
 594
 595
 596
 597
 598 TString* AliTreeDraw::FitPlane(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix, Int_t start, Int_t stop){
 599    //
 600    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 601    // returns chi2, fitParam and covMatrix
 602    // returns TString with fitted formula
 603    //
 604
 605    TString formulaStr(formula);
 606    TString drawStr(drawCommand);
 607    TString cutStr(cuts);
 608
 609    formulaStr.ReplaceAll("++", "~");
 610    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 611    Int_t dim = formulaTokens->GetEntriesFast();
 612
 613    fitParam.ResizeTo(dim);
 614    covMatrix.ResizeTo(dim,dim);
 615
 616    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 617    fitter->StoreData(kTRUE);
 618    fitter->ClearPoints();
 619
 620    Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 621    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 622    Double_t **values = new Double_t*[dim+1] ;
 623
 624    for (Int_t i = 0; i < dim + 1; i++){
 625       Int_t centries = 0;
 626       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 627       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 628
 629       if (entries != centries) return new TString("An ERROR has occured during fitting!");
 630       values[i] = new Double_t[entries];
 631       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t));
 632    }
 633
 634    // add points to the fitter
 635    for (Int_t i = 0; i < entries; i++){
 636       Double_t x[1000];
 637       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 638       fitter->AddPoint(x, values[dim][i], 1);
 639    }
 640
 641    fitter->Eval();
 642    fitter->GetParameters(fitParam);
 643    fitter->GetCovarianceMatrix(covMatrix);
 644    chi2 = fitter->GetChisquare();
 645    chi2 = chi2;
 646
 647    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 648
 649    for (Int_t iparam = 0; iparam < dim; iparam++) {
 650      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 651      if (iparam < dim-1) returnFormula.Append("+");
 652    }
 653    returnFormula.Append(" )");
 654    delete formulaTokens;
 655    delete fitter;
 656    delete[] values;
 657    return preturnFormula;
 658 }
 659