]>
Commit | Line | Data |
---|---|---|
07d62e30 | 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 | // Definitions the mathematical functions used in the DiHadronPID | |
18 | // analysis. | |
19 | // ----------------------------------------------------------------------- | |
20 | // Author: Misha Veldhoen (misha.veldhoen@cern.ch) | |
21 | ||
a5422983 | 22 | #include "AliFunctionsDiHadronPID.h" |
23 | ||
07d62e30 | 24 | #include <iostream> |
a5422983 | 25 | using namespace std; |
07d62e30 | 26 | |
27 | #include "AliExternalTrackParam.h" | |
07d62e30 | 28 | #include "TF1.h" |
29 | ||
07d62e30 | 30 | // ----------------------------------------------------------------------- |
31 | AliFunctionsDiHadronPID::AliFunctionsDiHadronPID() | |
32 | ||
33 | { | |
34 | ||
35 | // Constructor. | |
36 | ||
37 | } | |
38 | ||
39 | // ----------------------------------------------------------------------- | |
40 | AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID() | |
41 | ||
42 | { | |
43 | ||
44 | // Destructor. | |
45 | ||
46 | } | |
07d62e30 | 47 | |
01f5c9f4 LM |
48 | // ----------------------------------------------------------------------- |
49 | Int_t AliFunctionsDiHadronPID::Power(Int_t base, Int_t power) { | |
07d62e30 | 50 | |
01f5c9f4 | 51 | // Power function for integers (not available in TMath). |
07d62e30 | 52 | |
01f5c9f4 LM |
53 | if (power > 0) { |
54 | Int_t result = 1; | |
55 | for (Int_t ii = 0; ii < power; ++ii) {result *= base;} | |
56 | return result; | |
57 | } else { | |
58 | if (power == 0) {return 1;} | |
59 | else { | |
60 | cout << Form("%s::%s -> WARNING: Method doesn't work for negative powers.",__FILE__,__func__) << endl; | |
61 | return -999; | |
62 | } | |
63 | } | |
07d62e30 | 64 | |
65 | } | |
01f5c9f4 | 66 | |
07d62e30 | 67 | // ----------------------------------------------------------------------- |
fe463f34 | 68 | Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t binwidth) { |
07d62e30 | 69 | |
70 | // The other implementation should make use of this one. | |
71 | Double_t norm = (binwidth*integral)/(TMath::Sqrt(2.*TMath::Pi())*sigma); | |
72 | Double_t gaussian = TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma)); | |
73 | ||
74 | return (norm*gaussian); | |
75 | ||
76 | } | |
07d62e30 | 77 | |
07d62e30 | 78 | // ----------------------------------------------------------------------- |
fe463f34 | 79 | Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t tail, Double_t binwidth) { |
07d62e30 | 80 | |
81 | // Gaussian with exponential tail on the right, I is the integral. | |
82 | // For function definition see: FitFunctions.nb | |
83 | ||
84 | Double_t kappa = mu + tail; | |
85 | ||
86 | if (mu >= kappa) return 0.; // Function becomes ill-defined. | |
87 | ||
88 | Double_t beta = sigma*sigma/(kappa-mu); | |
89 | Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma*sigma) ); | |
90 | Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma*sigma) ); | |
91 | Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma) ); | |
92 | Double_t norm = norm1 + norm2; | |
93 | ||
94 | Double_t funcleft = binwidth * (integral/norm)*TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma)); | |
95 | Double_t funcright = binwidth * (integral/norm)*BB*TMath::Exp(-xx/beta); | |
96 | ||
97 | if (xx <= kappa) return funcleft; | |
98 | else return funcright; | |
99 | ||
100 | } | |
101 | ||
102 | // ----------------------------------------------------------------------- | |
fe463f34 | 103 | Double_t AliFunctionsDiHadronPID::Gaussian2D(Double_t xx, Double_t yy, Double_t integral, |
104 | Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay, | |
105 | Double_t binwidthx, Double_t binwidthy) { | |
07d62e30 | 106 | |
107 | // 2D Gaussian. | |
108 | Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx); | |
109 | Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy); | |
110 | ||
111 | return integral * GaussianX * GaussianY; | |
112 | ||
113 | } | |
114 | ||
115 | // ----------------------------------------------------------------------- | |
fe463f34 | 116 | Double_t AliFunctionsDiHadronPID::Gaussian2DTailX(Double_t xx, Double_t yy, Double_t integral, |
117 | Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay, | |
118 | Double_t tailx, Double_t binwidthx, Double_t binwidthy) { | |
07d62e30 | 119 | |
120 | // 2D Gaussian with exponential tail in X direction. | |
121 | Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx); | |
122 | Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy); | |
123 | ||
124 | return integral * GaussianTailX * GaussianY; | |
125 | ||
126 | } | |
127 | ||
128 | // ----------------------------------------------------------------------- | |
fe463f34 | 129 | Double_t AliFunctionsDiHadronPID::Gaussian2DTailY(Double_t xx, Double_t yy, Double_t integral, |
130 | Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay, | |
131 | Double_t taily, Double_t binwidthx, Double_t binwidthy) { | |
07d62e30 | 132 | |
133 | // 2D Gaussian with exponential tail in Y direction. | |
134 | Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx); | |
135 | Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy); | |
136 | ||
137 | return integral * GaussianX * GaussianTailY; | |
138 | ||
139 | } | |
140 | ||
141 | // ----------------------------------------------------------------------- | |
fe463f34 | 142 | Double_t AliFunctionsDiHadronPID::Gaussian2DTailXY(Double_t xx, Double_t yy, Double_t integral, |
143 | Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay, | |
144 | Double_t tailx, Double_t taily, Double_t binwidthx, Double_t binwidthy) { | |
07d62e30 | 145 | |
146 | // 2D Gaussian with exponential tail in X- and Y direction. | |
147 | Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx); | |
148 | Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy); | |
149 | ||
150 | return integral * GaussianTailX * GaussianTailY; | |
151 | ||
07d62e30 | 152 | } |
153 | ||
07d62e30 | 154 | // ----------------------------------------------------------------------- |
fe463f34 | 155 | Double_t AliFunctionsDiHadronPID::PolyPenalty(Double_t xx, Double_t center, Double_t flatwidth, const Int_t polyorder) { |
07d62e30 | 156 | |
157 | // Penalty function for a chi^2 fit. The function is defined as: | |
158 | // 1 for |xx - center| < flatwidth, | |
159 | // (|xx - center| - flatwidth) ^ polyorder for |xx - center| > flatwidth. | |
160 | ||
161 | Double_t fx = 1.; | |
162 | if (TMath::Abs(xx - center) > flatwidth) { | |
163 | fx = TMath::Power( (TMath::Abs(xx - center) - flatwidth), polyorder ) + 1.; | |
164 | } | |
165 | ||
166 | return fx; | |
167 | ||
168 | } | |
169 | ||
170 | // ----------------------------------------------------------------------- | |
fe463f34 | 171 | TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(Double_t range, Double_t center, Double_t flatwidth, const Int_t polyorder) { |
07d62e30 | 172 | |
173 | // Creates an example of the TestPolyPenalty function. | |
174 | TF1* tf = new TF1("tf",Form("AliFunctionsDiHadronPID::PolyPenalty(x,[0],[1],%i)",polyorder),-range,range); | |
175 | tf->SetParameters(center,flatwidth); | |
176 | TCanvas* cvs = TCanvas::MakeDefCanvas(); | |
177 | tf->Draw(); | |
178 | ||
179 | return cvs; | |
180 | ||
181 | } | |
182 | ||
07d62e30 | 183 | // ----------------------------------------------------------------------- |
fe463f34 | 184 | Double_t AliFunctionsDiHadronPID::TOFExpTime(Double_t pT, Double_t eta, Double_t mass) { |
07d62e30 | 185 | |
186 | // For description see ../Documents/TOFtime.tex | |
187 | ||
188 | Double_t AA = (2. * pT) / ( Charge() * BTPC() * GeVperkg() ); | |
189 | Double_t BB = TMath::ASin( (Charge() * BTPC() * 0.01 * RTOF() * GeVperkg() ) / (2. * pT * C()) ); | |
190 | Double_t CC = TMath::Sqrt( mass*mass/(pT*pT) + TMath::CosH(eta)*TMath::CosH(eta) ); | |
191 | ||
192 | return (1.e12*AA*BB*CC); // Time returned in ps. | |
193 | ||
194 | } | |
195 | ||
196 | // ----------------------------------------------------------------------- | |
fe463f34 | 197 | Double_t AliFunctionsDiHadronPID::TPCExpdEdX(Double_t pT, Double_t eta, Double_t mass) { |
07d62e30 | 198 | |
199 | // Not so neat solution, however the easiest for now. | |
200 | ||
201 | // Prameters taken from the constructor of AliTPCPIDResponse: | |
202 | Double_t MIP = 50.; | |
203 | Double_t Kp[5] = {0.0283086, 2.63394e+01, 5.04114e-11, 2.12543, 4.88663}; | |
204 | ||
205 | Double_t betaGamma = TMath::Abs( (pT * TMath::CosH(eta)) / mass ); | |
206 | ||
207 | // Implementation as in AliTPCPIDResponse. | |
208 | return MIP * AliExternalTrackParam::BetheBlochAleph(betaGamma,Kp[0],Kp[1],Kp[2],Kp[3],Kp[4]); | |
209 | ||
210 | } |