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 | } |