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