]>
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. // | |
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. // | |
27 | // // | |
dab811d0 | 28 | // Authors : // |
29 | // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)// | |
30 | // Alex Bercuci (a.bercuci@gsi.de) // | |
b52cb1b6 | 31 | // // |
32 | //////////////////////////////////////////////////////////////////////////// | |
b0f03c34 | 33 | |
b52cb1b6 | 34 | #include "AliLog.h" |
b0f03c34 | 35 | #include "AliESD.h" |
36 | #include "AliESDtrack.h" | |
c6011b06 | 37 | #include "AliTracker.h" |
b52cb1b6 | 38 | |
39 | #include "AliTRDpidESD.h" | |
40 | #include "AliTRDgeometry.h" | |
cc7cef99 | 41 | #include "AliTRDcalibDB.h" |
dab811d0 | 42 | #include "AliRun.h" |
43 | #include "AliTRDtrack.h" | |
7754cd1f | 44 | #include "Cal/AliTRDCalPIDLQ.h" |
b0f03c34 | 45 | |
c6011b06 | 46 | |
b0f03c34 | 47 | ClassImp(AliTRDpidESD) |
48 | ||
df34b21e | 49 | Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE; |
50 | Bool_t AliTRDpidESD::fCheckKinkStatus = kFALSE; | |
85815482 | 51 | Int_t AliTRDpidESD::fMinPlane = 0; |
df34b21e | 52 | |
b52cb1b6 | 53 | //_____________________________________________________________________________ |
dab811d0 | 54 | AliTRDpidESD::AliTRDpidESD() |
c6011b06 | 55 | :TObject(), fTrack(0x0) |
b0f03c34 | 56 | { |
57 | // | |
b52cb1b6 | 58 | // Default constructor |
b0f03c34 | 59 | // |
b52cb1b6 | 60 | |
b0f03c34 | 61 | } |
62 | ||
b52cb1b6 | 63 | //_____________________________________________________________________________ |
dab811d0 | 64 | AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) |
c6011b06 | 65 | :TObject(p), fTrack(0x0) |
eab5961e | 66 | { |
b0f03c34 | 67 | // |
b52cb1b6 | 68 | // AliTRDpidESD copy constructor |
b0f03c34 | 69 | // |
eab5961e | 70 | |
b52cb1b6 | 71 | ((AliTRDpidESD &) p).Copy(*this); |
72 | ||
73 | } | |
74 | ||
c6011b06 | 75 | //_____________________________________________________________________________ |
76 | AliTRDpidESD::~AliTRDpidESD() | |
77 | { | |
78 | // | |
79 | // Destructor | |
80 | // | |
81 | if(fTrack) delete fTrack; | |
82 | } | |
83 | ||
b52cb1b6 | 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 | ||
25ca55ce | 103 | ((AliTRDpidESD &) p).fCheckTrackStatus = fCheckTrackStatus; |
104 | ((AliTRDpidESD &) p).fCheckKinkStatus = fCheckKinkStatus; | |
105 | ((AliTRDpidESD &) p).fMinPlane = fMinPlane; | |
c6011b06 | 106 | ((AliTRDpidESD &) p).fTrack = 0x0; |
107 | ||
b0f03c34 | 108 | } |
109 | ||
b52cb1b6 | 110 | //_____________________________________________________________________________ |
b0f03c34 | 111 | Int_t AliTRDpidESD::MakePID(AliESD *event) |
bd50219c | 112 | { |
113 | // | |
b52cb1b6 | 114 | // This function calculates the PID probabilities based on TRD signals |
bd50219c | 115 | // |
dab811d0 | 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 AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual | |
120 | // implementation. | |
121 | // | |
122 | // Author | |
123 | // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 | |
124 | ||
125 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
126 | if (!calibration) { | |
127 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
128 | ,"No access to calibration data"); | |
129 | return -1; | |
130 | } | |
131 | ||
132 | // Retrieve the CDB container class with the probability distributions | |
133 | const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject(); | |
134 | if (!pd) { | |
135 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
136 | ,"No access to AliTRDCalPIDLQ"); | |
137 | return -1; | |
b52cb1b6 | 138 | } |
139 | ||
dab811d0 | 140 | // Loop through all ESD tracks |
141 | Double_t p[10]; | |
142 | AliESDtrack *t = 0x0; | |
c6011b06 | 143 | Float_t dedx[AliTRDtrack::kNslice], dEdx; |
144 | Int_t timebin; | |
145 | Float_t mom, length; | |
146 | Int_t nPlanePID; | |
dab811d0 | 147 | for (Int_t i=0; i<event->GetNumberOfTracks(); i++) { |
148 | t = event->GetTrack(i); | |
c6011b06 | 149 | |
dab811d0 | 150 | // Check track |
151 | if(!CheckTrack(t)) continue; | |
c6011b06 | 152 | |
dab811d0 | 153 | |
154 | // Skip tracks which have no TRD signal at all | |
155 | if (t->GetTRDsignal() == 0.) continue; | |
156 | ||
157 | // Loop over detector layers | |
c6011b06 | 158 | mom = 0.; |
159 | length = 0.; | |
dab811d0 | 160 | nPlanePID = 0; |
c6011b06 | 161 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1./AliPID::kSPECIES; |
dab811d0 | 162 | for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) { |
163 | // read data for track segment | |
164 | for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++) | |
165 | dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice); | |
166 | dEdx = t->GetTRDsignals(iPlan, -1); | |
167 | timebin = t->GetTRDTimBin(iPlan); | |
168 | ||
169 | // check data | |
170 | if ((dEdx <= 0.) || (timebin <= -1.)) continue; | |
171 | ||
172 | // retrive kinematic info for this track segment | |
c6011b06 | 173 | if(!RecalculateTrackSegmentKine(t, iPlan, mom, length)) continue; |
dab811d0 | 174 | |
175 | // this track segment has fulfilled all requierments | |
176 | nPlanePID++; | |
177 | ||
178 | // Get the probabilities for the different particle species | |
179 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { | |
180 | p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length); | |
c6011b06 | 181 | //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin); |
dab811d0 | 182 | } |
183 | } | |
c6011b06 | 184 | if(nPlanePID == 0) continue; |
185 | ||
dab811d0 | 186 | // normalize probabilities |
c6011b06 | 187 | Double_t probTotal = 0.; |
188 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies]; | |
189 | if(probTotal <= 0.){ | |
190 | AliWarningGeneral("AliTRDpidESD::MakePID()", | |
191 | Form("The total probability over all species <= 0 in ESD track %d. This may be caused by some error in reference data. Calculation continue but results might be corrupted.", i)); | |
192 | continue; | |
193 | } | |
194 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal; | |
dab811d0 | 195 | |
196 | // book PID to the track | |
197 | t->SetTRDpid(p); | |
c6011b06 | 198 | t->SetTRDpidQuality(nPlanePID); |
dab811d0 | 199 | } |
200 | ||
201 | return 0; | |
202 | } | |
b52cb1b6 | 203 | |
dab811d0 | 204 | //_____________________________________________________________________________ |
205 | Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t) | |
206 | { | |
207 | // | |
208 | // Check if track is eligible for PID calculations | |
209 | // | |
210 | ||
211 | // Check the ESD track status | |
212 | if (fCheckTrackStatus) { | |
213 | if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) && | |
214 | ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE; | |
215 | } | |
b52cb1b6 | 216 | |
dab811d0 | 217 | // Check for ESD kink tracks |
218 | if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE; | |
b52cb1b6 | 219 | |
dab811d0 | 220 | return kTRUE; |
221 | } | |
b52cb1b6 | 222 | |
dab811d0 | 223 | //_____________________________________________________________________________ |
c6011b06 | 224 | Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd, Int_t plan, Float_t &mom, Float_t &length) |
dab811d0 | 225 | { |
226 | // | |
227 | // Retrive momentum "mom" and track "length" in TRD chamber from plane | |
228 | // "plan" according to information stored in AliESDtrack "t". | |
c6011b06 | 229 | // |
230 | // Origin | |
231 | // Alex Bercuci (A.Bercuci@gsi.de) | |
dab811d0 | 232 | |
dab811d0 | 233 | const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.; |
f162af62 | 234 | const Float_t kDrWidth = AliTRDgeometry::DrThick(); |
c6011b06 | 235 | const Float_t kTime0 = AliTRDgeometry::GetTime0(plan); |
dab811d0 | 236 | |
c6011b06 | 237 | // set initial length value to chamber height |
238 | length = 2 * kAmHalfWidth + kDrWidth; | |
dab811d0 | 239 | |
c6011b06 | 240 | // retrive track's outer param |
241 | const AliExternalTrackParam *op = esd->GetOuterParam(); | |
242 | if(!op){ | |
243 | mom = esd->GetP(); | |
dab811d0 | 244 | return kFALSE; |
245 | } | |
dab811d0 | 246 | |
c6011b06 | 247 | AliExternalTrackParam *param = 0x0; |
248 | if(!fTrack){ | |
249 | fTrack = new AliExternalTrackParam(*op); | |
250 | param = fTrack; | |
251 | } else param = new(&fTrack) AliExternalTrackParam(*op); | |
252 | ||
253 | // retrive the magnetic field | |
254 | Double_t xyz0[3]; | |
255 | op->GetXYZ(xyz0); | |
256 | Float_t field = AliTracker::GetBz(xyz0); // Bz in kG at point xyz0 | |
257 | // propagate to chamber entrance | |
258 | if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){ | |
259 | mom = op->GetP(); | |
260 | return kFALSE; | |
261 | } | |
262 | mom = param->GetP(); | |
263 | Double_t s = param->GetSnp(); | |
264 | Double_t t = param->GetTgl(); | |
265 | length /= TMath::Sqrt((1. - s*s) / (1. - t*t)); | |
266 | ||
267 | // check if track is crossing tracking sector by propagating to chamber exit- maybe is too much :) | |
268 | Double_t alpha = param->GetAlpha(); | |
269 | if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE; | |
270 | ||
271 | // mark track segments which are crossing SM boundaries along chamber | |
272 | if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE; | |
273 | ||
dab811d0 | 274 | return kTRUE; |
b0f03c34 | 275 | } |
dab811d0 | 276 |