X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDpidESD.cxx;h=5f0739763aee3953fe304f8a21dac06b344c1059;hb=315b4daf2532e7b3902e8b9ee9bbc5f7d34cde58;hp=1b135254032b65d92c3e8c60863dabf03e814e70;hpb=eab5961ed9b69f5bb7036af12fd7b350e3a878f7;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDpidESD.cxx b/TRD/AliTRDpidESD.cxx index 1b135254032..5f0739763ae 100644 --- a/TRD/AliTRDpidESD.cxx +++ b/TRD/AliTRDpidESD.cxx @@ -13,102 +13,293 @@ * 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... -//----------------------------------------------------------------- +/* $Id$ */ -#include "AliTRDpidESD.h" -#include "AliESD.h" +//////////////////////////////////////////////////////////////////////////// +// // +// 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 "AliLog.h" +#include "AliESDEvent.h" #include "AliESDtrack.h" +#include "AliTracker.h" +#include "AliRun.h" + +#include "AliTRDReconstructor.h" +#include "AliTRDpidESD.h" +#include "AliTRDgeometry.h" +#include "AliTRDcalibDB.h" +#include "AliTRDtrack.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(0x0) +{ + // + // Default constructor + // + +} + +//_____________________________________________________________________________ +AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) + :TObject(p) + ,fTrack(0x0) { // - // The main constructor + // AliTRDpidESD copy constructor // - fMIP=param[0]; // MIP signal - fRes=param[1]; // relative resolution - fRange=param[2]; // PID "range" (in sigmas) + + ((AliTRDpidESD &) p).Copy(*this); + } -Double_t AliTRDpidESD::Bethe(Double_t bg) +//_____________________________________________________________________________ +AliTRDpidESD::~AliTRDpidESD() { // - // 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; - } + // 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 - // - static const Double_t masses[]={ - 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 - }; - 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=AliESDtrack::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; - } - p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; - } - t->SetTRDpid(p); - } - return 0; + // 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 = 0x0; + +} + +//_____________________________________________________________________________ +Int_t AliTRDpidESD::MakePID(AliESDEvent *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(AliTRDReconstructor::kLQPID/*rec->GetPIDMethod()*/); + if (!pd) { + AliErrorGeneral("AliTRDpidESD::MakePID()" + ,"No access to AliTRDCalPID"); + return -1; + } + + // Loop through all ESD tracks + Double_t p[10]; + AliESDtrack *t = 0x0; + Float_t dedx[AliTRDtrack::kNslice], 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() == 0.) 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); + } + } + 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->SetTRDpidQuality(nPlanePID); + } + + return 0; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *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 *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 = 0x0; + 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*s) / (1. + t*t)); + return kFALSE; + } + mom = param->GetP(); + s = param->GetSnp(); + t = param->GetTgl(); + if (s < 1.) length /= TMath::Sqrt((1. - s*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; + }