PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 //  Tools for drawing/ manipulating histograms. (NEEDS CLEANUP!)
  18 // -----------------------------------------------------------------------
  19 //  Author: Misha Veldhoen (misha.veldhoen@cern.ch)
  20
  21 #include "AliHistToolsDiHadronPID.h"
  22
  23 #include <iostream>
  24 using namespace std;
  25
  26 #include "TLegend.h"
  27 #include "TMath.h"
  28 #include "TF1.h"
  29 #include "THn.h"
  30 #include "TObjArray.h"
  31 #include "TStyle.h"
  32 #include "TArray.h"
  33
  34 // -----------------------------------------------------------------------
  35 TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(
  36         const TH1F* histIn, const Double_t* binsx, Int_t Nbinsx, Bool_t density) {
  37
  38         // Rebins a histogram (hin) with a variable binning to a histogram
  39         // with another variable binning (binsx). If the "density" flag is set,
  40         // then it is expected that the bins are per unit x-axis, otherwise an
  41         // absolute count is assumed.
  42         //
  43         // TODO: determine over/under-flow bins.
  44
  45         // Gather info from the original histogram.
  46         const TArrayD* binsxIn = histIn->GetXaxis()->GetXbins();
  47         TArrayD* binsxOut = new TArrayD(Nbinsx + 1,binsx);
  48         Int_t NbinsxIn = binsxIn->GetSize() - 1;
  49         const char* nameIn = histIn->GetName();
  50         const char* titleIn = histIn->GetTitle();
  51         const char* xAxisTitleIn = histIn->GetXaxis()->GetTitle();
  52         const char* yAxisTitleIn = histIn->GetYaxis()->GetTitle();
  53
  54         // Gather info from the output histogram and create it.
  55         Int_t NbinsxOut = binsxOut->GetSize() - 1;
  56         const char* nameOut = Form("%s (rebinned)",nameIn);
  57         const char* titleOut = Form("%s_rebinned;%s;%s",titleIn,xAxisTitleIn,yAxisTitleIn);
  58         TH1F* histOut = new TH1F(nameOut,titleOut,NbinsxOut,binsxOut->GetArray());
  59         //histOut->Sumw2();
  60
  61         // Filling the bins.
  62         for (Int_t iBinOut = 0; iBinOut < NbinsxOut; iBinOut++) {
  63
  64                 Double_t binOutLowEdge = binsxOut->At(iBinOut);
  65                 Double_t binOutHighEdge = binsxOut->At(iBinOut+1);
  66                 Double_t binOutWidth = binOutHighEdge - binOutLowEdge;
  67                 //Double_t binOutCenter = binOutHighEdge - binOutWidth/2.;
  68
  69                 // Setting all errors of the outgoing histogram to zero (needed later):
  70                 histOut->SetBinError(iBinOut+1,0.);
  71
  72                 for (Int_t iBinIn = 0; iBinIn < NbinsxIn; iBinIn++) {
  73
  74                         Double_t binInLowEdge = binsxIn->At(iBinIn);
  75                         Double_t binInHighEdge = binsxIn->At(iBinIn+1);
  76                         Double_t binInWidth = binInHighEdge - binInLowEdge;
  77                         //Double_t binInCenter = binInHighEdge - binInWidth/2.;
  78                         Double_t binInContent = histIn->GetBinContent(iBinIn+1);
  79                         Double_t binInError2 = histIn->GetBinError(iBinIn+1)*histIn->GetBinError(iBinIn+1);
  80
  81                         /* -- Determining what percentage of the in-bin is included
  82                            in the current out-bin -- */
  83                         Double_t PercentageIncluded = 0.;
  84
  85                         // Position of the low edge of the in-bin:
  86                         // In  1   2   3
  87                         // Out   |   |
  88                         Int_t binInLowEdgePos = 0;
  89                         if (binInLowEdge < binOutLowEdge) binInLowEdgePos = 1;
  90                         else if (binInLowEdge > binOutHighEdge) binInLowEdgePos = 3;
  91                         else binInLowEdgePos = 2;
  92
  93                         // Same for the high edge of the in-bin:
  94                         Int_t binInHighEdgePos = 0;
  95                         if (binInHighEdge < binOutLowEdge) binInHighEdgePos = 1;
  96                         else if (binInHighEdge > binOutHighEdge) binInHighEdgePos = 3;
  97                         else binInHighEdgePos = 2;
  98
  99                         // In-bin not included.
 100                         // In        | |     or     | |
 101                         // Out |   |                    |    |
 102                         if ( binInLowEdgePos == 3 || binInHighEdgePos == 1 ) {
 103                                 //cout<<"In-bin not included."<<endl;
 104                                 PercentageIncluded = 0.;
 105                         }
 106
 107                         // In-bin partially included (left side).
 108                         // In    |    |
 109                         // Out |  |
 110                         else if ( binInLowEdgePos == 2 && binInHighEdgePos == 3 ) {
 111                                 //cout<<"In-bin partially included (left side)."<<endl;
 112                                 PercentageIncluded = (binOutHighEdge - binInLowEdge)/(binInHighEdge - binInLowEdge);
 113                         }
 114
 115                         // In-bin partially included (middle).
 116                         // In    |     |
 117                         // Out     |  |
 118                         else if ( binInLowEdgePos == 1 && binInHighEdgePos == 3 ) {
 119                                 //cout<<"In-bin partially included (middle)."<<endl;
 120                                 PercentageIncluded = (binOutHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
 121                         }
 122
 123                         // In-bin partially included (right side).
 124                         // In    |     |
 125                         // Out        |  |
 126                         else if ( binInLowEdgePos == 1 && binInHighEdgePos == 2 ) {
 127                                 //cout<<"In-bin partially included (right side)."<<endl;
 128                                 PercentageIncluded = (binInHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
 129                         }
 130
 131                         // In-bin completely included.
 132                         // In    |  |
 133                         // Out  |    |
 134                         else if ( binInLowEdgePos == 2 && binInHighEdgePos == 2 ) {
 135                                 //cout<<"In-bin completely included."<<endl;
 136                                 PercentageIncluded = 1.;
 137                         }
 138
 139                         // Give some status:
 140                         //cout<<"Bin Out: "<<iBinOut+1<<" ("<<binOutLowEdge<<","<<binOutHighEdge<<")  Bin In: "<<iBinIn+1<<" ("<<binInLowEdge<<","<<binInHighEdge<<"), content: "<<binInContent<<" +/- "<<TMath::Sqrt(binInError2)<<". Bin in Pos: "<<binInLowEdgePos<<" "<<binInHighEdgePos<<" pct: "<<PercentageIncluded<<endl;
 141
 142                         //else cout<<"Something weird just happened."<<endl;
 143
 144                         Double_t binOutAddedError2 = 0.;
 145                         binOutAddedError2 = binInError2 * PercentageIncluded * PercentageIncluded;
 146                         Double_t binOutAddedContent = 0.;
 147                         binOutAddedContent = binInContent * PercentageIncluded;
 148
 149                         if (density) {
 150                                 binOutAddedContent*=binInWidth;
 151                         }
 152
 153                         //histOut->Fill(binOutCenter,binOutAddedContent);
 154                         histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent);
 155
 156                         //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 157                         //cout<<binOutAddedError2<<endl;
 158                         histOut->SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2);
 159                         //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 160                 }
 161
 162                 // Once a certain binOut has been filled all the way we have to scale
 163                 // it if the "density" flag is set.
 164
 165                 histOut->SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth);
 166                 histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1)));
 167                 //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 168         }
 169
 170         return histOut;
 171
 172 }
 173
 174 // -----------------------------------------------------------------------
 175 TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TH1F* histAxis, Bool_t density) {
 176
 177         // Rebins histogram histIn to the x-axis of histAxis
 178         const TAxis* xaxis = histAxis->GetXaxis();
 179         Int_t nbinsx = xaxis->GetNbins();
 180         const Double_t* binsx = (xaxis->GetXbins())->GetArray();
 181         return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);
 182
 183 }
 184
 185 // -----------------------------------------------------------------------
 186 TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TAxis* xaxis, Bool_t density) {
 187
 188         // Rebins histogram histIn to the x-axis of histAxis
 189         Int_t nbinsx = xaxis->GetNbins();
 190         const Double_t* binsx = (xaxis->GetXbins())->GetArray();
 191         return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);
 192
 193 }
 194
 195 // -----------------------------------------------------------------------
 196 TH1F* AliHistToolsDiHadronPID::TrimHisto(const TH1F* histo, Int_t firstbin, Int_t lastbin) {
 197
 198         const char* name = histo->GetName();
 199         const char* title = histo->GetTitle();
 200         if (firstbin < 0) return 0x0;
 201         if (lastbin > histo->GetNbinsX()) return 0x0;
 202         if (firstbin > lastbin) return 0x0;
 203
 204         Int_t Nbins = lastbin - firstbin + 1;
 205         cout<<firstbin<<" "<<lastbin<<" "<<Nbins<<endl;
 206         //const Int_t NbinsC = Nbins;
 207         Double_t newaxis[41];
 208
 209         for (Int_t ii = firstbin; ii < lastbin+1; ii++) {
 210                 newaxis[ii - firstbin] = histo->GetBinLowEdge(ii);
 211                 //cout<<ii-firstbin<<" "<<newaxis[ii - firstbin]<<endl;
 212         }
 213         newaxis[Nbins] = histo->GetBinLowEdge(lastbin+1);
 214
 215         TH1F* hout = new TH1F("hout","hout",Nbins,newaxis);
 216         hout->SetDirectory(0);
 217         hout->SetTitle(title);
 218         hout->SetName(name);
 219
 220         for (Int_t ii = 0; ii < Nbins+1; ii++) {
 221                 hout->SetBinContent(ii,histo->GetBinContent(firstbin+ii));
 222                 hout->SetBinError(ii,histo->GetBinError(firstbin+ii));
 223         }
 224
 225         return hout;
 226
 227 }
 228
 229 // -----------------------------------------------------------------------
 230 void AliHistToolsDiHadronPID::ConstMinusHist(TH1F* histo, Float_t cc) {
 231
 232         // h -> (c-h)
 233         Int_t nbins = histo->GetNbinsX();
 234         for (Int_t iBin = 0; iBin < (nbins + 1); iBin++) {
 235                 Float_t bincontent = histo->GetBinContent(iBin);
 236                 histo->SetBinContent(iBin,(cc - bincontent));
 237         }
 238
 239 }
 240
 241 // -----------------------------------------------------------------------
 242 TH3F* AliHistToolsDiHadronPID::MakeHist3D(const char* name, const char* title,
 243         Int_t nbinsX, Double_t minX, Double_t maxX,
 244         Int_t nbinsY, Double_t minY, Double_t maxY,
 245         Int_t nbinsZ, const Double_t* zaxis) {
 246
 247         const Double_t* xaxis = const_cast<Double_t*>(CreateAxis(nbinsX,minX,maxX));
 248         const Double_t* yaxis = const_cast<Double_t*>(CreateAxis(nbinsY,minY,maxY));
 249
 250         TH3F* hout = new TH3F(name,title,nbinsX,xaxis,nbinsY,yaxis,nbinsZ,zaxis);
 251
 252         return hout;
 253
 254 }
 255
 256 // -----------------------------------------------------------------------
 257 TH2F* AliHistToolsDiHadronPID::Function2DToHist2D(const TF2* function, const TH2* grid) {
 258
 259         // Creates a 2D histogram of "function" using the binning of
 260         // the histogram "grid".
 261         //  - We assume "grid" to have fixed binwidth.
 262         //  - Bins are only filled if they are within the domain of the function.
 263         //  - Histogram takes over the color of the function.
 264
 265         // Gathering info about the axes.
 266         const TAxis* Xaxis = grid->GetXaxis();
 267         Int_t NbinsX = Xaxis->GetNbins();
 268         const TAxis* Yaxis = grid->GetYaxis();
 269         Int_t NbinsY = Yaxis->GetNbins();
 270
 271         // Determining function range.
 272         Double_t Xmin = 0.;
 273         Double_t Xmax = 0.;
 274         Double_t Ymin = 0.;
 275         Double_t Ymax = 0.;
 276         function->GetRange(Xmin, Ymin, Xmax, Ymax);
 277
 278         // Creating the histogram.
 279         TH2F* hout = new TH2F(Form("%s_hist", function->GetName()),
 280                 Form("%s (Hist);%s;%s", function->GetTitle(), Xaxis->GetTitle(), Yaxis->GetTitle()),
 281                 NbinsX, Xaxis->GetBinLowEdge(1), Xaxis->GetBinUpEdge(NbinsX),
 282                 NbinsY, Yaxis->GetBinLowEdge(1), Yaxis->GetBinUpEdge(NbinsY));
 283         hout->SetLineColor(function->GetLineColor());
 284         hout->SetMarkerColor(function->GetLineColor());
 285         hout->SetDirectory(0);
 286
 287         // Filling the histogram.
 288         for (Int_t iBinX = 1; iBinX < (NbinsX + 1); iBinX++) {
 289                 for (Int_t iBinY = 1; iBinY < (NbinsY + 1); iBinY++) {
 290                         Double_t CoordX = Xaxis->GetBinCenter(iBinX);
 291                         Double_t CoordY = Yaxis->GetBinCenter(iBinY);
 292
 293                         if ( (CoordX < Xmin) || (CoordX > Xmax) || (CoordY < Ymin) || (CoordY > Ymax) ) {continue;}
 294
 295                         hout->SetBinContent(iBinX, iBinY, function->Eval(CoordX, CoordY));
 296                 }
 297         }
 298
 299         return hout;
 300
 301 }
 302
 303 // -----------------------------------------------------------------------
 304 TCanvas* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
 305         const char* title, const TH1F* h1, const TH1F* h2, Int_t markerstyle, Bool_t logy) {
 306
 307         // - Creates a window comparing two histograms h1, and h2.
 308         // - Returns an array of pointers to the objects created
 309         //   in this function.
 310
 311         //Int_t optstat = gStyle->GetOptStat();
 312         gStyle->SetOptStat(0);
 313         //c->UseCurrentStyle();
 314         //gStyle->SetOptStat(optstat);
 315
 316         TH1F* spectra[2];
 317         TH1F* division;
 318         Int_t spectracolor[2] = {1,2};
 319         Int_t divisioncolor = 4;
 320
 321         // Cloning and configuring spectra.
 322         spectra[0] = (TH1F*)h1->Clone();
 323         //spectra[0]->SetDirectory(0);
 324         spectra[1] = (TH1F*)h2->Clone();
 325         //spactra[1]->SetDirectory(0);
 326         for (Int_t iSpectra = 0; iSpectra < 2; iSpectra++) {
 327                 spectra[iSpectra]->Sumw2();
 328                 spectra[iSpectra]->SetMarkerStyle(markerstyle);
 329                 spectra[iSpectra]->SetLineColor(spectracolor[iSpectra]);
 330                 spectra[iSpectra]->SetMarkerColor(spectracolor[iSpectra]);
 331         }
 332
 333         // Creating the division.
 334         division = (TH1F*)spectra[0]->Clone("hDivision");
 335         division->Divide(spectra[1]);
 336         division->SetLineColor(divisioncolor);
 337         division->SetMarkerColor(divisioncolor);
 338         division->GetYaxis()->SetTitle("Ratio");
 339
 340         // Creating window
 341         //TCanvas* c = new TCanvas(name,title,0,0,400,400);
 342         TCanvas* c = TCanvas::MakeDefCanvas();
 343         c->SetName(name);
 344         c->SetTitle(title);
 345         TPad* p1 = new TPad("p1","pad1",0,0.25,1,1);
 346         p1->SetRightMargin(0.05);
 347         p1->SetBottomMargin(0.);
 348         TPad* p2 = new TPad("p2","pad2",0,0,1,0.25);
 349         p2->SetTopMargin(0.);
 350         p2->SetRightMargin(0.05);
 351         p2->SetBottomMargin(0.3);
 352         p2->SetFillStyle(4000);
 353         p1->Draw();
 354         p2->Draw();
 355
 356         // Determining plotting range.
 357         Double_t max = TMath::Max(spectra[0]->GetMaximum(),spectra[1]->GetMaximum());
 358         Double_t min = TMath::Min(spectra[0]->GetMinimum(),spectra[1]->GetMinimum());
 359         Double_t range = max-min;
 360         spectra[0]->SetMinimum(min-0.05*range);
 361         spectra[0]->SetMaximum(max+0.05*range);
 362
 363         // Drawing
 364         p1->cd();
 365         if (logy) {p1->SetLogy();}
 366         spectra[0]->Draw("e");
 367         spectra[0]->GetXaxis()->SetLabelSize(0);
 368         spectra[0]->GetYaxis()->SetLabelSize(0.05);
 369         spectra[0]->GetYaxis()->SetTitleSize(0.05);
 370         spectra[0]->GetYaxis()->SetTitleOffset(0.8);
 371         //Double_t labelsize = spectra[0]->GetYaxis()->GetLabelSize();
 372         spectra[1]->Draw("same e");
 373         p2->cd();
 374         division->SetTitle("");
 375         division->GetXaxis()->SetLabelSize(0.1);
 376         division->GetXaxis()->SetTitleSize(0.15);
 377         division->GetXaxis()->SetTitleOffset(0.6);
 378         division->GetYaxis()->SetLabelSize(0.1);
 379         division->GetYaxis()->SetTitleSize(0.12);
 380         division->GetYaxis()->SetTitleOffset(0.25);
 381         division->Draw("e");
 382
 383         //c->UseCurrentStyle();
 384
 385         //gStyle->SetOptStat(optstat);
 386
 387         TLegend* legend = new TLegend(0.75,0.8,0.95,0.95);
 388         legend->AddEntry(spectra[0],h1->GetTitle(),"lpe");
 389         legend->AddEntry(spectra[1],h2->GetTitle(),"lpe");
 390         p1->cd();
 391         legend->Draw("same");
 392
 393 /*
 394         // returning the created objects.
 395         TObjArray* returnobjects = new TObjArray(6);
 396         returnobjects->AddLast(c);
 397         returnobjects->AddLast(p1);
 398         returnobjects->AddLast(p2);
 399         returnobjects->AddLast(spectra[0]);
 400         returnobjects->AddLast(spectra[1]);
 401         returnobjects->AddLast(division);
 402 */
 403
 404         return c;
 405
 406 }
 407
 408 // -----------------------------------------------------------------------
 409 Double_t* AliHistToolsDiHadronPID::CreateAxis(Int_t nbins, Double_t min, Double_t max) {
 410
 411         if (nbins <= 0) return 0x0;
 412         if (max < min) return 0x0;
 413
 414         Double_t* axis = new Double_t[nbins + 1];
 415         Double_t binsize = (max - min)/((Double_t)nbins);
 416         for (Int_t iBin = 0; iBin < nbins + 1; iBin++) {
 417                 axis[iBin] = min + ((Double_t)iBin) * binsize;
 418         }
 419
 420         return axis;
 421 }
 422
 423
 424
 425
 426
 427