PWGPP/TPC/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   snprintf(cut,1000,"%s&&%s",selection,quality);
  98   char expression[1000];
  99   snprintf(expression,1000,"%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   snprintf(cut,1000,"%s&&%s",selection,quality);
 125   char expression[1000];
 126   snprintf(expression,1000,"%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   snprintf(inputGen,1000,"%s>>hGen", variable);
 151   fTree->Draw(inputGen, selection, "groff");
 152   char selectionRec[256];
 153   snprintf(selectionRec,256, "%s && %s", selection, quality);
 154   char inputRec[1000];
 155   snprintf(inputRec,1000,"%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   snprintf(inputGen,1000,"%s>>hGen", variable);
 180   fTree->Draw(inputGen, selection, "groff");
 181   char selectionRec[256];
 182   snprintf(selectionRec,256, "%s && %s", selection, quality);
 183   char inputRec[1000];
 184   snprintf(inputRec,1000,"%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         snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 329       } else if (bin == nBins+1) {
 330         snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 331       } else {
 332         snprintf(name,256, "%.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
 360   TH1F* hRes, *hMean;
 361   if (axis->GetXbins()->GetSize()){
 362     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 363     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 364   }
 365   else{
 366     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 367     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 368
 369   }
 370   hRes->SetStats(false);
 371   hRes->SetOption("E");
 372   hRes->SetMinimum(0.);
 373   //
 374   hMean->SetStats(false);
 375   hMean->SetOption("E");
 376
 377   // create the fit function
 378   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 379
 380   fitFunc->SetLineWidth(2);
 381   fitFunc->SetFillStyle(0);
 382   // create canvas for fits
 383   TCanvas* canBinFits = NULL;
 384   //Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 385   Int_t nPads = nBins;
 386   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 387   Int_t ny = (nPads-1) / nx + 1;
 388   if (drawBinFits) {
 389     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 390     if (canBinFits) delete canBinFits;
 391     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 392     canBinFits->Divide(nx, ny);
 393   }
 394
 395   // loop over x bins and fit projection
 396   //Int_t dBin = ((overflowBinFits) ? 1 : 0);
 397   Int_t dBin =  0;
 398   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 399     if (drawBinFits) canBinFits->cd(bin + dBin);
 400     Int_t bin0=TMath::Max(bin-integ,0);
 401     Int_t bin1=TMath::Min(bin+integ,nBins);
 402     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
 403     //
 404     if (hBin->GetEntries() > 5) {
 405       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 406       hBin->Fit(fitFunc,"s");
 407       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 408
 409       if (sigma > 0.){
 410         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 411         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 412       }
 413       else{
 414         hRes->SetBinContent(bin, 0.);
 415         hMean->SetBinContent(bin,0);
 416       }
 417       hRes->SetBinError(bin, fitFunc->GetParError(2));
 418       hMean->SetBinError(bin, fitFunc->GetParError(1));
 419
 420       //
 421       //
 422
 423     } else {
 424       hRes->SetBinContent(bin, 0.);
 425       hRes->SetBinError(bin, 0.);
 426       hMean->SetBinContent(bin, 0.);
 427       hMean->SetBinError(bin, 0.);
 428     }
 429
 430
 431     if (drawBinFits) {
 432       char name[256];
 433       if (bin == 0) {
 434         snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 435       } else if (bin == nBins+1) {
 436         snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 437       } else {
 438         snprintf(name,256, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 439                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 440       }
 441       canBinFits->cd(bin + dBin);
 442       hBin->SetTitle(name);
 443       hBin->SetStats(kTRUE);
 444       hBin->DrawCopy("E");
 445       canBinFits->Update();
 446       canBinFits->Modified();
 447       canBinFits->Update();
 448     }
 449
 450     delete hBin;
 451   }
 452
 453   delete fitFunc;
 454   currentPad->cd();
 455   *phMean = hMean;
 456   return hRes;
 457 }
 458
 459 TH1F* AliTreeDraw::CreateResHistoII(TH2F* hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits, Int_t cut)
 460 {
 461   TVirtualPad* currentPad = gPad;
 462   TAxis* axis = hRes2->GetXaxis();
 463   Int_t nBins = axis->GetNbins();
 464   //Bool_t overflowBinFits = kFALSE;
 465   TH1F* hRes, *hMean;
 466   if (axis->GetXbins()->GetSize()){
 467     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 468     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 469   }
 470   else{
 471     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 472     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 473
 474   }
 475   hRes->SetStats(false);
 476   hRes->SetOption("E");
 477   hRes->SetMinimum(0.);
 478   //
 479   hMean->SetStats(false);
 480   hMean->SetOption("E");
 481
 482   // create the fit function
 483   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 484
 485   fitFunc->SetLineWidth(2);
 486   fitFunc->SetFillStyle(0);
 487   // create canvas for fits
 488   TCanvas* canBinFits = NULL;
 489   //Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 490   Int_t nPads = nBins;
 491   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 492   Int_t ny = (nPads-1) / nx + 1;
 493   if (drawBinFits) {
 494     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 495     if (canBinFits) delete canBinFits;
 496     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 497     canBinFits->Divide(nx, ny);
 498   }
 499
 500   // loop over x bins and fit projection
 501   //Int_t dBin = ((overflowBinFits) ? 1 : 0);
 502   Int_t dBin = 0;
 503   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 504     if (drawBinFits) canBinFits->cd(bin + dBin);
 505     Int_t bin0=TMath::Max(bin-integ,0);
 506     Int_t bin1=TMath::Min(bin+integ,nBins);
 507     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
 508     //
 509     if (hBin->GetEntries() > cut) {
 510       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 511       hBin->Fit(fitFunc,"s");
 512       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 513
 514       if (sigma > 0.){
 515         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 516         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 517       }
 518       else{
 519         hRes->SetBinContent(bin, 0.);
 520         hMean->SetBinContent(bin,0);
 521       }
 522       hRes->SetBinError(bin, fitFunc->GetParError(2));
 523       hMean->SetBinError(bin, fitFunc->GetParError(1));
 524
 525       //
 526       //
 527
 528     } else {
 529       hRes->SetBinContent(bin, 0.);
 530       hRes->SetBinError(bin, 0.);
 531       hMean->SetBinContent(bin, 0.);
 532       hMean->SetBinError(bin, 0.);
 533     }
 534
 535
 536     if (drawBinFits) {
 537       char name[256];
 538       if (bin == 0) {
 539         snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 540       } else if (bin == nBins+1) {
 541         snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 542       } else {
 543         snprintf(name,256, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 544                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 545       }
 546       canBinFits->cd(bin + dBin);
 547       hBin->SetTitle(name);
 548       hBin->SetStats(kTRUE);
 549       hBin->DrawCopy("E");
 550       canBinFits->Update();
 551       canBinFits->Modified();
 552       canBinFits->Update();
 553     }
 554
 555     delete hBin;
 556   }
 557
 558   delete fitFunc;
 559   currentPad->cd();
 560   *phMean = hMean;
 561   return hRes;
 562 }
 563
 564
 565
 566
 567 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
 568                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
 569   //
 570   // load selected points from tree
 571   //
 572    if (!fPoints) fPoints = new TObjArray;
 573    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
 574    TBranch * br = tpoints->GetBranch(chpoints);
 575    if (!br) return;
 576    AliTrackPointArray * points = new AliTrackPointArray;
 577    br->SetAddress(&points);
 578
 579    Int_t npoints = fTree->Draw(label,selection);
 580    Float_t xyz[30000];
 581    rmin*=rmin;
 582    for (Int_t ii=0;ii<npoints;ii++){
 583      Int_t index = (Int_t)fTree->GetV1()[ii];
 584      tpoints->GetEntryWithIndex(index,index);
 585      if (points->GetNPoints()<2) continue;
 586      Int_t goodpoints=0;
 587      for (Int_t i=0;i<points->GetNPoints(); i++){
 588        xyz[goodpoints*3]   = points->GetX()[i];
 589        xyz[goodpoints*3+1] = points->GetY()[i];
 590        xyz[goodpoints*3+2] = points->GetZ()[i];
 591        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
 592      }
 593      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
 594      marker->SetMarkerColor(color);
 595      marker->SetMarkerStyle(1);
 596      fPoints->AddLast(marker);
 597    }
 598 }
 599
 600
 601
 602
 603 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){
 604    //
 605    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 606    // returns chi2, fitParam and covMatrix
 607    // returns TString with fitted formula
 608    //
 609
 610    TString formulaStr(formula);
 611    TString drawStr(drawCommand);
 612    TString cutStr(cuts);
 613
 614    formulaStr.ReplaceAll("++", "~");
 615    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 616    Int_t dim = formulaTokens->GetEntriesFast();
 617
 618    fitParam.ResizeTo(dim);
 619    covMatrix.ResizeTo(dim,dim);
 620
 621    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 622    fitter->StoreData(kTRUE);
 623    fitter->ClearPoints();
 624
 625    Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 626    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 627    Double_t **values = new Double_t*[dim+1] ;
 628
 629    for (Int_t i = 0; i < dim + 1; i++) {
 630       values[i] = 0;
 631    }
 632
 633    for (Int_t i = 0; i < dim + 1; i++) {
 634       Int_t centries = 0;
 635       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 636       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 637
 638       if (entries != centries) {
 639         for (Int_t j = 0; j < dim + 1; j++) {
 640           if(values[j]) delete values[j];
 641         }
 642         delete[] values;
 643         return new TString("An ERROR has occured during fitting!");
 644       }
 645       else {
 646         values[i] = new Double_t[entries];
 647         memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t));
 648       }
 649    }
 650
 651    // add points to the fitter
 652    for (Int_t i = 0; i < entries; i++) {
 653       Double_t x[1000];
 654       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 655       fitter->AddPoint(x, values[dim][i], 1);
 656    }
 657
 658    fitter->Eval();
 659    fitter->GetParameters(fitParam);
 660    fitter->GetCovarianceMatrix(covMatrix);
 661    chi2 = fitter->GetChisquare();
 662
 663    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 664
 665    for (Int_t iparam = 0; iparam < dim; iparam++) {
 666      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 667      if (iparam < dim-1) returnFormula.Append("+");
 668    }
 669    returnFormula.Append(" )");
 670    delete formulaTokens;
 671    delete fitter;
 672
 673    for (Int_t i = 0; i < dim + 1; i++) {
 674      delete[] values[i];
 675    }
 676    delete[] values;
 677    return preturnFormula;
 678 }
 679