]>
Commit | Line | Data |
---|---|---|
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 "Cal/AliTRDCalPID.h" | |
42 | ||
43 | 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 | :TObject() | |
52 | ,fTrack(NULL) | |
53 | { | |
54 | // | |
55 | // Default constructor | |
56 | // | |
57 | ||
58 | } | |
59 | ||
60 | //_____________________________________________________________________________ | |
61 | AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) | |
62 | :TObject(p) | |
63 | ,fTrack(NULL) | |
64 | { | |
65 | // | |
66 | // AliTRDpidESD copy constructor | |
67 | // | |
68 | ||
69 | ((AliTRDpidESD &) p).Copy(*this); | |
70 | ||
71 | } | |
72 | ||
73 | //_____________________________________________________________________________ | |
74 | AliTRDpidESD::~AliTRDpidESD() | |
75 | { | |
76 | // | |
77 | // Destructor | |
78 | // | |
79 | ||
80 | if(fTrack) delete fTrack; | |
81 | ||
82 | } | |
83 | ||
84 | //_____________________________________________________________________________ | |
85 | AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p) | |
86 | { | |
87 | // | |
88 | // Assignment operator | |
89 | // | |
90 | ||
91 | if (this != &p) ((AliTRDpidESD &) p).Copy(*this); | |
92 | return *this; | |
93 | ||
94 | } | |
95 | ||
96 | //_____________________________________________________________________________ | |
97 | void AliTRDpidESD::Copy(TObject &p) const | |
98 | { | |
99 | // | |
100 | // Copy function | |
101 | // | |
102 | ||
103 | ((AliTRDpidESD &) p).fgCheckTrackStatus = fgCheckTrackStatus; | |
104 | ((AliTRDpidESD &) p).fgCheckKinkStatus = fgCheckKinkStatus; | |
105 | ((AliTRDpidESD &) p).fgMinPlane = fgMinPlane; | |
106 | ((AliTRDpidESD &) p).fTrack = NULL; | |
107 | ||
108 | } | |
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 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
127 | if (!calibration) { | |
128 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
129 | ,"No access to calibration data"); | |
130 | 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 | const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ); | |
141 | if (!pd) { | |
142 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
143 | ,"No access to AliTRDCalPID"); | |
144 | return -1; | |
145 | } | |
146 | ||
147 | // Loop through all ESD tracks | |
148 | Double_t p[10]; | |
149 | AliESDtrack *t = NULL; | |
150 | Float_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx; | |
151 | Int_t timebin; | |
152 | Float_t mom, length; | |
153 | Int_t nPlanePID; | |
154 | for (Int_t i=0; i<event->GetNumberOfTracks(); i++) { | |
155 | t = event->GetTrack(i); | |
156 | ||
157 | // Check track | |
158 | if(!CheckTrack(t)) continue; | |
159 | ||
160 | ||
161 | // Skip tracks which have no TRD signal at all | |
162 | if (t->GetTRDsignal()<1.e-5) continue; | |
163 | ||
164 | // Loop over detector layers | |
165 | mom = 0.; | |
166 | length = 0.; | |
167 | nPlanePID = 0; | |
168 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) | |
169 | p[iSpecies] = 1./AliPID::kSPECIES; | |
170 | ||
171 | for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) { | |
172 | // read data for track segment | |
173 | for(int iSlice=0; iSlice<AliTRDCalPID::kNSlicesLQ; iSlice++) | |
174 | dedx[iSlice] = t->GetTRDslice(iLayer, iSlice); | |
175 | dEdx = t->GetTRDslice(iLayer, -1); | |
176 | timebin = t->GetTRDTimBin(iLayer); | |
177 | ||
178 | // check data | |
179 | if ((dEdx <= 0.) || (timebin <= -1.)) continue; | |
180 | ||
181 | // retrive kinematic info for this track segment | |
182 | 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 | nPlanePID++; | |
189 | ||
190 | // Get the probabilities for the different particle species | |
191 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { | |
192 | p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer); | |
193 | //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin); | |
194 | } | |
195 | } | |
196 | if(nPlanePID == 0) continue; | |
197 | ||
198 | // normalize probabilities | |
199 | Double_t probTotal = 0.; | |
200 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies]; | |
201 | if(probTotal <= 0.){ | |
202 | AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d." | |
203 | , probTotal, i)); | |
204 | AliWarning("This may be caused by some error in reference data."); | |
205 | AliWarning("Calculation continues but results might be corrupted."); | |
206 | continue; | |
207 | } | |
208 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal; | |
209 | ||
210 | // book PID to the track | |
211 | t->SetTRDpid(p); | |
212 | t->SetTRDntracklets(nPlanePID<<3); | |
213 | } | |
214 | ||
215 | return 0; | |
216 | ||
217 | } | |
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 | if (fgCheckTrackStatus) { | |
228 | if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) && | |
229 | ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE; | |
230 | } | |
231 | ||
232 | // Check for ESD kink tracks | |
233 | if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE; | |
234 | ||
235 | return kTRUE; | |
236 | ||
237 | } | |
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 | const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.; | |
254 | const Float_t kDrWidth = AliTRDgeometry::DrThick(); | |
255 | const Float_t kTime0 = AliTRDgeometry::GetTime0(layer); | |
256 | ||
257 | // set initial length value to chamber height | |
258 | length = 2 * kAmHalfWidth + kDrWidth; | |
259 | ||
260 | // retrive track's outer param | |
261 | const AliExternalTrackParam *op = esd->GetOuterParam(); | |
262 | if(!op){ | |
263 | mom = esd->GetP(); | |
264 | return kFALSE; | |
265 | } | |
266 | ||
267 | AliExternalTrackParam *param = NULL; | |
268 | if(!fTrack){ | |
269 | fTrack = new AliExternalTrackParam(*op); | |
270 | param = fTrack; | |
271 | } else param = new(fTrack) AliExternalTrackParam(*op); | |
272 | ||
273 | // retrive the magnetic field | |
274 | Double_t xyz0[3]; | |
275 | op->GetXYZ(xyz0); | |
276 | Float_t field = AliTracker::GetBz(xyz0); // Bz in kG at point xyz0 | |
277 | Double_t s, t; | |
278 | ||
279 | // propagate to chamber entrance | |
280 | if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){ | |
281 | mom = op->GetP(); | |
282 | s = op->GetSnp(); | |
283 | t = op->GetTgl(); | |
284 | if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t)); | |
285 | return kFALSE; | |
286 | } | |
287 | mom = param->GetP(); | |
288 | s = param->GetSnp(); | |
289 | t = param->GetTgl(); | |
290 | 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 | Double_t alpha = param->GetAlpha(); | |
294 | if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE; | |
295 | ||
296 | // mark track segments which are crossing SM boundaries along chamber | |
297 | if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE; | |
298 | ||
299 | return kTRUE; | |
300 | ||
301 | } |