]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/DPhi/DiHadronPID/AliFunctionsDiHadronPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
change binning, remove debug info
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / DiHadronPID / AliFunctionsDiHadronPID.cxx
CommitLineData
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 25using namespace std;
07d62e30 26
27#include "AliExternalTrackParam.h"
07d62e30 28#include "TF1.h"
29
07d62e30 30// -----------------------------------------------------------------------
31AliFunctionsDiHadronPID::AliFunctionsDiHadronPID()
32
33{
34
35 // Constructor.
36
37}
38
39// -----------------------------------------------------------------------
40AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID()
41
42{
43
44 // Destructor.
45
46}
69868b6b 47/*
07d62e30 48// -----------------------------------------------------------------------
fe463f34 49Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t *x, Double_t *par) {
07d62e30 50
51 // - Gaussian, I is the integral -
52 // f(x) = I/(Sqrt(2*pi)*sigma) * exp{-(x-mu)^2/2sigma^2}
53 // par = {I,mu,sigma}
54
55 Double_t norm = par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[2]);
56 Double_t gaussian = TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
57
58 return (norm*gaussian);
59
60}
69868b6b 61*/
07d62e30 62// -----------------------------------------------------------------------
fe463f34 63Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t binwidth) {
07d62e30 64
65 // The other implementation should make use of this one.
66 Double_t norm = (binwidth*integral)/(TMath::Sqrt(2.*TMath::Pi())*sigma);
67 Double_t gaussian = TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
68
69 return (norm*gaussian);
70
71}
69868b6b 72/*
07d62e30 73// -----------------------------------------------------------------------
fe463f34 74Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t *x,const Double_t *par) {
07d62e30 75
76 // Gaussian with exponential tail on the right, I is the integral.
77 // For function definition see: FitFunctions.nb
78
79 Double_t integral = par[0];
80 Double_t mu = par[1]; // Top of the gaussian.
81 Double_t kappa = par[1] + par[2]; // Point at which the gaussian becomes an exponential (w.r.t. to mu).
82 Double_t sigma_x = par[3];
83
84 if (mu >= kappa) return 0.; // Function becomes ill-defined.
85
86 Double_t beta = sigma_x*sigma_x/(kappa-mu);
87 Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma_x*sigma_x) );
88 Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma_x*sigma_x) );
89 Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma_x*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma_x) );
90 Double_t norm = norm1 + norm2;
91
92 Double_t funcleft = (integral/norm)*TMath::Exp(-(x[0]-mu)*(x[0]-mu)/(2.*sigma_x*sigma_x));
93 Double_t funcright = (integral/norm)*BB*TMath::Exp(-x[0]/beta);
94
95 if (x[0] <= kappa) return funcleft;
96 else return funcright;
97
98}
69868b6b 99*/
07d62e30 100// -----------------------------------------------------------------------
fe463f34 101Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t tail, Double_t binwidth) {
07d62e30 102
103 // Gaussian with exponential tail on the right, I is the integral.
104 // For function definition see: FitFunctions.nb
105
106 Double_t kappa = mu + tail;
107
108 if (mu >= kappa) return 0.; // Function becomes ill-defined.
109
110 Double_t beta = sigma*sigma/(kappa-mu);
111 Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma*sigma) );
112 Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma*sigma) );
113 Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma) );
114 Double_t norm = norm1 + norm2;
115
116 Double_t funcleft = binwidth * (integral/norm)*TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
117 Double_t funcright = binwidth * (integral/norm)*BB*TMath::Exp(-xx/beta);
118
119 if (xx <= kappa) return funcleft;
120 else return funcright;
121
122}
123
124// -----------------------------------------------------------------------
fe463f34 125Double_t AliFunctionsDiHadronPID::Gaussian2D(Double_t xx, Double_t yy, Double_t integral,
126 Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
127 Double_t binwidthx, Double_t binwidthy) {
07d62e30 128
129 // 2D Gaussian.
130 Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
131 Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
132
133 return integral * GaussianX * GaussianY;
134
135}
136
137// -----------------------------------------------------------------------
fe463f34 138Double_t AliFunctionsDiHadronPID::Gaussian2DTailX(Double_t xx, Double_t yy, Double_t integral,
139 Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
140 Double_t tailx, Double_t binwidthx, Double_t binwidthy) {
07d62e30 141
142 // 2D Gaussian with exponential tail in X direction.
143 Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
144 Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
145
146 return integral * GaussianTailX * GaussianY;
147
148}
149
150// -----------------------------------------------------------------------
fe463f34 151Double_t AliFunctionsDiHadronPID::Gaussian2DTailY(Double_t xx, Double_t yy, Double_t integral,
152 Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
153 Double_t taily, Double_t binwidthx, Double_t binwidthy) {
07d62e30 154
155 // 2D Gaussian with exponential tail in Y direction.
156 Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
157 Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
158
159 return integral * GaussianX * GaussianTailY;
160
161}
162
163// -----------------------------------------------------------------------
fe463f34 164Double_t AliFunctionsDiHadronPID::Gaussian2DTailXY(Double_t xx, Double_t yy, Double_t integral,
165 Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
166 Double_t tailx, Double_t taily, Double_t binwidthx, Double_t binwidthy) {
07d62e30 167
168 // 2D Gaussian with exponential tail in X- and Y direction.
169 Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
170 Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
171
172 return integral * GaussianTailX * GaussianTailY;
173
174}
69868b6b 175/*
07d62e30 176// -----------------------------------------------------------------------
fe463f34 177Double_t AliFunctionsDiHadronPID::Exponent(Double_t *x, Double_t *par) {
07d62e30 178
179 // f(x) = A * Exp(bx)
180 // par = {A,b}
181
182 return par[0]*TMath::Exp(par[1]*x[0]);
183
184}
185
186// -----------------------------------------------------------------------
187// COMBINED FUNCTIONS
188// -----------------------------------------------------------------------
fe463f34 189Double_t AliFunctionsDiHadronPID::SimpleTOFfit(Double_t *x, Double_t *par) {
07d62e30 190
191 // Signal fitted with a Gaussian, mismatches by an exponent.
192 return (Gaussian1D(x,&par[0]) + Exponent(x,&par[3]));
193
194}
195
196// -----------------------------------------------------------------------
fe463f34 197Double_t AliFunctionsDiHadronPID::SimpleTOFfitWithTail(Double_t *x, Double_t *par) {
07d62e30 198
199 // Signal fitted with a Gaussian with a tail, mismatches by an exponent.
200 return (Gaussian1D(x,&par[0]) + Exponent(x,&par[4]));
201
202}
69868b6b 203*/
07d62e30 204// -----------------------------------------------------------------------
205// PENALTY FUNCTIONS
206// -----------------------------------------------------------------------
fe463f34 207Double_t AliFunctionsDiHadronPID::PolyPenalty(Double_t xx, Double_t center, Double_t flatwidth, const Int_t polyorder) {
07d62e30 208
209 // Penalty function for a chi^2 fit. The function is defined as:
210 // 1 for |xx - center| < flatwidth,
211 // (|xx - center| - flatwidth) ^ polyorder for |xx - center| > flatwidth.
212
213 Double_t fx = 1.;
214 if (TMath::Abs(xx - center) > flatwidth) {
215 fx = TMath::Power( (TMath::Abs(xx - center) - flatwidth), polyorder ) + 1.;
216 }
217
218 return fx;
219
220}
221
222// -----------------------------------------------------------------------
fe463f34 223TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(Double_t range, Double_t center, Double_t flatwidth, const Int_t polyorder) {
07d62e30 224
225 // Creates an example of the TestPolyPenalty function.
226 TF1* tf = new TF1("tf",Form("AliFunctionsDiHadronPID::PolyPenalty(x,[0],[1],%i)",polyorder),-range,range);
227 tf->SetParameters(center,flatwidth);
228 TCanvas* cvs = TCanvas::MakeDefCanvas();
229 tf->Draw();
230
231 return cvs;
232
233}
234
235// -----------------------------------------------------------------------
236// PID Expected signal functions.
237// -----------------------------------------------------------------------
fe463f34 238Double_t AliFunctionsDiHadronPID::TOFExpTime(Double_t pT, Double_t eta, Double_t mass) {
07d62e30 239
240 // For description see ../Documents/TOFtime.tex
241
242 Double_t AA = (2. * pT) / ( Charge() * BTPC() * GeVperkg() );
243 Double_t BB = TMath::ASin( (Charge() * BTPC() * 0.01 * RTOF() * GeVperkg() ) / (2. * pT * C()) );
244 Double_t CC = TMath::Sqrt( mass*mass/(pT*pT) + TMath::CosH(eta)*TMath::CosH(eta) );
245
246 return (1.e12*AA*BB*CC); // Time returned in ps.
247
248}
249
250// -----------------------------------------------------------------------
fe463f34 251Double_t AliFunctionsDiHadronPID::TPCExpdEdX(Double_t pT, Double_t eta, Double_t mass) {
07d62e30 252
253 // Not so neat solution, however the easiest for now.
254
255 // Prameters taken from the constructor of AliTPCPIDResponse:
256 Double_t MIP = 50.;
257 Double_t Kp[5] = {0.0283086, 2.63394e+01, 5.04114e-11, 2.12543, 4.88663};
258
259 Double_t betaGamma = TMath::Abs( (pT * TMath::CosH(eta)) / mass );
260
261 // Implementation as in AliTPCPIDResponse.
262 return MIP * AliExternalTrackParam::BetheBlochAleph(betaGamma,Kp[0],Kp[1],Kp[2],Kp[3],Kp[4]);
263
264}