X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDpidESD.cxx;h=a7935c3ddb5e3fc415816ac47d94c0e1d4b3b884;hb=297c2e103544c15fe1d6721358c6a414a5b2b2cb;hp=eb4f888c967cedfe10b94286eaa53cfaad716fad;hpb=63a700c674cc00efd6fda0cb02d1b0eb717d7bc2;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDpidESD.cxx b/TRD/AliTRDpidESD.cxx index eb4f888c967..a7935c3ddb5 100644 --- a/TRD/AliTRDpidESD.cxx +++ b/TRD/AliTRDpidESD.cxx @@ -13,128 +13,286 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//----------------------------------------------------------------- -// Implementation of the TRD PID class -// Assigns the electron and pion liklihoods for each ESD track. -// The function MakePID(AliESD *event) calculates the probability -// of having dedx and the probability of having timbin at a given -// momentum (mom) and particle type k (0 for e) and (2 for pi) -// from the precalculated timbin distributions. -// Prashant Shukla -//----------------------------------------------------------------- +/* $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 "AliTRDpidESD.h" +#include "AliTRDgeometry.h" #include "AliTRDcalibDB.h" -#include "AliTRDCalPIDLQ.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) +{ + // + // AliTRDpidESD copy constructor + // + + ((AliTRDpidESD &) p).Copy(*this); + +} + +//_____________________________________________________________________________ +AliTRDpidESD::~AliTRDpidESD() +{ + // + // Destructor + // + + if(fTrack) delete fTrack; + +} + +//_____________________________________________________________________________ +AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p) +{ + // + // 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) { // - // The main constructor + // This function calculates the PID probabilities based on TRD signals // - fMIP=param[0]; // MIP signal - fRes=param[1]; // relative resolution - fRange=param[2]; // PID "range" (in sigmas) + // 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; + } -Double_t AliTRDpidESD::Bethe(Double_t bg) +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t) { // - // 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; - } + // 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; } -//_________________________________________________________________________ -Int_t AliTRDpidESD::MakePID(AliESD *event) +//_____________________________________________________________________________ +Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd + , Int_t plan + , Float_t &mom + , Float_t &length) { // - // This function calculates the "detector response" PID probabilities - // - - AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); - if (!calibration) - return -1; - - // The class AliTRDCalPIDLQ contains precalculated prob dis. - const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject(); - if (!pd) return -1; - - // Example to get mean for particle 2 (pi) and momentum number 4 (2 GeV) - // printf("%.2f \n", pd->GetMean(2, 4)); - // Example of use of Copy Constructor - // AliTRDCalPIDLQ *pd1 = new AliTRDCalPIDLQ(*pd); - - 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; - if(t->GetTRDsignal()==0) continue; - // Int_t ns=AliESDtrack::kSPECIES; - Int_t ns=AliPID::kSPECIES; - Double_t p[10]; - Double_t mom=t->GetP(); - Double_t probTotal=0.0; - for (Int_t j=0; jGetTRDsignals(ilayer); - Int_t timbin=t->GetTRDTimBin(ilayer); - p[j]*= pd->GetProbability(j,mom,dedx); - p[j]*= pd->GetProbabilityT(j,mom,timbin); - p[j]*= 100; - } // loop over layers - probTotal+=p[j]; - } //loop over particle species - // printf(" %f %d %f %f %f \n", mom, timbin, p[0], p[1], p[2]); - for (Int_t j=0; jSetTRDpid(p); - } //loop over tracks - return 0; + // 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; + } + + 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; + }