]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCtrackerMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug fix (bug in compiler)
[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
91162307 16/*
17$Log$
732b9e78 18Revision 1.13 2003/09/29 11:28:19 kowal2
19completly rewritten
20
91162307 21Revision 1.9.4.3 2003/06/23 14:47:10 hristov
22Minor fix
23
24Revision 1.9.4.2 2003/06/23 10:06:13 hristov
25Updated information about the overlapping clusters (M.Ivanov)
26
27Revision 1.9.4.1 2003/06/19 06:59:58 hristov
28Updated version of parallel tracking (M.Ivanov)
29
30Revision 1.9 2003/03/19 17:14:11 hristov
31Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
32
33Revision 1.8 2003/03/05 11:16:15 kowal2
34Logs added
35
36*/
37
38
39
40
41
42
1c53abe2 43
44/*
45 AliTPC parallel tracker -
46 How to use? -
47 run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
48 run AliTPCFindTracksMI.C macro - to find tracks
49 tracks are written to AliTPCtracks.root file
50 for comparison also seeds are written to the same file - to special branch
51*/
52
53//-------------------------------------------------------
54// Implementation of the TPC tracker
55//
56// Origin: Marian Ivanov Marian.Ivanov@cern.ch
57//
58//-------------------------------------------------------
1c53abe2 59#include <TObjArray.h>
60#include <TFile.h>
61#include <TTree.h>
cc5e9db0 62#include "Riostream.h"
1c53abe2 63
64#include "AliTPCtrackerMI.h"
65#include "AliTPCclusterMI.h"
66#include "AliTPCParam.h"
67#include "AliTPCClustersRow.h"
68#include "AliComplexCluster.h"
1627d1c4 69#include "AliTPCpolyTrack.h"
1c53abe2 70#include "TStopwatch.h"
91162307 71#include "AliESD.h"
72#include "AliHelix.h"
73#include "TGraph.h"
74//
75#include "AliRunLoader.h"
1c53abe2 76
77
c9427e08 78
91162307 79ClassImp(AliTPCseed)
80ClassImp(AliTPCtrackerMI)
c9427e08 81
82
91162307 83class TPCFastMath {
84public:
85 TPCFastMath();
86 static Double_t FastAsin(Double_t x);
87 private:
88 static Double_t fgFastAsin[20000];
89};
c9427e08 90
91162307 91Double_t TPCFastMath::fgFastAsin[20000];
c9427e08 92
91162307 93TPCFastMath::TPCFastMath(){
94 for (Int_t i=0;i<10000;i++){
95 fgFastAsin[2*i] = TMath::ASin(i/10000.);
96 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
97 }
c9427e08 98}
99
91162307 100Double_t TPCFastMath::FastAsin(Double_t x){
101 if (x>0){
102 Int_t index = int(x*10000);
103 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
104 }
105 x*=-1;
106 Int_t index = int(x*10000);
107 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
1c53abe2 108}
109
91162307 110TPCFastMath gTPCMath;
1c53abe2 111
112
113
91162307 114Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
1c53abe2 115
91162307 116 AliTPCclusterMI* c =track->fCurrentCluster;
117 if (accept>0) track->fCurrentClusterIndex1 |=0x8000; //sign not accepted clusters
c9427e08 118
91162307 119 UInt_t i = track->fCurrentClusterIndex1;
1c53abe2 120
121 Int_t sec=(i&0xff000000)>>24;
91162307 122 //Int_t row = (i&0x00ff0000)>>16;
1c53abe2 123 track->fRow=(i&0x00ff0000)>>16;
124 track->fSector = sec;
125 // Int_t index = i&0xFFFF;
126 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
91162307 127 track->SetClusterIndex2(track->fRow, i);
128
129 //track->fFirstPoint = row;
130 //if ( track->fLastPoint<row) track->fLastPoint =row;
131 // if (track->fRow<0 || track->fRow>160) {
132 // printf("problem\n");
133 //}
134 if (track->fFirstPoint>track->fRow)
135 track->fFirstPoint = track->fRow;
136 if (track->fLastPoint<track->fRow)
137 track->fLastPoint = track->fRow;
138
139
140 track->fClusterPointer[track->fRow] = c;
1c53abe2 141 //
142
1c53abe2 143 Float_t angle2 = track->GetSnp()*track->GetSnp();
144 angle2 = TMath::Sqrt(angle2/(1-angle2));
145 //
146 //SET NEW Track Point
147 //
91162307 148 // if (debug)
149 {
150 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow));
1c53abe2 151 //
91162307 152 point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
153 point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
154 point.SetErrY(sqrt(track->fErrorY2));
155 point.SetErrZ(sqrt(track->fErrorZ2));
1c53abe2 156 //
91162307 157 point.SetX(track->GetX());
158 point.SetY(track->GetY());
159 point.SetZ(track->GetZ());
160 point.SetAngleY(angle2);
161 point.SetAngleZ(track->GetTgl());
162 if (point.fIsShared){
163 track->fErrorY2 *= 4;
164 track->fErrorZ2 *= 4;
165 }
166 }
167
168 Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
169 //
170 track->fErrorY2 *= 1.3;
171 track->fErrorY2 += 0.01;
172 track->fErrorZ2 *= 1.3;
173 track->fErrorZ2 += 0.005;
174 //}
175 if (accept>0) return 0;
176 if (track->GetNumberOfClusters()%20==0){
177 // if (track->fHelixIn){
178 // TClonesArray & larr = *(track->fHelixIn);
179 // Int_t ihelix = larr.GetEntriesFast();
180 // new(larr[ihelix]) AliHelix(*track) ;
181 //}
1c53abe2 182 }
91162307 183 track->fNoCluster =0;
184 return track->Update(c,chi2,i);
185}
186
187
188
189Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
190 Float_t cory, Float_t corz)
191{
1c53abe2 192 //
91162307 193 // decide according desired precision to accept given
194 // cluster for tracking
195 Double_t sy2=ErrY2(seed,cluster)*cory;
196 Double_t sz2=ErrZ2(seed,cluster)*corz;
197 //sy2=ErrY2(seed,cluster)*cory;
198 //sz2=ErrZ2(seed,cluster)*cory;
1c53abe2 199
91162307 200 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
201 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
202
203 Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
204 (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
205 Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
206 (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
207
208 Double_t rdistance2 = rdistancey2+rdistancez2;
209 //Int_t accept =0;
1c53abe2 210
91162307 211 if (rdistance2>16) return 3;
212
213
214 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
215 return 2; //suspisiouce - will be changed
216
217 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
218 // strict cut on overlaped cluster
219 return 2; //suspisiouce - will be changed
220
221 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
222 && cluster->GetType()<0){
223 seed->fNFoundable--;
224 return 2;
1c53abe2 225 }
91162307 226 return 0;
227}
228
229
1c53abe2 230
1c53abe2 231
1c53abe2 232//_____________________________________________________________________________
f8aae377 233AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
1c53abe2 234AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
235{
236 //---------------------------------------------------------------------
237 // The main TPC tracker constructor
238 //---------------------------------------------------------------------
239 fInnerSec=new AliTPCSector[fkNIS];
240 fOuterSec=new AliTPCSector[fkNOS];
91162307 241
1c53abe2 242 Int_t i;
243 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
244 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
245
246 fN=0; fSectors=0;
247
1c53abe2 248 fSeeds=0;
249 fNtracks = 0;
91162307 250 fParam = par;
251 Int_t nrowlow = par->GetNRowLow();
252 Int_t nrowup = par->GetNRowUp();
253
254
255 for (Int_t i=0;i<nrowlow;i++){
256 fXRow[i] = par->GetPadRowRadiiLow(i);
257 fPadLength[i]= par->GetPadPitchLength(0,i);
258 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
259 }
260
261
262 for (Int_t i=0;i<nrowup;i++){
263 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
264 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
265 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
266 }
267 fSeeds=0;
268 //
269 fInput = 0;
270 fOutput = 0;
271 fSeedTree = 0;
272 fTreeDebug =0;
273 fNewIO =0;
274 fDebug =0;
275 fEvent =0;
1c53abe2 276}
277
278//_____________________________________________________________________________
279AliTPCtrackerMI::~AliTPCtrackerMI() {
280 //------------------------------------------------------------------
281 // TPC tracker destructor
282 //------------------------------------------------------------------
283 delete[] fInnerSec;
284 delete[] fOuterSec;
285 if (fSeeds) {
286 fSeeds->Delete();
287 delete fSeeds;
288 }
289}
290
91162307 291void AliTPCtrackerMI::SetIO()
292{
1c53abe2 293 //
91162307 294 fNewIO = kTRUE;
295 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::fgkDefaultEventFolderName);
296 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::fgkDefaultEventFolderName);
297 AliTPCtrack *iotrack= new AliTPCtrack;
298 // iotrack->fHelixIn = new TClonesArray("AliHelix");
299 //iotrack->fHelixOut = new TClonesArray("AliHelix");
300 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
301 delete iotrack;
302}
1c53abe2 303
91162307 304void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
305{
1c53abe2 306
91162307 307 // set input
308 fNewIO = kFALSE;
309 fInput = 0;
310 fOutput = 0;
311 fSeedTree = 0;
312 fTreeDebug =0;
313 fInput = input;
314 if (input==0){
315 return;
316 }
317 //set output
318 fOutput = output;
319 if (output){
320 AliTPCtrack *iotrack= new AliTPCtrack;
321 // iotrack->fHelixIn = new TClonesArray("AliHelix");
322 //iotrack->fHelixOut = new TClonesArray("AliHelix");
323 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
324 delete iotrack;
325 }
326 if (output && (fDebug&2)){
327 //write the full seed information if specified in debug mode
328 //
329 fSeedTree = new TTree("Seeds","Seeds");
330 AliTPCseed * vseed = new AliTPCseed;
331 //
332 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
333 arrtr->ExpandCreateFast(160);
334 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
335 //
336 vseed->fPoints = arrtr;
337 vseed->fEPoints = arre;
338 // vseed->fClusterPoints = arrcl;
339 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
340 delete arrtr;
341 delete arre;
342 fTreeDebug = new TTree("trackDebug","trackDebug");
343 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
344 fTreeDebug->Branch("debug",&arrd,32000,99);
1c53abe2 345 }
1c53abe2 346
1c53abe2 347
91162307 348 //set ESD event
349 fEvent = event;
350}
1c53abe2 351
91162307 352void AliTPCtrackerMI::WriteTracks()
353{
354 //
355 // write tracks to the given output tree -
356 // output specified with SetIO routine
357 if (!fSeeds) return;
358 if (fEvent){
359 // write tracks to the event
360 // store index of the track
361 Int_t nseed=fSeeds->GetEntriesFast();
362 for (Int_t i=0; i<nseed; i++) {
363 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
364 if (!pt) continue;
365 AliESDtrack iotrack;
366 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
367 //iotrack.SetTPCindex(i);
368 fEvent->AddTrack(&iotrack);
369 }
1c53abe2 370 }
371
1c53abe2 372
91162307 373 if (fOutput){
374 AliTPCtrack *iotrack= 0;
375 Int_t nseed=fSeeds->GetEntriesFast();
376 for (Int_t i=0; i<nseed; i++) {
377 iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
378 if (iotrack) break;
379 }
380
381 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
382 TBranch * br = fOutput->GetBranch("tracks");
383 br->SetAddress(&iotrack);
384 //
385 for (Int_t i=0; i<nseed; i++) {
386 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
387 if (!pt) continue;
388 iotrack = pt;
389 pt->fLab2 =i;
390 // br->SetAddress(&iotrack);
391 fOutput->Fill();
392 iotrack =0;
393 }
394 }
395 // delete iotrack;
396 //
397 if (fSeedTree){
398 //write the full seed information if specified in debug mode
399
400 AliTPCseed * vseed = new AliTPCseed;
401 //
402 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
403 arrtr->ExpandCreateFast(160);
404 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
405 //arrcl->ExpandCreateFast(160);
406 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
407 //
408 vseed->fPoints = arrtr;
409 vseed->fEPoints = arre;
410 // vseed->fClusterPoints = arrcl;
411 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
412 TBranch * brseed = fSeedTree->GetBranch("seeds");
413
414 Int_t nseed=fSeeds->GetEntriesFast();
415
416 for (Int_t i=0; i<nseed; i++) {
417 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
418 if (!pt) continue;
419 pt->fPoints = arrtr;
420 // pt->fClusterPoints = arrcl;
421 pt->fEPoints = arre;
422 pt->RebuildSeed();
423 vseed = pt;
424 brseed->SetAddress(&vseed);
425 fSeedTree->Fill();
426 pt->fPoints = 0;
427 pt->fEPoints = 0;
428 // pt->fClusterPoints = 0;
429 }
430 fSeedTree->Write();
431 if (fTreeDebug) fTreeDebug->Write();
432 }
1c53abe2 433
91162307 434}
1c53abe2 435
1c53abe2 436
1c53abe2 437
438
91162307 439Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
440 //
441 //
442 //seed->SetErrorY2(0.1);
443 //return 0.1;
444 //calculate look-up table at the beginning
445 static Bool_t ginit = kFALSE;
446 static Float_t gnoise1,gnoise2,gnoise3;
447 static Float_t ggg1[10000];
448 static Float_t ggg2[10000];
449 static Float_t ggg3[10000];
450 static Float_t glandau1[10000];
451 static Float_t glandau2[10000];
452 static Float_t glandau3[10000];
453 //
454 static Float_t gcor01[500];
455 static Float_t gcor02[500];
456 static Float_t gcorp[500];
457 //
1c53abe2 458
91162307 459 //
460 if (ginit==kFALSE){
461 for (Int_t i=1;i<500;i++){
462 Float_t rsigma = float(i)/100.;
463 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
464 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
465 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
466 }
1c53abe2 467
91162307 468 //
469 for (Int_t i=3;i<10000;i++){
470 //
471 //
472 // inner sector
473 Float_t amp = float(i);
474 Float_t padlength =0.75;
475 gnoise1 = 0.0004/padlength;
476 Float_t nel = 0.268*amp;
477 Float_t nprim = 0.155*amp;
478 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
479 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
480 if (glandau1[i]>1) glandau1[i]=1;
481 glandau1[i]*=padlength*padlength/12.;
482 //
483 // outer short
484 padlength =1.;
485 gnoise2 = 0.0004/padlength;
486 nel = 0.3*amp;
487 nprim = 0.133*amp;
488 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
489 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
490 if (glandau2[i]>1) glandau2[i]=1;
491 glandau2[i]*=padlength*padlength/12.;
492 //
493 //
494 // outer long
495 padlength =1.5;
496 gnoise3 = 0.0004/padlength;
497 nel = 0.3*amp;
498 nprim = 0.133*amp;
499 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
500 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
501 if (glandau3[i]>1) glandau3[i]=1;
502 glandau3[i]*=padlength*padlength/12.;
503 //
504 }
505 ginit = kTRUE;
506 }
1c53abe2 507 //
508 //
91162307 509 //
510 Int_t amp = int(TMath::Abs(cl->GetQ()));
511 if (amp>9999) {
512 seed->SetErrorY2(1.);
513 return 1.;
514 }
1c53abe2 515 Float_t snoise2;
c9427e08 516 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
91162307 517 Int_t ctype = cl->GetType();
518 Float_t padlength= GetPadPitchLength(seed->fRow);
519 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
520 angle2 = angle2/(1-angle2);
521 //
522 //cluster "quality"
523 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
524 Float_t res;
1c53abe2 525 //
1c53abe2 526 if (fSectors==fInnerSec){
91162307 527 snoise2 = gnoise1;
528 res = ggg1[amp]*z+glandau1[amp]*angle2;
529 if (ctype==0) res *= gcor01[rsigmay];
530 if ((ctype>0)){
531 res+=0.002;
532 res*= gcorp[rsigmay];
533 }
1c53abe2 534 }
535 else {
91162307 536 if (padlength<1.1){
537 snoise2 = gnoise2;
538 res = ggg2[amp]*z+glandau2[amp]*angle2;
539 if (ctype==0) res *= gcor02[rsigmay];
540 if ((ctype>0)){
541 res+=0.002;
542 res*= gcorp[rsigmay];
543 }
544 }
545 else{
546 snoise2 = gnoise3;
547 res = ggg3[amp]*z+glandau3[amp]*angle2;
548 if (ctype==0) res *= gcor02[rsigmay];
549 if ((ctype>0)){
550 res+=0.002;
551 res*= gcorp[rsigmay];
552 }
553 }
554 }
1c53abe2 555
91162307 556 if (ctype<0){
557 res+=0.005;
558 res*=2.4; // overestimate error 2 times
559 }
560 res+= snoise2;
561
1c53abe2 562 if (res<2*snoise2)
91162307 563 res = 2*snoise2;
564
565 seed->SetErrorY2(res);
1c53abe2 566 return res;
1c53abe2 567
568
91162307 569}
c9427e08 570
571
572
91162307 573Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
574 //
575 //
576 //seed->SetErrorY2(0.1);
577 //return 0.1;
578 //calculate look-up table at the beginning
579 static Bool_t ginit = kFALSE;
580 static Float_t gnoise1,gnoise2,gnoise3;
581 static Float_t ggg1[10000];
582 static Float_t ggg2[10000];
583 static Float_t ggg3[10000];
584 static Float_t glandau1[10000];
585 static Float_t glandau2[10000];
586 static Float_t glandau3[10000];
587 //
588 static Float_t gcor01[1000];
589 static Float_t gcor02[1000];
590 static Float_t gcorp[1000];
1627d1c4 591 //
1627d1c4 592
91162307 593 //
594 if (ginit==kFALSE){
595 for (Int_t i=1;i<1000;i++){
596 Float_t rsigma = float(i)/100.;
597 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
598 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
599 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
600 }
1c53abe2 601
91162307 602 //
603 for (Int_t i=3;i<10000;i++){
604 //
605 //
606 // inner sector
607 Float_t amp = float(i);
608 Float_t padlength =0.75;
609 gnoise1 = 0.0004/padlength;
610 Float_t nel = 0.268*amp;
611 Float_t nprim = 0.155*amp;
612 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
613 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
614 if (glandau1[i]>1) glandau1[i]=1;
615 glandau1[i]*=padlength*padlength/12.;
616 //
617 // outer short
618 padlength =1.;
619 gnoise2 = 0.0004/padlength;
620 nel = 0.3*amp;
621 nprim = 0.133*amp;
622 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
623 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
624 if (glandau2[i]>1) glandau2[i]=1;
625 glandau2[i]*=padlength*padlength/12.;
626 //
627 //
628 // outer long
629 padlength =1.5;
630 gnoise3 = 0.0004/padlength;
631 nel = 0.3*amp;
632 nprim = 0.133*amp;
633 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
634 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
635 if (glandau3[i]>1) glandau3[i]=1;
636 glandau3[i]*=padlength*padlength/12.;
637 //
638 }
639 ginit = kTRUE;
640 }
641 //
642 //
643 //
644 Int_t amp = int(TMath::Abs(cl->GetQ()));
645 if (amp>9999) {
646 seed->SetErrorY2(1.);
647 return 1.;
648 }
649 Float_t snoise2;
650 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
651 Int_t ctype = cl->GetType();
652 Float_t padlength= GetPadPitchLength(seed->fRow);
653 //
654 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
655 // if (angle2<0.6) angle2 = 0.6;
656 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
657 //
658 //cluster "quality"
659 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
660 Float_t res;
661 //
662 if (fSectors==fInnerSec){
663 snoise2 = gnoise1;
664 res = ggg1[amp]*z+glandau1[amp]*angle2;
665 if (ctype==0) res *= gcor01[rsigmaz];
666 if ((ctype>0)){
667 res+=0.002;
668 res*= gcorp[rsigmaz];
669 }
670 }
671 else {
672 if (padlength<1.1){
673 snoise2 = gnoise2;
674 res = ggg2[amp]*z+glandau2[amp]*angle2;
675 if (ctype==0) res *= gcor02[rsigmaz];
676 if ((ctype>0)){
677 res+=0.002;
678 res*= gcorp[rsigmaz];
679 }
680 }
681 else{
682 snoise2 = gnoise3;
683 res = ggg3[amp]*z+glandau3[amp]*angle2;
684 if (ctype==0) res *= gcor02[rsigmaz];
685 if ((ctype>0)){
686 res+=0.002;
687 res*= gcorp[rsigmaz];
688 }
689 }
690 }
691
692 if (ctype<0){
693 res+=0.002;
694 res*=1.3;
695 }
696 if ((ctype<0) &&amp<70){
697 res+=0.002;
698 res*=1.3;
699 }
700 res += snoise2;
701 if (res<2*snoise2)
702 res = 2*snoise2;
703 if (res>3) res =3;
704 seed->SetErrorZ2(res);
705 return res;
706}
707
708
709
710/*
711Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
712 //
713 //
714 //seed->SetErrorZ2(0.1);
715 //return 0.1;
716
717 Float_t snoise2;
718 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
719 //
720 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
721 Int_t ctype = cl->GetType();
722 Float_t amp = TMath::Abs(cl->GetQ());
723
724 Float_t nel;
725 Float_t nprim;
726 //
727 Float_t landau=2 ; //landau fluctuation part
728 Float_t gg=2; // gg fluctuation part
729 Float_t padlength= GetPadPitchLength(seed->GetX());
730
731 if (fSectors==fInnerSec){
732 snoise2 = 0.0004/padlength;
733 nel = 0.268*amp;
734 nprim = 0.155*amp;
735 gg = (2+0.001*nel/(padlength*padlength))/nel;
736 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
737 if (landau>1) landau=1;
738 }
739 else {
740 snoise2 = 0.0004/padlength;
741 nel = 0.3*amp;
742 nprim = 0.133*amp;
743 gg = (2+0.0008*nel/(padlength*padlength))/nel;
744 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
745 if (landau>1) landau=1;
746 }
747 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
748
749 //
750 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
751 angle2 = TMath::Sqrt((1-angle2));
752 if (angle2<0.6) angle2 = 0.6;
753 //angle2 = 1;
754
755 Float_t angle = seed->GetTgl()/angle2;
756 Float_t angular = landau*angle*angle*padlength*padlength/12.;
757 Float_t res = sdiff + angular;
758
759
760 if ((ctype==0) && (fSectors ==fOuterSec))
761 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
762
763 if ((ctype==0) && (fSectors ==fInnerSec))
764 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
765
766 if ((ctype>0)){
767 res+=0.005;
768 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
769 }
770 if (ctype<0){
771 res+=0.002;
772 res*=1.3;
773 }
774 if ((ctype<0) &&amp<70){
775 res+=0.002;
776 res*=1.3;
777 }
778 res += snoise2;
779 if (res<2*snoise2)
780 res = 2*snoise2;
781
782 seed->SetErrorZ2(res);
783 return res;
784}
785*/
786
787
788
789void AliTPCseed::Reset(Bool_t all)
790{
791 //
792 //
793 SetNumberOfClusters(0);
794 fNFoundable = 0;
795 SetChi2(0);
796 ResetCovariance();
797 /*
798 if (fTrackPoints){
799 for (Int_t i=0;i<8;i++){
800 delete [] fTrackPoints[i];
801 }
802 delete fTrackPoints;
803 fTrackPoints =0;
804 }
805 */
806
807 if (all){
808 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
809 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
810 }
811
812}
813
814
815void AliTPCseed::Modify(Double_t factor)
816{
817
818 //------------------------------------------------------------------
819 //This function makes a track forget its history :)
820 //------------------------------------------------------------------
821 if (factor<=0) {
822 ResetCovariance();
823 return;
824 }
825 fC00*=factor;
826 fC10*=factor; fC11*=factor;
827 fC20*=factor; fC21*=factor; fC22*=factor;
828 fC30*=factor; fC31*=factor; fC32*=factor; fC33*=factor;
829 fC40*=factor; fC41*=factor; fC42*=factor; fC43*=factor; fC44*=factor;
830 SetNumberOfClusters(0);
831 fNFoundable =0;
832 SetChi2(0);
833 fRemoval = 0;
834 fCurrentSigmaY2 = 0.000005;
835 fCurrentSigmaZ2 = 0.000005;
836 fNoCluster = 0;
837 //fFirstPoint = 160;
838 //fLastPoint = 0;
839}
840
841
842
843
844Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
845{
846 //-----------------------------------------------------------------
847 // This function find proloncation of a track to a reference plane x=xk.
848 // doesn't change internal state of the track
849 //-----------------------------------------------------------------
850
1c53abe2 851 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
91162307 852
853 if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
854 return 0;
855 }
856
1c53abe2 857 // Double_t y1=fP0, z1=fP1;
858 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
859 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
860
861 y = fP0;
862 z = fP1;
91162307 863 //y += dx*(c1+c2)/(r1+r2);
864 //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
865
866 Double_t dy = dx*(c1+c2)/(r1+r2);
867 Double_t dz = 0;
868 //
869 Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
870 /*
871 if (TMath::Abs(delta)>0.0001){
872 dz = fP3*TMath::ASin(delta)/fP4;
873 }else{
874 dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
875 }
876 */
877 dz = fP3*TPCFastMath::FastAsin(delta)/fP4;
878 //
879 y+=dy;
880 z+=dz;
881
882
883 return 1;
1c53abe2 884}
885
886
887//_____________________________________________________________________________
888Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
889{
890 //-----------------------------------------------------------------
891 // This function calculates a predicted chi2 increment.
892 //-----------------------------------------------------------------
893 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
894 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
895 r00+=fC00; r01+=fC10; r11+=fC11;
896
897 Double_t det=r00*r11 - r01*r01;
898 if (TMath::Abs(det) < 1.e-10) {
899 Int_t n=GetNumberOfClusters();
900 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
901 return 1e10;
902 }
903 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
904
905 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
906
907 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
908}
909
910
911//_________________________________________________________________________________________
912
913
914Int_t AliTPCseed::Compare(const TObject *o) const {
915 //-----------------------------------------------------------------
916 // This function compares tracks according to the sector - for given sector according z
917 //-----------------------------------------------------------------
918 AliTPCseed *t=(AliTPCseed*)o;
1c53abe2 919
91162307 920 if (fSort == 0){
921 if (t->fRelativeSector>fRelativeSector) return -1;
922 if (t->fRelativeSector<fRelativeSector) return 1;
923 Double_t z2 = t->GetZ();
924 Double_t z1 = GetZ();
925 if (z2>z1) return 1;
926 if (z2<z1) return -1;
927 return 0;
928 }
929 else {
930 Float_t f2 =1;
931 f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
932 if (t->fBConstrain) f2=1.2;
933
934 Float_t f1 =1;
935 f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
936
937 if (fBConstrain) f1=1.2;
938
939 if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
940 else return +1;
941 }
1c53abe2 942}
943
c9427e08 944void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
945{
946 //rotate to track "local coordinata
947 Float_t x = seed->GetX();
948 Float_t y = seed->GetY();
949 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
91162307 950
c9427e08 951 if (y > ymax) {
952 seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
953 if (!seed->Rotate(fSectors->GetAlpha()))
954 return;
955 } else if (y <-ymax) {
956 seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
957 if (!seed->Rotate(-fSectors->GetAlpha()))
958 return;
959 }
1c53abe2 960
c9427e08 961}
1c53abe2 962
963
964
965
966//_____________________________________________________________________________
967Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
968 //-----------------------------------------------------------------
969 // This function associates a cluster with this track.
970 //-----------------------------------------------------------------
1c53abe2 971 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
972
973 r00+=fC00; r01+=fC10; r11+=fC11;
974 Double_t det=r00*r11 - r01*r01;
975 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
976
977 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
978 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
979 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
980 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
981 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
982
983 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
984 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1627d1c4 985 if (TMath::Abs(cur*fX-eta) >= 0.9) {
1c53abe2 986 return 0;
987 }
988
989 fP0 += k00*dy + k01*dz;
990 fP1 += k10*dy + k11*dz;
991 fP2 = eta;
992 fP3 += k30*dy + k31*dz;
993 fP4 = cur;
994
995 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
996 Double_t c12=fC21, c13=fC31, c14=fC41;
997
998 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
999 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1000 fC40-=k00*c04+k01*c14;
1001
1002 fC11-=k10*c01+k11*fC11;
1003 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1004 fC41-=k10*c04+k11*c14;
1005
1006 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1007 fC42-=k20*c04+k21*c14;
1008
1009 fC33-=k30*c03+k31*c13;
1010 fC43-=k40*c03+k41*c13;
1011
1012 fC44-=k40*c04+k41*c14;
1013
1014 Int_t n=GetNumberOfClusters();
91162307 1015 // fIndex[n]=index;
1c53abe2 1016 SetNumberOfClusters(n+1);
1017 SetChi2(GetChi2()+chisq);
1018
1019 return 1;
1020}
1021
1022
1023
1024//_____________________________________________________________________________
91162307 1025Double_t AliTPCtrackerMI::f1old(Double_t x1,Double_t y1,
1c53abe2 1026 Double_t x2,Double_t y2,
1027 Double_t x3,Double_t y3)
1028{
1029 //-----------------------------------------------------------------
1030 // Initial approximation of the track curvature
1031 //-----------------------------------------------------------------
1032 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1033 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1034 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1035 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1036 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1037
1038 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
91162307 1039 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1c53abe2 1040 return -xr*yr/sqrt(xr*xr+yr*yr);
1041}
1042
1043
91162307 1044
1c53abe2 1045//_____________________________________________________________________________
91162307 1046Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
1047 Double_t x2,Double_t y2,
1048 Double_t x3,Double_t y3)
1049{
1050 //-----------------------------------------------------------------
1051 // Initial approximation of the track curvature
1052 //-----------------------------------------------------------------
1053 x3 -=x1;
1054 x2 -=x1;
1055 y3 -=y1;
1056 y2 -=y1;
1057 //
1058 Double_t det = x3*y2-x2*y3;
1059 if (det==0) {
1060 return 100;
1061 }
1062 //
1063 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1064 Double_t x0 = x3*0.5-y3*u;
1065 Double_t y0 = y3*0.5+x3*u;
1066 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1067 if (det<0) c2*=-1;
1068 return c2;
1069}
1070
1071
1c53abe2 1072Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
1073 Double_t x2,Double_t y2,
1074 Double_t x3,Double_t y3)
91162307 1075{
1076 //-----------------------------------------------------------------
1077 // Initial approximation of the track curvature
1078 //-----------------------------------------------------------------
1079 x3 -=x1;
1080 x2 -=x1;
1081 y3 -=y1;
1082 y2 -=y1;
1083 //
1084 Double_t det = x3*y2-x2*y3;
1085 if (det==0) {
1086 return 100;
1087 }
1088 //
1089 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1090 Double_t x0 = x3*0.5-y3*u;
1091 Double_t y0 = y3*0.5+x3*u;
1092 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1093 if (det<0) c2*=-1;
1094 x0+=x1;
1095 x0*=c2;
1096 return x0;
1097}
1098
1099
1100
1101//_____________________________________________________________________________
1102Double_t AliTPCtrackerMI::f2old(Double_t x1,Double_t y1,
1103 Double_t x2,Double_t y2,
1104 Double_t x3,Double_t y3)
1c53abe2 1105{
1106 //-----------------------------------------------------------------
1107 // Initial approximation of the track curvature times center of curvature
1108 //-----------------------------------------------------------------
1109 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1110 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1111 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1112 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1113 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1114
1115 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1116
1117 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1118}
1119
1120//_____________________________________________________________________________
1121Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
1122 Double_t x2,Double_t y2,
1123 Double_t z1,Double_t z2)
1124{
1125 //-----------------------------------------------------------------
1126 // Initial approximation of the tangent of the track dip angle
1127 //-----------------------------------------------------------------
1128 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1129}
1130
1131
91162307 1132Double_t AliTPCtrackerMI::f3n(Double_t x1,Double_t y1,
1133 Double_t x2,Double_t y2,
1134 Double_t z1,Double_t z2, Double_t c)
1c53abe2 1135{
91162307 1136 //-----------------------------------------------------------------
1137 // Initial approximation of the tangent of the track dip angle
1138 //-----------------------------------------------------------------
1139
1140 // Double_t angle1;
1141
1142 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1c53abe2 1143 //
91162307 1144 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1145 if (TMath::Abs(d*c*0.5)>1) return 0;
1146 // Double_t angle2 = TMath::ASin(d*c*0.5);
1147 // Double_t angle2 = TPCFastMath::FastAsin(d*c*0.5);
1148 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): TPCFastMath::FastAsin(d*c*0.5);
1149
1150 angle2 = (z1-z2)*c/(angle2*2.);
1151 return angle2;
1152}
1153
1154Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1155{//-----------------------------------------------------------------
1156 // This function find proloncation of a track to a reference plane x=x2.
1157 //-----------------------------------------------------------------
1158
1159 Double_t dx=x2-x1;
1160
1161 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1162 return kFALSE;
1c53abe2 1163 }
f8aae377 1164
91162307 1165 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1166 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1167 y = x[0];
1168 z = x[1];
1169
1170 Double_t dy = dx*(c1+c2)/(r1+r2);
1171 Double_t dz = 0;
1172 //
1173 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1174
1175 if (TMath::Abs(delta)>0.01){
1176 dz = x[3]*TMath::ASin(delta)/x[4];
1177 }else{
1178 dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
1179 }
1180
1181 //dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
f8aae377 1182
91162307 1183 y+=dy;
1184 z+=dz;
1185
1186 return kTRUE;
1c53abe2 1187}
1188
91162307 1189
1190
1191Int_t AliTPCtrackerMI::LoadClusters()
1c53abe2 1192{
1193 //
1194 // load clusters to the memory
91162307 1195 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1196 clrow->SetClass("AliTPCclusterMI");
1197 clrow->SetArray(0);
1198 clrow->GetArray()->ExpandCreateFast(10000);
1199 //
1200 // TTree * tree = fClustersArray.GetTree();
1201
1202 TTree * tree = fInput;
1203 TBranch * br = tree->GetBranch("Segment");
1204 br->SetAddress(&clrow);
1205 //
1206 Int_t j=Int_t(tree->GetEntries());
1c53abe2 1207 for (Int_t i=0; i<j; i++) {
91162307 1208 br->GetEntry(i);
1209 //
1210 Int_t sec,row;
1211 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1212 //
1213 AliTPCRow * tpcrow=0;
1214 Int_t left=0;
1215 if (sec<fkNIS*2){
1216 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1217 left = sec/fkNIS;
1218 }
1219 else{
1220 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1221 left = (sec-fkNIS*2)/fkNOS;
1222 }
1223 if (left ==0){
1224 tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1225 tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1226 for (Int_t i=0;i<tpcrow->fN1;i++)
1227 tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1228 }
1229 if (left ==1){
1230 tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1231 tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1232 for (Int_t i=0;i<tpcrow->fN2;i++)
1233 tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1234 }
1c53abe2 1235 }
91162307 1236 //
1237 delete clrow;
1238 LoadOuterSectors();
1239 LoadInnerSectors();
1240 return 0;
1c53abe2 1241}
1242
1243
91162307 1244void AliTPCtrackerMI::UnloadClusters()
1245{
1246 //
1247 // unload clusters from the memory
1248 //
1249 Int_t nrows = fOuterSec->GetNRows();
1250 for (Int_t sec = 0;sec<fkNOS;sec++)
1251 for (Int_t row = 0;row<nrows;row++){
1252 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1253 if (tpcrow){
1254 if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1255 if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1256 }
1c53abe2 1257 }
91162307 1258 //
1259 nrows = fInnerSec->GetNRows();
1260 for (Int_t sec = 0;sec<fkNIS;sec++)
1261 for (Int_t row = 0;row<nrows;row++){
1262 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1263 if (tpcrow){
1264 if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1265 if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1266 }
1267 }
1268
1269 return ;
1c53abe2 1270}
1271
1272
1273//_____________________________________________________________________________
91162307 1274Int_t AliTPCtrackerMI::LoadOuterSectors() {
1c53abe2 1275 //-----------------------------------------------------------------
91162307 1276 // This function fills outer TPC sectors with clusters.
1c53abe2 1277 //-----------------------------------------------------------------
91162307 1278 Int_t nrows = fOuterSec->GetNRows();
1279 UInt_t index=0;
1280 for (Int_t sec = 0;sec<fkNOS;sec++)
1281 for (Int_t row = 0;row<nrows;row++){
1282 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1283 Int_t sec2 = sec+2*fkNIS;
1284 //left
1285 Int_t ncl = tpcrow->fN1;
1286 while (ncl--) {
1287 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1288 index=(((sec2<<8)+row)<<16)+ncl;
1289 tpcrow->InsertCluster(c,index);
1290 }
1291 //right
1292 ncl = tpcrow->fN2;
1293 while (ncl--) {
1294 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1295 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1296 tpcrow->InsertCluster(c,index);
1297 }
1298 //
1299 // write indexes for fast acces
1300 //
1301 for (Int_t i=0;i<510;i++)
1302 tpcrow->fFastCluster[i]=-1;
1303 for (Int_t i=0;i<tpcrow->GetN();i++){
1304 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1305 tpcrow->fFastCluster[zi]=i; // write index
1306 }
1307 Int_t last = 0;
1308 for (Int_t i=0;i<510;i++){
1309 if (tpcrow->fFastCluster[i]<0)
1310 tpcrow->fFastCluster[i] = last;
1311 else
1312 last = tpcrow->fFastCluster[i];
1313 }
1314 }
1315 fN=fkNOS;
1316 fSectors=fOuterSec;
1317 return 0;
1318}
1319
1320
1321//_____________________________________________________________________________
1322Int_t AliTPCtrackerMI::LoadInnerSectors() {
1323 //-----------------------------------------------------------------
1324 // This function fills inner TPC sectors with clusters.
1325 //-----------------------------------------------------------------
1326 Int_t nrows = fInnerSec->GetNRows();
1327 UInt_t index=0;
1328 for (Int_t sec = 0;sec<fkNIS;sec++)
1329 for (Int_t row = 0;row<nrows;row++){
1330 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1331 //
1332 //left
1333 Int_t ncl = tpcrow->fN1;
1334 while (ncl--) {
1335 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1336 index=(((sec<<8)+row)<<16)+ncl;
1337 tpcrow->InsertCluster(c,index);
1338 }
1339 //right
1340 ncl = tpcrow->fN2;
1341 while (ncl--) {
1342 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1343 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1344 tpcrow->InsertCluster(c,index);
1345 }
1346 //
1347 // write indexes for fast acces
1348 //
1349 for (Int_t i=0;i<510;i++)
1350 tpcrow->fFastCluster[i]=-1;
1351 for (Int_t i=0;i<tpcrow->GetN();i++){
1352 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1353 tpcrow->fFastCluster[zi]=i; // write index
1354 }
1355 Int_t last = 0;
1356 for (Int_t i=0;i<510;i++){
1357 if (tpcrow->fFastCluster[i]<0)
1358 tpcrow->fFastCluster[i] = last;
1359 else
1360 last = tpcrow->fFastCluster[i];
1361 }
1c53abe2 1362
91162307 1363 }
1364
1c53abe2 1365 fN=fkNIS;
1366 fSectors=fInnerSec;
91162307 1367 return 0;
1368}
1369
1370
1371
1372//_________________________________________________________________________
1373AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1374 //--------------------------------------------------------------------
1375 // Return pointer to a given cluster
1376 //--------------------------------------------------------------------
1377 Int_t sec=(index&0xff000000)>>24;
1378 Int_t row=(index&0x00ff0000)>>16;
1379 Int_t ncl=(index&0x0000ffff)>>00;
1380
1381 const AliTPCRow * tpcrow=0;
1382 AliTPCclusterMI * clrow =0;
1383 if (sec<fkNIS*2){
1384 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1385 if (sec<fkNIS)
1386 clrow = tpcrow->fClusters1;
1387 else
1388 clrow = tpcrow->fClusters2;
1389 }
1390 else{
1391 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1392 if (sec-2*fkNIS<fkNOS)
1393 clrow = tpcrow->fClusters1;
1394 else
1395 clrow = tpcrow->fClusters2;
1396 }
1397 if (tpcrow==0) return 0;
1398 if (tpcrow->GetN()<=ncl) return 0;
1399 // return (AliTPCclusterMI*)(*tpcrow)[ncl];
1400 return &(clrow[ncl]);
1401
1c53abe2 1402}
1403
91162307 1404
1405
1c53abe2 1406Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1407 //-----------------------------------------------------------------
1408 // This function tries to find a track prolongation to next pad row
1409 //-----------------------------------------------------------------
1c53abe2 1410 //
91162307 1411 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1627d1c4 1412
91162307 1413 // if (t.GetRadius()>x+10 ) return 0;
1414 // t.PropagateTo(x+0.02);
1415 //t.PropagateTo(x+0.01);
1627d1c4 1416 if (!t.PropagateTo(x)) {
91162307 1417 t.fRemoval = 10;
1627d1c4 1418 return 0;
1419 }
91162307 1420 //
1421 Double_t y=t.GetY(), z=t.GetZ();
1422 if (TMath::Abs(y)>ymax){
1423 if (y > ymax) {
1424 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1425 if (!t.Rotate(fSectors->GetAlpha()))
1426 return 0;
1427 } else if (y <-ymax) {
1428 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1429 if (!t.Rotate(-fSectors->GetAlpha()))
1430 return 0;
1431 }
1432 if (!t.PropagateTo(x)) {
1433 return 0;
1434 }
1435 y=t.GetY();
1436 }
1437 //
1438 // update current shape info every 3 pad-row
1439 if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1440 //t.fCurrentSigmaY = GetSigmaY(&t);
1441 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1442 GetShape(&t,nr);
1443 }
1c53abe2 1444 //
1445 AliTPCclusterMI *cl=0;
1446 UInt_t index=0;
91162307 1447
1448
1449 //Int_t nr2 = nr;
1450 if (t.GetClusterIndex2(nr)>0){
1451 //
1452 //cl = GetClusterMI(t.GetClusterIndex2(nr));
1453 index = t.GetClusterIndex2(nr);
1454 cl = t.fClusterPointer[nr];
1455 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1456 t.fCurrentClusterIndex1 = index;
1457 }
1458
1459 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1460 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1461 Double_t roady =1.;
1462 Double_t roadz = 1.;
1463 //
1c53abe2 1464 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1465 t.fInDead = kTRUE;
91162307 1466 t.SetClusterIndex2(nr,-1);
1c53abe2 1467 return 0;
1468 }
1469 else
1470 {
1627d1c4 1471 if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
1472 else
1473 return 0;
1c53abe2 1474 }
1475 //calculate
91162307 1476 if (cl){
c9427e08 1477 t.fCurrentCluster = cl;
91162307 1478 t.fRow = nr;
1479 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1480 if (accept<3) {
1481 //if founded cluster is acceptible
1482 UpdateTrack(&t,accept);
1483 return 1;
1484 }
1485 }
c9427e08 1486
91162307 1487 if (krow) {
1488 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1489 cl = krow.FindNearest2(y,z,roady,roadz,index);
1490 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1491 }
1492 // t.fNoCluster++;
c9427e08 1493
91162307 1494 if (cl) {
1495 t.fCurrentCluster = cl;
1496 t.fRow = nr;
1497 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
c9427e08 1498
91162307 1499 if (t.fCurrentCluster->IsUsed(10)){
1500 //
1501 //
c9427e08 1502
91162307 1503 t.fNShared++;
1504 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1505 t.fRemoval =10;
1506 return 0;
1507 }
1508 }
1509
1510 if (accept<3) UpdateTrack(&t,accept);
c9427e08 1511
91162307 1512 } else {
1513 if (t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1514
1515 }
1516 return 1;
1517}
c9427e08 1518
91162307 1519Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1520 //-----------------------------------------------------------------
1521 // This function tries to find a track prolongation to next pad row
1522 //-----------------------------------------------------------------
1523 //
1524 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1525 Double_t y,z;
1526 if (!t.GetProlongation(x,y,z)) {
1527 t.fRemoval = 10;
1528 return 0;
1529 }
1530 //
1531 //
1532 if (TMath::Abs(y)>ymax){
1533 return 0;
1534
1c53abe2 1535 if (y > ymax) {
1536 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1537 if (!t.Rotate(fSectors->GetAlpha()))
1538 return 0;
1539 } else if (y <-ymax) {
1540 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1541 if (!t.Rotate(-fSectors->GetAlpha()))
1542 return 0;
91162307 1543 }
1544 if (!t.PropagateTo(x)) {
1545 return 0;
1546 }
1547 t.GetProlongation(x,y,z);
1548 }
1549 //
1550 // update current shape info every 3 pad-row
1551 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1552 // t.fCurrentSigmaY = GetSigmaY(&t);
1553 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1554 GetShape(&t,nr);
1555 }
1556 //
1557 AliTPCclusterMI *cl=0;
1558 UInt_t index=0;
1559
1560
1561 //Int_t nr2 = nr;
1562 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1563 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1564 Double_t roady =1.;
1565 Double_t roadz = 1.;
1566 //
1567 Int_t row = nr;
1568 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1569 t.fInDead = kTRUE;
1570 t.SetClusterIndex2(row,-1);
1571 return 0;
1572 }
1573 else
1574 {
1575 if (TMath::Abs(z)>(1.05*x+10)) t.SetClusterIndex2(row,-1);
1c53abe2 1576 }
91162307 1577 //calculate
1578
1579 if ((cl==0)&&(krow)) {
1580 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1581 cl = krow.FindNearest2(y,z,roady,roadz,index);
1582
1583 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1584 }
1585
1586 if (cl) {
1587 t.fCurrentCluster = cl;
1588 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1589 //if (accept<3){
1590 t.SetClusterIndex2(row,index);
1591 t.fClusterPointer[row] = cl;
1592 //}
1c53abe2 1593 }
1594 return 1;
1595}
1596
1597
91162307 1598
1599Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1c53abe2 1600 //-----------------------------------------------------------------
1601 // This function tries to find a track prolongation to next pad row
1602 //-----------------------------------------------------------------
1603 t.fCurrentCluster = 0;
1604 t.fCurrentClusterIndex1 = 0;
1c53abe2 1605
1606 Double_t xt=t.GetX();
91162307 1607 Int_t row = GetRowNumber(xt)-1;
1608 Double_t ymax= GetMaxY(nr);
1609
1c53abe2 1610 if (row < nr) return 1; // don't prolongate if not information until now -
91162307 1611 if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1612 t.fRemoval =10;
1613 return 0; // not prolongate strongly inclined tracks
1614 }
1615 if (TMath::Abs(t.GetSnp())>0.95) {
1616 t.fRemoval =10;
1617 return 0; // not prolongate strongly inclined tracks
1618 }
1619
1620 Double_t x= GetXrow(nr);
1621 Double_t y,z;
1622 //t.PropagateTo(x+0.02);
1623 //t.PropagateTo(x+0.01);
1627d1c4 1624 if (!t.PropagateTo(x)){
1627d1c4 1625 return 0;
1626 }
1c53abe2 1627 //
91162307 1628 y=t.GetY();
1629 z=t.GetZ();
1c53abe2 1630
91162307 1631 if (TMath::Abs(y)>ymax){
1632 if (y > ymax) {
1633 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1634 if (!t.Rotate(fSectors->GetAlpha()))
1635 return 0;
1636 } else if (y <-ymax) {
1637 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1638 if (!t.Rotate(-fSectors->GetAlpha()))
1639 return 0;
1640 }
1641 if (!t.PropagateTo(x)){
1642 return 0;
1c53abe2 1643 }
91162307 1644 y = t.GetY();
1c53abe2 1645 }
91162307 1646 //
1c53abe2 1647
91162307 1648 AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1c53abe2 1649
1650 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1651 t.fInDead = kTRUE;
91162307 1652 t.SetClusterIndex2(nr,-1);
1c53abe2 1653 return 0;
1654 }
1655 else
1656 {
1627d1c4 1657 if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
1658 else
1659 return 0;
1c53abe2 1660 }
1c53abe2 1661
91162307 1662 // update current
1663 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1664 // t.fCurrentSigmaY = GetSigmaY(&t);
1665 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1666 GetShape(&t,nr);
1667 }
1c53abe2 1668
91162307 1669 AliTPCclusterMI *cl=0;
1670 UInt_t index=0;
1671 //
1672 Double_t roady = 1.;
1673 Double_t roadz = 1.;
1674 //
1675 if (krow) {
1676 //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1677 cl = krow.FindNearest2(y,z,roady,roadz,index);
1678 }
1679 t.fCurrentCluster = cl;
1680 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1681 return 1;
1682}
1c53abe2 1683
1c53abe2 1684
91162307 1685Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1686 //-----------------------------------------------------------------
1687 // This function tries to find a track prolongation to next pad row
1688 //-----------------------------------------------------------------
1c53abe2 1689
91162307 1690 //update error according neighborhoud
1c53abe2 1691
91162307 1692 if (t.fCurrentCluster) {
1693 t.fRow = nr;
1694 Bool_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1695
1696 if (t.fCurrentCluster->IsUsed(10)){
1697 //
1698 //
1699 // t.fErrorZ2*=2;
1700 // t.fErrorY2*=2;
1701 t.fNShared++;
1702 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1703 t.fRemoval =10;
1704 return 0;
1705 }
1706 }
1707
1708 if (accept<3) UpdateTrack(&t,accept);
1c53abe2 1709
1710 } else {
91162307 1711 if (fIteration==0){
1712 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1713 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1714
1715 if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1c53abe2 1716 }
1717 }
1718 return 1;
1719}
1720
1721
1722
91162307 1723//_____________________________________________________________________________
1724Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1725 //-----------------------------------------------------------------
91162307 1726 // This function tries to find a track prolongation.
1c53abe2 1727 //-----------------------------------------------------------------
1728 Double_t xt=t.GetX();
1729 //
91162307 1730 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1c53abe2 1731 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1732 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1733 //
1734 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1735
91162307 1736 Int_t first = GetRowNumber(xt)-1;
1737 for (Int_t nr= first; nr>=rf; nr-=step) {
1738 if (FollowToNext(t,nr)==0)
1739 if (!t.IsActive()) return 0;
1740
1741 }
1742 return 1;
1743}
1744
1c53abe2 1745
1746//_____________________________________________________________________________
91162307 1747Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1748 //-----------------------------------------------------------------
1749 // This function tries to find a track prolongation.
1750 //-----------------------------------------------------------------
1751 Double_t xt=t.GetX();
1752 //
1753 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1754 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1755 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1756 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1757
1758 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1759
1760 if (FollowToNextFast(t,nr)==0)
1761 if (!t.IsActive()) return 0;
1c53abe2 1762
1c53abe2 1763 }
1764 return 1;
1765}
1766
1767
91162307 1768
1769
1770
1c53abe2 1771Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1772 //-----------------------------------------------------------------
1773 // This function tries to find a track prolongation.
1774 //-----------------------------------------------------------------
91162307 1775 // Double_t xt=t.GetX();
1c53abe2 1776 //
1777 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1778 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1779 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1780 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1781
91162307 1782 Int_t first = 0;
1783 first = t.fFirstPoint+3;
1784 //
1785 if (first<0) first=0;
1786 for (Int_t nr=first+1; nr<=rf; nr++) {
1787 //if ( (t.GetSnp()<0.9))
1788 FollowToNext(t,nr);
1c53abe2 1789 }
1790 return 1;
1791}
1792
1793
1794
1795
1796
1797Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1798{
1799 //
1800 //
1801 sum1=0;
1802 sum2=0;
1803 Int_t sum=0;
1c53abe2 1804 //
1805 Float_t dz2 =(s1->GetZ() - s2->GetZ());
c9427e08 1806 dz2*=dz2;
91162307 1807
c9427e08 1808 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1c53abe2 1809 dy2*=dy2;
1810 Float_t distance = TMath::Sqrt(dz2+dy2);
c9427e08 1811 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1c53abe2 1812
91162307 1813 // Int_t offset =0;
c9427e08 1814 Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1815 Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
c9427e08 1816 if (lastpoint>160)
1817 lastpoint =160;
1818 if (firstpoint<0)
1819 firstpoint = 0;
91162307 1820 if (firstpoint>lastpoint) {
1821 firstpoint =lastpoint;
1822 // lastpoint =160;
c9427e08 1823 }
1824
1825
91162307 1826 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1827 if (s1->GetClusterIndex2(i)>0) sum1++;
1828 if (s2->GetClusterIndex2(i)>0) sum2++;
1829 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1c53abe2 1830 sum++;
1831 }
1832 }
91162307 1833 if (sum<5) return 0;
1834
1627d1c4 1835 Float_t summin = TMath::Min(sum1+1,sum2+1);
1836 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 1837 return ratio;
1838}
1839
1840void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1841{
1842 //
1843 //
91162307 1844 if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
1845 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
1846
1c53abe2 1847 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1848 dz2*=dz2;
1849 Float_t dy2 =(s1->GetY() - s2->GetY());
1850 dy2*=dy2;
91162307 1851 Float_t distance = dz2+dy2;
1852 if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics
1853
1c53abe2 1854 //
91162307 1855 Int_t sumshared=0;
1856 //
1857 Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
1858 Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
1859 //
1860 if (firstpoint>=lastpoint-5) return;;
1c53abe2 1861
91162307 1862 for (Int_t i=firstpoint;i<lastpoint;i++){
1863 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1864 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1865 sumshared++;
1866 }
1867 }
1868 if (sumshared>4){
1869 // sign clusters
1870 //
1871 for (Int_t i=firstpoint;i<lastpoint;i++){
1872 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1873 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1874 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
1875 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
1876 if (s1->IsActive()&&s2->IsActive()){
1877 p1->fIsShared = kTRUE;
1878 p2->fIsShared = kTRUE;
1879 }
1880 }
1881 }
1882 }
1883 //
1884 if (sumshared>10){
1885 for (Int_t i=0;i<4;i++){
1886 if (s1->fOverlapLabels[3*i]==0){
1887 s1->fOverlapLabels[3*i] = s2->GetLabel();
1888 s1->fOverlapLabels[3*i+1] = sumshared;
1889 s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
1890 break;
1891 }
1892 }
1893 for (Int_t i=0;i<4;i++){
1894 if (s2->fOverlapLabels[3*i]==0){
1895 s2->fOverlapLabels[3*i] = s1->GetLabel();
1896 s2->fOverlapLabels[3*i+1] = sumshared;
1897 s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
1898 break;
1899 }
1900 }
1901 }
1c53abe2 1902
91162307 1903}
1c53abe2 1904
91162307 1905void AliTPCtrackerMI::SignShared(TObjArray * arr)
1906{
1c53abe2 1907 //
91162307 1908 //sort trackss according sectors
1909 //
c9427e08 1910 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1911 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1912 if (!pt) continue;
1913 //if (pt) RotateToLocal(pt);
1914 pt->fSort = 0;
c9427e08 1915 }
91162307 1916 arr->UnSort();
1c53abe2 1917 arr->Sort(); // sorting according z
1918 arr->Expand(arr->GetEntries());
91162307 1919 //
1920 //
1c53abe2 1921 Int_t nseed=arr->GetEntriesFast();
1c53abe2 1922 for (Int_t i=0; i<nseed; i++) {
1923 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1924 if (!pt) continue;
1925 for (Int_t j=0;j<=12;j++){
1926 pt->fOverlapLabels[j] =0;
1c53abe2 1927 }
91162307 1928 }
1929 for (Int_t i=0; i<nseed; i++) {
1930 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1931 if (!pt) continue;
1932 if (pt->fRemoval>10) continue;
1c53abe2 1933 for (Int_t j=i+1; j<nseed; j++){
1934 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
91162307 1935 // if (pt2){
1936 if (pt2->fRemoval<=10) {
1937 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1938 SignShared(pt,pt2);
c9427e08 1939 }
91162307 1940 }
1941 }
1942}
1943
1944void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
1945{
1946 //
1947 //sort trackss according sectors
1948 //
1949 if (fDebug&1) {
1950 printf("Number of tracks before double removal- %d\n",arr->GetEntries());
1951 }
1952 //
1953 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1954 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1955 if (!pt) continue;
1956 pt->fSort = 0;
1957 }
1958 arr->UnSort();
1959 arr->Sort(); // sorting according z
1960 arr->Expand(arr->GetEntries());
1961 //
1962 //reset overlap labels
1963 //
1964 Int_t nseed=arr->GetEntriesFast();
1965 for (Int_t i=0; i<nseed; i++) {
1966 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1967 if (!pt) continue;
1968 pt->SetUniqueID(i);
1969 for (Int_t j=0;j<=12;j++){
1970 pt->fOverlapLabels[j] =0;
1c53abe2 1971 }
1972 }
91162307 1973 //
1974 //sign shared tracks
1c53abe2 1975 for (Int_t i=0; i<nseed; i++) {
1976 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1977 if (!pt) continue;
1978 if (pt->fRemoval>10) continue;
1979 Float_t deltac = pt->GetC()*0.1;
1980 for (Int_t j=i+1; j<nseed; j++){
1981 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1982 // if (pt2){
1983 if (pt2->fRemoval<=10) {
1984 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1985 if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue;
1986 if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
1987 //
1988 SignShared(pt,pt2);
1989 }
1c53abe2 1990 }
1c53abe2 1991 }
91162307 1992 //
1993 // remove highly shared tracks
1994 for (Int_t i=0; i<nseed; i++) {
1995 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1996 if (!pt) continue;
1997 if (pt->fRemoval>10) continue;
1998 //
1999 Int_t sumshared =0;
2000 for (Int_t j=0;j<4;j++){
2001 sumshared = pt->fOverlapLabels[j*3+1];
2002 }
2003 Float_t factor = factor1;
2004 if (pt->fRemoval>0) factor = factor2;
2005 if (sumshared/pt->GetNumberOfClusters()>factor){
2006 for (Int_t j=0;j<4;j++){
2007 if (pt->fOverlapLabels[3*j]==0) continue;
2008 if (pt->fOverlapLabels[3*j+1]<5) continue;
2009 if (pt->fRemoval==removalindex) continue;
2010 AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2011 if (!pt2) continue;
2012 if (pt2->GetSigma2C()<pt->GetSigma2C()){
2013 // pt->fRemoval = removalindex;
2014 delete arr->RemoveAt(i);
2015 break;
1c53abe2 2016 }
91162307 2017 }
1c53abe2 2018 }
91162307 2019 }
2020 arr->Compress();
2021 if (fDebug&1) {
2022 printf("Number of tracks after double removal- %d\n",arr->GetEntries());
2023 }
2024}
2025
2026
2027
1c53abe2 2028
2029
91162307 2030
2031void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode)
2032{
2033 //
2034 //sort tracks in array according mode criteria
2035 Int_t nseed = arr->GetEntriesFast();
2036 for (Int_t i=0; i<nseed; i++) {
2037 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2038 if (!pt) {
2039 continue;
2040 }
2041 pt->fSort = mode;
2042 }
2043 arr->UnSort();
2044 arr->Sort();
1c53abe2 2045}
c9427e08 2046
91162307 2047void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
c9427e08 2048{
2049
2050 //Loop over all tracks and remove "overlaps"
2051 //
2052 //
2053 Int_t nseed = arr->GetEntriesFast();
2054 Int_t good =0;
91162307 2055
c9427e08 2056 for (Int_t i=0; i<nseed; i++) {
2057 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2058 if (!pt) {
91162307 2059 delete arr->RemoveAt(i);
c9427e08 2060 }
91162307 2061 else{
2062 pt->fSort =1;
2063 pt->fBSigned = kFALSE;
c9427e08 2064 }
91162307 2065 }
2066 arr->Compress();
2067 nseed = arr->GetEntriesFast();
2068 arr->UnSort();
2069 arr->Sort();
2070 //
2071 //unsign used
2072 UnsignClusters();
2073 //
2074 for (Int_t i=0; i<nseed; i++) {
2075 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2076 if (!pt) {
2077 continue;
2078 }
2079 Int_t found,foundable,shared;
2080 if (pt->IsActive())
2081 pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2082 else
2083 pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2084 //
2085 Double_t factor = factor2;
2086 if (pt->fBConstrain) factor = factor1;
2087
2088 if ((Float_t(shared)/Float_t(found))>factor){
c9427e08 2089 pt->Desactivate(removalindex);
91162307 2090 continue;
2091 }
2092
2093 good++;
2094 for (Int_t i=0; i<160; i++) {
2095 Int_t index=pt->GetClusterIndex2(i);
2096 if (index<0 || index&0x8000 ) continue;
2097 AliTPCclusterMI *c= pt->fClusterPointer[i];
2098 if (!c) continue;
2099 // if (!c->IsUsed(10)) c->Use(10);
2100 //if (pt->IsActive())
2101 c->Use(10);
2102 //else
2103 // c->Use(5);
c9427e08 2104 }
91162307 2105
c9427e08 2106 }
2107 fNtracks = good;
c9427e08 2108
91162307 2109 printf("\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
c9427e08 2110}
2111
91162307 2112void AliTPCtrackerMI::UnsignClusters()
1c53abe2 2113{
91162307 2114 //
2115 // loop over all clusters and unsign them
2116 //
2117
2118 for (Int_t sec=0;sec<fkNIS;sec++){
2119 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2120 AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2121 for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2122 // if (cl[icl].IsUsed(10))
2123 cl[icl].Use(-1);
2124 cl = fInnerSec[sec][row].fClusters2;
2125 for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2126 //if (cl[icl].IsUsed(10))
2127 cl[icl].Use(-1);
2128 }
2129 }
2130
2131 for (Int_t sec=0;sec<fkNOS;sec++){
2132 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2133 AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2134 for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2135 //if (cl[icl].IsUsed(10))
2136 cl[icl].Use(-1);
2137 cl = fOuterSec[sec][row].fClusters2;
2138 for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2139 //if (cl[icl].IsUsed(10))
2140 cl[icl].Use(-1);
2141 }
2142 }
2143
1c53abe2 2144}
2145
91162307 2146
2147
732b9e78 2148void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
1c53abe2 2149{
2150 //
91162307 2151 //sign clusters to be "used"
2152 //
2153 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2154 // loop over "primaries"
2155
2156 Float_t sumdens=0;
2157 Float_t sumdens2=0;
2158 Float_t sumn =0;
2159 Float_t sumn2 =0;
2160 Float_t sumchi =0;
2161 Float_t sumchi2 =0;
2162
2163 Float_t sum =0;
2164
1c53abe2 2165 TStopwatch timer;
91162307 2166 timer.Start();
1c53abe2 2167
91162307 2168 Int_t nseed = arr->GetEntriesFast();
2169 for (Int_t i=0; i<nseed; i++) {
2170 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2171 if (!pt) {
2172 continue;
2173 }
2174 if (!(pt->IsActive())) continue;
2175 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2176 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2177 sumdens += dens;
2178 sumdens2+= dens*dens;
2179 sumn += pt->GetNumberOfClusters();
2180 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2181 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2182 if (chi2>5) chi2=5;
2183 sumchi +=chi2;
2184 sumchi2 +=chi2*chi2;
2185 sum++;
2186 }
1627d1c4 2187 }
91162307 2188
2189 Float_t mdensity = 0.9;
2190 Float_t meann = 130;
2191 Float_t meanchi = 1;
2192 Float_t sdensity = 0.1;
2193 Float_t smeann = 10;
2194 Float_t smeanchi =0.4;
1627d1c4 2195
91162307 2196
2197 if (sum>20){
2198 mdensity = sumdens/sum;
2199 meann = sumn/sum;
2200 meanchi = sumchi/sum;
2201 //
2202 sdensity = sumdens2/sum-mdensity*mdensity;
2203 sdensity = TMath::Sqrt(sdensity);
2204 //
2205 smeann = sumn2/sum-meann*meann;
2206 smeann = TMath::Sqrt(smeann);
2207 //
2208 smeanchi = sumchi2/sum - meanchi*meanchi;
2209 smeanchi = TMath::Sqrt(smeanchi);
2210 }
2211
2212
2213 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2214 //
2215 for (Int_t i=0; i<nseed; i++) {
2216 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2217 if (!pt) {
2218 continue;
1c53abe2 2219 }
91162307 2220 if (pt->fBSigned) continue;
2221 if (pt->fBConstrain) continue;
2222 //if (!(pt->IsActive())) continue;
2223 /*
2224 Int_t found,foundable,shared;
2225 pt->GetClusterStatistic(0,160,found, foundable,shared);
2226 if (shared/float(found)>0.3) {
2227 if (shared/float(found)>0.9 ){
2228 //delete arr->RemoveAt(i);
2229 }
2230 continue;
c9427e08 2231 }
91162307 2232 */
2233 Bool_t isok =kFALSE;
2234 if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2235 isok = kTRUE;
2236 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2237 isok =kTRUE;
2238 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2239 isok =kTRUE;
2240 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2241 isok =kTRUE;
2242
2243 if (isok)
2244 for (Int_t i=0; i<160; i++) {
2245 Int_t index=pt->GetClusterIndex2(i);
2246 if (index<0) continue;
2247 AliTPCclusterMI *c= pt->fClusterPointer[i];
2248 if (!c) continue;
2249 //if (!(c->IsUsed(10))) c->Use();
2250 c->Use(10);
2251 }
2252 }
2253
c9427e08 2254
1c53abe2 2255 //
91162307 2256 Double_t maxchi = meanchi+2.*smeanchi;
2257
2258 for (Int_t i=0; i<nseed; i++) {
2259 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2260 if (!pt) {
1c53abe2 2261 continue;
91162307 2262 }
2263 //if (!(pt->IsActive())) continue;
2264 if (pt->fBSigned) continue;
2265 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2266 if (chi>maxchi) continue;
2267
2268 Float_t bfactor=1;
2269 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2270
2271 //sign only tracks with enoug big density at the beginning
2272
2273 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2274
2275
2276 Double_t mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2277 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2278
2279 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2280 if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2281 minn=0;
2282
2283 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2284 //Int_t noc=pt->GetNumberOfClusters();
2285 pt->fBSigned = kTRUE;
2286 for (Int_t i=0; i<160; i++) {
2287
2288 Int_t index=pt->GetClusterIndex2(i);
2289 if (index<0) continue;
2290 AliTPCclusterMI *c= pt->fClusterPointer[i];
2291 if (!c) continue;
2292 // if (!(c->IsUsed(10))) c->Use();
2293 c->Use(10);
2294 }
1c53abe2 2295 }
91162307 2296 }
2297 // gLastCheck = nseed;
2298 // arr->Compress();
2299 if (fDebug>0){
2300 timer.Print();
2301 }
1c53abe2 2302}
2303
2304
91162307 2305void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2)
2306{
2307 // stop not active tracks
2308 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2309 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2310 Int_t nseed = arr->GetEntriesFast();
2311 //
2312 for (Int_t i=0; i<nseed; i++) {
2313 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2314 if (!pt) {
2315 continue;
2316 }
2317 if (!(pt->IsActive())) continue;
2318 StopNotActive(pt,row0,th0, th1,th2);
2319 }
2320}
1c53abe2 2321
1c53abe2 2322
1c53abe2 2323
91162307 2324void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2325 Float_t th2)
2326{
2327 // stop not active tracks
2328 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2329 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2330 Int_t sumgood1 = 0;
2331 Int_t sumgood2 = 0;
2332 Int_t foundable = 0;
2333 Int_t maxindex = seed->fLastPoint; //last foundable row
2334 if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2335 seed->Desactivate(10) ;
2336 return;
2337 }
1c53abe2 2338
91162307 2339 for (Int_t i=row0; i++;i<maxindex){
2340 Int_t index = seed->GetClusterIndex2(i);
2341 if (index!=-1) foundable++;
2342 //if (!c) continue;
2343 if (foundable<=30) sumgood1++;
2344 if (foundable<=50) {
2345 sumgood2++;
2346 }
2347 else{
2348 break;
2349 }
2350 }
2351 if (foundable>=30.){
2352 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2353 }
2354 if (foundable>=50)
2355 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2356}
1c53abe2 2357
1c53abe2 2358
1c53abe2 2359
91162307 2360Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2361{
2362 //
2363 // back propagation of ESD tracks
2364 //
1c53abe2 2365
91162307 2366 fEvent = event;
2367 ReadSeeds(event);
2368 PropagateBack(fSeeds);
2369 Int_t nseed = fSeeds->GetEntriesFast();
2370 for (Int_t i=0;i<nseed;i++){
2371 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2372 AliESDtrack *esd=event->GetTrack(i);
2373 seed->CookdEdx(0.02,0.06);
2374 CookLabel(seed,0.1); //For comparison only
2375 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2376 }
2377 fEvent =0;
2378 WriteTracks();
2379 return 0;
2380}
2381
2382
2383void AliTPCtrackerMI::DeleteSeeds()
2384{
2385 Int_t nseed = fSeeds->GetEntriesFast();
2386 for (Int_t i=0;i<nseed;i++){
2387 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2388 if (seed) delete fSeeds->RemoveAt(i);
2389 }
2390 delete fSeeds;
2391 fSeeds =0;
2392}
2393
2394void AliTPCtrackerMI::ReadSeeds(AliESD *event)
2395{
2396 //
2397 //read seeds from the event
2398
2399 Int_t nentr=event->GetNumberOfTracks();
2400 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
2401 if (fSeeds)
2402 DeleteSeeds();
2403 if (!fSeeds){
2404 fSeeds = new TObjArray;
2405 }
2406
2407 // Int_t ntrk=0;
2408 for (Int_t i=0; i<nentr; i++) {
2409 AliESDtrack *esd=event->GetTrack(i);
2410 ULong_t status=esd->GetStatus();
2411 const AliTPCtrack t(*esd);
2412 AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2413 if (status==AliESDtrack::kTPCin) seed->Modify(0.8);
2414 //
2415 //
2416 // rotate to the local coordinate system
2417
2418 fSectors=fInnerSec; fN=fkNIS;
2419
2420 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2421 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2422 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2423 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2424 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2425 alpha-=seed->GetAlpha();
2426 if (!seed->Rotate(alpha)) continue;
2427 //
2428 seed->PropagateTo(fSectors->GetX(0));
2429 //
2430 // Int_t index = esd->GetTPCindex();
2431 //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2432
2433 fSeeds->AddLast(seed);
2434 }
2435}
2436
2437
2438
2439//_____________________________________________________________________________
2440void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2441 Float_t deltay, Int_t ddsec) {
2442 //-----------------------------------------------------------------
2443 // This function creates track seeds.
2444 // SEEDING WITH VERTEX CONSTRAIN
2445 //-----------------------------------------------------------------
2446 // cuts[0] - fP4 cut
2447 // cuts[1] - tan(phi) cut
2448 // cuts[2] - zvertex cut
2449 // cuts[3] - fP3 cut
2450 Int_t nin0 = 0;
2451 Int_t nin1 = 0;
2452 Int_t nin2 = 0;
2453 Int_t nin = 0;
2454 Int_t nout1 = 0;
2455 Int_t nout2 = 0;
2456
2457 Double_t x[5], c[15];
2458 // Int_t di = i1-i2;
2459 //
2460 AliTPCseed * seed = new AliTPCseed;
2461 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2462 Double_t cs=cos(alpha), sn=sin(alpha);
2463 //
2464 // Double_t x1 =fOuterSec->GetX(i1);
2465 //Double_t xx2=fOuterSec->GetX(i2);
2466
2467 Double_t x1 =GetXrow(i1);
2468 Double_t xx2=GetXrow(i2);
2469
2470 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2471
2472 Int_t imiddle = (i2+i1)/2; //middle pad row index
2473 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2474 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2475 //
2476 Int_t ns =sec;
2477
2478 const AliTPCRow& kr1=GetRow(ns,i1);
2479 Double_t ymax = GetMaxY(i1)-kr1.fDeadZone-1.5;
2480 Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2481
2482 //
2483 // change cut on curvature if it can't reach this layer
2484 // maximal curvature set to reach it
2485 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2486 if (dvertexmax*0.5*cuts[0]>0.85){
2487 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2488 }
2489 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2490
2491 // Int_t ddsec = 1;
2492 if (deltay>0) ddsec = 0;
2493 // loop over clusters
2494 for (Int_t is=0; is < kr1; is++) {
2495 //
2496 if (kr1[is]->IsUsed(10)) continue;
2497 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2498 //if (TMath::Abs(y1)>ymax) continue;
2499
2500 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2501
2502 // find possible directions
2503 Float_t anglez = (z1-z3)/(x1-x3);
2504 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2505 //
2506 //
2507 //find rotation angles relative to line given by vertex and point 1
2508 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2509 Double_t dvertex = TMath::Sqrt(dvertex2);
2510 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2511 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2512
2513 //
2514 // loop over 2 sectors
2515 Int_t dsec1=-ddsec;
2516 Int_t dsec2= ddsec;
2517 if (y1<0) dsec2= 0;
2518 if (y1>0) dsec1= 0;
2519
2520 Double_t dddz1=0; // direction of delta inclination in z axis
2521 Double_t dddz2=0;
2522 if ( (z1-z3)>0)
2523 dddz1 =1;
2524 else
2525 dddz2 =1;
2526 //
2527 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2528 Int_t sec2 = sec + dsec;
2529 //
2530 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2531 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2532 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2533 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2534 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2535 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2536
2537 // rotation angles to p1-p3
2538 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2539 Double_t x2, y2, z2;
2540 //
2541 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2542
2543 //
2544 Double_t dxx0 = (xx2-x3)*cs13r;
2545 Double_t dyy0 = (xx2-x3)*sn13r;
2546 for (Int_t js=index1; js < index2; js++) {
2547 const AliTPCclusterMI *kcl = kr2[js];
2548 if (kcl->IsUsed(10)) continue;
2549 //
2550 //calcutate parameters
2551 //
2552 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2553 // stright track
2554 if (TMath::Abs(yy0)<0.000001) continue;
2555 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2556 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2557 Double_t r02 = (0.25+y0*y0)*dvertex2;
2558 //curvature (radius) cut
2559 if (r02<r2min) continue;
2560
2561 nin0++;
2562 //
2563 Double_t c0 = 1/TMath::Sqrt(r02);
2564 if (yy0>0) c0*=-1.;
2565
2566
2567 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2568 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2569 Double_t dfi0 = 2.*TPCFastMath::FastAsin(dvertex*c0*0.5);
2570 Double_t dfi1 = 2.*TPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2571 //
2572 //
2573 Double_t z0 = kcl->GetZ();
2574 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2575 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2576 nin1++;
2577 //
2578 Double_t dip = (z1-z0)*c0/dfi1;
2579 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2580 //
2581 y2 = kcl->GetY();
2582 if (dsec==0){
2583 x2 = xx2;
2584 z2 = kcl->GetZ();
2585 }
2586 else
2587 {
2588 // rotation
2589 z2 = kcl->GetZ();
2590 x2= xx2*cs-y2*sn*dsec;
2591 y2=+xx2*sn*dsec+y2*cs;
2592 }
2593
2594 x[0] = y1;
2595 x[1] = z1;
2596 x[2] = x0;
2597 x[3] = dip;
2598 x[4] = c0;
2599 //
2600 //
2601 // do we have cluster at the middle ?
2602 Double_t ym,zm;
2603 GetProlongation(x1,xm,x,ym,zm);
2604 UInt_t dummy;
2605 AliTPCclusterMI * cm=0;
2606 if (TMath::Abs(ym)-ymaxm<0){
2607 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2608 if ((!cm) || (cm->IsUsed(10))) {
2609 continue;
2610 }
2611 }
2612 else{
2613 // rotate y1 to system 0
2614 // get state vector in rotated system
2615 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2616 Double_t xr2 = x0*cs+yr1*sn*dsec;
2617 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2618 //
2619 GetProlongation(xx2,xm,xr,ym,zm);
2620 if (TMath::Abs(ym)-ymaxm<0){
2621 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2622 if ((!cm) || (cm->IsUsed(10))) {
2623 continue;
2624 }
2625 }
2626 }
2627
2628
2629 Double_t dym = 0;
2630 Double_t dzm = 0;
2631 if (cm){
2632 dym = ym - cm->GetY();
2633 dzm = zm - cm->GetZ();
2634 }
2635 nin2++;
2636
2637
2638 //
2639 //
2640 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2641 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2642 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2643 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2644 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2645
2646 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2647 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2648 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2649 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2650 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1c53abe2 2651 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 2652
1c53abe2 2653 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2654 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2655 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2656 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
91162307 2657
1c53abe2 2658 c[0]=sy1;
2659 c[1]=0.; c[2]=sz1;
2660 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2661 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2662 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2663 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2664 c[13]=f30*sy1*f40+f32*sy2*f42;
2665 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
91162307 2666
2667 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2668
1c53abe2 2669 UInt_t index=kr1.GetIndex(is);
91162307 2670 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
2671
1c53abe2 2672 track->fIsSeeding = kTRUE;
91162307 2673 track->fSeed1 = i1;
2674 track->fSeed2 = i2;
2675 track->fSeedType=3;
c9427e08 2676
91162307 2677
2678 //if (dsec==0) {
2679 FollowProlongation(*track, (i1+i2)/2,1);
2680 Int_t foundable,found,shared;
2681 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2682 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2683 seed->Reset();
2684 seed->~AliTPCseed();
2685 continue;
2686 }
2687 //}
2688
2689 nin++;
2690 FollowProlongation(*track, i2,1);
2691
2692
2693 //Int_t rc = 1;
2694 track->fBConstrain =1;
2695 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2696 track->fLastPoint = i1; // first cluster in track position
2697 track->fFirstPoint = track->fLastPoint;
2698
2699 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2700 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
2701 track->fNShared>0.4*track->GetNumberOfClusters() ) {
2702 seed->Reset();
2703 seed->~AliTPCseed();
c9427e08 2704 continue;
2705 }
91162307 2706 nout1++;
2707 // Z VERTEX CONDITION
2708 Double_t zv;
2709 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2710 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2711 if (TMath::Abs(zv-z3)>cuts[2]) {
2712 FollowProlongation(*track, TMath::Max(i2-20,0));
2713 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2714 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2715 if (TMath::Abs(zv-z3)>cuts[2]){
2716 FollowProlongation(*track, TMath::Max(i2-40,0));
2717 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2718 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2719 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
2720 // make seed without constrain
2721 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
2722 FollowProlongation(*track2, i2,1);
2723 track2->fBConstrain = kFALSE;
2724 track2->fSeedType = 1;
2725 arr->AddLast(track2);
2726 seed->Reset();
2727 seed->~AliTPCseed();
2728 continue;
2729 }
2730 else{
2731 seed->Reset();
2732 seed->~AliTPCseed();
2733 continue;
2734
2735 }
2736 }
c9427e08 2737 }
1c53abe2 2738
91162307 2739 track->fSeedType =0;
2740 arr->AddLast(track);
2741 seed = new AliTPCseed;
2742 nout2++;
2743 // don't consider other combinations
2744 if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
2745 break;
1c53abe2 2746 }
2747 }
2748 }
91162307 2749 if (fDebug>1){
2750 // printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
2751 }
2752 delete seed;
1c53abe2 2753}
2754
1627d1c4 2755
91162307 2756void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2757 Float_t deltay) {
2758
2759
2760
1627d1c4 2761 //-----------------------------------------------------------------
91162307 2762 // This function creates track seeds.
1627d1c4 2763 //-----------------------------------------------------------------
91162307 2764 // cuts[0] - fP4 cut
2765 // cuts[1] - tan(phi) cut
2766 // cuts[2] - zvertex cut
2767 // cuts[3] - fP3 cut
2768
2769
2770 Int_t nin0 = 0;
2771 Int_t nin1 = 0;
2772 Int_t nin2 = 0;
2773 Int_t nin = 0;
2774 Int_t nout1 = 0;
2775 Int_t nout2 = 0;
2776 Int_t nout3 =0;
2777 Double_t x[5], c[15];
2778 //
2779 // make temporary seed
2780 AliTPCseed * seed = new AliTPCseed;
1627d1c4 2781 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
2782 // Double_t cs=cos(alpha), sn=sin(alpha);
91162307 2783 //
2784 //
1627d1c4 2785
91162307 2786 // first 3 padrows
2787 Double_t x1 = GetXrow(i1-1);
2788 const AliTPCRow& kr1=GetRow(sec,i1-1);
2789 Double_t y1max = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
2790 //
2791 Double_t x1p = GetXrow(i1);
2792 const AliTPCRow& kr1p=GetRow(sec,i1);
2793 //
2794 Double_t x1m = GetXrow(i1-2);
2795 const AliTPCRow& kr1m=GetRow(sec,i1-2);
1627d1c4 2796
91162307 2797 //
2798 //last 3 padrow for seeding
2799 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
2800 Double_t x3 = GetXrow(i1-7);
2801 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
2802 //
2803 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
2804 Double_t x3p = GetXrow(i1-6);
2805 //
2806 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
2807 Double_t x3m = GetXrow(i1-8);
1627d1c4 2808
91162307 2809 //
2810 //
2811 // middle padrow
2812 Int_t im = i1-4; //middle pad row index
2813 Double_t xm = GetXrow(im); // radius of middle pad-row
2814 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
2815 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
2816 //
2817 //
2818 Double_t deltax = x1-x3;
2819 Double_t dymax = deltax*cuts[1];
2820 Double_t dzmax = deltax*cuts[3];
2821 //
2822 // loop over clusters
2823 for (Int_t is=0; is < kr1; is++) {
1627d1c4 2824 //
91162307 2825 if (kr1[is]->IsUsed(10)) continue;
2826 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1627d1c4 2827 //
91162307 2828 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2829 //
2830 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
2831 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
2832 //
2833 Double_t y3, z3;
1627d1c4 2834 //
1627d1c4 2835 //
91162307 2836 UInt_t index;
2837 for (Int_t js=index1; js < index2; js++) {
2838 const AliTPCclusterMI *kcl = kr3[js];
2839 if (kcl->IsUsed(10)) continue;
2840 y3 = kcl->GetY();
2841 // apply angular cuts
2842 if (TMath::Abs(y1-y3)>dymax) continue;
2843 x3 = x3;
2844 z3 = kcl->GetZ();
2845 if (TMath::Abs(z1-z3)>dzmax) continue;
2846 //
2847 Double_t angley = (y1-y3)/(x1-x3);
2848 Double_t anglez = (z1-z3)/(x1-x3);
2849 //
2850 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
2851 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
2852 //
2853 Double_t yyym = angley*(xm-x1)+y1;
2854 Double_t zzzm = anglez*(xm-x1)+z1;
2855
2856 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
2857 if (!kcm) continue;
2858 if (kcm->IsUsed(10)) continue;
2859
2860 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
2861 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
2862 //
2863 //
2864 //
2865 Int_t used =0;
2866 Int_t found =0;
2867 //
2868 // look around first
2869 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
2870 anglez*(x1m-x1)+z1,
2871 erry,errz,index);
2872 //
2873 if (kc1m){
2874 found++;
2875 if (kc1m->IsUsed(10)) used++;
1627d1c4 2876 }
91162307 2877 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
2878 anglez*(x1p-x1)+z1,
2879 erry,errz,index);
1627d1c4 2880 //
91162307 2881 if (kc1p){
2882 found++;
2883 if (kc1p->IsUsed(10)) used++;
1627d1c4 2884 }
91162307 2885 if (used>1) continue;
2886 if (found<1) continue;
1627d1c4 2887
91162307 2888 //
2889 // look around last
2890 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
2891 anglez*(x3m-x3)+z3,
2892 erry,errz,index);
2893 //
2894 if (kc3m){
2895 found++;
2896 if (kc3m->IsUsed(10)) used++;
2897 }
2898 else
2899 continue;
2900 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
2901 anglez*(x3p-x3)+z3,
2902 erry,errz,index);
2903 //
2904 if (kc3p){
2905 found++;
2906 if (kc3p->IsUsed(10)) used++;
2907 }
2908 else
2909 continue;
2910 if (used>1) continue;
2911 if (found<3) continue;
2912 //
2913 Double_t x2,y2,z2;
2914 x2 = xm;
2915 y2 = kcm->GetY();
2916 z2 = kcm->GetZ();
2917 //
2918
1627d1c4 2919 x[0]=y1;
2920 x[1]=z1;
2921 x[4]=f1(x1,y1,x2,y2,x3,y3);
91162307 2922 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
2923 nin0++;
2924 //
1627d1c4 2925 x[2]=f2(x1,y1,x2,y2,x3,y3);
91162307 2926 nin1++;
2927 //
2928 x[3]=f3n(x1,y1,x2,y2,z1,z2,x[4]);
2929 //if (TMath::Abs(x[3]) > cuts[3]) continue;
2930 nin2++;
2931 //
2932 //
2933 Double_t sy1=0.1, sz1=0.1;
2934 Double_t sy2=0.1, sz2=0.1;
2935 Double_t sy3=0.1, sy=0.1, sz=0.1;
1627d1c4 2936
2937 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2938 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2939 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2940 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2941 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2942 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 2943
1627d1c4 2944 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2945 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2946 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2947 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2948
2949 c[0]=sy1;
91162307 2950 c[1]=0.; c[2]=sz1;
1627d1c4 2951 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2952 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2953 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2954 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2955 c[13]=f30*sy1*f40+f32*sy2*f42;
2956 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2957
91162307 2958 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2959
2960 UInt_t index=kr1.GetIndex(is);
2961 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
2962
2963 track->fIsSeeding = kTRUE;
2964
2965 nin++;
2966 FollowProlongation(*track, i1-7,1);
2967 if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
2968 track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
2969 seed->Reset();
2970 seed->~AliTPCseed();
2971 continue;
2972 }
2973 nout1++;
2974 nout2++;
2975 //Int_t rc = 1;
2976 FollowProlongation(*track, i2,1);
2977 track->fBConstrain =0;
2978 track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2979 track->fFirstPoint = track->fLastPoint;
2980
2981 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2982 track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
2983 track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
2984 seed->Reset();
2985 seed->~AliTPCseed();
2986 continue;
2987 }
2988
2989 {
2990 FollowProlongation(*track, TMath::Max(i2-10,0),1);
2991 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
2992 FollowProlongation(*track2, i2,1);
2993 track2->fBConstrain = kFALSE;
2994 track2->fSeedType = 4;
2995 arr->AddLast(track2);
2996 seed->Reset();
2997 seed->~AliTPCseed();
2998 }
2999
3000
3001 //arr->AddLast(track);
3002 //seed = new AliTPCseed;
3003 nout3++;
3004 }
3005 }
3006
3007 if (fDebug>1){
3008 // printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3009 }
3010 delete seed;
3011}
3012
3013
3014//_____________________________________________________________________________
3015void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3016 Float_t deltay, Bool_t bconstrain) {
3017 //-----------------------------------------------------------------
3018 // This function creates track seeds - without vertex constraint
3019 //-----------------------------------------------------------------
3020 // cuts[0] - fP4 cut - not applied
3021 // cuts[1] - tan(phi) cut
3022 // cuts[2] - zvertex cut - not applied
3023 // cuts[3] - fP3 cut
3024 Int_t nin0=0;
3025 Int_t nin1=0;
3026 Int_t nin2=0;
3027 Int_t nin3=0;
3028 // Int_t nin4=0;
3029 //Int_t nin5=0;
3030
3031
3032
3033 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3034 // Double_t cs=cos(alpha), sn=sin(alpha);
3035 Int_t row0 = (i1+i2)/2;
3036 Int_t drow = (i1-i2)/2;
3037 const AliTPCRow& kr0=fSectors[sec][row0];
3038 AliTPCRow * kr=0;
3039
3040 AliTPCpolyTrack polytrack;
3041 Int_t nclusters=fSectors[sec][row0];
3042 AliTPCseed * seed = new AliTPCseed;
3043
3044 Int_t sumused=0;
3045 Int_t cused=0;
3046 Int_t cnused=0;
3047 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3048 Int_t nfound =0;
3049 Int_t nfoundable =0;
3050 for (Int_t iter =1; iter<2; iter++){ //iterations
3051 const AliTPCRow& krm=fSectors[sec][row0-iter];
3052 const AliTPCRow& krp=fSectors[sec][row0+iter];
3053 const AliTPCclusterMI * cl= kr0[is];
3054
3055 if (cl->IsUsed(10)) {
3056 cused++;
3057 }
3058 else{
3059 cnused++;
3060 }
3061 Double_t x = kr0.GetX();
3062 // Initialization of the polytrack
3063 nfound =0;
3064 nfoundable =0;
3065 polytrack.Reset();
3066 //
3067 Double_t y0= cl->GetY();
3068 Double_t z0= cl->GetZ();
3069 Float_t erry = 0;
3070 Float_t errz = 0;
3071
3072 Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3073 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3074
3075 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3076 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3077 polytrack.AddPoint(x,y0,z0,erry, errz);
3078
3079 sumused=0;
3080 if (cl->IsUsed(10)) sumused++;
3081
3082
3083 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3084 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3085 //
3086 x = krm.GetX();
3087 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3088 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3089 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3090 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3091 if (cl1->IsUsed(10)) sumused++;
3092 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3093 }
3094 //
3095 x = krp.GetX();
3096 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3097 if (cl2) {
3098 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3099 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3100 if (cl2->IsUsed(10)) sumused++;
3101 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3102 }
3103 //
3104 if (sumused>0) continue;
3105 nin0++;
3106 polytrack.UpdateParameters();
3107 // follow polytrack
3108 roadz = 1.2;
3109 roady = 1.2;
3110 //
3111 Double_t yn,zn;
3112 nfoundable = polytrack.GetN();
3113 nfound = nfoundable;
3114 //
3115 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3116 Float_t maxdist = 0.8*(1.+3./(ddrow));
3117 for (Int_t delta = -1;delta<=1;delta+=2){
3118 Int_t row = row0+ddrow*delta;
3119 kr = &(fSectors[sec][row]);
3120 Double_t xn = kr->GetX();
3121 Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3122 polytrack.GetFitPoint(xn,yn,zn);
3123 if (TMath::Abs(yn)>ymax) continue;
3124 nfoundable++;
3125 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3126 if (cln) {
3127 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3128 if (dist<maxdist){
3129 /*
3130 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3131 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3132 if (cln->IsUsed(10)) {
3133 // printf("used\n");
3134 sumused++;
3135 erry*=2;
3136 errz*=2;
3137 }
3138 */
3139 erry=0.1;
3140 errz=0.1;
3141 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3142 nfound++;
3143 }
3144 }
3145 }
3146 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3147 polytrack.UpdateParameters();
3148 }
3149 }
3150 if ( (sumused>3) || (sumused>0.5*nfound)) {
3151 //printf("sumused %d\n",sumused);
3152 continue;
3153 }
3154 nin1++;
3155 Double_t dy,dz;
3156 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3157 AliTPCpolyTrack track2;
3158
3159 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3160 if (track2.GetN()<0.5*nfoundable) continue;
3161 nin2++;
3162
3163 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3164 //
3165 // test seed with and without constrain
3166 for (Int_t constrain=0; constrain<=0;constrain++){
3167 // add polytrack candidate
3168
3169 Double_t x[5], c[15];
3170 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3171 track2.GetBoundaries(x3,x1);
3172 x2 = (x1+x3)/2.;
3173 track2.GetFitPoint(x1,y1,z1);
3174 track2.GetFitPoint(x2,y2,z2);
3175 track2.GetFitPoint(x3,y3,z3);
3176 //
3177 //is track pointing to the vertex ?
3178 Double_t x0,y0,z0;
3179 x0=0;
3180 polytrack.GetFitPoint(x0,y0,z0);
3181
3182 if (constrain) {
3183 x2 = x3;
3184 y2 = y3;
3185 z2 = z3;
3186
3187 x3 = 0;
3188 y3 = 0;
3189 z3 = 0;
3190 }
3191 x[0]=y1;
3192 x[1]=z1;
3193 x[4]=f1(x1,y1,x2,y2,x3,y3);
3194
3195 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3196 x[2]=f2(x1,y1,x2,y2,x3,y3);
3197
3198 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3199 //x[3]=f3(x1,y1,x2,y2,z1,z2);
3200 x[3]=f3n(x1,y1,x3,y3,z1,z3,x[4]);
3201 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3202
3203
3204 Double_t sy =0.1, sz =0.1;
3205 Double_t sy1=0.02, sz1=0.02;
3206 Double_t sy2=0.02, sz2=0.02;
3207 Double_t sy3=0.02;
3208
3209 if (constrain){
3210 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3211 }
3212
3213 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3214 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3215 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3216 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3217 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3218 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3219
3220 Double_t f30=(f3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3221 Double_t f31=(f3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3222 Double_t f32=(f3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3223 Double_t f34=(f3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3224
3225
3226 c[0]=sy1;
3227 c[1]=0.; c[2]=sz1;
3228 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3229 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3230 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3231 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3232 c[13]=f30*sy1*f40+f32*sy2*f42;
3233 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3234
3235 //Int_t row1 = fSectors->GetRowNumber(x1);
3236 Int_t row1 = GetRowNumber(x1);
3237
3238 UInt_t index=0;
3239 //kr0.GetIndex(is);
3240 AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3241 track->fIsSeeding = kTRUE;
3242 Int_t rc=FollowProlongation(*track, i2);
3243 if (constrain) track->fBConstrain =1;
3244 else
3245 track->fBConstrain =0;
3246 track->fLastPoint = row1+fInnerSec->GetNRows(); // first cluster in track position
3247 track->fFirstPoint = track->fLastPoint;
3248
3249 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3250 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3251 track->fNShared>0.4*track->GetNumberOfClusters()) {
3252 //delete track;
3253 seed->Reset();
3254 seed->~AliTPCseed();
3255 }
3256 else {
3257 arr->AddLast(track);
3258 seed = new AliTPCseed;
3259 }
3260 nin3++;
3261 }
3262 } // if accepted seed
3263 }
3264 if (fDebug>1){
3265 printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3266 }
3267 delete seed;
3268}
3269
3270
3271AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3272{
3273 //
3274 //
3275 //reseed
3276 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3277 Int_t p1 = int(r1*track->GetNumberOfClusters());
3278 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3279 Int_t pp2;
3280 Double_t x0[3],x1[3],x2[3];
3281 x0[0]=-1;
3282 x0[0]=-1;
3283 x0[0]=-1;
3284
3285 // find track position at given ratio of the length
3286 Int_t sec0, sec1, sec2;
3287 Int_t index=-1;
3288 Int_t clindex;
3289 for (Int_t i=0;i<160;i++){
3290 if (track->fClusterPointer[i]){
3291 index++;
3292 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3293 if ( (index<p0) || x0[0]<0 ){
3294 if (trpoint->GetX()>1){
3295 clindex = track->GetClusterIndex2(i);
3296 if (clindex>0){
3297 x0[0] = trpoint->GetX();
3298 x0[1] = trpoint->GetY();
3299 x0[2] = trpoint->GetZ();
3300 sec0 = ((clindex&0xff000000)>>24)%18;
3301 }
3302 }
3303 }
3304
3305 if ( (index<p1) &&(trpoint->GetX()>1)){
3306 clindex = track->GetClusterIndex2(i);
3307 if (clindex>0){
3308 x1[0] = trpoint->GetX();
3309 x1[1] = trpoint->GetY();
3310 x1[2] = trpoint->GetZ();
3311 sec1 = ((clindex&0xff000000)>>24)%18;
3312 }
3313 }
3314 if ( (index<p2) &&(trpoint->GetX()>1)){
3315 clindex = track->GetClusterIndex2(i);
3316 if (clindex>0){
3317 x2[0] = trpoint->GetX();
3318 x2[1] = trpoint->GetY();
3319 x2[2] = trpoint->GetZ();
3320 sec2 = ((clindex&0xff000000)>>24)%18;
3321 pp2 = i;
3322 }
3323 }
3324 }
3325 }
3326
3327 Double_t alpha, cs,sn, xx2,yy2;
3328 //
3329 alpha = (sec1-sec2)*fSectors->GetAlpha();
3330 cs = TMath::Cos(alpha);
3331 sn = TMath::Sin(alpha);
3332 xx2= x1[0]*cs-x1[1]*sn;
3333 yy2= x1[0]*sn+x1[1]*cs;
3334 x1[0] = xx2;
3335 x1[1] = yy2;
3336 //
3337 alpha = (sec0-sec2)*fSectors->GetAlpha();
3338 cs = TMath::Cos(alpha);
3339 sn = TMath::Sin(alpha);
3340 xx2= x0[0]*cs-x0[1]*sn;
3341 yy2= x0[0]*sn+x0[1]*cs;
3342 x0[0] = xx2;
3343 x0[1] = yy2;
3344 //
3345 //
3346 //
3347 Double_t x[5],c[15];
3348 //
3349 x[0]=x2[1];
3350 x[1]=x2[2];
3351 x[4]=f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3352 // if (x[4]>1) return 0;
3353 x[2]=f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3354 x[3]=f3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3355 //if (TMath::Abs(x[3]) > 2.2) return 0;
3356 //if (TMath::Abs(x[2]) > 1.99) return 0;
3357 //
3358 Double_t sy =0.1, sz =0.1;
3359 //
3360 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3361 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3362 Double_t sy3=0.01+track->GetSigmaY2();
3363 //
3364 Double_t f40=(f1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3365 Double_t f42=(f1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3366 Double_t f43=(f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3367 Double_t f20=(f2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3368 Double_t f22=(f2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3369 Double_t f23=(f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3370 //
3371 Double_t f30=(f3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3372 Double_t f31=(f3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3373 Double_t f32=(f3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3374 Double_t f34=(f3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3375
3376
3377 c[0]=sy1;
3378 c[1]=0.; c[2]=sz1;
3379 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3380 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3381 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3382 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3383 c[13]=f30*sy1*f40+f32*sy2*f42;
3384 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3385
3386 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3387 AliTPCseed *seed=new AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3388 // Double_t y0,z0,y1,z1, y2,z2;
3389 //seed->GetProlongation(x0[0],y0,z0);
3390 // seed->GetProlongation(x1[0],y1,z1);
3391 //seed->GetProlongation(x2[0],y2,z2);
3392 // seed =0;
3393 seed->fLastPoint = pp2;
3394 seed->fFirstPoint = pp2;
3395
3396
3397 return seed;
3398}
3399
3400Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
3401{
3402 //
3403 //
3404 //
3405 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3406 //
3407 if (TMath::Abs(seed->GetC())>0.01) return 0;
3408 //
3409
3410 Float_t x[160], y[160], erry[160], z[160], errz[160];
3411 Int_t sec[160];
3412 Float_t xt[160], yt[160], zt[160];
3413 Int_t i1 = 200;
3414 Int_t i2 = 0;
3415 Int_t secm = -1;
3416 Int_t padm = -1;
3417 Int_t middle = seed->GetNumberOfClusters()/2;
3418 //
3419 //
3420 // find central sector, get local cooordinates
3421 Int_t count = 0;
3422 for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) {
3423 sec[i]= seed->GetClusterSector(i)%18;
3424 x[i] = GetXrow(i);
3425 if (sec[i]>=0) {
3426 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3427 // if (cl==0) cl = GetClusterMI(seed->GetClusterIndex2(i));
3428 if (cl==0) {
3429 sec[i] = -1;
3430 continue;
3431 }
3432 //
3433 //
3434 if (i>i2) i2 = i; //last point with cluster
3435 if (i2<i1) i1 = i; //first point with cluster
3436 y[i] = cl->GetY();
3437 z[i] = cl->GetZ();
3438 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3439 xt[i] = x[i];
3440 yt[i] = point->GetY();
3441 zt[i] = point->GetZ();
3442
3443 if (point->GetX()>0){
3444 erry[i] = point->GetErrY();
3445 errz[i] = point->GetErrZ();
3446 }
3447
3448 count++;
3449 if (count<middle) {
3450 secm = sec[i]; //central sector
3451 padm = i; //middle point with cluster
3452 }
3453 }
3454 }
3455 //
3456 // rotate position to global coordinate system connected to sector at last the point
3457 //
3458 for (Int_t i=i1;i<=i2;i++){
3459 //
3460 if (sec[i]<0) continue;
3461 Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha();
3462 Double_t cs = TMath::Cos(alpha);
3463 Double_t sn = TMath::Sin(alpha);
3464 Float_t xx2= x[i]*cs+y[i]*sn;
3465 Float_t yy2= -x[i]*sn+y[i]*cs;
3466 x[i] = xx2;
3467 y[i] = yy2;
3468 //
3469 xx2= xt[i]*cs+yt[i]*sn;
3470 yy2= -xt[i]*sn+yt[i]*cs;
3471 xt[i] = xx2;
3472 yt[i] = yy2;
3473
3474 }
3475 //get "state" vector
3476 Double_t xh[5],xm = x[padm];
3477 xh[0]=yt[i2];
3478 xh[1]=zt[i2];
3479 xh[4]=f1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3480 xh[2]=f2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3481 xh[3]=f3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]);
3482 //
3483 //
3484 for (Int_t i=i1;i<=i2;i++){
3485 Double_t yy,zz;
3486 if (sec[i]<0) continue;
3487 GetProlongation(x[i2], x[i],xh,yy,zz);
3488 if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){
3489 //Double_t xxh[5];
3490 //xxh[4]=f1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3491 //xxh[2]=f2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3492 printf("problem\n");
3493 }
3494 y[i] = y[i] - yy;
3495 z[i] = z[i] - zz;
3496 }
3497 Float_t dyup[160],dydown[160], dzup[160], dzdown[160];
3498 Float_t yup[160], ydown[160], zup[160], zdown[160];
3499
3500 AliTPCpolyTrack ptrack1,ptrack2;
3501 //
3502 // derivation up
3503 for (Int_t i=i1;i<=i2;i++){
3504 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3505 if (!cl) continue;
3506 if (cl->GetType()<0) continue;
3507 if (cl->GetType()>10) continue;
3508
3509 if (sec[i]>=0){
3510 ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3511 }
3512 if (ptrack1.GetN()>4.){
3513 ptrack1.UpdateParameters();
3514 Double_t ddy,ddz;
3515 ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz);
3516 Double_t yy,zz;
3517 ptrack1.GetFitPoint(x[i]-xm,yy,zz);
3518
3519 dyup[i] = ddy;
3520 dzup[i] = ddz;
3521 yup[i] = yy;
3522 zup[i] = zz;
3523
3524 }
3525 else{
3526 dyup[i]=0.; //not enough points
3527 }
3528 }
3529 //
3530 // derivation down
3531 for (Int_t i=i2;i>=i1;i--){
3532 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3533 if (!cl) continue;
3534 if (cl->GetType()<0) continue;
3535 if (cl->GetType()>10) continue;
3536 if (sec[i]>=0){
3537 ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3538 }
3539 if (ptrack2.GetN()>4){
3540 ptrack2.UpdateParameters();
3541 Double_t ddy,ddz;
3542 ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz);
3543 Double_t yy,zz;
3544 ptrack2.GetFitPoint(x[i]-xm,yy,zz);
3545
3546 dydown[i] = ddy;
3547 dzdown[i] = ddz;
3548 ydown[i] = yy;
3549 zdown[i] = zz;
3550 }
3551 else{
3552 dydown[i]=0.; //not enough points
3553 }
3554 }
3555 //
3556 //
3557 // find maximal difference of the derivation
3558 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3559
3560
3561 for (Int_t i=i1+10;i<i2-10;i++){
3562 if ( (TMath::Abs(dydown[i])<0.00000001) || (TMath::Abs(dyup[i])<0.00000001) ||i<30)continue;
3563 // printf("%f\t%f\t%f\t%f\t%f\n",x[i],dydown[i],dyup[i],dzdown[i],dzup[i]);
3564 //
3565 Float_t ddy = TMath::Abs(dydown[i]-dyup[i]);
3566 Float_t ddz = TMath::Abs(dzdown[i]-dzup[i]);
3567 if ( (ddy+ddz)> th){
3568 seed->fKinkPoint[0] = i;
3569 seed->fKinkPoint[1] = ddy;
3570 seed->fKinkPoint[2] = ddz;
3571 th = ddy+ddz;
3572 }
3573 }
3574
3575 if (fTreeDebug){
3576 //
3577 //write information to the debug tree
3578 TBranch * br = fTreeDebug->GetBranch("debug");
3579 TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
3580 arr->ExpandCreateFast(i2-i1);
3581 br->SetAddress(&arr);
3582 //
3583 AliTPCclusterMI cldummy;
3584 cldummy.SetQ(0);
3585 AliTPCTrackPoint2 pdummy;
3586 pdummy.GetTPoint().fIsShared = 10;
3587 //
3588 Double_t alpha = sec[i2]*fSectors->GetAlpha();
3589 Double_t cs = TMath::Cos(alpha);
3590 Double_t sn = TMath::Sin(alpha);
3591
3592 for (Int_t i=i1;i<i2;i++){
3593 AliTPCTrackPoint2 *trpoint = (AliTPCTrackPoint2*)arr->UncheckedAt(i-i1);
3594 //cluster info
3595 AliTPCclusterMI * cl0 = seed->fClusterPointer[i];
3596 //
3597 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3598
3599 if (cl0){
3600 Double_t x = GetXrow(i);
3601 trpoint->GetTPoint() = *point;
3602 trpoint->GetCPoint() = *cl0;
3603 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
3604 trpoint->fID = seed->GetUniqueID();
3605 trpoint->fLab = seed->GetLabel();
3606 //
3607 trpoint->fGX = cs *x + sn*point->GetY();
3608 trpoint->fGY = -sn *x + cs*point->GetY() ;
3609 trpoint->fGZ = point->GetZ();
3610 //
3611 trpoint->fDY = y[i];
3612 trpoint->fDZ = z[i];
3613 //
3614 trpoint->fDYU = dyup[i];
3615 trpoint->fDZU = dzup[i];
3616 //
3617 trpoint->fDYD = dydown[i];
3618 trpoint->fDZD = dzdown[i];
3619 //
3620 if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){
3621 trpoint->fDDY = dydown[i]-dyup[i];
3622 trpoint->fDDZ = dzdown[i]-dzup[i];
3623 }else{
3624 trpoint->fDDY = 0.;
3625 trpoint->fDDZ = 0.;
3626 }
3627 }
3628 else{
3629 *trpoint = pdummy;
3630 trpoint->GetCPoint()= cldummy;
3631 trpoint->fID = -1;
3632 }
3633 //
1627d1c4 3634 }
91162307 3635 fTreeDebug->Fill();
1627d1c4 3636 }
3637
3638
91162307 3639 return 0;
1627d1c4 3640
3641}
3642
3643
3644
3645
3646
91162307 3647AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
3648{
3649 //
3650 // reseed - refit - track
3651 //
3652 Int_t first = 0;
3653 // Int_t last = fSectors->GetNRows()-1;
3654 //
3655 if (fSectors == fOuterSec){
3656 first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
3657 //last =
3658 }
3659 else
3660 first = t->fFirstPoint;
3661 //
3662 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
3663 FollowBackProlongation(*t,fSectors->GetNRows()-1);
3664 t->Reset(kFALSE);
3665 FollowProlongation(*t,first);
3666 return seed;
3667}
3668
3669
3670
3671
3672
3673
3674
1c53abe2 3675//_____________________________________________________________________________
3676Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
3677 //-----------------------------------------------------------------
3678 // This function reades track seeds.
3679 //-----------------------------------------------------------------
3680 TDirectory *savedir=gDirectory;
3681
3682 TFile *in=(TFile*)inp;
3683 if (!in->IsOpen()) {
3684 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
3685 return 1;
3686 }
3687
3688 in->cd();
3689 TTree *seedTree=(TTree*)in->Get("Seeds");
3690 if (!seedTree) {
3691 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
3692 cerr<<"can't get a tree with track seeds !\n";
3693 return 2;
3694 }
3695 AliTPCtrack *seed=new AliTPCtrack;
3696 seedTree->SetBranchAddress("tracks",&seed);
3697
3698 if (fSeeds==0) fSeeds=new TObjArray(15000);
3699
3700 Int_t n=(Int_t)seedTree->GetEntries();
3701 for (Int_t i=0; i<n; i++) {
3702 seedTree->GetEvent(i);
3703 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
3704 }
3705
3706 delete seed;
3707 delete seedTree;
3708 savedir->cd();
3709 return 0;
3710}
3711
3712//_____________________________________________________________________________
f8aae377 3713Int_t AliTPCtrackerMI::Clusters2Tracks() {
1c53abe2 3714 //-----------------------------------------------------------------
3715 // This is a track finder.
3716 //-----------------------------------------------------------------
91162307 3717 TDirectory *savedir=gDirectory;
1c53abe2 3718 TStopwatch timer;
91162307 3719 //
3720 if (!fInput) SetIO(); //set default IO using loaders
3721 if (!fInput){
3722 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
3723 return 1;
3724 }
1c53abe2 3725 LoadClusters();
91162307 3726 //
3727 fIteration = 0;
3728 fSeeds = Tracking();
1c53abe2 3729
1c53abe2 3730
91162307 3731 printf("Time for tracking: \t");timer.Print();timer.Start();
3732
3733 //activate again some tracks
3734 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
3735 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3736 if (!pt) continue;
3737 Int_t nc=t.GetNumberOfClusters();
3738 if (nc<20) {
3739 delete fSeeds->RemoveAt(i);
3740 continue;
3741 }
3742 if (pt->fRemoval==10) {
3743 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
3744 pt->Desactivate(10); // make track again active
3745 else{
3746 pt->Desactivate(20);
3747 delete fSeeds->RemoveAt(i);
3748 }
3749 }
3750 }
3751 RemoveDouble(fSeeds,0.2,0.6,11);
3752 //RemoveUsed(fSeeds,0.9,0.9,6);
3753 //RemoveUsed(fSeeds,0.8,0.8,6);
3754 //RemoveUsed(fSeeds,0.7,0.7,6);
3755 RemoveUsed(fSeeds,0.5,0.5,6);
c9427e08 3756
1c53abe2 3757 //
91162307 3758 Int_t nseed=fSeeds->GetEntriesFast();
3759 Int_t found = 0;
3760 for (Int_t i=0; i<nseed; i++) {
3761 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3762 if (!pt) continue;
3763 Int_t nc=t.GetNumberOfClusters();
3764 if (nc<15) {
3765 delete fSeeds->RemoveAt(i);
3766 continue;
3767 }
3768 CookLabel(pt,0.1); //For comparison only
3769 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3770 if ((pt->IsActive() || (pt->fRemoval==10) )){
3771 cerr<<found++<<'\r';
3772 }
3773 else
3774 delete fSeeds->RemoveAt(i);
3775 pt->fLab2 = i;
3776 }
3777
3778
3779 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
3780 SignShared(fSeeds);
3781 //RemoveUsed(fSeeds,0.9,0.9,6);
1c53abe2 3782 //
91162307 3783 nseed=fSeeds->GetEntriesFast();
3784 found = 0;
3785 for (Int_t i=0; i<nseed; i++) {
3786 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3787 if (!pt) continue;
3788 Int_t nc=t.GetNumberOfClusters();
3789 if (nc<15) {
3790 delete fSeeds->RemoveAt(i);
3791 continue;
3792 }
3793 t.SetUniqueID(i);
3794 t.CookdEdx(0.02,0.6);
3795 // CheckKinkPoint(&t,0.05);
3796 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3797 if ((pt->IsActive() || (pt->fRemoval==10) )){
3798 cerr<<found++<<'\r';
3799 }
3800 else
3801 delete fSeeds->RemoveAt(i);
3802 pt->fLab2 = i;
3803 }
3804
3805 SortTracks(fSeeds, 1);
1c53abe2 3806
91162307 3807 /*
3808 fIteration = 1;
3809 PrepareForBackProlongation(fSeeds,0.5);
3810 PropagateBack(fSeeds);
3811 printf("Time for back propagation: \t");timer.Print();timer.Start();
3812
3813 fIteration = 2;
1c53abe2 3814
91162307 3815 PrepareForProlongation(fSeeds,1.);
3816 PropagateForward();
3817
3818 fSectors = fOuterSec;
3819 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3820 fSectors = fInnerSec;
3821 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3822 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
3823 // RemoveUsed(fSeeds,0.7,0.7,6);
3824 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
3825
1c53abe2 3826 nseed=fSeeds->GetEntriesFast();
91162307 3827 found = 0;
3828 for (Int_t i=0; i<nseed; i++) {
1c53abe2 3829 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3830 if (!pt) continue;
3831 Int_t nc=t.GetNumberOfClusters();
91162307 3832 if (nc<15) {
3833 delete fSeeds->RemoveAt(i);
3834 continue;
3835 }
1c53abe2 3836 t.CookdEdx(0.02,0.6);
91162307 3837 // CookLabel(pt,0.1); //For comparison only
3838 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3839 if ((pt->IsActive() || (pt->fRemoval==10) )){
3840 cerr<<found++<<'\r';
3841 }
3842 else
3843 delete fSeeds->RemoveAt(i);
3844 pt->fLab2 = i;
1c53abe2 3845 }
91162307 3846 */
3847
c9427e08 3848 // fNTracks = found;
91162307 3849 printf("Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
3850 //
3851 if (fOutput) {
3852 WriteTracks();
3853 }
3854 if (!fNewIO) fOutput->Write();
3855 else
3856 AliRunLoader::GetDetectorLoader("TPC",AliConfig::fgkDefaultEventFolderName)->WriteTracks("OVERWRITE");
1c53abe2 3857
1c53abe2 3858
91162307 3859 cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
3860 savedir->cd();
3861 //if (seedtree) delete seedtree;
3862 // UnloadClusters();
3863 //printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1c53abe2 3864
1c53abe2 3865 return 0;
3866}
3867
91162307 3868void AliTPCtrackerMI::Tracking(TObjArray * arr)
3869{
3870 //
3871 // tracking of the seeds
3872 //
3873
3874 fSectors = fOuterSec;
3875 ParallelTracking(arr,150,63);
3876 fSectors = fOuterSec;
3877 ParallelTracking(arr,63,0);
3878}
3879
3880TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
3881{
3882 //
3883 //
3884 //tracking routine
3885 TObjArray * arr = new TObjArray;
3886 //
3887 fSectors = fOuterSec;
3888 TStopwatch timer;
3889 timer.Start();
3890 for (Int_t sec=0;sec<fkNOS;sec++){
3891 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
3892 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
3893 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
3894 }
3895 if (fDebug>0){
3896 printf("\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
3897 timer.Print();
3898 timer.Start();
3899 }
3900 Tracking(arr);
3901 if (fDebug>0){
3902 timer.Print();
3903 }
3904
3905 return arr;
3906}
3907
3908TObjArray * AliTPCtrackerMI::Tracking()
3909{
3910 //
3911 //
3912 TStopwatch timer;
3913 timer.Start();
3914 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
3915
3916 TObjArray * seeds = new TObjArray;
3917 TObjArray * arr=0;
3918
3919 Int_t gap =20;
3920 Float_t cuts[4];
3921 cuts[0] = 0.002;
3922 cuts[1] = 1.5;
3923 cuts[2] = 3.;
3924 cuts[3] = 3.;
3925 Float_t fnumber = 3.0;
3926 Float_t fdensity = 3.0;
3927
3928 //
3929 //find primaries
3930 cuts[0]=0.0066;
3931 for (Int_t delta = 0; delta<18; delta+=6){
3932 //
3933 cuts[0]=0.0070;
3934 cuts[1] = 1.5;
3935 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
3936 SumTracks(seeds,arr);
3937 SignClusters(seeds,fnumber,fdensity);
3938 //
3939 for (Int_t i=2;i<6;i+=2){
3940 // seed high pt tracks
3941 cuts[0]=0.0022;
3942 cuts[1]=0.3;
3943 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
3944 SumTracks(seeds,arr);
3945 SignClusters(seeds,fnumber,fdensity);
3946 }
3947 }
3948 fnumber = 4;
3949 fdensity = 4.;
3950 // RemoveUsed(seeds,0.9,0.9,1);
3951 // UnsignClusters();
3952 // SignClusters(seeds,fnumber,fdensity);
3953
3954 //find primaries
3955 cuts[0]=0.0077;
3956 for (Int_t delta = 20; delta<120; delta+=10){
3957 //
3958 // seed high pt tracks
3959 cuts[0]=0.0060;
3960 cuts[1]=0.3;
3961 cuts[2]=6.;
3962 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
3963 SumTracks(seeds,arr);
3964 SignClusters(seeds,fnumber,fdensity);
3965
3966 cuts[0]=0.003;
3967 cuts[1]=0.3;
3968 cuts[2]=6.;
3969 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
3970 SumTracks(seeds,arr);
3971 SignClusters(seeds,fnumber,fdensity);
3972 }
3973
3974 cuts[0] = 0.01;
3975 cuts[1] = 2.0;
3976 cuts[2] = 3.;
3977 cuts[3] = 2.0;
3978 fnumber = 2.;
3979 fdensity = 2.;
3980
3981 if (fDebug>0){
3982 printf("\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
3983 timer.Print();
3984 timer.Start();
3985 }
3986 // RemoveUsed(seeds,0.75,0.75,1);
3987 //UnsignClusters();
3988 //SignClusters(seeds,fnumber,fdensity);
3989
3990 // find secondaries
3991
3992 cuts[0] = 0.3;
3993 cuts[1] = 1.5;
3994 cuts[2] = 3.;
3995 cuts[3] = 1.5;
3996
3997 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
3998 SumTracks(seeds,arr);
3999 SignClusters(seeds,fnumber,fdensity);
4000 //
4001 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
4002 SumTracks(seeds,arr);
4003 SignClusters(seeds,fnumber,fdensity);
4004 //
4005 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
4006 SumTracks(seeds,arr);
4007 SignClusters(seeds,fnumber,fdensity);
4008 //
4009
4010
4011 for (Int_t delta = 3; delta<30; delta+=5){
4012 //
4013 cuts[0] = 0.3;
4014 cuts[1] = 1.5;
4015 cuts[2] = 3.;
4016 cuts[3] = 1.5;
4017 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4018 SumTracks(seeds,arr);
4019 SignClusters(seeds,fnumber,fdensity);
4020 //
4021 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
4022 SumTracks(seeds,arr);
4023 SignClusters(seeds,fnumber,fdensity);
4024 //
4025 }
4026 fnumber = 1;
4027 fdensity = 1;
4028 //
4029 // change cuts
4030 fnumber = 2.;
4031 fdensity = 2.;
4032 cuts[0]=0.0080;
4033
4034 // find secondaries
4035 for (Int_t delta = 30; delta<70; delta+=10){
4036 //
4037 cuts[0] = 0.3;
4038 cuts[1] = 1.5;
4039 cuts[2] = 3.;
4040 cuts[3] = 1.5;
4041 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4042 SumTracks(seeds,arr);
4043 SignClusters(seeds,fnumber,fdensity);
4044 //
4045 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
4046 SumTracks(seeds,arr);
4047 SignClusters(seeds,fnumber,fdensity);
4048 }
4049
4050 if (fDebug>0){
4051 printf("\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
4052 timer.Print();
4053 timer.Start();
4054 }
4055
4056 return seeds;
4057 //
4058
4059}
4060
4061
4062void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2)
4063{
4064 //
4065 //sum tracks to common container
4066 //remove suspicious tracks
4067 Int_t nseed = arr2->GetEntriesFast();
4068 for (Int_t i=0;i<nseed;i++){
4069 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
4070 if (pt){
4071
4072 // NORMAL ACTIVE TRACK
4073 if (pt->IsActive()){
4074 arr1->AddLast(arr2->RemoveAt(i));
4075 continue;
4076 }
4077 //remove not usable tracks
4078 if (pt->fRemoval!=10){
4079 delete arr2->RemoveAt(i);
4080 continue;
4081 }
4082 // REMOVE VERY SHORT TRACKS
4083 if (pt->GetNumberOfClusters()<20){
4084 delete arr2->RemoveAt(i);
4085 continue;
4086 }
4087 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
4088 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4089 arr1->AddLast(arr2->RemoveAt(i));
4090 else{
4091 delete arr2->RemoveAt(i);
4092 }
4093 }
4094 }
4095 delete arr2;
4096}
4097
4098
1c53abe2 4099
91162307 4100void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
1c53abe2 4101{
4102 //
4103 // try to track in parralel
4104
91162307 4105 Int_t nseed=arr->GetEntriesFast();
1c53abe2 4106 //prepare seeds for tracking
4107 for (Int_t i=0; i<nseed; i++) {
91162307 4108 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1c53abe2 4109 if (!pt) continue;
4110 if (!t.IsActive()) continue;
4111 // follow prolongation to the first layer
91162307 4112 if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
c9427e08 4113 FollowProlongation(t, rfirst+1);
1c53abe2 4114 }
4115
4116
4117 //
4118 for (Int_t nr=rfirst; nr>=rlast; nr--){
4119 // make indexes with the cluster tracks for given
1c53abe2 4120
4121 // find nearest cluster
4122 for (Int_t i=0; i<nseed; i++) {
91162307 4123 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1c53abe2 4124 if (!pt) continue;
4125 if (!pt->IsActive()) continue;
91162307 4126 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
1627d1c4 4127 if (pt->fRelativeSector>17) {
4128 continue;
4129 }
91162307 4130 UpdateClusters(t,nr);
1c53abe2 4131 }
4132 // prolonagate to the nearest cluster - if founded
4133 for (Int_t i=0; i<nseed; i++) {
91162307 4134 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1c53abe2 4135 if (!pt) continue;
1627d1c4 4136 if (!pt->IsActive()) continue;
91162307 4137 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
1627d1c4 4138 if (pt->fRelativeSector>17) {
4139 continue;
4140 }
91162307 4141 FollowToNextCluster(*pt,nr);
1c53abe2 4142 }
91162307 4143 }
4144}
4145
4146void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac)
4147{
4148 //
4149 //
4150 // if we use TPC track itself we have to "update" covariance
4151 //
4152 Int_t nseed= arr->GetEntriesFast();
4153 for (Int_t i=0;i<nseed;i++){
4154 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4155 if (pt) {
4156 pt->Modify(fac);
4157 //
4158 //rotate to current local system at first accepted point
4159 Int_t index = pt->GetClusterIndex2(pt->fFirstPoint);
4160 Int_t sec = (index&0xff000000)>>24;
4161 sec = sec%18;
4162 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
4163 if (angle1>TMath::Pi())
4164 angle1-=2.*TMath::Pi();
4165 Float_t angle2 = pt->GetAlpha();
4166
4167 if (TMath::Abs(angle1-angle2)>0.001){
4168 pt->Rotate(angle1-angle2);
4169 //angle2 = pt->GetAlpha();
4170 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
4171 //if (pt->GetAlpha()<0)
4172 // pt->fRelativeSector+=18;
4173 //sec = pt->fRelativeSector;
4174 }
4175
4176 }
4177
4178 }
4179
4180
4181}
4182void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac)
4183{
4184 //
4185 //
4186 // if we use TPC track itself we have to "update" covariance
4187 //
4188 Int_t nseed= arr->GetEntriesFast();
4189 for (Int_t i=0;i<nseed;i++){
4190 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4191 if (pt) {
4192 pt->Modify(fac);
4193 pt->fFirstPoint = pt->fLastPoint;
4194 }
4195
4196 }
4197
4198
4199}
4200
4201Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
4202{
4203 //
4204 // make back propagation
4205 //
4206 Int_t nseed= arr->GetEntriesFast();
4207 for (Int_t i=0;i<nseed;i++){
4208 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4209 if (pt) {
4210 AliTPCseed *pt2 = new AliTPCseed(*pt);
4211 fSectors = fInnerSec;
4212 FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4213 fSectors = fOuterSec;
4214 FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4215 fSectors = fOuterSec;
4216 if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4217 printf("\n%d",pt->GetLabel());
4218 fSectors = fInnerSec;
4219 FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4220 fSectors = fOuterSec;
4221 FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4222 fSectors = fOuterSec;
4223 }
4224 }
4225 }
4226 return 0;
4227}
4228
4229
4230Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
4231{
4232 //
4233 // make forward propagation
4234 //
4235 Int_t nseed= arr->GetEntriesFast();
4236 for (Int_t i=0;i<nseed;i++){
4237 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4238 if (pt) {
4239 AliTPCseed *pt2 = new AliTPCseed(*pt);
4240 fSectors = fOuterSec;
4241 FollowProlongation(*pt,0);
4242 fSectors = fOuterSec;
4243 FollowProlongation(*pt,0);
4244 fSectors = fInnerSec;
4245 if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4246 printf("\n%d",pt->GetLabel());
4247 fSectors = fOuterSec;
4248 FollowProlongation(*pt2,0);
4249 fSectors = fOuterSec;
4250 FollowProlongation(*pt2,0);
4251 fSectors = fOuterSec;
4252 }
4253 }
4254 }
4255 return 0;
4256}
4257
4258
4259Int_t AliTPCtrackerMI::PropagateForward()
4260{
4261 fSectors = fOuterSec;
4262 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4263 fSectors = fInnerSec;
4264 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4265 //WriteTracks();
4266 return 1;
4267}
4268
4269
4270
4271
4272
4273
4274Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
4275{
4276 //
4277 // make back propagation, in between row0 and row1
4278 //
1c53abe2 4279
91162307 4280 if (pt) {
4281 fSectors = fInnerSec;
4282 Int_t r1;
4283 //
4284 if (row1<fSectors->GetNRows())
4285 r1 = row1;
4286 else
4287 r1 = fSectors->GetNRows()-1;
4288
4289 if (row0<fSectors->GetNRows()&& r1>0 )
4290 FollowBackProlongation(*pt,r1);
4291 if (row1<=fSectors->GetNRows())
4292 return 0;
4293 //
4294 r1 = row1 - fSectors->GetNRows();
4295 if (r1<=0) return 0;
4296 if (r1>=fOuterSec->GetNRows()) return 0;
4297 fSectors = fOuterSec;
4298 return FollowBackProlongation(*pt,r1);
4299 }
4300 return 0;
4301}
4302
4303
4304
4305
4306void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
4307{
4308 //
4309 //
4310 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4311 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4312 Float_t padlength = GetPadPitchLength(row);
4313 //
4314 Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4315 Float_t angulary = seed->GetSnp();
4316 angulary = angulary*angulary/(1-angulary*angulary);
4317 seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
4318 //
4319 Float_t sresz = fParam->GetZSigma();
4320 Float_t angularz = seed->GetTgl();
4321 seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
4322 /*
4323 Float_t wy = GetSigmaY(seed);
4324 Float_t wz = GetSigmaZ(seed);
4325 wy*=wy;
4326 wz*=wz;
4327 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
4328 printf("problem\n");
4329 }
4330 */
1c53abe2 4331}
4332
91162307 4333
1c53abe2 4334Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
4335{
4336 //
4337 //
c9427e08 4338 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1c53abe2 4339 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4340 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4341 Float_t angular = seed->GetSnp();
4342 angular = angular*angular/(1-angular*angular);
4343 // angular*=angular;
4344 //angular = TMath::Sqrt(angular/(1-angular));
4345 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
4346 return res;
4347}
4348Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
4349{
4350 //
4351 //
c9427e08 4352 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1c53abe2 4353 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4354 Float_t sres = fParam->GetZSigma();
4355 Float_t angular = seed->GetTgl();
4356 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
4357 return res;
4358}
4359
4360
4361
1c53abe2 4362
4363//__________________________________________________________________________
91162307 4364void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
1c53abe2 4365 //--------------------------------------------------------------------
4366 //This function "cooks" a track label. If label<0, this track is fake.
4367 //--------------------------------------------------------------------
4368 Int_t noc=t->GetNumberOfClusters();
91162307 4369 if (noc<10){
4370 printf("\nnot founded prolongation\n\n\n");
4371 t->Dump();
4372 return ;
4373 }
4374 Int_t lb[160];
4375 Int_t mx[160];
4376 AliTPCclusterMI *clusters[160];
4377 //
4378 for (Int_t i=0;i<160;i++) {
4379 clusters[i]=0;
4380 lb[i]=mx[i]=0;
4381 }
1c53abe2 4382
4383 Int_t i;
91162307 4384 Int_t current=0;
4385 for (i=0; i<160 && current<noc; i++) {
4386
4387 Int_t index=t->GetClusterIndex2(i);
4388 if (index<=0) continue;
4389 if (index&0x8000) continue;
4390 //
4391 //clusters[current]=GetClusterMI(index);
4392 if (t->fClusterPointer[i]){
4393 clusters[current]=t->fClusterPointer[i];
4394 current++;
4395 }
1c53abe2 4396 }
91162307 4397 noc = current;
1c53abe2 4398
4399 Int_t lab=123456789;
4400 for (i=0; i<noc; i++) {
4401 AliTPCclusterMI *c=clusters[i];
91162307 4402 if (!c) continue;
1c53abe2 4403 lab=TMath::Abs(c->GetLabel(0));
4404 Int_t j;
4405 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
4406 lb[j]=lab;
4407 (mx[j])++;
4408 }
4409
4410 Int_t max=0;
4411 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
4412
4413 for (i=0; i<noc; i++) {
9918f10a 4414 AliTPCclusterMI *c=clusters[i];
91162307 4415 if (!c) continue;
1c53abe2 4416 if (TMath::Abs(c->GetLabel(1)) == lab ||
4417 TMath::Abs(c->GetLabel(2)) == lab ) max++;
4418 }
4419
4420 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
4421
4422 else {
4423 Int_t tail=Int_t(0.10*noc);
4424 max=0;
91162307 4425 Int_t ind=0;
4426 for (i=1; i<=160&&ind<tail; i++) {
4427 // AliTPCclusterMI *c=clusters[noc-i];
4428 AliTPCclusterMI *c=clusters[i];
4429 if (!c) continue;
1c53abe2 4430 if (lab == TMath::Abs(c->GetLabel(0)) ||
4431 lab == TMath::Abs(c->GetLabel(1)) ||
4432 lab == TMath::Abs(c->GetLabel(2))) max++;
91162307 4433 ind++;
1c53abe2 4434 }
4435 if (max < Int_t(0.5*tail)) lab=-lab;
4436 }
4437
4438 t->SetLabel(lab);
4439
91162307 4440 // delete[] lb;
4441 //delete[] mx;
4442 //delete[] clusters;
1c53abe2 4443}
4444
4445//_________________________________________________________________________
4446void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
4447 //-----------------------------------------------------------------------
4448 // Setup inner sector
4449 //-----------------------------------------------------------------------
4450 if (f==0) {
4451 fAlpha=par->GetInnerAngle();
4452 fAlphaShift=par->GetInnerAngleShift();
4453 fPadPitchWidth=par->GetInnerPadPitchWidth();
4454 fPadPitchLength=par->GetInnerPadPitchLength();
4455 fN=par->GetNRowLow();
4456 fRow=new AliTPCRow[fN];
4457 for (Int_t i=0; i<fN; i++) {
4458 fRow[i].SetX(par->GetPadRowRadiiLow(i));
4459 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
4460 }
4461 } else {
4462 fAlpha=par->GetOuterAngle();
4463 fAlphaShift=par->GetOuterAngleShift();
4464 fPadPitchWidth = par->GetOuterPadPitchWidth();
4465 fPadPitchLength = par->GetOuter1PadPitchLength();
4466 f1PadPitchLength = par->GetOuter1PadPitchLength();
4467 f2PadPitchLength = par->GetOuter2PadPitchLength();
4468
4469 fN=par->GetNRowUp();
4470 fRow=new AliTPCRow[fN];
4471 for (Int_t i=0; i<fN; i++) {
4472 fRow[i].SetX(par->GetPadRowRadiiUp(i));
4473 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
4474 }
4475 }
4476}
4477
4478
4479AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
4480 //
1c53abe2 4481}
4482
4483
4484
1c53abe2 4485//_________________________________________________________________________
4486void
4487AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
4488 //-----------------------------------------------------------------------
4489 // Insert a cluster into this pad row in accordence with its y-coordinate
4490 //-----------------------------------------------------------------------
4491 if (fN==kMaxClusterPerRow) {
4492 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
4493 }
4494 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
4495 Int_t i=Find(c->GetZ());
4496 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
4497 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
4498 fIndex[i]=index; fClusters[i]=c; fN++;
4499}
4500
91162307 4501
1c53abe2 4502//___________________________________________________________________
4503Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
4504 //-----------------------------------------------------------------------
4505 // Return the index of the nearest cluster
4506 //-----------------------------------------------------------------------
4507 if (fN==0) return 0;
4508 if (z <= fClusters[0]->GetZ()) return 0;
4509 if (z > fClusters[fN-1]->GetZ()) return fN;
4510 Int_t b=0, e=fN-1, m=(b+e)/2;
4511 for (; b<e; m=(b+e)/2) {
4512 if (z > fClusters[m]->GetZ()) b=m+1;
4513 else e=m;
4514 }
4515 return m;
4516}
4517
4518
4519
1627d1c4 4520//___________________________________________________________________
4521AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
4522 //-----------------------------------------------------------------------
4523 // Return the index of the nearest cluster in z y
4524 //-----------------------------------------------------------------------
4525 Float_t maxdistance = roady*roady + roadz*roadz;
4526
4527 AliTPCclusterMI *cl =0;
4528 for (Int_t i=Find(z-roadz); i<fN; i++) {
4529 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4530 if (c->GetZ() > z+roadz) break;
4531 if ( (c->GetY()-y) > roady ) continue;
4532 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4533 if (maxdistance>distance) {
4534 maxdistance = distance;
4535 cl=c;
4536 }
4537 }
4538 return cl;
4539}
4540
91162307 4541AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4542{
4543 //-----------------------------------------------------------------------
4544 // Return the index of the nearest cluster in z y
4545 //-----------------------------------------------------------------------
4546 Float_t maxdistance = roady*roady + roadz*roadz;
4547 Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
4548 Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
4549
4550 AliTPCclusterMI *cl =0;
4551 //FindNearest3(y,z,roady,roadz,index);
4552 // for (Int_t i=Find(z-roadz); i<fN; i++) {
4553 for (Int_t i=iz1; i<iz2; i++) {
4554 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4555 if (c->GetZ() > z+roadz) break;
4556 if ( c->GetY()-y > roady ) continue;
4557 if ( y-c->GetY() > roady ) continue;
4558 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4559 if (maxdistance>distance) {
4560 maxdistance = distance;
4561 cl=c;
4562 index =i;
4563 //roady = TMath::Sqrt(maxdistance);
4564 }
4565 }
4566 return cl;
4567}
4568
4569
4570
4571AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4572{
4573 //-----------------------------------------------------------------------
4574 // Return the index of the nearest cluster in z y
4575 //-----------------------------------------------------------------------
4576 Float_t maxdistance = roady*roady + roadz*roadz;
4577 // Int_t iz = Int_t(z+255.);
4578 AliTPCclusterMI *cl =0;
4579 for (Int_t i=Find(z-roadz); i<fN; i++) {
4580 //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
4581 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4582 if (c->GetZ() > z+roadz) break;
4583 if ( c->GetY()-y > roady ) continue;
4584 if ( y-c->GetY() > roady ) continue;
4585 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4586 if (maxdistance>distance) {
4587 maxdistance = distance;
4588 cl=c;
4589 index =i;
4590 //roady = TMath::Sqrt(maxdistance);
4591 }
4592 }
4593 return cl;
4594}
1627d1c4 4595
4596
4597
4598
1c53abe2 4599AliTPCseed::AliTPCseed():AliTPCtrack(){
4600 //
4601 fRow=0;
4602 fRemoval =0;
91162307 4603 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4604 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
4605
1c53abe2 4606 fPoints = 0;
4607 fEPoints = 0;
4608 fNFoundable =0;
4609 fNShared =0;
91162307 4610 // fTrackPoints =0;
1c53abe2 4611 fRemoval = 0;
91162307 4612 fSort =0;
4613 fFirstPoint =0;
4614 fNoCluster =0;
4615 fBSigned = kFALSE;
4616 fSeed1 =-1;
4617 fSeed2 =-1;
1c53abe2 4618}
4619
4620AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
4621 fPoints = 0;
4622 fEPoints = 0;
4623 fNShared =0;
91162307 4624 // fTrackPoints =0;
1c53abe2 4625 fRemoval =0;
91162307 4626 fSort =0;
4627 for (Int_t i=0;i<160;i++) {
4628 fClusterPointer[i] = 0;
4629 Int_t index = t.GetClusterIndex(i);
4630 if (index>0) {
4631 SetClusterIndex2(i,index);
4632 }
4633 else{
4634 SetClusterIndex2(i,-3);
4635 }
4636 }
4637 fFirstPoint =0;
4638 fNoCluster =0;
4639 fBSigned = kFALSE;
4640 fSeed1 =-1;
4641 fSeed2 =-1;
1c53abe2 4642}
4643
4644AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
4645 fRow=0;
91162307 4646 for (Int_t i=0;i<160;i++) {
4647 fClusterPointer[i] = 0;
4648 Int_t index = t.GetClusterIndex(i);
4649 SetClusterIndex2(i,index);
4650 }
4651
1c53abe2 4652 fPoints = 0;
4653 fEPoints = 0;
4654 fNFoundable =0;
4655 fNShared =0;
91162307 4656 // fTrackPoints =0;
1c53abe2 4657 fRemoval =0;
91162307 4658 fSort = 0;
4659 fFirstPoint =0;
4660 fNoCluster =0;
4661 fBSigned = kFALSE;
4662 fSeed1 =-1;
4663 fSeed2 =-1;
1c53abe2 4664}
4665
4666AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
4667 Double_t xr, Double_t alpha):
4668 AliTPCtrack(index, xx, cc, xr, alpha) {
4669 //
4670 //
4671 fRow =0;
91162307 4672 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4673 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
1c53abe2 4674 fPoints = 0;
4675 fEPoints = 0;
4676 fNFoundable =0;
4677 fNShared = 0;
91162307 4678 // fTrackPoints =0;
1c53abe2 4679 fRemoval =0;
91162307 4680 fSort =0;
4681 fFirstPoint =0;
4682 // fHelixIn = new TClonesArray("AliHelix",0);
4683 //fHelixOut = new TClonesArray("AliHelix",0);
4684 fNoCluster =0;
4685 fBSigned = kFALSE;
4686 fSeed1 =-1;
4687 fSeed2 =-1;
1c53abe2 4688}
4689
4690AliTPCseed::~AliTPCseed(){
4691 if (fPoints) delete fPoints;
4692 fPoints =0;
91162307 4693 if (fEPoints) delete fEPoints;
1c53abe2 4694 fEPoints = 0;
91162307 4695 fNoCluster =0;
1c53abe2 4696}
4697
91162307 4698AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
1c53abe2 4699{
4700 //
4701 //
91162307 4702 return &fTrackPoints[i];
1c53abe2 4703}
4704
4705void AliTPCseed::RebuildSeed()
4706{
4707 //
4708 // rebuild seed to be ready for storing
91162307 4709 AliTPCclusterMI cldummy;
4710 cldummy.SetQ(0);
4711 AliTPCTrackPoint pdummy;
4712 pdummy.GetTPoint().fIsShared = 10;
1c53abe2 4713 for (Int_t i=0;i<160;i++){
91162307 4714 AliTPCclusterMI * cl0 = fClusterPointer[i];
4715 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
4716 if (cl0){
4717 trpoint->GetTPoint() = *(GetTrackPoint(i));
4718 trpoint->GetCPoint() = *cl0;
4719 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
4720 }
4721 else{
4722 *trpoint = pdummy;
4723 trpoint->GetCPoint()= cldummy;
4724 }
4725
4726 }
4727
4728}
4729
4730
4731Double_t AliTPCseed::GetDensityFirst(Int_t n)
4732{
4733 //
4734 //
4735 // return cluster for n rows bellow first point
4736 Int_t nfoundable = 1;
4737 Int_t nfound = 1;
4738 for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
4739 Int_t index = GetClusterIndex2(i);
4740 if (index!=-1) nfoundable++;
4741 if (index>0) nfound++;
1c53abe2 4742 }
91162307 4743 if (nfoundable<n) return 0;
4744 return Double_t(nfound)/Double_t(nfoundable);
4745
4746}
1c53abe2 4747
91162307 4748
4749void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
4750{
4751 // get cluster stat. on given region
4752 //
4753 found = 0;
4754 foundable = 0;
4755 shared =0;
4756 for (Int_t i=first;i<last; i++){
4757 Int_t index = GetClusterIndex2(i);
4758 if (index!=-1) foundable++;
4759 if (fClusterPointer[i]) {
4760 found++;
4761 }
4762 else
4763 continue;
4764
4765 if (fClusterPointer[i]->IsUsed(10)) {
4766 shared++;
4767 continue;
4768 }
4769 if (!plus2) continue; //take also neighborhoud
4770 //
4771 if ( (i>0) && fClusterPointer[i-1]){
4772 if (fClusterPointer[i-1]->IsUsed(10)) {
4773 shared++;
4774 continue;
4775 }
4776 }
4777 if ( fClusterPointer[i+1]){
4778 if (fClusterPointer[i+1]->IsUsed(10)) {
4779 shared++;
4780 continue;
4781 }
4782 }
4783
4784 }
4785 if (shared>found){
4786 printf("problem\n");
4787 }
1c53abe2 4788}
4789
4790//_____________________________________________________________________________
91162307 4791void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
1c53abe2 4792 //-----------------------------------------------------------------
4793 // This funtion calculates dE/dX within the "low" and "up" cuts.
4794 //-----------------------------------------------------------------
4795
4796 Float_t amp[200];
4797 Float_t angular[200];
4798 Float_t weight[200];
4799 Int_t index[200];
4800 //Int_t nc = 0;
4801 // TClonesArray & arr = *fPoints;
4802 Float_t meanlog = 100.;
4803
4804 Float_t mean[4] = {0,0,0,0};
4805 Float_t sigma[4] = {1000,1000,1000,1000};
4806 Int_t nc[4] = {0,0,0,0};
4807 Float_t norm[4] = {1000,1000,1000,1000};
4808 //
4809 //
4810 fNShared =0;
4811
4812 for (Int_t of =0; of<4; of++){
91162307 4813 for (Int_t i=of+i1;i<i2;i+=4)
1c53abe2 4814 {
91162307 4815 Int_t index = fIndex[i];
4816 if (index<0||index&0x8000) continue;
4817
1c53abe2 4818 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
91162307 4819 AliTPCTrackerPoint * point = GetTrackPoint(i);
4820 //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
4821 //AliTPCTrackerPoint * pointp = 0;
4822 //if (i<159) pointp = GetTrackPoint(i+1);
4823
1c53abe2 4824 if (point==0) continue;
91162307 4825 AliTPCclusterMI * cl = fClusterPointer[i];
4826 if (cl==0) continue;
4827 if (onlyused && (!cl->IsUsed(10))) continue;
4828 if (cl->IsUsed(11)) {
1c53abe2 4829 fNShared++;
4830 continue;
4831 }
91162307 4832 Int_t type = cl->GetType();
4833 //if (point->fIsShared){
4834 // fNShared++;
4835 // continue;
4836 //}
4837 //if (pointm)
4838 // if (pointm->fIsShared) continue;
4839 //if (pointp)
4840 // if (pointp->fIsShared) continue;
4841
4842 if (type<0) continue;
4843 //if (type>10) continue;
4844 //if (point->GetErrY()==0) continue;
4845 //if (point->GetErrZ()==0) continue;
4846
4847 //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
4848 //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
4849 //if ((ddy*ddy+ddz*ddz)>10) continue;
4850
4851
4852 // if (point->GetCPoint().GetMax()<5) continue;
4853 if (cl->GetMax()<5) continue;
4854 Float_t angley = point->GetAngleY();
4855 Float_t anglez = point->GetAngleZ();
4856
4857 Float_t rsigmay2 = point->GetSigmaY();
4858 Float_t rsigmaz2 = point->GetSigmaZ();
4859 /*
4860 Float_t ns = 1.;
4861 if (pointm){
4862 rsigmay += pointm->GetTPoint().GetSigmaY();
4863 rsigmaz += pointm->GetTPoint().GetSigmaZ();
4864 ns+=1.;
4865 }
4866 if (pointp){
4867 rsigmay += pointp->GetTPoint().GetSigmaY();
4868 rsigmaz += pointp->GetTPoint().GetSigmaZ();
4869 ns+=1.;
4870 }
4871 rsigmay/=ns;
4872 rsigmaz/=ns;
4873 */
4874
4875 Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
1c53abe2 4876
4877 Float_t ampc = 0; // normalization to the number of electrons
4878 if (i>64){
91162307 4879 // ampc = 1.*point->GetCPoint().GetMax();
4880 ampc = 1.*cl->GetMax();
1c53abe2 4881 //ampc = 1.*point->GetCPoint().GetQ();
4882 // AliTPCClusterPoint & p = point->GetCPoint();
4883 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
4884 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4885 //Float_t dz =
4886 // TMath::Abs( Int_t(iz) - iz + 0.5);
4887 //ampc *= 1.15*(1-0.3*dy);
4888 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 4889 // Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
4890 //ampc *=zfactor;
1c53abe2 4891 }
4892 else{
91162307 4893 //ampc = 1.0*point->GetCPoint().GetMax();
4894 ampc = 1.0*cl->GetMax();
1c53abe2 4895 //ampc = 1.0*point->GetCPoint().GetQ();
4896 //AliTPCClusterPoint & p = point->GetCPoint();
4897 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
4898 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4899 //Float_t dz =
4900 // TMath::Abs( Int_t(iz) - iz + 0.5);
4901
4902 //ampc *= 1.15*(1-0.3*dy);
4903 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 4904 // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
4905 //ampc *=zfactor;
1c53abe2 4906
4907 }
4908 ampc *= 2.0; // put mean value to channel 50
4909 //ampc *= 0.58; // put mean value to channel 50
4910 Float_t w = 1.;
4911 // if (type>0) w = 1./(type/2.-0.5);
91162307 4912 // Float_t z = TMath::Abs(cl->GetZ());
1c53abe2 4913 if (i<64) {
4914 ampc /= 0.6;
1627d1c4 4915 //ampc /= (1+0.0008*z);
4916 } else
4917 if (i>128){
4918 ampc /=1.5;
4919 //ampc /= (1+0.0008*z);
4920 }else{
4921 //ampc /= (1+0.0008*z);
4922 }
4923
1c53abe2 4924 if (type<0) { //amp at the border - lower weight
4925 // w*= 2.;
1627d1c4 4926
1c53abe2 4927 continue;
4928 }
4929 if (rsigma>1.5) ampc/=1.3; // if big backround
4930 amp[nc[of]] = ampc;
4931 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
4932 weight[nc[of]] = w;
4933 nc[of]++;
4934 }
4935
4936 TMath::Sort(nc[of],amp,index,kFALSE);
4937 Float_t sumamp=0;
4938 Float_t sumamp2=0;
4939 Float_t sumw=0;
4940 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
91162307 4941 meanlog = 50;
1c53abe2 4942 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
4943 Float_t ampl = amp[index[i]]/angular[index[i]];
4944 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
4945 //
4946 sumw += weight[index[i]];
4947 sumamp += weight[index[i]]*ampl;
4948 sumamp2 += weight[index[i]]*ampl*ampl;
4949 norm[of] += angular[index[i]]*weight[index[i]];
4950 }
4951 if (sumw<1){
4952 SetdEdx(0);
4953 }
4954 else {
4955 norm[of] /= sumw;
4956 mean[of] = sumamp/sumw;
4957 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
4958 if (sigma[of]>0.1)
4959 sigma[of] = TMath::Sqrt(sigma[of]);
4960 else
4961 sigma[of] = 1000;
4962
4963 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
4964 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
4965 //mean *=(1-0.1*(norm-1.));
4966 }
4967 }
4968
4969 Float_t dedx =0;
4970 fSdEdx =0;
4971 fMAngular =0;
4972 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
4973 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
4974
4975
4976 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
4977 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
4978
4979 Int_t norm2 = 0;
4980 Int_t norm3 = 0;
4981 for (Int_t i =0;i<4;i++){
4982 if (nc[i]>2&&nc[i]<1000){
4983 dedx += mean[i] *nc[i];
4984 fSdEdx += sigma[i]*(nc[i]-2);
4985 fMAngular += norm[i] *nc[i];
4986 norm2 += nc[i];
4987 norm3 += nc[i]-2;
4988 }
4989 fDEDX[i] = mean[i];
4990 fSDEDX[i] = sigma[i];
4991 fNCDEDX[i]= nc[i];
4992 }
4993
4994 if (norm3>0){
4995 dedx /=norm2;
4996 fSdEdx /=norm3;
4997 fMAngular/=norm2;
4998 }
4999 else{
5000 SetdEdx(0);
5001 return;
5002 }
5003 // Float_t dedx1 =dedx;
91162307 5004 /*
1c53abe2 5005 dedx =0;
5006 for (Int_t i =0;i<4;i++){
5007 if (nc[i]>2&&nc[i]<1000){
5008 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
5009 dedx += mean[i] *nc[i];
5010 }
5011 fDEDX[i] = mean[i];
5012 }
5013 dedx /= norm2;
91162307 5014 */
1c53abe2 5015
5016
5017 SetdEdx(dedx);
5018
5019 //mi deDX
5020
5021
5022
5023 //Very rough PID
5024 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
5025
5026 if (p<0.6) {
5027 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5028 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
5029 SetMass(0.93827); return;
5030 }
5031
5032 if (p<1.2) {
5033 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5034 SetMass(0.93827); return;
5035 }
5036
5037 SetMass(0.13957); return;
5038
5039}
5040
5041
91162307 5042
5043/*
5044
5045
5046
5047void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
5048 //-----------------------------------------------------------------
5049 // This funtion calculates dE/dX within the "low" and "up" cuts.
5050 //-----------------------------------------------------------------
5051
5052 Float_t amp[200];
5053 Float_t angular[200];
5054 Float_t weight[200];
5055 Int_t index[200];
5056 Bool_t inlimit[200];
5057 for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
5058 for (Int_t i=0;i<200;i++) amp[i]=10000;
5059 for (Int_t i=0;i<200;i++) angular[i]= 1;;
5060
5061
5062 //
5063 Float_t meanlog = 100.;
5064 Int_t indexde[4]={0,64,128,160};
5065
5066 Float_t amean =0;
5067 Float_t asigma =0;
5068 Float_t anc =0;
5069 Float_t anorm =0;
5070
5071 Float_t mean[4] = {0,0,0,0};
5072 Float_t sigma[4] = {1000,1000,1000,1000};
5073 Int_t nc[4] = {0,0,0,0};
5074 Float_t norm[4] = {1000,1000,1000,1000};
5075 //
5076 //
5077 fNShared =0;
5078
5079 // for (Int_t of =0; of<3; of++){
5080 // for (Int_t i=indexde[of];i<indexde[of+1];i++)
5081 for (Int_t i =0; i<160;i++)
5082 {
5083 AliTPCTrackPoint * point = GetTrackPoint(i);
5084 if (point==0) continue;
5085 if (point->fIsShared){
5086 fNShared++;
5087 continue;
5088 }
5089 Int_t type = point->GetCPoint().GetType();
5090 if (type<0) continue;
5091 if (point->GetCPoint().GetMax()<5) continue;
5092 Float_t angley = point->GetTPoint().GetAngleY();
5093 Float_t anglez = point->GetTPoint().GetAngleZ();
5094 Float_t rsigmay = point->GetCPoint().GetSigmaY();
5095 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
5096 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
5097
5098 Float_t ampc = 0; // normalization to the number of electrons
5099 if (i>64){
5100 ampc = point->GetCPoint().GetMax();
5101 }
5102 else{
5103 ampc = point->GetCPoint().GetMax();
5104 }
5105 ampc *= 2.0; // put mean value to channel 50
5106 // ampc *= 0.565; // put mean value to channel 50
5107
5108 Float_t w = 1.;
5109 Float_t z = TMath::Abs(point->GetCPoint().GetZ());
5110 if (i<64) {
5111 ampc /= 0.63;
5112 } else
5113 if (i>128){
5114 ampc /=1.51;
5115 }
5116 if (type<0) { //amp at the border - lower weight
5117 continue;
5118 }
5119 if (rsigma>1.5) ampc/=1.3; // if big backround
5120 angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5121 amp[i] = ampc/angular[i];
5122 weight[i] = w;
5123 anc++;
5124 }
5125
5126 TMath::Sort(159,amp,index,kFALSE);
5127 for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
5128 inlimit[index[i]] = kTRUE; // take all clusters
5129 }
5130
5131 // meanlog = amp[index[Int_t(anc*0.3)]];
5132 meanlog =10000.;
5133 for (Int_t of =0; of<3; of++){
5134 Float_t sumamp=0;
5135 Float_t sumamp2=0;
5136 Float_t sumw=0;
5137 for (Int_t i=indexde[of];i<indexde[of+1];i++)
5138 {
5139 if (inlimit[i]==kFALSE) continue;
5140 Float_t ampl = amp[i];
5141 ///angular[i];
5142 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
5143 //
5144 sumw += weight[i];
5145 sumamp += weight[i]*ampl;
5146 sumamp2 += weight[i]*ampl*ampl;
5147 norm[of] += angular[i]*weight[i];
5148 nc[of]++;
5149 }
5150 if (sumw<1){
5151 SetdEdx(0);
5152 }
5153 else {
5154 norm[of] /= sumw;
5155 mean[of] = sumamp/sumw;
5156 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5157 if (sigma[of]>0.1)
5158 sigma[of] = TMath::Sqrt(sigma[of]);
5159 else
5160 sigma[of] = 1000;
5161 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5162 }
5163 }
5164
5165 Float_t dedx =0;
5166 fSdEdx =0;
5167 fMAngular =0;
5168 //
5169 Int_t norm2 = 0;
5170 Int_t norm3 = 0;
5171 Float_t www[3] = {12.,14.,17.};
5172 //Float_t www[3] = {1.,1.,1.};
5173
5174 for (Int_t i =0;i<3;i++){
5175 if (nc[i]>2&&nc[i]<1000){
5176 dedx += mean[i] *nc[i]*www[i]/sigma[i];
5177 fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
5178 fMAngular += norm[i] *nc[i];
5179 norm2 += nc[i]*www[i]/sigma[i];
5180 norm3 += (nc[i]-2)*www[i]/sigma[i];
5181 }
5182 fDEDX[i] = mean[i];
5183 fSDEDX[i] = sigma[i];
5184 fNCDEDX[i]= nc[i];
5185 }
5186
5187 if (norm3>0){
5188 dedx /=norm2;
5189 fSdEdx /=norm3;
5190 fMAngular/=norm2;
5191 }
5192 else{
5193 SetdEdx(0);
5194 return;
5195 }
5196 // Float_t dedx1 =dedx;
5197
5198 dedx =0;
5199 Float_t norm4 = 0;
5200 for (Int_t i =0;i<3;i++){
5201 if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
5202 //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
5203 dedx += mean[i] *(nc[i])/(sigma[i]);
5204 norm4 += (nc[i])/(sigma[i]);
5205 }
5206 fDEDX[i] = mean[i];
5207 }
5208 if (norm4>0) dedx /= norm4;
5209
5210
5211
5212 SetdEdx(dedx);
5213
5214 //mi deDX
5215
5216}
5217
5218*/