git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Restoring some user friendly functions
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
CommitLineData
1c53abe2 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/*
22 AliTPC parallel tracker -
23 How to use? -
24 run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
25 run AliTPCFindTracksMI.C macro - to find tracks
26 tracks are written to AliTPCtracks.root file
27 for comparison also seeds are written to the same file - to special branch
28*/
29
30//-------------------------------------------------------
31// Implementation of the TPC tracker
32//
33// Origin: Marian Ivanov Marian.Ivanov@cern.ch
34//
35//-------------------------------------------------------
36
37#include <TObjArray.h>
38#include <TFile.h>
39#include <TTree.h>
cc5e9db0 40#include "Riostream.h"
1c53abe2 41
42#include "AliTPCtrackerMI.h"
43#include "AliTPCclusterMI.h"
44#include "AliTPCParam.h"
45#include "AliTPCClustersRow.h"
46#include "AliComplexCluster.h"
1627d1c4 47#include "AliTPCpolyTrack.h"
1c53abe2 48#include "TStopwatch.h"
49
50
51ClassImp(AliTPCseed)
52
53AliTPCclusterTracks::AliTPCclusterTracks(){
54 // class for storing overlaping info
55 fTrackIndex[0]=-1;
56 fTrackIndex[1]=-1;
57 fTrackIndex[2]=-1;
58 fDistance[0]=1000;
59 fDistance[1]=1000;
60 fDistance[2]=1000;
61}
62
63
64
65
66
67Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
68
69 Int_t sec=(i&0xff000000)>>24;
70 Int_t row = (i&0x00ff0000)>>16;
71 track->fRow=(i&0x00ff0000)>>16;
72 track->fSector = sec;
73 // Int_t index = i&0xFFFF;
74 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
75 track->fClusterIndex[track->fRow] = i;
76 track->fFirstPoint = row;
77 //
78
79 AliTPCTrackPoint *trpoint =track->GetTrackPoint(track->fRow);
80 Float_t angle2 = track->GetSnp()*track->GetSnp();
81 angle2 = TMath::Sqrt(angle2/(1-angle2));
82 //
83 //SET NEW Track Point
84 //
85 if (c!=0){
86 //if we have a cluster
87 trpoint->GetCPoint().SetY(c->GetY());
88 trpoint->GetCPoint().SetZ(c->GetZ());
89 //
90 trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
91 trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
92 //
93 trpoint->GetCPoint().SetType(c->GetType());
94 trpoint->GetCPoint().SetQ(c->GetQ());
95 trpoint->GetCPoint().SetMax(c->GetMax());
96 //
97 trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
98 trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
99 //
100 }
101 trpoint->GetTPoint().SetX(track->GetX());
102 trpoint->GetTPoint().SetY(track->GetY());
103 trpoint->GetTPoint().SetZ(track->GetZ());
104 //
105 trpoint->GetTPoint().SetAngleY(angle2);
106 trpoint->GetTPoint().SetAngleZ(track->GetTgl());
107
108
109
110 if (chi2>10){
111 // printf("suspicious chi2 %f\n",chi2);
112 }
113 // if (track->fIsSeeding){
114 track->fErrorY2 *= 1.2;
115 track->fErrorY2 += 0.0064;
116 track->fErrorZ2 *= 1.2;
117 track->fErrorY2 += 0.005;
118
119 //}
120
121 return track->Update(c,chi2,i);
122
123}
124//_____________________________________________________________________________
125AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn):
126AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
127{
128 //---------------------------------------------------------------------
129 // The main TPC tracker constructor
130 //---------------------------------------------------------------------
131 fInnerSec=new AliTPCSector[fkNIS];
132 fOuterSec=new AliTPCSector[fkNOS];
133
134 Int_t i;
135 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
136 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
137
138 fN=0; fSectors=0;
139
140 fClustersArray.Setup(par);
141 fClustersArray.SetClusterType("AliTPCclusterMI");
142
143 char cname[100];
144 if (eventn==-1) {
145 sprintf(cname,"TreeC_TPC");
146 }
147 else {
148 sprintf(cname,"TreeC_TPC_%d",eventn);
149 }
150
151 fClustersArray.ConnectTree(cname);
152
153 fEventN = eventn;
154 fSeeds=0;
155 fNtracks = 0;
156 fParam = par;
157}
158
159//_____________________________________________________________________________
160AliTPCtrackerMI::~AliTPCtrackerMI() {
161 //------------------------------------------------------------------
162 // TPC tracker destructor
163 //------------------------------------------------------------------
164 delete[] fInnerSec;
165 delete[] fOuterSec;
166 if (fSeeds) {
167 fSeeds->Delete();
168 delete fSeeds;
169 }
170}
171
172
173Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
174 //
175 //
176 Float_t snoise2;
177 Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
178
179 //cluster "quality"
180 Float_t rsigmay = 1;
181 Int_t ctype = 0;
182
183 //standard if we don't have cluster - take MIP
184 const Float_t chmip = 50.;
185 Float_t amp = chmip/0.3;
186 Float_t nel;
187 Float_t nprim;
188 if (cl){
189 amp = cl->GetQ();
190 rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
191 ctype = cl->GetType();
192 }
193
194
195 Float_t landau=2 ; //landau fluctuation part
196 Float_t gg=2; // gg fluctuation part
197 Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
198
199
200 if (fSectors==fInnerSec){
201 snoise2 = 0.0004;
202 nel = 0.268*amp;
203 nprim = 0.155*amp;
204 gg = (2+0.0002*amp)/nel;
205 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
206 if (landau>1) landau=1;
207 }
208 else {
209 snoise2 = 0.0004;
210 nel = 0.3*amp;
211 nprim = 0.133*amp;
212 gg = (2+0.0002*amp)/nel;
213 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
214 if (landau>1) landau=1;
215 }
216
217
218 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
219 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
220 angle2 = angle2/(1-angle2);
221 Float_t angular = landau*angle2*padlength*padlength/12.;
222 Float_t res = sdiff + angular;
223
224
225 if ((ctype==0) && (fSectors ==fOuterSec))
226 res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
227
228 if ((ctype==0) && (fSectors ==fInnerSec))
229 res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
230
231
232 if ((ctype>0))
233 res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
234
235 if (ctype<0)
236 res*=2.4; // overestimate error 2 times
237
238 res+= snoise2;
239
240 if (res<2*snoise2)
241 res = 2*snoise2;
242
243 seed->SetErrorY2(res);
244 return res;
245}
246
247
248
249
250Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
251 //
252 //
253 Float_t snoise2;
254 Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
255 //signal quality
256 Float_t rsigmaz=1;
257 Int_t ctype =0;
258
259 const Float_t chmip = 50.;
260 Float_t amp = chmip/0.3;
261 Float_t nel;
262 Float_t nprim;
263 if (cl){
264 amp = cl->GetQ();
265 rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
266 ctype = cl->GetType();
267 }
268
269 //
270 Float_t landau=2 ; //landau fluctuation part
271 Float_t gg=2; // gg fluctuation part
272 Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
273
274 if (fSectors==fInnerSec){
275 snoise2 = 0.0004;
276 nel = 0.268*amp;
277 nprim = 0.155*amp;
278 gg = (2+0.0002*amp)/nel;
279 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
280 if (landau>1) landau=1;
281 }
282 else {
283 snoise2 = 0.0004;
284 nel = 0.3*amp;
285 nprim = 0.133*amp;
286 gg = (2+0.0002*amp)/nel;
287 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
288 if (landau>1) landau=1;
289 }
290 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
291
292 Float_t angle = seed->GetTgl();
293 Float_t angular = landau*angle*angle*padlength*padlength/12.;
294 Float_t res = sdiff + angular;
295
296 if ((ctype==0) && (fSectors ==fOuterSec))
297 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
298
299 if ((ctype==0) && (fSectors ==fInnerSec))
300 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
301 if ((ctype>0))
302 res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064; //0.31+0.147*ctype;
303 if (ctype<0)
304 res*=1.3;
305 if ((ctype<0) &&amp<70)
306 res*=1.3;
307
308 res += snoise2;
309 if (res<2*snoise2)
310 res = 2*snoise2;
311
312 seed->SetErrorZ2(res);
313 return res;
314}
315
316
1627d1c4 317void AliTPCseed::Reset()
318{
319 //
320 //PH SetN(0);
321 fNFoundable = 0;
322 ResetCovariance();
323 SetChi2(0);
324 for (Int_t i=0;i<200;i++) fClusterIndex[i]=-1;
325}
326
1c53abe2 327
328Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
329{
330 //-----------------------------------------------------------------
331 // This function find proloncation of a track to a reference plane x=xk.
332 // doesn't change internal state of the track
333 //-----------------------------------------------------------------
334
335 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
336 // Double_t y1=fP0, z1=fP1;
337 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
338 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
339
340 y = fP0;
341 z = fP1;
342 y += dx*(c1+c2)/(r1+r2);
343 z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
344 return 0;
345}
346
347
348//_____________________________________________________________________________
349Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
350{
351 //-----------------------------------------------------------------
352 // This function calculates a predicted chi2 increment.
353 //-----------------------------------------------------------------
354 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
355 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
356 r00+=fC00; r01+=fC10; r11+=fC11;
357
358 Double_t det=r00*r11 - r01*r01;
359 if (TMath::Abs(det) < 1.e-10) {
360 Int_t n=GetNumberOfClusters();
361 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
362 return 1e10;
363 }
364 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
365
366 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
367
368 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
369}
370
371
372//_________________________________________________________________________________________
373
374
375Int_t AliTPCseed::Compare(const TObject *o) const {
376 //-----------------------------------------------------------------
377 // This function compares tracks according to the sector - for given sector according z
378 //-----------------------------------------------------------------
379 AliTPCseed *t=(AliTPCseed*)o;
380 if (t->fSector>fSector) return -1;
381 if (t->fSector<fSector) return 1;
382
383 Double_t z2 = t->GetZ();
384 Double_t z1 = GetZ();
385 if (z2>z1) return 1;
386 if (z2<z1) return -1;
387 return 0;
388}
389
390
391
392
393
394
395//_____________________________________________________________________________
396Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
397 //-----------------------------------------------------------------
398 // This function associates a cluster with this track.
399 //-----------------------------------------------------------------
400 // Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
401 //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
402 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
403
404 r00+=fC00; r01+=fC10; r11+=fC11;
405 Double_t det=r00*r11 - r01*r01;
406 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
407
408 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
409 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
410 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
411 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
412 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
413
414 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
415 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1627d1c4 416 if (TMath::Abs(cur*fX-eta) >= 0.9) {
417 // Int_t n=GetNumberOfClusters();
418 //if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
1c53abe2 419 return 0;
420 }
421
422 fP0 += k00*dy + k01*dz;
423 fP1 += k10*dy + k11*dz;
424 fP2 = eta;
425 fP3 += k30*dy + k31*dz;
426 fP4 = cur;
427
428 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
429 Double_t c12=fC21, c13=fC31, c14=fC41;
430
431 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
432 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
433 fC40-=k00*c04+k01*c14;
434
435 fC11-=k10*c01+k11*fC11;
436 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
437 fC41-=k10*c04+k11*c14;
438
439 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
440 fC42-=k20*c04+k21*c14;
441
442 fC33-=k30*c03+k31*c13;
443 fC43-=k40*c03+k41*c13;
444
445 fC44-=k40*c04+k41*c14;
446
447 Int_t n=GetNumberOfClusters();
448 fIndex[n]=index;
449 SetNumberOfClusters(n+1);
450 SetChi2(GetChi2()+chisq);
451
452 return 1;
453}
454
455
456
457//_____________________________________________________________________________
458Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
459 Double_t x2,Double_t y2,
460 Double_t x3,Double_t y3)
461{
462 //-----------------------------------------------------------------
463 // Initial approximation of the track curvature
464 //-----------------------------------------------------------------
465 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
466 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
467 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
468 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
469 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
470
471 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
472
473 return -xr*yr/sqrt(xr*xr+yr*yr);
474}
475
476
477//_____________________________________________________________________________
478Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
479 Double_t x2,Double_t y2,
480 Double_t x3,Double_t y3)
481{
482 //-----------------------------------------------------------------
483 // Initial approximation of the track curvature times center of curvature
484 //-----------------------------------------------------------------
485 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
486 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
487 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
488 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
489 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
490
491 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
492
493 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
494}
495
496//_____________________________________________________________________________
497Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
498 Double_t x2,Double_t y2,
499 Double_t z1,Double_t z2)
500{
501 //-----------------------------------------------------------------
502 // Initial approximation of the tangent of the track dip angle
503 //-----------------------------------------------------------------
504 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
505}
506
507
508void AliTPCtrackerMI::LoadClusters()
509{
510 //
511 // load clusters to the memory
512 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
513 for (Int_t i=0; i<j; i++) {
514 fClustersArray.LoadEntry(i);
515 }
516}
517
518void AliTPCtrackerMI::UnloadClusters()
519{
520 //
521 // load clusters to the memory
522 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
523 for (Int_t i=0; i<j; i++) {
524 fClustersArray.ClearSegment(i);
525 }
526}
527
528
529
530//_____________________________________________________________________________
531void AliTPCtrackerMI::LoadOuterSectors() {
532 //-----------------------------------------------------------------
533 // This function fills outer TPC sectors with clusters.
534 //-----------------------------------------------------------------
535 UInt_t index;
9918f10a 536 //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
537 Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
1c53abe2 538 for (Int_t i=0; i<j; i++) {
539 // AliSegmentID *s=fClustersArray.LoadEntry(i);
cc5e9db0 540 AliSegmentID *s= const_cast<AliSegmentID*>(fClustersArray.At(i));
9918f10a 541 if (!s) continue;
1c53abe2 542 Int_t sec,row;
543 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
544 par->AdjustSectorRow(s->GetID(),sec,row);
545 if (sec<fkNIS*2) continue;
546 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
547 Int_t ncl=clrow->GetArray()->GetEntriesFast();
548 while (ncl--) {
549 AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
550 index=(((sec<<8)+row)<<16)+ncl;
551 fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
552 }
553 }
554 fN=fkNOS;
555 fSectors=fOuterSec;
556}
557
558
559//_____________________________________________________________________________
560void AliTPCtrackerMI::LoadInnerSectors() {
561 //-----------------------------------------------------------------
562 // This function fills inner TPC sectors with clusters.
563 //-----------------------------------------------------------------
564 UInt_t index;
9918f10a 565 //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
566 Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
1c53abe2 567 for (Int_t i=0; i<j; i++) {
568 // AliSegmentID *s=fClustersArray.LoadEntry(i);
cc5e9db0 569 AliSegmentID *s=const_cast<AliSegmentID*>(fClustersArray.At(i));
9918f10a 570 if (!s) continue;
1c53abe2 571 Int_t sec,row;
572 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
573 par->AdjustSectorRow(s->GetID(),sec,row);
574 if (sec>=fkNIS*2) continue;
575 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
576 Int_t ncl=clrow->GetArray()->GetEntriesFast();
577 while (ncl--) {
578 AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
579 index=(((sec<<8)+row)<<16)+ncl;
580 fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
581 }
582 }
583
584 fN=fkNIS;
585 fSectors=fInnerSec;
586}
587
588Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
589 //-----------------------------------------------------------------
590 // This function tries to find a track prolongation to next pad row
591 //-----------------------------------------------------------------
592 // Double_t xt=t.GetX();
593 // Int_t row = fSectors->GetRowNumber(xt)-1;
594 // if (row < nr) return 1; // don't prolongate if not information until now -
595 //
596 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
1627d1c4 597 // if (t.GetRadius()>x+10 ) return 0;
598
599 if (!t.PropagateTo(x)) {
600 t.fStopped = kTRUE;
601 return 0;
602 }
1c53abe2 603 // update current
604 t.fCurrentSigmaY = GetSigmaY(&t);
605 t.fCurrentSigmaZ = GetSigmaZ(&t);
606 //
607 AliTPCclusterMI *cl=0;
608 UInt_t index=0;
609 const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
610 Double_t sy2=ErrY2(&t)*2;
611 Double_t sz2=ErrZ2(&t)*2;
612
613
614 Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
615 Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
616 Double_t y=t.GetY(), z=t.GetZ();
617
618 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
619 t.fInDead = kTRUE;
620 return 0;
621 }
622 else
623 {
1627d1c4 624 if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
625 else
626 return 0;
1c53abe2 627 }
628 //calculate
629 Float_t maxdistance = roady*roady + roadz*roadz;
630 if (krow) {
631 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
632 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
633 if (c->GetZ() > z+roadz) break;
634 if ( (c->GetY()-y) > roady ) continue;
635 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
636 if (maxdistance>distance) {
637 maxdistance = distance;
638 cl=c;
639 index=krow.GetIndex(i);
640 }
641 }
642 }
643 if (cl) {
644 //Double_t sy2=
645 ErrY2(&t,cl);
646 //Double_t sz2=
647 ErrZ2(&t,cl);
648 Double_t chi2= t.GetPredictedChi2(cl);
649 UpdateTrack(&t,cl,chi2,index);
650 } else {
651 if (y > ymax) {
652 t.fRelativeSector= (t.fRelativeSector+1) % fN;
653 if (!t.Rotate(fSectors->GetAlpha()))
654 return 0;
655 } else if (y <-ymax) {
656 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
657 if (!t.Rotate(-fSectors->GetAlpha()))
658 return 0;
659 }
660 }
661 return 1;
662}
663
664
665Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex, Int_t nr) {
666 //-----------------------------------------------------------------
667 // This function tries to find a track prolongation to next pad row
668 //-----------------------------------------------------------------
669 t.fCurrentCluster = 0;
670 t.fCurrentClusterIndex1 = 0;
671 t.fCurrentClusterIndex2 = 0;
672
673 Double_t xt=t.GetX();
674 Int_t row = fSectors->GetRowNumber(xt)-1;
675 if (row < nr) return 1; // don't prolongate if not information until now -
676 Double_t x=fSectors->GetX(nr);
1627d1c4 677 // if (t.fStopped) return 0;
678 // if (t.GetRadius()>x+10 ) return 0;
679 if (!t.PropagateTo(x)){
680 t.fStopped =kTRUE;
681 return 0;
682 }
1c53abe2 683 // update current
684 t.fCurrentSigmaY = GetSigmaY(&t);
685 t.fCurrentSigmaZ = GetSigmaZ(&t);
686
687 AliTPCclusterMI *cl=0;
688 UInt_t index=0;
689 AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
690 //
691 Double_t y=t.GetY(), z=t.GetZ();
692 Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
693 Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
694 //
695
696 Float_t maxdistance = 1000000;
697 if (krow) {
698 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
699 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
700 if (c->GetZ() > z+roadz) break;
701 if (TMath::Abs(c->GetY()-y)>roady) continue;
702
703 //krow.UpdateClusterTrack(i,trindex,&t);
704
705 Float_t dy2 = (c->GetY()- t.GetY());
706 dy2*=dy2;
707 Float_t dz2 = (c->GetZ()- t.GetZ());
708 dz2*=dz2;
709 //
710 Float_t distance = dy2+dz2;
711 //
712 if (distance > maxdistance) continue;
713 maxdistance = distance;
714 cl=c;
715 index=i;
716 }
717 }
718 t.fCurrentCluster = cl;
719 t.fCurrentClusterIndex1 = krow.GetIndex(index);
720 t.fCurrentClusterIndex2 = index;
721 return 1;
722}
723
724
725Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
726 //-----------------------------------------------------------------
727 // This function tries to find a track prolongation to next pad row
728 //-----------------------------------------------------------------
729 AliTPCseed & t = *((AliTPCseed*)(fSeeds->At(trindex)));
730 AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
731 // Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
732 Double_t y=t.GetY();
733 Double_t ymax=fSectors->GetMaxY(nr);
734
735 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
736 t.fInDead = kTRUE;
737 return 0;
738 }
739 else
740 {
1627d1c4 741 if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
742 else
743 return 0;
1c53abe2 744 }
745
746 if (t.fCurrentCluster) {
747 // Float_t l=fSectors->GetPadPitchLength();
748 // AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
749
750 Double_t sy2=ErrY2(&t,t.fCurrentCluster);
751 Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
752
753
754 Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
755 Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
756
757 Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
758 Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
759
760 Double_t rdistance = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
761
762
763 // printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
1627d1c4 764 if ( (rdistancey>1) || (rdistancez>1)) return 0;
1c53abe2 765 if (rdistance>4) return 0;
766
767 if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
768 return 0; //suspisiouce - will be changed
769
770 if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
771 // strict cut on overlaped cluster
772 return 0; //suspisiouce - will be changed
773
774 if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
775 && t.fCurrentCluster->GetType()<0)
776 return 0;
777
778 // t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
779 UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
780
781 } else {
782 if (y > ymax) {
783 t.fRelativeSector= (t.fRelativeSector+1) % fN;
784 if (!t.Rotate(fSectors->GetAlpha()))
785 return 0;
786 } else if (y <-ymax) {
787 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
788 if (!t.Rotate(-fSectors->GetAlpha()))
789 return 0;
790 }
791 }
792 return 1;
793}
794
795
796
797
798/*
799Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
800{
801 //-----------------------------------------------------------------
802 // fast prolongation mathod -
803 // don't update track only after step clusters
804 //-----------------------------------------------------------------
805 Double_t xt=t.GetX();
806 //
807 Double_t alpha=t.GetAlpha();
808 alpha =- fSectors->GetAlphaShift();
809 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
810 if (alpha < 0. ) alpha += 2.*TMath::Pi();
811 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
812 Int_t row0 = fSectors->GetRowNumber(xt);
813 Double_t x = fSectors->GetX(row0);
814 Double_t ymax = fSectors->GetMaxY(row0);
815 //
816 Double_t sy2=ErrY2(&t)*2;
817 Double_t sz2=ErrZ2(&t)*2;
818 Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
819 Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
820 Float_t maxdistance = roady*roady + roadz*roadz;
821 t.fCurrentSigmaY = GetSigmaY(&t);
822 t.fCurrentSigmaZ = GetSigmaZ(&t);
823 //
824 Int_t nclusters = 0;
825 Double_t y;
826 Double_t z;
827 Double_t yy[200]; //track prolongation
828 Double_t zz[200];
829 Double_t cy[200]; // founded cluster position
830 Double_t cz[200];
831 Double_t sy[200]; // founded cluster error
832 Double_t sz[200];
833 Bool_t hitted[200]; // indication of cluster presence
834 //
835
836 //
837 for (Int_t drow = step; drow>=0; drow--) {
838 Int_t row = row0-drow;
839 if (row<0) break;
840 Double_t x = fSectors->GetX(row);
841 Double_t ymax = fSectors->GetMaxY(row);
842 t.GetProlongation(x,y,z);
843 yy[drow] =y;
844 zz[drow] =z;
845 const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
846 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
847 t.fInDead = kTRUE;
848 break;
849 }
850 else
851 {
852 t.fNFoundable++;
853 }
854
855 //find nearest cluster
856 AliTPCclusterMI *cl= 0;
857 if (krow) {
858 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
859 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
860 if (c->GetZ() > z+roadz) break;
861 if ( (c->GetY()-y) > roady ) continue;
862 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
863 if (maxdistance>distance) {
864 maxdistance = distance;
865 cl=c;
866 // index=krow.GetIndex(i);
867 }
868 }
869 } // end of seearch
870 //update cluster information
871 if (cl){
872 cy[drow] = cl->GetY();
873 cz[drow] = cl->GetZ();
874 sy[drow] = ErrY2(&t,cl);
875 sz[drow] = ErrZ2(&t,cl);
876 hitted[drow] = kTRUE;
877 nclusters++;
878 }
879 else
880 hitted[drow] = kFALSE;
881 }
882 //if we have information - update track
883 if (nclusters>0){
884 Float_t sumyw0 = 0;
885 Float_t sumzw0 = 0;
886 Float_t sumyw = 0;
887 Float_t sumzw = 0;
888 Float_t sumrow = 0;
889 for (Int_t i=0;i<step;i++)
890 if (hitted[i]){
891 sumyw0+= 1/sy[i];
892 sumzw0+= 1/sz[i];
893 //
894 sumyw+= (cy[i]-yy[i])/sy[i];
895 sumzw+= (cz[i]-zz[i])/sz[i];
896 sumrow+= i;
897 }
898 Float_t dy = sumyw/sumyw0;
899 Float_t dz = sumzw/sumzw0;
900 Float_t mrow = sumrow/nclusters+row0;
901 Float_t x = fSectors->GetX(mrow);
902 t.PropagateTo(x);
903 AliTPCclusterMI cvirtual;
904 cvirtual.SetZ(dz+t.GetZ());
905 cvirtual.SetY(dy+t.GetY());
906 t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
907 t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
908 Float_t chi2 = t.GetPredictedChi2(&cvirtual);
909 t.Update(&cvirtual,chi2,0);
910 Int_t ncl = t.GetNumberOfClusters();
911 ncl = ncl-1+nclusters;
912 t.SetN(ncl);
913 }
914 return 1;
915}
916*/
917
918//_____________________________________________________________________________
919Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
920 //-----------------------------------------------------------------
921 // This function tries to find a track prolongation.
922 //-----------------------------------------------------------------
923 Double_t xt=t.GetX();
924 //
925 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
926 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
927 if (alpha < 0. ) alpha += 2.*TMath::Pi();
928 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
929
930 for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
931
932 if (FollowToNext(t,nr)==0) {
933 }
934 }
935 return 1;
936}
937
938
939//_____________________________________________________________________________
940Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
941 //-----------------------------------------------------------------
942 // This function tries to find a track prolongation.
943 //-----------------------------------------------------------------
944 Double_t xt=t.GetX();
945 //
946 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
947 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
948 if (alpha < 0. ) alpha += 2.*TMath::Pi();
949 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
950
951 for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
1627d1c4 952 if (t.GetSnp()<0.8)
1c53abe2 953 FollowToNext(t,nr);
954 }
955 return 1;
956}
957
958
959
960
961
962Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
963{
964 //
965 //
966 sum1=0;
967 sum2=0;
968 Int_t sum=0;
969 if (s1->fSector!=s2->fSector) return 0;
970 //
971 Float_t dz2 =(s1->GetZ() - s2->GetZ());
972 dz2*=dz2;
973 Float_t dy2 =(s1->GetY() - s2->GetY());
974 dy2*=dy2;
975 Float_t distance = TMath::Sqrt(dz2+dy2);
976 if (distance>5.) return 0; // if there are far away - not overlap - to reduce combinatorics
977
978 for (Int_t i=0;i<160;i++){
979 if (s1->fClusterIndex[i]>0) sum1++;
980 if (s2->fClusterIndex[i]>0) sum2++;
981 if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
982 sum++;
983 }
984 }
985
1627d1c4 986 Float_t summin = TMath::Min(sum1+1,sum2+1);
987 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 988 return ratio;
989}
990
991void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
992{
993 //
994 //
995 if (s1->fSector!=s2->fSector) return;
996 //
997 Float_t dz2 =(s1->GetZ() - s2->GetZ());
998 dz2*=dz2;
999 Float_t dy2 =(s1->GetY() - s2->GetY());
1000 dy2*=dy2;
1001 Float_t distance = TMath::Sqrt(dz2+dy2);
1002 if (distance>15.) return ; // if there are far away - not overlap - to reduce combinatorics
1003 //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
1004 // TClonesArray &pointarray1 = *(s1->fPoints);
1005 //TClonesArray &pointarray2 = *(s2->fPoints);
1006 //
1007 for (Int_t i=0;i<160;i++){
1008 if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
1009 // AliTPCTrackPoint *p1 = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
1010 //AliTPCTrackPoint *p2 = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i));
1011 AliTPCTrackPoint *p1 = s1->GetTrackPoint(i);
1012 AliTPCTrackPoint *p2 = s2->GetTrackPoint(i);;
1013 p1->fIsShared = kTRUE;
1014 p2->fIsShared = kTRUE;
1015 }
1016 }
1017}
1018
1019
1020
1021
cc5e9db0 1022void AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared){
1c53abe2 1023
1024
1025
1026 //
1027 // remove overlap - used removal factor - removal index stored in the track
1028 arr->Sort(); // sorting according z
1029 arr->Expand(arr->GetEntries());
1030 Int_t nseed=arr->GetEntriesFast();
1031 // printf("seeds \t%p \t%d\n",arr, nseed);
1032 // arr->Expand(arr->GetEntries()); //remove 0 pointers
1033 nseed = arr->GetEntriesFast();
1034 Int_t removed = 0;
1035 for (Int_t i=0; i<nseed; i++) {
1036 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1037 if (!pt) {
1038 continue;
1039 }
1040 if (!(pt->IsActive())) continue;
1041 for (Int_t j=i+1; j<nseed; j++){
1042 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1043 if ((pt2) && pt2->IsActive())
1044 if (pt->fSector == pt2->fSector)
1045 if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
1046 Int_t sum1,sum2;
1047 Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1627d1c4 1048 //if (sum1==0) {
1049 // pt->Desactivate(removalindex); // arr->RemoveAt(i);
1050 // break;
1051 //}
1c53abe2 1052 if (ratio>factor){
1053 // if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex); // arr->RemoveAt(j);
1054 Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1055 Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1056 removed++;
1057 if (TMath::Abs(ratio3)>0.025){ // if much more points
1058 if (sum1>sum2) pt2->Desactivate(removalindex);
1059 else {
1060 pt->Desactivate(removalindex); // arr->RemoveAt(i);
1061 break;
1062 }
1063 }
1064 else{ //decide on mean chi2
1065 if (ratio2<1)
1066 pt2->Desactivate(removalindex);
1067 else {
1068 pt->Desactivate(removalindex); // arr->RemoveAt(i);
1069 break;
1070 }
1071 }
1072
1073 } // if suspicious ratio
1074 }
1075 else
1076 break;
1077 }
1078 }
1079 // printf("removed\t%d\n",removed);
1080 Int_t good =0;
1081 for (Int_t i=0; i<nseed; i++) {
1082 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1083 if (!pt) break;
1084 if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1085 //desactivate tracks with small number of points
1086 // printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1087 pt->Desactivate(10); //desactivate - small muber of points
1088 }
1089 if (!(pt->IsActive())) continue;
1090 good++;
1091 }
1092
1093
1094 if (shared)
1095 for (Int_t i=0; i<nseed; i++) {
1096 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1097 if (!pt) continue;
1098 if (!(pt->IsActive())) continue;
1099 for (Int_t j=i+1; j<nseed; j++){
1100 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1101 if ((pt2) && pt2->IsActive()) {
1102 if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1103 SignShared(pt,pt2);
1104 }
1105 }
1106 }
1107 fNtracks = good;
1108 printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1109
1110
1111}
1112void AliTPCtrackerMI::MakeSeedsAll()
1113{
1114 if (fSeeds == 0) fSeeds = new TObjArray;
1115 TObjArray * arr;
1116 for (Int_t sec=0;sec<fkNOS;sec+=3){
1117 arr = MakeSeedsSectors(sec,sec+3);
1118 Int_t nseed = arr->GetEntriesFast();
1119 for (Int_t i=0;i<nseed;i++)
1120 fSeeds->AddLast(arr->RemoveAt(i));
1627d1c4 1121 }
1122 // fSeeds = MakeSeedsSectors(0,fkNOS);
1c53abe2 1123}
1124
1125TObjArray * AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1126{
1127 //
1128 // loop over all sectors and make seed
1129 //find track seeds
1130 TStopwatch timer;
1131 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1132 Int_t nrows=nlow+nup;
1133 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1134 // if (fSeeds==0) fSeeds = new TObjArray;
1135 TObjArray * arr = new TObjArray;
1136
1137 for (Int_t sec=sec1; sec<sec2;sec++){
1138 MakeSeeds(arr, sec, nup-1, nup-1-gap);
1139 MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
1140 }
1627d1c4 1141 gap = Int_t(0.3* nrows);
1142 for (Int_t sec=sec1; sec<sec2;sec++){
1143 //find secondaries
1144 MakeSeeds2(arr, sec, nup-1, nup-1-gap);
1145 MakeSeeds2(arr, sec, nup-1-shift, nup-1-shift-gap);
1146 MakeSeeds2(arr, sec, nup-1-2*shift, nup-1-2*shift-gap);
1147 //MakeSeeds2(arr, sec, nup-1-3*shift, nup-1-3*shift-gap);
1148 MakeSeeds2(arr, sec, 30, 0);
1149 }
1150
1c53abe2 1151 Int_t nseed=arr->GetEntriesFast();
1152 gap=Int_t(0.3*nrows);
1153 // continue seeds
1154 Int_t i;
1155 for (i=0; i<nseed; i++) {
1156 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1157 if (!pt) continue;
1158 if (FollowProlongation(t,nup-gap)) {
1159 pt->fIsSeeding =kFALSE;
1160 continue;
1161 }
1162 delete arr->RemoveAt(i);
1627d1c4 1163 }
1164
1c53abe2 1165 //
1166 //remove seeds which overlaps
1167 RemoveOverlap(arr,0.6,1);
1168 //delete seeds - which were sign
1169 nseed=arr->GetEntriesFast();
1170 for (i=0; i<nseed; i++) {
1171 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1172 //, &t=*pt;
1173 if (!pt) continue;
1174 if ((pt->IsActive()) && pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1627d1c4 1175 //pt->Reset();
1176 //FollowBackProlongation(*pt,nup-1);
1177 //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 )
1178 //delete arr->RemoveAt(i);
1179 //else
1180 // pt->Reset();
1c53abe2 1181 continue;
1182 }
1183 delete arr->RemoveAt(i);
1627d1c4 1184 }
1185 //RemoveOverlap(arr,0.6,1);
1c53abe2 1186 return arr;
1187}
1188
1189
1190
1191//_____________________________________________________________________________
1192void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1193 //-----------------------------------------------------------------
1194 // This function creates track seeds.
1195 //-----------------------------------------------------------------
1196 // if (fSeeds==0) fSeeds=new TObjArray(15000);
1197
1198 Double_t x[5], c[15];
1199
1200 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1201 Double_t cs=cos(alpha), sn=sin(alpha);
1202
1203 Double_t x1 =fOuterSec->GetX(i1);
1204 Double_t xx2=fOuterSec->GetX(i2);
1205 //
1206 // for (Int_t ns=0; ns<fkNOS; ns++)
1207 Int_t ns =sec;
1208 {
1209 Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1210 Int_t nm=fOuterSec[ns][i2];
1211 Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1212 const AliTPCRow& kr1=fOuterSec[ns][i1];
1213 AliTPCRow& kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1214 AliTPCRow& kr22 = fOuterSec[(ns)%fkNOS][i2];
1215 AliTPCRow& kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1216
1217 for (Int_t is=0; is < kr1; is++) {
1218 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1219 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1220
1221 Float_t anglez = (z1-z3)/(x1-x3);
1222 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
1223
1224 for (Int_t js=0; js < nl+nm+nu; js++) {
1225 const AliTPCclusterMI *kcl;
1226 Double_t x2, y2, z2;
1227 if (js<nl) {
1228 if (js==0) {
1229 js = kr21.Find(extraz-15.);
1230 if (js>=nl) continue;
1231 }
1232 kcl=kr21[js];
1233 z2=kcl->GetZ();
1234 if ((extraz-z2)>10) continue;
1235 if ((extraz-z2)<-10) {
ab9cc56f 1236 js = nl-1;
1c53abe2 1237 continue;
1238 }
1239 y2=kcl->GetY();
1240 x2= xx2*cs+y2*sn;
1241 y2=-xx2*sn+y2*cs;
1242 } else
1243 if (js<nl+nm) {
1244 if (js==nl) {
1245 js = nl+kr22.Find(extraz-15.);
1246 if (js>=nl+nm) continue;
1247 }
1248 kcl=kr22[js-nl];
1249 z2=kcl->GetZ();
1250 if ((extraz-z2)>10) continue;
1251 if ((extraz-z2)<-10) {
ab9cc56f 1252 js = nl+nm-1;
1c53abe2 1253 continue;
1254 }
1255 x2=xx2; y2=kcl->GetY();
1256 } else {
1257 //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
1258 if (js==nl+nm) {
1259 js = nl+nm+kr23.Find(extraz-15.);
1260 if (js>=nl+nm+nu) break;
1261 }
1262 kcl=kr23[js-nl-nm];
1263 z2=kcl->GetZ();
1264 if ((extraz-z2)>10) continue;
1265 if ((extraz-z2)<-10) {
1266 break;
1267 }
1268 y2=kcl->GetY();
1269 x2=xx2*cs-y2*sn;
1270 y2=xx2*sn+y2*cs;
1271 }
1272
1273 Double_t zz=z1 - anglez*(x1-x2);
1274 if (TMath::Abs(zz-z2)>10.) continue;
1275
1276 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1277 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1278
1279 x[0]=y1;
1280 x[1]=z1;
1281 x[4]=f1(x1,y1,x2,y2,x3,y3);
1282 if (TMath::Abs(x[4]) >= 0.0066) continue;
1283 x[2]=f2(x1,y1,x2,y2,x3,y3);
1284 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1285 x[3]=f3(x1,y1,x2,y2,z1,z2);
1286 if (TMath::Abs(x[3]) > 1.2) continue;
1287 Double_t a=asin(x[2]);
1288 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1289 if (TMath::Abs(zv-z3)>10.) continue;
1290
1291 Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1292 Double_t sy2=kcl->GetSigmaY2()*2, sz2=kcl->GetSigmaZ2()*4;
1293 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1294 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1295 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1296
1297 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1298 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1299 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1300 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1301 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1302 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1303 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1304 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1305 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1306 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1307
1308 c[0]=sy1;
1309 c[1]=0.; c[2]=sz1;
1310 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1311 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
1312 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1313 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1314 c[13]=f30*sy1*f40+f32*sy2*f42;
1315 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1316
1317 UInt_t index=kr1.GetIndex(is);
1318 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1319 track->fIsSeeding = kTRUE;
1320 Int_t rc=FollowProlongation(*track, i2);
1321 //FollowProlongationFast(*track, 5);
1322 //FollowProlongationFast(*track, 5);
1323 //FollowProlongationFast(*track, 5);
1324 //FollowProlongationFast(*track, 5);
1325 //Int_t rc = 1;
1326
1327 track->fLastPoint = i1; // first cluster in track position
1328 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1329 else arr->AddLast(track);
1330 }
1331 }
1332 }
1333}
1334
1627d1c4 1335
1336//_____________________________________________________________________________
1337void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1338 //-----------------------------------------------------------------
1339 // This function creates track seeds - without vertex constraint
1340 //-----------------------------------------------------------------
1341
1342 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1343 // Double_t cs=cos(alpha), sn=sin(alpha);
1344 Int_t row0 = (i1+i2)/2;
1345 Int_t drow = (i1-i2)/2;
1346 const AliTPCRow& kr0=fSectors[sec][row0];
1347 const AliTPCRow& krm=fSectors[sec][row0-1];
1348 const AliTPCRow& krp=fSectors[sec][row0+1];
1349 AliTPCRow * kr=0;
1350
1351 AliTPCpolyTrack polytrack;
1352 Int_t nclusters=fSectors[sec][row0];
1353
1354 for (Int_t is=0; is < nclusters; is++) {
1355 const AliTPCclusterMI * cl= kr0[is];
1356 Double_t x = kr0.GetX();
1357
1358 // Initialization of the polytrack
1359 polytrack.Reset();
1360
1361 Double_t y0= cl->GetY();
1362 Double_t z0= cl->GetZ();
1363 polytrack.AddPoint(x,y0,z0);
1364 Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2);
1365 Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2);
1366 //
1367 x = krm.GetX();
1368 cl = krm.FindNearest(y0,z0,roady,roadz);
1369 if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz);
1370 //
1371 x = krp.GetX();
1372 cl = krp.FindNearest(y0,z0,roady,roadz);
1373 if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05);
1374 //
1375 polytrack.UpdateParameters();
1376 // follow polytrack
1377 roadz = 0.6;
1378 roady = 0.6;
1379 //
1380 Double_t yn,zn;
1381 Int_t nfoundable = polytrack.GetN();
1382 Int_t nfound = nfoundable;
1383 for (Int_t ddrow = 2; ddrow<drow;ddrow++){
1384 for (Int_t delta = -1;delta<=1;delta+=2){
1385 Int_t row = row0+ddrow*delta;
1386 kr = &(fSectors[sec][row]);
1387 Double_t xn = kr->GetX();
1388 Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone;
1389 polytrack.GetFitPoint(xn,yn,zn);
1390 if (TMath::Abs(yn)>ymax) continue;
1391 nfoundable++;
1392 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
1393 if (cln) {
1394 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05);
1395 nfound++;
1396 }
1397 }
1398 polytrack.UpdateParameters();
1399 }
1400 if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
1401 // add polytrack candidate
1402 Double_t x[5], c[15];
1403 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
1404 polytrack.GetBoundaries(x3,x1);
1405 x2 = (x1+x3)/2.;
1406 polytrack.GetFitPoint(x1,y1,z1);
1407 polytrack.GetFitPoint(x2,y2,z2);
1408 polytrack.GetFitPoint(x3,y3,z3);
1409 //
1410 //is track pointing to the vertex ?
1411 Double_t x0,y0,z0;
1412 x0=0;
1413 polytrack.GetFitPoint(x0,y0,z0);
1414 if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint
1415 x3 = 0;
1416 y3 = GetY();
1417 z3 = GetZ();
1418 }
1419
1420 x[0]=y1;
1421 x[1]=z1;
1422 x[4]=f1(x1,y1,x2,y2,x3,y3);
1423 if (TMath::Abs(x[4]) >= 0.0066) continue;
1424 x[2]=f2(x1,y1,x2,y2,x3,y3);
1425 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1426 x[3]=f3(x1,y1,x2,y2,z1,z2);
1427 if (TMath::Abs(x[3]) > 1.2) continue;
1428 if (TMath::Abs(x[2]) > 0.99) continue;
1429 // Double_t a=asin(x[2]);
1430
1431
1432 Double_t sy=1.5, sz=1.5;
1433 Double_t sy1=1.5, sz1=1.5;
1434 Double_t sy2=1.3, sz2=1.3;
1435 Double_t sy3=1.5;
1436 //sz3=1.5;
1437 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1438 // Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1439 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1440
1441 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1442 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1443 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1444 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1445 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1446 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1447 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1448 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1449 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1450 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1451
1452 c[0]=sy1;
1453 c[1]=0.; c[2]=sz1;
1454 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1455 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
1456 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1457 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1458 c[13]=f30*sy1*f40+f32*sy2*f42;
1459 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1460
1461 UInt_t index=0;
1462 //kr0.GetIndex(is);
1463 AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift);
1464 track->fStopped =kFALSE;
1465 track->fIsSeeding = kTRUE;
1466 Int_t rc=FollowProlongation(*track, i2);
1467 track->fLastPoint = i1; // first cluster in track position
1468 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1469 else arr->AddLast(track);
1470 }
1471 }
1472
1473
1474
1475}
1476
1477
1478
1479
1480
1c53abe2 1481//_____________________________________________________________________________
1482Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1483 //-----------------------------------------------------------------
1484 // This function reades track seeds.
1485 //-----------------------------------------------------------------
1486 TDirectory *savedir=gDirectory;
1487
1488 TFile *in=(TFile*)inp;
1489 if (!in->IsOpen()) {
1490 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1491 return 1;
1492 }
1493
1494 in->cd();
1495 TTree *seedTree=(TTree*)in->Get("Seeds");
1496 if (!seedTree) {
1497 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1498 cerr<<"can't get a tree with track seeds !\n";
1499 return 2;
1500 }
1501 AliTPCtrack *seed=new AliTPCtrack;
1502 seedTree->SetBranchAddress("tracks",&seed);
1503
1504 if (fSeeds==0) fSeeds=new TObjArray(15000);
1505
1506 Int_t n=(Int_t)seedTree->GetEntries();
1507 for (Int_t i=0; i<n; i++) {
1508 seedTree->GetEvent(i);
1509 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1510 }
1511
1512 delete seed;
1513 delete seedTree;
1514 savedir->cd();
1515 return 0;
1516}
1517
1518//_____________________________________________________________________________
1519Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1520 //-----------------------------------------------------------------
1521 // This is a track finder.
1522 //-----------------------------------------------------------------
1523 TDirectory *savedir=gDirectory;
1524
1525 if (inp) {
1526 TFile *in=(TFile*)inp;
1527 if (!in->IsOpen()) {
1528 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1529 return 1;
1530 }
1531 }
1532
1533 if (!out->IsOpen()) {
1534 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1535 return 2;
1536 }
1537
1538 out->cd();
1539
1540 char tname[100];
1541 sprintf(tname,"TreeT_TPC_%d",fEventN);
1542 TTree tracktree(tname,"Tree with TPC tracks");
1543 TTree seedtree("Seeds","Seeds");
1544 AliTPCtrack *iotrack=0;
1545 AliTPCseed *ioseed=0;
1546 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1547 TStopwatch timer;
1548
1549 printf("Loading clusters \n");
1550 LoadClusters();
1551 printf("Time for loading clusters: \t");timer.Print();timer.Start();
1552
1553 printf("Loading outer sectors\n");
1554 LoadOuterSectors();
1555 printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1556
1557 printf("Loading inner sectors\n");
1558 LoadInnerSectors();
1559 printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1560 fSectors = fOuterSec;
1561 fN=fkNOS;
1562
1563
1564 //find track seeds
1565 MakeSeedsAll();
1566 printf("Time for seeding: \t"); timer.Print();timer.Start();
1567 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1568 Int_t nrows=nlow+nup;
1569
1570 Int_t gap=Int_t(0.3*nrows);
1571 Int_t i;
1627d1c4 1572 //RemoveOverlap(fSeeds,0.6,2);
1c53abe2 1573 Int_t nseed=fSeeds->GetEntriesFast();
1574 // outer sectors parallel tracking
1575 ParallelTracking(fSectors->GetNRows()-gap-1,0);
1627d1c4 1576 //ParallelTracking(fSectors->GetNRows()-1,0);
1577 //RemoveOverlap(fSeeds, 0.6,3);
1578 // ParallelTracking(49,0);
1c53abe2 1579 printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1580
1581 RemoveOverlap(fSeeds, 0.6,3);
1582 printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1583 //parallel tracking
1584 fSectors = fInnerSec;
1585 fN=fkNIS;
1586 ParallelTracking(fSectors->GetNRows()-1,0);
1587 printf("Number of tracks after inner tracking %d\n",fNtracks);
1588 printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1589 //
1590 RemoveOverlap(fSeeds,0.6,5,kTRUE); // remove overlap - shared points signed
1591 printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1592 //
1593 //
1594
1595 ioseed = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1596 AliTPCseed * vseed = new AliTPCseed;
1597 vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1598 vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1599 vseed->fPoints->ExpandCreateFast(2);
1600
1601 TBranch * seedbranch = seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1602 //delete vseed;
1603 nseed=fSeeds->GetEntriesFast();
1604
1605 Int_t found = 0;
1606 for (i=0; i<nseed; i++) {
1607 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1608 if (!pt) continue;
1609 Int_t nc=t.GetNumberOfClusters();
1627d1c4 1610 if (nc<20) continue;
1c53abe2 1611 t.CookdEdx(0.02,0.6);
1612 CookLabel(pt,0.1); //For comparison only
1613 // if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1627d1c4 1614 if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){
1c53abe2 1615 iotrack=pt;
1616 tracktree.Fill();
1617 cerr<<found++<<'\r';
1618 }
1619 else
1620 if ( (pt->IsActive())) fNtracks--;
1621 pt->RebuildSeed();
1622 seedbranch->SetAddress(&pt);
1623
1624 seedtree.Fill();
1625 for (Int_t j=0;j<160;j++){
1626 delete pt->fPoints->RemoveAt(j);
1627 }
1628 delete pt->fPoints;
1629 pt->fPoints =0;
1630 delete fSeeds->RemoveAt(i);
1631 }
1632 printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1633
1634 UnloadClusters();
1635 printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1636
1637 tracktree.Write();
1638 seedtree.Write();
1639 cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
1640
1641 savedir->cd();
1642
1643 return 0;
1644}
1645
1646
1647void AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1648{
1649 //
1650 // try to track in parralel
1651
1652 Int_t nseed=fSeeds->GetEntriesFast();
1653 //prepare seeds for tracking
1654 for (Int_t i=0; i<nseed; i++) {
1655 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1656 if (!pt) continue;
1657 if (!t.IsActive()) continue;
1658 // follow prolongation to the first layer
1659 FollowProlongation(t, rfirst+1);
1660 }
1661
1662
1663 //
1664 for (Int_t nr=rfirst; nr>=rlast; nr--){
1665 // make indexes with the cluster tracks for given
1666 // for (Int_t i = 0;i<fN;i++)
1667 // fSectors[i][nr].MakeClusterTracks();
1668
1669 // find nearest cluster
1670 for (Int_t i=0; i<nseed; i++) {
1671 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1672 if (!pt) continue;
1673 if (!pt->IsActive()) continue;
1627d1c4 1674 if (pt->fRelativeSector>17) {
1675 continue;
1676 }
1c53abe2 1677 UpdateClusters(t,i,nr);
1678 }
1679 // prolonagate to the nearest cluster - if founded
1680 for (Int_t i=0; i<nseed; i++) {
1681 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
1682 if (!pt) continue;
1627d1c4 1683 if (!pt->IsActive()) continue;
1684 if (pt->fRelativeSector>17) {
1685 continue;
1686 }
1c53abe2 1687 FollowToNextCluster(i,nr);
1688 }
1689 // for (Int_t i= 0;i<fN;i++)
1690 // fSectors[i][nr].ClearClusterTracks();
1691 }
1692
1693}
1694
1695Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1696{
1697 //
1698 //
1699 Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1700 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
1701 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1702 Float_t angular = seed->GetSnp();
1703 angular = angular*angular/(1-angular*angular);
1704 // angular*=angular;
1705 //angular = TMath::Sqrt(angular/(1-angular));
1706 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1707 return res;
1708}
1709Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1710{
1711 //
1712 //
1713 Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1714 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
1715 Float_t sres = fParam->GetZSigma();
1716 Float_t angular = seed->GetTgl();
1717 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1718 return res;
1719}
1720
1721
1722
1723//_________________________________________________________________________
1724AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1725 //--------------------------------------------------------------------
1726 // Return pointer to a given cluster
1727 //--------------------------------------------------------------------
1728 Int_t sec=(index&0xff000000)>>24;
1729 Int_t row=(index&0x00ff0000)>>16;
1730 Int_t ncl=(index&0x0000ffff)>>00;
1731
1732 AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
9918f10a 1733 if (!clrow) return 0;
1c53abe2 1734 return (AliTPCclusterMI*)(*clrow)[ncl];
1735}
1736
1737//__________________________________________________________________________
1738void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1739 //--------------------------------------------------------------------
1740 //This function "cooks" a track label. If label<0, this track is fake.
1741 //--------------------------------------------------------------------
1742 Int_t noc=t->GetNumberOfClusters();
1743 Int_t *lb=new Int_t[noc];
1744 Int_t *mx=new Int_t[noc];
1745 AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1746
1747 Int_t i;
1748 for (i=0; i<noc; i++) {
1749 lb[i]=mx[i]=0;
1750 Int_t index=t->GetClusterIndex(i);
1751 clusters[i]=GetClusterMI(index);
1752 }
1753
1754 Int_t lab=123456789;
1755 for (i=0; i<noc; i++) {
1756 AliTPCclusterMI *c=clusters[i];
9918f10a 1757 if (!clusters[i]) continue;
1c53abe2 1758 lab=TMath::Abs(c->GetLabel(0));
1759 Int_t j;
1760 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1761 lb[j]=lab;
1762 (mx[j])++;
1763 }
1764
1765 Int_t max=0;
1766 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1767
1768 for (i=0; i<noc; i++) {
9918f10a 1769 AliTPCclusterMI *c=clusters[i];
1770 if (!clusters[i]) continue;
1c53abe2 1771 if (TMath::Abs(c->GetLabel(1)) == lab ||
1772 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1773 }
1774
1775 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1776
1777 else {
1778 Int_t tail=Int_t(0.10*noc);
1779 max=0;
1780 for (i=1; i<=tail; i++) {
1781 AliTPCclusterMI *c=clusters[noc-i];
9918f10a 1782 if (!clusters[i]) continue;
1c53abe2 1783 if (lab == TMath::Abs(c->GetLabel(0)) ||
1784 lab == TMath::Abs(c->GetLabel(1)) ||
1785 lab == TMath::Abs(c->GetLabel(2))) max++;
1786 }
1787 if (max < Int_t(0.5*tail)) lab=-lab;
1788 }
1789
1790 t->SetLabel(lab);
1791
1792 delete[] lb;
1793 delete[] mx;
1794 delete[] clusters;
1795}
1796
1797//_________________________________________________________________________
1798void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1799 //-----------------------------------------------------------------------
1800 // Setup inner sector
1801 //-----------------------------------------------------------------------
1802 if (f==0) {
1803 fAlpha=par->GetInnerAngle();
1804 fAlphaShift=par->GetInnerAngleShift();
1805 fPadPitchWidth=par->GetInnerPadPitchWidth();
1806 fPadPitchLength=par->GetInnerPadPitchLength();
1807 fN=par->GetNRowLow();
1808 fRow=new AliTPCRow[fN];
1809 for (Int_t i=0; i<fN; i++) {
1810 fRow[i].SetX(par->GetPadRowRadiiLow(i));
1811 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
1812 }
1813 } else {
1814 fAlpha=par->GetOuterAngle();
1815 fAlphaShift=par->GetOuterAngleShift();
1816 fPadPitchWidth = par->GetOuterPadPitchWidth();
1817 fPadPitchLength = par->GetOuter1PadPitchLength();
1818 f1PadPitchLength = par->GetOuter1PadPitchLength();
1819 f2PadPitchLength = par->GetOuter2PadPitchLength();
1820
1821 fN=par->GetNRowUp();
1822 fRow=new AliTPCRow[fN];
1823 for (Int_t i=0; i<fN; i++) {
1824 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1825 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
1826 }
1827 }
1828}
1829
1830
1831AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
1832 //
1833 if (fClusterTracks) delete [] fClusterTracks;
1834 fClusterTracks = 0;
1835}
1836
1837void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
1838 //create cluster tracks
1839 if (fN>0)
1840 fClusterTracks = new AliTPCclusterTracks[fN];
1841}
1842
1843void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
1844 if (fClusterTracks) delete[] fClusterTracks;
1845 fClusterTracks =0;
1846}
1847
1848
1849
1850void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
1851 //
1852 //
1853 // update information of the cluster tracks - if track is nearer then other tracks to the
1854 // given track
1855 const AliTPCclusterMI * cl = (*this)[clindex];
1856 AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
1857 // find the distance of the cluster to the track
1858 Float_t dy2 = (cl->GetY()- seed->GetY());
1859 dy2*=dy2;
1860 Float_t dz2 = (cl->GetZ()- seed->GetZ());
1861 dz2*=dz2;
1862 //
1863 Float_t distance = TMath::Sqrt(dy2+dz2);
1864 if (distance> 3)
1865 return; // MI - to be changed - AliTPCtrackerParam
1866
1867 if ( distance < cltracks->fDistance[0]){
1868 cltracks->fDistance[2] =cltracks->fDistance[1];
1869 cltracks->fDistance[1] =cltracks->fDistance[0];
1870 cltracks->fDistance[0] =distance;
1871 cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1872 cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
1873 cltracks->fTrackIndex[0] =trindex;
1874 }
1875 else
1876 if ( distance < cltracks->fDistance[1]){
1877 cltracks->fDistance[2] =cltracks->fDistance[1];
1878 cltracks->fDistance[1] =distance;
1879 cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1880 cltracks->fTrackIndex[1] =trindex;
1881 } else
1882 if (distance < cltracks->fDistance[2]){
1883 cltracks->fDistance[2] =distance;
1884 cltracks->fTrackIndex[2] =trindex;
1885 }
1886}
1887
1888
1889//_________________________________________________________________________
1890void
1891AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
1892 //-----------------------------------------------------------------------
1893 // Insert a cluster into this pad row in accordence with its y-coordinate
1894 //-----------------------------------------------------------------------
1895 if (fN==kMaxClusterPerRow) {
1896 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1897 }
1898 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
1899 Int_t i=Find(c->GetZ());
1900 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
1901 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
1902 fIndex[i]=index; fClusters[i]=c; fN++;
1903}
1904
1905//___________________________________________________________________
1906Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
1907 //-----------------------------------------------------------------------
1908 // Return the index of the nearest cluster
1909 //-----------------------------------------------------------------------
1910 if (fN==0) return 0;
1911 if (z <= fClusters[0]->GetZ()) return 0;
1912 if (z > fClusters[fN-1]->GetZ()) return fN;
1913 Int_t b=0, e=fN-1, m=(b+e)/2;
1914 for (; b<e; m=(b+e)/2) {
1915 if (z > fClusters[m]->GetZ()) b=m+1;
1916 else e=m;
1917 }
1918 return m;
1919}
1920
1921
1922
1627d1c4 1923//___________________________________________________________________
1924AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
1925 //-----------------------------------------------------------------------
1926 // Return the index of the nearest cluster in z y
1927 //-----------------------------------------------------------------------
1928 Float_t maxdistance = roady*roady + roadz*roadz;
1929
1930 AliTPCclusterMI *cl =0;
1931 for (Int_t i=Find(z-roadz); i<fN; i++) {
1932 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
1933 if (c->GetZ() > z+roadz) break;
1934 if ( (c->GetY()-y) > roady ) continue;
1935 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
1936 if (maxdistance>distance) {
1937 maxdistance = distance;
1938 cl=c;
1939 }
1940 }
1941 return cl;
1942}
1943
1944
1945
1946
1947
1c53abe2 1948AliTPCseed::AliTPCseed():AliTPCtrack(){
1949 //
1950 fRow=0;
1951 fRemoval =0;
1952 memset(fClusterIndex,0,sizeof(Int_t)*200);
1953 fPoints = 0;
1954 fEPoints = 0;
1955 fNFoundable =0;
1956 fNShared =0;
1957 fTrackPoints =0;
1958 fRemoval = 0;
1959}
1960
1961AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
1962 fPoints = 0;
1963 fEPoints = 0;
1964 fNShared =0;
1965 fTrackPoints =0;
1966 fRemoval =0;
1967}
1968
1969AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
1970 fRow=0;
1971 memset(fClusterIndex,0,sizeof(Int_t)*200);
1972 fPoints = 0;
1973 fEPoints = 0;
1974 fNFoundable =0;
1975 fNShared =0;
1976 fTrackPoints =0;
1977 fRemoval =0;
1978}
1979
1980AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
1981 Double_t xr, Double_t alpha):
1982 AliTPCtrack(index, xx, cc, xr, alpha) {
1983 //
1984 //
1985 fRow =0;
1986 memset(fClusterIndex,0,sizeof(Int_t)*200);
1987 fPoints = 0;
1988 fEPoints = 0;
1989 fNFoundable =0;
1990 fNShared = 0;
1991 fTrackPoints =0;
1992 fRemoval =0;
1993}
1994
1995AliTPCseed::~AliTPCseed(){
1996 if (fPoints) delete fPoints;
1997 fPoints =0;
1998 fEPoints = 0;
1999 if (fTrackPoints){
2000 for (Int_t i=0;i<8;i++){
2001 delete [] fTrackPoints[i];
2002 }
2003 delete fTrackPoints;
2004 fTrackPoints =0;
2005 }
2006
2007}
2008
2009AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
2010{
2011 //
2012 //
2013 if (!fTrackPoints) {
2014 fTrackPoints = new AliTPCTrackPoint*[8];
2015 for ( Int_t i=0;i<8;i++)
2016 fTrackPoints[i]=0;
2017 }
2018 Int_t index1 = i/20;
2019 if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
2020 return &(fTrackPoints[index1][i%20]);
2021}
2022
2023void AliTPCseed::RebuildSeed()
2024{
2025 //
2026 // rebuild seed to be ready for storing
2027 fPoints = new TClonesArray("AliTPCTrackPoint",160);
2028 fPoints->ExpandCreateFast(160);
2029 fEPoints = new TClonesArray("AliTPCExactPoint",1);
2030 for (Int_t i=0;i<160;i++){
2031 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
2032 *trpoint = *(GetTrackPoint(i));
2033 }
2034
2035}
2036
2037//_____________________________________________________________________________
2038void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
2039 //-----------------------------------------------------------------
2040 // This funtion calculates dE/dX within the "low" and "up" cuts.
2041 //-----------------------------------------------------------------
2042
2043 Float_t amp[200];
2044 Float_t angular[200];
2045 Float_t weight[200];
2046 Int_t index[200];
2047 //Int_t nc = 0;
2048 // TClonesArray & arr = *fPoints;
2049 Float_t meanlog = 100.;
2050
2051 Float_t mean[4] = {0,0,0,0};
2052 Float_t sigma[4] = {1000,1000,1000,1000};
2053 Int_t nc[4] = {0,0,0,0};
2054 Float_t norm[4] = {1000,1000,1000,1000};
2055 //
2056 //
2057 fNShared =0;
2058
2059 for (Int_t of =0; of<4; of++){
2060 for (Int_t i=of;i<160;i+=4)
2061 {
2062 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
2063 AliTPCTrackPoint * point = GetTrackPoint(i);
2064 if (point==0) continue;
2065 if (point->fIsShared){
2066 fNShared++;
2067 continue;
2068 }
2069 if (point->GetCPoint().GetMax()<5) continue;
2070 Float_t angley = point->GetTPoint().GetAngleY();
2071 Float_t anglez = point->GetTPoint().GetAngleZ();
2072 Int_t type = point->GetCPoint().GetType();
2073 Float_t rsigmay = point->GetCPoint().GetSigmaY();
2074 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
2075 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
2076
2077 Float_t ampc = 0; // normalization to the number of electrons
2078 if (i>64){
2079 ampc = 1.*point->GetCPoint().GetMax();
2080 //ampc = 1.*point->GetCPoint().GetQ();
2081 // AliTPCClusterPoint & p = point->GetCPoint();
2082 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
2083 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2084 //Float_t dz =
2085 // TMath::Abs( Int_t(iz) - iz + 0.5);
2086 //ampc *= 1.15*(1-0.3*dy);
2087 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 2088 // Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
2089 //ampc *=zfactor;
1c53abe2 2090 }
2091 else{
2092 ampc = 1.0*point->GetCPoint().GetMax();
2093 //ampc = 1.0*point->GetCPoint().GetQ();
2094 //AliTPCClusterPoint & p = point->GetCPoint();
2095 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
2096 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2097 //Float_t dz =
2098 // TMath::Abs( Int_t(iz) - iz + 0.5);
2099
2100 //ampc *= 1.15*(1-0.3*dy);
2101 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 2102 // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
2103 //ampc *=zfactor;
1c53abe2 2104
2105 }
2106 ampc *= 2.0; // put mean value to channel 50
2107 //ampc *= 0.58; // put mean value to channel 50
2108 Float_t w = 1.;
2109 // if (type>0) w = 1./(type/2.-0.5);
1627d1c4 2110 Float_t z = TMath::Abs(point->GetCPoint().GetZ());
1c53abe2 2111 if (i<64) {
2112 ampc /= 0.6;
1627d1c4 2113 //ampc /= (1+0.0008*z);
2114 } else
2115 if (i>128){
2116 ampc /=1.5;
2117 //ampc /= (1+0.0008*z);
2118 }else{
2119 //ampc /= (1+0.0008*z);
2120 }
2121
1c53abe2 2122 if (type<0) { //amp at the border - lower weight
2123 // w*= 2.;
1627d1c4 2124
1c53abe2 2125 continue;
2126 }
2127 if (rsigma>1.5) ampc/=1.3; // if big backround
2128 amp[nc[of]] = ampc;
2129 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
2130 weight[nc[of]] = w;
2131 nc[of]++;
2132 }
2133
2134 TMath::Sort(nc[of],amp,index,kFALSE);
2135 Float_t sumamp=0;
2136 Float_t sumamp2=0;
2137 Float_t sumw=0;
2138 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
2139 meanlog = 200;
2140 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
2141 Float_t ampl = amp[index[i]]/angular[index[i]];
2142 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
2143 //
2144 sumw += weight[index[i]];
2145 sumamp += weight[index[i]]*ampl;
2146 sumamp2 += weight[index[i]]*ampl*ampl;
2147 norm[of] += angular[index[i]]*weight[index[i]];
2148 }
2149 if (sumw<1){
2150 SetdEdx(0);
2151 }
2152 else {
2153 norm[of] /= sumw;
2154 mean[of] = sumamp/sumw;
2155 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
2156 if (sigma[of]>0.1)
2157 sigma[of] = TMath::Sqrt(sigma[of]);
2158 else
2159 sigma[of] = 1000;
2160
2161 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
2162 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
2163 //mean *=(1-0.1*(norm-1.));
2164 }
2165 }
2166
2167 Float_t dedx =0;
2168 fSdEdx =0;
2169 fMAngular =0;
2170 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
2171 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
2172
2173
2174 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
2175 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
2176
2177 Int_t norm2 = 0;
2178 Int_t norm3 = 0;
2179 for (Int_t i =0;i<4;i++){
2180 if (nc[i]>2&&nc[i]<1000){
2181 dedx += mean[i] *nc[i];
2182 fSdEdx += sigma[i]*(nc[i]-2);
2183 fMAngular += norm[i] *nc[i];
2184 norm2 += nc[i];
2185 norm3 += nc[i]-2;
2186 }
2187 fDEDX[i] = mean[i];
2188 fSDEDX[i] = sigma[i];
2189 fNCDEDX[i]= nc[i];
2190 }
2191
2192 if (norm3>0){
2193 dedx /=norm2;
2194 fSdEdx /=norm3;
2195 fMAngular/=norm2;
2196 }
2197 else{
2198 SetdEdx(0);
2199 return;
2200 }
2201 // Float_t dedx1 =dedx;
2202
2203 dedx =0;
2204 for (Int_t i =0;i<4;i++){
2205 if (nc[i]>2&&nc[i]<1000){
2206 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
2207 dedx += mean[i] *nc[i];
2208 }
2209 fDEDX[i] = mean[i];
2210 }
2211 dedx /= norm2;
2212
2213
2214
2215 SetdEdx(dedx);
2216
2217 //mi deDX
2218
2219
2220
2221 //Very rough PID
2222 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
2223
2224 if (p<0.6) {
2225 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2226 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
2227 SetMass(0.93827); return;
2228 }
2229
2230 if (p<1.2) {
2231 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2232 SetMass(0.93827); return;
2233 }
2234
2235 SetMass(0.13957); return;
2236
2237}
2238
2239