Cosmetics
[u/mrichter/AliRoot.git] / 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 }