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