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