1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */
18 ///////////////////////////////////////////////////////////////////////////////
20 // class to determine centrality percentiles from 2D distributions //
22 ///////////////////////////////////////////////////////////////////////////////
34 #include <TGraphErrors.h>
38 #include "AliCentralityByFunction.h"
41 ClassImp(AliCentralityByFunction)
43 //______________________________________________________________________________
44 AliCentralityByFunction::AliCentralityByFunction() :
53 // standard constructor
54 fitter["fitf_pol2"] = new TF1("pol2",0,1,3);
55 fitter["fitf_pol2"]->SetLineColor(kRed);
57 fitter["fitf_pol3"] = new TF1("pol3",0,1,4);
58 fitter["fitf_pol3"]->SetLineColor(kRed);
60 fitter["fitf_pol4"] = new TF1("pol4",0,1,5);
61 fitter["fitf_pol4"]->SetLineColor(kRed);
63 fitter["fitf_pol6"] = new TF1("pol6",0,1,6);
64 fitter["fitf_pol6"]->SetLineColor(kRed);
67 void AliCentralityByFunction::SetFitFunction(TString distribution, TString func, Double_t xmin, Double_t xmax)
70 fitfunc[distribution] = func;
71 fitter[fitfunc[distribution]]->SetRange(xmin,xmax);
74 void AliCentralityByFunction::MakePercentiles(TString infilename)
76 // Make percentile bins
80 // open inrootfile, outrootfile
81 finrootfile = new TFile(infilename);
82 foutrootfile = new TFile(foutrootfilename,"RECREATE");
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);
91 // close inrootfile, outrootfile
93 foutrootfile->Close();
96 TH1D *AliCentralityByFunction::FitHisto(TString hdistributionName)
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");
106 fitter[fitfunc[hdistributionName]]->Write(hdistributionName.Append("_fit"));
108 return MakePercentHisto(hdist);
111 TH1D * AliCentralityByFunction::MakePercentHisto(TH2D *histo)
113 TH2D *htemp = (TH2D*)histo->Clone("htemp");
114 TString hdistributionName = htemp->GetTitle();
116 const int num_DIVISION=500;
118 TH1D *hpercent = new TH1D("","",num_DIVISION,0,htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX()));
121 double xpoint, ypoint, xcomp, ycomp, count;
123 double xmax = htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX());
126 double delta = (xmax - xmin) / num_DIVISION;
127 double slopePerp; // slope of the perpendicular line
128 double slopeFunction;
130 std::cout << "Start Percentile Histo for distribution " << hdistributionName << " with fit function " << fitfunc[hdistributionName] << std::endl;
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);
138 slopeFunction = fitter[fitfunc[hdistributionName]]->Derivative(xcomp,NULL,0.0001);
139 slopePerp = -1.0 /slopeFunction;
141 // equation of perpendicular line
142 // (y-Ycomp)/(x-Xcomp) = slopePerp
144 for (int ibiny = 1; ibiny < htemp->GetNbinsY(); ibiny++) {
145 ypoint = htemp->GetYaxis()->GetBinCenter(ibiny);
147 double xonLine =(ypoint - ycomp)/slopePerp + xcomp;
148 // double YonLine = slopePerp*(XonLine - Xcomp) +Ycomp;
150 for (int ibinx = 1; ibinx < htemp->GetNbinsX(); ibinx++) {
151 xpoint = htemp->GetXaxis()->GetBinCenter(ibinx);
153 // (XonLine,Ypoint) lies on the perpendicular
155 if ((xpoint > xonLine && xpoint > xcomp - 0.5)||
156 (xpoint > xcomp + 0.2)) {
158 count += htemp->GetBinContent(ibinx,ibiny);
162 count = count/htemp->Integral() * 100.0; // change to percentage
163 hpercent->SetBinContent(i, count);
166 hpercent->SetName(hdistributionName.Append("_percentile"));