]>
Commit | Line | Data |
---|---|---|
e6f3f2fe | 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 | //______________________________________________________________________________ | |
bdd08257 | 44 | AliCentralityByFunction::AliCentralityByFunction() : |
45 | finrootfile(0), | |
46 | foutrootfilename(0), | |
47 | foutrootfile(0), | |
48 | fhistnames(), | |
49 | fpercentXsec(0), | |
50 | fitfunc(), | |
51 | fitter() | |
52 | { | |
e6f3f2fe | 53 | // standard constructor |
87dea6ac | 54 | fitter["fitf_pol2"] = new TF1("pol2",0,1,3); |
e6f3f2fe | 55 | fitter["fitf_pol2"]->SetLineColor(kRed); |
56 | ||
87dea6ac | 57 | fitter["fitf_pol3"] = new TF1("pol3",0,1,4); |
e6f3f2fe | 58 | fitter["fitf_pol3"]->SetLineColor(kRed); |
59 | ||
87dea6ac | 60 | fitter["fitf_pol4"] = new TF1("pol4",0,1,5); |
e6f3f2fe | 61 | fitter["fitf_pol4"]->SetLineColor(kRed); |
62 | ||
87dea6ac | 63 | fitter["fitf_pol6"] = new TF1("pol6",0,1,6); |
e6f3f2fe | 64 | fitter["fitf_pol6"]->SetLineColor(kRed); |
65 | } | |
66 | ||
bdd08257 | 67 | void AliCentralityByFunction::SetFitFunction(TString distribution, TString func, Double_t xmin, Double_t xmax) |
68 | { | |
69 | // Set fit function | |
e6f3f2fe | 70 | fitfunc[distribution] = func; |
71 | fitter[fitfunc[distribution]]->SetRange(xmin,xmax); | |
72 | } | |
73 | ||
bdd08257 | 74 | void AliCentralityByFunction::MakePercentiles(TString infilename) |
75 | { | |
76 | // Make percentile bins | |
77 | ||
e6f3f2fe | 78 | TH1D *hpercentile; |
79 | ||
80 | // open inrootfile, outrootfile | |
bdd08257 | 81 | finrootfile = new TFile(infilename); |
82 | foutrootfile = new TFile(foutrootfilename,"RECREATE"); | |
e6f3f2fe | 83 | |
84 | // loop over all distribution names | |
bdd08257 | 85 | std::vector<TString>::const_iterator hni; |
86 | for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) { | |
e6f3f2fe | 87 | hpercentile = FitHisto(*hni); |
bdd08257 | 88 | foutrootfile->cd(); |
e6f3f2fe | 89 | hpercentile->Write(); |
90 | } | |
91 | // close inrootfile, outrootfile | |
bdd08257 | 92 | finrootfile->Close(); |
93 | foutrootfile->Close(); | |
e6f3f2fe | 94 | } |
95 | ||
bdd08257 | 96 | TH1D *AliCentralityByFunction::FitHisto(TString hdistributionName) |
97 | { | |
98 | // Fit histogram | |
99 | TH2D *hdist = (TH2D*) (finrootfile->Get(hdistributionName)); | |
e6f3f2fe | 100 | TProfile *profile =hdist->ProfileX(); |
101 | // fitter[fitfunc[hdistributionName]]->SetRange(0,profile->GetBinCenter(profile->GetNbinsX())); | |
102 | profile->Fit(fitter[fitfunc[hdistributionName]], "RNM"); | |
103 | ||
bdd08257 | 104 | foutrootfile->cd(); |
e6f3f2fe | 105 | profile->Write(); |
106 | fitter[fitfunc[hdistributionName]]->Write(hdistributionName.Append("_fit")); | |
107 | ||
108 | return MakePercentHisto(hdist); | |
109 | } | |
110 | ||
bdd08257 | 111 | TH1D * AliCentralityByFunction::MakePercentHisto(TH2D *histo) |
112 | { | |
e6f3f2fe | 113 | TH2D *htemp = (TH2D*)histo->Clone("htemp"); |
114 | TString hdistributionName = htemp->GetTitle(); | |
115 | ||
bdd08257 | 116 | const int num_DIVISION=500; |
e6f3f2fe | 117 | |
bdd08257 | 118 | TH1D *hpercent = new TH1D("","",num_DIVISION,0,htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX())); |
e6f3f2fe | 119 | hpercent->Reset(); |
120 | ||
bdd08257 | 121 | double xpoint, ypoint, xcomp, ycomp, count; |
e6f3f2fe | 122 | |
bdd08257 | 123 | double xmax = htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX()); |
124 | double xmin = 0; | |
e6f3f2fe | 125 | |
bdd08257 | 126 | double delta = (xmax - xmin) / num_DIVISION; |
e6f3f2fe | 127 | double slopePerp; // slope of the perpendicular line |
128 | double slopeFunction; | |
129 | ||
bdd08257 | 130 | std::cout << "Start Percentile Histo for distribution " << hdistributionName << " with fit function " << fitfunc[hdistributionName] << std::endl; |
e6f3f2fe | 131 | |
132 | // move perpendicular line along fitting curve | |
133 | // and count number of points on the right | |
bdd08257 | 134 | for (int i = 0; i <= num_DIVISION; i++) { |
135 | xcomp = xmin + i * delta; | |
136 | ycomp = fitter[fitfunc[hdistributionName]]->Eval(xcomp); | |
e6f3f2fe | 137 | count = 0.0; |
bdd08257 | 138 | slopeFunction = fitter[fitfunc[hdistributionName]]->Derivative(xcomp,NULL,0.0001); |
e6f3f2fe | 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++) { | |
bdd08257 | 145 | ypoint = htemp->GetYaxis()->GetBinCenter(ibiny); |
e6f3f2fe | 146 | |
bdd08257 | 147 | double xonLine =(ypoint - ycomp)/slopePerp + xcomp; |
e6f3f2fe | 148 | // double YonLine = slopePerp*(XonLine - Xcomp) +Ycomp; |
149 | ||
150 | for (int ibinx = 1; ibinx < htemp->GetNbinsX(); ibinx++) { | |
bdd08257 | 151 | xpoint = htemp->GetXaxis()->GetBinCenter(ibinx); |
e6f3f2fe | 152 | // |
153 | // (XonLine,Ypoint) lies on the perpendicular | |
154 | // | |
bdd08257 | 155 | if ((xpoint > xonLine && xpoint > xcomp - 0.5)|| |
156 | (xpoint > xcomp + 0.2)) { | |
e6f3f2fe | 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 | } |