Line data Source code
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 :
17 :
18 :
19 : //-----------------------------------------------------------------
20 : //
21 : // Implementation of the TPC seed class
22 : // This class is used by the AliTPCtracker class
23 : // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
24 : //-----------------------------------------------------------------
25 : #include <TVectorF.h>
26 : #include "TClonesArray.h"
27 : #include "TGraphErrors.h"
28 : #include "AliTPCseed.h"
29 : #include "AliTPCReconstructor.h"
30 : #include "AliTPCtracker.h"
31 : #include "AliTPCClusterParam.h"
32 : #include "AliTPCCalPad.h"
33 : #include "AliTPCCalROC.h"
34 : #include "AliTPCcalibDB.h"
35 : #include "AliTPCParam.h"
36 : #include "AliMathBase.h"
37 : #include "AliTPCTransform.h"
38 : #include "AliSplineFit.h"
39 : #include "AliCDBManager.h"
40 : #include "AliTPCcalibDButil.h"
41 : #include <AliCTPTimeParams.h>
42 :
43 :
44 16 : ClassImp(AliTPCseed)
45 :
46 :
47 :
48 14838 : AliTPCseed::AliTPCseed():
49 14838 : AliTPCtrack(),
50 14838 : fEsd(0x0),
51 14838 : fNClStore(0),
52 14838 : fClusterPointer(0),
53 14838 : fClusterOwner(kFALSE),
54 14838 : fRow(-1),
55 14838 : fSector(-1),
56 14838 : fRelativeSector(-1),
57 14838 : fCurrentSigmaY2(1e10),
58 14838 : fCurrentSigmaZ2(1e10),
59 14838 : fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30%
60 14838 : fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30%
61 14838 : fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2%
62 14838 : fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2%
63 : //
64 14838 : fErrorY2(1e10),
65 14838 : fErrorZ2(1e10),
66 14838 : fCurrentCluster(0x0),
67 14838 : fCurrentClusterIndex1(-1),
68 14838 : fNoCluster(0),
69 14838 : fSort(0),
70 14838 : fSeedType(0),
71 14838 : fSeed1(-1),
72 14838 : fSeed2(-1),
73 14838 : fCircular(0),
74 14838 : fMAngular(0),
75 14838 : fPoolID(-1),
76 14838 : fTrackPointsArr()
77 89028 : {
78 : //
79 4748160 : for (Int_t i=0;i<kMaxRow;i++) SetClusterIndex2(i,-3);
80 118704 : for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
81 178056 : for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0.2;
82 148380 : for (Int_t i=0;i<4;i++) {
83 59352 : fDEDX[i] = 0.;
84 59352 : fSDEDX[i] = 1e10;
85 59352 : fNCDEDX[i] = 0;
86 59352 : fNCDEDXInclThres[i] = 0;
87 : }
88 296760 : for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
89 385788 : for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
90 29676 : }
91 :
92 164 : AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner):
93 164 : AliTPCtrack(s),
94 164 : fEsd(0x0),
95 164 : fNClStore(0),
96 164 : fClusterPointer(0),
97 164 : fClusterOwner(clusterOwner),
98 164 : fRow(-1),
99 164 : fSector(-1),
100 164 : fRelativeSector(-1),
101 164 : fCurrentSigmaY2(-1),
102 164 : fCurrentSigmaZ2(-1),
103 164 : fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30%
104 164 : fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30%
105 164 : fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2%
106 164 : fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2%
107 164 : fErrorY2(1e10),
108 164 : fErrorZ2(1e10),
109 164 : fCurrentCluster(0x0),
110 164 : fCurrentClusterIndex1(-1),
111 164 : fNoCluster(0),
112 164 : fSort(0),
113 164 : fSeedType(0),
114 164 : fSeed1(-1),
115 164 : fSeed2(-1),
116 164 : fCircular(0),
117 164 : fMAngular(0),
118 164 : fPoolID(-1),
119 164 : fTrackPointsArr(s.fTrackPointsArr)
120 656 : {
121 : //---------------------
122 : // dummy copy constructor
123 : //-------------------------
124 164 : if (s.fClusterPointer) {
125 0 : fNClStore = kMaxRow;
126 0 : fClusterPointer = new AliTPCclusterMI*[kMaxRow];
127 0 : for (Int_t i=0;i<kMaxRow;i++) {
128 0 : fClusterPointer[i]=0;
129 0 : if (fClusterOwner){
130 0 : if (s.fClusterPointer[i])
131 0 : fClusterPointer[i] = new AliTPCclusterMI(*(s.fClusterPointer[i]));
132 : }else{
133 0 : fClusterPointer[i] = s.fClusterPointer[i];
134 : }
135 : }
136 0 : }
137 52480 : for (Int_t i=0;i<kMaxRow;i++) fIndex[i] = s.fIndex[i];
138 1968 : for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=s.fTPCr[i];
139 1640 : for (Int_t i=0;i<4;i++) {
140 656 : fDEDX[i] = s.fDEDX[i];
141 656 : fSDEDX[i] = s.fSDEDX[i];
142 656 : fNCDEDX[i] = s.fNCDEDX[i];
143 656 : fNCDEDXInclThres[i] = s.fNCDEDXInclThres[i];
144 : }
145 3280 : for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
146 :
147 4264 : for (Int_t i=0;i<12;i++) fOverlapLabels[i] = s.fOverlapLabels[i];
148 328 : }
149 :
150 :
151 272 : AliTPCseed::AliTPCseed(const AliTPCtrack &t):
152 272 : AliTPCtrack(t),
153 272 : fEsd(0x0),
154 272 : fNClStore(0),
155 272 : fClusterPointer(0),
156 272 : fClusterOwner(kFALSE),
157 272 : fRow(-1),
158 272 : fSector(-1),
159 272 : fRelativeSector(-1),
160 272 : fCurrentSigmaY2(-1),
161 272 : fCurrentSigmaZ2(-1),
162 272 : fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30%
163 272 : fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30%
164 272 : fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2%
165 272 : fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2%
166 272 : fErrorY2(1e10),
167 272 : fErrorZ2(1e10),
168 272 : fCurrentCluster(0x0),
169 272 : fCurrentClusterIndex1(-1),
170 272 : fNoCluster(0),
171 272 : fSort(0),
172 272 : fSeedType(0),
173 272 : fSeed1(-1),
174 272 : fSeed2(-1),
175 272 : fCircular(0),
176 272 : fMAngular(0),
177 272 : fPoolID(-1),
178 272 : fTrackPointsArr()
179 1632 : {
180 : //
181 : // Constructor from AliTPCtrack
182 : //
183 272 : fFirstPoint =0;
184 3264 : for (Int_t i=0;i<5;i++) fTPCr[i]=0.2;
185 87040 : for (Int_t i=0;i<kMaxRow;i++) {
186 43248 : Int_t index = t.GetClusterIndex(i);
187 43248 : if (index>=-1){
188 35490 : SetClusterIndex2(i,index);
189 35490 : }
190 : else{
191 7758 : SetClusterIndex2(i,-3);
192 : }
193 : }
194 2720 : for (Int_t i=0;i<4;i++) {
195 1088 : fDEDX[i] = 0.;
196 1088 : fSDEDX[i] = 1e10;
197 1088 : fNCDEDX[i] = 0;
198 1088 : fNCDEDXInclThres[i] = 0;
199 : }
200 5440 : for (Int_t i=0;i<9;i++) fDEDX[i] = fDEDX[i];
201 :
202 7072 : for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
203 544 : }
204 :
205 700 : AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5],
206 : const Double_t cc[15], Int_t index):
207 700 : AliTPCtrack(xr, alpha, xx, cc, index),
208 700 : fEsd(0x0),
209 700 : fNClStore(0),
210 700 : fClusterPointer(0),
211 700 : fClusterOwner(kFALSE),
212 700 : fRow(-1),
213 700 : fSector(-1),
214 700 : fRelativeSector(-1),
215 700 : fCurrentSigmaY2(-1),
216 700 : fCurrentSigmaZ2(-1),
217 700 : fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30%
218 700 : fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30%
219 700 : fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2%
220 700 : fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2%
221 700 : fErrorY2(1e10),
222 700 : fErrorZ2(1e10),
223 700 : fCurrentCluster(0x0),
224 700 : fCurrentClusterIndex1(-1),
225 700 : fNoCluster(0),
226 700 : fSort(0),
227 700 : fSeedType(0),
228 700 : fSeed1(-1),
229 700 : fSeed2(-1),
230 700 : fCircular(0),
231 700 : fMAngular(0),
232 700 : fPoolID(-1),
233 700 : fTrackPointsArr()
234 4200 : {
235 : //
236 : // Constructor
237 : //
238 700 : fFirstPoint =0;
239 224000 : for (Int_t i=0;i<kMaxRow;i++) SetClusterIndex2(i,-3);
240 8400 : for (Int_t i=0;i<5;i++) fTPCr[i]=0.2;
241 7000 : for (Int_t i=0;i<4;i++) {
242 2800 : fDEDX[i] = 0.;
243 2800 : fSDEDX[i] = 1e10;
244 2800 : fNCDEDX[i] = 0;
245 2800 : fNCDEDXInclThres[i] = 0;
246 : }
247 14000 : for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
248 :
249 18200 : for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
250 1400 : }
251 :
252 89596 : AliTPCseed::~AliTPCseed(){
253 : //
254 : // destructor
255 78160 : AliDebug(5,"Destruct AliTPCseed - Begin");
256 15632 : fNoCluster =0;
257 15632 : if (fClusterPointer) {
258 0 : if (fClusterOwner){
259 0 : for (Int_t icluster=0; icluster<kMaxRow; icluster++){
260 0 : if(fClusterPointer[icluster]) { // shared clusters might be already deleted!
261 0 : if (fClusterPointer[icluster]->TestBit(TObject::kNotDeleted)) delete fClusterPointer[icluster];
262 0 : fClusterPointer[icluster]=0;
263 0 : }
264 : }
265 0 : }
266 0 : delete[] fClusterPointer;
267 : }
268 78160 : AliDebug(5,"Destruct AliTPCseed - End");
269 44798 : }
270 : //_________________________________________________
271 : AliTPCseed & AliTPCseed::operator=(const AliTPCseed ¶m)
272 : {
273 : //
274 : // assignment operator
275 : // don't touch pool ID
276 : //
277 8472 : if(this!=¶m){
278 4236 : AliTPCtrack::operator=(param);
279 4236 : fTrackPointsArr = param.fTrackPointsArr;
280 4236 : fEsd =param.fEsd;
281 4236 : fClusterOwner = param.fClusterOwner;
282 4236 : if (param.fClusterPointer) {
283 0 : if (!fClusterPointer) {
284 0 : fClusterPointer = new AliTPCclusterMI*[kMaxRow];
285 0 : memset(fClusterPointer,0,kMaxRow*sizeof(AliTPCclusterMI*));
286 0 : fNClStore = kMaxRow;
287 0 : }
288 0 : if (!fClusterOwner) for(Int_t i = 0;i<kMaxRow;++i) fClusterPointer[i] = param.fClusterPointer[i];
289 0 : else for(Int_t i = 0;i<kMaxRow;++i) {
290 0 : delete fClusterPointer[i];
291 0 : if (param.fClusterPointer[i]) fClusterPointer[i] = new AliTPCclusterMI(*(param.fClusterPointer[i]));
292 0 : else fClusterPointer[i] = 0x0;
293 : }
294 : }
295 : // leave out fPoint, they are also not copied in the copy ctor...
296 : // but deleted in the dtor... strange...
297 4236 : fRow = param.fRow;
298 4236 : fSector = param.fSector;
299 4236 : fRelativeSector = param.fRelativeSector;
300 4236 : fCurrentSigmaY2 = param.fCurrentSigmaY2;
301 4236 : fCurrentSigmaZ2 = param.fCurrentSigmaZ2;
302 4236 : fCMeanSigmaY2p30 = param.fCMeanSigmaY2p30;
303 4236 : fCMeanSigmaZ2p30 = param.fCMeanSigmaZ2p30;
304 4236 : fCMeanSigmaY2p30R = param.fCMeanSigmaY2p30R;
305 4236 : fCMeanSigmaZ2p30R = param.fCMeanSigmaZ2p30R;
306 4236 : fErrorY2 = param.fErrorY2;
307 4236 : fErrorZ2 = param.fErrorZ2;
308 4236 : fCurrentCluster = param.fCurrentCluster; // this is not allocated by AliTPCSeed
309 4236 : fCurrentClusterIndex1 = param.fCurrentClusterIndex1;
310 4236 : fNoCluster = param.fNoCluster;
311 4236 : fSort = param.fSort;
312 42360 : for(Int_t i = 0;i<4;++i){
313 16944 : fSDEDX[i] = param.fSDEDX[i];
314 16944 : fNCDEDX[i] = param.fNCDEDX[i];
315 16944 : fNCDEDXInclThres[i] = param.fNCDEDXInclThres[i];
316 : }
317 :
318 50832 : for(Int_t i = 0;i<AliPID::kSPECIES;++i)fTPCr[i] = param.fTPCr[i];
319 :
320 4236 : fSeedType = param.fSeedType;
321 4236 : fSeed1 = param.fSeed1;
322 4236 : fSeed2 = param.fSeed2;
323 110136 : for(Int_t i = 0;i<12;++i)fOverlapLabels[i] = param.fOverlapLabels[i];
324 4236 : fMAngular = param.fMAngular;
325 4236 : fCircular = param.fCircular;
326 4236 : }
327 4236 : return (*this);
328 0 : }
329 :
330 : Double_t AliTPCseed::GetDensityFirst(Int_t n)
331 : {
332 : //
333 : //
334 : // return cluster for n rows bellow first point
335 : Int_t nfoundable = 1;
336 : Int_t nfound = 1;
337 110054 : for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
338 35588 : Int_t index = GetClusterIndex2(i);
339 67688 : if (index!=-1) nfoundable++;
340 64414 : if (index>0) nfound++;
341 : }
342 862 : if (nfoundable<n) return 0;
343 802 : return Double_t(nfound)/Double_t(nfoundable);
344 :
345 832 : }
346 :
347 :
348 : void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
349 : {
350 : // get cluster stat. on given region
351 : //
352 0 : found = 0;
353 0 : foundable = 0;
354 0 : shared =0;
355 0 : for (Int_t i=first;i<last; i++){
356 0 : Int_t index = GetClusterIndex2(i);
357 0 : if (index!=-1) foundable++;
358 0 : if (index&0x8000) continue;
359 0 : if (fClusterPointer && fClusterPointer[i]) {
360 0 : found++;
361 : }
362 : else
363 0 : continue;
364 :
365 0 : if (fClusterPointer && fClusterPointer[i]->IsUsed(10)) {
366 0 : shared++;
367 0 : continue;
368 : }
369 0 : if (!plus2) continue; //take also neighborhoud
370 : //
371 0 : if ( (i>0) && fClusterPointer && fClusterPointer[i-1]){
372 0 : if (fClusterPointer[i-1]->IsUsed(10)) {
373 0 : shared++;
374 0 : continue;
375 : }
376 : }
377 0 : if ( fClusterPointer && fClusterPointer[i+1]){
378 0 : if (fClusterPointer[i+1]->IsUsed(10)) {
379 0 : shared++;
380 0 : continue;
381 : }
382 : }
383 :
384 0 : }
385 : //if (shared>found){
386 : //Error("AliTPCseed::GetClusterStatistic","problem\n");
387 : //}
388 0 : }
389 :
390 :
391 : void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable)
392 : {
393 : // get cluster stat. on given region, faster than the version providing also shared info
394 : //
395 96 : found = foundable = 0;
396 6968 : for (Int_t i=first;i<last; i++){
397 3436 : Int_t index = GetClusterIndex2(i);
398 6730 : if (index!=-1) foundable++;
399 6926 : if (index<0 || index&0x8000) continue;
400 1700 : found++;
401 1700 : }
402 48 : }
403 :
404 : void AliTPCseed::Reset(Bool_t all)
405 : {
406 : //
407 : //
408 428 : SetNumberOfClusters(0);
409 428 : fNFoundable = 0;
410 428 : SetChi2(0);
411 428 : ResetCovariance(10.);
412 :
413 428 : if (all){
414 53760 : for (Int_t i=kMaxRow;i--;) SetClusterIndex2(i,-3);
415 168 : if (fClusterPointer) {
416 0 : if (!fClusterOwner) for (Int_t i=kMaxRow;i--;) fClusterPointer[i]=0;
417 0 : else for (Int_t i=kMaxRow;i--;) {delete fClusterPointer[i]; fClusterPointer[i]=0;}
418 : }
419 : }
420 :
421 428 : }
422 :
423 :
424 : void AliTPCseed::Modify(Double_t factor)
425 : {
426 :
427 : //------------------------------------------------------------------
428 : //This function makes a track forget its history :)
429 : //------------------------------------------------------------------
430 0 : if (factor<=0) {
431 0 : ResetCovariance(10.);
432 0 : return;
433 : }
434 0 : ResetCovariance(factor);
435 :
436 0 : SetNumberOfClusters(0);
437 0 : fNFoundable =0;
438 0 : SetChi2(0);
439 0 : fRemoval = 0;
440 0 : fCurrentSigmaY2 = 0.000005;
441 0 : fCurrentSigmaZ2 = 0.000005;
442 0 : fNoCluster = 0;
443 : //fFirstPoint = kMaxRow;
444 : //fLastPoint = 0;
445 0 : }
446 :
447 :
448 :
449 :
450 : Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
451 : {
452 : //-----------------------------------------------------------------
453 : // This function find proloncation of a track to a reference plane x=xk.
454 : // doesn't change internal state of the track
455 : //-----------------------------------------------------------------
456 :
457 141480 : Double_t x1=GetX(), x2=x1+(xk-x1), dx=x2-x1;
458 :
459 70740 : if (TMath::Abs(GetSnp()+GetC()*dx) >= AliTPCReconstructor::GetMaxSnpTrack()) {
460 36 : return 0;
461 : }
462 :
463 : // Double_t y1=fP0, z1=fP1;
464 70704 : Double_t c1=GetSnp(), r1=sqrt((1.-c1)*(1.+c1));
465 70704 : Double_t c2=c1 + GetC()*dx, r2=sqrt((1.-c2)*(1.+c2));
466 :
467 70704 : y = GetY();
468 70704 : z = GetZ();
469 : //y += dx*(c1+c2)/(r1+r2);
470 : //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
471 :
472 70704 : Double_t dy = dx*(c1+c2)/(r1+r2);
473 : Double_t dz = 0;
474 : //
475 70704 : Double_t delta = GetC()*dx*(c1+c2)/(c1*r2 + c2*r1);
476 : /*
477 : if (TMath::Abs(delta)>0.0001){
478 : dz = fP3*TMath::ASin(delta)/fP4;
479 : }else{
480 : dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
481 : }
482 : */
483 : // dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4;
484 70704 : dz = GetTgl()*TMath::ASin(delta)/GetC();
485 : //
486 70704 : y+=dy;
487 70704 : z+=dz;
488 :
489 :
490 : return 1;
491 70740 : }
492 :
493 :
494 : //_____________________________________________________________________________
495 : Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const
496 : {
497 : //-----------------------------------------------------------------
498 : // This function calculates a predicted chi2 increment.
499 : //-----------------------------------------------------------------
500 198496 : Double_t p[2]={c->GetY(), c->GetZ()};
501 99248 : Double_t cov[3]={fErrorY2, 0., fErrorZ2};
502 :
503 99248 : Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
504 99248 : if (TMath::Abs(dx)>0){
505 66766 : Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
506 66766 : Float_t dy = dx*ty;
507 66766 : Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
508 66766 : p[0] = c->GetY()-dy;
509 66766 : p[1] = c->GetZ()-dz;
510 66766 : }
511 198496 : return AliExternalTrackParam::GetPredictedChi2(p,cov);
512 99248 : }
513 :
514 : //_________________________________________________________________________________________
515 :
516 :
517 : Int_t AliTPCseed::Compare(const TObject *o) const {
518 : //-----------------------------------------------------------------
519 : // This function compares tracks according to the sector - for given sector according z
520 : //-----------------------------------------------------------------
521 3544 : AliTPCseed *t=(AliTPCseed*)o;
522 :
523 1772 : if (fSort == 0){
524 1610 : if (t->fRelativeSector>fRelativeSector) return -1;
525 1114 : if (t->fRelativeSector<fRelativeSector) return 1;
526 266 : Double_t z2 = t->GetZ();
527 266 : Double_t z1 = GetZ();
528 414 : if (z2>z1) return 1;
529 236 : if (z2<z1) return -1;
530 0 : return 0;
531 : }
532 : else {
533 : Float_t f2 =1;
534 622 : f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(t->OneOverPt()+0.0066);
535 622 : if (t->fBConstrain) f2=1.2;
536 :
537 : Float_t f1 =1;
538 622 : f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066);
539 :
540 622 : if (fBConstrain) f1=1.2;
541 :
542 930 : if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
543 314 : else return +1;
544 : }
545 1772 : }
546 :
547 :
548 :
549 :
550 : //_____________________________________________________________________________
551 : Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index)
552 : {
553 : //-----------------------------------------------------------------
554 : // This function associates a cluster with this track.
555 : //-----------------------------------------------------------------
556 198496 : Int_t n=GetNumberOfClusters();
557 99248 : Int_t idx=GetClusterIndex(n); // save the current cluster index
558 :
559 99248 : AliTPCclusterMI cl(*(AliTPCclusterMI*)c); cl.SetSigmaY2(fErrorY2); cl.SetSigmaZ2(fErrorZ2);
560 :
561 198496 : AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
562 :
563 297744 : Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
564 :
565 99248 : if( parcl ){
566 : Int_t padSize = 0; // short pads
567 99248 : if (cl.GetDetector() >= 36) {
568 : padSize = 1; // medium pads
569 83652 : if (cl.GetRow() > 63) padSize = 2; // long pads
570 : }
571 99248 : Float_t waveCorr = parcl->GetWaveCorrection( padSize, cl.GetZ(), cl.GetMax(),cl.GetPad(), ty );
572 99248 : cl.SetY( cl.GetY() - waveCorr );
573 99248 : }
574 :
575 198496 : Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
576 99248 : if (TMath::Abs(dx)>0){
577 66766 : Float_t dy = dx*ty;
578 133532 : Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
579 66766 : cl.SetY(cl.GetY()-dy);
580 66766 : cl.SetZ(cl.GetZ()-dz);
581 66766 : }
582 :
583 :
584 198498 : if (!AliTPCtrack::Update(&cl,chisq,index)) return kFALSE;
585 :
586 99246 : if (fCMeanSigmaY2p30<0){
587 1094 : fCMeanSigmaY2p30= c->GetSigmaY2(); //! current mean sigma Y2 - mean30%
588 1094 : fCMeanSigmaZ2p30= c->GetSigmaZ2(); //! current mean sigma Z2 - mean30%
589 1094 : fCMeanSigmaY2p30R = 1; //! current mean sigma Y2 - mean5%
590 1094 : fCMeanSigmaZ2p30R = 1; //! current mean sigma Z2 - mean5%
591 1094 : }
592 : //
593 99246 : fCMeanSigmaY2p30= 0.70*fCMeanSigmaY2p30 +0.30*c->GetSigmaY2();
594 99246 : fCMeanSigmaZ2p30= 0.70*fCMeanSigmaZ2p30 +0.30*c->GetSigmaZ2();
595 99246 : if (fCurrentSigmaY2>0){
596 99246 : fCMeanSigmaY2p30R = 0.7*fCMeanSigmaY2p30R +0.3*c->GetSigmaY2()/fCurrentSigmaY2;
597 99246 : fCMeanSigmaZ2p30R = 0.7*fCMeanSigmaZ2p30R +0.3*c->GetSigmaZ2()/fCurrentSigmaZ2;
598 99246 : }
599 :
600 :
601 99246 : SetClusterIndex(n,idx); // restore the current cluster index
602 99246 : return kTRUE;
603 99248 : }
604 :
605 :
606 :
607 : //_____________________________________________________________________________
608 : Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t /* onlyused */) {
609 : //-----------------------------------------------------------------
610 : // This funtion calculates dE/dX within the "low" and "up" cuts.
611 : //-----------------------------------------------------------------
612 : // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal)
613 1080 : AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
614 :
615 540 : Int_t row0 = param->GetNRowLow();
616 540 : Int_t row1 = row0+param->GetNRowUp1();
617 540 : Int_t row2 = row1+param->GetNRowUp2();
618 540 : const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
619 : Int_t useTot = 0;
620 1080 : if (recoParam) useTot = (recoParam->GetUseTotCharge())? 0:1;
621 : //
622 : //
623 540 : TVectorF i1i2;
624 540 : TVectorF irocTot;
625 540 : TVectorF oroc1Tot;
626 540 : TVectorF oroc2Tot;
627 540 : TVectorF forocTot;
628 : //
629 540 : TVectorF irocMax;
630 540 : TVectorF oroc1Max;
631 540 : TVectorF oroc2Max;
632 540 : TVectorF forocMax;
633 :
634 540 : CookdEdxAnalytical(low,up,useTot ,i1 ,i2, 0, 2, 0, &i1i2);
635 : //
636 540 : CookdEdxAnalytical(low,up,kTRUE ,0 ,row0, 0, 2, 0, &irocTot);
637 540 : CookdEdxAnalytical(low,up,kTRUE ,row0,row1, 0, 2, 0, &oroc1Tot);
638 540 : CookdEdxAnalytical(low,up,kTRUE ,row1,row2, 0, 2, 0, &oroc2Tot);
639 540 : CookdEdxAnalytical(low,up,kTRUE ,row0,row2, 0, 2, 0, &forocTot); // full OROC truncated mean
640 : //
641 540 : CookdEdxAnalytical(low,up,kFALSE ,0 ,row0, 0, 2, 0, &irocMax);
642 540 : CookdEdxAnalytical(low,up,kFALSE ,row0,row1, 0, 2, 0, &oroc1Max);
643 540 : CookdEdxAnalytical(low,up,kFALSE ,row1,row2, 0, 2, 0, &oroc2Max);
644 540 : CookdEdxAnalytical(low,up,kFALSE ,row0,row2, 0, 2, 0, &forocMax); // full OROC truncated mean
645 :
646 1080 : fDEDX[0] = i1i2(0);
647 : //
648 1080 : fDEDX[1] = irocTot(0);
649 1080 : fDEDX[2] = oroc1Tot(0);
650 1080 : fDEDX[3] = oroc2Tot(0);
651 1080 : fDEDX[4] = forocTot(0); // full OROC truncated mean
652 1080 : fDEDX[5] = irocMax(0);
653 1080 : fDEDX[6] = oroc1Max(0);
654 1080 : fDEDX[7] = oroc2Max(0);
655 1080 : fDEDX[8] = forocMax(0); // full OROC truncated mean
656 : //
657 1080 : fSDEDX[0] = i1i2(1);
658 1080 : fSDEDX[1] = irocTot(1);
659 1080 : fSDEDX[2] = oroc1Tot(1);
660 1080 : fSDEDX[3] = oroc2Tot(1);
661 : //
662 1080 : fNCDEDX[0] = TMath::Nint(i1i2(2));
663 :
664 1080 : fNCDEDX[1] = TMath::Nint( irocTot(2));
665 1080 : fNCDEDX[2] = TMath::Nint(oroc1Tot(2));
666 1080 : fNCDEDX[3] = TMath::Nint(oroc2Tot(2));
667 : //
668 1620 : fNCDEDXInclThres[0] = TMath::Nint(i1i2(2)+i1i2(9));
669 1620 : fNCDEDXInclThres[1] = TMath::Nint( irocTot(2)+ irocTot(9));
670 1620 : fNCDEDXInclThres[2] = TMath::Nint(oroc1Tot(2)+oroc1Tot(9));
671 1620 : fNCDEDXInclThres[3] = TMath::Nint(oroc2Tot(2)+oroc2Tot(9));
672 : //
673 540 : SetdEdx(fDEDX[0]);
674 540 : return fDEDX[0];
675 :
676 : // return CookdEdxNorm(low,up,0,i1,i2,1,0,2);
677 :
678 :
679 : // Float_t amp[200];
680 : // Float_t angular[200];
681 : // Float_t weight[200];
682 : // Int_t index[200];
683 : // //Int_t nc = 0;
684 : // Float_t meanlog = 100.;
685 :
686 : // Float_t mean[4] = {0,0,0,0};
687 : // Float_t sigma[4] = {1000,1000,1000,1000};
688 : // Int_t nc[4] = {0,0,0,0};
689 : // Float_t norm[4] = {1000,1000,1000,1000};
690 : // //
691 : // //
692 : // fNShared =0;
693 :
694 : // Float_t gainGG = 1;
695 : // if (AliTPCcalibDB::Instance()->GetParameters()){
696 : // gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000.; //relative gas gain
697 : // }
698 :
699 :
700 : // for (Int_t of =0; of<4; of++){
701 : // for (Int_t i=of+i1;i<i2;i+=4)
702 : // {
703 : // Int_t clindex = fIndex[i];
704 : // if (clindex<0||clindex&0x8000) continue;
705 :
706 : // //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
707 : // AliTPCTrackerPoint * point = GetTrackPoint(i);
708 : // //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
709 : // //AliTPCTrackerPoint * pointp = 0;
710 : // //if (i<kMaxRow) pointp = GetTrackPoint(i+1);
711 :
712 : // if (point==0) continue;
713 : // AliTPCclusterMI * cl = fClusterPointer[i];
714 : // if (cl==0) continue;
715 : // if (onlyused && (!cl->IsUsed(10))) continue;
716 : // if (cl->IsUsed(11)) {
717 : // fNShared++;
718 : // continue;
719 : // }
720 : // Int_t type = cl->GetType();
721 : // //if (point->fIsShared){
722 : // // fNShared++;
723 : // // continue;
724 : // //}
725 : // //if (pointm)
726 : // // if (pointm->fIsShared) continue;
727 : // //if (pointp)
728 : // // if (pointp->fIsShared) continue;
729 :
730 : // if (type<0) continue;
731 : // //if (type>10) continue;
732 : // //if (point->GetErrY()==0) continue;
733 : // //if (point->GetErrZ()==0) continue;
734 :
735 : // //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
736 : // //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
737 : // //if ((ddy*ddy+ddz*ddz)>10) continue;
738 :
739 :
740 : // // if (point->GetCPoint().GetMax()<5) continue;
741 : // if (cl->GetMax()<5) continue;
742 : // Float_t angley = point->GetAngleY();
743 : // Float_t anglez = point->GetAngleZ();
744 :
745 : // Float_t rsigmay2 = point->GetSigmaY();
746 : // Float_t rsigmaz2 = point->GetSigmaZ();
747 : // /*
748 : // Float_t ns = 1.;
749 : // if (pointm){
750 : // rsigmay += pointm->GetTPoint().GetSigmaY();
751 : // rsigmaz += pointm->GetTPoint().GetSigmaZ();
752 : // ns+=1.;
753 : // }
754 : // if (pointp){
755 : // rsigmay += pointp->GetTPoint().GetSigmaY();
756 : // rsigmaz += pointp->GetTPoint().GetSigmaZ();
757 : // ns+=1.;
758 : // }
759 : // rsigmay/=ns;
760 : // rsigmaz/=ns;
761 : // */
762 :
763 : // Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
764 :
765 : // Float_t ampc = 0; // normalization to the number of electrons
766 : // if (i>64){
767 : // // ampc = 1.*point->GetCPoint().GetMax();
768 : // ampc = 1.*cl->GetMax();
769 : // //ampc = 1.*point->GetCPoint().GetQ();
770 : // // AliTPCClusterPoint & p = point->GetCPoint();
771 : // // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
772 : // // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
773 : // //Float_t dz =
774 : // // TMath::Abs( Int_t(iz) - iz + 0.5);
775 : // //ampc *= 1.15*(1-0.3*dy);
776 : // //ampc *= 1.15*(1-0.3*dz);
777 : // // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
778 : // //ampc *=zfactor;
779 : // }
780 : // else{
781 : // //ampc = 1.0*point->GetCPoint().GetMax();
782 : // ampc = 1.0*cl->GetMax();
783 : // //ampc = 1.0*point->GetCPoint().GetQ();
784 : // //AliTPCClusterPoint & p = point->GetCPoint();
785 : // // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
786 : // //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
787 : // //Float_t dz =
788 : // // TMath::Abs( Int_t(iz) - iz + 0.5);
789 :
790 : // //ampc *= 1.15*(1-0.3*dy);
791 : // //ampc *= 1.15*(1-0.3*dz);
792 : // // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
793 : // //ampc *=zfactor;
794 :
795 : // }
796 : // ampc *= 2.0; // put mean value to channel 50
797 : // //ampc *= 0.58; // put mean value to channel 50
798 : // Float_t w = 1.;
799 : // // if (type>0) w = 1./(type/2.-0.5);
800 : // // Float_t z = TMath::Abs(cl->GetZ());
801 : // if (i<64) {
802 : // ampc /= 0.6;
803 : // //ampc /= (1+0.0008*z);
804 : // } else
805 : // if (i>128){
806 : // ampc /=1.5;
807 : // //ampc /= (1+0.0008*z);
808 : // }else{
809 : // //ampc /= (1+0.0008*z);
810 : // }
811 :
812 : // if (type<0) { //amp at the border - lower weight
813 : // // w*= 2.;
814 :
815 : // continue;
816 : // }
817 : // if (rsigma>1.5) ampc/=1.3; // if big backround
818 : // amp[nc[of]] = ampc;
819 : // amp[nc[of]] /=gainGG;
820 : // angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
821 : // weight[nc[of]] = w;
822 : // nc[of]++;
823 : // }
824 :
825 : // TMath::Sort(nc[of],amp,index,kFALSE);
826 : // Float_t sumamp=0;
827 : // Float_t sumamp2=0;
828 : // Float_t sumw=0;
829 : // //meanlog = amp[index[Int_t(nc[of]*0.33)]];
830 : // meanlog = 50;
831 : // for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
832 : // Float_t ampl = amp[index[i]]/angular[index[i]];
833 : // ampl = meanlog*TMath::Log(1.+ampl/meanlog);
834 : // //
835 : // sumw += weight[index[i]];
836 : // sumamp += weight[index[i]]*ampl;
837 : // sumamp2 += weight[index[i]]*ampl*ampl;
838 : // norm[of] += angular[index[i]]*weight[index[i]];
839 : // }
840 : // if (sumw<1){
841 : // SetdEdx(0);
842 : // }
843 : // else {
844 : // norm[of] /= sumw;
845 : // mean[of] = sumamp/sumw;
846 : // sigma[of] = sumamp2/sumw-mean[of]*mean[of];
847 : // if (sigma[of]>0.1)
848 : // sigma[of] = TMath::Sqrt(sigma[of]);
849 : // else
850 : // sigma[of] = 1000;
851 :
852 : // mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
853 : // //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
854 : // //mean *=(1-0.1*(norm-1.));
855 : // }
856 : // }
857 :
858 : // Float_t dedx =0;
859 : // fSdEdx =0;
860 : // fMAngular =0;
861 : // // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
862 : // // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
863 :
864 :
865 : // // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
866 : // // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
867 :
868 : // Int_t norm2 = 0;
869 : // Int_t norm3 = 0;
870 : // for (Int_t i =0;i<4;i++){
871 : // if (nc[i]>2&&nc[i]<1000){
872 : // dedx += mean[i] *nc[i];
873 : // fSdEdx += sigma[i]*(nc[i]-2);
874 : // fMAngular += norm[i] *nc[i];
875 : // norm2 += nc[i];
876 : // norm3 += nc[i]-2;
877 : // }
878 : // fDEDX[i] = mean[i];
879 : // fSDEDX[i] = sigma[i];
880 : // fNCDEDX[i]= nc[i];
881 : // }
882 :
883 : // if (norm3>0){
884 : // dedx /=norm2;
885 : // fSdEdx /=norm3;
886 : // fMAngular/=norm2;
887 : // }
888 : // else{
889 : // SetdEdx(0);
890 : // return 0;
891 : // }
892 : // // Float_t dedx1 =dedx;
893 : // /*
894 : // dedx =0;
895 : // for (Int_t i =0;i<4;i++){
896 : // if (nc[i]>2&&nc[i]<1000){
897 : // mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
898 : // dedx += mean[i] *nc[i];
899 : // }
900 : // fDEDX[i] = mean[i];
901 : // }
902 : // dedx /= norm2;
903 : // */
904 :
905 :
906 : // SetdEdx(dedx);
907 : // return dedx;
908 540 : }
909 :
910 : void AliTPCseed::CookPID()
911 : {
912 : //
913 : // cook PID information according dEdx
914 : //
915 : Double_t fRange = 10.;
916 : Double_t fRes = 0.1;
917 : Double_t fMIP = 47.;
918 : //
919 : Int_t ns=AliPID::kSPECIES;
920 : Double_t sumr =0;
921 1742 : for (Int_t j=0; j<ns; j++) {
922 670 : Double_t mass=AliPID::ParticleMass(j);
923 670 : Double_t mom=GetP();
924 670 : Double_t dedx=fdEdx/fMIP;
925 670 : Double_t bethe=AliMathBase::BetheBlochAleph(mom/mass);
926 670 : Double_t sigma=fRes*bethe;
927 670 : if (sigma>0.001){
928 670 : if (TMath::Abs(dedx-bethe) > fRange*sigma) {
929 0 : fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
930 0 : sumr+=fTPCr[j];
931 0 : continue;
932 : }
933 670 : fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
934 670 : sumr+=fTPCr[j];
935 670 : }
936 : else{
937 0 : fTPCr[j]=1.;
938 0 : sumr+=fTPCr[j];
939 : }
940 670 : }
941 1608 : for (Int_t j=0; j<ns; j++) {
942 670 : fTPCr[j]/=sumr; //normalize
943 : }
944 134 : }
945 :
946 : Double_t AliTPCseed::GetYat(Double_t xk) const {
947 : //-----------------------------------------------------------------
948 : // This function calculates the Y-coordinate of a track at the plane x=xk.
949 : //-----------------------------------------------------------------
950 0 : if (TMath::Abs(GetSnp())>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06
951 0 : Double_t c1=GetSnp(), r1=TMath::Sqrt((1.-c1)*(1.+c1));
952 0 : Double_t c2=c1+GetC()*(xk-GetX());
953 0 : if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0;
954 0 : Double_t r2=TMath::Sqrt((1.-c2)*(1.+c2));
955 0 : return GetY() + (xk-GetX())*(c1+c2)/(r1+r2);
956 0 : }
957 :
958 :
959 :
960 : Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Int_t posNorm, Int_t padNorm, Int_t returnVal){
961 :
962 : //
963 : // calculates dedx using the cluster
964 : // low - up specify trunc mean range - default form 0-0.7
965 : // type - 1 - max charge or 0- total charge in cluster
966 : // //2- max no corr 3- total+ correction
967 : // i1-i2 - the pad-row range used for calculation
968 : // shapeNorm - kTRUE -taken from OCDB
969 : //
970 : // posNorm - usage of pos normalization
971 : // padNorm - pad type normalization
972 : // returnVal - 0 return mean
973 : // - 1 return RMS
974 : // - 2 return number of clusters
975 : //
976 : // normalization parametrization taken from AliTPCClusterParam
977 : //
978 0 : AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
979 0 : AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
980 0 : if (!parcl) return 0;
981 0 : if (!param) return 0;
982 0 : Int_t row0 = param->GetNRowLow();
983 0 : Int_t row1 = row0+param->GetNRowUp1();
984 :
985 0 : Float_t amp[kMaxRow];
986 0 : Int_t indexes[kMaxRow];
987 : Int_t ncl=0;
988 : //
989 : //
990 : Float_t gainGG = 1; // gas gain factor -always enabled
991 : Float_t gainPad = 1; // gain map - used always
992 : Float_t corrShape = 1; // correction due angular effect, diffusion and electron attachment
993 : Float_t corrPos = 1; // local position correction - if posNorm enabled
994 : Float_t corrPadType = 1; // pad type correction - if padNorm enabled
995 : Float_t corrNorm = 1; // normalization factor - set Q to channel 50
996 : //
997 : //
998 : //
999 0 : if (AliTPCcalibDB::Instance()->GetParameters()){
1000 0 : gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain
1001 0 : gainGG *= AliTPCcalibDB::Instance()->GetParameters()->GetNtot()/36.82;//correction for the ionisation
1002 0 : }
1003 :
1004 0 : const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1005 : const Float_t kedgey =3.;
1006 : //
1007 : //
1008 0 : for (Int_t irow=i1; irow<i2; irow++){
1009 0 : AliTPCclusterMI* cluster = GetClusterPointer(irow);
1010 0 : if (!cluster) continue;
1011 0 : if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1012 0 : Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1013 : Int_t ipad= 0;
1014 0 : if (irow>=row0) ipad=1;
1015 0 : if (irow>=row1) ipad=2;
1016 : //
1017 : //
1018 : //
1019 0 : AliTPCCalPad * gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1020 0 : if (gainMap) {
1021 : //
1022 : // Get gainPad - pad by pad calibration
1023 : //
1024 : Float_t factor = 1;
1025 0 : AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1026 0 : if (irow < row0) { // IROC
1027 0 : factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1028 0 : } else { // OROC
1029 0 : factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1030 : }
1031 0 : if (factor>0.5) gainPad=factor;
1032 0 : }
1033 : //
1034 : //do position and angular normalization
1035 : //
1036 0 : if (shapeNorm){
1037 0 : if (type<=1){
1038 : //
1039 0 : const AliTPCTrackerPoints::Point * point = GetTrackPoint(irow);
1040 0 : Float_t ty = TMath::Abs(point->GetAngleY());
1041 0 : Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1042 :
1043 0 : Float_t dr = (250.-TMath::Abs(cluster->GetZ()))/250.;
1044 0 : corrShape = parcl->Qnorm(ipad,type,dr,ty,tz);
1045 0 : }
1046 : }
1047 :
1048 0 : if (posNorm>0){
1049 : //
1050 : // Do position normalization - relative distance to
1051 : // center of pad- time bin
1052 : // Work in progress
1053 : // corrPos = parcl->QnormPos(ipad,type, cluster->GetPad(),
1054 : // cluster->GetTimeBin(), cluster->GetZ(),
1055 : // cluster->GetSigmaY2(),cluster->GetSigmaZ2(),
1056 : // cluster->GetMax(),cluster->GetQ());
1057 : // scaled response function
1058 0 : Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1059 0 : Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1060 : //
1061 :
1062 0 : const AliTPCTrackerPoints::Point * point = GetTrackPoint(irow);
1063 0 : Float_t ty = TMath::Abs(point->GetAngleY());
1064 0 : Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1065 :
1066 0 : if (type==1) corrPos =
1067 0 : parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1068 0 : cluster->GetTimeBin(),ty,tz,yres0,zres0,0.4);
1069 0 : if (type==0) corrPos =
1070 0 : parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1071 0 : cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,0.4);
1072 0 : if (posNorm==3){
1073 0 : Float_t dr = (250.-TMath::Abs(cluster->GetZ()))/250.;
1074 0 : Double_t signtgl = (cluster->GetZ()*point->GetAngleZ()>0)? 1:-1;
1075 0 : Double_t p2 = TMath::Abs(TMath::Sin(TMath::ATan(ty)));
1076 0 : Float_t corrHis = parcl->QnormHis(ipad,type,dr,p2,TMath::Abs(point->GetAngleZ())*signtgl);
1077 0 : if (corrHis>0) corrPos*=corrHis;
1078 0 : }
1079 :
1080 0 : }
1081 :
1082 0 : if (padNorm==1){
1083 : //taken from OCDB
1084 0 : if (type==0 && parcl->QpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad];
1085 0 : if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[ipad];
1086 :
1087 : }
1088 0 : if (padNorm==2){
1089 0 : corrPadType =param->GetPadPitchLength(cluster->GetDetector(),cluster->GetRow());
1090 : //use hardwired - temp fix
1091 0 : if (type==0) corrNorm=3.;
1092 0 : if (type==1) corrNorm=1.;
1093 : }
1094 : //
1095 : //
1096 0 : amp[ncl]=charge;
1097 0 : amp[ncl]/=gainGG; // normalized gas gain
1098 0 : amp[ncl]/=gainPad; //
1099 0 : amp[ncl]/=corrShape;
1100 0 : amp[ncl]/=corrPadType;
1101 0 : amp[ncl]/=corrPos;
1102 0 : amp[ncl]/=corrNorm;
1103 : //
1104 0 : ncl++;
1105 0 : }
1106 :
1107 0 : if (type>3) return ncl;
1108 0 : TMath::Sort(ncl,amp, indexes, kFALSE);
1109 :
1110 0 : if (ncl<10) return 0;
1111 :
1112 : Float_t suma=0;
1113 : Float_t suma2=0;
1114 : Float_t sumn=0;
1115 0 : Int_t icl0=TMath::Nint(ncl*low);
1116 0 : Int_t icl1=TMath::Nint(ncl*up);
1117 0 : for (Int_t icl=icl0; icl<icl1;icl++){
1118 0 : suma+=amp[indexes[icl]];
1119 0 : suma2+=amp[indexes[icl]]*amp[indexes[icl]];
1120 0 : sumn++;
1121 : }
1122 0 : Float_t mean =suma/sumn;
1123 0 : Float_t rms =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1124 : //
1125 : // do time-dependent correction for pressure and temperature variations
1126 : UInt_t runNumber = 1;
1127 : Float_t corrTimeGain = 1;
1128 0 : AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1129 0 : const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1130 0 : if (trans && recoParam->GetUseGainCorrectionTime()>0) {
1131 0 : runNumber = trans->GetCurrentRunNumber();
1132 : //AliTPCcalibDB::Instance()->SetRun(runNumber);
1133 0 : TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1134 0 : if (timeGainSplines) {
1135 0 : UInt_t time = trans->GetCurrentTimeStamp();
1136 0 : AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1137 0 : AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1138 0 : if (fitMIP) {
1139 0 : corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
1140 0 : } else {
1141 0 : if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/
1142 : }
1143 0 : }
1144 0 : }
1145 0 : mean /= corrTimeGain;
1146 0 : rms /= corrTimeGain;
1147 : //
1148 0 : if (returnVal==1) return rms;
1149 0 : if (returnVal==2) return ncl;
1150 0 : return mean;
1151 0 : }
1152 :
1153 : Float_t AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres, Int_t mode, TVectorT<float> *returnVec){
1154 :
1155 : //
1156 : // calculates dedx using the cluster
1157 : // low - up specify trunc mean range - default form 0-0.7
1158 : // type - 1 - max charge or 0- total charge in cluster
1159 : // //2- max no corr 3- total+ correction
1160 : // i1-i2 - the pad-row range used for calculation
1161 : //
1162 : // posNorm - usage of pos normalization
1163 : // returnVal - 0 return mean
1164 : // - 1 return RMS
1165 : // - 2 return number of clusters
1166 : // - 3 ratio
1167 : // - 4 mean upper half
1168 : // - 5 mean - lower half
1169 : // - 6 third moment
1170 : // mode - 0 - linear
1171 : // - 1 - logatithmic
1172 : // rowThres - number of rows before and after given pad row to check for clusters below threshold
1173 : //
1174 : // normalization parametrization taken from AliTPCClusterParam
1175 : //
1176 9720 : if (returnVec) returnVec->ResizeTo(10);
1177 :
1178 4860 : AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
1179 4860 : AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
1180 4860 : AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1181 4860 : const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1182 4860 : if (!parcl) return 0;
1183 4860 : if (!param) return 0;
1184 4860 : Int_t row0 = param->GetNRowLow();
1185 4860 : Int_t row1 = row0+param->GetNRowUp1();
1186 :
1187 4860 : Float_t amp[kMaxRow];
1188 4860 : Int_t indexes[kMaxRow];
1189 : Int_t ncl=0;
1190 : Int_t nclBelowThr = 0; // counts number of clusters below threshold
1191 : //
1192 : //
1193 4860 : Float_t gainGG = 1; // gas gain factor -always enabled
1194 4860 : Float_t gainPad = 1; // gain map - used always
1195 4860 : Float_t corrPos = 1; // local position correction - if posNorm enabled
1196 : //
1197 : //
1198 : //
1199 4860 : if (AliTPCcalibDB::Instance()->GetParameters()){
1200 4860 : gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain
1201 4860 : gainGG *= AliTPCcalibDB::Instance()->GetParameters()->GetNtot()/36.82;//correction for the ionisation
1202 4860 : }
1203 : Double_t timeCut=0;
1204 4860 : if (AliTPCcalibDB::Instance()->IsTrgL0()){
1205 : // by defualt we assume L1 trigger is used - make a correction in case of L0
1206 0 : AliCTPTimeParams* ctp = AliTPCcalibDB::Instance()->GetCTPTimeParams();
1207 0 : Double_t delay = ctp->GetDelayL1L0()*0.000000025;
1208 0 : delay/=param->GetTSample();
1209 : timeCut=delay;
1210 0 : }
1211 4860 : timeCut += recoParam->GetSkipTimeBins();
1212 :
1213 : //
1214 : // extract time-dependent correction for pressure and temperature variations
1215 : //
1216 4860 : UInt_t runNumber = 1;
1217 : Float_t corrTimeGain = 1;
1218 : TObjArray * timeGainSplines = 0x0;
1219 : TGraphErrors * grPadEqual = 0x0;
1220 4860 : TGraphErrors* grChamberGain[4]={0x0,0x0,0x0,0x0};
1221 4860 : TF1* funDipAngle[4]={0x0,0x0,0x0,0x0};
1222 : //
1223 : //
1224 9720 : if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;
1225 4860 : UInt_t time = 1;//
1226 : //
1227 4860 : if (trans) {
1228 4860 : runNumber = trans->GetCurrentRunNumber();
1229 4860 : time = trans->GetCurrentTimeStamp();
1230 : //AliTPCcalibDB::Instance()->SetRun(runNumber);
1231 4860 : timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1232 4860 : if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
1233 0 : AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1234 0 : AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1235 0 : if (fitMIP) {
1236 0 : corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/
1237 0 : } else {
1238 0 : if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */
1239 : }
1240 : //
1241 0 : if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
1242 0 : if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
1243 0 : const char* names[4]={"SHORT","MEDIUM","LONG","ABSOLUTE"};
1244 0 : for (Int_t iPadRegion=0; iPadRegion<4; ++iPadRegion) {
1245 0 : grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion]));
1246 0 : if (type==1) funDipAngle[iPadRegion]=(TF1*)timeGainSplines->FindObject(Form("TF1_QMAX_DIPANGLE_%s_BEAM_ALL",names[iPadRegion]));
1247 0 : if (type==0) funDipAngle[iPadRegion]=(TF1*)timeGainSplines->FindObject(Form("TF1_QTOT_DIPANGLE_%s_BEAM_ALL",names[iPadRegion]));
1248 : }
1249 0 : }
1250 : }
1251 :
1252 : const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1253 4860 : const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1254 : const Float_t kedgey =3.;
1255 : //
1256 : //
1257 732240 : for (Int_t irow=i1; irow<i2; irow++){
1258 361260 : AliTPCclusterMI* cluster = GetClusterPointer(irow);
1259 361260 : if (!cluster && irow > 1 && irow < 157) {
1260 : Bool_t isClBefore = kFALSE;
1261 : Bool_t isClAfter = kFALSE;
1262 425072 : for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
1263 0 : AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1264 0 : if (clusterBefore) isClBefore = kTRUE;
1265 0 : AliTPCclusterMI * clusterAfter = GetClusterPointer(irow + ithres);
1266 0 : if (clusterAfter) isClAfter = kTRUE;
1267 : }
1268 212536 : if (isClBefore && isClAfter) nclBelowThr++;
1269 212536 : }
1270 579080 : if (!cluster) continue;
1271 143440 : if (cluster->GetTimeBin()<timeCut) continue; //reject clusters at the gating grid opening
1272 : //
1273 : //
1274 151972 : if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1275 : //
1276 134908 : const AliTPCTrackerPoints::Point * point = GetTrackPoint(irow);
1277 134908 : if (point==0) continue;
1278 134908 : Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1279 136760 : if (rsigmay > kClusterShapeCut) continue;
1280 : //
1281 133464 : if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1282 : //
1283 397944 : Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1284 132648 : Int_t ipad= 0;
1285 132648 : if (irow>=row0) ipad=1;
1286 132648 : if (irow>=row1) ipad=2;
1287 : //
1288 : //
1289 : //
1290 132648 : AliTPCCalPad * gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1291 132648 : if (gainMap) {
1292 : //
1293 : // Get gainPad - pad by pad calibration
1294 : //
1295 : Float_t factor = 1;
1296 132648 : AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1297 132648 : if (irow < row0) { // IROC
1298 36528 : factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1299 36528 : } else { // OROC
1300 96120 : factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1301 : }
1302 265296 : if (factor>0.3) gainPad=factor;
1303 132648 : }
1304 : //
1305 : // Do position normalization - relative distance to
1306 : // center of pad- time bin
1307 :
1308 132648 : Float_t ty = TMath::Abs(point->GetAngleY());
1309 132648 : Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1310 132648 : Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1311 132648 : Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1312 :
1313 132648 : yres0 *=parcl->GetQnormCorr(ipad, type,0);
1314 132648 : zres0 *=parcl->GetQnormCorr(ipad, type,1);
1315 132648 : Float_t effLength=parcl->GetQnormCorr(ipad, type,4)*0.5;
1316 132648 : Float_t effDiff =(parcl->GetQnormCorr(ipad, type,2)+parcl->GetQnormCorr(ipad, type,3))*0.5;
1317 132648 : Float_t corrThr=0;
1318 132648 : Float_t corrThrMax=0;
1319 : //
1320 132648 : if (type==1) {
1321 101248 : corrThr = parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1322 50624 : cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff);
1323 50624 : corrPos= parcl->GetQnormCorr(ipad, type,5)*corrThr;
1324 50624 : Float_t drm = 0.5-TMath::Abs(cluster->GetZ()/250.);
1325 50624 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1326 50624 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1327 50624 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1328 : //
1329 50624 : }
1330 132648 : if (type==0) {
1331 164048 : corrThr = parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1332 82024 : cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff);
1333 82024 : corrPos=parcl->GetQnormCorr(ipad, type,5)*corrThr;
1334 82024 : Float_t drm = 0.5-TMath::Abs(cluster->GetZ()/250.);
1335 82024 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1336 82024 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1337 82024 : corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1338 : //
1339 82024 : }
1340 : //
1341 : // pad region equalization outside of cluster param
1342 : //
1343 132648 : Float_t gainEqualPadRegion = 1;
1344 132648 : if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad);
1345 : //
1346 : // chamber-by-chamber equalization outside gain map
1347 : //
1348 132648 : Float_t gainChamber = 1;
1349 132648 : if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) {
1350 0 : gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector());
1351 0 : if (gainChamber==0) gainChamber=1; // in case old calibation was used before use no correction
1352 : }
1353 : //
1354 : // dip angle correction
1355 : //
1356 132648 : Float_t corrDipAngle = 1;
1357 132648 : Float_t corrDipAngleAbs = 1;
1358 : // if (grDipAngle[ipad]) corrDipAngle = grDipAngle[ipad]->Eval(GetTgl());
1359 132648 : Double_t tgl=GetTgl();
1360 132648 : if (funDipAngle[ipad]) corrDipAngle = funDipAngle[ipad]->Eval(tgl);
1361 132648 : if (funDipAngle[3]) corrDipAngleAbs = funDipAngle[3]->Eval(tgl);
1362 : //
1363 : // pressure temperature and high voltage correction
1364 : //
1365 132648 : Double_t correctionHVandPT = AliTPCcalibDB::Instance()->GetGainCorrectionHVandPT(time, runNumber,cluster->GetDetector(), 5 , recoParam->GetGainCorrectionHVandPTMode());
1366 : //
1367 132648 : if ((AliTPCReconstructor::StreamLevel()&AliTPCtracker::kStreamSeeddEdx)){ // this part of the code is for the test purposes only
1368 0 : TTreeSRedirector *pcstream = AliTPCReconstructor:: GetDebugStreamer();
1369 0 : TVectorF vecDEDX(9,fDEDX);
1370 0 : TVectorF vecSDEDX(4,fSDEDX);
1371 0 : TVectorF vecRDEDX(4);
1372 0 : corrThrMax = parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1373 0 : cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff);
1374 :
1375 0 : for (Int_t i=0; i<4; i++) vecRDEDX[i]=(fNCDEDXInclThres[i]>0)?Float_t(fNCDEDX[i])/Float_t(fNCDEDXInclThres[i]):0;
1376 0 : if (pcstream){
1377 0 : (*pcstream)<<"dEdxCorrDump"<< // streamer to check dEdx correction calibration
1378 : //
1379 0 : "cl.="<<cluster<<
1380 0 : "ipad="<<ipad<<
1381 0 : "time="<<time<<
1382 0 : "runNumber="<<runNumber<<
1383 0 : "vecDEDX.="<<&vecDEDX<<
1384 0 : "vecSDEDX.="<<&vecSDEDX<<
1385 0 : "vecRDEDX.="<<&vecRDEDX<<
1386 0 : "ty="<<ty<<
1387 0 : "tz="<<tz<<
1388 0 : "yres0="<<yres0<<
1389 0 : "zres0="<<zres0<<
1390 0 : "qpt="<<fP[4]<<
1391 0 : "type="<<type<<
1392 0 : "effLength="<<effLength<<
1393 0 : "effDiff="<<effDiff<<
1394 0 : "corrThr="<<corrThr<<
1395 0 : "corrThrMax="<<corrThrMax<<
1396 0 : "gainGG="<<gainGG<<
1397 0 : "correctionHVandPT="<<correctionHVandPT<<
1398 0 : "gainPad="<<gainPad<<
1399 0 : "corrPos="<<corrPos<<
1400 0 : "gainEqualPadRegion="<<gainEqualPadRegion<<
1401 0 : "gainChamber="<<gainChamber<<
1402 0 : "corrDipAngle="<<corrDipAngle<<
1403 0 : "corrDipAngleAbs="<<corrDipAngleAbs<<
1404 : "\n";
1405 : }
1406 0 : }
1407 132648 : amp[ncl]=charge;
1408 132648 : amp[ncl]/=gainGG; // nominal gas gain
1409 132648 : amp[ncl]/=correctionHVandPT; // correction for the HV and P/T - time dependent
1410 132648 : amp[ncl]/=gainPad; //
1411 132648 : amp[ncl]/=corrPos;
1412 132648 : amp[ncl]/=gainEqualPadRegion;
1413 132648 : amp[ncl]/=gainChamber;
1414 132648 : amp[ncl]/=corrDipAngle;
1415 132648 : amp[ncl]/=corrDipAngleAbs;
1416 : //
1417 132648 : ncl++;
1418 132648 : }
1419 :
1420 4860 : if (type==2) return ncl;
1421 4860 : TMath::Sort(ncl,amp, indexes, kFALSE);
1422 : //
1423 7490 : if (ncl<10) return 0;
1424 : //
1425 2230 : Double_t ampWithBelow[ncl + nclBelowThr];
1426 268860 : for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
1427 132200 : if (iCl < nclBelowThr) {
1428 0 : ampWithBelow[iCl] = amp[indexes[0]];
1429 0 : } else {
1430 132200 : ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
1431 : }
1432 : }
1433 : //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
1434 : //
1435 : Float_t suma=0;
1436 : Float_t suma2=0;
1437 : Float_t suma3=0;
1438 : Float_t sumaS=0;
1439 : Float_t sumn=0;
1440 : // upper,and lower part statistic
1441 : Float_t sumL=0, sumL2=0, sumLN=0;
1442 : Float_t sumD=0, sumD2=0, sumDN=0;
1443 :
1444 2230 : Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
1445 2230 : Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
1446 2230 : Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
1447 : //
1448 157648 : for (Int_t icl=icl0; icl<icl1;icl++){
1449 76594 : if (ampWithBelow[icl]<0.1) continue;
1450 76594 : Double_t camp=ampWithBelow[icl]/corrTimeGain;
1451 76594 : if (mode==1) camp= TMath::Log(camp);
1452 76594 : if (icl<icl1){
1453 76594 : suma+=camp;
1454 76594 : suma2+=camp*camp;
1455 76594 : suma3+=camp*camp*camp;
1456 76594 : sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
1457 76594 : sumn++;
1458 76594 : }
1459 76594 : if (icl>iclm){
1460 33410 : sumL+=camp;
1461 33410 : sumL2+=camp*camp;
1462 33410 : sumLN++;
1463 33410 : }
1464 76594 : if (icl<=iclm){
1465 43184 : sumD+=camp;
1466 43184 : sumD2+=camp*camp;
1467 43184 : sumDN++;
1468 43184 : }
1469 76594 : }
1470 : //
1471 : Float_t mean = 0;
1472 : Float_t meanL = 0;
1473 : Float_t meanD = 0; // lower half mean
1474 4460 : if (sumn > 1e-30) mean =suma/sumn;
1475 4460 : if (sumLN > 1e-30) meanL =sumL/sumLN;
1476 4460 : if (sumDN > 1e-30) meanD =(sumD/sumDN);
1477 : /*
1478 : Float_t mean =suma/sumn;
1479 : Float_t meanL = sumL/sumLN;
1480 : Float_t meanD =(sumD/sumDN); // lower half mean
1481 : */
1482 :
1483 : Float_t rms = 0;
1484 : Float_t mean2=0;
1485 : Float_t mean3=0;
1486 : Float_t meanS=0;
1487 :
1488 2230 : if(sumn>0){
1489 2230 : rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1490 : mean2=suma2/sumn;
1491 2230 : mean3=suma3/sumn;
1492 2230 : meanS=sumaS/sumn;
1493 2230 : }
1494 :
1495 4460 : if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
1496 4460 : if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
1497 4460 : if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
1498 : //
1499 2230 : if (mode==1) mean=TMath::Exp(mean);
1500 2230 : if (mode==1) meanL=TMath::Exp(meanL); // upper truncation
1501 2230 : if (mode==1) meanD=TMath::Exp(meanD); // lower truncation
1502 : //
1503 : //delete [] ampWithBelow; //return? // RS made on stack
1504 :
1505 :
1506 : //
1507 2230 : if(returnVec){
1508 2230 : (*returnVec)(0) = mean;
1509 2230 : (*returnVec)(1) = rms;
1510 2230 : (*returnVec)(2) = ncl;
1511 2230 : (*returnVec)(3) = Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1512 2230 : (*returnVec)(4) = meanL;
1513 2230 : (*returnVec)(5) = meanD;
1514 2230 : (*returnVec)(6) = mean2;
1515 2230 : (*returnVec)(7) = mean3;
1516 2230 : (*returnVec)(8) = meanS;
1517 2230 : (*returnVec)(9) = nclBelowThr;
1518 2230 : }
1519 :
1520 2230 : if (returnVal==1) return rms;
1521 2230 : if (returnVal==2) return ncl;
1522 2230 : if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1523 2230 : if (returnVal==4) return meanL;
1524 2230 : if (returnVal==5) return meanD;
1525 2230 : if (returnVal==6) return mean2;
1526 2230 : if (returnVal==7) return mean3;
1527 2230 : if (returnVal==8) return meanS;
1528 2230 : if (returnVal==9) return nclBelowThr;
1529 2230 : return mean;
1530 11950 : }
1531 :
1532 :
1533 :
1534 :
1535 : Float_t AliTPCseed::CookShape(Int_t type){
1536 : //
1537 : //
1538 : //
1539 : //-----------------------------------------------------------------
1540 : // This funtion calculates dE/dX within the "low" and "up" cuts.
1541 : //-----------------------------------------------------------------
1542 : Float_t means=0;
1543 : Float_t meanc=0;
1544 0 : for (Int_t i =0; i<kMaxRow;i++) {
1545 : //
1546 : //AliTPCclusterMI * cl = GetClusterPointer(i);
1547 : //if (cl==0) continue;
1548 0 : if (GetClusterIndex2(i)<0) continue;
1549 : //
1550 0 : const AliTPCTrackerPoints::Point * point = GetTrackPoint(i);
1551 0 : if (point==0) continue;
1552 0 : Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1553 0 : Float_t rsigmaz = TMath::Sqrt(point->GetSigmaZ());
1554 0 : Float_t rsigma = (rsigmay+rsigmaz)*0.5;
1555 0 : if (type==0) means+=rsigma;
1556 0 : if (type==1) means+=rsigmay;
1557 0 : if (type==2) means+=rsigmaz;
1558 0 : meanc++;
1559 0 : }
1560 0 : Float_t mean = (meanc>0)? means/meanc:0;
1561 0 : return mean;
1562 : }
1563 :
1564 :
1565 :
1566 : Int_t AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, AliExternalTrackParam * parout){
1567 : //
1568 : // Refit the track
1569 : // return value - number of used clusters
1570 : //
1571 : //
1572 : const Int_t kMinNcl =10;
1573 10 : static AliTPCseed trackTmp;
1574 : AliTPCseed* track = &trackTmp;
1575 2 : track->AliExternalTrackParam::operator=(*seed);
1576 : Int_t sector=-1;
1577 : // reset covariance
1578 : //
1579 2 : Double_t covar[15];
1580 64 : for (Int_t i=0;i<15;i++) covar[i]=0;
1581 2 : covar[0]=10.*10.;
1582 2 : covar[2]=10.*10.;
1583 2 : covar[5]=10.*10./(64.*64.);
1584 2 : covar[9]=10.*10./(64.*64.);
1585 2 : covar[14]=1*1;
1586 : //
1587 :
1588 : Float_t xmin=1000, xmax=-10000;
1589 : Int_t imin=158, imax=0;
1590 640 : for (Int_t i=0;i<kMaxRow;i++) {
1591 318 : AliTPCclusterMI *c=seed->GetClusterPointer(i);
1592 656 : if (!c || (seed->GetClusterIndex(i) & 0x8000)) continue;
1593 254 : if (sector<0) sector = c->GetDetector();
1594 254 : if (c->GetX()<xmin) xmin=c->GetX();
1595 504 : if (c->GetX()>xmax) xmax=c->GetX();
1596 254 : if (i<imin) imin=i;
1597 502 : if (i>imax) imax=i;
1598 252 : }
1599 2 : if(imax-imin<kMinNcl) {
1600 : //RS delete track;
1601 0 : return 0 ;
1602 : }
1603 : // Not succes to rotate
1604 2 : if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1605 : // delete track;
1606 0 : return 0;
1607 : }
1608 : //
1609 : //
1610 : // fit from inner to outer row
1611 : //
1612 2 : AliExternalTrackParam paramIn;
1613 2 : AliExternalTrackParam paramOut;
1614 : Bool_t isOK=kTRUE;
1615 : Int_t ncl=0;
1616 : //
1617 : //
1618 : //
1619 632 : for (Int_t i=imin; i<=imax; i++){
1620 314 : AliTPCclusterMI *c=seed->GetClusterPointer(i);
1621 920 : if (!c || (seed->GetClusterIndex(i) & 0x8000)) continue;
1622 : // if (RejectCluster(c,track)) continue;
1623 504 : sector = (c->GetDetector()%18);
1624 504 : if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1625 : //continue;
1626 : }
1627 252 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1628 252 : Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1629 504 : if (!track->PropagateTo(r[0])) {
1630 : isOK=kFALSE;
1631 0 : }
1632 504 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1633 252 : }
1634 : //RS if (!isOK) { delete track; return 0;}
1635 2 : if (!isOK) { return 0;}
1636 2 : track->AddCovariance(covar);
1637 : //
1638 : //
1639 : //
1640 632 : for (Int_t i=imax; i>=imin; i--){
1641 314 : AliTPCclusterMI *c=seed->GetClusterPointer(i);
1642 920 : if (!c || (seed->GetClusterIndex(i) & 0x8000)) continue;
1643 : //if (RejectCluster(c,track)) continue;
1644 504 : sector = (c->GetDetector()%18);
1645 504 : if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1646 : //continue;
1647 : }
1648 252 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1649 252 : Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1650 504 : if (!track->PropagateTo(r[0])) {
1651 : isOK=kFALSE;
1652 0 : }
1653 504 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1654 252 : }
1655 : //if (!isOK) { delete track; return 0;}
1656 2 : paramIn = *track;
1657 2 : track->AddCovariance(covar);
1658 : //
1659 : //
1660 632 : for (Int_t i=imin; i<=imax; i++){
1661 314 : AliTPCclusterMI *c=seed->GetClusterPointer(i);
1662 920 : if (!c || (seed->GetClusterIndex(i) & 0x8000)) continue;
1663 504 : sector = (c->GetDetector()%18);
1664 504 : if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1665 : //continue;
1666 : }
1667 252 : ncl++;
1668 : //if (RejectCluster(c,track)) continue;
1669 252 : Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1670 252 : Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1671 504 : if (!track->PropagateTo(r[0])) {
1672 : isOK=kFALSE;
1673 0 : }
1674 504 : if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1675 252 : }
1676 : //if (!isOK) { delete track; return 0;}
1677 2 : paramOut=*track;
1678 : //
1679 : //
1680 : //
1681 4 : if (parin) (*parin)=paramIn;
1682 4 : if (parout) (*parout)=paramOut;
1683 : //RS delete track;
1684 2 : return ncl;
1685 4 : }
1686 :
1687 :
1688 :
1689 : Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){
1690 : //
1691 : //
1692 : //
1693 0 : return kFALSE;
1694 : }
1695 :
1696 :
1697 :
1698 :
1699 :
1700 :
1701 : void AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1702 : Double_t& erry, Double_t &errz)
1703 : {
1704 : //
1705 : // Get cluster error at given position
1706 : //
1707 4792 : AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1708 : Double_t tany,tanz;
1709 2396 : Double_t snp1=param->GetSnp();
1710 2396 : tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1711 : //
1712 2396 : Double_t tgl1=param->GetTgl();
1713 2396 : tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1714 : //
1715 : Int_t padSize = 0; // short pads
1716 2396 : if (cluster->GetDetector() >= 36) {
1717 : padSize = 1; // medium pads
1718 1996 : if (cluster->GetRow() > 63) padSize = 2; // long pads
1719 : }
1720 :
1721 2396 : erry = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) );
1722 2396 : errz = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) );
1723 2396 : }
1724 :
1725 :
1726 : void AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1727 : Double_t& rmsy, Double_t &rmsz)
1728 : {
1729 : //
1730 : // Get cluster error at given position
1731 : //
1732 0 : AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1733 : Double_t tany,tanz;
1734 0 : Double_t snp1=param->GetSnp();
1735 0 : tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1736 : //
1737 0 : Double_t tgl1=param->GetTgl();
1738 0 : tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1739 : //
1740 : Int_t padSize = 0; // short pads
1741 0 : if (cluster->GetDetector() >= 36) {
1742 : padSize = 1; // medium pads
1743 0 : if (cluster->GetRow() > 63) padSize = 2; // long pads
1744 : }
1745 :
1746 0 : rmsy = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) );
1747 0 : rmsz = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax()));
1748 0 : }
1749 :
1750 :
1751 :
1752 : Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){
1753 : //Geoetrical
1754 : //ty - tangent in local y direction
1755 : //tz - tangent
1756 : //
1757 0 : Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz);
1758 0 : return norm;
1759 : }
1760 :
1761 : Double_t AliTPCseed::GetQCorrShape(Int_t ipad, Int_t type,Float_t z, Float_t ty, Float_t tz, Float_t /*q*/, Float_t /*thr*/){
1762 : //
1763 : // Q normalization
1764 : //
1765 : // return value = Q Normalization factor
1766 : // Normalization - 1 - shape factor part for full drift
1767 : // 1 - electron attachment for 0 drift
1768 :
1769 : // Input parameters:
1770 : //
1771 : // ipad - 0 short pad
1772 : // 1 medium pad
1773 : // 2 long pad
1774 : //
1775 : // type - 0 qmax
1776 : // - 1 qtot
1777 : //
1778 : //z - z position (-250,250 cm)
1779 : //ty - tangent in local y direction
1780 : //tz - tangent
1781 : //
1782 :
1783 0 : AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
1784 0 : AliTPCParam * paramTPC = AliTPCcalibDB::Instance()->GetParameters();
1785 :
1786 0 : if (!paramCl) return 1;
1787 : //
1788 0 : Double_t dr = 250.-TMath::Abs(z);
1789 0 : Double_t sy = paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty));
1790 0 : Double_t sy0= paramCl->GetRMS0(0,ipad, 250, 0);
1791 0 : Double_t sz = paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz));
1792 0 : Double_t sz0= paramCl->GetRMS0(1,ipad, 250, 0);
1793 :
1794 0 : Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz));
1795 :
1796 :
1797 0 : Double_t dt = 1000000*(dr/paramTPC->GetDriftV()); //time in microsecond
1798 0 : Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt);
1799 : //
1800 : //
1801 0 : if (type==0) return sfactorMax*attProb;
1802 :
1803 0 : return attProb;
1804 :
1805 :
1806 0 : }
1807 :
1808 :
1809 : /*
1810 : //_______________________________________________________________________
1811 : Float_t AliTPCseed::GetTPCClustInfo(Int_t nNeighbours, Int_t type, Int_t row0, Int_t row1, TVectorT<float> *returnVec)
1812 : {
1813 : //
1814 : // TPC cluster information
1815 : // type 0: get fraction of found/findable clusters with neighbourhood definition
1816 : // 1: found clusters
1817 : // 2: findable (number of clusters above and below threshold)
1818 : //
1819 : // definition of findable clusters:
1820 : // a cluster is defined as findable if there is another cluster
1821 : // within +- nNeighbours pad rows. The idea is to overcome threshold
1822 : // effects with a very simple algorithm.
1823 : //
1824 :
1825 : const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1826 : const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1827 : const Float_t kedgey =3.;
1828 :
1829 : Float_t ncl = 0;
1830 : Float_t nclBelowThr = 0; // counts number of clusters below threshold
1831 :
1832 : for (Int_t irow=row0; irow<row1; irow++){
1833 : AliTPCclusterMI* cluster = GetClusterPointer(irow);
1834 :
1835 : if (!cluster && irow > 1 && irow < 157) {
1836 : Bool_t isClBefore = kFALSE;
1837 : Bool_t isClAfter = kFALSE;
1838 : for(Int_t ithres = 1; ithres <= nNeighbours; ithres++) {
1839 : AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1840 : if (clusterBefore) isClBefore = kTRUE;
1841 : AliTPCclusterMI * clusterAfter = GetClusterPointer(irow + ithres);
1842 : if (clusterAfter) isClAfter = kTRUE;
1843 : }
1844 : if (isClBefore && isClAfter) nclBelowThr++;
1845 : }
1846 : if (!cluster) continue;
1847 : //
1848 : //
1849 : if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1850 : //
1851 : const AliTPCTrackerPoints::Point * point = GetTrackPoint(irow);
1852 : if (point==0) continue;
1853 : Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1854 : if (rsigmay > kClusterShapeCut) continue;
1855 : //
1856 : if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1857 : ncl++;
1858 : }
1859 : if(returnVec->GetNoElements != 3){
1860 : returnVec->ResizeTo(3);
1861 : }
1862 : Float_t nclAll = nclBelowThr+ncl;
1863 : returnVec(0) = nclAll>0?ncl/nclAll:0;
1864 : returnVec(1) = ncl;
1865 : returnVec(2) = nclAll;
1866 :
1867 : if(ncl<10)
1868 : return 0;
1869 : if(type==0)
1870 : if(nclAll>0)
1871 : return ncl/nclAll;
1872 : if(type==1)
1873 : return ncl;
1874 : if(type==2)
1875 : return nclAll;
1876 : return 0;
1877 : }
1878 : */
1879 : //_______________________________________________________________________
1880 : Int_t AliTPCseed::GetNumberOfClustersIndices() {
1881 : Int_t ncls = 0;
1882 0 : for (int i=0; i < kMaxRow; i++) {
1883 0 : if ((fIndex[i] & 0x8000) == 0)
1884 0 : ncls++;
1885 : }
1886 0 : return ncls;
1887 : }
1888 :
1889 : //_______________________________________________________________________
1890 : void AliTPCseed::Clear(Option_t*)
1891 : {
1892 : // formally seed may allocate memory for clusters (althought this should not happen for
1893 : // the seeds in the pool). Hence we need this method for fwd. compatibility
1894 0 : if (fClusterPointer) {
1895 0 : if (fClusterOwner) for (int i=kMaxRow;i--;) {delete fClusterPointer[i]; fClusterPointer[i] = 0;}
1896 0 : delete[] fClusterPointer;
1897 0 : fNClStore = 0;
1898 0 : }
1899 0 : fTrackPointsArr.Clear();
1900 0 : ResetBit(0xffffffff);
1901 0 : }
1902 :
1903 : //_______________________________________________________________________
1904 : void AliTPCseed::SetClusterPointer(Int_t irow, AliTPCclusterMI* cl)
1905 : {
1906 : // if needed, create array and set the pointer
1907 0 : if (!fClusterPointer) {
1908 0 : fClusterPointer = new AliTPCclusterMI*[kMaxRow];
1909 0 : fNClStore = kMaxRow;
1910 0 : memset(fClusterPointer,0,kMaxRow*sizeof(AliTPCclusterMI*));
1911 0 : }
1912 0 : fClusterPointer[irow]=cl;
1913 0 : }
1914 :
1915 :
1916 : TObject* AliTPCseed::Clone(const char* /*newname*/) const
1917 : {
1918 : // temporary override TObject::Clone to avoid crashes in reco
1919 : AliTPCseed* src = (AliTPCseed*)this;
1920 0 : AliTPCseed* dst = new AliTPCseed(*src,fClusterOwner);
1921 0 : return dst;
1922 0 : }
1923 :
1924 : //__________________________________________________
1925 0 : void AliTPCseed::TagSuppressSharedClusters()
1926 : {
1927 : // RS: special function to flag the cluster if not flagged yet or remove its pointer if flagged
1928 : // This is needed for deletion of the TPCseeds stored in the ESDfriendTracks in a safe way, avoiding
1929 : // deletion of shared clusters. DO not use this method except from
1930 : // AliESDfriend->DeleteTracksSafe->TagSuppressSharedObjectsBeforeDeletion(): the clusters are DISABLED
1931 : //
1932 0 : if (!fClusterOwner || !fClusterPointer) return;
1933 : AliTPCclusterMI* cl=0;
1934 0 : for (int i=kMaxRow;i--;) {
1935 0 : if (!(cl=fClusterPointer[i])) continue;
1936 0 : if (cl->IsDisabled()) fClusterPointer[i] = 0;
1937 0 : else cl->Disable();
1938 : }
1939 0 : }
|
