Modified for CAF
[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
1c53abe2 16
17//-------------------------------------------------------
18// Implementation of the TPC tracker
19//
20// Origin: Marian Ivanov Marian.Ivanov@cern.ch
21//
34acb742 22// AliTPC parallel tracker
1c53abe2 23//-------------------------------------------------------
47966a6d 24
25
26/* $Id$ */
27
cc5e9db0 28#include "Riostream.h"
6d171107 29#include <TClonesArray.h>
30#include <TFile.h>
31#include <TObjArray.h>
32#include <TTree.h>
a3232aae 33#include "AliLog.h"
1c53abe2 34
47966a6d 35#include "AliComplexCluster.h"
91162307 36#include "AliESD.h"
6c94f330 37#include "AliKink.h"
38#include "AliV0.h"
91162307 39#include "AliHelix.h"
91162307 40#include "AliRunLoader.h"
6d171107 41#include "AliTPCClustersRow.h"
42#include "AliTPCParam.h"
9996a03b 43#include "AliTPCReconstructor.h"
6d171107 44#include "AliTPCclusterMI.h"
45#include "AliTPCpolyTrack.h"
81e97e0d 46#include "AliTPCreco.h"
47#include "AliTPCseed.h"
b67e07dc 48#include "AliTPCtrackerMI.h"
6d171107 49#include "TStopwatch.h"
81e97e0d 50#include "AliTPCReconstructor.h"
81e97e0d 51#include "AliPID.h"
eea478d3 52#include "TTreeStream.h"
5d837844 53#include "AliAlignObj.h"
54#include "AliTrackPointArray.h"
55
6d171107 56//
c9427e08 57
91162307 58ClassImp(AliTPCtrackerMI)
c9427e08 59
60
b67e07dc 61class AliTPCFastMath {
91162307 62public:
b67e07dc 63 AliTPCFastMath();
91162307 64 static Double_t FastAsin(Double_t x);
65 private:
b67e07dc 66 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
91162307 67};
c9427e08 68
b67e07dc 69Double_t AliTPCFastMath::fgFastAsin[20000];
2274b54b 70AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
c9427e08 71
b67e07dc 72AliTPCFastMath::AliTPCFastMath(){
73 //
74 // initialized lookup table;
91162307 75 for (Int_t i=0;i<10000;i++){
76 fgFastAsin[2*i] = TMath::ASin(i/10000.);
77 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
78 }
c9427e08 79}
80
b67e07dc 81Double_t AliTPCFastMath::FastAsin(Double_t x){
82 //
83 // return asin using lookup table
91162307 84 if (x>0){
85 Int_t index = int(x*10000);
86 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
87 }
88 x*=-1;
89 Int_t index = int(x*10000);
90 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
1c53abe2 91}
92
93
94
95
91162307 96Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
b67e07dc 97 //
98 //update track information using current cluster - track->fCurrentCluster
99
1c53abe2 100
91162307 101 AliTPCclusterMI* c =track->fCurrentCluster;
102 if (accept>0) track->fCurrentClusterIndex1 |=0x8000; //sign not accepted clusters
c9427e08 103
91162307 104 UInt_t i = track->fCurrentClusterIndex1;
1c53abe2 105
106 Int_t sec=(i&0xff000000)>>24;
91162307 107 //Int_t row = (i&0x00ff0000)>>16;
1c53abe2 108 track->fRow=(i&0x00ff0000)>>16;
109 track->fSector = sec;
110 // Int_t index = i&0xFFFF;
111 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
d26d9159 112 track->SetClusterIndex2(track->fRow, i);
91162307 113 //track->fFirstPoint = row;
114 //if ( track->fLastPoint<row) track->fLastPoint =row;
115 // if (track->fRow<0 || track->fRow>160) {
116 // printf("problem\n");
117 //}
118 if (track->fFirstPoint>track->fRow)
119 track->fFirstPoint = track->fRow;
120 if (track->fLastPoint<track->fRow)
121 track->fLastPoint = track->fRow;
122
123
124 track->fClusterPointer[track->fRow] = c;
1c53abe2 125 //
126
3f82c4f2 127 Double_t angle2 = track->GetSnp()*track->GetSnp();
1c53abe2 128 //
129 //SET NEW Track Point
130 //
85e2b57d 131 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
91162307 132 {
85e2b57d 133 angle2 = TMath::Sqrt(angle2/(1-angle2));
91162307 134 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow));
1c53abe2 135 //
91162307 136 point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
137 point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
138 point.SetErrY(sqrt(track->fErrorY2));
139 point.SetErrZ(sqrt(track->fErrorZ2));
1c53abe2 140 //
91162307 141 point.SetX(track->GetX());
142 point.SetY(track->GetY());
143 point.SetZ(track->GetZ());
144 point.SetAngleY(angle2);
145 point.SetAngleZ(track->GetTgl());
146 if (point.fIsShared){
147 track->fErrorY2 *= 4;
148 track->fErrorZ2 *= 4;
149 }
150 }
151
152 Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
153 //
154 track->fErrorY2 *= 1.3;
155 track->fErrorY2 += 0.01;
156 track->fErrorZ2 *= 1.3;
157 track->fErrorZ2 += 0.005;
158 //}
159 if (accept>0) return 0;
160 if (track->GetNumberOfClusters()%20==0){
161 // if (track->fHelixIn){
162 // TClonesArray & larr = *(track->fHelixIn);
163 // Int_t ihelix = larr.GetEntriesFast();
164 // new(larr[ihelix]) AliHelix(*track) ;
165 //}
1c53abe2 166 }
91162307 167 track->fNoCluster =0;
168 return track->Update(c,chi2,i);
169}
170
171
172
173Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
174 Float_t cory, Float_t corz)
175{
1c53abe2 176 //
91162307 177 // decide according desired precision to accept given
178 // cluster for tracking
179 Double_t sy2=ErrY2(seed,cluster)*cory;
180 Double_t sz2=ErrZ2(seed,cluster)*corz;
181 //sy2=ErrY2(seed,cluster)*cory;
182 //sz2=ErrZ2(seed,cluster)*cory;
1c53abe2 183
91162307 184 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
185 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
186
187 Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
188 (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
189 Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
190 (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
191
192 Double_t rdistance2 = rdistancey2+rdistancez2;
193 //Int_t accept =0;
1c53abe2 194
91162307 195 if (rdistance2>16) return 3;
196
197
198 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
199 return 2; //suspisiouce - will be changed
200
201 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
202 // strict cut on overlaped cluster
203 return 2; //suspisiouce - will be changed
204
205 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
206 && cluster->GetType()<0){
207 seed->fNFoundable--;
208 return 2;
1c53abe2 209 }
91162307 210 return 0;
211}
212
213
1c53abe2 214
1c53abe2 215
1c53abe2 216//_____________________________________________________________________________
f8aae377 217AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
1c53abe2 218AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
219{
220 //---------------------------------------------------------------------
221 // The main TPC tracker constructor
222 //---------------------------------------------------------------------
223 fInnerSec=new AliTPCSector[fkNIS];
224 fOuterSec=new AliTPCSector[fkNOS];
91162307 225
1c53abe2 226 Int_t i;
227 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
228 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
229
230 fN=0; fSectors=0;
231
1c53abe2 232 fSeeds=0;
233 fNtracks = 0;
91162307 234 fParam = par;
235 Int_t nrowlow = par->GetNRowLow();
236 Int_t nrowup = par->GetNRowUp();
237
238
239 for (Int_t i=0;i<nrowlow;i++){
240 fXRow[i] = par->GetPadRowRadiiLow(i);
241 fPadLength[i]= par->GetPadPitchLength(0,i);
242 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
243 }
244
245
246 for (Int_t i=0;i<nrowup;i++){
247 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
248 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
249 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
250 }
251 fSeeds=0;
252 //
253 fInput = 0;
254 fOutput = 0;
255 fSeedTree = 0;
256 fTreeDebug =0;
257 fNewIO =0;
258 fDebug =0;
259 fEvent =0;
81e97e0d 260 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
1c53abe2 261}
2fc0c115 262//________________________________________________________________________
58251ea0 263AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
264 AliTracker(t),
265 fkNIS(t.fkNIS),
266 fkNOS(t.fkNOS)
267{
2fc0c115 268 //------------------------------------
269 // dummy copy constructor
270 //------------------------------------------------------------------
271}
272AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
273 //------------------------------
274 // dummy
275 //--------------------------------------------------------------
276 return *this;
277}
1c53abe2 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 }
81e97e0d 289 if (fDebugStreamer) delete fDebugStreamer;
1c53abe2 290}
291
91162307 292void AliTPCtrackerMI::SetIO()
293{
1c53abe2 294 //
91162307 295 fNewIO = kTRUE;
e191bb57 296 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
d26d9159 297
e191bb57 298 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
d26d9159 299 if (fOutput){
300 AliTPCtrack *iotrack= new AliTPCtrack;
301 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
302 delete iotrack;
303 }
91162307 304}
1c53abe2 305
d26d9159 306
91162307 307void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
308{
1c53abe2 309
91162307 310 // set input
311 fNewIO = kFALSE;
312 fInput = 0;
313 fOutput = 0;
314 fSeedTree = 0;
315 fTreeDebug =0;
316 fInput = input;
317 if (input==0){
318 return;
319 }
320 //set output
321 fOutput = output;
322 if (output){
323 AliTPCtrack *iotrack= new AliTPCtrack;
324 // iotrack->fHelixIn = new TClonesArray("AliHelix");
325 //iotrack->fHelixOut = new TClonesArray("AliHelix");
326 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
327 delete iotrack;
328 }
329 if (output && (fDebug&2)){
330 //write the full seed information if specified in debug mode
331 //
332 fSeedTree = new TTree("Seeds","Seeds");
333 AliTPCseed * vseed = new AliTPCseed;
334 //
335 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
336 arrtr->ExpandCreateFast(160);
337 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
338 //
339 vseed->fPoints = arrtr;
340 vseed->fEPoints = arre;
341 // vseed->fClusterPoints = arrcl;
342 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
343 delete arrtr;
344 delete arre;
345 fTreeDebug = new TTree("trackDebug","trackDebug");
346 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
347 fTreeDebug->Branch("debug",&arrd,32000,99);
1c53abe2 348 }
1c53abe2 349
1c53abe2 350
91162307 351 //set ESD event
352 fEvent = event;
353}
1c53abe2 354
d26d9159 355void AliTPCtrackerMI::FillESD(TObjArray* arr)
91162307 356{
47966a6d 357 //
358 //
359 //fill esds using updated tracks
91162307 360 if (fEvent){
361 // write tracks to the event
362 // store index of the track
d26d9159 363 Int_t nseed=arr->GetEntriesFast();
51ad6848 364 //FindKinks(arr,fEvent);
91162307 365 for (Int_t i=0; i<nseed; i++) {
d26d9159 366 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 367 if (!pt) continue;
51ad6848 368 pt->UpdatePoints();
eea478d3 369 // pt->PropagateTo(fParam->GetInnerRadiusLow());
370 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
371 pt->PropagateTo(fParam->GetInnerRadiusLow());
372 }
51ad6848 373
374 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
375 AliESDtrack iotrack;
376 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
377 iotrack.SetTPCPoints(pt->GetPoints());
378 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 379 iotrack.SetV0Indexes(pt->GetV0Indexes());
380 // iotrack.SetTPCpid(pt->fTPCr);
51ad6848 381 //iotrack.SetTPCindex(i);
382 fEvent->AddTrack(&iotrack);
383 continue;
384 }
385
4d158c36 386 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
d26d9159 387 AliESDtrack iotrack;
51ad6848 388 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
389 iotrack.SetTPCPoints(pt->GetPoints());
d26d9159 390 //iotrack.SetTPCindex(i);
51ad6848 391 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 392 iotrack.SetV0Indexes(pt->GetV0Indexes());
393 // iotrack.SetTPCpid(pt->fTPCr);
d26d9159 394 fEvent->AddTrack(&iotrack);
a42a6bae 395 continue;
396 }
51ad6848 397 //
398 // short tracks - maybe decays
399
a42a6bae 400 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
401 Int_t found,foundable,shared;
402 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
403 if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){
404 AliESDtrack iotrack;
405 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
406 //iotrack.SetTPCindex(i);
51ad6848 407 iotrack.SetTPCPoints(pt->GetPoints());
408 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 409 iotrack.SetV0Indexes(pt->GetV0Indexes());
410 //iotrack.SetTPCpid(pt->fTPCr);
a42a6bae 411 fEvent->AddTrack(&iotrack);
412 continue;
413 }
414 }
415
416 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) {
417 Int_t found,foundable,shared;
418 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
419 if (found<20) continue;
420 if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
421 //
422 AliESDtrack iotrack;
423 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
51ad6848 424 iotrack.SetTPCPoints(pt->GetPoints());
425 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 426 iotrack.SetV0Indexes(pt->GetV0Indexes());
427 //iotrack.SetTPCpid(pt->fTPCr);
51ad6848 428 //iotrack.SetTPCindex(i);
429 fEvent->AddTrack(&iotrack);
430 continue;
431 }
432 // short tracks - secondaties
433 //
434 if ( (pt->GetNumberOfClusters()>30) ) {
435 Int_t found,foundable,shared;
436 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
437 if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
438 AliESDtrack iotrack;
439 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
440 iotrack.SetTPCPoints(pt->GetPoints());
441 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 442 iotrack.SetV0Indexes(pt->GetV0Indexes());
443 //iotrack.SetTPCpid(pt->fTPCr);
51ad6848 444 //iotrack.SetTPCindex(i);
445 fEvent->AddTrack(&iotrack);
446 continue;
447 }
448 }
449
450 if ( (pt->GetNumberOfClusters()>15)) {
451 Int_t found,foundable,shared;
452 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
453 if (found<15) continue;
e7eb17e4 454 if (foundable<=0) continue;
51ad6848 455 if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
456 if (float(found)/float(foundable)<0.8) continue;
457 //
458 AliESDtrack iotrack;
459 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
460 iotrack.SetTPCPoints(pt->GetPoints());
461 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 462 iotrack.SetV0Indexes(pt->GetV0Indexes());
463 // iotrack.SetTPCpid(pt->fTPCr);
a42a6bae 464 //iotrack.SetTPCindex(i);
465 fEvent->AddTrack(&iotrack);
466 continue;
467 }
91162307 468 }
1c53abe2 469 }
51ad6848 470 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
d26d9159 471}
472
47966a6d 473void AliTPCtrackerMI::WriteTracks(TTree * tree)
474{
d26d9159 475 //
47966a6d 476 // write tracks from seed array to selected tree
d26d9159 477 //
478 fOutput = tree;
479 if (fOutput){
480 AliTPCtrack *iotrack= new AliTPCtrack;
481 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
482 }
483 WriteTracks();
484}
485
486void AliTPCtrackerMI::WriteTracks()
487{
488 //
489 // write tracks to the given output tree -
490 // output specified with SetIO routine
491 if (!fSeeds) return;
492 if (!fOutput){
493 SetIO();
494 }
1c53abe2 495
91162307 496 if (fOutput){
497 AliTPCtrack *iotrack= 0;
498 Int_t nseed=fSeeds->GetEntriesFast();
982aff31 499 //for (Int_t i=0; i<nseed; i++) {
500 // iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
501 // if (iotrack) break;
502 //}
91162307 503 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
504 TBranch * br = fOutput->GetBranch("tracks");
505 br->SetAddress(&iotrack);
506 //
507 for (Int_t i=0; i<nseed; i++) {
508 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
509 if (!pt) continue;
982aff31 510 AliTPCtrack * track = new AliTPCtrack(*pt);
511 iotrack = track;
91162307 512 pt->fLab2 =i;
513 // br->SetAddress(&iotrack);
514 fOutput->Fill();
982aff31 515 delete track;
91162307 516 iotrack =0;
517 }
982aff31 518 //fOutput->GetDirectory()->cd();
519 //fOutput->Write();
91162307 520 }
d26d9159 521 // delete iotrack;
522 //
91162307 523 if (fSeedTree){
524 //write the full seed information if specified in debug mode
525
526 AliTPCseed * vseed = new AliTPCseed;
527 //
528 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
529 arrtr->ExpandCreateFast(160);
530 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
531 //arrcl->ExpandCreateFast(160);
532 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
533 //
534 vseed->fPoints = arrtr;
535 vseed->fEPoints = arre;
536 // vseed->fClusterPoints = arrcl;
537 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
538 TBranch * brseed = fSeedTree->GetBranch("seeds");
539
540 Int_t nseed=fSeeds->GetEntriesFast();
541
542 for (Int_t i=0; i<nseed; i++) {
543 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
544 if (!pt) continue;
545 pt->fPoints = arrtr;
546 // pt->fClusterPoints = arrcl;
547 pt->fEPoints = arre;
548 pt->RebuildSeed();
549 vseed = pt;
550 brseed->SetAddress(&vseed);
551 fSeedTree->Fill();
552 pt->fPoints = 0;
553 pt->fEPoints = 0;
554 // pt->fClusterPoints = 0;
555 }
556 fSeedTree->Write();
557 if (fTreeDebug) fTreeDebug->Write();
558 }
1c53abe2 559
91162307 560}
1c53abe2 561
1c53abe2 562
1c53abe2 563
564
91162307 565Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
566 //
567 //
568 //seed->SetErrorY2(0.1);
569 //return 0.1;
570 //calculate look-up table at the beginning
571 static Bool_t ginit = kFALSE;
572 static Float_t gnoise1,gnoise2,gnoise3;
573 static Float_t ggg1[10000];
574 static Float_t ggg2[10000];
575 static Float_t ggg3[10000];
576 static Float_t glandau1[10000];
577 static Float_t glandau2[10000];
578 static Float_t glandau3[10000];
579 //
580 static Float_t gcor01[500];
581 static Float_t gcor02[500];
582 static Float_t gcorp[500];
583 //
1c53abe2 584
91162307 585 //
586 if (ginit==kFALSE){
587 for (Int_t i=1;i<500;i++){
588 Float_t rsigma = float(i)/100.;
589 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
590 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
591 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
592 }
1c53abe2 593
91162307 594 //
595 for (Int_t i=3;i<10000;i++){
596 //
597 //
598 // inner sector
599 Float_t amp = float(i);
600 Float_t padlength =0.75;
601 gnoise1 = 0.0004/padlength;
602 Float_t nel = 0.268*amp;
603 Float_t nprim = 0.155*amp;
604 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
605 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
606 if (glandau1[i]>1) glandau1[i]=1;
607 glandau1[i]*=padlength*padlength/12.;
608 //
609 // outer short
610 padlength =1.;
611 gnoise2 = 0.0004/padlength;
612 nel = 0.3*amp;
613 nprim = 0.133*amp;
614 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
615 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
616 if (glandau2[i]>1) glandau2[i]=1;
617 glandau2[i]*=padlength*padlength/12.;
618 //
619 //
620 // outer long
621 padlength =1.5;
622 gnoise3 = 0.0004/padlength;
623 nel = 0.3*amp;
624 nprim = 0.133*amp;
625 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
626 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
627 if (glandau3[i]>1) glandau3[i]=1;
628 glandau3[i]*=padlength*padlength/12.;
629 //
630 }
631 ginit = kTRUE;
632 }
1c53abe2 633 //
634 //
91162307 635 //
636 Int_t amp = int(TMath::Abs(cl->GetQ()));
637 if (amp>9999) {
638 seed->SetErrorY2(1.);
639 return 1.;
640 }
1c53abe2 641 Float_t snoise2;
c9427e08 642 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
91162307 643 Int_t ctype = cl->GetType();
644 Float_t padlength= GetPadPitchLength(seed->fRow);
3f82c4f2 645 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
91162307 646 angle2 = angle2/(1-angle2);
647 //
648 //cluster "quality"
649 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
650 Float_t res;
1c53abe2 651 //
1c53abe2 652 if (fSectors==fInnerSec){
91162307 653 snoise2 = gnoise1;
654 res = ggg1[amp]*z+glandau1[amp]*angle2;
655 if (ctype==0) res *= gcor01[rsigmay];
656 if ((ctype>0)){
657 res+=0.002;
658 res*= gcorp[rsigmay];
659 }
1c53abe2 660 }
661 else {
91162307 662 if (padlength<1.1){
663 snoise2 = gnoise2;
664 res = ggg2[amp]*z+glandau2[amp]*angle2;
665 if (ctype==0) res *= gcor02[rsigmay];
666 if ((ctype>0)){
667 res+=0.002;
668 res*= gcorp[rsigmay];
669 }
670 }
671 else{
672 snoise2 = gnoise3;
673 res = ggg3[amp]*z+glandau3[amp]*angle2;
674 if (ctype==0) res *= gcor02[rsigmay];
675 if ((ctype>0)){
676 res+=0.002;
677 res*= gcorp[rsigmay];
678 }
679 }
680 }
1c53abe2 681
91162307 682 if (ctype<0){
683 res+=0.005;
684 res*=2.4; // overestimate error 2 times
685 }
686 res+= snoise2;
687
1c53abe2 688 if (res<2*snoise2)
91162307 689 res = 2*snoise2;
690
691 seed->SetErrorY2(res);
1c53abe2 692 return res;
1c53abe2 693
694
91162307 695}
c9427e08 696
697
698
91162307 699Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
700 //
701 //
702 //seed->SetErrorY2(0.1);
703 //return 0.1;
704 //calculate look-up table at the beginning
705 static Bool_t ginit = kFALSE;
706 static Float_t gnoise1,gnoise2,gnoise3;
707 static Float_t ggg1[10000];
708 static Float_t ggg2[10000];
709 static Float_t ggg3[10000];
710 static Float_t glandau1[10000];
711 static Float_t glandau2[10000];
712 static Float_t glandau3[10000];
713 //
714 static Float_t gcor01[1000];
715 static Float_t gcor02[1000];
716 static Float_t gcorp[1000];
1627d1c4 717 //
1627d1c4 718
91162307 719 //
720 if (ginit==kFALSE){
721 for (Int_t i=1;i<1000;i++){
722 Float_t rsigma = float(i)/100.;
723 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
724 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
725 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
726 }
1c53abe2 727
91162307 728 //
729 for (Int_t i=3;i<10000;i++){
730 //
731 //
732 // inner sector
733 Float_t amp = float(i);
734 Float_t padlength =0.75;
735 gnoise1 = 0.0004/padlength;
736 Float_t nel = 0.268*amp;
737 Float_t nprim = 0.155*amp;
738 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
739 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
740 if (glandau1[i]>1) glandau1[i]=1;
741 glandau1[i]*=padlength*padlength/12.;
742 //
743 // outer short
744 padlength =1.;
745 gnoise2 = 0.0004/padlength;
746 nel = 0.3*amp;
747 nprim = 0.133*amp;
748 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
749 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
750 if (glandau2[i]>1) glandau2[i]=1;
751 glandau2[i]*=padlength*padlength/12.;
752 //
753 //
754 // outer long
755 padlength =1.5;
756 gnoise3 = 0.0004/padlength;
757 nel = 0.3*amp;
758 nprim = 0.133*amp;
759 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
760 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
761 if (glandau3[i]>1) glandau3[i]=1;
762 glandau3[i]*=padlength*padlength/12.;
763 //
764 }
765 ginit = kTRUE;
766 }
767 //
768 //
769 //
770 Int_t amp = int(TMath::Abs(cl->GetQ()));
771 if (amp>9999) {
772 seed->SetErrorY2(1.);
773 return 1.;
774 }
775 Float_t snoise2;
776 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
777 Int_t ctype = cl->GetType();
778 Float_t padlength= GetPadPitchLength(seed->fRow);
779 //
3f82c4f2 780 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
91162307 781 // if (angle2<0.6) angle2 = 0.6;
782 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
783 //
784 //cluster "quality"
785 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
786 Float_t res;
787 //
788 if (fSectors==fInnerSec){
789 snoise2 = gnoise1;
790 res = ggg1[amp]*z+glandau1[amp]*angle2;
791 if (ctype==0) res *= gcor01[rsigmaz];
792 if ((ctype>0)){
793 res+=0.002;
794 res*= gcorp[rsigmaz];
795 }
796 }
797 else {
798 if (padlength<1.1){
799 snoise2 = gnoise2;
800 res = ggg2[amp]*z+glandau2[amp]*angle2;
801 if (ctype==0) res *= gcor02[rsigmaz];
802 if ((ctype>0)){
803 res+=0.002;
804 res*= gcorp[rsigmaz];
805 }
806 }
807 else{
808 snoise2 = gnoise3;
809 res = ggg3[amp]*z+glandau3[amp]*angle2;
810 if (ctype==0) res *= gcor02[rsigmaz];
811 if ((ctype>0)){
812 res+=0.002;
813 res*= gcorp[rsigmaz];
814 }
815 }
816 }
817
818 if (ctype<0){
819 res+=0.002;
820 res*=1.3;
821 }
822 if ((ctype<0) &&amp<70){
823 res+=0.002;
824 res*=1.3;
825 }
826 res += snoise2;
827 if (res<2*snoise2)
828 res = 2*snoise2;
829 if (res>3) res =3;
830 seed->SetErrorZ2(res);
831 return res;
832}
833
834
835
836/*
837Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
838 //
839 //
840 //seed->SetErrorZ2(0.1);
841 //return 0.1;
842
843 Float_t snoise2;
844 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
845 //
846 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
847 Int_t ctype = cl->GetType();
848 Float_t amp = TMath::Abs(cl->GetQ());
849
850 Float_t nel;
851 Float_t nprim;
852 //
853 Float_t landau=2 ; //landau fluctuation part
854 Float_t gg=2; // gg fluctuation part
855 Float_t padlength= GetPadPitchLength(seed->GetX());
856
857 if (fSectors==fInnerSec){
858 snoise2 = 0.0004/padlength;
859 nel = 0.268*amp;
860 nprim = 0.155*amp;
861 gg = (2+0.001*nel/(padlength*padlength))/nel;
862 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
863 if (landau>1) landau=1;
864 }
865 else {
866 snoise2 = 0.0004/padlength;
867 nel = 0.3*amp;
868 nprim = 0.133*amp;
869 gg = (2+0.0008*nel/(padlength*padlength))/nel;
870 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
871 if (landau>1) landau=1;
872 }
873 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
874
875 //
876 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
877 angle2 = TMath::Sqrt((1-angle2));
878 if (angle2<0.6) angle2 = 0.6;
879 //angle2 = 1;
880
881 Float_t angle = seed->GetTgl()/angle2;
882 Float_t angular = landau*angle*angle*padlength*padlength/12.;
883 Float_t res = sdiff + angular;
884
885
886 if ((ctype==0) && (fSectors ==fOuterSec))
887 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
888
889 if ((ctype==0) && (fSectors ==fInnerSec))
890 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
891
892 if ((ctype>0)){
893 res+=0.005;
894 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
895 }
896 if (ctype<0){
897 res+=0.002;
898 res*=1.3;
899 }
900 if ((ctype<0) &&amp<70){
901 res+=0.002;
902 res*=1.3;
903 }
904 res += snoise2;
905 if (res<2*snoise2)
906 res = 2*snoise2;
907
908 seed->SetErrorZ2(res);
909 return res;
910}
911*/
912
913
914
1c53abe2 915
c9427e08 916void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
917{
918 //rotate to track "local coordinata
919 Float_t x = seed->GetX();
920 Float_t y = seed->GetY();
921 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
91162307 922
c9427e08 923 if (y > ymax) {
924 seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
925 if (!seed->Rotate(fSectors->GetAlpha()))
926 return;
927 } else if (y <-ymax) {
928 seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
929 if (!seed->Rotate(-fSectors->GetAlpha()))
930 return;
931 }
1c53abe2 932
c9427e08 933}
1c53abe2 934
935
936
1c53abe2 937//_____________________________________________________________________________
b67e07dc 938Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1c53abe2 939 Double_t x2,Double_t y2,
940 Double_t x3,Double_t y3)
941{
942 //-----------------------------------------------------------------
943 // Initial approximation of the track curvature
944 //-----------------------------------------------------------------
945 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
946 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
947 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
948 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
949 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
950
951 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
91162307 952 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1c53abe2 953 return -xr*yr/sqrt(xr*xr+yr*yr);
954}
955
956
91162307 957
1c53abe2 958//_____________________________________________________________________________
b67e07dc 959Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
91162307 960 Double_t x2,Double_t y2,
961 Double_t x3,Double_t y3)
962{
963 //-----------------------------------------------------------------
964 // Initial approximation of the track curvature
965 //-----------------------------------------------------------------
966 x3 -=x1;
967 x2 -=x1;
968 y3 -=y1;
969 y2 -=y1;
970 //
971 Double_t det = x3*y2-x2*y3;
972 if (det==0) {
973 return 100;
974 }
975 //
976 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
977 Double_t x0 = x3*0.5-y3*u;
978 Double_t y0 = y3*0.5+x3*u;
979 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
980 if (det<0) c2*=-1;
981 return c2;
982}
983
984
b67e07dc 985Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1c53abe2 986 Double_t x2,Double_t y2,
987 Double_t x3,Double_t y3)
91162307 988{
989 //-----------------------------------------------------------------
990 // Initial approximation of the track curvature
991 //-----------------------------------------------------------------
992 x3 -=x1;
993 x2 -=x1;
994 y3 -=y1;
995 y2 -=y1;
996 //
997 Double_t det = x3*y2-x2*y3;
998 if (det==0) {
999 return 100;
1000 }
1001 //
1002 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1003 Double_t x0 = x3*0.5-y3*u;
1004 Double_t y0 = y3*0.5+x3*u;
1005 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1006 if (det<0) c2*=-1;
1007 x0+=x1;
1008 x0*=c2;
1009 return x0;
1010}
1011
1012
1013
1014//_____________________________________________________________________________
b67e07dc 1015Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
91162307 1016 Double_t x2,Double_t y2,
1017 Double_t x3,Double_t y3)
1c53abe2 1018{
1019 //-----------------------------------------------------------------
1020 // Initial approximation of the track curvature times center of curvature
1021 //-----------------------------------------------------------------
1022 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1023 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1024 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1025 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1026 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1027
1028 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1029
1030 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1031}
1032
1033//_____________________________________________________________________________
b67e07dc 1034Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1c53abe2 1035 Double_t x2,Double_t y2,
1036 Double_t z1,Double_t z2)
1037{
1038 //-----------------------------------------------------------------
1039 // Initial approximation of the tangent of the track dip angle
1040 //-----------------------------------------------------------------
1041 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1042}
1043
1044
b67e07dc 1045Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
91162307 1046 Double_t x2,Double_t y2,
1047 Double_t z1,Double_t z2, Double_t c)
1c53abe2 1048{
91162307 1049 //-----------------------------------------------------------------
1050 // Initial approximation of the tangent of the track dip angle
1051 //-----------------------------------------------------------------
1052
1053 // Double_t angle1;
1054
1055 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1c53abe2 1056 //
91162307 1057 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1058 if (TMath::Abs(d*c*0.5)>1) return 0;
1059 // Double_t angle2 = TMath::ASin(d*c*0.5);
b67e07dc 1060 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1061 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
91162307 1062
1063 angle2 = (z1-z2)*c/(angle2*2.);
1064 return angle2;
1065}
1066
1067Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1068{//-----------------------------------------------------------------
1069 // This function find proloncation of a track to a reference plane x=x2.
1070 //-----------------------------------------------------------------
1071
1072 Double_t dx=x2-x1;
1073
1074 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1075 return kFALSE;
1c53abe2 1076 }
f8aae377 1077
91162307 1078 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1079 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1080 y = x[0];
1081 z = x[1];
1082
1083 Double_t dy = dx*(c1+c2)/(r1+r2);
1084 Double_t dz = 0;
1085 //
1086 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1087
1088 if (TMath::Abs(delta)>0.01){
1089 dz = x[3]*TMath::ASin(delta)/x[4];
1090 }else{
b67e07dc 1091 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
91162307 1092 }
1093
b67e07dc 1094 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
f8aae377 1095
91162307 1096 y+=dy;
1097 z+=dz;
1098
1099 return kTRUE;
1c53abe2 1100}
1101
d26d9159 1102Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1103{
1104 //
1105 //
1106 fInput = tree;
1107 return LoadClusters();
1108}
91162307 1109
1110Int_t AliTPCtrackerMI::LoadClusters()
1c53abe2 1111{
1112 //
1113 // load clusters to the memory
91162307 1114 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1115 clrow->SetClass("AliTPCclusterMI");
1116 clrow->SetArray(0);
1117 clrow->GetArray()->ExpandCreateFast(10000);
1118 //
1119 // TTree * tree = fClustersArray.GetTree();
1120
1121 TTree * tree = fInput;
1122 TBranch * br = tree->GetBranch("Segment");
1123 br->SetAddress(&clrow);
1124 //
1125 Int_t j=Int_t(tree->GetEntries());
1c53abe2 1126 for (Int_t i=0; i<j; i++) {
91162307 1127 br->GetEntry(i);
1128 //
1129 Int_t sec,row;
1130 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
0201b65c 1131 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1132 Transform((AliCluster*)(clrow->GetArray()->At(icl)));
1133 }
91162307 1134 //
1135 AliTPCRow * tpcrow=0;
1136 Int_t left=0;
1137 if (sec<fkNIS*2){
1138 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1139 left = sec/fkNIS;
1140 }
1141 else{
1142 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1143 left = (sec-fkNIS*2)/fkNOS;
1144 }
1145 if (left ==0){
1146 tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1147 tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1148 for (Int_t i=0;i<tpcrow->fN1;i++)
1149 tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1150 }
1151 if (left ==1){
1152 tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1153 tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1154 for (Int_t i=0;i<tpcrow->fN2;i++)
1155 tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1156 }
1c53abe2 1157 }
91162307 1158 //
1159 delete clrow;
1160 LoadOuterSectors();
1161 LoadInnerSectors();
1162 return 0;
1c53abe2 1163}
1164
1165
91162307 1166void AliTPCtrackerMI::UnloadClusters()
1167{
1168 //
1169 // unload clusters from the memory
1170 //
1171 Int_t nrows = fOuterSec->GetNRows();
1172 for (Int_t sec = 0;sec<fkNOS;sec++)
1173 for (Int_t row = 0;row<nrows;row++){
1174 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
982aff31 1175 // if (tpcrow){
1176 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1177 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1178 //}
1179 tpcrow->ResetClusters();
1c53abe2 1180 }
91162307 1181 //
1182 nrows = fInnerSec->GetNRows();
1183 for (Int_t sec = 0;sec<fkNIS;sec++)
1184 for (Int_t row = 0;row<nrows;row++){
1185 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
982aff31 1186 //if (tpcrow){
1187 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1188 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1189 //}
1190 tpcrow->ResetClusters();
91162307 1191 }
1192
1193 return ;
1c53abe2 1194}
1195
0201b65c 1196void AliTPCtrackerMI::Transform(AliCluster * cluster){
1197 //
1198 //
1199 //
1200 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1201 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
e7eb17e4 1202 //TGeoHMatrix mat;
0201b65c 1203 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1204 Double_t posC[3];
e7eb17e4 1205 //mat.LocalToMaster(pos,posC);
0201b65c 1206 mat->LocalToMaster(pos,posC);
1207 cluster->SetX(posC[0]);
1208 cluster->SetY(posC[1]);
1209 cluster->SetZ(posC[2]);
1210}
1c53abe2 1211
1212//_____________________________________________________________________________
91162307 1213Int_t AliTPCtrackerMI::LoadOuterSectors() {
1c53abe2 1214 //-----------------------------------------------------------------
91162307 1215 // This function fills outer TPC sectors with clusters.
1c53abe2 1216 //-----------------------------------------------------------------
91162307 1217 Int_t nrows = fOuterSec->GetNRows();
1218 UInt_t index=0;
1219 for (Int_t sec = 0;sec<fkNOS;sec++)
1220 for (Int_t row = 0;row<nrows;row++){
1221 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1222 Int_t sec2 = sec+2*fkNIS;
1223 //left
1224 Int_t ncl = tpcrow->fN1;
1225 while (ncl--) {
1226 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1227 index=(((sec2<<8)+row)<<16)+ncl;
1228 tpcrow->InsertCluster(c,index);
1229 }
1230 //right
1231 ncl = tpcrow->fN2;
1232 while (ncl--) {
1233 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1234 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1235 tpcrow->InsertCluster(c,index);
1236 }
1237 //
1238 // write indexes for fast acces
1239 //
1240 for (Int_t i=0;i<510;i++)
1241 tpcrow->fFastCluster[i]=-1;
1242 for (Int_t i=0;i<tpcrow->GetN();i++){
1243 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1244 tpcrow->fFastCluster[zi]=i; // write index
1245 }
1246 Int_t last = 0;
1247 for (Int_t i=0;i<510;i++){
1248 if (tpcrow->fFastCluster[i]<0)
1249 tpcrow->fFastCluster[i] = last;
1250 else
1251 last = tpcrow->fFastCluster[i];
1252 }
1253 }
1254 fN=fkNOS;
1255 fSectors=fOuterSec;
1256 return 0;
1257}
1258
1259
1260//_____________________________________________________________________________
1261Int_t AliTPCtrackerMI::LoadInnerSectors() {
1262 //-----------------------------------------------------------------
1263 // This function fills inner TPC sectors with clusters.
1264 //-----------------------------------------------------------------
1265 Int_t nrows = fInnerSec->GetNRows();
1266 UInt_t index=0;
1267 for (Int_t sec = 0;sec<fkNIS;sec++)
1268 for (Int_t row = 0;row<nrows;row++){
1269 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1270 //
1271 //left
1272 Int_t ncl = tpcrow->fN1;
1273 while (ncl--) {
1274 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1275 index=(((sec<<8)+row)<<16)+ncl;
1276 tpcrow->InsertCluster(c,index);
1277 }
1278 //right
1279 ncl = tpcrow->fN2;
1280 while (ncl--) {
1281 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1282 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1283 tpcrow->InsertCluster(c,index);
1284 }
1285 //
1286 // write indexes for fast acces
1287 //
1288 for (Int_t i=0;i<510;i++)
1289 tpcrow->fFastCluster[i]=-1;
1290 for (Int_t i=0;i<tpcrow->GetN();i++){
1291 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1292 tpcrow->fFastCluster[zi]=i; // write index
1293 }
1294 Int_t last = 0;
1295 for (Int_t i=0;i<510;i++){
1296 if (tpcrow->fFastCluster[i]<0)
1297 tpcrow->fFastCluster[i] = last;
1298 else
1299 last = tpcrow->fFastCluster[i];
1300 }
1c53abe2 1301
91162307 1302 }
1303
1c53abe2 1304 fN=fkNIS;
1305 fSectors=fInnerSec;
91162307 1306 return 0;
1307}
1308
1309
1310
1311//_________________________________________________________________________
1312AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1313 //--------------------------------------------------------------------
1314 // Return pointer to a given cluster
1315 //--------------------------------------------------------------------
e7eb17e4 1316 if (index<0) return 0; // no cluster
91162307 1317 Int_t sec=(index&0xff000000)>>24;
1318 Int_t row=(index&0x00ff0000)>>16;
d26d9159 1319 Int_t ncl=(index&0x00007fff)>>00;
91162307 1320
1321 const AliTPCRow * tpcrow=0;
1322 AliTPCclusterMI * clrow =0;
3da363a1 1323
ad23441a 1324 if (sec<0 || sec>=fkNIS*4) {
1325 AliWarning(Form("Wrong sector %d",sec));
1326 return 0x0;
1327 }
1328
91162307 1329 if (sec<fkNIS*2){
1330 tpcrow = &(fInnerSec[sec%fkNIS][row]);
3da363a1 1331 if (tpcrow==0) return 0;
1332
1333 if (sec<fkNIS) {
1334 if (tpcrow->fN1<=ncl) return 0;
91162307 1335 clrow = tpcrow->fClusters1;
3da363a1 1336 }
1337 else {
1338 if (tpcrow->fN2<=ncl) return 0;
91162307 1339 clrow = tpcrow->fClusters2;
3da363a1 1340 }
91162307 1341 }
3da363a1 1342 else {
91162307 1343 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
3da363a1 1344 if (tpcrow==0) return 0;
1345
1346 if (sec-2*fkNIS<fkNOS) {
1347 if (tpcrow->fN1<=ncl) return 0;
91162307 1348 clrow = tpcrow->fClusters1;
3da363a1 1349 }
1350 else {
1351 if (tpcrow->fN2<=ncl) return 0;
91162307 1352 clrow = tpcrow->fClusters2;
3da363a1 1353 }
91162307 1354 }
3da363a1 1355
91162307 1356 return &(clrow[ncl]);
1357
1c53abe2 1358}
1359
91162307 1360
1361
1c53abe2 1362Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1363 //-----------------------------------------------------------------
1364 // This function tries to find a track prolongation to next pad row
1365 //-----------------------------------------------------------------
1c53abe2 1366 //
91162307 1367 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
4d158c36 1368 AliTPCclusterMI *cl=0;
1369 Int_t tpcindex= t.GetClusterIndex2(nr);
1370 //
1371 // update current shape info every 5 pad-row
1372 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1373 GetShape(&t,nr);
1374 //}
1375 //
1376 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1377 //
1378 if (tpcindex==-1) return 0; //track in dead zone
1379 if (tpcindex>0){ //
1380 cl = t.fClusterPointer[nr];
e0656451 1381 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
4d158c36 1382 t.fCurrentClusterIndex1 = tpcindex;
1383 }
1384 if (cl){
1385 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1386 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1387 //
1388 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1389 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1390
1391 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1392 Double_t rotation = angle-t.GetAlpha();
1393 t.fRelativeSector= relativesector;
3f82c4f2 1394 if (!t.Rotate(rotation)) return 0;
4d158c36 1395 }
3f82c4f2 1396 if (!t.PropagateTo(x)) return 0;
4d158c36 1397 //
1398 t.fCurrentCluster = cl;
1399 t.fRow = nr;
1400 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1401 if ((tpcindex&0x8000)==0) accept =0;
1402 if (accept<3) {
1403 //if founded cluster is acceptible
1404 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1405 t.fErrorY2 += 0.03;
1406 t.fErrorZ2 += 0.03;
1407 t.fErrorY2 *= 3;
1408 t.fErrorZ2 *= 3;
1409 }
1410 t.fNFoundable++;
1411 UpdateTrack(&t,accept);
1412 return 1;
1413 }
1414 }
1627d1c4 1415 }
3f82c4f2 1416 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1417 if (fIteration>1){
1418 // not look for new cluster during refitting
1419 t.fNFoundable++;
1420 return 0;
1421 }
91162307 1422 //
4d158c36 1423 UInt_t index=0;
ca142b1f 1424 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
4d158c36 1425 Double_t y=t.GetYat(x);
91162307 1426 if (TMath::Abs(y)>ymax){
1427 if (y > ymax) {
1428 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1429 if (!t.Rotate(fSectors->GetAlpha()))
1430 return 0;
1431 } else if (y <-ymax) {
1432 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1433 if (!t.Rotate(-fSectors->GetAlpha()))
1434 return 0;
1435 }
4d158c36 1436 //return 1;
91162307 1437 }
1438 //
4d158c36 1439 if (!t.PropagateTo(x)) {
1440 if (fIteration==0) t.fRemoval = 10;
1441 return 0;
91162307 1442 }
4d158c36 1443 y=t.GetY();
1444 Double_t z=t.GetZ();
1445 //
a3232aae 1446
1447 if (!IsActive(t.fRelativeSector,nr)) {
1448 t.fInDead = kTRUE;
1449 t.SetClusterIndex2(nr,-1);
1450 return 0;
1451 }
1452 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1453 Bool_t isActive = IsActive(t.fRelativeSector,nr);
1454 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.fRelativeSector][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.fRelativeSector][nr].GetN()>0;
1455
1456 if (!isActive || !isActive2) return 0;
1457
91162307 1458 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1459 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1460 Double_t roady =1.;
1461 Double_t roadz = 1.;
1462 //
1c53abe2 1463 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1464 t.fInDead = kTRUE;
91162307 1465 t.SetClusterIndex2(nr,-1);
1c53abe2 1466 return 0;
1467 }
1468 else
1469 {
3f82c4f2 1470 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength() && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1471 t.fNFoundable++;
1627d1c4 1472 else
1473 return 0;
1c53abe2 1474 }
1475 //calculate
91162307 1476 if (krow) {
1477 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1478 cl = krow.FindNearest2(y,z,roady,roadz,index);
1479 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1480 }
91162307 1481 if (cl) {
1482 t.fCurrentCluster = cl;
1483 t.fRow = nr;
4d158c36 1484 if (fIteration==2&&cl->IsUsed(10)) return 0;
91162307 1485 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
4d158c36 1486 if (fIteration==2&&cl->IsUsed(11)) {
1487 t.fErrorY2 += 0.03;
1488 t.fErrorZ2 += 0.03;
1489 t.fErrorY2 *= 3;
1490 t.fErrorZ2 *= 3;
1491 }
d26d9159 1492 /*
91162307 1493 if (t.fCurrentCluster->IsUsed(10)){
1494 //
1495 //
c9427e08 1496
91162307 1497 t.fNShared++;
1498 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1499 t.fRemoval =10;
1500 return 0;
1501 }
1502 }
d26d9159 1503 */
91162307 1504 if (accept<3) UpdateTrack(&t,accept);
c9427e08 1505
91162307 1506 } else {
982aff31 1507 if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
91162307 1508
1509 }
1510 return 1;
1511}
c9427e08 1512
91162307 1513Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1514 //-----------------------------------------------------------------
1515 // This function tries to find a track prolongation to next pad row
1516 //-----------------------------------------------------------------
1517 //
1518 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1519 Double_t y,z;
1520 if (!t.GetProlongation(x,y,z)) {
1521 t.fRemoval = 10;
1522 return 0;
1523 }
1524 //
1525 //
1526 if (TMath::Abs(y)>ymax){
91162307 1527
1c53abe2 1528 if (y > ymax) {
1529 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1530 if (!t.Rotate(fSectors->GetAlpha()))
1531 return 0;
1532 } else if (y <-ymax) {
1533 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1534 if (!t.Rotate(-fSectors->GetAlpha()))
1535 return 0;
91162307 1536 }
1537 if (!t.PropagateTo(x)) {
1538 return 0;
1539 }
1540 t.GetProlongation(x,y,z);
1541 }
1542 //
1543 // update current shape info every 3 pad-row
1544 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1545 // t.fCurrentSigmaY = GetSigmaY(&t);
1546 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1547 GetShape(&t,nr);
1548 }
1549 //
1550 AliTPCclusterMI *cl=0;
1551 UInt_t index=0;
1552
1553
1554 //Int_t nr2 = nr;
1555 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1556 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1557 Double_t roady =1.;
1558 Double_t roadz = 1.;
1559 //
1560 Int_t row = nr;
1561 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1562 t.fInDead = kTRUE;
1563 t.SetClusterIndex2(row,-1);
1564 return 0;
1565 }
1566 else
1567 {
d7a11555 1568 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1c53abe2 1569 }
91162307 1570 //calculate
1571
1572 if ((cl==0)&&(krow)) {
1573 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1574 cl = krow.FindNearest2(y,z,roady,roadz,index);
1575
1576 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1577 }
1578
1579 if (cl) {
1580 t.fCurrentCluster = cl;
1581 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1582 //if (accept<3){
1583 t.SetClusterIndex2(row,index);
1584 t.fClusterPointer[row] = cl;
1585 //}
1c53abe2 1586 }
1587 return 1;
1588}
1589
1590
91162307 1591
5d837844 1592//_________________________________________________________________________
1593Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1594{
1595 // Get track space point by index
1596 // return false in case the cluster doesn't exist
1597 AliTPCclusterMI *cl = GetClusterMI(index);
1598 if (!cl) return kFALSE;
1599 Int_t sector = (index&0xff000000)>>24;
0201b65c 1600 // Int_t row = (index&0x00ff0000)>>16;
5d837844 1601 Float_t xyz[3];
0201b65c 1602 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1603 xyz[0] = cl->GetX();
5d837844 1604 xyz[1] = cl->GetY();
1605 xyz[2] = cl->GetZ();
1606 Float_t sin,cos;
1607 fParam->AdjustCosSin(sector,cos,sin);
1608 Float_t x = cos*xyz[0]-sin*xyz[1];
1609 Float_t y = cos*xyz[1]+sin*xyz[0];
1610 Float_t cov[6];
1611 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1612 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1613 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1614 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1615 cov[0] = sin*sin*sigmaY2;
1616 cov[1] = -sin*cos*sigmaY2;
1617 cov[2] = 0.;
1618 cov[3] = cos*cos*sigmaY2;
1619 cov[4] = 0.;
1620 cov[5] = sigmaZ2;
1621 p.SetXYZ(x,y,xyz[2],cov);
1622 AliAlignObj::ELayerID iLayer;
1623 Int_t idet;
1624 if (sector < fParam->GetNInnerSector()) {
1625 iLayer = AliAlignObj::kTPC1;
1626 idet = sector;
1627 }
1628 else {
1629 iLayer = AliAlignObj::kTPC2;
1630 idet = sector - fParam->GetNInnerSector();
1631 }
1632 UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,idet);
1633 p.SetVolumeID(volid);
1634 return kTRUE;
1635}
1636
1637
1638
91162307 1639Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1c53abe2 1640 //-----------------------------------------------------------------
1641 // This function tries to find a track prolongation to next pad row
1642 //-----------------------------------------------------------------
1643 t.fCurrentCluster = 0;
1644 t.fCurrentClusterIndex1 = 0;
1c53abe2 1645
1646 Double_t xt=t.GetX();
91162307 1647 Int_t row = GetRowNumber(xt)-1;
1648 Double_t ymax= GetMaxY(nr);
1649
1c53abe2 1650 if (row < nr) return 1; // don't prolongate if not information until now -
ca142b1f 1651// if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1652// t.fRemoval =10;
1653// return 0; // not prolongate strongly inclined tracks
1654// }
1655// if (TMath::Abs(t.GetSnp())>0.95) {
1656// t.fRemoval =10;
1657// return 0; // not prolongate strongly inclined tracks
1658// }// patch 28 fev 06
91162307 1659
1660 Double_t x= GetXrow(nr);
1661 Double_t y,z;
1662 //t.PropagateTo(x+0.02);
1663 //t.PropagateTo(x+0.01);
1627d1c4 1664 if (!t.PropagateTo(x)){
1627d1c4 1665 return 0;
1666 }
1c53abe2 1667 //
91162307 1668 y=t.GetY();
1669 z=t.GetZ();
1c53abe2 1670
91162307 1671 if (TMath::Abs(y)>ymax){
1672 if (y > ymax) {
1673 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1674 if (!t.Rotate(fSectors->GetAlpha()))
1675 return 0;
1676 } else if (y <-ymax) {
1677 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1678 if (!t.Rotate(-fSectors->GetAlpha()))
1679 return 0;
1680 }
982aff31 1681 // if (!t.PropagateTo(x)){
1682 // return 0;
1683 //}
1684 return 1;
1685 //y = t.GetY();
1c53abe2 1686 }
91162307 1687 //
3f82c4f2 1688 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
a3232aae 1689
1690 if (!IsActive(t.fRelativeSector,nr)) {
1691 t.fInDead = kTRUE;
1692 t.SetClusterIndex2(nr,-1);
1693 return 0;
1694 }
1695 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1696
91162307 1697 AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1c53abe2 1698
1699 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1700 t.fInDead = kTRUE;
91162307 1701 t.SetClusterIndex2(nr,-1);
1c53abe2 1702 return 0;
1703 }
1704 else
1705 {
3f82c4f2 1706 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.fNFoundable++;
1627d1c4 1707 else
1708 return 0;
1c53abe2 1709 }
1c53abe2 1710
91162307 1711 // update current
1712 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1713 // t.fCurrentSigmaY = GetSigmaY(&t);
1714 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1715 GetShape(&t,nr);
1716 }
1c53abe2 1717
91162307 1718 AliTPCclusterMI *cl=0;
d184f295 1719 Int_t index=0;
91162307 1720 //
1721 Double_t roady = 1.;
1722 Double_t roadz = 1.;
1723 //
d26d9159 1724
1725 if (!cl){
1726 index = t.GetClusterIndex2(nr);
1727 if ( (index>0) && (index&0x8000)==0){
1728 cl = t.fClusterPointer[nr];
1729 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1730 t.fCurrentClusterIndex1 = index;
1731 if (cl) {
1732 t.fCurrentCluster = cl;
1733 return 1;
1734 }
1735 }
1736 }
1737
3d172d79 1738 // if (index<0) return 0;
1739 UInt_t uindex = TMath::Abs(index);
d184f295 1740
91162307 1741 if (krow) {
d184f295 1742 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1743 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
91162307 1744 }
d26d9159 1745
d184f295 1746 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(uindex);
d26d9159 1747 t.fCurrentCluster = cl;
1748
91162307 1749 return 1;
1750}
1c53abe2 1751
1c53abe2 1752
91162307 1753Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1754 //-----------------------------------------------------------------
1755 // This function tries to find a track prolongation to next pad row
1756 //-----------------------------------------------------------------
1c53abe2 1757
91162307 1758 //update error according neighborhoud
1c53abe2 1759
91162307 1760 if (t.fCurrentCluster) {
1761 t.fRow = nr;
8c51a605 1762 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
91162307 1763
1764 if (t.fCurrentCluster->IsUsed(10)){
1765 //
1766 //
1767 // t.fErrorZ2*=2;
1768 // t.fErrorY2*=2;
1769 t.fNShared++;
1770 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1771 t.fRemoval =10;
1772 return 0;
1773 }
b364ca79 1774 }
d26d9159 1775 if (fIteration>0) accept = 0;
b364ca79 1776 if (accept<3) UpdateTrack(&t,accept);
1777
1c53abe2 1778 } else {
91162307 1779 if (fIteration==0){
1780 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1781 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1782
1783 if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1c53abe2 1784 }
1785 }
1786 return 1;
1787}
1788
1789
1790
91162307 1791//_____________________________________________________________________________
1792Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1793 //-----------------------------------------------------------------
91162307 1794 // This function tries to find a track prolongation.
1c53abe2 1795 //-----------------------------------------------------------------
1796 Double_t xt=t.GetX();
1797 //
91162307 1798 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1c53abe2 1799 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1800 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1801 //
1802 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1803
91162307 1804 Int_t first = GetRowNumber(xt)-1;
51ad6848 1805 for (Int_t nr= first; nr>=rf; nr-=step) {
1806 // update kink info
1807 if (t.GetKinkIndexes()[0]>0){
1808 for (Int_t i=0;i<3;i++){
1809 Int_t index = t.GetKinkIndexes()[i];
1810 if (index==0) break;
1811 if (index<0) continue;
1812 //
6c94f330 1813 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 1814 if (!kink){
1815 printf("PROBLEM\n");
1816 }
1817 else{
eea478d3 1818 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 1819 if (kinkrow==nr){
1820 AliExternalTrackParam paramd(t);
1821 kink->SetDaughter(paramd);
eea478d3 1822 kink->SetStatus(2,5);
51ad6848 1823 kink->Update();
51ad6848 1824 }
1825 }
1826 }
1827 }
1828
1829 if (nr==80) t.UpdateReference();
982aff31 1830 if (nr<fInnerSec->GetNRows())
1831 fSectors = fInnerSec;
1832 else
1833 fSectors = fOuterSec;
91162307 1834 if (FollowToNext(t,nr)==0)
4d158c36 1835 if (!t.IsActive())
1836 return 0;
91162307 1837
1838 }
1839 return 1;
1840}
1841
1c53abe2 1842
1843//_____________________________________________________________________________
91162307 1844Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1845 //-----------------------------------------------------------------
1846 // This function tries to find a track prolongation.
1847 //-----------------------------------------------------------------
1848 Double_t xt=t.GetX();
1849 //
1850 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1851 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1852 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1853 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1854
1855 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1856
1857 if (FollowToNextFast(t,nr)==0)
1858 if (!t.IsActive()) return 0;
1c53abe2 1859
1c53abe2 1860 }
1861 return 1;
1862}
1863
1864
91162307 1865
1866
1867
1c53abe2 1868Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1869 //-----------------------------------------------------------------
1870 // This function tries to find a track prolongation.
1871 //-----------------------------------------------------------------
1c53abe2 1872 //
eea478d3 1873 Double_t xt=t.GetX();
1c53abe2 1874 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1875 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1876 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1877 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1878
4d158c36 1879 Int_t first = t.fFirstPoint;
eea478d3 1880 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
91162307 1881 //
1882 if (first<0) first=0;
4d158c36 1883 for (Int_t nr=first; nr<=rf; nr++) {
ca142b1f 1884 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
51ad6848 1885 if (t.GetKinkIndexes()[0]<0){
1886 for (Int_t i=0;i<3;i++){
1887 Int_t index = t.GetKinkIndexes()[i];
1888 if (index==0) break;
1889 if (index>0) continue;
1890 index = TMath::Abs(index);
6c94f330 1891 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 1892 if (!kink){
1893 printf("PROBLEM\n");
1894 }
1895 else{
eea478d3 1896 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 1897 if (kinkrow==nr){
1898 AliExternalTrackParam paramm(t);
1899 kink->SetMother(paramm);
eea478d3 1900 kink->SetStatus(2,1);
51ad6848 1901 kink->Update();
51ad6848 1902 }
1903 }
eea478d3 1904 }
51ad6848 1905 }
eea478d3 1906 //
d26d9159 1907 if (nr<fInnerSec->GetNRows())
1908 fSectors = fInnerSec;
1909 else
1910 fSectors = fOuterSec;
c274e255 1911
d26d9159 1912 FollowToNext(t,nr);
1c53abe2 1913 }
1914 return 1;
1915}
1916
1917
1918
1919
1920
1921Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1922{
1923 //
1924 //
1925 sum1=0;
1926 sum2=0;
1927 Int_t sum=0;
1c53abe2 1928 //
1929 Float_t dz2 =(s1->GetZ() - s2->GetZ());
c9427e08 1930 dz2*=dz2;
91162307 1931
c9427e08 1932 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1c53abe2 1933 dy2*=dy2;
1934 Float_t distance = TMath::Sqrt(dz2+dy2);
c9427e08 1935 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1c53abe2 1936
91162307 1937 // Int_t offset =0;
c9427e08 1938 Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1939 Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
c9427e08 1940 if (lastpoint>160)
1941 lastpoint =160;
1942 if (firstpoint<0)
1943 firstpoint = 0;
91162307 1944 if (firstpoint>lastpoint) {
1945 firstpoint =lastpoint;
1946 // lastpoint =160;
c9427e08 1947 }
1948
1949
91162307 1950 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1951 if (s1->GetClusterIndex2(i)>0) sum1++;
1952 if (s2->GetClusterIndex2(i)>0) sum2++;
1953 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1c53abe2 1954 sum++;
1955 }
1956 }
91162307 1957 if (sum<5) return 0;
1958
1627d1c4 1959 Float_t summin = TMath::Min(sum1+1,sum2+1);
1960 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 1961 return ratio;
1962}
1963
1964void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1965{
1966 //
1967 //
91162307 1968 if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
1969 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
1970
1c53abe2 1971 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1972 dz2*=dz2;
1973 Float_t dy2 =(s1->GetY() - s2->GetY());
1974 dy2*=dy2;
91162307 1975 Float_t distance = dz2+dy2;
1976 if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics
1977
1c53abe2 1978 //
91162307 1979 Int_t sumshared=0;
1980 //
1981 Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
1982 Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
1983 //
1984 if (firstpoint>=lastpoint-5) return;;
1c53abe2 1985
91162307 1986 for (Int_t i=firstpoint;i<lastpoint;i++){
1987 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1988 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1989 sumshared++;
1990 }
1991 }
1992 if (sumshared>4){
1993 // sign clusters
1994 //
1995 for (Int_t i=firstpoint;i<lastpoint;i++){
1996 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1997 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1998 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
1999 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2000 if (s1->IsActive()&&s2->IsActive()){
2001 p1->fIsShared = kTRUE;
2002 p2->fIsShared = kTRUE;
2003 }
2004 }
2005 }
2006 }
2007 //
2008 if (sumshared>10){
2009 for (Int_t i=0;i<4;i++){
2010 if (s1->fOverlapLabels[3*i]==0){
2011 s1->fOverlapLabels[3*i] = s2->GetLabel();
2012 s1->fOverlapLabels[3*i+1] = sumshared;
2013 s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
2014 break;
2015 }
2016 }
2017 for (Int_t i=0;i<4;i++){
2018 if (s2->fOverlapLabels[3*i]==0){
2019 s2->fOverlapLabels[3*i] = s1->GetLabel();
2020 s2->fOverlapLabels[3*i+1] = sumshared;
2021 s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
2022 break;
2023 }
2024 }
2025 }
1c53abe2 2026
91162307 2027}
1c53abe2 2028
91162307 2029void AliTPCtrackerMI::SignShared(TObjArray * arr)
2030{
1c53abe2 2031 //
91162307 2032 //sort trackss according sectors
2033 //
c9427e08 2034 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2035 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2036 if (!pt) continue;
2037 //if (pt) RotateToLocal(pt);
2038 pt->fSort = 0;
c9427e08 2039 }
91162307 2040 arr->UnSort();
1c53abe2 2041 arr->Sort(); // sorting according z
2042 arr->Expand(arr->GetEntries());
91162307 2043 //
2044 //
1c53abe2 2045 Int_t nseed=arr->GetEntriesFast();
1c53abe2 2046 for (Int_t i=0; i<nseed; i++) {
2047 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2048 if (!pt) continue;
2049 for (Int_t j=0;j<=12;j++){
2050 pt->fOverlapLabels[j] =0;
1c53abe2 2051 }
91162307 2052 }
2053 for (Int_t i=0; i<nseed; i++) {
2054 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2055 if (!pt) continue;
2056 if (pt->fRemoval>10) continue;
1c53abe2 2057 for (Int_t j=i+1; j<nseed; j++){
2058 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
91162307 2059 // if (pt2){
2060 if (pt2->fRemoval<=10) {
2061 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2062 SignShared(pt,pt2);
c9427e08 2063 }
91162307 2064 }
2065 }
2066}
2067
2068void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
2069{
2070 //
2071 //sort trackss according sectors
2072 //
2073 if (fDebug&1) {
6bdc18d6 2074 Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
91162307 2075 }
2076 //
2077 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2078 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2079 if (!pt) continue;
2080 pt->fSort = 0;
2081 }
2082 arr->UnSort();
2083 arr->Sort(); // sorting according z
2084 arr->Expand(arr->GetEntries());
2085 //
2086 //reset overlap labels
2087 //
2088 Int_t nseed=arr->GetEntriesFast();
2089 for (Int_t i=0; i<nseed; i++) {
2090 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2091 if (!pt) continue;
2092 pt->SetUniqueID(i);
2093 for (Int_t j=0;j<=12;j++){
2094 pt->fOverlapLabels[j] =0;
1c53abe2 2095 }
2096 }
91162307 2097 //
2098 //sign shared tracks
1c53abe2 2099 for (Int_t i=0; i<nseed; i++) {
2100 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2101 if (!pt) continue;
2102 if (pt->fRemoval>10) continue;
2103 Float_t deltac = pt->GetC()*0.1;
2104 for (Int_t j=i+1; j<nseed; j++){
2105 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2106 // if (pt2){
2107 if (pt2->fRemoval<=10) {
2108 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2109 if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue;
2110 if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2111 //
2112 SignShared(pt,pt2);
2113 }
1c53abe2 2114 }
1c53abe2 2115 }
91162307 2116 //
2117 // remove highly shared tracks
2118 for (Int_t i=0; i<nseed; i++) {
2119 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2120 if (!pt) continue;
2121 if (pt->fRemoval>10) continue;
2122 //
2123 Int_t sumshared =0;
2124 for (Int_t j=0;j<4;j++){
2125 sumshared = pt->fOverlapLabels[j*3+1];
2126 }
2127 Float_t factor = factor1;
2128 if (pt->fRemoval>0) factor = factor2;
2129 if (sumshared/pt->GetNumberOfClusters()>factor){
2130 for (Int_t j=0;j<4;j++){
2131 if (pt->fOverlapLabels[3*j]==0) continue;
2132 if (pt->fOverlapLabels[3*j+1]<5) continue;
2133 if (pt->fRemoval==removalindex) continue;
2134 AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2135 if (!pt2) continue;
2136 if (pt2->GetSigma2C()<pt->GetSigma2C()){
2137 // pt->fRemoval = removalindex;
2138 delete arr->RemoveAt(i);
2139 break;
1c53abe2 2140 }
91162307 2141 }
1c53abe2 2142 }
91162307 2143 }
2144 arr->Compress();
2145 if (fDebug&1) {
6bdc18d6 2146 Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
91162307 2147 }
2148}
2149
2150
2151
1c53abe2 2152
2153
91162307 2154
47966a6d 2155void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
91162307 2156{
2157 //
2158 //sort tracks in array according mode criteria
2159 Int_t nseed = arr->GetEntriesFast();
2160 for (Int_t i=0; i<nseed; i++) {
2161 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2162 if (!pt) {
2163 continue;
2164 }
2165 pt->fSort = mode;
2166 }
2167 arr->UnSort();
2168 arr->Sort();
1c53abe2 2169}
c9427e08 2170
91162307 2171void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
c9427e08 2172{
2173
2174 //Loop over all tracks and remove "overlaps"
2175 //
2176 //
2177 Int_t nseed = arr->GetEntriesFast();
2178 Int_t good =0;
91162307 2179
c9427e08 2180 for (Int_t i=0; i<nseed; i++) {
2181 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2182 if (!pt) {
91162307 2183 delete arr->RemoveAt(i);
c9427e08 2184 }
91162307 2185 else{
2186 pt->fSort =1;
2187 pt->fBSigned = kFALSE;
c9427e08 2188 }
91162307 2189 }
2190 arr->Compress();
2191 nseed = arr->GetEntriesFast();
2192 arr->UnSort();
2193 arr->Sort();
2194 //
2195 //unsign used
2196 UnsignClusters();
2197 //
2198 for (Int_t i=0; i<nseed; i++) {
2199 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2200 if (!pt) {
2201 continue;
2202 }
2203 Int_t found,foundable,shared;
2204 if (pt->IsActive())
2205 pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2206 else
2207 pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2208 //
2209 Double_t factor = factor2;
2210 if (pt->fBConstrain) factor = factor1;
2211
2212 if ((Float_t(shared)/Float_t(found))>factor){
c9427e08 2213 pt->Desactivate(removalindex);
91162307 2214 continue;
2215 }
2216
2217 good++;
2218 for (Int_t i=0; i<160; i++) {
2219 Int_t index=pt->GetClusterIndex2(i);
2220 if (index<0 || index&0x8000 ) continue;
2221 AliTPCclusterMI *c= pt->fClusterPointer[i];
2222 if (!c) continue;
2223 // if (!c->IsUsed(10)) c->Use(10);
2224 //if (pt->IsActive())
2225 c->Use(10);
2226 //else
2227 // c->Use(5);
c9427e08 2228 }
91162307 2229
c9427e08 2230 }
2231 fNtracks = good;
6bdc18d6 2232 if (fDebug>0){
2233 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2234 }
c9427e08 2235}
2236
51ad6848 2237
2238void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2239{
2240
2241 //Loop over all tracks and remove "overlaps"
2242 //
2243 //
2244 UnsignClusters();
2245 //
2246 Int_t nseed = arr->GetEntriesFast();
2247 Float_t * quality = new Float_t[nseed];
2248 Int_t * indexes = new Int_t[nseed];
2249 Int_t good =0;
2250 //
2251 //
2252 for (Int_t i=0; i<nseed; i++) {
2253 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2254 if (!pt){
2255 quality[i]=-1;
2256 continue;
2257 }
2258 pt->UpdatePoints(); //select first last max dens points
2259 Float_t * points = pt->GetPoints();
2260 if (points[3]<0.8) quality[i] =-1;
2261 //
2262 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2263 }
2264 TMath::Sort(nseed,quality,indexes);
2265 //
2266 //
2267 for (Int_t itrack=0; itrack<nseed; itrack++) {
2268 Int_t trackindex = indexes[itrack];
2269 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2270 if (quality[trackindex]<0){
2271 if (pt) {
2272 delete arr->RemoveAt(trackindex);
2273 }
2274 else{
2275 arr->RemoveAt(trackindex);
2276 }
2277 continue;
2278 }
2279 //
2280 Int_t first = Int_t(pt->GetPoints()[0]);
2281 Int_t last = Int_t(pt->GetPoints()[2]);
2282 Double_t factor = (pt->fBConstrain) ? factor1: factor2;
2283 //
2284 Int_t found,foundable,shared;
2285 pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
eea478d3 2286 Float_t sharedfactor = Float_t(shared+1)/Float_t(found+1);
b406fdb0 2287 Bool_t itsgold =kFALSE;
2288 if (pt->fEsd){
ef7253ac 2289 Int_t dummy[12];
b406fdb0 2290 if (pt->fEsd->GetITSclusters(dummy)>4) itsgold= kTRUE;
51ad6848 2291 }
b406fdb0 2292 if (!itsgold){
2293 //
2294 if (Float_t(shared+1)/Float_t(found+1)>factor){
2295 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2296 delete arr->RemoveAt(trackindex);
2297 continue;
2298 }
2299 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2300 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2301 delete arr->RemoveAt(trackindex);
2302 continue;
2303 }
51ad6848 2304 }
2305
2306 good++;
2307 if (sharedfactor>0.4) continue;
b406fdb0 2308 if (pt->GetKinkIndexes()[0]>0) continue;
51ad6848 2309 for (Int_t i=first; i<last; i++) {
2310 Int_t index=pt->GetClusterIndex2(i);
2311 // if (index<0 || index&0x8000 ) continue;
2312 if (index<0 || index&0x8000 ) continue;
2313 AliTPCclusterMI *c= pt->fClusterPointer[i];
2314 if (!c) continue;
2315 c->Use(10);
2316 }
2317 }
2318 fNtracks = good;
2319 if (fDebug>0){
2320 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2321 }
2322 delete []quality;
2323 delete []indexes;
2324}
2325
91162307 2326void AliTPCtrackerMI::UnsignClusters()
1c53abe2 2327{
91162307 2328 //
2329 // loop over all clusters and unsign them
2330 //
2331
2332 for (Int_t sec=0;sec<fkNIS;sec++){
2333 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2334 AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2335 for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2336 // if (cl[icl].IsUsed(10))
2337 cl[icl].Use(-1);
2338 cl = fInnerSec[sec][row].fClusters2;
2339 for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2340 //if (cl[icl].IsUsed(10))
2341 cl[icl].Use(-1);
2342 }
2343 }
2344
2345 for (Int_t sec=0;sec<fkNOS;sec++){
2346 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2347 AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2348 for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2349 //if (cl[icl].IsUsed(10))
2350 cl[icl].Use(-1);
2351 cl = fOuterSec[sec][row].fClusters2;
2352 for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2353 //if (cl[icl].IsUsed(10))
2354 cl[icl].Use(-1);
2355 }
2356 }
2357
1c53abe2 2358}
2359
91162307 2360
2361
732b9e78 2362void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
1c53abe2 2363{
2364 //
91162307 2365 //sign clusters to be "used"
2366 //
2367 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2368 // loop over "primaries"
2369
2370 Float_t sumdens=0;
2371 Float_t sumdens2=0;
2372 Float_t sumn =0;
2373 Float_t sumn2 =0;
2374 Float_t sumchi =0;
2375 Float_t sumchi2 =0;
2376
2377 Float_t sum =0;
2378
1c53abe2 2379 TStopwatch timer;
91162307 2380 timer.Start();
1c53abe2 2381
91162307 2382 Int_t nseed = arr->GetEntriesFast();
2383 for (Int_t i=0; i<nseed; i++) {
2384 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2385 if (!pt) {
2386 continue;
2387 }
2388 if (!(pt->IsActive())) continue;
2389 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2390 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2391 sumdens += dens;
2392 sumdens2+= dens*dens;
2393 sumn += pt->GetNumberOfClusters();
2394 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2395 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2396 if (chi2>5) chi2=5;
2397 sumchi +=chi2;
2398 sumchi2 +=chi2*chi2;
2399 sum++;
2400 }
1627d1c4 2401 }
91162307 2402
2403 Float_t mdensity = 0.9;
2404 Float_t meann = 130;
2405 Float_t meanchi = 1;
2406 Float_t sdensity = 0.1;
2407 Float_t smeann = 10;
2408 Float_t smeanchi =0.4;
1627d1c4 2409
91162307 2410
2411 if (sum>20){
2412 mdensity = sumdens/sum;
2413 meann = sumn/sum;
2414 meanchi = sumchi/sum;
2415 //
2416 sdensity = sumdens2/sum-mdensity*mdensity;
2417 sdensity = TMath::Sqrt(sdensity);
2418 //
2419 smeann = sumn2/sum-meann*meann;
2420 smeann = TMath::Sqrt(smeann);
2421 //
2422 smeanchi = sumchi2/sum - meanchi*meanchi;
2423 smeanchi = TMath::Sqrt(smeanchi);
2424 }
2425
2426
2427 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2428 //
2429 for (Int_t i=0; i<nseed; i++) {
2430 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2431 if (!pt) {
2432 continue;
1c53abe2 2433 }
91162307 2434 if (pt->fBSigned) continue;
2435 if (pt->fBConstrain) continue;
2436 //if (!(pt->IsActive())) continue;
2437 /*
2438 Int_t found,foundable,shared;
2439 pt->GetClusterStatistic(0,160,found, foundable,shared);
2440 if (shared/float(found)>0.3) {
2441 if (shared/float(found)>0.9 ){
2442 //delete arr->RemoveAt(i);
2443 }
2444 continue;
c9427e08 2445 }
91162307 2446 */
2447 Bool_t isok =kFALSE;
2448 if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2449 isok = kTRUE;
2450 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2451 isok =kTRUE;
2452 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2453 isok =kTRUE;
2454 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2455 isok =kTRUE;
2456
2457 if (isok)
2458 for (Int_t i=0; i<160; i++) {
2459 Int_t index=pt->GetClusterIndex2(i);
2460 if (index<0) continue;
2461 AliTPCclusterMI *c= pt->fClusterPointer[i];
2462 if (!c) continue;
2463 //if (!(c->IsUsed(10))) c->Use();
2464 c->Use(10);
2465 }
2466 }
2467
c9427e08 2468
1c53abe2 2469 //
91162307 2470 Double_t maxchi = meanchi+2.*smeanchi;
2471
2472 for (Int_t i=0; i<nseed; i++) {
2473 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2474 if (!pt) {
1c53abe2 2475 continue;
91162307 2476 }
2477 //if (!(pt->IsActive())) continue;
2478 if (pt->fBSigned) continue;
2479 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2480 if (chi>maxchi) continue;
2481
2482 Float_t bfactor=1;
2483 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2484
2485 //sign only tracks with enoug big density at the beginning
2486
2487 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2488
2489
3e5d0aa2 2490 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
91162307 2491 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2492
2493 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2494 if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2495 minn=0;
2496
2497 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2498 //Int_t noc=pt->GetNumberOfClusters();
2499 pt->fBSigned = kTRUE;
2500 for (Int_t i=0; i<160; i++) {
2501
2502 Int_t index=pt->GetClusterIndex2(i);
2503 if (index<0) continue;
2504 AliTPCclusterMI *c= pt->fClusterPointer[i];
2505 if (!c) continue;
2506 // if (!(c->IsUsed(10))) c->Use();
2507 c->Use(10);
2508 }
1c53abe2 2509 }
91162307 2510 }
2511 // gLastCheck = nseed;
2512 // arr->Compress();
2513 if (fDebug>0){
2514 timer.Print();
2515 }
1c53abe2 2516}
2517
2518
47966a6d 2519void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
91162307 2520{
2521 // stop not active tracks
2522 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2523 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2524 Int_t nseed = arr->GetEntriesFast();
2525 //
2526 for (Int_t i=0; i<nseed; i++) {
2527 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2528 if (!pt) {
2529 continue;
2530 }
2531 if (!(pt->IsActive())) continue;
2532 StopNotActive(pt,row0,th0, th1,th2);
2533 }
2534}
1c53abe2 2535
1c53abe2 2536
1c53abe2 2537
91162307 2538void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
47966a6d 2539 Float_t th2) const
91162307 2540{
2541 // stop not active tracks
2542 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2543 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2544 Int_t sumgood1 = 0;
2545 Int_t sumgood2 = 0;
2546 Int_t foundable = 0;
2547 Int_t maxindex = seed->fLastPoint; //last foundable row
2548 if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2549 seed->Desactivate(10) ;
2550 return;
2551 }
1c53abe2 2552
3e5d0aa2 2553 for (Int_t i=row0; i<maxindex; i++){
91162307 2554 Int_t index = seed->GetClusterIndex2(i);
2555 if (index!=-1) foundable++;
2556 //if (!c) continue;
2557 if (foundable<=30) sumgood1++;
2558 if (foundable<=50) {
2559 sumgood2++;
2560 }
2561 else{
2562 break;
2563 }
2564 }
2565 if (foundable>=30.){
2566 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2567 }
2568 if (foundable>=50)
2569 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2570}
1c53abe2 2571
1c53abe2 2572
d26d9159 2573Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2574{
2575 //
2576 // back propagation of ESD tracks
2577 //
2578 //return 0;
2579 fEvent = event;
2580 ReadSeeds(event,2);
2581 fIteration=2;
982aff31 2582 //PrepareForProlongation(fSeeds,1);
2583 PropagateForward2(fSeeds);
af32720d 2584
4d158c36 2585 Int_t ntracks=0;
d26d9159 2586 Int_t nseed = fSeeds->GetEntriesFast();
2587 for (Int_t i=0;i<nseed;i++){
2588 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
4d158c36 2589 if (!seed) continue;
eea478d3 2590 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2591
d26d9159 2592 seed->PropagateTo(fParam->GetInnerRadiusLow());
51ad6848 2593 seed->UpdatePoints();
d26d9159 2594 AliESDtrack *esd=event->GetTrack(i);
2595 seed->CookdEdx(0.02,0.6);
2596 CookLabel(seed,0.1); //For comparison only
af32720d 2597 //
34acb742 2598 if (AliTPCReconstructor::StreamLevel()>0 && seed!=0&&esd!=0) {
af32720d 2599 TTreeSRedirector &cstream = *fDebugStreamer;
2600 cstream<<"Crefit"<<
2601 "Esd.="<<esd<<
2602 "Track.="<<seed<<
2603 "\n";
2604 }
c274e255 2605
51ad6848 2606 if (seed->GetNumberOfClusters()>15){
4d158c36 2607 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
51ad6848 2608 esd->SetTPCPoints(seed->GetPoints());
b406fdb0 2609 esd->SetTPCPointsF(seed->fNFoundable);
2610 Int_t ndedx = seed->fNCDEDX[0]+seed->fNCDEDX[1]+seed->fNCDEDX[2]+seed->fNCDEDX[3];
2611 Float_t sdedx = (seed->fSDEDX[0]+seed->fSDEDX[1]+seed->fSDEDX[2]+seed->fSDEDX[3])*0.25;
2612 Float_t dedx = seed->GetdEdx();
2613 esd->SetTPCsignal(dedx, sdedx, ndedx);
fdedfdec 2614 //
2615 // add seed to the esd track in Calib level
2616 //
2617 if (AliTPCReconstructor::StreamLevel()>0){
2618 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2619 esd->AddCalibObject(seedCopy);
2620 }
4d158c36 2621 ntracks++;
d26d9159 2622 }
2623 else{
2624 //printf("problem\n");
2625 }
2626 }
51ad6848 2627 //FindKinks(fSeeds,event);
4d158c36 2628 Info("RefitInward","Number of refitted tracks %d",ntracks);
d26d9159 2629 fEvent =0;
2630 //WriteTracks();
2631 return 0;
2632}
2633
1c53abe2 2634
91162307 2635Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2636{
2637 //
2638 // back propagation of ESD tracks
2639 //
1c53abe2 2640
91162307 2641 fEvent = event;
d26d9159 2642 fIteration = 1;
5d837844 2643 ReadSeeds(event,1);
b406fdb0 2644 PropagateBack(fSeeds);
2645 RemoveUsed2(fSeeds,0.4,0.4,20);
2646 //
91162307 2647 Int_t nseed = fSeeds->GetEntriesFast();
4d158c36 2648 Int_t ntracks=0;
91162307 2649 for (Int_t i=0;i<nseed;i++){
2650 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
d9e62e7e 2651 if (!seed) continue;
51ad6848 2652 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2653 seed->UpdatePoints();
91162307 2654 AliESDtrack *esd=event->GetTrack(i);
d26d9159 2655 seed->CookdEdx(0.02,0.6);
91162307 2656 CookLabel(seed,0.1); //For comparison only
51ad6848 2657 if (seed->GetNumberOfClusters()>15){
4d158c36 2658 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
51ad6848 2659 esd->SetTPCPoints(seed->GetPoints());
167c41ab 2660 esd->SetTPCPointsF(seed->fNFoundable);
2661 Int_t ndedx = seed->fNCDEDX[0]+seed->fNCDEDX[1]+seed->fNCDEDX[2]+seed->fNCDEDX[3];
2662 Float_t sdedx = (seed->fSDEDX[0]+seed->fSDEDX[1]+seed->fSDEDX[2]+seed->fSDEDX[3])*0.25;
2663 Float_t dedx = seed->GetdEdx();
2664 esd->SetTPCsignal(dedx, sdedx, ndedx);
4d158c36 2665 ntracks++;
ca142b1f 2666 Int_t eventnumber = event->GetEventNumber();// patch 28 fev 06
2667 (*fDebugStreamer)<<"Cback"<<
2668 "Tr0.="<<seed<<
2669 "EventNr="<<eventnumber<<
2670 "\n"; // patch 28 fev 06
4d158c36 2671 }
91162307 2672 }
51ad6848 2673 //FindKinks(fSeeds,event);
4d158c36 2674 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
91162307 2675 fEvent =0;
d26d9159 2676 //WriteTracks();
6c94f330 2677
91162307 2678 return 0;
2679}
2680
2681
2682void AliTPCtrackerMI::DeleteSeeds()
2683{
b67e07dc 2684 //
2685 //delete Seeds
6c94f330 2686
91162307 2687 Int_t nseed = fSeeds->GetEntriesFast();
2688 for (Int_t i=0;i<nseed;i++){
2689 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2690 if (seed) delete fSeeds->RemoveAt(i);
2691 }
2692 delete fSeeds;
6c94f330 2693
91162307 2694 fSeeds =0;
2695}
2696
d26d9159 2697void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
91162307 2698{
2699 //
2700 //read seeds from the event
2701
2702 Int_t nentr=event->GetNumberOfTracks();
6bdc18d6 2703 if (fDebug>0){
2704 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2705 }
91162307 2706 if (fSeeds)
2707 DeleteSeeds();
2708 if (!fSeeds){
4d158c36 2709 fSeeds = new TObjArray(nentr);
91162307 2710 }
4d158c36 2711 UnsignClusters();
2712 // Int_t ntrk=0;
91162307 2713 for (Int_t i=0; i<nentr; i++) {
2714 AliESDtrack *esd=event->GetTrack(i);
51ad6848 2715 ULong_t status=esd->GetStatus();
2716 if (!(status&AliESDtrack::kTPCin)) continue;
b67e07dc 2717 AliTPCtrack t(*esd);
5d837844 2718 t.SetNumberOfClusters(0);
eea478d3 2719 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2720 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2721 for (Int_t ikink=0;ikink<3;ikink++) {
2722 Int_t index = esd->GetKinkIndex(ikink);
2723 seed->GetKinkIndexes()[ikink] = index;
2724 if (index==0) continue;
2725 index = TMath::Abs(index);
2726 AliESDkink * kink = fEvent->GetKink(index-1);
2727 if (kink&&esd->GetKinkIndex(ikink)<0){
2728 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2729 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2730 }
2731 if (kink&&esd->GetKinkIndex(ikink)>0){
2732 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2733 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2734 }
2735
2736 }
6c94f330 2737 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2738 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
4d158c36 2739 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2740 fSeeds->AddAt(0,i);
2741 delete seed;
2742 continue;
2743 }
2744 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
c0b978f0 2745 Double_t par0[5],par1[5],alpha,x;
2746 esd->GetInnerExternalParameters(alpha,x,par0);
4d158c36 2747 esd->GetExternalParameters(x,par1);
2748 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2749 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
51ad6848 2750 Double_t trdchi2=0;
2751 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
4d158c36 2752 //reset covariance if suspicious
51ad6848 2753 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
6c94f330 2754 seed->ResetCovariance(10.);
4d158c36 2755 }
982aff31 2756
91162307 2757 //
2758 //
2759 // rotate to the local coordinate system
eea478d3 2760 //
2761 fSectors=fInnerSec; fN=fkNIS;
91162307 2762 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2763 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2764 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2765 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2766 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
4d158c36 2767 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2768 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
91162307 2769 alpha-=seed->GetAlpha();
d9b8978b 2770 if (!seed->Rotate(alpha)) {
2771 delete seed;
2772 continue;
2773 }
d26d9159 2774 seed->fEsd = esd;
4d158c36 2775 // sign clusters
b406fdb0 2776 if (esd->GetKinkIndex(0)<=0){
2777 for (Int_t irow=0;irow<160;irow++){
2778 Int_t index = seed->GetClusterIndex2(irow);
2779 if (index>0){
2780 //
2781 AliTPCclusterMI * cl = GetClusterMI(index);
2782 seed->fClusterPointer[irow] = cl;
2783 if (cl){
2784 if ((index & 0x8000)==0){
2785 cl->Use(10); // accepted cluster
2786 }else{
2787 cl->Use(6); // close cluster not accepted
2788 }
4d158c36 2789 }else{
b406fdb0 2790 Info("ReadSeeds","Not found cluster");
2791 }
4d158c36 2792 }
2793 }
2794 }
2795 fSeeds->AddAt(seed,i);
91162307 2796 }
2797}
2798
2799
2800
2801//_____________________________________________________________________________
2802void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2803 Float_t deltay, Int_t ddsec) {
2804 //-----------------------------------------------------------------
2805 // This function creates track seeds.
2806 // SEEDING WITH VERTEX CONSTRAIN
2807 //-----------------------------------------------------------------
2808 // cuts[0] - fP4 cut
2809 // cuts[1] - tan(phi) cut
2810 // cuts[2] - zvertex cut
2811 // cuts[3] - fP3 cut
2812 Int_t nin0 = 0;
2813 Int_t nin1 = 0;
2814 Int_t nin2 = 0;
2815 Int_t nin = 0;
2816 Int_t nout1 = 0;
2817 Int_t nout2 = 0;
2818
2819 Double_t x[5], c[15];
2820 // Int_t di = i1-i2;
2821 //
2822 AliTPCseed * seed = new AliTPCseed;
2823 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2824 Double_t cs=cos(alpha), sn=sin(alpha);
2825 //
2826 // Double_t x1 =fOuterSec->GetX(i1);
2827 //Double_t xx2=fOuterSec->GetX(i2);
2828
2829 Double_t x1 =GetXrow(i1);
2830 Double_t xx2=GetXrow(i2);
2831
2832 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2833
2834 Int_t imiddle = (i2+i1)/2; //middle pad row index
2835 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2836 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2837 //
2838 Int_t ns =sec;
2839
2840 const AliTPCRow& kr1=GetRow(ns,i1);
2841 Double_t ymax = GetMaxY(i1)-kr1.fDeadZone-1.5;
2842 Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2843
2844 //
2845 // change cut on curvature if it can't reach this layer
2846 // maximal curvature set to reach it
2847 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2848 if (dvertexmax*0.5*cuts[0]>0.85){
2849 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2850 }
2851 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2852
2853 // Int_t ddsec = 1;
2854 if (deltay>0) ddsec = 0;
2855 // loop over clusters
2856 for (Int_t is=0; is < kr1; is++) {
2857 //
2858 if (kr1[is]->IsUsed(10)) continue;
2859 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2860 //if (TMath::Abs(y1)>ymax) continue;
2861
2862 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2863
2864 // find possible directions
2865 Float_t anglez = (z1-z3)/(x1-x3);
2866 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2867 //
2868 //
2869 //find rotation angles relative to line given by vertex and point 1
2870 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2871 Double_t dvertex = TMath::Sqrt(dvertex2);
2872 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2873 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2874
2875 //
2876 // loop over 2 sectors
2877 Int_t dsec1=-ddsec;
2878 Int_t dsec2= ddsec;
2879 if (y1<0) dsec2= 0;
2880 if (y1>0) dsec1= 0;
2881
2882 Double_t dddz1=0; // direction of delta inclination in z axis
2883 Double_t dddz2=0;
2884 if ( (z1-z3)>0)
2885 dddz1 =1;
2886 else
2887 dddz2 =1;
2888 //
2889 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2890 Int_t sec2 = sec + dsec;
2891 //
2892 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2893 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2894 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2895 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2896 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2897 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2898
2899 // rotation angles to p1-p3
2900 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2901 Double_t x2, y2, z2;
2902 //
2903 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2904
2905 //
2906 Double_t dxx0 = (xx2-x3)*cs13r;
2907 Double_t dyy0 = (xx2-x3)*sn13r;
2908 for (Int_t js=index1; js < index2; js++) {
2909 const AliTPCclusterMI *kcl = kr2[js];
2910 if (kcl->IsUsed(10)) continue;
2911 //
2912 //calcutate parameters
2913 //
2914 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2915 // stright track
2916 if (TMath::Abs(yy0)<0.000001) continue;
2917 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2918 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2919 Double_t r02 = (0.25+y0*y0)*dvertex2;
2920 //curvature (radius) cut
2921 if (r02<r2min) continue;
2922
2923 nin0++;
2924 //
2925 Double_t c0 = 1/TMath::Sqrt(r02);
2926 if (yy0>0) c0*=-1.;
2927
2928
2929 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2930 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
b67e07dc 2931 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2932 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
91162307 2933 //
2934 //
2935 Double_t z0 = kcl->GetZ();
2936 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2937 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2938 nin1++;
2939 //
2940 Double_t dip = (z1-z0)*c0/dfi1;
2941 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2942 //
2943 y2 = kcl->GetY();
2944 if (dsec==0){
2945 x2 = xx2;
2946 z2 = kcl->GetZ();
2947 }
2948 else
2949 {
2950 // rotation
2951 z2 = kcl->GetZ();
2952 x2= xx2*cs-y2*sn*dsec;
2953 y2=+xx2*sn*dsec+y2*cs;
2954 }
2955
2956 x[0] = y1;
2957 x[1] = z1;
2958 x[2] = x0;
2959 x[3] = dip;
2960 x[4] = c0;
2961 //
2962 //
2963 // do we have cluster at the middle ?
2964 Double_t ym,zm;
2965 GetProlongation(x1,xm,x,ym,zm);
2966 UInt_t dummy;
2967 AliTPCclusterMI * cm=0;
2968 if (TMath::Abs(ym)-ymaxm<0){
2969 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2970 if ((!cm) || (cm->IsUsed(10))) {
2971 continue;
2972 }
2973 }
2974 else{
2975 // rotate y1 to system 0
2976 // get state vector in rotated system
2977 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2978 Double_t xr2 = x0*cs+yr1*sn*dsec;
2979 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2980 //
2981 GetProlongation(xx2,xm,xr,ym,zm);
2982 if (TMath::Abs(ym)-ymaxm<0){
2983 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2984 if ((!cm) || (cm->IsUsed(10))) {
2985 continue;
2986 }
2987 }
2988 }
2989
2990
2991 Double_t dym = 0;
2992 Double_t dzm = 0;
2993 if (cm){
2994 dym = ym - cm->GetY();
2995 dzm = zm - cm->GetZ();
2996 }
2997 nin2++;
2998
2999
3000 //
3001 //
3002 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3003 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3004 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3005 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3006 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3007
b67e07dc 3008 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3009 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3010 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3011 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3012 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3013 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 3014
b67e07dc 3015 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3016 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3017 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3018 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
91162307 3019
1c53abe2 3020 c[0]=sy1;
3021 c[1]=0.; c[2]=sz1;
3022 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3023 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3024 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3025 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3026 c[13]=f30*sy1*f40+f32*sy2*f42;
3027 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
91162307 3028
3029 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3030
1c53abe2 3031 UInt_t index=kr1.GetIndex(is);
6c94f330 3032 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
91162307 3033
1c53abe2 3034 track->fIsSeeding = kTRUE;
91162307 3035 track->fSeed1 = i1;
3036 track->fSeed2 = i2;
3037 track->fSeedType=3;
c9427e08 3038
91162307 3039
3040 //if (dsec==0) {
3041 FollowProlongation(*track, (i1+i2)/2,1);
3042 Int_t foundable,found,shared;
3043 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3044 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3045 seed->Reset();
3046 seed->~AliTPCseed();
3047 continue;
3048 }
3049 //}
3050
3051 nin++;
3052 FollowProlongation(*track, i2,1);
3053
3054
3055 //Int_t rc = 1;
3056 track->fBConstrain =1;
3057 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3058 track->fLastPoint = i1; // first cluster in track position
3059 track->fFirstPoint = track->fLastPoint;
3060
3061 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3062 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3063 track->fNShared>0.4*track->GetNumberOfClusters() ) {
3064 seed->Reset();
3065 seed->~AliTPCseed();
c9427e08 3066 continue;
3067 }
91162307 3068 nout1++;
3069 // Z VERTEX CONDITION
c274e255 3070 Double_t zv, bz=GetBz();
3071 if ( !track->GetZAt(0.,bz,zv) ) continue;
91162307 3072 if (TMath::Abs(zv-z3)>cuts[2]) {
3073 FollowProlongation(*track, TMath::Max(i2-20,0));
c274e255 3074 if ( !track->GetZAt(0.,bz,zv) ) continue;
91162307 3075 if (TMath::Abs(zv-z3)>cuts[2]){
3076 FollowProlongation(*track, TMath::Max(i2-40,0));
c274e255 3077 if ( !track->GetZAt(0.,bz,zv) ) continue;
91162307 3078 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
3079 // make seed without constrain
3080 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3081 FollowProlongation(*track2, i2,1);
3082 track2->fBConstrain = kFALSE;
3083 track2->fSeedType = 1;
3084 arr->AddLast(track2);
3085 seed->Reset();
3086 seed->~AliTPCseed();
3087 continue;
3088 }
3089 else{
3090 seed->Reset();
3091 seed->~AliTPCseed();
3092 continue;
3093
3094 }
3095 }
c9427e08 3096 }
1c53abe2 3097
91162307 3098 track->fSeedType =0;
3099 arr->AddLast(track);
3100 seed = new AliTPCseed;
3101 nout2++;
3102 // don't consider other combinations
3103 if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
3104 break;
1c53abe2 3105 }
3106 }
3107 }
6bdc18d6 3108 if (fDebug>3){
3109 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
91162307 3110 }
3111 delete seed;
1c53abe2 3112}
3113
1627d1c4 3114
91162307 3115void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3116 Float_t deltay) {
3117
3118
3119
1627d1c4 3120 //-----------------------------------------------------------------
91162307 3121 // This function creates track seeds.
1627d1c4 3122 //-----------------------------------------------------------------
91162307 3123 // cuts[0] - fP4 cut
3124 // cuts[1] - tan(phi) cut
3125 // cuts[2] - zvertex cut
3126 // cuts[3] - fP3 cut
3127
3128
3129 Int_t nin0 = 0;
3130 Int_t nin1 = 0;
3131 Int_t nin2 = 0;
3132 Int_t nin = 0;
3133 Int_t nout1 = 0;
3134 Int_t nout2 = 0;
3135 Int_t nout3 =0;
3136 Double_t x[5], c[15];
3137 //
3138 // make temporary seed
3139 AliTPCseed * seed = new AliTPCseed;
1627d1c4 3140 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3141 // Double_t cs=cos(alpha), sn=sin(alpha);
91162307 3142 //
3143 //
1627d1c4 3144
91162307 3145 // first 3 padrows
3146 Double_t x1 = GetXrow(i1-1);
3147 const AliTPCRow& kr1=GetRow(sec,i1-1);
3148 Double_t y1max = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
3149 //
3150 Double_t x1p = GetXrow(i1);
3151 const AliTPCRow& kr1p=GetRow(sec,i1);
3152 //
3153 Double_t x1m = GetXrow(i1-2);
3154 const AliTPCRow& kr1m=GetRow(sec,i1-2);
1627d1c4 3155
91162307 3156 //
3157 //last 3 padrow for seeding
3158 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3159 Double_t x3 = GetXrow(i1-7);
3160 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3161 //
3162 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3163 Double_t x3p = GetXrow(i1-6);
3164 //
3165 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3166 Double_t x3m = GetXrow(i1-8);
1627d1c4 3167
91162307 3168 //
3169 //
3170 // middle padrow
3171 Int_t im = i1-4; //middle pad row index
3172 Double_t xm = GetXrow(im); // radius of middle pad-row
3173 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
3174 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3175 //
3176 //
3177 Double_t deltax = x1-x3;
3178 Double_t dymax = deltax*cuts[1];
3179 Double_t dzmax = deltax*cuts[3];
3180 //
3181 // loop over clusters
3182 for (Int_t is=0; is < kr1; is++) {
1627d1c4 3183 //
91162307 3184 if (kr1[is]->IsUsed(10)) continue;
3185 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1627d1c4 3186 //
91162307 3187 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3188 //
3189 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3190 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3191 //
3192 Double_t y3, z3;
1627d1c4 3193 //
1627d1c4 3194 //
91162307 3195 UInt_t index;
3196 for (Int_t js=index1; js < index2; js++) {
3197 const AliTPCclusterMI *kcl = kr3[js];
3198 if (kcl->IsUsed(10)) continue;
3199 y3 = kcl->GetY();
3200 // apply angular cuts
3201 if (TMath::Abs(y1-y3)>dymax) continue;
3202 x3 = x3;
3203 z3 = kcl->GetZ();
3204 if (TMath::Abs(z1-z3)>dzmax) continue;
3205 //
3206 Double_t angley = (y1-y3)/(x1-x3);
3207 Double_t anglez = (z1-z3)/(x1-x3);
3208 //
3209 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3210 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3211 //
3212 Double_t yyym = angley*(xm-x1)+y1;
3213 Double_t zzzm = anglez*(xm-x1)+z1;
3214
3215 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3216 if (!kcm) continue;
3217 if (kcm->IsUsed(10)) continue;
3218
3219 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3220 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3221 //
3222 //
3223 //
3224 Int_t used =0;
3225 Int_t found =0;
3226 //
3227 // look around first
3228 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3229 anglez*(x1m-x1)+z1,
3230 erry,errz,index);
3231 //
3232 if (kc1m){
3233 found++;
3234 if (kc1m->IsUsed(10)) used++;
1627d1c4 3235 }
91162307 3236 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3237 anglez*(x1p-x1)+z1,
3238 erry,errz,index);
1627d1c4 3239 //
91162307 3240 if (kc1p){
3241 found++;
3242 if (kc1p->IsUsed(10)) used++;
1627d1c4 3243 }
91162307 3244 if (used>1) continue;
3245 if (found<1) continue;
1627d1c4 3246
91162307 3247 //
3248 // look around last
3249 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3250 anglez*(x3m-x3)+z3,
3251 erry,errz,index);
3252 //
3253 if (kc3m){
3254 found++;
3255 if (kc3m->IsUsed(10)) used++;
3256 }
3257 else
3258 continue;
3259 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3260 anglez*(x3p-x3)+z3,
3261 erry,errz,index);
3262 //
3263 if (kc3p){
3264 found++;
3265 if (kc3p->IsUsed(10)) used++;
3266 }
3267 else
3268 continue;
3269 if (used>1) continue;
3270 if (found<3) continue;
3271 //
3272 Double_t x2,y2,z2;
3273 x2 = xm;
3274 y2 = kcm->GetY();
3275 z2 = kcm->GetZ();
3276 //
3277
1627d1c4 3278 x[0]=y1;
3279 x[1]=z1;
b67e07dc 3280 x[4]=F1(x1,y1,x2,y2,x3,y3);
91162307 3281 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3282 nin0++;
3283 //
b67e07dc 3284 x[2]=F2(x1,y1,x2,y2,x3,y3);
91162307 3285 nin1++;
3286 //
b67e07dc 3287 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
91162307 3288 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3289 nin2++;
3290 //
3291 //
3292 Double_t sy1=0.1, sz1=0.1;
3293 Double_t sy2=0.1, sz2=0.1;
3294 Double_t sy3=0.1, sy=0.1, sz=0.1;
1627d1c4 3295
b67e07dc 3296 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3297 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3298 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3299 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3300 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3301 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 3302
b67e07dc 3303 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3304 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3305 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3306 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1627d1c4 3307
3308 c[0]=sy1;
91162307 3309 c[1]=0.; c[2]=sz1;
1627d1c4 3310 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3311 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3312 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3313 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3314 c[13]=f30*sy1*f40+f32*sy2*f42;
3315 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3316
91162307 3317 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3318
3319 UInt_t index=kr1.GetIndex(is);
6c94f330 3320 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
91162307 3321
3322 track->fIsSeeding = kTRUE;
3323
3324 nin++;
3325 FollowProlongation(*track, i1-7,1);
3326 if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
3327 track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3328 seed->Reset();
3329 seed->~AliTPCseed();
3330 continue;
3331 }
3332 nout1++;
3333 nout2++;
3334 //Int_t rc = 1;
3335 FollowProlongation(*track, i2,1);
3336 track->fBConstrain =0;
3337 track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3338 track->fFirstPoint = track->fLastPoint;
3339
3340 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3341 track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
3342 track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3343 seed->Reset();
3344 seed->~AliTPCseed();
3345 continue;
3346 }
3347
3348 {
3349 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3350 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3351 FollowProlongation(*track2, i2,1);
3352 track2->fBConstrain = kFALSE;
3353 track2->fSeedType = 4;
3354 arr->AddLast(track2);
3355 seed->Reset();
3356 seed->~AliTPCseed();
3357 }
3358
3359
3360 //arr->AddLast(track);
3361 //seed = new AliTPCseed;
3362 nout3++;
3363 }
3364 }
3365
6bdc18d6 3366 if (fDebug>3){
3367 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
91162307 3368 }
3369 delete seed;
3370}
3371
3372
3373//_____________________________________________________________________________
176aff27 3374void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3375 Float_t deltay, Bool_t /*bconstrain*/) {
91162307 3376 //-----------------------------------------------------------------
3377 // This function creates track seeds - without vertex constraint
3378 //-----------------------------------------------------------------
3379 // cuts[0] - fP4 cut - not applied
3380 // cuts[1] - tan(phi) cut
3381 // cuts[2] - zvertex cut - not applied
3382 // cuts[3] - fP3 cut
3383 Int_t nin0=0;
3384 Int_t nin1=0;
3385 Int_t nin2=0;
3386 Int_t nin3=0;
3387 // Int_t nin4=0;
3388 //Int_t nin5=0;
3389
3390
3391
3392 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3393 // Double_t cs=cos(alpha), sn=sin(alpha);
3394 Int_t row0 = (i1+i2)/2;
3395 Int_t drow = (i1-i2)/2;
3396 const AliTPCRow& kr0=fSectors[sec][row0];
3397 AliTPCRow * kr=0;
3398
3399 AliTPCpolyTrack polytrack;
3400 Int_t nclusters=fSectors[sec][row0];
3401 AliTPCseed * seed = new AliTPCseed;
3402
3403 Int_t sumused=0;
3404 Int_t cused=0;
3405 Int_t cnused=0;
3406 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3407 Int_t nfound =0;
3408 Int_t nfoundable =0;
3409 for (Int_t iter =1; iter<2; iter++){ //iterations
3410 const AliTPCRow& krm=fSectors[sec][row0-iter];
3411 const AliTPCRow& krp=fSectors[sec][row0+iter];
3412 const AliTPCclusterMI * cl= kr0[is];
3413
3414 if (cl->IsUsed(10)) {