PWGCF/Correlations/DPhi/DiHadronPID/AliFunctionsDiHadronPID.cxx

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