Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
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 : **************************************************************************/
15 :
16 : /* $Id$ */
17 :
18 : ////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // Implementation of the TRD PID class //
21 : // //
22 : // Assigns the electron and pion likelihoods to each ESD track. //
23 : // The function MakePID(AliESDEvent *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. //
27 : // //
28 : // Authors : //
29 : // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (orig. version) //
30 : // Alex Bercuci (a.bercuci@gsi.de) //
31 : // //
32 : ////////////////////////////////////////////////////////////////////////////
33 :
34 : #include "AliESDEvent.h"
35 : #include "AliESDtrack.h"
36 : #include "AliTracker.h"
37 :
38 : #include "AliTRDpidESD.h"
39 : #include "AliTRDgeometry.h"
40 : #include "AliTRDcalibDB.h"
41 : #include "AliTRDCalPID.h"
42 :
43 16 : ClassImp(AliTRDpidESD)
44 :
45 : Bool_t AliTRDpidESD::fgCheckTrackStatus = kTRUE;
46 : Bool_t AliTRDpidESD::fgCheckKinkStatus = kFALSE;
47 : Int_t AliTRDpidESD::fgMinPlane = 0;
48 :
49 : //_____________________________________________________________________________
50 : AliTRDpidESD::AliTRDpidESD()
51 0 : :TObject()
52 0 : ,fTrack(NULL)
53 0 : {
54 : //
55 : // Default constructor
56 : //
57 :
58 0 : }
59 :
60 : //_____________________________________________________________________________
61 : AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
62 0 : :TObject(p)
63 0 : ,fTrack(NULL)
64 0 : {
65 : //
66 : // AliTRDpidESD copy constructor
67 : //
68 :
69 0 : ((AliTRDpidESD &) p).Copy(*this);
70 :
71 0 : }
72 :
73 : //_____________________________________________________________________________
74 : AliTRDpidESD::~AliTRDpidESD()
75 0 : {
76 : //
77 : // Destructor
78 : //
79 :
80 0 : if(fTrack) delete fTrack;
81 :
82 0 : }
83 :
84 : //_____________________________________________________________________________
85 : AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
86 : {
87 : //
88 : // Assignment operator
89 : //
90 :
91 0 : if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
92 0 : return *this;
93 :
94 : }
95 :
96 : //_____________________________________________________________________________
97 : void AliTRDpidESD::Copy(TObject &p) const
98 : {
99 : //
100 : // Copy function
101 : //
102 :
103 0 : ((AliTRDpidESD &) p).fgCheckTrackStatus = fgCheckTrackStatus;
104 0 : ((AliTRDpidESD &) p).fgCheckKinkStatus = fgCheckKinkStatus;
105 0 : ((AliTRDpidESD &) p).fgMinPlane = fgMinPlane;
106 0 : ((AliTRDpidESD &) p).fTrack = NULL;
107 :
108 0 : }
109 :
110 : //_____________________________________________________________________________
111 : Int_t AliTRDpidESD::MakePID(AliESDEvent * const event)
112 : {
113 : //
114 : // This function calculates the PID probabilities based on TRD signals
115 : //
116 : // The method produces probabilities based on the charge
117 : // and the position of the maximum time bin in each layer.
118 : // The dE/dx information can be used as global charge or 2 to 3
119 : // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
120 : // implementation.
121 : //
122 : // Author
123 : // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
124 : //
125 :
126 0 : AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
127 0 : if (!calibration) {
128 0 : AliErrorGeneral("AliTRDpidESD::MakePID()"
129 : ,"No access to calibration data");
130 0 : return -1;
131 : }
132 :
133 : // AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
134 : // if (!rec) {
135 : // AliErrorGeneral("AliTRDpidESD::MakePID()", "No TRD reco param.");
136 : // return 0x0;
137 : // }
138 :
139 : // Retrieve the CDB container class with the probability distributions
140 0 : const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ);
141 0 : if (!pd) {
142 0 : AliErrorGeneral("AliTRDpidESD::MakePID()"
143 : ,"No access to AliTRDCalPID");
144 0 : return -1;
145 : }
146 :
147 : // Loop through all ESD tracks
148 0 : Double_t p[10];
149 : AliESDtrack *t = NULL;
150 0 : Float_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx;
151 : Int_t timebin;
152 0 : Float_t mom, length;
153 : Int_t nPlanePID;
154 0 : for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
155 0 : t = event->GetTrack(i);
156 :
157 : // Check track
158 0 : if(!CheckTrack(t)) continue;
159 :
160 :
161 : // Skip tracks which have no TRD signal at all
162 0 : if (t->GetTRDsignal()<1.e-5) continue;
163 :
164 : // Loop over detector layers
165 0 : mom = 0.;
166 0 : length = 0.;
167 : nPlanePID = 0;
168 0 : for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
169 0 : p[iSpecies] = 1./AliPID::kSPECIES;
170 :
171 0 : for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) {
172 : // read data for track segment
173 0 : for(int iSlice=0; iSlice<AliTRDCalPID::kNSlicesLQ; iSlice++)
174 0 : dedx[iSlice] = t->GetTRDslice(iLayer, iSlice);
175 0 : dEdx = t->GetTRDslice(iLayer, -1);
176 0 : timebin = t->GetTRDTimBin(iLayer);
177 :
178 : // check data
179 0 : if ((dEdx <= 0.) || (timebin <= -1.)) continue;
180 :
181 : // retrive kinematic info for this track segment
182 0 : if(!RecalculateTrackSegmentKine(t, iLayer, mom, length)){
183 : // information is not fully reliable especialy for length
184 : // estimation. To be used in the future.
185 : }
186 :
187 : // this track segment has fulfilled all requierments
188 0 : nPlanePID++;
189 :
190 : // Get the probabilities for the different particle species
191 0 : for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
192 0 : p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer);
193 : //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
194 : }
195 0 : }
196 0 : if(nPlanePID == 0) continue;
197 :
198 : // normalize probabilities
199 : Double_t probTotal = 0.;
200 0 : for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
201 0 : if(probTotal <= 0.){
202 0 : AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d."
203 : , probTotal, i));
204 0 : AliWarning("This may be caused by some error in reference data.");
205 0 : AliWarning("Calculation continues but results might be corrupted.");
206 0 : continue;
207 : }
208 0 : for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
209 :
210 : // book PID to the track
211 0 : t->SetTRDpid(p);
212 0 : t->SetTRDntracklets(nPlanePID<<3);
213 0 : }
214 :
215 : return 0;
216 :
217 0 : }
218 :
219 : //_____________________________________________________________________________
220 : Bool_t AliTRDpidESD::CheckTrack(AliESDtrack * const t)
221 : {
222 : //
223 : // Check if track is eligible for PID calculations
224 : //
225 :
226 : // Check the ESD track status
227 0 : if (fgCheckTrackStatus) {
228 0 : if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) &&
229 0 : ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
230 : }
231 :
232 : // Check for ESD kink tracks
233 0 : if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
234 :
235 0 : return kTRUE;
236 :
237 0 : }
238 :
239 : //_____________________________________________________________________________
240 : Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack * const esd
241 : , Int_t layer
242 : , Float_t &mom
243 : , Float_t &length)
244 : {
245 : //
246 : // Retrive momentum "mom" and track "length" in TRD chamber from plane
247 : // "plan" according to information stored in AliESDtrack "t".
248 : //
249 : // Origin
250 : // Alex Bercuci (A.Bercuci@gsi.de)
251 : //
252 :
253 0 : const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
254 0 : const Float_t kDrWidth = AliTRDgeometry::DrThick();
255 0 : const Float_t kTime0 = AliTRDgeometry::GetTime0(layer);
256 :
257 : // set initial length value to chamber height
258 0 : length = 2 * kAmHalfWidth + kDrWidth;
259 :
260 : // retrive track's outer param
261 0 : const AliExternalTrackParam *op = esd->GetOuterParam();
262 0 : if(!op){
263 0 : mom = esd->GetP();
264 0 : return kFALSE;
265 : }
266 :
267 : AliExternalTrackParam *param = NULL;
268 0 : if(!fTrack){
269 0 : fTrack = new AliExternalTrackParam(*op);
270 : param = fTrack;
271 0 : } else param = new(fTrack) AliExternalTrackParam(*op);
272 :
273 : // retrive the magnetic field
274 0 : Double_t xyz0[3];
275 0 : op->GetXYZ(xyz0);
276 0 : Float_t field = AliTracker::GetBz(xyz0); // Bz in kG at point xyz0
277 : Double_t s, t;
278 :
279 : // propagate to chamber entrance
280 0 : if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){
281 0 : mom = op->GetP();
282 0 : s = op->GetSnp();
283 0 : t = op->GetTgl();
284 0 : if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
285 0 : return kFALSE;
286 : }
287 0 : mom = param->GetP();
288 0 : s = param->GetSnp();
289 0 : t = param->GetTgl();
290 0 : if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
291 :
292 : // check if track is crossing tracking sector by propagating to chamber exit- maybe is too much :)
293 0 : Double_t alpha = param->GetAlpha();
294 0 : if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE;
295 :
296 : // mark track segments which are crossing SM boundaries along chamber
297 0 : if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE;
298 :
299 0 : return kTRUE;
300 :
301 0 : }
|
