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 "AliFunctionsDiHadronPID.h"
  23
  24 #include <iostream>
  25 using namespace std;
  26
  27 #include "AliExternalTrackParam.h"
  28 #include "TF1.h"
  29
  30 // -----------------------------------------------------------------------
  31 AliFunctionsDiHadronPID::AliFunctionsDiHadronPID()
  32
  33 {
  34
  35         // Constructor.
  36
  37 }
  38
  39 // -----------------------------------------------------------------------
  40 AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID()
  41
  42 {
  43
  44         // Destructor.
  45
  46 }
  47 /*
  48 // -----------------------------------------------------------------------
  49 Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t *x, Double_t *par) {
  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 }
  61 */
  62 // -----------------------------------------------------------------------
  63 Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t binwidth) {
  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 }
  72 /*
  73 // -----------------------------------------------------------------------
  74 Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t *x,const  Double_t *par) {
  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 }
  99 */
 100 // -----------------------------------------------------------------------
 101 Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t tail, Double_t binwidth) {
 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 // -----------------------------------------------------------------------
 125 Double_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) {
 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 // -----------------------------------------------------------------------
 138 Double_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) {
 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 // -----------------------------------------------------------------------
 151 Double_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) {
 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 // -----------------------------------------------------------------------
 164 Double_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) {
 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 }
 175 /*
 176 // -----------------------------------------------------------------------
 177 Double_t AliFunctionsDiHadronPID::Exponent(Double_t *x, Double_t *par) {
 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 // -----------------------------------------------------------------------
 189 Double_t AliFunctionsDiHadronPID::SimpleTOFfit(Double_t *x, Double_t *par) {
 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 // -----------------------------------------------------------------------
 197 Double_t AliFunctionsDiHadronPID::SimpleTOFfitWithTail(Double_t *x, Double_t *par) {
 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 }
 203 */
 204 // -----------------------------------------------------------------------
 205 //  PENALTY FUNCTIONS
 206 // -----------------------------------------------------------------------
 207 Double_t AliFunctionsDiHadronPID::PolyPenalty(Double_t xx, Double_t center, Double_t flatwidth, const Int_t polyorder) {
 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 // -----------------------------------------------------------------------
 223 TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(Double_t range, Double_t center, Double_t flatwidth, const Int_t polyorder) {
 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 // -----------------------------------------------------------------------
 238 Double_t AliFunctionsDiHadronPID::TOFExpTime(Double_t pT, Double_t eta, Double_t mass) {
 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 // -----------------------------------------------------------------------
 251 Double_t AliFunctionsDiHadronPID::TPCExpdEdX(Double_t pT, Double_t eta, Double_t mass) {
 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 }