+++ /dev/null
-/*************************************************************************
-* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
-* *
-* Author: The ALICE Off-line Project. *
-* Contributors are mentioned in the code where appropriate. *
-* *
-* Permission to use, copy, modify and distribute this software and its *
-* documentation strictly for non-commercial purposes is hereby granted *
-* without fee, provided that the above copyright notice appears in all *
-* copies and that both the copyright notice and this permission notice *
-* appear in the supporting documentation. The authors make no claims *
-* about the suitability of this software for any purpose. It is *
-* provided "as is" without express or implied warranty. *
-**************************************************************************/
-
-// -----------------------------------------------------------------------
-// Definitions the mathematical functions used in the DiHadronPID
-// analysis.
-// -----------------------------------------------------------------------
-// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
-
-#include "AliFunctionsDiHadronPID.h"
-
-#include <iostream>
-using namespace std;
-
-#include "AliExternalTrackParam.h"
-#include "TF1.h"
-
-// -----------------------------------------------------------------------
-AliFunctionsDiHadronPID::AliFunctionsDiHadronPID()
-
-{
-
- // Constructor.
-
-}
-
-// -----------------------------------------------------------------------
-AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID()
-
-{
-
- // Destructor.
-
-}
-
-// -----------------------------------------------------------------------
-Int_t AliFunctionsDiHadronPID::Power(Int_t base, Int_t power) {
-
- // Power function for integers (not available in TMath).
-
- if (power > 0) {
- Int_t result = 1;
- for (Int_t ii = 0; ii < power; ++ii) {result *= base;}
- return result;
- } else {
- if (power == 0) {return 1;}
- else {
- cout << Form("%s::%s -> WARNING: Method doesn't work for negative powers.",__FILE__,__func__) << endl;
- return -999;
- }
- }
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian1D(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t binwidth) {
-
- // The other implementation should make use of this one.
- Double_t norm = (binwidth*integral)/(TMath::Sqrt(2.*TMath::Pi())*sigma);
- Double_t gaussian = TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
-
- return (norm*gaussian);
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian1DTail(Double_t xx, Double_t integral, Double_t mu, Double_t sigma, Double_t tail, Double_t binwidth) {
-
- // Gaussian with exponential tail on the right, I is the integral.
- // For function definition see: FitFunctions.nb
-
- Double_t kappa = mu + tail;
-
- if (mu >= kappa) return 0.; // Function becomes ill-defined.
-
- Double_t beta = sigma*sigma/(kappa-mu);
- Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma*sigma) );
- Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma*sigma) );
- Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma) );
- Double_t norm = norm1 + norm2;
-
- Double_t funcleft = binwidth * (integral/norm)*TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
- Double_t funcright = binwidth * (integral/norm)*BB*TMath::Exp(-xx/beta);
-
- if (xx <= kappa) return funcleft;
- else return funcright;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian2D(Double_t xx, Double_t yy, Double_t integral,
- Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
- Double_t binwidthx, Double_t binwidthy) {
-
- // 2D Gaussian.
- Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
- Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
-
- return integral * GaussianX * GaussianY;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailX(Double_t xx, Double_t yy, Double_t integral,
- Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
- Double_t tailx, Double_t binwidthx, Double_t binwidthy) {
-
- // 2D Gaussian with exponential tail in X direction.
- Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
- Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
-
- return integral * GaussianTailX * GaussianY;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailY(Double_t xx, Double_t yy, Double_t integral,
- Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
- Double_t taily, Double_t binwidthx, Double_t binwidthy) {
-
- // 2D Gaussian with exponential tail in Y direction.
- Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
- Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
-
- return integral * GaussianX * GaussianTailY;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailXY(Double_t xx, Double_t yy, Double_t integral,
- Double_t mux, Double_t muy, Double_t sigmax, Double_t sigmay,
- Double_t tailx, Double_t taily, Double_t binwidthx, Double_t binwidthy) {
-
- // 2D Gaussian with exponential tail in X- and Y direction.
- Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
- Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
-
- return integral * GaussianTailX * GaussianTailY;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::PolyPenalty(Double_t xx, Double_t center, Double_t flatwidth, const Int_t polyorder) {
-
- // Penalty function for a chi^2 fit. The function is defined as:
- // 1 for |xx - center| < flatwidth,
- // (|xx - center| - flatwidth) ^ polyorder for |xx - center| > flatwidth.
-
- Double_t fx = 1.;
- if (TMath::Abs(xx - center) > flatwidth) {
- fx = TMath::Power( (TMath::Abs(xx - center) - flatwidth), polyorder ) + 1.;
- }
-
- return fx;
-
-}
-
-// -----------------------------------------------------------------------
-TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(Double_t range, Double_t center, Double_t flatwidth, const Int_t polyorder) {
-
- // Creates an example of the TestPolyPenalty function.
- TF1* tf = new TF1("tf",Form("AliFunctionsDiHadronPID::PolyPenalty(x,[0],[1],%i)",polyorder),-range,range);
- tf->SetParameters(center,flatwidth);
- TCanvas* cvs = TCanvas::MakeDefCanvas();
- tf->Draw();
-
- return cvs;
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::TOFExpTime(Double_t pT, Double_t eta, Double_t mass) {
-
- // For description see ../Documents/TOFtime.tex
-
- Double_t AA = (2. * pT) / ( Charge() * BTPC() * GeVperkg() );
- Double_t BB = TMath::ASin( (Charge() * BTPC() * 0.01 * RTOF() * GeVperkg() ) / (2. * pT * C()) );
- Double_t CC = TMath::Sqrt( mass*mass/(pT*pT) + TMath::CosH(eta)*TMath::CosH(eta) );
-
- return (1.e12*AA*BB*CC); // Time returned in ps.
-
-}
-
-// -----------------------------------------------------------------------
-Double_t AliFunctionsDiHadronPID::TPCExpdEdX(Double_t pT, Double_t eta, Double_t mass) {
-
- // Not so neat solution, however the easiest for now.
-
- // Prameters taken from the constructor of AliTPCPIDResponse:
- Double_t MIP = 50.;
- Double_t Kp[5] = {0.0283086, 2.63394e+01, 5.04114e-11, 2.12543, 4.88663};
-
- Double_t betaGamma = TMath::Abs( (pT * TMath::CosH(eta)) / mass );
-
- // Implementation as in AliTPCPIDResponse.
- return MIP * AliExternalTrackParam::BetheBlochAleph(betaGamma,Kp[0],Kp[1],Kp[2],Kp[3],Kp[4]);
-
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff