1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////
20 // Implementation of the TRD PID class //
22 // Assigns the electron and pion likelihoods to each ESD track. //
23 // The function MakePID(AliESD *event) calculates the probability //
24 // of having dedx and a maximum timbin at a given //
25 // momentum (mom) and particle type k //
26 // from the precalculated distributions. //
29 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)//
30 // Alex Bercuci (a.bercuci@gsi.de) //
32 ////////////////////////////////////////////////////////////////////////////
36 #include "AliESDtrack.h"
38 #include "AliTRDpidESD.h"
39 #include "AliTRDgeometry.h"
40 #include "AliTRDcalibDB.h"
42 #include "AliTRDtrack.h"
43 #include "Cal/AliTRDCalPIDLQ.h"
45 ClassImp(AliTRDpidESD)
47 Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
48 Bool_t AliTRDpidESD::fCheckKinkStatus = kFALSE;
49 Int_t AliTRDpidESD::fMinPlane = 0;
51 //_____________________________________________________________________________
52 AliTRDpidESD::AliTRDpidESD()
56 // Default constructor
61 //_____________________________________________________________________________
62 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
66 // AliTRDpidESD copy constructor
69 ((AliTRDpidESD &) p).Copy(*this);
73 //_____________________________________________________________________________
74 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
77 // Assignment operator
80 if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
85 //_____________________________________________________________________________
86 void AliTRDpidESD::Copy(TObject &p) const
92 ((AliTRDpidESD &) p).fCheckTrackStatus = fCheckTrackStatus;
93 ((AliTRDpidESD &) p).fCheckKinkStatus = fCheckKinkStatus;
94 ((AliTRDpidESD &) p).fMinPlane = fMinPlane;
98 //_____________________________________________________________________________
99 Int_t AliTRDpidESD::MakePID(AliESD *event)
102 // This function calculates the PID probabilities based on TRD signals
104 // The method produces probabilities based on the charge
105 // and the position of the maximum time bin in each layer.
106 // The dE/dx information can be used as global charge or 2 to 3
107 // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
111 // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
113 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
115 AliErrorGeneral("AliTRDpidESD::MakePID()"
116 ,"No access to calibration data");
120 // Retrieve the CDB container class with the probability distributions
121 const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
123 AliErrorGeneral("AliTRDpidESD::MakePID()"
124 ,"No access to AliTRDCalPIDLQ");
129 // Loop through all ESD tracks
131 AliESDtrack *t = 0x0;
132 Double_t dedx[AliTRDtrack::kNslice], dEdx;
134 Float_t mom, length, probTotal;
136 for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
137 t = event->GetTrack(i);
140 if(!CheckTrack(t)) continue;
142 // Skip tracks which have no TRD signal at all
143 if (t->GetTRDsignal() == 0.) continue;
145 // Loop over detector layers
146 mom = 0.; //t->GetP();
150 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1.;
151 for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
152 // read data for track segment
153 for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
154 dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice);
155 dEdx = t->GetTRDsignals(iPlan, -1);
156 timebin = t->GetTRDTimBin(iPlan);
159 if ((dEdx <= 0.) || (timebin <= -1.)) continue;
161 // retrive kinematic info for this track segment
162 if(!GetTrackSegmentKine(t, iPlan, mom, length)) continue;
164 // this track segment has fulfilled all requierments
167 // Get the probabilities for the different particle species
168 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
169 p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
170 p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
171 probTotal += p[iSpecies];
175 // normalize probabilities
177 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
178 if(nPlanePID > fMinPlane) p[iSpecies] /= probTotal;
179 else p[iSpecies] = 1.0;
182 // book PID to the track
189 //_____________________________________________________________________________
190 Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t)
193 // Check if track is eligible for PID calculations
196 // Check the ESD track status
197 if (fCheckTrackStatus) {
198 if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) &&
199 ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
202 // Check for ESD kink tracks
203 if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
208 //_____________________________________________________________________________
209 Bool_t AliTRDpidESD::GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length)
212 // Retrive momentum "mom" and track "length" in TRD chamber from plane
213 // "plan" according to information stored in AliESDtrack "t".
217 AliErrorGeneral("AliTRDpidESD::GetTrackSegmentKine()"
218 ,"No gAlice object to retrive TRDgeometry and Magnetic fied - this has to be removed in the future.");
222 // Retrieve TRD geometry -> Maybe there is a better way to do this
223 Bool_t kSelfGeom = kFALSE;
224 AliTRDgeometry *TRDgeom =0x0;
225 if(gAlice) TRDgeom = AliTRDgeometry::GetGeometry(gAlice->GetRunLoader());
227 AliWarningGeneral("AliTRDpidESD::GetTrackSegmentKine()", "Cannot load TRD geometry from gAlice! Build a new one.\n");
228 TRDgeom = new AliTRDgeometry();
231 const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
232 const Float_t kDrWidth = AliTRDgeometry::DrThick();
235 // retrive the magnetic field
236 Double_t xyz0[3] = { 0., 0., 0.}, xyz1[3];
237 Double_t b[3], alpha;
238 gAlice->Field(xyz0,b); // b[] is in kilo Gauss
239 Float_t field = b[2] * 0.1; // Tesla
242 // find momentum at chamber entrance and track length in chamber
243 AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam());
245 param->PropagateTo(TRDgeom->GetTime0(plan)+kAmHalfWidth, field);
247 alpha = param->GetAlpha();
248 param->PropagateTo(TRDgeom->GetTime0(plan)-kAmHalfWidth-kDrWidth, field);
249 // eliminate track segments which are crossing SM boundaries along chamber
250 if(TMath::Abs(alpha-param->GetAlpha())>.01){
252 if(kSelfGeom) delete TRDgeom;
257 (xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+
258 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])+
259 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2])
261 param->GetPxPyPz(xyz1);
262 mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]);
264 if(kSelfGeom) delete TRDgeom;