]>
Commit | Line | Data |
---|---|---|
e4f2f73d | 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 | // The TRD track seed // | |
21 | // // | |
22 | // Authors: // | |
23 | // Alex Bercuci <A.Bercuci@gsi.de> // | |
24 | // Markus Fasel <M.Fasel@gsi.de> // | |
25 | // // | |
26 | //////////////////////////////////////////////////////////////////////////// | |
27 | ||
28 | #include "TMath.h" | |
29 | #include "TLinearFitter.h" | |
30 | ||
31 | #include "AliLog.h" | |
32 | #include "AliMathBase.h" | |
33 | ||
34 | #include "AliTRDseedV1.h" | |
35 | #include "AliTRDcluster.h" | |
0906e73e | 36 | #include "AliTRDtrack.h" |
e4f2f73d | 37 | #include "AliTRDcalibDB.h" |
38 | #include "AliTRDstackLayer.h" | |
39 | #include "AliTRDrecoParam.h" | |
0906e73e | 40 | #include "AliTRDgeometry.h" |
41 | #include "Cal/AliTRDCalPID.h" | |
e4f2f73d | 42 | |
43 | #define SEED_DEBUG | |
44 | ||
45 | ClassImp(AliTRDseedV1) | |
46 | ||
47 | //____________________________________________________________________ | |
48 | AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) | |
49 | :AliTRDseed() | |
0906e73e | 50 | ,fPlane(layer) |
e4f2f73d | 51 | ,fOwner(kFALSE) |
0906e73e | 52 | ,fMom(0.) |
e4f2f73d | 53 | ,fRecoParam(p) |
54 | { | |
55 | // | |
56 | // Constructor | |
57 | // | |
0906e73e | 58 | for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.; |
59 | for(int itb=0; itb < knTimebins; itb++){ | |
60 | fdQdl[itb] = 0.; | |
61 | fdQ[itb] = 0.; | |
62 | } | |
63 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.; | |
e4f2f73d | 64 | } |
65 | ||
66 | //____________________________________________________________________ | |
0906e73e | 67 | AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref) |
e4f2f73d | 68 | :AliTRDseed((AliTRDseed&)ref) |
0906e73e | 69 | ,fPlane(ref.fPlane) |
e4f2f73d | 70 | ,fOwner(kFALSE) |
0906e73e | 71 | ,fMom(ref.fMom) |
e4f2f73d | 72 | ,fRecoParam(ref.fRecoParam) |
73 | { | |
74 | // | |
75 | // Copy Constructor performing a deep copy | |
76 | // | |
77 | ||
78 | //AliInfo(""); | |
0906e73e | 79 | if(ref.fOwner) SetOwner(); |
80 | for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice]; | |
81 | for(int itb=0; itb < knTimebins; itb++){ | |
82 | fdQdl[itb] = ref.fdQdl[itb]; | |
83 | fdQ[itb] = ref.fdQ[itb]; | |
fbb2ea06 | 84 | } |
0906e73e | 85 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec]; |
fbb2ea06 | 86 | } |
d9950a5a | 87 | |
0906e73e | 88 | |
e4f2f73d | 89 | //____________________________________________________________________ |
90 | AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref) | |
91 | { | |
92 | // | |
93 | // Assignment Operator using the copy function | |
94 | // | |
95 | ||
96 | //AliInfo(""); | |
97 | if(this != &ref){ | |
98 | ref.Copy(*this); | |
99 | } | |
100 | return *this; | |
101 | ||
102 | } | |
103 | ||
104 | //____________________________________________________________________ | |
105 | AliTRDseedV1::~AliTRDseedV1() | |
106 | { | |
107 | // | |
108 | // Destructor. The RecoParam object belongs to the underlying tracker. | |
109 | // | |
110 | ||
111 | //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO")); | |
112 | ||
0906e73e | 113 | if(fOwner) |
114 | for(int itb=0; itb<knTimebins; itb++){ | |
115 | if(!fClusters[itb]) continue; | |
116 | //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb)); | |
117 | delete fClusters[itb]; | |
118 | fClusters[itb] = 0x0; | |
119 | } | |
e4f2f73d | 120 | } |
121 | ||
122 | //____________________________________________________________________ | |
123 | void AliTRDseedV1::Copy(TObject &ref) const | |
124 | { | |
125 | // | |
126 | // Copy function | |
127 | // | |
128 | ||
129 | //AliInfo(""); | |
130 | AliTRDseedV1 &target = (AliTRDseedV1 &)ref; | |
131 | ||
0906e73e | 132 | target.fPlane = fPlane; |
133 | target.fMom = fMom; | |
134 | target.fRecoParam = fRecoParam; | |
135 | ||
136 | for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice]; | |
137 | for(int itb=0; itb < knTimebins; itb++){ | |
138 | target.fdQdl[itb] = fdQdl[itb]; | |
139 | target.fdQ[itb] = fdQ[itb]; | |
140 | } | |
141 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec]; | |
142 | ||
e4f2f73d | 143 | AliTRDseed::Copy(target); |
144 | } | |
145 | ||
0906e73e | 146 | |
147 | //____________________________________________________________ | |
148 | void AliTRDseedV1::Init(AliTRDtrack *track) | |
149 | { | |
150 | // Initialize this tracklet using the track information | |
151 | // | |
152 | // Parameters: | |
153 | // track - the TRD track used to initialize the tracklet | |
154 | // | |
155 | // Detailed description | |
156 | // The function sets the starting point and direction of the | |
157 | // tracklet according to the information from the TRD track. | |
158 | // | |
159 | // Caution | |
160 | // The TRD track has to be propagated to the beginning of the | |
161 | // chamber where the tracklet will be constructed | |
162 | // | |
163 | ||
164 | Double_t y, z; | |
165 | track->GetProlongation(fX0, y, z); | |
166 | fYref[0] = y; | |
33f721e9 | 167 | fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp()); |
0906e73e | 168 | fZref[0] = z; |
33f721e9 | 169 | fZref[1] = track->GetTgl(); |
0906e73e | 170 | |
171 | //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl()); | |
172 | } | |
173 | ||
174 | //____________________________________________________________________ | |
175 | Double_t* AliTRDseedV1::GetProbability() | |
176 | { | |
177 | // Fill probability array for tracklet from the DB. | |
178 | // | |
179 | // Parameters | |
180 | // | |
181 | // Output | |
182 | // returns pointer to the probability array and 0x0 if missing DB access | |
183 | // | |
184 | // Detailed description | |
185 | ||
186 | ||
187 | // retrive calibration db | |
188 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
189 | if (!calibration) { | |
190 | AliError("No access to calibration data"); | |
191 | return 0x0; | |
192 | } | |
193 | ||
194 | // Retrieve the CDB container class with the parametric detector response | |
195 | const AliTRDCalPID *pd = calibration->GetPIDObject(fRecoParam->GetPIDMethod()); | |
196 | if (!pd) { | |
197 | AliError("No access to AliTRDCalPID object"); | |
198 | return 0x0; | |
199 | } | |
200 | ||
201 | // calculate tracklet length TO DO | |
202 | Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()); | |
203 | /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane])); | |
204 | ||
205 | //calculate dE/dx | |
206 | CookdEdx(fRecoParam->GetNdEdxSlices()); | |
207 | ||
208 | // Sets the a priori probabilities | |
209 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) { | |
210 | fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane); | |
211 | } | |
212 | ||
213 | return &fProb[0]; | |
214 | } | |
215 | ||
216 | //____________________________________________________________________ | |
217 | void AliTRDseedV1::CookdEdx(Int_t nslices) | |
218 | { | |
219 | // Calculates average dE/dx for all slices and store them in the internal array fdEdx. | |
220 | // | |
221 | // Parameters: | |
222 | // nslices : number of slices for which dE/dx should be calculated | |
223 | // Output: | |
224 | // | |
225 | // Detailed description | |
226 | // Calculates average dE/dx for all slices. Depending on the PID methode | |
227 | // the number of slices can be 3 (LQ) or 8(NN). The calculation is based | |
228 | // on previously calculated quantities dQ/dl of each cluster. The | |
229 | // following effects are included in the calculation: | |
230 | // 1. calibration values for t0 and vdrift | |
231 | // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing()) | |
232 | // | |
233 | ||
234 | Int_t nclusters[knSlices]; | |
235 | for(int i=0; i<knSlices; i++){ | |
236 | fdEdx[i] = 0.; | |
237 | nclusters[i] = 0; | |
238 | } | |
239 | ||
240 | AliTRDcluster *cluster = 0x0; | |
241 | for(int ic=0; ic<fgTimeBins; ic++){ | |
242 | if(!(cluster = fClusters[ic])) continue; | |
243 | Int_t tb = cluster->GetLocalTimeBin(); | |
244 | ||
245 | // consider calibration effects | |
246 | if(tb < fTimeBin0 || tb >= fTimeBin0+fTimeBinsRange) continue; | |
247 | ||
248 | // consider cluster sharing ... TO DO | |
249 | //if(fRecoParam->GetClusterSharing() && cluster->GetSharing()) continue; | |
250 | ||
251 | Int_t slice = (tb-fTimeBin0)*nslices/fTimeBinsRange; | |
252 | fdEdx[slice] += fdQdl[ic]; | |
253 | nclusters[slice]++; | |
254 | } // End of loop over clusters | |
255 | ||
256 | // calculate mean charge per slice | |
257 | for(int is=0; is<nslices; is++) if(nclusters[is]) fdEdx[is] /= nclusters[is]; | |
258 | } | |
259 | ||
e4f2f73d | 260 | //____________________________________________________________________ |
261 | Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const | |
262 | { | |
263 | // | |
264 | // Returns a quality measurement of the current seed | |
265 | // | |
266 | ||
267 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
268 | return .5 * (18.0 - fN2) | |
269 | + 10.* TMath::Abs(fYfit[1] - fYref[1]) | |
270 | + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr) | |
271 | + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength; | |
272 | } | |
273 | ||
0906e73e | 274 | //____________________________________________________________________ |
275 | void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const | |
276 | { | |
277 | // Computes covariance in the y-z plane at radial point x | |
278 | ||
279 | const Float_t k0= .2; // to be checked in FindClusters | |
280 | Double_t sy20 = k0*TMath::Tan(fYfit[1]); sy20 *= sy20; | |
281 | ||
282 | Double_t sy2 = fSigmaY2*fSigmaY2 + sy20; | |
283 | Double_t sz2 = fPadLength/12.; | |
284 | ||
285 | //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2); | |
286 | ||
287 | cov[0] = sy2; | |
288 | cov[1] = fTilt*(sy2-sz2); | |
289 | cov[2] = sz2; | |
290 | } | |
291 | ||
292 | //____________________________________________________________________ | |
293 | void AliTRDseedV1::SetdQdl(Double_t length) | |
294 | { | |
295 | for(int ic=0; ic<fgTimeBins; ic++) fdQdl[ic] = fdQ[ic] *length; | |
296 | } | |
297 | ||
298 | //____________________________________________________________________ | |
299 | void AliTRDseedV1::SetOwner(Bool_t own) | |
300 | { | |
301 | //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO")); | |
302 | ||
303 | if(own){ | |
304 | for(int ic=0; ic<knTimebins; ic++){ | |
305 | if(!fClusters[ic]) continue; | |
306 | fClusters[ic] = new AliTRDcluster(*fClusters[ic]); | |
307 | } | |
308 | fOwner = kTRUE; | |
309 | } else { | |
310 | if(fOwner){ | |
311 | for(int ic=0; ic<knTimebins; ic++){ | |
312 | if(!fClusters[ic]) continue; | |
313 | delete fClusters[ic]; | |
314 | //fClusters[ic] = tracker->GetClusters(index) TODO | |
315 | } | |
316 | } | |
317 | fOwner = kFALSE; | |
318 | } | |
319 | } | |
320 | ||
e4f2f73d | 321 | //____________________________________________________________________ |
322 | Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer | |
323 | , Float_t quality | |
324 | , Bool_t kZcorr | |
325 | , AliTRDcluster *c) | |
326 | { | |
327 | // | |
328 | // Iterative process to register clusters to the seed. | |
329 | // In iteration 0 we try only one pad-row and if quality not | |
330 | // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows) | |
331 | // | |
332 | ||
333 | if(!fRecoParam){ | |
334 | AliError("Seed can not be used without a valid RecoParam."); | |
335 | return kFALSE; | |
336 | } | |
0906e73e | 337 | |
338 | //AliInfo(Form("TimeBins = %d TimeBinsRange = %d", fgTimeBins, fTimeBinsRange)); | |
339 | ||
e4f2f73d | 340 | Float_t tquality; |
341 | Double_t kroady = fRecoParam->GetRoad1y(); | |
342 | Double_t kroadz = fPadLength * .5 + 1.; | |
343 | ||
344 | // initialize configuration parameters | |
345 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
346 | Int_t niter = kZcorr ? 1 : 2; | |
347 | ||
348 | Double_t yexp, zexp; | |
349 | Int_t ncl = 0; | |
350 | // start seed update | |
351 | for (Int_t iter = 0; iter < niter; iter++) { | |
352 | //AliInfo(Form("iter = %i", iter)); | |
353 | ncl = 0; | |
0906e73e | 354 | for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) { |
e4f2f73d | 355 | // define searching configuration |
356 | Double_t dxlayer = layer[iTime].GetX() - fX0; | |
357 | if(c){ | |
358 | zexp = c->GetZ(); | |
359 | //Try 2 pad-rows in second iteration | |
360 | if (iter > 0) { | |
361 | zexp = fZref[0] + fZref[1] * dxlayer - zcorr; | |
362 | if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5; | |
363 | if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5; | |
364 | } | |
365 | } else zexp = fZref[0]; | |
366 | yexp = fYref[0] + fYref[1] * dxlayer - zcorr; | |
367 | // get cluster | |
368 | // printf("xexp = %3.3f ,yexp = %3.3f, zexp = %3.3f\n",layer[iTime].GetX(),yexp,zexp); | |
369 | // printf("layer[%i].GetNClusters() = %i\n", iTime, layer[iTime].GetNClusters()); | |
370 | Int_t index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz); | |
0906e73e | 371 | |
372 | // printf("%d[%d] x[%7.3f | %7.3f] y[%7.3f] z[%7.3f]\n", iTime, layer[iTime].GetNClusters(), dxlayer, layer[iTime].GetX(), yexp, zexp); | |
e4f2f73d | 373 | // for(Int_t iclk = 0; iclk < layer[iTime].GetNClusters(); iclk++){ |
374 | // AliTRDcluster *testcl = layer[iTime].GetCluster(iclk); | |
0906e73e | 375 | // printf("Cluster %i: %d x = %7.3f, y = %7.3f, z = %7.3f\n", iclk, testcl->GetLocalTimeBin(), testcl->GetX(), testcl->GetY(), testcl->GetZ()); |
e4f2f73d | 376 | // } |
377 | // printf("Index = %i\n",index); | |
0906e73e | 378 | |
e4f2f73d | 379 | if (index < 0) continue; |
380 | ||
381 | // Register cluster | |
382 | AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index); | |
383 | ||
384 | //printf("Cluster %i(0x%x): x = %3.3f, y = %3.3f, z = %3.3f\n", index, cl, cl->GetX(), cl->GetY(), cl->GetZ()); | |
385 | ||
d76231c8 | 386 | Int_t globalIndex = layer[iTime].GetGlobalIndex(index); |
387 | fIndexes[iTime] = globalIndex; | |
e4f2f73d | 388 | fClusters[iTime] = cl; |
389 | fX[iTime] = dxlayer; | |
390 | fY[iTime] = cl->GetY(); | |
391 | fZ[iTime] = cl->GetZ(); | |
0906e73e | 392 | fdQ[iTime] = cl->GetQ()/layer[iTime].GetdX(); |
393 | ||
e4f2f73d | 394 | // Debugging |
395 | ncl++; | |
396 | } | |
397 | ||
398 | #ifdef SEED_DEBUG | |
399 | // Int_t nclusters = 0; | |
400 | // Float_t fD[iter] = 0.; | |
0906e73e | 401 | // for(int ic=0; ic<fgTimeBins+1; ic++){ |
e4f2f73d | 402 | // AliTRDcluster *ci = fClusters[ic]; |
403 | // if(!ci) continue; | |
0906e73e | 404 | // for(int jc=ic+1; jc<fgTimeBins+1; jc++){ |
e4f2f73d | 405 | // AliTRDcluster *cj = fClusters[jc]; |
406 | // if(!cj) continue; | |
407 | // fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+ | |
408 | // (ci->GetZ()-cj->GetZ())*(ci->GetZ()-cj->GetZ())); | |
409 | // nclusters++; | |
410 | // } | |
411 | // } | |
412 | // if(nclusters) fD[iter] /= float(nclusters); | |
413 | #endif | |
414 | ||
415 | AliTRDseed::Update(); | |
416 | ||
417 | if(IsOK()){ | |
418 | tquality = GetQuality(kZcorr); | |
419 | if(tquality < quality) break; | |
420 | else quality = tquality; | |
421 | } | |
422 | kroadz *= 2.; | |
423 | } // Loop: iter | |
424 | if (!IsOK()) return kFALSE; | |
425 | ||
426 | CookLabels(); | |
427 | UpdateUsed(); | |
428 | return kTRUE; | |
429 | } | |
430 | ||
431 | //____________________________________________________________________ | |
0906e73e | 432 | Bool_t AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer |
433 | ,Bool_t kZcorr) | |
e4f2f73d | 434 | { |
435 | // | |
436 | // Projective algorithm to attach clusters to seeding tracklets | |
437 | // | |
438 | // Parameters | |
439 | // | |
440 | // Output | |
441 | // | |
442 | // Detailed description | |
443 | // 1. Collapse x coordinate for the full detector plane | |
444 | // 2. truncated mean on y (r-phi) direction | |
445 | // 3. purge clusters | |
446 | // 4. truncated mean on z direction | |
447 | // 5. purge clusters | |
448 | // 6. fit tracklet | |
449 | // | |
450 | ||
451 | if(!fRecoParam){ | |
452 | AliError("Seed can not be used without a valid RecoParam."); | |
453 | return kFALSE; | |
454 | } | |
455 | ||
0906e73e | 456 | const Int_t kClusterCandidates = 2 * knTimebins; |
e4f2f73d | 457 | |
458 | //define roads | |
459 | Double_t kroady = fRecoParam->GetRoad1y(); | |
460 | Double_t kroadz = fPadLength * 1.5 + 1.; | |
461 | // correction to y for the tilting angle | |
462 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
463 | ||
464 | // working variables | |
465 | AliTRDcluster *clusters[kClusterCandidates]; | |
0906e73e | 466 | Double_t cond[4], yexp[knTimebins], zexp[knTimebins], |
e4f2f73d | 467 | yres[kClusterCandidates], zres[kClusterCandidates]; |
0906e73e | 468 | Int_t ncl, *index = 0x0, tboundary[knTimebins]; |
e4f2f73d | 469 | |
470 | // Do cluster projection | |
471 | Int_t nYclusters = 0; Bool_t kEXIT = kFALSE; | |
0906e73e | 472 | for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) { |
e4f2f73d | 473 | fX[iTime] = layer[iTime].GetX() - fX0; |
474 | zexp[iTime] = fZref[0] + fZref[1] * fX[iTime]; | |
475 | yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr; | |
476 | ||
477 | // build condition and process clusters | |
478 | cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady; | |
479 | cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz; | |
480 | layer[iTime].GetClusters(cond, index, ncl); | |
481 | for(Int_t ic = 0; ic<ncl; ic++){ | |
0906e73e | 482 | AliTRDcluster *c = layer[iTime].GetCluster(index[ic]); |
e4f2f73d | 483 | clusters[nYclusters] = c; |
484 | yres[nYclusters++] = c->GetY() - yexp[iTime]; | |
485 | if(nYclusters >= kClusterCandidates) { | |
486 | AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates)); | |
487 | kEXIT = kTRUE; | |
488 | break; | |
489 | } | |
490 | } | |
491 | tboundary[iTime] = nYclusters; | |
492 | if(kEXIT) break; | |
493 | } | |
494 | ||
495 | // Evaluate truncated mean on the y direction | |
496 | Double_t mean, sigma; | |
497 | AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2); | |
498 | //purge cluster candidates | |
499 | Int_t nZclusters = 0; | |
500 | for(Int_t ic = 0; ic<nYclusters; ic++){ | |
501 | if(yres[ic] - mean > 4. * sigma){ | |
502 | clusters[ic] = 0x0; | |
503 | continue; | |
504 | } | |
505 | zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()]; | |
506 | } | |
507 | ||
508 | // Evaluate truncated mean on the z direction | |
509 | AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2); | |
510 | //purge cluster candidates | |
511 | for(Int_t ic = 0; ic<nZclusters; ic++){ | |
512 | if(zres[ic] - mean > 4. * sigma){ | |
513 | clusters[ic] = 0x0; | |
514 | continue; | |
515 | } | |
516 | } | |
517 | ||
518 | ||
519 | // Select only one cluster/TimeBin | |
520 | Int_t lastCluster = 0; | |
521 | fN2 = 0; | |
0906e73e | 522 | for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) { |
e4f2f73d | 523 | ncl = tboundary[iTime] - lastCluster; |
524 | if(!ncl) continue; | |
0906e73e | 525 | AliTRDcluster *c = 0x0; |
e4f2f73d | 526 | if(ncl == 1){ |
527 | c = clusters[lastCluster]; | |
528 | } else if(ncl > 1){ | |
529 | Float_t dold = 9999.; Int_t iptr = lastCluster; | |
530 | for(int ic=lastCluster; ic<tboundary[iTime]; ic++){ | |
531 | if(!clusters[ic]) continue; | |
532 | Float_t y = yexp[iTime] - clusters[ic]->GetY(); | |
533 | Float_t z = zexp[iTime] - clusters[ic]->GetZ(); | |
534 | Float_t d = y * y + z * z; | |
535 | if(d > dold) continue; | |
536 | dold = d; | |
537 | iptr = ic; | |
538 | } | |
539 | c = clusters[iptr]; | |
540 | } | |
0906e73e | 541 | //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index); |
542 | //fIndexes[iTime] = GlobalIndex; | |
e4f2f73d | 543 | fClusters[iTime] = c; |
544 | fY[iTime] = c->GetY(); | |
545 | fZ[iTime] = c->GetZ(); | |
0906e73e | 546 | fdQ[iTime] = c->GetQ()/layer[iTime].GetdX(); |
547 | lastCluster = tboundary[iTime]; | |
e4f2f73d | 548 | fN2++; |
549 | } | |
550 | ||
551 | // number of minimum numbers of clusters expected for the tracklet | |
0906e73e | 552 | Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins); |
e4f2f73d | 553 | if (fN2 < kClmin){ |
554 | AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin)); | |
555 | fN2 = 0; | |
556 | return kFALSE; | |
557 | } | |
0906e73e | 558 | |
559 | // update used clusters | |
560 | fNUsed = 0; | |
561 | for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) { | |
562 | if(!fClusters[iTime]) continue; | |
563 | if((fClusters[iTime]->IsUsed())) fNUsed++; | |
564 | } | |
565 | ||
566 | if (fN2-fNUsed < kClmin){ | |
567 | AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2)); | |
568 | fN2 = 0; | |
569 | return kFALSE; | |
570 | } | |
e4f2f73d | 571 | |
e4f2f73d | 572 | return kTRUE; |
573 | } | |
574 | ||
575 | //____________________________________________________________________ | |
0906e73e | 576 | Bool_t AliTRDseedV1::Fit() |
e4f2f73d | 577 | { |
578 | // | |
579 | // Linear fit of the tracklet | |
580 | // | |
581 | // Parameters : | |
582 | // | |
583 | // Output : | |
584 | // True if successful | |
585 | // | |
586 | // Detailed description | |
587 | // 2. Check if tracklet crosses pad row boundary | |
588 | // 1. Calculate residuals in the y (r-phi) direction | |
589 | // 3. Do a Least Square Fit to the data | |
590 | // | |
591 | ||
592 | //Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction | |
593 | Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing | |
594 | Float_t anglecor = fTilt * fZref[1]; // Correction to the angle | |
595 | ||
596 | // calculate residuals | |
0906e73e | 597 | Float_t yres[knTimebins]; // y (r-phi) residuals |
598 | Int_t zint[knTimebins], // Histograming of the z coordinate | |
599 | zout[2*knTimebins];// | |
e4f2f73d | 600 | |
601 | fN = 0; | |
0906e73e | 602 | for (Int_t iTime = 0; iTime < fTimeBinsRange; iTime++) { |
e4f2f73d | 603 | if (!fClusters[iTime]) continue; |
604 | yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime]; | |
0906e73e | 605 | zint[fN] = Int_t(fZ[iTime]); |
606 | fN++; | |
e4f2f73d | 607 | } |
608 | ||
609 | // calculate pad row boundary crosses | |
0906e73e | 610 | Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBinsRange); |
e4f2f73d | 611 | Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE); |
612 | fZProb = zout[0]; | |
613 | if(nz <= 1) zout[3] = 0; | |
614 | if(zout[1] + zout[3] < kClmin) { | |
615 | AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin)); | |
616 | return kFALSE; | |
617 | } | |
618 | // Z distance bigger than pad - length | |
619 | if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0; | |
620 | ||
621 | ||
622 | Double_t sumw = 0., | |
623 | sumwx = 0., | |
624 | sumwx2 = 0., | |
625 | sumwy = 0., | |
626 | sumwxy = 0., | |
627 | sumwz = 0., | |
628 | sumwxz = 0.; | |
629 | Int_t npads; | |
630 | fMPads = 0; | |
631 | fMeanz = 0.; | |
0906e73e | 632 | // we will use only the clusters which are in the detector range |
633 | for(int iTime=0; iTime<fTimeBinsRange; iTime++){ | |
e4f2f73d | 634 | fUsable[iTime] = kFALSE; |
635 | if (!fClusters[iTime]) continue; | |
636 | npads = fClusters[iTime]->GetNPads(); | |
637 | ||
638 | fUsable[iTime] = kTRUE; | |
639 | fN2++; | |
640 | fMPads += npads; | |
641 | Float_t weight = 1.0; | |
642 | if(npads > 5) weight = 0.2; | |
643 | else if(npads > 4) weight = 0.5; | |
644 | sumw += weight; | |
645 | sumwx += fX[iTime] * weight; | |
646 | sumwx2 += fX[iTime] * fX[iTime] * weight; | |
647 | sumwy += weight * yres[iTime]; | |
648 | sumwxy += weight * yres[iTime] * fX[iTime]; | |
649 | sumwz += weight * fZ[iTime]; | |
650 | sumwxz += weight * fZ[iTime] * fX[iTime]; | |
651 | } | |
652 | if (fN2 < kClmin){ | |
653 | AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin)); | |
654 | fN2 = 0; | |
655 | return kFALSE; | |
656 | } | |
657 | fMeanz = sumwz / sumw; | |
658 | fNChange = 0; | |
659 | ||
660 | // Tracklet on boundary | |
661 | Float_t correction = 0; | |
662 | if (fNChange > 0) { | |
663 | if (fMeanz < fZProb) correction = ycrosscor; | |
664 | if (fMeanz > fZProb) correction = -ycrosscor; | |
665 | } | |
666 | ||
667 | Double_t det = sumw * sumwx2 - sumwx * sumwx; | |
668 | fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det; | |
669 | fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det; | |
670 | ||
671 | fSigmaY2 = 0; | |
0906e73e | 672 | for (Int_t i = 0; i < fTimeBinsRange+1; i++) { |
e4f2f73d | 673 | if (!fUsable[i]) continue; |
674 | Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i]; | |
675 | fSigmaY2 += delta*delta; | |
676 | } | |
677 | fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2)); | |
678 | ||
679 | fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
680 | fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
681 | fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
682 | fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
683 | fYfitR[0] += fYref[0] + correction; | |
684 | fYfitR[1] += fYref[1]; | |
685 | fYfit[0] = fYfitR[0]; | |
686 | fYfit[1] = fYfitR[1]; | |
687 | ||
688 | return kTRUE; | |
689 | } | |
690 | ||
691 | //_____________________________________________________________________________ | |
692 | Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror) | |
693 | { | |
694 | // | |
695 | // Fit the Rieman tilt | |
696 | // | |
697 | ||
698 | // Fitting with tilting pads - kz not fixed | |
699 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
700 | Int_t nTimeBins = cal->GetNumberOfTimeBins(); | |
701 | TLinearFitter fitterT2(4,"hyp4"); | |
702 | fitterT2.StoreData(kTRUE); | |
703 | Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z | |
704 | ||
705 | Int_t npointsT = 0; | |
706 | fitterT2.ClearPoints(); | |
707 | ||
708 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
709 | // printf("\nLayer %d\n", iLayer); | |
710 | // cseed[iLayer].Print(); | |
711 | if (!cseed[iLayer].IsOK()) continue; | |
712 | Double_t tilt = cseed[iLayer].fTilt; | |
713 | ||
714 | for (Int_t itime = 0; itime < nTimeBins+1; itime++) { | |
715 | // printf("\ttime %d\n", itime); | |
716 | if (!cseed[iLayer].fUsable[itime]) continue; | |
717 | // x relative to the midle chamber | |
718 | Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2; | |
719 | Double_t y = cseed[iLayer].fY[itime]; | |
720 | Double_t z = cseed[iLayer].fZ[itime]; | |
721 | ||
722 | // | |
723 | // Tilted rieman | |
724 | // | |
725 | Double_t uvt[6]; | |
726 | Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x | |
727 | Double_t t = 1.0 / (x2*x2 + y*y); | |
728 | uvt[1] = t; | |
729 | uvt[0] = 2.0 * x2 * uvt[1]; | |
730 | uvt[2] = 2.0 * tilt * uvt[1]; | |
731 | uvt[3] = 2.0 * tilt *uvt[1] * x; | |
732 | uvt[4] = 2.0 * (y + tilt * z) * uvt[1]; | |
733 | ||
734 | Double_t error = 2.0 * uvt[1]; | |
735 | if (terror) { | |
736 | error *= cseed[iLayer].fSigmaY; | |
737 | } | |
738 | else { | |
739 | error *= 0.2; //Default error | |
740 | } | |
741 | // printf("\tadd point :\n"); | |
742 | // for(int i=0; i<5; i++) printf("%f ", uvt[i]); | |
743 | // printf("\n"); | |
744 | fitterT2.AddPoint(uvt,uvt[4],error); | |
745 | npointsT++; | |
746 | ||
747 | } | |
748 | ||
749 | } | |
750 | fitterT2.Eval(); | |
751 | Double_t rpolz0 = fitterT2.GetParameter(3); | |
752 | Double_t rpolz1 = fitterT2.GetParameter(4); | |
753 | ||
754 | // | |
755 | // Linear fitter - not possible to make boundaries | |
756 | // non accept non possible z and dzdx combination | |
757 | // | |
758 | Bool_t acceptablez = kTRUE; | |
759 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
760 | if (cseed[iLayer].IsOK()) { | |
761 | Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2); | |
762 | if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) { | |
763 | acceptablez = kFALSE; | |
764 | } | |
765 | } | |
766 | } | |
767 | if (!acceptablez) { | |
768 | Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0); | |
769 | Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5; | |
770 | fitterT2.FixParameter(3,zmf); | |
771 | fitterT2.FixParameter(4,dzmf); | |
772 | fitterT2.Eval(); | |
773 | fitterT2.ReleaseParameter(3); | |
774 | fitterT2.ReleaseParameter(4); | |
775 | rpolz0 = fitterT2.GetParameter(3); | |
776 | rpolz1 = fitterT2.GetParameter(4); | |
777 | } | |
778 | ||
779 | Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT); | |
780 | Double_t params[3]; | |
781 | params[0] = fitterT2.GetParameter(0); | |
782 | params[1] = fitterT2.GetParameter(1); | |
783 | params[2] = fitterT2.GetParameter(2); | |
784 | Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0]; | |
785 | ||
786 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
787 | ||
788 | Double_t x = cseed[iLayer].fX0; | |
789 | Double_t y = 0; | |
790 | Double_t dy = 0; | |
791 | Double_t z = 0; | |
792 | Double_t dz = 0; | |
793 | ||
794 | // y | |
795 | Double_t res2 = (x * params[0] + params[1]); | |
796 | res2 *= res2; | |
797 | res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2; | |
798 | if (res2 >= 0) { | |
799 | res2 = TMath::Sqrt(res2); | |
800 | y = (1.0 - res2) / params[0]; | |
801 | } | |
802 | ||
803 | //dy | |
804 | Double_t x0 = -params[1] / params[0]; | |
805 | if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) { | |
806 | Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); | |
807 | if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) { | |
808 | Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0)); | |
809 | if (params[0] < 0) res *= -1.0; | |
810 | dy = res; | |
811 | } | |
812 | } | |
813 | z = rpolz0 + rpolz1 * (x - xref2); | |
814 | dz = rpolz1; | |
815 | cseed[iLayer].fYref[0] = y; | |
816 | cseed[iLayer].fYref[1] = dy; | |
817 | cseed[iLayer].fZref[0] = z; | |
818 | cseed[iLayer].fZref[1] = dz; | |
819 | cseed[iLayer].fC = curvature; | |
820 | ||
821 | } | |
822 | ||
823 | return chi2TR; | |
824 | ||
825 | } | |
826 | ||
827 | //___________________________________________________________________ | |
828 | void AliTRDseedV1::Print() | |
829 | { | |
830 | // | |
831 | // Printing the seedstatus | |
832 | // | |
833 | ||
834 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
835 | Int_t nTimeBins = cal->GetNumberOfTimeBins(); | |
836 | ||
837 | printf("Seed status :\n"); | |
838 | printf(" fTilt = %f\n", fTilt); | |
839 | printf(" fPadLength = %f\n", fPadLength); | |
840 | printf(" fX0 = %f\n", fX0); | |
841 | for(int ic=0; ic<nTimeBins; ic++) { | |
842 | const Char_t *isUsable = fUsable[ic]?"Yes":"No"; | |
0906e73e | 843 | printf(" %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n" |
e4f2f73d | 844 | , ic |
845 | , fX[ic] | |
846 | , fY[ic] | |
847 | , fZ[ic] | |
848 | , fIndexes[ic] | |
0906e73e | 849 | , ((void*) fClusters[ic]) |
e4f2f73d | 850 | , isUsable); |
851 | } | |
852 | ||
853 | printf(" fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]); | |
854 | printf(" fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]); | |
855 | printf(" fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]); | |
856 | printf(" fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]); | |
857 | printf(" fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]); | |
858 | printf(" fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]); | |
859 | printf(" fSigmaY =%f\n", fSigmaY); | |
860 | printf(" fSigmaY2=%f\n", fSigmaY2); | |
861 | printf(" fMeanz =%f\n", fMeanz); | |
862 | printf(" fZProb =%f\n", fZProb); | |
863 | printf(" fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]); | |
864 | printf(" fN =%d\n", fN); | |
865 | printf(" fN2 =%d (>8 isOK)\n",fN2); | |
866 | printf(" fNUsed =%d\n", fNUsed); | |
867 | printf(" fFreq =%d\n", fFreq); | |
868 | printf(" fNChange=%d\n", fNChange); | |
869 | printf(" fMPads =%f\n", fMPads); | |
870 | ||
871 | printf(" fC =%f\n", fC); | |
872 | printf(" fCC =%f\n",fCC); | |
873 | printf(" fChi2 =%f\n", fChi2); | |
874 | printf(" fChi2Z =%f\n", fChi2Z); | |
875 | ||
876 | } |