]>
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" | |
b52cb1b6 | 37 | |
38 | #include "AliTRDpidESD.h" | |
39 | #include "AliTRDgeometry.h" | |
cc7cef99 | 40 | #include "AliTRDcalibDB.h" |
dab811d0 | 41 | #include "AliRun.h" |
42 | #include "AliTRDtrack.h" | |
7754cd1f | 43 | #include "Cal/AliTRDCalPIDLQ.h" |
b0f03c34 | 44 | |
45 | ClassImp(AliTRDpidESD) | |
46 | ||
df34b21e | 47 | Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE; |
48 | Bool_t AliTRDpidESD::fCheckKinkStatus = kFALSE; | |
85815482 | 49 | Int_t AliTRDpidESD::fMinPlane = 0; |
df34b21e | 50 | |
b52cb1b6 | 51 | //_____________________________________________________________________________ |
dab811d0 | 52 | AliTRDpidESD::AliTRDpidESD() |
53 | :TObject() | |
b0f03c34 | 54 | { |
55 | // | |
b52cb1b6 | 56 | // Default constructor |
b0f03c34 | 57 | // |
b52cb1b6 | 58 | |
b0f03c34 | 59 | } |
60 | ||
b52cb1b6 | 61 | //_____________________________________________________________________________ |
dab811d0 | 62 | AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p) |
63 | :TObject(p) | |
eab5961e | 64 | { |
b0f03c34 | 65 | // |
b52cb1b6 | 66 | // AliTRDpidESD copy constructor |
b0f03c34 | 67 | // |
eab5961e | 68 | |
b52cb1b6 | 69 | ((AliTRDpidESD &) p).Copy(*this); |
70 | ||
71 | } | |
72 | ||
73 | //_____________________________________________________________________________ | |
74 | AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p) | |
75 | { | |
76 | // | |
77 | // Assignment operator | |
78 | // | |
79 | ||
80 | if (this != &p) ((AliTRDpidESD &) p).Copy(*this); | |
81 | return *this; | |
82 | ||
83 | } | |
84 | ||
85 | //_____________________________________________________________________________ | |
86 | void AliTRDpidESD::Copy(TObject &p) const | |
87 | { | |
88 | // | |
89 | // Copy function | |
90 | // | |
91 | ||
25ca55ce | 92 | ((AliTRDpidESD &) p).fCheckTrackStatus = fCheckTrackStatus; |
93 | ((AliTRDpidESD &) p).fCheckKinkStatus = fCheckKinkStatus; | |
94 | ((AliTRDpidESD &) p).fMinPlane = fMinPlane; | |
eab5961e | 95 | |
b0f03c34 | 96 | } |
97 | ||
b52cb1b6 | 98 | //_____________________________________________________________________________ |
b0f03c34 | 99 | Int_t AliTRDpidESD::MakePID(AliESD *event) |
bd50219c | 100 | { |
101 | // | |
b52cb1b6 | 102 | // This function calculates the PID probabilities based on TRD signals |
bd50219c | 103 | // |
dab811d0 | 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 | |
108 | // implementation. | |
109 | // | |
110 | // Author | |
111 | // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 | |
112 | ||
113 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
114 | if (!calibration) { | |
115 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
116 | ,"No access to calibration data"); | |
117 | return -1; | |
118 | } | |
119 | ||
120 | // Retrieve the CDB container class with the probability distributions | |
121 | const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject(); | |
122 | if (!pd) { | |
123 | AliErrorGeneral("AliTRDpidESD::MakePID()" | |
124 | ,"No access to AliTRDCalPIDLQ"); | |
125 | return -1; | |
b52cb1b6 | 126 | } |
127 | ||
b52cb1b6 | 128 | |
dab811d0 | 129 | // Loop through all ESD tracks |
130 | Double_t p[10]; | |
131 | AliESDtrack *t = 0x0; | |
132 | Double_t dedx[AliTRDtrack::kNslice], dEdx; | |
133 | Int_t timebin; | |
134 | Float_t mom, length, probTotal; | |
135 | Int_t nPlanePID; | |
136 | for (Int_t i=0; i<event->GetNumberOfTracks(); i++) { | |
137 | t = event->GetTrack(i); | |
138 | ||
139 | // Check track | |
140 | if(!CheckTrack(t)) continue; | |
141 | ||
142 | // Skip tracks which have no TRD signal at all | |
143 | if (t->GetTRDsignal() == 0.) continue; | |
144 | ||
145 | // Loop over detector layers | |
146 | mom = 0.; //t->GetP(); | |
147 | length = 0.; | |
148 | probTotal = 0.; | |
149 | nPlanePID = 0; | |
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); | |
157 | ||
158 | // check data | |
159 | if ((dEdx <= 0.) || (timebin <= -1.)) continue; | |
160 | ||
161 | // retrive kinematic info for this track segment | |
162 | if(!GetTrackSegmentKine(t, iPlan, mom, length)) continue; | |
163 | ||
164 | // this track segment has fulfilled all requierments | |
165 | nPlanePID++; | |
166 | ||
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]; | |
172 | } | |
173 | } | |
174 | ||
175 | // normalize probabilities | |
176 | if(probTotal > 0.) | |
177 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) | |
178 | if(nPlanePID > fMinPlane) p[iSpecies] /= probTotal; | |
179 | else p[iSpecies] = 1.0; | |
180 | ||
181 | ||
182 | // book PID to the track | |
183 | t->SetTRDpid(p); | |
184 | } | |
185 | ||
186 | return 0; | |
187 | } | |
b52cb1b6 | 188 | |
dab811d0 | 189 | //_____________________________________________________________________________ |
190 | Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t) | |
191 | { | |
192 | // | |
193 | // Check if track is eligible for PID calculations | |
194 | // | |
195 | ||
196 | // Check the ESD track status | |
197 | if (fCheckTrackStatus) { | |
198 | if (((t->GetStatus() & AliESDtrack::kTRDout ) == 0) && | |
199 | ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE; | |
200 | } | |
b52cb1b6 | 201 | |
dab811d0 | 202 | // Check for ESD kink tracks |
203 | if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE; | |
b52cb1b6 | 204 | |
dab811d0 | 205 | return kTRUE; |
206 | } | |
b52cb1b6 | 207 | |
dab811d0 | 208 | //_____________________________________________________________________________ |
209 | Bool_t AliTRDpidESD::GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length) | |
210 | { | |
211 | // | |
212 | // Retrive momentum "mom" and track "length" in TRD chamber from plane | |
213 | // "plan" according to information stored in AliESDtrack "t". | |
214 | // | |
215 | ||
216 | if(!gAlice){ | |
217 | AliErrorGeneral("AliTRDpidESD::GetTrackSegmentKine()" | |
218 | ,"No gAlice object to retrive TRDgeometry and Magnetic fied - this has to be removed in the future."); | |
219 | return kFALSE; | |
220 | } | |
221 | ||
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()); | |
226 | if(!TRDgeom){ | |
227 | AliWarningGeneral("AliTRDpidESD::GetTrackSegmentKine()", "Cannot load TRD geometry from gAlice! Build a new one.\n"); | |
228 | TRDgeom = new AliTRDgeometry(); | |
229 | kSelfGeom = kTRUE; | |
230 | } | |
231 | const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.; | |
232 | const Float_t kDrWidth = AliTRDgeometry::DrThick(); | |
233 | ||
234 | ||
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 | |
240 | ||
241 | ||
242 | // find momentum at chamber entrance and track length in chamber | |
243 | AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam()); | |
244 | ||
245 | param->PropagateTo(TRDgeom->GetTime0(plan)+kAmHalfWidth, field); | |
246 | param->GetXYZ(xyz0); | |
247 | alpha = param->GetAlpha(); | |
248 | param->PropagateTo(TRDgeom->GetTime0(plan)-kAmHalfWidth-kDrWidth, field); | |
249 | // eliminate track segments which are crossing SM boundaries along chamber | |
aa612059 | 250 | if(TMath::Abs(alpha-param->GetAlpha())>.01){ |
dab811d0 | 251 | delete param; |
252 | if(kSelfGeom) delete TRDgeom; | |
253 | return kFALSE; | |
254 | } | |
255 | param->GetXYZ(xyz1); | |
256 | length = sqrt( | |
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]) | |
260 | ); | |
261 | param->GetPxPyPz(xyz1); | |
262 | mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]); | |
263 | delete param; | |
264 | if(kSelfGeom) delete TRDgeom; | |
265 | ||
266 | return kTRUE; | |
b0f03c34 | 267 | } |
dab811d0 | 268 |