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 "TH1F.h"
  26 #include "TH2F.h"
  27 #include "TH3F.h"
  28 #include "THn.h"
  29 #include "TObjArray.h"
  30 #include "TStyle.h"
  31 #include "TArray.h"
  32
  33 #include "AliHistToolsDiHadronPID.h"
  34
  35 using namespace std;
  36
  37 // -----------------------------------------------------------------------
  38 TObjArray* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
  39         const char* title, TH1F* h1, TH1F* h2, Int_t markerstyle, Bool_t logy) {
  40
  41         // - Creates a window comparing two histograms h1, and h2.
  42         // - Returns an array of pointers to the objects created
  43         //   in this function.
  44
  45         //Int_t optstat = gStyle->GetOptStat();
  46         gStyle->SetOptStat(0);
  47         //c->UseCurrentStyle();
  48         //gStyle->SetOptStat(optstat);
  49
  50         TH1F* spectra[2];
  51         TH1F* division;
  52         Int_t spectracolor[2] = {1,2};
  53         Int_t divisioncolor = 4;
  54
  55         // Cloning and configuring spectra.
  56         spectra[0] = (TH1F*)h1->Clone();
  57         //spectra[0]->SetDirectory(0);
  58         spectra[1] = (TH1F*)h2->Clone();
  59         //spactra[1]->SetDirectory(0);
  60         for (Int_t iSpectra = 0; iSpectra < 2; iSpectra++) {
  61                 spectra[iSpectra]->Sumw2();
  62                 spectra[iSpectra]->SetMarkerStyle(markerstyle);
  63                 spectra[iSpectra]->SetLineColor(spectracolor[iSpectra]);
  64                 spectra[iSpectra]->SetMarkerColor(spectracolor[iSpectra]);
  65         }
  66
  67         // Creating the division.
  68         division = (TH1F*)spectra[0]->Clone("hDivision");
  69         division->Divide(spectra[1]);
  70         division->SetLineColor(divisioncolor);
  71         division->SetMarkerColor(divisioncolor);
  72         division->GetYaxis()->SetTitle("Ratio");
  73
  74         // Creating window
  75         //TCanvas* c = new TCanvas(name,title,0,0,400,400);
  76         TCanvas* c = TCanvas::MakeDefCanvas();
  77         c->SetName(name);
  78         c->SetTitle(title);
  79         TPad* p1 = new TPad("p1","pad1",0,0.25,1,1);
  80         p1->SetRightMargin(0.05);
  81         p1->SetBottomMargin(0.);
  82         TPad* p2 = new TPad("p2","pad2",0,0,1,0.25);
  83         p2->SetTopMargin(0.);
  84         p2->SetRightMargin(0.05);
  85         p2->SetBottomMargin(0.3);
  86         p2->SetFillStyle(4000);
  87         p1->Draw();
  88         p2->Draw();
  89
  90         // Determining plotting range.
  91         Double_t max = TMath::Max(spectra[0]->GetMaximum(),spectra[1]->GetMaximum());
  92         Double_t min = TMath::Min(spectra[0]->GetMinimum(),spectra[1]->GetMinimum());
  93         Double_t range = max-min;
  94         spectra[0]->SetMinimum(min-0.05*range);
  95         spectra[0]->SetMaximum(max+0.05*range);
  96
  97         // Drawing
  98         p1->cd();
  99         if (logy) {p1->SetLogy();}
 100         spectra[0]->Draw("e");
 101         spectra[0]->GetXaxis()->SetLabelSize(0);
 102         spectra[0]->GetYaxis()->SetLabelSize(0.05);
 103         spectra[0]->GetYaxis()->SetTitleSize(0.05);
 104         spectra[0]->GetYaxis()->SetTitleOffset(1.1);
 105         //Double_t labelsize = spectra[0]->GetYaxis()->GetLabelSize();
 106         spectra[1]->Draw("same e");
 107         p2->cd();
 108         division->SetTitle("");
 109         division->GetXaxis()->SetLabelSize(0.1);
 110         division->GetXaxis()->SetTitleSize(0.15);
 111         division->GetXaxis()->SetTitleOffset(0.6);
 112         division->GetYaxis()->SetLabelSize(0.1);
 113         division->GetYaxis()->SetTitleSize(0.12);
 114         division->GetYaxis()->SetTitleOffset(0.3);
 115         division->Draw("e");
 116
 117         //c->UseCurrentStyle();
 118
 119         //gStyle->SetOptStat(optstat);
 120
 121         TLegend* legend = new TLegend(0.75,0.8,0.95,0.95);
 122         legend->AddEntry(spectra[0],h1->GetTitle(),"lpe");
 123         legend->AddEntry(spectra[1],h2->GetTitle(),"lpe");
 124         p1->cd();
 125         legend->Draw("same");
 126
 127         // returning the created objects.
 128         TObjArray* returnobjects = new TObjArray(6);
 129         returnobjects->AddLast(c);
 130         returnobjects->AddLast(p1);
 131         returnobjects->AddLast(p2);
 132         returnobjects->AddLast(spectra[0]);
 133         returnobjects->AddLast(spectra[1]);
 134         returnobjects->AddLast(division);
 135
 136         return returnobjects;
 137
 138 }
 139
 140 // -----------------------------------------------------------------------
 141 TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(
 142         TH1F* histIn, Double_t* binsx, Int_t Nbinsx, Bool_t density) {
 143
 144         // Rebins a histogram (hin) with a variable binning to a histogram
 145         // with another variable binning (binsx). If the "density" flag is set,
 146         // then it is expected that the bins are per unit x-axis, otherwise an
 147         // absolute count is assumed.
 148         //
 149         // TODO: determine over/under-flow bins.
 150
 151         // Gather info from the original histogram.
 152         const TArrayD* binsxIn = histIn->GetXaxis()->GetXbins();
 153         TArrayD* binsxOut = new TArrayD(Nbinsx + 1,binsx);
 154         Int_t NbinsxIn = binsxIn->GetSize() - 1;
 155         const char* nameIn = histIn->GetName();
 156         const char* titleIn = histIn->GetTitle();
 157         const char* xAxisTitleIn = histIn->GetXaxis()->GetTitle();
 158         const char* yAxisTitleIn = histIn->GetYaxis()->GetTitle();
 159
 160         // Gather info from the output histogram and create it.
 161         Int_t NbinsxOut = binsxOut->GetSize() - 1;
 162         const char* nameOut = Form("%s (rebinned)",nameIn);
 163         const char* titleOut = Form("%s_rebinned;%s;%s",titleIn,xAxisTitleIn,yAxisTitleIn);
 164         TH1F* histOut = new TH1F(nameOut,titleOut,NbinsxOut,binsxOut->GetArray());
 165         //histOut->Sumw2();
 166
 167         // Filling the bins.
 168         for (Int_t iBinOut = 0; iBinOut < NbinsxOut; iBinOut++) {
 169
 170                 Double_t binOutLowEdge = binsxOut->At(iBinOut);
 171                 Double_t binOutHighEdge = binsxOut->At(iBinOut+1);
 172                 Double_t binOutWidth = binOutHighEdge - binOutLowEdge;
 173                 //Double_t binOutCenter = binOutHighEdge - binOutWidth/2.;
 174
 175                 // Setting all errors of the outgoing histogram to zero (needed later):
 176                 histOut->SetBinError(iBinOut+1,0.);
 177
 178                 for (Int_t iBinIn = 0; iBinIn < NbinsxIn; iBinIn++) {
 179
 180                         Double_t binInLowEdge = binsxIn->At(iBinIn);
 181                         Double_t binInHighEdge = binsxIn->At(iBinIn+1);
 182                         Double_t binInWidth = binInHighEdge - binInLowEdge;
 183                         //Double_t binInCenter = binInHighEdge - binInWidth/2.;
 184                         Double_t binInContent = histIn->GetBinContent(iBinIn+1);
 185                         Double_t binInError2 = histIn->GetBinError(iBinIn+1)*histIn->GetBinError(iBinIn+1);
 186
 187                         /* -- Determining what percentage of the in-bin is included
 188                            in the current out-bin -- */
 189                         Double_t PercentageIncluded = 0.;
 190
 191                         // Position of the low edge of the in-bin:
 192                         // In  1   2   3
 193                         // Out   |   |
 194                         Int_t binInLowEdgePos = 0;
 195                         if (binInLowEdge < binOutLowEdge) binInLowEdgePos = 1;
 196                         else if (binInLowEdge > binOutHighEdge) binInLowEdgePos = 3;
 197                         else binInLowEdgePos = 2;
 198
 199                         // Same for the high edge of the in-bin:
 200                         Int_t binInHighEdgePos = 0;
 201                         if (binInHighEdge < binOutLowEdge) binInHighEdgePos = 1;
 202                         else if (binInHighEdge > binOutHighEdge) binInHighEdgePos = 3;
 203                         else binInHighEdgePos = 2;
 204
 205                         // In-bin not included.
 206                         // In        | |     or     | |
 207                         // Out |   |                    |    |
 208                         if ( binInLowEdgePos == 3 || binInHighEdgePos == 1 ) {
 209                                 //cout<<"In-bin not included."<<endl;
 210                                 PercentageIncluded = 0.;
 211                         }
 212
 213                         // In-bin partially included (left side).
 214                         // In    |    |
 215                         // Out |  |
 216                         else if ( binInLowEdgePos == 2 && binInHighEdgePos == 3 ) {
 217                                 //cout<<"In-bin partially included (left side)."<<endl;
 218                                 PercentageIncluded = (binOutHighEdge - binInLowEdge)/(binInHighEdge - binInLowEdge);
 219                         }
 220
 221                         // In-bin partially included (middle).
 222                         // In    |     |
 223                         // Out     |  |
 224                         else if ( binInLowEdgePos == 1 && binInHighEdgePos == 3 ) {
 225                                 //cout<<"In-bin partially included (middle)."<<endl;
 226                                 PercentageIncluded = (binOutHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
 227                         }
 228
 229                         // In-bin partially included (right side).
 230                         // In    |     |
 231                         // Out        |  |
 232                         else if ( binInLowEdgePos == 1 && binInHighEdgePos == 2 ) {
 233                                 //cout<<"In-bin partially included (right side)."<<endl;
 234                                 PercentageIncluded = (binInHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
 235                         }
 236
 237                         // In-bin completely included.
 238                         // In    |  |
 239                         // Out  |    |
 240                         else if ( binInLowEdgePos == 2 && binInHighEdgePos == 2 ) {
 241                                 //cout<<"In-bin completely included."<<endl;
 242                                 PercentageIncluded = 1.;
 243                         }
 244
 245                         // Give some status:
 246                         //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;
 247
 248                         //else cout<<"Something weird just happened."<<endl;
 249
 250                         Double_t binOutAddedError2 = 0.;
 251                         binOutAddedError2 = binInError2 * PercentageIncluded * PercentageIncluded;
 252                         Double_t binOutAddedContent = 0.;
 253                         binOutAddedContent = binInContent * PercentageIncluded;
 254
 255                         if (density) {
 256                                 binOutAddedContent*=binInWidth;
 257                         }
 258
 259                         //histOut->Fill(binOutCenter,binOutAddedContent);
 260                         histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent);
 261
 262                         //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 263                         //cout<<binOutAddedError2<<endl;
 264                         histOut->SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2);
 265                         //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 266                 }
 267
 268                 // Once a certain binOut has been filled all the way we have to scale
 269                 // it if the "density" flag is set.
 270
 271                 histOut->SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth);
 272                 histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1)));
 273                 //cout<<histOut->GetBinError(iBinOut+1)<<endl;
 274         }
 275
 276         return histOut;
 277
 278 }
 279
 280 // -----------------------------------------------------------------------
 281 TH1F* AliHistToolsDiHadronPID::TrimHisto(TH1F* histo, Int_t firstbin, Int_t lastbin) {
 282
 283         const char* name = histo->GetName();
 284         const char* title = histo->GetTitle();
 285         if (firstbin < 0) return 0x0;
 286         if (lastbin > histo->GetNbinsX()) return 0x0;
 287         if (firstbin > lastbin) return 0x0;
 288
 289         Int_t Nbins = lastbin - firstbin + 1;
 290         cout<<firstbin<<" "<<lastbin<<" "<<Nbins<<endl;
 291         //const Int_t NbinsC = Nbins;
 292         Double_t newaxis[41];
 293
 294         for (Int_t ii = firstbin; ii < lastbin+1; ii++) {
 295                 newaxis[ii - firstbin] = histo->GetBinLowEdge(ii);
 296                 //cout<<ii-firstbin<<" "<<newaxis[ii - firstbin]<<endl;
 297         }
 298         newaxis[Nbins] = histo->GetBinLowEdge(lastbin+1);
 299
 300         TH1F* hout = new TH1F("hout","hout",Nbins,newaxis);
 301         hout->SetDirectory(0);
 302         hout->SetTitle(title);
 303         hout->SetName(name);
 304
 305         for (Int_t ii = 0; ii < Nbins+1; ii++) {
 306                 hout->SetBinContent(ii,histo->GetBinContent(firstbin+ii));
 307                 hout->SetBinError(ii,histo->GetBinError(firstbin+ii));
 308         }
 309
 310         return hout;
 311
 312 }