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 "TGeometry.h"
  50 #include "TPolyLine3D.h"
  51 #include "TPolyMarker3D.h"
  52 #include "TObjString.h"
  53
  54
  55 //ALIROOT includes
  56 #include "AliTrackPointArray.h"
  57 #include "AliTreeDraw.h"
  58
  59 #endif
  60
  61 //
  62 //     Class for visualization and some statistacal analysis using tree
  63 //     To be used in comparisons
  64 //                and calib viewers based on tree
  65
  66
  67 ClassImp(AliTreeDraw)
  68
  69
  70 AliTreeDraw::AliTreeDraw():
  71   fTree(0),
  72   fRes(0),
  73   fMean(0),
  74   fPoints(0){
  75   //
  76   // default constructor
  77   //
  78 }
  79
  80 void  AliTreeDraw::ClearHisto(){
  81   //
  82   //
  83   delete fRes;
  84   delete fMean;
  85   fRes=0;
  86   fMean=0;
  87 }
  88
  89
  90
  91 TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection,
  92                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
  93 {
  94   //
  95   Double_t* bins = CreateLogBins(nbins, minx, maxx);
  96   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
  97   char cut[1000];
  98   sprintf(cut,"%s&&%s",selection,quality);
  99   char expression[1000];
 100   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 101   fTree->Draw(expression, cut, "groff");
 102   TH1F* hMean=0;
 103   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 104   AliLabelAxes(hRes, chx, chy);
 105   //
 106   delete hRes2;
 107   delete[] bins;
 108   ClearHisto();
 109   fRes  = hRes;
 110   fMean = hMean;
 111   return hRes;
 112 }
 113
 114
 115
 116 TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
 117                                     const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
 118 {
 119   //
 120   //
 121   //
 122   Double_t* bins = CreateLogBins(nbins, minx, maxx);
 123   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
 124   char cut[1000];
 125   sprintf(cut,"%s&&%s",selection,quality);
 126   char expression[1000];
 127   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 128   fTree->Draw(expression, cut, "groff");
 129   TH1F* hMean=0;
 130   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 131   AliLabelAxes(hRes, chx, chy);
 132   //
 133   delete hRes2;
 134   delete[] bins;
 135   ClearHisto();
 136   fRes  = hRes;
 137   fMean = hMean;
 138   return hRes;
 139 }
 140
 141 ///////////////////////////////////////////////////////////////////////////////////
 142 ///////////////////////////////////////////////////////////////////////////////////
 143 TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality,
 144                               Int_t nbins, Float_t min, Float_t max)
 145 {
 146   //
 147   //
 148   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
 149   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
 150   char inputGen[1000];
 151   sprintf(inputGen,"%s>>hGen", variable);
 152   fTree->Draw(inputGen, selection, "groff");
 153   char selectionRec[256];
 154   sprintf(selectionRec, "%s && %s", selection, quality);
 155   char inputRec[1000];
 156   sprintf(inputRec,"%s>>hRec", variable);
 157   fTree->Draw(inputRec, selectionRec, "groff");
 158   //
 159   TH1F* hEff = CreateEffHisto(hGen, hRec);
 160   AliLabelAxes(hEff, variable, "#epsilon [%]");
 161   fRes = hEff;
 162   delete hRec;
 163   delete hGen;
 164   return hEff;
 165 }
 166
 167
 168
 169 ///////////////////////////////////////////////////////////////////////////////////
 170 ///////////////////////////////////////////////////////////////////////////////////
 171 TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality,
 172                               Int_t nbins, Float_t min, Float_t max)
 173 {
 174   //
 175   //
 176   Double_t* bins = CreateLogBins(nbins, min, max);
 177   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
 178   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
 179   char inputGen[1000];
 180   sprintf(inputGen,"%s>>hGen", variable);
 181   fTree->Draw(inputGen, selection, "groff");
 182   char selectionRec[256];
 183   sprintf(selectionRec, "%s && %s", selection, quality);
 184   char inputRec[1000];
 185   sprintf(inputRec,"%s>>hRec", variable);
 186   fTree->Draw(inputRec, selectionRec, "groff");
 187   //
 188   TH1F* hEff = CreateEffHisto(hGen, hRec);
 189   AliLabelAxes(hEff, variable, "#epsilon [%]");
 190   fRes = hEff;
 191   delete hRec;
 192   delete hGen;
 193   delete[] bins;
 194   return hEff;
 195 }
 196
 197
 198 ///////////////////////////////////////////////////////////////////////////////////
 199 ///////////////////////////////////////////////////////////////////////////////////
 200
 201 Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
 202 {
 203   Double_t* bins = new Double_t[nBins+1];
 204   bins[0] = xMin;
 205   Double_t factor = pow(xMax/xMin, 1./nBins);
 206   for (Int_t i = 1; i <= nBins; i++)
 207     bins[i] = factor * bins[i-1];
 208   return bins;
 209 }
 210
 211
 212
 213
 214 void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
 215 {
 216   //
 217   histo->GetXaxis()->SetTitle(xAxisTitle);
 218   histo->GetYaxis()->SetTitle(yAxisTitle);
 219 }
 220
 221
 222 TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
 223 {
 224   //
 225   Int_t nBins = hGen->GetNbinsX();
 226   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 227   hEff->SetTitle("");
 228   hEff->SetStats(kFALSE);
 229   hEff->SetMinimum(0.);
 230   hEff->SetMaximum(110.);
 231   //
 232   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 233     Double_t nGen = hGen->GetBinContent(iBin);
 234     Double_t nRec = hRec->GetBinContent(iBin);
 235     if (nGen > 0) {
 236       Double_t eff = nRec/nGen;
 237       hEff->SetBinContent(iBin, 100. * eff);
 238       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
 239       //      if (error == 0) error = sqrt(0.1/nGen);
 240       //
 241       if (error == 0) error = 0.0001;
 242       hEff->SetBinError(iBin, 100. * error);
 243     } else {
 244       hEff->SetBinContent(iBin, 100. * 0.5);
 245       hEff->SetBinError(iBin, 100. * 0.5);
 246     }
 247   }
 248   return hEff;
 249 }
 250
 251
 252
 253 TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
 254                      Bool_t overflowBinFits)
 255 {
 256   TVirtualPad* currentPad = gPad;
 257   TAxis* axis = hRes2->GetXaxis();
 258   Int_t nBins = axis->GetNbins();
 259   TH1F* hRes, *hMean;
 260   if (axis->GetXbins()->GetSize()){
 261     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 262     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 263   }
 264   else{
 265     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 266     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 267
 268   }
 269   hRes->SetStats(false);
 270   hRes->SetOption("E");
 271   hRes->SetMinimum(0.);
 272   //
 273   hMean->SetStats(false);
 274   hMean->SetOption("E");
 275
 276   // create the fit function
 277   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 278
 279   fitFunc->SetLineWidth(2);
 280   fitFunc->SetFillStyle(0);
 281   // create canvas for fits
 282   TCanvas* canBinFits = NULL;
 283   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 284   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 285   Int_t ny = (nPads-1) / nx + 1;
 286   if (drawBinFits) {
 287     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 288     if (canBinFits) delete canBinFits;
 289     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 290     canBinFits->Divide(nx, ny);
 291   }
 292
 293   // loop over x bins and fit projection
 294   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 295   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 296     if (drawBinFits) canBinFits->cd(bin + dBin);
 297     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
 298     //
 299     if (hBin->GetEntries() > 5) {
 300       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 301       hBin->Fit(fitFunc,"s");
 302       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 303
 304       if (sigma > 0.){
 305         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 306         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 307       }
 308       else{
 309         hRes->SetBinContent(bin, 0.);
 310         hMean->SetBinContent(bin,0);
 311       }
 312       hRes->SetBinError(bin, fitFunc->GetParError(2));
 313       hMean->SetBinError(bin, fitFunc->GetParError(1));
 314
 315       //
 316       //
 317
 318     } else {
 319       hRes->SetBinContent(bin, 0.);
 320       hRes->SetBinError(bin, 0.);
 321       hMean->SetBinContent(bin, 0.);
 322       hMean->SetBinError(bin, 0.);
 323     }
 324
 325
 326     if (drawBinFits) {
 327       char name[256];
 328       if (bin == 0) {
 329         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 330       } else if (bin == nBins+1) {
 331         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 332       } else {
 333         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 334                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 335       }
 336       canBinFits->cd(bin + dBin);
 337       hBin->SetTitle(name);
 338       hBin->SetStats(kTRUE);
 339       hBin->DrawCopy("E");
 340       canBinFits->Update();
 341       canBinFits->Modified();
 342       canBinFits->Update();
 343     }
 344
 345     delete hBin;
 346   }
 347
 348   delete fitFunc;
 349   currentPad->cd();
 350   *phMean = hMean;
 351   return hRes;
 352 }
 353
 354 TH1F* AliTreeDraw::CreateResHistoI(TH2F* hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits)
 355 {
 356   TVirtualPad* currentPad = gPad;
 357   TAxis* axis = hRes2->GetXaxis();
 358   Int_t nBins = axis->GetNbins();
 359   Bool_t overflowBinFits = kFALSE;
 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 nx = Int_t(sqrt(nPads-1.));// + 1;
 386   Int_t ny = (nPads-1) / nx + 1;
 387   if (drawBinFits) {
 388     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 389     if (canBinFits) delete canBinFits;
 390     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 391     canBinFits->Divide(nx, ny);
 392   }
 393
 394   // loop over x bins and fit projection
 395   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 396   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 397     if (drawBinFits) canBinFits->cd(bin + dBin);
 398     Int_t bin0=TMath::Max(bin-integ,0);
 399     Int_t bin1=TMath::Min(bin+integ,nBins);
 400     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
 401     //
 402     if (hBin->GetEntries() > 5) {
 403       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 404       hBin->Fit(fitFunc,"s");
 405       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 406
 407       if (sigma > 0.){
 408         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 409         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 410       }
 411       else{
 412         hRes->SetBinContent(bin, 0.);
 413         hMean->SetBinContent(bin,0);
 414       }
 415       hRes->SetBinError(bin, fitFunc->GetParError(2));
 416       hMean->SetBinError(bin, fitFunc->GetParError(1));
 417
 418       //
 419       //
 420
 421     } else {
 422       hRes->SetBinContent(bin, 0.);
 423       hRes->SetBinError(bin, 0.);
 424       hMean->SetBinContent(bin, 0.);
 425       hMean->SetBinError(bin, 0.);
 426     }
 427
 428
 429     if (drawBinFits) {
 430       char name[256];
 431       if (bin == 0) {
 432         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 433       } else if (bin == nBins+1) {
 434         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 435       } else {
 436         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 437                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 438       }
 439       canBinFits->cd(bin + dBin);
 440       hBin->SetTitle(name);
 441       hBin->SetStats(kTRUE);
 442       hBin->DrawCopy("E");
 443       canBinFits->Update();
 444       canBinFits->Modified();
 445       canBinFits->Update();
 446     }
 447
 448     delete hBin;
 449   }
 450
 451   delete fitFunc;
 452   currentPad->cd();
 453   *phMean = hMean;
 454   return hRes;
 455 }
 456
 457
 458
 459
 460 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
 461                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
 462   //
 463   // load selected points from tree
 464   //
 465    if (!fPoints) fPoints = new TObjArray;
 466    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
 467    TBranch * br = tpoints->GetBranch(chpoints);
 468    if (!br) return;
 469    AliTrackPointArray * points = new AliTrackPointArray;
 470    br->SetAddress(&points);
 471
 472    Int_t npoints = fTree->Draw(label,selection);
 473    Float_t xyz[30000];
 474    rmin*=rmin;
 475    for (Int_t i=0;i<npoints;i++){
 476      Int_t index = (Int_t)fTree->GetV1()[i];
 477      tpoints->GetEntryWithIndex(index,index);
 478      if (points->GetNPoints()<2) continue;
 479      Int_t goodpoints=0;
 480      for (Int_t i=0;i<points->GetNPoints(); i++){
 481        xyz[goodpoints*3]   = points->GetX()[i];
 482        xyz[goodpoints*3+1] = points->GetY()[i];
 483        xyz[goodpoints*3+2] = points->GetZ()[i];
 484        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
 485      }
 486      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
 487      marker->SetMarkerColor(color);
 488      marker->SetMarkerStyle(1);
 489      fPoints->AddLast(marker);
 490    }
 491 }
 492
 493
 494
 495
 496 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){
 497    //
 498    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 499    // returns chi2, fitParam and covMatrix
 500    // returns TString with fitted formula
 501    //
 502
 503    TString formulaStr(formula);
 504    TString drawStr(drawCommand);
 505    TString cutStr(cuts);
 506
 507    formulaStr.ReplaceAll("++", "~");
 508    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 509    Int_t dim = formulaTokens->GetEntriesFast();
 510
 511    fitParam.ResizeTo(dim);
 512    covMatrix.ResizeTo(dim,dim);
 513
 514    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 515    fitter->StoreData(kTRUE);
 516    fitter->ClearPoints();
 517
 518    Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 519    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 520    Double_t **values = new Double_t*[dim+1] ;
 521
 522    for (Int_t i = 0; i < dim + 1; i++){
 523       Int_t centries = 0;
 524       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 525       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 526
 527       if (entries != centries) return new TString("An ERROR has occured during fitting!");
 528       values[i] = new Double_t[entries];
 529       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t));
 530    }
 531
 532    // add points to the fitter
 533    for (Int_t i = 0; i < entries; i++){
 534       Double_t x[1000];
 535       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 536       fitter->AddPoint(x, values[dim][i], 1);
 537    }
 538
 539    fitter->Eval();
 540    fitter->GetParameters(fitParam);
 541    fitter->GetCovarianceMatrix(covMatrix);
 542    chi2 = fitter->GetChisquare();
 543    chi2 = chi2;
 544
 545    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 546
 547    for (Int_t iparam = 0; iparam < dim; iparam++) {
 548      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 549      if (iparam < dim-1) returnFormula.Append("+");
 550    }
 551    returnFormula.Append(" )");
 552    delete formulaTokens;
 553    delete fitter;
 554    delete[] values;
 555    return preturnFormula;
 556 }
 557