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