X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDpidESD.cxx;h=a7935c3ddb5e3fc415816ac47d94c0e1d4b3b884;hb=fd3ef1361d8592c882ba4e000f6588d1acb17e76;hp=6b973081392d9be26400ab3910745fc97a8b0a28;hpb=b52cb1b6952f5882b016aa6abe43bd002afe7485;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDpidESD.cxx b/TRD/AliTRDpidESD.cxx index 6b973081392..a7935c3ddb5 100644 --- a/TRD/AliTRDpidESD.cxx +++ b/TRD/AliTRDpidESD.cxx @@ -13,47 +13,57 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + //////////////////////////////////////////////////////////////////////////// // // // Implementation of the TRD PID class // // // // Assigns the electron and pion likelihoods to each ESD track. // -// The function MakePID(AliESD *event) calculates the probability // +// 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. // // // -// Original version: // -// Prashant Shukla // +// Authors : // +// Prashant Shukla (orig. version) // +// Alex Bercuci (a.bercuci@gsi.de) // // // //////////////////////////////////////////////////////////////////////////// #include "AliLog.h" -#include "AliESD.h" +#include "AliESDEvent.h" #include "AliESDtrack.h" +#include "AliTracker.h" +#include "AliRun.h" #include "AliTRDpidESD.h" #include "AliTRDgeometry.h" #include "AliTRDcalibDB.h" -#include "Cal/AliTRDCalPIDLQ.h" +#include "AliTRDtrack.h" +#include "Cal/AliTRDCalPID.h" ClassImp(AliTRDpidESD) +Bool_t AliTRDpidESD::fgCheckTrackStatus = kTRUE; +Bool_t AliTRDpidESD::fgCheckKinkStatus = kFALSE; +Int_t AliTRDpidESD::fgMinPlane = 0; + //_____________________________________________________________________________ -AliTRDpidESD::AliTRDpidESD():TObject() +AliTRDpidESD::AliTRDpidESD() + :TObject() + ,fTrack(0x0) { // // Default constructor // - fCheckTrackStatus = kTRUE; - fCheckKinkStatus = kFALSE; - fMinPlane = 0; - } //_____________________________________________________________________________ -AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p) +AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) + :TObject(p) + ,fTrack(0x0) { // // AliTRDpidESD copy constructor @@ -63,6 +73,17 @@ AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p) } +//_____________________________________________________________________________ +AliTRDpidESD::~AliTRDpidESD() +{ + // + // Destructor + // + + if(fTrack) delete fTrack; + +} + //_____________________________________________________________________________ AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p) { @@ -82,116 +103,196 @@ void AliTRDpidESD::Copy(TObject &p) const // Copy function // - ((AliTRDpidESD &) p).fCheckTrackStatus = fCheckTrackStatus; - ((AliTRDpidESD &) p).fCheckKinkStatus = fCheckKinkStatus; - ((AliTRDpidESD &) p).fMinPlane = fMinPlane; - + ((AliTRDpidESD &) p).fgCheckTrackStatus = fgCheckTrackStatus; + ((AliTRDpidESD &) p).fgCheckKinkStatus = fgCheckKinkStatus; + ((AliTRDpidESD &) p).fgMinPlane = fgMinPlane; + ((AliTRDpidESD &) p).fTrack = 0x0; + } //_____________________________________________________________________________ -Int_t AliTRDpidESD::MakePID(AliESD *event) +Int_t AliTRDpidESD::MakePID(AliESDEvent *event) { // // This function calculates the PID probabilities based on TRD signals // - // So far this method produces probabilities based on the total charge - // in each layer and the position of the maximum time bin in each layer. - // In a final version this should also exploit the charge measurement in - // the different slices of a given layer. - // - - Double_t p[10]; - Int_t nSpecies = AliPID::kSPECIES; - Int_t nPlanePID = 0; - Double_t mom = 0.0; - Double_t probTotal = 0.0; - - AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); - if (!calibration) { - //AliError("No access to calibration data\n"); - return -1; - } - - // Retrieve the CDB container class with the probability distributions - const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject(); - if (!pd) { - //AliError("No access to AliTRDCalPIDLQ\n"); - return -1; - } - - // Loop through all ESD tracks - Int_t ntrk = event->GetNumberOfTracks(); - for (Int_t i = 0; i < ntrk; i++) { - - AliESDtrack *t = event->GetTrack(i); - - // Check the ESD track status - if (fCheckTrackStatus) { - if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) && - ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) { - continue; - } - } - - // Check for ESD kink tracks - if (fCheckKinkStatus) { - if (t->GetKinkIndex(0) != 0) { - continue; - } - } - - // Skip tracks that have no TRD signal at all - if (t->GetTRDsignal() == 0) { - continue; - } - - mom = t->GetP(); - probTotal = 0.0; - nPlanePID = 0; - for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - p[iSpecies] = 1.0; - } - - // Check the different detector layers - for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) { - - // Use the total charge in a given plane - Double_t dedx = t->GetTRDsignals(iPlan,-1); - Int_t timebin = t->GetTRDTimBin(iPlan); - if ((dedx > 0.0) && - (timebin > -1.0)) { - - nPlanePID++; - - // Get the probabilities for the different particle species - for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - - p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx); - p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin); - p[iSpecies] *= 100.0; // ?????????????? - - probTotal += p[iSpecies]; + // 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; + } + + // Retrieve the CDB container class with the probability distributions + const AliTRDCalPID *pd = calibration->GetPIDObject(1); + 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 iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) { + // read data for track segment + for(int iSlice=0; iSliceGetTRDsignals(iPlan, iSlice); + dEdx = t->GetTRDsignals(iPlan, -1); + timebin = t->GetTRDTimBin(iPlan); + + // check data + if ((dEdx <= 0.) || (timebin <= -1.)) continue; + + // retrive kinematic info for this track segment + if(!RecalculateTrackSegmentKine(t, iPlan, 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, iPlan); + //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; - } +} - for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - if ((probTotal > 0.0) && - (nPlanePID > fMinPlane)) { - p[iSpecies] /= probTotal; - } - else { - p[iSpecies] = 1.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; + +} - t->SetTRDpid(p); +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd + , Int_t plan + , 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(plan); + + // 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; + } - return 0; + 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; }