]>
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" | |
36 | #include "AliTRDcalibDB.h" | |
37 | #include "AliTRDstackLayer.h" | |
38 | #include "AliTRDrecoParam.h" | |
39 | ||
40 | #define SEED_DEBUG | |
41 | ||
42 | ClassImp(AliTRDseedV1) | |
43 | ||
44 | //____________________________________________________________________ | |
45 | AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) | |
46 | :AliTRDseed() | |
47 | ,fLayer(layer) | |
48 | ,fTimeBins(0) | |
49 | ,fOwner(kFALSE) | |
50 | ,fRecoParam(p) | |
51 | { | |
52 | // | |
53 | // Constructor | |
54 | // | |
55 | ||
56 | //AliInfo(""); | |
57 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
58 | fTimeBins = cal->GetNumberOfTimeBins(); | |
59 | ||
60 | } | |
61 | ||
62 | //____________________________________________________________________ | |
63 | AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner) | |
64 | :AliTRDseed((AliTRDseed&)ref) | |
65 | ,fLayer(ref.fLayer) | |
66 | ,fTimeBins(ref.fTimeBins) | |
67 | ,fOwner(kFALSE) | |
68 | ,fRecoParam(ref.fRecoParam) | |
69 | { | |
70 | // | |
71 | // Copy Constructor performing a deep copy | |
72 | // | |
73 | ||
74 | //AliInfo(""); | |
75 | ||
76 | if(owner){ | |
77 | for(int ic=0; ic<fTimeBins; ic++){ | |
78 | if(!fClusters[ic]) continue; | |
79 | fClusters[ic] = new AliTRDcluster(*fClusters[ic]); | |
80 | } | |
81 | fOwner = kTRUE; | |
82 | } | |
83 | ||
84 | } | |
85 | ||
86 | //____________________________________________________________________ | |
87 | AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref) | |
88 | { | |
89 | // | |
90 | // Assignment Operator using the copy function | |
91 | // | |
92 | ||
93 | //AliInfo(""); | |
94 | if(this != &ref){ | |
95 | ref.Copy(*this); | |
96 | } | |
97 | return *this; | |
98 | ||
99 | } | |
100 | ||
101 | //____________________________________________________________________ | |
102 | AliTRDseedV1::~AliTRDseedV1() | |
103 | { | |
104 | // | |
105 | // Destructor. The RecoParam object belongs to the underlying tracker. | |
106 | // | |
107 | ||
108 | //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO")); | |
109 | ||
110 | if(fOwner) delete [] fClusters; | |
111 | } | |
112 | ||
113 | //____________________________________________________________________ | |
114 | void AliTRDseedV1::Copy(TObject &ref) const | |
115 | { | |
116 | // | |
117 | // Copy function | |
118 | // | |
119 | ||
120 | //AliInfo(""); | |
121 | AliTRDseedV1 &target = (AliTRDseedV1 &)ref; | |
122 | ||
123 | target.fLayer = fLayer; | |
124 | target.fTimeBins = fTimeBins; | |
125 | target.fRecoParam = fRecoParam; | |
126 | AliTRDseed::Copy(target); | |
127 | } | |
128 | ||
129 | //____________________________________________________________________ | |
130 | Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const | |
131 | { | |
132 | // | |
133 | // Returns a quality measurement of the current seed | |
134 | // | |
135 | ||
136 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
137 | return .5 * (18.0 - fN2) | |
138 | + 10.* TMath::Abs(fYfit[1] - fYref[1]) | |
139 | + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr) | |
140 | + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength; | |
141 | } | |
142 | ||
143 | //____________________________________________________________________ | |
144 | Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer | |
145 | , Float_t quality | |
146 | , Bool_t kZcorr | |
147 | , AliTRDcluster *c) | |
148 | { | |
149 | // | |
150 | // Iterative process to register clusters to the seed. | |
151 | // In iteration 0 we try only one pad-row and if quality not | |
152 | // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows) | |
153 | // | |
154 | ||
155 | if(!fRecoParam){ | |
156 | AliError("Seed can not be used without a valid RecoParam."); | |
157 | return kFALSE; | |
158 | } | |
159 | ||
160 | Float_t tquality; | |
161 | Double_t kroady = fRecoParam->GetRoad1y(); | |
162 | Double_t kroadz = fPadLength * .5 + 1.; | |
163 | ||
164 | // initialize configuration parameters | |
165 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
166 | Int_t niter = kZcorr ? 1 : 2; | |
167 | ||
168 | Double_t yexp, zexp; | |
169 | Int_t ncl = 0; | |
170 | // start seed update | |
171 | for (Int_t iter = 0; iter < niter; iter++) { | |
172 | //AliInfo(Form("iter = %i", iter)); | |
173 | ncl = 0; | |
174 | for (Int_t iTime = 0; iTime < fTimeBins; iTime++) { | |
175 | // define searching configuration | |
176 | Double_t dxlayer = layer[iTime].GetX() - fX0; | |
177 | if(c){ | |
178 | zexp = c->GetZ(); | |
179 | //Try 2 pad-rows in second iteration | |
180 | if (iter > 0) { | |
181 | zexp = fZref[0] + fZref[1] * dxlayer - zcorr; | |
182 | if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5; | |
183 | if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5; | |
184 | } | |
185 | } else zexp = fZref[0]; | |
186 | yexp = fYref[0] + fYref[1] * dxlayer - zcorr; | |
187 | // get cluster | |
188 | // printf("xexp = %3.3f ,yexp = %3.3f, zexp = %3.3f\n",layer[iTime].GetX(),yexp,zexp); | |
189 | // printf("layer[%i].GetNClusters() = %i\n", iTime, layer[iTime].GetNClusters()); | |
190 | Int_t index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz); | |
191 | // for(Int_t iclk = 0; iclk < layer[iTime].GetNClusters(); iclk++){ | |
192 | // AliTRDcluster *testcl = layer[iTime].GetCluster(iclk); | |
193 | // printf("Cluster %i: x = %3.3f, y = %3.3f, z = %3.3f\n",iclk,testcl->GetX(), testcl->GetY(), testcl->GetZ()); | |
194 | // } | |
195 | // printf("Index = %i\n",index); | |
196 | if (index < 0) continue; | |
197 | ||
198 | // Register cluster | |
199 | AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index); | |
200 | ||
201 | //printf("Cluster %i(0x%x): x = %3.3f, y = %3.3f, z = %3.3f\n", index, cl, cl->GetX(), cl->GetY(), cl->GetZ()); | |
202 | ||
203 | Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index); | |
204 | fIndexes[iTime] = GlobalIndex; | |
205 | fClusters[iTime] = cl; | |
206 | fX[iTime] = dxlayer; | |
207 | fY[iTime] = cl->GetY(); | |
208 | fZ[iTime] = cl->GetZ(); | |
209 | ||
210 | // Debugging | |
211 | ncl++; | |
212 | } | |
213 | ||
214 | #ifdef SEED_DEBUG | |
215 | // Int_t nclusters = 0; | |
216 | // Float_t fD[iter] = 0.; | |
217 | // for(int ic=0; ic<fTimeBins+1; ic++){ | |
218 | // AliTRDcluster *ci = fClusters[ic]; | |
219 | // if(!ci) continue; | |
220 | // for(int jc=ic+1; jc<fTimeBins+1; jc++){ | |
221 | // AliTRDcluster *cj = fClusters[jc]; | |
222 | // if(!cj) continue; | |
223 | // fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+ | |
224 | // (ci->GetZ()-cj->GetZ())*(ci->GetZ()-cj->GetZ())); | |
225 | // nclusters++; | |
226 | // } | |
227 | // } | |
228 | // if(nclusters) fD[iter] /= float(nclusters); | |
229 | #endif | |
230 | ||
231 | AliTRDseed::Update(); | |
232 | ||
233 | if(IsOK()){ | |
234 | tquality = GetQuality(kZcorr); | |
235 | if(tquality < quality) break; | |
236 | else quality = tquality; | |
237 | } | |
238 | kroadz *= 2.; | |
239 | } // Loop: iter | |
240 | if (!IsOK()) return kFALSE; | |
241 | ||
242 | CookLabels(); | |
243 | UpdateUsed(); | |
244 | return kTRUE; | |
245 | } | |
246 | ||
247 | //____________________________________________________________________ | |
248 | Bool_t AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer | |
249 | , Float_t /*quality*/ | |
250 | , Bool_t kZcorr | |
251 | , AliTRDcluster *c) | |
252 | { | |
253 | // | |
254 | // Projective algorithm to attach clusters to seeding tracklets | |
255 | // | |
256 | // Parameters | |
257 | // | |
258 | // Output | |
259 | // | |
260 | // Detailed description | |
261 | // 1. Collapse x coordinate for the full detector plane | |
262 | // 2. truncated mean on y (r-phi) direction | |
263 | // 3. purge clusters | |
264 | // 4. truncated mean on z direction | |
265 | // 5. purge clusters | |
266 | // 6. fit tracklet | |
267 | // | |
268 | ||
269 | if(!fRecoParam){ | |
270 | AliError("Seed can not be used without a valid RecoParam."); | |
271 | return kFALSE; | |
272 | } | |
273 | ||
274 | const Int_t knTimeBins = 35; | |
275 | const Int_t kClusterCandidates = 2 * knTimeBins; | |
276 | ||
277 | //define roads | |
278 | Double_t kroady = fRecoParam->GetRoad1y(); | |
279 | Double_t kroadz = fPadLength * 1.5 + 1.; | |
280 | // correction to y for the tilting angle | |
281 | Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.; | |
282 | ||
283 | // working variables | |
284 | AliTRDcluster *clusters[kClusterCandidates]; | |
285 | Double_t cond[4], yexp[knTimeBins], zexp[knTimeBins], | |
286 | yres[kClusterCandidates], zres[kClusterCandidates]; | |
287 | Int_t ncl, *index = 0x0, tboundary[knTimeBins]; | |
288 | ||
289 | // Do cluster projection | |
290 | Int_t nYclusters = 0; Bool_t kEXIT = kFALSE; | |
291 | for (Int_t iTime = 0; iTime < fTimeBins; iTime++) { | |
292 | fX[iTime] = layer[iTime].GetX() - fX0; | |
293 | zexp[iTime] = fZref[0] + fZref[1] * fX[iTime]; | |
294 | yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr; | |
295 | ||
296 | // build condition and process clusters | |
297 | cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady; | |
298 | cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz; | |
299 | layer[iTime].GetClusters(cond, index, ncl); | |
300 | for(Int_t ic = 0; ic<ncl; ic++){ | |
301 | c = layer[iTime].GetCluster(index[ic]); | |
302 | clusters[nYclusters] = c; | |
303 | yres[nYclusters++] = c->GetY() - yexp[iTime]; | |
304 | if(nYclusters >= kClusterCandidates) { | |
305 | AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates)); | |
306 | kEXIT = kTRUE; | |
307 | break; | |
308 | } | |
309 | } | |
310 | tboundary[iTime] = nYclusters; | |
311 | if(kEXIT) break; | |
312 | } | |
313 | ||
314 | // Evaluate truncated mean on the y direction | |
315 | Double_t mean, sigma; | |
316 | AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2); | |
317 | //purge cluster candidates | |
318 | Int_t nZclusters = 0; | |
319 | for(Int_t ic = 0; ic<nYclusters; ic++){ | |
320 | if(yres[ic] - mean > 4. * sigma){ | |
321 | clusters[ic] = 0x0; | |
322 | continue; | |
323 | } | |
324 | zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()]; | |
325 | } | |
326 | ||
327 | // Evaluate truncated mean on the z direction | |
328 | AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2); | |
329 | //purge cluster candidates | |
330 | for(Int_t ic = 0; ic<nZclusters; ic++){ | |
331 | if(zres[ic] - mean > 4. * sigma){ | |
332 | clusters[ic] = 0x0; | |
333 | continue; | |
334 | } | |
335 | } | |
336 | ||
337 | ||
338 | // Select only one cluster/TimeBin | |
339 | Int_t lastCluster = 0; | |
340 | fN2 = 0; | |
341 | for (Int_t iTime = 0; iTime < fTimeBins; iTime++) { | |
342 | ncl = tboundary[iTime] - lastCluster; | |
343 | if(!ncl) continue; | |
344 | if(ncl == 1){ | |
345 | c = clusters[lastCluster]; | |
346 | } else if(ncl > 1){ | |
347 | Float_t dold = 9999.; Int_t iptr = lastCluster; | |
348 | for(int ic=lastCluster; ic<tboundary[iTime]; ic++){ | |
349 | if(!clusters[ic]) continue; | |
350 | Float_t y = yexp[iTime] - clusters[ic]->GetY(); | |
351 | Float_t z = zexp[iTime] - clusters[ic]->GetZ(); | |
352 | Float_t d = y * y + z * z; | |
353 | if(d > dold) continue; | |
354 | dold = d; | |
355 | iptr = ic; | |
356 | } | |
357 | c = clusters[iptr]; | |
358 | } | |
359 | //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index); | |
360 | //fIndexes[iTime] = GlobalIndex; | |
361 | fClusters[iTime] = c; | |
362 | fY[iTime] = c->GetY(); | |
363 | fZ[iTime] = c->GetZ(); | |
364 | lastCluster = tboundary[iTime]; | |
365 | fN2++; | |
366 | } | |
367 | ||
368 | // number of minimum numbers of clusters expected for the tracklet | |
369 | Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins); | |
370 | if (fN2 < kClmin){ | |
371 | AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin)); | |
372 | fN2 = 0; | |
373 | return kFALSE; | |
374 | } | |
375 | AliTRDseed::Update(); | |
376 | ||
377 | // // fit tracklet and update clusters | |
378 | // if(!FitTracklet()) return kFALSE; | |
379 | // UpdateUsed(); | |
380 | return kTRUE; | |
381 | } | |
382 | ||
383 | //____________________________________________________________________ | |
384 | Bool_t AliTRDseedV1::FitTracklet() | |
385 | { | |
386 | // | |
387 | // Linear fit of the tracklet | |
388 | // | |
389 | // Parameters : | |
390 | // | |
391 | // Output : | |
392 | // True if successful | |
393 | // | |
394 | // Detailed description | |
395 | // 2. Check if tracklet crosses pad row boundary | |
396 | // 1. Calculate residuals in the y (r-phi) direction | |
397 | // 3. Do a Least Square Fit to the data | |
398 | // | |
399 | ||
400 | //Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction | |
401 | Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing | |
402 | Float_t anglecor = fTilt * fZref[1]; // Correction to the angle | |
403 | ||
404 | // calculate residuals | |
405 | const Int_t knTimeBins = 35; | |
406 | Float_t yres[knTimeBins]; // y (r-phi) residuals | |
407 | Int_t zint[knTimeBins], // Histograming of the z coordinate | |
408 | zout[2*knTimeBins];// | |
409 | ||
410 | fN = 0; | |
411 | for (Int_t iTime = 0; iTime < fTimeBins; iTime++) { | |
412 | if (!fClusters[iTime]) continue; | |
413 | yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime]; | |
414 | zint[fN++] = Int_t(fZ[iTime]); | |
415 | } | |
416 | ||
417 | // calculate pad row boundary crosses | |
418 | Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins); | |
419 | Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE); | |
420 | fZProb = zout[0]; | |
421 | if(nz <= 1) zout[3] = 0; | |
422 | if(zout[1] + zout[3] < kClmin) { | |
423 | AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin)); | |
424 | return kFALSE; | |
425 | } | |
426 | // Z distance bigger than pad - length | |
427 | if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0; | |
428 | ||
429 | ||
430 | Double_t sumw = 0., | |
431 | sumwx = 0., | |
432 | sumwx2 = 0., | |
433 | sumwy = 0., | |
434 | sumwxy = 0., | |
435 | sumwz = 0., | |
436 | sumwxz = 0.; | |
437 | Int_t npads; | |
438 | fMPads = 0; | |
439 | fMeanz = 0.; | |
440 | for(int iTime=0; iTime<fTimeBins; iTime++){ | |
441 | fUsable[iTime] = kFALSE; | |
442 | if (!fClusters[iTime]) continue; | |
443 | npads = fClusters[iTime]->GetNPads(); | |
444 | ||
445 | fUsable[iTime] = kTRUE; | |
446 | fN2++; | |
447 | fMPads += npads; | |
448 | Float_t weight = 1.0; | |
449 | if(npads > 5) weight = 0.2; | |
450 | else if(npads > 4) weight = 0.5; | |
451 | sumw += weight; | |
452 | sumwx += fX[iTime] * weight; | |
453 | sumwx2 += fX[iTime] * fX[iTime] * weight; | |
454 | sumwy += weight * yres[iTime]; | |
455 | sumwxy += weight * yres[iTime] * fX[iTime]; | |
456 | sumwz += weight * fZ[iTime]; | |
457 | sumwxz += weight * fZ[iTime] * fX[iTime]; | |
458 | } | |
459 | if (fN2 < kClmin){ | |
460 | AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin)); | |
461 | fN2 = 0; | |
462 | return kFALSE; | |
463 | } | |
464 | fMeanz = sumwz / sumw; | |
465 | fNChange = 0; | |
466 | ||
467 | // Tracklet on boundary | |
468 | Float_t correction = 0; | |
469 | if (fNChange > 0) { | |
470 | if (fMeanz < fZProb) correction = ycrosscor; | |
471 | if (fMeanz > fZProb) correction = -ycrosscor; | |
472 | } | |
473 | ||
474 | Double_t det = sumw * sumwx2 - sumwx * sumwx; | |
475 | fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det; | |
476 | fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det; | |
477 | ||
478 | fSigmaY2 = 0; | |
479 | for (Int_t i = 0; i < fTimeBins+1; i++) { | |
480 | if (!fUsable[i]) continue; | |
481 | Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i]; | |
482 | fSigmaY2 += delta*delta; | |
483 | } | |
484 | fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2)); | |
485 | ||
486 | fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
487 | fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
488 | fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
489 | fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
490 | fYfitR[0] += fYref[0] + correction; | |
491 | fYfitR[1] += fYref[1]; | |
492 | fYfit[0] = fYfitR[0]; | |
493 | fYfit[1] = fYfitR[1]; | |
494 | ||
495 | return kTRUE; | |
496 | } | |
497 | ||
498 | //_____________________________________________________________________________ | |
499 | Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror) | |
500 | { | |
501 | // | |
502 | // Fit the Rieman tilt | |
503 | // | |
504 | ||
505 | // Fitting with tilting pads - kz not fixed | |
506 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
507 | Int_t nTimeBins = cal->GetNumberOfTimeBins(); | |
508 | TLinearFitter fitterT2(4,"hyp4"); | |
509 | fitterT2.StoreData(kTRUE); | |
510 | Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z | |
511 | ||
512 | Int_t npointsT = 0; | |
513 | fitterT2.ClearPoints(); | |
514 | ||
515 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
516 | // printf("\nLayer %d\n", iLayer); | |
517 | // cseed[iLayer].Print(); | |
518 | if (!cseed[iLayer].IsOK()) continue; | |
519 | Double_t tilt = cseed[iLayer].fTilt; | |
520 | ||
521 | for (Int_t itime = 0; itime < nTimeBins+1; itime++) { | |
522 | // printf("\ttime %d\n", itime); | |
523 | if (!cseed[iLayer].fUsable[itime]) continue; | |
524 | // x relative to the midle chamber | |
525 | Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2; | |
526 | Double_t y = cseed[iLayer].fY[itime]; | |
527 | Double_t z = cseed[iLayer].fZ[itime]; | |
528 | ||
529 | // | |
530 | // Tilted rieman | |
531 | // | |
532 | Double_t uvt[6]; | |
533 | Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x | |
534 | Double_t t = 1.0 / (x2*x2 + y*y); | |
535 | uvt[1] = t; | |
536 | uvt[0] = 2.0 * x2 * uvt[1]; | |
537 | uvt[2] = 2.0 * tilt * uvt[1]; | |
538 | uvt[3] = 2.0 * tilt *uvt[1] * x; | |
539 | uvt[4] = 2.0 * (y + tilt * z) * uvt[1]; | |
540 | ||
541 | Double_t error = 2.0 * uvt[1]; | |
542 | if (terror) { | |
543 | error *= cseed[iLayer].fSigmaY; | |
544 | } | |
545 | else { | |
546 | error *= 0.2; //Default error | |
547 | } | |
548 | // printf("\tadd point :\n"); | |
549 | // for(int i=0; i<5; i++) printf("%f ", uvt[i]); | |
550 | // printf("\n"); | |
551 | fitterT2.AddPoint(uvt,uvt[4],error); | |
552 | npointsT++; | |
553 | ||
554 | } | |
555 | ||
556 | } | |
557 | fitterT2.Eval(); | |
558 | Double_t rpolz0 = fitterT2.GetParameter(3); | |
559 | Double_t rpolz1 = fitterT2.GetParameter(4); | |
560 | ||
561 | // | |
562 | // Linear fitter - not possible to make boundaries | |
563 | // non accept non possible z and dzdx combination | |
564 | // | |
565 | Bool_t acceptablez = kTRUE; | |
566 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
567 | if (cseed[iLayer].IsOK()) { | |
568 | Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2); | |
569 | if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) { | |
570 | acceptablez = kFALSE; | |
571 | } | |
572 | } | |
573 | } | |
574 | if (!acceptablez) { | |
575 | Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0); | |
576 | Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5; | |
577 | fitterT2.FixParameter(3,zmf); | |
578 | fitterT2.FixParameter(4,dzmf); | |
579 | fitterT2.Eval(); | |
580 | fitterT2.ReleaseParameter(3); | |
581 | fitterT2.ReleaseParameter(4); | |
582 | rpolz0 = fitterT2.GetParameter(3); | |
583 | rpolz1 = fitterT2.GetParameter(4); | |
584 | } | |
585 | ||
586 | Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT); | |
587 | Double_t params[3]; | |
588 | params[0] = fitterT2.GetParameter(0); | |
589 | params[1] = fitterT2.GetParameter(1); | |
590 | params[2] = fitterT2.GetParameter(2); | |
591 | Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0]; | |
592 | ||
593 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
594 | ||
595 | Double_t x = cseed[iLayer].fX0; | |
596 | Double_t y = 0; | |
597 | Double_t dy = 0; | |
598 | Double_t z = 0; | |
599 | Double_t dz = 0; | |
600 | ||
601 | // y | |
602 | Double_t res2 = (x * params[0] + params[1]); | |
603 | res2 *= res2; | |
604 | res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2; | |
605 | if (res2 >= 0) { | |
606 | res2 = TMath::Sqrt(res2); | |
607 | y = (1.0 - res2) / params[0]; | |
608 | } | |
609 | ||
610 | //dy | |
611 | Double_t x0 = -params[1] / params[0]; | |
612 | if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) { | |
613 | Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); | |
614 | if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) { | |
615 | Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0)); | |
616 | if (params[0] < 0) res *= -1.0; | |
617 | dy = res; | |
618 | } | |
619 | } | |
620 | z = rpolz0 + rpolz1 * (x - xref2); | |
621 | dz = rpolz1; | |
622 | cseed[iLayer].fYref[0] = y; | |
623 | cseed[iLayer].fYref[1] = dy; | |
624 | cseed[iLayer].fZref[0] = z; | |
625 | cseed[iLayer].fZref[1] = dz; | |
626 | cseed[iLayer].fC = curvature; | |
627 | ||
628 | } | |
629 | ||
630 | return chi2TR; | |
631 | ||
632 | } | |
633 | ||
634 | //___________________________________________________________________ | |
635 | void AliTRDseedV1::Print() | |
636 | { | |
637 | // | |
638 | // Printing the seedstatus | |
639 | // | |
640 | ||
641 | AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); | |
642 | Int_t nTimeBins = cal->GetNumberOfTimeBins(); | |
643 | ||
644 | printf("Seed status :\n"); | |
645 | printf(" fTilt = %f\n", fTilt); | |
646 | printf(" fPadLength = %f\n", fPadLength); | |
647 | printf(" fX0 = %f\n", fX0); | |
648 | for(int ic=0; ic<nTimeBins; ic++) { | |
649 | const Char_t *isUsable = fUsable[ic]?"Yes":"No"; | |
650 | printf(" %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%#x] usable[%s]\n" | |
651 | , ic | |
652 | , fX[ic] | |
653 | , fY[ic] | |
654 | , fZ[ic] | |
655 | , fIndexes[ic] | |
656 | , ((Int_t) fClusters[ic]) | |
657 | , isUsable); | |
658 | } | |
659 | ||
660 | printf(" fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]); | |
661 | printf(" fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]); | |
662 | printf(" fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]); | |
663 | printf(" fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]); | |
664 | printf(" fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]); | |
665 | printf(" fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]); | |
666 | printf(" fSigmaY =%f\n", fSigmaY); | |
667 | printf(" fSigmaY2=%f\n", fSigmaY2); | |
668 | printf(" fMeanz =%f\n", fMeanz); | |
669 | printf(" fZProb =%f\n", fZProb); | |
670 | printf(" fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]); | |
671 | printf(" fN =%d\n", fN); | |
672 | printf(" fN2 =%d (>8 isOK)\n",fN2); | |
673 | printf(" fNUsed =%d\n", fNUsed); | |
674 | printf(" fFreq =%d\n", fFreq); | |
675 | printf(" fNChange=%d\n", fNChange); | |
676 | printf(" fMPads =%f\n", fMPads); | |
677 | ||
678 | printf(" fC =%f\n", fC); | |
679 | printf(" fCC =%f\n",fCC); | |
680 | printf(" fChi2 =%f\n", fChi2); | |
681 | printf(" fChi2Z =%f\n", fChi2Z); | |
682 | ||
683 | } |