X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TRD%2FAliTRDpidESD.cxx;h=705103a4d6fb34831dec6aaa47038aa3335d54b5;hp=bb017537b7f7e3728bb58c8991da8f732e4cfd42;hb=93daf3f7e3e7c2da5d535ac201dbdc6243d33467;hpb=304864ab927d252f20e224fae51913bbccd76283 diff --git a/TRD/AliTRDpidESD.cxx b/TRD/AliTRDpidESD.cxx index bb017537b7f..705103a4d6f 100644 --- a/TRD/AliTRDpidESD.cxx +++ b/TRD/AliTRDpidESD.cxx @@ -1,111 +1,301 @@ /************************************************************************** - * 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. * - **************************************************************************/ - -//----------------------------------------------------------------- -// Implementation of the TRD PID class -// Very naive one... And the implementation is even poorer... -// Should be made better by the detector experts... -//----------------------------------------------------------------- +* 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 "AliTRDpidESD.h" -#include "AliESD.h" +/* $Id$ */ + +//////////////////////////////////////////////////////////////////////////// +// // +// Implementation of the TRD PID class // +// // +// Assigns the electron and pion likelihoods to each ESD track. // +// The function MakePID(AliESDEvent *event) calculates the probability // +// of having dedx and a maximum timbin at a given // +// momentum (mom) and particle type k // +// from the precalculated distributions. // +// // +// Authors : // +// Prashant Shukla (orig. version) // +// Alex Bercuci (a.bercuci@gsi.de) // +// // +//////////////////////////////////////////////////////////////////////////// + +#include "AliESDEvent.h" #include "AliESDtrack.h" +#include "AliTracker.h" + +#include "AliTRDpidESD.h" +#include "AliTRDgeometry.h" +#include "AliTRDcalibDB.h" +#include "Cal/AliTRDCalPID.h" ClassImp(AliTRDpidESD) -//_________________________________________________________________________ -AliTRDpidESD::AliTRDpidESD(Double_t *param) +Bool_t AliTRDpidESD::fgCheckTrackStatus = kTRUE; +Bool_t AliTRDpidESD::fgCheckKinkStatus = kFALSE; +Int_t AliTRDpidESD::fgMinPlane = 0; + +//_____________________________________________________________________________ +AliTRDpidESD::AliTRDpidESD() + :TObject() + ,fTrack(NULL) { // - // The main constructor + // Default constructor // - fMIP=param[0]; // MIP signal - fRes=param[1]; // relative resolution - fRange=param[2]; // PID "range" (in sigmas) + } -Double_t AliTRDpidESD::Bethe(Double_t bg) +//_____________________________________________________________________________ +AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) + :TObject(p) + ,fTrack(NULL) { // - // Parametrization of the Bethe-Bloch-curve - // The parametrization is the same as for the TPC and is taken from Lehrhaus. - // - - // This parameters have been adjusted to averaged values from GEANT - const Double_t kP1 = 7.17960e-02; - const Double_t kP2 = 8.54196; - const Double_t kP3 = 1.38065e-06; - const Double_t kP4 = 5.30972; - const Double_t kP5 = 2.83798; - - // This parameters have been adjusted to Xe-data found in: - // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253 - //const Double_t kP1 = 0.76176E-1; - //const Double_t kP2 = 10.632; - //const Double_t kP3 = 3.17983E-6; - //const Double_t kP4 = 1.8631; - //const Double_t kP5 = 1.9479; - - // Lower cutoff of the Bethe-Bloch-curve to limit step sizes - const Double_t kBgMin = 0.8; - const Double_t kBBMax = 6.83298; - //const Double_t kBgMin = 0.6; - //const Double_t kBBMax = 17.2809; - //const Double_t kBgMin = 0.4; - //const Double_t kBBMax = 82.0; - - if (bg > kBgMin) { - Double_t yy = bg / TMath::Sqrt(1. + bg*bg); - Double_t aa = TMath::Power(yy,kP4); - Double_t bb = TMath::Power((1./bg),kP5); - bb = TMath::Log(kP3 + bb); - return ((kP2 - aa - bb)*kP1 / aa); - } - else { - return kBBMax; - } + // AliTRDpidESD copy constructor + // + + ((AliTRDpidESD &) p).Copy(*this); + +} + +//_____________________________________________________________________________ +AliTRDpidESD::~AliTRDpidESD() +{ + // + // Destructor + // + + if(fTrack) delete fTrack; } -//_________________________________________________________________________ -Int_t AliTRDpidESD::MakePID(AliESD *event) +//_____________________________________________________________________________ +AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p) { // - // This function calculates the "detector response" PID probabilities - // - Int_t ntrk=event->GetNumberOfTracks(); - for (Int_t i=0; iGetTrack(i); - if ((t->GetStatus()&AliESDtrack::kTRDin)==0) - if ((t->GetStatus()&AliESDtrack::kTRDout)==0) - if ((t->GetStatus()&AliESDtrack::kTRDrefit)==0) continue; - Int_t ns=AliPID::kSPECIES; - Double_t p[10]; - for (Int_t j=0; jGetP(); - Double_t dedx=t->GetTRDsignal()/fMIP; - Double_t bethe=Bethe(mom/mass); - Double_t sigma=fRes*bethe; - if (TMath::Abs(dedx-bethe) > fRange*sigma) { - p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; - continue; + // Assignment operator + // + + if (this != &p) ((AliTRDpidESD &) p).Copy(*this); + return *this; + +} + +//_____________________________________________________________________________ +void AliTRDpidESD::Copy(TObject &p) const +{ + // + // Copy function + // + + ((AliTRDpidESD &) p).fgCheckTrackStatus = fgCheckTrackStatus; + ((AliTRDpidESD &) p).fgCheckKinkStatus = fgCheckKinkStatus; + ((AliTRDpidESD &) p).fgMinPlane = fgMinPlane; + ((AliTRDpidESD &) p).fTrack = NULL; + +} + +//_____________________________________________________________________________ +Int_t AliTRDpidESD::MakePID(AliESDEvent * const event) +{ + // + // This function calculates the PID probabilities based on TRD signals + // + // The method produces probabilities based on the charge + // and the position of the maximum time bin in each layer. + // The dE/dx information can be used as global charge or 2 to 3 + // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual + // implementation. + // + // Author + // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 + // + + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliErrorGeneral("AliTRDpidESD::MakePID()" + ,"No access to calibration data"); + return -1; + } + +// AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam(); +// if (!rec) { +// AliErrorGeneral("AliTRDpidESD::MakePID()", "No TRD reco param."); +// return 0x0; +// } + + // Retrieve the CDB container class with the probability distributions + const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ); + if (!pd) { + AliErrorGeneral("AliTRDpidESD::MakePID()" + ,"No access to AliTRDCalPID"); + return -1; + } + + // Loop through all ESD tracks + Double_t p[10]; + AliESDtrack *t = NULL; + Float_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx; + Int_t timebin; + Float_t mom, length; + Int_t nPlanePID; + for (Int_t i=0; iGetNumberOfTracks(); i++) { + t = event->GetTrack(i); + + // Check track + if(!CheckTrack(t)) continue; + + + // Skip tracks which have no TRD signal at all + if (t->GetTRDsignal()<1.e-5) continue; + + // Loop over detector layers + mom = 0.; + length = 0.; + nPlanePID = 0; + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) + p[iSpecies] = 1./AliPID::kSPECIES; + + for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) { + // read data for track segment + for(int iSlice=0; iSliceGetTRDslice(iLayer, iSlice); + dEdx = t->GetTRDslice(iLayer, -1); + timebin = t->GetTRDTimBin(iLayer); + + // check data + if ((dEdx <= 0.) || (timebin <= -1.)) continue; + + // retrive kinematic info for this track segment + if(!RecalculateTrackSegmentKine(t, iLayer, mom, length)){ + // information is not fully reliable especialy for length + // estimation. To be used in the future. + } + + // this track segment has fulfilled all requierments + nPlanePID++; + + // Get the probabilities for the different particle species + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { + p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer); + //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin); } - p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; } + if(nPlanePID == 0) continue; + + // normalize probabilities + Double_t probTotal = 0.; + for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies]; + if(probTotal <= 0.){ + AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d." + , probTotal, i)); + AliWarning("This may be caused by some error in reference data."); + AliWarning("Calculation continues but results might be corrupted."); + continue; + } + for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal; + + // book PID to the track t->SetTRDpid(p); + t->SetTRDntracklets(nPlanePID<<3); } + return 0; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::CheckTrack(AliESDtrack * const t) +{ + // + // Check if track is eligible for PID calculations + // + + // Check the ESD track status + if (fgCheckTrackStatus) { + if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) && + ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE; + } + + // Check for ESD kink tracks + if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE; + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack * const esd + , Int_t layer + , Float_t &mom + , Float_t &length) +{ + // + // Retrive momentum "mom" and track "length" in TRD chamber from plane + // "plan" according to information stored in AliESDtrack "t". + // + // Origin + // Alex Bercuci (A.Bercuci@gsi.de) + // + + const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.; + const Float_t kDrWidth = AliTRDgeometry::DrThick(); + const Float_t kTime0 = AliTRDgeometry::GetTime0(layer); + + // set initial length value to chamber height + length = 2 * kAmHalfWidth + kDrWidth; + + // retrive track's outer param + const AliExternalTrackParam *op = esd->GetOuterParam(); + if(!op){ + mom = esd->GetP(); + return kFALSE; + } + + AliExternalTrackParam *param = NULL; + if(!fTrack){ + fTrack = new AliExternalTrackParam(*op); + param = fTrack; + } else param = new(fTrack) AliExternalTrackParam(*op); + + // retrive the magnetic field + Double_t xyz0[3]; + op->GetXYZ(xyz0); + Float_t field = AliTracker::GetBz(xyz0); // Bz in kG at point xyz0 + Double_t s, t; + + // propagate to chamber entrance + if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){ + mom = op->GetP(); + s = op->GetSnp(); + t = op->GetTgl(); + if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t)); + return kFALSE; + } + mom = param->GetP(); + s = param->GetSnp(); + t = param->GetTgl(); + if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t)); + + // check if track is crossing tracking sector by propagating to chamber exit- maybe is too much :) + Double_t alpha = param->GetAlpha(); + if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE; + + // mark track segments which are crossing SM boundaries along chamber + if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE; + + return kTRUE; + }