PWG2/EVCHAR/AliCentralityByFunction.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 /*   Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch                   */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  class to determine centrality percentiles from 2D distributions          //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TNamed.h>
  25 #include <TH1D.h>
  26 #include <TString.h>
  27 #include <TFile.h>
  28 #include <TMath.h>
  29 #include <TROOT.h>
  30 #include <TH2F.h>
  31 #include <TProfile.h>
  32 #include <TF1.h>
  33 #include <TStyle.h>
  34 #include <TGraphErrors.h>
  35 #include <vector>
  36 #include <TMap.h>
  37 #include <map>
  38 #include "AliCentralityByFunction.h"
  39
  40
  41 ClassImp(AliCentralityByFunction)
  42
  43 //______________________________________________________________________________
  44 AliCentralityByFunction::AliCentralityByFunction() :
  45   finrootfile(0),
  46   foutrootfilename(0),
  47   foutrootfile(0),
  48   fhistnames(),
  49   fpercentXsec(0),
  50   fitfunc(),
  51   fitter()
  52 {
  53   // standard constructor
  54   fitter["fitf_pol2"] = new TF1("pol2",0,1,3);
  55   fitter["fitf_pol2"]->SetLineColor(kRed);
  56
  57   fitter["fitf_pol3"] = new TF1("pol3",0,1,4);
  58   fitter["fitf_pol3"]->SetLineColor(kRed);
  59
  60   fitter["fitf_pol4"] = new TF1("pol4",0,1,5);
  61   fitter["fitf_pol4"]->SetLineColor(kRed);
  62
  63   fitter["fitf_pol6"] = new TF1("pol6",0,1,6);
  64   fitter["fitf_pol6"]->SetLineColor(kRed);
  65 }
  66
  67 void AliCentralityByFunction::SetFitFunction(TString distribution, TString func, Double_t xmin, Double_t xmax)
  68 {
  69   // Set fit function
  70   fitfunc[distribution] = func;
  71   fitter[fitfunc[distribution]]->SetRange(xmin,xmax);
  72 }
  73
  74 void AliCentralityByFunction::MakePercentiles(TString infilename)
  75 {
  76   // Make percentile bins
  77
  78   TH1D *hpercentile;
  79
  80   // open inrootfile, outrootfile
  81   finrootfile  = new TFile(infilename);
  82   foutrootfile = new TFile(foutrootfilename,"RECREATE");
  83
  84   // loop over all distribution names
  85   std::vector<TString>::const_iterator hni;
  86   for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
  87     hpercentile = FitHisto(*hni);
  88     foutrootfile->cd();
  89     hpercentile->Write();
  90   }
  91   // close inrootfile, outrootfile
  92   finrootfile->Close();
  93   foutrootfile->Close();
  94 }
  95
  96 TH1D *AliCentralityByFunction::FitHisto(TString hdistributionName)
  97 {
  98   // Fit histogram
  99   TH2D *hdist  = (TH2D*) (finrootfile->Get(hdistributionName));
 100   TProfile *profile =hdist->ProfileX();
 101   //  fitter[fitfunc[hdistributionName]]->SetRange(0,profile->GetBinCenter(profile->GetNbinsX()));
 102   profile->Fit(fitter[fitfunc[hdistributionName]], "RNM");
 103
 104   foutrootfile->cd();
 105   profile->Write();
 106   fitter[fitfunc[hdistributionName]]->Write(hdistributionName.Append("_fit"));
 107
 108   return MakePercentHisto(hdist);
 109 }
 110
 111 TH1D * AliCentralityByFunction::MakePercentHisto(TH2D *histo)
 112 {
 113   TH2D *htemp = (TH2D*)histo->Clone("htemp");
 114   TString hdistributionName = htemp->GetTitle();
 115
 116   const int num_DIVISION=500;
 117
 118   TH1D *hpercent  = new TH1D("","",num_DIVISION,0,htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX()));
 119   hpercent->Reset();
 120
 121   double xpoint, ypoint, xcomp, ycomp, count;
 122
 123   double xmax = htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX());
 124   double xmin = 0;
 125
 126   double delta = (xmax - xmin) / num_DIVISION;
 127   double slopePerp;   // slope of the perpendicular line
 128   double slopeFunction;
 129
 130   std::cout << "Start Percentile Histo for distribution " << hdistributionName << " with fit function " << fitfunc[hdistributionName] << std::endl;
 131
 132   // move perpendicular line along fitting curve
 133   // and count number of points on the right
 134   for (int i = 0; i <= num_DIVISION; i++) {
 135     xcomp = xmin + i * delta;
 136     ycomp = fitter[fitfunc[hdistributionName]]->Eval(xcomp);
 137     count = 0.0;
 138     slopeFunction = fitter[fitfunc[hdistributionName]]->Derivative(xcomp,NULL,0.0001);
 139     slopePerp = -1.0 /slopeFunction;
 140     //
 141     // equation of perpendicular line
 142     // (y-Ycomp)/(x-Xcomp) = slopePerp
 143     //
 144     for (int ibiny = 1; ibiny < htemp->GetNbinsY(); ibiny++) {
 145       ypoint = htemp->GetYaxis()->GetBinCenter(ibiny);
 146
 147       double xonLine =(ypoint - ycomp)/slopePerp + xcomp;
 148       //      double YonLine = slopePerp*(XonLine - Xcomp) +Ycomp;
 149
 150       for (int ibinx = 1; ibinx < htemp->GetNbinsX(); ibinx++) {
 151         xpoint = htemp->GetXaxis()->GetBinCenter(ibinx);
 152         //
 153         // (XonLine,Ypoint) lies on the perpendicular
 154         //
 155          if ((xpoint > xonLine && xpoint > xcomp - 0.5)||
 156              (xpoint > xcomp + 0.2)) {
 157
 158           count  += htemp->GetBinContent(ibinx,ibiny);
 159         }
 160       }
 161     }
 162     count = count/htemp->Integral() * 100.0;   // change to percentage
 163     hpercent->SetBinContent(i, count);
 164   }
 165
 166   hpercent->SetName(hdistributionName.Append("_percentile"));
 167
 168   return hpercent;
 169 }