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
  53 //ALIROOT includes
  54 #include "AliTrackPointArray.h"
  55 #include "AliTreeDraw.h"
  56
  57 #endif
  58
  59 //
  60 //     Class for visualization and some statistacal analysis using tree
  61 //     To be used in comparisons
  62 //                and calib viewers based on tree
  63
  64
  65 ClassImp(AliTreeDraw)
  66
  67
  68 AliTreeDraw::AliTreeDraw():
  69   fTree(0),
  70   fRes(0),
  71   fMean(0),
  72   fPoints(0){
  73   //
  74   // default constructor
  75   //
  76 }
  77
  78 void  AliTreeDraw::ClearHisto(){
  79   //
  80   //
  81   delete fRes;
  82   delete fMean;
  83   fRes=0;
  84   fMean=0;
  85 }
  86
  87
  88
  89 TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection,
  90                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
  91 {
  92   //
  93   Double_t* bins = CreateLogBins(nbins, minx, maxx);
  94   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
  95   char cut[1000];
  96   sprintf(cut,"%s&&%s",selection,quality);
  97   char expression[1000];
  98   sprintf(expression,"%s:%s>>hRes2",chy,chx);
  99   fTree->Draw(expression, cut, "groff");
 100   TH1F* hMean=0;
 101   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 102   AliLabelAxes(hRes, chx, chy);
 103   //
 104   delete hRes2;
 105   delete[] bins;
 106   ClearHisto();
 107   fRes  = hRes;
 108   fMean = hMean;
 109   return hRes;
 110 }
 111
 112
 113
 114 TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
 115                                     const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
 116 {
 117   //
 118   //
 119   //
 120   Double_t* bins = CreateLogBins(nbins, minx, maxx);
 121   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
 122   char cut[1000];
 123   sprintf(cut,"%s&&%s",selection,quality);
 124   char expression[1000];
 125   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 126   fTree->Draw(expression, cut, "groff");
 127   TH1F* hMean=0;
 128   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 129   AliLabelAxes(hRes, chx, chy);
 130   //
 131   delete hRes2;
 132   delete[] bins;
 133   ClearHisto();
 134   fRes  = hRes;
 135   fMean = hMean;
 136   return hRes;
 137 }
 138
 139 ///////////////////////////////////////////////////////////////////////////////////
 140 ///////////////////////////////////////////////////////////////////////////////////
 141 TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality,
 142                               Int_t nbins, Float_t min, Float_t max)
 143 {
 144   //
 145   //
 146   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
 147   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
 148   char inputGen[1000];
 149   sprintf(inputGen,"%s>>hGen", variable);
 150   fTree->Draw(inputGen, selection, "groff");
 151   char selectionRec[256];
 152   sprintf(selectionRec, "%s && %s", selection, quality);
 153   char inputRec[1000];
 154   sprintf(inputRec,"%s>>hRec", variable);
 155   fTree->Draw(inputRec, selectionRec, "groff");
 156   //
 157   TH1F* hEff = CreateEffHisto(hGen, hRec);
 158   AliLabelAxes(hEff, variable, "#epsilon [%]");
 159   fRes = hEff;
 160   delete hRec;
 161   delete hGen;
 162   return hEff;
 163 }
 164
 165
 166
 167 ///////////////////////////////////////////////////////////////////////////////////
 168 ///////////////////////////////////////////////////////////////////////////////////
 169 TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality,
 170                               Int_t nbins, Float_t min, Float_t max)
 171 {
 172   //
 173   //
 174   Double_t* bins = CreateLogBins(nbins, min, max);
 175   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
 176   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
 177   char inputGen[1000];
 178   sprintf(inputGen,"%s>>hGen", variable);
 179   fTree->Draw(inputGen, selection, "groff");
 180   char selectionRec[256];
 181   sprintf(selectionRec, "%s && %s", selection, quality);
 182   char inputRec[1000];
 183   sprintf(inputRec,"%s>>hRec", variable);
 184   fTree->Draw(inputRec, selectionRec, "groff");
 185   //
 186   TH1F* hEff = CreateEffHisto(hGen, hRec);
 187   AliLabelAxes(hEff, variable, "#epsilon [%]");
 188   fRes = hEff;
 189   delete hRec;
 190   delete hGen;
 191   delete[] bins;
 192   return hEff;
 193 }
 194
 195
 196 ///////////////////////////////////////////////////////////////////////////////////
 197 ///////////////////////////////////////////////////////////////////////////////////
 198
 199 Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
 200 {
 201   Double_t* bins = new Double_t[nBins+1];
 202   bins[0] = xMin;
 203   Double_t factor = pow(xMax/xMin, 1./nBins);
 204   for (Int_t i = 1; i <= nBins; i++)
 205     bins[i] = factor * bins[i-1];
 206   return bins;
 207 }
 208
 209
 210
 211
 212 void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
 213 {
 214   //
 215   histo->GetXaxis()->SetTitle(xAxisTitle);
 216   histo->GetYaxis()->SetTitle(yAxisTitle);
 217 }
 218
 219
 220 TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
 221 {
 222   //
 223   Int_t nBins = hGen->GetNbinsX();
 224   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 225   hEff->SetTitle("");
 226   hEff->SetStats(kFALSE);
 227   hEff->SetMinimum(0.);
 228   hEff->SetMaximum(110.);
 229   //
 230   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 231     Double_t nGen = hGen->GetBinContent(iBin);
 232     Double_t nRec = hRec->GetBinContent(iBin);
 233     if (nGen > 0) {
 234       Double_t eff = nRec/nGen;
 235       hEff->SetBinContent(iBin, 100. * eff);
 236       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
 237       //      if (error == 0) error = sqrt(0.1/nGen);
 238       //
 239       if (error == 0) error = 0.0001;
 240       hEff->SetBinError(iBin, 100. * error);
 241     } else {
 242       hEff->SetBinContent(iBin, 100. * 0.5);
 243       hEff->SetBinError(iBin, 100. * 0.5);
 244     }
 245   }
 246   return hEff;
 247 }
 248
 249
 250
 251 TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
 252                      Bool_t overflowBinFits)
 253 {
 254   TVirtualPad* currentPad = gPad;
 255   TAxis* axis = hRes2->GetXaxis();
 256   Int_t nBins = axis->GetNbins();
 257   TH1F* hRes, *hMean;
 258   if (axis->GetXbins()->GetSize()){
 259     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 260     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 261   }
 262   else{
 263     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 264     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 265
 266   }
 267   hRes->SetStats(false);
 268   hRes->SetOption("E");
 269   hRes->SetMinimum(0.);
 270   //
 271   hMean->SetStats(false);
 272   hMean->SetOption("E");
 273
 274   // create the fit function
 275   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 276
 277   fitFunc->SetLineWidth(2);
 278   fitFunc->SetFillStyle(0);
 279   // create canvas for fits
 280   TCanvas* canBinFits = NULL;
 281   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 282   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 283   Int_t ny = (nPads-1) / nx + 1;
 284   if (drawBinFits) {
 285     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 286     if (canBinFits) delete canBinFits;
 287     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 288     canBinFits->Divide(nx, ny);
 289   }
 290
 291   // loop over x bins and fit projection
 292   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 293   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 294     if (drawBinFits) canBinFits->cd(bin + dBin);
 295     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
 296     //
 297     if (hBin->GetEntries() > 5) {
 298       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 299       hBin->Fit(fitFunc,"s");
 300       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 301
 302       if (sigma > 0.){
 303         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 304         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 305       }
 306       else{
 307         hRes->SetBinContent(bin, 0.);
 308         hMean->SetBinContent(bin,0);
 309       }
 310       hRes->SetBinError(bin, fitFunc->GetParError(2));
 311       hMean->SetBinError(bin, fitFunc->GetParError(1));
 312
 313       //
 314       //
 315
 316     } else {
 317       hRes->SetBinContent(bin, 0.);
 318       hRes->SetBinError(bin, 0.);
 319       hMean->SetBinContent(bin, 0.);
 320       hMean->SetBinError(bin, 0.);
 321     }
 322
 323
 324     if (drawBinFits) {
 325       char name[256];
 326       if (bin == 0) {
 327         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 328       } else if (bin == nBins+1) {
 329         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 330       } else {
 331         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 332                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 333       }
 334       canBinFits->cd(bin + dBin);
 335       hBin->SetTitle(name);
 336       hBin->SetStats(kTRUE);
 337       hBin->DrawCopy("E");
 338       canBinFits->Update();
 339       canBinFits->Modified();
 340       canBinFits->Update();
 341     }
 342
 343     delete hBin;
 344   }
 345
 346   delete fitFunc;
 347   currentPad->cd();
 348   *phMean = hMean;
 349   return hRes;
 350 }
 351
 352
 353
 354
 355 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
 356                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
 357   //
 358   // load selected points from tree
 359   //
 360    if (!fPoints) fPoints = new TObjArray;
 361    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
 362    TBranch * br = tpoints->GetBranch(chpoints);
 363    if (!br) return;
 364    AliTrackPointArray * points = new AliTrackPointArray;
 365    br->SetAddress(&points);
 366
 367    Int_t npoints = fTree->Draw(label,selection);
 368    Float_t xyz[30000];
 369    rmin*=rmin;
 370    for (Int_t i=0;i<npoints;i++){
 371      Int_t index = (Int_t)fTree->GetV1()[i];
 372      tpoints->GetEntryWithIndex(index,index);
 373      if (points->GetNPoints()<2) continue;
 374      Int_t goodpoints=0;
 375      for (Int_t i=0;i<points->GetNPoints(); i++){
 376        xyz[goodpoints*3]   = points->GetX()[i];
 377        xyz[goodpoints*3+1] = points->GetY()[i];
 378        xyz[goodpoints*3+2] = points->GetZ()[i];
 379        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
 380      }
 381      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
 382      marker->SetMarkerColor(color);
 383      marker->SetMarkerStyle(1);
 384      fPoints->AddLast(marker);
 385    }
 386 }
 387