+/**************************************************************************
+ * Copyright(c) 1998-1999, 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. *
+ **************************************************************************/
+
+#include "AliTRDPartID.h"
+#include "AliESDtrack.h"
+#include "TProfile.h"
+#include "TF1.h"
+
+#include <TMath.h>
+
+ClassImp(AliTRDPartID)
+
+
+AliTRDPartID::AliTRDPartID()
+{
+// create a TRD particle identifier
+
+ fBetheBloch = NULL;
+ fRes = 0.2;
+ fRange = 3.;
+}
+
+AliTRDPartID::AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range)
+{
+// create a TRD particle identifier with custom settings
+
+ fBetheBloch = betheBloch;
+ fRes = res;
+ fRange = range;
+}
+
+AliTRDPartID::~AliTRDPartID()
+{
+ if (fBetheBloch) delete fBetheBloch;
+}
+
+
+Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
+{
+// This function calculates the "detector response" PID probabilities
+
+ static const Double_t masses[]={
+ 0.000511, 0.105658, 0.139570, 0.493677, 0.938272
+ };
+
+ if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
+ ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
+ Double_t momentum = track->GetP();
+ if (momentum < 0.001) return kFALSE;
+
+ // get the probability densities
+ Double_t pSum = 0;
+ for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
+ Double_t expectedSignal = fBetheBloch->Eval(momentum/masses[iSpecies]);
+ Double_t expectedError = fRes * expectedSignal;
+ Double_t measuredSignal = track->GetTRDsignal();
+ if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
+ track->SetTRDpid(iSpecies, 0.);
+ } else {
+ Double_t p = TMath::Gaus(measuredSignal, expectedSignal, expectedError);
+ pSum += p;
+ track->SetTRDpid(iSpecies, p);
+ }
+ }
+
+ // "normalize" the probability densities
+ if (pSum <= 0) return kFALSE;
+ for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
+ track->SetTRDpid(iSpecies, track->GetTRDpid(iSpecies) / pSum);
+ }
+
+ return kTRUE;
+}
+
+
+void AliTRDPartID::FitBetheBloch(TProfile* dEdxVsBetaGamma)
+{
+// determine the parameters of the bethe bloch function
+
+ if (fBetheBloch) delete fBetheBloch;
+ fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch,
+ 0.001, 100000., 5);
+ fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
+ fBetheBloch->SetParLimits(2, 0., 0.01);
+ fBetheBloch->SetParLimits(3, 0., 10.);
+ fBetheBloch->SetParLimits(4, 0., 10.);
+ fBetheBloch->SetFillStyle(0);
+ fBetheBloch->SetLineColor(kRed);
+ dEdxVsBetaGamma->Fit(fBetheBloch, "N");
+}
+
+TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
+{
+// create a function for expected dE/dx vs p
+
+ TF1* result = new TF1("betheBlochMass", fcnBetheBlochMass,
+ 0.001, 100000., 6);
+ result->SetParameter(0, mass);
+ if (fBetheBloch) {
+ for (Int_t iPar = 0; iPar < 5; iPar++) {
+ result->SetParameter(iPar+1, fBetheBloch->GetParameter(iPar));
+ result->SetParError(iPar+1, fBetheBloch->GetParError(iPar));
+ }
+ }
+ result->SetFillStyle(0);
+ return result;
+}
+
+Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
+{
+// parametrized bethe bloch function (xx[0] = beta*gamma = p/m):
+//
+// p0/beta^p3 * [ p1 - log(p2 + 1/(beta*gamma)^p4) - beta^p3 ]
+//
+
+ Double_t betaGamma2 = xx[0] * xx[0];
+ Double_t beta2 = betaGamma2 / (1. + betaGamma2);
+ Double_t betaPar3 = TMath::Power(beta2, par[3]/2.);
+ return par[0]/betaPar3 * (par[1] - TMath::Log(TMath::Abs(par[2] + TMath::Power(betaGamma2, -par[4]/2.))) - betaPar3);
+}
+
+Double_t AliTRDPartID::fcnBetheBlochMass(Double_t* xx, Double_t* par)
+{
+// parametrized bethe bloch function: xx[0] = p, p0 = mass
+
+ Double_t betaGamma = xx[0] / par[0];
+ return fcnBetheBloch(&betaGamma, &par[1]);
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff